Leonardo Robol

Leonardo Robol, Raf Vandebril, Paul Van Dooren

We present a framework for the construction of linearizations for scalar and matrix polynomials based on dual bases which, in the case of orthogonal polynomials, can be described by the associated recurrence relations. The framework provides an extension of the classical linearization theory for polynomials expressed in non-monomial bases and allows to represent polynomials expressed in product families, that is as a linear combination of elements of the form $ϕ_i(λ) ψ_j(λ)$, where $\{ ϕ_i(λ) \}$ and $\{ ψ_j(λ) \}$ can either be polynomial bases or polynomial families which satisfy some mild assumptions. We show that this general construction can be used for many different purposes. Among them, we show how to linearize sums of polynomials and rational functions expressed in different bases. As an example, this allows to look for intersections of functions interpolated on different nodes without converting them to the same basis. We then provide some constructions for structured linearizations for $\star$-even and $\star$-palindromic matrix polynomials. The extensions of these constructions to $\star$-odd and $\star$-antipalindromic of odd degree is discussed and follows immediately from the previous results.

Daniel Kressner, Stefano Massei, Leonardo Robol

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.

Dario A. Bini, Stefano Massei, Leonardo Robol

A Quasi Toeplitz (QT) matrix is a semi-infinite matrix of the kind $A=T(a)+E$ where $T(a)=(a_{j-i})_{i,j\in\mathbb Z^+}$, $E=(e_{i,j})_{i,j\in\mathbb Z^+}$ is compact and the norms $\lVert a\rVert_{\mathcal W} = \sum_{i\in\mathbb Z}|a_i|$ and $\lVert E \rVert_2$ are finite. These properties allow to approximate any QT-matrix, within any given precision, by means of a finite number of parameters. QT-matrices, equipped with the norm $\lVert A \rVert_{\mathcal QT}=α\lVert a\rVert_{\mathcal{W}} \lVert E \rVert_2$, for $α= (1+\sqrt 5)/2$, are a Banach algebra with the standard arithmetic operations. We provide an algorithmic description of these operations on the finite parametrization of QT-matrices, and we develop a MATLAB toolbox implementing them in a transparent way. The toolbox is then extended to perform arithmetic operations on matrices of finite size that have a Toeplitz plus low-rank structure. This enables the development of algorithms for Toeplitz and quasi-Toeplitz matrices whose cost does not necessarily increase with the dimension of the problem. Some examples of applications to computing matrix functions and to solving matrix equations are presented, and confirm the effectiveness of the approach.

Maria Girardi, Cristina Padovani, Daniele Pellegrini et al.

