diff --git a/python/atoms/add.h b/python/atoms/add.h
index 4b86b79..2f62a7a 100644
--- a/python/atoms/add.h
+++ b/python/atoms/add.h
@@ -24,6 +24,7 @@ static PyObject *py_make_add(PyObject *self, PyObject *args)
         PyErr_SetString(PyExc_RuntimeError, "failed to create add node");
         return NULL;
     }
+    expr_retain(node); /* Capsule owns a reference */
     return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
 }
 
diff --git a/python/atoms/constant.h b/python/atoms/constant.h
index da2af09..d5fcba3 100644
--- a/python/atoms/constant.h
+++ b/python/atoms/constant.h
@@ -28,6 +28,7 @@ static PyObject *py_make_constant(PyObject *self, PyObject *args)
         PyErr_SetString(PyExc_RuntimeError, "failed to create constant node");
         return NULL;
     }
+    expr_retain(node); /* Capsule owns a reference */
     return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
 }
 
diff --git a/python/atoms/exp.h b/python/atoms/exp.h
index 9505daf..9ac7fc5 100644
--- a/python/atoms/exp.h
+++ b/python/atoms/exp.h
@@ -23,6 +23,7 @@ static PyObject *py_make_exp(PyObject *self, PyObject *args)
         PyErr_SetString(PyExc_RuntimeError, "failed to create exp node");
         return NULL;
     }
+    expr_retain(node); /* Capsule owns a reference */
     return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
 }
 
diff --git a/python/atoms/linear.h b/python/atoms/linear.h
index 0f85b91..c7afcae 100644
--- a/python/atoms/linear.h
+++ b/python/atoms/linear.h
@@ -54,6 +54,7 @@ static PyObject *py_make_linear(PyObject *self, PyObject *args)
         PyErr_SetString(PyExc_RuntimeError, "failed to create linear node");
         return NULL;
     }
+    expr_retain(node); /* Capsule owns a reference */
     return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
 }
 
diff --git a/python/atoms/log.h b/python/atoms/log.h
index 3d0314e..738b4e5 100644
--- a/python/atoms/log.h
+++ b/python/atoms/log.h
@@ -23,6 +23,7 @@ static PyObject *py_make_log(PyObject *self, PyObject *args)
         PyErr_SetString(PyExc_RuntimeError, "failed to create log node");
         return NULL;
     }
+    expr_retain(node); /* Capsule owns a reference */
     return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
 }
 
diff --git a/python/atoms/neg.h b/python/atoms/neg.h
index bda69d0..3927b0c 100644
--- a/python/atoms/neg.h
+++ b/python/atoms/neg.h
@@ -23,6 +23,7 @@ static PyObject *py_make_neg(PyObject *self, PyObject *args)
         PyErr_SetString(PyExc_RuntimeError, "failed to create neg node");
         return NULL;
     }
+    expr_retain(node); /* Capsule owns a reference */
     return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
 }
 
diff --git a/python/atoms/promote.h b/python/atoms/promote.h
index fb258b4..97e405c 100644
--- a/python/atoms/promote.h
+++ b/python/atoms/promote.h
@@ -24,6 +24,7 @@ static PyObject *py_make_promote(PyObject *self, PyObject *args)
         PyErr_SetString(PyExc_RuntimeError, "failed to create promote node");
         return NULL;
     }
+    expr_retain(node); /* Capsule owns a reference */
     return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
 }
 
diff --git a/python/atoms/sum.h b/python/atoms/sum.h
index 8a8480e..8d01604 100644
--- a/python/atoms/sum.h
+++ b/python/atoms/sum.h
@@ -24,6 +24,7 @@ static PyObject *py_make_sum(PyObject *self, PyObject *args)
         PyErr_SetString(PyExc_RuntimeError, "failed to create sum node");
         return NULL;
     }
+    expr_retain(node); /* Capsule owns a reference */
     return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
 }
 
diff --git a/python/atoms/variable.h b/python/atoms/variable.h
index 5438c44..2dc1972 100644
--- a/python/atoms/variable.h
+++ b/python/atoms/variable.h
@@ -17,6 +17,7 @@ static PyObject *py_make_variable(PyObject *self, PyObject *args)
         PyErr_SetString(PyExc_RuntimeError, "failed to create variable node");
         return NULL;
     }
+    expr_retain(node); /* Capsule owns a reference */
     return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
 }
 
diff --git a/python/bindings.c b/python/bindings.c
index 197390f..6e202d3 100644
--- a/python/bindings.c
+++ b/python/bindings.c
@@ -16,6 +16,7 @@
 /* Include problem bindings */
 #include "problem/constraint_forward.h"
 #include "problem/gradient.h"
