File Exchange

Inverse and determinant of square matrix

version 1.5 (1.44 KB) by

Inverse and determinant of a square matrix are determined using only simple matrix multiplication

Updated

The inverse (AI) and determinant (det) of a given square matrix (AO) may be directly found by
[AI,det] = inv1(AO)
It uses automatic pivoting scheme. All computations involves only simple matrix multiplication.

The direct result without pivoting may also be found by
AI = inv0(AO)
The sourse code is only 4 statement lines. Yet it works for AO = randn(n), even n = 1000.
However, it fails for some simple peculiar matrix, such as AO = [0 1; 1 0].

hüseyin çilsalar

Feng Cheng Chang

Feng Cheng Chang (view profile)

Thank you for all of your comments, John and Jan. Afterall I have decided to continue instead of giving up. Here, I want to say thank again to the encouragement from Jan.

The presentation in File Exchange is now updated to include two files: inv0.m and inv0det_pvt.m.

(1) AI = inv0(AO) -- Find AI(inverse) directly from AO(given matrix) without any pivoting. The code is only 4 statement lines total.

(2) [AI,det] = inv0det_pvt(AO,pvt) -- Find both AI and det(determinant) from AO, with or without pivoting. The "write-in" is used for the pivoting (pvt = 2). It is only good for AO of small order and not practical at all. However it may resolve the problem when the divisor 'd' in the k-loop is zero or close to zero. Such as d=AO(1,1)=0.

The option of automatic pivoting with the preset critique is the way to proceed. I think the critique for pivoting is to select the pivoting element such that the divisor d is the absolute maximum among all computed d's. If the value of d, resulted from the selected pivoting element, is zero or close to zero, then the given matrix AO must be singular.

Anyway the derived algorithm for matrix inversion is very pratical compare to some other classical approaches, such as Gauss-Jordan Elimination. I hope our FEX users will like to use it.

I am not majored in mathematics. I am retired now and soon approach to 80. I need some help regarding the option for automatic pivoting.

Jan Simon

Jan Simon (view profile)

@John: Your point has become clear. I'm convinced, that Feng did not want to insult you, me or any other FEX user. The quality of his linear algebra function is low and you rated it three times with 1 star. I'm sure that users, who are interested in linear algebra methods, will remember, that they have heard different strategies for pivoting in their numeric lessons.
According to his profile Feng Cheng Chang is a retired professor and has a PhD in Electrical Engineering. I'm sure, that he has reasons to post this submission in its current form. His calm and friendly reactions let me think, that these reasons are neither the will to provoke nor to seduce beginners to apply poor LA methods.

In my opinion, there are too few ratings and comments in the FEX. You, John, are very active here and your ratings increase the value of the complete platform. Beside others, I personally would benefit, if you use the limited time to analyze and comment other submission instead of this one.

John D'Errico

John D'Errico (view profile)

There are several reasons why you might post code for others. One is for them to use. Tools that expand MATLAB are always good, as long as they are well done. This tool has incredibly serious problems for anyone who might want to use it, and the description of this code strongly implies that it is for use! That is how I reviewed it.

There is also a teaching use of the FEX. There others can look at your code and learn from it, learning how did you solve the problem. As a learning tool, comments are extremely valuable. Code for this purpose is useless IF it claims to teach yet has no comments that explain why you have done something and what you are doing. Good teaching code should be overflowing with helpful comments. It should use descriptive variable names.

In fact, the first case I mentioned also needs internal comments, as debugging is a forever task. Things break. Changes to MATLAB happen. You cannot simply write code and forget it in any language.

The question of pivots is an important one. When I mentioned the lack of pivots, the only thing the author did was to add the options of random pivots and user supplied pivots. Both of these options are foolish, and are an insult to the user, an insult to the student who would learn from you, and an insult to me. No rational person will ever want to use them for such a code, and teaching a student that random pivoting is a good idea is crazy.

Had the author chosen to add various other pivoting options that ARE designed to improve the stability of the problem while also including the poor schemes along with clear and exhaustive explanations of why they are terribly poor choices, I'd have had no problem of their inclusion in a teaching tool. But no such thing was done.

Jan Simon

Jan Simon (view profile)

