Experiment: Hessian of the lagrangian (#31)

* attempting outline for hessian

* first draft at a hessian

---------

Co-authored-by: William Zijie Zhang <[email protected]>
Co-authored-by: Daniel Cederberg <[email protected]>

Author: Transurgeon

cvxpy/atoms/atom.py

@@ -439,6 +439,16 @@ def grad(self):
 
         return result
 
+    def hess(self, duals):
+        """
+        Compute the hessian-vector product of a scalar function with dual
+        values.
+        Here the dual values should correspond to the dual values of the
+        constraint function or None if we are taking the hessian of the objective.
+        We must first slice the duals array to get the relevant values.
+        """
+        pass
+
     @abc.abstractmethod
     def _grad(self, values):
         """Gives the (sub/super)gradient of the atom w.r.t. each argument.

cvxpy/reductions/solvers/nlp_solvers/ipopt_nlpif.py

@@ -213,7 +213,6 @@ def __init__(self, problem, inital_point):
 
         def objective(self, x):
             """Returns the scalar value of the objective given x."""
-            # Set the variable value
             offset = 0
             for var in self.main_var:
                 size = var.size
@@ -227,7 +226,6 @@ def gradient(self, x):
             #import pdb 
             #pdb.set_trace()
             """Returns the gradient of the objective with respect to x."""
-            # compute the gradient using _grad
             offset = 0
             for var in self.main_var:
                 size = var.size
@@ -248,15 +246,13 @@ def gradient(self, x):
 
         def constraints(self, x):
             """Returns the constraint values."""
-            # Set the variable value
             offset = 0
             #import pdb 
             #pdb.set_trace()
             for var in self.main_var:
                 size = var.size
                 var.value = x[offset:offset+size].reshape(var.shape, order='F')
                 offset += size
-            
             # Evaluate all constraints
             constraint_values = []
             for constraint in self.problem.constraints:
@@ -265,15 +261,12 @@ def constraints(self, x):
 
         def jacobian(self, x):
             """Returns only the non-zero values of the Jacobian."""
-            #import pdb 
-            #pdb.set_trace()
             # Set variable values
             offset = 0
             for var in self.main_var:
                 size = var.size
                 var.value = x[offset:offset+size].reshape(var.shape, order='F')
                 offset += size
-            
             values = []
             for constraint in self.problem.constraints:
                 # get the jacobian of the constraint
@@ -329,9 +322,44 @@ def jacobianstructure(self):
                     col_offset += var.size
                 row_offset += constraint.size
                 self.jacobian_idxs[constraint] = constraint_jac
-
             return (np.array(rows), np.array(cols))
 
+        def hessian(self, x, duals, obj_factor):
+            offset = 0
+            for var in self.main_var:
+                size = var.size
+                var.value = x[offset:offset+size].reshape(var.shape, order='F')
+                offset += size
+            hess = np.zeros((x.size, x.size), dtype=np.float64)
+            # hess_dict = self.problem.objective.expr.hess(obj_factor)
+            # if we specify the problem in graph form (i.e. t=obj),
+            # the objective hessian will always be zero.
+            # To compute the hessian of the constraints, we need to gather the
+            # second derivatives from each constraint.
+            # This is done by looping through each constraint and each
+            # pair of variables and summing up the hessian contributions.
+            constr_offset = 0
+            for constraint in self.problem.constraints:
+                constr_hess = constraint.expr.hess()
+                # we have to make sure that the dual variables correspond
+                # to the constraints in the same order
+                constr_dual = duals[constr_offset:constr_offset + constraint.size]
+                row_offset = 0
+                for var1 in self.main_var:
+                    col_offset = 0
+                    for var2 in self.main_var:
+                        hess[var1.index * row_offset, var2.index * col_offset] += (
+                            constr_dual * constr_hess.get((var1, var2), 0).toarray()
+                        )
+                        col_offset += var2.size
+                    row_offset += var1.size
+                constr_offset += constraint.size
+            return hess
+
+        def hessianstructure(self):
+            pass
+
+
     class Bounds():
         def __init__(self, problem):
             self.problem = problem