Karl Meerbergen

Karl Meerbergen, Zhijun Wang

A popular method for solving large sparse regular eigenvalue problem is the shift-and-invert Arnoldi method. This paper aims to use the method for large sparse singular pencils. In three recent papers, {\em Hochstenbach, Mehl, and Plestenjak, 2019, 2023, and 2024}, propose regularization of the singular pencil, using randomly chosen regularization matrices. We propose sparse regularization matrices obtained from the pivoting sequence of a sparse LU factorization. As a side effect, the LU factorization often is rank revealing, which facilitates finding a regularization. Numerical examples illustrate that the LU factorization mostly detects the normal rank and finds a suitable sparse regularization. A rank correction method is proposed for the cases where the normal rank is not determined correctly. For full rank rectangular eigenvalue problems, the pivoting sequence of existing sparse direct system solvers can be used. We compare with randomized regularization methods: preservation of sparsity is beneficial for performance, and often, the accuracy of the eigenvalue solver.

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.

Pieter Lietaert, Javier Pérez, Bart Vandereycken et al.

We present a method for solving nonlinear eigenvalue problems using rational approximation. The method uses the AAA method by Nakatsukasa, Sète, and Trefethen to approximate the nonlinear eigenvalue problem by a rational eigenvalue problem and is embedded in the state space representation of a rational polynomial by Su and Bai. The advantage of the method, compared to related techniques such as NLEIGS and infinite Arnoldi, is the efficient computation by an automatic procedure. In addition, a set-valued approach is developed that allows building a low degree rational approximation of a nonlinear eigenvalue problem. The method perfectly fits the framework of the Compact rational Krylov methods (CORK and TS-CORK), allowing to efficiently solve large scale nonlinear eigenvalue problems. Numerical examples show that the presented framework is competitive with NLEIGS and usually produces smaller linearizations with the same accuracy but with less effort for the user.

Michiel E. Hochstenbach, Karl Meerbergen, Emre Mengi et al.

We propose subspace methods for 3-parameter eigenvalue problems. Such problems arise when separation of variables is applied to separable boundary value problems; a particular example is the Helmholtz equation in ellipsoidal and paraboloidal coordinates. While several subspace methods for 2-parameter eigenvalue problems exist, their extensions to three parameter setting seem to be challenging. An inherent difficulty is that, while for 2-parameter eigenvalue problems we can exploit a relation to Sylvester equations to obtain a fast Arnoldi type method, such a relation does not seem to exist when there are three or more parameters. Instead, we introduce a subspace iteration method with projections onto generalized Krylov subspaces that are constructed from scratch at every iteration using certain Ritz vectors as the initial vectors. Another possibility is a Jacobi--Davidson type method for three or more parameters, which we generalize from its 2-parameter counterpart. For both approaches, we introduce a selection criterion for deflation that is based on the angles between left and right eigenvectors. The Jacobi--Davidson approach is devised to locate eigenvalues close to a prescribed target, yet it often also performs well when eigenvalues are sought based on the proximity of one of the components to a prescribed target. The subspace iteration method is devised specifically for the latter task. The proposed approaches are suitable especially for problems where the computation of several eigenvalues is required with high accuracy. Matlab implementations of both methods have been made available in the package MultiParEig.

Joris Tavernier, Jaak Simm, Karl Meerbergen et al.

Solving linear systems is often the computational bottleneck in real-life problems. Iterative solvers are the only option due to the complexity of direct algorithms or because the system matrix is not explicitly known. Here, we develop a two-level preconditioner for regularized least squares linear systems involving a feature or data matrix. Variants of this linear system may appear in machine learning applications, such as ridge regression, logistic regression, support vector machines and Bayesian regression. We use clustering algorithms to create a coarser level that preserves the principal components of the covariance or Gram matrix. This coarser level approximates the dominant eigenvectors and is used to build a subspace preconditioner accelerating the Conjugate Gradient method. We observed speed-ups for artificial and real-life data.

Daan Camps, Karl Meerbergen, Raf Vandebril

We propose a rational QZ method for the solution of the dense, unsymmetric generalized eigenvalue problem. This generalization of the classical QZ method operates implicitly on a Hessenberg, Hessenberg pencil instead of on a Hessenberg, triangular pencil. Whereas the QZ method performs nested subspace iteration driven by a polynomial, the rational QZ method allows for nested subspace iteration driven by a rational function, this creates the additional freedom of selecting poles. In this article we study Hessenberg, Hessenberg pencils, link them to rational Krylov subspaces, propose a direct reduction method to such a pencil, and introduce the implicit rational QZ step. The link with rational Krylov subspaces allows us to prove essential uniqueness (implicit Q theorem) of the rational QZ iterates as well as convergence of the proposed method. In the proofs, we operate directly on the pencil instead of rephrasing it all in terms of a single matrix. Numerical experiments are included to illustrate competitiveness in terms of speed and accuracy with the classical approach. Two other types of experiments exemplify new possibilities. First we illustrate that good pole selection can be used to deflate the original problem during the reduction phase, and second we use the rational QZ method to implicitly filter a rational Krylov subspace in an iterative method.

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.

Luis Fernando Cusicanqui Lopez, Ramadan Krasniqi, Florian Feppon et al.

