PaperlayerView exact revision on arXiv

Source-linked AI summary

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

Noah Amsel, Yves Baumann, Paul Beckman, Peter Bürgisser, Chris Camaño, Tyler Chen, Edmond Chow, Anil Damle, Michal Derezinski, Mark Embree, Ethan N. Epperly, Robert Falgout, Mark Fornace, Anne Greenbaum, Chen Greif, Diana Halikias, Zhen Huang, Elias Jarlebring, Yiannis Koutis, Daniel Kressner, Rasmus Kyng, Jörg Liesen, Jackie Lok, Raphael A. Meyer, Yuji Nakatsukasa, Kate Pearce, Richard Peng, David Persson, Eliza Rebrova, Ryan Schneider, Rikhav Shah, Edgar Solomonik, Nikhil Srivastava, Alex Townsend, Robert J. Webber, Jess Williams

arXiv:2602.05394v3math.NAcs.DS

TL;DR

The paper asks how to formulate and analyze open problems in numerical linear algebra at the intersection of theoretical computer science and numerical analysis. It organizes workshop-generated questions across linear solvers, eigenvalues, low-rank approximation, sketching, and related topics, while highlighting unresolved limitations and recent progress.

  • Problem

    The workshop sought well-defined research questions about the complexity and numerical solution of linear algebra problems across theoretical computer science and numerical analysis.

  • Method

    Working groups from diverse research areas organized and analyzed open questions spanning iterative solvers, eigenvalue computation, low-rank approximation, randomized sketching, and broader matrix-computation topics.

  • Results

    The document records progress on sparse-embedding Problem 5.4, including a full two-sided OSE with ζ = O((log r)/ε) and k = O(r/ε2) for constant failure probability.

  • Takeaways & Limitations

    The collection identifies concrete directions for extending scalable algorithms and clarifying guarantees across numerical linear algebra, while documenting open bounds and derandomization questions.

  • Takeaways & Limitations

    For nonlinear eigenvalue problems, optimization methods have limited study of their complexity and convergence properties, while the generalization of randomized MR3 guarantees to singular values remains unresolved.

Abstract

from arXiv · show

This document presents a series of open questions arising in matrix computations, i.e., the numerical solution of linear algebra problems. It is a result of working groups at the workshop Linear Systems and Eigenvalue Problems, which was organized at the Simons Institute for the Theory of Computing program on Complexity and Linear Algebra in Fall 2025. The complexity and numerical solution of linear algebra problems is a crosscutting area between theoretical computer science and numerical analysis. The value of the particular problem formulations here is that they were produced via discussions between researchers from both groups. The open questions are organized in five categories: iterative solvers for linear systems, eigenvalue computation, low-rank approximation, randomized sketching, and other areas including tensors, quantum systems, and matrix functions. (Updated to reflect the status of the open problems as of August 20, 2026.)

1 Introduction

The workshop assembled well-defined open questions at the intersection of theoretical computer science and numerical analysis, spanning iterative solvers, eigenvalue computation, low-rank approximation, randomized sketching, and broader matrix-computation topics.

  • Workshop scope: The working groups refined 55 proposed questions into well-defined open questions through discussions among researchers from diverse research areas.The workshop took place October 6–10, 2025, at the Simons Institute.
  • Iterative solvers: Open questions on iterative linear solvers target benchmark problems and the numerical behavior and convergence of multigrid, conjugate gradient, and GMRES methods.The focus is large-scale simulation systems, including discretized partial differential equations.
  • Eigenvalue computation: Eigenvalue questions include derandomizing a provably convergent randomized QR variant, analyzing tridiagonal and bidiagonal problems, and understanding Krylov-method convergence through Ritz-value distributions.The randomized QR approach addresses a longstanding open problem for general matrices.
  • Low-rank approximation: Low-rank approximation questions examine greedy and randomized methods, structured matrix classes, column subset selection, and hierarchical or block approximations.The discussion includes Laplacians and diagonally dominant matrices, with tensor approximations revisited later.
  • Randomized sketching: Randomized-sketching questions study injections, sparse embeddings, and structured embeddings, including their sample complexity and approximation properties.Injections relax standard subspace-embedding requirements while supporting effective approximations for problems such as linear least squares.
  • Other areas: Additional questions concern tensor decompositions, eigenvector-dependent nonlinear eigenvalue problems in quantum chemistry, and matrix functions computable with a fixed number of matrix-matrix products.The matrix-function questions also examine approximation properties of this practically relevant class.

2.1 Scalable iterative solvers for a set of parameterized benchmark problems

