4 ms·

> The 'C' way is very efficient for computing Ax which presumes vector x to be a Column vector. When you want to compute xA instead which presumes vector x to b

> The 'C' way is very efficient for computing Ax which presumes vector x to be a Column vector. When you want to compute xA instead which presumes vector x to be a row vector, you want column-major ordering for efficiency. You have this the wrong way around. The C row-major ordering is better for A^Tx, the Fortran column-major format is better for Ax. If you do Ax in row-major, you end up multiplying a row of A by x as your simd vector op. This then leaves you needing to horizontally reduce the result, which you avoid by using a column-major layout that instead accumulates multiple different results.

bruce343434 4y ago
wait, so a row major stored matrix should be multiplied by extracting the columns first? Or what? Is there some sort of animation that illustrates Ax and xA for both row major storage and column major storage?
- kergonath 4y ago
  > wait, so a row major stored matrix should be multiplied by extracting the columns first? Or what? No, the other way around. In row-major storage, all the elements of the same row are contiguous, so it is faster to process the array row by row. > Is there some sort of animation that illustrates Ax and xA for both row major storage and column major storage? The explanations here are decent enough. All you need is the figure at the top, really: https://en.m.wikipedia.org/wiki/Row-_and_column-major_order https://en.m.wikipedia.org/wiki/Row-_and_column-major_order .
- adrian_b 4y ago
  The schoolbook formulae for most linear algebra operations are written using scalar products, e.g. for Ax as scalar products of rows of A with x. On a computer, such operations should almost never be implemented with scalar products, because the speed of a scalar product is limited by the latency of the FMA operation, not by its throughput. It is possible to accelerate the computation of a scalar product by ordering the elementary operations in a tree, but that complicates the program and it still does not reach the maximum speed of the hardware. Except for vector-vector operations (i.e. for level 2 BLAS or greater levels), it is possible to change the order of the loops so that the scalar products are replaced by AXPY operations (whose result vectors should be kept in registers, to not be limited by the memory transfer throughput), which are not limited by the latency of the FMA, like the scalar products. Except for vector-vector operations and matrix-vector operations (i.e. for level 3 BLAS or greater levels), it is possible to change the order of the loops so that the scalar products are replaced by rank-one matrix updates. Most of the computation of a rank-one matrix update consists of the tensor product of 2 vectors. Therefore, a matrix-matrix multiplication can be computed either by scalar products of the rows of the 1st matrix with the columns of the 2nd matrix, or by the tensor products of the columns of the 1st matrix with the rows of the 2nd matrix. The second method is much faster. The result of the tensor product must be kept in registers for the duration of the computation, so large matrices are partitioned in small blocks, i.e. submatrices, which are multiplied directly (an additional complication that increases the number of nested loops is that the small blocks that can be multiplied directly must be grouped in larger blocks that can be kept in the level-2 cache memory and reused for computations). For matrix-matrix multiplications it is less important whether they are stored in column major order or in row major order, because when submatrices are brought into the L2 cache, the matrix elements will be gathered or scattered to be in the best order, before being loaded into registers repeatedly. For matrix-vector multiplications, it is better to have the matrix in column major order, in order to compute the product by AXPY operations between columns and elements of the vector, and not by scalar products between rows and the vector (computing AXPY operations with shorter parts of a column at a time, so that the result vector can be kept in registers).
  - rocqua 4y ago
    Thanks for this excellent explanation. I had questions, but they were answered by a closer reading of your explanation. I recently did some work where, interestingly enough, this did not apply. I was doing linear algebra modulo 2 (i.e. binary entries, adding is xor, multiplication is and). In order to have any kind of memory (throughput) efficiency, you want to pack your matrix entries into bits of integers so a single 64 int can represent 64 entries. This packing (when done row-wise or column-wise) gives a very quick primitive for computing inner products. <a, b> is just (in pseudo-code) sum(popcnt(a[i] & b[i]) for i in range len(a)) % 2 Hence in matrix multiplication over binary matrix you very much do want the l.h.s. to be row-major and the r.h.s. to be column major. And, in general, you very much want to read and write matrices in a contiguous manner. Because otherwise, depending on alignment, you are accessing 8, 32, or 64 bits of memory per single bit entry of matrix your are reading or writing. I wonder if a similar very fast scalar product primitive for floats (or ints?) could work. If not, I imagine that is due to FMA having some inherent efficiencies that cannot be gained by a very quick accumulation step. One of the beauties in the binary case is the popcnt instruction doing a 64-input accumulation very quickly.
