381 - A18 - Matrix Multiplication Six Ways

340 days ago by Professor381

John Travis

Mississippi College

In this worksheet, the various ways to create the product of two matrices $A$ and $B$ to get $C = AB$ are demonstrated.  Notice that the major change in each is to simply cycle through the 6 different ways that three loops can be arranged.

The first two approaches create $C1$ and $C2$ and generally follow the standard approach taught for matrix multiplication.

The second two approaches create $C3$ and $C4$ and follow the linear combination of columns (or rows) approach for matrix multiplication.  This is similar to how Lay's Linear Algebra first introduces matrix-vector and then matrix-matrix multiplication. This approach might be useful if one had matrices stored in memory by column/column or by row/row.

The final two approaches create $C5$ and $C6$ which essentially do all of the accumulation parts of the standard approach but only one term at a time.  That is, normally you'd have something like

$c_{i,j} = \Sigma_{k=1}^n a_{i,k}b_{k,j}$

$ = a_{i,1}b_{1,j} +a_{i,2}b_{2,j} + ... + a_{i,n}b_{n,j}$

and one might do that entire sum in one fell swoop and then procees to the next i-j entry.  Instead, a "rank-one" update will do only one of these terms in the sum but for every entry in the product matric C.  That is, for a fixed k, only do $a_{i,k}b_{k,j}$ for EVERY i-j.

In pondering each of these approaches, one might take into consideration how data is stored in memory and how each approach accesses that data.  Indeed, a matrix is likely stored in memory as one long string of numbers perhaps by one row followed by another row and followed again by others.  Or perhaps stored by one column and then another and followed again by others.  It may be more efficient to only fetch rows of both A and B or only fetch columns of A and B rather than rows of A and columns of B. 

Also, with the use of parallel processing, it might be more efficient to assign tasks to various processors using one of these approaches versions versus another.

# Let's set up some initial matrices to play with... A = [ 1 2 6 -9; -4 3 2 7; -5 0 1 -1; 2 8 -3 5] B = [1 2 3 4; 5 6 7 8; -1 2 3 -4; -3 1 6 2] C = A*B m = size(A)(1); n = size(B)(2); Bnumrows = size(B)(1); Anumcols = size(A)(2); if Anumcols != Bnumrows disp("Dimensions don't agree. Can't multiply AB") else disp("Dimensions agree. Proceed to multiplication") end 
       
A =

 1 2 6 -9
 -4 3 2 7
 -5 0 1 -1
 2 8 -3 5

B =

 1 2 3 4
 5 6 7 8
 -1 2 3 -4
 -3 1 6 2

C =

 32 17 -19 -22
 -12 21 57 14
 -3 -9 -18 -26
 30 51 83 94

Dimensions agree. Proceed to multiplication
A =

 1 2 6 -9
 -4 3 2 7
 -5 0 1 -1
 2 8 -3 5

B =

 1 2 3 4
 5 6 7 8
 -1 2 3 -4
 -3 1 6 2

C =

 32 17 -19 -22
 -12 21 57 14
 -3 -9 -18 -26
 30 51 83 94

Dimensions agree. Proceed to multiplication
disp("Do the multiplication ourselves the standard way by combining A's row by B's column") C1 = zeros(m,m); for row = 1:m for column = 1:n for k = 1:Anumcols # Create C by accessing A by row and B by column C1(row,column) += A(row,k)*B(k,column); end end end C1 disp("Taking C - this approach gives") C-C1 
       
Do the multiplication ourselves the standard way by combining A's row by
B's column
C1 =

 32 17 -19 -22 0
 -12 21 57 14 0
 -3 -9 -18 -26 0
 2 11 19 10 0
 30 51 83 94 0

Taking C - this approach gives
error:
/home/sageserver/.sage/temp/travis-PowerEdge-T430/317822/interface/tmp31\
7829: operator -: nonconformant arguments (op1 is 5x4, op2 is 5x5)
error: called from
   
