Making a Matrix Strictly Diagonally-Dominant
조회 수: 234 (최근 30일)
이전 댓글 표시
Eligijus Rupsys
2020년 3월 20일
편집: Khandoker Mohammad Faisal Karim
2021년 8월 21일
Hello everyone ! Hope everyone is safe and healthy in light of the recent developments.
I'm trying to create a matlab code that takes a given matrix, firstly tests if the matrix is diagonally-dominant, if it is not, then the matrix rows are randomly swapped and the test is carried out again until the matrix is diagonally dominant. I know that this is definitaly not the most efficient way to convert a matrix to be diagonally dominant, however it is the best approach i could come up with the MATLAB knowledge that i know.
This is a script that tests if the matrix is diagonally dominant;
function [isdom] = IsDiagDom( A )
isdom = true;
for r = 1:size(A,1)
rowdom = 2 * abs(A(r,r)) > sum(abs(A(r,:)));
isdom = isdom && rowdom;
end
if isdom == 0
disp (['Matrix A is not diagonally-dominant']);
elseif isdom == 1
disp (['Matrix A is diagonally-dominant']);
end
end
And this is the script that im trying to make work that if the matrix is not diagonally dominat, the rows are randomly swapped and tested till it becomes diagonally dominant;
When i try to run the script, it says ' Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check
for mismatched delimiters. ' in the 'for ...' line
function [ A ] = DiagDom
A = [ 4 -28 -7 1; 4 -1 10 -1; -4 0 -3 11; 19.375 5 8 -3 ];
for IsDiagDom( A );
if isdom == 1
disp (['Matrix A is diagonally-dominant']);
elseif A = A(randperm(size(A, 1)), :); % Randomly swaps rows
return
end
end
end
Thank you a lot !
댓글 수: 1
John D'Errico
2020년 3월 20일
편집: John D'Errico
2020년 3월 20일
Please see my answer.
You should understand why it is that the use of random permutations is a bad idea. That is so because if the matrix is even remotely large, and here a 15 by 15 matrix is essentially huge, then the number of permutations will be immense.
The number of permutations of N numbers is factorial(N). If N is 15, then we see
factorial(15)
ans =
1307674368000
So over 1 TRILLION permutations are possible. A matrix with 20 rows would have
factorial(20)
ans =
2.43290200817664e+18
In words, that is...
two quintillion, four hundred thirty two quadrillion, nine hundred two trillion, eight billion, one hundred seventy six million, six hundred forty thousand
row permutations possible for a matrix with 20 rows.
So why are random row permutations a bad idea? Because there is such a simple non-random solution possible. It takes little more than a call to the function max to find that permutation, and to see if a permutation does exist at all.
채택된 답변
Sriram Tadavarty
2020년 3월 20일
HI Eligijus,
The way the for loop is used here caused the issue. Update the second part of code as below and it works:
function [ A ] = DiagDom
A = [ 4 -28 -7 1; 4 -1 10 -1; -4 0 -3 11; 19.375 5 8 -3 ];
while(1) % Perform infinite loop, till you find the diagonally dominant matrix
if IsDiagDom (A) % If this is diagonally dominant, disp and break the loop
disp (['Matrix A is diagonally-dominant']);
break;
else
A = A(randperm(size(A, 1)), :); % Randomly swaps rows
end
end
end
Hope this helps.
Regards,
Sriram
댓글 수: 2
K Bayoud
2021년 6월 22일
Hello!
There is a case where the matrix cannot be diagonally dominant even though we swap its rows.
It will be an infinit loop.
추가 답변 (2개)
John D'Errico
2020년 3월 20일
편집: John D'Errico
2020년 3월 20일
Can you solve this? Yes, sometimes, and there is no need for random permutations of the matrix. In fact, that is a poor solution, since there is indeed a simple solution that has no need for random swaps. But first...
A serious flaw in your problem is there are some matrices (easy to construct) that can NEVER be made diagonally dominant using simply row exchanges.
I'll paste in the important wording here:
"a square matrix is said to be diagonally dominant if, for every row of the matrix, the magnitude of the diagonal entry in a row is larger than or equal to the sum of the magnitudes of all the other (non-diagonal) entries in that row. "
So it is clearly true that there can easily be rows that can never satisfy that requirement. For example, consider the row vector:
Row = [10 10 1 1 1 1 10];
Suppose we made this to be the first row of the matrix? Well, then we must have 10 (the first element) being larger than the sum of the magnitudes of the other elements.
Likewise, if we made it the second row, or the last row, then we still have the same problem.
Row = [10 10 1 1 1 1 10];
sum(abs(Row(2:end)))
ans =
24
As long as that row is in the matrix, there is NO possible re-ordering that will make the matrix diagonally dominant. It simply cannot happen, because no matter which row you swap it to, it will always fail the requirement.
In fact, I could have made it even simpler. How about this row vector?
Row = ones(1,10);
Where would you swap that row to, such that the matrix will now be diagonally dominant? You cannot ever find a solution, even disregarding all other rows of the matrix. If your matrix has such a row, then you can never succeed.
Even more interesting though, is we can show that any row can only ever live in ONE position, IF the matrix is to be strictly diagonally dominant. That is because we need only find the largest element in any row in abolute magnitude. The position of that element tell you which row it needs to be in. Case closed.
Now, having said that, why did I say that it is possible to find a non-random solution SOME of the time? In fact, it is simple to derive such an algorithm. Consder ANY row. Find the maximum absolute value of that element. If that value exceeds the absolute sum of the remainder of the row elements then that row is POTENTIALLY a candidate for being in a diagonally dominant matrix.
Is there a problem here? Well yes. suppose that two rows must both be row 1? Consider these two rows:
Row1 = [10,rand(1,5)];
Row2 = [11,rand(1,5)];
There is only one position for either of those rows to live in, IF the corresponding matrix will be DD. If your matrix has both of those rows, then you are stuck, up a creek without a paddle. There would be no solution.
A = randi(5,[5,5]) + eye(5)*100;
A = A(randperm(5),:)
A =
4 2 1 101 1
5 104 3 4 1
105 5 2 2 3
5 4 105 4 4
2 4 4 1 101
If we consider the matrix A, as I created it there is CLEARLY a permutation that will yield a diagonally dominant matrix as a solution. What is it? All we need is ONE simple call to the function max do most of the work.
[maxrow,maxind] = max(abs(A),[],2)
maxrow =
101
104
105
105
101
maxind =
4
2
1
3
5
Now, CAN the matrix be made to be diagonally dominant? there are two tests necessary. First, we need for this to be true:
all(maxrow > (sum(abs(A),2) - maxrow))
ans =
logical
1
Think about why it is necessary. In order for the matrix to be STRICTLY diagonally dominant, we need that strict inequality too. A simpler >= will not suffice.
Next, we need for the vector maxind to be a permutation of the numbers 1:5. We might write it like this:
isequal(sort(maxind),(1:numel(maxind))')
ans =
logical
1
There are other ways I could have written that test, but it is sufficient and necessary.
Regardless, now what is the solution? SIMPLE!
A(maxind,:) = A
A =
105 5 2 2 3
5 104 3 4 1
5 4 105 4 4
4 2 1 101 1
2 4 4 1 101
As such, the code to perform what you asked for is both trivial to write and fast to execute.
function A = makeDD(A)
% takes a square matrix A and permutes the rows if possible so that A is diagonally dominant
[n,m] = size(A);
if n~= m
error('A is not square')
end
% test to see if a valid permutation exists
[maxrow,maxind] = max(abs(A),[],2);
if all(maxrow > (sum(abs(A),2) - maxrow)) && isequal(sort(maxind),(1:numel(maxind))')
% success is both possible and easy to achieve
A(maxind,:) = A;
else
disp('Sorry, but this matrix can never be made to be diagonally dominant')
A = [];
end
end
Let me test it now.
A = magic(3)
A =
8 1 6
3 5 7
4 9 2
makeDD(magic(3))
Sorry, but this matrix can never be made to be diagonally dominant
ans =
[]
As you can see, even though A has distinct maximal elements which are larger than the rest in that row, AND they fall in distinct columns, it still fails the other test, that for the second row of A, we must have had 7 > (3+5).
Change A just a tiny bit by changing one element, we can succeed however.
A(2,3) = 12
A =
8 1 6
3 5 12
4 9 2
makeDD(A)
ans =
8 1 6
4 9 2
3 5 12
In all of this you need to see the solution is always trivial to find, IF one exists, and that it requires no random permutations, Finally, see that the solution, if it DOES exist, is unique.
As I said, the code I wrote is blazingly fast, even for huge matrices. Consider this case for a 100x100 row-randomized matrix. Again, I'll construct it where the matrix is known to have a solution.
A = randi(5,[100,100]) + eye(100)*1000;
A = A(randperm(100),:);
tic,Ahat = makeDD(A);toc
Elapsed time is 0.002210 seconds.
So 0.002 seconds to solve a problem that if we used random permutations would take the lifetime of the universe to solve, even using a computer the size of the entire universe.
factorial(100)
ans =
9.33262154439441e+157
댓글 수: 6
Khandoker Mohammad Faisal Karim
2021년 8월 21일
편집: Khandoker Mohammad Faisal Karim
2021년 8월 21일
"Diagonally dominant: The coefficient on the diagonal must be at least equal to the sum of the other coefficients in that row and at least one row with a diagonal coefficient greater than the sum of the other coefficients in that row."
To implement this:
if all(maxrow >= (sum(abs(A),2) - maxrow)) && any(maxrow > (sum(abs(A),2) - maxrow))...
&& isequal(sort(maxind),(1:numel(maxind))')
A(maxind,:) = A;
else
disp('Sorry, but this matrix can never be made to be diagonally dominant')
A = [];
end
Anita Mundra
2020년 4월 23일
hell Eligijus Rupsys
i am also looking for such loop code, but unable to trace out.
did you get solution for your problem.
as the code taht is mentioned is not running.
it is giving an error
"Error in DiagDom (line 5)
if IsDiagDom (A) % If this is diagonally dominant, disp and break the loop"
if you can please share the code with me.
my email id is: anitamundra26@gmail.com
댓글 수: 0
참고 항목
카테고리
Help Center 및 File Exchange에서 Creating and Concatenating Matrices에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!