Solve the following problem in MATLAB: Write a routine that uses a one-sided pla

Question

Solve the following problem in MATLAB:

Write a routine that uses a one-sided plane rotation to symmetrize an arbitrary 2 times 2 matrix. That is, given a 2 times 2 matrix A, choose c and s so that [c -s -s c] [a_11 a_12 a_21 a_22] = [b_11 b_12 b_12 b_22] is symmetric. Write a routine that uses a two-sided plane rotation to annihilate the off-diagonal entries of an arbitrary 2 times 2 symmetric matrix. That is, given a symmetric 2 times 2 matrix B, choose c and s so that [c -s -s c] [b_11 b_12 b_12 b_22] [c s -s c] = [d_11 0 0 d_22] is diagonal. Combine the two routines developed in parts a and b to obtain a routine for computing the singular value decomposition A = U sigma V^T of an arbitrary 2 times 2 matrix A. Note that U will be a product to two plane rotations, whereas V will be a single plane rotation. Test your routine on a few randomly chosen 2 times 2 matrices and compare the results with those for a library SVD routine. By systematically solving successive 2 times 2 sub-problems, the module you have just developed can be used to compute the SVD of an arbitrary m times n matrix in a manner analogous to the Jacobi method for the symmetric eigenvalue problem.

Explanation / Answer

rng(0); % ensure same results m=200; % number of data rows r=10000; % number of bootstrap repeats A=randn(m,4); A(:,2)=A(:,1)+A(:,2); b=A(:,1)+0.5*A(:,2); A(:,5)=A(:,1); b1=b+0.01*randn(m,1); %% Simple regression case % Columns 1 and 5 are the same giving redundent variable columns, % and the regression matrix is singular. % Columns 3 and 4 are not related to the dependent variable b. % % The result includes both variables 1 and 5 in equal amounts: disp(' ') disp('Test case 1') [x1,vx1,info1]=linfitregsel(A,b1) %% Regression case with specified variables % Variables 1 and 2 are available for selection, % and variable 3 is forced into the regression. % % Variable 3 is now not set to zero, but receives a small coefficient: disp(' ') disp('Test case 2') [x2,vx2,info2]=linfitregsel(A,b1,struct('sel',[2 2 1 0 0])) %% Bootstrap example % r values of the regression coefficients are collected. % The variance of these is compared with the average calculated variance. % % The t-value for rejection is set high to exclude the occasional % coefficient for columns 3 and 4 that gets a larger t-value. % % Mean values of coefficient are close to values used to generate data. % Bootstrap and calculated values of standard deviation and % also agree to expecteed level: tic disp(' ') disp('Bootstrap test case (about 8 seconds)') x=zeros(r,5); var1=zeros(5,5); for i=1:r [xx,vmx,info]=linfitregsel(A,b+0.01*randn(m,1),struct('tval',6)); var1=var1+vmx; x(i,:)=xx'; if(mod(i,10000)==0) fprintf('%i/%i %#6.1f seconds ',i,r,toc) end end var1=var1/r; toc disp(' ') disp('Mean value of x') disp(mean(x)) disp(' ') disp('Standard deviation of x values (from bootstrap)') disp(std(x)); disp(' ') disp('Standard deviation from linfitregsel variance matrix') disp(sqrt(diag(var1)')) disp(' ') disp(' ') disp('Variance of x values (from bootstrap)') xx=bsxfun(@minus,x,mean(x)); disp(xx'*xx/r) disp(' ') disp('Calculated variance from linfitregsel variance matrix') disp(var1) %% Comparison with ordinary linear regression % Generate some data rng(0) m=10; n=4; a=randn(m,n); a1=a+0.01*randn(m,n); a2=a+0.01*randn(m,n); aa=[a1,a2]; b=sum(a,2); %% Solve least squares fits % First using linfitregsel and second using % Solution from linfitregsel is close to expected values, % while gives much different values disp(' ') disp('Solution from linfitregsel') [x1,~,info]=linfitregsel(aa,b) disp(' ') disp('Solution from ') x2=aa %% Fit to a new set of data % Compare accuracy of predictions using new data % linfitregsel gives a more accuray fit to new data % new random pertubation into data a1a=a+0.01*randn(m,n); a2a=a+0.01*randn(m,n); aaa=[a1a,a2a]; disp(' ') disp('Root mean square of residuals using linfitregsel solution') t1=aaa*x1-b; disp(sqrt(t1'*t1/length(t1))) disp('Root mean square of residuals using solution') t2=aaa*x2-b; disp(sqrt(t2'*t2/length(t2)))

Related Questions

Navigate

• Browse (All)
• Browse S
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.