MATH 612 Computational methods for equation solving and function - PowerPoint PPT Presentation

MATH 612 Computational methods for equation solving and function minimization – Week # 3 Instructor: Francisco-Javier ‘Pancho’ Sayas Spring 2014 – University of Delaware FJS MATH 612 1 / 28

Plan for this week Discuss any problems you couldn’t solve of Lectures 3 to 5 (Lectures are the chapters of the book) Read Lectures 6 and 7 Homework is due Wednesday at class time The first coding assignment will be posted at the end of the week Friday is a work’n’code day. You’ll be given some problems to work on and a small coding assignment. Warning Typically, I’ll keep on updating, correcting, and modifying the slides until the end of each week. FJS MATH 612 2 / 28

ORTHOGONAL PROJECTIONS FJS MATH 612 3 / 28

Orthogonal decomposition Consider an orthonormal basis of C m (or R m ) q 1 , . . . , q n , q n + 1 , . . . , q m . ( S 1 ⊥ S 2 ) � �� � � �� � basis of S 1 basis of S 2 Given a vector x the vector   n n � � ( q ∗  q j q ∗  x = Px j x ) q j = j j = 1 j = 1 is in S 1 and x − Px is in the orthogonal complement of S 1 . In other words, the decomposition x = Px + ( I − P ) x is orthogonal. FJS MATH 612 4 / 28

The orthogonal projector Starting again: given orthonormal vectors q 1 , . . . , q n in C m , the matrix n � q j q ∗ P = j j = 1 gives the component in � q 1 , . . . , q n � of the decomposition of x as a sum of a vector in this subspace plus a vector orthogonal to it. This matrix is called the orthogonal projector onto the subspace S 1 = � q 1 , . . . , q n � . FJS MATH 612 5 / 28

Easy properties n � q j q ∗ P = q 1 , . . . , q n orthonormal j j = 1 P is hermitian ( P ∗ = P ) P 2 = P If n = m , then P = I Px = x for all x ∈ range ( P ) = � q 1 , . . . , q n � Px = 0 if and only if x ⊥ range ( P ) The eigenvalues of P are 1 (multiplicity n ) and 0 (multiplicity n − m ) I − P is also an orthogonal projector (what is its range?) FJS MATH 612 6 / 28

The SVD of an orthogonal projector n � q j q ∗ P = q 1 , . . . , q n orthonormal j j = 1 If � Q is the matrix whose columns are q 1 , . . . , q n , then P = � Q � Q ∗ . Q ∗ � (Note that � Q = I n though.) If Q is unitary and contains q 1 , . . . , q n as its first columns, then � I � O P = Q Σ Q ∗ , Σ = O O FJS MATH 612 7 / 28

Orthogonal projectors without orthonormal bases Let S 1 = � a 1 , . . . , a n � be a subspace of C m , given as the span of n linearly independent vectors. The orthogonal projector onto S 1 is given by the matrix P = A ( A ∗ A ) − 1 A ∗ . First of all, P ∗ = P and P 2 = P . (This is easy to verify. Do it!) Next, note how ( Px ) ∗ (( I − P ) y ) x ∗ P ( I − P ) y = x ∗ P ( I − P ) y = x ∗ ( P − P 2 ) y = 0 . = This means that the decomposition x = Px + ( I − P ) x is orthogonal and therefore P is an orthogonal projector onto its own range. FJS MATH 612 8 / 28

Orthogonal projectors without orthonormal bases (2) We only need to check that the range of P is the range of A , that is, S 1 . Two things to prove: if x ∈ range ( P ) , then x = A ( A ∗ A ) − 1 A ∗ y ∈ range ( A ) . � �� � z if x ∈ range ( A ) , then x = Ay and Px = A ( A ∗ A ) − 1 A ∗ A y = Ay = x ∈ range ( P ) � �� � I FJS MATH 612 9 / 28

How to compute an orthogonal projection We will find better ways in the next chapters. Here’s the most direct method You do not construct the matrix P . If you only have a basis of the subspace, construct the matrix A whose columns are the elements of the basis. Then Multiply y = A ∗ x Solve the system ( A ∗ A ) z = y = A ∗ x Multiply Az The result is the projection of x onto the space spanned by the columns of A . FJS MATH 612 10 / 28

Oblique projectors Any square matrix P satisfying P 2 = P is called a projector (an oblique projector). Two important results on oblique projectors: If P 2 = P , then 1 ( I − P ) 2 = I − 2 P + P 2 = I − P , so I − P is a projector too. It is called the complementary projector If P is a projector, P is an orthogonal projector if and only if 2 P ∗ = P . And... if P is a projector, I − 2 P is invertible. (Why? Show that ( I − 2 P ) 2 = I )) FJS MATH 612 11 / 28

Oblique projectors (2) Consider an oblique projector P . Let S 1 = range ( P ) . By definition of projector x ∈ range ( P ) ⇐ ⇒ Px = x . (Prove this!) Also, x ∈ range ( I − P ) ⇐ ⇒ Px = 0 , and we denote S 2 = range ( I − P ) = null ( P ) . The projector P projects onto S 1 along S 2 . We then have a decomposition for every vector x = Px + ( I − P ) x . ���� � �� � ∈ S 1 ∈ S 2 Geometrically speaking, to know the projection P , you need to know its range S 1 and the range of its complementary projection S 2 . With orthogonal projections, hey are orthogonal to each other and only one of them is needed. FJS MATH 612 12 / 28

