Converting a piece of octave (mrdivide) code to numpy

Question

I am trying to convert a piece of octave code to numpy but I got different results from octave and numpy. Here are my Data (in fact they are far larger than given data below):

A =

 Columns 1 through 6:

   1.000000000000000e+00   9.954090291002812e-01   9.820064469806473e-01   9.908807141567176e-01   9.954090291002812e-01   9.908807141567179e-01
   9.954090291002812e-01   1.000000000000000e+00   9.954090291002812e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01
   9.820064469806473e-01   9.954090291002812e-01   1.000000000000000e+00   9.567491788333946e-01   9.776578008378289e-01   9.908807141567179e-01
   9.908807141567176e-01   9.776578008378289e-01   9.567491788333946e-01   1.000000000000000e+00   9.954090291002812e-01   9.820064469806473e-01
   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01   9.954090291002812e-01   1.000000000000000e+00   9.954090291002812e-01
   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01   9.820064469806473e-01   9.954090291002812e-01   1.000000000000000e+00
   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01   9.608235911808946e-01   9.820064469806473e-01   9.954090291002812e-01
   9.776578008378289e-01   9.649505047327671e-01   9.448300707857176e-01   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01
   9.820064469806473e-01   9.776578008378289e-01   9.649505047327671e-01   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01
   9.776578008378289e-01   9.820064469806473e-01   9.776578008378289e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01
   9.649505047327671e-01   9.776578008378289e-01   9.820064469806473e-01   9.567491788333946e-01   9.776578008378289e-01   9.908807141567179e-01
   9.567491788333946e-01   9.448300707857176e-01   9.259179407344914e-01   9.820064469806473e-01   9.776578008378289e-01   9.649505047327671e-01
   9.608235911808946e-01   9.567491788333946e-01   9.448300707857176e-01   9.776578008378289e-01   9.820064469806473e-01   9.776578008378289e-01
   9.567491788333946e-01   9.608235911808946e-01   9.567491788333946e-01   9.649505047327671e-01   9.776578008378289e-01   9.820064469806473e-01
   9.448300707857176e-01   9.567491788333946e-01   9.608235911808946e-01   9.448300707857176e-01   9.649505047327673e-01   9.776578008378289e-01

 Columns 7 through 12:

   9.776578008378289e-01   9.776578008378289e-01   9.820064469806473e-01   9.776578008378289e-01   9.649505047327671e-01   9.567491788333946e-01
   9.908807141567179e-01   9.649505047327671e-01   9.776578008378289e-01   9.820064469806473e-01   9.776578008378289e-01   9.448300707857176e-01
   9.954090291002812e-01   9.448300707857176e-01   9.649505047327671e-01   9.776578008378289e-01   9.820064469806473e-01   9.259179407344914e-01
   9.608235911808946e-01   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01   9.567491788333946e-01   9.820064469806473e-01
   9.820064469806473e-01   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01   9.776578008378289e-01
   9.954090291002812e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01   9.649505047327671e-01
   1.000000000000000e+00   9.567491788333946e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01   9.448300707857176e-01
   9.567491788333946e-01   1.000000000000000e+00   9.954090291002812e-01   9.820064469806473e-01   9.608235911808946e-01   9.954090291002812e-01
   9.776578008378289e-01   9.954090291002812e-01   1.000000000000000e+00   9.954090291002812e-01   9.820064469806473e-01   9.908807141567179e-01
   9.908807141567179e-01   9.820064469806473e-01   9.954090291002812e-01   1.000000000000000e+00   9.954090291002812e-01   9.776578008378289e-01
   9.954090291002812e-01   9.608235911808946e-01   9.820064469806473e-01   9.954090291002812e-01   1.000000000000000e+00   9.567491788333946e-01
   9.448300707857176e-01   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01   9.567491788333946e-01   1.000000000000000e+00
   9.649505047327673e-01   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01   9.954090291002812e-01
   9.776578008378289e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01   9.820064469806473e-01
   9.820064469806473e-01   9.567491788333946e-01   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01   9.608235911808946e-01

 Columns 13 through 15:

   9.608235911808946e-01   9.567491788333946e-01   9.448300707857176e-01
   9.567491788333946e-01   9.608235911808946e-01   9.567491788333946e-01
   9.448300707857176e-01   9.567491788333946e-01   9.608235911808946e-01
   9.776578008378289e-01   9.649505047327671e-01   9.448300707857176e-01
   9.820064469806473e-01   9.776578008378289e-01   9.649505047327673e-01
   9.776578008378289e-01   9.820064469806473e-01   9.776578008378289e-01
   9.649505047327673e-01   9.776578008378289e-01   9.820064469806473e-01
   9.908807141567179e-01   9.776578008378289e-01   9.567491788333946e-01
   9.954090291002812e-01   9.908807141567179e-01   9.776578008378289e-01
   9.908807141567179e-01   9.954090291002812e-01   9.908807141567179e-01
   9.776578008378289e-01   9.908807141567179e-01   9.954090291002812e-01
   9.954090291002812e-01   9.820064469806473e-01   9.608235911808946e-01
   1.000000000000000e+00   9.954090291002812e-01   9.820064469806473e-01
   9.954090291002812e-01   1.000000000000000e+00   9.954090291002812e-01
   9.820064469806473e-01   9.954090291002812e-01   1.000000000000000e+00

