Arnoldi Method Research Articles

A preconditioned block Arnoldi method for large scale Lyapunov and algebraic Riccati equations

In the present paper, we propose a preconditioned Newton---Block Arnoldi method for solving large continuous time algebraic Riccati equations. Such equations appear in control theory, model reduction, circuit simulation amongst other problems. At each step of the Newton process, we solve a large Lyapunov matrix equation with a low rank right hand side. These equations are solved by using the block Arnoldi process associated with a preconditioner based on the alternating direction implicit iteration method. We give some theoretical results and report numerical tests to show the effectiveness of the proposed approach.

An inexact shift‐and‐invert Arnoldi algorithm for Toeplitz matrix exponential

SummaryWe revisit the shift‐and‐invert Arnoldi method proposed in [S. Lee, H. Pang, and H. Sun.Shift‐invert Arnoldi approximation to the Toeplitz matrix exponential, SIAM J. Sci. Comput., 32: 774–792, 2010] for numerical approximation to the product of Toeplitz matrix exponential with a vector. In this approach, one has to solve two large‐scale Toeplitz linear systems in advance. However, if the desired accuracy is high, the cost will be prohibitive. Therefore, it is interesting to investigate how to solve the Toeplitz systems inexactly in this method. The contribution of this paper is in three regards. First, we give a new stability analysis on the Gohberg–Semencul formula (GSF) and define the GSF condition number of a Toeplitz matrix. It is shown that when the size of the Toeplitz matrix is large, our result is sharper than the one given in [M. Gutknecht and M. Hochbruck.The stability of inversion formulas for Toeplitz matrices, Linear Algebra Appl., 223/224: 307–324, 1995]. Second, we establish a relation between the error of Toeplitz systems and the residual of Toeplitz matrix exponential. We show that if the GSF condition number of the Toeplitz matrix is medium‐sized, then the Toeplitz systems can be solved in a low accuracy. Third, based on this relationship, we present a practical stopping criterion for relaxing the accuracy of the Toeplitz systems and propose an inexact shift‐and‐invert Arnoldi algorithm for the Toeplitz matrix exponential problem. Numerical experiments illustrate the numerical behavior of the new algorithm and show the effectiveness of our theoretical results. Copyright © 2015 John Wiley & Sons, Ltd.

An Inexact Shift-and-Invert Arnoldi Algorithm for Large Non-Hermitian Generalised Toeplitz Eigenproblems

AbstractThe shift-and-invert Arnoldi method is a most effective approach to compute a few eigenpairs of a large non-Hermitian Toeplitz matrix pencil, where the Gohberg-Semencul formula can be used to obtain the Toeplitz inverse. However, two large non-Hermitian Toeplitz systems must be solved in the first step of this method, and the cost becomes prohibitive if the desired accuracy for this step is high — especially for some ill-conditioned problems. To overcome this difficulty, we establish a relationship between the errors in solving these systems and the residual of the Toeplitz eigenproblem. We consequently present a practical stopping criterion for their numerical solution, and propose an inexact shift-and-invert Arnoldi algorithm for the generalised Toeplitz eigenproblem. Numerical experiments illustrate our theoretical results and demonstrate the efficiency of the new algorithm.

Reduction of large-scale dynamical systems by the Krylov subspaces method: Analysis of approaches

Approaches and numerical algorithms for reducing the order of mathematical models of multidimensional dynamical systems that are based on the Krylov's subspaces method are described. To calculate matrices representing the reduced models in the state space, Lanczos' and Arnoldi's methods are used. Practical examples of the reduction of large-scale systems are presented.

IMPLICITLY RESTARTED GENERALIZED SECOND-ORDER ARNOLDI TYPE ALGORITHMS FOR THE QUADRATIC EIGENVALUE PROBLEM

We investigate the generalized second-order Arnoldi (GSOAR) method, a generalization of the SOAR method proposed by Bai and Su [ SIAM J. Matrix Anal. Appl. , 26 (2005): 640--659.], and the Refined GSOAR (RGSOAR) method for the quadratic eigenvalue problem (QEP). The two methods use the GSOAR procedure to generate an orthonormal basis of a given generalized second-order Krylov subspace, and with such basis they project the QEP onto the subspace and compute the Ritz pairs and the refined Ritz pairs, respectively. We develop implicitly restarted GSOAR and RGSOAR algorithms, in which we propose certain exact and refined shifts for respective use within the two algorithms. Numerical experiments on real-world problems illustrate the efficiency of the restarted algorithms and the superiority of the restarted RGSOAR to the restarted GSOAR. The experiments also demonstrate that both IGSOAR and IRGSOAR generally perform much better than the implicitly restarted Arnoldi method applied to the corresponding linearization problems, in terms of the accuracy and the computational efficiency.

