Best way of finding KKT points for a Sympy polynomial

Question

I'm experimenting with running Gradient Descent (GD) on polynomials of some low degree - say 3 - with n variables, on a domain constrained by a hypercube [-1,1]^n. I want to compare the termination point of GD with the actual KKT point calculated with Sympy. For example, I have a Sympy expression with 3 variables a, b and c (that is n=3) like this:

0.2*a**3 - 0.1*a**2*b - 0.9*a**2*c + 0.5*a**2 - 0.1*a*b**2 - 0.4*a*b*c + 0.5*a*b - 0.6*a*c**2 - 0.5*a*c - 0.6*a - 0.1*b**3 + 0.4*b**2*c - 0.4*b**2 + 0.7*b*c**2 - 0.4*b*c + 0.9*b + 0.09*c**3 - 0.6*c**2 + 0.4*c + 0.22

What is the best, quickest way of finding all KKT points for this expressions with the constraints:

-1 <= a <= 1
-1 <= b <= 1
-1 <= c <= 1

(ideally the search would be limited only to the vicinity of the GD termination point to save computation time but I'm not sure how to do that rigorously)

Answer 1

Pre-compile the expression and its Jacobian to function objects, then pass them into minimize. Note that it attempts to be somewhat smarter than your intended gradient descent which may or may not be helpful; this defaults to Kraft's sequential least-squares programming method.

import string

import numpy as np
import scipy.optimize
import sympy

expr_str = (
    '0.2*a**3 - 0.1*a**2*b - 0.9*a**2*c + 0.5*a**2 - 0.1*a*b**2 - 0.4*a*b*c '
    '+ 0.5*a*b - 0.6*a*c**2 - 0.5*a*c - 0.6*a - 0.1*b**3 + 0.4*b**2*c '
    '- 0.4*b**2 + 0.7*b*c**2 - 0.4*b*c + 0.9*b + 0.09*c**3 - 0.6*c**2 + 0.4*c + 0.22'
)
symbols = {
    c: sympy.Symbol(name=c, real=True, finite=True)
    for c in string.ascii_lowercase
    if c in expr_str
}
expr = sympy.parse_expr(s=expr_str, local_dict=symbols)
expr_callable = sympy.lambdify(
    args=[symbols.values()], expr=expr,
)

jac = [
    sympy.diff(expr, sym)
    for sym in symbols.values()
]
jac_callable = sympy.lambdify(
    args=[symbols.values()], expr=jac,
)

result = scipy.optimize.minimize(
    fun=expr_callable,
    jac=jac_callable,
    x0=np.zeros(len(symbols)),
    bounds=scipy.optimize.Bounds(lb=-1, ub=1),
)
assert result.success
print(result)

CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
  message: CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL
  success: True
   status: 0
      fun: -4.020346688237437
        x: [ 5.139e-01 -1.000e+00 -1.000e+00]
      nit: 5
      jac: [-1.222e-08  3.839e+00  4.398e+00]
     nfev: 7
     njev: 7
 hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64>

If you really want a symbolic solution to a where b and c are equal to -1, you can then work backward, substituting b and c into the first Jacobian component and solving for a.

jac_rhs = jac[0].subs({
    symbols['b']: result.x[1],
    symbols['c']: result.x[2],
})
print(jac_rhs)
soln = sympy.solve(jac_rhs, symbols['a'],
                   dict=True, numerical=False, simplify=False, rational=True)
pprint(soln)

0.6*a**2 + 3.0*a - 1.7
[{a: -5/2 + sqrt(327)/6}, {a: -sqrt(327)/6 - 5/2}]