sparse solve returns wrong (and non-deterministic) results for symmetric indefinite KKT systems

[binxuwang2002-lang]

cuDSS issue draft — "sparse solve returns wrong (and non-deterministic) results for symmetric indefinite KKT systems"

Status: ready to file. Target repo: NVIDIA/cudss (issues). Filing checklist at the bottom.
All numbers below are reproducible with the attached MRE (code/cudss_mre.c) and the exported
matrices (code/export_kkt_mm.py + mre_out/iter_*.txt). Raw cross-check logs:
results/gpu01/xcheck.json (jax/spineax path, RTX 3080, cu12) and
results/gpu01/xcheck_a100_patched.log (jax/spineax path, A100, official cu13 wheel).

---

Title

cudssSolve returns incorrect and non-deterministic results for symmetric indefinite systems (0.8.0.10, both cu12 and cu13 builds)

Summary

A minimal C program that calls the plain cuDSS flow (cudssCreate → cudssExecute(ANALYSIS) →
cudssExecute(FACTORIZATION) → cudssExecute(SOLVE), CUDSS_MTYPE_SYMMETRIC /
CUDSS_MVIEW_UPPER, CSR, double) gets residuals of 0.5–7 in ||Ax-b||_inf on 8 of 15 real
symmetric indefinite KKT matrices of size n=1136 / nnz=5571. LAPACK (numpy.linalg.solve on the
same matrix in the same precision) returns 1e-11 on all 15. Repeated runs of the same binary on
the same input return different residuals every time (e.g. 3.2, 52.3, 43.3, 26.6, 98.1), so the
result is not just inaccurate, it is non-deterministic.

The matrices are the per-iteration KKT systems of an interior-point solver
(z = [x; s; y_c; y_d], symmetric indefinite, inertia roughly [+626, −510]). Some of them are
near-singular (one has a dense-solve residual of ~1e-4 already), but the majority are perfectly
solvable in double precision, and for those the dense reference is at 1e-11.

Environment (both stacks independently reproduce)

	machine A	machine B
GPU	A100 80GB PCIe (CC 8.0)	GeForce RTX 3080 20GB
driver	595.91.07	535.309.01
cuDSS	0.8.0.10 (nvidia-cudss-cu13 pip wheel, linked against nvidia/cu13/lib/libcudss.so.0)	0.8.0.10 (nvidia-cudss-cu12)
CUDA runtime	13.4.2 (toolkit 12.4 present, unused by the MRE)	12.6
host	x86_64 Ubuntu 22.04, gcc 11.4	x86_64 WSL2, gcc 11.4

The reproducer does not use any language binding — it links cudss.h + libcudss.so.0 directly
and manages all buffers manually, so nothing on our side (JAX, sympy codegen, our own CSR
assembly) is in the loop.

Repro

# build (libs from the nvidia pip wheels; any cuDSS 0.8.0.10 install works)
gcc -O2 -o cudss_mre cudss_mre.c \
    -I$VENV/lib/python3.12/site-packages/nvidia/cu13/include \
    -I/usr/local/cuda/include \
    -L$VENV/lib/python3.12/site-packages/nvidia/cu13/lib \
    -l:libcudss.so.0 -l:libcudart.so.13 -lm \
    -Wl,-rpath,$VENV/lib/python3.12/site-packages/nvidia/cu13/lib

# one failing system (n=1136, nnz=5571):
./cudss_mre iter_10.txt
# [mre] n=1136 nnz=5571 max_col=1135 b_norm=6.260191e+02
# RESULT n=1136 nnz=5571 ||Ax-b||_inf=1.807062e+01 ||b||_inf=6.145458e+02 rel=2.940484e-02 <-- FAILURE

# determinism: same input, five runs
for i in 1 2 3 4 5; do ./cudss_mre iter_10.txt | grep RESULT; done
# 3.2e+00 / 5.2e+01 / 4.3e+01 / 2.7e+01 / 9.8e+01   (all different)

# reference: dense LAPACK on the same matrix and rhs
python3 -c "import numpy as np; ... "   # ||Ax-b||_inf = 4.5e-12

