From 3fba5813651cbeb2e44063c9194a89ccdffd27d8 Mon Sep 17 00:00:00 2001 From: Daniel Date: Fri, 9 Jan 2026 14:48:42 -0800 Subject: [PATCH 1/3] efficient assembly of problem hessian --- docs/problem_struct_design.md | 12 +-- include/problem.h | 6 +- include/utils/CSR_Matrix.h | 3 +- python/problem/jacobian.h | 2 +- src/affine/sum.c | 1 + src/affine/variable.c | 15 +++ src/problem.c | 177 +++++++++++++++++++++++++++++++--- src/utils/CSR_Matrix.c | 2 +- tests/all_tests.c | 1 + tests/problem/test_problem.h | 83 +++++++++++++++- 10 files changed, 276 insertions(+), 26 deletions(-) diff --git a/docs/problem_struct_design.md b/docs/problem_struct_design.md index 4c9144e..77ea2d3 100644 --- a/docs/problem_struct_design.md +++ b/docs/problem_struct_design.md @@ -37,7 +37,7 @@ typedef struct problem /* Pre-allocated storage */ double *constraint_values; double *gradient_values; /* Dense gradient array */ - CSR_Matrix *stacked_jac; + CSR_Matrix *jacobian; } problem; problem *new_problem(expr *objective, expr **constraints, int n_constraints); @@ -85,16 +85,16 @@ void problem_allocate(problem *prob, const double *u) } /* Allocate stacked jacobian with total_constraint_size rows */ - prob->stacked_jac = alloc_csr(prob->total_constraint_size, prob->n_vars, total_nnz); + prob->jacobian = alloc_csr(prob->total_constraint_size, prob->n_vars, total_nnz); /* Note: The actual nnz may be smaller after evaluation due to - * cancellations. Update stacked_jac->nnz after problem_jacobian(). */ + * cancellations. Update jacobian->nnz after problem_jacobian(). */ } ``` ### `free_problem` - Call `free_expr` on objective and all constraints (decrements refcount) -- Free allocated arrays and stacked_jac +- Free allocated arrays and jacobian ### `problem_forward` ```c @@ -127,7 +127,7 @@ double problem_forward(problem *prob, const double *u) - Stack CSR matrices vertically: - Total rows = `total_constraint_size` - Copy row pointers with offset, copy col indices and values -- Lazy allocate/reallocate `stacked_jac` based on total nnz +- Lazy allocate/reallocate `jacobian` based on total nnz --- @@ -236,6 +236,6 @@ class Problem: - **Two-phase init**: `new_problem` creates struct, `problem_allocate` allocates arrays - Constraint values array: size = `total_constraint_size` - Jacobian: initialize all constraint jacobians first, count total nnz, allocate CSR - - The allocated nnz may be a slight overestimate; update `stacked_jac->nnz` after evaluation + - The allocated nnz may be a slight overestimate; update `jacobian->nnz` after evaluation - **Hessian**: Deferred - not allocated in this design (to be added later) - **API**: Returns internal pointers (caller should NOT free) diff --git a/include/problem.h b/include/problem.h index 3159938..8f685b5 100644 --- a/include/problem.h +++ b/include/problem.h @@ -17,7 +17,10 @@ typedef struct problem double *gradient_values; /* Allocated by problem_init_derivatives */ - CSR_Matrix *stacked_jac; + CSR_Matrix *jacobian; + CSR_Matrix *lagrange_hessian; + int *hess_idx_map; /* Maps all wsum_hess nnz to lagrange_hessian (obj + + constraints) */ } problem; /* Retains objective and constraints (shared ownership with caller) */ @@ -29,5 +32,6 @@ double problem_objective_forward(problem *prob, const double *u); void problem_constraint_forward(problem *prob, const double *u); void problem_gradient(problem *prob); void problem_jacobian(problem *prob); +void problem_hessian(problem *prob, double obj_w, const double *w); #endif diff --git a/include/utils/CSR_Matrix.h b/include/utils/CSR_Matrix.h index aa4d7ab..6dbeb24 100644 --- a/include/utils/CSR_Matrix.h +++ b/include/utils/CSR_Matrix.h @@ -93,7 +93,8 @@ void sum_all_rows_csr_fill_sparsity_and_idx_map(const CSR_Matrix *A, CSR_Matrix // void sum_all_rows_csr_fill_values(const CSR_Matrix *A, CSR_Matrix *C, // const int *idx_map); -/* Fill accumulator for summing rows using precomputed idx_map for each nnz of A*/ +/* Fill accumulator for summing rows using precomputed idx_map for each nnz of A. + Must memset accumulator to zero before calling. */ void idx_map_accumulator(const CSR_Matrix *A, const int *idx_map, double *accumulator); diff --git a/python/problem/jacobian.h b/python/problem/jacobian.h index 69f240d..77fe795 100644 --- a/python/problem/jacobian.h +++ b/python/problem/jacobian.h @@ -33,7 +33,7 @@ static PyObject *py_problem_jacobian(PyObject *self, PyObject *args) problem_jacobian(prob); - CSR_Matrix *jac = prob->stacked_jac; + CSR_Matrix *jac = prob->jacobian; npy_intp nnz = jac->nnz; npy_intp m_plus_1 = jac->m + 1; diff --git a/src/affine/sum.c b/src/affine/sum.c index ac40679..ba3193b 100644 --- a/src/affine/sum.c +++ b/src/affine/sum.c @@ -100,6 +100,7 @@ static void eval_jacobian(expr *node) /* we have precomputed an idx map between the nonzeros of the child's jacobian and this node's jacobian, so we just accumulate accordingly */ + memset(node->jacobian->x, 0, node->jacobian->nnz * sizeof(double)); idx_map_accumulator(x->jacobian, ((sum_expr *) node)->idx_map, node->jacobian->x); } diff --git a/src/affine/variable.c b/src/affine/variable.c index 8ce8c6e..3e44dca 100644 --- a/src/affine/variable.c +++ b/src/affine/variable.c @@ -25,6 +25,19 @@ static void eval_jacobian(expr *node) (void) node; } +static void wsum_hess_init(expr *node) +{ + /* Variables have zero Hessian */ + node->wsum_hess = new_csr_matrix(node->n_vars, node->n_vars, 0); +} + +static void wsum_hess_eval(expr *node, const double *w) +{ + /* Variables have zero Hessian */ + (void) node; + (void) w; +} + static bool is_affine(const expr *node) { (void) node; @@ -39,6 +52,8 @@ expr *new_variable(int d1, int d2, int var_id, int n_vars) 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 = wsum_hess_eval; return node; } diff --git a/src/problem.c b/src/problem.c index 945d0a2..a297f02 100644 --- a/src/problem.c +++ b/src/problem.c @@ -2,6 +2,20 @@ #include #include +// TODO: move this into common file +#define MAX(a, b) ((a) > (b) ? (a) : (b)) + +/* forward declaration */ +static void problem_lagrange_hess_fill_sparsity(problem *prob, int *iwork); + +/* Helper function to compare integers for qsort */ +static int compare_int_asc(const void *a, const void *b) +{ + int ia = *((const int *) a); + int ib = *((const int *) b); + return (ia > ib) - (ia < ib); +} + problem *new_problem(expr *objective, expr **constraints, int n_constraints) { problem *prob = (problem *) calloc(1, sizeof(problem)); @@ -49,35 +63,129 @@ problem *new_problem(expr *objective, expr **constraints, int n_constraints) prob->gradient_values = (double *) calloc(prob->n_vars, sizeof(double)); /* Derivative structures allocated by problem_init_derivatives */ - prob->stacked_jac = NULL; + prob->jacobian = NULL; return prob; } void problem_init_derivatives(problem *prob) { - /* 1. Initialize objective jacobian */ + // ------------------------------------------------------------------------------- + // Jacobian structure + // ------------------------------------------------------------------------------- prob->objective->jacobian_init(prob->objective); - - /* 2. Initialize constraint jacobians and count total nnz */ - int total_nnz = 0; + int nnz = 0; for (int i = 0; i < prob->n_constraints; i++) { expr *c = prob->constraints[i]; c->jacobian_init(c); - total_nnz += c->jacobian->nnz; + nnz += c->jacobian->nnz; } - /* 3. Allocate stacked jacobian */ - if (prob->total_constraint_size > 0) + prob->jacobian = new_csr_matrix(prob->total_constraint_size, prob->n_vars, nnz); + + // ------------------------------------------------------------------------------- + // Lagrange Hessian structure + // ------------------------------------------------------------------------------- + prob->objective->wsum_hess_init(prob->objective); + nnz = prob->objective->wsum_hess->nnz; + + for (int i = 0; i < prob->n_constraints; i++) + { + prob->constraints[i]->wsum_hess_init(prob->constraints[i]); + nnz += prob->constraints[i]->wsum_hess->nnz; + } + + prob->lagrange_hessian = new_csr_matrix(prob->n_vars, prob->n_vars, nnz); + prob->hess_idx_map = (int *) malloc(nnz * sizeof(int)); + int *iwork = (int *) malloc(MAX(nnz, prob->n_vars) * sizeof(int)); + problem_lagrange_hess_fill_sparsity(prob, iwork); + free(iwork); +} + +static void problem_lagrange_hess_fill_sparsity(problem *prob, int *iwork) +{ + expr **constrs = prob->constraints; + int *cols = iwork; + int *col_to_pos = iwork; /* reused after qsort */ + int nnz = 0; + CSR_Matrix *H_obj = prob->objective->wsum_hess; + CSR_Matrix *H_c; + CSR_Matrix *H = prob->lagrange_hessian; + H->p[0] = 0; + + // ---------------------------------------------------------------------- + // Fill sparsity pattern + // ---------------------------------------------------------------------- + for (int row = 0; row < H->m; row++) { - prob->stacked_jac = - new_csr_matrix(prob->total_constraint_size, prob->n_vars, total_nnz); + /* gather columns from objective hessian */ + int count = H_obj->p[row + 1] - H_obj->p[row]; + memcpy(cols, H_obj->i + H_obj->p[row], count * sizeof(int)); + + /* gather columns from constraint hessians */ + for (int c_idx = 0; c_idx < prob->n_constraints; c_idx++) + { + H_c = constrs[c_idx]->wsum_hess; + int c_len = H_c->p[row + 1] - H_c->p[row]; + memcpy(cols + count, H_c->i + H_c->p[row], c_len * sizeof(int)); + count += c_len; + } + + /* find unique columns */ + qsort(cols, count, sizeof(int), compare_int_asc); + int prev_col = -1; + for (int j = 0; j < count; j++) + { + if (cols[j] != prev_col) + { + H->i[nnz] = cols[j]; + nnz++; + prev_col = cols[j]; + } + } + + H->p[row + 1] = nnz; } - /* TODO: 4. Initialize objective wsum_hess */ + H->nnz = nnz; + + // ---------------------------------------------------------------------- + // Build idx map + // ---------------------------------------------------------------------- + int idx_offset = 0; + + /* map objective hessian entries */ + for (int row = 0; row < H->m; row++) + { + for (int idx = H->p[row]; idx < H->p[row + 1]; idx++) + { + col_to_pos[H->i[idx]] = idx; + } + + for (int j = H_obj->p[row]; j < H_obj->p[row + 1]; j++) + { + prob->hess_idx_map[idx_offset++] = col_to_pos[H_obj->i[j]]; + } + } - /* TODO: 5. Initialize constraint wsum_hess */ + /* map constraint hessian entries */ + for (int c_idx = 0; c_idx < prob->n_constraints; c_idx++) + { + H_c = constrs[c_idx]->wsum_hess; + for (int row = 0; row < H->m; row++) + { + for (int idx = H->p[row]; idx < H->p[row + 1]; idx++) + { + col_to_pos[H->i[idx]] = idx; + } + + for (int j = H_c->p[row]; j < H_c->p[row + 1]; j++) + { + prob->hess_idx_map[idx_offset++] = col_to_pos[H_c->i[j]]; + } + } + } } void free_problem(problem *prob) @@ -87,7 +195,9 @@ void free_problem(problem *prob) /* Free allocated arrays */ free(prob->constraint_values); free(prob->gradient_values); - free_csr_matrix(prob->stacked_jac); + free_csr_matrix(prob->jacobian); + free_csr_matrix(prob->lagrange_hessian); + free(prob->hess_idx_map); /* Release expression references (decrements refcount) */ free_expr(prob->objective); @@ -141,7 +251,7 @@ void problem_gradient(problem *prob) void problem_jacobian(problem *prob) { - CSR_Matrix *stacked = prob->stacked_jac; + CSR_Matrix *stacked = prob->jacobian; /* Initialize row pointers */ stacked->p[0] = 0; @@ -149,6 +259,7 @@ void problem_jacobian(problem *prob) int row_offset = 0; int nnz_offset = 0; + /* TODO: here we eventually want to evaluate affine jacobians only once */ for (int i = 0; i < prob->n_constraints; i++) { expr *c = prob->constraints[i]; @@ -176,3 +287,41 @@ void problem_jacobian(problem *prob) /* Update actual nnz (may be less than allocated) */ stacked->nnz = nnz_offset; } + +void problem_hessian(problem *prob, double obj_w, const double *w) +{ + // ------------------------------------------------------------------------ + // evaluate hessian of objective and constraints + // ------------------------------------------------------------------------ + expr *obj = prob->objective; + expr **constrs = prob->constraints; + + obj->eval_wsum_hess(obj, &obj_w); + + int offset = 0; + for (int i = 0; i < prob->n_constraints; i++) + { + constrs[i]->eval_wsum_hess(constrs[i], w + offset); + offset += constrs[i]->size; + } + + // ------------------------------------------------------------------------ + // assemble Lagrange hessian using index map + // ------------------------------------------------------------------------ + CSR_Matrix *H = prob->lagrange_hessian; + int *idx_map = prob->hess_idx_map; + + /* zero out hessian before adding contribution from obj and constraints */ + memset(H->x, 0, H->nnz * sizeof(double)); + + /* accumulate objective function */ + idx_map_accumulator(obj->wsum_hess, idx_map, H->x); + offset = obj->wsum_hess->nnz; + + /* accumulate constraint functions */ + for (int i = 0; i < prob->n_constraints; i++) + { + idx_map_accumulator(constrs[i]->wsum_hess, idx_map + offset, H->x); + offset += constrs[i]->wsum_hess->nnz; + } +} diff --git a/src/utils/CSR_Matrix.c b/src/utils/CSR_Matrix.c index 5e9aed8..6b63c7e 100644 --- a/src/utils/CSR_Matrix.c +++ b/src/utils/CSR_Matrix.c @@ -675,7 +675,7 @@ void sum_evenly_spaced_rows_csr_fill_sparsity_and_idx_map(const CSR_Matrix *A, void idx_map_accumulator(const CSR_Matrix *A, const int *idx_map, double *accumulator) { - memset(accumulator, 0, A->nnz * sizeof(double)); + // memset(accumulator, 0, A->nnz * sizeof(double)); for (int row = 0; row < A->m; row++) { diff --git a/tests/all_tests.c b/tests/all_tests.c index ccfbd56..f59d647 100644 --- a/tests/all_tests.c +++ b/tests/all_tests.c @@ -172,6 +172,7 @@ int main(void) mu_run_test(test_problem_jacobian, tests_run); mu_run_test(test_problem_jacobian_multi, tests_run); mu_run_test(test_problem_constraint_forward, tests_run); + mu_run_test(test_problem_hessian, tests_run); printf("\n=== All %d tests passed ===\n", tests_run); diff --git a/tests/problem/test_problem.h b/tests/problem/test_problem.h index 7ea502b..d2b713e 100644 --- a/tests/problem/test_problem.h +++ b/tests/problem/test_problem.h @@ -174,7 +174,7 @@ const char *test_problem_jacobian(void) problem_constraint_forward(prob, u); problem_jacobian(prob); - CSR_Matrix *jac = prob->stacked_jac; + CSR_Matrix *jac = prob->jacobian; /* Check dimensions */ mu_assert("jac rows wrong", jac->m == 2); @@ -238,7 +238,7 @@ const char *test_problem_jacobian_multi(void) problem_constraint_forward(prob, u); problem_jacobian(prob); - CSR_Matrix *jac = prob->stacked_jac; + CSR_Matrix *jac = prob->jacobian; /* Check dimensions: 4 rows (2 + 2), 2 cols */ mu_assert("jac rows wrong", jac->m == 4); @@ -267,4 +267,83 @@ const char *test_problem_jacobian_multi(void) return 0; } +/* + * Test problem_hessian: Lagrange Hessian of: + * Objective: sum(log(x)) where x is 3x1 + * Constraint 1: exp(x) + * Constraint 2: sin(x) + * + * Evaluate at x = [1, 2, 3] with w = [1, 2, 3, 4, 5, 6], w_obj = 2 + * + * Lagrange function: + * L = 2*sum(log(x)) + w[0:3]^T exp(x) + w[3:6]^T sin(x) + * + * Hessian at x = [1, 2, 3]: + * H_obj = diag([-1/1^2, -1/2^2, -1/3^2]) = diag([-1, -0.25, -0.111111]) + * H_c1 = diag([exp(1), exp(2), exp(3)]) (element i has weight w[i]) + * H_c2 = diag([-sin(1), -sin(2), -sin(3)]) (element i has weight w[i+3]) + * + * Lagrange Hessian diagonal elements: + * H[0,0] = 2*(-1) + 1*exp(1) + 4*(-sin(1)) + * H[1,1] = 2*(-0.25) + 2*exp(2) + 5*(-sin(2)) + * H[2,2] = 2*(-0.111111) + 3*exp(3) + 6*(-sin(3)) + */ +const char *test_problem_hessian(void) +{ + int n_vars = 3; + + /* Shared variable */ + expr *x = new_variable(3, 1, 0, n_vars); + + /* Objective: sum(log(x)) */ + expr *log_obj = new_log(x); + expr *objective = new_sum(log_obj, -1); + + /* Constraint 1: exp(x) */ + expr *exp_c1 = new_exp(x); + + /* Constraint 2: sin(x) */ + expr *sin_c2 = new_sin(x); + + expr *constraints[2] = {exp_c1, sin_c2}; + + problem *prob = new_problem(objective, constraints, 2); + + double u[3] = {1.0, 2.0, 3.0}; + problem_init_derivatives(prob); + + /* Forward pass */ + problem_objective_forward(prob, u); + problem_constraint_forward(prob, u); + + /* Evaluate Lagrange Hessian */ + double w[6] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0}; + double w_obj = 2.0; + problem_hessian(prob, w_obj, w); + + CSR_Matrix *H = prob->lagrange_hessian; + + /* Check dimensions: 3x3 symmetric */ + mu_assert("H rows wrong", H->m == 3); + mu_assert("H cols wrong", H->n == 3); + + /* Compute expected diagonal values */ + double expected_H00 = 2.0 * (-1.0) + 1.0 * exp(1.0) + 4.0 * (-sin(1.0)); + double expected_H11 = 2.0 * (-0.25) + 2.0 * exp(2.0) + 5.0 * (-sin(2.0)); + double expected_H22 = 2.0 * (-1.0 / 9.0) + 3.0 * exp(3.0) + 6.0 * (-sin(3.0)); + + /* Since Hessian is diagonal, check the diagonal entries + * Row pointers should be [0, 1, 2, 3] */ + int expected_p[4] = {0, 1, 2, 3}; + int expected_i[3] = {0, 1, 2}; + double expected_x[3] = {expected_H00, expected_H11, expected_H22}; + mu_assert("H->p wrong", cmp_int_array(H->p, expected_p, 4)); + mu_assert("H->i wrong", cmp_int_array(H->i, expected_i, 3)); + mu_assert("H->x wrong", cmp_double_array(H->x, expected_x, 3)); + + free_problem(prob); + + return 0; +} + #endif /* TEST_PROBLEM_H */ From 2eb7554689e9e07900a12b2492901ea723ef8b1b Mon Sep 17 00:00:00 2001 From: William Zijie Zhang Date: Sat, 10 Jan 2026 02:18:35 -0500 Subject: [PATCH 2/3] Add Python Hessian bindings and fix memory management bugs - Add expr_retain() in all Python atom bindings to properly increment refcount when capsules are created. Without this, expression trees were getting double-freed when Python garbage collected capsules. - Add wsum_hess_init and eval_wsum_hess to constant.c. Constants were missing Hessian functions, causing NULL pointer crash when constraint expressions contained constants (e.g., "0 + -log(x)"). - Add Python bindings for problem_hessian (python/problem/hessian.h) - Add hessian() method to C_problem class in convert.py - Add comprehensive Hessian tests to test_constrained.py Co-Authored-By: Claude Opus 4.5 --- python/atoms/add.h | 1 + python/atoms/constant.h | 1 + python/atoms/exp.h | 1 + python/atoms/linear.h | 1 + python/atoms/log.h | 1 + python/atoms/neg.h | 1 + python/atoms/promote.h | 1 + python/atoms/sum.h | 1 + python/atoms/variable.h | 1 + python/bindings.c | 3 + python/convert.py | 19 ++++ python/problem/hessian.h | 76 +++++++++++++++ python/tests/test_constrained.py | 154 +++++++++++++++++++++++++++++++ src/affine/constant.c | 15 +++ 14 files changed, 276 insertions(+) create mode 100644 python/problem/hessian.h 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; } From 4adf4a86e7bc0bae16244f4cc54dea0b592864b1 Mon Sep 17 00:00:00 2001 From: William Zijie Zhang Date: Sat, 10 Jan 2026 02:20:49 -0500 Subject: [PATCH 3/3] Revert "Add Python Hessian bindings and fix memory management bugs" This reverts commit 2eb7554689e9e07900a12b2492901ea723ef8b1b. --- python/atoms/add.h | 1 - python/atoms/constant.h | 1 - python/atoms/exp.h | 1 - python/atoms/linear.h | 1 - python/atoms/log.h | 1 - python/atoms/neg.h | 1 - python/atoms/promote.h | 1 - python/atoms/sum.h | 1 - python/atoms/variable.h | 1 - python/bindings.c | 3 - python/convert.py | 19 ---- python/problem/hessian.h | 76 --------------- python/tests/test_constrained.py | 154 ------------------------------- src/affine/constant.c | 15 --- 14 files changed, 276 deletions(-) delete mode 100644 python/problem/hessian.h diff --git a/python/atoms/add.h b/python/atoms/add.h index 2f62a7a..4b86b79 100644 --- a/python/atoms/add.h +++ b/python/atoms/add.h @@ -24,7 +24,6 @@ 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 d5fcba3..da2af09 100644 --- a/python/atoms/constant.h +++ b/python/atoms/constant.h @@ -28,7 +28,6 @@ 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 9ac7fc5..9505daf 100644 --- a/python/atoms/exp.h +++ b/python/atoms/exp.h @@ -23,7 +23,6 @@ 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 c7afcae..0f85b91 100644 --- a/python/atoms/linear.h +++ b/python/atoms/linear.h @@ -54,7 +54,6 @@ 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 738b4e5..3d0314e 100644 --- a/python/atoms/log.h +++ b/python/atoms/log.h @@ -23,7 +23,6 @@ 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 3927b0c..bda69d0 100644 --- a/python/atoms/neg.h +++ b/python/atoms/neg.h @@ -23,7 +23,6 @@ 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 97e405c..fb258b4 100644 --- a/python/atoms/promote.h +++ b/python/atoms/promote.h @@ -24,7 +24,6 @@ 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 8d01604..8a8480e 100644 --- a/python/atoms/sum.h +++ b/python/atoms/sum.h @@ -24,7 +24,6 @@ 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 2dc1972..5438c44 100644 --- a/python/atoms/variable.h +++ b/python/atoms/variable.h @@ -17,7 +17,6 @@ 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 6e202d3..197390f 100644 --- a/python/bindings.c +++ b/python/bindings.c @@ -16,7 +16,6 @@ /* 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" @@ -54,8 +53,6 @@ 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 669f269..92920cb 100644 --- a/python/convert.py +++ b/python/convert.py @@ -159,22 +159,3 @@ 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 deleted file mode 100644 index 41c2788..0000000 --- a/python/problem/hessian.h +++ /dev/null @@ -1,76 +0,0 @@ -#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 4409b13..9574d3d 100644 --- a/python/tests/test_constrained.py +++ b/python/tests/test_constrained.py @@ -186,147 +186,6 @@ 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!") @@ -340,17 +199,4 @@ def test_hessian_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 145a075..577029b 100644 --- a/src/affine/constant.c +++ b/src/affine/constant.c @@ -21,19 +21,6 @@ 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; @@ -48,8 +35,6 @@ 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; }