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

3 次查看（过去 30 天）
rengell
rengell 2015-3-20
编辑： Greg Dionne 2015-3-23
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!
  2 个评论
rengell
rengell 2015-3-20
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!

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

John D'Errico
John D'Errico 2015-3-20
编辑：John D'Errico 2015-3-20
"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.)
  3 个评论
显示 1更早的评论
John D'Errico
John D'Errico 2015-3-20
编辑：John D'Errico 2015-3-20
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
rengell 2015-3-20
编辑：rengell 2015-3-20
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;];

请先登录，再进行评论。

Greg Dionne
Greg Dionne 2015-3-20
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.
  2 个评论
rengell
rengell 2015-3-20
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
Greg Dionne 2015-3-23
编辑：Greg Dionne 2015-3-23
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...

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by