This section proposes parameterized PDE discretizations as benchmark linear systems for connecting theoretical and numerical linear algebra, then identifies solver challenges across diffusion, convection-diffusion, and Helmholtz problems.

  • Benchmark design: Parameterized PDE discretizations are proposed as benchmark systems that are relevant in practice, suitable for theoretical analysis, and connected to graph-Laplacian research.The proposed benchmark pipeline is intended to broaden the problem classes addressed by theoretical computer science.
  • Benchmark design: The two central problems are to generate accepted parameterized benchmark systems and to develop theoretically and practically fast algorithms for each class.The benchmark classes should support analysis across varying problem parameters.
  • Diffusion problems: Diffusion benchmarks become challenging when the diffusion tensor has large jumps or when higher-order discretizations produce positive off-diagonal entries.These matrix properties can challenge existing graph-based approaches.
  • Convection-diffusion problems: Convection terms create nonsymmetric systems, and large convection coefficients can produce positive off-diagonal entries and loss of diagonal dominance.Upwind differences can restore diagonal dominance by choosing the one-sided direction according to the convection signs.
  • Convection-diffusion problems: For nonsymmetric systems, GMRES convergence depends on eigenvector conditioning, while efficient short-recurrence methods such as conjugate gradient are unavailable.The desired solution accuracy should be comparable to the discretization error, O(h^2), in the same norm.
  • Helmholtz problems: Helmholtz discretizations can produce indefinite matrices with several negative eigenvalues, making iterative solution challenging; shifted-Laplacian preconditioning is a developed approach.The shifted Laplacian moves the preconditioned spectrum away from zero using real parameters with a complex shift.

2.2 Correctness of multigrid methods beyond the standard model problems

The section asks whether multigrid guarantees and near-linear complexity can extend beyond standard model problems, including polynomially bounded SWCDDM matrices and broader elliptic discretizations.

  • Algebraic multigrid: The multigrid extension problem seeks correctness guarantees beyond standard model problems for polynomially bounded symmetric weakly-chained diagonally dominant matrices.The matrix class is defined using diagonal dominance, non-positive off-diagonals, and chains from weakly to strictly dominant rows.
  • Algebraic multigrid: The target is an implicitly represented symmetric operator Z computable in O~(nnz(M)) time, applicable in O(nnz(M)) time, and spectrally comparable to M^-1 with size-independent constants.The required comparison is Ω(1)M^-1 ⪯ Z ⪯ O(1)M^-1.
  • Algebraic multigrid: Variants ask for the same result on planar matrices, with faster construction, or with slightly slower application of the implicit operator.The proposed construction and application times vary through poly(log log(nnz(M))) factors.
  • Standard model problems: For standard elliptic model problems, smoothing and approximation properties define the two-grid framework, whose multilevel extension achieves optimal O(n) complexity.The convergence rate is uniform with respect to mesh size h when the smoothing-step count is independent of h.

2.3 Conjugate gradient versus sketch-and-project

This section asks when sketch-and-project methods can outperform conjugate gradient, using epoch-based comparisons and a stylized random positive-definite model while leaving the broader characterization open.

  • Comparison framework: The motivating question is which eigenvalue conditions let sketch-and-project reach ε accuracy in fewer epochs than conjugate gradient.An epoch processes randomly selected rows for sketch-and-project and all rows through one matrix-vector multiplication for a Krylov method.
  • Randomized coordinate descent: Randomized coordinate descent selects coordinate i_t with probability A(i_t,i_t)/tr(A) and updates that coordinate to satisfy its row equation exactly.This defines the coordinate-update method used in the comparison setup.
  • Stylized model: The stylized open problem compares the asymptotic stopping times of conjugate gradient and randomized coordinate descent for random positive-definite matrices with polynomially decaying eigenvalues.The matrices use Haar-unitary eigenvectors and eigenvalues 1^-p, 2^-p, ..., n^-p, with right-hand sides generated from uniformly spherical vectors.
  • Status update: An August 20, 2026 update records a reference claiming to solve the stylized problem, while the broader motivating question remains open.The update identifies the reference as Iteris: Agentic Research Loops for Computational Mathematics.
  • Comparison framework: The broader open problem covers general spectra, right-hand sides, preconditioning, and finite precision beyond the stylized comparison.The section also notes acceleration mechanisms such as batching rows and momentum.

2.4 The effect of outlying singular values (or eigenvalues)

