How can I derive inverse of the matrix with infinite determinant?

Hi guys, there is a matrix A which has finite dimension and finite-valued elements.

However, because of its large size and values, the determinant of A becomes infinite on MATLAB.

load("myMatrix.mat","A")
size(A)

ans = 1×2

320 320

max(max(abs(A)))

ans = 2.2913e+04

det(A)

ans = Inf

Therefore, it cannot compute the inverse of the matrix precisely.

B = inv(A);

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.336932e-19.

isdiag(A*B)

ans = logical
   0

Is there way to compute the inverse of the matrix A?

Cheers,

Answer 1

[Bruno Luong]

Welcome to the world of numerical calculation. Every conclusion you made is wrong.

"det(A) = Inf Therefore, it cannot compute the inverse of the matrix precisely. ".

Wrong

A=eye(200)*1024;
det(A)

ans = Inf

B=inv(A);
all(A*B == eye(size(A)), 'all') % A*B is exatly I, despite det(A) is Inf due to overflow

ans = logical
   1

Check B is inverse of A using isdiag(A*B).

This is a bad numerically criteria as showed in the false negative answer example

A=[1 2; 3 4]

A = 2×2

1 2 3 4

B = inv(A)

B = 2×2

-2.0000 1.0000 1.5000 -0.5000

AB = A*B;
isdiag(AB) % Suppose to be TRUE

ans = logical
   0

AB(2,1) % Suppose to be 0, but not exactly

ans = 8.8818e-16

[Junho Kweon]

You are right. I didn't consider the computational limitation. Haha :)

[Bruno Luong]

The funny thing is the inconsistency of linear algebra - mostly det - in floating point world.

The matrix is numerical singular, so the determinant should be 0, yet it is Inf.

[Junho Kweon]

@Bruno Luong Yeah, it is quite interesting that detmerninant is not zero though it is mathematically 0. May be in the mathematical computation of determinant, (very large number)- (very large number), and it is over the computational limitation of MATLAB. Therefore, it might be Inf comes out, and MATLAB just conclude that the determinant is infinite.

Answer 2

[Torsten]

What do you get if you type

rank(A)

?

If it's smaller than 320, no inverse exists (or is almost impossible to calculate).

[Junho Kweon]

I was silly not to check that one. Thanks

Answer 3

[Bruno Luong]

"Is there way to compute the inverse of the matrix A?"

No.

The sum of all columns of your matrix A is numerically 0, therefore your matrix is not invertible.

[Steven Lord]

The first of Cleve's Golden Rules of computation applies here. "The hardest things to compute are things that do not exist."

[Bruno Luong]

Hmm read these two rules from Cleve I wonder this: things that do not exist are they unique?

It sounds like a non-sense question, but why it is non-sense?

If it does make sense I have no preference for the answer.

Answer 4

[Walter Roberson]

In theory if you have the symbolic toolbox, you could

digits(16);
sA = sym(A, 'd');
digits(32);
Ainv = inv(sA);

... But test it with a much smaller matrix first, as it might take very very long.

