Linear Systems and Eigenvalue Problems: Open Questions from a Simons Workshop

Scope

• Across these topics, the document connects classical numerical analysis with randomized algorithms, finite-precision and bit-complexity concerns, and matrix problems arising in scientific computing and machine learning.The introduction identifies randomized matrix sketching as a major recent influence on basic numerical linear algebra.Source: §1 · p. 3

Workshop framework

• The workshop refined 55 proposed questions through discussions among diverse groups of 5–10 researchers into agreed-upon open-question formulations.The workshop was held October 6–10, 2025, at the Simons Institute for the Theory of Computing.Source: §1 · p. 3

Iterative solvers for linear systems

• Section 2 targets iterative solvers by seeking stronger benchmark problems and improved understanding of multigrid, conjugate-gradient, and GMRES convergence and numerical behavior.These solvers are widely used for large-scale simulations, including discretized fluid-flow problems, but important theoretical gaps remain.Source: §1 · p. 3

Eigenvalue problems

• Section 3 asks whether randomized QR convergence can be derandomized and pursues stable, efficient eigenvalue, singular-value, and Krylov-method analyses.The questions include deterministic perturbations that spread eigenvalue spectra, tridiagonal and bidiagonal computations, and Ritz-value distributions.Source: §1 · p. 3

Low-rank approximation and randomized sketching

• Sections 4 and 5 examine low-rank approximation and randomized embeddings, emphasizing structured matrix classes, column subset selection, injection-based sketches, sparse embeddings, and subsampled randomized Hadamard transforms.The low-rank questions cover greedy and randomized methods, hierarchical and block approximations, and matrices such as Laplacians and diagonally dominant matrices; injections can reduce sketch size while supporting effective least-squares approximation.Source: §1 · pp. 3–4

Other areas

• Section 6 extends the agenda to tensor decompositions with bounded error, eigenvector-dependent nonlinear eigenvalue problems, and practically computable matrix-function classes.The nonlinear eigenvalue model is based on the 2-local Hamiltonian problem, while the matrix-function questions concern fixed numbers of matrix-matrix products and approximation properties.Source: §1 · p. 4

2.1 Scalable iterative solvers for a set of parameterized benchmark problems

• Parameterized PDE discretizations are proposed as benchmarks because they combine broad scientific and engineering relevance with connections to graph-Laplacian research and systematic class generation.Unlike existing matrix collections, the proposed pipeline would generate samples across parameter ranges, exposing solver behavior across varied problem characteristics.Source: §2.1 · p. 5
• The agenda is to develop accepted benchmark-generation software and theoretically and practically fast algorithms for every parameterized linear-system class.The target is scalable convergence, ideally linear or nearly linear and independent of problem size, while extending beyond SDDM and SDD systems.Source: §2.1 · pp. 4–5
• The proposed PDE classes are intended to be tractable extensions of existing TCS techniques, not the most difficult multiscale or multiphysics PDE workloads.Some classes may nevertheless challenge NLA methods when robust, linear-time performance is required across wide parameter ranges.Source: §2.1 · p. 5
• Diffusion benchmarks vary the diffusion tensor, with large coefficient jumps and higher-order discretizations producing challenging matrices, including positive off-diagonal entries.Finite-difference and finite-element discretizations produce the linear systems; reference methods include Boomer-AMG and GENEO.Source: §2.1 · p. 5
• Convection-diffusion benchmarks generate nonsymmetric systems whose difficulty increases when convection causes positive off-diagonal entries and loss of diagonal dominance.Centered differences can become unstable for large convection coefficients, whereas one-sided upwind differences restore diagonal dominance; irregular geometries may require finite elements.Source: §2.1 · p. 6
• For convection-diffusion systems, GMRES convergence depends on eigenvalues and eigenvector conditioning, limiting direct use of Krylov-based TCS techniques and motivating stationary-iteration analyses.The target solution accuracy should match the discretization error, O(h2), in the same norm.Source: §2.1 · p. 7
• Helmholtz discretizations become indefinite when solving A − k2I for wave numbers whose matrices have several negative eigenvalues, making iterative solution challenging.The standard approach combines an iterative method with a shifted-Laplacian preconditioner whose real parameters move the preconditioned spectrum away from zero.Source: §2.1 · pp. 7–8