@Feng: I understand John's arguments. I miss comments, which explain the algorithm you are using. Without them, a user has to guess what you are doing and it is nearly impossible to find bugs.
A random pivoting is a surprising technique. Actually the pivoting should increase the stability of the algorithm and avoid errors caused by zeros on the diagonal. But random pivoting will let the output change (slightly or massively depending on the condition of the matrix) in a not reproducible manner.
If you are interested in the source of INV and DET, look in the corresponding BLAS and LAPACK functions, which can be downloaded at www.netlib.org. Even TMW itself did not dare to implement such fundamental functions, but existing and well established code has been chosen. This should let you and other developers of alternative INV versions think twice.

@John: Thanks for your comments and the exhaustive testing. I'm sure that Feng did not want to occupy your time.

John D'Errico

John D'Errico (view profile)

Further testing shows the latest version of the code fails to run due to a syntax error. Must I debug your code for you? Next time test your code BEFORE you post it!

So if pvt = 1 is the method we should use for normal use (can anybody actually think that is a good idea, to use random pivots?) then surely this tool can solve for the inverse of an identity matrix?

Or, maybe not.

>> [I,D] = inv_det_pvt(eye(3),1)
Warning: Input arguments must be scalar.
> In inv_det_pvt at 23
Warning: Input arguments must be scalar.
> In inv_det_pvt at 23
ijp =
2 2
1 3
3 1
AOp =
1 0 0
0 0 1
0 1 0
AIp =
NaN NaN NaN
NaN NaN NaN
NaN NaN NaN
detp =
NaN
I =
NaN NaN NaN
NaN NaN NaN
NaN NaN NaN
D =
NaN

This is not worth any more of my time.

John D'Errico

John D'Errico (view profile)

But this code does NOT show how inv is computed in MATLAB! It does not show how det is computed! It shows how they MIGHT be computed using sub-standard algorithms for both purposes.

This is my point. By giving this tool as an example, it might convince somebody to use it, thinking you have provided something useful, when in fact they would be wildly wrong in that assumption. There will surely be someone out there who would think that your code is efficient and use it for actual work. Hey, it does TWO things at once! Even worse, because you return the determinant, they might think the determinant will actually tell them some useful information about their matrix! It is not to be trusted, as the determinant can always be confused with certain problems.

I go into great depth to explain why the determinant is problematic in the link I supplied before. Even when computed in an intelligent means, as done in MATLAB (and not so by your code) the determinant is simply a bad thing to use. By connecting it to a code that returns the inverse, somebody might assume that your code returns something useful to say about the matrix inverse.

Or, somebody might assume that just because your code claims that it does pivoting, that the pivoting was done intelligently!

Offering a choice of random pivoting is simply insane. There is no mathematical reason why anyone would use that method. To assert that its use fixes a problem is simply a lucky circumstance, as the next time it is used, it may not do so.

As far as offering user provided pivoting, that is just as insane! This forces the choice of pivots onto a user who has no knowledge of what will happen deep inside the algorithm. I note that the choice of pivots chosen by an intelligent algorithmic decision was never provided. Pivoting is applied for a purpose, at least that is so in mathematics. There is a reason why pivoting is done, to make the result stable. Your code is a farce in this respect, because rather than do what makes sense with pivoting, no attempt was made by you to fix the problem.

Finally telling people to develop their own code for these purposes is itself often a poor idea. Helping them to understand how those tools work MIGHT be meaningful, IF you did so in good ways. But you don't do so. Your code has no comments internally, helping them to understand what you are doing. Teaching code should have that. It should be incredibly well commented. Yours is not.

Feng Cheng Chang

Feng Cheng Chang (view profile)

Thank you again for your comments. Afterward I do not expect people would want to use my code at all, since the MATLAB functions 'inv' and 'det' have already been here to be used conveniently. However, it does not show how these built-in functions are developed, such as by co-factors, eigenvalues, Guanssian elimination, etc. Perpaps, this is why people like to develop their own codes.

Sometimes, people want to derive his own routine that may be better than the built-in one in some extent. For example, the routine 'polyroots' presented in Mathworks Files Exchange is
more suitable than the built-in function 'roots' for any given polynomials with many high multiple roots.

If anyone like to develop his own codes about the inverse and determinant of a given square matrix, he may do it by applying the approach used in my code. It is very siple and straightforward, and only six statement lines.

Finally, the errors have been found during repeatly running the code
with pvt = 1 option. It will yield the determinant of correct
absolute value but with different sign, either + and - . It needs to be corrected.

I may withdraw this contribution if I feel no good.

John D'Errico

