Add h_sparse_triplets method for sparse Hessian

include/hypergraph/hypergraph.h

@@ -20,6 +20,7 @@
 #include <cstddef>
 #include <memory>
 #include <stdexcept>
+#include <tuple>
 
 namespace hypergraph {
 
@@ -654,6 +655,84 @@ class HyperGraph {
         }
     }
 
+    index num_variables() const
+    {
+        return length(m_variables);
+    }
+
+    std::tuple<std::vector<index>, std::vector<index>, std::vector<T>>
+    h_sparse_triplets(const bool full = false) const
+    {
+        const index n = length(m_variables);
+
+        std::vector<index> rows;
+        std::vector<index> cols;
+        std::vector<T> values;
+
+        // Build a map from vertex_id → variable_index for efficient lookup
+        tsl::robin_map<index, index> vid_to_var;
+        vid_to_var.reserve(n);
+        for (index i = 0; i < n; i++) {
+            vid_to_var[m_variables[i].id()] = i;
+        }
+
+        // Diagonal entries from m_self_second_order_edges
+        for (index i = 0; i < n; i++) {
+            const T val = m_self_second_order_edges[m_variables[i].id()];
+            if (val != T{}) {
+                rows.push_back(i);
+                cols.push_back(i);
+                values.push_back(double(val));
+            }
+        }
+
+        // Off-diagonal entries from m_second_order_edges
+        // The structure stores entries with key = min(vid_a, vid_b) inside
+        // m_second_order_edges[max(vid_a, vid_b)].
+        for (index outer_vid = 0; outer_vid < length(m_second_order_edges); outer_vid++) {
+            const auto& btree = m_second_order_edges[outer_vid];
+            if (btree.empty()) {
+                continue;
+            }
+
+            // Check if the outer vertex_id corresponds to a variable
+            auto outer_it = vid_to_var.find(outer_vid);
+            if (outer_it == vid_to_var.end()) {
+                continue;
+            }
+            const index outer_var_idx = outer_it->second;
+
+            for (const auto& [inner_vid, val] : btree) {
+                if (val == T{}) {
+                    continue;
+                }
+
+                auto inner_it = vid_to_var.find(inner_vid);
+                if (inner_it == vid_to_var.end()) {
+                    continue;
+                }
+                const index inner_var_idx = inner_it->second;
+
+                // inner_vid < outer_vid by construction (minmax in second_order_edge)
+                // So inner_var_idx may or may not be < outer_var_idx
+                const auto [row, col] = hypergraph::minmax(inner_var_idx, outer_var_idx);
+
+                // Upper triangle: row <= col
+                rows.push_back(row);
+                cols.push_back(col);
+                values.push_back(double(val));
+
+                if (full && row != col) {
+                    rows.push_back(col);
+                    cols.push_back(row);
+                    values.push_back(double(val));
+                }
+            }
+        }
+
+        return {std::move(rows), std::move(cols), std::move(values)};
+    }
+
 public: // python
     template <typename TModule>
     static void register_python(TModule& m)
@@ -666,14 +745,46 @@ class HyperGraph {
         nb::class_<Type>(m, name.c_str())
             // constructors
             .def(nb::init<>())
+            // properties
+            .def_prop_ro("num_variables", &Type::num_variables)
             // methods
             .def("new_variable", &Type::new_variable, "value"_a)
             .def("new_variables", &Type::new_variables, "values"_a)
             .def("compute", &Type::compute, "expression"_a)
             .def("g", nb::overload_cast<>(&Type::g, nb::const_))
             .def("g", nb::overload_cast<Eigen::Ref<Eigen::VectorXd>>(&Type::g, nb::const_), "out"_a)
             .def("h", nb::overload_cast<const bool>(&Type::h, nb::const_), "full"_a = false)
-            .def("h", nb::overload_cast<Eigen::Ref<Eigen::MatrixXd>, const bool>(&Type::h, nb::const_), "out"_a, "full"_a = false);
+            .def("h", nb::overload_cast<Eigen::Ref<Eigen::MatrixXd>, const bool>(&Type::h, nb::const_), "out"_a, "full"_a = false)
+            .def("h_sparse_triplets", &Type::h_sparse_triplets, "full"_a = false)
+            .def("h_sparse", [](const Type& self, const std::string& format, const bool full) -> nb::object {
+                auto [rows, cols, values] = self.h_sparse_triplets(full);
+                const auto n = self.num_variables();
+
+                nb::module_ scipy_sparse = nb::module_::import_("scipy.sparse");
+
+                // Build numpy arrays from the triplet vectors
+                nb::module_ np = nb::module_::import_("numpy");
+                nb::object np_array = np.attr("array");
+
+                nb::object py_rows = np_array(nb::cast(rows), "dtype"_a = np.attr("int64"));
+                nb::object py_cols = np_array(nb::cast(cols), "dtype"_a = np.attr("int64"));
+                nb::object py_vals = np_array(nb::cast(values), "dtype"_a = np.attr("float64"));
+
+                nb::tuple data_ij = nb::make_tuple(py_vals, nb::make_tuple(py_rows, py_cols));
+                nb::tuple shape = nb::make_tuple(n, n);
+
+                nb::object coo = scipy_sparse.attr("coo_matrix")(data_ij, "shape"_a = shape);
+
+                if (format == "coo") {
+                    return coo;
+                } else if (format == "csc") {
+                    return coo.attr("tocsc")();
+                } else if (format == "csr") {
+                    return coo.attr("tocsr")();
+                } else {
+                    throw std::invalid_argument("h_sparse: format must be 'coo', 'csc', or 'csr'");
+                }
+            }, "format"_a = "csc", "full"_a = false);
     }
 };
 