Motivation

• The motivation is to extend multigrid’s optimal O(n) solver complexity from specialized boundary-value discretizations to general algebraically specified systems resembling elliptic PDE discretizations.A proposed route is to establish smoothing and approximation properties using suitably chosen smoothers and prolongators.Source: §2.2 · pp. 8–9

Matrix class

• SWCDDM matrices are SDDM matrices whose weakly diagonally dominant rows connect through nonzero A-entries to strictly diagonally dominant rows.SDDM additionally requires symmetry, diagonal dominance, and non-positive off-diagonals.Source: §2.2 · p. 8

Algebraic multigrid analysis

• The central open problem is to compute, for polynomially bounded SWCDDM matrices, an implicit symmetric Z satisfying Ω(1)M^-1 ⪯ Z ⪯ O(1)M^-1 in ˜O(nnz(M)) time.The operator should be applicable in O(nnz(M)) time, with constants independent of n.Source: §2.2 · pp. 8–9

Variants

• Variants restrict the problem to planar matrices or trade construction time against operator-application time using poly(log log(nnz(M))) factors.The stated variants consider faster implicit construction or slightly slower application while retaining the core approximation goal.Source: §2.2 · p. 9

Convergence framework

• The standard convergence framework requires level-wise smoothers, prolongation, restriction, and coarse-grid correction satisfying smoothing and approximation properties.For elliptic model problems, these conditions yield convergence uniform in mesh size and extend to multilevel V-cycles with O(n) complexity.Source: §2.2 · pp. 9–10

Conjugate gradient versus sketch-and-project

• The central open question is which eigenvalue conditions let sketch-and-project reach ε accuracy in fewer epochs than conjugate gradient [GR15, HS52].A sketch-and-project epoch processes n randomly selected rows, whereas a Krylov epoch performs one matrix–vector multiplication over all n rows.Source: §2.3 · p. 10
• The simplest proposed comparison is between unpreconditioned conjugate gradient and randomized coordinate descent, despite possible preconditioning, batching, and momentum accelerations [LR25] [GHRS18].RCD applies to a strictly positive definite system, samples index i with probability A(i,i)/tr(A), and enforces the selected row exactly.Source: §2.3 · p. 10
• The unresolved stylized problem asks how CG and RCD stopping times compare asymptotically for randomly rotated positive-definite matrices with eigenvalues 1^-p, 2^-p, ..., n^-p.The right-hand side is generated as b_n=A_n z_n, with z_n uniform on the sphere, and the limit is taken as n→∞.Source: §2.3 · p. 11
• An August 20, 2026 update says Iteris claims to solve Problem 2.4, but the broader motivating question remains open.That broader question covers general spectra, right-hand sides, preconditioning, and finite precision.Source: §2.3 · p. 11, §4.1.2 · p. 32

Outlying singular values

• Clustered spectra with outlying singular values can make the overall condition number a poor complexity indicator, motivating the (n, k, l, κ) model that separates top, bottom, and clustered singular values.The model permits one-sided outliers by setting k or l to zero and generalizes earlier Krylov analyses such as.Source: §2.4 · pp. 11–12

Large outlying singular values

• For large outliers, direct-access randomized sketching and matrix-vector Krylov methods exhibit a complexity separation, while sparse systems motivate seeking a ˜O(kω + nnz(A)κ) solver.Krylov methods require ˜O(k + κ) matrix-vector queries, with an Ω(k) lower bound, whereas sparse matrix-vector costs are ˜O(nnz(A) · k); the proposed direct-access target is Problem 2.5.Source: §2.4 · p. 12

Averaged condition numbers

• Averaged condition numbers may substantially reduce complexity—for example, σ_i = 1/i gives κ = n but ¯κ = O(log n)—yet the smallest usable exponent p remains open across dense and sparse access models.[LS13a] replaces κ with ¯κ for psd systems; [DS25] achieves p > 1/2 generally and p > 1/4 for psd systems when T_A = n^2, while finite-p guarantees for T_A = nnz(A) remain unclear.Source: §2.4 · p. 13

