Copyright (c) 2014, Feng Cheng Chang
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
Feng Cheng Chang (view profile)
Thank you, John, for your comment and advice.
In order not to use 'det' to verify matrix singular, the following simple code is developed to deal with this matter.
function p = M_p(A)
% F C Chang 11/20/14
M = A; p = [ ];
for k = 1:size(A,1),
n = size(M,1);
[mp,rcm] = max(abs(M(:)));
[r,c] = ind2sub([n,n],rcm);
rx = [1:r-1,r+1:n];
cx = [1:c-1,c+1:n];
p = [p,M(r,c)*(-1)^(r+c)];
if mp < 1.e-10, p; return, end;
M = M(rx,cx)-M(rx,c)*M(r,cx)/M(r,c);
end;
Given a matrix A, the total or partial elements of an array p are generated. It follows:
(1) The determinant of A is found to be equal to the product of total elements of p, if A is non-singular.
(2) The last element of p is found dip down steeply toward zero, then the matrix A is said to be singular.
To show the merit of the derived code, some examples using this code are shown below.
It is noted that even A is regard as singular, yet det(A) is found to be very hugh.
>> A=magic(7); p=M_p(A); p, prod(p), det(A),
p =
-49.0000 47.1837 -42.0484 41.5685 -40.9333 -47.2149 -44.5643
ans =
-3.4805e+011
ans =
-3.4805e+011
>> A=magic(8); p=M_p(A); p, prod(p), det(A),
p =
64.0000 -62.1250 -12.8169 0.0000
ans =
1.8105e-010
ans =
0
>>M=randn(100,99); M(:,100)=M(:,99)+eps; M([94:100],[94:100])
>>ans =
-0.4511 -0.5694 -1.3372 0.3640 -0.7417 0.6269 0.6269
-0.5323 -0.7567 -0.0995 -0.8734 0.0907 -0.1194 -0.1194
0.4189 1.9004 -1.1543 0.1443 1.4559 0.0251 0.0251
-0.6081 -0.1819 1.3628 0.7203 -0.1243 -1.0554 -1.0554
2.5348 0.5820 -0.3999 0.5663 0.1373 -1.0932 -1.0932
0.7985 -0.9134 1.5652 -2.1671 -1.6688 1.1840 1.1840
0.8397 -0.1540 1.6039 -0.5653 0.2192 -0.7909 -0.7909
>>A=M; p=M_p(A); p([94:100]),,prod(p), det(A)
ans =
-5.0236 4.7189 4.3518 -5.4064 -3.5973 -1.7538 0.0000
ans =
-5.1420e+063
ans =
-7.5012e+063
>>A=10*M; p=M_p(A); p([94:100]), prod(p), det(A)
ans =
-50.2364 47.1890 43.5181 -54.0645 -35.9727 -17.5377 -0.0000
ans =
5.2564e+163
ans =
-1.6784e+163
>>A=0.1*M; p=M_p(A); p([94:100]), prod(p), det(A)
ans =
-0.5024 0.4719 0.4352 -0.5406 -0.3597 -0.1754 0.0000
ans =
-4.8905e-038
ans =
-7.4221e-038
John D'Errico (view profile)
I never said you could use norm to test for singularity of a matrix!
I did say you could use rank, cond, or svd to determine that status, AND be far ahead of the use of det.
I already showed how poor det is, using several examples. It is simply a terrible idea to use that tool. Yes, I know that your teacher told you about it, or some book you read. The fact is, they probably learned it from someone before them, and they are just propagating the use of a terrible tool for the purpose.
When you use det in floating point arithmetic, is is simply a bad thing to do. This is what those people who taught you did not understand. They thought that because something is valid in a mathematical sense, that it is a good thing to do on a computer. Sorry, but computer programming is NOT mathematics. Floating point arithmetic is a different animal, and is only vaguely related to mathematics.
The simple answer is to just use rank. If the matrix is less than full rank, then it is singular.
If you want, you can use cond. Large condition numbers are bad. A condition number that is roughly 1/eps indicates singularity. The problem with cond is it takes some (minor) judgment on your part. How close does a condition number need to be to 1/eps for the matrix to be numerically singular? At the same time, use of cond gives you some information that rank does not, since it tells you how close to singularity is the matrix.
And of course, you can use svd, which tells you more yet. Since the condition number is the ratio of the largest to smallest singular values, you can see if the matrix is singular easily enough. But several singular values that are relatively zero compared to the largest tells you much about the matrix too.
Again, rank is simplest here, FAR better than det. Avoid det, like you would avoid the bubonic plague. I showed several cases where det simply fails to yield an answer that is even close to truth.
Feng Cheng Chang (view profile)
Thanks for your length comment. I accept wholly your criticism.
The main object of this code 'M_singular' was to get computed square matrix M singular, and in fact it did. There is no need to prove that the matrix M is
singular by saying abs(det(M)) < 1e-8, therefore the last two line of the code should be avoided and removed.
It is pretty poor to show that using 'det(A)=0' or 'inv(A)=inf' to see whether the computed matrix A is singular or not. However, I still like to see if there are any other simple ways to get it.
Anyway, the original code is slightly revised by adding related output arguments:
function [M,A,H,s,D] = M_singular
% Making a non-singular square matrix singular
% by perturbing with prescribed distribution
% F C Chang 04/21/2013 updated 11/17/14
A = input('Enter A = ');
DetA = det(A);
AI=inv(A);
n = size(A,1);
r = input('Enter r = ');
c = input('Enter c = ');
H = r'*c;
s = -1/sum(sum(H.*AI'));
D = H*s;
M = A+D;
% DetM = det(M); % DetM,
% if abs(DetM) < 1.e-8, disp("The perturbed matrix M=A+D is singular !"), end;
For your reference, the last example A=rand(5) in your comment is applyed here to run the routine as follows:
>> [M,A,H,s,D] = M_singular;
Enter A = rand(5)
Enter r = ones(1,5)
Enter c = ones(1,5)
>> A, detA=det(A)
A =
0.98081 0.31531 0.11638 0.81645 0.99564
0.51656 0.3867 0.38323 0.74401 0.41795
0.25996 0.31241 0.33713 0.7494 0.061339
0.10318 0.58642 0.098538 0.94941 0.53959
0.53426 0.4308 0.025338 0.68107 0.81786
detA = -0.0034162
>> M, detM=det(M)
M =
0.51217 -0.15333 -0.35227 0.34781 0.527
0.047915 -0.081951 -0.08542 0.27536 -0.050693
-0.20869 -0.15624 -0.13151 0.28075 -0.40731
-0.36547 0.11778 -0.37011 0.48076 0.070939
0.065616 -0.037845 -0.44331 0.21242 0.34921
detM = -3.2035e-017
The above is for the final check, poorly apply the MATLAB build-in file 'det'.
Anyway, I still do not know how to check the matrix singular by applying 'norm' or 'cond', etc. I would like to know further.
John D'Errico (view profile)
Why do people use det to determine if a matrix is singular? For example, in this code, we see the terribly poor code...
DetA = det(A);
if abs(DetA) < 1.e-8, disp("*** The original matrix [A] is singular ! ***"), return; end;
So why is it poor?
Consider the matrix
A = eye(10);
Clearly, it is not singular. How about the related matrix B though?
B = 0.1*A;
Does everyone agree that B is not singular? Yet see that D has a small determinant!
det(A)
ans =
1
det(B)
ans =
1e-10
Or, how about this simple matrix?
A = randn(100,99);
A(:,100) = A(:,99) + eps;
It seems pretty clear that A is singular. But does det tell us that?
det(A)
ans =
6.1242e+62
Well, if the test in the code by this author is any indication, that matrix is definitely NOT singular! Of course, we can use the same scaling trick to decide it is singular.
det(A*0.1)
ans =
-1.0936e-38
The point is, det is arguably a terrible way to determine if a matrix is singular!
Instead, use rank, or cond or svd to determine that characteristic of a matrix. Of course, cond and svd tell you more about the matrix, but you also need to understand how to use those tools to infer if the matrix is singular.
The rest of the code itself is also not terribly good. It is a function, yet it uses input to pass in the matrix! I shiver at the idea, as this is a terrible user interface. A function has arguments. Learn to use them!
As I look at the code itself, I see that it will fail to execute. It fails to use valid MATLAB syntax. Here is the line I saw:
disp("*** The original matrix [A] is singular ! ***")
See that MATLAB strings are delimited by single quotes, i.e.: '
Thus we define a string as 'This is a string'
However, "This is not a string"
Next, the code says that it can be used to modify a single element, a row, a column, or the entire array, yet no explanation of how one would do any of those things. (Yes, I know after reading through the code, one might infer that these things are controlled by the vectors r and c. However, unless you understand the linear algebra involved, then one would never know what is meant there.)
An interesting question is if the perturbation generated is a minimum norm perturbation. The answer is it is not so at all.
For example, if I generate a random 5x5 matrix A as
A = rand(5);
It would seem that to choose a modification of the matrix A that allows all elements to change, then I want to specify r and c as:
r = ones(1,5);
c = ones(1,5);
Then the matrix perturbation provided by this code is given in the code as D. How much of a perturbation is D?
norm(D)
ans =
2.1399
It is pretty big actually. Instead, we can compare it to a minimum norm perturbation based on the SVD.
n = size(A,1);
[U,S,V] = svd(A,0);
Dmin = U(:,n)*S(n,n)*V(:,n)';
norm(Dmin)
ans =
0.049
And in fact the svd based solution has a larger condition number here, so it is closer to being numerically singular than that provided by the code from the author.
cond(A - Dmin)
ans =
3.5015e+16
cond(A - D)
ans =
8.2915e+15
When I originally looked at this code, I was hoping to see some nice, well written, numerically stable code to solve a moderately interesting problem. I did not find that here. I'm sorry, but there are too many problems to call this code good given all of these issues.