Elias Jarlebring

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.

Elias Jarlebring

In this work we present a new method to compute the delays of delay differential equations (DDEs), such that the DDE has a purely imaginary eigenvalue. For delay differential equations with multiple delays, the critical curves or critical surfaces in delay space (that is, the set of delays where the DDE has a purely imaginary eigenvalue) are parameterized. We show how the method is related to other works in the field by treating the case where the delays are integer multiples of some delay value, i.e., commensurate delays. The parametrization is done by solving a {\em quadratic eigenvalue problem} which is constructed from the vectorization of a matrix equation and hence typically of large size. For commensurate delay differential equations, the corresponding equation is a polynomial eigenvalue problem. As a special case of the proposed method, we find a closed form for a parameterization of the critical surface for the scalar case. We provide several examples with visualizations where the computation is done with some exploitation of the structure of eigenvalue problems.

Elias Jarlebring, Simen Kvaal, Wim Michiels

Consider a symmetric matrix $A(v)\in\RR^{n\times n}$ depending on a vector $v\in\RR^n$ and satisfying the property $A(αv)=A(v)$ for any $α\in\RR\backslash{0}$. We will here study the problem of finding $(λ,v)\in\RR\times \RR^n\backslash\{0\}$ such that $(λ,v)$ is an eigenpair of the matrix $A(v)$ and we propose a generalization of inverse iteration for eigenvalue problems with this type of eigenvector nonlinearity. The convergence of the proposed method is studied and several convergence properties are shown to be analogous to inverse iteration for standard eigenvalue problems, including local convergence properties. The algorithm is also shown to be equivalent to a particular discretization of an associated ordinary differential equation, if the shift is chosen in a particular way. The algorithm is adapted to a variant of the Schrödinger equation known as the Gross-Pitaevskii equation. We use numerical simulations toillustrate the convergence properties, as well as the efficiency of the algorithm and the adaption.

Elias Jarlebring, Karl Meerbergen, Wim Michiels

The partial Schur factorization can be used to represent several eigenpairs of a matrix in a numerically robust way. Different adaptions of the Arnoldi method are often used to compute partial Schur factorizations. We propose here a technique to compute a partial Schur factorization of a nonlinear eigenvalue problem (NEP). The technique is inspired by the algorithm in [8], now called the infinite Arnoldi method. The infinite Arnoldi method is a method designed for NEPs, and can be interpreted as Arnoldi's method applied to a linear infinite-dimensional operator, whose reciprocal eigenvalues are the solutions to the NEP. As a first result we show that the invariant pairs of the operator are equivalent to invariant pairs of the NEP. We characterize the structure of the invariant pairs of the operator and show how one can carry out a modification of the infinite Arnoldi method by respecting the structure. This also allows us to naturally add the feature known as locking. We nest this algorithm with an outer iteration, where the infinite Arnoldi method for a particular type of structured functions is appropriately restarted. The restarting exploits the structure and is inspired by the well-known implicitly restarted Arnoldi method for standard eigenvalue problems. The final algorithm is applied to examples from a benchmark collection, showing that both processing time and memory consumption can be considerably reduced with the restarting technique.

Elias Jarlebring, Michiel E. Hochstenbach

Several recent methods used to analyze asymptotic stability of delay-differential equations (DDEs) involve determining the eigenvalues of a matrix, a matrix pencil or a matrix polynomial constructed by Kronecker products. Despite some similarities between the different types of these so-called matrix pencil methods, the general ideas used as well as the proofs differ considerably. Moreover, the available theory hardly reveals the relations between the different methods. In this work, a different derivation of various matrix pencil methods is presented using a unifying framework of a new type of eigenvalue problem: the polynomial two-parameter eigenvalue problem, of which the quadratic two-parameter eigenvalue problem is a special case. This framework makes it possible to establish relations between various seemingly different methods and provides further insight in the theory of matrix pencil methods. We also recognize a few new matrix pencil variants to determine DDE stability. Finally, the recognition of the new types of eigenvalue problem opens a door to efficient computation of DDE stability.

Parikshit Upadhyaya, Elias Jarlebring, Emanuel H. Rubensson

