Numpys pinv seem slow?
I was playing around with numpys pinv, and it seems like it is quite slow. Compared to a homebrewed function:
def moore_penrose_inverse(data: np.ndarray) -> np.ndarray: if data.shape > data.shape: results = np.linalg.inv(np.dot(data.T, data)) @ data.T elif data.shape < data.shape: results = data.T @ np.linalg.inv(np.dot(data, data.T)) else: results = np.linalg.inv(data) return results
It seem to be way slower:
data = np.random.random((500,1000)) %timeit np.linalg.pinv(data) 335 ms ± 57.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) %timeit moore_penrose_inverse(data) 42.6 ms ± 1.64 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
An they do return the same result:
np.allclose(np.linalg.pinv(data), moore_penrose_inverse(data)) True np.allclose(np.linalg.pinv(data.T), moore_penrose_inverse(data.T)) True
For larger datasets, it also seemed that the pinv was around a factor 10 slower than the homebrewed.
I know numpy uses svd to calculate the pseudo inverse, but what are the benefits of this and is it worth compared to the loss of speed?
Solution – 1
Based on my understanding of the Moore–Penrose inverse, the code of Numpy of both
pinv, and profiling results, I can say that the performance gap comes from the computation of the singular value decomposition which takes about 96% of the time (at least on my machine).
I do not have a strong mathematical background on this, so please correct me if I am wrong: AFAIK, this SVD seems required for the pinv to provide correct/numerically-stable results when it is useful to compute this generalized inverse. AFAIK, the random matrix used in the example is not a pathological case where the generalized inverse is useful and a basic inverse can be used instead. I think it is related to the rank of the input matrix and random input matrices have the property of being full-rank while this method is useful when they are not.