inv_det_p(M)
Finding inverse and determinant of matrix by order expansion and condensation
Author: Feng Cheng Chang

The routine [iM,dM,p]=inv_det_p(M) is very compact, and involves only simple arithmatic multiplication/division operations. The total number of operations is about N^3, which is just required for the product of any two NxN matrices. The given matrix M may be real or complex, and it works without failure even for very large matrices, such as M=magic(55)+1i*ones(55,55).

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

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)

In the routine [iA,dA,p] = inv_det_0(A), the locations of array p are fixed along the diagonal of the given matrix A. This routine will fail if A is a singular matrix or P(k)->0 at the k-th iteration step.
Therefore, if p(k) is picked at absolute maximum among possible locations, the routine will not fail. I will present it soon.

The routine [iM,dM,p]=inv_det_p(M) is very compact, and involves only simple arithmatic multiplication/division operations. The total number of operations is about N^3, which is just required for the product of any two NxN matrices. The given matrix M may be real or complex, and it works without failure even for very large matrices, such as M=magic(55)+1i*ones(55,55).

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

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.

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)

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.

The routine [iM,dM,p]=inv_det_p(M) is very compact, and involves only simple arithmatic multiplication/division operations. The total number of operations is about N^3, which is just required for the product of any two NxN matrices. The given matrix M may be real or complex, and it works without failure even for very large matrices, such as M=magic(55)+1i*ones(55,55).

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

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.

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.

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.

Comment only