diff --git a/SparseDiffEngine b/SparseDiffEngine
index 7e7875d..4353001 160000
--- a/SparseDiffEngine
+++ b/SparseDiffEngine
@@ -1 +1 @@
-Subproject commit 7e7875df92f92b0ef2ecd431d77174768aa22e01
+Subproject commit 435300196a87485222b3fd028dc84702eb78765e
diff --git a/sparsediffpy/_bindings/atoms/asinh.h b/sparsediffpy/_bindings/atoms/asinh.h
index a09934e..09724c8 100644
--- a/sparsediffpy/_bindings/atoms/asinh.h
+++ b/sparsediffpy/_bindings/atoms/asinh.h
@@ -2,7 +2,7 @@
 #define ATOM_ASINH_H
 
 #include "common.h"
-#include "elementwise_univariate.h"
+#include "elementwise_full_dom.h"
 
 static PyObject *py_make_asinh(PyObject *self, PyObject *args)
 {
diff --git a/sparsediffpy/_bindings/atoms/atanh.h b/sparsediffpy/_bindings/atoms/atanh.h
index fd0568d..86519b9 100644
--- a/sparsediffpy/_bindings/atoms/atanh.h
+++ b/sparsediffpy/_bindings/atoms/atanh.h
@@ -2,7 +2,7 @@
 #define ATOM_ATANH_H
 
 #include "common.h"
-#include "elementwise_univariate.h"
+#include "elementwise_restricted_dom.h"
 
 static PyObject *py_make_atanh(PyObject *self, PyObject *args)
 {
diff --git a/sparsediffpy/_bindings/atoms/common.h b/sparsediffpy/_bindings/atoms/common.h
index 83101dd..72cffdb 100644
--- a/sparsediffpy/_bindings/atoms/common.h
+++ b/sparsediffpy/_bindings/atoms/common.h
@@ -6,7 +6,8 @@
 #include <numpy/arrayobject.h>
 
 #include "affine.h"
-#include "elementwise_univariate.h"
+#include "elementwise_full_dom.h"
+#include "elementwise_restricted_dom.h"
 #include "expr.h"
 
 #define EXPR_CAPSULE_NAME "DNLP_EXPR"
diff --git a/sparsediffpy/_bindings/atoms/const_scalar_mult.h b/sparsediffpy/_bindings/atoms/const_scalar_mult.h
index 1a83e5a..365ef10 100644
--- a/sparsediffpy/_bindings/atoms/const_scalar_mult.h
+++ b/sparsediffpy/_bindings/atoms/const_scalar_mult.h
@@ -22,9 +22,18 @@ static PyObject *py_make_const_scalar_mult(PyObject *self, PyObject *args)
         return NULL;
     }
 
-    expr *node = new_const_scalar_mult(a, child);
+    expr *a_node = new_parameter(1, 1, PARAM_FIXED, child->n_vars, &a);
+    if (!a_node)
+    {
+        PyErr_SetString(PyExc_RuntimeError,
+                        "failed to create parameter node for scalar");
+        return NULL;
+    }
+
+    expr *node = new_scalar_mult(a_node, child);
     if (!node)
     {
+        free_expr(a_node);
         PyErr_SetString(PyExc_RuntimeError,
                         "failed to create const_scalar_mult node");
         return NULL;
diff --git a/sparsediffpy/_bindings/atoms/const_vector_mult.h b/sparsediffpy/_bindings/atoms/const_vector_mult.h
index aa1f83e..8b29ae7 100644
--- a/sparsediffpy/_bindings/atoms/const_vector_mult.h
+++ b/sparsediffpy/_bindings/atoms/const_vector_mult.h
@@ -42,12 +42,22 @@ static PyObject *py_make_const_vector_mult(PyObject *self, PyObject *args)
 
     double *a_data = (double *) PyArray_DATA(a_array);
 
-    expr *node = new_const_vector_mult(a_data, child);
+    expr *a_node =
+        new_parameter(a_size, 1, PARAM_FIXED, child->n_vars, a_data);
 
     Py_DECREF(a_array);
 
+    if (!a_node)
+    {
+        PyErr_SetString(PyExc_RuntimeError,
+                        "failed to create parameter node for vector");
+        return NULL;
+    }
+
+    expr *node = new_vector_mult(a_node, child);
     if (!node)
     {
+        free_expr(a_node);
         PyErr_SetString(PyExc_RuntimeError,
                         "failed to create const_vector_mult node");
         return NULL;
diff --git a/sparsediffpy/_bindings/atoms/constant.h b/sparsediffpy/_bindings/atoms/constant.h
index d5fcba3..d820f09 100644
--- a/sparsediffpy/_bindings/atoms/constant.h
+++ b/sparsediffpy/_bindings/atoms/constant.h
@@ -19,8 +19,8 @@ static PyObject *py_make_constant(PyObject *self, PyObject *args)
         return NULL;
     }
 
-    expr *node =
-        new_constant(d1, d2, n_vars, (const double *) PyArray_DATA(values_array));
+    expr *node = new_parameter(d1, d2, PARAM_FIXED, n_vars,
+                               (const double *) PyArray_DATA(values_array));
     Py_DECREF(values_array);
 
     if (!node)
diff --git a/sparsediffpy/_bindings/atoms/cos.h b/sparsediffpy/_bindings/atoms/cos.h
index 3f798be..152400d 100644
--- a/sparsediffpy/_bindings/atoms/cos.h
+++ b/sparsediffpy/_bindings/atoms/cos.h
@@ -2,7 +2,7 @@
 #define ATOM_COS_H
 
 #include "common.h"
-#include "elementwise_univariate.h"
+#include "elementwise_full_dom.h"
 
 static PyObject *py_make_cos(PyObject *self, PyObject *args)
 {
diff --git a/sparsediffpy/_bindings/atoms/dense_matmul.h b/sparsediffpy/_bindings/atoms/dense_matmul.h
index d8aef8f..b3fab51 100644
--- a/sparsediffpy/_bindings/atoms/dense_matmul.h
+++ b/sparsediffpy/_bindings/atoms/dense_matmul.h
@@ -7,17 +7,20 @@
 /* Dense left matrix multiplication: A @ f(x) where A is a dense matrix.
  *
  * Python signature:
- *   make_dense_left_matmul(child, A_data_flat, m, n)
+ *   make_dense_left_matmul(param_or_none, child, A_data_flat, m, n)
  *
+ * - param_or_none: None for constant matrix, or a parameter capsule.
  * - child: the child expression capsule f(x).
  * - A_data_flat: contiguous row-major numpy float64 array of size m*n.
  * - m, n: dimensions of matrix A. */
 static PyObject *py_make_dense_left_matmul(PyObject *self, PyObject *args)
 {
+    PyObject *param_obj;
     PyObject *child_capsule;
     PyObject *data_obj;
     int m, n;
-    if (!PyArg_ParseTuple(args, "OOii", &child_capsule, &data_obj, &m, &n))
+    if (!PyArg_ParseTuple(args, "OOOii", &param_obj, &child_capsule,
+                          &data_obj, &m, &n))
     {
         return NULL;
     }
@@ -38,11 +41,38 @@ static PyObject *py_make_dense_left_matmul(PyObject *self, PyObject *args)
 
     double *A_data = (double *) PyArray_DATA(data_array);
 
-    expr *node = new_left_matmul_dense(child, m, n, A_data);
+    /* Build the parameter node: use provided capsule or create PARAM_FIXED */
+    expr *param_node = NULL;
+    if (param_obj == Py_None)
+    {
+        param_node =
+            new_parameter(m * n, 1, PARAM_FIXED, child->n_vars, A_data);
+        if (!param_node)
+        {
+            Py_DECREF(data_array);
+            PyErr_SetString(PyExc_RuntimeError,
+                            "failed to create parameter node for dense matrix");
+            return NULL;
+        }
+    }
+    else
+    {
+        param_node =
+            (expr *) PyCapsule_GetPointer(param_obj, EXPR_CAPSULE_NAME);
+        if (!param_node)
+        {
+            Py_DECREF(data_array);
+            PyErr_SetString(PyExc_ValueError, "invalid parameter capsule");
+            return NULL;
+        }
+    }
+
+    expr *node = new_left_matmul_dense(param_node, child, m, n, A_data);
     Py_DECREF(data_array);
 
     if (!node)
     {
+        if (param_obj == Py_None) free_expr(param_node);
         PyErr_SetString(PyExc_RuntimeError,
                         "failed to create dense_left_matmul node");
         return NULL;
@@ -54,17 +84,20 @@ static PyObject *py_make_dense_left_matmul(PyObject *self, PyObject *args)
 /* Dense right matrix multiplication: f(x) @ A where A is a dense matrix.
  *
  * Python signature:
- *   make_dense_right_matmul(child, A_data_flat, m, n)
+ *   make_dense_right_matmul(param_or_none, child, A_data_flat, m, n)
  *
+ * - param_or_none: None for constant matrix, or a parameter capsule.
  * - child: the child expression capsule f(x).
  * - A_data_flat: contiguous row-major numpy float64 array of size m*n.
  * - m, n: dimensions of matrix A. */
 static PyObject *py_make_dense_right_matmul(PyObject *self, PyObject *args)
 {
+    PyObject *param_obj;
     PyObject *child_capsule;
     PyObject *data_obj;
     int m, n;
-    if (!PyArg_ParseTuple(args, "OOii", &child_capsule, &data_obj, &m, &n))
+    if (!PyArg_ParseTuple(args, "OOOii", &param_obj, &child_capsule,
+                          &data_obj, &m, &n))
     {
         return NULL;
     }
@@ -85,11 +118,38 @@ static PyObject *py_make_dense_right_matmul(PyObject *self, PyObject *args)
 
     double *A_data = (double *) PyArray_DATA(data_array);
 
-    expr *node = new_right_matmul_dense(child, m, n, A_data);
+    /* Build the parameter node: use provided capsule or create PARAM_FIXED */
+    expr *param_node = NULL;
+    if (param_obj == Py_None)
+    {
+        param_node =
+            new_parameter(m * n, 1, PARAM_FIXED, child->n_vars, A_data);
+        if (!param_node)
+        {
+            Py_DECREF(data_array);
+            PyErr_SetString(PyExc_RuntimeError,
+                            "failed to create parameter node for dense matrix");
+            return NULL;
+        }
+    }
+    else
+    {
+        param_node =
+            (expr *) PyCapsule_GetPointer(param_obj, EXPR_CAPSULE_NAME);
+        if (!param_node)
+        {
+            Py_DECREF(data_array);
+            PyErr_SetString(PyExc_ValueError, "invalid parameter capsule");
+            return NULL;
+        }
+    }
+
+    expr *node = new_right_matmul_dense(param_node, child, m, n, A_data);
     Py_DECREF(data_array);
 
     if (!node)
     {
+        if (param_obj == Py_None) free_expr(param_node);
         PyErr_SetString(PyExc_RuntimeError,
                         "failed to create dense_right_matmul node");
         return NULL;
@@ -98,4 +158,4 @@ static PyObject *py_make_dense_right_matmul(PyObject *self, PyObject *args)
     return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
 }
 
-#endif /* ATOM_DENSE_MATMUL_H */
\ No newline at end of file
+#endif /* ATOM_DENSE_MATMUL_H */
diff --git a/sparsediffpy/_bindings/atoms/entr.h b/sparsediffpy/_bindings/atoms/entr.h
index 9a37668..899ca00 100644
--- a/sparsediffpy/_bindings/atoms/entr.h
+++ b/sparsediffpy/_bindings/atoms/entr.h
@@ -2,7 +2,7 @@
 #define ATOM_ENTR_H
 
 #include "common.h"
-#include "elementwise_univariate.h"
+#include "elementwise_restricted_dom.h"
 
 static PyObject *py_make_entr(PyObject *self, PyObject *args)
 {
diff --git a/sparsediffpy/_bindings/atoms/left_matmul.h b/sparsediffpy/_bindings/atoms/left_matmul.h
index 09fa985..8de5ebd 100644
--- a/sparsediffpy/_bindings/atoms/left_matmul.h
+++ b/sparsediffpy/_bindings/atoms/left_matmul.h
@@ -4,14 +4,24 @@
 #include "bivariate.h"
 #include "common.h"
 
-/* Left matrix multiplication: A @ f(x) where A is a constant matrix */
+/* Left matrix multiplication: A @ f(x) where A is a constant or parameter
+ * sparse matrix.
+ *
+ * Python signature:
+ *   make_sparse_left_matmul(param_or_none, child, data, indices, indptr, m, n)
+ *
+ * - param_or_none: None for constant matrix, or a parameter capsule.
+ * - child: the child expression capsule f(x).
+ * - data, indices, indptr: CSR arrays for matrix A.
+ * - m, n: dimensions of matrix A. */
 static PyObject *py_make_sparse_left_matmul(PyObject *self, PyObject *args)
 {
+    PyObject *param_obj;
     PyObject *child_capsule;
     PyObject *data_obj, *indices_obj, *indptr_obj;
     int m, n;
-    if (!PyArg_ParseTuple(args, "OOOOii", &child_capsule, &data_obj, &indices_obj,
-                          &indptr_obj, &m, &n))
+    if (!PyArg_ParseTuple(args, "OOOOOii", &param_obj, &child_capsule,
+                          &data_obj, &indices_obj, &indptr_obj, &m, &n))
     {
         return NULL;
     }
@@ -39,6 +49,37 @@ static PyObject *py_make_sparse_left_matmul(PyObject *self, PyObject *args)
     }
 
     int nnz = (int) PyArray_SIZE(data_array);
+
+    /* Build the parameter node: use provided capsule or create PARAM_FIXED */
+    expr *param_node = NULL;
+    if (param_obj == Py_None)
+    {
+        param_node = new_parameter(nnz, 1, PARAM_FIXED, child->n_vars,
+                                   (const double *) PyArray_DATA(data_array));
+        if (!param_node)
+        {
+            Py_DECREF(data_array);
+            Py_DECREF(indices_array);
+            Py_DECREF(indptr_array);
+            PyErr_SetString(PyExc_RuntimeError,
+                            "failed to create parameter node for matrix");
+            return NULL;
+        }
+    }
+    else
+    {
+        param_node =
+            (expr *) PyCapsule_GetPointer(param_obj, EXPR_CAPSULE_NAME);
+        if (!param_node)
+        {
+            Py_DECREF(data_array);
+            Py_DECREF(indices_array);
+            Py_DECREF(indptr_array);
+            PyErr_SetString(PyExc_ValueError, "invalid parameter capsule");
+            return NULL;
+        }
+    }
+
     CSR_Matrix *A = new_csr_matrix(m, n, nnz);
     memcpy(A->x, PyArray_DATA(data_array), nnz * sizeof(double));
     memcpy(A->i, PyArray_DATA(indices_array), nnz * sizeof(int));
@@ -48,11 +89,12 @@ static PyObject *py_make_sparse_left_matmul(PyObject *self, PyObject *args)
     Py_DECREF(indices_array);
     Py_DECREF(indptr_array);
 
-    expr *node = new_left_matmul(child, A);
+    expr *node = new_left_matmul(param_node, child, A);
     free_csr_matrix(A);
 
     if (!node)
     {
+        if (param_obj == Py_None) free_expr(param_node);
         PyErr_SetString(PyExc_RuntimeError, "failed to create left_matmul node");
         return NULL;
     }
diff --git a/sparsediffpy/_bindings/atoms/logistic.h b/sparsediffpy/_bindings/atoms/logistic.h
index e903661..86fe1c8 100644
--- a/sparsediffpy/_bindings/atoms/logistic.h
+++ b/sparsediffpy/_bindings/atoms/logistic.h
@@ -2,7 +2,7 @@
 #define ATOM_LOGISTIC_H
 
 #include "common.h"
-#include "elementwise_univariate.h"
+#include "elementwise_full_dom.h"
 
 static PyObject *py_make_logistic(PyObject *self, PyObject *args)
 {
diff --git a/sparsediffpy/_bindings/atoms/parameter.h b/sparsediffpy/_bindings/atoms/parameter.h
new file mode 100644
index 0000000..940bf59
--- /dev/null
+++ b/sparsediffpy/_bindings/atoms/parameter.h
@@ -0,0 +1,24 @@
+#ifndef ATOM_PARAMETER_H
+#define ATOM_PARAMETER_H
+
+#include "common.h"
+
+static PyObject *py_make_parameter(PyObject *self, PyObject *args)
+{
+    int d1, d2, param_id, n_vars;
+    if (!PyArg_ParseTuple(args, "iiii", &d1, &d2, &param_id, &n_vars))
+    {
+        return NULL;
+    }
+
+    expr *node = new_parameter(d1, d2, param_id, n_vars, NULL);
+    if (!node)
+    {
+        PyErr_SetString(PyExc_RuntimeError, "failed to create parameter node");
+        return NULL;
+    }
+    expr_retain(node); /* Capsule owns a reference */
+    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
+}
+
+#endif /* ATOM_PARAMETER_H */
diff --git a/sparsediffpy/_bindings/atoms/power.h b/sparsediffpy/_bindings/atoms/power.h
index da9b7bf..a9fe96b 100644
--- a/sparsediffpy/_bindings/atoms/power.h
+++ b/sparsediffpy/_bindings/atoms/power.h
@@ -2,7 +2,7 @@
 #define ATOM_POWER_H
 
 #include "common.h"
-#include "elementwise_univariate.h"
+#include "elementwise_full_dom.h"
 
 static PyObject *py_make_power(PyObject *self, PyObject *args)
 {
diff --git a/sparsediffpy/_bindings/atoms/right_matmul.h b/sparsediffpy/_bindings/atoms/right_matmul.h
index f808058..4338b0c 100644
--- a/sparsediffpy/_bindings/atoms/right_matmul.h
+++ b/sparsediffpy/_bindings/atoms/right_matmul.h
@@ -4,14 +4,24 @@
 #include "bivariate.h"
 #include "common.h"
 
-/* Right matrix multiplication: f(x) @ A where A is a constant matrix */
+/* Right matrix multiplication: f(x) @ A where A is a constant or parameter
+ * sparse matrix.
+ *
+ * Python signature:
+ *   make_sparse_right_matmul(param_or_none, child, data, indices, indptr, m, n)
+ *
+ * - param_or_none: None for constant matrix, or a parameter capsule.
+ * - child: the child expression capsule f(x).
+ * - data, indices, indptr: CSR arrays for matrix A.
+ * - m, n: dimensions of matrix A. */
 static PyObject *py_make_sparse_right_matmul(PyObject *self, PyObject *args)
 {
+    PyObject *param_obj;
     PyObject *child_capsule;
     PyObject *data_obj, *indices_obj, *indptr_obj;
     int m, n;
-    if (!PyArg_ParseTuple(args, "OOOOii", &child_capsule, &data_obj, &indices_obj,
-                          &indptr_obj, &m, &n))
+    if (!PyArg_ParseTuple(args, "OOOOOii", &param_obj, &child_capsule,
+                          &data_obj, &indices_obj, &indptr_obj, &m, &n))
     {
         return NULL;
     }
@@ -39,6 +49,37 @@ static PyObject *py_make_sparse_right_matmul(PyObject *self, PyObject *args)
     }
 
     int nnz = (int) PyArray_SIZE(data_array);
+
+    /* Build the parameter node: use provided capsule or create PARAM_FIXED */
+    expr *param_node = NULL;
+    if (param_obj == Py_None)
+    {
+        param_node = new_parameter(nnz, 1, PARAM_FIXED, child->n_vars,
+                                   (const double *) PyArray_DATA(data_array));
+        if (!param_node)
+        {
+            Py_DECREF(data_array);
+            Py_DECREF(indices_array);
+            Py_DECREF(indptr_array);
+            PyErr_SetString(PyExc_RuntimeError,
+                            "failed to create parameter node for matrix");
+            return NULL;
+        }
+    }
+    else
+    {
+        param_node =
+            (expr *) PyCapsule_GetPointer(param_obj, EXPR_CAPSULE_NAME);
+        if (!param_node)
+        {
+            Py_DECREF(data_array);
+            Py_DECREF(indices_array);
+            Py_DECREF(indptr_array);
+            PyErr_SetString(PyExc_ValueError, "invalid parameter capsule");
+            return NULL;
+        }
+    }
+
     CSR_Matrix *A = new_csr_matrix(m, n, nnz);
     memcpy(A->x, PyArray_DATA(data_array), nnz * sizeof(double));
     memcpy(A->i, PyArray_DATA(indices_array), nnz * sizeof(int));
@@ -48,12 +89,14 @@ static PyObject *py_make_sparse_right_matmul(PyObject *self, PyObject *args)
     Py_DECREF(indices_array);
     Py_DECREF(indptr_array);
 
-    expr *node = new_right_matmul(child, A);
+    expr *node = new_right_matmul(param_node, child, A);
     free_csr_matrix(A);
 
     if (!node)
     {
-        PyErr_SetString(PyExc_RuntimeError, "failed to create right_matmul node");
+        if (param_obj == Py_None) free_expr(param_node);
+        PyErr_SetString(PyExc_RuntimeError,
+                        "failed to create right_matmul node");
         return NULL;
     }
     expr_retain(node); /* Capsule owns a reference */
diff --git a/sparsediffpy/_bindings/atoms/scalar_mult.h b/sparsediffpy/_bindings/atoms/scalar_mult.h
new file mode 100644
index 0000000..133ca11
--- /dev/null
+++ b/sparsediffpy/_bindings/atoms/scalar_mult.h
@@ -0,0 +1,45 @@
+#ifndef ATOM_SCALAR_MULT_H
+#define ATOM_SCALAR_MULT_H
+
+#include "common.h"
+
+/* Parameter scalar multiplication: a * f(x) where a comes from a parameter node
+ */
+static PyObject *py_make_param_scalar_mult(PyObject *self, PyObject *args)
+{
+    PyObject *param_capsule;
+    PyObject *child_capsule;
+
+    if (!PyArg_ParseTuple(args, "OO", &param_capsule, &child_capsule))
+    {
+        return NULL;
+    }
+
+    expr *param_node =
+        (expr *) PyCapsule_GetPointer(param_capsule, EXPR_CAPSULE_NAME);
+    if (!param_node)
+    {
+        PyErr_SetString(PyExc_ValueError, "invalid parameter capsule");
+        return NULL;
+    }
+
+    expr *child =
+        (expr *) PyCapsule_GetPointer(child_capsule, EXPR_CAPSULE_NAME);
+    if (!child)
+    {
+        PyErr_SetString(PyExc_ValueError, "invalid child capsule");
+        return NULL;
+    }
+
+    expr *node = new_scalar_mult(param_node, child);
+    if (!node)
+    {
+        PyErr_SetString(PyExc_RuntimeError,
+                        "failed to create scalar_mult node");
+        return NULL;
+    }
+    expr_retain(node); /* Capsule owns a reference */
+    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
+}
+
+#endif /* ATOM_SCALAR_MULT_H */
diff --git a/sparsediffpy/_bindings/atoms/sin.h b/sparsediffpy/_bindings/atoms/sin.h
index 426d302..89ca790 100644
--- a/sparsediffpy/_bindings/atoms/sin.h
+++ b/sparsediffpy/_bindings/atoms/sin.h
@@ -2,7 +2,7 @@
 #define ATOM_SIN_H
 
 #include "common.h"
-#include "elementwise_univariate.h"
+#include "elementwise_full_dom.h"
 
 static PyObject *py_make_sin(PyObject *self, PyObject *args)
 {
diff --git a/sparsediffpy/_bindings/atoms/sinh.h b/sparsediffpy/_bindings/atoms/sinh.h
index be86cf5..d3a1abb 100644
--- a/sparsediffpy/_bindings/atoms/sinh.h
+++ b/sparsediffpy/_bindings/atoms/sinh.h
@@ -2,7 +2,7 @@
 #define ATOM_SINH_H
 
 #include "common.h"
-#include "elementwise_univariate.h"
+#include "elementwise_full_dom.h"
 
 static PyObject *py_make_sinh(PyObject *self, PyObject *args)
 {
diff --git a/sparsediffpy/_bindings/atoms/tan.h b/sparsediffpy/_bindings/atoms/tan.h
index 3ec83b7..837128f 100644
--- a/sparsediffpy/_bindings/atoms/tan.h
+++ b/sparsediffpy/_bindings/atoms/tan.h
@@ -2,7 +2,7 @@
 #define ATOM_TAN_H
 
 #include "common.h"
-#include "elementwise_univariate.h"
+#include "elementwise_restricted_dom.h"
 
 static PyObject *py_make_tan(PyObject *self, PyObject *args)
 {
diff --git a/sparsediffpy/_bindings/atoms/tanh.h b/sparsediffpy/_bindings/atoms/tanh.h
index e60d420..6013173 100644
--- a/sparsediffpy/_bindings/atoms/tanh.h
+++ b/sparsediffpy/_bindings/atoms/tanh.h
@@ -2,7 +2,7 @@
 #define ATOM_TANH_H
 
 #include "common.h"
-#include "elementwise_univariate.h"
+#include "elementwise_full_dom.h"
 
 static PyObject *py_make_tanh(PyObject *self, PyObject *args)
 {
diff --git a/sparsediffpy/_bindings/atoms/vector_mult.h b/sparsediffpy/_bindings/atoms/vector_mult.h
new file mode 100644
index 0000000..a34d831
--- /dev/null
+++ b/sparsediffpy/_bindings/atoms/vector_mult.h
@@ -0,0 +1,45 @@
+#ifndef ATOM_VECTOR_MULT_H
+#define ATOM_VECTOR_MULT_H
+
+#include "common.h"
+
+/* Parameter vector elementwise multiplication: a . f(x) where a comes from a
+ * parameter node */
+static PyObject *py_make_param_vector_mult(PyObject *self, PyObject *args)
+{
+    PyObject *param_capsule;
+    PyObject *child_capsule;
+
+    if (!PyArg_ParseTuple(args, "OO", &param_capsule, &child_capsule))
+    {
+        return NULL;
+    }
+
+    expr *param_node =
+        (expr *) PyCapsule_GetPointer(param_capsule, EXPR_CAPSULE_NAME);
+    if (!param_node)
+    {
+        PyErr_SetString(PyExc_ValueError, "invalid parameter capsule");
+        return NULL;
+    }
+
+    expr *child =
+        (expr *) PyCapsule_GetPointer(child_capsule, EXPR_CAPSULE_NAME);
+    if (!child)
+    {
+        PyErr_SetString(PyExc_ValueError, "invalid child capsule");
+        return NULL;
+    }
+
+    expr *node = new_vector_mult(param_node, child);
+    if (!node)
+    {
+        PyErr_SetString(PyExc_RuntimeError,
+                        "failed to create vector_mult node");
+        return NULL;
+    }
+    expr_retain(node); /* Capsule owns a reference */
+    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
+}
+
+#endif /* ATOM_VECTOR_MULT_H */
diff --git a/sparsediffpy/_bindings/atoms/xexp.h b/sparsediffpy/_bindings/atoms/xexp.h
index 1dda312..324a9b7 100644
--- a/sparsediffpy/_bindings/atoms/xexp.h
+++ b/sparsediffpy/_bindings/atoms/xexp.h
@@ -2,7 +2,7 @@
 #define ATOM_XEXP_H
 
 #include "common.h"
-#include "elementwise_univariate.h"
+#include "elementwise_full_dom.h"
 
 static PyObject *py_make_xexp(PyObject *self, PyObject *args)
 {
diff --git a/sparsediffpy/_bindings/bindings.c b/sparsediffpy/_bindings/bindings.c
index cbda8e0..79edfb9 100644
--- a/sparsediffpy/_bindings/bindings.c
+++ b/sparsediffpy/_bindings/bindings.c
@@ -25,6 +25,7 @@
 #include "atoms/matmul.h"
 #include "atoms/multiply.h"
 #include "atoms/neg.h"
+#include "atoms/parameter.h"
 #include "atoms/normal_cdf.h"
 #include "atoms/power.h"
 #include "atoms/prod.h"
@@ -38,6 +39,7 @@
 #include "atoms/rel_entr_vector_scalar.h"
 #include "atoms/reshape.h"
 #include "atoms/right_matmul.h"
+#include "atoms/scalar_mult.h"
 #include "atoms/sin.h"
 #include "atoms/sinh.h"
 #include "atoms/sum.h"
@@ -46,6 +48,7 @@
 #include "atoms/trace.h"
 #include "atoms/transpose.h"
 #include "atoms/variable.h"
+#include "atoms/vector_mult.h"
 #include "atoms/xexp.h"
 
 /* Include problem bindings */
@@ -58,6 +61,8 @@
 #include "problem/jacobian.h"
 #include "problem/make_problem.h"
 #include "problem/objective_forward.h"
+#include "problem/register_params.h"
+#include "problem/update_params.h"
 
 static int numpy_initialized = 0;
 
@@ -71,6 +76,8 @@ static int ensure_numpy(void)
 
 static PyMethodDef DNLPMethods[] = {
     {"make_variable", py_make_variable, METH_VARARGS, "Create variable node"},
+    {"make_parameter", py_make_parameter, METH_VARARGS,
+     "Create updatable parameter node"},
     {"make_constant", py_make_constant, METH_VARARGS, "Create constant node"},
     {"make_linear", py_make_linear, METH_VARARGS, "Create linear op node"},
     {"make_log", py_make_log, METH_VARARGS, "Create log node"},
@@ -94,6 +101,10 @@ static PyMethodDef DNLPMethods[] = {
      "Create constant scalar multiplication node (a * f(x))"},
     {"make_const_vector_mult", py_make_const_vector_mult, METH_VARARGS,
      "Create constant vector multiplication node (a ∘ f(x))"},
+    {"make_param_scalar_mult", py_make_param_scalar_mult, METH_VARARGS,
+     "Create parameter scalar multiplication node (a * f(x))"},
+    {"make_param_vector_mult", py_make_param_vector_mult, METH_VARARGS,
+     "Create parameter vector elementwise multiplication node (a ∘ f(x))"},
     {"make_power", py_make_power, METH_VARARGS, "Create power node"},
     {"make_prod", py_make_prod, METH_VARARGS, "Create prod node"},
     {"make_prod_axis_zero", py_make_prod_axis_zero, METH_VARARGS,
@@ -170,6 +181,10 @@ static PyMethodDef DNLPMethods[] = {
      "Get Hessian sparsity in COO format (lower triangular)"},
     {"problem_eval_hessian_vals_coo", py_problem_eval_hessian_vals_coo, METH_VARARGS,
      "Evaluate Hessian and return COO values array"},
+    {"problem_register_params", py_problem_register_params, METH_VARARGS,
+     "Register parameter nodes with a problem"},
+    {"problem_update_params", py_problem_update_params, METH_VARARGS,
+     "Update parameter values from theta array"},
     {NULL, NULL, 0, NULL}};
 
 static struct PyModuleDef sparsediffpy_module = {
diff --git a/sparsediffpy/_bindings/problem/register_params.h b/sparsediffpy/_bindings/problem/register_params.h
new file mode 100644
index 0000000..346eb33
--- /dev/null
+++ b/sparsediffpy/_bindings/problem/register_params.h
@@ -0,0 +1,61 @@
+#ifndef PROBLEM_REGISTER_PARAMS_H
+#define PROBLEM_REGISTER_PARAMS_H
+
+#include "common.h"
+
+static PyObject *py_problem_register_params(PyObject *self, PyObject *args)
+{
+    PyObject *problem_capsule;
+    PyObject *param_list;
+
+    if (!PyArg_ParseTuple(args, "OO", &problem_capsule, &param_list))
+    {
+        return NULL;
+    }
+
+    problem *prob = (problem *) PyCapsule_GetPointer(problem_capsule,
+                                                     PROBLEM_CAPSULE_NAME);
+    if (!prob)
+    {
+        PyErr_SetString(PyExc_ValueError, "invalid problem capsule");
+        return NULL;
+    }
+
+    if (!PyList_Check(param_list))
+    {
+        PyErr_SetString(PyExc_TypeError,
+                        "second argument must be a list of parameter capsules");
+        return NULL;
+    }
+
+    Py_ssize_t n = PyList_Size(param_list);
+    expr **param_nodes = (expr **) malloc(n * sizeof(expr *));
+    if (!param_nodes)
+    {
+        PyErr_SetString(PyExc_MemoryError,
+                        "failed to allocate parameter node array");
+        return NULL;
+    }
+
+    for (Py_ssize_t i = 0; i < n; i++)
+    {
+        PyObject *capsule = PyList_GetItem(param_list, i);
+        expr *node =
+            (expr *) PyCapsule_GetPointer(capsule, EXPR_CAPSULE_NAME);
+        if (!node)
+        {
+            free(param_nodes);
+            PyErr_SetString(PyExc_ValueError,
+                            "invalid parameter capsule in list");
+            return NULL;
+        }
+        param_nodes[i] = node;
+    }
+
+    problem_register_params(prob, param_nodes, (int) n);
+    free(param_nodes);
+
+    Py_RETURN_NONE;
+}
+
+#endif /* PROBLEM_REGISTER_PARAMS_H */
diff --git a/sparsediffpy/_bindings/problem/update_params.h b/sparsediffpy/_bindings/problem/update_params.h
new file mode 100644
index 0000000..547de0c
--- /dev/null
+++ b/sparsediffpy/_bindings/problem/update_params.h
@@ -0,0 +1,37 @@
+#ifndef PROBLEM_UPDATE_PARAMS_H
+#define PROBLEM_UPDATE_PARAMS_H
+
+#include "common.h"
+
+static PyObject *py_problem_update_params(PyObject *self, PyObject *args)
+{
+    PyObject *problem_capsule;
+    PyObject *theta_obj;
+
+    if (!PyArg_ParseTuple(args, "OO", &problem_capsule, &theta_obj))
+    {
+        return NULL;
+    }
+
+    problem *prob = (problem *) PyCapsule_GetPointer(problem_capsule,
+                                                     PROBLEM_CAPSULE_NAME);
+    if (!prob)
+    {
+        PyErr_SetString(PyExc_ValueError, "invalid problem capsule");
+        return NULL;
+    }
+
+    PyArrayObject *theta_array = (PyArrayObject *) PyArray_FROM_OTF(
+        theta_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
+    if (!theta_array)
+    {
+        return NULL;
+    }
+
+    problem_update_params(prob, (const double *) PyArray_DATA(theta_array));
+    Py_DECREF(theta_array);
+
+    Py_RETURN_NONE;
+}
+
+#endif /* PROBLEM_UPDATE_PARAMS_H */