Arnoldi Algorithm Research Articles

Structure-preserving model reduction based on multiorder Arnoldi for discrete-time second-order systems

This paper investigates the new structure-preserving model reduction methods based on multiorder Arnoldi for discrete-time second-order systems by using discrete Laguerre polynomials. The expansion coefficients of discrete-time second-order system and its equivalent first-order system satisfy the implicit recursive expressions in space spanned by discrete Laguerre polynomials, and then multiorder Arnoldi algorithm is applied to the implicit recursive expressions to produce the column-orthogonal matrices, where the orthogonal projection matrix for the equivalent first-order system is constructed by partitioning the resulting column-orthogonal matrix appropriately. The reduced-order systems can preserve the special structures of the original systems as well as match a desired number of expansion coefficients of the original outputs. The error estimate of output variables is also given for discrete-time second-order system. Two benchmark examples are simulated to illustrate the effectiveness of the proposed methods.

Mechanism of vortex oscillation around a hemisphere–cylinder body

Previous studies have shown that low-frequency vortex oscillations occur around a hemisphere–cylinder body at different angles of attack, but the underlying mechanism is still unclear. In this study, we examine the origin of the vortex oscillation using numerical simulations and global linear stability analysis. The vortex oscillation is reproduced using numerical simulations, and the oscillatory modes are computed through dynamic mode decomposition (DMD). We obtain the base flow through a selective frequency damping method, which exhibits a pair of steady leeward vortices over the body. The four unstable modes are computed using a modified Arnoldi iteration. The antisymmetric mode with a Strouhal number of 0.105 is discovered to be responsible for the alternate oscillation of the vortex pair, and the mode with a Strouhal number of 0.220 corresponds to the in-phase vortex oscillation. Their frequencies have good agreement with the modes of DMD. The other two unstable modes with higher frequencies, one antisymmetric and one symmetric, are harmonic frequencies of the above two modes. The study conclusively verifies that the vortex oscillation over a hemisphere–cylinder body originates from a global flow instability.

3D modelling of ground magnetic-source airborne TEM data for anisotropic media based on the rational Krylov method

With the rapid development of the transient electromagnetic method (TEM), the numerical simulations of three-dimensional anisotropic media have garnered attention. This study used the rational Krylov method to realise the three-dimensional modelling of the ground magnetic-source airborne transient electromagnetic data for arbitrary anisotropic media. First, based on Maxwell's equations, a time-domain electromagnetic equation was derived for anisotropic media. The finite difference method was employed to discretise the equation, and the electromagnetic field was then expressed as the product of an exponential matrix and a vector related to the source. Subsequently, the rational Arnoldi algorithm was adopted to construct the rational Krylov subspace orthogonal basis vectors, and the electromagnetic fields at different times were obtained based on a model reduction strategy. The accuracy of the proposed algorithm was verified through the comparisons of its results with those of other software packages. We calculated the response for different anisotropic parameters of the typical three-dimensional anisotropic models, and analysed the influence of the anisotropic media on the ground magnetic-source airborne transient electromagnetic data. The proposed three-dimensional forward TEM algorithm provides important technical support for identifying the anisotropic geological bodies and lays the foundation for further inversion research.

Broadband topology optimization of three-dimensional structural-acoustic interaction with reduced order isogeometric FEM/BEM

This paper presents a model order reduction method to accelerate broadband topology optimization of structural-acoustic interaction systems by coupling Finite Element Methods and Boundary Element Methods. The finite element method is used for simulating thin-shell vibration and the boundary element method for exterior acoustic fields. Moreover, the finite element and boundary element methods are implemented in the context of isogeometric analysis, whereby the geometric accuracy and high order continuity of Kirchhoff-Love shells can be guaranteed and meantime no meshing is necessary. The topology optimization method takes continuous material interpolation functions in the density and bulk modulus, and adopts adjoint variable methods for sensitivity analysis. The reduced order model is constructed based on second-order Arnoldi algorithm combined with Taylor's expansions which eliminate the frequency dependence of the system matrices. Numerical results show that the proposed algorithm can significantly improve the efficiency of broadband topology optimization analysis.

Fast Subsurface EM Simulation Based on the Mixed Finite-Element Method and Second-Order Arnoldi Algorithm

Fast Subsurface EM Simulation Based on the Mixed Finite-Element Method and Second-Order Arnoldi Algorithm

Padé via Arnoldi with single‐size matrix simplification for electromagnetic fast frequency sweep

AbstractThis paper proposes a new Padé via Arnoldi algorithm with single‐size matrix simplification for electromagnetic (EM) fast frequency sweep. New equations are derived to reduce the double‐size system matrix to single‐size system matrix. We also propose a systematic algorithm to calculate S‐parameters using the simplified single‐size system matrix. Using the proposed algorithm, the EM responses can be obtained with the same accuracy while consuming much less time compared with that using the existing double‐size matrix Padé via Lanczos. The proposed algorithm is demonstrated by two microwave examples.