Input format (written by export_kkt_mm.py from the solver's own CSR arrays): n nnz, then
nnz lines of row col val (0-based, upper triangle including diagonal, column indices are
not required to be sorted — see question 3 below), then n right-hand-side values.

Observed (15 exported systems, one solve each)

iter_1  2.268977e+00   FAIL      iter_10 (5 runs)  3.2 / 52.3 / 43.3 / 26.6 / 98.1   FAIL
iter_2  9.426486e-10   ok        iter_11          3.310241e-11                     ok
iter_3  8.441980e-11   ok        iter_12          6.407230e-11                     ok
iter_4  5.016689e-10   ok        iter_13          5.983779e-01                     FAIL
iter_5  1.006816e+00   FAIL      iter_14          7.709061e-01                     FAIL
iter_6  8.308368e-01   FAIL      iter_15          4.935396e-01                     FAIL
iter_7  1.006560e+00   FAIL
iter_8  7.450581e-09   ok
iter_9  7.250556e-11   ok

Pass rate 7/15. Failing residuals are 12 orders of magnitude above the dense reference.

Additional observation: a single manual iterative-refinement step (compute the residual on the
host, re-solve for the correction with the same factorization) reduces the residual to ~1e-4
for most failing cases (iter_1: 1.9 → 8.2e-5, iter_5: 0.77 → 1.0e-4, iter_13: 0.70 → 2.7e-4),
but makes one case worse (iter_10: 13.6 → 38.3). CUDSS_PHASE_SOLVE includes
CUDSS_PHASE_SOLVE_REFINEMENT in 0.8.0.10, so I would expect the library's own refinement to
cover this; it evidently does not.

Expected

||Ax-b||_inf on the order of eps * cond(A) * ||b||, i.e. what a dense LU gives on the same
inputs (1e-11 … 1e-13 here), and a deterministic answer for a deterministic input.

What would help

1. Confirm/deny on your side with the attached MRE + matrices (iter_*.txt, ~1.5 MB total).
   The matrix set covers both "ordinary" and near-singular members of the same family.
2. Is there a documented path by which cudssSolve can silently produce a wrong answer for
   CUDSS_MTYPE_SYMMETRIC + CUDSS_MVIEW_UPPER when (a) the matrix is indefinite,
   (b) column indices within a row are not sorted, and/or (c) the matrix is near-singular?
   Any of these three is present in our set; I have not isolated which one is the trigger.
3. The run-to-run variation on identical input suggests a race (library-internal threading?).
   Is there a supported way to force single-threaded/deterministic execution to check, other
   than rebuilding with a different threading layer?
4. Is iterative refinement supposed to run as part of CUDSS_PHASE_SOLVE, and under what
   conditions is it skipped?

Not ruled out on our side

• Our CSR arrays are produced by a code generator (jax2sympy) and are passed exactly as
  generated; rowPtr/colInd are validated (0 ≤ idx < n, rowPtr monotone, nnz consistent)
  inside the MRE, and the host residual uses a matrix rebuilt from the same arrays, so an
  input-construction error would show up as a large residual for the dense solve too — it
  does not.
• We have not tried other reorderings (CUDSS_CONFIG_REORDERING) or the hybrid memory mode.
• Two independent cuDSS builds (cu12 and cu13 wheels) and two different GPUs/drivers show the
  same behaviour, which is why I believe this is worth your time rather than a local install
  problem.

Attachments

• cudss_mre.c — self-contained reproducer (no dependencies beyond cuDSS + CUDA runtime).
• iter_1..15.txt — the 15 systems in the format described above.
• xcheck.json / xcheck_a100_patched.log — the same failure observed through a Python
  binding (for cross-reference only; the C path is the authoritative one).

Filing checklist

• attach cudss_mre.c and the 15 .txt matrices (or a tarball mre.tar.gz).
• state versions exactly as in the table above.
• include the determinism loop output verbatim.
• mention that CUDSS_PHASE_SOLVE includes the refinement phase and it did not help.
• do not claim "this is a bug in NVIDIA's library"; claim "the plain documented flow
  returns wrong/non-deterministic results on these inputs" — the numbers speak for themselves.

cudss_mre_pkg.tar.gz

[kvoronin]

Hi @binxuwang2002-lang!

Thanks for reporting the issue with all details, it made investigation very easy!

Our analysis:

• the issue is confirmed with the reproducer for cudss 0.8.0
• primary reason for poor accuracy is a large number of perturbed pivots due to a (typical for KKT) large zero block with the cudss settings used.
• for cudss 0.8.0, the way to restore the accuracy is enabling matching (also described in https://docs.nvidia.com/cuda/cudss/tips_and_tricks.html#accuracy-issues). However, it is sub-optimal because in cudss 0.8.0 matching is non-symmetric and thus makes cudss effectively use non-symmetric code path internally (~up to 1.5-2x more memory and less perf potentially).
• when we introduce the new features, the right settings would be to use symmetric matching and Bunch-Kaufman pivoting + potentially IR if there are still perturbed pivots reported.

Let us know if this solves your issue or you have questions remaining.

[binxuwang2002-lang]

Hi @kvoronin,

Thank you for the quick and thorough analysis — knowing the root cause (pivot perturbation
on the KKT zero block) rather than just "works in a newer version" is exactly what we needed,
and it made our own follow-up testing much easier to design.

We re-ran the reproducer on a third independent stack to see how far the 0.8.0.10 settings
get us (RTX 4070 Laptop, driver 566.26, nvidia-cudss-cu12 0.8.0.10, CUDA 12.9; same 15
systems, median of 3 runs, host-side residual code unchanged):

configuration                                  failing systems   notes
matching OFF (default)                              8/15          residual spread across repeats up to 80
+ CUDSS_MATCHING_ALG_MAX_DIAG_COUNT                  8/15          no change
+ CUDSS_MATCHING_ALG_MAX_MIN_DIAG                    6/15          fixes 5 of the 8, but 2 previously-passing systems now fail
+ IR_N_STEPS=2                                       4/15          median residual ~5e-12
+ MAX_MIN_DIAG + IR_N_STEPS=2                        3/8*          one stubborn case gets worse (2.1e+03)
* of the originally failing 8

The honest summary: we could not find a configuration in 0.8.0.10 that fixes all systems,
so we cannot recommend matching as a mitigation to our users in good conscience. That said,
your root-cause analysis is corroborated by exactly this table: the failing set moves with
the pivot/matching choices and is never empty — which is what one expects if perturbed pivots
on the zero block are the mechanism. (If any of this is a misconfiguration on our side, we
would genuinely appreciate being corrected — we are driving cuDSS through the spineax
binding, and we may well be missing something.)

Three small questions, only if convenient:

1. For 0.8.x users, is the correct guidance simply "no reliable workaround until the
   symmetric-matching release"? We want to document this accurately for downstream users.
2. Is there a public tracking item for symmetric matching + Bunch-Kaufman + IR, so we can
   watch for it and re-validate?
3. For the cases where iterative refinement makes things worse (iteration 6 above): is
   CUDSS_CONFIG_PIVOT_EPSILON the intended knob in 0.8.0.10, or is it disabled on the
   symmetric path?

The 15 matrices and the reproducer remain available, and we are happy to re-run them
against any new build, or generate additional systems (different sizes and inertias) if
that would be useful for your regression testing.

Thanks again for the investigation — reports with a confirmed root cause are rare and very
much appreciated.

Best regards,
binxuwang2002-lang

cudss_mre_configs_pkg.tar.gz

[kvoronin]

Hi @binxuwang2002-lang!

Could you check the most robust matching algorithm which is CUDSS_MATCHING_ALG_MAX_DIAG_PRODUCT? I think the matching algorithms which you have used are often worse (one of the good reasons for that is that the option above also applies scaling).

Answering your questions:

1. For 0.8.x users, is the correct guidance simply "no reliable workaround until the
   symmetric-matching release"? We want to document this accurately for downstream users.

The correct guidance is to enable the default matching algorithm which is CUDSS_MATCHING_ALG_AUTO (for 0.8.0 equivalent to CUDSS_MATCHING_ALG_MAX_DIAG_PRODUCT) and (if applicable) tolerance-based iterative refinement (which should step back if an iteration makes the accuracy worse, but costs some extra time for checking the residual). After the factorization, one should also check the number of perturbed pivots via cudssDataGet() with CUDSS_DATA_NPIVOTS and if the returned number is greater than zero, expect some iterative refinement iterations to be beneficial.

2. Is there a public tracking item for symmetric matching + Bunch-Kaufman + IR, so we can
   watch for it and re-validate?

As far as I know, no, there is no such public tracking place and I am limited in sharing exact information about next releases but we acknowledge it is a very important feature. When we publish our new versions, we update our release notes at (almost) the very same moment and release notes will reflect the changes made in a new release.

3. For the cases where iterative refinement makes things worse (iteration 6 above): is
   CUDSS_CONFIG_PIVOT_EPSILON the intended knob in 0.8.0.10, or is it disabled on the
   symmetric path?

The math for iterative refinement is relatively simple and unfortunately, impractical (= does not help in practice): without having spectral information (which is much harder to get than solving a system), it is very hard to say whether iterative refinement will make accuracy worse or not. The most robust way is to check the tolerance after solving the system (without IR). If it is not good enough, do iterative refinement and see if accuracy gets better. Once it stops getting better, refinement loop should be exited.

Changing CUDSS_CONFIG_PIVOT_EPSILON (which then requires re-doing factorization) can and usually affects the IR convergence but the presence of any perturbed pivots implies that making any guarantees is hard/impossible.

Also, potentially, if you know more about your systems, you might want to use a more complicated refinement process (than the so called "simple iteration") but as any iterative method convergence is tricky to achieve in a robust way.

Best regards,
Kirill