Sparse interval-MDP abstraction (memory priority, ISSUE-0006)

Build the sparse solve::IMDPModel directly (only nonzero kernel-window successors),
the same sparse representation as IntervalMDP.jl => synthesis + storage are O(nnz),
not O(states^2).

src/abstraction.{h,cpp}: interval transition-bound kernel transitionInterval1D
(closed form: mass of N(mu,sigma^2) in [a,b] is unimodal in mu) + transitionIntervalBox
(axis-decoupled product) + a 1-D sparse reach abstraction for affine diagonal-Gaussian
systems (absorbing target/sink, feasible bounds).

Verified (tests/unit/test_abstraction.cpp): kernel vs brute-force mean sampling; box ==
product; lossless pruning (sparse synthesis == dense-window synthesis); target cell = 1;
O(N) sparsity (grow domain at fixed step/noise -> bounded nnz/cell).

Scaling (benchmarks/sparse_scaling.cpp): 125k states ~255 MB + 0.42s vs dense ~250 GB.
Honest caveat: O(N) holds with kernel-in-cells fixed (grow domain); refining eta at
fixed sigma is a constant-factor saving (window ~ 6 sigma / eta). v1 dense path unchanged.

Suite: 38/38 cases, 54375 assertions, 0 failures.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Author: Kiguli

benchmarks/RESULTS.md

@@ -22,7 +22,28 @@ report before publication):
 > co-safe spec collapsed into separate reach-avoid targets. Phase 2 replaces this
 > with a real DFA product (no reduction).
 
-## Baseline (v1, branch main) — TODO: fill via Docker/CI
+## v1 baseline note (ISSUE-0006)
+The v1 dense path OOMs even on the smallest shipped example: 2D-robot-RA builds a
+dense (state×input)×state matrix of 698103×1583 ≈ **8.84 GB** → killed under a 4 GB
+cap. Confirms the need for sparse storage before larger benchmarks.
+
+## Sparse abstraction scaling (v2.0, `benchmarks/sparse_scaling.cpp`)
+1-D reach IMDP, domain grown with grid step (eta=0.1) and noise (sigma=0.3) FIXED, so
+the kernel spans a constant #cells ⇒ nnz is O(N) while a dense (min+max) matrix is
+O(N²). Synthesis = optimistic VI (eps=1e-6).
+
+| domain | cells | nnz | nnz/cell | sparse (MB) | dense (MB) | synth (s) |
+|---|---|---|---|---|---|---|
+| 20 | 200 | 18,005 | 89 | 0.36 | 0.6 | 0.01 |
+| 100 | 1,000 | 99,980 | 100 | 2.0 | 16 | 0.00 |
+| 500 | 5,000 | 507,980 | 102 | 10.2 | 400 | 0.01 |
+| 2,500 | 25,000 | 2,547,980 | 102 | 51 | 10,000 | 0.08 |
+| 12,500 | 125,000 | 12,747,980 | 102 | **255** | **250,000** | 0.42 |
+
+At 125k states the sparse model is ~255 MB and synthesizes in 0.42 s, where the dense
+matrix would be ~250 GB (infeasible). This is the IntervalMDP.jl-style memory win.
+
+## Baseline (v1, branch main) — TODO: fill via Docker/CI on a tiny-enough model
 ```
 # paste run_archcomp.py output here
 ```

benchmarks/sparse_scaling.cpp