On the Forsythe conjecture

Forsythe formulated a conjecture about the asymptotic behavior of the restarted conjugate gradient method in 1968. We translate several of his results into modern terms, and pose an analogous version of the conjecture (originally formulated only for symmetric positive definite matrices) for symmetric and nonsymmetric matrices. Our version of the conjecture uses a two-sided or cross iteration with the given matrix and its transpose, which is based on the projection process used in the Arnoldi (or for symmetric matrices the Lanczos) algorithm. We prove several new results about the limiting behavior of this iteration, but the conjecture still remains largely open. We hope that our paper motivates further research that eventually leads to a proof of the conjecture.

Rapid Prediction of Compressor Rotating Stall Inception Using Arnoldi Eigenvalue Algorithm

This paper presents a stability model that can make a rapid prediction of the rotating stall inception in turbomachinery and provide the spatial distribution of the corresponding instability mode. In addition, this model can take the three-dimensional geometry of blades and complex flow details in the compressor into consideration, and the solution of the development process of small perturbations can be converted to a nonlinear eigenvalue problem. We propose a solution method by converting the nonlinear eigenvalue problem into a generalized one; then, it can be solved by the Arnoldi algorithm. The proposed method can shorten the elapsed time from hundreds of hours to a few minutes, as compared with the methods adopted in previous works, substantially reducing the computational cost. Furthermore, the spatial distribution of eigenvectors can be obtained to investigate the characteristics of the perturbation mode, which can be applied as a foundation to set the inlet/outlet boundary conditions and select the eigenvalue representing the rotating stall inception. In the cases of a transonic isolated rotor and a subsonic one-stage compressor, the results are in accordance with those measured in experiments, verifying the accuracy and effectiveness of the stability model. Therefore, the model can be applied to evaluate the flow stability in the design stage of compressors with low computational cost.

Speed control strategy for power line inspection robot servo system considering time-varying parameters

Most servo systems in power line inspection robots consist of a motor, an independent joint, and a load. In the process of crossing obstacles, the parameters in the servo systems have conspicuous time-varying properties due to the posture changes. The time-varying properties of dynamic parameters and the flexibility of the load would cause the rotation speed of the inspection robot to fluctuate, thereby affecting the motion accuracy. In this paper, the pole placement strategy was proposed to optimize the parameters in the proportional integral (PI) controller. The optimal controller parameters were selected in different postures to ensure steady speed output in the inspection robot servo system. First, the dynamic equations of the inspection robot servo system were established. Both joint flexibility and load flexibility were considered in the modeling process. Then, the Arnoldi algorithm was used to reduce the order of the servo system, and the transfer function from the speed to the drive torque was obtained. Next, the controller parameters were optimized using the pole placement method. By reasonably selecting the pole damping coefficient, the inspection robot could obtain a stable speed output. Finally, the numerical analysis and speed control of the inspection robot in different postures were analyzed. The results showed that the control strategy of pole placement could achieve a stable rotation speed for the inspection robot.

Complex moment-based eigensolver coupled with two Krylov subspaces

Complex moment-based eigensolvers have been well studied for solving interior eigenvalue problems because of their high parallel efficiency. Recently, as a time-efficient complex moment-based eigensolver, the block SS–CAA method was proposed that is based on the communication-avoiding Arnoldi procedure regarding the high-order complex moments. These complex moment-based eigensolvers, including the block SS–CAA method, usually use with a subspace iteration technique for improving the accuracy. In this paper, we aim to improve the convergence behavior of the block SS–CAA method using a block Arnoldi iteration instead of the subspace iteration. The proposed method is based on a special subspace coupled with two Krylov subspaces: one is for the high-order complex moments, and the other is for the iteration technique. Numerical experiments indicate that the proposed method has a higher convergence rate than the block SS–CAA method and the FEAST eigensolver.

A Structure Preserving Shift-Invert Infinite Arnoldi Algorithm for a Class of Delay Eigenvalue Problems with Hamiltonian Symmetry