+#include "problem/hessian.h"
 #include "problem/init_derivatives.h"
 #include "problem/jacobian.h"
 #include "problem/make_problem.h"
@@ -53,6 +54,8 @@ static PyMethodDef DNLPMethods[] = {
      "Compute objective gradient"},
     {"problem_jacobian", py_problem_jacobian, METH_VARARGS,
      "Compute constraint jacobian"},
+    {"problem_hessian", py_problem_hessian, METH_VARARGS,
+     "Compute Lagrangian Hessian"},
     {NULL, NULL, 0, NULL}};
 
 static struct PyModuleDef dnlp_module = {PyModuleDef_HEAD_INIT, "DNLP_diff_engine",
diff --git a/python/convert.py b/python/convert.py
index 92920cb..669f269 100644
--- a/python/convert.py
+++ b/python/convert.py
@@ -159,3 +159,22 @@ def jacobian(self) -> sparse.csr_matrix:
         """Compute jacobian of constraints. Call constraint_forward first. Returns scipy CSR matrix."""
         data, indices, indptr, shape = diffengine.problem_jacobian(self._capsule)
         return sparse.csr_matrix((data, indices, indptr), shape=shape)
+
+    def hessian(self, obj_factor: float, lagrange: np.ndarray) -> sparse.csr_matrix:
+        """Compute Lagrangian Hessian.
+
+        Computes: obj_factor * H_obj + sum(lagrange_i * H_constraint_i)
+
+        Call objective_forward and constraint_forward before this.
+
+        Args:
+            obj_factor: Weight for objective Hessian
+            lagrange: Array of Lagrange multipliers (length = total_constraint_size)
+
+        Returns:
+            scipy CSR matrix of shape (n_vars, n_vars)
+        """
+        data, indices, indptr, shape = diffengine.problem_hessian(
+            self._capsule, obj_factor, lagrange
+        )
+        return sparse.csr_matrix((data, indices, indptr), shape=shape)
diff --git a/python/problem/hessian.h b/python/problem/hessian.h
new file mode 100644
index 0000000..41c2788
--- /dev/null
+++ b/python/problem/hessian.h
@@ -0,0 +1,76 @@
+#ifndef PROBLEM_HESSIAN_H
+#define PROBLEM_HESSIAN_H
+
+#include "common.h"
+
+/*
+ * py_problem_hessian: Compute Lagrangian Hessian
+ *
+ * Args:
+ *   prob_capsule: PyCapsule containing problem pointer
+ *   obj_factor: Scaling factor for objective Hessian (double)
+ *   lagrange: Array of Lagrange multipliers (numpy array, length = total_constraint_size)
+ *
+ * Returns:
+ *   Tuple of (data, indices, indptr, (m, n)) for scipy.sparse.csr_matrix
+ */
+static PyObject *py_problem_hessian(PyObject *self, PyObject *args)
+{
+    PyObject *prob_capsule;
+    double obj_factor;
+    PyObject *lagrange_obj;
+
+    if (!PyArg_ParseTuple(args, "OdO", &prob_capsule, &obj_factor, &lagrange_obj))
+    {
+        return NULL;
+    }
+
+    problem *prob =
+        (problem *) PyCapsule_GetPointer(prob_capsule, PROBLEM_CAPSULE_NAME);
+    if (!prob)
+    {
+        PyErr_SetString(PyExc_ValueError, "invalid problem capsule");
+        return NULL;
+    }
+
+    /* Convert lagrange to contiguous C array */
+    PyArrayObject *lagrange_arr = (PyArrayObject *) PyArray_FROM_OTF(
+        lagrange_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
+    if (!lagrange_arr)
+    {
+        return NULL;
+    }
+
+    double *lagrange = (double *) PyArray_DATA(lagrange_arr);
+
+    /* Compute Hessian */
+    problem_hessian(prob, obj_factor, lagrange);
+
+    Py_DECREF(lagrange_arr);
+
+    /* Extract CSR components and return as tuple */
+    CSR_Matrix *H = prob->lagrange_hessian;
+    npy_intp nnz = H->nnz;
+    npy_intp n_plus_1 = H->n + 1;
+
+    PyObject *data = PyArray_SimpleNew(1, &nnz, NPY_DOUBLE);
+    PyObject *indices = PyArray_SimpleNew(1, &nnz, NPY_INT32);
+    PyObject *indptr = PyArray_SimpleNew(1, &n_plus_1, NPY_INT32);
+
+    if (!data || !indices || !indptr)
+    {
+        Py_XDECREF(data);
+        Py_XDECREF(indices);
+        Py_XDECREF(indptr);
+        return NULL;
+    }
+
+    /* Copy CSR data using memcpy for efficiency */
+    memcpy(PyArray_DATA((PyArrayObject *) data), H->x, nnz * sizeof(double));
+    memcpy(PyArray_DATA((PyArrayObject *) indices), H->i, nnz * sizeof(int));
+    memcpy(PyArray_DATA((PyArrayObject *) indptr), H->p, n_plus_1 * sizeof(int));
+
+    return Py_BuildValue("(OOO(ii))", data, indices, indptr, H->m, H->n);
+}
+
+#endif /* PROBLEM_HESSIAN_H */
diff --git a/python/tests/test_constrained.py b/python/tests/test_constrained.py
index 9574d3d..4409b13 100644
--- a/python/tests/test_constrained.py
+++ b/python/tests/test_constrained.py
@@ -186,6 +186,147 @@ def test_repeated_evaluations():
     assert np.allclose(jac2.toarray(), np.diag(-np.exp(u2)))
 
 
+def test_hessian_objective_only():
+    """Test Hessian of sum(log(x)) with no constraints.
+
+    For f(x) = sum(log(x)), Hessian is diagonal: d²f/dx_i² = -1/x_i²
+    """
+    x = cp.Variable(3)
+    obj = cp.sum(cp.log(x))
+    problem = cp.Problem(cp.Minimize(obj))
+    prob = C_problem(problem)
+
+    u = np.array([1.0, 2.0, 4.0])
+    prob.init_derivatives()
+
+    prob.objective_forward(u)
+    # No constraints, so pass empty lagrange array
+    H = prob.hessian(obj_factor=1.0, lagrange=np.array([]))
+
+    # Expected: diag([-1/1², -1/2², -1/4²]) = diag([-1, -0.25, -0.0625])
+    expected_H = np.diag([-1.0, -0.25, -0.0625])
+    assert np.allclose(H.toarray(), expected_H)
+
+
+def test_hessian_with_constraint():
+    """Test Lagrangian Hessian: sum(log(x)) + lagrange * (-exp(x)).
+
+    Objective Hessian: diag([-1/x²])
+    Constraint (-exp(x)) Hessian: diag([-exp(x)])
+
+    Lagrangian = obj + lagrange * constraint
+    """
+    x = cp.Variable(2)
+    obj = cp.sum(cp.log(x))
+    constraints = [cp.exp(x) >= 0]  # becomes -exp(x) <= 0
+
+    problem = cp.Problem(cp.Minimize(obj), constraints)
+    prob = C_problem(problem)
+
+    u = np.array([1.0, 2.0])
+    prob.init_derivatives()
+
+    prob.objective_forward(u)
+    prob.constraint_forward(u)
+
+    # Lagrange multipliers for the constraint (size 2)
+    lagrange = np.array([0.5, 0.5])
+    H = prob.hessian(obj_factor=1.0, lagrange=lagrange)
+
+    # Constraint expr is -exp(x), so its Hessian is -exp(x)
+    # H[0,0] = -1/1² + 0.5*(-exp(1))
+    # H[1,1] = -1/4 + 0.5*(-exp(2))
+    expected_00 = -1.0 + 0.5 * (-np.exp(1.0))
+    expected_11 = -0.25 + 0.5 * (-np.exp(2.0))
+    expected_H = np.diag([expected_00, expected_11])
+
+    assert np.allclose(H.toarray(), expected_H)
+
+
+def test_hessian_obj_factor_zero():
+    """Test Hessian with obj_factor=0 (constraint Hessian only)."""
+    x = cp.Variable(2)
+    obj = cp.sum(cp.log(x))
+    constraints = [cp.exp(x) >= 0]  # becomes -exp(x) <= 0
+
+    problem = cp.Problem(cp.Minimize(obj), constraints)
+    prob = C_problem(problem)
+
+    u = np.array([1.0, 2.0])
+    prob.init_derivatives()
+
+    prob.objective_forward(u)
+    prob.constraint_forward(u)
+
+    # With obj_factor=0, only constraint Hessian contributes
+    # Constraint is -exp(x), so Hessian is -exp(x)
+    lagrange = np.array([1.0, 1.0])
+    H = prob.hessian(obj_factor=0.0, lagrange=lagrange)
+
+    expected_H = np.diag([-np.exp(1.0), -np.exp(2.0)])
+    assert np.allclose(H.toarray(), expected_H)
+
+
+def test_hessian_two_constraints():
+    """Test Hessian with two constraints."""
+    x = cp.Variable(2)
+    obj = cp.sum(cp.log(x))
+    constraints = [
+        cp.log(x) >= 0,  # becomes -log(x) <= 0, Hessian: 1/x²
+        cp.exp(x) >= 0,  # becomes -exp(x) <= 0, Hessian: -exp(x)
+    ]
+
+    problem = cp.Problem(cp.Minimize(obj), constraints)
+    prob = C_problem(problem)
+
+    u = np.array([2.0, 3.0])
+    prob.init_derivatives()
+
+    prob.objective_forward(u)
+    prob.constraint_forward(u)
+
+    # Lagrange: [lambda1_1, lambda1_2, lambda2_1, lambda2_2]
+    # First constraint (-log): 2 elements, second constraint (-exp): 2 elements
+    lagrange = np.array([1.0, 1.0, 0.5, 0.5])
+    H = prob.hessian(obj_factor=1.0, lagrange=lagrange)
+
+    # Objective Hessian: diag([-1/x²])
+    # Constraint 1 (-log) Hessian: diag([1/x²]) (second derivative of -log is 1/x²)
+    # Constraint 2 (-exp) Hessian: diag([-exp(x)])
+    # Total: (-1/x²) + 1.0*(1/x²) + 0.5*(-exp(x)) = -0.5*exp(x)
+    expected_00 = -1.0/(2.0**2) + 1.0*(1.0/(2.0**2)) + 0.5*(-np.exp(2.0))
+    expected_11 = -1.0/(3.0**2) + 1.0*(1.0/(3.0**2)) + 0.5*(-np.exp(3.0))
+    expected_H = np.diag([expected_00, expected_11])
+
+    assert np.allclose(H.toarray(), expected_H)
+
+
+def test_hessian_repeated_evaluations():
+    """Test Hessian with repeated evaluations at different points."""
+    x = cp.Variable(2)
+    obj = cp.sum(cp.log(x))
+    problem = cp.Problem(cp.Minimize(obj))
+    prob = C_problem(problem)
+
+    prob.init_derivatives()
+
+    # First evaluation
+    u1 = np.array([1.0, 2.0])
+    prob.objective_forward(u1)
+    H1 = prob.hessian(obj_factor=1.0, lagrange=np.array([]))
+
+    # Second evaluation at different point
+    u2 = np.array([2.0, 4.0])
+    prob.objective_forward(u2)
+    H2 = prob.hessian(obj_factor=1.0, lagrange=np.array([]))
+
+    expected_H1 = np.diag([-1.0, -0.25])
+    expected_H2 = np.diag([-0.25, -0.0625])
+
+    assert np.allclose(H1.toarray(), expected_H1)
+    assert np.allclose(H2.toarray(), expected_H2)
+
+
 if __name__ == "__main__":
     test_single_constraint()
     print("test_single_constraint passed!")
@@ -199,4 +340,17 @@ def test_repeated_evaluations():
     print("test_larger_scale passed!")
     test_repeated_evaluations()
     print("test_repeated_evaluations passed!")
+
+    # Hessian tests
+    test_hessian_objective_only()
+    print("test_hessian_objective_only passed!")
+    test_hessian_with_constraint()
+    print("test_hessian_with_constraint passed!")
+    test_hessian_obj_factor_zero()
+    print("test_hessian_obj_factor_zero passed!")
+    test_hessian_two_constraints()
+    print("test_hessian_two_constraints passed!")
+    test_hessian_repeated_evaluations()
+    print("test_hessian_repeated_evaluations passed!")
+
     print("\nAll constrained tests passed!")
diff --git a/src/affine/constant.c b/src/affine/constant.c
index 577029b..145a075 100644
--- a/src/affine/constant.c
+++ b/src/affine/constant.c
@@ -21,6 +21,19 @@ static void eval_jacobian(expr *node)
     (void) node;
 }
 
+static void wsum_hess_init(expr *node)
+{
+    /* Constant Hessian is all zeros: n_vars x n_vars with 0 nonzeros. */
+    node->wsum_hess = new_csr_matrix(node->n_vars, node->n_vars, 0);
+}
+
+static void eval_wsum_hess(expr *node, const double *w)
+{
+    /* Constant Hessian is always zero - nothing to compute */
+    (void) node;
+    (void) w;
+}
+
 static bool is_affine(const expr *node)
 {
     (void) node;
@@ -35,6 +48,8 @@ expr *new_constant(int d1, int d2, int n_vars, const double *values)
     node->is_affine = is_affine;
     node->jacobian_init = jacobian_init;
     node->eval_jacobian = eval_jacobian;
+    node->wsum_hess_init = wsum_hess_init;
+    node->eval_wsum_hess = eval_wsum_hess;
 
     return node;
 }