Fast solver for linear system Ax=b when only x(1) is needed

Hello,

Since i need to solve the following problem many times I am looking for the fastest way for the solution of Ax=b, where:

• A(nxn): complex , positive definite , symmetric , poorly sparse(approx 1/3 of elements are 0)
• b(n): [ 1 0 0 0 ... 0]
• n= 6 or 8 ;

Plus,i just need to compute the solution for the first term x(1) ONLY . What is in your opinion the best method to achieve that?

Thank you very much in advance!

[Andrei Bobrov]

xout = A\b

[rengell]

thanks but I know how to solve the system in general. my concern here is performance.I just need something fast considering that i need only one term of the solution and not the entire vector!

Answer 1

[John D'Errico]

"Poorly" sparse, thus only 33% zero, is NOT sparse in any meaningful sense of the word. In fact, it may take more memory to store that array than a dense array. Similarly, your solve time will probably be no faster for this array than a dense one, and very likely slower.

Solving for x(1) only is also not going to be anything useful in terms of a speedup. Almost all of the time is invested in the factorization. The actual solve is fast. Yes, you CAN solve for x(1) only, and you MAY gain a fractional speed up in theory.

However, in order to realize these theoretical gains, you would need to use a cholesky decomposition, returning that factor to you. Then you need to do a back solve (using \), and finally solve for that one element that you need. Again, my guess is the extra overhead of you calling each of these tools successively, plus explicitly returning a cholesky factor, may all be more than a direct call, just

x = A\b;
x = x(1);

The only place you might gain is that \ needs to determine that your matrix is actually SPD, whereas you already know it is.

So perhaps a slight gain may come from using linsolve instead of backslash, since here you can tell the code that A is SPD.

A quick test might help. First, is a 30% sparse matrix even a space savings? No.

A = sprand(5000,5000,.7);
FA = full(A);

whos
  Name         Size                  Bytes  Class     Attributes
  A         5000x5000            201372440  double    sparse    
  FA        5000x5000            200000000  double

A 30% sparse matrix takes MORE space to store, as I suggested. Don't even bother with sparse.

I won't even test the speed for a sparse call there, since creating a known SPD matrix that is 30% zero is a complete waste of time. (Sorry, but it is.) Just make that matrix dense. People think that because they have a few zeros in their matrix, that it is sparse. NO. NO. NO. Don't waste your time with sparse here.

n = 8;
b = [1;zeros(n-1,1)];
A = rand(n,n);
Aspd = A'*A; % clearly dense SPD matrix

opts.SYM = true;opts.POSDEF = true;

timeit(@() linsolve(Aspd,b,opts))
ans =
   4.4486e-05

timeit(@() Aspd\b)
ans =
   1.3111e-05

In another test with n = 5000, linsolve was actually 10% faster than backslash. But for small matrices like an 8x8 matrix, backslash was seriously faster as you can see here. The presence of a few zeros would not be significant.

How about other tools, like lsqr? NO. Don't bother. lsqr is useful for large sparse problems, where your solver will need to wander into virtual memory space. For a 6x6 or 8x8, DON"T BOTHER. lsqr will be significantly slower. Iterative methods are NOT good for small problems. (A quick test for that same 8x8 problem as above took 100 times as long with lsqr.)

Just use A\b, then take the first element. Again, remember that the time spent here is mainly in the factorization, even for a small matrix like an 8x8 system. Factoring the matrix is an O(n^3) operation. The actual solve is an O(n^2) operation. That suggests that the factorization part is roughly 87% of the total time needed. Therefore the time savings from taking only one element of that result is still quite small, roughly 50% of that 13%.

(Were you to write high quality compiled C code for this, you might gain some speed. That is something I cannot compare, but even so, unless you seriously know what you are doing, you won't beat MATLAB here either.)

[rengell]

Hi, Thank you deeply for your very detailed answer. Very clear! Thruth is that I did run some tests in the while. As you pointed out sparse matrix is making no sense in this context. For now I find out that the fastest way is through LU decomposition, considering that for complex cases Cholesky cannot be adopted.

[John D'Errico]

Where did you say it was complex? Oh, I guess you did mention that.

Why do you think that chol cannot be used for complex matrices though?

A = rand(4) + rand(4)*sqrt(-1);
A = A'*A;
chol(A)
ans =
       1.7294 +          0i       1.0432 -    0.26648i      0.94045 -    0.45878i       1.4522 -    0.53865i
            0 +          0i        1.352 +          0i      0.16539 +    0.22396i      0.24492 +    0.16164i
            0 +          0i            0 +          0i      0.70945 +          0i     -0.15539 -    0.27396i
            0 +          0i            0 +          0i            0 +          0i      0.21471 +          0i

Chol has no problem at all.

My guess is your matrix was not created as truly SPD. This is a common problem that people have. They think their matrix is, but it is not so.

[rengell]

The matrix A = A'*A is different from the initial matrix. You will notice that by finidng different solutions for the problem.

In fact cholesky factorization requires matrix to be hermitian.My Symmetric complex matrix is not hermitian (A is different from A') because diagonal terms are complex too. That's why i think it's not possible to use cholesky factorization in this case.

My matrix is truly SPD:

A=[2469175345.45771 + 246917534.545771i 5937774.84578943 + 593777.484578943i -2469105733.36118 - 246910573.336118i -14844528.7413451 - 1484452.87413451i 0.00000000000000 + 0.00000000000000i 0.00000000000000 + 0.00000000000000i;5937774.84578943 + 593777.484578943i 10699569438.2229 + 1069956943.82229i 14844528.7413451 + 1484452.87413451i -10699553373.8215 - 1069955337.38215i 0.00000000000000 + 0.00000000000000i 0.00000000000000 + 0.00000000000000i; -2469105733.36118 - 246910573.336118i 14844528.7413451 + 1484452.87413451i 2765484822.65343 + 258769913.633600i -4750243.35624848 - 546276.224997304i -296286273.357373 - 11851450.9342949i -2968879.62714829 - 118755.185085932i; -14844528.7413451 - 1484452.87413451i -10699553373.8215 - 1069955337.38215i -4750243.35624848 - 546276.224997304i 11983513709.2159 + 1121314714.66201i 2968879.62714829 + 118755.185085932i -1283938916.19484 - 51357556.6477936i; 0.00000000000000 + 0.00000000000000i 0.00000000000000 + 0.00000000000000i -296286273.357373 - 11851450.9342949i 2968879.62714829 + 118755.185085932i 305035715.550811 + 12725002.9233380i -599150.333205392 + 11336.8560519182i; 0.00000000000000 + 0.00000000000000i 0.00000000000000 + 0.00000000000000i -2968879.62714829 - 118755.185085932i -1283938916.19484 - 51357556.6477936i -599150.333205392 + 11336.8560519182i 1321214194.17106 + 55084763.1575303i;];

Answer 2

[Greg Dionne]

This seems to work o.k.-ish.

det(A(2:end, 2:end)) / det(A)

I have no idea how well it compares against the backslash operator though.

[rengell]

hey, no cramer's rule is faaar to slow,i was not even considering it. It's almost 40 times slower tha LU. Thanks anyway!

[Greg Dionne]

Well... someone had to say it! Although I was secretly hoping the recursive block matrix form might have worked ok (e.g. det([A B; C D]) = det(A)*det(D-C*(A\B)) if A is invertible...