Ubah urutan matriks dan pertahankan tata letak memori kolom-utama

Dec 10 2020

Ilustrasi masalah: norma baris dari sebuah matriks

Pertimbangkan contoh mainan ini di mana saya menghitung norma semua kolom dari matriks acak 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)

Kemudian hitung norma baris

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)

Faktor antara dua eksekusi ini karena (menurut saya) tata letak memori dari matriks (kolom-utama). Memang perhitungan norma baris adalah loop pada data yang tidak bersebelahan, yang mengarah ke kode non-vektorisasi dengan cache miss. Saya ingin waktu eksekusi yang sama untuk kedua perhitungan norma.

Apakah mungkin untuk mengubah tata letak Mmenjadi baris mayor untuk mendapatkan kecepatan yang sama saat menghitung norma baris?

Apa yang saya coba

Saya mencoba dengan transposedan permutedimstanpa hasil, tampaknya ketika menggunakan fungsi-fungsi ini memori sekarang berada di baris-mayor (jadi kolom utama dari matriks asli).

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)

Saya digunakan di copysini untuk mencoba memaksa tata letak baru.

Bagaimana cara memaksa tata letak utama kolom dari transposisi, atau tata letak baris-utama dari matriks asli?

EDIT

Seperti yang ditunjukkan oleh @mcabbott dan @ przemyslaw-szufel, ada kesalahan dalam kode terakhir yang saya tunjukkan, saya menghitung norma baris Mtalih-alih norma kolom.

Ujilah norma kolom Mtmemberi sebagai gantinya:

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)

Jadi pada akhirnya tampaknya permutedimstoko - toko di kolom-mayor, seperti yang diharapkan. Faktanya, array Julia selalu disimpan di kolom-mayor. transposeadalah semacam pengecualian karena merupakan baris utama viewdari matriks tersimpan utama kolom.

Jawaban

3 PrzemyslawSzufel Dec 11 2020 at 02:39

Ada beberapa masalah di sini:

  • Anda salah mengukur kode Anda - kemungkinan besar menguji kode yang dikompilasi pada proses pertama dan kode yang tidak dikompilasi (dan karenanya mengukur waktu kompilasi) pada proses kedua. Anda harus selalu menjalankan @timedua kali atau menggunakan BenchmarkTools sebagai gantinya
  • kode Anda tidak efisien - melakukan penyalinan memori yang tidak perlu
  • tipe M tidak stabil dan oleh karena itu pengukuran mencakup waktu untuk mengetahui jenisnya yang tidak terjadi ketika Anda biasanya menjalankan fungsi Julia
  • Anda tidak perlu memiliki lambda - Anda bisa mengurai fungsi secara langsung.
  • seperti yang disebutkan oleh @mcabbott kode Anda mengandung bug dan Anda mengukur dua kali hal yang sama. Setelah membersihkan kode Anda terlihat seperti ini:
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)

Sekarang pertanyaannya tentang tata letak data. Julia menggunakan tata letak memori kolom-utama. Karenanya operasi yang bekerja pada kolom akan lebih cepat daripada yang bekerja pada baris. Salah satu solusi yang mungkin adalah memiliki salinan yang diubah urutannya dari M:

const Mᵀ = collect(M')

Ini membutuhkan beberapa waktu untuk menyalin tetapi memungkinkan Anda nanti untuk menyesuaikan kinerja:

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

Anda membuang banyak waktu untuk membuat salinan setiap kolom / baris saat menghitung norma. Gunakan views sebagai gantinya, atau lebih baik lagi, eachcol/ eachrow, yang juga tidak mengalokasikan:

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)