This work considers a class of delay eigenvalue problems that admit a spectrum similar to that of a Hamiltonian matrix, in the sense that the spectrum is symmetric with respect to both the real and imaginary axes. We propose a method to iteratively approximate the eigenvalues closest to a given purely real or imaginary shift, while preserving the symmetries of the spectrum. To this end, our method exploits the equivalence between the considered delay eigenvalue problem and the eigenvalue problem associated with a linear but infinite-dimensional operator. To compute the eigenvalues closest to the given shift, we apply a specifically chosen shift-invert transformation to this linear operator and compute the eigenvalues with the largest modulus of the new shifted and inverted operator using an (infinite) Arnoldi procedure. The advantage of the chosen shift-invert transformation is that the spectrum of the transformed operator has a “real skew-Hamiltonian”-like structure. Furthermore, it is proven that the Krylov subspace constructed by applying this operator satisfies an orthogonality property in terms of a specifically chosen bilinear form. By taking this property into account during the orthogonalization process, it is ensured that, even in the presence of rounding errors, the obtained approximation for, e.g., a simple, purely imaginary eigenvalue is simple and purely imaginary. The presented work can thus be seen as an extension of [V. Mehrmann and D. Watkins, SIAM J. Sci. Comput., 22 (2001), pp. 1905–1925] to the considered class of delay eigenvalue problems. Although the presented method is initially defined on function spaces, it can be implemented using finite-dimensional linear algebra operations. The performance of the resulting numerical algorithm is verified for two example problems: the first example illustrates the advantage of the proposed approach in preserving purely imaginary eigenvalues when working in finite precision, while the second one demonstrates its applicability to a large scale problem.

Operator growth in open quantum systems: lessons from the dissipative SYK

We study the operator growth in open quantum systems with dephasing dissipation terms, extending the Krylov complexity formalism of [1]. Our results are based on the study of the dissipative q-body Sachdev-Ye-Kitaev (SYKq) model, governed by the Markovian dynamics. We introduce a notion of “operator size concentration” which allows a diagrammatic and combinatorial proof of the asymptotic linear behavior of the two sets of Lanczos coefficients (an and bn) in the large q limit. Our results corroborate with the semi-analytics in finite q in the large N limit, and the numerical Arnoldi iteration in finite q and finite N limit. As a result, Krylov complexity exhibits exponential growth following a saturation at a time that grows logarithmically with the inverse dissipation strength. The growth of complexity is suppressed compared to the closed system results, yet it upper bounds the growth of the normalized out-of-time-ordered correlator (OTOC). We provide a plausible explanation of the results from the dual gravitational side.

A Spurious-Free Domain Decomposition Method for 3-D Maxwell’s Eigenvalue Problems

A spurious-free domain decomposition method (DDM) based on a nonconforming tetrahedron–hexahedron mesh is proposed to solve 3-D Maxwell’s eigenvalue problems with several boundary conditions. In this DDM, the nonzero spurious modes and zero (dc) spurious modes are eliminated by locally modifying the tetrahedral mesh at the subdomain interface and implementing a tree–cotree technique. This method is a generalization of the Nitsche’s method, in which the mutually conjugate transposed terms present are omitted to reduce the computational costs. Moreover, we derive the constrained equations for the quartic Maxwell’s eigenvalue problem and impose them on Ritz vectors in the Arnoldi algorithm to remove the dc spurious modes. Numerical experiments demonstrate that the DDM is about 2.6 times more efficient than the commercial software COMSOL.

Operator growth and Krylov construction in dissipative open quantum systems

Inspired by the universal operator growth hypothesis, we extend the formalism of Krylov construction in dissipative open quantum systems connected to a Markovian bath. Our construction is based upon the modification of the Liouvillian superoperator by the appropriate Lindbladian, thereby following the vectorized Lanczos algorithm and the Arnoldi iteration. This is well justified due to the incorporation of non-Hermitian effects due to the environment. We study the growth of Lanczos coefficients in the transverse field Ising model (integrable and chaotic limits) for boundary amplitude damping and bulk dephasing. Although the direct implementation of the Lanczos algorithm fails to give physically meaningful results, the Arnoldi iteration retains the generic nature of the integrability and chaos as well as the signature of non-Hermiticity through separate sets of coefficients (Arnoldi coefficients) even after including the dissipative environment. Our results suggest that the Arnoldi iteration is meaningful and more appropriate in dealing with open systems.

Optimal computational parameters for maximum accuracy and minimum cost of Arnoldi-based time-stepping methods for flow global stability analysis

Global instability analysis of flows is often performed via time-stepping methods, based on the Arnoldi algorithm. When setting up these methods, several computational parameters must be chosen, which affect intrinsic errors of the procedure, such as the truncation errors, the discretization error of the flow solver, the error associated with the nonlinear terms of the Navier-Stokes equations and the error associated with the limited size of the approximation of the Jacobian matrix. This paper develops theoretical equations for the estimation of optimal balance between accuracy and cost for each case. The 2D open cavity flow is used both for explaining the effect of the parameters on the accuracy and the cost of the solution, and for verifying the quality of the predictions. The equations demonstrate the impact of each parameter on the quality of the solution. For example, if higher order methods are used for approaching a Fr\'echet derivative in the procedure, it is shown that the solution deteriorates more rapidly for larger grids or less accurate flow solvers. On the other hand, lower order approximations are more sensitive to the initial disturbance magnitude. Nevertheless, for accurate flow solvers and moderate grid dimensions, first order Fr\'echet derivative approximation with optimal computational parameters can provide 5 decimal place accurate eigenvalues. It is further shown that optimal parameters based on accuracy tend to also lead to the most cost-effective solution. The predictive equations, guidelines and conclusions are general, and, in principle, applicable to any flow, including 3D ones.