rocqua 4y ago
edit: My argument was based on the naive implementation of matrix-vector multiplication. But, as explained here [1] there are more efficient algorithms for small-ish matrices that flip this around. For bigger matrices it remains better to do the efficient algorithm block-wise. There is an argument that, for 'obviousness', the naive method matter more than the efficient one. But that gets tenuous. Especially because with some exceptions (finite fields, ugh) people are probably much better of using libraries for matrix operations, rather than writing their own. [1] https://news.ycombinator.com/item?id=33359654 https://news.ycombinator.com/item?id=33359654

bruce343434 4y ago

wait, so a row major stored matrix should be multiplied by extracting the columns first? Or what? Is there some sort of animation that illustrates Ax and xA for both row major storage and column major storage?

kergonath 4y ago
> wait, so a row major stored matrix should be multiplied by extracting the columns first? Or what? No, the other way around. In row-major storage, all the elements of the same row are contiguous, so it is faster to process the array row by row. > Is there some sort of animation that illustrates Ax and xA for both row major storage and column major storage? The explanations here are decent enough. All you need is the figure at the top, really: https://en.m.wikipedia.org/wiki/Row-_and_column-major_order https://en.m.wikipedia.org/wiki/Row-_and_column-major_order .
adrian_b 4y ago
The schoolbook formulae for most linear algebra operations are written using scalar products, e.g. for Ax as scalar products of rows of A with x. On a computer, such operations should almost never be implemented with scalar products, because the speed of a scalar product is limited by the latency of the FMA operation, not by its throughput. It is possible to accelerate the computation of a scalar product by ordering the elementary operations in a tree, but that complicates the program and it still does not reach the maximum speed of the hardware. Except for vector-vector operations (i.e. for level 2 BLAS or greater levels), it is possible to change the order of the loops so that the scalar products are replaced by AXPY operations (whose result vectors should be kept in registers, to not be limited by the memory transfer throughput), which are not limited by the latency of the FMA, like the scalar products. Except for vector-vector operations and matrix-vector operations (i.e. for level 3 BLAS or greater levels), it is possible to change the order of the loops so that the scalar products are replaced by rank-one matrix updates. Most of the computation of a rank-one matrix update consists of the tensor product of 2 vectors. Therefore, a matrix-matrix multiplication can be computed either by scalar products of the rows of the 1st matrix with the columns of the 2nd matrix, or by the tensor products of the columns of the 1st matrix with the rows of the 2nd matrix. The second method is much faster. The result of the tensor product must be kept in registers for the duration of the computation, so large matrices are partitioned in small blocks, i.e. submatrices, which are multiplied directly (an additional complication that increases the number of nested loops is that the small blocks that can be multiplied directly must be grouped in larger blocks that can be kept in the level-2 cache memory and reused for computations). For matrix-matrix multiplications it is less important whether they are stored in column major order or in row major order, because when submatrices are brought into the L2 cache, the matrix elements will be gathered or scattered to be in the best order, before being loaded into registers repeatedly. For matrix-vector multiplications, it is better to have the matrix in column major order, in order to compute the product by AXPY operations between columns and elements of the vector, and not by scalar products between rows and the vector (computing AXPY operations with shorter parts of a column at a time, so that the result vector can be kept in registers).
- rocqua 4y ago
  Thanks for this excellent explanation. I had questions, but they were answered by a closer reading of your explanation. I recently did some work where, interestingly enough, this did not apply. I was doing linear algebra modulo 2 (i.e. binary entries, adding is xor, multiplication is and). In order to have any kind of memory (throughput) efficiency, you want to pack your matrix entries into bits of integers so a single 64 int can represent 64 entries. This packing (when done row-wise or column-wise) gives a very quick primitive for computing inner products. <a, b> is just (in pseudo-code) sum(popcnt(a[i] & b[i]) for i in range len(a)) % 2 Hence in matrix multiplication over binary matrix you very much do want the l.h.s. to be row-major and the r.h.s. to be column major. And, in general, you very much want to read and write matrices in a contiguous manner. Because otherwise, depending on alignment, you are accessing 8, 32, or 64 bits of memory per single bit entry of matrix your are reading or writing. I wonder if a similar very fast scalar product primitive for floats (or ints?) could work. If not, I imagine that is due to FMA having some inherent efficiencies that cannot be gained by a very quick accumulation step. One of the beauties in the binary case is the popcnt instruction doing a 64-input accumulation very quickly.

rocqua 4y ago

edit: My argument was based on the naive implementation of matrix-vector multiplication. But, as explained here [1] there are more efficient algorithms for small-ish matrices that flip this around. For bigger matrices it remains better to do the efficient algorithm block-wise. There is an argument that, for 'obviousness', the naive method matter more than the efficient one. But that gets tenuous. Especially because with some exceptions (finite fields, ugh) people are probably much better of using libraries for matrix operations, rather than writing their own. [1] https://news.ycombinator.com/item?id=33359654 https://news.ycombinator.com/item?id=33359654