PaperlayerView exact revision on arXiv

Source-linked AI summary

DBDSVDMR3: An MRRR-Based Bidiagonal SVD Solver Built from DSTEGR

James Demmel, Inderjit Dhillon, Sai Surya Duvvuri, Julien Langou, Osni Marques, Ryan Schneider

arXiv:2610.30227v1math.NAcs.NA

TL;DR

MRRR is optimal for symmetric tridiagonal eigenproblems, but reliable adaptation to bidiagonal SVD remained unresolved after earlier implementations were abandoned. This paper builds DBDSVDMR3 from DSTEGR to implement the MR3-GK approach. It passes LAPACK tests and performs well on most test matrices, while retaining accuracy and complexity limitations.

  • Problem

    Earlier MRRR-based bidiagonal-SVD implementations were abandoned, leaving open whether the MR3-GK approach had a theoretical flaw or the implementations had bugs.

    Source: §1 · p. 3

  • Method

    DBDSVDMR3 implements Willems and Lang’s MR3-GK bidiagonal-SVD algorithm using code from LAPACK 3.0 DSTEGR.

    Source: §1 · p. 3

  • Results

    DBDSVDMR3 passes LAPACK tests and keeps error ratios below 100 on 88% of STCollection matrices, while outperforming DBDSVDX and DBDSQR by one to two orders of magnitude.

    Source: §1 · p. 3, §4.4.1 · p. 18

  • Takeaways & Limitations

    The results show that MRRR can compute the bidiagonal SVD via the Golub–Kahan matrix with minimal modification.

    Source: §4.4.1 · p. 18

    On tightly clustered eigenvalue cases, unguarded shifts can cause accuracy loss proportional to the depth of the representation tree.

    Source: §4.3 · p. 15

Abstract

from arXiv · show

The optimal algorithm for the symmetric tridiagonal eigenvalue problem is Multiple Relatively Robust Representations (MRRR or MR$^3$), which is capable of computing $k$ eigenvalue/eigenvector pairs in only $O(kn)$ floating-point operations. While MRRR is nearly 30 years old, it has never successfully been adapted to the bidiagonal SVD. This is not for lack of trying; two previous attempts produced implementations that were considered for LAPACK but eventually abandoned, a result of the fact that neither was reliable enough when run on test problems. The question therefore remained: is there a problem with the theory of MRRR for the bidiagonal SVD or were there simply bugs in the previous implementations that went undetected? This paper resolves this question by presenting DBDSVDMR3 -- a new implementation of the MRRR-based algorithm for the bidiagonal SVD originally proposed by Willems and Lang [ETNA 2012], which, for the first time, performs well on a variety of test problems while achieving significant speedups over both LAPACK's bisection-based solver and, in certain cases, divide-and-conquer for the bidiagonal SVD (thereby indicating that the theory is sound). DBDSVDMR3 will be available in the next release of LAPACK. In this paper, we discuss the design of our implementation, which was obtained with the help of Claude Code, and exhibit its performance on a diverse set of test problems.

1. Introduction.

The paper asks whether MRRR can reliably compute the bidiagonal SVD through the Golub–Kahan matrix, or whether earlier implementations failed because of bugs. It presents a new implementation and reports that it passes LAPACK tests and performs well on many test matrices, while identifying accuracy and complexity boundaries.

1. Introduction.

  • MRRR can compute k eigenpairs of a symmetric tridiagonal matrix in O(kn) flops, using accurate eigenvalues and relative gaps to obtain mutually orthogonal inverse-iteration eigenvectors by default.A bidiagonal SVD can be obtained by diagonalizing B^T B, BB^T, or the Golub–Kahan tridiagonal matrix, motivating an MRRR-based approach.Source: §1 · pp. 1–2
  • The open question is whether MRRR has a fundamental flaw for the bidiagonal SVD via MR3-GK, or whether earlier implementations failed because of code bugs.Two proposed implementations were abandoned, and the later MR3-GK implementation could fail even on simple examples.Source: §1 · p. 3
  • Error ratios stay below 100 on 110 of 125 STCollection matrices (88%), although the remaining clustered cases lose accuracy.The implementation passes LAPACK tests; its median ratios are about ten times those of existing LAPACK bidiagonal-SVD routines, while DBDSDC and DBDSQR are accurate throughout these tests and DBDSVDX fails twice.Source: §1 · p. 3
  • The new DBDSVDMR3 implementation builds Willems and Lang’s MR3-GK algorithm from LAPACK 3.0 DSTEGR code with help from Claude Code.It was merged into the reference LAPACK repository and is slated for the next release.Source: §1 · p. 3
  • DBDSVDMR3 does not yet provide high relative accuracy for singular values or O(kn) work for requested subsets, since it computes all n singular vectors.The implementation applies an absolute perturbation when preprocessing entries, and its subset computation therefore remains O(n^2).Source: §1 · p. 4

