WWW home page: http://www.ba.cnr.it/~irmanm21/Welcome.html Abstract. Classically, several numeric... more WWW home page: http://www.ba.cnr.it/~irmanm21/Welcome.html Abstract. Classically, several numerical linear algebra problems are solved by means of tridiagonal (symmetric) and Hessenberg (unsymmet- ric) matrices. In this talk, it will be shown that a similar role can be played by semiseparable matrices. A matrix is called semiseparable if all submatrices that can be taken out of the lower triangular part (the main
A code implementing the ST-HOSVD, using the Tensor Toolbox. See our paper on "A new truncati... more A code implementing the ST-HOSVD, using the Tensor Toolbox. See our paper on "A new truncation strategy for the higher-order singular value decomposition" for details.
Electronic transactions on numerical analysis ETNA
In this paper we present a survey of various algorithms for computing matrix geometric means and ... more In this paper we present a survey of various algorithms for computing matrix geometric means and derive new second-order optimization algorithms to compute the Karcher mean. These new algorithms are constructed using the standard definition of the Riemannian Hessian. The survey includes the ALM list of desired properties for a geometric mean, the analytical expression for the mean of two matrices, algorithms based on the centroid computation in Euclidean (flat) space, and Riemannian optimization techniques to compute the Karcher mean (preceded by a short introduction into differential geometry). A change of metric is considered in the optimization techniques to reduce the complexity of the structures used in these algorithms. Numerical experiments are presented to compare the existing and the newly developed algorithms. We conclude that currently first-order algorithms are best suited for this optimization problem as the size and/or number of the matrices increase.
We present a numerical example illustrating that the deflation procedure in Francis's implicitly ... more We present a numerical example illustrating that the deflation procedure in Francis's implicitly shifted QR algorithm can be improved by a deflation criteria based on the QR decomposition of the upper Hessenberg matrix.
SIAM Journal on Matrix Analysis and Applications, 2015
A stable algorithm to compute the roots of polynomials is presented. The roots are found by compu... more A stable algorithm to compute the roots of polynomials is presented. The roots are found by computing the eigenvalues of the associated companion matrix by Francis's implicitly-shifted QR algorithm. A companion matrix is an upper Hessenberg matrix that is unitary-plus-rank-one, that is, it is the sum of a unitary matrix and a rank-one matrix. These properties are preserved by iterations of Francis's algorithm, and it is these properties that are exploited here. The matrix is represented as a product of 3n − 1 Givens rotators plus the rank-one part, so only O(n) storage space is required. In fact, the information about the rank-one part is also encoded in the rotators, so it is not necessary to store the rank-one part explicitly. Francis's algorithm implemented on this representation requires only O(n) flops per iteration and thus O(n 2) flops overall. The algorithm is described, backward stability is proved under certain conditions on the polynomial coefficients, and an extensive set of numerical experiments is presented. The algorithm is shown to be about as accurate as the (slow) Francis QR algorithm applied to the companion matrix without exploiting the structure. It is faster than other fast methods that have been proposed, and its accuracy is comparable or better.
... 3. William L. Briggs , Van Emden Henson , Steve F. McCormick, A multigrid tutorial: second ed... more ... 3. William L. Briggs , Van Emden Henson , Steve F. McCormick, A multigrid tutorial: second edition, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2000. 4. [4] O. Buneman. ... 32. [32] WM Pickering and PJ Harley FFT solution of the Robbins problem, IMAJ. Num. ...
There are many classical results in which orthogonal vectors stemming from Krylov subspaces are l... more There are many classical results in which orthogonal vectors stemming from Krylov subspaces are linked to short recurrence relations, e.g., three-terms recurrences for Hermitian and short rational recurrences for unitary matrices. These recurrence coefficients can be captured in a Hessenberg matrix, whose structure reflects the relation between the spectrum of the original matrix and the recurrences. The easier the recurrences, the faster the orthogonal vectors can be computed possibly resulting in computational savings in the design of, e.g., iterative solvers. In this article we focus on multiple recursions, i.e., the (j + 1)st orthogonal vector satisfies q j+1 = j i=j−m ρ j,i Aq i − j i=j− γ j,i q i , with ρ j,i , γ j,i scalars, and A the matrix defining the Krylov space. Though many compelling results are around, the theoretical proof of the existence of this structure, and the associated construction of algorithms to fastly retrieve the orthogonal vectors are rapidly becoming technically involved. In this article we first review classical results on short multiple recurrences for normal matrices whose Hermitian conjugate can be written as a 'low degree' rational function of itself. Instead of considering innerproduct relations, we reformulate this theory to a matrix setting. Moreover, the matrix building blocks allow us to also derive multiple recursions for B-normal matrices, normal matrices whose eigenvalues lie on the union of curves in the plane, normal matrices perturbed by a low rank, and normal matrices A satisfying a 'low degree' relation s(A H) = p(A)/q(A). The theoretical results on the structure available in the Hessenberg matrix lead to a new manner to compute the orthogonal vectors. The numerical experiments illustrate, however, that straightforward algorithms exploiting this structure are very sensitive, and more research is required to develop robust algorithms.
In this manuscript a new type of QR-iteration will be presented. Each step of this new iteration ... more In this manuscript a new type of QR-iteration will be presented. Each step of this new iteration consists of two substeps. In the explicit version, rst an RQ-factorization of the initial matrix A κI = RQ will be computed, followed by a QR-factorization of the matrix (A σI)Q H. Applying the unitary similarity transformation dened by the QR-factorization of the transformed matrix (A σI)Q H , will yield interesting convergence properties. It will be shown that the con- vergence behavior is related to a subspace iteration based on a rational function in A namely (A σI)(A κI) 1 . Convergence properties of this new iteration will be investigated and exam- ples will be presented, illustrating the effectiveness of this approach with respect to some specic classes of matrices.
SIAM Journal on Matrix Analysis and Applications, 2014
The Schmidt-Eckart-Young theorem for matrices states that the optimal rankr approximation of a ma... more The Schmidt-Eckart-Young theorem for matrices states that the optimal rankr approximation of a matrix is obtained by retaining the first r terms from the singular value decomposition of that matrix. This paper considers a generalization of this optimal truncation property to the rank decomposition (Candecomp/Parafac) of tensors and establishes a necessary orthogonality condition. We prove that this condition is not satisfied at least by an open set of positive Lebesgue measure in complex tensor spaces. It is proved, moreover, that for complex tensors of small rank this condition can be satisfied only by a set of tensors of Lebesgue measure zero. Finally, we demonstrate that generic tensors in cubic tensor spaces are not optimally truncatable.
Recent progress in signal processing and estimation has gener- ated considerable interest in the ... more Recent progress in signal processing and estimation has gener- ated considerable interest in the problem of computing the smallest eigenvalue of symmetric positive definite Toeplitz matrices. Several algorithms have been proposed in the literature. Many of them com- pute the smallest eigenvalue in an iterative fashion, relying on the Levinson-Durbin solution of sequences of Yule-Walker systems. Exploiting the properties of
Various generalizations of companion matrices to companion pencils are presented. Companion matri... more Various generalizations of companion matrices to companion pencils are presented. Companion matrices link to monic polynomials, whereas companion pencils do not require monicity of the corresponding polynomial. In the classical companion pencil case (A, B) only the coefficient of the highest degree appears in B's lower right corner. We will show, however, that all coefficients of the polynomial can be distributed over both A and B creating additional flexibility. Companion matrices admit a so-called Fiedler factorization into essentially 2 × 2 matrices. These Fiedler factors can be reordered without affecting the eigenvalues (which equal the polynomial roots) of the assembled matrix. We will propose a generalization of this factorization and the reshuffling property for companion pencils. Special examples of the factorizations and extensions to matrix polynomials and product eigenvalue problems are included.
In this talk we will discuss numerical techniques for computing a matrix geometric mean, arising ... more In this talk we will discuss numerical techniques for computing a matrix geometric mean, arising in quite some applications. The generalization of the scalar geometric mean to matrices has proven to be non-trivial, due to noncommutativity, and as a consequence lacks a unique definition. We will first analyse the properties of some existing algorithms for calculating geometric matrix means. These algorithms were constructed from geometric and axiomatic points of view to satisfy some necessary and desirable properties of the averaging of symmetric positive definite matrices. It is shown that these properties do not specify a unique definition for the geometric matrix mean, resulting in different valid algorithms of varying performance. Afterwards we compare these results to several manifold optimization techniqus. These means are constructed from a center-of-mass-principle by optimizing on the manifold of positive definite matrices, utilizing thereby differential geometric tools.
It has been shown that approximate extended Krylov subspaces can be computed-under certain assump... more It has been shown that approximate extended Krylov subspaces can be computed-under certain assumptions-without any explicit inversion or system solves. Instead the necessary products A −1 v are obtained in an implicit way retrieved from an enlarged Krylov subspace. In this paper this approach is generalized to rational Krylov subspaces, which contain besides poles at infinite and zero also finite non-zero poles. Also an adaption of the algorithm to the block and the symmetric case is presented. For all variants of the algorithm numerical experiments underpin the power of the new approach. Rational Krylov subspaces can be used, e.g., to approximate matrix functions or the solutions of matrix equations.
WWW home page: http://www.ba.cnr.it/~irmanm21/Welcome.html Abstract. Classically, several numeric... more WWW home page: http://www.ba.cnr.it/~irmanm21/Welcome.html Abstract. Classically, several numerical linear algebra problems are solved by means of tridiagonal (symmetric) and Hessenberg (unsymmet- ric) matrices. In this talk, it will be shown that a similar role can be played by semiseparable matrices. A matrix is called semiseparable if all submatrices that can be taken out of the lower triangular part (the main
A code implementing the ST-HOSVD, using the Tensor Toolbox. See our paper on "A new truncati... more A code implementing the ST-HOSVD, using the Tensor Toolbox. See our paper on "A new truncation strategy for the higher-order singular value decomposition" for details.
Electronic transactions on numerical analysis ETNA
In this paper we present a survey of various algorithms for computing matrix geometric means and ... more In this paper we present a survey of various algorithms for computing matrix geometric means and derive new second-order optimization algorithms to compute the Karcher mean. These new algorithms are constructed using the standard definition of the Riemannian Hessian. The survey includes the ALM list of desired properties for a geometric mean, the analytical expression for the mean of two matrices, algorithms based on the centroid computation in Euclidean (flat) space, and Riemannian optimization techniques to compute the Karcher mean (preceded by a short introduction into differential geometry). A change of metric is considered in the optimization techniques to reduce the complexity of the structures used in these algorithms. Numerical experiments are presented to compare the existing and the newly developed algorithms. We conclude that currently first-order algorithms are best suited for this optimization problem as the size and/or number of the matrices increase.
We present a numerical example illustrating that the deflation procedure in Francis's implicitly ... more We present a numerical example illustrating that the deflation procedure in Francis's implicitly shifted QR algorithm can be improved by a deflation criteria based on the QR decomposition of the upper Hessenberg matrix.
SIAM Journal on Matrix Analysis and Applications, 2015
A stable algorithm to compute the roots of polynomials is presented. The roots are found by compu... more A stable algorithm to compute the roots of polynomials is presented. The roots are found by computing the eigenvalues of the associated companion matrix by Francis's implicitly-shifted QR algorithm. A companion matrix is an upper Hessenberg matrix that is unitary-plus-rank-one, that is, it is the sum of a unitary matrix and a rank-one matrix. These properties are preserved by iterations of Francis's algorithm, and it is these properties that are exploited here. The matrix is represented as a product of 3n − 1 Givens rotators plus the rank-one part, so only O(n) storage space is required. In fact, the information about the rank-one part is also encoded in the rotators, so it is not necessary to store the rank-one part explicitly. Francis's algorithm implemented on this representation requires only O(n) flops per iteration and thus O(n 2) flops overall. The algorithm is described, backward stability is proved under certain conditions on the polynomial coefficients, and an extensive set of numerical experiments is presented. The algorithm is shown to be about as accurate as the (slow) Francis QR algorithm applied to the companion matrix without exploiting the structure. It is faster than other fast methods that have been proposed, and its accuracy is comparable or better.
... 3. William L. Briggs , Van Emden Henson , Steve F. McCormick, A multigrid tutorial: second ed... more ... 3. William L. Briggs , Van Emden Henson , Steve F. McCormick, A multigrid tutorial: second edition, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2000. 4. [4] O. Buneman. ... 32. [32] WM Pickering and PJ Harley FFT solution of the Robbins problem, IMAJ. Num. ...
There are many classical results in which orthogonal vectors stemming from Krylov subspaces are l... more There are many classical results in which orthogonal vectors stemming from Krylov subspaces are linked to short recurrence relations, e.g., three-terms recurrences for Hermitian and short rational recurrences for unitary matrices. These recurrence coefficients can be captured in a Hessenberg matrix, whose structure reflects the relation between the spectrum of the original matrix and the recurrences. The easier the recurrences, the faster the orthogonal vectors can be computed possibly resulting in computational savings in the design of, e.g., iterative solvers. In this article we focus on multiple recursions, i.e., the (j + 1)st orthogonal vector satisfies q j+1 = j i=j−m ρ j,i Aq i − j i=j− γ j,i q i , with ρ j,i , γ j,i scalars, and A the matrix defining the Krylov space. Though many compelling results are around, the theoretical proof of the existence of this structure, and the associated construction of algorithms to fastly retrieve the orthogonal vectors are rapidly becoming technically involved. In this article we first review classical results on short multiple recurrences for normal matrices whose Hermitian conjugate can be written as a 'low degree' rational function of itself. Instead of considering innerproduct relations, we reformulate this theory to a matrix setting. Moreover, the matrix building blocks allow us to also derive multiple recursions for B-normal matrices, normal matrices whose eigenvalues lie on the union of curves in the plane, normal matrices perturbed by a low rank, and normal matrices A satisfying a 'low degree' relation s(A H) = p(A)/q(A). The theoretical results on the structure available in the Hessenberg matrix lead to a new manner to compute the orthogonal vectors. The numerical experiments illustrate, however, that straightforward algorithms exploiting this structure are very sensitive, and more research is required to develop robust algorithms.
In this manuscript a new type of QR-iteration will be presented. Each step of this new iteration ... more In this manuscript a new type of QR-iteration will be presented. Each step of this new iteration consists of two substeps. In the explicit version, rst an RQ-factorization of the initial matrix A κI = RQ will be computed, followed by a QR-factorization of the matrix (A σI)Q H. Applying the unitary similarity transformation dened by the QR-factorization of the transformed matrix (A σI)Q H , will yield interesting convergence properties. It will be shown that the con- vergence behavior is related to a subspace iteration based on a rational function in A namely (A σI)(A κI) 1 . Convergence properties of this new iteration will be investigated and exam- ples will be presented, illustrating the effectiveness of this approach with respect to some specic classes of matrices.
SIAM Journal on Matrix Analysis and Applications, 2014
The Schmidt-Eckart-Young theorem for matrices states that the optimal rankr approximation of a ma... more The Schmidt-Eckart-Young theorem for matrices states that the optimal rankr approximation of a matrix is obtained by retaining the first r terms from the singular value decomposition of that matrix. This paper considers a generalization of this optimal truncation property to the rank decomposition (Candecomp/Parafac) of tensors and establishes a necessary orthogonality condition. We prove that this condition is not satisfied at least by an open set of positive Lebesgue measure in complex tensor spaces. It is proved, moreover, that for complex tensors of small rank this condition can be satisfied only by a set of tensors of Lebesgue measure zero. Finally, we demonstrate that generic tensors in cubic tensor spaces are not optimally truncatable.
Recent progress in signal processing and estimation has gener- ated considerable interest in the ... more Recent progress in signal processing and estimation has gener- ated considerable interest in the problem of computing the smallest eigenvalue of symmetric positive definite Toeplitz matrices. Several algorithms have been proposed in the literature. Many of them com- pute the smallest eigenvalue in an iterative fashion, relying on the Levinson-Durbin solution of sequences of Yule-Walker systems. Exploiting the properties of
Various generalizations of companion matrices to companion pencils are presented. Companion matri... more Various generalizations of companion matrices to companion pencils are presented. Companion matrices link to monic polynomials, whereas companion pencils do not require monicity of the corresponding polynomial. In the classical companion pencil case (A, B) only the coefficient of the highest degree appears in B's lower right corner. We will show, however, that all coefficients of the polynomial can be distributed over both A and B creating additional flexibility. Companion matrices admit a so-called Fiedler factorization into essentially 2 × 2 matrices. These Fiedler factors can be reordered without affecting the eigenvalues (which equal the polynomial roots) of the assembled matrix. We will propose a generalization of this factorization and the reshuffling property for companion pencils. Special examples of the factorizations and extensions to matrix polynomials and product eigenvalue problems are included.
In this talk we will discuss numerical techniques for computing a matrix geometric mean, arising ... more In this talk we will discuss numerical techniques for computing a matrix geometric mean, arising in quite some applications. The generalization of the scalar geometric mean to matrices has proven to be non-trivial, due to noncommutativity, and as a consequence lacks a unique definition. We will first analyse the properties of some existing algorithms for calculating geometric matrix means. These algorithms were constructed from geometric and axiomatic points of view to satisfy some necessary and desirable properties of the averaging of symmetric positive definite matrices. It is shown that these properties do not specify a unique definition for the geometric matrix mean, resulting in different valid algorithms of varying performance. Afterwards we compare these results to several manifold optimization techniqus. These means are constructed from a center-of-mass-principle by optimizing on the manifold of positive definite matrices, utilizing thereby differential geometric tools.
It has been shown that approximate extended Krylov subspaces can be computed-under certain assump... more It has been shown that approximate extended Krylov subspaces can be computed-under certain assumptions-without any explicit inversion or system solves. Instead the necessary products A −1 v are obtained in an implicit way retrieved from an enlarged Krylov subspace. In this paper this approach is generalized to rational Krylov subspaces, which contain besides poles at infinite and zero also finite non-zero poles. Also an adaption of the algorithm to the block and the symmetric case is presented. For all variants of the algorithm numerical experiments underpin the power of the new approach. Rational Krylov subspaces can be used, e.g., to approximate matrix functions or the solutions of matrix equations.
Uploads
Papers by Raf Vandebril