# matrix multiplication speed

Consider the codes below:

dim = 500;

rd = RandomReal[1, {dim, dim}];

W = Table[{rd[[i]]}\[Transpose].{rd[[i]]}, {i, 1, dim}];

A = RandomReal[1, {dim, dim}];

x = Table[A.W[[i]], {i, 1, dim}]; // AbsoluteTiming

y = Table[A.{rd[[i]]}\[Transpose].{rd[[i]]}, {i, 1, dim}]; // AbsoluteTiming


W[[i]]is equal to {rd[[i]]}\[Transpose].{rd[[i]]}so basically the x multiplication is the same as y one. However y is around 3 times faster than x. Does anybody know why y is faster? I though because I had stored matrix W[[i]] then multiplication x would be faster, specially that y needs to compute W[[i]]first but x has it already.

It seems if we have a matrix in this form W = a.b in which a is a column matrix and b is a row matrix then W.a.b is faster than W.W:

In[1]:= dim = 2000;

In[2]:= a = RandomReal[1, {dim, 1}];

In[3]:= b = RandomReal[1, {1, dim}];

In[4]:= W = a.b;

In[5]:= W.(a.b);(*n^2+n^3 operations*)// AbsoluteTiming

Out[5]= {0.229023, Null}

In[6]:= W.W;(*n^3 operations*)// AbsoluteTiming

Out[6]= {0.183018, Null}

In[7]:= (W.a).b;(*2 n^2 operations*)// AbsoluteTiming

Out[7]= {0.021002, Null}


Does anybody know if it is possible to decompose every arbitrary matrix like this: W =a.b?

Edit

I added the number of multiplication operations for each calculation after reading the answer of bill s. Now, its clear why the last one is the fastest. So, the order by which we do the multiplication matters, at least sometimes.

Here is one way to see what is happening:

z1 = Table[A.({rd[[i]]}\[Transpose].{rd[[i]]}), {i, 1,  dim}]; // AbsoluteTiming
z2 = Table[(A.{rd[[i]]}\[Transpose]).{rd[[i]]}, {i, 1, dim}]; // AbsoluteTiming


which explicitly groups the calculations using parenthesis. As you found, z1 is slower than z2, which makes sense because z1 is essentially the product of two matrices, while z2 is the product of a matrix and a vector, followed by the product of a vector and vector.

For example, say you are multiplying as in z1. This takes n^3 multiplications where the matrices are size n (500 in the above example). On the other hand, A.rd takes n^2 multiplications and results in a vector. Then the outer product of the two vectors requires another n^2 multiplications, so z2 takes a total of 2 n^2 multiplications.

It is certainly not possible to reduce every matrix to a product of two vectors. One way to see this is that the rank of a product of two vectors is at most 1, so any W with rank greater than 1 cannot be decomposed this way.

• I didn't get. Isn't the number of arithmetic calculations the same? I guess there are more efficient algorithms for the product of a matrix and a vector than that of a matrix and a matrix. – MOON Nov 10 '14 at 17:02
• @yashar, Matrix-matrix and matrix-vertor multiplications have a sufficiently different number of operations. – ybeltukov Nov 10 '14 at 17:36
• @yashar Have a look at an example in David Wagner's book named "Matrix chain multiplication". He uses dynamic programming to optimize the multiplication strategy, and discusses this issue at length. – Leonid Shifrin Nov 10 '14 at 17:37
• @ybeltukov z1 is a matrix matrix product and z2 is a matrix vector and then a vector vector product. The number of arithmetic calculations must be the same. – MOON Nov 10 '14 at 19:39
• @LeonidShifrin the only book I found was this mathematica.stackexchange.com/questions/16485/… – MOON Nov 10 '14 at 19:50