Search

preprint at https://arxiv.org/abs/1811.05500 ispartof: ELECTRONIC TRANSACTIONS ON NUMERICAL ANALYSIS vol:51 pages:240-261 status: Published online

In recent years, several new approaches for solving the large-scale continuous-time algebraic Riccati equation have appeared in the literature. Amodei and Buchot suggest computing a low-dimensional invariant subspace of the associated Hamiltonian matrix. Simoncini and Lin also target the Hamiltonian matrix, but in a different way: they iterate on the Cayley-transformed matrix with various shifts. Wong and Balakrishnan directly generalize the Lyapunov ADI-method to the Riccati equation. In this paper we introduce another method, inspired by the Cholesky-factored variant of the Lyapunov ADI-method. The advantage of the new algorithm is in its immediate and efficient low-rank formulation, and a simpler implementation compared to the three algorithms mentioned above. We discuss the theoretical properties of the new method, as well as various shift selection strategies. Finally, we show that all of the seemingly different methods listed above in fact produce exactly the same iterates when used with the same parameters: they are algorithmically different descriptions of the same approximation sequence to the Riccati solution.

Recovering brain connectivity from tract tracing data is an important computational problem in the neurosciences. Mesoscopic connectome reconstruction was previously formulated as a structured matrix regression problem (Harris et al. in Neural Information Processing Systems, 2016), but existing techniques do not scale to the whole-brain setting. The corresponding matrix equation is challenging to solve due to large scale, ill-conditioning, and a general form that lacks a convergent splitting. We propose a greedy low-rank algorithm for the connectome reconstruction problem in very high dimensions. The algorithm approximates the solution by a sequence of rank-one updates which exploit the sparse and positive definite problem structure. This algorithm was described previously (Kressner and Sirković in Numer Lin Alg Appl 22(3):564–583, 2015) but never implemented for this connectome problem, leading to a number of challenges. We have had to design judicious stopping criteria and employ efficient solvers for the three main sub-problems of the algorithm, including an efficient GPU implementation that alleviates the main bottleneck for large datasets. The performance of the method is evaluated on three examples: an artificial “toy” dataset and two whole-cortex instances using data from the Allen Mouse Brain Connectivity Atlas. We find that the method is significantly faster than previous methods and that moderate ranks offer a good approximation. This speedup allows for the estimation of increasingly large-scale connectomes across taxa as these data become available from tracing experiments. The data and code are available online. ispartof: Journal Of Mathematical Neuroscience vol:9 issue:1 ispartof: location:Germany status: Published online

Low-rank Krylov methods are one of the few options available in the literature to address the numerical solution of large-scale general linear matrix equations. These routines amount to well-known Krylov schemes that have been equipped with a couple of low-rank truncations to maintain a feasible storage demand in the overall solution procedure. However, such truncations may affect the convergence properties of the adopted Krylov method. In this paper we show how the truncation steps have to be performed in order to maintain the convergence of the Krylov routine. Several numerical experiments validate our theoretical findings.

In this work, we consider two types of large-scale quadratic matrix equations: continuous-time algebraic Riccati equations, which play a central role in optimal and robust control, and unilateral quadratic matrix equations, which arise from stochastic processes on 2D lattices and vibrating systems. We propose a simple and fast way to update the solution to such matrix equations under low-rank modifications of the coefficients. Based on this procedure, we develop a divide-and-conquer method for quadratic matrix equations with coefficients that feature a specific type of hierarchical low-rank structure, which includes banded matrices. This generalizes earlier work on linear matrix equations. Numerical experiments indicate the advantages of our newly proposed method versus iterative schemes combined with hierarchical low-rank arithmetic.

In this paper, balancing based model order reduction (MOR) for large-scale linear discrete-time time-invariant systems in prescribed finite time intervals is studied. The first main topic is the development of error bounds regarding the approximated output vector within the time limits. The influence of different components in the established bounds will be highlighted. After that, the second part of the article proposes strategies that enable an efficient numerical execution of time-limited balanced truncation for large-scale systems. Numerical experiments illustrate the performance of the proposed techniques., Comment: 23 pages, 4 figures

