matlab code to transform linear systems to strictly diagonally dominant matrix
3 views (last 30 days)
Show older comments
segun egbekunle
on 18 Mar 2016
Edited: John D'Errico
on 16 Apr 2016
I need matlab syntax to transform a linear system Ax=b to strictly diagonally dominant matrix. For example given A=[6 5 7; 4 3 5; 2 3 4] b=[18 12 9]' I want to transform the coefficient matrix A to another matrix B such that matrix B is strictly diagonally dominant and b to another vector d
2 Comments
Accepted Answer
John D'Errico
on 16 Apr 2016
Edited: John D'Errico
on 16 Apr 2016
I'm sorry, but this is not possible to do. In fact, you can prove the impossibility, by a simple counterexample.
Let the matrix A be ones(3,3). This matrix is singular, worse, it has a rank of 1. No linear transformation that you can apply to A is sufficient to make A STRICTLY diagonally dominant, since a strictly diagonally dominant matrix would be NON-SINGULAR.
Sorry, but mathematics (in this case, linear algebra) is simple, and unrelenting. It does not allow you to do just anything.
In your case, we have a non-singular matrix A, so there are many (all trivial) solutions.
A=[6 5 7; 4 3 5; 2 3 4];
rank(A)
ans =
3
As I said, there are trivial solutions, and many of them, IF A is non-singular. For example, all you need do is multiply both sides by the inverse of A. Thus, if we start with the linear system:
A*x = b
then as long as A is not singular (and square) then we can write it as
B*x = eye(size(A))*x = inv(A)*b
So, your new coefficient matrix is
B = eye(size(A))
and the transformed right hand side (b) is now
inv(A)*b
That satisfies the requirements of this question, as long as A is non-singular. Since you cannot get any more diagonally dominant than an identity matrix, this is the answer, and no answer can be better.
Another simple answer is to use pinv, which for your purposes is again only valid if A is non-singular, if the created matrix B is to be strictly diagonally dominant. As I said, no linear transformation of A will suffice for the bad case. So, again we have
B = eye(size(A));
bhat = pinv(A)*b;
so the new linear system is essentially
eye(size(A))*x = pinv(A)*b
As it turns out, for a non-singular matrix A, pinv(A) is mathematically equivalent to inv(A). pinv is arguably a little better behaved for some nearly singular matrices, but if the matrix is nearly singular, you are in deep trouble anyway with any approach.
Alternatively, one can use a QR factorization of A to do the transformation. It will take slightly more effort to do (but really only a few extra characters.) If you want the system to be stable, then a pivoted QR would be a better choice than a simple QR, but a pivoted QR will implicitly re-sequence the unknown vector x, in order to get the form you seem to desire.
If you do not re-sequence the vector x, then the solution may be less numerically stable for SOME problems. You can resolve the re-sequencing simply enough. For example...
[Q,R,E] = qr(A);
Here, E is a permutation matrix. The necessary linear transformation is now
B = eye(size(A));
bhat = E'*(R\Q')*b;
Or you could use an LU decomposition. Again, a pivoted LU is a better choice for stability.
You could also use an eigenvalue decomposition in a similar way, as long as A has a complete set of eigenvalues and eigenvectors.
Do I need to go on writing what could be a basic linear algebra text?
0 Comments
More Answers (1)
Ced
on 16 Apr 2016
Without giving away the whole solution, one way would be to diagonalize your linear system. What you need to do is thus to find the transformation matrices P such that
B = P*A*P.'
is diagonal (and apply the same transformation to b), i.e. using eig.
If it is not diagonalizable, you can use svd to transform it in such a way, or see here. Note that all this is only possible if A is non-singular (otherwise, you can reshape A with linear row operations such that a row of zeros appears, which is obviously not strictly diagonally dominant).
0 Comments
See Also
Categories
Find more on Operating on Diagonal 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!