Транспонировать матрицу и сохранить макет памяти по основным столбцам
Иллюстрация проблемы: строковые нормы матрицы
Рассмотрим этот игрушечный пример, в котором я вычисляю нормы всех столбцов случайной матрицы 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
хранимой матрицы основного столбца.
Ответы
Здесь есть несколько проблем:
- вы неправильно используете свой код - скорее всего, вы тестируете скомпилированный код при первом запуске и некомпилированный код (и, следовательно, измеряете время компиляции) во втором. Всегда следует запускать
@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
Вы тратите много времени на создание копий каждого столбца / строки при расчете норм. 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)