Add support for objective and constraint hessians and jac=True option for constraints in the scipy interface
#121
Conversation
…vectar-valued constraint functions. The latter is needed when evaluating the hessian.
jac=True option for constraints in the scipy interface
|
@jhelgert Is this ready to be reviewed? Whenever you think it is done, let us know. |
|
@moorepants Yes, I think it's ready. However, there are some unintentional formatting changes because my code formatter seems to ignore the settings :/. |
Co-authored-by: Jason K. Moore <moorepants@gmail.com>
…t/cyipopt into scipy_interface_hessians
|
@moorepants I think it's ready to be reviewed again. Thanks! |
|
@jhelgert Can you merge in master so the tests at least run on linux? |
Co-authored-by: Sam Brockie <sgb39@cam.ac.uk>
|
I guess appveyor is still running, just not webhooked into here anymore. Appveyor also failed: https://ci.appveyor.com/project/moorepants/cyipopt/builds/40165736 |
|
Thanks to both of you for the reviews so far! PS: Appveyor fails with This was already fixed with the last commit. |
|
@moorepants All tests pass now on Linux. However, there's still one Appveyor job that fails: https://ci.appveyor.com/project/moorepants/cyipopt/builds/40167487 |
Looking at the Appveyor readout it is failing on |
|
@brocksam Thanks for the discussion! Can you please have another look? I think it should be much clearer now thanks to Edit: Here's a small performance comparison (using import numpy as np
from scipy.optimize import rosen, rosen_der
from cyipopt import minimize_ipopt
# Tiny example
M = 2000
N = 3000
A = np.random.randn(M, N)
b = np.random.randn(M)
# Precalculate the matrix-vector products
AtA = A.T @ A
btA = b.T @ A
btb = b.T @ b
Atb = btA.T
# return the constraint value and the jacobian
def con_and_jac(x, rhs):
AtAx = AtA @ x
con_value = x.T @ AtAx - 2*btA @ x + btb
jac = 2*AtAx - 2*Atb
return con_value - rhs, jac
## Variant 1 (use jac=True in order to cache AtA @ x)
# x.T @ AtA @ x - 2*btA @ x + btb - 60000 >= 0
constr1 = {"type": "ineq", "fun": lambda x: con_and_jac(x, 60000), "jac": True}
## Variant 2 (no caching of intermediate results)
# x.T @ AtA @ x - 2*btA @ x + btb - 60000 >= 0
constr2 = {"type": "ineq", "fun": lambda x: x.T @ AtA @ x - 2*btA @ x + btb - 60000, "jac": lambda x: 2*AtA @ x - 2*Atb}Timing res = minimize_ipopt(rosen, x0=2.0*np.ones(N), constraints=constr1, options={'disp': 5, 'maxiter': 100})
res = minimize_ipopt(rosen, x0=2.0*np.ones(N), constraints=constr2, options={'disp': 5, 'maxiter': 100})yields 11.2s and 15.1s on my machine. |
Appveryor still failing; Environment: PYTHON=C:\Miniconda37-x64, PYTHON_VERSION=3.8, PYTHON_ARCH=64, CONDA_PY=38, IPOPT=, IPOPTWINDIR=C:\projects\cyipopt\Ipopt-3.13.3-win64-msvs2019-md
@jhelgert this is still a problem. |
|
I think we are getting closer. I removed the moody # macOS: 1.087746, 1.143973, 1.306229, 1.710728, 2.880089
# Win: 1.087746, 1.143977, 1.306237, 1.710748, 2.880159
# Linux: 1.087745, 1.143968, 1.306218, 1.710698, 2.879984Both |
To make sure I've understood this correctly, you're still operating under the assumption that there is a corner case where the different |
Yes, exactly.
Before writing it in C++, I have done it quickly with CasADi which computes all gradients, jacobians and Hessians by algorithmic differentiation and interfaces different solvers including Ipopt. Here, I can reproduce the problem: import casadi as cas
opti = cas.Opti()
x = opti.variable(5)
rosen = sum(100*(x[i+1]-x[i]**2)**2 + (1-x[i])**2 for i in range(4))
opti.minimize(rosen)
opti.subject_to(rosen - 1.0 >= 0)
opti.solver("ipopt", {}, {"mu_strategy": "adaptive"})
opti.set_initial(x, [1.3, 0.7, 0.8, 1.9, 1.2])
sol = opti.solve()
print(sol.value(x))
print(sol.value(rosen))which gives me Using the default 'mu_strategy' the objective value is again 1.0000000578559098 as expected. Therefore, it should indeed be a quirk of Ipopt. Should I still double-check it again with Ipopt's C++ API? |
|
Excellent idea to use CasADi! As you've confirmed it's an Ipopt quirk independent of CyIpopt, that's definitely sufficient. |
There was a problem hiding this comment.
Other than the docstring on the new test, I'm very happy with this. Thanks so much for all of your work! I'll defer to @moorepants to give final approval and make the merge.
|
@moorepants Any objections from your side? |
|
I'm on vacation right now and haven't had time to do a final check. If Sam has time for a final review and merge, he can do that. Otherwise, I'll try to review in the next week if I can find a bit of time. |
Implements #99 and #122.
hessfield in the constraint dictionary. The hessian must have the signaturehess(x,v) -> np.ndarrayand returnv[0] * Hg1_(x) + ... + v[m-1] * Hg_m(x). Here,vis the np.ndarray with shapemcontaining the lagrangian multipliers for the constraint andHg_idenotes the hessian matrix of the constraint componentg_i : R^N -> R.Current limitations:
hessp.