GR decompositions and their relations to Cholesky-like factorizations
GGR decompositions and their relations toCholesky-like factorizations
Peter Benner * Carolin Penke ∗ † June 12, 2020
Abstract
For a given matrix, we are interested in computing GR decompositions A = GR ,where G is an isometry with respect to given scalar products. The orthogonal QRdecomposition is the representative for the Euclidian scalar product. For a signaturematrix, a respective factorization is given as the hyperbolic QR decomposition. Con-sidering a skew-symmetric matrix leads to the symplectic QR decomposition. Thestandard approach for computing GR decompositions is based on the successive elim-ination of subdiagonal matrix entries. For the hyperbolic and symplectic case, thisapproach does in general not lead to a satisfying numerical accuracy. An alternativeapproach computes the QR decomposition via a Cholesky factorization, but also hasbad stability. It is improved by repeating the procedure a second time. In the same way,the hyperbolic and the symplectic QR decomposition are related to the LDL T and askew-symmetric Cholesky-like factorization. We show that methods exploiting thisconnection can provide better numerical stability than elimination-based approaches. Bilinear forms on R m with respect to a nonsingular matrix M ∈ R m × m are defined as h x, y i M = x T M y [1]. The adjoint of a matrix A ∈ R m × m is given as A ? M ∈ R m × m and is uniquely defined by h Ax, y i M = h x, A ? M y i M for all x, y ∈ R m . The adjointgeneralizes the transpose . T . It holds A ? M = M − A T M . Similarly, the adjoint of arectangular matrix A ∈ R m × n is given with respect to two bilinear forms induced bymatrices M ∈ R m × m and N ∈ R n × n as A ? M,N [2]. It is defined by satisfying the iden-tity h Ax, y i M = h x, A ? M,N y i N and it holds A ? M,N = N − A T M . We are interestedin computing decompositions A = GR ∈ R m × n , m ≥ n , where G ∈ R m × n is an ( M, N ) -isometry, i.e. G ? M,N G = I n and R ∈ R n × n [3]. ( M, N ) -isometries are use-ful for devising structure-preserving methods, for example in the context of eigenvaluecomputations [4]. This work considers bilinear forms in real space but the theory is eas-ily extended to complex space or sesquilinear forms. The most well known represen-tative of this class of decompositions is the (thin) QR decomposition. Here M = I m , N = I n and R is upper triangular. With respect to these matrices, an isometry is a * Computational Methods in Systems and Control Theory, Max Planck Institute forDynamics of Complex Technical Systems, Sandtorstr. 1, 39106 Magdeburg, Germany † Corresponding author: [email protected] a r X i v : . [ m a t h . NA ] J un atrix with orthonormal columns. Typically, the QR decomposition is computed ina stable fashion by successively eliminating subdiagonal entries of the matrix usingorthogonal transformations. The decomposition has a well known connection to theCholesky factorization. Let A have full column rank. It holds that A = QR is a thinQR decomposition if and only if R defines a Cholesky decomposition R T R = A T A .Computing Q := AR − provides an alternative to the column elimination approach.For tall and skinny matrices, this method has a much lower computational effort but isknown to be unstable. However, the stability can be drastically improved by doing asecond repetition, i.e. compute the QR decomposition of Q [5]. This is also done in thecontext of H -matrices [6]. In this work, we investigate whether this observation alsoholds for other QR-like decompositions. We now consider a scalar product induced by a signature matrix Σ m = diag( σ , . . . , σ m ) ,where σ , . . . , σ m ∈ { +1 , − } . A (Σ m , Σ n ) -isometry H is called hyperbolic and ful-fills the property H T Σ m H = Σ n . For a given Σ m , the hyperbolic QR decomposition A = HR , where H ∈ R m × n , R ∈ R n × n upper triangular, exists if all principal subma-trices of A are nonsingular [3]. It can be computed via successive column elimination,similar to the orthogonal case [7]. The diagonal values of Σ n are determined by theused transformations and are a subset of the diagonal values of Σ m . The role of theCholesky factorization is now played by the LDL T factorization. A = HR is a hyper-bolic QR decomposition with respect to Σ m and Σ n if and only if R T Σ n R = A T Σ m A gives a scaled LDL T factorization. As the computation of the LDL T factorizationcan be unstable, one typically relies on the slightly altered Bunch–Kaufman factoriza-tion [8]. Here, D is allowed to have × diagonal blocks and pivoting is introducedin form of a permutation P . Using this factorization as a starting point, we arrive at adifferent variant of the HR decomposition, called indefinite QR decomposition in [9].It can be computed from the Bunch-Kaufman (BK) factorization in the following way.1. Compute BK factorization A T Σ A =: P LDL T P T . 2. Diagonalize D =: V Λ V T , Σ n := sign(Λ) R := | Λ | V T L T P T , H := AR − In the resulting QR-like decomposition, R is no longer upper triangular, but blockupper-triangular with permuted columns. In Figure 1 we see how applying this methoda second time to H affects the numerical accuracy.We now consider scalar products induced by J m = (cid:20) I m − I m (cid:21) . A ( J m , J n ) -isometry S ∈ R m × n fulfills the property S T J m S = J n and is called symplectic.The symplectic QR decomposition A = SR can again be computed by the successiveintroduction of zeros in the columns using symplectic transformations. In contrastto the hyperbolic and orthogonal QR decomposition, R is not upper triangular butof the form R = " (cid:64)(cid:64)(cid:64) (cid:64)(cid:64) (cid:64)(cid:64)(cid:64)(cid:64) (cid:64)(cid:64) (cid:64)(cid:64)(cid:64) = P T s ˆ RP s , where ˆ R is upper triangular with × diagonal blocks on the diagonal and P s = [ e , e , . . . , e n − , e , e . . . , e n ] is theperfect shuffle. This variant of the symplectic QR decomposition corresponds to theskew-symmetric Cholesky factorization described in [10, 11]. A = SR is a symplecticQR decomposition if and only if ˆ R T ( P s J n P T s ) ˆ R is a Cholesky-like decomposition of2he skew-symmetric matrix P s A T J m AP T s . Similar to the hyperbolic case, pivotingin form of a permutation matrix P can be introduced to increase the stability of theCholesky-like decomposition, leading to an altered SR decomposition AP T P s = SR .A Cholesky-based computation takes the form1. Compute Cholesky-like factorization A T J m A =: P ˆ R T P s J n P T s ˆ RP . 2. R := P T s ˆ RP s . 3. S := AP T P S R − .Again there exists the possibility to repeat the procedure for the computed symplecticfactor S . − − − cond(A) HR residual k A − HR k F / k A k F − − − cond(A) k H T Σ H − ˆΣ k F Column elimination 1 x
LDL T LDL T − − − cond(A) SR residual k A − SR k F / k A k F − − cond(A) k S T JS − J k F Column Elimination 1x Chol-like, no pivoting2x Chol-like, no pivoting 1x Chol-like with pivoting2x Chol-like with pivoting in first passFigure 1: Numerical results for computing decompositions of randomly generated ma-trices of size × (HR) or × (SR).Figure 1 shows how the accuracy of HR and SR decompositions can be improvedby computing them via the LDL T and the Cholesky-like decomposition. Computingthe decompositions via factorizations employing pivoting with two iterations leads toaccuracies that do not derail for badly conditioned matrices. References [1] N. Higham, D. Mackey, N. Mackey, and F. Tisseur. Functions preserving matrixgroups and iterations for the matrix square root.
SIAM Journal on Matrix Analysisand Applications , 26(3):849–877, 2005.32] N. Higham, C. Mehl, and F. Tisseur. The canonical generalized polar decom-position.
SIAM Journal on Matrix Analysis and Applications , 31(4):2163–2180,2010.[3] D. Watkins.
The Matrix Eigenvalue Problem . Society for Industrial and AppliedMathematics, 2007.[4] Angelika Bunse-Gerstner. Matrix factorizations for symplectic QR -like methods. Linear Algebra and its Applications , 83:49–77, 1986.[5] Yamamoto Y., Y. Nakatsukasa, Y. Yanagisawa, and T. Fukaya. Roundoff erroranalysis of the Cholesky QR2 algorithm.
Electronic Transactions on NumericalAnalysis , 44:306–326, 2015.[6] M. Lintner. The eigenvalue problem for the 2D Laplacian in H -matrix arithmeticand application to the heat and wave equation. Computing. Archives for ScientificComputing , 72(3-4):293–323, 2004.[7] Ivo Houtzager. JQR/JRQ/JQL/JLQ factorizations.
MATLAB Central File Ex-change , July 2015. Retrieved February 12, 2020.[8] C. Ashcraft, R. G. Grimes, and J. G. Lewis. Accurate symmetric indefinite linearequation solvers.
SIAM Journal on Matrix Analysis and Applications , 20(2):513–561, 1999.[9] S. Singer. Indefinite QR factorization.
BIT. Numerical Mathematics , 46(1):141–161, 2006.[10] P. Benner, R. Byers, H. Fassbender, V. Mehrmann, and D. Watkins. Cholesky-likefactorizations of skew-symmetric matrices.
Electronic Transactions on Numeri-cal Analysis , 11:85–93, 2000.[11] James R. Bunch. A note on the stable decompostion of skew-symmetric matrices.