Search code examples
pythonsympysymbolic-mathequation-solvingderivative

Set derivatives of a function to 0 in Python


I am trying to simultaneously solve a system of equations. The equations themselves are found by calculating the gradient of a function with some variables. I am using sympy and here is the code:

from sympy import *
m = Matrix(symbols('a b c', positive = True))

y = 4*log(m[0]) + 4*log(m[1]) + 4*log(m[2]) - 2*log(m[1] + m[2]) \
    - 2*log(m[0] + m[2]) - 2*log(m[0] + m[1]) - 6*log(m[0] + m[1] + m[2])


s = [diff(y, i) for i in m]
solve(s,m)

However I am getting the following error: "raise NotImplementedError('could not solve %s' % eq2)"

Can somebody please help me in solving this. Or there some other way in which I can calculate a bunch of gradients and then solve the system of equations obtained? I am fine with getting a numerical approximation and if multiple solutions are present, even one solution will be sufficient.

EDIT I understand that he objective that I have in the code shown above will have symmetric gradients. So here I am looking for a solution like (1,1,1) or (2,2,2). But in actual implementation, my objective function will have gradients that are not symmetric. So in a way I need to know the ratio between them.


Solution

  • import sympy as sp
    
    # define a vector of variables
    vm = sp.var('m0:3', real = True)
    
    y = 4*sp.log(vm[0]) + 4*sp.log(vm[1]) + 4*sp.log(vm[2]) - 2*sp.log(vm[1] + vm[2]) \
        - 2*sp.log(vm[0] + vm[2]) - 2*sp.log(vm[0] + vm[1]) - 6*sp.log(vm[0] + vm[1] + vm[2])
    

    The gradient w.r.t. to vm can be obtained as

    grad = [sp.diff(y, i) for i in vm]
    

    However, the result is a set of complicated rational polynomials which cannot be handled by sp.solve(grad, vm). We can help solve by performing some preprocessing on the equations, namely, factorizing and considering only the numerator:

    grad_numerators = [sp.numer(sp.diff(y, i).together()).factor() for i in vm]
    

    Now a call to

    sp.solve(grad_numerators,vm)
    

    gives a set of possible solutions.

    [{m1: -m2, m0: 0},
     {m0: -m1, m2: 0},
     {m1: m2, m0: -4*m2/3},
     {m1: 0, m0: -m2},
     {m1: -3*m2/4, m0: -3*m2/4},
     {m1: -4*m2/3, m0: m2}]
    

    Note that some of them may be invalid in the sense that they may correspond to a zero denominator for (some of) the grad elements, which were ignored in this derivation.