# more numerically accurate inverse matrix

I encountered the following matrix

mat = {{2,
2.161209223472559 + 1.682941969615793 I}, {2.161209223472559 -
1.682941969615793 I, 2}}


and Inverse[mat] will give

{{-0.57092 - 1.06364*10^-16 I,
0.616939 + 0.480412 I}, {0.616939 - 0.480412 I, -0.57092 -
1.11022*10^-16 I}}


Notice there is small imaginary part and they ought to be zero. as you can see from general symbolic calculation.

In:= Inverse[{{a, b + I c}, {b - I c, a}}]

Out= {{1/(1 - b^2 - c^2), (-b - I c)/(
1 - b^2 - c^2)}, {(-b + I c)/(1 - b^2 - c^2), 1/(1 - b^2 - c^2)}}


Normally, those small part won't bother. But I am doing iteration calculation right now, in each iteration step, there is inversion process and I found those small part in every step will greatly affect the result after only 30 iteration.

According to this link and this link, inverse matrix should never be performed and they recommend linear equation solving. So for example

mat1={{I, -1}, {-1, -I}}
mat2={{-I, -1}, {-1, I}}


and we want to calculate

mat1.Inverse[mat].mat2


we could do it without inverse like this

mat1.LinearSolve[mat,mat2]


But I found this give exactly the same result as directly calculating Inverse

The way I can think of is to Chop matrix at every step, but I don't know whether it will accumulate other kind of error in the iteration process.

So what is the correct way to deal with inverse matrix in this case?

PS: I also found Python's numpy gives more accurate inverse than Mathematica. for example, numpy.linalg.inv(mat) gives imaginary part 6x10^-17, and np.dot(np.dot(mat1,mat),mat2) will not present imaginary part.

Update

first

As Karsten 7. has pointed out, Method -> "CofactorExpansion" will give correct result. And I think this method is actually using direct formula like the following for 2x2 matrix

$$\mathbf A^{-1}=\frac1{\det \mathbf A}\begin{bmatrix}d&-b\\-c&a\end{bmatrix}=\frac1{ad-bc}\begin{bmatrix}d&-b\\-c&a\end{bmatrix}$$

Though CofactorExpansion is slower than direct Inverse, but for 2x2 matrix, it is acceptable efficient. But I have no idea about the numerical stability.

Second

The thing I don't understand right now is that even using Mathematica's LinearAlgebraLapack(see here) still can't get the correct answer, but simple fortran coding using lapack do give the correct answer !!!!

According to Lapack,

"getrf" Computes the LU factorization of a general m-by-n matrix.

"getri" Computes the inverse of an LU-factored general matrix.

so we have the following code

tmp = mat;
ipiv = ConstantArray[1, Length@tmp];
LinearAlgebraLAPACKGETRF[tmp, ipiv, info];
LinearAlgebraLAPACKGETRI[tmp, ipiv, info];
tmp


But this gives result

{{-0.57092 - 1.06364*10^-16 I,
0.616939 + 0.480412 I}, {0.616939 - 0.480412 I, -0.57092 -
1.11022*10^-16 I}}


exactly the same wrong result as direct Inverse!!

But I have tried fortran coding below

program testinversion
use lapack95
use f95_precision
implicit none
complex*16,dimension(2,2)::a
integer,dimension(2)::ipiv
integer::info
a=reshape((/(2.,0.),(2.161209223472559,-1.682941969615793),(2.161209223472559,1.682941969615793),(2.,0.)/),(/2,2/))
call zgetrf(2, 2, a, 2, ipiv, info)
if(info==0) then
call getri(a,ipiv,info)
print*,a
else
print*,"error"
endif
end program testinversion


Compile it with

ifort testinversion.f90 -mkl=sequential -lmkl_blas95_lp64 -lmkl_lapack95_lp64


it gives the correct result! What is wrong with Mathematica's LAPACK?

update2

Solving inverse of 2x2 mat is actually equivalent to solving $$mat \cdot x = \left( {\begin{array}{*{20}{c}} 1&0\\ 0&1 \end{array}} \right)$$

But LinearSolve[mat,{1,0}] gives

{-0.57092 - 5.91873*10^-17 I, 0.616939 - 0.480412 I}


looking at the document of LinearSolve, I found more things: "CofactorExpansion","DivisionFreeRowReduction","OneStepRowReduction" are actually designed for symbolic solving. For numerical solving, There are "Banded" "Cholesky","Krylov","Multifrontal".

For this paticular simple problm, mathematica choosed a wrong method??!!

As I tested, both "Krylov" and "Multifrontal" can give correct answer.

So now, I don't know whether Inverse is using LinearSolve, but apparently, automatic method in LinearSolve is also bugged.

• (1) Version 10.3 on Linux here is what I see for those products. In:= mat1.LinearSolve[mat, mat2] Out= {{-0.181014708395 + 0. I, 0. + 0.181014708395 I}, {0. - 0.181014708395 I, -0.181014708395 + 0. I}} In:= mat1.Inverse[mat].mat2 Out= {{-0.181014708395 + 0. I, 0. + 0.181014708395 I}, {0. - 0.181014708395 I, -0.181014708395 + 0. I}}. For the inverse itself, I get In:= Inverse[mat] Out= {{-0.570919803469 + 6.01693250938*10^-17 I, 0.61693857256 + 0.480412449271 I}, {0.61693857256 - 0.480412449271 I, -0.570919803469 + 5.55111512313*10^-17 I}}. Dec 17 '15 at 16:58
• (2) This is almost certainly BLAS and OS dependent. Dec 17 '15 at 16:58
• @DanielLichtblau You are right. The result of imaginary part is OS dependent. But that is not the whole story, see my update! Dec 17 '15 at 17:04
• Just to be sure. You aren't calculating the Inverse in order to solve a linear system, are you? Dec 17 '15 at 20:44
• "I have to multiply this inverse matrix with other matrix" - then yes, you are solving a linear system. As to why LinearSolve[] is acting like this for a manifestly Hermitian system, I don't yet have any ideas. Dec 18 '15 at 10:38

