Dario A. Bini

Dario A. Bini, Massimiliano Fasi, Bruno Iannazzo

We consider the numerical evaluation of the quantity $Af(A^{-1}B)$, where $A$ is Hermitian positive definite, $B$ is Hermitian, and $f$ is a function defined on the spectrum of $A^{-1}B$. This problem is related to the Hermitian-definite matrix pencil $B-λA$. We study the conditioning of the problem, and we introduce several algorithms that combine the Schur decomposition with either the matrix square root or the Cholesky factorization. We study the numerical behavior of these algorithms in floating-point arithmetic, assess their computational costs, and compare their numerical performance. Our analysis suggests that the algorithms based on the Cholesky factorization will be more accurate and efficient than those based on the matrix square root. This is confirmed by our numerical experiments.

Dario A. Bini, Vanni Noferini, Meisam Sharify

Some known results for locating the roots of polynomials are extended to the case of matrix polynomials. In particular, a theorem by A.E. Pellet [Bulletin des Sciences Mathématiques, (2), vol 5 (1881), pp.393-395], some results of D.A. Bini [Numer. Algorithms 13:179-200, 1996] based on the Newton polygon technique, and recent results of M. Akian, S. Gaubert and M. Sharify (see in particular [LNCIS, 389, Springer p.p.291-303] and [M. Sharify, Ph.D. thesis, École Polytechnique, ParisTech, 2011]). These extensions are applied for determining effective initial approximations for the numerical computation of the eigenvalues of matrix polynomials by means of simultaneous iterations, like the Ehrlich-Aberth method. Numerical experiments that show the computational advantage of these results are presented.

Dario A. Bini, V. Noferini

Given the $n\times n$ matrix polynomial $P(x)=\sum_{i=0}^kP_i x^i$, we consider the associated polynomial eigenvalue problem. This problem, viewed in terms of computing the roots of the scalar polynomial $\det P(x)$, is treated in polynomial form rather than in matrix form by means of the Ehrlich-Aberth iteration. The main computational issues are discussed, namely, the choice of the starting approximations needed to start the Ehrlich-Aberth iteration, the computation of the Newton correction, the halting criterion, and the treatment of eigenvalues at infinity. We arrive at an effective implementation which provides more accurate approximations to the eigenvalues with respect to the methods based on the QZ algorithm. The case of polynomials having special structures, like palindromic, Hamiltonian, symplectic, etc., where the eigenvalues have special symmetries in the complex plane, is considered. A general way to adapt the Ehrlich-Aberth iteration to structured matrix polynomial is introduced. Numerical experiments which confirm the effectiveness of this approach are reported.

Dario A. Bini, Stefano Massei, Beatrice Meini

Denote by $\mathcal{W}_1$ the set of complex valued functions of the form $a(z)=\sum_{i=-\infty}^{+\infty}a_iz^i$ which are continuous on the unit circle, and such that $\sum_{i=-\infty}^{+\infty}|ia_i|<\infty$. We call CQT matrix a quasi-Toeplitz matrix $A$, associated with a continuous symbol $a(z)\in\mathcal W_1$, of the form $A=T(a)+E$, where $T(a)=(t_{i,j})_{i,j\in\mathbb{Z}^+}$ is the semi-infinite Toeplitz matrix such that $t_{i,j}=a_{j-i}$, for $i,j\in\mathbb Z^+$, and $E=(e_{i,j})_{i,j\in\mathbb{Z}^+}$ is a semi-infinite matrix such that $\sum_{i,j=1}^{+\infty}|e_{i,j}|$ is finite. We prove that the class of CQT matrices is a Banach algebra with a suitable sub-multiplicative matrix norm $\|\cdot\|$. We introduce a finite representation of CQT matrices together with algorithms which implement elementary matrix operations. An application to solving quadratic matrix equations of the kind $AX^2+BX+C=0$, encountered in the solution of Quasi-Birth and Death (QBD) stochastic processes with a denumerable set of phases, is presented where $A,B,C$ are CQT matrices.

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.

Dario A. Bini, Stefano Massei, Beatrice Meini

Let $a(z)=\sum_{i\in\mathbb Z}a_iz^i$ be a complex valued continuous function, defined for $|z|=1$, such that $\sum_{i=-\infty}^{+\infty}|ia_i|<\infty$. Consider the semi-infinite Toeplitz matrix $T(a)=(t_{i,j})_{i,j\in\mathbb Z^+}$ associated with the symbol $a(z)$ such that $t_{i,j}=a_{j-i}$. A quasi-Toeplitz matrix associated with the continuous symbol $a(z)$ is a matrix of the form $A=T(a)+E$ where $E=(e_{i,j})$, $\sum_{i,j\in\mathbb Z^+}|e_{i,j}|<\infty$, and is called a CQT-matrix. Given a function $f(x)$ and a CQT matrix $M$, we provide conditions under which $f(M)$ is well defined and is a CQT matrix. Moreover, we introduce a parametrization of CQT matrices and algorithms for the computation of $f(M)$. We treat the case where $f(x)$ is assigned in terms of power series and the case where $f(x)$ is defined in terms of a Cauchy integral. This analysis is applied also to finite matrices which can be written as the sum of a Toeplitz matrix and of a low rank correction.

Dario A. Bini, Beatrice Meini