In this paper, we present a local convergence analysis of the self-consistent field (SCF) iteration using the density matrix as the state of a fixed-point iteration. Sufficient and almost necessary conditions for local convergence are formulated in terms of the spectral radius of the Jacobian of a fixed-point map. The relationship between convergence and certain properties of the problem is explored by deriving upper bounds expressed in terms of higher gaps. This gives more information regarding how the gaps between eigenvalues of the problem affect the convergence, and hence these bounds are more insightful on the convergence behaviour than standard convergence results. We also provide a detailed analysis to describe the difference between the bounds and the exact convergence factor for an illustrative example. Finally we present numerical examples and compare the exact value of the convergence factor with the observed behaviour of SCF, along with our new bounds and the characterization using the higher gaps. We provide heuristic convergence factor estimates in situations where the bounds fail to well capture the convergence.

Sarah W. Gaaf, Elias Jarlebring

We propose a two-sided Lanczos method for the nonlinear eigenvalue problem (NEP). This two-sided approach provides approximations to both the right and left eigenvectors of the eigenvalues of interest. The method implicitly works with matrices and vectors with infinite size, but because particular (starting) vectors are used, all computations can be carried out efficiently with finite matrices and vectors. We specifically introduce a new way to represent infinite vectors that span the subspace corresponding to the conjugate transpose operation for approximating the left eigenvectors. Furthermore, we show that also in this infinite-dimensional interpretation the short recurrences inherent to the Lanczos procedure offer an efficient algorithm regarding both the computational cost and the storage.

Juan Carlos Araujo-Cabarcas, Christian Engstrom, Elias Jarlebring

We present an efficient procedure for computing resonances and resonant modes of Helmholtz problems posed in exterior domains. The problem is formulated as a nonlinear eigenvalue problem (NEP), where the nonlinearity arises from the use of a Dirichlet-to-Neumann map, which accounts for modeling unbounded domains. We consider a variational formulation and show that the spectrum consists of isolated eigenvalues of finite multiplicity that only can accumulate at infinity. The proposed method is based on a high order finite element discretization combined with a specialization of the Tensor Infinite Arnoldi method. Using Toeplitz matrices, we show how to specialize this method to our specific structure. In particular we introduce a pole cancellation technique in order to increase the radius of convergence for computation of eigenvalues that lie close to the poles of the matrix-valued function. The solution scheme can be applied to multiple resonators with a varying refractive index that is not necessarily piecewise constant. We present two test cases to show stability, performance and numerical accuracy of the method. In particular the use of a high order finite element discretization together with TIAR results in an efficient and reliable method to compute resonances.

Giampaolo Mele, Elias Jarlebring

An efficient and robust restart strategy is important for any Krylov-based method for eigenvalue problems. The tensor infinite Arnoldi method (TIAR) is a Krylov-based method for solving nonlinear eigenvalue problems (NEPs). This method can be interpreted as an Arnoldi method applied to a linear and infinite dimensional eigenvalue problem where the Krylov basis consists of polynomials. We propose new restart techniques for TIAR and analyze efficiency and robustness. More precisely, we consider an extension of TIAR which corresponds to generating the Krylov space using not only polynomials but also structured functions that are sums of exponentials and polynomials, while maintaining a memory efficient tensor representation. We propose two restarting strategies, both derived from the specific structure of the infinite dimensional Arnoldi factorization. One restarting strategy, which we call semi-explicit TIAR restart, provides the possibility to carry out locking in a compact way. The other strategy, which we call implicit TIAR restart, is based on the Krylov-Schur restart method for linear eigenvalue problem and preserves its robustness. Both restarting strategies involve approximations of the tensor structured factorization in order to reduce complexity and required memory resources. We bound the error in the infinite dimensional Arnoldi factorization showing that the approximation does not substantially influence the robustness of the restart approach. We illustrate the approaches by applying them to solve large scale NEPs that arise from a delay differential equation and a wave propagation problem. The advantages in comparison to other restart methods are also illustrated.

Elias Jarlebring, Antti Koskela, Giampaolo Mele

In this paper we take a quasi-Newton approach to nonlinear eigenvalue problems (NEPs) of the type $M(λ)v=0$, where $M:\mathbb{C}\rightarrow\mathbb{C}^{n\times n}$ is a holomorphic function. We investigate which types of approximations of the Jacobian matrix lead to competitive algorithms, and provide convergence theory. The convergence analysis is based on theory for quasi-Newton methods and Keldysh's theorem for NEPs. We derive new algorithms and also show that several well-established methods for NEPs can be interpreted as quasi-Newton methods, and thereby provide insight to their convergence behavior. In particular, we establish quasi-Newton interpretations of Neumaier's residual inverse iteration and Ruhe's method of successive linear problems.