Small outlying singular values

• Small outlying singular values appear harder to exploit: Krylov methods need ˜O(lκ) matrix-vector products, and the direct-access complexity for (n, 0, l, κ) systems remains unknown.This regime is motivated by PDE spectra and certain random matrices, whereas no improvement over the stated Krylov bound is known for dense or sparse direct access.Source: §2.4 · p. 13

Solving a random linear system

• For random Gaussian, Bernoulli, or Rademacher systems, excluding the bottom l singular values improves the condition number to O(n/l), but it remains open whether high-accuracy solution is possible faster than O(nω).The related matrix-vector-product question asks whether Ω(n) queries are necessary, even under classical matrix multiplication with ω = 3.Source: §2.4 · p. 14

Stability and finite-precision arithmetic

• Finite-precision stability and bit complexity remain unresolved: Krylov guarantees may require extremely high precision, while a 2026 result nearly attains the κ = O(1) target only for dense normal equations with explicit B access.The result gives bit complexity n^o(1) ˜O(k^(ω+η) + n^2) for every fixed η > 0; more general systems remain open, and sparse Word-RAM Krylov bounds are still worse than Real-RAM bounds [PV21, Nie22].Source: §2.4 · pp. 14–15

Preconditioning framework

• Preconditioning accelerates CG and GMRES by replacing Ax = b with a modified system using an easily computable approximate inverse M −1.Effective preconditioners should approximate A−1 while keeping the action of M −1 inexpensive; left, right, and split forms have different effects.Source: §2.5 · p. 15

CG convergence and spectral localization

• For CG, reducing the condition number does not necessarily accelerate convergence because the method adapts nonlinearly to the entire eigenvalue spectrum.For elliptic PDE model problems, [GMNS19] [GNS20] obtain exact eigenvalue localization under Laplacian preconditioning; extending this approach to other PDEs or preconditioners is open.Source: §2.5 · pp. 15–16

Nonnormal GMRES

• GMRES convergence for nonnormal systems may depend on eigenvector conditioning and transient behavior, motivating an open search for descriptive early-iteration convergence bounds.The analysis may require ϵ-pseudospectra or the field of values rather than eigenvalues alone.Source: §2.5 · p. 16

Eigenvalue clustering

• An open problem is to identify nonsymmetric model problems where m eigenvalue clusters in M −1A rigorously correspond to acceptable GMRES solutions in O(m) steps.The goal is to determine which problem properties make eigenvalue-clustering predictions match actual convergence.Source: §2.5 · p. 16

Convection-diffusion GMRES

• For SUPG convection-diffusion discretizations, the post-transient GMRES acceleration remains unexplained, including whether it follows from eigenvalue locations or eigenvector conditioning.Field-of-values analysis can establish scalability in some cases but typically does not provide tight iteration-count estimates; explains how boundary conditions affect the initial slow phase.Source: §2.5 · p. 17

Scalability

• A central open direction is a fully scalable preconditioner whose iteration count is nearly dimension-independent and almost insensitive to physical parameters such as small viscosity coefficients.Practical solvers should use few iterations relative to matrix size and achieve a reasonable, not machine-roundoff-level, error or residual norm.Source: §2.5 · p. 17

• Finite-precision analyses must determine convergence and attainable accuracy for short-recurrence Krylov methods applied to SPD systems or matrix-function evaluations [Che24, Chapter 4].Prior results include.Source: §2.6 · pp. 17–18
• Because recursively updated residuals are not true residuals and may shrink below machine precision, Problem 2.15 asks when this occurs for CG or steepest descent.The problem also seeks counterexamples where the updated residual norms do not drop below machine precision.Source: §2.6 · p. 18
• Problem 2.16 asks whether standard good implementations of CG and Lanczos behave similarly in finite precision for SPD systems, despite their unexplained implementation-level relationship.In exact arithmetic, CG can be viewed as computing the Cholesky factorization of Lanczos’s tridiagonal matrix.Source: §2.6 · p. 18
• Problem 2.17 asks how many precision bits guarantee normwise backward error at most ϵ for CG solving an n by n SPD system within n steps.Source: §2.6 · p. 18
• Problem 2.18 seeks Lanczos analyses for inexact matrix-vector products and preconditioner applications that generalize Paige’s results [Pai76, Pai80] and clarify algorithm behavior.Existing analysis does not cover preconditioners stored in memory or applied as implicit, inexact operators; generalized models use perturbations such as ∆A and ∆M.Source: §2.6 · pp. 18–19
• Problems 2.19 and related work seek deeper finite-precision understanding through orthogonal-polynomial connections and a Paige-style analysis of block Lanczos [Pai76, Pai80].The block case must address rank-deficient blocks and deflation, while the orthogonal-polynomial view links forward instability to an ill-conditioned moments-to-tridiagonal-matrix bijection [Kni96, CT24, Che24].Source: §2.6 · p. 19

