Ubah urutan matriks dan pertahankan tata letak memori kolom-utama
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 M
menjadi baris mayor untuk mendapatkan kecepatan yang sama saat menghitung norma baris?
Apa yang saya coba
Saya mencoba dengan transpose
dan permutedims
tanpa 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 copy
sini 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 Mt
alih-alih norma kolom.
Ujilah norma kolom Mt
memberi 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 permutedims
toko - toko di kolom-mayor, seperti yang diharapkan. Faktanya, array Julia selalu disimpan di kolom-mayor. transpose
adalah semacam pengecualian karena merupakan baris utama view
dari matriks tersimpan utama kolom.
Jawaban
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
@time
dua 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
Anda membuang banyak waktu untuk membuat salinan setiap kolom / baris saat menghitung norma. Gunakan view
s 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)