The numerical treatment of large-scale, nonsymmetric algebraic Riccati equations (NAREs) by a low-rank variant of Newton׳s method is considered. We discuss a method to compute approximations to the solution of the NARE in a factorized form of low rank. The occurring large-scale Sylvester equations are dealt with using the factored alternating direction implicit iteration (fADI). Several performance enhancing strategies available for the factored ADI as well as the related Newton-ADI for symmetric algebraic Riccati equations are generalized to this combination. This includes the efficient computation of the norm of the residual matrix, adapted shift parameter strategies for fADI, and an acceleration of Newton׳s scheme by means of a Galerkin projection. Numerical experiments illustrate the capabilities of the proposed method to solve high-dimensional NAREs.

In this paper, we discuss numerical methods for solving large-scale continuous-time algebraic Riccati equations. These methods have been the focus of intensive research in recent years, and significant progress has been made in both the theoretical understanding and efficient implementation of various competing algorithms. There are several goals of this manuscript. The first is to gather in one place an overview of different approaches for solving large-scale Riccati equations, and to point to the recent advances in each of them. The second goal is to analyze and compare the main computational ingredients of these algorithms and to detect their strong points and their potential bottlenecks. Finally, we want to compare the effective implementations of all methods on a set of relevant benchmark examples, giving an indication of their relative performance.

Model order reduction (MOR) is often applied to spatially-discretized partial differential equations to reduce their order and hence decrease computational complexity. A reduced system can be obtained, e.g., by time-limited balanced truncation, a method that aims to construct an accurate reduced order model on a given finite time interval [ 0 , T ] . This particular balancing related MOR technique is studied in this paper. An L T 2 -error bound based on the truncated time-limited singular values is proved and is the main result of this paper. The derived error bound converges (as T → ∞ ) to the well-known H ∞ -error bound of unrestricted balanced truncation, a scheme that is used to construct a good reduced system on the entire time line. The techniques within the proofs of this paper can also be applied to unrestricted balanced truncation so that a relatively short time domain proof of the H ∞ -error bound is found here.

The rational Krylov subspace method (RKSM) and the low-rank alternating directions implicit (LR-ADI) iteration are established numerical tools for computing low-rank solution factors of large-scale Lyapunov equations. In order to generate the basis vectors for the RKSM, or extend the low-rank factors within the LR-ADI method the repeated solution to a shifted linear system is necessary. For very large systems this solve is usually implemented using iterative methods, leading to inexact solves within this inner iteration. We derive theory for a relaxation strategy within these inexact solves, both for the RKSM and the LR-ADI method. Practical choices for relaxing the solution tolerance within the inner linear system are then provided. The theory is supported by several numerical examples., 31 pages, 4 figures, 2 tables

The convergence of GMRES for solving linear systems can be influenced heavily by the structure of the right hand side. Within the solution of eigenvalue problems via inverse iteration or subspace iteration, the right hand side is generally related to an approximate invariant subspace of the linear system. We give detailed and new bounds on (block) GMRES that take the special behavior of the right hand side into account and explain the initial sharp decrease of the GMRES residual. The bounds give rise to adapted preconditioners applied to the eigenvalue problems, e.g. tuned and polynomial preconditioners. The numerical results show that the new (block) GMRES bounds are much sharper than conventional bounds and that preconditioned subspace iteration with either a tuned or polynomial preconditioner should be used in practice., second revised version

