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?

Answer 1

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.

Answer 2

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