John D'Errico (view profile)

Ok, so "pivoting" is now in there. But the main point of my comments is why would ANYBODY use this code? If you really desperately want the inverse of a matrix, MATLAB already has inv, an immensely faster tool than this. It already has det, again faster than this.

I suppose your argument is this gives BOTH the inverse and the determinant. But why is that important? Surely NOT to test for singularity? That is perhaps the single worst use of a matrix determinant that anyone has ever invented. I go into some detail about why the determinant is so useless here:

http://stackoverflow.com/questions/13145948/how-to-find-out-if-a-matrix-is-singular/13146750#13146750

As for computing and returning both results, this function does it far better and far faster, even though it technically does more work than necessary:

Ainv = inv(A);

Feel free to add your own help, since the function is of no value anyway. But that function will be far faster than the code provided by the author. The author brags about his short code, but mine is only 2 lines of MATLAB code.

As far as the new options go, random pivoting? Is it really a good idea to provide a function that will randomly fail or not on some problems? Or the second option, wherein the user must supply the pivot sequence? Give me a break!

I think the author does not understand the idea of why pivoting is important in these problems. Go back and check your linear algebra 101 notes.

Pivoting makes the solution more stable, IF you use intelligent pivots. But making the user supply the pivots is incredibly insane, and random pivots are a complete waste of time.

So still I see no value in the tool the author has provided. In fact, the changes have made the code worse, because someone might believe that what was done was good. Just because some arbitrary pivoting was done does not make it well written.

Feng Cheng Chang

Feng Cheng Chang (view profile)

I do appreciate John's comments, so that I am to force myself to do some further homework. Hope it will give him some answers but not expect all.

The code is now updated to include the pivoting scheme. There are three options to choose:
(1) pvt=0 ---> no pivoting,
(2) pvt=1 ---> pivoting elements randomly selected,
(3) pvt=2 ---> pivoting elements by writing-in.

If the very first element of the given square matrix is zero, it does surely fail for option(1), but it will be OK by sucessively running either option(2)or(3) for any non-singular matrix.

The code derived is very short (10 lines for the original and less than 30 for the updated). Yet it will give both inverse and determinant of a given square matrix.

After all, nothing is completely perfect.

John D'Errico

John D'Errico (view profile)

This is BAD for many reasons. For example, it completely fails on the very simple 2x2 matrix,

A = [0 1;1 0];
[I,D] = inv_det(A)
I =
NaN NaN
NaN 0
D =
NaN

That seems funny, since I could have sworn that A has an inverse. Oh, yeah. A is its own inverse. So I wonder why inv_det fails here? POOR CODE.

MATLAB has no problem with this nice matrix.

inv(A)
ans =
0 1
1 0

Is it faster than other tools in MATLAB? NO!

A = randn(1000);
tic,[I,D] = inv_det(A);toc
Elapsed time is 13.865905 seconds.

Admittedly, the determinant of such a large matrix is meaningless.

tic,I = inv(A);D = det(A);toc
Elapsed time is 0.123640 seconds.

Yeah, MATLAB itself does a better job. It is faster. It is more accurate. Why do I say more accurate? Because inv_det makes no attempt to do pivoting. This is why the simple 2x2 example above failed miserably.

So is the help any good? No.

help inv_det
Inverse and determinant of square matrix
F C Chang 10/29/2012

It tells me nothing that I need to know. Which order do the outputs arrive in? Must the input array be double? May it be single? Is int8 acceptable?

This code is basically just a homework assignment, and one that would not receive an A either. Maybe a C.

John D'Errico

Feng Cheng Chang

Feng Cheng Chang (view profile)

If the run fails due to det = 0 at any stage, please try the following added routine:

function [AI,det,er] = inv_det_pij(AO,pij)
if nargin < 2, AO = AO(pij(:,1),pij(:,2)); end;
[AO,det,er] = inv_det(AO);
if nargin < 2, AO(pij(:,2),pij(:,1)) = AO; end;

where the input 'pij' is nx2 array to modify the given matrix such that the diagonal elements are all replaced.

For example,

>> A=[0 1 2; 3 4 5; 6 7 8], pij=[2 1; 3 2; 1 3], Ap=A(pij(:,1),pij(:,2)),

A =
0 1 2
3 4 5
6 7 8
pij =
2 1
3 2
1 3
Ap =
3 4 5
6 7 8
0 1 2