diff --git a/include/expr.h b/include/expr.h
index fbe5a675..85e91c89 100644
--- a/include/expr.h
+++ b/include/expr.h
@@ -23,6 +23,7 @@
 #include "utils/matrix.h"
 #include <stdbool.h>
 #include <stddef.h>
+#include <stdint.h>
 #include <string.h>
 
 #define JAC_IDXS_NOT_SET -1
@@ -39,6 +40,7 @@ typedef void (*local_jacobian_fn)(struct expr *node, double *out);
 typedef void (*local_wsum_hess_fn)(struct expr *node, double *out, const double *w);
 typedef bool (*is_affine_fn)(const struct expr *node);
 typedef void (*free_type_data_fn)(struct expr *node);
+typedef void (*set_needs_refresh_children_fn)(struct expr *node);
 
 /* Workspace for derivative computation */
 typedef struct
@@ -48,10 +50,19 @@ typedef struct
     CSC_matrix *jacobian_csc;
     int *csc_work; /* for CSR_matrix-CSC_matrix conversion */
 
-    /* jacobian_csc_filled is only used for affine functions to avoid redundant
-       conversions. Could become relevant for non-affine functions if we start
-       supporting common subexpressions on the Python side. */
-    bool jacobian_csc_filled;
+    /* jacobian->values_version that the jacobian_csc mirror reflects;
+       expr_refresh_jacobian_csc refills iff it differs. */
+    uint64_t jacobian_csc_seen;
+
+    /* node->is_affine(node), computed once by jacobian_init (affinity is
+       structural, and the recursive is_affine is too costly per eval). */
+    bool is_affine_cached;
+
+    /* True once eval_jacobian has run this parameter epoch; cleared by
+       expr_set_needs_refresh. Only consulted for affine nodes, where it
+       lets the eval_jacobian wrapper skip the values_version bump (same
+       role the old jacobian_csc_filled latch played). */
+    bool jacobian_evaluated;
     double *local_jac_diag; /* cached f'(g(x)) diagonal */
     matrix *hess_term1;     /* Jg^T D Jg workspace */
     matrix *hess_term2;     /* child wsum_hess workspace */
@@ -76,8 +87,8 @@ typedef struct expr
     forward_fn forward;
     jacobian_init_fn jacobian_init_impl;
     wsum_hess_init_fn wsum_hess_init_impl;
-    eval_jacobian_fn eval_jacobian;
-    wsum_hess_fn eval_wsum_hess;
+    eval_jacobian_fn eval_jacobian_impl;
+    wsum_hess_fn eval_wsum_hess_impl;
 
     // ------------------------------------------------------------------------
     //                      other things
@@ -86,7 +97,11 @@ typedef struct expr
     local_jacobian_fn local_jacobian;   /* used by elementwise univariate atoms*/
     local_wsum_hess_fn local_wsum_hess; /* used by elementwise univariate atoms*/
     free_type_data_fn free_type_data;   /* Cleanup for type-specific fields */
-    Expr_Work *work;                    /* derivative workspace */
+    /* Recursion hook for expr_set_needs_refresh: atoms holding children
+       outside left/right (hstack's args[]) set this so the parameter-refresh
+       walk reaches them. NULL for binary/unary atoms. */
+    set_needs_refresh_children_fn set_needs_refresh_children;
+    Expr_Work *work; /* derivative workspace */
     /* Set to true on all nodes by problem_update_params() via
        expr_set_needs_refresh(). Atoms that cache parameter data
        (e.g. left_matmul_dense) check this flag before their forward
@@ -111,6 +126,15 @@ void free_expr(expr *node);
 void jacobian_init(expr *node);
 void wsum_hess_init(expr *node);
 
+/* Eval wrappers: run the atom's eval_*_impl and bump the output matrix's
+ * values_version so version-guarded caches (CSC mirrors, spd CSR views)
+ * refresh. Always call these instead of the impl slots. */
+void eval_jacobian(expr *node);
+void eval_wsum_hess(expr *node, const double *w);
+
+/* Refresh work->jacobian_csc from node->jacobian iff its values changed. */
+void expr_refresh_jacobian_csc(expr *node);
+
 /* Initialize CSC_matrix form of the Jacobian from the CSR_matrix Jacobian.
  * Must be called after jacobian_init. */
 void jacobian_csc_init(expr *node);
diff --git a/include/utils/matrix.h b/include/utils/matrix.h
index 55d99a1c..1aad9322 100644
--- a/include/utils/matrix.h
+++ b/include/utils/matrix.h
@@ -21,6 +21,7 @@
 #include "CSC_matrix.h"
 #include "CSR_matrix.h"
 #include <stdbool.h>
+#include <stdint.h>
 
 /* Broadcast shape used by the broadcast atom and its vtable methods. */
 typedef enum
@@ -77,7 +78,8 @@ typedef void (*matrix_transpose_fill_values_fn)(const matrix *A, matrix *AT);
 typedef CSR_matrix *(*matrix_to_csr_fn)(matrix *A);
 
 /* Refresh any internal caches (e.g. a CSC_matrix mirror) so subsequent ATA /
-   ATDA calls reflect the current values. */
+   ATDA calls reflect the current values. Version-guarded: a no-op when the
+   cache already matches values_version, so it is cheap to call when fresh. */
 typedef void (*matrix_refresh_csc_values_fn)(matrix *A);
 
 /* Allocate C = A[indices, :] */
@@ -128,6 +130,14 @@ struct matrix
     bool is_permuted_dense;
     bool is_stacked_pd;
 
+    /* Monotone counter bumped whenever the matrix's values change. Consumers
+       that mirror the values into a cache (CSC mirror, CSR view, ...) record
+       the version they last saw and refresh iff it differs. Code that writes
+       x directly must call matrix_values_changed on the OWNER of the buffer —
+       aliased children (spd blocks, cache views) have no version of their
+       own. */
+    uint64_t values_version;
+
     /* Operator ops */
     matrix_block_left_mult_vec_fn block_left_mult_vec;
     matrix_block_left_mult_sparsity_fn block_left_mult_sparsity;
@@ -160,6 +170,12 @@ struct matrix
     matrix_free_fn free_fn;
 };
 
+/* Notify the library after writing A->x directly. */
+static inline void matrix_values_changed(matrix *A)
+{
+    A->values_version++;
+}
+
 /* Free helper */
 static inline void free_matrix(matrix *m)
 {
diff --git a/include/utils/sparse_matrix.h b/include/utils/sparse_matrix.h
index 5ef7788e..85ee3f35 100644
--- a/include/utils/sparse_matrix.h
+++ b/include/utils/sparse_matrix.h
@@ -28,6 +28,7 @@ typedef struct sparse_matrix
     matrix base;
     CSR_matrix *csr;
     CSC_matrix *csc_cache;
+    uint64_t csc_seen; /* base.values_version the csc_cache values reflect */
     int *csc_iwork;
     int *transpose_iwork; /* sized csr->n; allocated by sparse_transpose_alloc
                              on the output sm and reused by
diff --git a/include/utils/stacked_pd.h b/include/utils/stacked_pd.h
index 568e3824..37ca3a16 100644
--- a/include/utils/stacked_pd.h
+++ b/include/utils/stacked_pd.h
@@ -48,6 +48,7 @@ typedef struct stacked_pd
 
     /* lazily built CSR view */
     CSR_matrix *csr_cache;
+    uint64_t csr_seen; /* base.values_version the csr_cache values reflect */
 
     /* Private permuted_dense scratch owned by the kernel that produced
        this spd. Allocated by the producing _alloc, used (without
diff --git a/src/atoms/affine/add.c b/src/atoms/affine/add.c
index 7f434e0b..0e02a60a 100644
--- a/src/atoms/affine/add.c
+++ b/src/atoms/affine/add.c
@@ -50,11 +50,11 @@ static void jacobian_init_impl(expr *node)
     sum_matrices_alloc(node->left->jacobian, node->right->jacobian, node->jacobian);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     /* evaluate children's jacobians */
-    node->left->eval_jacobian(node->left);
-    node->right->eval_jacobian(node->right);
+    eval_jacobian(node->left);
+    eval_jacobian(node->right);
 
     /* sum children's jacobians */
     sum_matrices_fill_values(node->left->jacobian, node->right->jacobian,
@@ -76,11 +76,11 @@ static void wsum_hess_init_impl(expr *node)
                        node->wsum_hess);
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     /* evaluate children's wsum_hess */
-    node->left->eval_wsum_hess(node->left, w);
-    node->right->eval_wsum_hess(node->right, w);
+    eval_wsum_hess(node->left, w);
+    eval_wsum_hess(node->right, w);
 
     /* sum children's wsum_hess */
     sum_matrices_fill_values(node->left->wsum_hess, node->right->wsum_hess,
@@ -97,7 +97,8 @@ expr *new_add(expr *left, expr *right)
     assert(left->d1 == right->d1 && left->d2 == right->d2);
     expr *node = (expr *) sp_calloc(1, sizeof(expr));
     init_expr(node, left->d1, left->d2, left->n_vars, forward, jacobian_init_impl,
-              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL);
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, NULL);
     node->left = left;
     node->right = right;
     expr_retain(left);
diff --git a/src/atoms/affine/broadcast.c b/src/atoms/affine/broadcast.c
index f5792f83..1a1cccc2 100644
--- a/src/atoms/affine/broadcast.c
+++ b/src/atoms/affine/broadcast.c
@@ -77,9 +77,9 @@ static void jacobian_init_impl(expr *node)
         x->jacobian->broadcast_alloc(x->jacobian, bcast->type, node->d1, node->d2);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
-    node->left->eval_jacobian(node->left);
+    eval_jacobian(node->left);
 
     /* fill values into the preallocated output. */
     broadcast_expr *bcast = (broadcast_expr *) node;
@@ -99,7 +99,7 @@ static void wsum_hess_init_impl(expr *node)
     node->work->dwork = sp_malloc(node->size * sizeof(double));
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     broadcast_expr *bcast = (broadcast_expr *) node;
     expr *x = node->left;
@@ -139,7 +139,7 @@ static void eval_wsum_hess(expr *node, const double *w)
         }
     }
 
-    x->eval_wsum_hess(x, node->work->dwork);
+    eval_wsum_hess(x, node->work->dwork);
     memcpy(node->wsum_hess->x, x->wsum_hess->x,
            node->wsum_hess->nnz * sizeof(double));
 }
@@ -183,7 +183,8 @@ expr *new_broadcast(expr *child, int d1, int d2)
     //                  initialize the rest of the expression
     // --------------------------------------------------------------------------
     init_expr(node, d1, d2, child->n_vars, forward, jacobian_init_impl,
-              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL);
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, NULL);
     node->left = child;
     expr_retain(child);
     bcast->type = type;
diff --git a/src/atoms/affine/convolve.c b/src/atoms/affine/convolve.c
index 0ce44e49..c9ed047d 100644
--- a/src/atoms/affine/convolve.c
+++ b/src/atoms/affine/convolve.c
@@ -91,12 +91,12 @@ static void jacobian_init_impl(expr *node)
         new_sparse_matrix(csr_csc_matmul_alloc(cnode->T, cnode->Jchild_CSC));
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     expr *child = node->left;
     convolve_expr *cnode = (convolve_expr *) node;
 
-    child->eval_jacobian(child);
+    eval_jacobian(child);
 
     /* J = T @ J_child */
     csr_to_csc_fill_values(child->jacobian->to_csr(child->jacobian),
@@ -115,7 +115,7 @@ static void wsum_hess_init_impl(expr *node)
     node->work->dwork = (double *) sp_malloc(cnode->n * sizeof(double));
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     expr *child = node->left;
     convolve_expr *cnode = (convolve_expr *) node;
@@ -133,7 +133,7 @@ static void eval_wsum_hess(expr *node, const double *w)
         w_prime[j] = sum;
     }
 
-    child->eval_wsum_hess(child, w_prime);
+    eval_wsum_hess(child, w_prime);
     memcpy(node->wsum_hess->x, child->wsum_hess->x,
            node->wsum_hess->nnz * sizeof(double));
 }
@@ -181,8 +181,8 @@ expr *new_convolve(expr *param_node, expr *child)
     convolve_expr *cnode = (convolve_expr *) sp_calloc(1, sizeof(convolve_expr));
     expr *node = &cnode->base;
     init_expr(node, d1, d2, child->n_vars, forward, jacobian_init_impl,
-              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess,
-              free_type_data);
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, free_type_data);
     node->left = child;
     expr_retain(child);
 
diff --git a/src/atoms/affine/diag_vec.c b/src/atoms/affine/diag_vec.c
index d5d74a04..955ca88c 100644
--- a/src/atoms/affine/diag_vec.c
+++ b/src/atoms/affine/diag_vec.c
@@ -54,9 +54,9 @@ static void jacobian_init_impl(expr *node)
     node->jacobian = x->jacobian->diag_vec_alloc(x->jacobian);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
-    node->left->eval_jacobian(node->left);
+    eval_jacobian(node->left);
 
     /* fill the diagonal rows of the preallocated output. */
     node->left->jacobian->diag_vec_fill_values(node->left->jacobian, node->jacobian);
@@ -77,7 +77,7 @@ static void wsum_hess_init_impl(expr *node)
     node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     expr *x = node->left;
     int n = x->size;
@@ -89,7 +89,7 @@ static void eval_wsum_hess(expr *node, const double *w)
     }
 
     /* Evaluate child's Hessian with extracted weights */
-    x->eval_wsum_hess(x, node->work->dwork);
+    eval_wsum_hess(x, node->work->dwork);
     memcpy(node->wsum_hess->x, x->wsum_hess->x,
            node->wsum_hess->nnz * sizeof(double));
 }
@@ -107,8 +107,9 @@ expr *new_diag_vec(expr *child)
     /* n is the number of elements (works for both row and column vectors) */
     int n = child->size;
     expr *node = (expr *) sp_calloc(1, sizeof(expr));
-    init_expr(node, n, n, child->n_vars, forward, jacobian_init_impl, eval_jacobian,
-              is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL);
+    init_expr(node, n, n, child->n_vars, forward, jacobian_init_impl,
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, NULL);
     node->left = child;
     expr_retain(child);
 
diff --git a/src/atoms/affine/hstack.c b/src/atoms/affine/hstack.c
index 9c4795a4..71554cf7 100644
--- a/src/atoms/affine/hstack.c
+++ b/src/atoms/affine/hstack.c
@@ -87,7 +87,7 @@ static void jacobian_init_impl(expr *node)
     node->jacobian = new_sparse_matrix(A);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     hstack_expr *hnode = (hstack_expr *) node;
     int cursor = 0;
@@ -95,7 +95,7 @@ static void eval_jacobian(expr *node)
     for (int i = 0; i < hnode->n_args; i++)
     {
         expr *child = hnode->args[i];
-        child->eval_jacobian(child);
+        eval_jacobian(child);
         /* to_csr needed for stacked_pd */
         CSR_matrix *child_csr = child->jacobian->to_csr(child->jacobian);
         memcpy(node->jacobian->x + cursor, child_csr->x,
@@ -148,7 +148,7 @@ static void wsum_hess_eval(expr *node, const double *w)
     for (int i = 0; i < hnode->n_args; i++)
     {
         expr *child = hnode->args[i];
-        child->eval_wsum_hess(child, w + row_offset);
+        eval_wsum_hess(child, w + row_offset);
         copy_CSR_matrix(H, hnode->CSR_work);
         sum_csr_fill_values(hnode->CSR_work,
                             child->wsum_hess->to_csr(child->wsum_hess), H);
@@ -170,6 +170,17 @@ static bool is_affine(const expr *node)
     return true;
 }
 
+/* Children live in args[], not left/right, so the parameter-refresh walk
+   needs this hook to reach them. */
+static void set_needs_refresh_children(expr *node)
+{
+    hstack_expr *hnode = (hstack_expr *) node;
+    for (int i = 0; i < hnode->n_args; i++)
+    {
+        expr_set_needs_refresh(hnode->args[i]);
+    }
+}
+
 static void free_type_data(expr *node)
 {
     hstack_expr *hnode = (hstack_expr *) node;
@@ -199,9 +210,11 @@ expr *new_hstack(expr **args, int n_args, int n_vars)
     hstack_expr *hnode = (hstack_expr *) sp_calloc(1, sizeof(hstack_expr));
     expr *node = &hnode->base;
     init_expr(node, args[0]->d1, d2, n_vars, forward, jacobian_init_impl,
-              eval_jacobian, is_affine, wsum_hess_init_impl, wsum_hess_eval,
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl, wsum_hess_eval,
               free_type_data);
 
+    node->set_needs_refresh_children = set_needs_refresh_children;
+
     /* Set type-specific fields (deep copy args array) */
     hnode->args = (expr **) sp_calloc(n_args, sizeof(expr *));
     hnode->n_args = n_args;
diff --git a/src/atoms/affine/index.c b/src/atoms/affine/index.c
index c0de149f..28baa1d0 100644
--- a/src/atoms/affine/index.c
+++ b/src/atoms/affine/index.c
@@ -70,11 +70,11 @@ static void jacobian_init_impl(expr *node)
         x->jacobian->index_alloc(x->jacobian, idx->indices, idx->n_idxs);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     expr *x = node->left;
     index_expr *idx = (index_expr *) node;
-    x->eval_jacobian(x);
+    eval_jacobian(x);
 
     /* copy values of the selected rows into the preallocated output */
     x->jacobian->index_fill_values(x->jacobian, idx->indices, idx->n_idxs,
@@ -100,7 +100,7 @@ static void wsum_hess_init_impl(expr *node)
     node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     expr *x = node->left;
     index_expr *idx = (index_expr *) node;
@@ -124,7 +124,7 @@ static void eval_wsum_hess(expr *node, const double *w)
     }
 
     /* evalute hessian of child */
-    x->eval_wsum_hess(x, node->work->dwork);
+    eval_wsum_hess(x, node->work->dwork);
     memcpy(node->wsum_hess->x, x->wsum_hess->x,
            node->wsum_hess->nnz * sizeof(double));
 }
@@ -152,8 +152,8 @@ expr *new_index(expr *child, int d1, int d2, const int *indices, int n_idxs)
     expr *node = &idx->base;
 
     init_expr(node, d1, d2, child->n_vars, forward, jacobian_init_impl,
-              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess,
-              free_type_data);
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, free_type_data);
 
     node->left = child;
     expr_retain(child);
diff --git a/src/atoms/affine/kron.c b/src/atoms/affine/kron.c
index d3e97940..b0dc5e17 100644
--- a/src/atoms/affine/kron.c
+++ b/src/atoms/affine/kron.c
@@ -107,12 +107,12 @@ static void jacobian_init_impl(expr *node)
     node->jacobian = new_sparse_matrix(Jk);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     expr *child = node->left;
     kron_expr *knode = (kron_expr *) node;
 
-    child->eval_jacobian(child);
+    eval_jacobian(child);
 
     /* Sparsity is fixed after jacobian_init, so the row offsets still align;
        refill active rows as scale * child-row-values. */
@@ -146,7 +146,7 @@ static void wsum_hess_init_impl(expr *node)
     node->work->dwork = (double *) sp_malloc(child->size * sizeof(double));
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     expr *child = node->left;
     kron_expr *knode = (kron_expr *) node;
@@ -166,7 +166,7 @@ static void eval_wsum_hess(expr *node, const double *w)
         }
     }
 
-    child->eval_wsum_hess(child, w_prime);
+    eval_wsum_hess(child, w_prime);
     memcpy(node->wsum_hess->x, child->wsum_hess->x,
            node->wsum_hess->nnz * sizeof(double));
 }
@@ -198,8 +198,8 @@ static kron_expr *new_kron_common(expr *param_node, expr *child, int p, int q, i
     kron_expr *knode = (kron_expr *) sp_calloc(1, sizeof(kron_expr));
     expr *node = &knode->base;
     init_expr(node, p * r, q * s, child->n_vars, forward, jacobian_init_impl,
-              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess,
-              free_type_data);
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, free_type_data);
     node->left = child;
     expr_retain(child);
 
diff --git a/src/atoms/affine/left_matmul.c b/src/atoms/affine/left_matmul.c
index a983577a..e4680f6f 100644
--- a/src/atoms/affine/left_matmul.c
+++ b/src/atoms/affine/left_matmul.c
@@ -139,7 +139,7 @@ static void eval_jacobian_dense(expr *node)
     /* evaluate jacobian of child */
     left_matmul_expr *lnode = (left_matmul_expr *) node;
     expr *x = node->left;
-    x->eval_jacobian(x);
+    eval_jacobian(x);
 
     /* must refresh CSC cache if x->jacobian is sparse_matrix */
     x->jacobian->refresh_csc_values(x->jacobian);
@@ -180,7 +180,7 @@ static void eval_jacobian_sparse(expr *node)
     /* evaluate jacobian of child */
     left_matmul_expr *lnode = (left_matmul_expr *) node;
     expr *x = node->left;
-    x->eval_jacobian(x);
+    eval_jacobian(x);
 
     /* evaluate this node's jacobian */
     CSC_matrix *Jchild_CSC = lnode->Jchild_CSC;
@@ -207,7 +207,7 @@ static void wsum_hess_init_impl(expr *node)
     node->work->dwork = (double *) sp_malloc(dim * sizeof(double));
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     left_matmul_expr *lnode = (left_matmul_expr *) node;
 
@@ -216,7 +216,7 @@ static void eval_wsum_hess(expr *node, const double *w)
     int n_blocks = lnode->n_blocks;
     AT->block_left_mult_vec(AT, w, node->work->dwork, n_blocks);
 
