Benner, P.

Penke, C.

English

MPG.PuRe

GR decompositions and their relations toCholesky-like factorizationsPeter Benner* Carolin Penke∗ †June 12, 2020AbstractFor 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 LDLT and askew-symmetric Cholesky-like factorization. We show that methods exploiting thisconnection can provide better numerical stability than elimination-based approaches.1 Introduction and PreliminariesBilinear forms on Rm with respect to a nonsingular matrix M ∈ Rm×m are defined as〈x, y〉M = xTMy [1]. The adjoint of a matrix A ∈ Rm×m is given as A?M ∈ Rm×mand is uniquely defined by 〈Ax, y〉M = 〈x,A?M y〉M for all x, y ∈ Rm. The adjointgeneralizes the transpose .T. It holds A?M = M−1ATM . Similarly, the adjoint of arectangular matrix A ∈ Rm×n is given with respect to two bilinear forms induced bymatricesM ∈ Rm×m andN ∈ Rn×n asA?M,N [2]. It is defined by satisfying the iden-tity 〈Ax, y〉M = 〈x,A?M,N y〉N and it holds A?M,N = N−1ATM . We are interestedin computing decompositions A = GR ∈ Rm×n, m ≥ n, where G ∈ Rm×n is an(M,N)-isometry, i.e. G?M,NG = In and R ∈ Rn×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 = Im,N = In 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: penke@mpi-magdeburg.mpg.de1arXiv:2006.06558v1  [math.NA]  11 Jun 2020matrix 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 RTR = ATA.Computing Q := AR−1 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 ofQ [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.2 GR decompositions and Cholesky-like factorizationsWe now consider a scalar product induced by a signature matrix Σm = diag(σ1, . . . , σm),where σ1, . . . , σm ∈ {+1,−1}. A (Σm,Σn)-isometry H is called hyperbolic and ful-fills the property HTΣmH = Σn. For a given Σm, the hyperbolic QR decompositionA = HR, whereH ∈ Rm×n,R ∈ Rn×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 LDLT factorization. A = HR is a hyper-bolic QR decomposition with respect to Σm and Σn if and only if RTΣnR = ATΣmAgives a scaled LDLT factorization. As the computation of the LDLT factorizationcan be unstable, one typically relies on the slightly altered Bunch–Kaufman factoriza-tion [8]. Here, D is allowed to have 2 × 2 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 ATΣA =: PLDLTPT. 2. Diagonalize D =:V ΛV T , Σn := sign(Λ)3. R := |Λ| 12V TLTPT, H := AR−1In 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 Jm =[0 Im−Im 0]. A (Jm, Jn)-isometry S ∈ R2m×2n fulfills the property STJmS = Jn 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 =[@@@@00@@@00@@@]= PTs RˆPs, where Rˆ is upper triangular with 2 × 2diagonal blocks on the diagonal and Ps = [e1, e3, . . . , e2n−1, e2, e4 . . . , e2n] 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(PsJnPTs )Rˆ is a Cholesky-like decomposition of2the skew-symmetric matrix PsATJmAPTs . 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 APTPs = SR.A Cholesky-based computation takes the form1. Compute Cholesky-like factorization ATJmA =: PRˆTPsJnPTs RˆP . 2. R :=PTs RˆPs. 3. S := APTPSR−1.Again there exists the possibility to repeat the procedure for the computed symplecticfactor S.102 104 106 10810−1510−1210−9cond(A)HR residual ‖A−HR‖F /‖A‖F102 104 106 10810−1510−810−1cond(A)‖HTΣH − Σˆ‖FColumn elimination 1 x LDLT 2 x LDLT102 104 106 10810−1410−1210−10cond(A)SR residual ‖A− SR‖F /‖A‖F102 104 106 10810−1510−5105cond(A)‖STJS − J‖FColumn 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 500× 500 (HR) or 1000× 1000 (SR).Figure 1 shows how the accuracy of HR and SR decompositions can be improvedby computing them via the LDLT 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.3[2] 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 symplecticQR-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 inH-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.Mathematics of Computation, 38(158):475–479, 1982.4

GR Decompositions and their Relations to Cholesky-like Factorizations

Abstract

Similar works

Full text

Available Versions

MPG.PuRe