Arnoldi methods for image deblurring with anti-reflective boundary conditions

Image deblurring with anti-reflective boundary conditions and a non-symmetric point spread function is considered. Several iterative methods based on Krylov subspace projections, as well as Arnoldi–Tikhonov regularization methods, with reblurring right or left preconditioners are compared. The aim of the preconditioner is not to accelerate the convergence, but to improve the quality of the computed solution and to increase the robustness of the regularization method. Right preconditioning in conjunction with methods based on the Arnoldi process are found to be robust and give high-quality restorations. In particular, when the observed image is contaminated by motion blur, our new method is much more effective than other approaches described in the literature, such as range restricted Arnoldi methods and the conjugate gradient method applied to the normal equations (implemented with the reblurring approach).

A key to choose subspace size in implicitly restarted Arnoldi method

The implicitly restarted Arnoldi method (IRAM) computes some eigenpairs of large sparse non Hermitian matrices. However, the size of the subspace in this method is chosen empirically. A poor choice of this size could lead to the non-convergence of the method. In this paper we propose a technique to improve the choice of the size of subspace. This approach, called multiple implicitly restarted Arnoldi method with nested subspaces (MIRAMns) is based on the projection of the problem on several nested subspaces instead of a single one. Thus, it takes advantage of several different sized subspaces. MIRAMns updates the restarting vector of an IRAM by taking the eigen-information of interest obtained in all subspaces into account. With almost the same complexity as IRAM, according to our experiments, MIRAMns improves the convergence of IRAM.

Singular Value Decompositions for Single-Curl Operators in Three-Dimensional Maxwell's Equations for Complex Media

This article focuses on solving the generalized eigenvalue problems (GEP) arising in the source-free Maxwell equation with magnetoelectric coupling effects that models three-dimensional complex media. The goal is to compute the smallest positive eigenvalues, and the main challenge is that the coefficient matrix in the discrete Maxwell equation is indefinite and degenerate. To overcome this difficulty, we derive a singular value decomposition (SVD) of the discrete single-curl operator and then explicitly express the basis of the invariant subspace corresponding to the nonzero eigenvalues of the GEP. Consequently, we reduce the GEP to a null space free standard eigenvalue problem (NFSEP) that contains only the nonzero (complex) eigenvalues of the GEP and can be solved by the shift-and-invert Arnoldi method without being disturbed by the null space. Furthermore, the basis of the eigendecomposition is chosen carefully so that we can apply fast Fourier transformation (FFT-) based matrix vector multiplication t...

Efficient solution of the simplified [formula omitted] equations

In this paper we show new solver strategies for the multigroup SPN equations for nuclear reactor analysis. By forming the complete matrix over space, moments, and energy, a robust set of solution strategies may be applied. Power iteration, shifted power iteration, Rayleigh quotient iteration, Arnoldi's method, and a generalized Davidson method, each using algebraic and physics-based multigrid preconditioners, have been compared on the C5G7 MOX test problem as well as an operational pressurized water reactor model. Our results show that the most efficient approach is the generalized Davidson method, which is 30–40 times faster than traditional power iteration and 6–10 times faster than Arnoldi's method.

A REDUCED ORDER MODEL BASED ON BLOCK ARNOLDI METHOD FOR AEROELASTIC SYSTEM

A reduced-order model (ROM) based on block Arnoldi algorithm to quickly predict flutter boundary of aeroelastic system is investigated. First, a mass–damper–spring dynamic system is tested, which shows that the low dimension system produced by the block Arnoldi method can keep a good dynamic property with the original system in low and high frequencies. Then a two-degree of freedom transonic nonlinear aerofoil aeroelastic system is used to validate the suitability of the block Arnoldi method in flutter prediction analysis. In the aerofoil case, the ROM based on a linearized model is obtained through a high-fidelity nonlinear computational fluid dynamics (CFD) calculation. The order of the reduced model is only 8 while it still has nearly the same accuracy as the full 9600-order model. Compared with the proper orthogonal decomposition (POD) method, the results show that, without snapshots the block Arnoldi/ROM has a unique superiority by maintaining the system stability aspect. The flutter boundary of the aeroelastic system predicted by the block Arnoldi/ROM agrees well with the CFD and reference results. The Arnoldi/ROM provides an efficient and convenient tool to quick analyze the system stability of nonlinear transonic aeroelastic systems.