-    node->left->eval_wsum_hess(node->left, node->work->dwork);
+    eval_wsum_hess(node->left, node->work->dwork);
     memcpy(node->wsum_hess->x, node->left->wsum_hess->x,
            node->wsum_hess->nnz * sizeof(double));
 }
@@ -270,8 +270,8 @@ expr *new_left_matmul(expr *param_node, expr *u, const CSR_matrix *A)
     expr *node = &lnode->base;
     /* Sparse A — always the general CSC-mirror path. */
     init_expr(node, d1, d2, u->n_vars, forward, jacobian_init_sparse,
-              eval_jacobian_sparse, is_affine, wsum_hess_init_impl, eval_wsum_hess,
-              free_type_data);
+              eval_jacobian_sparse, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, free_type_data);
     node->left = u;
     expr_retain(u);
 
@@ -310,8 +310,8 @@ expr *new_left_matmul_dense(expr *param_node, expr *u, int m, int n,
         (left_matmul_expr *) sp_calloc(1, sizeof(left_matmul_expr));
     expr *node = &lnode->base;
     init_expr(node, d1, d2, u->n_vars, forward, jacobian_init_dense,
-              eval_jacobian_dense, is_affine, wsum_hess_init_impl, eval_wsum_hess,
-              free_type_data);
+              eval_jacobian_dense, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, free_type_data);
     node->left = u;
     expr_retain(u);
 
diff --git a/src/atoms/affine/neg.c b/src/atoms/affine/neg.c
index a42c305f..fe4bd396 100644
--- a/src/atoms/affine/neg.c
+++ b/src/atoms/affine/neg.c
@@ -43,10 +43,10 @@ static void jacobian_init_impl(expr *node)
     node->jacobian = x->jacobian->copy_sparsity(x->jacobian);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     /* evaluate child's jacobian */
-    node->left->eval_jacobian(node->left);
+    eval_jacobian(node->left);
 
     /* negate values only (sparsity pattern set in jacobian_init_impl) */
     for (int k = 0; k < node->left->jacobian->nnz; k++)
@@ -66,10 +66,10 @@ static void wsum_hess_init_impl(expr *node)
     node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     /* For f(x) = -g(x): d²f/dx² = -d²g/dx² */
-    node->left->eval_wsum_hess(node->left, w);
+    eval_wsum_hess(node->left, w);
 
     /* negate values (sparsity pattern set in wsum_hess_init_impl) */
     for (int k = 0; k < node->left->wsum_hess->nnz; k++)
@@ -87,7 +87,8 @@ expr *new_neg(expr *child)
 {
     expr *node = (expr *) sp_calloc(1, sizeof(expr));
     init_expr(node, child->d1, child->d2, child->n_vars, forward, jacobian_init_impl,
-              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL);
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, NULL);
     node->left = child;
     expr_retain(child);
 
diff --git a/src/atoms/affine/parameter.c b/src/atoms/affine/parameter.c
index a757edc9..39b010e1 100644
--- a/src/atoms/affine/parameter.c
+++ b/src/atoms/affine/parameter.c
@@ -36,7 +36,7 @@ static void jacobian_init_impl(expr *node)
     node->jacobian = new_sparse_matrix_alloc(node->size, node->n_vars, 0);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     (void) node;
 }
@@ -47,7 +47,7 @@ static void wsum_hess_init_impl(expr *node)
     node->wsum_hess = new_sparse_matrix_alloc(node->n_vars, node->n_vars, 0);
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     (void) node;
     (void) w;
@@ -63,8 +63,8 @@ expr *new_parameter(int d1, int d2, int param_id, int n_vars, const double *valu
 {
     parameter_expr *pnode = (parameter_expr *) sp_calloc(1, sizeof(parameter_expr));
     expr *node = &pnode->base;
-    init_expr(node, d1, d2, n_vars, forward, jacobian_init_impl, eval_jacobian,
-              is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL);
+    init_expr(node, d1, d2, n_vars, forward, jacobian_init_impl, eval_jacobian_impl,
+              is_affine, wsum_hess_init_impl, eval_wsum_hess_impl, NULL);
 
     // TODO we should assert that the values array has the correct size.
     pnode->param_id = param_id;
diff --git a/src/atoms/affine/promote.c b/src/atoms/affine/promote.c
index 529f846e..7fea9140 100644
--- a/src/atoms/affine/promote.c
+++ b/src/atoms/affine/promote.c
@@ -45,9 +45,9 @@ static void jacobian_init_impl(expr *node)
     node->jacobian = x->jacobian->promote_alloc(x->jacobian, node->size);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
-    node->left->eval_jacobian(node->left);
+    eval_jacobian(node->left);
 
     /* tile the child's single row into the preallocated output. */
     node->left->jacobian->promote_fill_values(node->left->jacobian, node->jacobian);
@@ -59,7 +59,7 @@ static void wsum_hess_init_impl(expr *node)
     node->wsum_hess = node->left->wsum_hess->copy_sparsity(node->left->wsum_hess);
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     /* Sum all weights (they all correspond to the same scalar child) */
     double sum_w = 0.0;
@@ -69,7 +69,7 @@ static void eval_wsum_hess(expr *node, const double *w)
     }
 
     /* evaluate child's wsum_hess with summed weight */
-    node->left->eval_wsum_hess(node->left, &sum_w);
+    eval_wsum_hess(node->left, &sum_w);
 
     /* copy values */
     memcpy(node->wsum_hess->x, node->left->wsum_hess->x,
@@ -86,7 +86,8 @@ expr *new_promote(expr *child, int d1, int d2)
     assert(child->size == 1);
     expr *node = (expr *) sp_calloc(1, sizeof(expr));
     init_expr(node, d1, d2, child->n_vars, forward, jacobian_init_impl,
-              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL);
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, NULL);
     node->left = child;
     expr_retain(child);
 
diff --git a/src/atoms/affine/reshape.c b/src/atoms/affine/reshape.c
index 55171125..68625f61 100644
--- a/src/atoms/affine/reshape.c
+++ b/src/atoms/affine/reshape.c
@@ -39,10 +39,10 @@ static void jacobian_init_impl(expr *node)
     node->jacobian = x->jacobian->copy_sparsity(x->jacobian);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     expr *x = node->left;
-    x->eval_jacobian(x);
+    eval_jacobian(x);
     memcpy(node->jacobian->x, x->jacobian->x, x->jacobian->nnz * sizeof(double));
 }
 
@@ -53,10 +53,10 @@ static void wsum_hess_init_impl(expr *node)
     node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     expr *x = node->left;
-    x->eval_wsum_hess(x, w);
+    eval_wsum_hess(x, w);
     memcpy(node->wsum_hess->x, x->wsum_hess->x,
            node->wsum_hess->nnz * sizeof(double));
 }
@@ -71,7 +71,8 @@ expr *new_reshape(expr *child, int d1, int d2)
     assert(d1 * d2 == child->size);
     expr *node = (expr *) sp_calloc(1, sizeof(expr));
     init_expr(node, d1, d2, child->n_vars, forward, jacobian_init_impl,
-              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL);
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, NULL);
     node->left = child;
     expr_retain(child);
 
diff --git a/src/atoms/affine/scalar_mult.c b/src/atoms/affine/scalar_mult.c
index adb0c6a2..e2886bb5 100644
--- a/src/atoms/affine/scalar_mult.c
+++ b/src/atoms/affine/scalar_mult.c
@@ -65,13 +65,13 @@ static void jacobian_init_impl(expr *node)
     node->jacobian = x->jacobian->copy_sparsity(x->jacobian);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     expr *child = node->left;
     double a = ((scalar_mult_expr *) node)->param_source->value[0];
 
     /* evaluate child */
-    child->eval_jacobian(child);
+    eval_jacobian(child);
 
     /* scale child's jacobian */
     for (int j = 0; j < child->jacobian->nnz; j++)
@@ -91,10 +91,10 @@ static void wsum_hess_init_impl(expr *node)
     node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     expr *x = node->left;
-    x->eval_wsum_hess(x, w);
+    eval_wsum_hess(x, w);
 
     double a = ((scalar_mult_expr *) node)->param_source->value[0];
     for (int j = 0; j < x->wsum_hess->nnz; j++)
@@ -125,8 +125,8 @@ expr *new_scalar_mult(expr *param_node, expr *child)
     expr *node = &mult_node->base;
 
     init_expr(node, child->d1, child->d2, child->n_vars, forward, jacobian_init_impl,
-              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess,
-              free_type_data);
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, free_type_data);
     node->left = child;
     expr_retain(child);
 
diff --git a/src/atoms/affine/sum.c b/src/atoms/affine/sum.c
index edbb2571..42aa2ecd 100644
--- a/src/atoms/affine/sum.c
+++ b/src/atoms/affine/sum.c
@@ -94,12 +94,12 @@ static void jacobian_init_impl(expr *node)
                                                           x->d1, snode->idx_map);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     expr *child = node->left;
 
     /* evaluate child's jacobian */
-    child->eval_jacobian(child);
+    eval_jacobian(child);
 
     /* we have precomputed an idx map between the nonzeros of the child's jacobian
        and this node's jacobian, so we just accumulate accordingly */
@@ -119,7 +119,7 @@ static void wsum_hess_init_impl(expr *node)
     node->work->dwork = sp_malloc(child->size * sizeof(double));
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     expr *child = node->left;
     sum_expr *snode = (sum_expr *) node;
@@ -138,7 +138,7 @@ static void eval_wsum_hess(expr *node, const double *w)
         tile_double(node->work->dwork, w, child->d1, child->d2);
     }
 
-    child->eval_wsum_hess(child, node->work->dwork);
+    eval_wsum_hess(child, node->work->dwork);
 
     memcpy(node->wsum_hess->x, child->wsum_hess->x,
            node->wsum_hess->nnz * sizeof(double));
@@ -181,8 +181,9 @@ expr *new_sum(expr *child, int axis)
 
     /* to be consistent with CVXPY and NumPy we treat the result from
        sum with an axis argument as a row vector */
-    init_expr(node, 1, d2, child->n_vars, forward, jacobian_init_impl, eval_jacobian,
-              is_affine, wsum_hess_init_impl, eval_wsum_hess, free_type_data);
+    init_expr(node, 1, d2, child->n_vars, forward, jacobian_init_impl,
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, free_type_data);
     node->left = child;
     expr_retain(child);
 
diff --git a/src/atoms/affine/trace.c b/src/atoms/affine/trace.c
index d6477a50..80e16650 100644
--- a/src/atoms/affine/trace.c
+++ b/src/atoms/affine/trace.c
@@ -86,13 +86,13 @@ static void jacobian_init_impl(expr *node)
     node->jacobian = new_sparse_matrix(jac);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     expr *x = node->left;
     trace_expr *tnode = (trace_expr *) node;
 
     /* evaluate child's jacobian */
-    x->eval_jacobian(x);
+    eval_jacobian(x);
 
     /* local jacobian */
     memset(node->jacobian->x, 0, node->jacobian->nnz * sizeof(double));
@@ -117,7 +117,7 @@ static void wsum_hess_init_impl(expr *node)
     node->wsum_hess = x->wsum_hess->copy_sparsity(x->wsum_hess);
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     expr *x = node->left;
 
@@ -127,7 +127,7 @@ static void eval_wsum_hess(expr *node, const double *w)
         node->work->dwork[i] = w[0];
     }
 
-    x->eval_wsum_hess(x, node->work->dwork);
+    eval_wsum_hess(x, node->work->dwork);
 
     memcpy(node->wsum_hess->x, x->wsum_hess->x,
            node->wsum_hess->nnz * sizeof(double));
@@ -151,8 +151,9 @@ expr *new_trace(expr *child)
 {
     trace_expr *tnode = (trace_expr *) sp_calloc(1, sizeof(trace_expr));
     expr *node = &tnode->base;
-    init_expr(node, 1, 1, child->n_vars, forward, jacobian_init_impl, eval_jacobian,
-              is_affine, wsum_hess_init_impl, eval_wsum_hess, free_type_data);
+    init_expr(node, 1, 1, child->n_vars, forward, jacobian_init_impl,
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, free_type_data);
     node->left = child;
     expr_retain(child);
 
diff --git a/src/atoms/affine/transpose.c b/src/atoms/affine/transpose.c
index 47603761..40bf61de 100644
--- a/src/atoms/affine/transpose.c
+++ b/src/atoms/affine/transpose.c
@@ -56,10 +56,10 @@ static void jacobian_init_impl(expr *node)
     node->work->iwork = indices;
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     expr *child = node->left;
-    child->eval_jacobian(child);
+    eval_jacobian(child);
     child->jacobian->index_fill_values(child->jacobian, node->work->iwork,
                                        node->size, node->jacobian);
 }
@@ -76,7 +76,7 @@ static void wsum_hess_init_impl(expr *node)
     /* for computing Kw where K is the commutation matrix */
     node->work->dwork = (double *) sp_malloc(node->size * sizeof(double));
 }
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     int d2 = node->d2;
     int d1 = node->d1;
@@ -90,7 +90,7 @@ static void eval_wsum_hess(expr *node, const double *w)
         }
     }
 
-    node->left->eval_wsum_hess(node->left, node->work->dwork);
+    eval_wsum_hess(node->left, node->work->dwork);
 
     /* copy to this node's hessian */
     memcpy(node->wsum_hess->x, node->left->wsum_hess->x,
@@ -106,7 +106,8 @@ expr *new_transpose(expr *child)
 {
     expr *node = (expr *) sp_calloc(1, sizeof(expr));
     init_expr(node, child->d2, child->d1, child->n_vars, forward, jacobian_init_impl,
-              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL);
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, NULL);
     node->left = child;
     expr_retain(child);
 
diff --git a/src/atoms/affine/variable.c b/src/atoms/affine/variable.c
index 3291829f..f3cc7d09 100644
--- a/src/atoms/affine/variable.c
+++ b/src/atoms/affine/variable.c
@@ -39,7 +39,7 @@ static void jacobian_init_impl(expr *node)
     node->jacobian = new_sparse_matrix(jac);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     /* Variable jacobian never changes - nothing to evaluate */
     (void) node;
@@ -67,7 +67,7 @@ static bool is_affine(const expr *node)
 expr *new_variable(int d1, int d2, int var_id, int n_vars)
 {
     expr *node = (expr *) sp_calloc(1, sizeof(expr));
-    init_expr(node, d1, d2, n_vars, forward, jacobian_init_impl, eval_jacobian,
+    init_expr(node, d1, d2, n_vars, forward, jacobian_init_impl, eval_jacobian_impl,
               is_affine, wsum_hess_init_impl, wsum_hess_eval, NULL);
     node->var_id = var_id;
 
diff --git a/src/atoms/affine/vector_mult.c b/src/atoms/affine/vector_mult.c
index e5cc5f3e..c4635748 100644
--- a/src/atoms/affine/vector_mult.c
+++ b/src/atoms/affine/vector_mult.c
@@ -65,13 +65,13 @@ static void jacobian_init_impl(expr *node)
     node->jacobian = x->jacobian->copy_sparsity(x->jacobian);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     expr *x = node->left;
     const double *a = ((vector_mult_expr *) node)->param_source->value;
 
     /* evaluate jacobian of child */
-    x->eval_jacobian(x);
+    eval_jacobian(x);
 
     /* row-wise scale child's jacobian: diag(a) @ Jx */
     x->jacobian->DA_fill_values(a, x->jacobian, node->jacobian);
@@ -91,7 +91,7 @@ static void wsum_hess_init_impl(expr *node)
     node->work->dwork = (double *) sp_malloc(node->size * sizeof(double));
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     expr *x = node->left;
     const double *a = ((vector_mult_expr *) node)->param_source->value;
@@ -102,7 +102,7 @@ static void eval_wsum_hess(expr *node, const double *w)
         node->work->dwork[i] = a[i] * w[i];
     }
 
-    x->eval_wsum_hess(x, node->work->dwork);
+    eval_wsum_hess(x, node->work->dwork);
 
     /* copy values from child to this node */
     memcpy(node->wsum_hess->x, x->wsum_hess->x,
@@ -131,8 +131,8 @@ expr *new_vector_mult(expr *param_node, expr *child)
     expr *node = &vnode->base;
 
     init_expr(node, child->d1, child->d2, child->n_vars, forward, jacobian_init_impl,
-              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess,
-              free_type_data);
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, free_type_data);
     node->left = child;
     expr_retain(child);
 
diff --git a/src/atoms/bivariate_full_dom/matmul.c b/src/atoms/bivariate_full_dom/matmul.c
index c85fde8f..39202f19 100644
--- a/src/atoms/bivariate_full_dom/matmul.c
+++ b/src/atoms/bivariate_full_dom/matmul.c
@@ -257,13 +257,11 @@ static void eval_jacobian_chain_rule(expr *node)
     int k = f->d2;
     int n = g->d2;
 
-    /* evaluate Jacobians of children */
-    f->eval_jacobian(f);
-    g->eval_jacobian(g);
-    csr_to_csc_fill_values(f->jacobian->to_csr(f->jacobian), f->work->jacobian_csc,
-                           f->work->csc_work);
-    csr_to_csc_fill_values(g->jacobian->to_csr(g->jacobian), g->work->jacobian_csc,
-                           g->work->csc_work);
+    /* evaluate Jacobians of children and refresh their CSC mirrors */
+    eval_jacobian(f);
+    eval_jacobian(g);
+    expr_refresh_jacobian_csc(f);
+    expr_refresh_jacobian_csc(g);
 
     /* evaluate term1, term2, and their sum */
     YT_kron_I_fill_values(m, k, n, g->value, f->work->jacobian_csc,
@@ -494,27 +492,11 @@ static void eval_wsum_hess_chain_rule(expr *node, const double *w)
     CSC_matrix *Jf = f->work->jacobian_csc;
     CSC_matrix *Jg = g->work->jacobian_csc;
 
-    /* refresh child Jacobian CSC_matrix values (cache if affine) */
-    if (!f->work->jacobian_csc_filled)
-    {
-        csr_to_csc_fill_values(f->jacobian->to_csr(f->jacobian), Jf,
-                               f->work->csc_work);
-        if (is_f_affine)
-        {
-            f->work->jacobian_csc_filled = true;
-        }
-    }
-
-    /* refresh child Jacobian CSC_matrix values (cache if affine) */
-    if (!g->work->jacobian_csc_filled)
-    {
-        csr_to_csc_fill_values(g->jacobian->to_csr(g->jacobian), Jg,
-                               g->work->csc_work);
-        if (is_g_affine)
-        {
-            g->work->jacobian_csc_filled = true;
-        }
-    }
+    /* refresh child Jacobian CSC mirrors (version-guarded: no-ops when the
+       jacobian pass already refreshed them this eval, or when an affine
+       child's values are unchanged) */
+    expr_refresh_jacobian_csc(f);
+    expr_refresh_jacobian_csc(g);
 
     /* compute C = J_f^T @ B(w) @ J_g */
     fill_cross_hessian_values(m, k, n, w, mnode->B);
@@ -529,14 +511,14 @@ static void eval_wsum_hess_chain_rule(expr *node, const double *w)
     if (!is_f_affine)
     {
         Y_kron_I_vec(m, k, n, g->value, w, node->work->dwork);
-        f->eval_wsum_hess(f, node->work->dwork);
+        eval_wsum_hess(f, node->work->dwork);
     }
 
     /* compute Hessian of g */
     if (!is_g_affine)
     {
         I_kron_XT_vec(m, k, n, f->value, w, node->work->dwork);
-        g->eval_wsum_hess(g, node->work->dwork);
+        eval_wsum_hess(g, node->work->dwork);
     }
 
     /* accumulate H = C + C^T + H_f + H_g */
diff --git a/src/atoms/bivariate_full_dom/multiply.c b/src/atoms/bivariate_full_dom/multiply.c
index 1b183d85..dbfc115d 100644
--- a/src/atoms/bivariate_full_dom/multiply.c
+++ b/src/atoms/bivariate_full_dom/multiply.c
@@ -61,13 +61,13 @@ static void jacobian_init_impl(expr *node)
     sum_matrices_alloc(node->left->jacobian, node->right->jacobian, node->jacobian);
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     expr *x = node->left;
     expr *y = node->right;
 
-    x->eval_jacobian(x);
-    y->eval_jacobian(y);
+    eval_jacobian(x);
+    eval_jacobian(y);
 
     /* chain rule: the jacobian of h(x) = f(g1(x), g2(x))) is Jh = J_{f, 1} J_{g1} +
      * J_{f, 2} J_{g2} */
@@ -216,7 +216,7 @@ static void wsum_hess_init_impl(expr *node)
     }
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     expr *x = node->left;
     expr *y = node->right;
@@ -235,31 +235,18 @@ static void eval_wsum_hess(expr *node, const double *w)
         // ----------------------------------------------------------------------
         //  Refresh each operand's CSC_matrix cache as needed for the (Sparse,
         //  Sparse) dispatch path. For PD operands, refresh_csc_values is a no-op.
-        //  The jacobian_csc_filled flag preserves the affine optimization: we only
-        //  refresh on the first eval for affine children.
+        //  The refresh is version-guarded, so an affine child's cache is only
+        //  filled once per parameter epoch.
         // ----------------------------------------------------------------------
-        if (!x->work->jacobian_csc_filled)
-        {
-            x->jacobian->refresh_csc_values(x->jacobian);
-            if (is_x_affine)
-            {
-                x->work->jacobian_csc_filled = true;
-            }
-        }
-        if (!y->work->jacobian_csc_filled)
-        {
-            y->jacobian->refresh_csc_values(y->jacobian);
-            if (is_y_affine)
-            {
-                y->work->jacobian_csc_filled = true;
-            }
-        }
+        x->jacobian->refresh_csc_values(x->jacobian);
+        y->jacobian->refresh_csc_values(y->jacobian);
 
         // ---------------------------------------------------------------
         //                    compute C and CT
         // ---------------------------------------------------------------
         elementwise_mult_expr *mul_node = (elementwise_mult_expr *) node;
         BTDA_matrices_fill_values(x->jacobian, w, y->jacobian, mul_node->C);
+        matrix_values_changed(mul_node->C);
         mul_node->C->transpose_fill_values(mul_node->C, mul_node->CT);
 
         // ---------------------------------------------------------------
@@ -271,7 +258,7 @@ static void eval_wsum_hess(expr *node, const double *w)
             {
                 node->work->dwork[i] = w[i] * y->value[i];
             }
-            x->eval_wsum_hess(x, node->work->dwork);
+            eval_wsum_hess(x, node->work->dwork);
         }
 
         if (!is_y_affine)
@@ -280,7 +267,7 @@ static void eval_wsum_hess(expr *node, const double *w)
             {
                 node->work->dwork[i] = w[i] * x->value[i];
             }
-            y->eval_wsum_hess(y, node->work->dwork);
+            eval_wsum_hess(y, node->work->dwork);
         }
 
         // ---------------------------------------------------------------
