16
$\begingroup$

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[494]:= Inverse[{{a, b + I c}, {b - I c, a}}]

Out[494]= {{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 LinearAlgebra``Lapack(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];
LinearAlgebra`LAPACK`GETRF[tmp, ipiv, info];
LinearAlgebra`LAPACK`GETRI[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!

enter image description here

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.

Improve this question
$\endgroup$
14
  • 1
    $\begingroup$ (1) Version 10.3 on Linux here is what I see for those products. In[59]:= mat1.LinearSolve[mat, mat2] Out[59]= {{-0.181014708395 + 0. I, 0. + 0.181014708395 I}, {0. - 0.181014708395 I, -0.181014708395 + 0. I}} In[60]:= mat1.Inverse[mat].mat2 Out[60]= {{-0.181014708395 + 0. I, 0. + 0.181014708395 I}, {0. - 0.181014708395 I, -0.181014708395 + 0. I}}. For the inverse itself, I get In[37]:= Inverse[mat] Out[37]= {{-0.570919803469 + 6.01693250938*10^-17 I, 0.61693857256 + 0.480412449271 I}, {0.61693857256 - 0.480412449271 I, -0.570919803469 + 5.55111512313*10^-17 I}}. $\endgroup$
    – Daniel Lichtblau
    Commented Dec 17, 2015 at 16:58
  • 1
    $\begingroup$ (2) This is almost certainly BLAS and OS dependent. $\endgroup$
    – Daniel Lichtblau
    Commented Dec 17, 2015 at 16:58
  • 1
    $\begingroup$ @DanielLichtblau You are right. The result of imaginary part is OS dependent. But that is not the whole story, see my update! $\endgroup$
    – matheorem
    Commented Dec 17, 2015 at 17:04
  • 1
    $\begingroup$ Just to be sure. You aren't calculating the Inverse in order to solve a linear system, are you? $\endgroup$
    – Karsten7
    Commented Dec 17, 2015 at 20:44
  • 2
    $\begingroup$ "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. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Dec 18, 2015 at 10:38

2 Answers 2

Reset to default
13
$\begingroup$

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]]

out

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, ∞]]

out

N @ %

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

Improve this answer
$\endgroup$
5
  • $\begingroup$ 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? $\endgroup$
    – matheorem
    Commented Dec 17, 2015 at 12:20
  • 1
    $\begingroup$ 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 $\endgroup$
    – matheorem
    Commented Dec 17, 2015 at 12:47
  • $\begingroup$ @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. $\endgroup$
    – Karsten7
    Commented Dec 17, 2015 at 14:57
  • 1
    $\begingroup$ @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. $\endgroup$
    – Karsten7
    Commented Dec 17, 2015 at 15:04
  • $\begingroup$ what does maximum possible precision mean? btw, i updated my post $\endgroup$
    – matheorem
    Commented Dec 17, 2015 at 17:20
7
$\begingroup$

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]

Mathematica graphics

It is an exact zero for the complex part:

Mathematica graphics

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]

Mathematica graphics

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]

Mathematica graphics

Here is Maple 2015 result also:

mat:=<<2|2.161209223472559 + 1.682941969615793*I>,
      <2.161209223472559 - 1.682941969615793*I|2>>;
LinearAlgebra[MatrixInverse](mat);

Mathematica graphics

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ 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? $\endgroup$
    – matheorem
    Commented Dec 17, 2015 at 12:25
  • $\begingroup$ 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[68]:= Inverse[mat] // InputForm Out[68]//InputForm= {{-0.5709198034691264 + 6.016932509384263*^-17*I, 0.6169385725603085 + 0.4804124492714965*I}, {0.6169385725603084 - 0.4804124492714966*I, -0.5709198034691264 + 5.551115123125783*^-17*I}} $\endgroup$
    – Daniel Lichtblau
    Commented Dec 17, 2015 at 17:10
  • $\begingroup$ i updated my post. hope you have a look $\endgroup$
    – matheorem
    Commented Dec 17, 2015 at 17:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.