How can I approximate this large matrix?

I have a sparse symmetric matrix M, roughly 10^6 x 10^6, but only about 100 values in each column (or row) are non-zero, so they'll all fit memory / disk space. Values are complex, and all magnitudes are < 10^-1 I need an approximation of

M * inv(eye() - M) = M + M^2 + M^3 + M^4 + ...

The inverse in there won't be sparse, but the produc hopefully will be.

So, if I have the nonzero values and their locations tabled up, what's the best way to go about this calculation in MATLAB?

Are there library tools avaiable for large sparse matrices?

And should I try to approximate the LHS or the RHS? I know the RHS can be approximated by calculating M^2, M^4, M^8,... and then multiplying M * (I+M) * (I+M^2) * (I+M^4) * (I+M^8) = M + M^2 + ... + M^16. But I'm hoping there might be some iterative way to approximate the LHS. And in either case I'd still need to know how to set M up as a sparse matrix and take advantage of M's sparseness.

Answer 1

[Bruno Luong]

What about

T = M;
niter = 5:
for k=1:niter
    T = M*(speye(size(M)) + T);
end

This returns

T = M + M^2 + .... + M^6

[Jerry Guern]

Yeah, I thought about doing something like that...

It would be faster to calculate M^2, M^4, M^8, then just multiply

M * (I+M) * (I+M^2) * (I+M^4) * (I+M^8) = M + M^2 + ... + M^16

but that code is just doing straight matrix operations on M, not taking any particular advantage of its spareness. Do tools like that exist in MATLAB? This operation will take forever if we don't do that.

[Bruno Luong]

"but that code is just doing straight matrix operations on M, not taking any particular advantage of its spareness"

Of course it takes the advantage; The operation * and + DO take the advantage of the manipulated matrix is sparse.

The approximation formula you want to use has nothing specific to sparse matrix, it can be applied on full matrix as well provide the norm(M) < 1.

[Jerry Guern]

But storing M as a matrix in the first place is going to be a problem because of its size. How would I set it up as explicitly sparse in the first place, analogous to how speye(n) works? This is new to me on MATLAB.

[Jerry Guern]

I understand that these formula have nothing to do with sparseness and can be applied to full matrices. But I can't do any operations on a matrix that's too big for memory. So I also need tools for handling the matrix explicitly as sparse, so that only the non-zeros are stored.

[Bruno Luong]

But MATLAB sparse stores only non-zeros elements.

If your matrix are big and not sparse enough to be store as sparse in the memory of your PC, then you are simply stuck.

I don't know what you meant by tool, MATLAB matrix library comprises a comprehensive array of functions that can deal with sparse matrix.

[Jerry Guern]

Ah, thank you, this was helpful on the sparse matrix handling part. I realize now I probably should have posted this as two separate questions, one on the sparse matrix tools and one on the approximation.