python - Compute sum of power of large sparse matrix

Question

Given a query vector (one-hot-vector) q with size of 50000x1 and a large sparse matrix A with size of 50000 x 50000 and nnz of A is 0.3 billion, I want to compute r=(A + A^2 + ... + A^S)q (usually 4 <= S <=6).

I can above equation iteratively using loop

r = np.zeros((50000,1))
for i in range(S):
   q = A.dot(q)
   r += q

but I want to more fast method.

First thought was A can be symmetric, so eigen decomposition would help for compute power of A. But since A is large sparse matrix, decomposition makes dense matrix with same size as A which leads to performance degradation (in aspect of memory and speed).

Low Rank Approximation was also considered. But A is large and sparse, so not sure which rank r is appropriate.

It is totally fine to pre-compute something, like A + A^2 + ... + A^S = B. But I hope last computation will be fast: compute Bq less than 40ms.

Is there any reference or paper or tricks for that?

question from:https://stackoverflow.com/questions/66056756/compute-sum-of-power-of-large-sparse-matrix

Answer 1

Even if the matrix is not sparse this the iterative method is the way to go.

Multiplying A.dot(q) has complexity O(N^2), while computing A.dot(A^i) has complexity O(N^3).

The fact that q is sparse (indeed much more sparse than A) may help. For the first iteration A*q can be computed as A[q_hot_index,:].T.

For the second iteration the expected density of A @ q has the same expected density as A for A (about 10%) so it is still good to do it sparsely.

For the third iteration afterwards the A^i @ q will be dense.

Since you are accumulating the result it is good that your r is not sparse, it prevents index manipulation.

Eigen decomposition is good when you need to compute a P(A), to compute a P(A)*q the eigen decomposition becomes advantageous only when P(A) has degree of the order of the size of A. Eigen-decomposition has complexity O(N^3), a matrix-vector product has complexity O(N^2), the evaluation of a polynomial P(A) of degree D using the eigen decomposition can be achieved in O(N^3 + N*D).