Computing the Matrix Square Root in a Matrix Group Nick Higham - PowerPoint PPT Presentation

Computing the Matrix Square Root in a Matrix Group Nick Higham Department of Mathematics The Victoria University of Manchester higham@ma.man.ac.uk http://www.ma.man.ac.uk/~higham/ Joint work with Niloufer Mackey, D. Steven Mackey, and Françoise Tisseur. Matrix Group Iterations – p. 1/20

Matrix Square Root X is a square root of A ∈ C n × n ⇐ ⇒ X 2 = A . Number of square roots may be zero, finite or infinite. Principal Square Root For A with no eigenvalues on R − = { x ∈ R : x ≤ 0 } is the unique square root with spectrum in the open right half-plane. Denoted by A 1 / 2 . Matrix Group Iterations – p. 2/20

Newton’s Method A = ( X + E ) 2 : drop second order term and solve for E . X 0 given, Solve X k E k + E k X k = A − X 2 � k k = 0 , 1 , 2 , . . . X k +1 = X k + E k Prohibitively expensive. Matrix Group Iterations – p. 3/20

Newton’s Method A = ( X + E ) 2 : drop second order term and solve for E . X 0 given, Solve X k E k + E k X k = A − X 2 � k k = 0 , 1 , 2 , . . . X k +1 = X k + E k Prohibitively expensive. But observe that X 0 A = AX 0 ⇒ X k A = AX k for all k . Matrix Group Iterations – p. 3/20

Newton Iteration 2 ( X k + X − 1 X k +1 = 1 k A ) , X 0 = A . If A ∈ C n × n has no eigenvalues on R − then Iterates X k nonsingular, converge quadratically to A 1 / 2 . Related to Newton sign function iterates S k +1 = 1 S 0 = A 1 / 2 2( S k + S − 1 k ) , by X k ≡ A 1 / 2 S k . Matrix Group Iterations – p. 4/20

Newton Iteration 2 ( X k + X − 1 X k +1 = 1 k A ) , X 0 = A . If A ∈ C n × n has no eigenvalues on R − then Iterates X k nonsingular, converge quadratically to A 1 / 2 . Related to Newton sign function iterates S k +1 = 1 S 0 = A 1 / 2 2( S k + S − 1 k ) , by X k ≡ A 1 / 2 S k . Problem : numerical instability (Laasonen, 1958; H, 1996). Matrix Group Iterations – p. 4/20

Newton Variants X k +1 = 1 X k + Y − 1 � � , X 0 = A, DB iteration : k 2 Y k +1 = 1 [Denman–Beavers, 1976] Y k + X − 1 � � , Y 0 = I. k 2 I + M k + M − 1 M k +1 = 1 � � k , M 0 = A, 2 2 Product form DB : X k +1 = 1 2 X k ( I + M − 1 [Cheng-Higham- k ) , X 0 = A, Kenney-Laub, 2001] Y k +1 = 1 2 Y k ( I + M − 1 k ) , Y 0 = I. Y k +1 = − Y k Z − 1 k Y k , Y 0 = I − A, CR iteration : [Meini, 2004] Z k +1 = Z k + 2 Y k +1 , Z 0 = 2( I + A ) . Matrix Group Iterations – p. 5/20

Content ⋆ Characterizations of functions f that preserve automorphism group structure. ⋆ New, numerically stable square root iterations. ⋆ Unified stability analysis of square root iterations based on Fréchet derivatives. Matrix Group Iterations – p. 6/20

Group Background Given nonsingular M and K = R or C , � x T My, real or complex bilinear forms, � x, y � M = sesquilinear forms. x ∗ My, Define automorphism group G = { A ∈ K n × n : � Ax, Ay � M = � x, y � M , ∀ x, y ∈ K n } . Adjoint A⋆ of A ∈ K n × n wrt �· , ·� M defined by � Ax, y � M = � x, A⋆y � M ∀ x, y ∈ K n . � M − 1 A T M, for bilinear forms, A⋆ = Can show: M − 1 A ∗ M, for sesquilinear forms. G = { A ∈ K n × n : A⋆ = A − 1 } . Matrix Group Iterations – p. 7/20

