Alright, so I have 3 numpy matrices :
m1 = [[ 3 2 2 ... 2 2 3]
[ 3 2 2 ... 3 3 2]
[500 501 502 ... 625 626 627]
...
[623 624 625 ... 748 749 750]
[624 625 626 ... 749 750 751]
[625 626 627 ... 750 751 752]]
m2 = [[ 3 2 500 ... 623 624 625]
[ 3 2 500 ... 623 624 625]
[ 2 3 500 ... 623 624 625]
...
[ 2 2 500 ... 623 624 625]
[ 2 2 500 ... 623 624 625]
[ 3 2 500 ... 623 624 625]]
m3 = [[ 813 827 160500 ... 199983 200304 200625]
[ 830 843 164000 ... 204344 204672 205000]
[ 181317 185400 36064000 ... 44935744 45007872 45080000]
...
[ 221046 225867 43936000 ... 54744256 54832128 54920000]
[ 221369 226196 44000000 ... 54824000 54912000 55000000]
[ 221692 226525 44064000 ... 54903744 54991872 55080000]]
m1, m2 and m3 are very large square matrices (those examples are 128x128, but they can go up to 2048x2048). Also m1*m2=m3.
My goal is to obtain m2 by using only m1 and m3. Someone told me this was possible, as m1*m2=m3 implies that (m1**-1) * m3 = m2 (I believe it was that, please correct me if I'm wrong) ; so I calculated the inverse of m1 :
m1**-1 = [[ 7.70884284e-01 -8.13188394e-01 -1.65131146e+13 ... -2.49697170e+12
-7.70160676e+12 -4.13395320e+13]
[-3.38144598e-01 2.54532610e-01 1.01286404e+13 ... -3.64296085e+11
2.60327813e+12 2.41783491e+13]
[ 1.77721050e-01 -3.54566231e-01 -5.00564604e+12 ... 5.82415184e+10
-5.98354744e+11 -1.29817153e+13]
...
[-6.56772812e-02 1.54498025e-01 3.21826474e+12 ... 2.61432526e+11
1.14203762e+12 3.61036457e+12]
[ 5.82732587e-03 -3.44252762e-02 -4.79430664e+11 ... 5.10855381e+11
-1.07679881e+11 -1.71485373e+12]
[ 6.55360708e-02 -8.24446025e-02 -1.19618881e+12 ... 4.45713678e+11
-3.48073716e+11 -4.89344092e+12]]
The result looked rather messy so I ran a test and multiplied m1**-1 and m1 to see if it worked :
(m1**-1)*m1 = [[-125.296875 , -117.34375 , -117.390625 , ..., -139.15625 ,
-155.203125 , -147.25 ],
[ 483.1640625 , 483.953125 , 482.7421875 , ..., 603.796875 ,
590.5859375 , 593.375 ],
[-523.22851562, -522.36328125, -523.49804688, ..., -633.07421875,
-635.20898438, -637.34375 ],
...,
[ 10.58691406, 11.68945312, 10.29199219, ..., 14.40429688,
13.00683594, 11.609375 ],
[ -5.32177734, -5.47949219, -4.63720703, ..., -5.28613281,
-5.31884766, -5.6015625 ],
[ -4.93554688, -3.58984375, -3.24414062, ..., -8.72265625,
-5.37695312, -8.03125 ]]
The result is different from the one expected (identity matrix). My guess is that m1 is too big, causing numerical imprecision. But if that previous calculation to get an identity matrix doesn't work properly, then (m1**-1)*m3 surely won't (and it doesn't). But I really can't decrease the matrix sizes for m1, m2 and m3 and in fact I'd like it to work with even bigger sizes (as said before, max size would be 2048x2048).
Would there be any way to be more precise with such calculations ? Is there an alternative that could work for bigger matrices ?
You are right, inverting a large matrix can be inefficient and numerically unstable. Luckily, there are methods in linear algebra that solve this problem without requiring an inverse.
In this case, m2 = np.linalg.solve(m1, m3) works.