The Forsythe conjecture

• The conjecture remains largely open despite extensive use of steepest descent, gradient descent, and CG; [FLT23] surveys results through 2023 and extends the framework to nonsymmetric and cross iterations.For certain orthogonal matrices with s = 1, the cross iteration exhibits similar limiting behavior.Source: Scribe: J¨org Liesen · pp. 20–21
• The restarted method converges to the exact solution under the stated nontriviality assumptions 1 ≤s < d(A) and d(A, r0) ≥s + 1.The grade condition ensures the next residual is nonzero and preserves the corresponding lower-bound condition for subsequent iterations.Source: Scribe: J¨org Liesen · p. 20
• For 2 ≤s < d(A), each even and odd subsequence of normalized residuals is conjectured to have a single limit vector.The iteration repeatedly minimizes the A-norm error over an s-dimensional Krylov subspace, equivalently performing restarted CG in exact arithmetic.Source: Scribe: J¨org Liesen · p. 20
• If true, the conjecture would make even and odd residuals asymptotically collinear, so residual or error norms could converge to zero at best linearly.Forsythe established the two-direction cycling for s = 1, while numerical evidence suggests analogous behavior for larger restart lengths.Source: Scribe: J¨org Liesen · pp. 19, 21
• Three 2026 references claim to have solved Problem 2.20 for s = 2, including work by Colbrook, Stepaniants, and Townsend, Liesen, and Peng.The supplied update lists these references but does not provide their proofs or detailed results.Source: Scribe: J¨org Liesen · p. 21

3.1 Deterministic pseudospectral shattering / Minami bound

• A deterministic perturbation algorithm remains open for ensuring polynomially bounded diagonalization conditioning within O(n^3 log^c(n/δ)) operations and comparable precision.Random perturbations achieve polynomial or weaker bounds with high probability, while deterministic regularization would replace the randomness used by current provable algorithms.Source: §3.1 · pp. 21–22
• A related open problem seeks an efficient deterministic diagonal perturbation giving Minami-type eigenvalue gaps for Hermitian tridiagonal matrices, including faster bounds for Krylov methods.The target is O(n^3 log^c(n/δ)) for diagonalization and ambitiously O(nnz(A) log^c(n/δ)) or O(nnz(A) log^c(n)/δ^c) for Krylov methods.Source: §3.1 · p. 23

3.2 Accuracy vs. precision in efficient diagonalization algorithms

• Problem 3.3 asks for randomized diagonalization or Schur-form algorithms using O(n^3 log^c(n/δ)) time and only O(log(n/δ)) bits of precision.Existing methods either require more precision, as in and [BGVS22], or use the target precision with substantially worse runtime, as in [ABB+18]; inverse-free sign-function methods may achieve the target precision, but end-to-end floating-point analysis remains open.Source: §3.2 · pp. 23–24

3.3 Ritz values and approximate invariant subspaces

• Ritz-value theory asks when augmenting an approximate invariant subspace with arbitrary extra vectors preserves eigenvalue approximations, and how Ritz values distribute over W(A) for non-Hermitian matrices.Problem 3.4 concerns Q = bV + W and eigenvalues of Q*AQ; Problem 3.6 asks deterministic or probabilistic distribution questions for arbitrary, Krylov, and random subspaces.Source: §3.3 · pp. 24–25
• A 2026 construction claims negative answers to both questions in Problem 3.5, ruling out uniform poly(n)ε^β bounds for every β > 1, while normal or polynomially conditioned cases remain open.The construction uses real diagonalizable matrices with exponentially growing κV(A), so it does not settle the question under normality or polynomial eigenvector conditioning.Source: §3.3 · p. 25