Let $a(z)=\sum_{i\in\mathbb Z}a_iz^i$ be a complex valued function defined for $|z|=1$, such that $\sum_{i\in\mathbb Z}|ia_i|<\infty$, and let $E=(e_{i,j})_{i,j\in\mathbb {Z}^+}$ be such that $\sum_{i,j\in\mathbb{Z}^+}|e_{i,j}|<\infty$. A semi-infinite quasi-Toeplitz matrix is a matrix of the kind $A=T(a)+E$, where $T(a)=(t_{i,j})_{i,j\in\mathbb{Z}^+}$ is the semi-infinite Toeplitz matrix associated with the symbol $a(z)$, that is, $t_{i,j}=a_{j-i}$ for $i,j\in\mathbb Z^+$. We analyze theoretical and computational properties of the exponential of $A$. More specifically, it is shown that $\exp(A)=T(\exp(a))+F$ where $F=(f_{i,j})_{i,j\in\mathbb{Z}^+}$ is such that $\sum_{i,j\in\mathbb{Z}^+}|f_{i,j}|$ is finite, i.e., $\exp(A)$ is a semi-infinite quasi-Toeplitz matrix as well, and an effective algorithm for its computation is given. These results can be extended from the function $\exp(z)$ to any function $f(z)$ satisfying mild conditions, and can be applied to finite quasi-Toeplitz matrices.

Dario A. Bini, Beatrice Meini, Federico Poloni

We present some advances, both from a theoretical and from a computational point of view, on a quadratic vector equation (QVE) arising in Markovian Binary Trees. Concerning the theoretical advances, some irreducibility assumptions are relaxed, and the minimality of the solution of the QVE is expressed in terms of properties of the Jacobian of a suitable function. From the computational point of view, we elaborate on the Perron vector-based iteration proposed in [http://arxiv.org/abs/1006.0577]. In particular we provide a condition which ensures that the Perron iteration converges to the sought solution of the QVE. Moreover we introduce a variant of the algorithm which consists in applying the Newton method instead of a fixed-point iteration. This method has the same convergence behaviour as the Perron iteration, since its convergence speed increases for close-to-critical problems. Moreover, unlike the Perron iteration, the method has a quadratic convergence. Finally, we show that it is possible to alter the bilinear form defining the QVE in several ways without changing the solution. This modification has an impact on convergence speed of the algorithms.

Dario A. Bini, Sarah Dendievel, Guy Latouche et al.

We consider the Poisson equation $(I-P)\boldsymbol{u}=\boldsymbol{g}$, where $P$ is the transition matrix of a Quasi-Birth-and-Death (QBD) process with infinitely many levels, $\bm g$ is a given infinite dimensional vector and $\bm u$ is the unknown. Our main result is to provide the general solution of this equation. To this purpose we use the block tridiagonal and block Toeplitz structure of the matrix $P$ to obtain a set of matrix difference equations, which are solved by constructing suitable resolvent triples.

Dario A. Bini, Guy Latouche, Beatrice Meini

We revisit the shift technique applied to Quasi-Birth and Death (QBD) processes (He, Meini, Rhee, SIAM J. Matrix Anal. Appl., 2001) by bringing the attention to the existence and properties of canonical factorizations. To this regard, we prove new results concerning the solutions of the quadratic matrix equations associated with the QBD. These results find applications to the solution of the Poisson equation for QBDs.

Dario A. Bini, Beatrice Meini

Given a square matrix $A$, Brauer's theorem [Duke Math. J. 19 (1952), 75--91] shows how to modify one single eigenvalue of $A$ via a rank-one perturbation, without changing any of the remaining eigenvalues. We reformulate Brauer's theorem in functional form and provide extensions to matrix polynomials and to matrix Laurent series $A(z)$ together with generalizations to shifting a set of eigenvalues. We provide conditions under which the modified function $\widetilde A(z)$ has a canonical factorization $\widetilde A(z)=\widetilde U(z)\widetilde L(z^{-1})$ and we provide explicit expressions of the factors $\widetilde U(z)$ and $\widetilde L(z)$. Similar conditions and expressions are given for the factorization of $\widetilde A(z^{-1})$. Some applications are discussed.

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

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$.

Dario A. Bini, Bruno Iannazzo, Beatrice Meini et al.

Component-wise accurate algorithms for computing the principal square root of an M-matrix are designed in terms of triplet representations. A triplet representation of an M-matrix $A$ is the triple $(P, {\bf u},{\bf v})$, where the matrix $P$ is such that $p_{ij}=-a_{ij}$ for $i\ne j$, $p_{ii}=0$, and ${\bf u}>0$, ${\bf v}\ge 0$ are two vectors such that $A{\bf u}={\bf v}$. It is shown that if $A$ is an M-matrix representable by a triplet, then its principal square root exists and is an M-matrix represented by a triplet as well. New versions of the Cyclic Reduction and the Incremental Newton iterations are provided in terms of triplets, to compute the principal matrix square root of $A$. It is shown that these algorithms are component-wise numerically stable independently of the singularity of $A$ and of its condition number. Numerical experiments are shown to confirm the component-wise stability.

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

Dario A. Bini, Federico Poloni

We review some known inclusion results for the roots of a polynomial, and adapt them to a conjecture recently presented by S. A. Vavasis. In particular, we provide strict upper and lower bounds to the distance of the closest root of a polynomial p(z) from a given root of p'(z).