/home/sageserver/.sage/temp/travis-PowerEdge-T430/317822/interface/tmp31\
7829 at line 14 column 2
Do the multiplication ourselves the standard way by combining A's row by B's column
C1 =

 32 17 -19 -22 0
 -12 21 57 14 0
 -3 -9 -18 -26 0
 2 11 19 10 0
 30 51 83 94 0

Taking C - this approach gives
error: /home/sageserver/.sage/temp/travis-PowerEdge-T430/317822/interface/tmp317829: operator -: nonconformant arguments (op1 is 5x4, op2 is 5x5)
error: called from
    /home/sageserver/.sage/temp/travis-PowerEdge-T430/317822/interface/tmp317829 at line 14 column 2
disp("Do the multiplication by combining B's column by A's row") C2 = zeros(m,m); for column = 1:n for row = 1:m for k = 1:Anumcols # Create C by accessing A by row and B by column C2(row,column) += A(row,k)*B(k,column); end end end C2 disp("Taking C - this approach gives") C-C2 
       
Do the multiplication by combining B's column by A's row
C2 =

 32 17 -19 -22
 -12 21 57 14
 -3 -9 -18 -26
 30 51 83 94

Taking C - this approach gives
ans =

 0 0 0 0
 0 0 0 0
 0 0 0 0
 0 0 0 0
Do the multiplication by combining B's column by A's row
C2 =

 32 17 -19 -22
 -12 21 57 14
 -3 -9 -18 -26
 30 51 83 94

Taking C - this approach gives
ans =

 0 0 0 0
 0 0 0 0
 0 0 0 0
 0 0 0 0
disp("Do the multiplication by accumulating all columns with a particular row") C3 = zeros(m,m); for row = 1:m for k = 1:Anumcols for column = 1:n # Update C using each row of B by one entry in A # That is both A and B are accessed by ROWS C3(row,column) += A(row,k)*B(k,column); end end end C3 disp("Taking C - this approach gives") C-C3 
       
Do the multiplication by accumulating all columns with a particular row
Jake!
C3 =

 32 17 -19 -22
 -12 21 57 14
 -3 -9 -18 -26
 30 51 83 94

Taking C - this approach gives
ans =

 0 0 0 0
 0 0 0 0
 0 0 0 0
 0 0 0 0
Do the multiplication by accumulating all columns with a particular row
Jake!
C3 =

 32 17 -19 -22
 -12 21 57 14
 -3 -9 -18 -26
 30 51 83 94

Taking C - this approach gives
ans =

 0 0 0 0
 0 0 0 0
 0 0 0 0
 0 0 0 0
disp("Do the multiplication by accumulating all rows with a particular column") C4 = zeros(m,m); for column = 1:n for k = 1:Anumcols for row = 1:m # Update C using each column of A by one entry in B # That is both A and B are accessed by COLUMNS C4(row,column) += A(row,k)*B(k,column); end end end C4 disp("Taking C - this approach gives") C-C4 
       
Do the multiplication by accumulating all rows with a particular column
C4 =

 32 17 -19 -22
 -12 21 57 14
 -3 -9 -18 -26
 30 51 83 94

Taking C - this approach gives
ans =

 0 0 0 0
 0 0 0 0
 0 0 0 0
 0 0 0 0
Do the multiplication by accumulating all rows with a particular column
C4 =

 32 17 -19 -22
 -12 21 57 14
 -3 -9 -18 -26
 30 51 83 94

Taking C - this approach gives
ans =

 0 0 0 0
 0 0 0 0
 0 0 0 0
 0 0 0 0
disp("Do the multiplication by 'rank-one updates' by doing only one of the multiplications for each C_ij ") disp("Proceed through C by going across columns from left to right") A B C5 = zeros(m,m); for k = 1:Anumcols for column = 1:n for row = 1:m # Update C by doing only one entry in the normal multiplication but for every term in C # This is what is known as a rank-one update C5(row,column) += A(row,k)*B(k,column); end end C5 end disp("Taking C - this approach gives") C-C5 
       