A novel method for performing model updating on finite element models is presented. The approach is particularly tailored to modal analyses of buildings, by which the lowest frequencies, obtained by using sensors and system identification approaches, need to be matched to the numerical ones predicted by the model. This is done by optimizing some unknown material parameters (such as mass density and Young's modulus) of the materials and/or the boundary conditions, which are often known only approximately. In particular, this is the case when considering historical buildings. The straightforward application of a general-purpose optimizer can be impractical, given the large size of the model involved. In the paper, we show that, by slightly modifying the projection scheme used to compute the eigenvalues at the lowest end of the spectrum one can obtain local parametric reduced order models that, embedded in a trust-region scheme, form the basis for a reliable and efficient specialized algorithm. We describe an optimization strategy based on this approach, and we provide numerical experiments that confirm its effectiveness and accuracy.

Dario A. Bini, Leonardo Robol

We present a novel algorithm to perform the Hessenberg reduction of an $n\times n$ matrix $A$ of the form $A = D + UV^*$ where $D$ is diagonal with real entries and $U$ and $V$ are $n\times k$ matrices with $k\le n$. The algorithm has a cost of $O(n^2k)$ arithmetic operations and is based on the quasiseparable matrix technology. Applications are shown to solving polynomial eigenvalue problems and some numerical experiments are reported in order to analyze the stability of the approach

Stefano Massei, Davide Palitta, Leonardo Robol

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.

Luca Gemignani, Leonardo Robol

We develop two fast algorithms for Hessenberg reduction of a structured matrix $A = D + UV^H$ where $D$ is a real or unitary $n \times n$ diagonal matrix and $U, V \in\mathbb{C}^{n \times k}$. The proposed algorithm for the real case exploits a two--stage approach by first reducing the matrix to a generalized Hessenberg form and then completing the reduction by annihilation of the unwanted sub-diagonals. It is shown that the novel method requires $O(n^2k)$ arithmetic operations and it is significantly faster than other reduction algorithms for rank structured matrices. The method is then extended to the unitary plus low rank case by using a block analogue of the CMV form of unitary matrices. It is shown that a block Lanczos-type procedure for the block tridiagonalization of $\Re(D)$ induces a structured reduction on $A$ in a block staircase CMV--type shape. Then, we present a numerically stable method for performing this reduction using unitary transformations and we show how to generalize the sub-diagonal elimination to this shape, while still being able to provide a condensed representation for the reduced matrix. In this way the complexity still remains linear in $k$ and, moreover, the resulting algorithm can be adapted to deal efficiently with block companion matrices.

Jared L. Aurentz, Thomas Mach, Leonardo Robol et al.

This work is a continuation of "Fast and backward stable computation of roots of polynomials" by J.L. Aurentz, T. Mach, R. Vandebril, and D.S. Watkins, SIAM Journal on Matrix Analysis and Applications, 36(3): 942--973, 2015. In that paper we introduced a companion QR algorithm that finds the roots of a polynomial by computing the eigenvalues of the companion matrix in $O(n^{2})$ time using $O(n)$ memory. We proved that the method is backward stable. Here we introduce, as an alternative, a companion QZ algorithm that solves a generalized eigenvalue problem for a companion pencil. More importantly, we provide an improved backward error analysis that takes advantage of the special structure of the problem. The improvement is also due, in part, to an improvement in the accuracy (in both theory and practice) of the turnover operation, which is the key component of our algorithms. We prove that for the companion QR algorithm, the backward error on the polynomial coefficients varies linearly with the norm of the polynomial's vector of coefficients. Thus the companion QR algorithm has a smaller backward error than the unstructured QR algorithm (used by MATLAB's \texttt{roots} command, for example), for which the backward error on the polynomial coefficients grows quadratically with the norm of the coefficient vector. The companion QZ algorithm has the same favorable backward error as companion QR, provided that the polynomial coefficients are properly scaled.

Stefano Massei, Leonardo Robol

Given matrices $A$ and $B$ such that $B=f(A)$, where $f(z)$ is a holomorphic function, we analyze the relation between the singular values of the off-diagonal submatrices of $A$ and $B$. We provide family of bounds which depend on the interplay between the spectrum of the argument $A$ and the singularities of the function. In particular, these bounds guarantee the numerical preservation of quasiseparable structures under mild hypotheses. We extend the Dunford-Cauchy integral formula to the case in which some poles are contained inside the contour of integration. We use this tool together with the technology of hierarchical matrices ($\mathcal H$-matrices) for the effective computation of matrix functions with quasiseparable arguments.

Dario A. Bini, Stefano Massei, Leonardo Robol

We provide effective algorithms for solving block tridiagonal block Toeplitz systems with $m\times m$ quasiseparable blocks, as well as quadratic matrix equations with $m\times m$ quasiseparable coefficients, based on cyclic reduction and on the technology of rank-structured matrices. The algorithms rely on the exponential decay of the singular values of the off-diagonal submatrices generated by cyclic reduction. We provide a formal proof of this decay in the Markovian framework. The results of the numerical experiments that we report confirm a significant speed up over the general algorithms, already starting with the moderately small size $m\approx 10^2$.

Fernando De Terán, Bruno Iannazzo, Federico Poloni et al.

We provide necessary and sufficient conditions for the generalized $\star$-Sylvester matrix equation, $AXB + CX^\star D = E$, to have exactly one solution for any right-hand side E. These conditions are given for arbitrary coefficient matrices $A, B, C, D$ (either square or rectangular) and generalize existing results for the same equation with square coefficients. We also review the known results regarding the existence and uniqueness of solution for generalized Sylvester and $\star$-Sylvester equations.

Alberto Bucci, Silvia Gazzola, Leonardo Robol

