Julia中3D数组与矩阵向量存储方式下矩阵乘向量的性能差异：为何向量存储版本更快？

Great question! The performance difference you're seeing is completely normal, and it boils down to how Julia handles memory layout and array access patterns. Let's break this down step by step:

1. Memory Contiguity is the Biggest Factor #

Julia uses column-major (Fortran-style) memory layout for arrays. This means elements are stored column-by-column, then row-by-row, then higher dimensions.

• For your 3D array X (shape nid × npar × 3), memory is ordered first by nid, then npar, then the third dimension. When you slice X[:,:,i], you're accessing non-consecutive chunks of memory—each column of the sliced matrix is separated by nid × npar elements in memory. This leads to terrible cache utilization, since the CPU can't efficiently prefetch data for the matrix multiplication.
• For your vector of matrices X1, each individual matrix is a contiguous 2D array. When you run X[i] * beta, the CPU can load entire chunks of the matrix into cache at once, letting the BLAS-backed matrix-vector multiplication run at full speed.

2. Slicing 3D Arrays Introduces Hidden Overhead #

Even though Julia creates views (not copies) for slices like X[:,:,i], the non-contiguous memory access still adds overhead. BLAS routines (which power Julia's core linear algebra operations) are hyper-optimized for contiguous memory, so they perform best when working with uninterrupted blocks of data.

Can We Optimize the 3D Array Version? #

Absolutely! You can rewrite the 3D array function to avoid non-contiguous slicing and leverage contiguous memory operations. Here's an optimized version:

function f_optimized(X::Array{Int,3}, beta::Vector{Float64})::Array{Float64}
    # Reshape beta to broadcast with the second dimension of X
    beta_broadcast = reshape(beta, 1, :, 1)
    # Element-wise multiplication + sum over the parameter dimension
    intermediate = sum(X .* beta_broadcast, dims=2)
    # Reshape to match the output format of f/g
    return hcat(intermediate[:, 1, :]...)
end

This version uses broadcasting and summation over the npar dimension, which operates entirely on contiguous memory. When you benchmark this, you'll see performance close to (or even matching) the vector-of-matrices version, since it eliminates the non-contiguous slicing overhead.

Summary #

Your initial results are totally expected: the vector-of-matrices approach uses contiguous memory blocks that play perfectly with Julia's (and BLAS's) optimization strategies. The 3D array version is slower because of non-contiguous memory access, but with a few tweaks to how you structure the computation, you can get comparable performance.

内容的提问来源于stack exchange，提问作者djourd1