Matrix fill-in when solving sparse linear systems

I am solving the sparse linear system y=c'\b where c (10'000x10'000 density 0.0076) is the result of a Cholesky factorization while b (10'000x23'000 density 0.0027) is a generic sparse matrix. The result y is 10'000x23'000 and it has a density equal to 0.0711. The problem is that many of the non-zero elements are between 1e-15 and 1e-50 while (I suppose) they should be zero.

Things are even worse when I subsequently solve x=c\y in which case the density of x is 0.77 but again most of the elements are lower than eps.

Any idea why this happens?

Indeed, if I use an iterative method (such as the Biconjugate gradients method) and solve c\b(:,i) I get a vector where those near-zero elements are actually zero but I cannot iterate the iterative method over the 23'000 columns of b as it takes ages.

Thanks for your help

Answer 1

[Richard Brown]

Firstly, is there any reason to believe that x should be sparse? Just because your matrix A = C*C' and your right hand sides b are both sparse, doesn't mean that A^(-1)*b should be ... in fact I'd expect that x would not be sparse unless your matrix A has some kind of special structure.

If x should in fact be sparse, then you could solve your problem in batches, maybe doing 100, or 500 or 1000 columns at a time, and thresholding the solution x before storing it

[Francesco Finazzi]

What I am solving is something like W-Z*A\Z' where W and A are symmetric and s.p.d. sparse matrices. I believe that Z*A\Z' has the same sparseness structure of W. In fact the zero elements of W are around 1e-15 in the result of Z*A\Z'. I will set these elements to zero but I am quite sure that the fill-in during A\Z' increases the computational time unnecessarily.

Thanks for your answer.

[Richard Brown]

Can you be a little more precise? Do you mean you are trying to evaluate W - ZA^-1Z^ T and that the properties of your problem should guarantee that both terms have the same sparsity? Is there anything special about Z? Are you able to make your matrices available?