Parallel model order reduction based on block discrete Fourier transform and Krylov subspace for parametric systems

This paper explores a time-domain parallel parametric model order reduction (PMOR) method for parametric systems based on the block discrete Fourier transform (DFT) and Krylov subspace. The proposed method is suitable for parametric systems with non-affine parametric dependence. With Taylor expansion, the expansion coefficients of the state variable are first obtained. Then, we show that the subspace spanned by the expansion coefficients belongs to a Krylov subspace. To speed up the PMOR process, a parallel strategy based on the block DFT and the structured matrices is proposed to compute the matrices involved in the Krylov subspace. This can avoid directly computing the inverse of the large-scale matrix. After that, the reduced parametric systems are constructed with the projection matrix obtained by the Arnoldi algorithm and orthonormalisation. Furthermore, we analyse the invertibility and the error estimations to guarantee the feasibility of the proposed PMOR method. Finally, the numerical experiments are given to demonstrate the effectiveness of the proposed method.

On Some Numerical Methods for Solving Large Differential Nonsymmetric Stein Matrix Equations

In this paper, we propose a new numerical method based on the extended block Arnoldi algorithm for solving large-scale differential nonsymmetric Stein matrix equations with low-rank right-hand sides. This algorithm is based on projecting the initial problem on the extended block Krylov subspace to obtain a low-dimensional differential Stein matrix equation. The obtained reduced-order problem is solved by the backward differentiation formula (BDF) method or the Rosenbrock (Ros) method, the obtained solution is used to build the low-rank approximate solution of the original problem. We give some theoretical results and report some numerical experiments.

Scaling Up Influence Functions

We address efficient calculation of influence functions for tracking predictions back to the training data. We propose and analyze a new approach to speeding up the inverse Hessian calculation based on Arnoldi iteration. With this improvement, we achieve, to the best of our knowledge, the first successful implementation of influence functions that scales to full-size (language and vision) Transformer models with several hundreds of millions of parameters. We evaluate our approach in image classification and sequence-to-sequence tasks with tens to a hundred of millions of training examples. Our code is available at https://github.com/google-research/jax-influence.

Randomized Gram--Schmidt Process with Application to GMRES

A randomized Gram--Schmidt algorithm is developed for orthonormalization of high-dimensional vectors or QR factorization. The proposed process can be less computationally expensive than the classical Gram--Schmidt process while being at least as numerically stable as the modified Gram--Schmidt process. Our approach is based on random sketching, which is a dimension reduction technique consisting in estimation of inner products of high-dimensional vectors by inner products of their small efficiently computable random images, so-called sketches. In this way, an approximate orthogonality of the full vectors can be obtained by orthogonalization of their sketches. The proposed Gram--Schmidt algorithm can provide computational cost reduction in any architecture. The benefit of random sketching can be amplified by performing the nondominant operations in higher precision. In this case the numerical stability can be guaranteed with a working unit roundoff independent of the dimension of the problem. The proposed Gram--Schmidt process can be applied to Arnoldi iteration and results in new Krylov subspace methods for solving high-dimensional systems of equations or eigenvalue problems. Among them we chose the randomized GMRES method as a practical application of the methodology.

Low-synch Gram–Schmidt with delayed reorthogonalization for Krylov solvers

The parallel strong-scaling of iterative methods is often determined by the number of global reductions at each iteration. Low-synch Gram–Schmidt algorithms are applied here to the Arnoldi algorithm to reduce the number of global reductions and therefore to improve the parallel strong-scaling of iterative solvers for nonsymmetric matrices such as the GMRES and the Krylov–Schur iterative methods. In the Arnoldi context, the QR factorization is “left-looking” and processes one column at a time. Among the methods for generating an orthogonal basis for the Arnoldi algorithm, the classical Gram–Schmidt algorithm, with reorthogonalization (CGS2) requires three global reductions per iteration. A new variant of CGS2 that requires only one reduction per iteration is presented and applied to the Arnoldi algorithm. Delayed CGS2 (DCGS2) employs the minimum number of global reductions per iteration (one) for a one-column at-a-time algorithm. The main idea behind the new algorithm is to group global reductions by rearranging the order of operations. DCGS2 must be carefully integrated into an Arnoldi expansion or a GMRES solver. Numerical stability experiments assess robustness for Krylov–Schur eigenvalue computations. Performance experiments on the ORNL Summit supercomputer then establish the superiority of DCGS2 over CGS2.