# Exact Hessian for IPOPT

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

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.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."""

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
constraints=fun_g,  # to evaluate the constraints
jacobian=jac_g,  # to evaluate the constraint Jacobian
hessian = Hes,  # to evaluate the reduced functional Hessian
),
)

if isinstance(self.problem, MaximizationProblem):

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?