Exact Hessian for IPOPT


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


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)


    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
                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!