pyproject.toml

@@ -22,6 +22,7 @@ classifiers = [
 ]
 dependencies = [
     "numpy",
+    "scipy",
 ]
 
 [project.urls]

src/Module.cpp

@@ -2,6 +2,7 @@
 #include <nanobind/eigen/dense.h>
 #include <nanobind/operators.h>
 #include <nanobind/stl/string.h>
+#include <nanobind/stl/tuple.h>
 #include <nanobind/stl/vector.h>
 
 #include <hypergraph/hypergraph.h>

tests/TestHyperGraph.py

@@ -2,6 +2,7 @@
 import math
 import hypergraph as hg
 import numpy as np
+import scipy.sparse
 from numpy.testing import assert_equal, assert_array_equal, assert_almost_equal, assert_array_almost_equal
 
 
@@ -1154,5 +1155,204 @@ def test_hypot_pythagorean(self):
         assert_almost_equal(result.value, 13)
 
 
+    # sparse Hessian
+
+    def test_h_sparse_triplets_basic(self):
+        """h_sparse_triplets returns correct (rows, cols, values) for upper triangle"""
+        graph = hg.HyperGraph()
+
+        a, b = graph.new_variables([5, 6])
+
+        result = a * b
+        graph.compute(result)
+
+        rows, cols, values = graph.h_sparse_triplets()
+        # a*b has H = [[0, 1], [0, 0]] (upper triangle only)
+        # Only one non-zero entry: (0, 1) = 1
+        assert_equal(len(rows), 1)
+        assert_equal(len(cols), 1)
+        assert_equal(len(values), 1)
+        assert_equal(rows[0], 0)
+        assert_equal(cols[0], 1)
+        assert_almost_equal(values[0], 1.0)
+
+    def test_h_sparse_triplets_full(self):
+        """h_sparse_triplets with full=True returns both triangles"""
+        graph = hg.HyperGraph()
+
+        a, b = graph.new_variables([5, 6])
+
+        result = a * b
+        graph.compute(result)
+
+        rows, cols, values = graph.h_sparse_triplets(full=True)
+        # full=True: (0,1)=1 and (1,0)=1
+        assert_equal(len(rows), 2)
+        # Convert to a set of (row, col, val) for order-independent check
+        entries = set(zip(rows, cols, values))
+        self.assertIn((0, 1, 1.0), entries)
+        self.assertIn((1, 0, 1.0), entries)
+
+    def test_h_sparse_triplets_diagonal(self):
+        """h_sparse_triplets includes diagonal entries"""
+        graph = hg.HyperGraph()
+
+        a, b = graph.new_variables([3, 5])
+
+        result = a.square()
+        graph.compute(result)
+
+        rows, cols, values = graph.h_sparse_triplets()
+        # d²(x²)/dx² = 2, only diagonal entry (0,0)
+        assert_equal(len(rows), 1)
+        assert_equal(rows[0], 0)
+        assert_equal(cols[0], 0)
+        assert_almost_equal(values[0], 2.0)
+
+    def test_h_sparse_triplets_empty(self):
+        """h_sparse_triplets returns empty lists when Hessian is zero"""
+        graph = hg.HyperGraph()
+
+        a, b = graph.new_variables([5, 6])
+
+        result = a + b
+        graph.compute(result)
+
+        rows, cols, values = graph.h_sparse_triplets()
+        assert_equal(len(rows), 0)
+        assert_equal(len(cols), 0)
+        assert_equal(len(values), 0)
+
+    def test_h_sparse_matches_dense(self):
+        """h_sparse().toarray() matches h(full=True) for a complex expression"""
+        graph = hg.HyperGraph()
+
+        ax, ay, az, bx, by, bz = graph.new_variables([1, 2, 3, 4, 5, 6])
+
+        a = np.array([ax, ay, az])
+        b = np.array([bx, by, bz])
+
+        result = np.linalg.norm(np.cross(a, b))
+
+        graph.compute(result)
+
+        h_dense = graph.h(full=True)
+        h_sparse = graph.h_sparse(full=True)
+
+        assert_array_almost_equal(h_sparse.toarray(), h_dense)
+
+    def test_h_sparse_upper_triangle_matches(self):
+        """h_sparse() default (upper triangle) matches h() upper triangle"""
+        graph = hg.HyperGraph()
+
+        ax, ay, az, bx, by, bz = graph.new_variables([1, 2, 3, 4, 5, 6])
+
+        a = np.array([ax, ay, az])
+        b = np.array([bx, by, bz])
+
+        result = np.linalg.norm(np.cross(a, b))
+
+        graph.compute(result)
+
+        h_dense_upper = graph.h(full=False)
+        h_sparse_upper = graph.h_sparse(full=False)
+
+        # The sparse version should have the upper triangle entries
+        # When converted to dense, lower triangle should be zero
+        sparse_dense = h_sparse_upper.toarray()
+        assert_array_almost_equal(sparse_dense, h_dense_upper)
+
+    def test_h_sparse_coo_format(self):
+        """h_sparse(format='coo') returns COO matrix"""
+        graph = hg.HyperGraph()
+
+        a, b = graph.new_variables([5, 6])
+        result = a * b
+        graph.compute(result)
+
+        h_sparse = graph.h_sparse(format='coo')
+        self.assertIsInstance(h_sparse, scipy.sparse.coo_matrix)
+
+    def test_h_sparse_csc_format(self):
+        """h_sparse(format='csc') returns CSC matrix (default)"""
+        graph = hg.HyperGraph()
+
+        a, b = graph.new_variables([5, 6])
+        result = a * b
+        graph.compute(result)
+
+        h_sparse = graph.h_sparse(format='csc')
+        self.assertIsInstance(h_sparse, scipy.sparse.csc_matrix)
+
+    def test_h_sparse_csr_format(self):
+        """h_sparse(format='csr') returns CSR matrix"""
+        graph = hg.HyperGraph()
+
+        a, b = graph.new_variables([5, 6])
+        result = a * b
+        graph.compute(result)
+
+        h_sparse = graph.h_sparse(format='csr')
+        self.assertIsInstance(h_sparse, scipy.sparse.csr_matrix)
+
+    def test_h_sparse_default_is_csc(self):
+        """h_sparse() default format is CSC"""
+        graph = hg.HyperGraph()
+
+        a, b = graph.new_variables([5, 6])
+        result = a * b
+        graph.compute(result)
+
+        h_sparse = graph.h_sparse()
+        self.assertIsInstance(h_sparse, scipy.sparse.csc_matrix)
+
+    def test_h_sparse_invalid_format_raises(self):
+        """h_sparse with invalid format raises ValueError"""
+        graph = hg.HyperGraph()
+
+        a, b = graph.new_variables([5, 6])
+        result = a * b
+        graph.compute(result)
+
+        with self.assertRaises(Exception):
+            graph.h_sparse(format='invalid')
+
+    def test_h_sparse_shape(self):
+        """h_sparse returns matrix with correct shape"""
+        graph = hg.HyperGraph()
+
+        a, b, c = graph.new_variables([1, 2, 3])
+        result = a * b + b * c
+        graph.compute(result)
+
+        h_sparse = graph.h_sparse()
+        assert_equal(h_sparse.shape, (3, 3))
+
+    def test_h_sparse_multiplication(self):
+        """h_sparse values match dense h for multiplication"""
+        graph = hg.HyperGraph()
+
+        a, b = graph.new_variables([5, 6])
+        result = a * b
+        graph.compute(result)
+
+        h_dense = graph.h(full=True)
+        h_sparse = graph.h_sparse(full=True)
+
+        assert_array_almost_equal(h_sparse.toarray(), h_dense)
+
+    def test_num_variables(self):
+        """num_variables property returns correct count"""
+        graph = hg.HyperGraph()
+
+        assert_equal(graph.num_variables, 0)
+
+        graph.new_variable(1)
+        assert_equal(graph.num_variables, 1)
+
+        graph.new_variables([2, 3, 4])
+        assert_equal(graph.num_variables, 4)
+
+
 if __name__ == '__main__':
     unittest.main()