Transponer la matriz y mantener el diseño de la memoria de la columna principal

Question

Ilustración del problema: las normas de fila de una matriz

Considere este ejemplo de juguete donde calculo las normas de todas las columnas de una matriz aleatoria 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)

Luego calcule las normas de las filas

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)

El factor entre las dos ejecuciones se debe (creo) al diseño de memoria de la matriz (columna principal). De hecho, el cálculo de las normas de fila es un bucle sobre datos no contiguos, lo que conduce a un código no vectorizado con falta de caché. Me gustaría tener el mismo tiempo de ejecución para los cálculos de ambas normas.

¿Es posible convertir el diseño de Muna fila mayor para obtener la misma velocidad al calcular las normas de las filas?

Que intenté

Lo intenté con transposey permutedimssin éxito, parece que al usar estas funciones, la memoria ahora está en la fila principal (por lo que las columnas principales de la matriz original).

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)

Usé copyaquí para intentar forzar el nuevo diseño.

¿Cómo puedo forzar el diseño de columna principal de la transposición o el diseño de fila principal de la matriz original?

EDITAR

Como lo señalaron @mcabbott y @ przemyslaw-szufel, hubo un error en el último código que mostré, calculé las normas de las filas de en Mtlugar de las normas de las columnas.

La prueba sobre las normas de las columnas de Mtdar en cambio:

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)

Así que al final parece que las permutedimstiendas en columna-mayor, como era de esperar. De hecho, las matrices de Julia siempre se almacenan en column-major. transposees una especie de excepción porque es una fila principal viewde una matriz almacenada de columna principal.

Answer 1

Hay varios problemas aquí:

• está probando incorrectamente su código; lo más probable es que pruebe el código compilado en la primera ejecución y el código no compilado (y, por lo tanto, mida los tiempos de compilación) en la segunda ejecución. Siempre debe ejecutar @timedos veces o usar BenchmarkTools en su lugar
• su código es ineficiente - ¿copia de memoria innecesaria
• el tipo de M es inestable y, por lo tanto, la medición incluye el tiempo para averiguar su tipo, lo cual no es un caso cuando normalmente se ejecuta una función de Julia
• no es necesario tener una lambda; simplemente puede analizar la función directamente.
• como lo menciona @mcabbott, su código contiene un error y está midiendo el doble de lo mismo. Después de limpiar su código se ve así:

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)

Ahora la pregunta sobre el diseño de datos. Julia está usando un diseño de memoria de columna principal. Por lo tanto, las operaciones que funcionan en columnas serán más rápidas que las que funcionan en filas. Una posible solución alternativa sería tener una copia transpuesta de M:

const Mᵀ = collect(M')

Esto requiere algo de tiempo para copiar, pero luego le permite igualar el rendimiento:

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

Answer 2

Está perdiendo mucho tiempo creando copias de cada columna / fila al calcular las normas. Use views en su lugar, o mejor aún, eachcol/ eachrow, que tampoco asigna:

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)