Some Automorphism Groups A⋆ Space M Automorphism group , G Groups corresponding to a bilinear form R n A T I Real orthogonals C n A T I Complex orthogonals Σ p,q A T Σ p,q R n Σ p,q Pseudo-orthogonals RA T R R n Real perplectics R − JA T J R 2 n Real symplectics J − JA T J C 2 n Complex symplectics J Groups corresponding to a sesquilinear form C n Unitaries I A ∗ C n Pseudo-unitaries Σ p,q Σ p,q A ∗ Σ p,q C 2 n − JA ∗ J Conjugate symplectics J � � ...     1 0 In Ip 0 R = , J =  , Σp,q =     − In    0 0 − Iq 1 Matrix Group Iterations – p. 8/20

Bilinear Forms Theorem 1 (a) For any f and A ∈ K n × n , f ( A⋆ ) = f ( A ) ⋆ . (b) For A ∈ G , f ( A ) ∈ G iff f ( A − 1 ) = f ( A ) − 1 . Proof . (a) We have f ( A⋆ ) = f ( M − 1 A T M ) = M − 1 f ( A T ) M = M − 1 f ( A ) T M = f ( A ) ⋆. (b) For A ∈ G , consider f ( A ) ⋆ f ( A⋆ ) = � f ( A − 1 ) Matrix Group Iterations – p. 9/20

Bilinear Forms Theorem 1 (a) For any f and A ∈ K n × n , f ( A⋆ ) = f ( A ) ⋆ . (b) For A ∈ G , f ( A ) ∈ G iff f ( A − 1 ) = f ( A ) − 1 . Proof . (a) We have f ( A⋆ ) = f ( M − 1 A T M ) = M − 1 f ( A T ) M = M − 1 f ( A ) T M = f ( A ) ⋆. (b) For A ∈ G , consider f ( A ) ⋆ f ( A⋆ ) = � � f ( A ) − 1 f ( A − 1 ) = Matrix Group Iterations – p. 9/20

Implications For bilinear forms, f preserves group structure of A when f ( A − 1 ) = f ( A ) − 1 . This condition holds for all A for Matrix sign function , sign( A ) = A ( A 2 ) − 1 / 2 . Any matrix power A α , subject to suitable choice of branches. In particular, the principal matrix p th root A 1 /p ( p ∈ Z + , Λ ( A ) ∩ R − = ∅ ): unique X such that 1. X p = A . 2. − π/p < arg( λ ( X )) < π/p . Matrix Group Iterations – p. 10/20

Group Newton Iteration Theorem 2 Let A ∈ G ( any group ) , Λ ( A ) ∩ R − = ∅ , and Y k +1 = 1 2( Y k + Y − ⋆ ) k = 1 Y 1 = 1 2( Y k + M − 1 Y − T M ) , 2( I + A ) . k Then Y k → A 1 / 2 quadratically and Y k ≡ X k ( k ≥ 1 ) , where X k are the Newton iterates: X k +1 = 1 2 ( X k + X − 1 k A ) , X 0 = A . Matrix Group Iterations – p. 11/20

Group Newton Iteration Theorem 2 Let A ∈ G ( any group ) , Λ ( A ) ∩ R − = ∅ , and Y k +1 = 1 2( Y k + Y − ⋆ ) k = 1 Y 1 = 1 2( Y k + M − 1 Y − T M ) , 2( I + A ) . k Then Y k → A 1 / 2 quadratically and Y k ≡ X k ( k ≥ 1 ) , where X k are the Newton iterates: X k +1 = 1 2 ( X k + X − 1 k A ) , X 0 = A . Cf. Cardoso, Kenney & Silva Leite (2003, App. Num. Math.): bilinear forms with M T = ± M , M T M = I . H (2003, SIREV): M = Σ p,q . Matrix Group Iterations – p. 11/20

Generalized Polar Decomposition Theorem 2 For any group G , A ∈ K n × n has a unique generalized polar decomposition A = WS where S⋆ = S, ( i.e., W⋆ = W − 1 ) , W ∈ G and Λ ( S ) ∈ open right half-plane ( i.e., sign( S ) = I ) iff ( A⋆ ) ⋆ = A and Λ ( A⋆A ) ∩ R − = ∅ . Note ( A⋆ ) ⋆ = A holds for all G in the earlier table. Other gpd’s exist with different conditions on Λ ( S ) (Rodman & co-authors). Matrix Group Iterations – p. 12/20