[John D'Errico]

The symbolic inverse of a 320x320 matrix? Unless the matrix is REALLY, simple, that could take the lifetime of the universe to compute.

[Bruno Luong]

I don't know symbolic at all but is there any approximation going on as well (digits=32)? Why it does give out Ainv for the matrix that has numerical singular?

How to exploit the bunch of fractions in the result and be sure it is correct with still a finite truncation (digits=32)?

It doesn't look like it is any better than pure numerical method.

A = rand(8,7)*rand(7,8)

A = 8×8

1.1932 1.4076 1.2416 0.7862 1.4254 1.3012 1.7597 1.5144 1.6909 2.2753 2.1156 1.4951 2.5238 2.3885 2.7331 2.1692 1.1620 1.6231 1.5313 1.2976 1.9202 1.8231 2.1256 1.5481 1.1399 1.2803 1.2986 0.9485 1.6928 1.6402 2.0115 1.3554 1.7376 2.4798 2.3184 1.5994 2.6161 2.4092 2.8091 2.1660 1.3393 1.5057 1.3630 0.9368 1.4926 1.5373 2.0491 2.0210 1.3573 2.2611 1.8506 1.2858 2.4179 2.2090 2.5930 2.0655 1.8750 1.6381 2.0915 1.2237 2.0587 2.2006 2.3607 1.8280

B = A\eye(size(A)) % for those who think \ is better than B=inv(A)

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.262494e-18.

B = 8×8

1.0e+15 * 0.3034 -0.4816 0.3936 -0.2782 -0.1810 -0.1330 0.3088 0.2056 0.7213 -1.1449 0.9358 -0.6614 -0.4303 -0.3162 0.7342 0.4889 -0.5342 0.8479 -0.6930 0.4898 0.3187 0.2342 -0.5437 -0.3620 0.0684 -0.1086 0.0888 -0.0628 -0.0408 -0.0300 0.0697 0.0464 -1.0208 1.6202 -1.3242 0.9360 0.6090 0.4475 -1.0390 -0.6918 0.8036 -1.2754 1.0425 -0.7368 -0.4794 -0.3523 0.8179 0.5446 0.2064 -0.3277 0.2678 -0.1893 -0.1232 -0.0905 0.2101 0.1399 -0.4766 0.7565 -0.6183 0.4370 0.2843 0.2090 -0.4852 -0.3230

Inum = A*B % Garbage

Inum = 8×8

0.8339 0.1022 -0.4563 -0.0586 0.1096 0.0810 -0.1324 -0.0255 -0.2712 1.1957 -0.2699 0.0979 0.1207 0.1741 0.0739 -0.0912 0.1829 -0.1298 0.8032 -0.1076 0.0034 0.0307 -0.1420 0.1355 -0.2903 0.1039 -0.6101 0.9326 0.0593 0.0702 -0.0003 -0.0164 -0.0857 -0.2977 0.0205 -0.0615 1.0466 0.1174 0.2904 0.2597 -0.0443 0.1953 -0.3567 0.0948 0.0285 1.1345 0.1275 0.0570 -0.0902 -0.1424 0.0200 -0.3266 -0.1365 0.1368 1.1092 0.1673 0.0828 -0.2810 -0.1661 -0.1480 0.0341 -0.0215 0.1781 1.1653

norm(Inum-eye(size(A)))

ans = 1.0418

digits(16);
sA = sym(A, 'd');
digits(32);
Ainv = inv(sA)

Ainv = 

Isym = A*double(Ainv) % Still garbage, is it a correct test?

Isym = 8×8

1.1474 -0.4078 0.4030 -0.2153 -0.2193 -0.1139 0.1195 0.1647 -0.2932 1.1517 0.0988 -0.0115 0.0175 0.0423 -0.1806 -0.0345 -0.7159 1.2818 0.1346 0.6660 0.4150 0.2834 -0.8288 -0.4359 -0.0244 -0.1726 0.1465 0.9150 -0.0716 -0.0160 0.0375 -0.0374 1.0789 -1.5676 1.3428 -0.9192 0.3049 -0.4418 0.9138 0.6292 0.4085 -0.6508 0.6104 -0.2542 -0.2783 0.8002 0.5759 0.3289 -0.3197 0.3399 -0.5910 0.4714 0.3066 0.2587 0.4518 -0.2483 -0.6429 0.9187 -0.6401 0.6024 0.2725 0.2275 -0.5178 0.5798

norm(Isym-eye(size(A)))

ans = 4.2935

[Walter Roberson]

sorry should have been 'd' flag so that it creates symbolic floating point values instead of symbolic fractions.

It is not documented whether the symbolic engine is able use an algorithm similar to numeric inv() in the case where all of the entries are symbolic floating point. If it has to resort to minors then certainly it would take aeons for a matrix that long.

[Bruno Luong]

OK I make change accordingly to 'd' flag in my test code, and as you can see the symbolic is no bettrer than numerical method with "\", at least in this toy example.

Even if it use minor formula, there are a lot of cancellation going on in determinant so the summation order and tree (parenthesis) matter. The symbolic probably does the sumation in the implicit order of the implementation as it progres and no care about cancellation.

I'm not sure if the symbolic solution is any better than numerical method for precision point of view, let alone the computation time that goes as factorial of the dimesion.

[Walter Roberson]

The point about this approach was that it should avoid the intermediate inf