The Sherman–Morrison–Woodbury formula for generalized linear matrix equations and applications

We discuss the use of a matrix-oriented approach for numerically solving the dense matrix equation AX + XAT + M1XN1 + … + MℓXNℓ = F, with ℓ ≥ 1, and Mi, Ni, i = 1, …, ℓ of low rank. The approach relies on the Sherman–Morrison–Woodbury formula formally defined in the vectorized form of the problem, but applied in the matrix setting. This allows one to solve medium size dense problems with computational costs and memory requirements dramatically lower than with a Kronecker formulation. Application problems leading to medium size equations of this form are illustrated and the performance of the matrix-oriented method is reported. The application of the procedure as the core step in the solution of the large-scale problem is also shown. In addition, a new explicit method for linear tensor equations is proposed, that uses the discussed matrix equation procedure as a key building block

Low-rank updates and a divide-and-conquer method for linear matrix equations

Linear matrix equations, such as the Sylvester and Lyapunov equations, play an important role in various applications, including the stability analysis and dimensionality reduction of linear dynamical control systems and the solution of partial differential equations. In this work, we present and analyze a new algorithm, based on tensorized Krylov subspaces, for quickly updating the solution of such a matrix equation when its coefficients undergo low-rank changes. We demonstrate how our algorithm can be utilized to accelerate the Newton method for solving continuous-time algebraic Riccati equations. Our algorithm also forms the basis of a new divide-and-conquer approach for linear matrix equations with coefficients that feature hierarchical low-rank structure, such as HODLR, HSS, and banded matrices. Numerical experiments demonstrate the advantages of divide-and-conquer over existing approaches, in terms of computational time and memory consumption

On the numerical solution of a T-Sylvester type matrix equation arising in the control of stochastic partial differential equations

We outline a derivation of a nonlinear system of equations, which finds the entries of an m×N matrix K, given the eigenvalues of a matrix D, a diagonal N ×N matrix A and an N ×m matrix B. These matrices are related through the matrix equation D = 2A + BK + K tB t , which is sometimes called a t-Sylvester equation. The need to prescribe the eigenvalues of the matrix D is motivated by the control of the surface roughness of certain nonlinear SPDEs (e.g., the stochastic Kuramoto-Sivashinsky equation) using nontrivial controls. We implement the methodology to solve numerically the nonlinear system for various test cases, including matrices related to the control of the stochastic Kuramoto-Sivashinsky equation and for randomly generated matrices. We study the effect of increasing the dimensions of the system and changing the size of the matrices B and K (which correspond to using more or less controls) and find good convergence of the solutions

Interpolatory methods for H∞\mathcal{H}_\infty model reduction of multi-input/multi-output systems

We develop here a computationally effective approach for producing high-quality $\mathcal{H}_\infty$-approximations to large scale linear dynamical systems having multiple inputs and multiple outputs (MIMO). We extend an approach for $\mathcal{H}_\infty$ model reduction introduced by Flagg, Beattie, and Gugercin for the single-input/single-output (SISO) setting, which combined ideas originating in interpolatory $\mathcal{H}_2$-optimal model reduction with complex Chebyshev approximation. Retaining this framework, our approach to the MIMO problem has its principal computational cost dominated by (sparse) linear solves, and so it can remain an effective strategy in many large-scale settings. We are able to avoid computationally demanding $\mathcal{H}_\infty$ norm calculations that are normally required to monitor progress within each optimization cycle through the use of "data-driven" rational approximations that are built upon previously computed function samples. Numerical examples are included that illustrate our approach. We produce high fidelity reduced models having consistently better $\mathcal{H}_\infty$ performance than models produced via balanced truncation; these models often are as good as (and occasionally better than) models produced using optimal Hankel norm approximation as well. In all cases considered, the method described here produces reduced models at far lower cost than is possible with either balanced truncation or optimal Hankel norm approximation

Inverting a complex matrix

We analyze a complex matrix inversion algorithm proposed by Frobenius, which we call the Frobenius inversion. We show that the Frobenius inversion uses the least number of real matrix multiplications and inversions among all complex matrix inversion algorithms. We also analyze numerical properties of the Frobenius inversion. We prove that the Frobenius inversion runs faster than the widely used method based on LU decomposition if and only if the ratio of the running time of the real matrix inversion to that of the real matrix multiplication is greater than $5/4$. We corroborate this theoretical result by numerical experiments. Moreover, we apply the Frobenius inversion to matrix sign function, Sylvester equation, and polar decomposition. In each of these examples, the Frobenius inversion is more efficient than inversion via LU-decomposition

Approximate solutions to large nonsymmetric differential Riccati problems with applications to transport theory

In the present paper, we consider large scale nonsymmetric differential matrix Riccati equations with low rank right hand sides. These matrix equations appear in many applications such as control theory, transport theory, applied probability and others. We show how to apply Krylov-type methods such as the extended block Arnoldi algorithm to get low rank approximate solutions. The initial problem is projected onto small subspaces to get low dimensional nonsymmetric differential equations that are solved using the exponential approximation or via other integration schemes such as Backward Differentiation Formula (BDF) or Rosenbrok method. We also show how these technique could be easily used to solve some problems from the well known transport equation. Some numerical experiments are given to illustrate the application of the proposed methods to large-scale problem

Numerical solution of a class of quasi-linear matrix equations

Given the matrix equation ${\bf A X} + {\bf X B} + f({\bf X }) {\bf C} ={\bf D}$ in the unknown $n\times m$ matrix ${\bf X }$, we analyze existence and uniqueness conditions, together with computational solution strategies for $f \,: \mathbb{R}^{n \times m} \to \mathbb{R}$ being a linear or nonlinear function. We characterize different properties of the matrix equation and of its solution, depending on the considered classes of functions $f$. Our analysis mainly concerns small dimensional problems, though several considerations also apply to large scale matrix equations

Solving rank structured Sylvester and Lyapunov equations

We consider the problem of efficiently solving Sylvester and Lyapunov equations of medium and large scale, in case of rank-structured data, i.e., when the coefficient matrices and the right-hand side have low-rank off-diagonal blocks. This comprises problems with banded data, recently studied by Haber and Verhaegen in "Sparse solution of the Lyapunov equation for large-scale interconnected systems", Automatica, 2016, and by Palitta and Simoncini in "Numerical methods for large-scale Lyapunov equations with symmetric banded data", SISC, 2018, which often arise in the discretization of elliptic PDEs. We show that, under suitable assumptions, the quasiseparable structure is guaranteed to be numerically present in the solution, and explicit novel estimates of the numerical rank of the off-diagonal blocks are provided. Efficient solution schemes that rely on the technology of hierarchical matrices are described, and several numerical experiments confirm the applicability and efficiency of the approaches. We develop a MATLAB toolbox that allows easy replication of the experiments and a ready-to-use interface for the solvers. The performances of the different approaches are compared, and we show that the new methods described are efficient on several classes of relevant problems