# Scalability of factor-solve vs. pseudoinverse + product

I’ve always read advice and warnings about the poor scalability and time spent trying to invert a matrix and why it’s better to solve a system of equations whenever the inverted matrix will be used with a product of something, e.g. (A^-1)b vs. solving Ax=b.

However, I’ve been playing with matrix factorization through alternating least squares and it seems to run reasonably faster at different problem sizes if I solve it by using inverses than by solving systems of linear equations (which is how this is implemented in most software), which is what everyone says shouldn’t happen.

For example, let's say I generate a random matrix X and factorize it as the product of two lower dimensional matrices A and B (initialized randomly), so as to minimize sum( (X-AB)^2 ):

import numpy as np

m1=10000
m2=5000
k=300

X=np.random.uniform(low=1,high=5,size=(m1,m2))
A=np.random.normal(size=(m1,k))
B=np.random.normal(size=(k,m2))


Then I solve this problem using factor-solve:

%%timeit
for i in range(10):
A[:,:]=np.linalg.solve(B.dot(B.T),B.dot(X.T)).T
B[:,:]=np.linalg.solve(A.T.dot(A), A.T.dot(X))

1 loop, best of 3: 37.6 s per loop


And solve the same problem by inverting the first term:

%%timeit
for i in range(10):
A[:,:]=(np.linalg.pinv(B.dot(B.T)).dot(B.dot(X.T))).T
B[:,:]=np.linalg.pinv(A.T.dot(A)).dot(A.T.dot(X))
1 loop, best of 3: 36.2 s per loop


Am I misunderstanding something? Why is the second option faster than the first one? I think this problem should be large enough for the algorithmic complexity to play a big role on running times, but it seems to bring no advantage.