Setting the Method option to "CofactorExpansion" results in the correct output.

mat = {{2, 2.161209223472559 + 1.682941969615793 I},
{2.161209223472559 - 1.682941969615793 I, 2}}

Inverse[mat, Method -> "CofactorExpansion"]


$\$ {{-0.57092 + 0. I, 0.616939 + 0.480412 I}, {0.616939 - 0.480412 I, -0.57092 + 0. I}}

As you want to perform iterative calculations, it might be in general a good idea to Rationalize mat.

Inverse[Rationalize[mat, 0]] N @ %


$\$ {{-0.57092, 0.616939 + 0.480412 I}, {0.616939 - 0.480412 I, -0.57092}}

One can also use SetPrecision as suggested in a comment by J. M.♦, but needs to set it to Infinity.

Inverse[SetPrecision[mat, ∞]] N @ %


$\$ {{-0.57092, 0.616939 + 0.480412 I}, {0.616939 - 0.480412 I, -0.57092}}

• Thank you. Rationalize is not suitable for me. Because I have to repeat iteration for thousands of times and rationalize is too slow. And for Method -> "CofactorExpansion", it indeed works. But on what circumstances that we should use "CofactorExpansion"? How can we be sure that it will be suitable in later iteration process?Why mma doesn't automatically choose the right way? Dec 17 '15 at 12:20
• Oh, I just read wikiwand.com/en/Invertible_matrix#/Eigen_decomposition . It seems that cofactor expansion is the same thing as the symbolic inversion formula and should be effective only for small matrix. But numerical instability is not mentioned in wikipedia Dec 17 '15 at 12:47
• @matheorem With using Rationalize I didn't mean to Rationalize mat every time you use Inverse, but to Rationalize it just one. If that is possible and makes sense depends on the rest of your iteration process. Dec 17 '15 at 14:57
• @matheorem With respect to why Mma doesn't choose "CofactorExpansion", I can only cite the documentation: "For matrices with approximate real or complex numbers, the inverse is generated to the maximum possible precision given the input." "The default setting of Automatic switches among these methods depending on the matrix given." That means that the default isn't designed for your special situation, where you want to have a result with a higher precision than the input. Another reason not to use "CofactorExpansion" is that it is slower than Automatic for bigger matrices. Dec 17 '15 at 15:04
• what does maximum possible precision mean? btw, i updated my post Dec 17 '15 at 17:20

I agree the imaginary parts should be zero. I do not know why they are not zero. But this is what I found, too small to put in comment.

First, Matlab does give zero for the exact same input:

format long g
mat = [2, 2.161209223472559 + 1.682941969615793*1j;
2.161209223472559 - 1.682941969615793*1j, 2]

inv(mat)
-0.570919803469126 + 0i                   0.616938572560308 + 0.480412449271497i
0.616938572560308 - 0.480412449271497i    -0.570919803469126+0i


You can get same output in Mathematica by doing the direct computation itself without calling Inverse

foo[x11_, x12_, x21_, x22_] :=
Module[{}, {{x22, -x12}, {-x21, x11}}/(x11*x22 - x21*x12)]
foo[2, 2.161209223472559 + 1.682941969615793 I,
2.161209223472559 - 1.682941969615793 I, 2] It is an exact zero for the complex part: And can see it match the Matlab output. Even compiled version did not resolve the issue (even though told it to run in hardware floating point)

cf = Compile[{{x11, _Complex},{x12, _Complex},{x21, _Complex},
{x22,_Complex}},
Inverse[{{x22, -x12}, {-x21, x11}}], RuntimeOptions -> "Speed"
];

cf[2, 2.161209223472559^1 + 1.682941969615793^1 I,
2.161209223472559^1 - 1.682941969615793^1 I, 2] If we do break the inverse to 2 parts, and do one by 'hand' and then use Det only, then now the accuracy improves a little, and now it is of order 10^-17

foo1[x11_, x12_, x21_, x22_] :=
Module[{}, {{x22, -x12}, {-x21, x11}}/Det[{{x11, x12}, {x21, x22}}]]
foo1[2, 2.161209223472559 + 1.682941969615793 I,
2.161209223472559 - 1.682941969615793 I, 2] Here is Maple 2015 result also:

mat:=<<2|2.161209223472559 + 1.682941969615793*I>,
<2.161209223472559 - 1.682941969615793*I|2>>;
LinearAlgebra[MatrixInverse](mat); • So matlab and maple, no matter how they done it internally, they done it right. Can we say that Inverse in mma has some defects? But I think for numerical inversion, they should be all implementing BLAS or Lapack, right? Why such a difference? Dec 17 '15 at 12:25
• Possibly Matlab and Mathematica are getting the same result. It seems Mathematica shows one more decimal place (I don't know what the Matlab result would look like if it showed one more place), In:= Inverse[mat] // InputForm Out//InputForm= {{-0.5709198034691264 + 6.016932509384263*^-17*I, 0.6169385725603085 + 0.4804124492714965*I}, {0.6169385725603084 - 0.4804124492714966*I, -0.5709198034691264 + 5.551115123125783*^-17*I}} Dec 17 '15 at 17:10
• i updated my post. hope you have a look Dec 17 '15 at 17:21