This section asks how singular-value distributions beyond the overall condition number affect linear-system complexity, especially when spectra contain outliers. It contrasts direct-access and matrix-vector-product models and identifies unresolved complexity and stability questions.

  • Open complexity questions: The open problems seek natural characterizations of linear-system complexity in terms of eigenvalue or singular-value distributions.This includes systems with large or small outlying singular values and random-matrix test cases.
  • Spectral structure: Outlying singular values can make the overall condition number a poor complexity indicator, while iterative methods may remain effective for large-condition-number systems.The section models spectra as top and bottom outliers surrounding a cluster with condition number κ.
  • Computational models: Direct-access algorithms and matrix-vector-product methods can have different complexities for systems with large outlying singular values.In the matrix-vector-product model, Krylov methods require ˜O(k + κ) queries, while lower bounds give Ω(k) queries; dense systems can incur Ω(n2k) operations.
  • Computational models: For sparse matrices, the Krylov cost is ˜O(nnz(A) · k), motivating whether direct-access randomized methods can achieve comparable input-sparsity dependence.The associated open question asks for ˜O(kω + nnz(A)κ) arithmetic operations.
  • Averaged condition numbers: The smallest p enabling complexity ˜O(TA · ¯κp) remains open, although known results achieve p > 1/2 for general systems and p > 1/4 for positive semidefinite systems when TA = n2.For TA = nnz(A), the existence of such a result for finite p is unclear.
  • Stability and precision: Finite-precision stability and bit complexity remain concerns for Krylov and randomized sketching methods, especially when guarantees rely on outlying singular values.A randomized result nearly attains the target for dense normal equations with κ = O(1), while the general case remains open.

2.5 Preconditioning and convergence of iterative solvers

This section studies preconditioning and convergence for iterative solvers, emphasizing scalable methods whose iteration counts remain stable across dimensions and physical parameters. It also highlights that eigenvalue clustering alone does not generally explain nonsymmetric solver behavior.

  • Preconditioning: Preconditioning aims to accelerate iterative methods by replacing Ax = b with a system whose inverse preconditioner is both effective and easy to apply.The section distinguishes left, right, and split preconditioning, while noting that their differences are outside its scope.
  • Nonnormal convergence: For nonnormal systems, GMRES convergence can depend on eigenvectors and quantities such as pseudospectra or the field of values, not only eigenvalue locations.The proposed problems seek descriptive transient bounds and model cases where eigenvalue clustering corresponds to actual convergence.
  • Nonnormal convergence: Open questions ask whether eigeninformation explains convergence acceleration after the transient phase in convection-diffusion problems.The requested eigeninformation includes eigenvalue location and eigenvector conditioning.
  • Scalable preconditioning: A central goal is a fully scalable preconditioned solver whose iteration count is nearly independent of matrix dimensions and problem parameters.The target concerns systems from PDE discretizations and practical error or residual tolerances.
  • Scalable preconditioning: One proposed direction is a problem-dependent preconditioner that is nearly insensitive to changes in physical parameters such as small viscosity coefficients.This is posed as a research direction for robust scalable solution procedures.

2.6 Symmetric Krylov methods in finite precision arithmetic

This section examines how symmetric Krylov methods behave in finite precision, where exact-arithmetic equivalences and residual recurrences may fail. It proposes analyses for inexact operations, block methods, and the relationship between CG and Lanczos.

  • Finite-precision behavior: Finite-precision behavior of short-recurrence Krylov methods can differ substantially from exact-arithmetic behavior, leaving convergence and attainable accuracy unresolved.The questions concern SPD solves and matrix-function-vector products.
  • Residuals: Computed residual updates may not equal the true residuals b − Axk in finite precision, motivating conditions under which their norms fall below machine precision.The associated problem asks for guarantees or counterexamples.
  • CG versus Lanczos: Although CG and Lanczos are equivalent in exact arithmetic, their relationship and similarity of behavior in finite precision remain unexplained.The proposed comparison depends on using standard good implementations of both methods.
  • Precision requirements: A proposed problem asks how many precision bits guarantee normwise backward error at most ϵ for CG on an n by n SPD system within n steps.This question directly targets precision requirements for finite-precision CG.
  • Lanczos analysis: Existing Lanczos finite-precision analysis does not account for preconditioners or more general inexact matrix-vector products.The proposed extension models operations through perturbations such as (A + ∆A)x and inexact solutions involving M + ∆M.
  • Block Krylov methods: Extending Paige-style analysis to block Lanczos is difficult because blocks may become rank deficient and require deflation.The open problem allows assumptions about within-block orthogonalization or analysis of a standard accurate implementation.

2.7 The Forsythe conjecture

