Add parameter support to C diff engine

Add parameter node type and parameter-aware variants of scalar mult,
vector mult, and left matmul. Parameters store an offset into a global
theta vector and can be updated via problem_update_params without
rebuilding the expression tree.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: Transurgeon

include/affine.h

@@ -34,6 +34,7 @@ expr *new_trace(expr *child);
 
 expr *new_constant(int d1, int d2, int n_vars, const double *values);
 expr *new_variable(int d1, int d2, int var_id, int n_vars);
+expr *new_parameter(int d1, int d2, int param_id, int n_vars);
 
 expr *new_index(expr *child, int d1, int d2, const int *indices, int n_idxs);
 expr *new_reshape(expr *child, int d1, int d2);

include/bivariate.h

@@ -42,4 +42,13 @@ expr *new_const_scalar_mult(double a, expr *child);
 /* Constant vector elementwise multiplication: a ∘ f(x) where a is constant */
 expr *new_const_vector_mult(const double *a, expr *child);
 
+/* Left matrix multiplication with parameter source: P @ f(x) where P is a parameter */
+expr *new_left_param_matmul(expr *param_node, expr *u, int A_m, int A_n);
+
+/* Parameter scalar multiplication: p * f(x) where p is a parameter */
+expr *new_param_scalar_mult(expr *param_node, expr *child);
+
+/* Parameter vector elementwise multiplication: p ∘ f(x) where p is a parameter */
+expr *new_param_vector_mult(expr *param_node, expr *child);
+
 #endif /* BIVARIATE_H */

include/problem.h

@@ -59,6 +59,11 @@ typedef struct problem
      * hessian are called */
     bool jacobian_called;
 
+    /* Parameter tracking for fast parameter updates */
+    expr **param_nodes;     /* weak references to parameter nodes in tree */
+    int n_param_nodes;
+    int n_params;           /* total scalar parameters */
+
     /* Statistics for performance measurement */
     Diff_engine_stats stats;
     bool verbose;
@@ -78,4 +83,9 @@ void problem_gradient(problem *prob);
 void problem_jacobian(problem *prob);
 void problem_hessian(problem *prob, double obj_w, const double *w);
 
+/* Parameter support */
+void problem_register_params(problem *prob, expr **param_nodes,
+                             int n_param_nodes, int n_params);
+void problem_update_params(problem *prob, const double *theta);
+
 #endif

include/subexpr.h

@@ -25,6 +25,13 @@
 /* Forward declaration */
 struct int_double_pair;
 
+/* Parameter node: like constant but with updatable values via problem_update_params */
+typedef struct parameter_expr
+{
+    expr base;
+    int param_id; /* offset into global theta vector */
+} parameter_expr;
+
 /* Type-specific expression structures that "inherit" from expr */
 
 /* Linear operator: y = A * x + b */
@@ -110,6 +117,8 @@ typedef struct left_matmul_expr
     CSR_Matrix *A;
     CSR_Matrix *AT;
     CSC_Matrix *CSC_work;
+    expr *param_source; /* if non-NULL, refresh A/AT values from param_source->value */
+    int src_m, src_n;   /* original (non-block-diag) matrix dimensions */
 } left_matmul_expr;
 
 /* Right matrix multiplication: y = f(x) * A where f(x) is an expression.
@@ -128,13 +137,15 @@ typedef struct const_scalar_mult_expr
 {
     expr base;
     double a;
+    expr *param_source; /* if non-NULL, read a from param_source->value[0] */
 } const_scalar_mult_expr;
 
 /* Constant vector elementwise multiplication: y = a \circ child for constant a */
 typedef struct const_vector_mult_expr
 {
     expr base;
-    double *a; /* length equals node->size */
+    double *a;          /* length equals node->size */
+    expr *param_source; /* if non-NULL, use param_source->value instead of a */
 } const_vector_mult_expr;
 
 /* Index/slicing: y = child[indices] where indices is a list of flat positions */

src/affine/parameter.c

@@ -0,0 +1,74 @@
+/*
+ * Copyright 2026 Daniel Cederberg and William Zhang
+ *
+ * This file is part of the DNLP-differentiation-engine project.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/* Parameter leaf node: behaviorally identical to constant (zero derivatives
+   w.r.t. variables), but its values are updatable via problem_update_params.
+   This allows re-solving with different parameter values without rebuilding
+   the expression tree. */
+
+#include "affine.h"
+#include "subexpr.h"
+#include <stdlib.h>
+#include <string.h>
+
+static void forward(expr *node, const double *u)
+{
+    /* Values are set by problem_update_params, not by forward pass */
+    (void)node;
+    (void)u;
+}
+
+static void jacobian_init(expr *node)
+{
+    /* Parameter jacobian is all zeros: size x n_vars with 0 nonzeros */
+    node->jacobian = new_csr_matrix(node->size, node->n_vars, 0);
+}
+
+static void eval_jacobian(expr *node)
+{
+    /* Parameter jacobian never changes */
+    (void)node;
+}
+
+static void wsum_hess_init(expr *node)
+{
+    /* Parameter Hessian is all zeros */
+    node->wsum_hess = new_csr_matrix(node->n_vars, node->n_vars, 0);
+}
+
+static void eval_wsum_hess(expr *node, const double *w)
+{
+    (void)node;
+    (void)w;
+}
+
+static bool is_affine(const expr *node)
+{
+    (void)node;
+    return true;
+}
+
+expr *new_parameter(int d1, int d2, int param_id, int n_vars)
+{
+    parameter_expr *pnode = (parameter_expr *)calloc(1, sizeof(parameter_expr));
+    init_expr(&pnode->base, d1, d2, n_vars, forward, jacobian_init, eval_jacobian,
+              is_affine, wsum_hess_init, eval_wsum_hess, NULL);
+    pnode->param_id = param_id;
+    /* values will be populated by problem_update_params */
+    return &pnode->base;
+}