We consider the problem of determining an optimal semi-active damping of vibrating systems. For this damping optimization we use a minimization criterion based on the impulse response energy of the system. The optimization approach yields a large number of Lyapunov equations which have to be solved. In this work, we propose an optimization approach that works with reduced systems which are generated using the parametric dominant pole algorithm. This optimization process is accelerated with a modal approach while the initial parameters for the parametric dominant pole algorithm are chosen in advance using residual bounds. Our approach calculates a satisfactory approximation of the impulse response energy while providing a significant acceleration of the optimization process. Numerical results illustrate the effectiveness of the proposed algorithm.

In this article we investigate model order reduction of large-scale systems using time-limited balanced truncation, which restricts the well known balanced truncation framework to prescribed finite time intervals. The main emphasis is on the efficient numerical realization of this model reduction approach in case of large system dimensions. We discuss numerical methods to deal with the resulting matrix exponential functions and Lyapunov equations which are solved for low-rank approximations. Our main tool for this purpose are rational Krylov subspace methods. We also discuss the eigenvalue decay and numerical rank of the solutions of the Lyapunov equations. These results, and also numerical experiments, will show that depending on the final time horizon, the numerical rank of the Lyapunov solutions in time-limited balanced truncation can be smaller compared to standard balanced truncation. In numerical experiments we test the approaches for computing low-rank factors of the involved Lyapunov solutions and illustrate that time-limited balanced truncation can generate reduced order models having a higher accuracy in the considered time region.

SUMMARY Convergence results are provided for inexact two-sided inverse and Rayleigh quotient iteration, which extend the previously established results to the generalized non-Hermitian eigenproblem and inexact solves with a decreasing solve tolerance. Moreover, the simultaneous solution of the forward and adjoint problem arising in two-sided methods is considered, and the successful tuning strategy for preconditioners is extended to two-sided methods, creating a novel way of preconditioning two-sided algorithms. Furthermore, it is shown that inexact two-sided Rayleigh quotient iteration and the inexact two-sided Jacobi-Davidson method (without subspace expansion) applied to the generalized preconditioned eigenvalue problem are equivalent when a certain number of steps of a Petrov–Galerkin–Krylov method is used and when this specific tuning strategy is applied. Copyright © 2014 John Wiley & Sons, Ltd.

We investigate the factored alternating directions implicit (ADI) iteration for large and sparse Sylvester equations. A novel low-rank expression for the associated Sylvester residual is established which enables cheap computations of the residual norm along the iteration, and which yields a reformulated factored ADI iteration. The application to generalized Sylvester equations is considered as well. We also discuss the efficient handling of complex shift parameters and reveal interconnections between the ADI iterates w.r.t. those complex shifts. This yields a further modification of the factored ADI iteration which employs only an absolutely necessary amount of complex arithmetic operations and storage, and which produces low-rank solution factors consisting of entirely real data. Certain linear matrix equations, such as, e.g., cross Gramian Sylvester and Stein equations, are in fact special cases of generalized Sylvester equations and we show how specially tailored low-rank ADI iterations can be deduced from the generalized factored ADI iteration.

We consider balanced truncation model order reduction for symmetric second-order systems. The occurring large-scale generalized and structured Lyapunov equations are solved with a specially adapted low-rank alternating directions implicit (ADI) type method. Stopping criteria for this iteration are investigated, and a new result concerning the Lyapunov residual within the low-rank ADI method is established. We also propose a goal-oriented stopping criterion which tries to incorporate the balanced truncation approach already during the ADI iteration. The model reduction approach using the ADI method with different stopping criteria is evaluated on several test systems.