@@ -322,8 +309,8 @@ expr *new_elementwise_mult(expr *left, expr *right)
     expr *node = &mul_node->base;
 
     init_expr(node, left->d1, left->d2, left->n_vars, forward, jacobian_init_impl,
-              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess,
-              free_type_data);
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, free_type_data);
     node->left = left;
     node->right = right;
     expr_retain(left);
diff --git a/src/atoms/bivariate_restricted_dom/quad_over_lin.c b/src/atoms/bivariate_restricted_dom/quad_over_lin.c
index b4460f91..36bec791 100644
--- a/src/atoms/bivariate_restricted_dom/quad_over_lin.c
+++ b/src/atoms/bivariate_restricted_dom/quad_over_lin.c
@@ -129,13 +129,15 @@ static void jacobian_init_impl(expr *node)
 
         /* prepare CSC_matrix form of child jacobian for chain rule.
          * For a linear operator the values are constant, so fill
-         * them once here. */
+         * them once here and record the version so the counter stays
+         * truthful. */
         jacobian_csc_init(x);
         csr_to_csc_fill_values(Jx, x->work->jacobian_csc, x->work->csc_work);
+        x->work->jacobian_csc_seen = x->jacobian->values_version;
     }
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     expr *x = node->left;
     expr *y = node->right;
@@ -266,7 +268,7 @@ static void wsum_hess_init_impl(expr *node)
     }
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     double *x = node->left->value;
     double y = node->right->value[0];
@@ -342,8 +344,9 @@ expr *new_quad_over_lin(expr *left, expr *right)
     }
 
     expr *node = (expr *) sp_calloc(1, sizeof(expr));
-    init_expr(node, 1, 1, left->n_vars, forward, jacobian_init_impl, eval_jacobian,
-              is_affine, wsum_hess_init_impl, eval_wsum_hess, NULL);
+    init_expr(node, 1, 1, left->n_vars, forward, jacobian_init_impl,
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, NULL);
     node->left = left;
     node->right = right;
     expr_retain(left);
diff --git a/src/atoms/elementwise_full_dom/common.c b/src/atoms/elementwise_full_dom/common.c
index c834b737..0775382e 100644
--- a/src/atoms/elementwise_full_dom/common.c
+++ b/src/atoms/elementwise_full_dom/common.c
@@ -64,7 +64,7 @@ void eval_jacobian_elementwise(expr *node)
     else
     {
         /* jacobian of h(x) = f(g(x)) is Jf @ Jg, and here Jf is diagonal */
-        child->eval_jacobian(child);
+        eval_jacobian(child);
         node->local_jacobian(node, node->work->local_jac_diag);
         memcpy(node->work->dwork, node->work->local_jac_diag,
                node->size * sizeof(double));
@@ -139,13 +139,9 @@ void eval_wsum_hess_elementwise(expr *node, const double *w)
     {
         if (child->is_affine(child))
         {
-            /* Refresh the child Jacobian's CSC_matrix mirror once; subsequent calls
-               skip since the affine child's values don't change. */
-            if (!child->work->jacobian_csc_filled)
-            {
-                child->jacobian->refresh_csc_values(child->jacobian);
-                child->work->jacobian_csc_filled = true;
-            }
+            /* Version-guarded: fills the affine child's CSC_matrix mirror on the
+               first call, no-ops after. */
+            child->jacobian->refresh_csc_values(child->jacobian);
 
             node->local_wsum_hess(node, node->work->dwork, w);
             child->jacobian->ATDA_fill_values(child->jacobian, node->work->dwork,
@@ -169,9 +165,10 @@ void eval_wsum_hess_elementwise(expr *node, const double *w)
                 node->work->dwork[k] *= w[k];
             }
 
-            child->eval_wsum_hess(child, node->work->dwork);
+            eval_wsum_hess(child, node->work->dwork);
             memcpy(node->work->hess_term2->x, child->wsum_hess->x,
                    child->wsum_hess->nnz * sizeof(double));
+            matrix_values_changed(node->work->hess_term2);
 
             /* wsum_hess = term1 + term2 */
             sum_matrices_fill_values(node->work->hess_term1, node->work->hess_term2,
diff --git a/src/atoms/elementwise_restricted_dom/atanh.c b/src/atoms/elementwise_restricted_dom/atanh.c
index c2c71c2f..48fd6bf1 100644
--- a/src/atoms/elementwise_restricted_dom/atanh.c
+++ b/src/atoms/elementwise_restricted_dom/atanh.c
@@ -52,7 +52,7 @@ expr *new_atanh(expr *child)
 {
     expr *node = new_restricted(child);
     node->forward = atanh_forward;
-    node->eval_jacobian = atanh_eval_jacobian;
-    node->eval_wsum_hess = atanh_eval_wsum_hess;
+    node->eval_jacobian_impl = atanh_eval_jacobian;
+    node->eval_wsum_hess_impl = atanh_eval_wsum_hess;
     return node;
 }
diff --git a/src/atoms/elementwise_restricted_dom/entr.c b/src/atoms/elementwise_restricted_dom/entr.c
index c3d29cbd..b9c65c61 100644
--- a/src/atoms/elementwise_restricted_dom/entr.c
+++ b/src/atoms/elementwise_restricted_dom/entr.c
@@ -53,7 +53,7 @@ expr *new_entr(expr *child)
 {
     expr *node = new_restricted(child);
     node->forward = entr_forward;
-    node->eval_jacobian = entr_eval_jacobian;
-    node->eval_wsum_hess = entr_eval_wsum_hess;
+    node->eval_jacobian_impl = entr_eval_jacobian;
+    node->eval_wsum_hess_impl = entr_eval_wsum_hess;
     return node;
 }
diff --git a/src/atoms/elementwise_restricted_dom/log.c b/src/atoms/elementwise_restricted_dom/log.c
index d9753e13..929b0f47 100644
--- a/src/atoms/elementwise_restricted_dom/log.c
+++ b/src/atoms/elementwise_restricted_dom/log.c
@@ -52,7 +52,7 @@ expr *new_log(expr *child)
 {
     expr *node = new_restricted(child);
     node->forward = log_forward;
-    node->eval_jacobian = log_eval_jacobian;
-    node->eval_wsum_hess = log_eval_wsum_hess;
+    node->eval_jacobian_impl = log_eval_jacobian;
+    node->eval_wsum_hess_impl = log_eval_wsum_hess;
     return node;
 }
diff --git a/src/atoms/elementwise_restricted_dom/tan.c b/src/atoms/elementwise_restricted_dom/tan.c
index 3a4e4dc7..0e9fc503 100644
--- a/src/atoms/elementwise_restricted_dom/tan.c
+++ b/src/atoms/elementwise_restricted_dom/tan.c
@@ -53,7 +53,7 @@ expr *new_tan(expr *child)
 {
     expr *node = new_restricted(child);
     node->forward = tan_forward;
-    node->eval_jacobian = tan_eval_jacobian;
-    node->eval_wsum_hess = tan_eval_wsum_hess;
+    node->eval_jacobian_impl = tan_eval_jacobian;
+    node->eval_wsum_hess_impl = tan_eval_wsum_hess;
     return node;
 }
diff --git a/src/atoms/other/prod.c b/src/atoms/other/prod.c
index 3feab4a5..5d1078e2 100644
--- a/src/atoms/other/prod.c
+++ b/src/atoms/other/prod.c
@@ -89,7 +89,7 @@ static void jacobian_init_impl(expr *node)
     }
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     expr *x = node->left;
     prod_expr *pnode = (prod_expr *) node;
@@ -167,7 +167,7 @@ static void wsum_hess_init_impl(expr *node)
     }
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     expr *x = node->left;
     int num_of_zeros = ((prod_expr *) node)->num_of_zeros;
@@ -219,8 +219,9 @@ expr *new_prod(expr *child)
     /* Output is scalar: 1 x 1 */
     prod_expr *pnode = (prod_expr *) sp_calloc(1, sizeof(prod_expr));
     expr *node = &pnode->base;
-    init_expr(node, 1, 1, child->n_vars, forward, jacobian_init_impl, eval_jacobian,
-              is_affine, wsum_hess_init_impl, eval_wsum_hess, free_type_data);
+    init_expr(node, 1, 1, child->n_vars, forward, jacobian_init_impl,
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, free_type_data);
     node->left = child;
     expr_retain(child);
     return node;
diff --git a/src/atoms/other/prod_axis_one.c b/src/atoms/other/prod_axis_one.c
index 591657b1..6e74cc30 100644
--- a/src/atoms/other/prod_axis_one.c
+++ b/src/atoms/other/prod_axis_one.c
@@ -110,7 +110,7 @@ static void jacobian_init_impl(expr *node)
     }
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     expr *x = node->left;
     prod_axis *pnode = (prod_axis *) node;
@@ -344,7 +344,7 @@ static inline void wsum_hess_row_many_zeros(expr *node, int row, int d2)
     }
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     expr *x = node->left;
     prod_axis *pnode = (prod_axis *) node;
@@ -403,8 +403,8 @@ expr *new_prod_axis_one(expr *child)
 
     /* output is always a row vector 1 x d1 (one product per row) */
     init_expr(node, 1, child->d1, child->n_vars, forward, jacobian_init_impl,
-              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess,
-              free_type_data);
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, free_type_data);
 
     /* allocate arrays to store per-row statistics */
     pnode->num_of_zeros = (int *) sp_calloc(child->d1, sizeof(int));
diff --git a/src/atoms/other/prod_axis_zero.c b/src/atoms/other/prod_axis_zero.c
index b3987ac7..888e8f15 100644
--- a/src/atoms/other/prod_axis_zero.c
+++ b/src/atoms/other/prod_axis_zero.c
@@ -105,7 +105,7 @@ static void jacobian_init_impl(expr *node)
     }
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     expr *x = node->left;
     prod_axis *pnode = (prod_axis *) node;
@@ -299,7 +299,7 @@ static inline void wsum_hess_column_many_zeros(expr *node, const double *w, int
     (void) w;
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     expr *x = node->left;
     prod_axis *pnode = (prod_axis *) node;
@@ -363,8 +363,8 @@ expr *new_prod_axis_zero(expr *child)
 
     /* output is always a row vector 1 x d2 - TODO: is that correct? */
     init_expr(node, 1, child->d2, child->n_vars, forward, jacobian_init_impl,
-              eval_jacobian, is_affine, wsum_hess_init_impl, eval_wsum_hess,
-              free_type_data);
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, free_type_data);
 
     /* allocate arrays to store per-column statistics */
     pnode->num_of_zeros = (int *) sp_calloc(child->d2, sizeof(int));
diff --git a/src/atoms/other/quad_form.c b/src/atoms/other/quad_form.c
index 07141aaa..5ba79ae6 100644
--- a/src/atoms/other/quad_form.c
+++ b/src/atoms/other/quad_form.c
@@ -113,7 +113,7 @@ static void jacobian_init_impl(expr *node)
     }
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     quad_form_expr *qnode = (quad_form_expr *) node;
     expr *x = node->left;
@@ -129,18 +129,8 @@ static void eval_jacobian(expr *node)
     else
     {
         /* jacobian = 2 * (Q @ f(x))^T @ J_f */
-        x->eval_jacobian(x);
-
-        if (!x->work->jacobian_csc_filled)
-        {
-            csr_to_csc_fill_values(x->jacobian->to_csr(x->jacobian),
-                                   x->work->jacobian_csc, x->work->csc_work);
-
-            if (x->is_affine(x))
-            {
-                x->work->jacobian_csc_filled = true;
-            }
-        }
+        eval_jacobian(x);
+        expr_refresh_jacobian_csc(x);
 
         /* The jacobian has same values as the gradient, which is
            J_f^T (Q @ f(x)). Here, dwork stores Q @ f(x) from forward */
@@ -225,18 +215,11 @@ static void eval_wsum_hess_sparse(expr *node, const double *w)
     }
     else
     {
-        /* fill the CSC_matrix representation of the Jacobian of the child */
+        /* refresh the CSC_matrix mirror of the child's Jacobian (version-guarded:
+           no-ops when the jacobian pass already refreshed it this eval, or when
+           an affine child's values are unchanged) */
         CSC_matrix *Jf = x->work->jacobian_csc;
-        if (!x->work->jacobian_csc_filled)
-        {
-            csr_to_csc_fill_values(x->jacobian->to_csr(x->jacobian), Jf,
-                                   x->work->csc_work);
-
-            if (x->is_affine(x))
-            {
-                x->work->jacobian_csc_filled = true;
-            }
-        }
+        expr_refresh_jacobian_csc(x);
 
         CSC_matrix *QJf = qnode->QJf;
         CSR_matrix *term1 = node->work->hess_term1->to_csr(node->work->hess_term1);
@@ -246,7 +229,7 @@ static void eval_wsum_hess_sparse(expr *node, const double *w)
         BTDA_fill_values(Jf, QJf, NULL, term1);
 
         /* term2 */
-        x->eval_wsum_hess(x, node->work->dwork);
+        eval_wsum_hess(x, node->work->dwork);
         memcpy(node->work->hess_term2->x, x->wsum_hess->x,
                x->wsum_hess->nnz * sizeof(double));
 
@@ -255,6 +238,7 @@ static void eval_wsum_hess_sparse(expr *node, const double *w)
                     1);
         cblas_dscal(node->work->hess_term2->nnz, two_w, node->work->hess_term2->x,
                     1);
+        matrix_values_changed(node->work->hess_term2);
 
         /* sum the two terms */
         sum_matrices_fill_values(node->work->hess_term1, node->work->hess_term2,
@@ -347,11 +331,12 @@ static void eval_wsum_hess_dense(expr *node, const double *w)
                                   node->work->hess_term1);
 
         /* term2 = 2w sum_i (Q f(x))_i nabla^2 f_i (dwork = Q f(x) from forward) */
-        x->eval_wsum_hess(x, node->work->dwork);
+        eval_wsum_hess(x, node->work->dwork);
         memcpy(node->work->hess_term2->x, x->wsum_hess->x,
                x->wsum_hess->nnz * sizeof(double));
         cblas_dscal(node->work->hess_term2->nnz, two_w, node->work->hess_term2->x,
                     1);
+        matrix_values_changed(node->work->hess_term2);
 
         sum_matrices_fill_values(node->work->hess_term1, node->work->hess_term2,
                                  node->wsum_hess);
@@ -395,9 +380,9 @@ expr *new_quad_form_sparse(expr *left, CSR_matrix *Q)
     quad_form_expr *qnode = (quad_form_expr *) sp_calloc(1, sizeof(quad_form_expr));
     expr *node = &qnode->base;
 
-    init_expr(node, 1, 1, left->n_vars, forward, jacobian_init_impl, eval_jacobian,
-              is_affine, wsum_hess_init_sparse, eval_wsum_hess_sparse,
-              free_type_data);
+    init_expr(node, 1, 1, left->n_vars, forward, jacobian_init_impl,
+              eval_jacobian_impl, is_affine, wsum_hess_init_sparse,
+              eval_wsum_hess_sparse, free_type_data);
     node->left = left;
     expr_retain(left);
 
@@ -419,8 +404,9 @@ expr *new_quad_form_dense(expr *child, int n, const double *P_data,
     quad_form_expr *qnode = (quad_form_expr *) sp_calloc(1, sizeof(quad_form_expr));
     expr *node = &qnode->base;
 
-    init_expr(node, 1, 1, child->n_vars, forward, jacobian_init_impl, eval_jacobian,
-              is_affine, wsum_hess_init_dense, eval_wsum_hess_dense, free_type_data);
+    init_expr(node, 1, 1, child->n_vars, forward, jacobian_init_impl,
+              eval_jacobian_impl, is_affine, wsum_hess_init_dense,
+              eval_wsum_hess_dense, free_type_data);
     node->left = child;
     expr_retain(child);
 
diff --git a/src/expr.c b/src/expr.c
index 30d5b08a..8e0ef2c1 100644
--- a/src/expr.c
+++ b/src/expr.c
@@ -36,10 +36,10 @@ void init_expr(expr *node, int d1, int d2, int n_vars, forward_fn forward,
     node->var_id = NOT_A_VARIABLE;
     node->forward = forward;
     node->jacobian_init_impl = jacobian_init;
-    node->eval_jacobian = eval_jacobian;
+    node->eval_jacobian_impl = eval_jacobian;
     node->is_affine = is_affine;
     node->wsum_hess_init_impl = wsum_hess_init;
-    node->eval_wsum_hess = eval_wsum_hess;
+    node->eval_wsum_hess_impl = eval_wsum_hess;
     node->free_type_data = free_type_data;
     node->work = (Expr_Work *) sp_calloc(1, sizeof(Expr_Work));
 }
@@ -53,6 +53,19 @@ void jacobian_csc_init(expr *node)
     node->work->csc_work = (int *) sp_malloc(node->n_vars * sizeof(int));
     node->work->jacobian_csc = csr_to_csc_alloc(
         node->jacobian->to_csr(node->jacobian), node->work->csc_work);
+
+    /* Deliberately stale: the mirror is built during the symbolic phase,
+       before values are meaningful, so the first refresh must fill. */
+    node->work->jacobian_csc_seen = node->jacobian->values_version - 1;
+}
+
+/* Refresh work->jacobian_csc from node->jacobian iff its values changed. */
+void expr_refresh_jacobian_csc(expr *node)
+{
+    if (node->work->jacobian_csc_seen == node->jacobian->values_version) return;
+    csr_to_csc_fill_values(node->jacobian->to_csr(node->jacobian),
+                           node->work->jacobian_csc, node->work->csc_work);
+    node->work->jacobian_csc_seen = node->jacobian->values_version;
 }
 
 void free_expr(expr *node)
@@ -98,7 +111,14 @@ void free_expr(expr *node)
 
 void jacobian_init(expr *node)
 {
-    if (node == NULL || node->jacobian != NULL) return;
+    if (node == NULL) return;
+
+    /* Affinity is structural; cache it so the eval_jacobian wrapper doesn't
+       pay the O(subtree) recursive is_affine per eval. Set before the init
+       guard so nodes that preset their jacobian at construction get it too. */
+    node->work->is_affine_cached = node->is_affine(node);
+
+    if (node->jacobian != NULL) return;
     node->jacobian_init_impl(node);
 }
 
@@ -108,13 +128,44 @@ void wsum_hess_init(expr *node)
     node->wsum_hess_init_impl(node);
 }
 
+/* Runs the atom's eval and bumps the Jacobian's values_version — except
+   for an affine node already evaluated this parameter epoch: its values
+   are constant, so consumers correctly see "unchanged". (The impl still
+   runs; skipping it entirely is a planned follow-up.) */
+void eval_jacobian(expr *node)
+{
+    if (node == NULL) return;
+    node->eval_jacobian_impl(node);
+    if (!(node->work->is_affine_cached && node->work->jacobian_evaluated))
+    {
+        matrix_values_changed(node->jacobian);
+    }
+    node->work->jacobian_evaluated = true;
+}
+
+/* Runs the atom's eval and bumps the Hessian's values_version (no affine
+   handling — the weights change every call). */
+void eval_wsum_hess(expr *node, const double *w)
+{
+    if (node == NULL) return;
+    node->eval_wsum_hess_impl(node, w);
+    matrix_values_changed(node->wsum_hess);
+}
+
 void expr_set_needs_refresh(expr *node)
 {
     if (node == NULL) return;
     node->needs_parameter_refresh = true;
-    node->work->jacobian_csc_filled = false;
+
+    /* Re-arm the eval_jacobian wrapper's values_version bump: the next eval
+       after a parameter update may change even an affine node's values. */
+    node->work->jacobian_evaluated = false;
     expr_set_needs_refresh(node->left);
     expr_set_needs_refresh(node->right);
+    if (node->set_needs_refresh_children != NULL)
+    {
+        node->set_needs_refresh_children(node);
+    }
 }
 
 void expr_retain(expr *node)