Density-based topology optimization has become a powerful method for automatically generating optimized designs in a wide variety of applications. However, it comes with a large computational cost when solving the physical model requires large-scale simulations. Here, we investigate the use of model order reduction (MOR) techniques to accelerate the simulations in the context of thermal design applications. We project the governing and the adjoint equations onto a low-dimensional subspace by constructing two distinct reduced bases -- one for the forward state and one for the adjoint system -- using solution snapshots from previous design iterations. These snapshots are generated using either the high-fidelity solver or inaccurate fast solvers, such as the one-shot method \citep{amir2024one}. Additionally, we demonstrate that properly selecting the stopping criterion for the iterative linear solver is crucial for the effective use of reduced models. In our 3D example, the proposed framework reduces the overall total simulation time relative to the high-fidelity workflow by a factor up to $3$ when combined with high-fidelity solves and a factor up to $16$ when combined with the one-shot method. Moreover, we find that the reduced order model approach is able to achieve a speed up of $1.54$ with respect to the one-shot method.

Joris Tavernier, Jaak Simm, Adam Arany et al.

Bayesian regression remains a simple but effective tool based on Bayesian inference techniques. For large-scale applications, with complicated posterior distributions, Markov Chain Monte Carlo methods are applied. To improve the well-known computational burden of Markov Chain Monte Carlo approach for Bayesian regression, we developed a multilevel Gibbs sampler for Bayesian regression of linear mixed models. The level hierarchy of data matrices is created by clustering the features and/or samples of data matrices. Additionally, the use of correlated samples is investigated for variance reduction to improve the convergence of the Markov Chain. Testing on a diverse set of data sets, speed-up is achieved for almost all of them without significant loss in predictive performance.

Koen Ruymbeek, Karl Meerbergen, Wim Michiels

In applications of linear algebra including nuclear physics and structural dynamics, there is a need to deal with uncertainty in the matrices. We focus on matrices that depend on a set of parameters $ω$ and we are interested in the minimal eigenvalue of a matrix pencil $(A,B)$ with $A,B$ symmetric and $B$ positive definite. If $ω$ can be interpreted as the realisation of random variables, one may be interested in statistical moments of the minimal eigenvalue. In order to obtain statistical moments, we need a fast evaluation of the eigenvalue as a function of $ω$. Since this is costly for large matrices,we are looking for a small parametrized eigenvalue problem whose minimal eigenvalue makes a small error with the minimal eigenvalue of the large eigenvalue problem. The advantage, in comparison with a global polynomial approximation (on which, e.g., the polynomial chaos approximation relies), is that we do not suffer from the possible non-smoothness of the minimal eigenvalue. The small scale eigenvalue problem is obtained by projection of the large scale problem. Our main contribution is that for constructing the subspace we use multiple eigenvectors as well as derivatives of eigenvectors.We provide theoretical results and document numerical experiments regarding the beneficial effect of adding multiple eigenvectors and derivatives.

Joris Tavernier, Jaak Simm, Karl Meerbergen et al.

High-dimensional data requires scalable algorithms. We propose and analyze three scalable and related algorithms for semi-supervised discriminant analysis (SDA). These methods are based on Krylov subspace methods which exploit the data sparsity and the shift-invariance of Krylov subspaces. In addition, the problem definition was improved by adding centralization to the semi-supervised setting. The proposed methods are evaluated on a industry-scale data set from a pharmaceutical company to predict compound activity on target proteins. The results show that SDA achieves good predictive performance and our methods only require a few seconds, significantly improving computation time on previous state of the art.

Karl Meerbergen, Javier Pérez

We revisit the numerical stability of the two-level orthogonal Arnoldi (TOAR) method for computing an orthonormal basis of a second--order Krylov subspace associated with two given matrices. We show that the computed basis is close (on certain subspace metric sense) to a basis for a second-order Krylov subspace associated with nearby coefficient matrices, provided that the norms of the given matrices are not too large or too small. Thus, the results in this work provide for the first time conditions that guarantee the numerical stability of the TOAR method in computing orthonormal bases of second-order Krylov subspaces. We also study scaling the quadratic problem for improving the numerical stability of the TOAR procedure when the norms of the matrices are too large or too small. We show that in many cases the TOAR procedure applied to scaled matrices is numerically stable when the scaling introduced by Fan, Lin and Van Dooren is used.

Fatih Kangal, Karl Meerbergen, Emre Mengi et al.

We consider the minimization or maximization of the $J$th largest eigenvalue of an analytic and Hermitian matrix-valued function, and build on Mengi et al. (2014, SIAM J. Matrix Anal. Appl., 35, 699-724). This work addresses the setting when the matrix-valued function involved is very large. We describe subspace procedures that convert the original problem into a small-scale one by means of orthogonal projections and restrictions to certain subspaces, and that gradually expand these subspaces based on the optimal solutions of small-scale problems. Global convergence and superlinear rate-of-convergence results with respect to the dimensions of the subspaces are presented in the infinite dimensional setting, where the matrix-valued function is replaced by a compact operator depending on parameters. In practice, it suffices to solve eigenvalue optimization problems involving matrices with sizes on the scale of tens, instead of the original problem involving matrices with sizes on the scale of thousands.