Implement spin_op::to_sparse_matrix().

python/runtime/cudaq/spin/py_matrix.cpp

@@ -52,8 +52,18 @@ void bindComplexMatrix(py::module &mod) {
            }),
            "Create a :class:`ComplexMatrix` from a buffer of data, such as a "
            "numpy.ndarray.")
-      .def("__getitem__", &complex_matrix::operator(),
-           "Return the matrix element at i, j.")
+      .def(
+          "__getitem__",
+          [](complex_matrix &m, std::size_t i, std::size_t j) {
+            return m(i, j);
+          },
+          "Return the matrix element at i, j.")
+      .def(
+          "__getitem__",
+          [](complex_matrix &m, std::tuple<std::size_t, std::size_t> rowCol) {
+            return m(std::get<0>(rowCol), std::get<1>(rowCol));
+          },
+          "Return the matrix element at i, j.")
       .def("minimal_eigenvalue", &complex_matrix::minimal_eigenvalue,
            "Return the lowest eigenvalue for this :class:`ComplexMatrix`.")
       .def(

python/runtime/cudaq/spin/py_spin_op.cpp

@@ -139,7 +139,11 @@ void bindSpinOperator(py::module &mod) {
            "represented as a double.")
       .def("to_matrix", &spin_op::to_matrix,
            "Return `self` as a :class:`ComplexMatrix`.")
-
+      .def("to_sparse_matrix", &spin_op::to_sparse_matrix,
+           "Return `self` as a sparse matrix representation. This "
+           "representation is a `Tuple[list[complex],list[int], list[int]]`, "
+           "encoding the non-zero values, rows, and columns of the matrix. "
+           "This format is supported by `scipy.sparse.csr_array`.")
       .def(
           "__iter__",
           [](spin_op &self) {

python/tests/unittests/test_SpinOperator.py

@@ -16,7 +16,6 @@
 from cudaq import spin
 import numpy as np
 
-
 def assert_close(want, got, tolerance=1.e-5) -> bool:
     return abs(want - got) < tolerance
 
@@ -335,6 +334,26 @@ def test_spin_op_iter():
         count += 1
     assert count == 5
 
+def test_spin_op_sparse_matrix():
+    """
+    Test that the `cudaq.SpinOperator` can produce its sparse matrix representation 
+    and that we can use that matrix with standard python packages like numpy.
+    """
+    hamiltonian = 5.907 - 2.1433 * spin.x(0) * spin.x(1) - 2.1433 * spin.y(
+        0) * spin.y(1) + .21829 * spin.z(0) - 6.125 * spin.z(1)
+    numQubits = hamiltonian.get_qubit_count()
+    mat = hamiltonian.to_matrix()
+    data, rows, cols = hamiltonian.to_sparse_matrix()
+    for i, value in enumerate(data):
+        print(rows[i],cols[i],value)
+        assert np.isclose(mat[rows[i],cols[i]], value)
+
+    # can use scipy 
+    # scipyM = scipy.sparse.csr_array((data, (rows, cols)), shape=(2**numQubits,2**numQubits))
+    # E, ev = scipy.sparse.linalg.eigsh(scipyM, k=1, which='SA')
+    # assert np.isclose(E[0], -1.7488, 1e-2)
+
+
 
 def test_spin_op_from_word():
     s = cudaq.SpinOperator.from_word("ZZZ")

runtime/common/EigenSparse.h

@@ -0,0 +1,25 @@
+/****************************************************************-*- C++ -*-****
+ * Copyright (c) 2022 - 2023 NVIDIA Corporation & Affiliates.                  *
+ * All rights reserved.                                                        *
+ *                                                                             *
+ * This source code and the accompanying materials are made available under    *
+ * the terms of the Apache License 2.0 which accompanies this distribution.    *
+ ******************************************************************************/
+
+#ifdef __clang__
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wdeprecated-anon-enum-enum-conversion"
+#pragma clang diagnostic ignored "-Wcovered-switch-default"
+#endif
+#if (defined(__GNUC__) && !defined(__clang__) && !defined(__INTEL_COMPILER))
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wdeprecated-enum-enum-conversion"
+#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
+#endif
+#include <Eigen/Sparse>
+#if (defined(__GNUC__) && !defined(__clang__) && !defined(__INTEL_COMPILER))
+#pragma GCC diagnostic pop
+#endif
+#ifdef __clang__
+#pragma clang diagnostic pop
+#endif

runtime/cudaq/spin/spin_op.cpp

@@ -7,13 +7,14 @@
  ******************************************************************************/
 
 #include "common/EigenDense.h"
+#include "common/EigenSparse.h"
 #include "common/FmtCore.h"
 #include <cudaq/spin_op.h>
 #include <stdint.h>
+#include <unsupported/Eigen/KroneckerProduct>
 #ifdef CUDAQ_HAS_OPENMP
 #include <omp.h>
 #endif
-
 #include <algorithm>
 #include <array>
 #include <cassert>
@@ -199,6 +200,59 @@ complex_matrix spin_op::to_matrix() const {
   return A;
 }
 
+spin_op::csr_spmatrix spin_op::to_sparse_matrix() const {
+  auto n = num_qubits();
+  auto dim = 1UL << n;
+  using Triplet = Eigen::Triplet<std::complex<double>>;
+  using SpMat = Eigen::SparseMatrix<std::complex<double>>;
+  std::vector<Triplet> xT{Triplet{0, 1, 1}, Triplet{1, 0, 1}},
+      iT{Triplet{0, 0, 1}, Triplet{1, 1, 1}},
+      yT{{0, 1, std::complex<double>{0, -1}},
+         {1, 0, std::complex<double>{0, 1}}},
+      zT{Triplet{0, 0, 1}, Triplet{1, 1, -1}};
+  SpMat x(2, 2), y(2, 2), z(2, 2), i(2, 2), mat(dim, dim);
+  x.setFromTriplets(xT.begin(), xT.end());
+  y.setFromTriplets(yT.begin(), yT.end());
+  z.setFromTriplets(zT.begin(), zT.end());
+  i.setFromTriplets(iT.begin(), iT.end());
+
+  auto kronProd = [](const std::vector<SpMat> &ops) -> SpMat {
+    SpMat ret = ops[0];
+    for (std::size_t k = 1; k < ops.size(); ++k)
+      ret = Eigen::kroneckerProduct(ret, ops[k]).eval();
+    return ret;
+  };
+
+  for_each_term([&](spin_op &term) {
+    auto termStr = term.to_string(false);
+    std::vector<SpMat> operations;
+    for (std::size_t k = 0; k < termStr.length(); ++k) {
+      auto pauli = termStr[k];
+      if (pauli == 'X')
+        operations.emplace_back(x);
+      else if (pauli == 'Y')
+        operations.emplace_back(y);
+      else if (pauli == 'Z')
+        operations.emplace_back(z);
+      else
+        operations.emplace_back(i);
+    }
+
+    mat += term.get_coefficient() * kronProd(operations);
+  });
+
+  std::vector<std::complex<double>> values;
+  std::vector<std::size_t> rows, cols;
+  for (int k = 0; k < mat.outerSize(); ++k)
+    for (SpMat::InnerIterator it(mat, k); it; ++it) {
+      values.emplace_back(it.value());
+      rows.emplace_back(it.row());
+      cols.emplace_back(it.col());
+    }
+
+  return std::make_tuple(values, rows, cols);
+}
+
 std::complex<double> spin_op::get_coefficient() const {
   if (terms.size() != 1)
     throw std::runtime_error(

runtime/cudaq/spin_op.h

@@ -349,6 +349,15 @@ class spin_op {
   /// @brief Return a dense matrix representation of this
   /// spin_op.
   complex_matrix to_matrix() const;
+
+  /// @brief Typedef for a vector of non-zero sparse matrix elements.
+  using csr_spmatrix =
+      std::tuple<std::vector<std::complex<double>>, std::vector<std::size_t>,
+                 std::vector<std::size_t>>;
+
+  /// @brief Return a sparse matrix representation of this `spin_op`. The
+  /// return type encodes all non-zero `(row, col, value)` elements.
+  csr_spmatrix to_sparse_matrix() const;
 };
 
 /// @brief Add a double and a spin_op

unittests/spin_op/SpinOpTester.cpp

@@ -149,6 +149,19 @@ TEST(SpinOpTester, canBuildDeuteron) {
   EXPECT_EQ(2, H.num_qubits());
 }
 
+TEST(SpinOpTester, checkGetSparseMatrix) {
+  auto H = 5.907 - 2.1433 * x(0) * x(1) - 2.1433 * y(0) * y(1) + .21829 * z(0) -
+           6.125 * z(1);
+  auto matrix = H.to_matrix();
+  matrix.dump();
+  auto [values, rows, cols] = H.to_sparse_matrix();
+  for (std::size_t i = 0; auto &el : values) {
+    std::cout << rows[i] << ", " << cols[i] << ", " << el << "\n";
+    EXPECT_NEAR(matrix(rows[i], cols[i]).real(), el.real(), 1e-3);
+    i++;
+  }
+}
+
 TEST(SpinOpTester, checkGetMatrix) {
   auto H = 5.907 - 2.1433 * x(0) * x(1) - 2.1433 * y(0) * y(1) + .21829 * z(0) -
            6.125 * z(1);