2. Background.

The section reviews MRRR for symmetric tridiagonals and its extension to bidiagonal SVD through the Golub–Kahan matrix. Under stated assumptions, the extension yields numerically orthogonal singular vectors and small coupled residuals.

2.1. MRRR for symmetric tridiagonals.

  • MRRR traverses a representation tree, refining eigenvalues and computing vectors from relatively robust representations; clustered eigenvalues trigger a shift and a new representation.For well-separated eigenvalues, inverse iteration on the representation computes vectors; differential qd routines support eigenvalue computation and shifting.Source: §2.1 · pp. 4, 6
  • The eigenvector error bound depends on the residual divided by the distance to other eigenvalues, explaining why accurate approximations and relative separation yield near-aligned vectors.Proposition 2.1 bounds sin(θ_i) by the residual norm divided by the minimum eigenvalue separation.Source: §2.1 · p. 5

2.2. MR3-GK for bidiagonals.

  • MR3-GK runs MRRR on the Golub–Kahan matrix using it as the root representation, then de-interleaves its positive-eigenvalue vectors into singular-vector pairs.The required nearly constant diagonal condition keeps representations close to ones associated with constant-diagonal tridiagonals, whose invariant subspaces have GK structure.Source: §2.1 · p. 6, §2.2 · pp. 6–7

2.3. Unified theory of correctness.

  • Under Willems and Lang’s assumptions plus the nearly constant diagonal condition, MR3-GK yields numerically orthogonal singular vectors and coupled residuals of order O(εn∥B∥2/gaptol).The requirements are listed in Table 2.1; the theory extends the assumptions for tridiagonal MRRR to the bidiagonal algorithm.Source: §2.3 · p. 7, Table 2.1 · p. 7

3. Adapting DSTEGR to the Bidiagonal SVD.

DBDSVDMR3 adapts DSTEGR to the bidiagonal SVD by preprocessing the Golub–Kahan matrix and preserving its structure through MRRR. The design avoids a shift that would merge positive and negative eigenvalue copies, and disables aggressive eigenvector truncation.

3. Adapting DSTEGR to the Bidiagonal SVD.

  • DBDSVDMR3 implements MR3-GK on DSTEGR, with changes to preprocessing, root selection, and eigenvector handling.The implementation was developed with an AI coding assistant under the authors’ direction.Source: §3 · pp. 7–8

3.1. Preprocessing.

  • Scaling protects dqds and twisted factorizations from overflow or underflow; zeroing entries at most ε∥B∥ enables splitting TGK into independent subproblems.This perturbation can forfeit high relative accuracy in singular values; when no entry falls below ε∥B∥, the split is skipped and singular values retain full relative accuracy.Source: §3.1 · pp. 8, 10
  • Zeroing a diagonal entry can produce odd-dimensional Golub–Kahan blocks, which remain valid root representations but require combining zero eigenvectors from separate blocks to recover both singular vectors.The zeroing step may sacrifice relative accuracy; the authors compare DBDSVDMR3 and DBDSQR’s accurate tiny singular values with DBDSDC’s absolute-level result for the stated examples.Source: §3.1 · pp. 9–10

3.2. TGK as root representation.

  • The usual positive-definite shift would merge positive and negative copies of each singular value into a cluster, so DBDSVDMR3 skips it and uses TGK directly as the root representation.The authors expect omitted element-growth monitoring to be a remedy for accuracy losses on heavily clustered spectra; it remains future work.Source: §3.2 · p. 10

3.3. Full-length eigenvectors.

  • DSTEGR’s early truncation can stop too soon because TGK eigenvectors interleave left and right singular-vector entries, so DBDSVDMR3 computes full-length eigenvectors.The paper notes a less aggressive truncation strategy might save work, but disables the optimization to avoid potential failures.Source: §3.3 · pp. 11–12

4. Numerical Results.

Tests compare DBDSVDMR3 with three LAPACK solvers and unmodified MR3GK on LAPACK-style and practical/stress-test matrices. DBDSVDMR3 is faster across nontrivial spectra but is less accurate on practical matrices and can fail on tightly clustered spectra.

4. Numerical Results.

  • The evaluation uses 16 classes of LAPACK-style tests and STCollection practical matrices and stress tests, comparing DBDSVDMR3 with three LAPACK solvers and MR3GK.MR3GK fails even on simple cases, while DBDSVDMR3 passes LAPACK’s test code.Source: §4 · p. 12
  • At n = 2000, DBDSVDMR3 is 96× faster than DBDSVDX on arithmetic spectra and 116× faster on the 32-decade log-uniform spectrum; it is 12–13× faster than DBDSDC on arithmetic spectra.Across 25 large practical matrices, median speedup over DBDSDC is 2.9×; with optimized BLAS, DBDSDC is faster for full SVD except deflation-free spectra, while DBDSVDMR3 remains faster for subset requests.Source: §4 · p. 12
  • On STCollection, DBDSVDMR3’s median error ratios are 4.2 and 6.6 versus 0.3–0.7 for LAPACK solvers, and it exceeds 100 on 15 of 125 matrices with tightly clustered eigenvalues.On the LAPACK test suite, median residual and orthogonality error ratios are 0.03 and 0.18, matching the LAPACK solvers.Source: §4 · p. 12