Elias Jarlebring

Broyden's method is a general method commonly used for nonlinear systems of equations, when very little information is available about the problem. We develop an approach based on Broyden's method for nonlinear eigenvalue problems. Our approach is designed for problems where the evaluation of a matrix vector product is computationally expensive, essentially as expensive as solving the corresponding linear system of equations. We show how the structure of the Jacobian matrix can be incorporated into the algorithm to improve convergence. The algorithm exhibits local superlinear convergence for simple eigenvalues, and we characterize the convergence. We show how deflation can be integrated and combined such that the method can be used to compute several eigenvalues. A specific problem in machine tool milling, coupled with a PDE is used to illustrate the approach. The simulations are done in the julia programming language, and are provided as publicly available module for reproducability.

Antti Koskela, Elias Jarlebring

Exponential integrators that use Krylov approximations of matrix functions have turned out to be efficient for the time-integration of certain ordinary differential equations (ODEs). This holds in particular for linear homogeneous ODEs, where the exponential integrator is equivalent to approximating the product of the matrix exponential and a vector. In this paper, we consider linear inhomogeneous ODEs, $y'(t)=Ay(t)+g(t)$, where the function $g(t)$ is assumed to satisfy certain regularity conditions. We derive an algorithm for this problem which is equivalent to approximating the product of the matrix exponential and a vector using Arnoldi's method. The construction is based on expressing the function $g(t)$ as a linear combination of given basis functions $[ϕ_i]_{i=0}^\infty$ with particular properties. The properties are such that the inhomogeneous ODE can be restated as an infinite-dimensional linear homogeneous ODE. Moreover, the linear homogeneous infinite-dimensional ODE has properties that directly allow us to extend a Krylov method for finite-dimensional linear ODEs. Although the construction is based on an infinite-dimensional operator, the algorithm can be carried out with operations involving matrices and vectors of finite size. This type of construction resembles in many ways the infinite Arnoldi method for nonlinear eigenvalue problems. We prove convergence of the algorithm under certain natural conditions, and illustrate properties of the algorithm with examples stemming from the discretization of partial differential equations.

Elias Jarlebring, Emanuel H. Rubensson

Consider a matrix function f defined for Hermitian matrices. The purpose of this paper is two-fold. We derive new results for the absolute structured condition number of the matrix function and we derive new bounds for the perturbation ||f(A+E)-f(A)|| expressed in terms of eigenvalues of A and A+E. The results are general and under certain conditions hold for an arbitrary unitarily invariant matrix norm ||\cdot||. We illustrate the use of the formulas with an application to the error analysis of calculations in electronic structure theory.

Vilhelm Peterson Lithell, Victor Janssens, Elias Jarlebring et al.

We consider eigenvalue condition numbers and backward errors for a class of symmetric nonlinear eigenvalue problems with eigenvector nonlinearities. For both of these quantities, we derive explicit and computable expressions that can be evaluated with little computational effort for a given eigenpair, assuming the matrix perturbations are measured by the spectral or Frobenius norm. We also show how symmetric perturbations can be exploited in the analysis. By means of two numerical experiments we demonstrate that problems incorporating eigenvector nonlinearities potentially need to be treated with additional care, when compared to the linear or eigenvalue-nonlinear theory.

Antti Koskela, Elias Jarlebring

The Bessel-Neumann expansion (of integer order) of a function $g:\mathbb{C}\rightarrow\mathbb{C}$ corresponds to representing $g$ as a linear combination of basis functions $ϕ_0,ϕ_1,\ldots$, i.e., $g(z)=\sum_{\ell = 0}^\infty w_\ell ϕ_\ell(s)$, where $ϕ_i(z)=J_i(z)$, $i=0,\ldots$, are the Bessel functions. In this work, we study an expansion for a more general class of basis functions. More precisely, we assume that the basis functions satisfy an infinite dimensional linear ordinary differential equation associated with a Hessenberg matrix, motivated by the fact that these basis functions occur in certain iterative methods. A procedure to compute the basis functions as well as the coefficients is proposed. Theoretical properties of the expansion are studied. We illustrate that non-standard basis functions can give faster convergence than the Bessel functions.