@@ -0,0 +1,41 @@
+// Sparse-abstraction scaling demo (ISSUE-0006): grow the DOMAIN with the grid step
+// (eta) and noise (sigma) FIXED, so the Gaussian kernel spans a constant number of
+// cells. Then stored nonzeros are O(N) while a dense (min+max) matrix is O(N^2):
+// large state counts that would OOM densely run in megabytes.
+// (Refining eta with fixed sigma instead would widen the kernel-in-cells and is NOT
+// asymptotically sparse — the win is fixed kernel width, growing N.)
+//
+// Build:  c++ -std=c++17 -O2 benchmarks/sparse_scaling.cpp \
+//             src/abstraction.cpp src/solve.cpp src/omaximization.cpp src/graph_utils.cpp -o /tmp/scaling
+// Run:    /tmp/scaling
+#include "../src/abstraction.h"
+#include "../src/solve.h"
+#include <cstdio>
+#include <chrono>
+
+using namespace impact;
+
+int main() {
+    printf("%10s %10s %12s %10s %12s %12s %9s\n",
+           "domain", "cells", "nnz", "nnz/cell", "sparse(MB)", "dense(MB)", "synth(s)");
+    for (double half : {10.0, 50.0, 250.0, 1250.0, 6250.0}) {   // grow domain, fixed eta+sigma
+        abstraction::System1D s;
+        s.a = 0.9; s.b = 1.0; s.sigma = 0.3;
+        s.xlb = -half; s.xub = half; s.eta = 0.1;   // fixed step => kernel ~ constant #cells
+        s.ulb = -1; s.uub = 1; s.ueta = 1.0;        // 3 inputs
+        s.tlo = half - 2; s.thi = half;
+
+        auto ab = abstraction::buildSparseReach1D(s, 1e-7);
+        auto t0 = std::chrono::steady_clock::now();
+        auto r = solve::maxReachOptimistic(ab.model, ab.targets, 1e-6);
+        double secs = std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count();
+
+        double sparseMB = ab.nnz * (double)(sizeof(int) + 2 * sizeof(double)) / 1e6;
+        double denseMB  = (double)ab.nCells * ab.nCells * 2 * sizeof(double) / 1e6; // min+max N x N
+        printf("%10.0f %10d %12lld %10.1f %12.2f %12.1f %9.2f\n",
+               2 * half, ab.nCells, ab.nnz, (double)ab.nnz / (ab.nCells + 2),
+               sparseMB, denseMB, secs);
+        (void)r;
+    }
+    return 0;
+}

issues/0006-dense-transition-matrix-oom.md

@@ -1,7 +1,7 @@
 ---
 id: ISSUE-0006
 title: Dense transition matrix causes memory overflow even on "small" cases
-status: open
+status: in-progress
 severity: high
 labels: scalability, performance, tool-v1, memory
 created: 2026-06-25
@@ -30,10 +30,19 @@ abstraction vectors computed (is/ss/ts/min*/max* h5 produced), then:
 - Confirms the plan's dense→sparse concern; products (IMDP×automaton, Phase 2/3)
   multiply state count and will hit this much sooner.
 
-## Resolution / plan
-- Short term: generate goldens from a *genuinely tiny* model (few states/inputs)
-  or raise the memory cap only where the host allows.
-- Core: **sparse transition representation** (CSR-style) for transition matrices
-  and products — scheduled as the Phase 2/3 cross-cutting "sparse" task in the plan.
-  This is the right place to also reduce the (state×input)×state blow-up.
-- Keep "very small case studies only" for local/CI runs until sparse lands.
+## Resolution / progress
+- **DONE — sparse synthesis path.** The v2.0 solver works on a sparse model
+  (`solve::IMDPModel` = per-state lists of only nonzero successor intervals), so
+  synthesis memory is O(nnz), not O(states²) — same model as IntervalMDP.jl.
+- **DONE — sparse abstraction (`src/abstraction.{h,cpp}`).** Builds the sparse IMDP
+  directly for affine, axis-decoupled, diagonal-Gaussian systems (1-D end-to-end),
+  storing only kernel-window successors. Verified: kernel vs brute-force; lossless
+  pruning (sparse synthesis == dense-window synthesis); O(N) scaling
+  (`benchmarks/sparse_scaling.cpp`): 125k states ≈ 255 MB + 0.42 s vs dense ≈ 250 GB.
+  **Caveat:** the O(N) win holds with grid step + noise FIXED and the domain growing
+  (constant kernel-in-cells). Refining eta at fixed sigma widens the kernel-in-cells
+  (window ≈ 6σ/eta), giving a constant-factor (not asymptotic) saving — documented.
+- **TODO:** n-D (diagonal) and coupled dynamics (mixed-monotone / NLopt bounds, as v1);
+  sparse for IMDP×automaton products; and either refactor v1 `IMDP.cpp` to sparse or
+  drive ARCH benchmarks through the new sparse pipeline. v1's dense path is unchanged
+  for now (still OOMs); use the sparse pipeline for large runs.

