Hello,
I’d like to provide an exact hessian to the IPOPT solver. If I understand correctly, the hessian of the reduced functional should be sufficient. Method hessian(self, m_array, m_dot_array) in the source returns hessian matrix action in the given direction. So I thought that to obtain hessian itself I just need to specify the identity matrix as m_dot_array. But apparently, it’s wrong, as I’m getting hessian related error during optimization.
As MWE I took poisson mother demo for IPOPT simplified by removing subdomain. I modified a couple of classes, defining methods that provide exact hessian according to this cyipopt documentation (page 16). I didn’t provide hessian structure, since I found it to be dense.
import numpy as np
from dolfin import *
from dolfin_adjoint import *
from pyadjoint.reduced_functional_numpy import ReducedFunctionalNumPy
from pyadjoint.optimization.ipopt_solver import _IPOptProblem
from pyadjoint.optimization.optimization_solver import OptimizationSolver
from pyadjoint.optimization.optimization_problem import MaximizationProblem
from functools import partial
set_log_level(LogLevel.ERROR)
n = 64
mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, "CG", 1)
f = interpolate(Expression("x[0]+x[1]", degree=1), V)
u = Function(V, name='State')
v = TestFunction(V)
F = (inner(grad(u), grad(v)) - f * v) * dx
bc = DirichletBC(V, 0.0, "on_boundary")
solve(F == 0, u, bc)
x = SpatialCoordinate(mesh)
w = Expression("sin(pi*x[0])*sin(pi*x[1])", degree=3)
d = 1 / (2 * pi ** 2)
d = Expression("d*w", d=d, w=w, degree=3)
alpha = Constant(1e-6)
J = assemble((0.5 * inner(u - d, u - d)) * dx + alpha / 2 * f ** 2 * dx)
control = Control(f)
rf = ReducedFunctional(J, control)
rf_np = ReducedFunctionalNumPy(rf)
class _MyIPOptProblem(_IPOptProblem):
"""API used by cyipopt for wrapping the problem"""
def __init__(self, objective, gradient, constraints, jacobian, hessian):
self.objective = objective
self.gradient = gradient
self.constraints = constraints
self.jacobian = jacobian
self.hessian = hessian
class MyIPOPTSolver(IPOPTSolver):
def __init__(self, problem, parameters=None):
OptimizationSolver.__init__(self, problem, parameters)
self.__build_ipopt_problem()
self._IPOPTSolver__set_parameters()
def __build_ipopt_problem(self):
"""Build the ipopt problem from the OptimizationProblem instance."""
from pyadjoint.ipopt import cyipopt
self.rfn = ReducedFunctionalNumPy(self.problem.reduced_functional)
(lb, ub) = self._IPOPTSolver__get_bounds()
(nconstraints, fun_g, jac_g, clb, cub) = self._IPOPTSolver__get_constraints()
# A callback that evaluates the functional and derivative.
J = self.rfn.__call__
dJ = partial(self.rfn.derivative, forget=False)
IdentityMatrix = np.identity(n + 1).flatten()
Hes = lambda x, lambd, objective_factor: objective_factor*self.rfn.hessian(x,IdentityMatrix)
nlp = cyipopt.problem(
n=len(ub), # length of control vector
lb=lb, # lower bounds on control vector
ub=ub, # upper bounds on control vector
m=nconstraints, # number of constraints
cl=clb, # lower bounds on constraints
cu=cub, # upper bounds on constraints
problem_obj=_MyIPOptProblem(
objective=J, # to evaluate the functional
gradient=dJ, # to evaluate the gradient
constraints=fun_g, # to evaluate the constraints
jacobian=jac_g, # to evaluate the constraint Jacobian
hessian = Hes, # to evaluate the reduced functional Hessian
),
)
nlp.addOption("print_level", 6)
if isinstance(self.problem, MaximizationProblem):
nlp.addOption('obj_scaling_factor', -1.0)
self.ipopt_problem = nlp
problem = MinimizationProblem(rf)
parameters = {"tolerance": 1e-12, } # 'print_level' : 5
solver = MyIPOPTSolver(problem, parameters=parameters)
f_opt = solver.solve()
This is the error I’m getting.
EXIT: Invalid number in NLP function or derivative detected.
IndexError: Out of bounds on buffer access (axis 0)
Exception ignored in: 'ipopt_wrapper.hessian_cb'
Traceback (most recent call last):
File "/home/ruslan/.local/lib/python3.8/site-packages/pyadjoint/optimization/ipopt_solver.py", line 199, in solve
results = self.ipopt_problem.solve(guess)
IndexError: Out of bounds on buffer access (axis 0)
Could you explain and show how I should use hessian() method of ReducedFunctionalNumPy class, if this is my mistake?
Thank you in advance!