Emanuel H. Rubensson, Elias Jarlebring, Gustaf Lorentzon

We consider the problem of efficiently computing the matrix step function of a large dense symmetric matrix. To this end, we introduce a recursive polynomial expansion method in which a composite polynomial of high degree is built recursively from component polynomials of degree eight. The component polynomial used in each iteration is designed to achieve strong amplification of the spectral gap across the step while favorably positioning the updated gap for subsequent iterations. A key ingredient is a novel evaluation scheme for arbitrary matrix polynomials of degree exactly eight requiring only three matrix-matrix multiplications and three matrices in memory. This scheme makes available a substantially larger class of component polynomials than previously possible within a three-multiplication budget, thereby expanding the class of composite polynomials that can be generated. Together with our polynomial selection strategy, this leads to a significant and consistent reduction in the number of matrix-matrix multiplications required to compute the matrix step function compared to existing recursive expansion methods.

Ricardo Vinuesa, Paola Cinnella, Jean Rabault et al. · uw

As modern scientific instruments generate vast amounts of data and the volume of information in the scientific literature continues to grow, machine learning (ML) has become an essential tool for organising, analysing, and interpreting these complex datasets. This paper explores the transformative role of ML in accelerating breakthroughs across a range of scientific disciplines. By presenting key examples -- such as brain mapping and exoplanet detection -- we demonstrate how ML is reshaping scientific research. We also explore different scenarios where different levels of knowledge of the underlying phenomenon are available, identifying strategies to overcome limitations and unlock the full potential of ML. Despite its advances, the growing reliance on ML poses challenges for research applications and rigorous validation of discoveries. We argue that even with these challenges, ML is poised to disrupt traditional methodologies and advance the boundaries of knowledge by enabling researchers to tackle increasingly complex problems. Thus, the scientific community can move beyond the necessary traditional oversimplifications to embrace the full complexity of natural systems, ultimately paving the way for interdisciplinary breakthroughs and innovative solutions to humanity's most pressing challenges.

Antti Koskela, Elias Jarlebring, Michiel E. Hochstenbach

We propose a new numerical method to solve linear ordinary differential equations of the type $\frac{\partial u}{\partial t}(t,\varepsilon) = A(\varepsilon) \, u(t,\varepsilon)$, where $A:\mathbb{C}\rightarrow\mathbb{C}^{n\times n}$ is a matrix polynomial with large and sparse matrix coefficients. The algorithm computes an explicit parameterization of approximations of $u(t,\varepsilon)$ such that approximations for many different values of $\varepsilon$ and $t$ can be obtained with a very small additional computational effort. The derivation of the algorithm is based on a reformulation of the parameterization as a linear parameter-free ordinary differential equation and on approximating the product of the matrix exponential and a vector with a Krylov method. The Krylov approximation is generated with Arnoldi's method and the structure of the coefficient matrix turns out to have an independence on the truncation parameter so that it can also be interpreted as Arnoldi's method applied to an infinite dimensional matrix. We prove the superlinear convergence of the algorithm and provide a posteriori error estimates to be used as termination criteria. The behavior of the algorithm is illustrated with examples stemming from spatial discretizations of partial differential equations.

Elias Jarlebring, Giampaolo Mele, Olof Runborg

We present a new computational approach for a class of large-scale nonlinear eigenvalue problems (NEPs) that are nonlinear in the eigenvalue. The contribution of this paper is two-fold. We derive a new iterative algorithm for NEPs, the tensor infinite Arnoldi method (TIAR), which is applicable to a general class of NEPs, and we show how to specialize the algorithm to a specific NEP: the waveguide eigenvalue problem. The waveguide eigenvalue problem arises from a finite-element discretization of a partial differential equation (PDE) used in the study waves propagating in a periodic medium. The algorithm is successfully applied to accurately solve benchmark problems as well as complicated waveguides. We study the complexity of the specialized algorithm with respect to the number of iterations m and the size of the problem n, both from a theoretical perspective and in practice. For the waveguide eigenvalue problem, we establish that the computationally dominating part of the algorithm has complexity O(nm^2 + sqrt(n) m^3). Hence, the asymptotic complexity of TIAR applied to the waveguide eigenvalue problem, for n that goes to infinity, is the same as for Arnoldi's method for standard eigenvalue problems.