diff --git a/cmake/fix_stubgen.py b/cmake/fix_stubgen.py
new file mode 100755
index 00000000..1554b722
--- /dev/null
+++ b/cmake/fix_stubgen.py
@@ -0,0 +1,43 @@
+#!/usr/bin/env python3
+
+"""
+Fixes errors in .pyi files generated by stubgen.
+"""
+
+import re
+import os
+import sys
+
+
+def main():
+    # sys.argv[1] should be a directory in which to search for .pyi files
+    filenames = [
+        os.path.join(dp, f)
+        for dp, dn, fn in os.walk(sys.argv[1])
+        for f in fn
+        if f.endswith(".pyi")
+    ]
+    for filename in filenames:
+        with open(filename) as f:
+            content = f.read()
+
+        # Fix errors
+        content = content.replace(
+            "<ExpressionType.NONE: 0>", "ExpressionType.NONE"
+        ).replace("<SolverExitCondition.SUCCESS: 0>", "SolverExitCondition.SUCCESS")
+
+        # Convert parameter names from camel case to snake case
+        new_content = ""
+        extract_location = 0
+        for match in re.finditer(r"(?<=Parameter ``)(.*?)(?=``:)", content):
+            new_content += content[extract_location : match.start()]
+            new_content += re.sub(r"(?<!^)(?=[A-Z])", "_", match.group()).lower()
+            extract_location = match.end()
+        content = new_content + content[extract_location:]
+
+        with open(filename, mode="w") as f:
+            f.write(content)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/cmake/modules/Pybind11Stubgen.cmake b/cmake/modules/Pybind11Stubgen.cmake
index 8813806c..ca4cae5a 100644
--- a/cmake/modules/Pybind11Stubgen.cmake
+++ b/cmake/modules/Pybind11Stubgen.cmake
@@ -8,6 +8,10 @@ function(pybind11_stubgen target)
             --print-invalid-expressions-as-is --exit-code
             $<TARGET_FILE_BASE_NAME:${target}> -o
             $<TARGET_FILE_DIR:${target}>-stubs
+        COMMAND
+            ${Python3_EXECUTABLE}
+            ${CMAKE_CURRENT_SOURCE_DIR}/cmake/fix_stubgen.py
+            $<TARGET_FILE_DIR:${target}>-stubs
         WORKING_DIRECTORY $<TARGET_FILE_DIR:${target}>
         USES_TERMINAL
     )