The Forsythe conjecture concerns whether normalized residual directions in restarted conjugate-gradient iterations settle into two limiting vectors for restart lengths 2 ≤ s < d(A). If true, the restarted method's residual norms can converge at best linearly, and the conjecture was largely open in general despite numerical evidence and reported progress for s = 2.

  • Iteration: Forsythe analyzed the optimum s-gradient method, which in modern terms is CG restarted every s steps for a symmetric positive definite system.Each restarted iterate is the A-norm best approximation over a Krylov subspace.
  • Iteration: The restarted method converges to x = A−1b under the stated setup, while the conjecture concerns the limiting directions of its normalized residuals.The iteration repeatedly restarts after an s-step Krylov update.
  • Conjecture: The Forsythe conjecture states that, for 2 ≤ s < d(A), the even and odd normalized-residual subsequences each have a single limit vector.The conjecture concerns limiting directions rather than convergence speed.
  • Consequence: If the conjecture holds, even-indexed and odd-indexed residual directions become separately asymptotically collinear, so residual norms can converge to zero at best linearly.This consequence applies to the restarted iteration's error or residual norms.
  • Status: The conjecture remained largely open for general restart lengths despite the popularity of steepest descent, gradient descent, and CG.Numerical evidence suggested it was true, and later references reported claims addressing s = 2.

3 Eigenvalue problems

This section frames eigenvalue computation around randomized perturbations that regularize diagonalization, then asks whether comparable guarantees, runtimes, and precision can be achieved deterministically or more efficiently. It also identifies open difficulties for invariant-subspace Ritz approximations and SVD-related eigensolver adaptations.

  • Pseudospectral shattering: Random perturbations make the eigenvalue condition number polynomially bounded with high probability, enabling provable diagonalization and Schur-form algorithms.Replacing A by A + E is called pseudospectral shattering; existing provable algorithms rely on this regularization.
  • Pseudospectral shattering: Deterministic pseudospectral shattering asks for an efficient perturbation E with ∥E∥≤δ and κeig(A+E)≤C(δ/n)^c.The target runtime is O(n^3 log^c(n/δ)) arithmetic operations, with an ambitious finite-precision version using O(log^c(n/δ)) bits.
  • Pseudospectral shattering: Sparse complex Gaussian perturbations yield the weaker bound κeig(A+E)≤exp(log^2(n/δ)) with high probability, which remains useful for algorithms depending logarithmically on κeig.Dense complex Gaussian perturbations are known to obtain the stronger bound with c≤5 with high probability.
  • Precision and derandomization: A separate open problem seeks randomized diagonalization or Schur-form algorithms with O(n^3 log^c(n/δ)) runtime and only O(log(n/δ)) bits of precision.Existing approaches use more precision, while inverse-free replacements may recover the desired precision bound.
  • Ritz values and SVD: The section also proposes open questions on Ritz values from approximate invariant subspaces and notes that independent MR3 computations for left and right singular vectors can create residual-error difficulties.Modified MR3 variants for bidiagonal SVD have been proposed, but reported implementations encountered inaccuracies or were abandoned.

4 Low-rank approximation and index selection

This section surveys open questions on low-rank approximation and index selection, focusing on when greedy and randomized methods can achieve quasi-optimal guarantees. It also discusses kernel approximation, rank-structured matrices, and related subset-selection questions.

  • Column subset selection: The column subset selection problem seeks k columns whose projection gives a good low-rank approximation, but finding the optimum requires exhaustive search.Efficient methods instead seek quasi-optimal choices with slowly growing approximation factors.
  • Rank-revealing methods: Rank-revealing methods must balance approximation quality with computational cost, interpretability, and preservation of matrix structure.The truncated SVD achieves the optimal approximation factor, but practical algorithms may be preferred for other reasons.
  • Greedy and randomized selection: Worst-case bounds for CPQR and complete-pivoting LU can grow exponentially with k, motivating structural assumptions under which greedy or randomized methods become quasi-optimal.The stated open problem asks for efficient algorithms satisfying CSSP or RRQR guarantees under such assumptions.
  • Kernel approximation: Empirical GECP performance is closer to the sharper low-rank rate than existing theoretical bounds, especially when Λ may be as large as 10^5 or 10^6.This highlights a substantial gap between observed behavior and current guarantees.
  • Kernel approximation: For the fermionic kernel, GECP achieves k = O(log(Λ) log(1/ε)) while guaranteeing ∥K − bK∥∞ ≤ ε, matching the rank scaling known for existence of a low-rank approximation.The same scaling holds for finite matrices sampled at distinct time and frequency nodes.
  • Structured matrix classes: For positive-definite SDD matrices, submodularity has a negative answer via a 3×3 parametrized counterexample, whereas the corresponding SDDM case has a positive answer.The broader behavior for related matrix classes remains an open direction.

5 Sketching