A positivity preserving inexact Noda iteration for computing the smallest eigenpair of a large irreducible $$M$$ M -matrix

In this paper, based on the Noda iteration, we present inexact Noda iterations (INI), to find the smallest eigenvalue and the associated positive eigenvector of a large irreducible nonsingular $$M$$M-matrix. The positivity of approximations is critical in applications, and if the approximations lose the positivity then they may be meaningless and could not be interpreted. We propose two different inner tolerance strategies for solving the inner linear systems involved, and prove that the convergence of resulting INI algorithms is globally linear and superlinear with the convergence order $$\frac{1+\sqrt{5}}{2}$$1+52, respectively. The proposed INI algorithms are structure preserving and maintains the positivity of approximate eigenvectors. We also revisit the exact Noda iteration and establish a new quadratic convergence result. All the above is first done for the problem of computing the Perron root and the positive Perron vector of an irreducible nonnegative matrix and is then adapted to computing the smallest eigenpair of the irreducible nonsingular $$M$$M-matrix. Numerical examples illustrate that the proposed INI algorithms are practical, and they always preserve the positivity of approximate eigenvectors. We compare them with the Jacobi---Davidson method, the implicitly restarted Arnoldi method and the explicitly restarted Krylov---Schur method, all of which cannot guarantee the positivity of approximate eigenvectors, and illustrate that the overall efficiency of the INI algorithms is competitive with and can be considerably higher than the latter three methods.

PageRank Computation Using a Multiple Implicitly Restarted Arnoldi Method for Modeling Epidemic Spread

A parallel implementation based on implicitly restarted Arnoldi method (MIRAM) is proposed for calculating dominant eigenpair of stochastic matrices derived from very large real networks. Their high damping factor makes many existing algorithms less efficient, while MIRAM could be promising. Also, we apply this method in an epidemic application. We describe in this paper a stochastic model based on PageRank to simulate the epidemic spread, where a PageRank-like infection vector is calculated by MIRAM to help establish efficient vaccination strategy. MIRAM is implemented within the framework of Trilinos, targeting big data and sparse matrices representing scale-free networks, also known as power law networks. Hypergraph partitioning approach is employed to minimize the communication overhead. The algorithm is tested on a nation wide cluster of clusters Grid5000. Experiments on very large networks such as twitter and yahoo with over 1 billion nodes are conducted. With our parallel implementation, a speedup of $$27\times $$ is met compared to the sequential solver.

Boundary fidelity and entanglement in the symmetry protected topological phase of the SSH model

We present a detailed study of the fidelity, the entanglement entropy and the entanglement spectrum, for a dimerized chain of spinless fermions—a simplified Su–Schrieffer–Heeger (SSH) model—with open boundary conditions which is a well-known example for a model supporting a symmetry protected topological (SPT) phase. In the non-interacting case the Hamiltonian matrix is tridiagonal and the eigenvalues and vectors can be given explicitly as a function of a single parameter which is known analytically for odd chain lengths and can be determined numerically in the even length case. From a scaling analysis of these data for essentially semi-infinite chains we obtain the fidelity susceptibility and show that it contains a boundary contribution which is different in the topologically ordered than in the topologically trivial phase. For the entanglement spectrum and entropy we confirm predictions from massive field theory for a block in the middle of an infinite chain but also consider blocks containing the edge of the chain. For the latter case we show that in the SPT phase additional entanglement—as compared to the trivial phase—is present which is localized at the boundary. Finally, we extend our study to the dimerized chain with a nearest-neighbour interaction using exact diagonalization, Arnoldi and density-matrix renormalization group methods and show that a phase transition into a topologically trivial charge-density wave phase occurs.

An adaptive dynamic Implicitly Restarted Arnoldi method for the small signal stability eigen analysis of large power systems

This paper deals with small signal stability analysis for large scale power systems using Krylov subspace techniques, we develop an algorithm based on the Implicitly Restarted method coupled with a dynamic switching approach in such manner to adapt the well known algorithm to our study issue which is converging to the few sought eigenvalues with poor damping.