3.4 Revisiting MR3

• MR3 lacks a fully rigorous, readily verifiable floating-point correctness proof for arbitrary symmetric tridiagonal matrices, and adapting it reliably to bidiagonal SVD remains unresolved.Independent MR3 applications can produce large residuals, while Golub–Kahan formulations can lose orthogonality; prior modified implementations were abandoned or deemed unreliable [GL03, WLV06, WL12a, WL12b, WL13, MDV20].Source: §3.4 · pp. 26–27

3.5 Largest / rightmost eigenvalue

• For a general nonnormal A, the section asks whether poly(log(nκV), 1/ε) matrix-vector products suffice to approximate every eigenvalue near the spectral radius, with a rightmost-eigenvalue variant.The power method supplies a mixed forward-backward approximation in O(log(nκV)/ε) iterations, but converting this into the requested eigenvalue guarantees remains open.Source: §3.5 · p. 28

4 Low-rank approximation and index selection

• The broader framework defines CSSP as selecting k columns for a good projection approximation, while rank-revealing factorizations balance approximation quality with computational cost, interpretability, and structural preservation.The truncated SVD is optimal with µ(m,n,k) = 1, but GE- and QR-based methods can be cheaper, more interpretable, and better at retaining matrix structure [MG03, GE96].Source: §4 · pp. 28–29

4.1 When do greedy/randomized column selection algorithms work?

• Greedy and randomized methods often perform far better than worst-case bounds suggest, motivating structural assumptions under which they can provably achieve quasi-optimal CSSP or RRQR guarantees.CPQR has a worst-case factor 2^k√(n−k), while related exponential factors become polynomial for positive definite Lipschitz-smooth kernels or doubly diagonally-dominant matrices [GE96, JT25, CKM20].Source: §4.1 · pp. 29–30

4.1.1 The discrete Lehmann representation

• For the discrete Lehmann representation, GECP empirically follows a sharper low-rank rate than existing theory, especially when Λ is as large as 10^5 or 10^6, motivating stronger kernel-structure-based bounds.The approximation is accurate when selected time and frequency nodes form a good cross-approximation to the fermionic kernel; Problem 4.2 asks for stronger GECP guarantees.Source: §4.1.1 · p. 31
• GECP for the fermionic kernel is claimed to achieve k = O(log(Λ) log(1/ε)) and ||K − K̂||∞ ≤ ε, matching the optimal rank scaling up to universal constants.The same bound is stated for finite matrices obtained by sampling the kernel at distinct time and frequency nodes.Source: §4.1.1 · p. 32

4.1.2 Orthogonal rows

• For orthonormal-row matrices, standard QRCP is posed as an open question because existing bounds are exponential in k, while a related Bischof–Stewart pivoting strategy achieves the desired bound.The update reports that Theorem 2 of Chen et al. claims no polynomial bound for standard QRCP, whereas Damle reports the positive result for the alternative pivoting strategy.Source: §4.1.2 · p. 32

4.1.3 Incoherence and random row subsampling

• For kernel matrices with many more rows than columns, the open problem is to identify efficiently verifiable conditions on f under which O(k) random row samples preserve effective column selection.The motivation is reducing direct O(nmk) interpolative-decomposition cost to O(k^2m) when |I| = 5k, without requiring computationally prohibitive sparsity or incoherence verification.Source: §4.1.3 · p. 33

4.1.4 Block-column subset selection for rank-structured matrices

• Block-column subset selection seeks quasi-optimal hierarchical approximations while adaptively choosing partitions whose structure reflects compressibility, including anisotropic transport geometry in hyperbolic PDEs.The open question asks which PDE classes admit practical efficient algorithms and the optimal dependence on partition parameters; fixed distance-based hierarchies can fail for characteristic-manifold structure.Source: §4.1.4 · p. 34

4.1.5 Inverses of Laplacians and related matrices