System analysis and optimization of combustion engines and engine components are increasingly supported by digital simulations. In the simulation process of combustion engines multi physics simulations are used. As an example, in the simulation of a crank drive the mechanical subsystem is coupled to a hydrodynamic subsystem. As far as the modeling of the mechanical subsystems is concerned, elastic multibody systems are frequently used. During the simulation many equations must be solved simultaneously, the hydrodynamic equations as well as the equations of motion of each body in the elastic multibody system. Since the discretization of the elastic bodies, e.g with the help of the finite element method, introduces a large number of elastic degrees of freedom, an efficient simulation of the system becomes difficult. The linear model reduction of the elastic degrees of freedom is a key step for using flexible bodies in multibody systems and turning simulations more efficient from a computational point of view. In recent years, a variety of new reduction methods alongside the traditional techniques were developed in applied mathematics. Some of these methods are reviewed and compared for reducing the equations of motion of an elastic body used in multibody systems. The special focus of this work is on balanced truncation model order reduction, which is a singular value based reduction technique using the Gramian matrices of the system. We investigate a version of this method that is adapted to the structure of a special class of second order dynamical systems which is important for the particular application discussed here. The simulation of a crank drive with a flexible crankshaft is taken as technically relevant example. The results are compared to other methods like Krylov approaches or modal reduction.

The solution of large-scale Lyapunov equations is a crucial problem for several fields of modern applied mathematics. The low-rank Cholesky factor version of the alternating directions implicit method (LRCF-ADI) is one iterative algorithm that computes approximate low-rank factors of the solution. In order to achieve fast convergence it requires adequate shift parameters, which can be complex if the matrices defining the Lyapunov equation are unsymmetric. This will require complex arithmetic computations as well as storage of complex data and thus, increase the overall complexity and memory requirements of the method. In this article we propose a novel reformulation of LRCF-ADI which generates real low-rank factors by carefully exploiting the dependencies of the iterates with respect to pairs of complex conjugate shift parameters. It significantly reduces the amount of complex arithmetic calculations and requirements for complex storage. It is hence often superior in terms of efficiency compared to other real formulations.

In this contribution we propose a novel dual Lyapunov solver that computes low rank factors of the reachability and observability Gramians of a control system simultaneously. This is especially helpful in balanced truncation model order reduction applications, where the singular values of the product of these Gramian factors are the basis of the truncation process. For stopping the dual iteration scheme we therefore propose a tailored stopping criterion aiming at the accurate computation of these singular values.

We investigate second order balanced truncation model order reduction for symmetric linear time invariant second order control systems. This special structure decreases the required computational effort significantly. Moreover, we show how stability of the original model can be preserved for such systems. We briefly discuss the numerical solution of the occurring large-scale Lyapunov equations with a modified low-rank ADI method. The approach is tested on finite element models of mechanical structures.

Based on our recent findings [1, 2] regarding the low-rank ADI iteration for large-scale Lyapunov equations, we propose a mathematically equivalent formulation of the LR-ADI whose iteration is directly connected to low rank factors of the Lyapunov residual. If complex shift parameters occur and are handled efficiently [1], this reformulation is slightly more efficient from a computational point of view, since it guarantees that the right hand sides of all linear systems that need to be solved are real. (© 2013 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)

We present a new reformulation of the low-rank ADI method for solving large-scale Lyapunov equations which uses only real arithmetic operation and storage in the presence of complex shift parameters. This makes the method applicable on computing environments where complex computations and storage are not supported or not efficiently available. For generalized Lyapunov equations it is significantly more efficient than the older completely real formulation of the low-rank ADI as confirmed by numerical examples. (© 2012 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)

Technische Universiteit Eindhoven: Centre for Analysis, Scientific Computing and Applications, Den Dolech 2, NL-5612AZ EindhovenOne crucial step of the solution of large-scale generalized eigenvalue problems with iterative subspace methods, e.g. Arnoldi,Jacobi-Davidson, is a projection of the original large-scale problem onto a low dimensional subspaces. Here we investigatetwo-sided methods, where approximate eigenvalues together with their right and left eigenvectors of the full-size problem areextracted from the resulting small eigenproblem. The two-sided Ritz-Galerkin projection can be seen as the most basic formof this approach. It usually provides a good convergence towards the extremal eigenvalues of the spectrum. For improvingthe convergence towards interior eigenvalues, we investigate two approaches based on harmonic subspace extractions for thegeneralized eigenvalue problem.