diff --git a/src/old-code/linear_op.c b/src/old-code/linear_op.c
index 739f38dc..58f5dd0d 100644
--- a/src/old-code/linear_op.c
+++ b/src/old-code/linear_op.c
@@ -66,7 +66,7 @@ static void jacobian_init_impl(expr *node)
     (void) node;
 }
 
-static void eval_jacobian(expr *node)
+static void eval_jacobian_impl(expr *node)
 {
     /* Linear operator jacobian never changes - nothing to evaluate */
     (void) node;
@@ -79,7 +79,7 @@ static void wsum_hess_init_impl(expr *node)
         new_sparse_matrix(new_CSR_matrix(node->n_vars, node->n_vars, 0));
 }
 
-static void eval_wsum_hess(expr *node, const double *w)
+static void eval_wsum_hess_impl(expr *node, const double *w)
 {
     /* Linear operator Hessian is always zero - nothing to evaluate */
     (void) node;
@@ -93,8 +93,9 @@ expr *new_linear(expr *u, const CSR_matrix *A, const double *b)
     linear_op_expr *lin_node =
         (linear_op_expr *) sp_calloc(1, sizeof(linear_op_expr));
     expr *node = &lin_node->base;
-    init_expr(node, A->m, 1, u->n_vars, forward, jacobian_init_impl, eval_jacobian,
-              is_affine, wsum_hess_init_impl, eval_wsum_hess, free_type_data);
+    init_expr(node, A->m, 1, u->n_vars, forward, jacobian_init_impl,
+              eval_jacobian_impl, is_affine, wsum_hess_init_impl,
+              eval_wsum_hess_impl, free_type_data);
     node->left = u;
     expr_retain(u);
 
diff --git a/src/problem.c b/src/problem.c
index fe89ff00..11b49f00 100644
--- a/src/problem.c
+++ b/src/problem.c
@@ -495,7 +495,7 @@ void problem_gradient(problem *prob)
     clock_gettime(CLOCK_MONOTONIC, &timer.start);
 
     /* evaluate jacobian of objective */
-    prob->objective->eval_jacobian(prob->objective);
+    eval_jacobian(prob->objective);
 
     /* copy sparse jacobian to dense gradient */
     memset(prob->gradient_values, 0, prob->n_vars * sizeof(double));
@@ -528,7 +528,7 @@ void problem_jacobian(problem *prob)
             continue;
         }
 
-        c->eval_jacobian(c);
+        eval_jacobian(c);
         memcpy(J->x + nnz_offset, c->jacobian->x, c->jacobian->nnz * sizeof(double));
         nnz_offset += c->jacobian->nnz;
     }
@@ -548,7 +548,7 @@ void problem_hessian(problem *prob, double obj_w, const double *w)
     //             evaluate hessian of objective and constraints
     // ------------------------------------------------------------------------
     expr *obj = prob->objective;
-    obj->eval_wsum_hess(obj, &obj_w);
+    eval_wsum_hess(obj, &obj_w);
 
     int offset = 0;
     expr **constrs = prob->constraints;
@@ -560,7 +560,7 @@ void problem_hessian(problem *prob, double obj_w, const double *w)
             offset += constrs[i]->size;
             continue;
         }
-        constrs[i]->eval_wsum_hess(constrs[i], w + offset);
+        eval_wsum_hess(constrs[i], w + offset);
         offset += constrs[i]->size;
     }
 
diff --git a/src/utils/sparse_matrix.c b/src/utils/sparse_matrix.c
index 31f191cc..8a7162f2 100644
--- a/src/utils/sparse_matrix.c
+++ b/src/utils/sparse_matrix.c
@@ -62,14 +62,19 @@ static void sparse_free(matrix *self)
 /* Forward decl: ctor is referenced by copy_sparsity below. */
 matrix *new_sparse_matrix(CSR_matrix *A);
 
-/* Build the CSC_matrix cache structure if absent. Values are NOT filled here; caller
-   must call refresh_csc_values before consuming. ATA_alloc only needs structure,
-   so it's safe to call without a subsequent refresh. */
+/* Build the CSC_matrix cache structure if absent. Values are NOT filled here; call
+   refresh_csc_values before consuming (it no-ops when the cache is already
+   fresh). ATA_alloc only needs structure, so it's safe to call without a
+   subsequent refresh. */
 void sparse_matrix_ensure_csc_cache(sparse_matrix *sm)
 {
     if (sm->csc_cache != NULL) return;
     sm->csc_iwork = (int *) sp_malloc(sm->csr->n * sizeof(int));
     sm->csc_cache = csr_to_csc_alloc(sm->csr, sm->csc_iwork);
+
+    /* Deliberately stale: the structure is built during the symbolic phase,
+       before values are meaningful, so the first refresh must fill. */
+    sm->csc_seen = sm->base.values_version - 1;
 }
 
 static matrix *sparse_copy_sparsity(const matrix *self)
@@ -92,7 +97,8 @@ static matrix *sparse_ATA_alloc(matrix *self)
     return new_sparse_matrix(ATA_alloc(sm->csc_cache));
 }
 
-/* Caller must have called refresh_csc_values since the last change to csr->x. */
+/* Call refresh_csc_values before consuming; it no-ops when the cache is
+   already fresh. */
 static void sparse_ATDA_fill_values(const matrix *self, const double *d, matrix *out)
 {
     const sparse_matrix *sm = (const sparse_matrix *) self;
@@ -321,13 +327,15 @@ static void sparse_diag_vec_fill_values(matrix *self, matrix *out)
     }
 }
 
