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.

Elias Jarlebring, Giampaolo Mele, Davide Palitta et al.

We consider generalizations of the Sylvester matrix equation, consisting of the sum of a Sylvester operator and a linear operator $Π$ with a particular structure. More precisely, the commutator of the matrix coefficients of the operator $Π$ and the Sylvester operator coefficients are assumed to be matrices with low rank. We show (under certain additional conditions) low-rank approximability of this problem, i.e., the solution to this matrix equation can be approximated with a low-rank matrix. Projection methods have successfully been used to solve other matrix equations with low-rank approximability. We propose a new projection method for this class of matrix equations. The choice of subspace is a crucial ingredient for any projection method for matrix equations. Our method is based on an adaption and extension of the extended Krylov subspace method for Sylvester equations. A constructive choice of the starting vector/block is derived from the low-rank commutators. We illustrate the effectiveness of our method by solving large-scale matrix equations arising from applications in control theory and the discretization of PDEs. The advantages of our approach in comparison to other methods are also illustrated.

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.

Emil Ringh, Giampaolo Mele, Johan Karlsson et al.

We consider a nonlinear eigenvalue problem (NEP) arising from absorbing boundary conditions in the study of a partial differential equation (PDE) describing a waveguide. We propose a new computational approach for this large scale NEP based on residual inverse iteration (Resinv) with preconditioned iterative solves. Similar to many preconditioned iterative methods for discretized PDEs, this approach requires the construction of an accurate and efficient preconditioner. For the waveguide eigenvalue problem, the associated linear system can be formulated as a generalized Sylvester equation. The equation is approximated by a low-rank correction of a Sylvester equation, which we use as a preconditioner. The action of the preconditioner is efficiently computed using the matrix equation version of the Sherman-Morrison-Woodbury (SMW) formula. We show how the preconditioner can be integrated into Resinv. The results are illustrated by applying the method to large-scale problems.

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.

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.

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.

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.

Felix Nannesson Meli, Johan Tell, Shirwan Piroti et al.

Mobile networks consist of interconnected radio nodes strategically positioned across various geographical regions to provide connectivity services. The set of relations between these radio nodes, referred to as the \emph{mobile network topology}, is vital in the construction of the networking infrastructure. Typically, the connections between radio nodes and their associated cells are defined by software features that establish mobility relations (referred to as \emph{edges} in this paper) within the mobile network graph through heuristic methods. Although these approaches are efficient, they encounter significant limitations, particularly since edges can only be established prior to the installation of physical hardware. In this work, we use graph-based deep learning methods to determine mobility relations (edges), trained on radio node configuration data and reliable mobility relations set by Automatic Neighbor Relations (ANR) in stable networks. This paper focuses on measuring the accuracy and precision of different graph-based deep learning approaches applied to real-world mobile networks. We evaluated two deep learning models. Our comprehensive experiments on Telecom datasets obtained from operational Telecom Networks demonstrate the effectiveness of the graph neural network (GNN) model and multilayer perceptron. Our evaluation showed that considering graph structure improves results, which motivates the use of GNNs. Additionally, we investigated the use of heuristics to reduce the training time based on the distance between radio nodes to eliminate irrelevant cases. Our investigation showed that the use of these heuristics improved precision and accuracy considerably.

Elias Jarlebring, Federico Poloni

The delay Lyapunov equation is an important matrix boundary-value problem which arises as an analogue of the Lyapunov equation in the study of time-delay systems $\dot{x}(t) = A_0x(t)+A_1x(t-τ)+B_0u(t)$. We propose a new algorithm for the solution of the delay Lyapunov equation. Our method is based on the fact that the delay Lyapunov equation can be expressed as a linear system of equations, whose unknown is the value $U(τ/2)\in\mathbb{R}^{n\times n}$, i.e., the delay Lyapunov matrix at time $τ/2$. This linear matrix equation with $n^2$ unknowns is solved by adapting a preconditioned iterative method such as GMRES. The action of the $n^2\times n^2$ matrix associated to this linear system can be computed by solving a coupled matrix initial-value problem. A preconditioner for the iterative method is proposed based on solving a T-Sylvester equation $MX+X^TN=C$, for which there are methods available in the literature. We prove that the preconditioner is effective under certain assumptions. The efficiency of the approach is illustrated by applying it to a time-delay system stemmingfrom the discretization of a partial differential equation with delay. Approximate solutions to this problem can be obtained for problems of size up to $n\approx 1000$, i.e., a linear system with $n^2\approx 10^6$ unknowns, a dimension which is outside of the capabilities of the other existing methods for the delay Lyapunov equation.

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.