Scientific Computing Language Homework 3
Problem 1 (Least squares, 20 pts). Define the following data in Matlab:
1 2
t=(0:.5:10)’; y = tanh(t);
1. Fit the data to a cubic polynomial (using least squares) and plot the data together with the polynomial fit.
2. Fit the data to the function c1 + c2z + c3z2 + c4z3, where z = t2/(1 + t2). Plot the data together with the fit. What feature of z makes this fit much better than the original cubic?
Problem 2 (Householder, 10 pts). Let P be a Householder reflection matrix.
1. Find a vector u such that Pu = −u.
2. What algebraic condition is necessary and sufficient for a vector w to satisfy Pw = w? In n dimensions, how many such linear independent vectors are there?
Problem 3 (QR, 10 pts). Let A = QR be the thin QR factorization where Q ∈ Rm×n and R ∈ Rn×n. The matrix P = QQT has some interesting and important properties:
1. Show that P = AA+.
2. Prove that P2 = P. Moreover, any vector x can be written as x = u+v where u = Px and
v = (I − P )x. Prove that u and v are orthogonal.
Problem 4 (Condition number, 15 pts). Suppose A ∈ Rm×n has full column rank with m > n. Condition number is defined similar to that of square matrix by replacing the inverse with pseudo inverse:
κ2(A) = ∥A∥2 · ∥A+∥2.
Show that κ2(A) = σ1 where σ1 and σn are the largest and smallest singular values of A respectively.
σn
Use the definition to show that κ2(A) = κ2(R). Moreover, show that κ2(AT A) = κ2(A)2.
Problem 5 (Power method, 15 pts). Suppose A ∈ Rm×n is nonsingular and that Q ∈ Rn×p has orthonormal columns. The following iteration is referred to as inverse orthogonal iteration.
for k = 1,2,… do
Solve AZk = Qk−1 for Zk ∈ Rn×p; Apply QR factorization: Zk = QkRk;
end
Explain why this iteration can usually be used to compute the p smallest eigenvalues of A in absolute value. Note that to implement this iteration it is necessary to be able to solve linear systems that involve A. If p = 1, the method is referred to as the inverse power method.
1
Problem 6 (SVD, 30pts). The idea and data for this lab come from Stephens et al., “Dimensionality and Dynamics in the Behavior of C. elegans,” PLoS Comput Biol, 2008. In this work the authors captured video of worms moving as they were subjected to stimuli (from “standard” to “painful”). Using image processing, they found 100 points representing a path along the back of single worm within each frame of the video and computed a representation independent of rotation and translation using the tangents to the path. The result is a 100 × n matrix of tangent angles for n frames of video. They used n = 56200.
The authors then noted that dimension reduction by the SVD is extraordinarily successful for this data set; they showed that the data are well characterized by a small number of “eigenworms”.1 This makes some sense, as the motions are constrained a great deal by anatomy and kinematics.
Suppose T is the matrix of angles; i.e., column tj is the vector of angles in the jth video frame. SupposewehavethefullSVDT =USVT. ThiswouldmakeV ann×nmatrix,whichwouldnotfit in memory, so we have to use a thin SVD. In the text we only did this for an m×n matrix with m > n, which is not the case for T. However, TT = VSTU does have a thin form in which TT = VSTU, where V is only n × 100 and S is 100 × 100. Finally, this gives us
the thin SVD we need. Observe that
T =USVT, (1)
100
T =σkukvkT. (2) k=1
Because the singular values are always in decreasing order, we may approximate the original matrix by
r
Tr =σkukvkT, (3)
k=1
for some rank r ≪ 100. The range of this matrix is spanned by u1 , . . . , ur , which are the eigenworms.
A good way to express the proportion of T that is captured by Tr is the ratio
τr = ∥sr∥2 , (4)
∥s100 ∥2
wheresk isthevectorσ1,…,σk.
One use of the eigenworms is to create a compact representation of the data. A column tj of the
original data can be expressed in terms of its closest approximation as a linear combination of the eigenworms:
tj ≈c1u1 +···+crur =Urc,
where Ur is 100 × r. This is a least squares problem, but by orthogonality its solution is c = UrT tj . The r values in this vector give the components of tj in the principal eigenworm directions and could be used as a low-dimensional representation for further analysis.
Goals
Given the data matrix, you will perform the SVD analysis, find a reasonable value for the cutoff rank r using the coefficients τk, and extract the eigenworms.
1One often sees the “eigen-” prefix used in this context, but the SVD is a more fundamental description of the mathematics.
2
Procedure
1. Load the shapes.mat file from the assignment site. It has the matrix T.
2. On one graph, plot the first three columns of T . (These are tangent angles of the worm’s body
as a function of arc length.)
3. Using svd, compute the three matrices in the thin SVD (1).
4. Let s be the vector of singular values. Using it, plot 1 − τr versus r on a semi-log scale, for r = 1, . . . , 100. From this plot it should be clear that r = 4 is a compelling choice. Compute and print out the value of τ4.
5. In a 2-by-2 subplot grid, plot the first 4 eigenworms.
6. For the first three columns of {T}, plot the best approximation of each column by the leading 4 eigenworms. (The results will be much smoother curves than in step 1.)
3