paper/IMPaCT-v2.md

@@ -160,6 +160,20 @@ Format per entry: `date — feature | decision + rationale | algorithm | refs |
   addition to LTL. Phase 3 pipelines chosen from verified literature (Dutreix-Coogan
   Rabin+permanent-WC primary; Weininger game-reduction; Asadi qualitative).
 
+- **Sparse abstraction (priority: memory / larger benchmarks, ISSUE-0006).**
+  `src/abstraction.{h,cpp}`: builds the sparse `solve::IMDPModel` directly for affine,
+  axis-decoupled, diagonal-Gaussian systems (storing only kernel-window successors) —
+  same sparse model as IntervalMDP.jl, so memory is O(nnz). Verified: the interval
+  kernel `transitionInterval1D` vs brute-force mean sampling; box == product of 1-D;
+  lossless pruning (sparse synthesis == dense-window synthesis); target-cell = 1.
+  Scaling (`benchmarks/sparse_scaling.cpp`, domain grown at fixed step/noise):
+  **125,000 states ≈ 255 MB + 0.42 s** vs a dense (min+max) matrix ≈ **250 GB** — the
+  large runs that v1's dense path OOMs on (8.84 GB on the *smallest* 2D-robot case).
+  Honest caveat recorded: O(N) holds for fixed kernel-in-cells (grow domain); refining
+  eta at fixed sigma is only a constant-factor saving. Suite: 38/38, 54375 assertions.
+  Next: n-D/coupled abstraction (mixed-monotone/NLopt bounds), sparse products, and
+  driving ARCH benchmarks through the sparse pipeline. Branch pushed to origin.
+
 <!-- add new dated entries above this line -->
 
 ---

src/abstraction.cpp

@@ -0,0 +1,108 @@
+#include "abstraction.h"
+
+#include <cmath>
+#include <algorithm>
+
+namespace impact {
+namespace abstraction {
+
+namespace { constexpr double INV_SQRT2 = 0.70710678118654752440; }
+
+double normalCdf(double z) { return 0.5 * std::erfc(-z * INV_SQRT2); }
+
+double massInInterval(double mu, double sigma, double a, double b) {
+    return normalCdf((b - mu) / sigma) - normalCdf((a - mu) / sigma);
+}
+
+Bound transitionInterval1D(double muLo, double muHi, double sigma, double a, double b) {
+    if (muHi < muLo) std::swap(muLo, muHi);
+    const double center = 0.5 * (a + b);
+    const double muMax = std::min(std::max(center, muLo), muHi);   // closest to centre
+    double hi = massInInterval(muMax, sigma, a, b);
+    double lo = std::min(massInInterval(muLo, sigma, a, b),
+                         massInInterval(muHi, sigma, a, b));         // farther endpoint
+    lo = std::max(0.0, std::min(1.0, lo));
+    hi = std::max(0.0, std::min(1.0, hi));
+    if (lo > hi) lo = hi;
+    return {lo, hi};
+}
+
+Bound transitionIntervalBox(const std::vector<double>& muLo,
+                            const std::vector<double>& muHi,
+                            const std::vector<double>& sigma,
+                            const std::vector<double>& aLo,
+                            const std::vector<double>& aHi) {
+    double lo = 1.0, hi = 1.0;
+    for (size_t d = 0; d < sigma.size(); ++d) {
+        Bound b = transitionInterval1D(muLo[d], muHi[d], sigma[d], aLo[d], aHi[d]);
+        lo *= b.lo; hi *= b.hi;
+    }
+    return {lo, hi};
+}
+
+SparseReach buildSparseReach1D(const System1D& sys, double prune) {
+    const int N = std::max(1, (int)std::llround((sys.xub - sys.xlb) / sys.eta));
+    const int TARGET = N, SINK = N + 1;
+    const int M = std::max(0, (int)std::llround((sys.uub - sys.ulb) / sys.ueta));  // input grid points = M+1
+
+    SparseReach out;
+    out.nCells = N;
+    out.nnz = 0;
+    out.model.assign(N + 2, {});
+    out.targets.insert(TARGET);
+
+    auto cellLo = [&](int j) { return sys.xlb + j * sys.eta; };
+    auto isTargetCell = [&](int j) {
+        double lo = cellLo(j), hi = lo + sys.eta;
+        return lo >= sys.tlo - 1e-12 && hi <= sys.thi + 1e-12;
+    };
+
+    const double W = 6.0 * sys.sigma;   // truncation window (~all mass within 6 sigma)
+
+    for (int i = 0; i < N; ++i) {
+        if (isTargetCell(i)) {                              // absorbing target cell
+            out.model[i].push_back({ {TARGET, 1.0, 1.0} });
+            out.nnz += 1;
+            continue;
+        }
+        const double xl = cellLo(i), xr = xl + sys.eta;
+        for (int k = 0; k <= M; ++k) {
+            const double u = sys.ulb + k * sys.ueta;
+            // affine mean range over the source cell (handle sign of a)
+            double m1 = sys.a * xl + sys.b * u, m2 = sys.a * xr + sys.b * u;
+            double muLo = std::min(m1, m2), muHi = std::max(m1, m2);
+
+            solve::ActionDist row;
+            // mass into the target region [tlo,thi]
+            Bound bT = transitionInterval1D(muLo, muHi, sys.sigma, sys.tlo, sys.thi);
+            if (bT.hi > prune) row.push_back({TARGET, bT.lo, bT.hi});
+
+            // mass into non-target grid cells within the support window
+            int jmin = (int)std::floor((muLo - W - sys.xlb) / sys.eta);
+            int jmax = (int)std::floor((muHi + W - sys.xlb) / sys.eta);
+            jmin = std::max(0, jmin); jmax = std::min(N - 1, jmax);
+            for (int j = jmin; j <= jmax; ++j) {
+                if (isTargetCell(j)) continue;              // counted in the TARGET aggregate
+                double cl = cellLo(j), cr = cl + sys.eta;
+                Bound b = transitionInterval1D(muLo, muHi, sys.sigma, cl, cr);
+                if (b.hi > prune) row.push_back({j, b.lo, b.hi});
+            }
+
+            // remaining mass (outside grid / pruned) -> SINK, with feasible bounds
+            double sumLo = 0, sumHi = 0;
+            for (const auto& iv : row) { sumLo += iv.lo; sumHi += iv.hi; }
+            double sinkLo = std::max(0.0, 1.0 - sumHi);
+            double sinkHi = std::min(1.0, 1.0 - sumLo);
+            row.push_back({SINK, sinkLo, sinkHi});
+
+            out.nnz += (long long)row.size();
+            out.model[i].push_back(std::move(row));
+        }
+    }
+    out.model[TARGET].push_back({ {TARGET, 1.0, 1.0} });    // absorbing
+    out.model[SINK].push_back({ {SINK, 1.0, 1.0} });
+    return out;
+}
+
+} // namespace abstraction
+} // namespace impact

src/abstraction.h

@@ -0,0 +1,75 @@
+#ifndef IMPACT_ABSTRACTION_H
+#define IMPACT_ABSTRACTION_H
+
+// ============================================================================
+// Sparse interval-MDP abstraction of continuous-state stochastic systems
+// (Phase: sparse / scalability — ISSUE-0006). Produces a SPARSE solve::IMDPModel
+// (per-state lists of only the nonzero successor probability intervals), the same
+// memory model as IntervalMDP.jl, so memory is O(nnz) not O(states^2).
+//
+// This header provides the correctness-critical kernel: the transition-probability
+// INTERVAL [lo,hi] from a SOURCE CELL (a box of states) to a TARGET box, under a
+// Gaussian noise model. For affine, axis-decoupled dynamics the per-dimension mean
+// ranges are independent, so the box probability is the product of 1-D factors and
+// its bounds are the product of the per-dimension bounds.
+//
+// Refs: interval-MDP abstraction of stochastic systems — Lahijanian-Andersson-Belta
+// (IEEE TAC 2015); FAUST^2 (TACAS 2015); AMYTISS (CAV 2020); IntervalMDP.jl (arXiv
+// 2401.04068) for the sparse representation. Verified against brute-force numerics
+// in tests/unit/test_abstraction.cpp.
+// ============================================================================
+
+#include <vector>
+#include "solve.h"
+
+namespace impact {
+namespace abstraction {
+
+    // Standard normal CDF.
+    double normalCdf(double z);
+
+    // Probability mass of N(mu, sigma^2) in the interval [a, b].
+    double massInInterval(double mu, double sigma, double a, double b);
+
+    struct Bound { double lo; double hi; };
+
+    // The 1-D transition probability mass of N(mu, sigma^2) in [a,b], minimized and
+    // maximized over the mean mu ranging in [muLo, muHi]. Closed form: the mass is
+    // unimodal in mu with peak at the box centre (a+b)/2 — so the max is at mu
+    // clamped to the centre and the min at the farther endpoint. Requires muLo<=muHi,
+    // a<=b, sigma>0.
+    Bound transitionInterval1D(double muLo, double muHi, double sigma, double a, double b);
+
+    // Axis-decoupled n-D box bound = product of per-dimension 1-D bounds.
+    Bound transitionIntervalBox(const std::vector<double>& muLo,
+                                const std::vector<double>& muHi,
+                                const std::vector<double>& sigma,
+                                const std::vector<double>& aLo,
+                                const std::vector<double>& aHi);
+
+    // --- 1-D sparse reach abstraction (the first end-to-end, verifiable case) -----
+    // System: x' = a*x + b*u + N(0, sigma^2), x on grid [xlb,xub] step eta, inputs
+    // on grid [ulb,uub] step ueta. Target region [tlo,thi] (absorbing). Mass leaving
+    // the grid (or into the target) is handled by absorbing TARGET / SINK states.
+    struct System1D {
+        double a, b, sigma;
+        double xlb, xub, eta;
+        double ulb, uub, ueta;
+        double tlo, thi;
+    };
+
+    struct SparseReach {
+        solve::IMDPModel model;     // states: 0..N-1 cells, N = TARGET, N+1 = SINK
+        std::set<int> targets;      // { N }
+        int nCells;
+        long long nnz;              // total stored successor intervals (sparsity metric)
+    };
+
+    // Build the sparse IMDP. `prune` drops successor cells whose upper bound <= prune
+    // (mass folded into SINK); prune=0 keeps every cell (a dense-equivalent build).
+    SparseReach buildSparseReach1D(const System1D& sys, double prune);
+
+} // namespace abstraction
+} // namespace impact
+
+#endif // IMPACT_ABSTRACTION_H

tests/Makefile

@@ -22,14 +22,16 @@ IMPL_SRCS := ../src/omaximization.cpp \
              ../src/graph_utils.cpp \
              ../src/solve.cpp \
              ../src/ltl.cpp \
-             ../src/product.cpp
+             ../src/product.cpp \
+             ../src/abstraction.cpp
 
 TEST_SRCS := unit/test_main.cpp \
              unit/test_omaximization.cpp \
              unit/test_graph.cpp \
              unit/test_interval_iteration.cpp \
              unit/test_ltl_dfa.cpp \
-             unit/test_product.cpp
+             unit/test_product.cpp \
+             unit/test_abstraction.cpp
 
 BIN := run_tests
 

tests/contracts/contracts.h

@@ -39,4 +39,7 @@
 // Phase 2 part 3 — IMDP x DFA product + co-safe reachability.
 #include "../../src/product.h"
 
+// Sparse interval-MDP abstraction of continuous systems (scalability / ISSUE-0006).
+#include "../../src/abstraction.h"
+
 #endif // IMPACT_TEST_CONTRACTS_H