GPD Iteration & Square Root Theorem 3 Suppose the iteration X k +1 = X k h ( X 2 k ) , X 0 = A converges to sign( A ) with order m . If A has the generalized polar decomposition A = WS w.r.t. a scalar product then Y k +1 = Y k h ( Y ⋆ k Y k ) , Y 0 = A converges to W with order of convergence m . Matrix Group Iterations – p. 13/20

GPD Iteration & Square Root Theorem 3 Suppose the iteration X k +1 = X k h ( X 2 k ) , X 0 = A converges to sign( A ) with order m . If A has the generalized polar decomposition A = WS w.r.t. a scalar product then Y k +1 = Y k h ( Y ⋆ k Y k ) , Y 0 = A converges to W with order of convergence m . Theorem 4 Let G be any automorphism group and A ∈ G . If Λ ( A ) ∩ R − = ∅ then I + A = WS is a generalized polar decomposition with W = A 1 / 2 and S = A − 1 / 2 + A 1 / 2 . Matrix Group Iterations – p. 13/20

Application Newton Sign : X k +1 = X k · 1 2 ( I + ( X 2 k ) − 1 ) ≡ X k h ( X 2 k ) , X 0 = A . Group sqrt : Y k +1 = 1 k Y k ) − 1 ) = 1 2 Y k ( I + ( Y ⋆ 2( Y k + Y − ⋆ ) , Y 0 = I + A. k Schulz Sign : X k +1 = X k · 1 2 (3 I − X 2 k ) ≡ X k h ( X 2 k ) , X 0 = A . Group Schulz : Y k +1 = 1 2 Y k (3 I − Y ⋆ k Y k ) , Y 0 = I + A. Matrix Group Iterations – p. 14/20

Class of Square Root Iterations Theorem 5 Suppose the iteration X k +1 = X k h ( X 2 k ) , X 0 = A converges to sign( A ) with order m . If Λ ( A ) ∩ R − = ∅ and Y k +1 = Y k h ( Z k Y k ) , Y 0 = A, Z k +1 = h ( Z k Y k ) Z k , Z 0 = I, then Y k → A 1 / 2 and Z k → A − 1 / 2 as k → ∞ , both with order m , and Y k = AZ k for all k . Moreover, if X ∈ G implies Xh ( X 2 ) ∈ G then A ∈ G implies Y k ∈ G and Z k ∈ G for all k . �� �� � A 1 / 2 � 0 A 0 Proof makes use of sign = . A − 1 / 2 I 0 0 Newton sign leads to DB iteration. Matrix Group Iterations – p. 15/20

Padé Square Root Iterations Example: structure-preserving cubic: Y k +1 = Y k (3 I + Z k Y k )( I + 3 Z k Y k ) − 1 , Y 0 = A, Z k +1 = (3 I + Z k Y k )( I + 3 Z k Y k ) − 1 Z k , Z 0 = I. If Λ ( A ) ∩ R − = ∅ then • Y k → A 1 / 2 and Z k → A − 1 / 2 cubically, • A ∈ G ⇒ X k ∈ G , Y k ∈ G for all k . Matrix Group Iterations – p. 16/20

Stability Define X k +1 = g ( X k ) to be stable in nbhd of fixed point X = g ( X ) if for X 0 := X + H 0 , with arbitrary error H 0 , the H k := X k − X satisfy H k +1 = L X ( H k ) + O ( � H k � 2 ) , where L X is a linear operator with bounded powers . Theorem 6 Consider the mapping � Y h ( ZY ) � G ( Y, Z ) = , h ( ZY ) Z where X k +1 = X k h ( X 2 k ) is any superlinear matrix sign iteration. Any P = ( B, B − 1 ) is a fixed point for G , and E − BFB dG P ( E, F ) = 1 � � . F − B − 1 EB − 1 2 dG P is idempotent , that is, dG P ◦ dG P = dG P . Matrix Group Iterations – p. 17/20