This section develops open questions about when weaker sketching guarantees suffice for accurate numerical linear algebra, especially least squares, low-rank approximation, and sparse embeddings. It also records several 2026 updates that answer some questions negatively or establish conjectured parameter scalings.

  • 5 Sketching: Sketching reduces large linear algebra problems through random linear maps, but basic theoretical and practical questions about effective sketching matrices remain.The section focuses on identifying sufficient properties of sketching matrices for particular algorithms.
  • 5.1 Injections vs. embeddings: The meta-problem asks for the weakest assumption on a sketching matrix that guarantees a desired accuracy level for a given algorithm.Candidate assumptions include symmetric or asymmetric OSE guarantees and OSI guarantees, while candidate algorithms include sketch-and-solve, sketch-and-precondition, and randomized SVD.
  • 5.1 Injections vs. embeddings: An August 2026 update reports a negative answer for sketch-and-solve relative error: an OSI with injectivity 1−ε need not provide a 1+O(ε) guarantee controlled solely by its failure parameter.The counterexample uses a 2 × 1 matrix A, and the weakest convenient additional condition remains open.
  • 5.1 Injections vs. embeddings: A separate August 2026 update reports a 2 × 2 counterexample showing that an OSI does not necessarily guarantee relative-error randomized SVD.For least squares with p=1, an OSI with failure probability δ=O(ε) is known to suffice, motivating questions about whether smaller failure probabilities yield broader guarantees.
  • 5.1 Injections vs. embeddings: For ℓp regression, a p-isotropic ℓp-OSI property is reported to yield a constant-factor sketch-and-solve guarantee for 1 ≤ p < ∞.The corresponding open problem asks for an OSI generalization that ensures constant-factor accuracy.
  • 5.2 Sparse dimensionality reduction maps: Sparse embedding questions examine whether very sparse maps can satisfy OSI or OSE guarantees with near-linear sketch dimensions.Updates report full two-sided OSE bounds for SparseStack at ζ=O((log r)/ε) and k=O(r/ε^2), while another result rules out ζ=4 and k=2r for constant injectivity in a broader setting.
  • 5.2 Sparse dimensionality reduction maps: A traditional SRHT has constant-distortion OSE dimension O((r + log n) log r), but some r-dimensional subspaces are not injected when k=o(r log r).The rerandomized SRHT with k=O(r/ε^2) remains posed as an open OSE question.

6 Tensor, quantum, and other problems

This section poses open problems in tensor approximation, nonlinear eigenvalue problems arising from Hamiltonians, and matrix-function approximation. The questions seek stronger approximation guarantees, efficient convergent algorithms, and optimal computation under restricted matrix operations.

  • 6.1 Optimal bounds for Tensor Train decomposition: Tree tensor decompositions include tensor trains and Tucker decompositions, but their best general error bounds arise from sequential low-rank matrix approximations.These algorithms perform m−1 SVDs following the tree structure.
  • 6.1 Optimal bounds for Tensor Train decomposition: For Tucker decomposition, the sequential HOSVD error bound is tight for tensors of certain order and size, and no polynomial-time general improvement is known for tree decompositions.The resulting open problem asks for a polynomial-time approximation algorithm with a strictly better guarantee.
  • 6.2 Eigenvalue systems arising from Hamiltonians: The Hamiltonian problem minimizes a tensor-product state energy and reformulates the objective as a nonlinear eigenvalue problem of dimension nd.The construction satisfies z^Hz=f(x1,...,xd).
  • 6.2 Eigenvalue systems arising from Hamiltonians: The open algorithmic goal is an efficient solver for this nonlinear eigenvalue problem with superlinear local convergence and iteration complexity linear in d.The formulation is motivated by numerical solvers for quantum-chemical models and related self-consistent-field procedures.
  • 6.2 Eigenvalue systems arising from Hamiltonians: A general quadratically convergent method for these nonlinear eigenvalue problems requires an nd×nd Hessian, while the complexity and convergence of methods for this family remain limited in study.This motivates seeking algorithms that avoid that large Hessian requirement.
  • 6.3 Matrix functions: One open question asks whether the ℓ2-optimal polynomial approximation to sign satisfies equioscillation, which has been observed empirically but is established here only for ℓ∞ approximation.The restricted composition class consists of T polynomials pt(x)=atx+btx^3, with each factor computable using two matrix-matrix multiplications.
  • 6.3 Matrix functions: Matrix-function questions study polynomial approximations to sign(A) using matrix multiplications, additions, and scalings rather than eigendecomposition.The approximation targets the spectral set [−1,−δ] ∪ [δ,1].
Loading 2602.05394v3…