Description
Tasks
1. Open an M-file called m551lab11.m and add the following commands
t=linspace(0,2*pi,100);
X=[cos(t);sin(t)];
figure(1);
plot(X(1,:),X(2,:),’b’);
axis equal
this will create a 2 × 100 matrix X = [x1 x2 · · · x100] whose columns are unit vectors pointing in various
directions. The plot will show a blue circle corresponding to the endpoints of these vectors.
Variables: X
2. Now in the M-file, define a variable A holding the matrix
A =
2 1
−1 1
and a variable AX holding the product A*X. Note that this product equals
AX = [Ax1 Ax2 · · · Ax100],
so the columns of AX show where the matrix A maps each of the columns of X. Add the following lines
to your M-file after the definition of AX.
1 Math 551
figure(2);
plot(AX(1,:),AX(2,:),’b’);
axis equal
Now when you run the M-file you should see a blue circle in Figure 1 and a blue ellipse in Figure 2. This
shows that A maps vectors that end on the unit circle to vectors ending on an ellipse.
Variables: A, AX
3. Add the following line to your M-file and re-run the file.
[U,S,V] = svd(A)
Your output should look like this
U =
-0.9571 0.2898
0.2898 0.9571
S =
2.3028 0
0 1.3028
V =
-0.9571 -0.2898
-0.2898 0.9571
The function svd computes the singular value decomposition of A:
A = USVT
where U = [u1 u2] and V = [v1 v2] are orthogonal matrices and
S =
σ1 0
0 σ2
is a non-negative diagonal matrix. Moreover
[Av1 Av2] = AV = USVT V = US = [u1 u2]
σ1 0
0 σ2
= [σ1u1 σ2u2],
showing that
Av1 = σ1u1 and Av2 = σ2v2.
You can check this fact in the command window as follows.
>> A – U*S*V’
ans =
1.0e-15 *
0.8882 0.4441
-0.4441 0.2220
>> U’*U
ans =
2
1.0000 0
0 1.0000
>> V’*V
ans =
1 0
0 1
Variables: U, S, V
4. Plot the columns of V on Figure 1 by adding the following to your M-file.
figure(1);
hold on;
quiver(0,0,V(1,1),V(2,1),0,’r’);
quiver(0,0,V(1,2),V(2,2),0,’g’);
hold off;
The vector corresponding to the first column will be plotted in red and the second column in green.
5. Plot the columns of US on Figure 2 by adding the following to your M-file.
figure(2);
hold on;
quiver(0,0,S(1,1)*U(1,1),S(1,1)*U(2,1),0,’r’);
quiver(0,0,S(2,2)*U(1,2),S(2,2)*U(2,2),0,’g’);
hold off;
The vector corresponding to the first column will be plotted in red and the second column in green.
6. Now we’ll work on a bigger matrix. Add the following commands to your M-file and run it. (If you want,
you can create a new cell by entering two percent signs (%%) at the beginning of a line and just run the
cell from now on.)
%% (start a new cell)
Img = double(imread(’TarantulaNebula.png’));
figure(3);
image(Img);
colormap(gray(256));
This code loads an image as a real matrix and displays it on the screen. Each entry in the matrix
corresponds to a pixel on the screen and takes a value somewhere between 0 (black) and 255 (white).
Variables: Img
7. Perform the singular value decomposition of Img, save the output in matrices UImg, SImg, and VImg.
Variables: UImg, SImg, VImg
8. Add the following code and run the M-file.
figure(4);
semilogy(diag(SImg));
This shows the singular values (the diagonal entries of the S matrix) for our image matrix (the scale of
the y axis is logarithmic). Notice that the diagonal entries of S are ordered so that σ1 ≥ σ2 ≥ σ3 ≥ · · · .
One way of writing the SVD is
A = σ1u1v
T
1 + σ2u2v
T
2 + · · · + σrurv
T
r
3 Math 551
where r is the rank of A and ui and vi are the ith columns of U and V respectively.
If for some k < r, σk+1 is very small compared to σ1, we should expect
A ≈ σ1u1v
T
1 + σ2u2v
T
2 + · · · + σkukv
T
k
to be a good approximation of A. (This is called a truncated SVD.) This idea can be used for image
compression as follows. Instead of storing the whole m × n matrix A, we instead store the m × k and
n × k matrices
C = [u1 u2 · · · uk] and D = [σ1v1 σ2v2 · · · σkvk].
If k(m + n) is much smaller than mn, then storing C and D will take much less space than storing A.
Moreover, if we wish to display the image, we can reconstruct A (approximately) as A ≈ CDT
.
9. Add the following code to your M-file.
%% (start a new cell)
k = 10;
% compressed image
C = UImg(:,1:k);
D = VImg(:,1:k)*SImg(1:k,1:k);
% compression percentage
pct = 1 – (numel(C)+numel(D))/numel(Img);
figure(5);
image(C*D’);
colormap(gray(256));
title( sprintf( ’%.2f%% compression’, 100*pct ) );
figure(6);
image(Img);
colormap(gray(256));
title( ’Original’ );
This code compresses the image as described above using k = 10 singular values, then displays the
reconstructed image along with the original. It also shows the percent decrease in storage size required
for the compressed image. Try several values of k to see how the quality and compression percentage vary
with k. Set k so that the compression percentage is 80.23% before turning in your M-file.
Variables: k, C, D
4 Math 551
