Mark Embree

Noah Amsel, Yves Baumann, Paul Beckman et al. · berkeley

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.

Sean Reiter, Mark Embree, Serkan Gugercin et al.

This work investigates the reduction of phasor measurement unit (PMU) data through low-rank matrix approximations. To reconstruct a PMU data matrix from fewer measurements, we propose the framework of interpolatory matrix decompositions (IDs). In contrast to methods relying on principal component analysis or singular value decomposition, IDs recover the complete data matrix using only a few of its rows (PMU datastreams) and/or a few of its columns (snapshots in time). This row-/column-based compression enables real-time monitoring of power transmission systems using measurements from a smaller subset of pilot datastreams, thereby minimizing communication bandwidth. The ID perspective gives a rigorous error bound on the quality of the data compression. We propose selecting the pilot measurements used in an ID via the discrete empirical interpolation method (DEIM), a greedy algorithm that aims to control the error bound. This bound yields a computable estimate of the reconstruction error during online operations. A violation of this estimate suggests a change in the system's operating conditions and thus serves as a tool for fault detection. Following a disturbance, DEIM can be used to localize the event source across all buses with high accuracy. Numerical tests on synthetic PMU data demonstrate DEIM's excellent performance in data compression and validate the proposed DEIM-based fault-detection and localization method.

David Damanik, Mark Embree, Anton Gorodetski

We survey results that have been obtained for self-adjoint operators, and especially Schrödinger operators, associated with mathematical models of quasicrystals. After presenting general results that hold in arbitrary dimensions, we focus our attention on the one-dimensional case, and in particular on several key examples. The most prominent of these is the Fibonacci Hamiltonian, for which much is known by now and to which an entire section is devoted here. Other examples that are discussed in detail are given by the more general class of Schrödinger operators with Sturmian potentials. We put some emphasis on the methods that have been introduced quite recently in the study of these operators, many of them coming from hyperbolic dynamics. We conclude with a multitude of numerical calculations that illustrate the validity of the known rigorous results and suggest conjectures for further exploration.

Jonathan Baker, Mark Embree, John Sabino

Lyapunov equations with low-rank right-hand sides often have solutions whose singular values decay rapidly, enabling iterative methods that produce low-rank approximate solutions. All previously known bounds on this decay involve quantities that depend quadratically on the departure of the coefficient matrix from normality: these bounds suggest that the larger the departure from normality, the slower the singular values will decay. We show this is only true up to a threshold, beyond which a larger departure from normality can actually correspond to faster decay of singular values: if the singular values decay slowly, the numerical range cannot extend far into the right-half plane.

Mark Embree, Ronald B. Morgan, Huy V. Nguyen

The convergence of the restarted GMRES method can be significantly improved, for some problems, by using a weighted inner product that changes at each restart. How does this weighting affect convergence, and when is it useful? We show that weighted inner products can help in two distinct ways: when the coefficient matrix has localized eigenvectors, weighting can allow restarted GMRES to focus on eigenvalues that otherwise slow convergence; for general problems, weighting can break the cyclic convergence pattern into which restarted GMRES often settles. The eigenvectors of matrices derived from differential equations are often not localized, thus limiting the impact of weighting. For such problems, incorporating the discrete cosine transform into the inner product can significantly improve GMRES convergence, giving a method we call W-GMRES-DCT. Integrating weighting with eigenvalue deflation via GMRES-DR also can give effective solutions.

Mark Embree, Jennifer A. Loe, Ronald B. Morgan

Polynomial preconditioning can improve the convergence of the Arnoldi method for computing eigenvalues. Such preconditioning significantly reduces the cost of orthogonalization; for difficult problems, it can also reduce the number of matrix-vector products. Parallel computations can particularly benefit from the reduction of communication-intensive operations. The GMRES algorithm provides a simple and effective way of generating the preconditioning polynomial. For some problems high degree polynomials are especially effective, but they can lead to stability problems that must be mitigated. A two-level "double polynomial preconditioning" strategy provides an effective way to generate high-degree preconditioners.

Mark Embree

Projection methods provide an appealing way to construct reduced-order models of large-scale linear dynamical systems: they are intuitively motivated and fairly easy to compute. Unfortunately, the resulting reduced models need not inherit the stability of the original system. How many unstable modes can these reduced models have? This note investigates this question, using theory originally motivated by iterative methods for linear algebraic systems and eigenvalue problems, and illustrating the theory with a number of small examples. From these results follow rigorous upper bounds on the number of unstable modes in reduced models generated via orthogonal projection, for both continuous- and discrete-time systems. Can anything be learned from the unstable modes in reduced-order models? Several examples illustrate how such instability can helpfully signal transient growth in the original system.

Mark Embree, Blake Keeler

To understand the solution of a linear, time-invariant differential-algebraic equation, one must analyze a matrix pencil (A,E) with singular E. Even when this pencil is stable (all its finite eigenvalues fall in the left-half plane), the solution can exhibit transient growth before its inevitable decay. When the equation results from the linearization of a nonlinear system, this transient growth gives a mechanism that can promote nonlinear instability. One might hope to enrich the conventional large-scale eigenvalue calculation used for linear stability analysis to signal the potential for such transient growth. Toward this end, we introduce a new definition of the pseudospectrum of a matrix pencil, use it to bound transient growth, explain how to incorporate a physically-relevant norm, and derive approximate pseudospectra using the invariant subspace computed in conventional linear stability analysis. We apply these tools to several canonical test problems in fluid mechanics, an important source of differential-algebraic equations.

Charles Puelz, Mark Embree, Jake Fillman

Quasiperiodic Jacobi operators arise as mathematical models of quasicrystals and in more general studies of structures exhibiting aperiodic order. The spectra of these self-adjoint operators can be quite exotic, such as Cantor sets, and their fine properties yield insight into associated dynamical systems. Quasiperiodic operators can be approximated by periodic ones, the spectra of which can be computed via two finite dimensional eigenvalue problems. Since long periods are necessary to get detailed approximations, both computational efficiency and numerical accuracy become a concern. We describe a simple method for numerically computing the spectrum of a period-$K$ Jacobi operator in $O(K^2)$ operations, and use it to investigate the spectra of Schrödinger operators with Fibonacci, period doubling, and Thue-Morse potentials.