• For inverse Laplacians, nuclear Nyström error is submodular over nonempty index sets, yielding greedy guarantees, while its extension to positive-definite SDDM and SDD matrices remains to be proved or disproved.The SDDM extension is expected to be straightforward, whereas the positive-definite SDD case is less clear and may use reductions from SDD systems to Laplacians.Source: §4.1.5 · p. 35
• Problem 4.6 receives a positive answer for (a) and a negative answer for (b), via a 3×3 parametrized counterexample for positive definite SDD matrices.Source: §4.1.5 · p. 35

4.2 Clarifying relations between volume sampling and optimal subset selection

• Problem 4.7 asks for tight bounds comparing worst-case optimal column subset-selection trace error with the objective produced by volume sampling, for matrices with prescribed eigenvalues λ.Source: §4.2 · p. 36
• A simple tightness bound would clarify limits of spectrum-only column subset-selection bounds and strengthen interpretation of the relation between greedy optimization and volume sampling established in.Source: §4.2 · p. 36
• For k = n −1, Schur complements and the Schur-Horn theorem show that xn−1(λ) = yn−1(λ) for a symmetric matrix K whose inverse has equal diagonal entries.Source: §4.2 · p. 36
• For other k, the relationship between the two quantities remains unclear, although limited empirical investigation finds them usually quite close.The variance and other higher-order statistics of volume sampling are also identified as related open questions.Source: §4.2 · p. 36

5.1 Injections vs. embeddings

• The section asks whether weaker oblivious subspace injection assumptions can guarantee accurate sketch-and-solve algorithms, including ℓp regression, rather than requiring restrictive subspace embeddings.The meta-problem targets the weakest assumption ensuring a desired accuracy level; a p-isotropic ℓp-OSI is reported to yield constant-factor guarantees for 1 ≤p < ∞.Source: §5.1 · pp. 37, 39
• Updates report negative answers for relative-error guarantees from oblivious subspace injections in sketch-and-solve least squares and randomized SVD, while the weakest additional condition remains open.Townsend and Wang [CEMT25?] construct a 2 × 1 least-squares counterexample and a 2 × 2 randomized-SVD counterexample.Source: §5.1 · p. 38

5.2 Sparse dimensionality reduction maps

• Sparse-map questions examine whether SparseStack achieves OSE or OSI guarantees with near-linear output dimension and low sparsity, including ζ = 4 and k = 2r.The supplied account claims the Nelson–Nguyen OSE conjecture is settled affirmatively with ζ = O((log r)/ε) and k = O(r/ε2) for constant failure probability, while a broader result rules out ζ = 4, k = 2r, and constant α for OSI.Source: §5.2 · pp. 40–41

5.3 Subsampled randomized trigonometric transforms

• Traditional SRHTs require embedding dimension O((r + log n) log r), motivating the open question whether rerandomized SRHTs achieve ε-accurate OSE with k = O(r/ε2).The theorem also gives an r-dimensional subspace where traditional SRHTs fail to be injections when k = o(r log r), despite empirical O(r) scaling for rerandomized variants.Source: §5.3 · p. 42

6.1 Optimal bounds for Tensor Train decomposition

• Tree tensor decompositions currently inherit an (m−1)-factor error bound from sequential SVD approximations, motivating algorithms that improve it or establish its optimality.The bound is tight for HOSVD in some Tucker cases, and no polynomial-time general improvement is known.Source: §6.1 · p. 43

6.2 Eigenvalue systems arising from Hamiltonians

• Hamiltonian tensor-product minimization can be reformulated as a dimension-nd nonlinear eigenvalue problem, but efficient methods with superlinear local convergence and iteration complexity linear in d remain open.Existing general quadratically convergent NEPv methods require an nd×nd Hessian, while the minimizer corresponds to the NEPv solution.Source: §6.2 · pp. 44–45

6.3 Matrix sign function

• Matrix-sign approximation targets GPU-friendly computation using matrix multiplications, additions, and scalings for spectra in [−1,−δ]∪[δ,1], with applications including eigensolvers, matrix approximation, and electronic-structure purification.For unrestricted polynomials, degree and multiplication-cost bounds are known, but exploiting lower multiplication complexity for some polynomials could improve the estimate.Source: §6.3 · pp. 46–47