Working with very large matrices in numpy

Question

I have a transition matrix for which I want to calculate a steady state vector. The code I'm using is adapted from this question, and it works well for matrices of normal size:

def steady_state(matrix):
    dim = matrix.shape[0]
    q = (matrix - np.eye(dim))
    ones = np.ones(dim)
    q = np.c_[q, ones]
    qtq = np.dot(q, q.T)
    bqt = np.ones(dim)
    return np.linalg.solve(qtq, bqt)

However, the matrix I'm working with has about 1.5 million rows and columns. It isn't a sparse matrix either; most entries are small but non-zero. Of course, just trying to build that matrix throws a memory error.

How can I modify the above code to work with huge matrices? I've heard of solutions like PyTables, but I'm not sure how to apply them, and I don't know if they would work for tasks like np.linalg.solve.

Being very new to numpy and very inexperienced with linear algebra, I'd very much appreciate an example of what to do in my case. I'm open to using something other than numpy, and even something other than Python if needed.

https://stackoverflow.com/questions/14351255/techniques-for-working-with-large-numpy-arrays — AMC, Mar 10 '20 at 00:30
Does this answer your question? [Very large matrices using Python and NumPy](https://stackoverflow.com/questions/1053928/very-large-matrices-using-python-and-numpy) — AMC, Mar 10 '20 at 00:30
Not sure this will completely solve the question, but consider checking out [dask library](https://dask.org/) — FBruzzesi, Mar 10 '20 at 08:08

score 0 · Answer 1 · answered Mar 14 '20 at 06:42

Here's some ideas to start with:

We can use the fact that any initial probability vector will converge on the steady state under time evolution (assuming it's ergodic, aperiodic, regular, etc).

For small matrices we could use

def steady_state(matrix):
    dim = matrix.shape[0]
    prob = np.ones(dim) / dim
    other = np.zeros(dim)
    while np.linalg.norm(prob - other) > 1e-3:
        other = prob.copy()
        prob = other @ matrix
    return prob

(I think the conventions assumed by the function in the question is that distributions go in rows).

Now we can use the fact that matrix multiplication and norm can be done chunk by chunk:

def steady_state_chunk(matrix, block_in=100, block_out=10):
    dim = matrix.shape[0]
    prob = np.ones(dim) / dim
    error = 1.
    while error > 1e-3:
        error = 0.
        other = prob.copy()
        for i in range(0, dim, block_out):
            outs = np.s_[i:i+block_out]
            vec_out = np.zeros(block_out)
            for j in range(0, dim, block_in):
                ins = np.s_[j:j+block_in]
                vec_out += other[ins] @ matrix[ins, outs]
            error += np.linalg.norm(vec_out - prob[outs])**2
            prob[outs] = vec_out
        error = np.sqrt(error)
    return prob

This should use less memory for temporaries, thought you could do better by using the out parameter of np.matmul. I should add something to deal with the last slice in each loop, in case dim isn't divisible by block_*, but I hope you get the idea.

For arrays that don't fit in memory to start with, you can apply the tools from the links in the comments above.

Working with very large matrices in numpy

1 Answers1