In this laboratory session we will learn how to1. Find the SVD decomposition of a matrix using MATLAB2. Use the SVD to perform Image Compression.IntroductionLet A be an m n matrix (m n). Then A can be factored asA = USV T =. . .u1 u2 . . . um. . .mm1 0 . . . 00 2 . . . 00 0. . . 00 0 . . . n. . . 00 0 . . . 0mn. . .v1 v2 . . . vn. . .Tnn(L8.1)By convention 1 2 . . . n 0. U and V are orthogonal matrices (square matrices withorthonormal column vectors).Note that the matrix S is usually denoted by , however we will use S to be consistent with MATLABnotation.One useful application of the SVD is the approximation of a large matrix by another matrix of lowerrank. The lower-rank matrix contains less data and so can be represented more compactly than theoriginal one. This is the idea behind one form of signal compression.The Frobenius Norm and Lower-Rank ApproximationThe Frobenius norm is one way to measure the size of a matrix
a bc d
F=
a2 + b2 + c2 + d2and in general, AF =
i,j a2ij1/2, that is, the square root of the sum of squares of the entries of A.If A is an approximation to A, then we can quantify the goodness of fit by A AF. This is just aleast-squares measure.The SVD has the property that, if you want to approximate a matrix A by a matrix A of lower rank,then the matrix that minimizes A AF among all rank-1 matrices is the matrix A = u1 m1111v1T1n= 1u1vT1In general, the best rank-r approximation to A is given by A = . . .u1 u2 . . . ur. . .mr1 0 . . . 00 2 . . . 00 0. . . 00 0 . . . rrr. . .v1 v2 . . . vr. . .Trn= 1u1vT1 + 2u2vT2 + . . . + rurvTrc2011 Stefania Tracogna, SoMSS, ASU 1MATLAB sessions: Laboratory 8and it can be shown thatA AF =2r+1 + 2r+2 + . . . + 2n(L8.2)The MATLAB command[U, S, V] = svd(A)returns the SVD decomposition of the matrix A, that is, it returns matrices U, S and V such thatA = USV T .ExampleThe SVD of the following matrix A is:A =2 8 2014 19 102 2 1 =3545 04535 00 0 130 0 00 15 00 0 3132323231323232313(a) Enter the matrix A and compute U, S and V using the svd comand.Verify that A = USV T .Answer:>> A=[-2,8,20;14,19,10;2, -2, 1];>> [U,S,V]=svd(A)U =0.6000 -0.8000 0.00000.8000 0.6000 0.00000.0000 -0.0000 1.0000S =30.0000 0 00 15.0000 00 0 3.0000V =0.3333 0.6667 0.66670.6667 0.3333 -0.66670.6667 -0.6667 0.3333We can easily verify that U*S*V returns the matrix A.(b) Use MATLAB to find the best rank-1 approximation to A (with respect to the Frobenius norm),that is, A1 = 1u1vT1 .Use the command rank to verify that the matrix A1 has indeed rank one.Evaluate A A1F by typing norm(A A1,fro) and verify Eq. (L8.2).Answer:>> A1 = S(1,1)*U(:,1)*V(:,1)A1 =6.0000 12.0000 12.00008.0000 16.0000 16.00000.0000 0.0000 0.0000>> rank(A1)ans =1c2011 Stefania Tracogna, SoMSS, ASU 2MATLAB sessions: Laboratory 8>> norm(A-A1,fro)ans =15.2971>> sqrt(S(2,2)^2+S(3,3)^2)ans =15.2971Note that, although A1 is 33 and therefore it contains nine entries, we only need seven values togenerate it: (one singular value)+ (one column of U) +( one column of V ) = 1 + 3 + 3 = 7.This is less than the number of entries in the matrix A.(c) Find the best rank-2 approximation to A. Check the rank and the Frobenius norm of AA2 usingMATLAB and verify (L8.2).Answer:>> A2 = A1+S(2,2)*U(:,2)*V(:,2)A2 =-2.0000 8.0000 20.000014.0000 19.0000 10.0000-0.0000 0.0000 0.0000>> rank(A2)ans =2>> norm(A-A2,fro)ans =3.0000>> S(3,3)ans =3Note that we need 14 entries to generate A2:(2 singular values) + (2 columns of U) + (2 columns of V ) = 2 + 2 3 + 2 3 = 14.This is more than the number of entries in the original matrix A.Image Compression ExercisesInstructions: The following problems can be done interactively or by writing the commands in an Mfile(or by a combination of the two). In either case, record all MATLAB input commands and output ina text document and edit it according to the instructions of LAB 1 and LAB 2. For problem 2, includea picture of the rank-1 approximation. For problem 3, include a picture of the rank-10 approximationand for problem 4, include a picture of the lower rank approximation that, in your opinion, gives anacceptable approximation to the original picture. Crop and resize the pictures so that they do not takeup too much space and paste them into your lab report in the appropriate order.Step 1. Download the file cauchybw.jpg and save it to your working MATLAB directory. Thenload it into MATLAB with the commandA = imread(cauchybw.jpg); note the semicolonThe semicolon is necessary so that MATLAB does not print out many screenfuls of data. Theresult is a matrix of grayscale values corresponding to a black and white picture of a dog. (Thematrix has 104, 780 entries). Now, A is actually 3103383. To create an image from the matrixA, MATLAB uses the three values A(i, j, 1 : 3) as the RGB values to use to color the ijth pixel.We have a black and white picture so A(:,:,1) = A(:,:,2) = A(:,:,3) and we only need towork with A(:,:,1).c2011 Stefania Tracogna, SoMSS, ASU 3MATLAB sessions: Laboratory 8Step 2. We need to do some initial processing. TypeB = double(A(:,:,1)) + 1; dont forget the semicolonwhich converts A into the double-precision format that is needed for the singular value decomposition.Now typeB = B/256; semicolon![U S V] = svd(B); semicolon!This decomposition is just Eq. (L8.1).The gray scale goes from 0 to 256 in a black- and-white JPEG image. We divide B by 256 toobtain values between 0 and 1, which is required for MATLABs image routines, which we will uselater.PROBLEM 1. What are the dimensions of U, S, and V ? (Find out by typing size(U) withoutthe semicolon and likewise for the others.)Here S has more columns than rows; in fact, columns 311 to 338 are all zeros (When A has morecolumns than rows, we pad S on the right with zero columns to turn S into an m n matrix ).Otherwise, with this modification, the SVD is just like Eq. (L8.1).PROBLEM 2. Compute the best rank-1 approximation to B and store it in the matrix rank1(Use the commands given in the Example).Step 3. Lets visualize rank1. To do that, first createC = zeros(size(A)); semicolon!This creates an array of zeros, C, of the same dimension as the original image matrix A.Step 4. Copy the rank-1 image into C as follows:C(:,:,1) = rank1;C(:,:,2) = rank1;C(:,:,3) = rank1;Step 5: We are almost done, except for one hitch. MATLAB does all its scaling using valuesfrom 0 to 1 (and maps them to values between 0 and 256 for the graphics hardware). Lower-rankapproximations to the actual image can have values that are slightly less than 0 and greater than1. So we will truncate them to fit, as follows:C = max(0,min(1,C));Step 6. View the resulting image:image(C) no semicolonPROBLEM 3. Create and view a rank-10 approximation to the original picture (Use Steps 4-6but with rank10 instead of rank1. If you mess up for example, you get an all-black picture then start over from Step 3.) It is convenient to create an M-file with a for loop to evaluate1u1vT1 + . . . + 10u10vT10PROBLEM 4. Repeat with rank-20, 30 and 40 approximations (and any other ranks that youdlike to experiment with). What is the smallest rank that, in your opinion, gives an acceptableapproximation to the original picture?In your lab write-up, include the code that gives an acceptable approximation and the correspondingpicture.PROBLEM 5. What rank-r approximation exactly reproduces the original picture?c2011 Stefania Tracogna, SoMSS, ASU 4MATLAB sessions: Laboratory 8PROBLEM 6. Suppose that a rank-k approximation, for some k, gives an acceptable approximation.How much data is needed to represent the rank-k approximation? (Hint: you need k columnsof U, k columns of V , and k singular values of S.) The ratio of the amount of data used for theapproximation and the amount of data of the (original format of the) picture is the compressionrate. Find the compression rate for the value of the rank you determined in Problem 4. What doesthe compression rate represent?PROBLEM 7. If we use a high rank approximation, then the amount of data needed for theapproximation may exceed the amount of data used in the original representation of the picture.Find the smallest value of k such that the rank-k approximation of the matrix uses the same ormore amount of data as the original picture. Approximations with ranks higher than this k cannotbe used for image compression since they defeat the purpose of using less data than in the originalrepresentation.c2011 Stefania Tracogna, SoMSS, ASU 5
Reviews
There are no reviews yet.