Very different matrix inverse result in complex ring CC and in real ring RR while entry only contain real elements

asked 2020-12-10 01:11:16 +0100

I would use an dramatized example here but it was encountered during a project as well.

Suppose two matrix consisted of the same entry of only real elements. However, one was in the real domain RR, and one was in the complex domain CC. The inverse of the complex matrix and the real matrix different significantly.

Example:

test_S=matrix(RR,5)
for ix in range(0,5):
    for iy in range(0,5):
        test_S[ix,iy]=sin(ix+iy)*sin(iy+1);

and

test_S_CC=matrix(CC,5)
for ix in range(0,5):
    for iy in range(0,5):
        test_S_CC[ix,iy]=sin(ix+iy)*sin(iy+1);

which returned

test_S

[  0.000000000000000   0.765147401234293   0.128320060202457  -0.106799974237582   0.725716283876408]
[  0.708073418273571   0.826821810431806  0.0199148566748170   0.572750016904307   0.919535764538226]
[  0.765147401234293   0.128320060202457  -0.106799974237582   0.725716283876408   0.267938303940044]
[  0.118748392158235  -0.688158561598754  -0.135323401369264   0.211462346264655  -0.630000397639817]
[ -0.636827341031836  -0.871947375471875 -0.0394311173578842  -0.497209097294248  -0.948719639025321]

test_S_CC

[  0.000000000000000   0.765147401234293   0.128320060202457  -0.106799974237582   0.725716283876408]
[  0.708073418273571   0.826821810431806  0.0199148566748170   0.572750016904307   0.919535764538226]
[  0.765147401234293   0.128320060202457  -0.106799974237582   0.725716283876408   0.267938303940044]
[  0.118748392158235  -0.688158561598754  -0.135323401369264   0.211462346264655  -0.630000397639817]
[ -0.636827341031836  -0.871947375471875 -0.0394311173578842  -0.497209097294248  -0.948719639025321]

However, their inverse

(test_S)^-1

[ 2.87043424444896e15 -9.36229301025494e15  1.13569572695379e16 -8.12065238284797e15  1.72141024999095e15]
[-1.10320559226948e15  3.19497062929143e15  5.27553667803348e15 -7.22878751888619e15  8.54302211670640e15]
[ 1.94664421070014e16 -3.02112860544102e15    0.000000000000000  1.80143985094820e16    0.000000000000000]
[ 8.29739198582516e14  9.39732200089164e15 -1.02939420054183e16  1.02939420054183e16    0.000000000000000]
[-2.15677123538261e15 -1.45141888820365e15 -7.07708512872506e15  5.95118522188244e15 -9.00719925474099e15]

(test_S_CC)^-1

[-1.00798918109020e16  1.87532050063097e15  1.07471168476593e16 -2.09935440833922e16  1.10831351979113e16]
[ 2.64758058163658e15  7.44366790239271e15 -8.04492396917722e15  9.82621377301499e15  4.42753908914828e14]
[ 2.35701470124050e16 -4.71646526274252e15  1.02939420054183e16  1.08658276723860e16  9.15017067148291e15]
[ 1.53627372015977e16 -1.32579194813951e15 -1.02939420054183e16  2.51629693465780e16 -9.15017067148291e15]
[-4.69821640166124e15 -7.20926442861200e15  5.14697100270914e15 -8.57828500451523e15 -3.43131400180609e15]

where

(test_S)^(-1)-test_S_CC^(-1)

[ 1.29503260553510e16 -1.12376135108859e16  6.09840421878568e14  1.28728917005442e16 -9.36172494792037e15]
[-3.75078617390606e15 -4.24869727310128e15  1.33204606472107e16 -1.70550012919012e16  8.10026820779157e15]
[-4.10370490540364e15  1.69533665730150e15 -1.02939420054183e16  7.14857083709603e15 -9.15017067148291e15]
[-1.45329980030152e16  1.07231139490312e16    0.000000000000000 -1.48690273411597e16  9.15017067148291e15]
[ 2.54144516627863e15  5.75784554040835e15 -1.22240561314342e16  1.45294702263977e16 -5.57588525293490e15]

Is this a bug? How to fix the issue?

edit retag flag offensive close merge delete

add a comment

3 Answers

Sort by » oldest newest most voted

answered 2020-12-10 01:56:51 +0100

Florentin Jaffredo
636 ●2 ●13 ●29

updated 2020-12-10 02:00:31 +0100

That's not a bug, your matrix is simply not invertible. You can see it by switching the ring to SR to perform exact computations:

