Anull =
Why is my sparse complex lower triangular matrix singular when the full matrix is not?
Mostrar comentarios más antiguos
I have a sparse, complex-valued symmetric matrix, Afull, where I have stored only the lower triangular component (and diagonal) as matrix A.
spy(A) looks like:
mldivide of A is poorly conditioned:
>> u=A\b;
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 7.011677e-19.
I get the same error when I specify that the matrix is lower-triangular in LINSOLVE:
>> opts.LT=true;u=linsolve(full(A),b,opts);
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.
RCOND = 7.011677e-19.
When I re-build the full matrix from the LT(+diagonal) component A:
>> Afull=A+A';
>> for ii=1:size(A,1)
Afull(ii,ii)=A(ii,ii);
end
spy(Afull) looks like:
>> ufull=Afull\b;
>>
The matrix condition is now fine.
Is it possible to get Matlab to solve the système with ONLY the lower component (A)?
2 comentarios
I've often (actually, not that often) wondered if Matlab should offer a hermitian/symmetric class to minimize storage of lots of large, hermitian/symmetric, 2D arrays. You can always submit an enhancement request for such a feature.
This comment states: "Many matrix solvers let you pass just half of the symmetric matrix for solution, I was a bit surprised to find that I can't seem to do this with Matlab." Do you have any links to such solvers with code available to the public? I checked the BLAS routines, which I think Matlab uses under the hood, and I didn't see a matrix multiply routine with this capability.
The original LAPACK functions only expect the lower or upper triangular part of the symmetric matrix, depending on a flag "UPLO" passed to the solver. E.g. for positive-definite matrices
I guess MATLAB will half the (full) user matrix accordingly.
Respuestas (3)
Consider A = [0 0;1 0]. Then A is lower triangular and singular. But A + A' = [0 1;1 0] is regular. So why do you think Afull should have the same rank property as A ?
Zeros on the main diagonal of A immediately make A singular while this is not true for A + A'.
2 comentarios
Elijah
hace 25 minutos
What you solve in your code is A*u = b and [(A+A')-diag(diag(A))]*u = b. There are cases where A is singular or badly conditionned (if the diagonal of A has zero or small entries) whereas (A+A')-diag(diag(A)) is regular (see my example). Thus the MATLAB results from above are no surprise for me.
John D'Errico
hace alrededor de 3 horas
Editada: John D'Errico
hace alrededor de 2 horas
You don't show your matrix. Only a picture of it, and a picture is not quite worth a thousand words all of the time. ;-)
But consider this matrix:
format long g
A = [1 0 0;1 1e-5 0;1 1 1e-16]
cond(A)
Is A singular, or nearly so? Yes, it definitely is at least numerically singular. It does not have szeros on the diagonal, but because A is only a 3x3, I had to push things a bit. Had I created a much larger matrix, I could have been far more creative with diagonal elements that were not nearly so nasty.
But now I'll create a new matrix from A, that is now symmetric using the same computation you did. Surely by your logic, B must also be singular?
B = A + A'
cond(B)
B is not even remotely singular!
What happened? What happened is that your conclusion does not follow! I could have done as easily with a 2x2 matrix.
A = [1 0;1 1e-8]
B = A + A'
[cond(A),cond(B)]
A similar example could be contrived with a much larger martrix, since the diagonal elements do not need to be at all small for a large matrix.
Or, I could have contrived an example with a 3x3 matrix which is exactly provably singular. For example...
A = [4 5 2;6 8 3;6 7 3]
rank(A)
In fact, A is actually truly singular, not just numerically singular. We can see that because there exists a simple vector Anull, that completely kills A.
Anull = null(sym(A))
A*Anull
The product is EXACTLY zero. And if you look at the 1st and third columns of A, you will see that must be true.
However, consider the symmetrized version B.
B = A + A';
cond(B)
B is not remotely singular. Again, your conclusion is completely false. Adding a matrix to its transpose does not maintain the rank of the matrix. The resulting sum will often (but not always) have larger rank than the original. Even in a trivial case...
A = [1 1 1;0 0 0;0 0 0] % A has a rank of exactly 1
rank(A)
B = A + A'
rank(B)
However, B has a rank of 2.
You should be able to trivially construct examples where the rank does not change though. Consider
A = [1 0; 0 0]
B = A + A'
the rank of A is 1, but the rank of B is clearly identical to the rank of A. Of more interest would be to find an example where the rank of the symmetrixed sum is less than the original. I'm not sure I can find a simple example of a decrease in rank, though I'm not willing to claim it impossible without some thought. That is, is it always true that
rank(A) <= rank(A+A')
My gut says this may be true, but It is time for breakfast, and my gut may be a liar. (Note that what you did on the diagonal changes nothing about my statements here.)
3 comentarios
Elijah
hace 41 minutos
Elijah
hace 40 minutos
John D'Errico
hace 21 minutos
Editada: John D'Errico
hace 19 minutos
No. You still do not understand!
It is not just that you wish to store only the lower triangle of a symmetric matrix. The fundamental problem is that the rank of the lower triangular part is NOT the same as the fully symmetric version. So trying to use the solve A\y is meaningless, since the lower triangular part does not have full rank.
So the comparison you were trying to make is essentially meaningless.
The real question appears to be what you added at the end, after you modified your question. That is...
Can you solve the linear system (normally done as B\y) if ONLY the lower triangle of B is given and stored in memory?
The answer is itself trivial. Yes, you can, by forming the complete symmetric form of the matrix. However, you need not make the matrix a full matrix. It can remain in sparse form. For example...
n = 1000;
A = tril(sprand(n,n,0.01));
spy(A)
A is sparse and not even remotely full rank.
condest(A)
rank(full(A))
rhs = randn(n,1);
x = A\rhs;
This result was meaningless garbage, since A was singular! However, if I do this:
B = A + A' - diag(diag(A));
whos B
As you can see, B is still stored in sparse form. There is no need to make B a full matrix!
spy(B)
As you can see, I formed B as the symmetrized version. (I did not have any need to use a loop however.)
condest(B)
rank(full(B))
So even while A was rank deficient, B is fine.
x = B\rhs;
As you can see, there was no warning in the solve. backslash is perfectly happy.
Steven Lord
hace 36 minutos
0 votos
The backslash operator (mldivide, \) does not have a way to specify that the coefficient matrix you pass to it actually only represents half the system. Nor does the linsolve function (specifying the LT field in the options structure indicates that the matrix is lower triangular, not that MATLAB should use only the lower triangular part.)
Since you have a sparse coefficient matrix, you might want to use one of the iterative solvers. All of these allow you to pass in a function handle that returns a function of the coefficient matrix A and a vector x (often A*x or A'*x.) See the "Using Linear Operators Instead of Matrices" section on the documentation page linked above. If your actualCoefficientMatrix is A+A', you could pass a function handle that computes A*x+A'*x in to compute actualCoefficientMatrix*x without explicitly creating actualCoefficientMatrix, just using A.
1 comentario
John D'Errico
hace 16 minutos
Actually, the desired symmetrized matrix is not A+A', but
A + A' - diag(diag(A))
A+A' produces a result where the diagonal elements are multiplied by 2.
hace alrededor de 8 horas
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
