Precision loss numpy - mpmath

Question

I use numpy and mpmath in my Python programm. I use numpy, because it allows an easy access to many linear algebra operations. But because numpy's solver for linear equations is not that exact, i use mpmath for more precision operations. After i compute the solution of a System:

solution = mpmath.lu_solve(A,b)

i want the solution as an array. So i use

array = np.zeros(m)

and then do a loop for setting the values:

for i in range(m):
    array[i] = solution[i]

or

for i in range(m):
    array.put([i],solution[i])

but with both ways i get again numerical instabilities like:

solution[0] = 12.375
array[0] = 12.37500000000000177636

Is there a way to avoid these errors?

that's not a 'numerical instability', its the limited precision of floating point numbers. either don't use floating point numbers, or reconsider the question if this limited precision is infact a problem of any significance for your application. — Eelco Hoogendoorn, Aug 12 '14 at 14:10

score 2 · Accepted Answer · answered Jan 16 '14 at 15:30

2

numpy ndarrays have homogeneous type. When you make array, the default dtype will be some type of float, which doesn't have as much precision as you want:

>>> array = np.zeros(3)
>>> array
array([ 0.,  0.,  0.])
>>> array.dtype
dtype('float64')

You can get around this by using dtype=object:

>>> mp.mp.prec = 65
>>> mp.mpf("12.37500000000000177636")
mpf('12.37500000000000177636')
>>> array = np.zeros(3, dtype=object)
>>> array[0] = 12.375
>>> array[1] = mp.mpf("12.37500000000000177636")
>>> array
array([12.375, mpf('12.37500000000000177636'), 0], dtype=object)

but note that there's a significant performance hit when you do this.

answered Jan 16 '14 at 15:30

DSM

342,061
65
592
494

Can i go on with this array and do linear algebra? For example, use your last array as a right side of a linear System? And what type is the solution then? – dalvo Jan 16 '14 at 15:38
You can perform user-level manipulations like base arithmetic. Whether any given builtin function will coerce down to float depends upon how it was implemented, so I'm not sure in general. `.dot` seems to work okay. – DSM Jan 16 '14 at 15:42
Ok, i will try it. But maybe there is another way. Thanks so far. – dalvo Jan 16 '14 at 15:44
2

I'm almost certain that most `linalg` functions, e.g. `solve`, will cast the data to `np.double` before solving. – Jaime Jan 16 '14 at 16:46
But i use the lu_solve() from the mpmath package, not from the numpy. And i think that the mpmath precision for solving the System is the same as definied by mp.prec. But the question is more about copy the mpmath-solution stable to numpy arrays. – dalvo Jan 16 '14 at 16:52

score 1 · Answer 2 · answered Aug 12 '14 at 11:44

For the completeness, and for people like me who stumbled upon this question because numpy's linear solver is not exact enough (it seems to be able to handle 64bit numbers, only), there is also sympy.

The API is somewhat similar to numpy, but needs a few tweaks every now and then.

In [104]: A = Matrix([
[17928014155669123106522437234449354737723367262236489360399559922258155650097260907649387867023242961198972825743674594974017771680414642705007756271459833, 13639120912900071306285490050678803027789658562874829601953000463099941578381997997439951418291413106684405816668933580642992992427754602071359317117391198,  2921704428390104906296926198429197524950528698980675801502622843572749765891484935643316840553487900050392953088680445022408396921815210925936936841894852],
[14748352608418286197003564525494635018601621545162877692512866070970871867506554630832144013042243382377181384934564249544095078709598306314885920519882886,  2008780320611667023380867301336185953729900109553256489538663036530355388609791926150229595099056264556936500639831205368501493630132784265435798020329993,  6522019637107271075642013448499575736343559556365957230686263307525076970365959710482607736811185215265865108154015798779344536283754814484220624537361159],
[ 5150176345214410132408539250659057272148578629292610140888144535962281110335200803482349370429701981480772441369390017612518504366585966665444365683628345,  1682449005116141683379601801172780644784704357790687066410713584101285844416803438769144460036425678359908733332182367776587521824356306545308780262109501, 16960598957857158004200152340723768697140517883876375860074812414430009210110270596775612236591317858945274366804448872120458103728483749408926203642159476]])

In [105]: B = Matrix([
   .....: [13229751631544149067279482127723938103350882358472000559554652108051830355519740001369711685002280481793927699976894894174915494730969365689796995942384549941729746359],
   .....: [ 6297029075285965452192058994038386805630769499325569222070251145048742874865001443012028449109256920653330699534131011838924934761256065730590598587324702855905568792],
   .....: [ 2716399059127712137195258185543449422406957647413815998750448343033195453621025416723402896107144066438026581899370740803775830118300462801794954824323092548703851334]])

In [106]: A.solve(B)
Out[106]: 
Matrix([
[358183301733],
[498758543457],
[  1919512167]])

In [107]:

Precision loss numpy - mpmath

2 Answers2

Linked