LSQR and LSMR are iterative methods, based on the Golub-Kahan bidiagonalization algorithm, widely used for large-scale linear least squares problems. FLSQR and FLSMR are flexible variants of LSQR and LSMR, respectively, based on a flexible Golub-Kahan (Arnoldi-like) factorization algorithm, which naturally allow modifications of the solution approximation subspace and/or handling inexact matrix-vector multiplications with the (transpose of the) coefficient matrix, thereby enabling to enforce prior information into the computed solution. The goal of this paper is to introduce sFLSQR and sFLSMR, i.e., sketched variants of FLSQR and FLSMR, respectively, where randomization becomes particularly effective, as it allows to recover short recurrences for the solution approximation. In particular, this paper explores applications to large-scale inverse problems, showing the ability of the new randomized solvers to alleviate computational bottlenecks while preserving reconstruction quality. A theoretical analysis of sFLSQR and sFLSMR is provided, and their performance is validated through numerical experiments.

Stefano Massei, Leonardo Robol

This paper investigates explicit expressions for the error associated with the block rational Krylov approximation of matrix functions. Two formulas are proposed, both derived from characterizations of the block FOM residual. The first formula employs a block generalization of the residual polynomial, while the second leverages the block collinearity of the residuals. A posteriori error bounds based on the knowledge of spectral information of the argument are derived and tested on a set of examples. Notably, both error formulas and their corresponding upper bounds do not require the use of quadratures for their practical evaluation.

Jared Aurentz, Thomas Mach, Leonardo Robol et al.

In the last decade matrix polynomials have been investigated with the primary focus on adequate linearizations and good scaling techniques for computing their eigenvalues and eigenvectors. In this article we propose a new method for computing a factored Schur form of the associated companion pencil. The algorithm has a quadratic cost in the degree of the polynomial and a cubic one in the size of the coefficient matrices. Also the eigenvectors can be computed at the same cost. The algorithm is a variant of Francis's implicitly shifted QR algorithm applied on the companion pencil. A preprocessing unitary equivalence is executed on the matrix polynomial to simultaneously bring the leading matrix coefficient and the constant matrix term to triangular form before forming the companion pencil. The resulting structure allows us to stably factor each matrix of the pencil as a product of $k$ matrices of unitary-plus-rank-one form, admitting cheap and numerically reliable storage. The problem is then solved as a product core chasing eigenvalue problem. A backward error analysis is included, implying normwise backward stability after a proper scaling. Computing the eigenvectors via reordering the Schur form is discussed as well. Numerical experiments illustrate stability and efficiency of the proposed methods.

Dario Andrea Bini, Stefano Massei, Leonardo Robol

It was recently observed that the singular values of the off-diagonal blocks of the matrix sequences generated by the Cyclic Reduction algorithm decay exponentially. This property was used to solve, with a higher efficiency, certain quadratic matrix equations encountered in the analysis of queueing models. In this paper, we provide a sharp theoretical bound to the basis of this exponential decay together with a tool for its estimation based on a rational interpolation problem. Applications to solving $n\times n$ block tridiagonal block Toeplitz systems with $n\times n$ semiseparable blocks and certain generalized Sylvester equations in $O(n^2\log n)$ arithmetic operations are shown.

Dario A. Bini, Leonardo Robol

A new class of linearizations and $\ell$-ifications for $m\times m$ matrix polynomials $P(x)$ of degree $n$ is proposed. The $\ell$-ifications in this class have the form $A(x) = D(x) + (e\otimes I_m) W(x)$ where $D$ is a block diagonal matrix polynomial with blocks $B_i(x)$ of size $m$, $W$ is an $m\times qm$ matrix polynomial and $e=(1,\ldots,1)^t\in\mathbb C^q$, for a suitable integer $q$. The blocks $B_i(x)$ can be chosen a priori, subjected to some restrictions. Under additional assumptions on the blocks $B_i(x)$ the matrix polynomial $A(x)$ is a strong $\ell$-ification, i.e., the reversed polynomial of $A(x)$ defined by $A^\#(x) := x^{\mathrm{deg} A(x)} A(x^{-1})$ is an $\ell$-ification of $P^\#(x)$. The eigenvectors of the matrix polynomials $P(x)$ and $A(x)$ are related by means of explicit formulas. Some practical examples of $\ell$-ifications are provided. A strategy for choosing $B_i(x)$ in such a way that $A(x)$ is a well conditioned linearization of $P(x)$ is proposed. Some numerical experiments that validate the theoretical results are reported