Determinant in sympy

Question

Hello I have this code:

import numpy as np
from sympy import *
import sympy
 
x = Symbol('x')
x0 = Symbol('x0')
B = Symbol('B')
E = Symbol('E')
linewidth = Symbol('linewidth')

f = 1/pi*(linewidth/2)/((x - (E+2*B))**2 + (linewidth/2)**2)

def makeOneMeasurement(pCov,f,sigma,pvec,pmean):
    Iprior = np.linalg.inv(pCov)
    derivative_f = [f.diff(E),f.diff(B)]
    F = []
    for i in range(len(derivative_f)):
        for j in range(len(derivative_f)):
            g = derivative_f[i]*derivative_f[j]
            for k in range(len(pvec)):
                g = g.subs(pvec[k], pmean[k])
            F = F + [g]
    F = np.reshape(F,(2,2))/sigma**2
    F = sympy.Matrix(F)
    Iprior = sympy.Matrix(Iprior)
    
    return (Iprior + F).det()

If I call this function like this:

makeOneMeasurement([[1,0],[0,1]],f,0.1,["linewidth","E","B"],[0.1,990,0.11]).subs(x,990).evalf()

I get the output 1 (exactly 1). However, if instead of returning return (Iprior + F).det() I return just return (Iprior + F), and I call:

makeOneMeasurement([[1,0],[0,1]],f,0.1,["linewidth","E","B"],[0.1,990,0.11]).subs(x,990).evalf().det()

(basically exactly the same thing, but the det is called outside instead of inside the function), I get:

3653.9564196053

What am I doing wrong? I want to be able to return the determinant as a function of x (so I would like the use the first version). Why do I get 1 all the time instead of a function of x? Thank you!

Answer 1

Change your code to avoid all use of floats which includes avoiding all use of numpy routines like inv and reshape (use the sympy equivalents):

In [39]: import numpy as np
    ...: from sympy import *
    ...: import sympy
    ...: 
    ...: x = Symbol('x')
    ...: x0 = Symbol('x0')
    ...: B = Symbol('B')
    ...: E = Symbol('E')
    ...: linewidth = Symbol('linewidth')
    ...: 
    ...: f = 1/pi*(linewidth/2)/((x - (E+2*B))**2 + (linewidth/2)**2)
    ...: 
    ...: def makeOneMeasurement(pCov,f,sigma,pvec,pmean):
    ...:     Iprior = Matrix(pCov).inv()
    ...:     derivative_f = [f.diff(E),f.diff(B)]
    ...:     F = []
    ...:     for i in range(len(derivative_f)):
    ...:         for j in range(len(derivative_f)):
    ...:             g = derivative_f[i]*derivative_f[j]
    ...:             for k in range(len(pvec)):
    ...:                 g = g.subs(pvec[k], pmean[k])
    ...:             F = F + [g]
    ...:     F = Matrix(F).reshape(2,2)/sigma**2
    ...:     F = sympy.Matrix(F)
    ...:     Iprior = sympy.Matrix(Iprior)
    ...: 
    ...:     return (Iprior + F)
    ...: 

In [40]: M = makeOneMeasurement([[1,0],[0,1]],f,Rational('0.1'),["linewidth","E","B"],[Rational('0.1'
    ...: ),S(990),Rational('0.11')])

Then the exact matrix and determinant is:

In [41]: M
Out[41]: 
⎡                   2                                          ⎤
⎢      ⎛      49511⎞              ⎛      49511⎞ ⎛      99022⎞  ⎥
⎢      ⎜2⋅x - ─────⎟              ⎜2⋅x - ─────⎟⋅⎜4⋅x - ─────⎟  ⎥
⎢      ⎝        25 ⎠              ⎝        25 ⎠ ⎝        25 ⎠  ⎥
⎢────────────────────────── + 1   ───────────────────────────  ⎥
⎢                         4                                 4  ⎥
⎢     ⎛           2      ⎞              ⎛           2      ⎞   ⎥
⎢   2 ⎜⎛    49511⎞     1 ⎟            2 ⎜⎛    49511⎞     1 ⎟   ⎥
⎢4⋅π ⋅⎜⎜x - ─────⎟  + ───⎟         4⋅π ⋅⎜⎜x - ─────⎟  + ───⎟   ⎥
⎢     ⎝⎝      50 ⎠    400⎠              ⎝⎝      50 ⎠    400⎠   ⎥
⎢                                                              ⎥
⎢                                                   2          ⎥
⎢ ⎛      49511⎞ ⎛      99022⎞          ⎛      99022⎞           ⎥
⎢ ⎜2⋅x - ─────⎟⋅⎜4⋅x - ─────⎟          ⎜4⋅x - ─────⎟           ⎥
⎢ ⎝        25 ⎠ ⎝        25 ⎠          ⎝        25 ⎠           ⎥
⎢ ───────────────────────────    ────────────────────────── + 1⎥
⎢                           4                             4    ⎥
⎢       ⎛           2      ⎞          ⎛           2      ⎞     ⎥
⎢     2 ⎜⎛    49511⎞     1 ⎟        2 ⎜⎛    49511⎞     1 ⎟     ⎥
⎢  4⋅π ⋅⎜⎜x - ─────⎟  + ───⎟     4⋅π ⋅⎜⎜x - ─────⎟  + ───⎟     ⎥
⎣       ⎝⎝      50 ⎠    400⎠          ⎝⎝      50 ⎠    400⎠     ⎦

In [42]: M.det()
Out[42]: 
                   2  8                         2  7                             2  6              
10000000000000000⋅π ⋅x  - 79217600000000000000⋅π ⋅x  + 274549981652000000000000⋅π ⋅x  - 54372976605
───────────────────────────────────────────────────────────────────────────────────────────────────
                                                        2  8                         2  7          
                                     10000000000000000⋅π ⋅x  - 79217600000000000000⋅π ⋅x  + 2745499

                  2  5                                   2  4                                      
8974880000000000⋅π ⋅x  + 673015111919049368767000000000⋅π ⋅x  - 533146420076341661747549392000000⋅π
───────────────────────────────────────────────────────────────────────────────────────────────────
                   2  6                                2  5                                   2  4 
81652000000000000⋅π ⋅x  - 543729766058974880000000000⋅π ⋅x  + 673015111919049368767000000000⋅π ⋅x  

2  3                      2                                         2  2                           
 ⋅x  + 50000000000000000⋅x  + 263966124524722600510236863978120000⋅π ⋅x  - 746812961137426061526108
───────────────────────────────────────────────────────────────────────────────────────────────────
                                     2  3                                         2  2             
- 533146420076341661747549392000000⋅π ⋅x  + 263966124524722600510236863978120000⋅π ⋅x  - 7468129611

                2                                                                                  
47700028830400⋅π ⋅x - 99022000000000000000⋅x + 49026782420000000000000 + 92438641532871794599827086
───────────────────────────────────────────────────────────────────────────────────────────────────
                              2                                               2                    
3742606152610847700028830400⋅π ⋅x + 9243864153287179459982708676348253060561⋅π                     

                2
76348253060561⋅π 
─────────────────
                 

As you can see there are huge numbers in there. You can use an exact substitution to get the formula for the determinant and then evaluate it:

In [43]: M.det().subs(x, 990)
Out[43]: 
             2                   
67122964561⋅π  + 2420000000000000
─────────────────────────────────
                       2         
          67122964561⋅π          

In [44]: M.det().subs(x, 990).n()
Out[44]: 3653.95642136943

It is not reasonable to expect to be able to substitute low precision floats into an expression like this without a significant loss of accuracy so use exact calculations until the final call to evalf.