and

b =

  -1.024208397018539
  -1.055718555015945
  -1.066560607689640
  -1.187368188387253
  -1.258866007703282
  -1.305258462589997
  -1.321354530870290
  -1.333661132027602
  -1.421384660329320
  -1.478743779481671
  -1.498725636719488
  -1.385960967135295
  -1.479663779776475
  -1.541078471216082
  -1.562500000000000

In octave I have x = b'/A. I could not find corresponding numpy function for /. So far I tried x = np.dot(b.T,np.linalg.inv(A)) from numpy but the results are not same as octave.

The results of octave for x = b'/A are:

x =

 Columns 1 through 6:

  -5.642309525591432e+00   7.813412870218545e+00  -1.559855155426489e+02  -3.597241224212262e+01   2.201914551287831e+02  -3.100354445411479e+02

 Columns 7 through 12:

   7.253956488595386e+02   7.369595892720794e+01  -4.313273469816049e+02   6.064725968037579e+02  -9.855235323530542e+02  -4.111380598448122e+01

 Columns 13 through 15:

   2.334109900297194e+02  -3.269547704109582e+02   4.254317069619117e+02

and the results of numpy are

x=np.array([[-5.642310487492068,    7.813414778371225,  -155.9855165364133, -35.9724138597885,  220.1914623928024,  -310.0354544342263, 725.3956532399461,  73.69596218669903,  -431.3273588509765, 606.472611254314,   -985.5235383154941, -41.11380770278629, 233.4109958125337,  -326.9547770833597, 425.4317096135928]])

I would appreciate if someone can help me to find same results as in octave by numpy. Or closer than current the results to Octave results.

I just ran both arrays through `numpy.allclose()`. They are the same. — Nils Werner, Aug 30 '19 at 13:20
Explicitly computing the inverse leads to poor results. Use [scypy.linalg.solve](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html) instead (NumPy also has one). You will have to re-write your problem from `xA=b` to `Ax=b`. — Cris Luengo, Aug 30 '19 at 13:20
@KingStone, actually the code is long. But you can save data in octave by `save -6 matrix.mat A b` and then import to numpy. — For Bonder, Aug 30 '19 at 14:28

score 0 · Answer 1 · answered Aug 30 '19 at 13:15

It looks like Matlab and numpy use different precisions for their calculations. I'd try to force both to use the some floating point precision.

By default matlab uses a 16-bit precision for numerical calculation (see https://www.mathworks.com/help/symbolic/increase-precision-of-numeric-calculations.html). On this site it is also describe how to change the precision.

numpy, as far as I am concerned uses 32-bit precision (You might find this article useful: how set numpy floating point accuracy?).

score 0 · Answer 2 · answered Aug 30 '19 at 13:25

0

Take both arrays (copy and pasted from your post)

x = [-5.642309525591432e+00, 7.813412870218545e+00, -1.559855155426489e+02, -3.597241224212262e+01, 2.201914551287831e+02, -3.100354445411479e+02, 7.253956488595386e+02, 7.369595892720794e+01, -4.313273469816049e+02, 6.064725968037579e+02, -9.855235323530542e+02, -4.111380598448122e+01, 2.334109900297194e+02, -3.269547704109582e+02, 4.254317069619117e+02]
y = [-5.642310487492068,    7.813414778371225,  -155.9855165364133, -35.9724138597885,  220.1914623928024,  -310.0354544342263, 725.3956532399461,  73.69596218669903,  -431.3273588509765, 606.472611254314,   -985.5235383154941, -41.11380770278629, 233.4109958125337,  -326.9547770833597, 425.4317096135928]

and run them through allclose

import numpy as np
np.allclose(x, y) # True

They are (within a some tolerance) the same.

answered Aug 30 '19 at 13:25

Nils Werner

34,832
7
76
98

I need more accuracy – For Bonder Aug 30 '19 at 13:27
More than 1e-5? Why? And you know that float numbers are inherently of limited precision? – Nils Werner Aug 30 '19 at 13:29
I know, but I need more accurate(closer to octave) results if possible. Thank you. – For Bonder Aug 30 '19 at 13:38
Also keep in mind that the values displayed to you might be only an approximation of the value saved – Nils Werner Aug 30 '19 at 13:38

Converting a piece of octave (mrdivide) code to numpy

2 Answers2