-/* Build CSC_matrix structure on first call; refill values from csr->x on every call.
- */
+/* Build CSC_matrix structure on first call; refill values from csr->x iff they
+   changed since the last refresh (tracked via base.values_version). */
 static void sparse_refresh_csc_values(matrix *self)
 {
     sparse_matrix *sm = (sparse_matrix *) self;
     sparse_matrix_ensure_csc_cache(sm);
+    if (sm->csc_seen == sm->base.values_version) return;
     csr_to_csc_fill_values(sm->csr, sm->csc_cache, sm->csc_iwork);
+    sm->csc_seen = sm->base.values_version;
 }
 
 static matrix *sparse_sum_row_partition_alloc(matrix *self, int axis, int d1,
diff --git a/src/utils/stacked_pd.c b/src/utils/stacked_pd.c
index caea15e9..f29bcc2f 100644
--- a/src/utils/stacked_pd.c
+++ b/src/utils/stacked_pd.c
@@ -69,6 +69,7 @@ static void stacked_pd_vtable_DA_fill_values(const double *d, const matrix *self
                                              matrix *out)
 {
     DA_spd_fill_values(d, (const stacked_pd *) self, (stacked_pd *) out);
+    matrix_values_changed(out);
 }
 
 static matrix *stacked_pd_vtable_ATA_alloc(matrix *self)
@@ -80,6 +81,7 @@ static void stacked_pd_vtable_ATDA_fill_values(const matrix *self, const double
                                                matrix *out)
 {
     ATDA_spd_fill_values((const stacked_pd *) self, d, (stacked_pd *) out);
+    matrix_values_changed(out);
 }
 
 static matrix *stacked_pd_vtable_transpose_alloc(const matrix *self)
@@ -90,6 +92,7 @@ static matrix *stacked_pd_vtable_transpose_alloc(const matrix *self)
 static void stacked_pd_vtable_transpose_fill_values(const matrix *self, matrix *out)
 {
     transpose_spd_fill_values((const stacked_pd *) self, (stacked_pd *) out);
+    matrix_values_changed(out);
 }
 
 static void stacked_pd_vtable_refresh_csc_values(matrix *self)
@@ -134,9 +137,17 @@ static CSR_matrix *stacked_pd_to_csr(matrix *self)
     if (spd->csr_cache == NULL)
     {
         spd->csr_cache = stacked_pd_to_csr_alloc(spd);
+
+        /* Deliberately stale: force the first call to copy the values. */
+        spd->csr_seen = spd->base.values_version - 1;
     }
 
-    /* refresh values every call (block X buffers may have changed). */
+    /* refresh values iff the block X buffers changed since the last call */
+    if (spd->csr_seen == spd->base.values_version)
+    {
+        return spd->csr_cache;
+    }
+    spd->csr_seen = spd->base.values_version;
     int *p_arr = spd->csr_cache->p;
     for (int k = 0; k < spd->n_blocks; k++)
     {
diff --git a/tests/all_tests.c b/tests/all_tests.c
index f5ddbed6..c6f137dd 100644
--- a/tests/all_tests.c
+++ b/tests/all_tests.c
@@ -56,6 +56,7 @@
 #include "jacobian_tests/other/test_prod_axis_one.h"
 #include "jacobian_tests/other/test_prod_axis_zero.h"
 #include "jacobian_tests/other/test_quad_form.h"
+#include "jacobian_tests/test_values_version.h"
 #include "numerical_diff/test_numerical_diff.h"
 #include "old-code/test_old_permuted_dense.h"
 #include "problem/test_param_broadcast.h"
@@ -209,6 +210,12 @@ int main(void)
     mu_run_test(test_quad_over_lin4, tests_run);
     mu_run_test(test_quad_over_lin5, tests_run);
     mu_run_test(test_quad_form, tests_run);
+    mu_run_test(test_values_version_non_affine, tests_run);
+    mu_run_test(test_values_version_affine, tests_run);
+    mu_run_test(test_values_version_csc_mirror_dedup, tests_run);
+    mu_run_test(test_values_version_stacked_pd_to_csr, tests_run);
+    mu_run_test(test_values_version_param_under_hstack, tests_run);
+    mu_run_test(test_values_version_spd_hess_terms, tests_run);
     /* commented out - see test_quad_form.h */
     // mu_run_test(test_quad_form2, tests_run);
     mu_run_test(test_jacobian_prod_no_zero, tests_run);
diff --git a/tests/jacobian_tests/affine/test_broadcast.h b/tests/jacobian_tests/affine/test_broadcast.h
index a1f79edc..32cfd242 100644
--- a/tests/jacobian_tests/affine/test_broadcast.h
+++ b/tests/jacobian_tests/affine/test_broadcast.h
@@ -30,7 +30,7 @@ const char *test_broadcast_row_jacobian(void)
     expr *bcast = new_broadcast(var, 2, 3);
     bcast->forward(bcast, u);
     jacobian_init(bcast);
-    bcast->eval_jacobian(bcast);
+    eval_jacobian(bcast);
 
     /* Each variable affects 2 elements (m times) */
     double expected_x[6] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
@@ -72,7 +72,7 @@ const char *test_broadcast_col_jacobian(void)
     expr *bcast = new_broadcast(var, 3, 2);
     bcast->forward(bcast, u);
     jacobian_init(bcast);
-    bcast->eval_jacobian(bcast);
+    eval_jacobian(bcast);
 
     /* Each variable affects 2 elements (n times) */
     double expected_x[6] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
@@ -110,7 +110,7 @@ const char *test_broadcast_scalar_to_matrix_jacobian(void)
     expr *bcast = new_broadcast(var, 2, 3);
     bcast->forward(bcast, u);
     jacobian_init(bcast);
-    bcast->eval_jacobian(bcast);
+    eval_jacobian(bcast);
 
     /* All 6 elements depend on the single input variable */
     double expected_x[6] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
@@ -139,7 +139,7 @@ const char *test_double_broadcast(void)
 
     sum->forward(sum, x_vals);
     jacobian_init(sum);
-    sum->eval_jacobian(sum);
+    eval_jacobian(sum);
 
     /* TODO: what is this test? */
 
diff --git a/tests/jacobian_tests/affine/test_convolve.h b/tests/jacobian_tests/affine/test_convolve.h
index 0eb4db1f..3d9151f9 100644
--- a/tests/jacobian_tests/affine/test_convolve.h
+++ b/tests/jacobian_tests/affine/test_convolve.h
@@ -27,7 +27,7 @@ const char *test_jacobian_convolve(void)
     double u[3] = {1.0, 2.0, 3.0};
     y->forward(y, u);
     jacobian_init(y);
-    y->eval_jacobian(y);
+    eval_jacobian(y);
 
     mu_assert("Jacobian should have 5 rows", y->jacobian->m == 5);
     mu_assert("Jacobian should have 3 columns", y->jacobian->n == 3);
diff --git a/tests/jacobian_tests/affine/test_diag_mat.h b/tests/jacobian_tests/affine/test_diag_mat.h
index 9b669a07..9c377a18 100644
--- a/tests/jacobian_tests/affine/test_diag_mat.h
+++ b/tests/jacobian_tests/affine/test_diag_mat.h
@@ -17,7 +17,7 @@ const char *test_diag_mat_jacobian_variable(void)
 
     dm->forward(dm, u);
     jacobian_init(dm);
-    dm->eval_jacobian(dm);
+    eval_jacobian(dm);
 
     double expected_x[2] = {1.0, 1.0};
     int expected_p[3] = {0, 1, 2};
@@ -46,7 +46,7 @@ const char *test_diag_mat_jacobian_of_log(void)
 
     dm->forward(dm, u);
     jacobian_init(dm);
-    dm->eval_jacobian(dm);
+    eval_jacobian(dm);
 
     double expected_x[2] = {1.0, 0.25};
     int expected_p[3] = {0, 1, 2};
diff --git a/tests/jacobian_tests/affine/test_hstack.h b/tests/jacobian_tests/affine/test_hstack.h
index 017730ee..bd61df82 100644
--- a/tests/jacobian_tests/affine/test_hstack.h
+++ b/tests/jacobian_tests/affine/test_hstack.h
@@ -32,7 +32,7 @@ const char *test_jacobian_hstack_vectors(void)
 
     stack->forward(stack, u);
     jacobian_init(stack);
-    stack->eval_jacobian(stack);
+    eval_jacobian(stack);
 
     /* Expected jacobian: 9x3 with 9 nonzeros (diagonal blocks) */
     double expected_Ax[9] = {1.0,      0.5,      1.0 / 3.0, /* d(log)/dx */
@@ -74,7 +74,7 @@ const char *test_jacobian_hstack_matrix(void)
 
     stack->forward(stack, u);
     jacobian_init(stack);
-    stack->eval_jacobian(stack);
+    eval_jacobian(stack);
 
     /* Expected jacobian: 18x6 with 18 nonzeros (diagonal) */
     double expected_Ax[18] = {
diff --git a/tests/jacobian_tests/affine/test_index.h b/tests/jacobian_tests/affine/test_index.h
index 1d759140..0b71454d 100644
--- a/tests/jacobian_tests/affine/test_index.h
+++ b/tests/jacobian_tests/affine/test_index.h
@@ -49,7 +49,7 @@ const char *test_index_jacobian_of_variable(void)
     expr *idx = new_index(var, 1, 2, indices, 2);
     idx->forward(idx, u);
     jacobian_init(idx);
-    idx->eval_jacobian(idx);
+    eval_jacobian(idx);
 
     /* Jacobian is 2x3 with pattern: row 0 selects col 0, row 1 selects col 2 */
     double expected_x[2] = {1.0, 1.0};
@@ -74,7 +74,7 @@ const char *test_index_jacobian_of_log(void)
     expr *idx = new_index(log_node, 1, 2, indices, 2);
     idx->forward(idx, u);
     jacobian_init(idx);
-    idx->eval_jacobian(idx);
+    eval_jacobian(idx);
 
     /* d/dx log(x) = diag(1/x), then select rows 0 and 2.
      * Each selected row has exactly one nonzero (the diagonal entry).
@@ -101,7 +101,7 @@ const char *test_index_jacobian_repeated(void)
     expr *idx = new_index(var, 1, 2, indices, 2);
     idx->forward(idx, u);
     jacobian_init(idx);
-    idx->eval_jacobian(idx);
+    eval_jacobian(idx);
 
     /* Both rows select column 0 */
     double expected_x[2] = {1.0, 1.0};
@@ -128,7 +128,7 @@ const char *test_index_jacobian_duplicates_exceed_source_nnz(void)
     expr *idx = new_index(var, 1, 6, indices, 6);
     idx->forward(idx, u);
     jacobian_init(idx);
-    idx->eval_jacobian(idx);
+    eval_jacobian(idx);
 
     double expected_x[6] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
     int expected_p[7] = {0, 1, 2, 3, 4, 5, 6};
@@ -155,7 +155,7 @@ const char *test_sum_of_index(void)
 
     s->forward(s, u);
     jacobian_init(s);
-    s->eval_jacobian(s);
+    eval_jacobian(s);
 
     /* Gradient: [1, 0, 1] in sparse form. Single output row holds both nnz. */
     double expected_x[2] = {1.0, 1.0};
diff --git a/tests/jacobian_tests/affine/test_kron.h b/tests/jacobian_tests/affine/test_kron.h
index 22308616..59a7c936 100644
--- a/tests/jacobian_tests/affine/test_kron.h
+++ b/tests/jacobian_tests/affine/test_kron.h
@@ -22,7 +22,7 @@ const char *test_jacobian_kron_const_left(void)
     double u[4] = {5.0, 7.0, 6.0, 8.0};
     Z->forward(Z, u);
     jacobian_init(Z);
-    Z->eval_jacobian(Z);
+    eval_jacobian(Z);
 
     mu_assert("kron J rows", Z->jacobian->m == 16);
     mu_assert("kron J cols", Z->jacobian->n == 4);
@@ -52,7 +52,7 @@ const char *test_jacobian_kron_const_right(void)
     double u[4] = {5.0, 7.0, 6.0, 8.0};
     Z->forward(Z, u);
     jacobian_init(Z);
-    Z->eval_jacobian(Z);
+    eval_jacobian(Z);
 
     mu_assert("kron J rows", Z->jacobian->m == 16);
     mu_assert("kron J nnz", Z->jacobian->nnz == 16);
@@ -83,7 +83,7 @@ const char *test_jacobian_kron_sparse(void)
     double u[4] = {1.0, 3.0, 2.0, 4.0};
     Z->forward(Z, u);
     jacobian_init(Z);
-    Z->eval_jacobian(Z);
+    eval_jacobian(Z);
 
     mu_assert("kron sparse J rows", Z->jacobian->m == 36);
     mu_assert("kron sparse J pruned to 12 nnz", Z->jacobian->nnz == 12);
diff --git a/tests/jacobian_tests/affine/test_left_matmul.h b/tests/jacobian_tests/affine/test_left_matmul.h
index 18a27c87..3e667cb7 100644
--- a/tests/jacobian_tests/affine/test_left_matmul.h
+++ b/tests/jacobian_tests/affine/test_left_matmul.h
@@ -47,7 +47,7 @@ const char *test_jacobian_left_matmul_log(void)
 
     A_log_x->forward(A_log_x, x_vals);
     jacobian_init(A_log_x);
-    A_log_x->eval_jacobian(A_log_x);
+    eval_jacobian(A_log_x);
 
     /* Expected jacobian values: A @ diag(1/x) */
     double expected_Ax[7] = {
@@ -91,7 +91,7 @@ const char *test_jacobian_left_matmul_log_matrix(void)
 
     A_log_x->forward(A_log_x, x_vals);
     jacobian_init(A_log_x);
-    A_log_x->eval_jacobian(A_log_x);
+    eval_jacobian(A_log_x);
 
     /* Expected Jacobian: block-diagonal repeat of A scaled by diag(1./x) */
     double expected_Ax[14] = {/* first column block (x = [1, 2, 3]) */
@@ -169,7 +169,7 @@ const char *test_jacobian_left_matmul_pd_from_composite_child(void)
     double x_vals[2] = {0.5, -1.5};
     A2_A1_x->forward(A2_A1_x, x_vals);
     jacobian_init(A2_A1_x);
-    A2_A1_x->eval_jacobian(A2_A1_x);
+    eval_jacobian(A2_A1_x);
 
     /* Structural: outer's Jacobian must be PD (produced by the
        jacobian_init_pd path via BA_pd_matrices_alloc). */
@@ -220,7 +220,7 @@ const char *test_jacobian_left_matmul_pd_param(void)
     double x_vals[2] = {0.5, -1.5};
     A_x->forward(A_x, x_vals);
     jacobian_init(A_x);
-    A_x->eval_jacobian(A_x);
+    eval_jacobian(A_x);
 
     /* Structural: Jacobian must be PD. */
     mu_assert("Jacobian should be PD", A_x->jacobian->is_permuted_dense);
@@ -248,7 +248,7 @@ const char *test_jacobian_left_matmul_pd_param(void)
     expr_set_needs_refresh(A_x);
 
     A_x->forward(A_x, x_vals);
-    A_x->eval_jacobian(A_x);
+    eval_jacobian(A_x);
 
     mu_assert("Jacobian still PD after refresh", A_x->jacobian->is_permuted_dense);
     double expected_X2[6] = {7.0, 8.0, 9.0, 10.0, 11.0, 12.0};
diff --git a/tests/jacobian_tests/affine/test_neg.h b/tests/jacobian_tests/affine/test_neg.h
index f4d859cd..7064f94a 100644
--- a/tests/jacobian_tests/affine/test_neg.h
+++ b/tests/jacobian_tests/affine/test_neg.h
@@ -12,7 +12,7 @@ const char *test_neg_jacobian(void)
     expr *neg_node = new_neg(var);
     neg_node->forward(neg_node, u);
     jacobian_init(neg_node);
-    neg_node->eval_jacobian(neg_node);
+    eval_jacobian(neg_node);
 
     /* Jacobian of neg(x) is -I (diagonal with -1) */
     double expected_x[3] = {-1.0, -1.0, -1.0};
@@ -40,7 +40,7 @@ const char *test_neg_chain(void)
     mu_assert("neg chain forward failed", cmp_double_array(neg2->value, u, 3));
 
     jacobian_init(neg2);
-    neg2->eval_jacobian(neg2);
+    eval_jacobian(neg2);
 
     /* Jacobian of neg(neg(x)) is (-1)*(-1)*I = I */
     double expected_x[3] = {1.0, 1.0, 1.0};
diff --git a/tests/jacobian_tests/affine/test_promote.h b/tests/jacobian_tests/affine/test_promote.h
index c0e31ed9..47466453 100644
--- a/tests/jacobian_tests/affine/test_promote.h
+++ b/tests/jacobian_tests/affine/test_promote.h
@@ -15,7 +15,7 @@ const char *test_promote_scalar_jacobian(void)
     expr *promote_node = new_promote(var, 3, 1);
     promote_node->forward(promote_node, u);
     jacobian_init(promote_node);
-    promote_node->eval_jacobian(promote_node);
+    eval_jacobian(promote_node);
 
     /* Jacobian is 3x1 with all 1s (each output depends on same input) */
     double expected_x[3] = {1.0, 1.0, 1.0};
@@ -44,7 +44,7 @@ const char *test_promote_scalar_to_matrix_jacobian(void)
               cmp_double_array(promote_node->value, expected_val, 6));
 
     jacobian_init(promote_node);
-    promote_node->eval_jacobian(promote_node);
+    eval_jacobian(promote_node);
 
     /* Jacobian is 6x1 with all 1s (each output depends on same scalar input) */
     double expected_x[6] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
diff --git a/tests/jacobian_tests/affine/test_right_matmul.h b/tests/jacobian_tests/affine/test_right_matmul.h
index 42591722..dc571e9e 100644
--- a/tests/jacobian_tests/affine/test_right_matmul.h
+++ b/tests/jacobian_tests/affine/test_right_matmul.h
@@ -31,7 +31,7 @@ const char *test_jacobian_right_matmul_log(void)
 
     log_x_A->forward(log_x_A, x_vals);
     jacobian_init(log_x_A);
-    log_x_A->eval_jacobian(log_x_A);
+    eval_jacobian(log_x_A);
 
     /* Expected jacobian values */
     double expected_Ax[8] = {
@@ -80,7 +80,7 @@ const char *test_jacobian_right_matmul_log_vector(void)
 
     log_x_A->forward(log_x_A, x_vals);
     jacobian_init(log_x_A);
-    log_x_A->eval_jacobian(log_x_A);
+    eval_jacobian(log_x_A);
 
     /* Expected jacobian values: A^T @ diag(1/x) */
     double expected_Ax[4] = {
diff --git a/tests/jacobian_tests/affine/test_scalar_mult.h b/tests/jacobian_tests/affine/test_scalar_mult.h
index afdcdd4f..19fe8061 100644
--- a/tests/jacobian_tests/affine/test_scalar_mult.h
+++ b/tests/jacobian_tests/affine/test_scalar_mult.h
@@ -29,7 +29,7 @@ const char *test_jacobian_scalar_mult_log_vector(void)
 
     /* Initialize and evaluate jacobian */
     jacobian_init(y);
-    y->eval_jacobian(y);
+    eval_jacobian(y);
 
     /* Expected jacobian: 2.5 * [1/1, 1/2, 1/4] = [2.5, 1.25, 0.625] */
     double expected_x[3] = {2.5, 1.25, 0.625};
@@ -63,7 +63,7 @@ const char *test_jacobian_scalar_mult_log_matrix(void)
 
     /* Initialize and evaluate jacobian */
     jacobian_init(y);
-    y->eval_jacobian(y);
+    eval_jacobian(y);
 
     /* Expected jacobian: 3.0 * [1/1, 1/2, 1/4, 1/8] = [3.0, 1.5, 0.75, 0.375] */
     double expected_x[4] = {3.0, 1.5, 0.75, 0.375};
diff --git a/tests/jacobian_tests/affine/test_sum.h b/tests/jacobian_tests/affine/test_sum.h
index c2b61bd8..12cd22b4 100644
--- a/tests/jacobian_tests/affine/test_sum.h
+++ b/tests/jacobian_tests/affine/test_sum.h
@@ -21,7 +21,7 @@ const char *test_jacobian_sum_log(void)
     expr *sum_node = new_sum(log_node, -1);
     sum_node->forward(sum_node, u_vals);
     jacobian_init(sum_node);
-    sum_node->eval_jacobian(sum_node);
+    eval_jacobian(sum_node);
     double expected_Ax[3] = {1.0, 0.5, 1.0 / 3.0};
     int expected_Ap[2] = {0, 3};
     int expected_Ai[3] = {2, 3, 4};
@@ -54,7 +54,7 @@ const char *test_jacobian_sum_mult(void)
 
     sum_node->forward(sum_node, u_vals);
     jacobian_init(sum_node);
-    sum_node->eval_jacobian(sum_node);
+    eval_jacobian(sum_node);
 
     double expected_Ax[6] = {2, 3, 4, 1, 2, 3};
     int expected_Ap[2] = {0, 6}; /* 1x10 matrix: row 0 spans all 6 nonzeros */
@@ -94,7 +94,7 @@ const char *test_jacobian_sum_log_axis_0(void)
     expr *sum_node = new_sum(log_node, 0);
     sum_node->forward(sum_node, u_vals);
     jacobian_init(sum_node);
-    sum_node->eval_jacobian(sum_node);
+    eval_jacobian(sum_node);
 
     double expected_Ax[6] = {1.0, 1.0 / 3.0, 1.0 / 5.0, 0.5, 0.25, 1.0 / 6.0};
     int expected_Ap[3] = {0, 3, 6};
@@ -136,7 +136,7 @@ const char *test_jacobian_sum_add_log_axis_0(void)
 
     sum_node->forward(sum_node, u_vals);
     jacobian_init(sum_node);
-    sum_node->eval_jacobian(sum_node);
+    eval_jacobian(sum_node);
 
     /* Expected jacobian values for both rows */
     double expected_Ax[12] = {1.0, 1.0 / 3.0, 1.0 / 5.0,  /* d(sum[0])/dx */
@@ -183,7 +183,7 @@ const char *test_jacobian_sum_log_axis_1(void)
     expr *sum_node = new_sum(log_node, 1);
     sum_node->forward(sum_node, u_vals);
     jacobian_init(sum_node);
-    sum_node->eval_jacobian(sum_node);
+    eval_jacobian(sum_node);
 
     double expected_Ax[6] = {1.0, 0.5, 1.0 / 3.0, 0.25, 1.0 / 5.0, 1.0 / 6.0};
     int expected_Ap[4] = {0, 2, 4, 6};
diff --git a/tests/jacobian_tests/affine/test_trace.h b/tests/jacobian_tests/affine/test_trace.h
index 808d6343..1391e73b 100644
--- a/tests/jacobian_tests/affine/test_trace.h
+++ b/tests/jacobian_tests/affine/test_trace.h
@@ -33,7 +33,7 @@ const char *test_jacobian_trace_variable(void)
 
     trace_node->forward(trace_node, u_vals);
     jacobian_init(trace_node);
-    trace_node->eval_jacobian(trace_node);
+    eval_jacobian(trace_node);
 
     double expected_Ax[3] = {1.0, 1.0, 1.0};
     int expected_Ap[2] = {0, 3};
@@ -78,7 +78,7 @@ const char *test_jacobian_trace_composite(void)
 
     jacobian_init(trace_node);
     trace_node->forward(trace_node, u_vals);
-    trace_node->eval_jacobian(trace_node);
+    eval_jacobian(trace_node);
 
     /* Expected values: d(log(x_ii) + exp(x_ii))/dx_ii = 1/x_ii + exp(x_ii)
      * At x_00 = 1: 1/1 + exp(1) = 1 + 2.718281828...
diff --git a/tests/jacobian_tests/affine/test_transpose.h b/tests/jacobian_tests/affine/test_transpose.h
index 3302545b..40b1b763 100644
--- a/tests/jacobian_tests/affine/test_transpose.h
+++ b/tests/jacobian_tests/affine/test_transpose.h
@@ -27,7 +27,7 @@ const char *test_jacobian_transpose(void)
     double u[4] = {1, 3, 2, 4};
     transpose_AX->forward(transpose_AX, u);
     jacobian_init(transpose_AX);
-    transpose_AX->eval_jacobian(transpose_AX);
+    eval_jacobian(transpose_AX);
 
     // Jacobian of transpose_AX
     double expected_x[8] = {1, 2, 1, 2, 3, 4, 3, 4};
@@ -63,7 +63,7 @@ const char *test_jacobian_transpose_pd_preserved(void)
     double u_vals[2] = {0.5, -1.5};
     T->forward(T, u_vals);
     jacobian_init(T);
-    T->eval_jacobian(T);
+    eval_jacobian(T);
 
     /* Structural: output Jacobian must be a PD. */
     mu_assert("transpose Jacobian should be PD", T->jacobian->is_permuted_dense);
diff --git a/tests/jacobian_tests/affine/test_upper_tri.h b/tests/jacobian_tests/affine/test_upper_tri.h
index 82d74195..a0b580b7 100644
--- a/tests/jacobian_tests/affine/test_upper_tri.h
+++ b/tests/jacobian_tests/affine/test_upper_tri.h
@@ -21,7 +21,7 @@ const char *test_upper_tri_jacobian_variable(void)
 
     ut->forward(ut, u);
     jacobian_init(ut);
-    ut->eval_jacobian(ut);
+    eval_jacobian(ut);
 
     double expected_x[6] = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
     int expected_p[7] = {0, 1, 2, 3, 4, 5, 6};
@@ -53,7 +53,7 @@ const char *test_upper_tri_jacobian_of_log(void)
 
     ut->forward(ut, u);
     jacobian_init(ut);
-    ut->eval_jacobian(ut);
+    eval_jacobian(ut);
 
     double expected_x[6] = {0.2, 1.0 / 9.0, 1.0 / 13.0, 0.1, 1.0 / 14.0, 1.0 / 15.0};
     int expected_p[7] = {0, 1, 2, 3, 4, 5, 6};
diff --git a/tests/jacobian_tests/affine/test_vector_mult.h b/tests/jacobian_tests/affine/test_vector_mult.h
index 8196f9e0..0ffa348c 100644
--- a/tests/jacobian_tests/affine/test_vector_mult.h
+++ b/tests/jacobian_tests/affine/test_vector_mult.h
@@ -29,7 +29,7 @@ const char *test_jacobian_vector_mult_log_vector(void)
 
     /* Initialize and evaluate jacobian */
     jacobian_init(y);
-    y->eval_jacobian(y);
+    eval_jacobian(y);
 
     double expected_x[3] = {2.0, 1.5, 1.0};
     int expected_p[4] = {0, 1, 2, 3};
@@ -62,7 +62,7 @@ const char *test_jacobian_vector_mult_log_matrix(void)
 
     /* Initialize and evaluate jacobian */
     jacobian_init(y);
-    y->eval_jacobian(y);
+    eval_jacobian(y);
 
     double expected_x[4] = {1.5, 1.25, 0.875, 0.5625};
     int expected_p[5] = {0, 1, 2, 3, 4};
diff --git a/tests/jacobian_tests/affine/test_vstack.h b/tests/jacobian_tests/affine/test_vstack.h
index 86f4ac8d..7240082c 100644
--- a/tests/jacobian_tests/affine/test_vstack.h
+++ b/tests/jacobian_tests/affine/test_vstack.h
@@ -32,7 +32,7 @@ const char *test_jacobian_vstack_vectors(void)
 
     stack->forward(stack, u);
     jacobian_init(stack);
-    stack->eval_jacobian(stack);
+    eval_jacobian(stack);
 
     double expected_x[6] = {1.0, 0.5, 1.0 / 3.0, exp(1.0), exp(2.0), exp(3.0)};
     int expected_i[6] = {0, 1, 2, 0, 1, 2};
@@ -78,7 +78,7 @@ const char *test_jacobian_vstack_matrix(void)
 
     stack->forward(stack, u);
     jacobian_init(stack);
-    stack->eval_jacobian(stack);
+    eval_jacobian(stack);
 
     double expected_x[9] = {1.0,      0.5, exp(7.0),  1.0 / 3.0, 0.25,
                             exp(8.0), 0.2, 1.0 / 6.0, exp(9.0)};
diff --git a/tests/jacobian_tests/bivariate_full_dom/test_elementwise_mult.h b/tests/jacobian_tests/bivariate_full_dom/test_elementwise_mult.h
index b42dbe84..7fbb0a1f 100644
--- a/tests/jacobian_tests/bivariate_full_dom/test_elementwise_mult.h
+++ b/tests/jacobian_tests/bivariate_full_dom/test_elementwise_mult.h
@@ -17,7 +17,7 @@ const char *test_jacobian_elementwise_mult_1(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     double vals[6] = {y->value[0], x->value[0], y->value[1],
                       x->value[1], y->value[2], x->value[2]};
@@ -41,7 +41,7 @@ const char *test_jacobian_elementwise_mult_2(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     double vals[6] = {x->value[0], y->value[0], x->value[1],
                       y->value[1], x->value[2], y->value[2]};
@@ -97,7 +97,7 @@ const char *test_jacobian_elementwise_mult_3(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     /* Correct answer:
       0 0 19 38 0  0 0 5  15  0
@@ -136,7 +136,7 @@ const char *test_jacobian_elementwise_mult_duplicate_gathers(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     double fwd[6] = {5.0, 5.0, 12.0, 21.0, 21.0, 32.0};
     mu_assert("forward fail", cmp_double_array(node->value, fwd, 6));
@@ -178,7 +178,7 @@ const char *test_jacobian_elementwise_mult_4(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     /* Correct answer: 0 0 10 20 0  0 0 0 0 0
                        0 0 24 24 72 0 0 0 0 0
diff --git a/tests/jacobian_tests/bivariate_full_dom/test_matmul.h b/tests/jacobian_tests/bivariate_full_dom/test_matmul.h
index 6155252e..5312ea1e 100644
--- a/tests/jacobian_tests/bivariate_full_dom/test_matmul.h
+++ b/tests/jacobian_tests/bivariate_full_dom/test_matmul.h
@@ -46,7 +46,7 @@ const char *test_jacobian_matmul(void)
     /* Forward pass and jacobian initialization */
     Z->forward(Z, u_vals);
     jacobian_init(Z);
-    Z->eval_jacobian(Z);
+    eval_jacobian(Z);
 
     /* Verify sparsity pattern */
     mu_assert("Jacobian should have 18 columns", Z->jacobian->n == n_vars);
diff --git a/tests/jacobian_tests/bivariate_restricted_dom/test_quad_over_lin.h b/tests/jacobian_tests/bivariate_restricted_dom/test_quad_over_lin.h
index 1a44469e..22f09599 100644
--- a/tests/jacobian_tests/bivariate_restricted_dom/test_quad_over_lin.h
+++ b/tests/jacobian_tests/bivariate_restricted_dom/test_quad_over_lin.h
@@ -19,7 +19,7 @@ const char *test_quad_over_lin1(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     double expected_Ax[4] = {2.0 / 4.0, 4.0 / 4.0, 6.0 / 4.0, -14.0 / 16.0};
     int expected_Ap[2] = {0, 4};
@@ -43,7 +43,7 @@ const char *test_quad_over_lin2(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     double expected_Ax[4] = {-14.0 / 16.0, 2.0 / 4.0, 4.0 / 4.0, 6.0 / 4.0};
     int expected_Ap[2] = {0, 4};
@@ -80,7 +80,7 @@ const char *test_quad_over_lin3(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     double expected_vals[4] = {71.0, 94.0, 117.0, -76.25};
     int expected_Ap[2] = {0, 4};
@@ -120,7 +120,7 @@ const char *test_quad_over_lin4(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     double expected_vals[4] = {-76.25, 71.0, 94.0, 117.0};
     int expected_Ap[2] = {0, 4};
@@ -160,7 +160,7 @@ const char *test_quad_over_lin5(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     double expected_vals[7] = {12, 36, -117, 36, 84, 114, 144};
     int expected_Ap[2] = {0, 7};
diff --git a/tests/jacobian_tests/bivariate_restricted_dom/test_rel_entr.h b/tests/jacobian_tests/bivariate_restricted_dom/test_rel_entr.h
index e9998fad..1bc91cde 100644
--- a/tests/jacobian_tests/bivariate_restricted_dom/test_rel_entr.h
+++ b/tests/jacobian_tests/bivariate_restricted_dom/test_rel_entr.h
@@ -17,7 +17,7 @@ const char *test_jacobian_rel_entr_vector_args_1(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     double a = log(1.0 / 4.0) + 1.0;
     double b = log(2.0 / 5.0) + 1.0;
@@ -48,7 +48,7 @@ const char *test_jacobian_rel_entr_vector_args_2(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     double a = log(1.0 / 4.0) + 1.0;
     double b = log(2.0 / 5.0) + 1.0;
@@ -82,7 +82,7 @@ const char *test_jacobian_rel_entr_matrix_args(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     double dx0 = log(1.0 / 6.0) + 1.0;
     double dx1 = log(2.0 / 5.0) + 1.0;
diff --git a/tests/jacobian_tests/bivariate_restricted_dom/test_rel_entr_scalar_vector.h b/tests/jacobian_tests/bivariate_restricted_dom/test_rel_entr_scalar_vector.h
index 5236d3e7..4dcf6afa 100644
--- a/tests/jacobian_tests/bivariate_restricted_dom/test_rel_entr_scalar_vector.h
+++ b/tests/jacobian_tests/bivariate_restricted_dom/test_rel_entr_scalar_vector.h
@@ -15,7 +15,7 @@ const char *test_jacobian_rel_entr_scalar_vector(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     double a = log(1.0 / 1.0) + 1.0;
     double b = log(1.0 / 2.0) + 1.0;
diff --git a/tests/jacobian_tests/bivariate_restricted_dom/test_rel_entr_vector_scalar.h b/tests/jacobian_tests/bivariate_restricted_dom/test_rel_entr_vector_scalar.h
index d2f594f4..3fef0535 100644
--- a/tests/jacobian_tests/bivariate_restricted_dom/test_rel_entr_vector_scalar.h
+++ b/tests/jacobian_tests/bivariate_restricted_dom/test_rel_entr_vector_scalar.h
@@ -15,7 +15,7 @@ const char *test_jacobian_rel_entr_vector_scalar(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     double a = log(1.0 / 4.0) + 1.0;
     double b = log(2.0 / 4.0) + 1.0;
diff --git a/tests/jacobian_tests/composite/test_composite_exp.h b/tests/jacobian_tests/composite/test_composite_exp.h
index 73a2e95a..48b36434 100644
--- a/tests/jacobian_tests/composite/test_composite_exp.h
+++ b/tests/jacobian_tests/composite/test_composite_exp.h
@@ -24,7 +24,7 @@ const char *test_jacobian_composite_exp(void)
     expr *exp_node = new_exp(Au);
     jacobian_init(exp_node);
     exp_node->forward(exp_node, u_vals);
-    exp_node->eval_jacobian(exp_node);
+    eval_jacobian(exp_node);
 
     /* A*u = [10, 7], so exp(A*u) = [exp(10), exp(7)]
      * J = diag(exp(A*u)) * A */
diff --git a/tests/jacobian_tests/elementwise_restricted_dom/test_log.h b/tests/jacobian_tests/elementwise_restricted_dom/test_log.h
index b671910b..ce123706 100644
--- a/tests/jacobian_tests/elementwise_restricted_dom/test_log.h
+++ b/tests/jacobian_tests/elementwise_restricted_dom/test_log.h
@@ -17,7 +17,7 @@ const char *test_jacobian_log(void)
     expr *log_node = new_log(u);
     log_node->forward(log_node, u_vals);
     jacobian_init(log_node);
-    log_node->eval_jacobian(log_node);
+    eval_jacobian(log_node);
     mu_assert("vals fail", cmp_values(log_node->jacobian, expected_Ax, 3));
     mu_assert("sparsity fail",
               cmp_sparsity(log_node->jacobian, expected_Ap, expected_Ai, 3, 3));
@@ -35,7 +35,7 @@ const char *test_jacobian_log_matrix(void)
     expr *log_node = new_log(u);
     log_node->forward(log_node, u_vals);
     jacobian_init(log_node);
-    log_node->eval_jacobian(log_node);
+    eval_jacobian(log_node);
     mu_assert("vals fail", cmp_values(log_node->jacobian, expected_Ax, 4));
     mu_assert("sparsity fail",
               cmp_sparsity(log_node->jacobian, expected_Ap, expected_Ai, 4, 4));
diff --git a/tests/jacobian_tests/other/test_prod.h b/tests/jacobian_tests/other/test_prod.h
index 7875437a..bd7e9b95 100644
--- a/tests/jacobian_tests/other/test_prod.h
+++ b/tests/jacobian_tests/other/test_prod.h
@@ -19,7 +19,7 @@ const char *test_jacobian_prod_no_zero(void)
 
     p->forward(p, u_vals);
     jacobian_init(p);
-    p->eval_jacobian(p);
+    eval_jacobian(p);
 
     double expected_Ax[4] = {24.0, 12.0, 8.0, 6.0};
     int expected_Ap[2] = {0, 4};
@@ -44,7 +44,7 @@ const char *test_jacobian_prod_one_zero(void)
 
     p->forward(p, u_vals);
     jacobian_init(p);
-    p->eval_jacobian(p);
+    eval_jacobian(p);
 
     double expected_Ax[4] = {0.0, 12.0, 0.0, 0.0};
     int expected_Ap[2] = {0, 4};
@@ -67,7 +67,7 @@ const char *test_jacobian_prod_two_zeros(void)
 
     p->forward(p, u_vals);
     jacobian_init(p);
-    p->eval_jacobian(p);
+    eval_jacobian(p);
 
     double expected_Ax[4] = {0.0, 0.0, 0.0, 0.0};
     int expected_Ap[2] = {0, 4};
diff --git a/tests/jacobian_tests/other/test_prod_axis_one.h b/tests/jacobian_tests/other/test_prod_axis_one.h
index 8310f863..afcc5ded 100644
--- a/tests/jacobian_tests/other/test_prod_axis_one.h
+++ b/tests/jacobian_tests/other/test_prod_axis_one.h
@@ -35,7 +35,7 @@ const char *test_jacobian_prod_axis_one(void)
 
     p->forward(p, u_vals);
     jacobian_init(p);
-    p->eval_jacobian(p);
+    eval_jacobian(p);
 
     /* CSR_matrix format for 3x10 Jacobian with row-strided structure */
     double expected_Ax[9] = {28.0, 7.0, 4.0, 40.0, 16.0, 10.0, 54.0, 27.0, 18.0};
@@ -78,7 +78,7 @@ const char *test_jacobian_prod_axis_one_one_zero(void)
 
     p->forward(p, u_vals);
     jacobian_init(p);
-    p->eval_jacobian(p);
+    eval_jacobian(p);
 
     /* CSR_matrix format for 3x10 Jacobian with row-strided structure */
     double expected_Ax[9] = {28.0, 7.0, 4.0, 0.0, 16.0, 0.0, 54.0, 27.0, 18.0};
diff --git a/tests/jacobian_tests/other/test_prod_axis_zero.h b/tests/jacobian_tests/other/test_prod_axis_zero.h
index 519a40ca..ccc45ba9 100644
--- a/tests/jacobian_tests/other/test_prod_axis_zero.h
+++ b/tests/jacobian_tests/other/test_prod_axis_zero.h
@@ -29,7 +29,7 @@ const char *test_jacobian_prod_axis_zero(void)
 
     p->forward(p, u_vals);
     jacobian_init(p);
-    p->eval_jacobian(p);
+    eval_jacobian(p);
 
     /* CSR_matrix format for 3x8 Jacobian with block diagonal structure */
     double expected_Ax[6] = {2.0, 1.0, 4.0, 3.0, 6.0, 5.0};
diff --git a/tests/jacobian_tests/other/test_quad_form.h b/tests/jacobian_tests/other/test_quad_form.h
index f4c1b7e5..6eb0c705 100644
--- a/tests/jacobian_tests/other/test_quad_form.h
+++ b/tests/jacobian_tests/other/test_quad_form.h
@@ -25,7 +25,7 @@ const char *test_quad_form(void)
 
     jacobian_init(node);
     node->forward(node, u_vals);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     double expected_Ax[3] = {10.0, 16.0, 24.0};
     int expected_Ap[2] = {0, 3};
@@ -71,7 +71,7 @@ expr *node = new_quad_form_sparse(Au, Q);
 
 jacobian_init(node);
 node->forward(node, u_vals);
-node->eval_jacobian(node);
+eval_jacobian(node);
 
 double expected_Ax[5] = {422, 2222, 3244, 3786, 120};
 int expected_Ap[2] = {0, 5};
diff --git a/tests/jacobian_tests/test_values_version.h b/tests/jacobian_tests/test_values_version.h
new file mode 100644
index 00000000..74ffbba6
--- /dev/null
+++ b/tests/jacobian_tests/test_values_version.h
@@ -0,0 +1,229 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "atoms/affine.h"
+#include "atoms/elementwise_full_dom.h"
+#include "expr.h"
+#include "minunit.h"
+#include "numerical_diff.h"
+#include "test_helpers.h"
+#include "utils/sparse_matrix.h"
+
+/* Non-affine node: every eval_jacobian call bumps the jacobian's
+ * values_version. */
+const char *test_values_version_non_affine(void)
+{
+    double u[3] = {0.1, 0.2, 0.3};
+    expr *x = new_variable(3, 1, 0, 3);
+    expr *e = new_exp(x);
+
+    jacobian_init(e);
+    e->forward(e, u);
+
+    uint64_t v0 = e->jacobian->values_version;
+    eval_jacobian(e);
+    mu_assert("first eval must bump", e->jacobian->values_version == v0 + 1);
+    eval_jacobian(e);
+    mu_assert("second eval must bump", e->jacobian->values_version == v0 + 2);
+
+    free_expr(e);
+    return 0;
+}
+
+/* Affine node: the bump is skipped once the node has been evaluated this
+ * parameter epoch (values provably identical), and re-armed by
+ * expr_set_needs_refresh. */
+const char *test_values_version_affine(void)
+{
+    double u[3] = {0.1, 0.2, 0.3};
+    expr *x = new_variable(3, 1, 0, 3);
+    expr *m = new_neg(x);
+
+    jacobian_init(m);
+    m->forward(m, u);
+
+    uint64_t v0 = m->jacobian->values_version;
+    eval_jacobian(m);
+    mu_assert("first eval must bump", m->jacobian->values_version == v0 + 1);
+    eval_jacobian(m);
+    mu_assert("re-eval of affine node must not bump",
+              m->jacobian->values_version == v0 + 1);
+
+    expr_set_needs_refresh(m);
+    eval_jacobian(m);
+    mu_assert("eval after parameter refresh must bump",
+              m->jacobian->values_version == v0 + 2);
+
+    free_expr(m);
+    return 0;
+}
+
+/* sparse_matrix CSC mirror: refresh_csc_values fills once per write and
+ * dedupes repeated calls with no write in between. */
+const char *test_values_version_csc_mirror_dedup(void)
+{
+    double u[3] = {0.1, 0.2, 0.3};
+    expr *x = new_variable(3, 1, 0, 3);
+    expr *e = new_exp(x);
+
+    jacobian_init(e);
+    e->forward(e, u);
+    eval_jacobian(e);
+
+    /* jacobian of exp over a leaf variable is a diagonal sparse_matrix, so
+       the CSC mirror holds the same values in the same order */
+    double expected[3] = {exp(u[0]), exp(u[1]), exp(u[2])};
+    sparse_matrix *sm = (sparse_matrix *) e->jacobian;
+
+    e->jacobian->refresh_csc_values(e->jacobian);
+    mu_assert("csc_seen must catch up to values_version",
+              sm->csc_seen == sm->base.values_version);
+    mu_assert("csc values after first refresh",
+              cmp_double_array(sm->csc_cache->x, expected, 3));
+
+    /* no write in between: second refresh must be a no-op and leave the
+       values correct */
+    e->jacobian->refresh_csc_values(e->jacobian);
+    mu_assert("csc values after deduped refresh",
+              cmp_double_array(sm->csc_cache->x, expected, 3));
+
+    free_expr(e);
+    return 0;
+}
+
+/* stacked_pd CSR view: to_csr returns a stable pointer, skips the value
+ * copy when nothing was written, and refreshes after a write + bump. */
+const char *test_values_version_stacked_pd_to_csr(void)
+{
+    /* exp(A @ X) with X a 3x2 matrix variable: left_matmul_dense has
+       n_blocks = 2, so the child jacobian is stacked_pd and the (non-affine)
+       exp node's jacobian mirrors that type. */
+    double A[6] = {1.0, -0.5, 2.0, 0.5, 1.5, -1.0};
+    double u1[6] = {0.1, 0.2, 0.3, -0.1, -0.2, -0.3};
+    double u2[6] = {0.4, -0.3, 0.2, 0.1, 0.3, -0.4};
+
+    expr *x = new_variable(3, 2, 0, 6);
+    expr *L = new_left_matmul_dense(NULL, x, 2, 3, A);
+    expr *e = new_exp(L);
+
+    jacobian_init(e);
+    e->forward(e, u1);
+    eval_jacobian(e);
+    mu_assert("jacobian must be stacked_pd", e->jacobian->is_stacked_pd);
+
+    CSR_matrix *view = e->jacobian->to_csr(e->jacobian);
+    double *vals1 = (double *) malloc(view->nnz * sizeof(double));
+    memcpy(vals1, view->x, view->nnz * sizeof(double));
+
+    /* no write in between: same view pointer, same values */
+    mu_assert("view pointer must be stable",
+              e->jacobian->to_csr(e->jacobian) == view);
+    mu_assert("values unchanged without a write",
+              cmp_double_array(view->x, vals1, view->nnz));
+
+    /* write + bump: the view must reflect the new values */
+    e->forward(e, u2);
+    eval_jacobian(e);
+    mu_assert("view pointer must be stable after re-eval",
+              e->jacobian->to_csr(e->jacobian) == view);
+    mu_assert("values must change after a write",
+              memcmp(view->x, vals1, view->nnz * sizeof(double)) != 0);
+    mu_assert("jacobian values fresh after write+bump",
+              check_jacobian_num(e, u2, NUMERICAL_DIFF_DEFAULT_H));
+
+    free(vals1);
+    free_expr(e);
+    return 0;
+}
+
+/* Regression test for n-ary parameter-refresh propagation:
+ * expr_set_needs_refresh must reach hstack args[] children. A
+ * parameter-dependent affine child with a stacked_pd Jacobian would
+ * otherwise keep its bump-skip armed across a parameter update, and
+ * hstack's to_csr read of the child would serve the previous epoch's
+ * values. */
+const char *test_values_version_param_under_hstack(void)
+{
+    int n_vars = 4;
+    double u[4] = {1.0, 2.0, 3.0, 4.0};
+    expr *X = new_variable(2, 2, 0, n_vars);
+
+    /* A @ X with dense 2x2 A: n_blocks = 2, so the jacobian is stacked_pd */
+    double A[4] = {1.0, 2.0, 3.0, 4.0};
+    expr *AX = new_left_matmul_dense(NULL, X, 2, 2, A);
+
+    /* p * (A @ X): affine, parameter-dependent */
+    double p0 = 2.0;
+    expr *p = new_parameter(1, 1, 0, n_vars, &p0);
+    expr *pAX = new_scalar_mult(p, AX);
+
+    expr *args[2] = {pAX, X};
+    expr *h = new_hstack(args, 2, n_vars);
+
+    jacobian_init(h);
+    h->forward(h, u);
+    eval_jacobian(h);
+
+    /* the first hstack block holds pAX's jacobian, which is linear in p */
+    int nnz1 = pAX->jacobian->nnz;
+    double *expected = (double *) malloc(nnz1 * sizeof(double));
+    memcpy(expected, h->jacobian->x, nnz1 * sizeof(double));
+
+    /* problem_update_params equivalent: p = 2.0 -> 3.0 */
+    double p1 = 3.0;
+    memcpy(p->value, &p1, sizeof(double));
+    expr_set_needs_refresh(h);
+
+    h->forward(h, u);
+    eval_jacobian(h);
+
+    for (int k = 0; k < nnz1; k++)
+    {
+        expected[k] *= p1 / p0;
+    }
+    mu_assert("hstack jacobian must pick up the new parameter value",
+              cmp_double_array(h->jacobian->x, expected, nnz1));
+
+    free(expected);
+    free_expr(h);
+    return 0;
+}
+
+/* Regression test for stale spd Hessian work matrices: two chained
+ * non-affine elementwise atoms over a dense left_matmul with n_blocks > 1
+ * make hess_term1/hess_term2 stacked_pd. They are written outside the eval
+ * wrappers and read via to_csr in sum_matrices_fill_values, so a missing
+ * values_version bump would leave the second Hessian eval reading stale
+ * CSR caches. */
+const char *test_values_version_spd_hess_terms(void)
+{
+    double A[6] = {0.8, -0.4, 0.6, 0.3, 0.9, -0.7};
+    double u1[6] = {0.1, 0.2, 0.3, -0.1, -0.2, -0.3};
+    double u2[6] = {0.3, -0.2, 0.1, 0.2, -0.3, 0.4};
+    double w1[4] = {1.0, -2.0, 0.5, 1.5};
+    double w2[4] = {-0.5, 1.0, 2.0, -1.5};
+
+    expr *x = new_variable(3, 2, 0, 6);
+    expr *L = new_left_matmul_dense(NULL, x, 2, 3, A);
+    expr *inner = new_exp(L);
+    expr *outer = new_exp(inner);
+
+    /* make sure this test exercises the intended path: both Hessian work
+       matrices of the outer node must be stacked_pd */
+    jacobian_init(outer);
+    wsum_hess_init(outer);
+    mu_assert("hess_term1 must be stacked_pd",
+              outer->work->hess_term1->is_stacked_pd);
+    mu_assert("hess_term2 must be stacked_pd",
+              outer->work->hess_term2->is_stacked_pd);
+
+    mu_assert("first hessian eval",
+              check_wsum_hess(outer, u1, w1, NUMERICAL_DIFF_DEFAULT_H));
+    mu_assert("second hessian eval with different input",
+              check_wsum_hess(outer, u2, w2, NUMERICAL_DIFF_DEFAULT_H));
+
+    free_expr(outer);
+    return 0;
+}
diff --git a/tests/numerical_diff.c b/tests/numerical_diff.c
index 75faf8c6..68ca576c 100644
--- a/tests/numerical_diff.c
+++ b/tests/numerical_diff.c
@@ -68,7 +68,7 @@ int check_jacobian_num(expr *node, const double *u, double h)
 
     jacobian_init(node);
     node->forward(node, u);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
 
     double *J_num = numerical_jacobian(node, u, h);
 
@@ -131,7 +131,8 @@ double *numerical_wsum_hess(expr *node, const double *u, const double *w, double
        and permuted_dense the cached csr->x aliases the underlying value
        buffer (so the same call returns the same valid view), but for
        stacked_pd the csr_cache->x is a separate buffer that to_csr
-       memcpy-refreshes from block X each call. */
+       refreshes from block X whenever the jacobian's values_version moved
+       since the previous call. */
     CSR_matrix *jac;
 
     for (int j = 0; j < n; j++)
@@ -139,7 +140,7 @@ double *numerical_wsum_hess(expr *node, const double *u, const double *w, double
         /* g(u + h*e_j) */
         u_work[j] = u[j] + h;
         node->forward(node, u_work);
-        node->eval_jacobian(node);
+        eval_jacobian(node);
         jac = node->jacobian->to_csr(node->jacobian);
         memset(g_plus, 0, n * sizeof(double));
         csr_transpose_mult_vec(jac, w, g_plus);
@@ -147,7 +148,7 @@ double *numerical_wsum_hess(expr *node, const double *u, const double *w, double
         /* g(u - h*e_j) */
         u_work[j] = u[j] - h;
         node->forward(node, u_work);
-        node->eval_jacobian(node);
+        eval_jacobian(node);
         jac = node->jacobian->to_csr(node->jacobian);
         memset(g_minus, 0, n * sizeof(double));
         csr_transpose_mult_vec(jac, w, g_minus);
@@ -176,8 +177,8 @@ int check_wsum_hess(expr *node, const double *u, const double *w, double h)
     /* Now compute analytical (reuses jacobian from numerical) */
     wsum_hess_init(node);
     node->forward(node, u);
-    node->eval_jacobian(node);
-    node->eval_wsum_hess(node, w);
+    eval_jacobian(node);
+    eval_wsum_hess(node, w);
 
     double *H_ana = calloc((size_t) n * n, sizeof(double));
     csr_to_dense(node->wsum_hess->to_csr(node->wsum_hess), H_ana);
diff --git a/tests/profiling/profile_hessian_exp_AX.h b/tests/profiling/profile_hessian_exp_AX.h
index a5cdafc2..a426e1ba 100644
--- a/tests/profiling/profile_hessian_exp_AX.h
+++ b/tests/profiling/profile_hessian_exp_AX.h
@@ -74,8 +74,8 @@ const char *profile_hessian_exp_AX(void)
     Timer t1f;
     clock_gettime(CLOCK_MONOTONIC, &t1f.start);
     node->forward(node, X_vals);
-    node->eval_jacobian(node);
-    node->eval_wsum_hess(node, w);
+    eval_jacobian(node);
+    eval_wsum_hess(node, w);
     clock_gettime(CLOCK_MONOTONIC, &t1f.end);
 
     matrix *H1 = node->wsum_hess;
diff --git a/tests/profiling/profile_left_matmul.h b/tests/profiling/profile_left_matmul.h
index db6b8328..a206818a 100644
--- a/tests/profiling/profile_left_matmul.h
+++ b/tests/profiling/profile_left_matmul.h
@@ -51,7 +51,7 @@ const char *profile_left_matmul(void)
     printf("left_matmul jacobian init time: %8.3f seconds\n",
            GET_ELAPSED_SECONDS(timer));
     clock_gettime(CLOCK_MONOTONIC, &timer.start);
-    AX->eval_jacobian(AX);
+    eval_jacobian(AX);
     clock_gettime(CLOCK_MONOTONIC, &timer.end);
     printf("left_matmul jacobian eval time: %8.3f seconds\n",
            GET_ELAPSED_SECONDS(timer));
diff --git a/tests/profiling/profile_log_reg.h b/tests/profiling/profile_log_reg.h
index 1bba52b7..7fcc69e6 100644
--- a/tests/profiling/profile_log_reg.h
+++ b/tests/profiling/profile_log_reg.h
@@ -56,10 +56,10 @@ const char *profile_log_reg(void)
     Timer t_a_jac, t_a_hess;
     double w_one = 1.0;
     clock_gettime(CLOCK_MONOTONIC, &t_a_jac.start);
-    obj->eval_jacobian(obj);
+    eval_jacobian(obj);
     clock_gettime(CLOCK_MONOTONIC, &t_a_jac.end);
     clock_gettime(CLOCK_MONOTONIC, &t_a_hess.start);
-    obj->eval_wsum_hess(obj, &w_one);
+    eval_wsum_hess(obj, &w_one);
     clock_gettime(CLOCK_MONOTONIC, &t_a_hess.end);
     double sec_a_jac = GET_ELAPSED_SECONDS(t_a_jac);
     double sec_a_hess = GET_ELAPSED_SECONDS(t_a_hess);
diff --git a/tests/profiling/profile_trimmed_log_reg.h b/tests/profiling/profile_trimmed_log_reg.h
index 3fe2ae60..3583905e 100644
--- a/tests/profiling/profile_trimmed_log_reg.h
+++ b/tests/profiling/profile_trimmed_log_reg.h
@@ -67,14 +67,14 @@ const char *profile_trimmed_log_reg(void)
     double w_one = 1.0;
     Timer t_jac, t_hess;
     clock_gettime(CLOCK_MONOTONIC, &t_jac.start);
-    obj->eval_jacobian(obj);
+    eval_jacobian(obj);
     clock_gettime(CLOCK_MONOTONIC, &t_jac.end);
 
-    obj->eval_wsum_hess(obj, &w_one); /* warm-up */
+    eval_wsum_hess(obj, &w_one); /* warm-up */
     clock_gettime(CLOCK_MONOTONIC, &t_hess.start);
     for (int it = 0; it < N_HESS_ITERS; it++)
     {
-        obj->eval_wsum_hess(obj, &w_one);
+        eval_wsum_hess(obj, &w_one);
     }
     clock_gettime(CLOCK_MONOTONIC, &t_hess.end);
 
diff --git a/tests/wsum_hess/affine/test_broadcast.h b/tests/wsum_hess/affine/test_broadcast.h
index 6635a7de..2428e1c0 100644
--- a/tests/wsum_hess/affine/test_broadcast.h
+++ b/tests/wsum_hess/affine/test_broadcast.h
@@ -39,7 +39,7 @@ const char *test_wsum_hess_broadcast_row(void)
      *      col0      col1         col2
      */
     double w[6] = {1.0, 0.5, 2.0, 1.0, 0.25, 0.125};
-    bcast->eval_wsum_hess(bcast, w);
+    eval_wsum_hess(bcast, w);
 
     /* For broadcast_row, weights are summed across the m replicas:
      * Accumulated weights for log(x):
@@ -91,7 +91,7 @@ const char *test_wsum_hess_broadcast_col(void)
      *      col0            col1
      */
     double w[6] = {1.0, 0.5, 0.25, 2.0, 1.0, 0.5};
-    bcast->eval_wsum_hess(bcast, w);
+    eval_wsum_hess(bcast, w);
 
     /* For broadcast_col, weights are summed across the n replicas:
      * Accumulated weights for log(x):
@@ -141,7 +141,7 @@ const char *test_wsum_hess_broadcast_scalar_to_matrix(void)
      * w = [1.0, 0.5, 2.0, 1.0, 0.25, 0.125]
      */
     double w[6] = {1.0, 0.5, 2.0, 1.0, 0.25, 0.125};
-    bcast->eval_wsum_hess(bcast, w);
+    eval_wsum_hess(bcast, w);
 
     /* For broadcast_scalar, all weights are summed:
      * w_acc[0] = sum(w) = 1.0 + 0.5 + 2.0 + 1.0 + 0.25 + 0.125 = 4.875
diff --git a/tests/wsum_hess/affine/test_convolve.h b/tests/wsum_hess/affine/test_convolve.h
index 98dd7ab0..496e569f 100644
--- a/tests/wsum_hess/affine/test_convolve.h
+++ b/tests/wsum_hess/affine/test_convolve.h
@@ -24,7 +24,7 @@ const char *test_wsum_hess_convolve(void)
     y->forward(y, u);
     jacobian_init(y);
     wsum_hess_init(y);
-    y->eval_wsum_hess(y, w);
+    eval_wsum_hess(y, w);
 
     mu_assert("Convolve wsum_hess should be 3x3", y->wsum_hess->m == 3);
     mu_assert("Convolve wsum_hess should be square", y->wsum_hess->n == 3);
diff --git a/tests/wsum_hess/affine/test_diag_mat.h b/tests/wsum_hess/affine/test_diag_mat.h
index 2774d5f7..58fa5eb6 100644
--- a/tests/wsum_hess/affine/test_diag_mat.h
+++ b/tests/wsum_hess/affine/test_diag_mat.h
@@ -30,7 +30,7 @@ const char *test_wsum_hess_diag_mat_log(void)
     dm->forward(dm, u);
     jacobian_init(dm);
     wsum_hess_init(dm);
-    dm->eval_wsum_hess(dm, w);
+    eval_wsum_hess(dm, w);
 
     double expected_x[4] = {-1.0, 0.0, 0.0, -0.0625};
     int expected_p[5] = {0, 1, 2, 3, 4};
diff --git a/tests/wsum_hess/affine/test_hstack.h b/tests/wsum_hess/affine/test_hstack.h
index 49809a9e..2b6666bd 100644
--- a/tests/wsum_hess/affine/test_hstack.h
+++ b/tests/wsum_hess/affine/test_hstack.h
@@ -38,7 +38,7 @@ const char *test_wsum_hess_hstack(void)
     hstack_node->forward(hstack_node, u_vals);
     jacobian_init(hstack_node);
     wsum_hess_init(hstack_node);
-    hstack_node->eval_wsum_hess(hstack_node, w);
+    eval_wsum_hess(hstack_node, w);
 
     /* Expected Hessian:
      * log(x): d²/dx² = -1/x²
@@ -134,7 +134,7 @@ const char *test_wsum_hess_hstack_matrix(void)
     hstack_node->forward(hstack_node, u_vals);
     jacobian_init(hstack_node);
     wsum_hess_init(hstack_node);
-    hstack_node->eval_wsum_hess(hstack_node, w);
+    eval_wsum_hess(hstack_node, w);
 
     /* Expected Hessian (diagonal):
      * log(x): w[0:5] * (-1/x[0:5]²) at indices 0-5
diff --git a/tests/wsum_hess/affine/test_index.h b/tests/wsum_hess/affine/test_index.h
index 6d77a394..0685597b 100644
--- a/tests/wsum_hess/affine/test_index.h
+++ b/tests/wsum_hess/affine/test_index.h
@@ -27,7 +27,7 @@ const char *test_wsum_hess_index_log(void)
     idx->forward(idx, u);
     jacobian_init(idx);
     wsum_hess_init(idx);
-    idx->eval_wsum_hess(idx, w);
+    eval_wsum_hess(idx, w);
 
     /* Expected diagonal values:
      * H[0,0] = -1 * 1/1^2 = -1.0
@@ -60,7 +60,7 @@ const char *test_wsum_hess_index_repeated(void)
     idx->forward(idx, u);
     jacobian_init(idx);
     wsum_hess_init(idx);
-    idx->eval_wsum_hess(idx, w);
+    eval_wsum_hess(idx, w);
 
     /* Hessian of log at x=2 is -1/4
      * weighted by 3 (accumulated) -> -3/4 = -0.75
@@ -95,7 +95,7 @@ const char *test_wsum_hess_sum_index_log(void)
     sum_node->forward(sum_node, u);
     jacobian_init(sum_node);
     wsum_hess_init(sum_node);
-    sum_node->eval_wsum_hess(sum_node, &w);
+    eval_wsum_hess(sum_node, &w);
 
     /* Expected diagonal values:
      * H[0,0] = -1 * 1/1^2 = -1.0
diff --git a/tests/wsum_hess/affine/test_kron.h b/tests/wsum_hess/affine/test_kron.h
index 801db65c..8e7fb9e3 100644
--- a/tests/wsum_hess/affine/test_kron.h
+++ b/tests/wsum_hess/affine/test_kron.h
@@ -24,7 +24,7 @@ const char *test_wsum_hess_kron(void)
     Z->forward(Z, u);
     jacobian_init(Z);
     wsum_hess_init(Z);
-    Z->eval_wsum_hess(Z, w);
+    eval_wsum_hess(Z, w);
 
     mu_assert("kron wsum_hess square", Z->wsum_hess->m == 4 && Z->wsum_hess->n == 4);
     mu_assert("kron wsum_hess zero for linear arg", Z->wsum_hess->nnz == 0);
diff --git a/tests/wsum_hess/affine/test_left_matmul.h b/tests/wsum_hess/affine/test_left_matmul.h
index 4869dbad..56286bf2 100644
--- a/tests/wsum_hess/affine/test_left_matmul.h
+++ b/tests/wsum_hess/affine/test_left_matmul.h
@@ -68,7 +68,7 @@ const char *test_wsum_hess_left_matmul(void)
     A_log_x->forward(A_log_x, x_vals);
     jacobian_init(A_log_x);
     wsum_hess_init(A_log_x);
-    A_log_x->eval_wsum_hess(A_log_x, w);
+    eval_wsum_hess(A_log_x, w);
 
     /* Expected wsum_hess: diagonal matrix with all 3 entries
      * (sparsity matches child's diagonal Hessian) */
@@ -174,7 +174,7 @@ const char *test_wsum_hess_left_matmul_matrix(void)
     A_log_x->forward(A_log_x, x_vals);
     jacobian_init(A_log_x);
     wsum_hess_init(A_log_x);
-    A_log_x->eval_wsum_hess(A_log_x, w);
+    eval_wsum_hess(A_log_x, w);
 
     /* Expected wsum_hess: 6x6 diagonal matrix with all 6 entries */
     double expected_x[6] = {
diff --git a/tests/wsum_hess/affine/test_right_matmul.h b/tests/wsum_hess/affine/test_right_matmul.h
index 661e8e0f..bc66b5fa 100644
--- a/tests/wsum_hess/affine/test_right_matmul.h
+++ b/tests/wsum_hess/affine/test_right_matmul.h
@@ -38,7 +38,7 @@ const char *test_wsum_hess_right_matmul(void)
     log_x_A->forward(log_x_A, x_vals);
     jacobian_init(log_x_A);
     wsum_hess_init(log_x_A);
-    log_x_A->eval_wsum_hess(log_x_A, w);
+    eval_wsum_hess(log_x_A, w);
 
     /* Expected wsum_hess: diagonal matrix with 4 entries */
     double expected_x[4] = {
@@ -87,7 +87,7 @@ const char *test_wsum_hess_right_matmul_vector(void)
     log_x_A->forward(log_x_A, x_vals);
     jacobian_init(log_x_A);
     wsum_hess_init(log_x_A);
-    log_x_A->eval_wsum_hess(log_x_A, w);
+    eval_wsum_hess(log_x_A, w);
 
     /* Expected wsum_hess: diagonal matrix with 3 entries */
     double expected_x[3] = {
diff --git a/tests/wsum_hess/affine/test_scalar_mult.h b/tests/wsum_hess/affine/test_scalar_mult.h
index d08e1343..2903f7f1 100644
--- a/tests/wsum_hess/affine/test_scalar_mult.h
+++ b/tests/wsum_hess/affine/test_scalar_mult.h
@@ -32,7 +32,7 @@ const char *test_wsum_hess_scalar_mult_log_vector(void)
     jacobian_init(y);
     wsum_hess_init(y);
     double w[3] = {1.0, 0.5, 0.25};
-    y->eval_wsum_hess(y, w);
+    eval_wsum_hess(y, w);
 
     double expected_x[3] = {-2.5, -0.3125, -0.0390625};
     int expected_p[4] = {0, 1, 2, 3};
@@ -67,7 +67,7 @@ const char *test_wsum_hess_scalar_mult_log_matrix(void)
     jacobian_init(y);
     wsum_hess_init(y);
     double w[4] = {1.0, 1.0, 1.0, 1.0};
-    y->eval_wsum_hess(y, w);
+    eval_wsum_hess(y, w);
 
     double expected_x[4] = {-3.0, -0.75, -0.1875, -0.046875};
     int expected_p[5] = {0, 1, 2, 3, 4};
diff --git a/tests/wsum_hess/affine/test_sum.h b/tests/wsum_hess/affine/test_sum.h
index 2c739c94..b6004e5f 100644
--- a/tests/wsum_hess/affine/test_sum.h
+++ b/tests/wsum_hess/affine/test_sum.h
@@ -54,7 +54,7 @@ const char *test_wsum_hess_sum_log_axis0(void)
     sum_node->forward(sum_node, x);
     jacobian_init(sum_node);
     wsum_hess_init(sum_node);
-    sum_node->eval_wsum_hess(sum_node, w);
+    eval_wsum_hess(sum_node, w);
 
     /* Expected diagonal values */
     double expected_x[6] = {-w[0] / (x[0] * x[0]), -w[0] / (x[1] * x[1]),
@@ -89,7 +89,7 @@ const char *test_wsum_hess_sum_log_axis1(void)
     sum_node->forward(sum_node, x);
     jacobian_init(sum_node);
     wsum_hess_init(sum_node);
-    sum_node->eval_wsum_hess(sum_node, w);
+    eval_wsum_hess(sum_node, w);
 
     /* Expected diagonal values */
     double expected_x[6] = {-w[0] / (x[0] * x[0]), -w[1] / (x[1] * x[1]),
diff --git a/tests/wsum_hess/affine/test_trace.h b/tests/wsum_hess/affine/test_trace.h
index bfbb5e07..75a4510b 100644
--- a/tests/wsum_hess/affine/test_trace.h
+++ b/tests/wsum_hess/affine/test_trace.h
@@ -32,7 +32,7 @@ const char *test_wsum_hess_trace_variable(void)
     trace_node->forward(trace_node, u_vals);
     jacobian_init(trace_node);
     wsum_hess_init(trace_node);
-    trace_node->eval_wsum_hess(trace_node, &w);
+    eval_wsum_hess(trace_node, &w);
 
     /* For a linear operation (variable), Hessian is zero */
     mu_assert("wsum_hess should be empty", trace_node->wsum_hess->nnz == 0);
@@ -64,7 +64,7 @@ const char *test_wsum_hess_trace_log_variable(void)
     trace_node->forward(trace_node, u_vals);
     jacobian_init(trace_node);
     wsum_hess_init(trace_node);
-    trace_node->eval_wsum_hess(trace_node, &w);
+    eval_wsum_hess(trace_node, &w);
 
     double expected_Ax[9] = {-2.0, 0, 0, 0, -0.08, 0, 0, 0, -0.024691358024691357};
     int expected_Ap[14] = {0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 9, 9, 9};
@@ -111,7 +111,7 @@ const char *test_wsum_hess_trace_composite(void)
     trace_node->forward(trace_node, u_vals);
     jacobian_init(trace_node);
     wsum_hess_init(trace_node);
-    trace_node->eval_wsum_hess(trace_node, &w);
+    eval_wsum_hess(trace_node, &w);
 
     /* Expected diagonal Hessian values at indices [1,1], [5,5], [9,9]
      * d²(log(x_ii) + exp(x_ii))/dx_ii² = -1/x_ii² + exp(x_ii)
diff --git a/tests/wsum_hess/affine/test_transpose.h b/tests/wsum_hess/affine/test_transpose.h
index ad4e93f1..2f4c1b51 100644
--- a/tests/wsum_hess/affine/test_transpose.h
+++ b/tests/wsum_hess/affine/test_transpose.h
@@ -21,7 +21,7 @@ const char *test_wsum_hess_transpose(void)
     jacobian_init(XYT);
     wsum_hess_init(XYT);
     double w[4] = {1, 2, 3, 4};
-    XYT->eval_wsum_hess(XYT, w);
+    eval_wsum_hess(XYT, w);
 
     double expected_x[16] = {1, 2, 3, 4, 1, 2, 3, 4, 1, 3, 1, 3, 2, 4, 2, 4};
     int expected_p[9] = {0, 2, 4, 6, 8, 10, 12, 14, 16};
diff --git a/tests/wsum_hess/affine/test_upper_tri.h b/tests/wsum_hess/affine/test_upper_tri.h
index 0cff5b6a..2b9710e2 100644
--- a/tests/wsum_hess/affine/test_upper_tri.h
+++ b/tests/wsum_hess/affine/test_upper_tri.h
@@ -40,7 +40,7 @@ const char *test_wsum_hess_upper_tri_log(void)
     ut->forward(ut, u);
     jacobian_init(ut);
     wsum_hess_init(ut);
-    ut->eval_wsum_hess(ut, w);
+    eval_wsum_hess(ut, w);
 
     double expected_x[16] = {0.0,          0.0,          0.0,          0.0,
                              -1.0 / 25.0,  0.0,          0.0,          0.0,
diff --git a/tests/wsum_hess/affine/test_vector_mult.h b/tests/wsum_hess/affine/test_vector_mult.h
index 5d0d1404..e5e7bbb5 100644
--- a/tests/wsum_hess/affine/test_vector_mult.h
+++ b/tests/wsum_hess/affine/test_vector_mult.h
@@ -32,7 +32,7 @@ const char *test_wsum_hess_vector_mult_log_vector(void)
     jacobian_init(y);
     wsum_hess_init(y);
     double w[3] = {1.0, 0.5, 0.25};
-    y->eval_wsum_hess(y, w);
+    eval_wsum_hess(y, w);
 
     double expected_x[3] = {-2.0, -0.375, -0.0625};
     int expected_p[4] = {0, 1, 2, 3};
@@ -67,7 +67,7 @@ const char *test_wsum_hess_vector_mult_log_matrix(void)
     jacobian_init(y);
     wsum_hess_init(y);
     double w[4] = {1.0, 1.0, 1.0, 1.0};
-    y->eval_wsum_hess(y, w);
+    eval_wsum_hess(y, w);
 
     double expected_x[4] = {-1.5, -0.625, -0.21875, -0.0703125};
     int expected_p[5] = {0, 1, 2, 3, 4};
diff --git a/tests/wsum_hess/affine/test_vstack.h b/tests/wsum_hess/affine/test_vstack.h
index cb473a19..ee1ae1cf 100644
--- a/tests/wsum_hess/affine/test_vstack.h
+++ b/tests/wsum_hess/affine/test_vstack.h
@@ -32,7 +32,7 @@ const char *test_wsum_hess_vstack_vectors(void)
     stack->forward(stack, u);
     jacobian_init(stack);
     wsum_hess_init(stack);
-    stack->eval_wsum_hess(stack, w);
+    eval_wsum_hess(stack, w);
 
     double expected_x[3] = {-1.0 + 4.0 * exp(1.0), -0.5 + 5.0 * exp(2.0),
                             -1.0 / 3.0 + 6.0 * exp(3.0)};
@@ -85,7 +85,7 @@ const char *test_wsum_hess_vstack_matrix(void)
     stack->forward(stack, u);
     jacobian_init(stack);
     wsum_hess_init(stack);
-    stack->eval_wsum_hess(stack, w);
+    eval_wsum_hess(stack, w);
 
     double expected_x[9] = {-1.0,            /* x0: w[0]*(-1/1) */
                             -0.5,            /* x1: w[1]*(-1/4) */
diff --git a/tests/wsum_hess/bivariate_full_dom/test_matmul.h b/tests/wsum_hess/bivariate_full_dom/test_matmul.h
index 9241609d..2eaaf386 100644
--- a/tests/wsum_hess/bivariate_full_dom/test_matmul.h
+++ b/tests/wsum_hess/bivariate_full_dom/test_matmul.h
@@ -46,7 +46,7 @@ const char *test_wsum_hess_matmul(void)
     Z->forward(Z, u_vals);
     jacobian_init(Z);
     wsum_hess_init(Z);
-    Z->eval_wsum_hess(Z, w);
+    eval_wsum_hess(Z, w);
 
     /* Verify Hessian dimensions and sparsity */
     mu_assert("Hessian should be 18 cols", Z->wsum_hess->n == n_vars);
@@ -140,7 +140,7 @@ const char *test_wsum_hess_matmul_yx(void)
     Z->forward(Z, u_vals);
     jacobian_init(Z);
     wsum_hess_init(Z);
-    Z->eval_wsum_hess(Z, w);
+    eval_wsum_hess(Z, w);
 
     /* Verify Hessian dimensions and sparsity */
     mu_assert("Hessian should be 18 cols", Z->wsum_hess->n == n_vars);
diff --git a/tests/wsum_hess/bivariate_full_dom/test_multiply.h b/tests/wsum_hess/bivariate_full_dom/test_multiply.h
index 36786403..b180c3d7 100644
--- a/tests/wsum_hess/bivariate_full_dom/test_multiply.h
+++ b/tests/wsum_hess/bivariate_full_dom/test_multiply.h
@@ -24,7 +24,7 @@ const char *test_wsum_hess_multiply_1(void)
     node->forward(node, u_vals);
     jacobian_init(node);
     wsum_hess_init(node);
-    node->eval_wsum_hess(node, w);
+    eval_wsum_hess(node, w);
 
     int expected_p[13] = {0, 0, 0, 0, 1, 2, 3, 3, 3, 4, 5, 6, 6};
     int expected_i[6] = {8, 9, 10, 3, 4, 5};
@@ -83,7 +83,7 @@ const char *test_wsum_hess_multiply_sparse_random(void)
     jacobian_init(mult_node);
     wsum_hess_init(mult_node);
     double w[5] = {0.50646339, 0.44756224, 0.67295241, 0.16424956, 0.03031469};
-    mult_node->eval_wsum_hess(mult_node, w);
+    eval_wsum_hess(mult_node, w);
 
     /* Expected Hessian in CSR_matrix format (10x10) */
     int expected_p[11] = {0, 6, 9, 13, 18, 19, 20, 20, 22, 25, 29};
@@ -164,7 +164,7 @@ const char *test_wsum_hess_multiply_linear_ops(void)
 
     /* Evaluate Hessian with weights */
     double w[4] = {1.0, 2.0, 3.0, 4.0};
-    mult_node->eval_wsum_hess(mult_node, w);
+    eval_wsum_hess(mult_node, w);
 
     /* Check sparsity pattern and values */
     /* Expected CSR_matrix format:
@@ -206,9 +206,9 @@ const char *test_wsum_hess_multiply_duplicate_gathers(void)
 
     node->forward(node, u_vals);
     jacobian_init(node);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
     wsum_hess_init(node);
-    node->eval_wsum_hess(node, w);
+    eval_wsum_hess(node, w);
 
     int expected_p[9] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
     int expected_i[8] = {4, 5, 6, 7, 0, 1, 2, 3};
@@ -241,7 +241,7 @@ const char *test_wsum_hess_multiply_2(void)
     jacobian_init(node);
     wsum_hess_init(node);
     node->forward(node, u_vals);
-    node->eval_wsum_hess(node, w);
+    eval_wsum_hess(node, w);
 
     int expected_p[13] = {0, 0, 0, 0, 1, 2, 3, 3, 3, 4, 5, 6, 6};
     int expected_i[6] = {8, 9, 10, 3, 4, 5};
diff --git a/tests/wsum_hess/bivariate_restricted_dom/test_quad_over_lin.h b/tests/wsum_hess/bivariate_restricted_dom/test_quad_over_lin.h
index 2cb5dc44..cac8cda0 100644
--- a/tests/wsum_hess/bivariate_restricted_dom/test_quad_over_lin.h
+++ b/tests/wsum_hess/bivariate_restricted_dom/test_quad_over_lin.h
@@ -19,7 +19,7 @@ const char *test_wsum_hess_quad_over_lin_xy(void)
     node->forward(node, u_vals);
     jacobian_init(node);
     wsum_hess_init(node);
-    node->eval_wsum_hess(node, &w);
+    eval_wsum_hess(node, &w);
 
     int expected_p[10] = {0, 0, 0, 2, 4, 6, 6, 6, 10, 10};
     int expected_i[10] = {2, 7, 3, 7, 4, 7, 2, 3, 4, 7};
@@ -49,7 +49,7 @@ const char *test_wsum_hess_quad_over_lin_yx(void)
     node->forward(node, u_vals);
     jacobian_init(node);
     wsum_hess_init(node);
-    node->eval_wsum_hess(node, &w);
+    eval_wsum_hess(node, &w);
 
     int expected_p[10] = {0, 0, 0, 4, 4, 4, 6, 8, 10, 10};
     int expected_i[10] = {2, 5, 6, 7, 2, 5, 2, 6, 2, 7};
diff --git a/tests/wsum_hess/bivariate_restricted_dom/test_rel_entr.h b/tests/wsum_hess/bivariate_restricted_dom/test_rel_entr.h
index 5c3addeb..0c6a344e 100644
--- a/tests/wsum_hess/bivariate_restricted_dom/test_rel_entr.h
+++ b/tests/wsum_hess/bivariate_restricted_dom/test_rel_entr.h
@@ -23,7 +23,7 @@ const char *test_wsum_hess_rel_entr_1(void)
     node->forward(node, u_vals);
     jacobian_init(node);
     wsum_hess_init(node);
-    node->eval_wsum_hess(node, w);
+    eval_wsum_hess(node, w);
 
     int expected_p[11] = {0, 0, 2, 4, 6, 6, 6, 8, 10, 12, 12};
     int expected_i[12] = {1, 6, 2, 7, 3, 8, 1, 6, 2, 7, 3, 8};
@@ -55,7 +55,7 @@ const char *test_wsum_hess_rel_entr_2(void)
     node->forward(node, u_vals);
     jacobian_init(node);
     wsum_hess_init(node);
-    node->eval_wsum_hess(node, w);
+    eval_wsum_hess(node, w);
 
     int expected_p[11] = {0, 0, 2, 4, 6, 6, 6, 8, 10, 12, 12};
     int expected_i[12] = {1, 6, 2, 7, 3, 8, 1, 6, 2, 7, 3, 8};
@@ -87,7 +87,7 @@ const char *test_wsum_hess_rel_entr_matrix(void)
     node->forward(node, u_vals);
     jacobian_init(node);
     wsum_hess_init(node);
-    node->eval_wsum_hess(node, w);
+    eval_wsum_hess(node, w);
 
     int expected_p[13] = {0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24};
     int expected_i[24] = {0, 6, 1, 7, 2, 8, 3, 9, 4, 10, 5, 11,
diff --git a/tests/wsum_hess/bivariate_restricted_dom/test_rel_entr_scalar_vector.h b/tests/wsum_hess/bivariate_restricted_dom/test_rel_entr_scalar_vector.h
index 74ba8b79..118b9fdb 100644
--- a/tests/wsum_hess/bivariate_restricted_dom/test_rel_entr_scalar_vector.h
+++ b/tests/wsum_hess/bivariate_restricted_dom/test_rel_entr_scalar_vector.h
@@ -18,7 +18,7 @@ const char *test_wsum_hess_rel_entr_scalar_vector(void)
     node->forward(node, u_vals);
     jacobian_init(node);
     wsum_hess_init(node);
-    node->eval_wsum_hess(node, w);
+    eval_wsum_hess(node, w);
 
     int expected_p[5] = {0, 4, 6, 8, 10};
     int expected_i[10] = {0, 1, 2, 3, 0, 1, 0, 2, 0, 3};
diff --git a/tests/wsum_hess/bivariate_restricted_dom/test_rel_entr_vector_scalar.h b/tests/wsum_hess/bivariate_restricted_dom/test_rel_entr_vector_scalar.h
index d74e92b7..12f89e3b 100644
--- a/tests/wsum_hess/bivariate_restricted_dom/test_rel_entr_vector_scalar.h
+++ b/tests/wsum_hess/bivariate_restricted_dom/test_rel_entr_vector_scalar.h
@@ -18,7 +18,7 @@ const char *test_wsum_hess_rel_entr_vector_scalar(void)
     node->forward(node, u_vals);
     jacobian_init(node);
     wsum_hess_init(node);
-    node->eval_wsum_hess(node, w);
+    eval_wsum_hess(node, w);
 
     int expected_p[5] = {0, 2, 4, 6, 10};
     int expected_i[10] = {0, 3, 1, 3, 2, 3, 0, 1, 2, 3};
diff --git a/tests/wsum_hess/composite/test_chain_rule_wsum_hess.h b/tests/wsum_hess/composite/test_chain_rule_wsum_hess.h
index ed2c405a..3c95daa5 100644
--- a/tests/wsum_hess/composite/test_chain_rule_wsum_hess.h
+++ b/tests/wsum_hess/composite/test_chain_rule_wsum_hess.h
@@ -93,8 +93,8 @@ const char *test_wsum_hess_sum_exp_left_matmul_dense_transpose(void)
     jacobian_init(node);
     wsum_hess_init(node);
     node->forward(node, u_vals);
-    node->eval_jacobian(node);
-    node->eval_wsum_hess(node, &w);
+    eval_jacobian(node);
+    eval_wsum_hess(node, &w);
 
     CSR_matrix *H = node->wsum_hess->to_csr(node->wsum_hess);
     double dense[16] = {0};
diff --git a/tests/wsum_hess/elementwise_full_dom/test_exp.h b/tests/wsum_hess/elementwise_full_dom/test_exp.h
index 891eae0e..05b97f14 100644
--- a/tests/wsum_hess/elementwise_full_dom/test_exp.h
+++ b/tests/wsum_hess/elementwise_full_dom/test_exp.h
@@ -20,7 +20,7 @@ const char *test_wsum_hess_exp(void)
     exp_node->forward(exp_node, u_vals);
     jacobian_init(exp_node);
     wsum_hess_init(exp_node);
-    exp_node->eval_wsum_hess(exp_node, w);
+    eval_wsum_hess(exp_node, w);
 
     /* Expected values on the diagonal: w_i * exp(x_i) */
     double expected_x[3] = {1.0 * exp(1.0), 2.0 * exp(2.0), 3.0 * exp(3.0)};
diff --git a/tests/wsum_hess/elementwise_full_dom/test_hyperbolic.h b/tests/wsum_hess/elementwise_full_dom/test_hyperbolic.h
index cf543a71..98a006db 100644
--- a/tests/wsum_hess/elementwise_full_dom/test_hyperbolic.h
+++ b/tests/wsum_hess/elementwise_full_dom/test_hyperbolic.h
@@ -27,7 +27,7 @@ const char *test_wsum_hess_sinh(void)
     sinh_node->forward(sinh_node, u_vals);
     jacobian_init(sinh_node);
     wsum_hess_init(sinh_node);
-    sinh_node->eval_wsum_hess(sinh_node, w);
+    eval_wsum_hess(sinh_node, w);
 
     /* Expected values on the diagonal: w_i * sinh(x_i) */
     double expected_x[3] = {1.0 * sinh(1.0), 2.0 * sinh(2.0), 3.0 * sinh(3.0)};
@@ -60,7 +60,7 @@ const char *test_wsum_hess_tanh(void)
     tanh_node->forward(tanh_node, u_vals);
     jacobian_init(tanh_node);
     wsum_hess_init(tanh_node);
-    tanh_node->eval_wsum_hess(tanh_node, w);
+    eval_wsum_hess(tanh_node, w);
 
     /* Expected values on the diagonal: w_i * (-2*tanh(x_i)/cosh^2(x_i)) */
     double expected_x[3] = {1.0 * (-2.0 * tanh(1.0) / pow(cosh(1.0), 2)),
@@ -95,7 +95,7 @@ const char *test_wsum_hess_asinh(void)
     asinh_node->forward(asinh_node, u_vals);
     jacobian_init(asinh_node);
     wsum_hess_init(asinh_node);
-    asinh_node->eval_wsum_hess(asinh_node, w);
+    eval_wsum_hess(asinh_node, w);
 
     /* Expected values on the diagonal: w_i * (-x_i/(1+x_i^2)^(3/2)) */
     double expected_x[3] = {1.0 * (-1.0 / pow(1.0 + 1.0 * 1.0, 1.5)),
@@ -131,7 +131,7 @@ const char *test_wsum_hess_atanh(void)
     atanh_node->forward(atanh_node, u_vals);
     jacobian_init(atanh_node);
     wsum_hess_init(atanh_node);
-    atanh_node->eval_wsum_hess(atanh_node, w);
+    eval_wsum_hess(atanh_node, w);
 
     /* Expected values on the diagonal: w_i * (2*x_i/(1-x_i^2)^2) */
     double expected_x[3] = {1.0 * (2.0 * 0.1 / pow(1.0 - 0.1 * 0.1, 2)),
diff --git a/tests/wsum_hess/elementwise_full_dom/test_logistic.h b/tests/wsum_hess/elementwise_full_dom/test_logistic.h
index 025ec8cb..097aa6c9 100644
--- a/tests/wsum_hess/elementwise_full_dom/test_logistic.h
+++ b/tests/wsum_hess/elementwise_full_dom/test_logistic.h
@@ -27,10 +27,10 @@ const char *test_wsum_hess_logistic(void)
     expr *logistic_node = new_logistic(x);
     logistic_node->forward(logistic_node, u_vals);
     jacobian_init(logistic_node);
-    logistic_node->eval_jacobian(logistic_node);
+    eval_jacobian(logistic_node);
     jacobian_init(logistic_node);
     wsum_hess_init(logistic_node);
-    logistic_node->eval_wsum_hess(logistic_node, w);
+    eval_wsum_hess(logistic_node, w);
 
     /* Expected values on the diagonal: w_i * σ(x_i) * (1 - σ(x_i)) */
     double sigma1 = 1.0 / (1.0 + exp(-1.0));
diff --git a/tests/wsum_hess/elementwise_full_dom/test_power.h b/tests/wsum_hess/elementwise_full_dom/test_power.h
index 1456bef7..f841c93a 100644
--- a/tests/wsum_hess/elementwise_full_dom/test_power.h
+++ b/tests/wsum_hess/elementwise_full_dom/test_power.h
@@ -20,7 +20,7 @@ const char *test_wsum_hess_power(void)
     power_node->forward(power_node, u_vals);
     jacobian_init(power_node);
     wsum_hess_init(power_node);
-    power_node->eval_wsum_hess(power_node, w);
+    eval_wsum_hess(power_node, w);
 
     /* Expected values on the diagonal: w_i * 6 * x_i */
     double expected_x[3] = {6.0 * 1.0, 6.0 * 4.0, 6.0 * 9.0};
diff --git a/tests/wsum_hess/elementwise_full_dom/test_trig.h b/tests/wsum_hess/elementwise_full_dom/test_trig.h
index 550c2727..6218109c 100644
--- a/tests/wsum_hess/elementwise_full_dom/test_trig.h
+++ b/tests/wsum_hess/elementwise_full_dom/test_trig.h
@@ -20,7 +20,7 @@ const char *test_wsum_hess_sin(void)
     sin_node->forward(sin_node, u_vals);
     jacobian_init(sin_node);
     wsum_hess_init(sin_node);
-    sin_node->eval_wsum_hess(sin_node, w);
+    eval_wsum_hess(sin_node, w);
 
     /* Expected values on the diagonal: -w_i * sin(x_i) */
     double expected_x[3] = {-1.0 * sin(1.0), -2.0 * sin(2.0), -3.0 * sin(3.0)};
@@ -46,7 +46,7 @@ const char *test_wsum_hess_cos(void)
     cos_node->forward(cos_node, u_vals);
     jacobian_init(cos_node);
     wsum_hess_init(cos_node);
-    cos_node->eval_wsum_hess(cos_node, w);
+    eval_wsum_hess(cos_node, w);
 
     /* Expected values on the diagonal: -w_i * cos(x_i) */
     double expected_x[3] = {-1.0 * cos(1.0), -2.0 * cos(2.0), -3.0 * cos(3.0)};
@@ -72,7 +72,7 @@ const char *test_wsum_hess_tan(void)
     tan_node->forward(tan_node, u_vals);
     jacobian_init(tan_node);
     wsum_hess_init(tan_node);
-    tan_node->eval_wsum_hess(tan_node, w);
+    eval_wsum_hess(tan_node, w);
 
     /* Expected values on the diagonal: w_i * 2 * sin(x_i) / cos^3(x_i) */
     double expected_x[3] = {1.0 * 2.0 * sin(1.0) / pow(cos(1.0), 3),
diff --git a/tests/wsum_hess/elementwise_full_dom/test_xexp.h b/tests/wsum_hess/elementwise_full_dom/test_xexp.h
index 55e80ef5..547014fb 100644
--- a/tests/wsum_hess/elementwise_full_dom/test_xexp.h
+++ b/tests/wsum_hess/elementwise_full_dom/test_xexp.h
@@ -20,7 +20,7 @@ const char *test_wsum_hess_xexp(void)
     xexp_node->forward(xexp_node, u_vals);
     jacobian_init(xexp_node);
     wsum_hess_init(xexp_node);
-    xexp_node->eval_wsum_hess(xexp_node, w);
+    eval_wsum_hess(xexp_node, w);
 
     /* Expected values on the diagonal: w_i * (2+x_i) * exp(x_i) */
     double expected_x[3] = {1.0 * 3.0 * exp(1.0), 2.0 * 4.0 * exp(2.0),
diff --git a/tests/wsum_hess/elementwise_restricted_dom/test_entr.h b/tests/wsum_hess/elementwise_restricted_dom/test_entr.h
index f7c83ca6..5858751d 100644
--- a/tests/wsum_hess/elementwise_restricted_dom/test_entr.h
+++ b/tests/wsum_hess/elementwise_restricted_dom/test_entr.h
@@ -20,7 +20,7 @@ const char *test_wsum_hess_entr(void)
     entr_node->forward(entr_node, u_vals);
     jacobian_init(entr_node);
     wsum_hess_init(entr_node);
-    entr_node->eval_wsum_hess(entr_node, w);
+    eval_wsum_hess(entr_node, w);
 
     /* Expected values on the diagonal: -w_i/x_i */
     double expected_x[3] = {-1.0, -1.0, -1.0};
diff --git a/tests/wsum_hess/elementwise_restricted_dom/test_log.h b/tests/wsum_hess/elementwise_restricted_dom/test_log.h
index c8acff59..c44b865c 100644
--- a/tests/wsum_hess/elementwise_restricted_dom/test_log.h
+++ b/tests/wsum_hess/elementwise_restricted_dom/test_log.h
@@ -32,7 +32,7 @@ const char *test_wsum_hess_log(void)
     log_node->forward(log_node, u_vals);
     jacobian_init(log_node);
     wsum_hess_init(log_node);
-    log_node->eval_wsum_hess(log_node, w);
+    eval_wsum_hess(log_node, w);
 
     /* Expected values on the diagonal: -w_i/x_i^2 */
     double expected_x[3] = {-1.0, -0.5, -1.0 / 3.0};
diff --git a/tests/wsum_hess/other/test_prod.h b/tests/wsum_hess/other/test_prod.h
index 28184c18..0da9ea96 100644
--- a/tests/wsum_hess/other/test_prod.h
+++ b/tests/wsum_hess/other/test_prod.h
@@ -19,7 +19,7 @@ const char *test_wsum_hess_prod_no_zero(void)
     p->forward(p, u_vals);
     jacobian_init(p);
     wsum_hess_init(p);
-    p->eval_wsum_hess(p, &w);
+    eval_wsum_hess(p, &w);
 
     /* Row-major over dense 4x4 block */
     double expected_x[16] = {0.0, 12.0, 8.0, 6.0, 12.0, 0.0, 4.0, 3.0,
@@ -47,7 +47,7 @@ const char *test_wsum_hess_prod_one_zero(void)
     p->forward(p, u_vals);
     jacobian_init(p);
     wsum_hess_init(p);
-    p->eval_wsum_hess(p, &w);
+    eval_wsum_hess(p, &w);
 
     double expected_x[16];
     memset(expected_x, 0, sizeof(expected_x));
@@ -81,7 +81,7 @@ const char *test_wsum_hess_prod_two_zeros(void)
     p->forward(p, u_vals);
     jacobian_init(p);
     wsum_hess_init(p);
-    p->eval_wsum_hess(p, &w);
+    eval_wsum_hess(p, &w);
 
     double expected_x[16];
     memset(expected_x, 0, sizeof(expected_x));
@@ -110,7 +110,7 @@ const char *test_wsum_hess_prod_many_zeros(void)
     p->forward(p, u_vals);
     jacobian_init(p);
     wsum_hess_init(p);
-    p->eval_wsum_hess(p, &w);
+    eval_wsum_hess(p, &w);
 
     double expected_x[16];
     memset(expected_x, 0, sizeof(expected_x));
diff --git a/tests/wsum_hess/other/test_prod_axis_one.h b/tests/wsum_hess/other/test_prod_axis_one.h
index 21b810da..ad68487d 100644
--- a/tests/wsum_hess/other/test_prod_axis_one.h
+++ b/tests/wsum_hess/other/test_prod_axis_one.h
@@ -32,7 +32,7 @@ const char *test_wsum_hess_prod_axis_one_no_zeros(void)
     p->forward(p, u_vals);
     jacobian_init(p);
     wsum_hess_init(p);
-    p->eval_wsum_hess(p, w_vals);
+    eval_wsum_hess(p, w_vals);
 
     double expected_x[12] = {/* Var 1 (row 0, col 0): [5, 3] (excludes col 0) */
                              5.0, 3.0,
@@ -100,7 +100,7 @@ const char *test_wsum_hess_prod_axis_one_one_zero(void)
     p->forward(p, u_vals);
     jacobian_init(p);
     wsum_hess_init(p);
-    p->eval_wsum_hess(p, w_vals);
+    eval_wsum_hess(p, w_vals);
 
     double expected_x[18];
     memset(expected_x, 0, sizeof(expected_x));
@@ -204,7 +204,7 @@ const char *test_wsum_hess_prod_axis_one_mixed_zeros(void)
     p->forward(p, u_vals);
     jacobian_init(p);
     wsum_hess_init(p);
-    p->eval_wsum_hess(p, w_vals);
+    eval_wsum_hess(p, w_vals);
 
     double expected_x[30];
     memset(expected_x, 0, sizeof(expected_x));
@@ -355,7 +355,7 @@ const char *test_wsum_hess_prod_axis_one_2x2(void)
     p->forward(p, u_vals);
     jacobian_init(p);
     wsum_hess_init(p);
-    p->eval_wsum_hess(p, w_vals);
+    eval_wsum_hess(p, w_vals);
 
     /* Expected sparse structure (nnz = 4, each row has 1 nnz) */
     double expected_x[4] = {1.0,  /* Var 0 (excludes col 0) */
diff --git a/tests/wsum_hess/other/test_prod_axis_zero.h b/tests/wsum_hess/other/test_prod_axis_zero.h
index e057e47b..13b74ea4 100644
--- a/tests/wsum_hess/other/test_prod_axis_zero.h
+++ b/tests/wsum_hess/other/test_prod_axis_zero.h
@@ -34,7 +34,7 @@ const char *test_wsum_hess_prod_axis_zero_no_zeros(void)
     p->forward(p, u_vals);
     jacobian_init(p);
     wsum_hess_init(p);
-    p->eval_wsum_hess(p, w_vals);
+    eval_wsum_hess(p, w_vals);
 
     /* Block diagonal structure: 3 blocks of 2x2 = 6 nnz total
      * Each block has 4 entries (2x2 dense)
@@ -110,7 +110,7 @@ const char *test_wsum_hess_prod_axis_zero_mixed_zeros(void)
     p->forward(p, u_vals);
     jacobian_init(p);
     wsum_hess_init(p);
-    p->eval_wsum_hess(p, w_vals);
+    eval_wsum_hess(p, w_vals);
 
     /* Block 0: 5x5 all off-diagonal = 1.0, total = 25 entries (indices 0-24)
      * Block 1: 5x5 with 1 zero, total = 25 entries (indices 25-49)
@@ -213,7 +213,7 @@ const char *test_wsum_hess_prod_axis_zero_one_zero(void)
     p->forward(p, u_vals);
     jacobian_init(p);
     wsum_hess_init(p);
-    p->eval_wsum_hess(p, w_vals);
+    eval_wsum_hess(p, w_vals);
 
     /* Block 0 (no zeros): w[0]*f[0] = 1
      *   (0,1) = 1/(1*1) = 1, (1,0) = 1
diff --git a/tests/wsum_hess/other/test_quad_form.h b/tests/wsum_hess/other/test_quad_form.h
index fa8956f8..3c2104f6 100644
--- a/tests/wsum_hess/other/test_quad_form.h
+++ b/tests/wsum_hess/other/test_quad_form.h
@@ -32,7 +32,7 @@ const char *test_wsum_hess_quad_form(void)
     node->forward(node, u_vals);
     jacobian_init(node);
     wsum_hess_init(node);
-    node->eval_wsum_hess(node, &w);
+    eval_wsum_hess(node, &w);
 
     int expected_p[11] = {0, 0, 0, 0, 2, 5, 8, 10, 10, 10, 10};
     int expected_i[10] = {3, 4, 3, 4, 5, 4, 5, 6, 5, 6};
diff --git a/tests/wsum_hess/other/test_quad_form_dense.h b/tests/wsum_hess/other/test_quad_form_dense.h
index becaa4f8..3c99d2a0 100644
--- a/tests/wsum_hess/other/test_quad_form_dense.h
+++ b/tests/wsum_hess/other/test_quad_form_dense.h
@@ -28,9 +28,9 @@ const char *test_wsum_hess_quad_form_dense(void)
 
     jacobian_init(node);
     node->forward(node, u_vals);
-    node->eval_jacobian(node);
+    eval_jacobian(node);
     wsum_hess_init(node);
-    node->eval_wsum_hess(node, &w);
+    eval_wsum_hess(node, &w);
 
     /* forward value */
     mu_assert("dense quad_form value fail", fabs(node->value[0] - 57.0) < 1e-9);
@@ -120,8 +120,8 @@ const char *test_wsum_hess_quad_form_dense_param(void)
     jacobian_init(node);
     wsum_hess_init(node);
     node->forward(node, u_vals);
-    node->eval_jacobian(node);
-    node->eval_wsum_hess(node, &w);
+    eval_jacobian(node);
+    eval_wsum_hess(node, &w);
 
     /* parameter P1: check values and sparsity */
     int expected_jp[2] = {0, 3};
@@ -144,8 +144,8 @@ const char *test_wsum_hess_quad_form_dense_param(void)
     memcpy(param_P->value, P2, 9 * sizeof(double));
     expr_set_needs_refresh(node);
     node->forward(node, u_vals);
-    node->eval_jacobian(node);
-    node->eval_wsum_hess(node, &w);
+    eval_jacobian(node);
+    eval_wsum_hess(node, &w);
 
     /* parameter P2: only the values change, sparsity is unchanged */
     double grad2[3] = {8.0, 18.0, 30.0};