src/bivariate/const_scalar_mult.c

@@ -24,6 +24,11 @@
 
 /* Constant scalar multiplication: y = a * child where a is a constant double */
 
+static inline double get_scalar(const const_scalar_mult_expr *sn)
+{
+    return sn->param_source ? sn->param_source->value[0] : sn->a;
+}
+
 static void forward(expr *node, const double *u)
 {
     expr *child = node->left;
@@ -32,7 +37,7 @@ static void forward(expr *node, const double *u)
     child->forward(child, u);
 
     /* local forward pass: multiply each element by scalar a */
-    double a = ((const_scalar_mult_expr *) node)->a;
+    double a = get_scalar((const_scalar_mult_expr *) node);
     for (int i = 0; i < node->size; i++)
     {
         node->value[i] = a * child->value[i];
@@ -55,7 +60,7 @@ static void jacobian_init(expr *node)
 static void eval_jacobian(expr *node)
 {
     expr *child = node->left;
-    double a = ((const_scalar_mult_expr *) node)->a;
+    double a = get_scalar((const_scalar_mult_expr *) node);
 
     /* evaluate child */
     child->eval_jacobian(child);
@@ -85,7 +90,7 @@ static void eval_wsum_hess(expr *node, const double *w)
     expr *x = node->left;
     x->eval_wsum_hess(x, w);
 
-    double a = ((const_scalar_mult_expr *) node)->a;
+    double a = get_scalar((const_scalar_mult_expr *) node);
     for (int j = 0; j < x->wsum_hess->nnz; j++)
     {
         node->wsum_hess->x[j] = a * x->wsum_hess->x[j];
@@ -108,7 +113,25 @@ expr *new_const_scalar_mult(double a, expr *child)
               eval_jacobian, is_affine, wsum_hess_init, eval_wsum_hess, NULL);
     node->left = child;
     mult_node->a = a;
+    mult_node->param_source = NULL;
+    expr_retain(child);
+
+    return node;
+}
+
+expr *new_param_scalar_mult(expr *param_node, expr *child)
+{
+    const_scalar_mult_expr *mult_node =
+        (const_scalar_mult_expr *) calloc(1, sizeof(const_scalar_mult_expr));
+    expr *node = &mult_node->base;
+
+    init_expr(node, child->d1, child->d2, child->n_vars, forward, jacobian_init,
+              eval_jacobian, is_affine, wsum_hess_init, eval_wsum_hess, NULL);
+    node->left = child;
+    mult_node->a = param_node->value[0]; /* initial value */
+    mult_node->param_source = param_node;
     expr_retain(child);
+    expr_retain(param_node);
 
     return node;
 }

src/bivariate/const_vector_mult.c

@@ -23,10 +23,15 @@
 
 /* Constant vector elementwise multiplication: y = a \circ child */
 
+static inline const double *get_vector(const const_vector_mult_expr *vn)
+{
+    return vn->param_source ? vn->param_source->value : vn->a;
+}
+
 static void forward(expr *node, const double *u)
 {
     expr *child = node->left;
-    const double *a = ((const_vector_mult_expr *) node)->a;
+    const double *a = get_vector((const_vector_mult_expr *) node);
 
     /* child's forward pass */
     child->forward(child, u);
@@ -54,7 +59,7 @@ static void jacobian_init(expr *node)
 static void eval_jacobian(expr *node)
 {
     expr *x = node->left;
-    const double *a = ((const_vector_mult_expr *) node)->a;
+    const double *a = get_vector((const_vector_mult_expr *) node);
 
     /* evaluate x */
     x->eval_jacobian(x);
@@ -87,7 +92,7 @@ static void wsum_hess_init(expr *node)
 static void eval_wsum_hess(expr *node, const double *w)
 {
     expr *x = node->left;
-    const double *a = ((const_vector_mult_expr *) node)->a;
+    const double *a = get_vector((const_vector_mult_expr *) node);
 
     /* scale weights w by a */
     for (int i = 0; i < node->size; i++)
@@ -128,6 +133,39 @@ expr *new_const_vector_mult(const double *a, expr *child)
     /* copy a vector */
     vnode->a = (double *) malloc(child->size * sizeof(double));
     memcpy(vnode->a, a, child->size * sizeof(double));
+    vnode->param_source = NULL;
+
+    return node;
+}
+
+static void free_param_type_data(expr *node)
+{
+    const_vector_mult_expr *vnode = (const_vector_mult_expr *) node;
+    /* a is not owned when param_source is set */
+    free(vnode->a);
+    if (vnode->param_source)
+    {
+        free_expr(vnode->param_source);
+    }
+}
+
+expr *new_param_vector_mult(expr *param_node, expr *child)
+{
+    const_vector_mult_expr *vnode =
+        (const_vector_mult_expr *) calloc(1, sizeof(const_vector_mult_expr));
+    expr *node = &vnode->base;
+
+    init_expr(node, child->d1, child->d2, child->n_vars, forward, jacobian_init,
+              eval_jacobian, is_affine, wsum_hess_init, eval_wsum_hess,
+              free_param_type_data);
+    node->left = child;
+    expr_retain(child);
+
+    /* Still allocate a copy for initial values (used before first update_params) */
+    vnode->a = (double *) malloc(child->size * sizeof(double));
+    memcpy(vnode->a, param_node->value, child->size * sizeof(double));
+    vnode->param_source = param_node;
+    expr_retain(param_node);
 
     return node;
 }