A particular case 1 If S 1 = � p � = � q � where q = � p � p , then the orthogonal projector onto S 1 is given by 1 P = qq ∗ = � p � 2 pp ∗ , its complementary projection is 1 I − P = I − qq ∗ = I − � p � 2 pp ∗ . Finally I − 2 P is an operator for symmetry with respect to the hyperplane S ⊥ 1 . Note that I − 2 P is a unitary matrix, because P 2 = P = P ∗ . FJS MATH 612 13 / 28

Reduced QR decompositions FJS MATH 612 14 / 28

Reminder Let q 1 , . . . , q j be an orthonormal collection of vectors. The orthogonal projection of v onto � q 1 , . . . , q j � is the vector j � ( q ∗ j v ) q j . i = 1 FJS MATH 612 15 / 28

Reduced QR decomposition For the moment being, let A be an m × n matrix with full rank by columns. (Therefore m ≥ n .) A reduced QR decomposition is a factorization A = � Q � R , where: � Q has orthonormal columns, that is, Q ∗ � � Q = I n � R is upper triangular (with non-zero diagonal elements – why?) in some cases, it is also required that r ii > 0 for all i . FJS MATH 612 16 / 28

What does the decomposition mean? This is the matrix form:       r 11 . . . r 1 n   . ...  a 1  =  q 1  . · · · a n · · · q n   . r nn Now with vectors: a 1 = r 11 q 1 a 2 = r 12 q 1 + r 22 q 2 . . . a n = r 1 n q 1 + . . . + r nn q n FJS MATH 612 17 / 28

What does the decomposition mean? (2) With vectors, everything in one expression j � a j = r 1 j q 1 + . . . + r jj q j = r ij q i , j = 1 , . . . , n . i = 1 Something equivalent but written in a funny way (thanks to r jj � = 0 for all j ) � j − 1 � � 1 a j − r ij q i = q j , j = 1 , . . . , n . r jj i = 1 Can you see how � a 1 , . . . , a j � = � q 1 , . . . , q j � j = 1 , . . . , n ? Finding a reduced QR decomposition is equivalent to finding orthonormal bases for the growing column subspaces. FJS MATH 612 18 / 28

Existence, uniqueness, computation We are going to focus our attention in the case of an m × 3 matrix with full rank by columns. We want to find { q 1 , q 2 , q 3 } orthonormal, and the elements of an upper triangular matris r ij with positive diagonal 1 , satisfying a 1 = r 11 q 1 , = r 12 q 1 + r 22 q 2 , a 2 a 3 = r 13 q 1 + r 23 q 2 + r 33 q 3 . Since � q 1 � 2 = 1 and r 11 > 0, there are no many options for the first equation q 1 = 1 r 11 = � a 1 � 2 , v 1 . r 11 1 Otherwise, there’s no uniqueness. We’ll see how Matlab always computes another QR decomposition FJS MATH 612 19 / 28

Existence, uniqueness, computation (2) What we have so far... a 1 = r 11 q 1 , a 2 = r 12 q 1 + r 22 q 2 , = r 13 q 1 + r 23 q 2 + r 33 q 3 . a 3 Looking at the second equation, we use that q ∗ 1 q 2 = 0 and q ∗ 1 q 1 = 1 to obtain r 12 = q ∗ 1 a 2 . Then a 2 − ( q ∗ r 22 q 2 = 1 a 2 ) q 1 = a 2 − ( orth.proj. of a 2 onto � q 1 � ) = v 2 leading to q 2 = 1 r 22 = � v 2 � 2 , v 2 . r 22 FJS MATH 612 20 / 28

Existence, uniqueness, computation (3) What we have so far... a 1 = r 11 q 1 , a 2 = r 12 q 1 + r 22 q 2 , a 3 = r 13 q 1 + r 23 q 2 + r 33 q 3 . With the same idea (use that q ∗ i q j = δ ij ) r 13 = q ∗ r 23 = q ∗ 1 a 3 , 2 a 3 , and a 3 − ( q ∗ 1 a 3 ) q 1 − ( q ∗ r 33 q 3 = 2 a 3 ) q 2 = a 3 − ( orth.proj. of a 3 onto � q 1 , q 2 � ) = v 3 , leading to q 3 = 1 r 33 = � v 3 � 2 , v 3 . r 33 FJS MATH 612 21 / 28

Existence, uniqueness, computation (4) We can continue with as many vectors as we have... j − 1 � a j = r ij q i + r jj q j . i = 1 We first compute the coefficients r ij = q ∗ i a j i = 1 , . . . , j − 1 (upper part of the j -th column of � R ), substract the result j − 1 � r jj q j = v j = a j − r ij q i = a j − ( orth.proj. of a j onto � q 1 , . . . , q j − 1 � ) i = 1 and finally choose r jj = � v j � 2 and q j = 1 r jj v j . FJS MATH 612 22 / 28

Sounds familiar? This is the algorithm in algebraic form: for j = 1 , . . . , n , compute j − 1 � r ij = q ∗ i a j ( i = 1 , . . . , j − 1 ) , v j = a j − r ij q i , i = 1 and then q j = 1 r jj = � v j � 2 , v j . r jj (This is the well-known Gram-Schmidt orthogonalization method.) FJS MATH 612 23 / 28