4.1. Experimental setup.

The experiments compare solver accuracy and efficiency using normalized residual and orthogonality errors across LAPACK-generated and STCollection matrices. Timing focuses on solver-call ratios under specified single-threaded setups.

4.1. Experimental setup.

  • The section establishes the experimental setup for comparing accuracy and efficiency across the solver test suites.The authors begin by summarizing the test setup.Source: §4.1 · p. 12

4.1.1. Notation and metrics.

  • Errors are measured by normalized residual and U- and V-orthogonality ratios; any ratio above 100 counts as a failure.Ratios near one indicate accuracy at floating-point roundoff level; DBDSVDMR3 is not required to deliver high relative accuracy for singular values.Source: §4.1.1 · pp. 12–13

4.1.2. Software versions.

  • Solver versions are drawn from the LAPACK development tree following release 3.12.1, with MR3GK compiled unmodified.Unless stated otherwise, all five solvers use the reference BLAS and run single-threaded.Source: §4.1.2 · p. 13

4.1.3. Test suite.

  • Tests use DCHKBD matrices and STCollection problems; accuracy comparisons use 122 shared STCollection matrices, while efficiency tests use 35 larger ones.DCHKBD provides 16 matrix classes; the STCollection accuracy set excludes three largest cases for DBDSQR because its vector computation did not finish within the time budget.Source: §4.1.3 · p. 13

4.1.4. Timing.

  • Timing measures solver calls after a warm-up and reports mean times, emphasizing solver ratios because absolute times drift on the virtualized host.The host uses one unpinned Intel Xeon core, reference BLAS, and gfortran 12 at -O2.Source: §4.1.4 · p. 13

4.2. Accuracy on LAPACK test matrices.

On the LAPACK test matrices, DBDSVDMR3 generally matches the LAPACK solvers, but loses accuracy on a tightly clustered Cholesky-derived variety. The black-box MR3GK and DBDSVDX also show distinct failure cases.

4.2. Accuracy on LAPACK test matrices.

  • DBDSVDMR3 matches DBDSDC and DBDSQR on 15 varieties, with median residual 0.03 and median orthogonality 0.18.Its distributions coincide with those solvers over most of the range.Source: §4.2 · p. 14, Figure 4.1 · p. 14
  • Ratios reach 20 for residual and 31 for orthogonality on the clustered Cholesky-derived geometric variety.Two-thirds of its singular values lie in one cluster with no relative gap above gaptol; long sequences of unguarded shifted child representations are associated with the loss of accuracy.Source: §4.2 · p. 14
  • MR3GK reaches orthogonality errors of 2 × 10^15 on varieties with small singular values, while DBDSVDX has a 9 × 10^13 orthogonality error on the log-uniform variety.MR3GK also does not terminate on the zero matrix; DBDSVDX’s large value occurs where half the singular values fall below ulp · ∥B∥.Source: §4.2 · p. 14

4.3. Accuracy on STCollection.

On STCollection, DBDSVDMR3 is less accurate in median ratios than the LAPACK solvers, and its largest errors occur on deeply clustered spectra. DBDSDC and DBDSQR have no failures in either suite, while DBDSVDX has specific failure modes.

4.3. Accuracy on STCollection.

  • DBDSDC has median residual and orthogonality ratios of 0.41 and 0.29; DBDSVDMR3’s are 4.2 and 6.6, while MR3GK’s are 5.7 and 16.DBDSQR has no ratio above 4.4 on any matrix, and DBDSVDMR3 does not reach MR3GK’s 10^13 tail on 10% of matrices.Source: §4.3 · p. 14, Figure 4.2 · p. 16
  • Table 4.1 reports worst ratios by matrix variety across five sizes, separating DBDSVDMR3 residual and U/V-orthogonality ratios from the other solvers’ largest ratio.Variety 7 is the Cholesky-derived geometric spectrum.Source: Table 4.1 · p. 15
  • On glued Wilkinson matrix T W21 g 1e+02, DBDSVDMR3’s residual and orthogonality ratios reach 8 × 10^5 and 9 × 10^5, while DBDSDC and DBDSQR stay below 3.The matrix repeats the 21 eigenvalues of W21 a hundred times, with coupling between 10^-14 and 10^14; deep shift sequences on clustered spectra are the solver’s main limitation, and it returns info = 0 without warning.Source: §4.3 · pp. 14–15
  • DBDSDC and DBDSQR fail on no matrix in either suite, whereas DBDSVDX fails on B bug414 and produces non-orthogonal vectors for singular values it cannot separate.On B bug414, underflow leads to info = 2n + 1; DBDSVDX also has orthogonality ratio 1.5 × 10^4 on B 16 and the log-uniform class.Source: §4.3 · p. 15