Do the multiplication by 'rank-one updates' by doing only one of the
multiplications for each C_ij 
Proceed through C by going across columns from left to right
A =

 1 2 6 -9
 -4 3 2 7
 -5 0 1 -1
 2 8 -3 5

B =

 1 2 3 4
 5 6 7 8
 -1 2 3 -4
 -3 1 6 2

C5 =

 1 2 3 4
 -4 -8 -12 -16
 -5 -10 -15 -20
 2 4 6 8

C5 =

 11 14 17 20
 11 10 9 8
 -5 -10 -15 -20
 42 52 62 72

C5 =

 5 26 35 -4
 9 14 15 0
 -6 -8 -12 -24
 45 46 53 84

C5 =

 32 17 -19 -22
 -12 21 57 14
 -3 -9 -18 -26
 30 51 83 94

Taking C - this approach gives
ans =

 0 0 0 0
 0 0 0 0
 0 0 0 0
 0 0 0 0
Do the multiplication by 'rank-one updates' by doing only one of the multiplications for each C_ij 
Proceed through C by going across columns from left to right
A =

 1 2 6 -9
 -4 3 2 7
 -5 0 1 -1
 2 8 -3 5

B =

 1 2 3 4
 5 6 7 8
 -1 2 3 -4
 -3 1 6 2

C5 =

 1 2 3 4
 -4 -8 -12 -16
 -5 -10 -15 -20
 2 4 6 8

C5 =

 11 14 17 20
 11 10 9 8
 -5 -10 -15 -20
 42 52 62 72

C5 =

 5 26 35 -4
 9 14 15 0
 -6 -8 -12 -24
 45 46 53 84

C5 =

 32 17 -19 -22
 -12 21 57 14
 -3 -9 -18 -26
 30 51 83 94

Taking C - this approach gives
ans =

 0 0 0 0
 0 0 0 0
 0 0 0 0
 0 0 0 0
disp("Do the multiplication by 'rank-one updates' by doing only one of the multiplications for each C_ij ") disp("Proceed through C by going across rows from top to bottom") C6 = zeros(m,m); for k = 1:Anumcols for row = 1:m for column = 1:n # Update C by doing only one entry in the normal multiplication but for every term in C # This is again what is known as a rank-one update C6(row,column) += A(row,k)*B(k,column); end end C6 end disp("Taking C - this approach gives") C-C6 
       
Do the multiplication by 'rank-one updates' by doing only one of the
multiplications for each C_ij 
Proceed through C by going across rows from top to bottom
C6 =

 1 2 3 4
 -4 -8 -12 -16
 -5 -10 -15 -20
 2 4 6 8

C6 =

 11 14 17 20
 11 10 9 8
 -5 -10 -15 -20
 42 52 62 72

C6 =

 5 26 35 -4
 9 14 15 0
 -6 -8 -12 -24
 45 46 53 84

C6 =

 32 17 -19 -22
 -12 21 57 14
 -3 -9 -18 -26
 30 51 83 94

Taking C - this approach gives
ans =

 0 0 0 0
 0 0 0 0
 0 0 0 0
 0 0 0 0
Do the multiplication by 'rank-one updates' by doing only one of the multiplications for each C_ij 
Proceed through C by going across rows from top to bottom
C6 =

 1 2 3 4
 -4 -8 -12 -16
 -5 -10 -15 -20
 2 4 6 8

C6 =

 11 14 17 20
 11 10 9 8
 -5 -10 -15 -20
 42 52 62 72

C6 =

 5 26 35 -4
 9 14 15 0
 -6 -8 -12 -24
 45 46 53 84

C6 =

 32 17 -19 -22
 -12 21 57 14
 -3 -9 -18 -26
 30 51 83 94

Taking C - this approach gives
ans =

 0 0 0 0
 0 0 0 0
 0 0 0 0
 0 0 0 0