Anull =
Can I solve a symetric system using only the lower triangular portion of the matrix?
Show older comments
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)?
7 Comments
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.
Paul
15 minutes ago
The dposv documenation states:
" A is DOUBLE PRECISION array, dimension (LDA,N)
!> On entry, the symmetric matrix A. If UPLO = 'U', the leading
!> N-by-N upper triangular part of A contains the upper
!> triangular part of the matrix A, and the strictly lower
!> triangular part of A is not referenced. If UPLO = 'L', the
!> leading N-by-N lower triangular part of A contains the lower
!> triangular part of the matrix A, and the strictly upper
!> triangular part of A is not referenced."
My reading is that on input A is 2D array with storage for the entire matrix. Even though only the upper or lower triangular portion of A is used based on UPLO, and therefore only the triangular portion has to be properly populated, there is no storage savings at the caller.
Upon further reflection, I think I misunderstood what the OP wants. I was focused on storage of A at the caller, though I don't think that's what the OP is concerned about.
Ah, so this question really is (or at least partially) about storage. Thanks for the links.
Again, you can submit an enhancement request for a hermitian/symmetric matrix class. That would help with storage in the workspace. However, I believe that the update to take advantage of such a class under the hood with the actual routines that do the linear algebra would be a significant effort. Maybe I'm wrong about that. No idea what the implications would be be on functions like eig, svd, etc., e.g., is there an eigenvalue solver for a symmetric matrix that takes an input in half-storage format?
You are right - the link I gave was wrong. This is the function where only the upper or lower triangular part of the matrix needs to be supplied (as a 1d-array of size n*(n+1)/2) :
Asking AI, LAPACK supports full, packed and band storage schemes.
Paul
5 minutes ago
Thanks for the update. AFAICT after a quick check, LAPACK supports those storage schemes for the solvers for Ax = B, for eigenvalue solvers, and Level-2 BLAS (matrix-vector). But I didn't see anything other than full storage supported for Level-3 BLAS (matrix-matrix). I did not review all LAPACK functionality, just a bit that is of interest to me.
Accepted Answer
Steven Lord
35 minutes ago
0 votes
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 Comment
John D'Errico
15 minutes ago
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.
More Answers (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 Comments
Elijah
1 minute ago
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
about 11 hours ago
Edited: John D'Errico
about 1 hour ago
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. And we need not have the matrix be triangular. 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.)
Edit: Does there exist a matrix A, such that rank(A) > rank(A+A'), and the answer is in fact trival given some thought. Consider any skew symmetric matrix for A. A skew-symmetric matrix has the property that A' = -A. For example...
A = [0 -1;1 0]
Then A+A' will be identically zero.
A + A'
And so we can see a decrease in rank from the sum. Next, suppose A is lower triangular. That would seem a little more difficult. But still possible.
A = [-1 0 0;-1 -1 0;-1 1 -1]
A is lower triangular. And with clearly non-zero diagonal elements, it will have full rank.
rank(A)
B = A+A'
rank(B)
In general however, it appears to be far more common for the sum of A+A' to have a higher rank than the original, at least based on my tests.
Finally, does there exist a lower triangular matrix A, such that rank(A) > rank(A+A'-diag(diag(A)))? This seems more difficult yet, as long as the matrix has non-zero diagonal elements. I don't believe that can happen, but a counter-example may exist.
3 Comments
Elijah
11 minutes ago
Elijah
11 minutes ago
John D'Errico
21 minutes ago
Edited: John D'Errico
19 minutes ago
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.
Christine Tobler
8 minutes ago
The mldivide and linsolve functions will interpret a triangular matrix as being the full matrix - there is no way to make them interpret it as defining a symmetric matrix.
The decomposition object allows one to specify that a symmetric matrix is defined by the upper or lower triangular part of a matrix:
A = sparse(tril(randn(100)));
b = randn(100, 1);
dA = decomposition(A, 'ldl', 'lower'); % same is possible with 'chol' for s.p.d. matrices
x = dA\b;
This requires one to specify the use of a symmetric decomposition: 'chol' for Cholesky in the case of a symmetric positive definite matrix, or 'ldl' for a general symmetric matrix.
By the way: You mention a symmetric complex-valued matrix. Do you mean symmetric (transpose(A) == A) or Hermitian (conj(transpose(A)) == A)? MATLAB's 'chol' and 'ldl' in the complex case only work for Hermitian matrices.
Categories
Find more on Sparse Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
