Sparse parity-check matrix support for decoders (#550)

Continuation of [#440](#440).
Closes [#379](#379).

Introduces `cudaq::qec::sparse_binary_matrix` (CSC/CSR, `uint32_t`
indices) as the decoder PCM type, so large PCMs no longer have to be
materialized as a dense `cudaqx::tensor`. The sparse path is wired
through PCM generation, round slicing, decoder construction, and Python
bindings. `sliding_window` canonicalizes once and reads `H.ptr() /
H.indices()` directly per window.

## Breaking changes

- **Decoder constructor** now takes `const sparse_binary_matrix &H`.
Dense `cudaqx::tensor<uint8_t>` callers still compile via implicit
conversion, but the extension-point template changed — **out-of-tree
plugins must be recompiled**.
- **`generate_random_pcm` throws above `k_max_dense_pcm_elements`
(400M)**. Use `generate_random_pcm_sparse` for larger matrices.
- **`single_error_lut_example`** now uses column `qErr` of H directly
instead of XOR-of-other-columns (the old form was only correct when
every row of H had even weight). Example plugin only; production
`single_error_lut` was already correct.

## Follow-up PRs

- Sparse-aware PCM utilities: migrate sort_pcm_columns, simplify_pcm,
reorder_pcm_columns, shuffle_pcm_columns, pcm_extend_to_n_rounds,
pcm_to_sparse_*, pcm_from_sparse_*, and
detector_error_model::canonicalize_for_rounds to sparse_binary_matrix so
the full #379 scale pipeline avoids any dense cudaqx::tensor<uint8_t>
allocation.
- Docs: Sphinx entries for the new type/functions and the 400M dense
cap.

Note to reviewers: Quick context on why it's larger than the ~1000-line
guideline: it's a foundational type-system change rather than an
additive feature. The new sparse_binary_matrix type, the decoder
constructor breaking change, all the plugin updates, the sparse
generator, and the Python bindings all have to land together for the
#379 fix to actually be usable end-to-end. Splitting it would create
intermediate states where the breaking change has landed but users still
can't construct a >50k-mechanism PCM which would mean shipping the API
churn without the benefit.

---------

Signed-off-by: Ben Howe <bhowe@nvidia.com>
Signed-off-by: vedika-saravanan <vsaravanan@nvidia.com>
Signed-off-by: Angela Burton <angelab@nvidia.com>
Co-authored-by: Ben Howe <bhowe@nvidia.com>
Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
Co-authored-by: GitHub Action <action@github.com>
Co-authored-by: Angela Burton <angelab@nvidia.com>

Author: 5 people

libs/qec/include/cudaq/qec/decoder.h

@@ -1,5 +1,5 @@
 /****************************************************************-*- C++ -*-****
- * Copyright (c) 2024 - 2025 NVIDIA Corporation & Affiliates.                  *
+ * Copyright (c) 2024 - 2026 NVIDIA Corporation & Affiliates.                  *
  * All rights reserved.                                                        *
  *                                                                             *
  * This source code and the accompanying materials are made available under    *
@@ -11,6 +11,7 @@
 #include "cuda-qx/core/extension_point.h"
 #include "cuda-qx/core/heterogeneous_map.h"
 #include "cuda-qx/core/tensor.h"
+#include "sparse_binary_matrix.h"
 #include <future>
 #include <optional>
 #include <vector>
@@ -121,7 +122,8 @@ class async_decoder_result {
 /// arbitrary constructor parameters that can be unique to each specific
 /// decoder.
 class decoder
-    : public cudaqx::extension_point<decoder, const cudaqx::tensor<uint8_t> &,
+    : public cudaqx::extension_point<decoder,
+                                     const cudaq::qec::sparse_binary_matrix &,
                                      const cudaqx::heterogeneous_map &> {
 private:
   struct rt_impl;
@@ -134,11 +136,9 @@ class decoder
   decoder() = delete;
 
   /// @brief Constructor
-  /// @param H Decoder's parity check matrix represented as a tensor. The tensor
-  /// is required be rank 2 and must be of dimensions \p syndrome_size x
-  /// \p block_size.
-  /// will use the same \p H.
-  decoder(const cudaqx::tensor<uint8_t> &H);
+  /// @param H Decoder's parity check matrix. Taken by value so rvalue
+  /// arguments are moved into the base member.
+  decoder(cudaq::qec::sparse_binary_matrix H);
 
   /// @brief Decode a single syndrome
   /// @param syndrome A vector of syndrome measurements where the floating point
@@ -174,7 +174,7 @@ class decoder
 
   /// @brief This `get` overload supports default values.
   static std::unique_ptr<decoder>
-  get(const std::string &name, const cudaqx::tensor<uint8_t> &H,
+  get(const std::string &name, const cudaq::qec::sparse_binary_matrix &H,
       const cudaqx::heterogeneous_map &param_map = cudaqx::heterogeneous_map());
 
   std::size_t get_block_size() { return block_size; }
@@ -273,7 +273,7 @@ class decoder
   std::size_t syndrome_size = 0;
 
   /// @brief The decoder's parity check matrix
-  cudaqx::tensor<uint8_t> H;
+  sparse_binary_matrix H;
 
   /// @brief The decoder's observable matrix in sparse format
   std::vector<std::vector<uint32_t>> O_sparse;
@@ -425,6 +425,6 @@ inline void convert_vec_hard_to_soft(const std::vector<std::vector<t_hard>> &in,
 }
 
 std::unique_ptr<decoder>
-get_decoder(const std::string &name, const cudaqx::tensor<uint8_t> &H,
+get_decoder(const std::string &name, const cudaq::qec::sparse_binary_matrix &H,
             const cudaqx::heterogeneous_map options = {});
 } // namespace cudaq::qec

libs/qec/include/cudaq/qec/pcm_utils.h

@@ -1,5 +1,5 @@
 /****************************************************************-*- C++ -*-****
- * Copyright (c) 2025 NVIDIA Corporation & Affiliates.                         *
+ * Copyright (c) 2025 - 2026 NVIDIA Corporation & Affiliates.                  *
  * All rights reserved.                                                        *
  *                                                                             *
  * This source code and the accompanying materials are made available under    *
@@ -8,6 +8,7 @@
 #pragma once
 
 #include "cuda-qx/core/tensor.h"
+#include "sparse_binary_matrix.h"
 #include <limits>
 #include <random>
 
@@ -111,6 +112,14 @@ get_sorted_pcm_column_indices(const cudaqx::tensor<uint8_t> &pcm,
 bool pcm_is_sorted(const cudaqx::tensor<uint8_t> &pcm,
                    std::uint32_t num_syndromes_per_round = 0);
 
+/// @brief Check if a PCM is sorted.
+/// @param sparse_pcm The sparse PCM to check (in the same format as
+/// dense_to_sparse())
+/// @param num_syndromes_per_round The number of syndromes per round.
+/// @return True if the PCM is sorted, false otherwise.
+bool pcm_is_sorted(const std::vector<std::vector<std::uint32_t>> &sparse_pcm,
+                   std::uint32_t num_syndromes_per_round = 0);
+
 /// @brief Reorder the columns of a PCM according to the given column order.
 /// Note: this may return a subset of the columns in the original PCM if the
 /// \p column_order does not contain all of the columns in the original PCM.
@@ -165,7 +174,53 @@ get_pcm_for_rounds(const cudaqx::tensor<uint8_t> &pcm,
                    bool straddle_start_round = false,
                    bool straddle_end_round = false);
 
+/// @brief Same semantics as the overload taking a dense tensor \p pcm, but
+/// reads from \p pcm as ``sparse_binary_matrix`` so the full dense PCM is not
+/// required (only the returned sub-matrix is dense). Parameter meanings match
+/// the dense overload.
+///
+/// @param pcm_is_canonical If true, the caller asserts \p pcm has
+/// sorted-unique per-group indices (i.e. is the output of `canonicalize_pcm`
+/// or was constructed canonically). In that case the per-call
+/// canonicalization step is skipped — useful for callers like
+/// `sliding_window` that canonicalize once at construction and then call
+/// this function in a per-window loop on the same matrix. Default `false`
+/// preserves the original "canonicalize on entry" behavior.
+///
+/// @warning When \p pcm_is_canonical is `true` the precondition is *not*
+/// checked. `select_pcm_columns_for_round_range` reads `.front()` / `.back()`
+/// of each CSC column to derive first/last round; those are only the true
+/// min/max row if the column list is sorted-unique. Passing a non-canonical
+/// PCM with this flag (e.g. a raw DEM decomposition, a hand-built sparse
+/// matrix with duplicate or unsorted indices, or a canonical CSR that the
+/// caller forgot to re-canonicalize after conversion) silently produces
+/// wrong round assignments. If unsure, leave the flag `false`.
+/// @return A tuple with the new PCM with the columns in the range [start_round,
+/// end_round], the first column included, and the last column included.
+std::tuple<cudaqx::tensor<uint8_t>, std::uint32_t, std::uint32_t>
+get_pcm_for_rounds(const sparse_binary_matrix &pcm,
+                   std::uint32_t num_syndromes_per_round,
+                   std::uint32_t start_round, std::uint32_t end_round,
+                   bool straddle_start_round = false,
+                   bool straddle_end_round = false,
+                   bool pcm_is_canonical = false);
+
+/// @brief Upper bound on \f$\texttt{rows} \times \texttt{cols}\f$ for the
+/// dense `generate_random_pcm` path. Above this, the call throws and the
+/// caller should switch to `generate_random_pcm_sparse`. Exposed as a named
+/// symbol so tests and users can reference it without grepping the source.
+inline constexpr std::size_t k_max_dense_pcm_elements =
+    static_cast<std::size_t>(400u) * 1024u * 1024u;
+
 /// @brief Generate a random PCM with the given parameters.
+///
+/// The PCM has shape `(n_rounds * n_syndromes_per_round) × (n_rounds *
+/// n_errs_per_round)`. If the product
+/// \f$n\_rounds^2 \times n\_syndromes\_per\_round \times n\_errs\_per\_round\f$
+/// exceeds `k_max_dense_pcm_elements`
+/// (\f$400 \times 1024 \times 1024 \approx 4.19 \times 10^8\f$ entries),
+/// throws `std::invalid_argument`; use `generate_random_pcm_sparse` for
+/// matrices that large without a dense allocation.
 /// @param n_rounds The number of rounds in the PCM.
 /// @param n_errs_per_round The number of errors per round in the PCM.
 /// @param n_syndromes_per_round The number of syndromes per round in the PCM.
@@ -178,6 +233,27 @@ cudaqx::tensor<uint8_t> generate_random_pcm(std::size_t n_rounds,
                                             std::size_t n_syndromes_per_round,
                                             int weight, std::mt19937_64 &&rng);
 
+/// @brief Same distribution as generate_random_pcm, but constructs a CSC
+/// sparse_binary_matrix directly without allocating a dense rank-2 tensor.
+/// Intended for large PCMs whose dense form would exceed
+/// k_max_dense_pcm_elements.
+sparse_binary_matrix
+generate_random_pcm_sparse(std::size_t n_rounds, std::size_t n_errs_per_round,
+                           std::size_t n_syndromes_per_round, int weight,
+                           std::mt19937_64 &&rng);
+
+/// @brief Return a GF(2)-canonical copy of \p pcm: each column (CSC) or row
+/// (CSR) has its indices sorted ascending, and duplicate indices are XOR-
+/// merged — an index that appears \p k times is kept iff \p k is odd. The
+/// output has the same layout as the input and is idempotent under further
+/// `canonicalize_pcm` calls.
+///
+/// Use this when a PCM source, such as a DEM decomposition, may legitimately
+/// emit duplicate indices within a column/row and the caller wants to apply
+/// GF(2) duplicate-collapse semantics before passing the matrix to consumers
+/// that require at-most-one entry per row per column.
+sparse_binary_matrix canonicalize_pcm(const sparse_binary_matrix &pcm);
+
 /// @brief Randomly permute the columns of a PCM.
 /// @param pcm The PCM to permute.
 /// @param rng The random number generator to use (e.g.

libs/qec/include/cudaq/qec/sparse_binary_matrix.h

@@ -0,0 +1,127 @@
+/****************************************************************-*- C++ -*-****
+ * Copyright (c) 2026 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.    *
+ ******************************************************************************/
+#pragma once
+
+#include "cuda-qx/core/tensor.h"
+#include <cstdint>
+#include <vector>
+
+namespace cudaq::qec {
+
+/// @brief Storage layout for the sparse PCM: Compressed Sparse Column (CSC)
+/// or Compressed Sparse Row (CSR). All non-zero entries are assumed to be 1;
+/// values are not stored.
+enum class sparse_binary_matrix_layout { csc, csr };
+
+/// @brief Sparse parity-check matrix in either CSC or CSR form.
+///
+/// Input index lists are stored as given: not required to be sorted or
+/// GF(2)-unique. Consumers that require cuSPARSE-style compressed groups can
+/// call `validate_sorted_unique_indices`; consumers that need GF(2)-collapsed
+/// per-group indices can call `cudaq::qec::canonicalize_pcm` on entry.
+///
+/// `index_type` is `uint32_t`, so each dimension and `nnz` must fit in
+/// `~4×10^9`.
+class sparse_binary_matrix {
+public:
+  using index_type = std::uint32_t;
+
+  /// @brief Construct a sparse PCM in CSC form.
+  /// @param num_rows Number of rows.
+  /// @param num_cols Number of columns.
+  /// @param col_ptrs Column pointer array (length num_cols + 1); column j has
+  /// indices in \p row_indices[col_ptrs[j] .. col_ptrs[j+1]-1].
+  /// @param row_indices Row indices of non-zeros (length nnz).
+  static sparse_binary_matrix from_csc(index_type num_rows, index_type num_cols,
+                                       std::vector<index_type> col_ptrs,
+                                       std::vector<index_type> row_indices);
+
+  /// @brief Construct a sparse PCM in CSR form.
+  /// @param num_rows Number of rows.
+  /// @param num_cols Number of columns.
+  /// @param row_ptrs Row pointer array (length num_rows + 1); row i has
+  /// indices in \p col_indices[row_ptrs[i] .. row_ptrs[i+1]-1].
+  /// @param col_indices Column indices of non-zeros (length nnz).
+  static sparse_binary_matrix from_csr(index_type num_rows, index_type num_cols,
+                                       std::vector<index_type> row_ptrs,
+                                       std::vector<index_type> col_indices);
+
+  /// @brief Construct from nested CSC: \p nested[j] is the list of row indices
+  /// for column j; \p nested.size() must equal \p num_cols.
+  static sparse_binary_matrix
+  from_nested_csc(index_type num_rows, index_type num_cols,
+                  const std::vector<std::vector<index_type>> &nested);
+
+  /// @brief Construct from nested CSR: \p nested[i] is the list of column
+  /// indices for row i; \p nested.size() must equal \p num_rows.
+  static sparse_binary_matrix
+  from_nested_csr(index_type num_rows, index_type num_cols,
+                  const std::vector<std::vector<index_type>> &nested);
+
+  /// @brief Construct from a rank-2 dense PCM (any non-zero treated as 1).
+  /// Intentionally not `explicit` so call sites that take
+  /// `sparse_binary_matrix` accept a dense `cudaqx::tensor` unchanged.
+  sparse_binary_matrix(
+      const cudaqx::tensor<std::uint8_t> &dense,
+      sparse_binary_matrix_layout layout = sparse_binary_matrix_layout::csc);
+
+  sparse_binary_matrix() = default;
+  sparse_binary_matrix(const sparse_binary_matrix &) = default;
+  sparse_binary_matrix(sparse_binary_matrix &&) noexcept = default;
+  sparse_binary_matrix &operator=(const sparse_binary_matrix &) = default;
+  sparse_binary_matrix &operator=(sparse_binary_matrix &&) noexcept = default;
+
+  sparse_binary_matrix_layout layout() const { return layout_; }
+  index_type num_rows() const { return num_rows_; }
+  index_type num_cols() const { return num_cols_; }
+  index_type num_nnz() const {
+    return static_cast<index_type>(indices_.size());
+  }
+
+  /// @brief For CSC: ptr has length num_cols+1; for CSR: ptr has length
+  /// num_rows+1.
+  const std::vector<index_type> &ptr() const { return ptr_; }
+  /// @brief For CSC: row indices; for CSR: column indices.
+  const std::vector<index_type> &indices() const { return indices_; }
+
+  /// @brief Throw if each compressed column/row does not have strictly
+  /// increasing indices. This rejects duplicate entries in the stored layout.
+  void validate_sorted_unique_indices(
+      const char *context = "sparse_binary_matrix") const;
+
+  /// @brief Return a copy of this matrix in CSC layout. No-op if already CSC.
+  sparse_binary_matrix to_csc() const;
+
+  /// @brief Return a copy of this matrix in CSR layout. No-op if already CSR.
+  sparse_binary_matrix to_csr() const;
+
+  /// @brief Convert to a dense PCM tensor (rows x columns). Non-zero entries
+  /// are set to 1.
+  cudaqx::tensor<std::uint8_t> to_dense() const;
+
+  /// @brief Nested CSC: outer vector has size num_cols; inner vector for
+  /// column j lists row indices of non-zeros in that column.
+  std::vector<std::vector<index_type>> to_nested_csc() const;
+
+  /// @brief Nested CSR: outer vector has size num_rows; inner vector for row i
+  /// lists column indices of non-zeros in that row.
+  std::vector<std::vector<index_type>> to_nested_csr() const;
+
+private:
+  sparse_binary_matrix(sparse_binary_matrix_layout layout, index_type num_rows,
+                       index_type num_cols, std::vector<index_type> ptr,
+                       std::vector<index_type> indices);
+
+  sparse_binary_matrix_layout layout_ = sparse_binary_matrix_layout::csc;
+  index_type num_rows_ = 0;
+  index_type num_cols_ = 0;
+  std::vector<index_type> ptr_;
+  std::vector<index_type> indices_;
+};
+
+} // namespace cudaq::qec

libs/qec/lib/CMakeLists.txt

@@ -1,5 +1,5 @@
 # ============================================================================ #
-# Copyright (c) 2024 - 2025 NVIDIA Corporation & Affiliates.                   #
+# Copyright (c) 2024 - 2026 NVIDIA Corporation & Affiliates.                   #
 # All rights reserved.                                                         #
 #                                                                              #
 # This source code and the accompanying materials are made available under     #
@@ -43,6 +43,7 @@ set(QEC_SOURCES
   experiments.cpp
   pcm_utils.cpp
   plugin_loader.cpp
+  sparse_binary_matrix.cpp
   stabilizer_utils.cpp 
   decoders/lut.cpp
   decoders/sliding_window.cpp

libs/qec/lib/decoder.cpp

@@ -1,5 +1,5 @@
 /*******************************************************************************
- * Copyright (c) 2022 - 2025 NVIDIA Corporation & Affiliates.                  *
+ * Copyright (c) 2022 - 2026 NVIDIA Corporation & Affiliates.                  *
  * All rights reserved.                                                        *
  *                                                                             *
  * This source code and the accompanying materials are made available under    *
@@ -17,8 +17,10 @@
 #include <filesystem>
 #include <vector>
 
-INSTANTIATE_REGISTRY(cudaq::qec::decoder, const cudaqx::tensor<uint8_t> &)
-INSTANTIATE_REGISTRY(cudaq::qec::decoder, const cudaqx::tensor<uint8_t> &,
+INSTANTIATE_REGISTRY(cudaq::qec::decoder,
+                     const cudaq::qec::sparse_binary_matrix &)
+INSTANTIATE_REGISTRY(cudaq::qec::decoder,
+                     const cudaq::qec::sparse_binary_matrix &,
                      const cudaqx::heterogeneous_map &)
 
 // Include decoder implementations AFTER registry instantiation
@@ -70,12 +72,11 @@ struct decoder::rt_impl {
 
 void decoder::rt_impl_deleter::operator()(rt_impl *p) const { delete p; }
 
-decoder::decoder(const cudaqx::tensor<uint8_t> &H)
-    : H(H), pimpl(std::unique_ptr<rt_impl, rt_impl_deleter>(new rt_impl())) {
-  const auto H_shape = H.shape();
-  assert(H_shape.size() == 2 && "H tensor must be of rank 2");
-  syndrome_size = H_shape[0];
-  block_size = H_shape[1];
+decoder::decoder(cudaq::qec::sparse_binary_matrix H)
+    : H(std::move(H)),
+      pimpl(std::unique_ptr<rt_impl, rt_impl_deleter>(new rt_impl())) {
+  syndrome_size = this->H.num_rows();
+  block_size = this->H.num_cols();
   reset_decoder();
   pimpl->persistent_detector_buffer.resize(this->syndrome_size);
   pimpl->persistent_soft_detector_buffer.resize(this->syndrome_size);
@@ -130,7 +131,7 @@ decoder::decode_async(const std::vector<float_t> &syndrome) {
 }
 
 std::unique_ptr<decoder>
-decoder::get(const std::string &name, const cudaqx::tensor<uint8_t> &H,
+decoder::get(const std::string &name, const cudaq::qec::sparse_binary_matrix &H,
              const cudaqx::heterogeneous_map &param_map) {
   auto [mutex, registry] = get_registry();
   std::lock_guard<std::recursive_mutex> lock(mutex);
@@ -479,7 +480,7 @@ void decoder::reset_decoder() {
 }
 
 std::unique_ptr<decoder> get_decoder(const std::string &name,
-                                     const cudaqx::tensor<uint8_t> &H,
+                                     const cudaq::qec::sparse_binary_matrix &H,
                                      const cudaqx::heterogeneous_map options) {
   return decoder::get(name, H, options);
 }