diff --git a/include/sleipnir/autodiff/VariableBlock.hpp b/include/sleipnir/autodiff/VariableBlock.hpp
index 44c48fb6..3e06601a 100644
--- a/include/sleipnir/autodiff/VariableBlock.hpp
+++ b/include/sleipnir/autodiff/VariableBlock.hpp
@@ -26,6 +26,8 @@ class VariableBlock {
 
   /**
    * Assigns a VariableBlock to the block.
+   *
+   * @param values VariableBlock of values.
    */
   VariableBlock<Mat>& operator=(const VariableBlock<Mat>& values) {
     if (this == &values) {
@@ -56,6 +58,8 @@ class VariableBlock {
 
   /**
    * Assigns a VariableBlock to the block.
+   *
+   * @param values VariableBlock of values.
    */
   VariableBlock<Mat>& operator=(VariableBlock<Mat>&& values) {
     if (this == &values) {
@@ -124,6 +128,8 @@ class VariableBlock {
    * Assigns a double to the block.
    *
    * This only works for blocks with one row and one column.
+   *
+   * @param value Value to assign.
    */
   VariableBlock<Mat>& SetValue(double value) {
     Assert(Rows() == 1 && Cols() == 1);
@@ -135,6 +141,8 @@ class VariableBlock {
 
   /**
    * Assigns an Eigen matrix to the block.
+   *
+   * @param values Eigen matrix of values to assign.
    */
   template <typename Derived>
   VariableBlock<Mat>& operator=(const Eigen::MatrixBase<Derived>& values) {
@@ -152,6 +160,8 @@ class VariableBlock {
 
   /**
    * Sets block's internal values.
+   *
+   * @param values Eigen matrix of values.
    */
   template <typename Derived>
     requires std::same_as<typename Derived::Scalar, double>
@@ -170,6 +180,8 @@ class VariableBlock {
 
   /**
    * Assigns a VariableMatrix to the block.
+   *
+   * @param values VariableMatrix of values.
    */
   VariableBlock<Mat>& operator=(const Mat& values) {
     Assert(Rows() == values.Rows());
@@ -185,6 +197,8 @@ class VariableBlock {
 
   /**
    * Assigns a VariableMatrix to the block.
+   *
+   * @param values VariableMatrix of values.
    */
   VariableBlock<Mat>& operator=(Mat&& values) {
     Assert(Rows() == values.Rows());
diff --git a/include/sleipnir/autodiff/VariableMatrix.hpp b/include/sleipnir/autodiff/VariableMatrix.hpp
index 68c71a5a..ecd7a82d 100644
--- a/include/sleipnir/autodiff/VariableMatrix.hpp
+++ b/include/sleipnir/autodiff/VariableMatrix.hpp
@@ -137,6 +137,8 @@ class SLEIPNIR_DLLEXPORT VariableMatrix {
 
   /**
    * Constructs a VariableMatrix from an Eigen matrix.
+   *
+   * @param values Eigen matrix of values.
    */
   template <typename Derived>
   VariableMatrix(const Eigen::MatrixBase<Derived>& values)  // NOLINT
@@ -152,6 +154,8 @@ class SLEIPNIR_DLLEXPORT VariableMatrix {
 
   /**
    * Constructs a VariableMatrix from an Eigen diagonal matrix.
+   *
+   * @param values Diagonal matrix of values.
    */
   template <typename Derived>
   VariableMatrix(const Eigen::DiagonalBase<Derived>& values)  // NOLINT
@@ -171,6 +175,8 @@ class SLEIPNIR_DLLEXPORT VariableMatrix {
 
   /**
    * Assigns an Eigen matrix to a VariableMatrix.
+   *
+   * @param values Eigen matrix of values.
    */
   template <typename Derived>
   VariableMatrix& operator=(const Eigen::MatrixBase<Derived>& values) {
@@ -188,6 +194,8 @@ class SLEIPNIR_DLLEXPORT VariableMatrix {
 
   /**
    * Sets the VariableMatrix's internal values.
+   *
+   * @param values Eigen matrix of values.
    */
   template <typename Derived>
     requires std::same_as<typename Derived::Scalar, double>
@@ -206,6 +214,8 @@ class SLEIPNIR_DLLEXPORT VariableMatrix {
 
   /**
    * Constructs a scalar VariableMatrix from a Variable.
+   *
+   * @param variable Variable.
    */
   VariableMatrix(const Variable& variable)  // NOLINT
       : m_rows{1}, m_cols{1} {
@@ -214,6 +224,8 @@ class SLEIPNIR_DLLEXPORT VariableMatrix {
 
   /**
    * Constructs a scalar VariableMatrix from a Variable.
+   *
+   * @param variable Variable.
    */
   VariableMatrix(Variable&& variable) : m_rows{1}, m_cols{1} {  // NOLINT
     m_storage.emplace_back(std::move(variable));
@@ -221,6 +233,8 @@ class SLEIPNIR_DLLEXPORT VariableMatrix {
 
   /**
    * Constructs a VariableMatrix from a VariableBlock.
+   *
+   * @param values VariableBlock of values.
    */
   VariableMatrix(const VariableBlock<VariableMatrix>& values)  // NOLINT
       : m_rows{values.Rows()}, m_cols{values.Cols()} {
@@ -233,6 +247,8 @@ class SLEIPNIR_DLLEXPORT VariableMatrix {
 
   /**
    * Constructs a VariableMatrix from a VariableBlock.
+   *
+   * @param values VariableBlock of values.
    */
   VariableMatrix(const VariableBlock<const VariableMatrix>& values)  // NOLINT
       : m_rows{values.Rows()}, m_cols{values.Cols()} {
diff --git a/jormungandr/cpp/Binders.hpp b/jormungandr/cpp/Binders.hpp
index 390637eb..4b9802f5 100644
--- a/jormungandr/cpp/Binders.hpp
+++ b/jormungandr/cpp/Binders.hpp
@@ -3,23 +3,41 @@
 #pragma once
 
 #include <pybind11/pybind11.h>
+#include <sleipnir/autodiff/ExpressionType.hpp>
+#include <sleipnir/autodiff/Gradient.hpp>
+#include <sleipnir/autodiff/Hessian.hpp>
+#include <sleipnir/autodiff/Jacobian.hpp>
+#include <sleipnir/autodiff/Variable.hpp>
+#include <sleipnir/autodiff/VariableBlock.hpp>
+#include <sleipnir/autodiff/VariableMatrix.hpp>
+#include <sleipnir/optimization/Constraints.hpp>
+#include <sleipnir/optimization/OptimizationProblem.hpp>
+#include <sleipnir/optimization/SolverExitCondition.hpp>
+#include <sleipnir/optimization/SolverIterationInfo.hpp>
+#include <sleipnir/optimization/SolverStatus.hpp>
 
 namespace py = pybind11;
 
 namespace sleipnir {
 
-void BindExpressionType(py::module_& autodiff);
-void BindGradient(py::module_& autodiff);
-void BindHessian(py::module_& autodiff);
-void BindJacobian(py::module_& autodiff);
-void BindVariable(py::module_& autodiff);
-void BindVariableBlock(py::module_& autodiff);
-void BindVariableMatrix(py::module_& autodiff);
-
-void BindConstraints(py::module_& optimization);
-void BindOptimizationProblem(py::module_& optimization);
-void BindSolverExitCondition(py::module_& optimization);
-void BindSolverIterationInfo(py::module_& optimization);
-void BindSolverStatus(py::module_& optimization);
+void BindExpressionType(py::enum_<ExpressionType>& e);
+
+void BindVariable(py::module_& autodiff, py::class_<Variable>& cls);
+void BindVariableMatrix(py::module_& autodiff, py::class_<VariableMatrix>& cls);
+void BindVariableBlock(py::module_& autodiff,
+                       py::class_<VariableBlock<VariableMatrix>>& cls);
+
+void BindGradient(py::class_<Gradient>& cls);
+void BindHessian(py::class_<Hessian>& cls);
+void BindJacobian(py::class_<Jacobian>& cls);
+
+void BindEqualityConstraints(py::class_<EqualityConstraints>& cls);
+void BindInequalityConstraints(py::class_<InequalityConstraints>& cls);
+
+void BindSolverExitCondition(py::enum_<SolverExitCondition>& e);
+void BindSolverIterationInfo(py::class_<SolverIterationInfo>& cls);
+void BindSolverStatus(py::class_<SolverStatus>& cls);
+
+void BindOptimizationProblem(py::class_<OptimizationProblem>& cls);
 
 }  // namespace sleipnir
diff --git a/jormungandr/cpp/Docstrings.hpp b/jormungandr/cpp/Docstrings.hpp
index dedd62f0..9b02d663 100644
--- a/jormungandr/cpp/Docstrings.hpp
+++ b/jormungandr/cpp/Docstrings.hpp
@@ -1085,9 +1085,16 @@ static const char *__doc_sleipnir_VariableBlock_Rows = R"doc(Returns number of r
 static const char *__doc_sleipnir_VariableBlock_SetValue =
 R"doc(Assigns a double to the block.
 
-This only works for blocks with one row and one column.)doc";
+This only works for blocks with one row and one column.
+
+Parameter ``value``:
+    Value to assign.)doc";
+
+static const char *__doc_sleipnir_VariableBlock_SetValue_2 =
+R"doc(Sets block's internal values.
 
-static const char *__doc_sleipnir_VariableBlock_SetValue_2 = R"doc(Sets block's internal values.)doc";
+Parameter ``values``:
+    Eigen matrix of values.)doc";
 
 static const char *__doc_sleipnir_VariableBlock_T = R"doc(Returns the transpose of the variable matrix.)doc";
 
@@ -1194,20 +1201,40 @@ static const char *__doc_sleipnir_VariableBlock_m_mat = R"doc()doc";
 
 static const char *__doc_sleipnir_VariableBlock_m_rowOffset = R"doc()doc";
 
-static const char *__doc_sleipnir_VariableBlock_operator_assign = R"doc(Assigns a VariableBlock to the block.)doc";
+static const char *__doc_sleipnir_VariableBlock_operator_assign =
+R"doc(Assigns a VariableBlock to the block.
+
+Parameter ``values``:
+    VariableBlock of values.)doc";
+
+static const char *__doc_sleipnir_VariableBlock_operator_assign_2 =
+R"doc(Assigns a VariableBlock to the block.
 
-static const char *__doc_sleipnir_VariableBlock_operator_assign_2 = R"doc(Assigns a VariableBlock to the block.)doc";
+Parameter ``values``:
+    VariableBlock of values.)doc";
 
 static const char *__doc_sleipnir_VariableBlock_operator_assign_3 =
 R"doc(Assigns a double to the block.
 
 This only works for blocks with one row and one column.)doc";
 
-static const char *__doc_sleipnir_VariableBlock_operator_assign_4 = R"doc(Assigns an Eigen matrix to the block.)doc";
+static const char *__doc_sleipnir_VariableBlock_operator_assign_4 =
+R"doc(Assigns an Eigen matrix to the block.
+
+Parameter ``values``:
+    Eigen matrix of values to assign.)doc";
+
+static const char *__doc_sleipnir_VariableBlock_operator_assign_5 =
+R"doc(Assigns a VariableMatrix to the block.
 
-static const char *__doc_sleipnir_VariableBlock_operator_assign_5 = R"doc(Assigns a VariableMatrix to the block.)doc";
+Parameter ``values``:
+    VariableMatrix of values.)doc";
+
+static const char *__doc_sleipnir_VariableBlock_operator_assign_6 =
+R"doc(Assigns a VariableMatrix to the block.
 
-static const char *__doc_sleipnir_VariableBlock_operator_assign_6 = R"doc(Assigns a VariableMatrix to the block.)doc";
+Parameter ``values``:
+    VariableMatrix of values.)doc";
 
 static const char *__doc_sleipnir_VariableBlock_operator_call =
 R"doc(Returns a scalar subblock at the given row and column.
@@ -1373,7 +1400,11 @@ Parameter ``offset``:
 Parameter ``length``:
     The length of the segment.)doc";
 
-static const char *__doc_sleipnir_VariableMatrix_SetValue = R"doc(Sets the VariableMatrix's internal values.)doc";
+static const char *__doc_sleipnir_VariableMatrix_SetValue =
+R"doc(Sets the VariableMatrix's internal values.
+
+Parameter ``values``:
+    Eigen matrix of values.)doc";
 
 static const char *__doc_sleipnir_VariableMatrix_T = R"doc(Returns the transpose of the variable matrix.)doc";
 
@@ -1433,17 +1464,41 @@ This overload is for Python bindings only.
 Parameter ``list``:
     The nested list of Variables.)doc";
 
-static const char *__doc_sleipnir_VariableMatrix_VariableMatrix_7 = R"doc(Constructs a VariableMatrix from an Eigen matrix.)doc";
+static const char *__doc_sleipnir_VariableMatrix_VariableMatrix_7 =
+R"doc(Constructs a VariableMatrix from an Eigen matrix.
+
+Parameter ``values``:
+    Eigen matrix of values.)doc";
+
+static const char *__doc_sleipnir_VariableMatrix_VariableMatrix_8 =
+R"doc(Constructs a VariableMatrix from an Eigen diagonal matrix.
 
-static const char *__doc_sleipnir_VariableMatrix_VariableMatrix_8 = R"doc(Constructs a VariableMatrix from an Eigen diagonal matrix.)doc";
+Parameter ``values``:
+    Diagonal matrix of values.)doc";
+
+static const char *__doc_sleipnir_VariableMatrix_VariableMatrix_9 =
+R"doc(Constructs a scalar VariableMatrix from a Variable.
 
-static const char *__doc_sleipnir_VariableMatrix_VariableMatrix_9 = R"doc(Constructs a scalar VariableMatrix from a Variable.)doc";
+Parameter ``variable``:
+    Variable.)doc";
 
-static const char *__doc_sleipnir_VariableMatrix_VariableMatrix_10 = R"doc(Constructs a scalar VariableMatrix from a Variable.)doc";
+static const char *__doc_sleipnir_VariableMatrix_VariableMatrix_10 =
+R"doc(Constructs a scalar VariableMatrix from a Variable.
 
-static const char *__doc_sleipnir_VariableMatrix_VariableMatrix_11 = R"doc(Constructs a VariableMatrix from a VariableBlock.)doc";
+Parameter ``variable``:
+    Variable.)doc";
 
-static const char *__doc_sleipnir_VariableMatrix_VariableMatrix_12 = R"doc(Constructs a VariableMatrix from a VariableBlock.)doc";
+static const char *__doc_sleipnir_VariableMatrix_VariableMatrix_11 =
+R"doc(Constructs a VariableMatrix from a VariableBlock.
+
+Parameter ``values``:
+    VariableBlock of values.)doc";
+
+static const char *__doc_sleipnir_VariableMatrix_VariableMatrix_12 =
+R"doc(Constructs a VariableMatrix from a VariableBlock.
+
+Parameter ``values``:
+    VariableBlock of values.)doc";
 
 static const char *__doc_sleipnir_VariableMatrix_VariableMatrix_13 =
 R"doc(Constructs a column vector wrapper around a Variable array.
@@ -1528,7 +1583,11 @@ static const char *__doc_sleipnir_VariableMatrix_m_storage = R"doc()doc";
 
 static const char *__doc_sleipnir_VariableMatrix_operator_Variable = R"doc(Implicit conversion operator from 1x1 VariableMatrix to Variable.)doc";
 
-static const char *__doc_sleipnir_VariableMatrix_operator_assign = R"doc(Assigns an Eigen matrix to a VariableMatrix.)doc";
+static const char *__doc_sleipnir_VariableMatrix_operator_assign =
+R"doc(Assigns an Eigen matrix to a VariableMatrix.
+
+Parameter ``values``:
+    Eigen matrix of values.)doc";
 
 static const char *__doc_sleipnir_VariableMatrix_operator_call =
 R"doc(Returns a block pointing to the given row and column.
diff --git a/jormungandr/cpp/Main.cpp b/jormungandr/cpp/Main.cpp
index 12bbdc34..ed7b845c 100644
--- a/jormungandr/cpp/Main.cpp
+++ b/jormungandr/cpp/Main.cpp
@@ -1,8 +1,17 @@
 // Copyright (c) Sleipnir contributors
 
 #include <pybind11/pybind11.h>
+#include <sleipnir/autodiff/Variable.hpp>
+#include <sleipnir/autodiff/VariableBlock.hpp>
+#include <sleipnir/autodiff/VariableMatrix.hpp>
+#include <sleipnir/optimization/Constraints.hpp>
+#include <sleipnir/optimization/OptimizationProblem.hpp>
+#include <sleipnir/optimization/SolverExitCondition.hpp>
+#include <sleipnir/optimization/SolverIterationInfo.hpp>
+#include <sleipnir/optimization/SolverStatus.hpp>
 
 #include "Binders.hpp"
+#include "Docstrings.hpp"
 
 namespace py = pybind11;
 
@@ -14,20 +23,55 @@ PYBIND11_MODULE(_jormungandr, m) {
       "problem solver that uses the interior-point method.";
 
   py::module_ autodiff = m.def_submodule("autodiff");
-  BindExpressionType(autodiff);
-  BindGradient(autodiff);
-  BindHessian(autodiff);
-  BindJacobian(autodiff);
-  BindVariable(autodiff);
-  BindVariableMatrix(autodiff);
-  BindVariableBlock(autodiff);  // Must be bound after VariableMatrix
-
   py::module_ optimization = m.def_submodule("optimization");
-  BindConstraints(optimization);
-  BindSolverExitCondition(optimization);
-  BindSolverIterationInfo(optimization);
-  BindSolverStatus(optimization);
-  BindOptimizationProblem(optimization);
+
+  py::enum_<ExpressionType> expression_type{autodiff, "ExpressionType",
+                                            DOC(sleipnir, ExpressionType)};
+
+  py::class_<Variable> variable{autodiff, "Variable", DOC(sleipnir, Variable)};
+  py::class_<VariableMatrix> variable_matrix{autodiff, "VariableMatrix",
+                                             DOC(sleipnir, VariableMatrix)};
+  py::class_<VariableBlock<VariableMatrix>> variable_block{
+      autodiff, "VariableBlock", DOC(sleipnir, VariableBlock)};
+
+  py::class_<Gradient> gradient{autodiff, "Gradient", DOC(sleipnir, Gradient)};
+  py::class_<Hessian> hessian{autodiff, "Hessian", DOC(sleipnir, Hessian)};
+  py::class_<Jacobian> jacobian{autodiff, "Jacobian", DOC(sleipnir, Jacobian)};
+
+  py::class_<EqualityConstraints> equality_constraints{
+      optimization, "EqualityConstraints", DOC(sleipnir, EqualityConstraints)};
+  py::class_<InequalityConstraints> inequality_constraints{
+      optimization, "InequalityConstraints",
+      DOC(sleipnir_InequalityConstraints)};
+
+  py::enum_<SolverExitCondition> solver_exit_condition{
+      optimization, "SolverExitCondition", DOC(sleipnir, SolverExitCondition)};
+  py::class_<SolverIterationInfo> solver_iteration_info{
+      optimization, "SolverIterationInfo", DOC(sleipnir, SolverIterationInfo)};
+  py::class_<SolverStatus> solver_status{optimization, "SolverStatus",
+                                         DOC(sleipnir, SolverStatus)};
+
+  py::class_<OptimizationProblem> optimization_problem{
+      optimization, "OptimizationProblem", DOC(sleipnir, OptimizationProblem)};
+
+  BindExpressionType(expression_type);
+
+  BindVariable(autodiff, variable);
+  BindVariableMatrix(autodiff, variable_matrix);
+  BindVariableBlock(autodiff, variable_block);
+
+  BindGradient(gradient);
+  BindHessian(hessian);
+  BindJacobian(jacobian);
+
+  BindEqualityConstraints(equality_constraints);
+  BindInequalityConstraints(inequality_constraints);
+
+  BindSolverExitCondition(solver_exit_condition);
+  BindSolverIterationInfo(solver_iteration_info);
+  BindSolverStatus(solver_status);
+
+  BindOptimizationProblem(optimization_problem);
 }
 
 }  // namespace sleipnir
diff --git a/jormungandr/cpp/autodiff/BindExpressionType.cpp b/jormungandr/cpp/autodiff/BindExpressionType.cpp
index c049e14e..3d277979 100644
--- a/jormungandr/cpp/autodiff/BindExpressionType.cpp
+++ b/jormungandr/cpp/autodiff/BindExpressionType.cpp
@@ -9,20 +9,16 @@ namespace py = pybind11;
 
 namespace sleipnir {
 
-void BindExpressionType(py::module_& autodiff) {
-  py::enum_<ExpressionType> expressionType{autodiff, "ExpressionType",
-                                           DOC(sleipnir, ExpressionType)};
-  expressionType
-      .value("NONE", ExpressionType::kNone,
-             DOC(sleipnir, ExpressionType, kNone))
-      .value("CONSTANT", ExpressionType::kConstant,
-             DOC(sleipnir, ExpressionType, kConstant))
-      .value("LINEAR", ExpressionType::kLinear,
-             DOC(sleipnir, ExpressionType, kLinear))
-      .value("QUADRATIC", ExpressionType::kQuadratic,
-             DOC(sleipnir, ExpressionType, kQuadratic))
-      .value("NONLINEAR", ExpressionType::kNonlinear,
-             DOC(sleipnir, ExpressionType, kNonlinear));
+void BindExpressionType(py::enum_<ExpressionType>& e) {
+  e.value("NONE", ExpressionType::kNone, DOC(sleipnir, ExpressionType, kNone));
+  e.value("CONSTANT", ExpressionType::kConstant,
+          DOC(sleipnir, ExpressionType, kConstant));
+  e.value("LINEAR", ExpressionType::kLinear,
+          DOC(sleipnir, ExpressionType, kLinear));
+  e.value("QUADRATIC", ExpressionType::kQuadratic,
+          DOC(sleipnir, ExpressionType, kQuadratic));
+  e.value("NONLINEAR", ExpressionType::kNonlinear,
+          DOC(sleipnir, ExpressionType, kNonlinear));
 }
 
 }  // namespace sleipnir
diff --git a/jormungandr/cpp/autodiff/BindGradient.cpp b/jormungandr/cpp/autodiff/BindGradient.cpp
index 7b0f07a3..6a5af676 100644
--- a/jormungandr/cpp/autodiff/BindGradient.cpp
+++ b/jormungandr/cpp/autodiff/BindGradient.cpp
@@ -10,19 +10,19 @@ namespace py = pybind11;
 
 namespace sleipnir {
 
-void BindGradient(py::module_& autodiff) {
-  py::class_<Gradient> cls{autodiff, "Gradient", DOC(sleipnir, Gradient)};
-  cls.def(py::init<Variable, Variable>(), DOC(sleipnir, Gradient, Gradient))
-      .def(py::init<Variable, VariableMatrix>(),
-           DOC(sleipnir, Gradient, Gradient, 2))
-      .def("get", &Gradient::Get, DOC(sleipnir, Gradient, Get))
-      .def(
-          "value",
-          [](Gradient& self) {
-            return Eigen::SparseMatrix<double>{self.Value()};
-          },
-          DOC(sleipnir, Gradient, Value))
-      .def("update", &Gradient::Update, DOC(sleipnir, Gradient, Update));
+void BindGradient(py::class_<Gradient>& cls) {
+  using namespace pybind11::literals;
+
+  cls.def(py::init<Variable, Variable>(), "variable"_a, "wrt"_a,
+          DOC(sleipnir, Gradient, Gradient));
+  cls.def(py::init<Variable, VariableMatrix>(), "variable"_a, "wrt"_a,
+          DOC(sleipnir, Gradient, Gradient, 2));
+  cls.def("get", &Gradient::Get, DOC(sleipnir, Gradient, Get));
+  cls.def(
+      "value",
+      [](Gradient& self) { return Eigen::SparseMatrix<double>{self.Value()}; },
+      DOC(sleipnir, Gradient, Value));
+  cls.def("update", &Gradient::Update, DOC(sleipnir, Gradient, Update));
 }
 
 }  // namespace sleipnir
diff --git a/jormungandr/cpp/autodiff/BindHessian.cpp b/jormungandr/cpp/autodiff/BindHessian.cpp
index 4d7da4f2..7514cedb 100644
--- a/jormungandr/cpp/autodiff/BindHessian.cpp
+++ b/jormungandr/cpp/autodiff/BindHessian.cpp
@@ -10,12 +10,14 @@ namespace py = pybind11;
 
 namespace sleipnir {
 
-void BindHessian(py::module_& autodiff) {
-  py::class_<Hessian> cls{autodiff, "Hessian", DOC(sleipnir, Hessian)};
-  cls.def(py::init<Variable, VariableMatrix>(), DOC(sleipnir, Hessian, Hessian))
-      .def("get", &Hessian::Get, DOC(sleipnir, Hessian, Get))
-      .def("value", &Hessian::Value, DOC(sleipnir, Hessian, Value))
-      .def("update", &Hessian::Update, DOC(sleipnir, Hessian, Update));
+void BindHessian(py::class_<Hessian>& cls) {
+  using namespace pybind11::literals;
+
+  cls.def(py::init<Variable, VariableMatrix>(), "variable"_a, "wrt"_a,
+          DOC(sleipnir, Hessian, Hessian));
+  cls.def("get", &Hessian::Get, DOC(sleipnir, Hessian, Get));
+  cls.def("value", &Hessian::Value, DOC(sleipnir, Hessian, Value));
+  cls.def("update", &Hessian::Update, DOC(sleipnir, Hessian, Update));
 }
 
 }  // namespace sleipnir
diff --git a/jormungandr/cpp/autodiff/BindJacobian.cpp b/jormungandr/cpp/autodiff/BindJacobian.cpp
index 39b2e731..6ce7f782 100644
--- a/jormungandr/cpp/autodiff/BindJacobian.cpp
+++ b/jormungandr/cpp/autodiff/BindJacobian.cpp
@@ -10,13 +10,14 @@ namespace py = pybind11;
 
 namespace sleipnir {
 
-void BindJacobian(py::module_& autodiff) {
-  py::class_<Jacobian> cls{autodiff, "Jacobian", DOC(sleipnir, Jacobian)};
-  cls.def(py::init<VariableMatrix, VariableMatrix>(),
-          DOC(sleipnir, Jacobian, Jacobian))
-      .def("get", &Jacobian::Get, DOC(sleipnir, Jacobian, Get))
-      .def("value", &Jacobian::Value, DOC(sleipnir, Jacobian, Value))
-      .def("update", &Jacobian::Update, DOC(sleipnir, Jacobian, Update));
+void BindJacobian(py::class_<Jacobian>& cls) {
+  using namespace pybind11::literals;
+
+  cls.def(py::init<VariableMatrix, VariableMatrix>(), "variables"_a, "wrt"_a,
+          DOC(sleipnir, Jacobian, Jacobian));
+  cls.def("get", &Jacobian::Get, DOC(sleipnir, Jacobian, Get));
+  cls.def("value", &Jacobian::Value, DOC(sleipnir, Jacobian, Value));
+  cls.def("update", &Jacobian::Update, DOC(sleipnir, Jacobian, Update));
 }
 
 }  // namespace sleipnir
diff --git a/jormungandr/cpp/autodiff/BindVariable.cpp b/jormungandr/cpp/autodiff/BindVariable.cpp
index f64e3ca4..9093c5d8 100644
--- a/jormungandr/cpp/autodiff/BindVariable.cpp
+++ b/jormungandr/cpp/autodiff/BindVariable.cpp
@@ -22,186 +22,195 @@ PYBIND11_WARNING_DISABLE_CLANG("-Wself-assign-overloaded")
 
 namespace sleipnir {
 
-void BindVariable(py::module_& autodiff) {
-  py::class_<Variable> variable{autodiff, "Variable", DOC(sleipnir, Variable)};
-  variable.def(py::init<>(), DOC(sleipnir, Variable, Variable));
-  variable.def(py::init<double>(), DOC(sleipnir, Variable, Variable, 2));
-  variable.def(py::init<int>(), DOC(sleipnir, Variable, Variable, 3));
-  variable.def("set_value", py::overload_cast<double>(&Variable::SetValue),
-               DOC(sleipnir, Variable, SetValue));
-  variable.def("set_value", py::overload_cast<int>(&Variable::SetValue),
-               DOC(sleipnir, Variable, SetValue, 2));
-  variable.def(double() * py::self);
-  variable.def(py::self * double());
-  variable.def(py::self * py::self);
-  variable.def(py::self *= double(), DOC(sleipnir, Variable, operator, imul));
-  variable.def(py::self *= py::self, DOC(sleipnir, Variable, operator, imul));
-  variable.def(double() / py::self);
-  variable.def(py::self / double());
-  variable.def(py::self / py::self);
-  variable.def(py::self /= double(), DOC(sleipnir, Variable, operator, idiv));
-  variable.def(py::self /= py::self, DOC(sleipnir, Variable, operator, idiv));
-  variable.def(double() + py::self);
-  variable.def(py::self + double());
-  variable.def(py::self + py::self);
-  variable.def(py::self += double(), DOC(sleipnir, Variable, operator, iadd));
-  variable.def(py::self += py::self, DOC(sleipnir, Variable, operator, iadd));
-  variable.def(double() - py::self);
-  variable.def(py::self - double());
-  variable.def(py::self - py::self);
-  variable.def(py::self -= double(), DOC(sleipnir, Variable, operator, isub));
-  variable.def(py::self -= py::self, DOC(sleipnir, Variable, operator, isub));
-  variable.def(
+void BindVariable(py::module_& autodiff, py::class_<Variable>& cls) {
+  using namespace pybind11::literals;
+
+  cls.def(py::init<>(), DOC(sleipnir, Variable, Variable));
+  cls.def(py::init<double>(), "value"_a, DOC(sleipnir, Variable, Variable, 2));
+  cls.def(py::init<int>(), "value"_a, DOC(sleipnir, Variable, Variable, 3));
+  cls.def("set_value", py::overload_cast<double>(&Variable::SetValue),
+          "value"_a, DOC(sleipnir, Variable, SetValue));
+  cls.def("set_value", py::overload_cast<int>(&Variable::SetValue), "value"_a,
+          DOC(sleipnir, Variable, SetValue, 2));
+  cls.def(double() * py::self, "lhs"_a);
+  cls.def(py::self * double(), "rhs"_a);
+  cls.def(py::self * py::self, "rhs"_a);
+  cls.def(py::self *= double(), "rhs"_a,
+          DOC(sleipnir, Variable, operator, imul));
+  cls.def(py::self *= py::self, "rhs"_a,
+          DOC(sleipnir, Variable, operator, imul));
+  cls.def(double() / py::self, "lhs"_a);
+  cls.def(py::self / double(), "rhs"_a);
+  cls.def(py::self / py::self, "rhs"_a);
+  cls.def(py::self /= double(), "rhs"_a,
+          DOC(sleipnir, Variable, operator, idiv));
+  cls.def(py::self /= py::self, DOC(sleipnir, Variable, operator, idiv),
+          "rhs"_a);
+  cls.def(double() + py::self, "lhs"_a);
+  cls.def(py::self + double(), "rhs"_a);
+  cls.def(py::self + py::self, "rhs"_a);
+  cls.def(py::self += double(), "rhs"_a,
+          DOC(sleipnir, Variable, operator, iadd));
+  cls.def(py::self += py::self, "rhs"_a,
+          DOC(sleipnir, Variable, operator, iadd));
+  cls.def(double() - py::self, "lhs"_a);
+  cls.def(py::self - double(), "rhs"_a);
+  cls.def(py::self - py::self, "rhs"_a);
+  cls.def(py::self -= double(), "rhs"_a,
+          DOC(sleipnir, Variable, operator, isub));
+  cls.def(py::self -= py::self, "rhs"_a,
+          DOC(sleipnir, Variable, operator, isub));
+  cls.def(
       "__pow__",
       [](const Variable& self, int power) {
         return sleipnir::pow(self, power);
       },
-      py::is_operator());
-  variable.def(-py::self);
-  variable.def(+py::self);
-  variable.def("value", &Variable::Value, DOC(sleipnir, Variable, Value));
-  variable.def("type", &Variable::Type, DOC(sleipnir, Variable, Type));
-  variable.def("update", &Variable::Update, DOC(sleipnir, Variable, Update));
-  variable.def(py::self == py::self, DOC(sleipnir, operator, eq));
-  variable.def(py::self < py::self, DOC(sleipnir, operator, lt));
-  variable.def(py::self <= py::self, DOC(sleipnir, operator, le));
-  variable.def(py::self > py::self, DOC(sleipnir, operator, gt));
-  variable.def(py::self >= py::self, DOC(sleipnir, operator, ge));
-  variable.def(py::self == double(), DOC(sleipnir, operator, eq));
-  variable.def(py::self < double(), DOC(sleipnir, operator, lt));
-  variable.def(py::self <= double(), DOC(sleipnir, operator, le));
-  variable.def(py::self > double(), DOC(sleipnir, operator, gt));
-  variable.def(py::self >= double(), DOC(sleipnir, operator, ge));
-  variable.def(double() == py::self, DOC(sleipnir, operator, eq));
-  variable.def(double() < py::self, DOC(sleipnir, operator, lt));
-  variable.def(double() <= py::self, DOC(sleipnir, operator, le));
-  variable.def(double() > py::self, DOC(sleipnir, operator, gt));
-  variable.def(double() >= py::self, DOC(sleipnir, operator, ge));
+      py::is_operator(), "power"_a);
+  cls.def(-py::self);
+  cls.def(+py::self);
+  cls.def("value", &Variable::Value, DOC(sleipnir, Variable, Value));
+  cls.def("type", &Variable::Type, DOC(sleipnir, Variable, Type));
+  cls.def("update", &Variable::Update, DOC(sleipnir, Variable, Update));
+  cls.def(py::self == py::self, "rhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(py::self < py::self, "rhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(py::self <= py::self, "rhs"_a, DOC(sleipnir, operator, le));
+  cls.def(py::self > py::self, "rhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(py::self >= py::self, "rhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(py::self == double(), "rhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(py::self < double(), "rhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(py::self <= double(), "rhs"_a, DOC(sleipnir, operator, le));
+  cls.def(py::self > double(), "rhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(py::self >= double(), "rhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(double() == py::self, "lhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(double() < py::self, "lhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(double() <= py::self, "lhs"_a, DOC(sleipnir, operator, le));
+  cls.def(double() > py::self, "lhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(double() >= py::self, "lhs"_a, DOC(sleipnir, operator, ge));
 
   autodiff.def(
-      "abs", [](double x) { return sleipnir::abs(Variable{x}); },
+      "abs", [](double x) { return sleipnir::abs(Variable{x}); }, "x"_a,
       DOC(sleipnir, abs));
-  autodiff.def("abs", static_cast<Variable (*)(const Variable&)>(&abs),
+  autodiff.def("abs", static_cast<Variable (*)(const Variable&)>(&abs), "x"_a,
                DOC(sleipnir, abs));
   autodiff.def(
-      "acos", [](double x) { return sleipnir::acos(Variable{x}); },
+      "acos", [](double x) { return sleipnir::acos(Variable{x}); }, "x"_a,
       DOC(sleipnir, acos));
-  autodiff.def("acos", static_cast<Variable (*)(const Variable&)>(&acos),
+  autodiff.def("acos", static_cast<Variable (*)(const Variable&)>(&acos), "x"_a,
                DOC(sleipnir, acos));
   autodiff.def(
-      "asin", [](double x) { return sleipnir::asin(Variable{x}); },
+      "asin", [](double x) { return sleipnir::asin(Variable{x}); }, "x"_a,
       DOC(sleipnir, asin));
-  autodiff.def("asin", static_cast<Variable (*)(const Variable&)>(&asin),
+  autodiff.def("asin", static_cast<Variable (*)(const Variable&)>(&asin), "x"_a,
                DOC(sleipnir, asin));
   autodiff.def(
-      "atan", [](double x) { return sleipnir::atan(Variable{x}); },
+      "atan", [](double x) { return sleipnir::atan(Variable{x}); }, "x"_a,
       DOC(sleipnir, atan));
-  autodiff.def("atan", static_cast<Variable (*)(const Variable&)>(&atan),
+  autodiff.def("atan", static_cast<Variable (*)(const Variable&)>(&atan), "x"_a,
                DOC(sleipnir, atan));
   autodiff.def(
       "atan2",
-      [](double y, const Variable& x) { return sleipnir::atan2(y, x); },
-      DOC(sleipnir, atan2));
+      [](double y, const Variable& x) { return sleipnir::atan2(y, x); }, "y"_a,
+      "x"_a, DOC(sleipnir, atan2));
   autodiff.def(
       "atan2",
-      [](const Variable& y, double x) { return sleipnir::atan2(y, x); },
-      DOC(sleipnir, atan2));
+      [](const Variable& y, double x) { return sleipnir::atan2(y, x); }, "y"_a,
+      "x"_a, DOC(sleipnir, atan2));
   autodiff.def(
       "atan2",
       [](const Variable& y, const Variable& x) {
         return sleipnir::atan2(y, x);
       },
-      DOC(sleipnir, atan2));
+      "y"_a, "x"_a, DOC(sleipnir, atan2));
   autodiff.def(
-      "cos", [](double x) { return sleipnir::cos(Variable{x}); },
+      "cos", [](double x) { return sleipnir::cos(Variable{x}); }, "x"_a,
       DOC(sleipnir, cos));
-  autodiff.def("cos", static_cast<Variable (*)(const Variable&)>(&cos),
+  autodiff.def("cos", static_cast<Variable (*)(const Variable&)>(&cos), "x"_a,
                DOC(sleipnir, cos));
   autodiff.def(
-      "cosh", [](double x) { return sleipnir::cosh(Variable{x}); },
+      "cosh", [](double x) { return sleipnir::cosh(Variable{x}); }, "x"_a,
       DOC(sleipnir, cosh));
-  autodiff.def("cosh", static_cast<Variable (*)(const Variable&)>(&cosh),
+  autodiff.def("cosh", static_cast<Variable (*)(const Variable&)>(&cosh), "x"_a,
                DOC(sleipnir, cosh));
   autodiff.def(
-      "erf", [](double x) { return sleipnir::erf(Variable{x}); },
+      "erf", [](double x) { return sleipnir::erf(Variable{x}); }, "x"_a,
       DOC(sleipnir, erf));
-  autodiff.def("erf", static_cast<Variable (*)(const Variable&)>(&erf),
+  autodiff.def("erf", static_cast<Variable (*)(const Variable&)>(&erf), "x"_a,
                DOC(sleipnir, erf));
   autodiff.def(
-      "exp", [](double x) { return sleipnir::exp(Variable{x}); },
+      "exp", [](double x) { return sleipnir::exp(Variable{x}); }, "x"_a,
       DOC(sleipnir, exp));
-  autodiff.def("exp", static_cast<Variable (*)(const Variable&)>(&exp),
+  autodiff.def("exp", static_cast<Variable (*)(const Variable&)>(&exp), "x"_a,
                DOC(sleipnir, exp));
   autodiff.def(
       "hypot",
-      [](double x, const Variable& y) { return sleipnir::hypot(x, y); },
-      DOC(sleipnir, hypot));
+      [](double x, const Variable& y) { return sleipnir::hypot(x, y); }, "x"_a,
+      "y"_a, DOC(sleipnir, hypot));
   autodiff.def(
       "hypot",
-      [](const Variable& x, double y) { return sleipnir::hypot(x, y); },
-      DOC(sleipnir, hypot));
+      [](const Variable& x, double y) { return sleipnir::hypot(x, y); }, "x"_a,
+      "y"_a, DOC(sleipnir, hypot));
   autodiff.def(
       "hypot",
       static_cast<Variable (*)(const Variable&, const Variable&)>(&hypot),
-      DOC(sleipnir, hypot));
+      "x"_a, "y"_a, DOC(sleipnir, hypot));
   autodiff.def("hypot",
                static_cast<Variable (*)(const Variable&, const Variable&,
                                         const Variable&)>(&hypot),
-               DOC(sleipnir, hypot, 2));
+               "x"_a, "y"_a, "z"_a, DOC(sleipnir, hypot, 2));
   autodiff.def(
-      "log", [](double x) { return sleipnir::log(Variable{x}); },
+      "log", [](double x) { return sleipnir::log(Variable{x}); }, "x"_a,
       DOC(sleipnir, log));
-  autodiff.def("log", static_cast<Variable (*)(const Variable&)>(&log),
+  autodiff.def("log", static_cast<Variable (*)(const Variable&)>(&log), "x"_a,
                DOC(sleipnir, log));
   autodiff.def(
-      "log10", [](double x) { return sleipnir::log10(Variable{x}); },
+      "log10", [](double x) { return sleipnir::log10(Variable{x}); }, "x"_a,
       DOC(sleipnir, log10));
   autodiff.def("log10", static_cast<Variable (*)(const Variable&)>(&log10),
-               DOC(sleipnir, log10));
+               "x"_a, DOC(sleipnir, log10));
   autodiff.def(
       "pow",
       [](double base, const Variable& power) {
         return sleipnir::pow(base, power);
       },
-      DOC(sleipnir, pow));
+      "base"_a, "power"_a, DOC(sleipnir, pow));
   autodiff.def(
       "pow",
       [](const Variable& base, double power) {
         return sleipnir::pow(base, power);
       },
-      DOC(sleipnir, pow));
+      "base"_a, "power"_a, DOC(sleipnir, pow));
   autodiff.def(
       "pow", static_cast<Variable (*)(const Variable&, const Variable&)>(&pow),
-      DOC(sleipnir, pow));
+      "base"_a, "power"_a, DOC(sleipnir, pow));
   autodiff.def(
-      "sign", [](double x) { return sleipnir::sign(Variable{x}); },
+      "sign", [](double x) { return sleipnir::sign(Variable{x}); }, "x"_a,
       DOC(sleipnir, sign));
-  autodiff.def("sign", static_cast<Variable (*)(const Variable&)>(&sign),
+  autodiff.def("sign", static_cast<Variable (*)(const Variable&)>(&sign), "x"_a,
                DOC(sleipnir, sign));
   autodiff.def(
-      "sin", [](double x) { return sleipnir::sin(Variable{x}); },
+      "sin", [](double x) { return sleipnir::sin(Variable{x}); }, "x"_a,
       DOC(sleipnir, sin));
-  autodiff.def("sin", static_cast<Variable (*)(const Variable&)>(&sin),
+  autodiff.def("sin", static_cast<Variable (*)(const Variable&)>(&sin), "x"_a,
                DOC(sleipnir, sin));
   autodiff.def(
-      "sinh", [](double x) { return sleipnir::sinh(Variable{x}); },
+      "sinh", [](double x) { return sleipnir::sinh(Variable{x}); }, "x"_a,
       DOC(sleipnir, sinh));
-  autodiff.def("sinh", static_cast<Variable (*)(const Variable&)>(&sinh),
+  autodiff.def("sinh", static_cast<Variable (*)(const Variable&)>(&sinh), "x"_a,
                DOC(sleipnir, sinh));
   autodiff.def(
-      "sqrt", [](double x) { return sleipnir::sqrt(Variable{x}); },
+      "sqrt", [](double x) { return sleipnir::sqrt(Variable{x}); }, "x"_a,
       DOC(sleipnir, sqrt));
-  autodiff.def("sqrt", static_cast<Variable (*)(const Variable&)>(&sqrt),
+  autodiff.def("sqrt", static_cast<Variable (*)(const Variable&)>(&sqrt), "x"_a,
                DOC(sleipnir, sqrt));
   autodiff.def(
-      "tan", [](double x) { return sleipnir::tan(Variable{x}); },
+      "tan", [](double x) { return sleipnir::tan(Variable{x}); }, "x"_a,
       DOC(sleipnir, tan));
-  autodiff.def("tan", static_cast<Variable (*)(const Variable&)>(&tan),
+  autodiff.def("tan", static_cast<Variable (*)(const Variable&)>(&tan), "x"_a,
                DOC(sleipnir, tan));
   autodiff.def(
-      "tanh", [](double x) { return sleipnir::tanh(Variable{x}); },
+      "tanh", [](double x) { return sleipnir::tanh(Variable{x}); }, "x"_a,
       DOC(sleipnir, tanh));
-  autodiff.def("tanh", static_cast<Variable (*)(const Variable&)>(&tanh),
+  autodiff.def("tanh", static_cast<Variable (*)(const Variable&)>(&tanh), "x"_a,
                DOC(sleipnir, tanh));
 }
 
diff --git a/jormungandr/cpp/autodiff/BindVariableBlock.cpp b/jormungandr/cpp/autodiff/BindVariableBlock.cpp
index 19f90bcc..bc4c33ab 100644
--- a/jormungandr/cpp/autodiff/BindVariableBlock.cpp
+++ b/jormungandr/cpp/autodiff/BindVariableBlock.cpp
@@ -17,109 +17,118 @@ namespace py = pybind11;
 
 namespace sleipnir {
 
-void BindVariableBlock(py::module_& autodiff) {
-  py::class_<VariableBlock<VariableMatrix>> variable_block{
-      autodiff, "VariableBlock", DOC(sleipnir, VariableBlock)};
+void BindVariableBlock(py::module_& autodiff,
+                       py::class_<VariableBlock<VariableMatrix>>& cls) {
+  using namespace pybind11::literals;
 
   // VariableBlock-VariableMatrix overloads
-  variable_block.def(py::self * VariableMatrix());
-  variable_block.def(py::self + VariableMatrix());
-  variable_block.def(py::self - VariableMatrix());
-  variable_block.def(py::self == VariableMatrix(), DOC(sleipnir, operator, eq));
-  variable_block.def(py::self < VariableMatrix(), DOC(sleipnir, operator, lt));
-  variable_block.def(py::self <= VariableMatrix(), DOC(sleipnir, operator, le));
-  variable_block.def(py::self > VariableMatrix(), DOC(sleipnir, operator, gt));
-  variable_block.def(py::self >= VariableMatrix(), DOC(sleipnir, operator, ge));
+  cls.def(py::self * VariableMatrix(), "rhs"_a);
+  cls.def(py::self + VariableMatrix(), "rhs"_a);
+  cls.def(py::self - VariableMatrix(), "rhs"_a);
+  cls.def(py::self == VariableMatrix(), "rhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(py::self < VariableMatrix(), "rhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(py::self <= VariableMatrix(), "rhs"_a, DOC(sleipnir, operator, le));
+  cls.def(py::self > VariableMatrix(), "rhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(py::self >= VariableMatrix(), "rhs"_a, DOC(sleipnir, operator, ge));
 
-  variable_block.def(py::init<VariableMatrix&>(),
-                     DOC(sleipnir, VariableBlock, VariableBlock, 3));
-  variable_block.def(py::init<VariableMatrix&, int, int, int, int>(),
-                     DOC(sleipnir, VariableBlock, VariableBlock, 4));
-  variable_block.def(
+  cls.def(py::init<VariableMatrix&>(), "mat"_a,
+          DOC(sleipnir, VariableBlock, VariableBlock, 3));
+  cls.def(py::init<VariableMatrix&, int, int, int, int>(), "mat"_a,
+          "row_offset"_a, "col_offset"_a, "block_rows"_a, "block_cols"_a,
+          DOC(sleipnir, VariableBlock, VariableBlock, 4));
+  cls.def(
       "set_value",
       [](VariableBlock<VariableMatrix>& self, double value) {
         self.SetValue(value);
       },
-      DOC(sleipnir, VariableBlock, SetValue));
-  variable_block.def(
+      "value"_a, DOC(sleipnir, VariableBlock, SetValue));
+  cls.def(
       "set_value",
       [](VariableBlock<VariableMatrix>& self, const Eigen::MatrixXd& values) {
         self.SetValue(values);
       },
-      DOC(sleipnir, VariableBlock, SetValue, 2));
-  variable_block.def("__setitem__",
-                     [](VariableBlock<VariableMatrix>& self, int row,
-                        const Variable& value) { return self(row) = value; });
+      "values"_a, DOC(sleipnir, VariableBlock, SetValue, 2));
+  cls.def(
+      "__setitem__",
+      [](VariableBlock<VariableMatrix>& self, int row, const Variable& value) {
+        return self(row) = value;
+      },
+      "row"_a, "value"_a);
   // TODO: Support slice stride other than 1
-  variable_block.def("__setitem__", [](VariableBlock<VariableMatrix>& self,
-                                       py::tuple slices, py::object value) {
-    if (slices.size() != 2) {
-      throw py::index_error(
-          std::format("Expected 2 slices, got {}.", slices.size()));
-    }
+  cls.def(
+      "__setitem__",
+      [](VariableBlock<VariableMatrix>& self, py::tuple slices,
+         py::object value) {
+        if (slices.size() != 2) {
+          throw py::index_error(
+              std::format("Expected 2 slices, got {}.", slices.size()));
+        }
 
-    int rowOffset = 0;
-    int colOffset = 0;
-    int blockRows = self.Rows();
-    int blockCols = self.Cols();
+        int rowOffset = 0;
+        int colOffset = 0;
+        int blockRows = self.Rows();
+        int blockCols = self.Cols();
 
-    size_t start;
-    size_t stop;
-    size_t step;
-    size_t sliceLength;
+        size_t start;
+        size_t stop;
+        size_t step;
+        size_t sliceLength;
 
-    // Row slice
-    const auto& rowElem = slices[0];
-    if (py::isinstance<py::slice>(rowElem)) {
-      const auto& rowSlice = rowElem.cast<py::slice>();
-      if (!rowSlice.compute(self.Rows(), &start, &stop, &step, &sliceLength)) {
-        throw py::error_already_set();
-      }
-      rowOffset = start;
-      blockRows = stop - start;
-    } else {
-      rowOffset = rowElem.cast<int>();
-      blockRows = 1;
-    }
+        // Row slice
+        const auto& rowElem = slices[0];
+        if (py::isinstance<py::slice>(rowElem)) {
+          const auto& rowSlice = rowElem.cast<py::slice>();
+          if (!rowSlice.compute(self.Rows(), &start, &stop, &step,
+                                &sliceLength)) {
+            throw py::error_already_set();
+          }
+          rowOffset = start;
+          blockRows = stop - start;
+        } else {
+          rowOffset = rowElem.cast<int>();
+          blockRows = 1;
+        }
 
-    // Column slice
-    const auto& colElem = slices[1];
-    if (py::isinstance<py::slice>(colElem)) {
-      const auto& colSlice = colElem.cast<py::slice>();
-      if (!colSlice.compute(self.Cols(), &start, &stop, &step, &sliceLength)) {
-        throw py::error_already_set();
-      }
-      colOffset = start;
-      blockCols = stop - start;
-    } else {
-      colOffset = colElem.cast<int>();
-      blockCols = 1;
-    }
+        // Column slice
+        const auto& colElem = slices[1];
+        if (py::isinstance<py::slice>(colElem)) {
+          const auto& colSlice = colElem.cast<py::slice>();
+          if (!colSlice.compute(self.Cols(), &start, &stop, &step,
+                                &sliceLength)) {
+            throw py::error_already_set();
+          }
+          colOffset = start;
+          blockCols = stop - start;
+        } else {
+          colOffset = colElem.cast<int>();
+          blockCols = 1;
+        }
 
-    if (py::isinstance<VariableMatrix>(value)) {
-      self.Block(rowOffset, colOffset, blockRows, blockCols) =
-          value.cast<VariableMatrix>();
-    } else if (py::isinstance<VariableBlock<VariableMatrix>>(value)) {
-      self.Block(rowOffset, colOffset, blockRows, blockCols) =
-          value.cast<VariableBlock<VariableMatrix>>();
-    } else if (IsNumPyArithmeticArray(value)) {
-      self.Block(rowOffset, colOffset, blockRows, blockCols) =
-          value.cast<Eigen::MatrixXd>();
-    } else if (py::isinstance<Variable>(value)) {
-      self.Block(rowOffset, colOffset, blockRows, blockCols) =
-          value.cast<Variable>();
-    } else if (py::isinstance<py::float_>(value)) {
-      self.Block(rowOffset, colOffset, blockRows, blockCols) =
-          value.cast<double>();
-    } else if (py::isinstance<py::int_>(value)) {
-      self.Block(rowOffset, colOffset, blockRows, blockCols) =
-          value.cast<int>();
-    } else {
-      throw py::value_error(
-          "VariableBlock.__setitem__ not implemented for value");
-    }
-  });
-  variable_block.def(
+        if (py::isinstance<VariableMatrix>(value)) {
+          self.Block(rowOffset, colOffset, blockRows, blockCols) =
+              value.cast<VariableMatrix>();
+        } else if (py::isinstance<VariableBlock<VariableMatrix>>(value)) {
+          self.Block(rowOffset, colOffset, blockRows, blockCols) =
+              value.cast<VariableBlock<VariableMatrix>>();
+        } else if (IsNumPyArithmeticArray(value)) {
+          self.Block(rowOffset, colOffset, blockRows, blockCols) =
+              value.cast<Eigen::MatrixXd>();
+        } else if (py::isinstance<Variable>(value)) {
+          self.Block(rowOffset, colOffset, blockRows, blockCols) =
+              value.cast<Variable>();
+        } else if (py::isinstance<py::float_>(value)) {
+          self.Block(rowOffset, colOffset, blockRows, blockCols) =
+              value.cast<double>();
+        } else if (py::isinstance<py::int_>(value)) {
+          self.Block(rowOffset, colOffset, blockRows, blockCols) =
+              value.cast<int>();
+        } else {
+          throw py::value_error(
+              "VariableBlock.__setitem__ not implemented for value");
+        }
+      },
+      "slices"_a, "value"_a);
+  cls.def(
       "__getitem__",
       [](VariableBlock<VariableMatrix>& self, int row) -> Variable& {
         if (row < 0) {
@@ -127,9 +136,10 @@ void BindVariableBlock(py::module_& autodiff) {
         }
         return self(row);
       },
-      py::keep_alive<0, 1>(), DOC(sleipnir, VariableBlock, operator, call, 3));
+      py::keep_alive<0, 1>(), "row"_a,
+      DOC(sleipnir, VariableBlock, operator, call, 3));
   // TODO: Support slice stride other than 1
-  variable_block.def(
+  cls.def(
       "__getitem__",
       [](VariableBlock<VariableMatrix>& self, py::tuple slices) -> py::object {
         if (slices.size() != 2) {
@@ -205,31 +215,30 @@ void BindVariableBlock(py::module_& autodiff) {
 
         return py::cast(self.Block(rowOffset, colOffset, blockRows, blockCols));
       },
-      py::keep_alive<0, 1>(), DOC(sleipnir, VariableBlock, operator, call));
-  variable_block.def(
-      "row", py::overload_cast<int>(&VariableBlock<VariableMatrix>::Row),
-      DOC(sleipnir, VariableBlock, Row));
-  variable_block.def(
-      "col", py::overload_cast<int>(&VariableBlock<VariableMatrix>::Col),
-      DOC(sleipnir, VariableBlock, Col));
-  variable_block.def(
+      py::keep_alive<0, 1>(), DOC(sleipnir, VariableBlock, operator, call),
+      "slices"_a);
+  cls.def("row", py::overload_cast<int>(&VariableBlock<VariableMatrix>::Row),
+          "row"_a, DOC(sleipnir, VariableBlock, Row));
+  cls.def("col", py::overload_cast<int>(&VariableBlock<VariableMatrix>::Col),
+          "col"_a, DOC(sleipnir, VariableBlock, Col));
+  cls.def(
       "__mul__",
       [](const VariableBlock<VariableMatrix>& lhs,
          const VariableBlock<VariableMatrix>& rhs) { return lhs * rhs; },
-      py::is_operator());
-  variable_block.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(
       "__rmul__",
       [](const VariableBlock<VariableMatrix>& rhs,
          const VariableBlock<VariableMatrix>& lhs) { return lhs * rhs; },
-      py::is_operator());
-  variable_block.def(
+      py::is_operator(), "lhs"_a);
+  cls.def(
       "__matmul__",
       [](const VariableBlock<VariableMatrix>& lhs,
          const VariableBlock<VariableMatrix>& rhs) { return lhs * rhs; },
-      py::is_operator());
+      py::is_operator(), "rhs"_a);
 
   // https://numpy.org/doc/stable/user/basics.dispatch.html
-  variable_block.def(
+  cls.def(
       "__array_ufunc__",
       [](VariableBlock<VariableMatrix>& self, py::object ufunc, py::str method,
          py::args inputs, const py::kwargs& kwargs) -> py::object {
@@ -295,158 +304,158 @@ void BindVariableBlock(py::module_& autodiff) {
                         "implemented for ({}, {})",
                         method_name, ufunc_name, input1_name, input2_name));
         return py::cast(VariableMatrix{self});
-      });
+      },
+      "ufunc"_a, "method"_a);
 
-  variable_block.def(py::self * Variable());
-  variable_block.def(py::self * double());
-  variable_block.def(Variable() * py::self);
-  variable_block.def(double() * py::self);
-  variable_block.def(py::self / Variable());
-  variable_block.def(py::self / double());
-  variable_block.def(py::self + py::self);
-  variable_block.def(
+  cls.def(py::self * Variable(), "rhs"_a);
+  cls.def(py::self * double(), "rhs"_a);
+  cls.def(Variable() * py::self, "lhs"_a);
+  cls.def(double() * py::self, "lhs"_a);
+  cls.def(py::self / Variable(), "rhs"_a);
+  cls.def(py::self / double(), "rhs"_a);
+  cls.def(py::self + py::self, "rhs"_a);
+  cls.def(
       "__add__",
       [](const VariableBlock<VariableMatrix>& lhs,
          const Eigen::Ref<const Eigen::MatrixXd>& rhs) -> VariableMatrix {
         return lhs + rhs;
       },
-      py::is_operator());
-  variable_block.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(
       "__add__",
       [](const Eigen::Ref<const Eigen::MatrixXd>& lhs,
          const VariableBlock<VariableMatrix>& rhs) -> VariableMatrix {
         return lhs + rhs;
       },
-      py::is_operator());
-  variable_block.def(py::self - py::self);
-  variable_block.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(py::self - py::self, "rhs"_a);
+  cls.def(
       "__sub__",
       [](const VariableBlock<VariableMatrix>& lhs,
          const Eigen::Ref<const Eigen::MatrixXd>& rhs) -> VariableMatrix {
         return lhs - rhs;
       },
-      py::is_operator());
-  variable_block.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(
       "__sub__",
       [](const Eigen::Ref<const Eigen::MatrixXd>& lhs,
          const VariableBlock<VariableMatrix>& rhs) -> VariableMatrix {
         return lhs - rhs;
       },
-      py::is_operator());
-  variable_block.def(-py::self);
-  variable_block.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(-py::self);
+  cls.def(
       "__pow__",
       [](const VariableBlock<VariableMatrix>& self, int power) {
         return sleipnir::pow(VariableMatrix{self}, power);
       },
-      py::is_operator());
-  variable_block.def_property_readonly("T", &VariableBlock<VariableMatrix>::T,
-                                       DOC(sleipnir, VariableBlock, T));
-  variable_block.def("rows", &VariableBlock<VariableMatrix>::Rows,
-                     DOC(sleipnir, VariableBlock, Rows));
-  variable_block.def("cols", &VariableBlock<VariableMatrix>::Cols,
-                     DOC(sleipnir, VariableBlock, Cols));
-  variable_block.def_property_readonly(
-      "shape", [](const VariableBlock<VariableMatrix>& self) {
-        return py::make_tuple(self.Rows(), self.Cols());
-      });
-  variable_block.def(
+      py::is_operator(), "power"_a);
+  cls.def_property_readonly("T", &VariableBlock<VariableMatrix>::T,
+                            DOC(sleipnir, VariableBlock, T));
+  cls.def("rows", &VariableBlock<VariableMatrix>::Rows,
+          DOC(sleipnir, VariableBlock, Rows));
+  cls.def("cols", &VariableBlock<VariableMatrix>::Cols,
+          DOC(sleipnir, VariableBlock, Cols));
+  cls.def_property_readonly("shape",
+                            [](const VariableBlock<VariableMatrix>& self) {
+                              return py::make_tuple(self.Rows(), self.Cols());
+                            });
+  cls.def(
       "value",
       static_cast<double (VariableBlock<VariableMatrix>::*)(int, int) const>(
           &VariableBlock<VariableMatrix>::Value),
-      DOC(sleipnir, VariableBlock, Value));
-  variable_block.def(
-      "value",
-      static_cast<double (VariableBlock<VariableMatrix>::*)(int) const>(
-          &VariableBlock<VariableMatrix>::Value),
-      DOC(sleipnir, VariableBlock, Value, 2));
-  variable_block.def(
+      "row"_a, "col"_a, DOC(sleipnir, VariableBlock, Value));
+  cls.def("value",
+          static_cast<double (VariableBlock<VariableMatrix>::*)(int) const>(
+              &VariableBlock<VariableMatrix>::Value),
+          "index"_a, DOC(sleipnir, VariableBlock, Value, 2));
+  cls.def(
       "value",
       static_cast<Eigen::MatrixXd (VariableBlock<VariableMatrix>::*)() const>(
           &VariableBlock<VariableMatrix>::Value),
       DOC(sleipnir, VariableBlock, Value, 3));
-  variable_block.def(
+  cls.def(
       "cwise_transform",
       [](const VariableBlock<VariableMatrix>& self,
          const std::function<Variable(const Variable&)>& func) {
         return self.CwiseTransform(func);
       },
-      DOC(sleipnir, VariableBlock, CwiseTransform));
-  variable_block.def(py::self == py::self, DOC(sleipnir, operator, eq));
-  variable_block.def(py::self == Variable(), DOC(sleipnir, operator, eq));
-  variable_block.def(py::self == double(), DOC(sleipnir, operator, eq));
-  variable_block.def(py::self == int(), DOC(sleipnir, operator, eq));
-  variable_block.def(Variable() == py::self, DOC(sleipnir, operator, eq));
-  variable_block.def(double() == py::self, DOC(sleipnir, operator, eq));
-  variable_block.def(int() == py::self, DOC(sleipnir, operator, eq));
-  variable_block.def(py::self < py::self, DOC(sleipnir, operator, lt));
-  variable_block.def(py::self < Variable(), DOC(sleipnir, operator, lt));
-  variable_block.def(py::self < double(), DOC(sleipnir, operator, lt));
-  variable_block.def(py::self < int(), DOC(sleipnir, operator, lt));
-  variable_block.def(Variable() < py::self, DOC(sleipnir, operator, lt));
-  variable_block.def(double() < py::self, DOC(sleipnir, operator, lt));
-  variable_block.def(int() < py::self, DOC(sleipnir, operator, lt));
-  variable_block.def(py::self <= py::self, DOC(sleipnir, operator, le));
-  variable_block.def(py::self <= Variable(), DOC(sleipnir, operator, le));
-  variable_block.def(py::self <= double(), DOC(sleipnir, operator, le));
-  variable_block.def(py::self <= int(), DOC(sleipnir, operator, le));
-  variable_block.def(Variable() <= py::self, DOC(sleipnir, operator, le));
-  variable_block.def(double() <= py::self, DOC(sleipnir, operator, le));
-  variable_block.def(int() <= py::self, DOC(sleipnir, operator, le));
-  variable_block.def(py::self > py::self, DOC(sleipnir, operator, gt));
-  variable_block.def(py::self > Variable(), DOC(sleipnir, operator, gt));
-  variable_block.def(py::self > double(), DOC(sleipnir, operator, gt));
-  variable_block.def(py::self > int(), DOC(sleipnir, operator, gt));
-  variable_block.def(Variable() > py::self, DOC(sleipnir, operator, gt));
-  variable_block.def(double() > py::self, DOC(sleipnir, operator, gt));
-  variable_block.def(int() > py::self, DOC(sleipnir, operator, gt));
-  variable_block.def(py::self >= py::self, DOC(sleipnir, operator, ge));
-  variable_block.def(py::self >= Variable(), DOC(sleipnir, operator, ge));
-  variable_block.def(py::self >= double(), DOC(sleipnir, operator, ge));
-  variable_block.def(py::self >= int(), DOC(sleipnir, operator, ge));
-  variable_block.def(Variable() >= py::self, DOC(sleipnir, operator, ge));
-  variable_block.def(double() >= py::self, DOC(sleipnir, operator, ge));
-  variable_block.def(int() >= py::self, DOC(sleipnir, operator, ge));
-  variable_block.def(
+      "func"_a, DOC(sleipnir, VariableBlock, CwiseTransform));
+  cls.def(py::self == py::self, "rhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(py::self == Variable(), "rhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(py::self == double(), "rhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(py::self == int(), "rhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(Variable() == py::self, "lhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(double() == py::self, "lhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(int() == py::self, "lhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(py::self < py::self, "rhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(py::self < Variable(), "rhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(py::self < double(), "rhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(py::self < int(), "rhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(Variable() < py::self, "lhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(double() < py::self, "lhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(int() < py::self, "lhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(py::self <= py::self, "rhs"_a, DOC(sleipnir, operator, le));
+  cls.def(py::self <= Variable(), "rhs"_a, DOC(sleipnir, operator, le));
+  cls.def(py::self <= double(), "rhs"_a, DOC(sleipnir, operator, le));
+  cls.def(py::self <= int(), "rhs"_a, DOC(sleipnir, operator, le));
+  cls.def(Variable() <= py::self, "lhs"_a, DOC(sleipnir, operator, le));
+  cls.def(double() <= py::self, "lhs"_a, DOC(sleipnir, operator, le));
+  cls.def(int() <= py::self, "lhs"_a, DOC(sleipnir, operator, le));
+  cls.def(py::self > py::self, "rhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(py::self > Variable(), "rhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(py::self > double(), "rhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(py::self > int(), "rhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(Variable() > py::self, "lhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(double() > py::self, "lhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(int() > py::self, "lhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(py::self >= py::self, "rhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(py::self >= Variable(), "rhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(py::self >= double(), "rhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(py::self >= int(), "rhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(Variable() >= py::self, "lhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(double() >= py::self, "lhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(int() >= py::self, "lhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(
       "__eq__",
       [](const VariableBlock<VariableMatrix>& lhs,
          const py::array_t<double>& rhs) {
         return lhs == rhs.cast<Eigen::MatrixXd>();
       },
-      py::is_operator(), DOC(sleipnir, operator, eq));
-  variable_block.def(
+      py::is_operator(), "rhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(
       "__lt__",
       [](const VariableBlock<VariableMatrix>& lhs,
          const py::array_t<double>& rhs) {
         return lhs < rhs.cast<Eigen::MatrixXd>();
       },
-      py::is_operator(), DOC(sleipnir, operator, lt));
-  variable_block.def(
+      py::is_operator(), "rhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(
       "__le__",
       [](const VariableBlock<VariableMatrix>& lhs,
          const py::array_t<double>& rhs) {
         return lhs <= rhs.cast<Eigen::MatrixXd>();
       },
-      py::is_operator(), DOC(sleipnir, operator, le));
-  variable_block.def(
+      py::is_operator(), "rhs"_a, DOC(sleipnir, operator, le));
+  cls.def(
       "__gt__",
       [](const VariableBlock<VariableMatrix>& lhs,
          const py::array_t<double>& rhs) {
         return lhs > rhs.cast<Eigen::MatrixXd>();
       },
-      py::is_operator(), DOC(sleipnir, operator, gt));
-  variable_block.def(
+      py::is_operator(), "rhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(
       "__ge__",
       [](const VariableBlock<VariableMatrix>& lhs,
          const py::array_t<double>& rhs) {
         return lhs >= rhs.cast<Eigen::MatrixXd>();
       },
-      py::is_operator(), DOC(sleipnir, operator, ge));
+      py::is_operator(), "rhs"_a, DOC(sleipnir, operator, ge));
 
-  variable_block.def("__len__", &VariableBlock<VariableMatrix>::Rows,
-                     DOC(sleipnir, VariableBlock, Rows));
+  cls.def("__len__", &VariableBlock<VariableMatrix>::Rows,
+          DOC(sleipnir, VariableBlock, Rows));
 
-  variable_block.def(
+  cls.def(
       "__iter__",
       [](const VariableBlock<VariableMatrix>& self) {
         return py::make_iterator(self.begin(), self.end());
@@ -462,7 +471,7 @@ void BindVariableBlock(py::module_& autodiff) {
          const std::function<Variable(const Variable&, const Variable&)> func) {
         return CwiseReduce(lhs, rhs, func);
       },
-      DOC(sleipnir, CwiseReduce));
+      "lhs"_a, "rhs"_a, "func"_a, DOC(sleipnir, CwiseReduce));
 }
 
 }  // namespace sleipnir
diff --git a/jormungandr/cpp/autodiff/BindVariableMatrix.cpp b/jormungandr/cpp/autodiff/BindVariableMatrix.cpp
index bbdc8b69..4791e7f6 100644
--- a/jormungandr/cpp/autodiff/BindVariableMatrix.cpp
+++ b/jormungandr/cpp/autodiff/BindVariableMatrix.cpp
@@ -17,104 +17,109 @@ namespace py = pybind11;
 
 namespace sleipnir {
 
-void BindVariableMatrix(py::module_& autodiff) {
-  py::class_<VariableMatrix> variable_matrix{autodiff, "VariableMatrix",
-                                             DOC(sleipnir, VariableMatrix)};
-
-  variable_matrix.def(py::init<>(),
-                      DOC(sleipnir, VariableMatrix, VariableMatrix));
-  variable_matrix.def(py::init<int>(),
-                      DOC(sleipnir, VariableMatrix, VariableMatrix, 2));
-  variable_matrix.def(py::init<int, int>(),
-                      DOC(sleipnir, VariableMatrix, VariableMatrix, 3));
-  variable_matrix.def(py::init<std::vector<std::vector<double>>>(),
-                      DOC(sleipnir, VariableMatrix, VariableMatrix, 5));
-  variable_matrix.def(py::init<std::vector<std::vector<Variable>>>(),
-                      DOC(sleipnir, VariableMatrix, VariableMatrix, 6));
-  variable_matrix.def(py::init<const Variable&>(),
-                      DOC(sleipnir, VariableMatrix, VariableMatrix, 9));
-  variable_matrix.def(py::init<const VariableBlock<VariableMatrix>&>(),
-                      DOC(sleipnir, VariableMatrix, VariableMatrix, 11));
-  variable_matrix.def(
+void BindVariableMatrix(py::module_& autodiff,
+                        py::class_<VariableMatrix>& cls) {
+  using namespace pybind11::literals;
+
+  cls.def(py::init<>(), DOC(sleipnir, VariableMatrix, VariableMatrix));
+  cls.def(py::init<int>(), "rows"_a,
+          DOC(sleipnir, VariableMatrix, VariableMatrix, 2));
+  cls.def(py::init<int, int>(), "rows"_a, "cols"_a,
+          DOC(sleipnir, VariableMatrix, VariableMatrix, 3));
+  cls.def(py::init<std::vector<std::vector<double>>>(), "list"_a,
+          DOC(sleipnir, VariableMatrix, VariableMatrix, 5));
+  cls.def(py::init<std::vector<std::vector<Variable>>>(), "list"_a,
+          DOC(sleipnir, VariableMatrix, VariableMatrix, 6));
+  cls.def(py::init<const Variable&>(), "values"_a,
+          DOC(sleipnir, VariableMatrix, VariableMatrix, 9));
+  cls.def(py::init<const VariableBlock<VariableMatrix>&>(), "values"_a,
+          DOC(sleipnir, VariableMatrix, VariableMatrix, 11));
+  cls.def(
       "set_value",
       [](VariableMatrix& self, const Eigen::MatrixXd& values) {
         self.SetValue(values);
       },
-      DOC(sleipnir, VariableMatrix, SetValue));
-  variable_matrix.def("__setitem__",
-                      [](VariableMatrix& self, int row, const Variable& value) {
-                        return self(row) = value;
-                      });
+      "values"_a, DOC(sleipnir, VariableMatrix, SetValue));
+  cls.def(
+      "__setitem__",
+      [](VariableMatrix& self, int row, const Variable& value) {
+        return self(row) = value;
+      },
+      "row"_a, "value"_a);
   // TODO: Support slice stride other than 1
-  variable_matrix.def("__setitem__", [](VariableMatrix& self, py::tuple slices,
-                                        py::object value) {
-    if (slices.size() != 2) {
-      throw py::index_error(
-          std::format("Expected 2 slices, got {}.", slices.size()));
-    }
-
-    int rowOffset = 0;
-    int colOffset = 0;
-    int blockRows = self.Rows();
-    int blockCols = self.Cols();
-
-    size_t start;
-    size_t stop;
-    size_t step;
-    size_t sliceLength;
-
-    // Row slice
-    const auto& rowElem = slices[0];
-    if (py::isinstance<py::slice>(rowElem)) {
-      const auto& rowSlice = rowElem.cast<py::slice>();
-      if (!rowSlice.compute(self.Rows(), &start, &stop, &step, &sliceLength)) {
-        throw py::error_already_set();
-      }
-      rowOffset = start;
-      blockRows = stop - start;
-    } else {
-      rowOffset = rowElem.cast<int>();
-      blockRows = 1;
-    }
-
-    // Column slice
-    const auto& colElem = slices[1];
-    if (py::isinstance<py::slice>(colElem)) {
-      const auto& colSlice = colElem.cast<py::slice>();
-      if (!colSlice.compute(self.Cols(), &start, &stop, &step, &sliceLength)) {
-        throw py::error_already_set();
-      }
-      colOffset = start;
-      blockCols = stop - start;
-    } else {
-      colOffset = colElem.cast<int>();
-      blockCols = 1;
-    }
-
-    if (py::isinstance<VariableMatrix>(value)) {
-      self.Block(rowOffset, colOffset, blockRows, blockCols) =
-          value.cast<VariableMatrix>();
-    } else if (py::isinstance<VariableBlock<VariableMatrix>>(value)) {
-      self.Block(rowOffset, colOffset, blockRows, blockCols) =
-          value.cast<VariableBlock<VariableMatrix>>();
-    } else if (IsNumPyArithmeticArray(value)) {
-      self.Block(rowOffset, colOffset, blockRows, blockCols) =
-          value.cast<Eigen::MatrixXd>();
-    } else if (py::isinstance<Variable>(value)) {
-      self.Block(rowOffset, colOffset, blockRows, blockCols) =
-          value.cast<Variable>();
-    } else if (py::isinstance<py::float_>(value)) {
-      self.Block(rowOffset, colOffset, blockRows, blockCols) =
-          value.cast<double>();
-    } else if (py::isinstance<py::int_>(value)) {
-      self.Block(rowOffset, colOffset, blockRows, blockCols) =
-          value.cast<int>();
-    } else {
-      throw py::value_error(
-          "VariableMatrix.__setitem__ not implemented for value");
-    }
-  });
-  variable_matrix.def(
+  cls.def(
+      "__setitem__",
+      [](VariableMatrix& self, py::tuple slices, py::object value) {
+        if (slices.size() != 2) {
+          throw py::index_error(
+              std::format("Expected 2 slices, got {}.", slices.size()));
+        }
+
+        int rowOffset = 0;
+        int colOffset = 0;
+        int blockRows = self.Rows();
+        int blockCols = self.Cols();
+
+        size_t start;
+        size_t stop;
+        size_t step;
+        size_t sliceLength;
+
+        // Row slice
+        const auto& rowElem = slices[0];
+        if (py::isinstance<py::slice>(rowElem)) {
+          const auto& rowSlice = rowElem.cast<py::slice>();
+          if (!rowSlice.compute(self.Rows(), &start, &stop, &step,
+                                &sliceLength)) {
+            throw py::error_already_set();
+          }
+          rowOffset = start;
+          blockRows = stop - start;
+        } else {
+          rowOffset = rowElem.cast<int>();
+          blockRows = 1;
+        }
+
+        // Column slice
+        const auto& colElem = slices[1];
+        if (py::isinstance<py::slice>(colElem)) {
+          const auto& colSlice = colElem.cast<py::slice>();
+          if (!colSlice.compute(self.Cols(), &start, &stop, &step,
+                                &sliceLength)) {
+            throw py::error_already_set();
+          }
+          colOffset = start;
+          blockCols = stop - start;
+        } else {
+          colOffset = colElem.cast<int>();
+          blockCols = 1;
+        }
+
+        if (py::isinstance<VariableMatrix>(value)) {
+          self.Block(rowOffset, colOffset, blockRows, blockCols) =
+              value.cast<VariableMatrix>();
+        } else if (py::isinstance<VariableBlock<VariableMatrix>>(value)) {
+          self.Block(rowOffset, colOffset, blockRows, blockCols) =
+              value.cast<VariableBlock<VariableMatrix>>();
+        } else if (IsNumPyArithmeticArray(value)) {
+          self.Block(rowOffset, colOffset, blockRows, blockCols) =
+              value.cast<Eigen::MatrixXd>();
+        } else if (py::isinstance<Variable>(value)) {
+          self.Block(rowOffset, colOffset, blockRows, blockCols) =
+              value.cast<Variable>();
+        } else if (py::isinstance<py::float_>(value)) {
+          self.Block(rowOffset, colOffset, blockRows, blockCols) =
+              value.cast<double>();
+        } else if (py::isinstance<py::int_>(value)) {
+          self.Block(rowOffset, colOffset, blockRows, blockCols) =
+              value.cast<int>();
+        } else {
+          throw py::value_error(
+              "VariableMatrix.__setitem__ not implemented for value");
+        }
+      },
+      "slices"_a, "value"_a);
+  cls.def(
       "__getitem__",
       [](VariableMatrix& self, int row) -> Variable& {
         if (row < 0) {
@@ -122,9 +127,10 @@ void BindVariableMatrix(py::module_& autodiff) {
         }
         return self(row);
       },
-      py::keep_alive<0, 1>(), DOC(sleipnir, VariableMatrix, operator, call, 3));
+      py::keep_alive<0, 1>(), "row"_a,
+      DOC(sleipnir, VariableMatrix, operator, call, 3));
   // TODO: Support slice stride other than 1
-  variable_matrix.def(
+  cls.def(
       "__getitem__",
       [](VariableMatrix& self, py::tuple slices) -> py::object {
         if (slices.size() != 2) {
@@ -200,32 +206,33 @@ void BindVariableMatrix(py::module_& autodiff) {
 
         return py::cast(self.Block(rowOffset, colOffset, blockRows, blockCols));
       },
-      py::keep_alive<0, 1>(), DOC(sleipnir, VariableMatrix, operator, call));
-  variable_matrix.def("row", py::overload_cast<int>(&VariableMatrix::Row),
-                      DOC(sleipnir, VariableMatrix, Row));
-  variable_matrix.def("col", py::overload_cast<int>(&VariableMatrix::Col),
-                      DOC(sleipnir, VariableMatrix, Col));
-  variable_matrix.def(
+      py::keep_alive<0, 1>(), "slices"_a,
+      DOC(sleipnir, VariableMatrix, operator, call));
+  cls.def("row", py::overload_cast<int>(&VariableMatrix::Row), "row"_a,
+          DOC(sleipnir, VariableMatrix, Row));
+  cls.def("col", py::overload_cast<int>(&VariableMatrix::Col), "col"_a,
+          DOC(sleipnir, VariableMatrix, Col));
+  cls.def(
       "__mul__",
       [](const VariableMatrix& lhs, const VariableMatrix& rhs) {
         return lhs * rhs;
       },
-      py::is_operator());
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(
       "__rmul__",
       [](const VariableMatrix& rhs, const VariableMatrix& lhs) {
         return lhs * rhs;
       },
-      py::is_operator());
-  variable_matrix.def(
+      py::is_operator(), "lhs"_a);
+  cls.def(
       "__matmul__",
       [](const VariableMatrix& lhs, const VariableMatrix& rhs) {
         return lhs * rhs;
       },
-      py::is_operator());
+      py::is_operator(), "rhs"_a);
 
   // https://numpy.org/doc/stable/user/basics.dispatch.html
-  variable_matrix.def(
+  cls.def(
       "__array_ufunc__",
       [](VariableMatrix& self, py::object ufunc, py::str method,
          py::args inputs, const py::kwargs& kwargs) -> py::object {
@@ -291,261 +298,259 @@ void BindVariableMatrix(py::module_& autodiff) {
                         "implemented for ({}, {})",
                         method_name, ufunc_name, input1_name, input2_name));
         return py::cast(VariableMatrix{self});
-      });
-
-  variable_matrix.def(py::self * Variable());
-  variable_matrix.def(py::self * double());
-  variable_matrix.def(Variable() * py::self);
-  variable_matrix.def(double() * py::self);
-
-  variable_matrix.def(py::self / Variable());
-  variable_matrix.def(py::self / double());
-  variable_matrix.def(py::self /= Variable(),
-                      DOC(sleipnir, VariableMatrix, operator, idiv));
-  variable_matrix.def(py::self /= double(),
-                      DOC(sleipnir, VariableMatrix, operator, idiv));
-
-  variable_matrix.def(py::self + py::self);
-  variable_matrix.def(
+      },
+      "ufunc"_a, "method"_a);
+
+  cls.def(py::self * Variable(), "rhs"_a);
+  cls.def(py::self * double(), "rhs"_a);
+  cls.def(Variable() * py::self, "lhs"_a);
+  cls.def(double() * py::self, "lhs"_a);
+
+  cls.def(py::self / Variable(), "rhs"_a);
+  cls.def(py::self / double(), "rhs"_a);
+  cls.def(py::self /= Variable(), "rhs"_a,
+          DOC(sleipnir, VariableMatrix, operator, idiv));
+  cls.def(py::self /= double(), "rhs"_a,
+          DOC(sleipnir, VariableMatrix, operator, idiv));
+
+  cls.def(py::self + py::self, "rhs"_a);
+  cls.def(
       "__add__",
       [](const VariableMatrix& lhs, const Variable& rhs) {
         return lhs + VariableMatrix{rhs};
       },
-      py::is_operator());
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(
       "__radd__",
       [](const VariableMatrix& rhs, const Variable& lhs) {
         return VariableMatrix{lhs} + rhs;
       },
-      py::is_operator());
-  variable_matrix.def(
+      py::is_operator(), "lhs"_a);
+  cls.def(
       "__add__",
       [](double lhs, const VariableMatrix& rhs) {
         return VariableMatrix{Variable{lhs}} + rhs;
       },
-      py::is_operator());
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(
       "__radd__",
       [](const VariableMatrix& rhs, double lhs) {
         return VariableMatrix{Variable{lhs}} + rhs;
       },
-      py::is_operator());
-  variable_matrix.def(
+      py::is_operator(), "lhs"_a);
+  cls.def(
       "__add__",
       [](const VariableMatrix& lhs,
          const Eigen::Ref<const Eigen::MatrixXd>& rhs) -> VariableMatrix {
         return lhs + rhs;
       },
-      py::is_operator());
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(
       "__radd__",
       [](const VariableMatrix& rhs,
          const Eigen::Ref<const Eigen::MatrixXd>& lhs) -> VariableMatrix {
         return lhs + rhs;
       },
-      py::is_operator());
+      py::is_operator(), "lhs"_a);
 
-  variable_matrix.def(py::self - py::self);
-  variable_matrix.def(
+  cls.def(py::self - py::self, "rhs"_a);
+  cls.def(
       "__sub__",
       [](const VariableMatrix& lhs, const Variable& rhs) {
         return lhs - VariableMatrix{rhs};
       },
-      py::is_operator());
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(
       "__rsub__",
       [](const VariableMatrix& rhs, const Variable& lhs) {
         return VariableMatrix{lhs} - rhs;
       },
-      py::is_operator());
-  variable_matrix.def(
+      py::is_operator(), "lhs"_a);
+  cls.def(
       "__sub__",
       [](const VariableMatrix& lhs,
          const Eigen::Ref<const Eigen::MatrixXd>& rhs) -> VariableMatrix {
         return lhs - rhs;
       },
-      py::is_operator());
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(
       "__rsub__",
       [](const VariableMatrix& rhs,
          const Eigen::Ref<const Eigen::MatrixXd>& lhs) -> VariableMatrix {
         return lhs - rhs;
       },
-      py::is_operator());
+      py::is_operator(), "lhs"_a);
 
-  variable_matrix.def(-py::self);
-  variable_matrix.def(
+  cls.def(-py::self);
+  cls.def(
       "__pow__",
       [](const VariableMatrix& self, int power) {
         return sleipnir::pow(self, power);
       },
-      py::is_operator());
-  variable_matrix.def_property_readonly("T", &VariableMatrix::T,
-                                        DOC(sleipnir, VariableMatrix, T));
-  variable_matrix.def("rows", &VariableMatrix::Rows,
-                      DOC(sleipnir, VariableMatrix, Rows));
-  variable_matrix.def("cols", &VariableMatrix::Cols,
-                      DOC(sleipnir, VariableMatrix, Cols));
-  variable_matrix.def_property_readonly(
-      "shape", [](const VariableMatrix& self) {
-        return py::make_tuple(self.Rows(), self.Cols());
-      });
-  variable_matrix.def("value",
-                      static_cast<double (VariableMatrix::*)(int, int) const>(
-                          &VariableMatrix::Value),
-                      DOC(sleipnir, VariableMatrix, Value));
-  variable_matrix.def("value",
-                      static_cast<double (VariableMatrix::*)(int) const>(
-                          &VariableMatrix::Value),
-                      DOC(sleipnir, VariableMatrix, Value, 2));
-  variable_matrix.def("value",
-                      static_cast<Eigen::MatrixXd (VariableMatrix::*)() const>(
-                          &VariableMatrix::Value),
-                      DOC(sleipnir, VariableMatrix, Value, 3));
-  variable_matrix.def(
+      py::is_operator(), "power"_a);
+  cls.def_property_readonly("T", &VariableMatrix::T,
+                            DOC(sleipnir, VariableMatrix, T));
+  cls.def("rows", &VariableMatrix::Rows, DOC(sleipnir, VariableMatrix, Rows));
+  cls.def("cols", &VariableMatrix::Cols, DOC(sleipnir, VariableMatrix, Cols));
+  cls.def_property_readonly("shape", [](const VariableMatrix& self) {
+    return py::make_tuple(self.Rows(), self.Cols());
+  });
+  cls.def("value",
+          static_cast<double (VariableMatrix::*)(int, int) const>(
+              &VariableMatrix::Value),
+          "row"_a, "col"_a, DOC(sleipnir, VariableMatrix, Value));
+  cls.def("value",
+          static_cast<double (VariableMatrix::*)(int) const>(
+              &VariableMatrix::Value),
+          "index"_a, DOC(sleipnir, VariableMatrix, Value, 2));
+  cls.def("value",
+          static_cast<Eigen::MatrixXd (VariableMatrix::*)() const>(
+              &VariableMatrix::Value),
+          DOC(sleipnir, VariableMatrix, Value, 3));
+  cls.def(
       "cwise_transform",
       [](const VariableMatrix& self,
          const std::function<Variable(const Variable&)>& func) {
         return self.CwiseTransform(func);
       },
-      DOC(sleipnir, VariableMatrix, CwiseTransform));
-  variable_matrix.def_static("zero", &VariableMatrix::Zero,
-                             DOC(sleipnir, VariableMatrix, Zero));
-  variable_matrix.def_static("ones", &VariableMatrix::Ones,
-                             DOC(sleipnir, VariableMatrix, Ones));
-  variable_matrix.def(py::self == py::self, DOC(sleipnir, operator, eq));
-  variable_matrix.def(py::self == Variable(), DOC(sleipnir, operator, eq));
-  variable_matrix.def(py::self == double(), DOC(sleipnir, operator, eq));
-  variable_matrix.def(py::self == int(), DOC(sleipnir, operator, eq));
-  variable_matrix.def(Variable() == py::self, DOC(sleipnir, operator, eq));
-  variable_matrix.def(double() == py::self, DOC(sleipnir, operator, eq));
-  variable_matrix.def(int() == py::self, DOC(sleipnir, operator, eq));
-  variable_matrix.def(py::self < py::self, DOC(sleipnir, operator, lt));
-  variable_matrix.def(py::self < Variable(), DOC(sleipnir, operator, lt));
-  variable_matrix.def(py::self < double(), DOC(sleipnir, operator, lt));
-  variable_matrix.def(py::self < int(), DOC(sleipnir, operator, lt));
-  variable_matrix.def(Variable() < py::self, DOC(sleipnir, operator, lt));
-  variable_matrix.def(double() < py::self, DOC(sleipnir, operator, lt));
-  variable_matrix.def(int() < py::self, DOC(sleipnir, operator, lt));
-  variable_matrix.def(py::self <= py::self, DOC(sleipnir, operator, le));
-  variable_matrix.def(py::self <= Variable(), DOC(sleipnir, operator, le));
-  variable_matrix.def(py::self <= double(), DOC(sleipnir, operator, le));
-  variable_matrix.def(py::self <= int(), DOC(sleipnir, operator, le));
-  variable_matrix.def(Variable() <= py::self, DOC(sleipnir, operator, le));
-  variable_matrix.def(double() <= py::self, DOC(sleipnir, operator, le));
-  variable_matrix.def(int() <= py::self, DOC(sleipnir, operator, le));
-  variable_matrix.def(py::self > py::self, DOC(sleipnir, operator, gt));
-  variable_matrix.def(py::self > Variable(), DOC(sleipnir, operator, gt));
-  variable_matrix.def(py::self > double(), DOC(sleipnir, operator, gt));
-  variable_matrix.def(py::self > int(), DOC(sleipnir, operator, gt));
-  variable_matrix.def(Variable() > py::self, DOC(sleipnir, operator, gt));
-  variable_matrix.def(double() > py::self, DOC(sleipnir, operator, gt));
-  variable_matrix.def(int() > py::self, DOC(sleipnir, operator, gt));
-  variable_matrix.def(py::self >= py::self, DOC(sleipnir, operator, ge));
-  variable_matrix.def(py::self >= Variable(), DOC(sleipnir, operator, ge));
-  variable_matrix.def(py::self >= double(), DOC(sleipnir, operator, ge));
-  variable_matrix.def(py::self >= int(), DOC(sleipnir, operator, ge));
-  variable_matrix.def(Variable() >= py::self, DOC(sleipnir, operator, ge));
-  variable_matrix.def(double() >= py::self, DOC(sleipnir, operator, ge));
-  variable_matrix.def(int() >= py::self, DOC(sleipnir, operator, ge));
+      "func"_a, DOC(sleipnir, VariableMatrix, CwiseTransform));
+  cls.def_static("zero", &VariableMatrix::Zero, "rows"_a, "cols"_a,
+                 DOC(sleipnir, VariableMatrix, Zero));
+  cls.def_static("ones", &VariableMatrix::Ones, "rows"_a, "cols"_a,
+                 DOC(sleipnir, VariableMatrix, Ones));
+  cls.def(py::self == py::self, "rhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(py::self == Variable(), "rhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(py::self == double(), "rhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(py::self == int(), "rhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(Variable() == py::self, "lhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(double() == py::self, "lhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(int() == py::self, "lhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(py::self < py::self, "rhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(py::self < Variable(), "rhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(py::self < double(), "rhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(py::self < int(), "rhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(Variable() < py::self, "lhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(double() < py::self, "lhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(int() < py::self, "lhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(py::self <= py::self, "rhs"_a, DOC(sleipnir, operator, le));
+  cls.def(py::self <= Variable(), "rhs"_a, DOC(sleipnir, operator, le));
+  cls.def(py::self <= double(), "rhs"_a, DOC(sleipnir, operator, le));
+  cls.def(py::self <= int(), "rhs"_a, DOC(sleipnir, operator, le));
+  cls.def(Variable() <= py::self, "lhs"_a, DOC(sleipnir, operator, le));
+  cls.def(double() <= py::self, "lhs"_a, DOC(sleipnir, operator, le));
+  cls.def(int() <= py::self, "lhs"_a, DOC(sleipnir, operator, le));
+  cls.def(py::self > py::self, "rhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(py::self > Variable(), "rhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(py::self > double(), "rhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(py::self > int(), "rhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(Variable() > py::self, "lhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(double() > py::self, "lhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(int() > py::self, "lhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(py::self >= py::self, "rhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(py::self >= Variable(), "rhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(py::self >= double(), "rhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(py::self >= int(), "rhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(Variable() >= py::self, "lhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(double() >= py::self, "lhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(int() >= py::self, "lhs"_a, DOC(sleipnir, operator, ge));
   py::implicitly_convertible<VariableMatrix, Variable>();
 
   // VariableMatrix-VariableBlock overloads
-  variable_matrix.def(
+  cls.def(
       "__mul__",
       [](const VariableMatrix& lhs, const VariableBlock<VariableMatrix>& rhs) {
         return lhs * rhs;
       },
-      py::is_operator());
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(
       "__rmul__",
       [](const VariableMatrix& rhs, const VariableBlock<VariableMatrix>& lhs) {
         return lhs * rhs;
       },
-      py::is_operator());
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(
       "__matmul__",
       [](const VariableMatrix& lhs, const VariableBlock<VariableMatrix>& rhs) {
         return lhs * rhs;
       },
-      py::is_operator());
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(
       "__add__",
       [](const VariableMatrix& lhs, const VariableBlock<VariableMatrix>& rhs) {
         return lhs + rhs;
       },
-      py::is_operator());
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(
       "__sub__",
       [](const VariableMatrix& lhs, const VariableBlock<VariableMatrix>& rhs) {
         return lhs - rhs;
       },
-      py::is_operator());
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a);
+  cls.def(
       "__eq__",
       [](const VariableMatrix& lhs, const VariableBlock<VariableMatrix>& rhs) {
         return lhs == rhs;
       },
-      py::is_operator(), DOC(sleipnir, operator, eq));
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(
       "__lt__",
       [](const VariableMatrix& lhs, const VariableBlock<VariableMatrix>& rhs) {
         return lhs < rhs;
       },
-      py::is_operator(), DOC(sleipnir, operator, lt));
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(
       "__le__",
       [](const VariableMatrix& lhs, const VariableBlock<VariableMatrix>& rhs) {
         return lhs <= rhs;
       },
-      py::is_operator(), DOC(sleipnir, operator, le));
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a, DOC(sleipnir, operator, le));
+  cls.def(
       "__gt__",
       [](const VariableMatrix& lhs, const VariableBlock<VariableMatrix>& rhs) {
         return lhs > rhs;
       },
-      py::is_operator(), DOC(sleipnir, operator, gt));
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(
       "__ge__",
       [](const VariableMatrix& lhs, const VariableBlock<VariableMatrix>& rhs) {
         return lhs >= rhs;
       },
-      py::is_operator(), DOC(sleipnir, operator, ge));
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a, DOC(sleipnir, operator, ge));
+  cls.def(
       "__eq__",
       [](const VariableMatrix& lhs, const py::array_t<double>& rhs) {
         return lhs == rhs.cast<Eigen::MatrixXd>();
       },
-      py::is_operator(), DOC(sleipnir, operator, eq));
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a, DOC(sleipnir, operator, eq));
+  cls.def(
       "__lt__",
       [](const VariableMatrix& lhs, const py::array_t<double>& rhs) {
         return lhs < rhs.cast<Eigen::MatrixXd>();
       },
-      py::is_operator(), DOC(sleipnir, operator, lt));
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a, DOC(sleipnir, operator, lt));
+  cls.def(
       "__le__",
       [](const VariableMatrix& lhs, const py::array_t<double>& rhs) {
         return lhs <= rhs.cast<Eigen::MatrixXd>();
       },
-      py::is_operator(), DOC(sleipnir, operator, le));
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a, DOC(sleipnir, operator, le));
+  cls.def(
       "__gt__",
       [](const VariableMatrix& lhs, const py::array_t<double>& rhs) {
         return lhs > rhs.cast<Eigen::MatrixXd>();
       },
-      py::is_operator(), DOC(sleipnir, operator, gt));
-  variable_matrix.def(
+      py::is_operator(), "rhs"_a, DOC(sleipnir, operator, gt));
+  cls.def(
       "__ge__",
       [](const VariableMatrix& lhs, const py::array_t<double>& rhs) {
         return lhs >= rhs.cast<Eigen::MatrixXd>();
       },
-      py::is_operator(), DOC(sleipnir, operator, ge));
+      py::is_operator(), "rhs"_a, DOC(sleipnir, operator, ge));
 
-  variable_matrix.def("__len__", &VariableMatrix::Rows,
-                      DOC(sleipnir, VariableMatrix, Rows));
+  cls.def("__len__", &VariableMatrix::Rows,
+          DOC(sleipnir, VariableMatrix, Rows));
 
-  variable_matrix.def(
+  cls.def(
       "__iter__",
       [](const VariableMatrix& self) {
         return py::make_iterator(self.begin(), self.end());
@@ -558,13 +563,13 @@ void BindVariableMatrix(py::module_& autodiff) {
          const std::function<Variable(const Variable&, const Variable&)> func) {
         return CwiseReduce(lhs, rhs, func);
       },
-      DOC(sleipnir, CwiseReduce));
+      "lhs"_a, "rhs"_a, "func"_a, DOC(sleipnir, CwiseReduce));
 
   autodiff.def(
       "block",
       static_cast<VariableMatrix (*)(std::vector<std::vector<VariableMatrix>>)>(
           &Block),
-      DOC(sleipnir, Block));
+      "list"_a, DOC(sleipnir, Block));
 }  // NOLINT(readability/fn_size)
 
 }  // namespace sleipnir
diff --git a/jormungandr/cpp/optimization/BindConstraints.cpp b/jormungandr/cpp/optimization/BindConstraints.cpp
deleted file mode 100644
index 3837bc9b..00000000
--- a/jormungandr/cpp/optimization/BindConstraints.cpp
+++ /dev/null
@@ -1,29 +0,0 @@
-// Copyright (c) Sleipnir contributors
-
-#include <pybind11/operators.h>
-#include <pybind11/pybind11.h>
-#include <sleipnir/optimization/Constraints.hpp>
-
-#include "Docstrings.hpp"
-
-namespace py = pybind11;
-
-namespace sleipnir {
-
-void BindConstraints(py::module_& optimization) {
-  py::class_<EqualityConstraints> equalityConstraints{
-      optimization, "EqualityConstraints", DOC(sleipnir, EqualityConstraints)};
-  equalityConstraints.def(
-      "__bool__", [](const EqualityConstraints& self) -> bool { return self; },
-      py::is_operator(), DOC(sleipnir, EqualityConstraints, operator, bool));
-
-  py::class_<InequalityConstraints> inequalityConstraints{
-      optimization, "InequalityConstraints",
-      DOC(sleipnir_InequalityConstraints)};
-  inequalityConstraints.def(
-      "__bool__",
-      [](const InequalityConstraints& self) -> bool { return self; },
-      py::is_operator(), DOC(sleipnir, InequalityConstraints, operator, bool));
-}
-
-}  // namespace sleipnir
diff --git a/jormungandr/cpp/optimization/BindEqualityConstraints.cpp b/jormungandr/cpp/optimization/BindEqualityConstraints.cpp
new file mode 100644
index 00000000..e4052f71
--- /dev/null
+++ b/jormungandr/cpp/optimization/BindEqualityConstraints.cpp
@@ -0,0 +1,19 @@
+// Copyright (c) Sleipnir contributors
+
+#include <pybind11/operators.h>
+#include <pybind11/pybind11.h>
+#include <sleipnir/optimization/Constraints.hpp>
+
+#include "Docstrings.hpp"
+
+namespace py = pybind11;
+
+namespace sleipnir {
+
+void BindEqualityConstraints(py::class_<EqualityConstraints>& cls) {
+  cls.def(
+      "__bool__", [](const EqualityConstraints& self) -> bool { return self; },
+      py::is_operator(), DOC(sleipnir, EqualityConstraints, operator, bool));
+}
+
+}  // namespace sleipnir
diff --git a/jormungandr/cpp/optimization/BindInequalityConstraints.cpp b/jormungandr/cpp/optimization/BindInequalityConstraints.cpp
new file mode 100644
index 00000000..f23a338d
--- /dev/null
+++ b/jormungandr/cpp/optimization/BindInequalityConstraints.cpp
@@ -0,0 +1,20 @@
+// Copyright (c) Sleipnir contributors
+
+#include <pybind11/operators.h>
+#include <pybind11/pybind11.h>
+#include <sleipnir/optimization/Constraints.hpp>
+
+#include "Docstrings.hpp"
+
+namespace py = pybind11;
+
+namespace sleipnir {
+
+void BindInequalityConstraints(py::class_<InequalityConstraints>& cls) {
+  cls.def(
+      "__bool__",
+      [](const InequalityConstraints& self) -> bool { return self; },
+      py::is_operator(), DOC(sleipnir, InequalityConstraints, operator, bool));
+}
+
+}  // namespace sleipnir
diff --git a/jormungandr/cpp/optimization/BindOptimizationProblem.cpp b/jormungandr/cpp/optimization/BindOptimizationProblem.cpp
index 81875a6e..a22207dd 100644
--- a/jormungandr/cpp/optimization/BindOptimizationProblem.cpp
+++ b/jormungandr/cpp/optimization/BindOptimizationProblem.cpp
@@ -14,11 +14,9 @@ namespace py = pybind11;
 
 namespace sleipnir {
 
-void BindOptimizationProblem(py::module_& optimization) {
+void BindOptimizationProblem(py::class_<OptimizationProblem>& cls) {
   using namespace pybind11::literals;
 
-  py::class_<OptimizationProblem> cls{optimization, "OptimizationProblem",
-                                      DOC(sleipnir, OptimizationProblem)};
   cls.def(py::init<>(),
           DOC(sleipnir, OptimizationProblem, OptimizationProblem));
   cls.def("decision_variable",
@@ -33,40 +31,40 @@ void BindOptimizationProblem(py::module_& optimization) {
           DOC(sleipnir, OptimizationProblem, SymmetricDecisionVariable));
   cls.def("minimize",
           py::overload_cast<const Variable&>(&OptimizationProblem::Minimize),
-          DOC(sleipnir, OptimizationProblem, Minimize));
+          "cost"_a, DOC(sleipnir, OptimizationProblem, Minimize));
   cls.def(
       "minimize",
       [](OptimizationProblem& self, const VariableMatrix& cost) {
         self.Minimize(cost);
       },
-      DOC(sleipnir, OptimizationProblem, Minimize));
+      "cost"_a, DOC(sleipnir, OptimizationProblem, Minimize));
   cls.def(
       "minimize",
       [](OptimizationProblem& self, double cost) { self.Minimize(cost); },
-      DOC(sleipnir, OptimizationProblem, Minimize));
+      "cost"_a, DOC(sleipnir, OptimizationProblem, Minimize));
   cls.def("maximize",
           py::overload_cast<const Variable&>(&OptimizationProblem::Maximize),
-          DOC(sleipnir, OptimizationProblem, Maximize));
+          "objective"_a, DOC(sleipnir, OptimizationProblem, Maximize));
   cls.def(
       "maximize",
       [](OptimizationProblem& self, const VariableMatrix& objective) {
         self.Maximize(objective);
       },
-      DOC(sleipnir, OptimizationProblem, Maximize));
+      "objective"_a, DOC(sleipnir, OptimizationProblem, Maximize));
   cls.def(
       "maximize",
       [](OptimizationProblem& self, double objective) {
         self.Maximize(objective);
       },
-      DOC(sleipnir, OptimizationProblem, Maximize));
+      "objective"_a, DOC(sleipnir, OptimizationProblem, Maximize));
   cls.def("subject_to",
           py::overload_cast<const EqualityConstraints&>(
               &OptimizationProblem::SubjectTo),
-          DOC(sleipnir, OptimizationProblem, SubjectTo));
+          "constraint"_a, DOC(sleipnir, OptimizationProblem, SubjectTo));
   cls.def("subject_to",
           py::overload_cast<const InequalityConstraints&>(
               &OptimizationProblem::SubjectTo),
-          DOC(sleipnir, OptimizationProblem, SubjectTo, 3));
+          "constraint"_a, DOC(sleipnir, OptimizationProblem, SubjectTo, 3));
   cls.def(
       "solve",
       [](OptimizationProblem& self, const py::kwargs& kwargs) {
@@ -99,7 +97,7 @@ void BindOptimizationProblem(py::module_& optimization) {
          std::function<bool(const SolverIterationInfo&)> callback) {
         self.Callback(std::move(callback));
       },
-      DOC(sleipnir, OptimizationProblem, Callback, 2));
+      "callback"_a, DOC(sleipnir, OptimizationProblem, Callback, 2));
 }
 
 }  // namespace sleipnir
diff --git a/jormungandr/cpp/optimization/BindSolverExitCondition.cpp b/jormungandr/cpp/optimization/BindSolverExitCondition.cpp
index 00ea1a15..90bec286 100644
--- a/jormungandr/cpp/optimization/BindSolverExitCondition.cpp
+++ b/jormungandr/cpp/optimization/BindSolverExitCondition.cpp
@@ -9,9 +9,7 @@ namespace py = pybind11;
 
 namespace sleipnir {
 
-void BindSolverExitCondition(py::module_& optimization) {
-  py::enum_<SolverExitCondition> e{optimization, "SolverExitCondition",
-                                   DOC(sleipnir, SolverExitCondition)};
+void BindSolverExitCondition(py::enum_<SolverExitCondition>& e) {
   e.value("SUCCESS", SolverExitCondition::kSuccess,
           DOC(sleipnir, SolverExitCondition, kSuccess));
   e.value("SOLVED_TO_ACCEPTABLE_TOLERANCE",
diff --git a/jormungandr/cpp/optimization/BindSolverIterationInfo.cpp b/jormungandr/cpp/optimization/BindSolverIterationInfo.cpp
index 671a9422..7f3fb495 100644
--- a/jormungandr/cpp/optimization/BindSolverIterationInfo.cpp
+++ b/jormungandr/cpp/optimization/BindSolverIterationInfo.cpp
@@ -10,9 +10,7 @@ namespace py = pybind11;
 
 namespace sleipnir {
 
-void BindSolverIterationInfo(py::module_& optimization) {
-  py::class_<SolverIterationInfo> cls{optimization, "SolverIterationInfo",
-                                      DOC(sleipnir, SolverIterationInfo)};
+void BindSolverIterationInfo(py::class_<SolverIterationInfo>& cls) {
   cls.def_readonly("iteration", &SolverIterationInfo::iteration,
                    DOC(sleipnir, SolverIterationInfo, iteration));
   cls.def_property_readonly(
diff --git a/jormungandr/cpp/optimization/BindSolverStatus.cpp b/jormungandr/cpp/optimization/BindSolverStatus.cpp
index 547f8641..1b09dfdb 100644
--- a/jormungandr/cpp/optimization/BindSolverStatus.cpp
+++ b/jormungandr/cpp/optimization/BindSolverStatus.cpp
@@ -9,11 +9,9 @@ namespace py = pybind11;
 
 namespace sleipnir {
 
-void BindSolverStatus(py::module_& optimization) {
+void BindSolverStatus(py::class_<SolverStatus>& cls) {
   using namespace pybind11::literals;
 
-  py::class_<SolverStatus> cls{optimization, "SolverStatus",
-                               DOC(sleipnir, SolverStatus)};
   cls.def(py::init<>());
   cls.def(py::init<ExpressionType, ExpressionType, ExpressionType,
                    SolverExitCondition, double>(),
diff --git a/pyproject.toml b/pyproject.toml
index 141638db..3348f35a 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -26,6 +26,7 @@ include = [
   "SleipnirConfig.cmake.in",
   "cmake/modules/*.cmake",
   "cmake/*.patch",
+  "cmake/*.py",
   "jormungandr/*",
   "include/*",
   "src/*"
diff --git a/tools/pydoc.py b/tools/pydoc.py
index c3ad7828..a7472b0c 100755
--- a/tools/pydoc.py
+++ b/tools/pydoc.py
@@ -10,6 +10,15 @@
 def main():
     # Clear workspace
     shutil.rmtree(".py-build-cmake_cache", ignore_errors=True)
+    subprocess.run(
+        [
+            "git",
+            "restore",
+            "jormungandr/autodiff/__init__.py",
+            "jormungandr/optimization/__init__.py",
+        ],
+        check=True,
+    )
 
     # Generate .pyi files
     subprocess.run(
@@ -28,40 +37,27 @@ def main():
     for package in ["autodiff", "optimization"]:
         # Read .pyi
         with open(os.path.join(PYI_PATH, package + ".pyi")) as f:
-            pyi_content = f.read()
+            package_content = f.read()
 
-        # Fix up contents
+        # Remove redundant prefixes for documentation
         if package == "autodiff":
-            pyi_content = pyi_content.replace("sleipnir::", "").replace(
-                "VariableBlock<VariableMatrix>", "VariableBlock"
-            )
+            package_content = package_content.replace(
+                "import _jormungandr.optimization\n", ""
+            ).replace("_jormungandr.optimization.", "")
         elif package == "optimization":
-            pyi_content = (
-                pyi_content.replace("_jormungandr.autodiff.", "")
-                .replace("import _jormungandr.autodiff\n", "")
-                .replace("<ExpressionType.NONE: 0>", "ExpressionType.NONE")
-                .replace(
-                    "<SolverExitCondition.SUCCESS: 0>", "SolverExitCondition.SUCCESS"
-                )
-            )
+            package_content = package_content.replace(
+                "import _jormungandr.autodiff\n", ""
+            ).replace("_jormungandr.autodiff.", "")
 
         # Replace _jormungandr.package import with contents
-        with (
-            open(os.path.join("jormungandr", package, "__init__.py")) as input,
-            open(
-                os.path.join("jormungandr", package, "__init__.py.new"), mode="w"
-            ) as output,
-        ):
-            package_content = input.read()
-            output.write(
-                package_content.replace(
-                    f"from .._jormungandr.{package} import *", pyi_content
+        with open(os.path.join("jormungandr", package, "__init__.py")) as f:
+            init_content = f.read()
+        with open(os.path.join("jormungandr", package, "__init__.py"), mode="w") as f:
+            f.write(
+                init_content.replace(
+                    f"from .._jormungandr.{package} import *", package_content
                 )
             )
-        shutil.move(
-            os.path.join("jormungandr", package, "__init__.py.new"),
-            os.path.join("jormungandr", package, "__init__.py"),
-        )
 
 
 if __name__ == "__main__":