3D stability analysis of Rayleigh–Bénard convection of a liquid metal layer in the presence of a magnetic field—effect of wall electrical conductivity

Rayleigh–Bénard stability of a liquid metal layer of rectangular cross section is examined in the presence of a strong magnetic field that is aligned with the horizontal direction of the cross section. The latter is much longer than the vertical direction and the cross section assumes a large aspect ratio. The side walls are treated as highly conducting. Linear stability analysis is performed allowing for three-dimensional instabilities that develop along the longitudinal direction. The finite element methodology is employed for the discretization of the stability analysis formulation while accounting for the electrical conductivity of the cavity walls. The Arnoldi method provides the dominant eigenvalues and eigenvectors of the problem. In order to facilitate parallel implementation of the numerical solution at large Hartmann numbers, Ha, domain decomposition is employed along the horizontal direction of the cross section. As the Hartmann number increases a real eigenvalue emerges as the dominant unstable eigenmode, signifying the onset of thermal convection, whose major vorticity component in the core of the layer is aligned with the direction of the magnetic field. Its wavelength along the longitudinal direction of the layer is on the order of twice its height and increases as Ha increases. The critical Grashof was obtained for large Ha and it was seen to scale like Ha2 signifying the balance between buoyancy and Lorentz forces. For well conducting side walls, the nature of the emerging flow pattern is determined by the combined conductivity of Hartmann walls and Hartmann layers, cH + Ha−1. When poor conducting Hartmann walls are considered, cH ≪ 1, the critical eigensolution is characterized by well defined Hartmann and side layers. The side layers are characterized by fast fluid motion in the magnetic field direction as a result of the electromagnetic pumping in the vicinity of the Hartmann walls. Increasing the electrical conductivity of the Hartmann walls was seen to delay the onset of thermal convection, while retaining the above scaling at criticality. Furthermore, for both conducting and insulating Hartmann walls and the entire range of Ha numbers that was examined, there was no tendency for a well defined quasi two-dimensional structure to develop owing to the convective motion in the core. A connection is made between the above findings and previous experimental investigations indicating the onset of standing waves followed by travelling waves as Gr is further increased beyond its critical value.

Preface to BIT 54:2

Now we are back with a set of contributions from the regular stream of incoming BIT manuscripts. We note that several contributions deal with algorithms for time propagation, where the system is linear and the linear algebra is crucial for the success of a numerical algorithm for propagation in time, see the papers by Arnold, Bodendiek and Freitag. The last paper by Yan, deals with propagation governed by a fractional derivative. The papers byNovati andSundealwith regularization of ill posed problems, while those of Goda, Kressner and Ling treat empirical approximation problems. We have the two papers by Juntunen and Lamichhane that handle finite elements, and the one by Barreras on high relative accuracy issues in numerical linear algebra. These are the papers: Andreas Arnold andTobias Jahnke study a systemof ordinary differential equations, where the unknown is amultidimensional tensor, represented in the hierarchical Tucker format. This way, the time propagation can be described in fewer dimensions than for general tensors. It opens up for handling some practically interestingmodels in physics chemistry and economics. Alvaro Barreras and Juan Manuel Pena describe a subtraction free algorithm for LDU decomposition of an almost diagonally dominant Z-matrix. It is marginally more expensive than Gaussian elimination and expands the class of linear systems that can be solved with high relative accuracy. Andre Bodendiek and Matthias Bollhofer describe an adaptive-order rational Arnoldi method for moment matchingmodel reduction of a descriptor system, coming from time dependent Maxwell equations in computational electromagnetism. Crucial

Solving leaky modes on a dielectric slab waveguide involving materials of arbitrary dielectric anisotropy with a finite-element formulation

A planar optical waveguide mode solver is established based on a finite-element (FE) formulation for determining the guided and leaky modes that exist on waveguides made of anisotropic materials with an arbitrary permittivity tensor, for example, with arbitrary optic-axis orientation in the uniaxially anisotropic material case. Correct numerical determination of the complex effective index, especially its imaginary part that gives the modal leakage, is particularly emphasized referring to available analytical solutions. For the situation when the optic axis changes its direction only in the plane parallel to the waveguide interface planes, analytical characteristic equations for solving purely guided and leaky modes are separately derived in a more systematic manner compared with prior analytical formulae reported more than three decades ago, with the obtained complex effective indices agreeing excellently with FE solutions. It is found that in the FE analysis of leaky modes, the thickness of the perfectly matched layer (PML) and the PML theoretical reflection coefficient should be properly chosen. The FE formulation is based on either the three electric-field components or the three magnetic-field components using quadratic nodal bases, resulting in a quadratic eigenvalue equation that is then solved by the shift-and-invert Arnoldi method.

