APPM 4720/5720 Take-Home Midterm Solutions This exam is due

APPM 4720/5720
Take-Home Midterm Solutions
This exam is due at the start of class on Friday, March 14th. Minimal credit will be given for incomplete
solutions or solutions that do not provide details on how the solution is found. Problems marked with a ?
are required only for students enrolled in APPM 5720. You may consult textbooks and other sources but
collaboration with other people (e.g. classmates) is expressly prohibited.
1. Suppose that H = A + iB is Hermitian and positive definite with A, B ∈ Rn×n , i.e. H ∗ = H and
z∗ Hz > 0 whenever z 6= 0 for z ∈ Cn . Here H ∗ indicates the complex conjugate transpose of H. The
real matrix
C=
A
B
−B
A
is called the equivalent real form of H. Equivalent real forms are often used to perform computations
on complex systems without defined complex data types.
(a) Let z = x + iy ∈ Cn for x, y ∈ Rn and H be as described above. Show how you can use C to
compute Hz using only real matrix arithmetic.
Solution:
Notice that in complex arithmetic we have
Hz = (A + iB) (x + iy) = Ax − By + i (Ay + Bx)
Now, if we stack the real and imaginary components of z in one vector, and multiply by C, we
have
−B
A
A
B
x
y
=
Ax − By
Ay + Bx
So, the first half of the entries in the vector correspond to the real part of the result, and the
second half correspond to the imaginary part.
(b) Prove that C is symmetric positive definite.
Solutions: First let’s deal with symmetry. H hermitian means that H ∗ = H. Then we have
H∗ = H
⇒
∗
(A + iB) = (A + iB)
⇒
AT − iB T = A + iB
Equating the real and imaginary parts we have that
A = AT
and B T = −B.
In other words, A is symmetric and B is skew-symmetric. Using this, we have
T
C =
AT
−B T
BT
AT
A
B
−B
A
=C
so C is symmetric. To show that H positive definite implies C positive definite we need to prove
a little lemma about skew-symmetric matrices.
Lemma: Let B ∈ Rn×n be skew-symmetric. Then xT Bx = 0 for all x.
Proof: Since xT Bx is just a scalar we must have that
xT Bx = xT Bx
T
= xT B T x = −xT Bx
Since we’re working with real numbers here, the only way for the above to be true is if xT Bx = 0
Now, let’s see what H positive definite means in terms of A and B. Let z = (x + iy). Then
z∗ Hx
=
=
∗
(x + iy) (A + iB) (x + iy)
xT − iyT (A + iB) (x + iy)
!
0 T *0
* T = x Ax − x By + y Ax + y By + i x Ay + x Bx − y
Ax + y
By
T
T
T
T
*
T
:
T
= xT Ax − xT By + yT Ax + yT By > 0
In the third step the second and fourth terms in the imaginary part are 0 by our lemma, and the
first and third terms cancel out because A is symmetric. The inequality in the last step follows
from the assumption that H is positive definite.
Now let’s investigate the positive definiteness of C. For arbitrary nonzero vectors x and y ∈ Rn
we have
x
T
y
T
A −B
B
A
x
y
=
T
x
y
T
Ax − By
Bx + Ay
= xT Ax − xT By + yT Bx + yT Ay > 0
where here the inequality follows from the positive definiteness of H. Thus, C is positive definite.
(c) Formulate an algorithm for solving (A + iB)(x + iy) = (b + ic), where b, c, x, y are in Rn . It
should involve 8n3 /3 flops. How much storage is required?
Solution:
From the analysis above, it is clear that we can solve the complex system (A + iB) (x + iy) =
(b + ic) by solving the equivalent real form version
A
B
−B
A
x
y
=
b
c
This is an SPD linear system of size 2n × 2n which we can solve by performing Cholesky Factorization of C and then doing a forward and a backward solve. For an n × n SPD matrix the Cholesky
Factor can be computed in approximately n3 /3 flops. Since here we have a 2n × 2n system the
3
computation takes (2n) /3 = 8n3 /3 flops. The forward
and back solves cost approximately 8n2
3
flops, which does not contribute to the leading O n term in the asymptotic cost.
Since C is symmetric we only need to need to store the upper triangular part of C which uses
2n2 memory. We also need 4n storage for the vector of unknowns and the right-hand side vector.
Since efficient versions of the Cholesky algorithm overwrite C with R, this does not require any
additional storage. Note that it might look like we could save some storage by only storing A
once, but since we will eventually need that storage for R this doesn’t really help.
2. Let x and y in Rn be such that kxk2 = kyk2 .
(a) Show how to construct an orthogonal matrix Q using reflectors such that QT x = y.
Solution: Consider the following picture showing x and y with equal magnitude
x
y
Since kxk2 = kyk2 we should be able to find an orthogonal matrix Q that reflects x onto y by
QT x = y. Recall the derivation of the QR decomposition using reflectors from class. In that
case we wanted to build reflector Q that reflected x onto the x-axis by QT x = kxk2 e1 . The only
difference now is that we’re reflecting onto y. Consider the following diagram
x
u=y−x
y
Recalling the definition of the reflector, it is clear that Q takes on the form
Q = I − γuuT
where
u=y−x
and γ =
2
kuk22
(b) Repeat part (a) but with rotators instead of reflectors.
Solution: There are several solutions to this problem. The first approach is to use the same
orthogonal matrices based on rotators that we used in the QR algorithm that rotate x and y onto
the first canonical basis vector, and then combine them in a clever way to get the desired result.
Consider the following 2D diagram:
e2
x
y
e1
Let Qy be the orthogonal matrix based on rotators, computed in the usual way, such that QTy
maps y onto e1 . That is, QTy y = kyk2 e1 . Note that since Qy is an orthogonal matrix we have
Qy QTy y = y. In other words, QTy rotates y onto the first canonical basis vector, and Qy performs
the opposite rotation. Similarly, let Qx be the orthogonal matrix that rotates x onto e1 , that is,
QTx x = kxk2 e1 .
Then, the orthogonal matrix Q that rotates x onto y can be constructed by applying the rotation
QTx to rotate x onto e1 , and then applying the reverse rotation Qy to rotate the result onto y.
The process in the 2D case looks like
e2
e2
x
Qy QTx x = y
e1
e1
kxk2 e1 = QTx x
QTx x
The orthogonal matrix that takes x onto y via QT x = y is then
Q = Qy QTx
T
= Qx QTy
If we do not explicitly form Q then this solution is approximately twice as expensive as rotating
a vector onto the first canonical basis vector as seen in the QR algorithm with rotators.
The second approach is to construct the rotators component-wise as in the usual QR algorithm
based on rotators, but instead of rotating x onto e1 we perform plane-rotations and that rotate
the components of x directly onto y. The development of the algorithm is as follows:
Suppose we want to find a sequence of plane rotators that uses the first entry in x and y to rotate
the k th component of x onto the k th component of y. That is, we want to pick c and s such that
c
−s
s
c
x1
xk
r1
yk
=
T
Note that we might be tempted to try to rotate x1 xk
directly onto y1 ykT . Of course
we cannot accomplish this with an orthogonal matrix in general because even though x and y
T
have the same length, it is highly unlikely that this is true of the sub-vectors x1 xk
and
T
y1 yk
. The value r1 is chosen such that the length of the rotated sub-vector stays the same.
In this case we take
r1 = sign(y1 )
q
x21 + x2k − yk2
Note that the sign term is needed to ensure that in the last step r1 = y1 is satisfied automatically.
The values of c and s must be chosen such that
cx1 + sxk
=
r1
−sx1 + cxk
=
yk
2
=
1
2
c +s
Solving this system for c and s yields
c=
xk yk + x1 r1
x21 + x2k
and s =
xk r1 − x1 yk
x21 + x2k
We then proceed exactly as in the case of QR with rotators. Pseudocode for the proposed algorithm
is given below which computes QT x and also constructs the orthogonal matrix Q.
Initialize Q to In
for k = 2, . . . , n do
p
r1 = sign(y1 ) x21 + x2k − yk2
c = (xk yk + x1 r1 ) / x21 + x2k
s = (xk r1 − x1 yk ) / x21 + x2k
x1 = cx1 + sxk (= r1 )
xk = −sx1 + cxk (= yk )
c s
Q1:n,[1 k] = Q1:n,[1 k]
−s c
end
where here Q1:n,[1 k] denotes the tall skinny matrix created by extracting the 1st and k th columns
of Q.
There is of course an inconvenient flaw in the algorithm as we have described it here. It is possible
that yk might be large enough that r1 will be imaginary. In practice, if you have access to complex
data types, then the algorithm as it’s presented will achieve the desired result, but the added cost
of using complex data types is significant. We can fix the algorithm so that imaginary numbers
are not encountered by being careful with the order in which we alter elements of the vector x.
There is no reason that we need to use x1 as the pivot in each operation. Instead it would be
better to pick the largest entry in x as the pivot. This makes it less likely that a negative number
will be encountered under the square root. We can also choose to loop over the entries in y such
a way that we encounter smaller yk ’s first. This also helps ensure that the quantity under the
square root is nonnegative. An algorithm that works in entirely real arithmetic is given below
Initialize Q to In
Set kmax to smallest k s.t. |xk | = max {|x1 | , |x2 | , . . . , |xn |}
Set ind to vector of indices k = 1, . . . , n sorted by size of |yk | and excluding kmax
for k in ind do
rkmax = sign(ykmax )
p
x2kmax + x2k − yk2
c = (xk yk + xkmax rkmax ) / x2kmax + x2k
s = (xk rkmax − xkmax yk ) / x2kmax + x2k
xkmax = cxkmax + sxk (= rkmax )
xk = −sxkmax + cxk (= yk )
Q1:n,[kmax k] = Q1:n,[kmax k]
c
−s
s
c
end
Note that due to the additional flops in computing c and s this algorithm is again about twice as
expensive as the usual QR algorithm based on rotators and the same cost as alternate solution to
this problem. There is a storage savings of about a factor of 2 compared to the alternate algorithm
because we only need to store one set of rotators for each entry instead of 2.
3. In exact arithmetic it is well known that the matrix A ∈ Rn×n with full rank has both a right inverse
and a left inverse, and they are the same, i.e. there exists B ∈ Rn×n such that
AB = BA = I
However, numerically the left inverse (i.e the matrix BL that satisfies BL A = I) is not necessarily a
good right inverse (i.e. the matrix BR that satisfies ABR = I) and vice versa.
(a) Let A = U ΣV T be the Singular Value Decomposition of A where U = [u1 , . . . , un ] and
V = [v1 , . . . , vn ] are orthogonal, and Σ = diag(σ1 , . . . , σn ) with σ1 ≥ . . . ≥ σn > 0. Show that A
is invertible, and give the explicit form of A−1 .
Solution: Since A is square and all of it’s singular values are nonnegative, A is invertible. To find
the explicit form of A−1 we note that
V Σ−1 U T
U ΣV T
V Σ−1 U T = U ΣV T V Σ−1 U T = U ΣΣ−1 U T = U U T = I
U ΣV T = V Σ−1 U T U ΣV T = V Σ−1 ΣV T = V V T = I
So A−1 = V Σ−1 U T .
(b) Let X = A−1 + vn uT1 . Express AX and XA as rank-one updates of the identity matrix.
Solution: We have
AX
= A A−1 + vn uT1
= I + U ΣV T vn uT1
= I + U Σen uT1
= I + σn U en uT1
= I + σn un uT1
and similarly
XA
A−1 + vn uT1 A
=
= I + vn uT1 U ΣV T
= I + vn eT1 ΣV T
= I + σ1 vn eT1 V T
= I + σ1 vn v1T
(c) Prove that if x and y are vectors in Rn , then kxyT k2 = kxk2 kyk2 . Use this to compute kAX − Ik2
and kXA − Ik2 . What can you say about the accuracy of X as a left and right inverse when κ2(A)
is large?
Solution: We have
kxyT k2 = max kxyT zk2
kzk=1
Since yT z is a scalar it is clear that the product will be maximized when yT z is as large as possible.
This is achieved when z = y/kyk2 . We then have
kxyT k2
=
max kxyT zk2
kzk=1
= kxyT yk2 /kyk2
= kyk22 kxk2 /kyk2
= kxk2 kyk2
Then
kAX − Ik2
= kσn un uT1 k2
= σn kun uT1 k2
= σn kun k2 ku1 k2
= σn
and similarly
kXA − Ik2
= kσ1 vn v1T k2
= σ1 kvn v1T k2
= σ1 kvn k2 kv1 k2
= σ1
Since κ2(A) = σ1 /σn , A ill-conditioned implies that σ1 σn . From the expression above we see
that if A is ill-conditioned then X is a much better right inverse than it is a left inverse since σ1
is much larger than σn .
4. ? We discussed the solution of the least squares problem at great length in class. Let A ∈ Rm×n and
C ∈ Rn×p with m ≥ n ≥ p and consider instead the constrained least squares problem
min kb − Axk2
s.t.
CT x = d
CT
A
C1T
A1
C2T
A2
x
where A and C have full rank. Let
=
where C1 ∈ Rp×p has rank p.
u
(a) Show that if x =
, then
v
b−A
u
v
ˆ − A2 − A1 C −T C T v,
=b
2
1
ˆ = b − A1 C −T d.
where b
1
Solution:
T
C x=
C1T
C2T
u
v
=d
⇒
C1T u + C2T v = d
u = C1−T d − C2T v
⇒
where here we know that C1−T exists because C1 is p × p and full rank. Then
b−A
u
v
C1−T d − C2T v
b−A
v
−T
b − A1 C1
d − C2T v − A2 v
=
=
=
=
b − A1 C1−T d − A2 − A1 C1−T C2T v
ˆ − A2 − A1 C −T C T v
b
2
1
(b) Write down the decomposition
Elim T in block form that results from applying p stepsof Gaussian
C
CT
ination to the matrix
. To make life easier, you may assume that
is column
A
A
diagonally dominant so you don’t have to worry about pivoting.
T C
Solution: Applying p steps of Gaussian Elimination without pivoting to
yields the folA
lowing partial LU Decomposition:
CT
A
=
C1T
A1
C2T
A2
=
L11
L22
0
I
U11
0
U12
U22
where L11 is unit lower triangular, U11 is upper triangular, and L22 , U12 and U22 are dense.
Furthermore, since C1 is nonsingular, so are L11 and U11 .
(c) Argue that if x satisfies the constraint C T x = d and minimizes the following
T d
C
min −
x
.
b
A
x
2
then x is a solution to the original constrained least squares problem.
Solution: Squaring the term we’re trying to minimize we have
T 2
d
C
kd − C T xk22 + kb − Axk22 .
min −
x
= min
b
A
x
x
2
Assuming x satisfies the constraint C T x = d the first term in the functional is 0 and we have
T 2
d
C
kb − Axk22
−
x
min = min
b
A
x
x
2
subject to C T x = d
which is the original constrained minimization problem.
(d) Show that the solution to the constrained least squares problem can be found by solving an
unconstrained least squares problem of the form
ˆ − Xvk2
min kb
v
ˆ and X depend on the blocks of the factorization you defined in part (b).
where b
Solution: From parts (d) and (a) we have
T 2
d
C
min −
x
b
A
x
2
=
min kb − Axk22
=
ˆ − A2 − A1 C −T C T vk2
min kb
2
2
1
x
v
Let’s see if we can simplify the matrix above. By equating blocks in the partial LU decomposition
in (b) we have
A1
= L22 U11
A2
= L22 U12 + U22
C1T
= L11 U11
C2T
= L11 U12
Then
A2 − A1 C1−T C2T
= L22 U12 + U22 − L22 U11 U11 −1 L11 −1 L11 U12
= L22 U12 + U22 − L22 U12
= U22
So the unconstrained minimization problem becomes
ˆ − U22 vk2 .
min kb
v
ˆ which is given by
We of course need to compute b
ˆ
b
= b − A1 C1−T d
= b − L22 U11 U11 −1 L11 −1 d
= b − L22 L11 −1 d.
ˆ is computed we can solve the unconstrained minimization problem via our favorite least
Once b
squares solver, such as QR. Then, to get u we have
u = C1−T d − C2T v
u = U11 −1 L11 −1 d − C2T v
which can be solved using one step of forward and back substitution.
(e) Write down a step-by-step algorithm, based on the steps above, for solving the constrained least
squares problem. I want to see a general outline. There is no need to give pseudocode.
i. Do p steps of Gaussian elimination to obtain
T T
L11
C
C1 C2T
=
=
L22
A
A1 A2
0
I
U11
0
U12
U22
ˆ = b − L22 L11 −1 d using forward substitution to apply the inverse.
ii. Compute b
iii. Solve the unconstrained minimization problem for v.
iv. Solve for u using a forward and a back substitution.

Download Report