orthogonalization L. Olson October 27, 2015 Department of Computer - PowerPoint PPT Presentation

orthogonalization L. Olson October 27, 2015 Department of Computer Science University of Illinois at Urbana-Champaign 1

objectives • Revisit SVD and Orthogonal Matrices • Create orthogonal vectors • Outline the Gram-Schmidt algorithm for orthogonalization 2

normal equations: conditioning The normal equations tend to worsen the condition of the matrix. Theorem cond ( A T A ) = ( cond ( A )) 2 1 A = np.random.rand(10,10) 2 print(np.linalg.cond(A)) 3 print(np.linalg.cond(A.T.dot(A))) 4 50.0972712517 5 2509.73658686 6 3

other approaches • QR factorization. • For A ∈ R m × n , factor A = QR where • Q is an m × m orthogonal matrix • R is an m × n upper triangular matrix (since R is an m × n upper � � R ′ triangular matrix we can write R = where R is n × n upper 0 triangular and 0 is the ( m − n ) × n matrix of zeros) • SVD - singular value decomposition • For A ∈ R m × n , factor A = USV T where • U is an m × m orthogonal matrix • V is an n × n orthogonal matrix • S is an m × n diagonal matrix whose elements are the singular values. 4

orthogonal matrices Definition A matrix Q is orthogonal if Q T Q = QQ T = I Orthogonal matrices preserve the Euclidean norm of any vector v , || Qv || 2 2 = ( Qv ) T ( Qv ) = v T Q T Qv = v T v = || v || 2 2 . 5

gram-schmidt orthogonalization One way to obtain the QR factorization of a matrix A is by Gram-Schmidt orthogonalization. We are looking for a set of orthogonal vectors q that span the range of A . For the simple case of 2 vectors { a 1 , a 2 } , first normalize a 1 and obtain a 1 q 1 = || a 1 || . Now we need q 2 such that q T 1 q 2 = 0 and q 2 = a 2 + cq 1 . That is, R ( q 1 , q 2 ) = R ( a 1 , a 2 ) Enforcing orthogonality gives: q T 1 q 2 = 0 = q T 1 a 2 + cq T 1 q 1 6

gram-schmidt orthogonalization q T 1 q 2 = 0 = q T 1 a 2 + cq T 1 q 1 Solving for the constant c. c = − q T 1 a 2 q T 1 q 1 reformulating q 2 gives. q 2 = a 2 − q T 1 a 2 q 1 q T 1 q 1 Adding another vector a 3 and we have for q 3 , q 3 = a 3 − q T q 2 − q T 2 a 3 1 a 3 q 1 q T q T 2 q 2 1 q 1 Repeating this idea for n columns gives us Gram-Schmidt orthogonalization. 7

gram-schmidt orthogonalization Since R is upper triangular and A = QR we have a 1 q 1 r 11 = a 2 = q 1 r 12 + q 2 r 22 . . = . . . . a n = q 1 r 1 n + q 2 r 2 n + ... + q n r nn From this we see that r ij = q T i a j i q i , j > i q T 8

orthogonal projection The orthogonal projector onto the range of q 1 can be written: q 1 q T 1 q T 1 q 1 . Application of this operator to a vector a orthogonally projects a onto q 1 . If we subtract the result from a we are left with a vector that is orthogonal to q 1 . 1 ( I − q 1 q T q T 1 ) a = 0 q T 1 q 1 9

gram-schmidt orthogonalization 1 def qr(A): 2 Q = np.zeros(A.shape) 3 4 for k in range(A.shape[1]): 5 avec = A[:, k] 6 q = avec 7 for j in range(k): 8 q = q - np.dot(avec, Q[:,j])*Q[:,j] 9 10