test_S_SR = matrix(SR, 5)
for ix in range(0, 5):
    for iy in range(0, 5):
        test_S_SR[ix, iy] = sin(ix+iy)*sin(iy+1)

sage: test_S_SR

[            0 sin(2)*sin(1) sin(3)*sin(2) sin(4)*sin(3) sin(5)*sin(4)]
[     sin(1)^2      sin(2)^2      sin(3)^2      sin(4)^2      sin(5)^2]
[sin(2)*sin(1) sin(3)*sin(2) sin(4)*sin(3) sin(5)*sin(4) sin(6)*sin(5)]
[sin(3)*sin(1) sin(4)*sin(2) sin(5)*sin(3) sin(6)*sin(4) sin(7)*sin(5)]
[sin(4)*sin(1) sin(5)*sin(2) sin(6)*sin(3) sin(7)*sin(4) sin(8)*sin(5)]

sage: test_S_SR.determinant().trig_reduce() # can take a few seconds
0

The float version looks invertible but that's only due to machine precision. Since the internal mechanism are probably different for RR and CC, it's normal that you get different answers.

edit flag offensive delete link

Comments

But this happened to invertible matrix as well. The above result was only for demonstration, since it was calculated with float precision. I encountered this in a project where an invertible matrix calculated in CC had a very different inverse than the one calculated with RR. And a final result ought to be 0 but returned 0.0354399342695375 - 0.109072902259946*I.

One observation was that the invertible matrix calculated in CC was almost 1e2 times less preciece than the standard numpy matrix. After I used RR and complete the matrix inverse by directly calling the inverse function in numpy package, the result then returned correctly. There was something very wrong with the precision control of matrix inverse in CC.

ShoutOutAndCalculate ( 2020-12-10 04:27:07 +0100 )edit

In CC

from numpy.linalg import inv as np_inv
matrix(np_inv(rho_num_S))*rho_num_S

returned 0.9999999999999999 and 0.9999999999999999 for one of the entry where

rho_num_S^-1*rho_num_S

returned 1.07087986853908 and 0.723606797749979 for two of the entry.

where in RR both expression returned 1.00000... I don't have enough space to copy the function but you can try couple of square matrix of higher order.

For the interest in our project, default matrix inverse operation in CC worked for square matrix of 3 by 3, but failed for square matrix of 6 by 6, and default matrix inverse operation in RR worked for square matrix 6 by 6 but failed for square matrix of 10 by 10. But so far, the inverse function provided by numpy package had been working for all those cases.

ShoutOutAndCalculate ( 2020-12-10 04:34:57 +0100 )edit

add a comment

answered 2020-12-10 12:46:43 +0100

Emmanuel Charpentier
7774 ●6 ●51 ●154

WorksForMe (on 9.3.beta3) :

sage: MR=MatrixSpace(RR,3).random_element()
sage: MR
[-0.0413379232879041  -0.383726107921675   0.402484964228335]
[ 0.0264553329794861   0.345098500636929   0.264737280303141]
[ -0.141479831307881  -0.998284842827106   0.178457190865074]
sage: MR.det()
0.0117348614307324
sage: MR^-1
[  27.7692689578674  -28.4039106845136  -20.4930893556867]
[ -3.59408846998909   4.22386369616762   1.83995041270540]
[  1.91007962570815   1.10972616436864 -0.350583891628722]
sage: MR*MR^-1
[     1.00000000000000 -4.99600361081320e-16  1.11022302462516e-16]
[-2.22044604925031e-16      1.00000000000000 -8.32667268468867e-17]
[-2.77555756156289e-16 -2.77555756156289e-17      1.00000000000000]
sage: MC=MR.change_ring(CC)
sage: MC
[-0.0413379232879041  -0.383726107921675   0.402484964228335]
[ 0.0264553329794861   0.345098500636929   0.264737280303141]
[ -0.141479831307881  -0.998284842827106   0.178457190865074]
sage: MC.det()
0.0117348614307324
sage: MC^-1
[  27.7692689578674  -28.4039106845136  -20.4930893556867]
[ -3.59408846998909   4.22386369616762   1.83995041270540]
[  1.91007962570815   1.10972616436864 -0.350583891628722]
sage: (MC^-1)-(MR^-1)
[0.000000000000000 0.000000000000000 0.000000000000000]
[0.000000000000000 0.000000000000000 0.000000000000000]
[0.000000000000000 0.000000000000000 0.000000000000000]

Of course, that's not a proof, just an example. If you have a counter-example, I'd be glad to see it...

IMHO, Florentin is right : any deviation has probably to be attributed to imprecisions in floating point approximate computations.

HTH,