On expansion of search subspaces for large non-Hermitian eigenproblems

How to expand a given subspace efficiently is an important task for large sparse eigenvalue problems. In this paper, we are interested in the following question: which vector from a given subspace, after multiplied by the matrix A, will give a better expanded search subspace. In (2008) [33], Ye investigated this problem and determined how to choose the optimal expansion vector theoretically. Unfortunately, the result is not computable since it involves some unknown information on the desired eigenvector. When A is symmetric, it was suggested that the Ritz vector may lead to a good candidate for subspace expansion. However, when A is non-Hermitian, the Ritz vector may not be a satisfactory choice. The contribution of this paper is twofold. First, we suggest to use the refined Ritz vector to expand the search subspace, and propose a residual expansion Arnoldi method for subspace expansion. Theoretical results justify the use of the refined Ritz vector. Second, we prove that the elements of the primitive refined Ritz vector have a decreasing pattern going to zero. We then show that the decreasing pattern exists in an arbitrary primitive approximate eigenvector, and derive an inexact residual expansion Arnoldi method for subspace expansion. Numerical examples show the effectiveness of our theoretical results and illustrate the numerical behavior of the new method for subspace expansion.

Inverse eigenvalue problems for extended Hessenberg and extended tridiagonal matrices

In inverse eigenvalue problems one tries to reconstruct a matrix, satisfying some constraints, given some spectral information. Here, two inverse eigenvalue problems are solved.First, given the eigenvalues and the first components of the associated eigenvectors (called the weight vector) an extended Hessenberg matrix with prescribed poles is computed possessing these eigenvalues and satisfying the eigenvector constraints. The extended Hessenberg matrix is retrieved by executing particularly designed unitary similarity transformations on the diagonal matrix containing the eigenvalues. This inverse problem closely links to orthogonal rational functions: the extended Hessenberg matrix contains the recurrence coefficients given the nodes (eigenvalues), poles (poles of the extended Hessenberg matrix), and a weight vector (first eigenvector components) determining the discrete inner product. Moreover, it is also sort of the inverse of the (rational) Arnoldi algorithm: instead of using the (rational) Arnoldi method to compute a Krylov basis to approximate the spectrum, we will reconstruct the orthogonal Krylov basis given the spectral info.In the second inverse eigenvalue problem, we do the same, but refrain from unitarity. As a result we execute possibly non-unitary similarity transformations on the diagonal matrix of eigenvalues to retrieve a (non)-symmetric extended tridiagonal matrix. The algorithm will be less stable, but it will be faster, as the extended tridiagonal matrix admits a low cost factorization of O(n) (n equals the number of eigenvalues), whereas the extended Hessenberg matrix does not. Again there is a close link with orthogonal function theory, the extended tridiagonal matrix captures the recurrence coefficients of bi-orthogonal rational functions. Moreover, it is again sort of inverse of the nonsymmetric Lanczos algorithm: given spectral properties, we reconstruct the two basis Krylov matrices linked to the nonsymmetric Lanczos algorithm.

Eigenvalue solvers for three dimensional photonic crystals with face-centered cubic lattice

To numerically determine the band structure of three-dimensional photonic crystals with face-centered cubic lattices, we study how the associated large-scale generalized eigenvalue problem (GEP) can be solved efficiently. The main computational challenge is due to the complexity of the coefficient matrix and the fact that the desired eigenvalues are interior. For solving the GEP by the shift-and-invert Lanczos method, we propose a preconditioning for the associated linear system therein. Recently, a way to reformat the GEP to the null space free eigenvalue problem (NFEP) is proposed. For solving the NFEP, we analyze potential advantages and disadvantages of the null space free inverse Lanczos method, the shift–invert residual Arnoldi method, and the Jacobi–Davidson method from theoretical viewpoints. These four approaches are compared numerically to find out their properties. The numerical results suggest that the shift–invert residual Arnoldi method with an initialization scheme is the fastest and the most robust eigensolver for the target eigenvalue problems. Our findings promise to play an essential role in simulating photonic crystals.