Dynamic Systems and Control, Assignment 2 -Solutions

MASSACHUSETTS INSTITUTE OF TECHNOLOGY Department of Electrical Engineering and Computer Science 6.241: Dynamic Systems—Fall 2007 Homework 2 Solutions Exercise 1.4 a) First define all the spaces: R(A) = {y ∈ Cm | ∃x ∈ Cn such that y = Ax} R⊥ (A) = {z ∈ Cm | y � z = z � y = 0, ∀y ∈ R(A)} R(A� ) = {p ∈ Cn | ∃v ∈ Cm such that p = A� v} N (A) = {x ∈ Cn | Ax = 0} N (A� ) = {q ∈ Cm | A� q = 0} i) Prove that R⊥ (A) = N (A� ). Proof: Let z ∈ R⊥ (A) → y � z = 0 ∀y ∈ R(A) → x� A� z = 0 ∀x ∈ Cn → A� z = 0 → z ∈ N (A� ) → R⊥ (A) ⊂ N (A� ). Now let q ∈ N (A� ) → A� q = 0 → x� A� q = 0 ∀x ∈ Cn → y � q = 0 ∀y ∈ R(A) → q ∈ R⊥ (A) → N (A� ) ⊂ R⊥ (A). Therefore R⊥ (A) = N (A� ). ii) Prove that N ⊥ (A) = R(A� ). Proof: From i) we know that N (A) = R⊥ (A� ) by switching A with A� . That implies that N ⊥ (A) = {R⊥ (A� )}⊥ = R(A� ). b) Show that rank(A) + rank(B) − n ≤ rank(AB) ≤ min{rank(A), rank(B)}. Proof: i) Show that rank(AB) ≤ min{rank(A), rank(B)}. It can be proved as follows: Each column of AB is a combination of the columns of A, which implies that R(AB) ⊆ R(A). Hence, dim(R(AB)) ≤ dim(R(A)), or equivalently, rank(AB) ≤ rank(A). Each row of AB is a combination of the rows of B → rowspace (AB) ⊆ rowspace (B), but the dimension of rowspace = dimension of column space = rank, so that rank(AB) ≤ rank(B). Therefore, 1 rank(AB) ≤ min{rank(A), rank(B)}. ii) Show that rank(A) + rank(B) − n ≤ rank(AB). Let rB = rank(B) rA = rank(A) where A ∈ Cm×n , B ∈ Cn×p . Now, let {v1 , · · · , vrB } be a basis set of R(B), and add n − rB linearly independent vectors {w1 , · · · , wn−rB } to this basis to span all of Cn , {v1 , v2 , · · · , vn , w1 , · · · , wn−rB }. Let M= � v1 | v2 · · · vrB | w1 | · · · wn−rB � = � V W � . Suppose x ∈ Cn , then x = M α for some α ∈ Cn . 1. R(A) = R(AM ) = R([AV |AW ]). Proof: i) Let x ∈ R(A). Then Ay = x for some y ∈ Cn . But y can be written as a linear combination of the basis vectors of Cn , so y = M α for some α ∈ Cn . Then, Ay = AM α = x → x ∈ R(AM ) → R(A) ⊂ R(AM ). ii) Let x ∈ R(AM ). Then AM y = x for some y ∈ Cn . But M y = z ∈ Cn → Az = x → x ∈ R(A) → R(AM ) ⊂ R(A). Therefore, R(A) = R(AM ) = R([AV |AW ]). 2. R(AB) = R(AV ). Proof: i) Let x ∈ R(AV ). Then AV y = x for some y ∈ CrB . Yet, V y = Bα for some α ∈ Cp since the columns of V and B span the same space. That implies that AV y = ABα = x → x ∈ R(AB) → R(AV ) ⊂ R(AB). ii) Let x ∈ R(AB). Then (AB)y = x for some y ∈ Cp . Yet, again By = V θ for some θ ∈ CrB → ABy = AV θ = x → x ∈ R(AV ) → R(AB) ⊂ R(AV ). Therefore, R(AV ) = R(AB). Using fact 1, we see that the number of linearly independent columns of A is less than or equal to the number of linearly independent columns of AV + the number of linearly independent columns of AW , which means that rank(A) ≤ rank(AV ) + rank(AW ). Using fact 2, we see that rank(AV ) = rank(AB) → rank(A) ≤ rank(AB) + rank(AW ), yet, there re only n − rB columns in AW . Thus, → rank(AW ) ≤ n − rB → rank(A) − rank(AB) ≤ rank(AW ) ≤ n − rB → rA − (n − rB ) ≤ rAB . 2 This completes the proof. Exercise 2.2 (a) For the 2nd order polynomial p2 (t) = a0 + a1 t + a2 t2 , we have f (ti ) = p2 (ti ) + ei i = 1, . . . , 16, and ti ∈ T. We can express the relationship between yi and the polynomial as follows; ⎤ ⎡ ⎤ ⎡ ⎤ 1 t1 t21 ⎤ ⎡ ⎡ y1 e 1 a ⎢ . . .. ⎥ 0 ⎢ .. ⎥ ⎢ .. .. . ⎥ . ⎥ ⎣ a1 ⎦ + ⎢ ⎥ ⎣ .. ⎦ ⎣ . ⎦=⎢ 2 ⎣ 1 t16 t16 ⎦ a2 e16 y16 The coefficients a0 , a1 , and a2 are determined by the least � �� squares solution to this (overconstrained) problem, a = (A� A)−1 A� y, where aLS = a0 a1 a2 . Numerically, the values of the coefficients are: ⎡ ⎤ 0.5296 aLS = ⎣ 0.2061 ⎦ 0.375 For the 15th order polynomial, by a similar reasoning we can express the relation between data points yi and the polynomial as follows: ⎤⎡ ⎤ ⎤ ⎡ ⎤ ⎡ ⎡ 1 t1 t21 · · · t115 a0 e1 y1 ⎥ ⎢ .. ⎥ ⎢ .. ⎥ ⎢ . ⎥ ⎢ . . .. .. ... ⎦ ⎣ . ⎦+⎣ . ⎦ ⎣ .. ⎦ = ⎣ .. .. . . 2 15 y16 a15 e16 1 t16 t16 · · · t16 This can be rewritten as y = Aa + e. Observe that matrix A is invertible for distinct t�i s. So � �� the coefficients ai of the polynomial are aexact = A−1 y, where aexact = a0 a1 · · · a15 . The resulting error in fitting the data is e = 0, thus we have a perfect fit at these particular time instants. Numerically, the values of the coefficients of are: ⎡ ⎤ 0.49999998876521 ⎢ 0.39999826604650 ⎥ ⎢ ⎥ ⎢ 0.16013119161635 ⎥ ⎢ ⎥ ⎢ 0.04457531385982 ⎥ ⎢ ⎥ ⎢ 0.00699544100513 ⎥ ⎢ ⎥ ⎢ −0.00976690595462 ⎥ ⎢ ⎥ ⎢ −0.02110628552919 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ 0.02986537283027 ⎥ aexact = ⎢ ⎥ ⎢ −0.03799813521505 ⎥ ⎢ ⎥ ⎢ 0.00337725219202 ⎥ ⎢ ⎥ ⎢ −0.00252507772183 ⎥ ⎢ ⎥ ⎢ 0.00072658523695 ⎥ ⎢ ⎥ ⎢ −0.00021752221402 ⎥ ⎥ ⎢ ⎢ −0.00009045014791 ⎥ ⎥ ⎢ ⎣ −0.00015170733465 ⎦ −0.00001343734075 3 Part A Dashed: 2nd−degree approximation ; Dash−dotted: 15th−degree approximation 2.5 2 1.5 1 0.5 0 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Figure 2.2a The function f (t) as well as the approximating polynomials p15 (t) and p2 (t) are plotted in Figure 2.2b. Note that while both polynomials are a good fit, the fifteenth order polynomial is a better approximation, as expected. b) Now we have measurements affected by some noise. The corrupted data is ỹi = f (ti ) + e(ti ) i = 1, . . . , 16 ti ∈ T where the noise e(ti ) is generated by a command “randn” in Matlab. Following the reasoning in part (a), we can express the relation between the noisy data points ỹi and the polynomial as follows: ỹ = Aa + ẽ The solution procedure is the same as in part (a), with y replaced by ỹ. Numerically, the values of the coefficients are: ⎡ 0.00001497214861 ⎢ 0.00089442543781 ⎢ ⎢ −0.01844588716755 ⎢ ⎢ 0.14764397515270 ⎢ ⎢ −0.63231582484352 ⎢ ⎢ 1.62190727992829 ⎢ ⎢ −2.61484909708492 ⎢ ⎢ ⎢ 2.67459894145774 aexact = ⎢ ⎢ −1.67594757924772 ⎢ ⎢ 0.56666848864500 ⎢ ⎢ −0.06211921500456 ⎢ ⎢ 0.00219622725954 ⎢ ⎢ −0.01911248745682 ⎢ ⎢ 0.01085690854235 ⎢ ⎣ −0.00207893294346 0.00010788458590 4 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ∗ 105 ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ Part B Dashed: 2nd−degree approximation ; Dash−dotted: 15th−degree approximation 7 6 5 4 3 2 1 0 −1 −2 −3 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Figure 2.2b and ⎡ aLS ⎤ 1.2239 = ⎣ −0.1089 ⎦ 0.3219 The function f (t) as well as the approximating polynomials p15 (t) and p2 (t) are plotted in Figure 2.2b. The second order polynomial does much better in this case as the fifteenth order polynomial ends up fitting the noise. Overfitting is a common problem encountered when trying to fit a finite data set corrupted by noise using a class of models that is too rich. Additional Comments A stochastic derivation shows that the “minimum variance unbiased estimator” for a is â = argmin�ỹ − Aa�2W where W = Rn−1 , and Rn is the covariance matrix of the random variable e. So, â = (A� W A)−1 A� W ỹ. Roughly speaking, this is saying that measurements with more noise are given less weight in the estimate of a. In our problem, Rn = I because the e�i s are independent, zero mean and have unit variance. That is, each of the measurments is “equally noisy” or treated as equally reliable. c) p2 (t) can be written as p2 (t) = a0 + a1 t + a2 t2 . In order to minimize the approximation error in least square sense, the optimal p̂2 (t) must be such that the error, f − pˆ2 , is orthogonal to the span of {1, t, t2 }: < f − pˆ2 , 1 >= 0 →< f, 1 >=< pˆ2 , 1 > < f − pˆ2 , t >= 0 →< f, t >=< pˆ2 , t > < f − pˆ2 , t2 >= 0 →< f, t2 >=< pˆ2 , t2 > . 5 Part C Dashed: 2nd−degree approximation 2.5 2 1.5 1 0.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Figure 2.2c We have that f = 12 e0.8t for t ∈ [0, 2], So, � 2 < f, 1 >= 0 2 � < f, t >= 0 < f, t2 >= � 0 1 0.8t 5 8 5 e dt = e 5 − 2 8 8 t 0.8t 15 8 25 e dt = e 5 + 2 32 32 2 2 t 2 e0.8t dt = 85 8 125 e5 − . 64 64 And, 8 < pˆ2 , 1 >= 2a0 + 2a1 + a2 3 8 < pˆ2 , t >= 2a0 + a1 + 4a2 3 8 32 < pˆ2 , t2 >= a0 + 4a1 + a2 3 5 Therefore the problem reduces to solving another set of linear equations: ⎡ ⎤⎡ ⎤ ⎡ ⎤ 2 2 38 a0 < f, 1 > ⎣ 2 8 4 ⎦ ⎣ a1 ⎦ = ⎣ < f, t > ⎦ . 3 8 32 4 a2 < f, t2 > 3 5 Numerically, the values of the coefficients are: ⎡ ⎤ 0.5353 a = ⎣ 0.2032 ⎦ 0.3727 The function f (t) and the approximating polynomial p2 (t) are plotted in Figure 2.2c. Here we use a different notion for the closeness of the approximating polynomial, p̂2 (t), to the original function, f . Roughly speaking, in parts (a) and (b), the optimal polynomial will be the one for 6 which there is smallest discrepancy between f (ti ) and p2 (ti ) for all ti , i.e., the polynomial that will come closest to passing through all the sample points, f (ti ). All that matters is the 16 sample points, f (ti ). In this part however, all the points of f matter. � Exercise 2.3 Let y = y1 y2 � � � � � � � C1 e1 S1 0 ,A= ,e= and S = . C2 e2 0 S2 Note that A has full column rank because C1 has full column rank. Also note that S is symmetric positive definite since both S1 and S2 are symmetric positive definite. Therefore, we know that x̂ = argmin e� Se exists and is unique and is given by x̂ = (A� SA)−1 A� Sy. Thus by direct substitution of terms, we have: x̂ = (C1� S1 C1 + C2� S2 C2 )−1 (C1� S1 y1 + C2� S2 y2 ) Recall that x̂1 = (C1� S1 C1 )−1 C1� S1 y1 and that x̂2 = (C2� S2 C2 )−1 C2� S2 y2 . Hence x̂ can be re-written as: x̂ = (Q1 + Q2 )−1 (Q1 x̂1 + Q2 x̂2 ) Exercise 2.8 We can think of the two data sets as sequentially available data sets. x̂ is the least squares solution to y ≈ Ax �corresponding � � �to minimizing the euclidean norm of e1 = y − Ax. y A x̄ is the least squares solution to ≈ x corresponding to minimizing e�1 e1 + e�2 Se2 where z D e2 = z − Dx and S is a symmetric (hermitian) positive definite matrix of weights. By the recursion formula, we have: x̄ = x̂ + (A� A + D� SD)−1 D� S(z − Dx̂) This can be re-written as: x̄ = x̂ + (I + (A� A)−1 D� SD)−1 (A� A)−1 D� S(z − Dx̂) = x̂ + (A� A)−1 D� (I + SD(A� A)−1 D� )−1 S(z − Dx̂) This step follows from the result in Problem 1.3 (b). Hence x̄ = x̂ + (A� A)−1 D� (SS −1 + SD(A� A)−1 D� )−1 S(z − Dx̂) = x̂ + (A� A)−1 D� (S −1 + D(A� A)−1 D� )−1 S −1 S(z − Dx̂) = x̂ + (A� A)−1 D� (S −1 + D(A� A)−1 D� )−1 (z − Dx̂) In order to ensure that the constraint z = Dx is satisfied exactly, we need to penalize the corresponding error term heavily (S → ∞). Since D has full row rank, we know there exists at least one value of x that satisfies equation z = Dx exactly. Hence the optimization problem we are setting up does indeed have a solution. Taking the limiting case as S → ∞, hence as S −1 → 0, we get the desired expression: x̄ = x̂ + (A� A)−1 D� (D(A� A)−1 D� )−1 (z − Dx̂) 7 In the ’trivial’ case where D is a square (hence non-singular) matrix, the set of values of x over which we seek to minimize the cost function consists of a single element, D−1 z. Thus, x̄ in this case is simply x̄ = D−1 z. It is easy to verify that the expression we obtained does in fact reduce to this when D is invertible. Exercise 3.1 The first and the third facts given in the problem are the keys to solve this problem, in addition to the fact that: � � R UA = . 0 Here note that R is a nonsingular, upper-triangular matrix so that it can be inverted. Now the problem reduces to show that x̂ = arg min �y − Ax�22 = arg min(y − Ax)� (y − Ax) x x is indeed equal to x̂ = R−1 y1 . Let’s transform the problem into the familiar form. We introduce an error e such that y = Ax + e, and we would like to minimize �e�2 which is equivalent to minimizing �y−Ax�2 . Using the property of an orthogonal matrix, we have that �e�2 = �U e�2 . Thus with e = y − Ax, we have �e�22 = �U e�22 = e� U � U e = (U (y − Ax))� (U (y − Ax)) = �U y − U Ax�22 �� � � � �2 � y1 � R � � � = � − x� � = (y1 − Rx) (y1 − Rx) + y2 y2 . y2 0 2 Since �y2 �22 = y2� y2 is just a constant, it does not play any role in this minimization. Thus we would lik to have y1 − Rx̂ = 0 and because R is an invertible matrix, x̂ = R−1 y1 . 8