4.4. Timing.

DBDSVDMR3’s timing advantage depends on the spectrum, whether the full SVD or a subset is requested, and the BLAS implementation. It has the lowest worst-case normalized time on the LAPACK suite and performs especially well on large practical matrices and subset requests.

4.4. Timing.

  • On the LAPACK suite, DBDSVDMR3 was 12–13× faster than DBDSDC on arithmetic varieties at n=2000, but 1.0–4.1× slower on geometric and clustered varieties.It was 96× faster than DBDSVDX on the arithmetic variety at n=2000 and 116× faster on the 32-decade log-uniform variety.Source: §4.4 · p. 16
  • DBDSVDMR3 had the lowest maximum normalized time on the LAPACK suite, while DBDSDC had the lowest median and DBDSVDX lay an order of magnitude to the right.DBDSVDMR3’s maximum was 266 ns against DBDSDC’s 981 ns, and its total time at n=2000 was 4.2 s against 14.6 s.Source: §4.4 · p. 16, Figure 4.3 · p. 17

4.4.1. Optimized BLAS and subsets of the spectrum.

  • For the largest 10% of triplets, DBDSVDMR3 was faster than DBDSDC on 25 of 35 matrices and about as fast as DBDSVDX, with a median time ratio of 1.06.The subset cost was 4.7× lower than the full SVD; DBDSDC still won on the glued-Wilkinson family, which it deflated to a full SVD in 0.2 s.Source: §4.4 · p. 17, Figure 4.6 · p. 19
  • With single-threaded OpenBLAS, DBDSDC won on 31 of 35 practical full-SVD cases, while DBDSVDMR3 stayed ahead on spectra without deflation.DBDSDC became five times faster; the comparison used reference-BLAS timings as the control.Source: §4.4.1 · p. 17

4.4. Timing.

  • 199 ns was DBDSVDMR3’s median normalized time on 35 practical matrices, versus 711 ns for DBDSDC; it was faster on 25 matrices, with a 2.9× median factor.Across all 35, total times were 89 s for DBDSVDMR3 and 618 s for DBDSDC; the largest margins were on structural-engineering matrices, while DBDSDC won on matrices whose spectra deflate well.Source: §4.4 · pp. 16–17, Figure 4.4 · p. 18

4.4.1. Optimized BLAS and subsets of the spectrum.

  • The reference-BLAS results support MRRR’s advantage on spectra without deflation and for subset requests, while divide-and-conquer can win with substantial deflation on full spectra.The paper states the subset advantage assumes the subset path is merged; the merged implementation still computes all vectors regardless of RANGE.Source: §4.4 · p. 17, §4.4.1 · pp. 18–19

5. Recommendations for AI-assisted code development.

The authors describe an evaluator-guided workflow for adapting established numerical software with AI assistance. Human mathematical direction and expert review were central to the process.

5. Recommendations for AI-assisted code development.

  • The human supplied DSTEGR’s starting code, mathematical instructions, test matrices, and acceptance criteria; Claude Code wrote the Fortran, test harness, and test-suite extensions.AI-proposed changes were accepted only after accuracy and timing experiments were rerun on the same matrices.Source: §5 · p. 20
  • The human prescribed using TGK as the root representation, while the assistant proposed preprocessing, splitting, odd-block handling, null-vector handling, and removal of support-window truncation.The human checked the accepted modifications against the paper’s theory; a proposed random perturbation was rejected based on evidence.Source: §5 · p. 20
  • The authors attribute part of the effort’s success to an author’s familiarity with the original DSTEGR code and recommend AI-assisted coding for well-understood implementations.They chose DSTEGR rather than attempting to fix the Willems–Lang code or use DSTEMR.Source: §5 · p. 20

6. Availability and reproduction.

The solver and related software are in the reference LAPACK repository, while reproduction materials are available in the xmr-evolve repository.

6. Availability and reproduction.

  • DBDSVDMR3, its single-precision version, LAPACKE bindings, and dense [SDCZ]GESVDMR3 drivers are part of reference LAPACK.The subset path is the subject of a follow-up pull request.Source: §6 · p. 20
  • The xmr-evolve repository contains the test harness, matrix lists, raw logs, and scripts to regenerate the paper’s figures and tables.The repository is pinned by commit.Source: §6 · p. 20
Loading 2610.30227v1…