Транспонировать матрицу и сохранить макет памяти по основным столбцам

Dec 10 2020

Иллюстрация проблемы: строковые нормы матрицы

Рассмотрим этот игрушечный пример, в котором я вычисляю нормы всех столбцов случайной матрицы M

julia> M = rand(Float64, 10000, 10000);

julia> @time map(x -> norm(x), M[:,j] for j in 1:size(M)[2]);
  0.363795 seconds (166.70 k allocations: 770.086 MiB, 27.78% gc time)

Затем вычислите нормы строки

julia> @time map(x -> norm(x), M[:,i] for i in 1:size(M)[1]);
  1.288872 seconds (176.19 k allocations: 770.232 MiB, 0.37% gc time)

Фактор между двумя выполнениями связан (я думаю) с расположением памяти матрицы (по столбцам). Действительно, вычисление норм строк - это цикл для несмежных данных, который приводит к невекторизованному коду с промахами в кэше. Я хотел бы иметь одинаковое время выполнения для обоих вычислений норм.

Можно ли преобразовать раскладку из Mстрочки в большую, чтобы получить ту же скорость при расчете норм рядов?

Что я пробовал

Я пробовал с успехом transposeи permutedimsбезуспешно, кажется, что при использовании этих функций память теперь находится в строчном порядке (поэтому столбцы в исходной матрице).

julia> Mt = copy(transpose(M));

julia> @time map(x -> norm(x), Mt[j,:] for j in 1:size(M)[2]);
  1.581778 seconds (176.19 k allocations: 770.230 MiB)

julia> Mt = copy(permutedims(M,[2,1]));

julia> @time map(x -> norm(x), Mt[j,:] for j in 1:size(M)[2]);
  1.454153 seconds (176.19 k allocations: 770.236 MiB, 9.98% gc time)

Я использовал copyздесь, чтобы попытаться заставить новый макет.

Как я могу заставить макет транспонирования по основным столбцам или макет по строкам исходной матрицы?

РЕДАКТИРОВАТЬ

Как указали @mcabbott и @ przemyslaw-szufel, в последнем показанном мной коде была ошибка: я вычислил нормы строк Mtвместо норм столбцов.

Тест на норму столбцов Mtвыдаем взамен:

julia> Mt = transpose(M);
julia> @time map(x -> norm(x), M[:,j] for j in 1:size(M)[2]);
  1.307777 seconds (204.52 k allocations: 772.032 MiB, 0.45% gc time)

julia> Mt = permutedims(M)
julia> @time map(x -> norm(x), M[:,j] for j in 1:size(M)[2]); 
  0.334047 seconds (166.53 k allocations: 770.079 MiB, 1.42% gc time)

Так что, в конце концов, кажется, что permutedimsмагазины в столбцах, как и следовало ожидать. Фактически, массивы Julia всегда хранятся по столбцам. transposeявляется своего рода исключением, потому что это основная строка viewхранимой матрицы основного столбца.

Ответы

3 PrzemyslawSzufel Dec 11 2020 at 02:39

Здесь есть несколько проблем:

  • вы неправильно используете свой код - скорее всего, вы тестируете скомпилированный код при первом запуске и некомпилированный код (и, следовательно, измеряете время компиляции) во втором. Всегда следует запускать @timeдважды или использовать вместо этого BenchmarkTools
  • ваш код неэффективен - делает ненужное копирование памяти
  • тип M нестабилен, и, следовательно, измерение включает время, необходимое для определения его типа, что не является случаем, когда вы обычно выполняете функцию Julia
  • вам не нужно иметь лямбда - вы можете просто проанализировать функцию напрямую.
  • как упоминалось @mcabbott, ваш код содержит ошибку, и вы измеряете дважды одно и то же. После очистки ваш код выглядит так:
julia> using LinearAlgebra, BenchmarkTools

julia> const M = rand(10000, 10000);

julia> @btime map(norm, @view M[:,j] for j in 1:size(M)[2]);
  49.378 ms (2 allocations: 78.20 KiB)

julia> @btime map(norm, @view M[i, :] for i in 1:size(M)[1]);
  1.013 s (2 allocations: 78.20 KiB)

Теперь вопрос о расположении данных. Джулия использует макет памяти по основным столбцам. Следовательно, операции со столбцами будут выполняться быстрее, чем операции со строками. Одним из возможных решений может быть транспонированная копия M:

const Mᵀ = collect(M')

Это требует некоторого времени на копирование, но позволяет позже добиться соответствия производительности:

julia> @btime map(norm, @view Mᵀ[:,j] for j in 1:size(M)[2]);
  48.455 ms (2 allocations: 78.20 KiB)

julia> map(norm, Mᵀ[:,j] for j in 1:size(M)[2]) == map(norm, M[i,:] for i in 1:size(M)[1])
true
2 DNF Dec 11 2020 at 18:00

Вы тратите много времени на создание копий каждого столбца / строки при расчете норм. viewВместо этого используйте s или, еще лучше, eachcol/ eachrow, которые также не выделяют:

julia> M = rand(1000, 1000);

julia> @btime map(norm, $M[:,j] for j in 1:size($M, 2));  # slow and ugly
  946.301 μs (1001 allocations: 7.76 MiB)

julia> @btime map(norm, eachcol($M)); # fast and nice 223.199 μs (1 allocation: 7.94 KiB) julia> @btime norm.(eachcol($M));  # even nicer, but allocates more for some reason.
  227.701 μs (3 allocations: 47.08 KiB)