2071 Lab07 Summary M. M. Sussman $Id: summary.txt,v 1.12 2017/03/14 19:05:06 mike Exp $ {@(#) Tue Jan 31 11:08:45 2017 } EXERCISE 1 1. unstable_gs.m written 2. Q=identity 3. Q= 5.7735e-01 4.0825e-01 7.0711e-01 5.7735e-01 4.0825e-01 -7.0711e-01 5.7735e-01 -8.1650e-01 0 4. Yes, got (was in web page) 0.89443 0.35857 0.19518 -0.44721 0.71714 0.39036 0.00000 -0.59761 0.58554 0.00000 0.00000 -0.68313 verify: Q'*Q=identity (1-liner: norm(Q'*Q-eye(3),'fro')=0) 5. Q is not orthogonal because Q*Q' is not the identity. surprised? (or because Q is not square, so cannot be invertible) 6. Q not orthog: norm(Q1'*Q1-eye(10),'fro')=3.4632 ~=0 DESCRIBE TESTING 7. B1=Q1'*Q1: B1(1:3,1:3) = identity up to roundoff B1(8:10,8:10)= 1.00000 -0.99996 -0.99991 -0.99996 1.00000 0.99999 -0.99991 0.99999 1.00000 EXERCISE 2 1. modified_gs.m 2. compare against unstable_gs for at least 3 3X3 matrices: DESCRIBE FULLY [A=rand(3,3);Qm=modified_gs(A);Qu=unstable_gs(A);norm(Qm-Qu,'fro')] 3. norm(B1-eye(10),'fro')=3.4619, norm(B2-eye(10),'fro')= 3.5403e-4 EXERCISE 3 1. gs_factor.m 2. compare Q from gs_factor with Q from gram_schmidt.m: same also check that A=Q*R Q = [0.57735 0.40825 0.70711 0.57735 0.40825 -0.70711 0.57735 -0.81650 0.00000] R = [1.73205 1.15470 0.57735 0.00000 0.81650 0.40825 0.00000 0.00000 0.70711] 3. norm(Q'*Q-eye(100),'fro')=roundoff norm(Q*Q'-eye(100),'fro')=roundoff => Q is orthogonal R is upper triangular: norm(R-triu(R),'fro') is zero also norm(tril(R,-1),'fro') is zero A=Q*R: norm(Q*R-A,'fro') is roundoff 4. Q = 0 0 1.0000 0.0463 0.9707 0.0000 0.2314 -0.2400 -0.0000 0.9718 0.0109 -0.0000 R = 21.6102 34.9835 56.5937 0 0.3927 0.3927 0 0 1.0000 Q'*Q is identity Q*Q' is not identity Q is not orthogonal R is upper triangular A=Q*R norm(Q*R-A,'fro') is roundoff EXERCISE 4 1-3 householder.m, comments explain k 4. Check orthogonality of H both times k=1 k=4 -1.9621e+01 1.0000e+01 -4.0419e-16 9.0000e+00 1.0113e-15 8.0000e+00 -7.2858e-17 -1.1832e+01 1.8562e-16 -1.6653e-16 2.7756e-16 -3.4694e-18 3.5302e-16 4.9266e-16 -6.8522e-17 3.7470e-16 1.5743e-16 2.8796e-16 1.5959e-16 1.1102e-16 EXERCISE 5 1. A=magic(5) 2. H1 = -0.523387 -0.708111 -0.123150 -0.307875 -0.338662 -0.708111 0.670851 -0.057243 -0.143108 -0.157419 -0.123150 -0.057243 0.990045 -0.024888 -0.027377 -0.307875 -0.143108 -0.024888 0.937779 -0.068443 -0.338662 -0.157419 -0.027377 -0.068443 0.924713 A1 = -3.2481e+01 -2.6631e+01 -2.1397e+01 -2.3706e+01 -2.5861e+01 0 -1.8535e+01 -3.4109e+00 -7.3797e-01 -2.9935e+00 0 1.9070e+00 1.1189e+01 1.7437e+01 1.8697e+01 0 1.7675e+00 1.4474e+01 1.4592e+01 -5.2580e+00 0 6.7443e+00 2.0021e+01 -5.0486e+00 -8.3855e-02 3. H2 = 1.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.93166 0.09586 0.08885 0.33901 0.00000 0.09586 0.99524 -0.00441 -0.01682 0.00000 0.08885 -0.00441 0.99591 -0.01559 0.00000 0.33901 -0.01682 -0.01559 0.94050 A2 = -3.2481e+01 -2.6631e+01 -2.1397e+01 -2.3706e+01 -2.5861e+01 0 1.9894e+01 1.2323e+01 1.9439e+00 4.0856e+00 0 0 1.0409e+01 1.7304e+01 1.8345e+01 0 0 1.3750e+01 1.4469e+01 -5.5836e+00 0 0 1.7260e+01 -5.5193e+00 -1.3262e+00 EXERCISE 6 1. h_factor.m 2. Q = -0 1 0 -1 -0 0 0 0 -1 R = -1 -1 -1 -0 1 1 0 0 -1 3. h_factor(magic(5)) Q = -0.523387 0.505754 0.673470 -0.121542 0.044104 -0.708111 -0.696573 -0.017727 0.081542 0.080002 -0.123150 0.136742 -0.355751 -0.630744 0.664635 -0.307875 0.191057 -0.412231 -0.424672 -0.720021 -0.338662 0.451439 -0.499631 0.632767 0.177441 R = -3.2481e+01 -2.6631e+01 -2.1397e+01 -2.3706e+01 -2.5861e+01 -1.2948e-15 1.9894e+01 1.2323e+01 1.9439e+00 4.0856e+00 8.0813e-16 8.9824e-16 -2.4399e+01 -1.1632e+01 -3.7415e+00 3.6963e-16 -5.2526e-16 1.9380e-15 -2.0098e+01 -9.9739e+00 1.5466e-16 -4.2955e-16 -7.5355e-16 -2.1185e-15 1.6000e+01 Agrees with gs_factor except cols 1, 3, 4 of Q are mult by -1 as are rows 1, 3, 4 of R 4. Agrees with qr, last column Q mult by -1, last row of R mult by -1 5. norm(B3-eye(10),'fro')=roundoff, much smaller than B2 or B1 from ex2 EXERCISE 7 h_solve.m (x=u_solve(R,Q'*b);) x=x2=[1 2 3 4 5]', norm of difference is roundoff EXERCISE 8 1. L1 is not the Cholesky factor because it has some neg. diag. 2. cholesky.m 3. A=eye(6)+pascal(6), check L*L'=A 4. Choose matrix, compare cholesky with chol: norm(L'-chol(A),'fro')=0 Describe tests! 5. Download l_solve.m 6. Use given code to generate SPD matrix A and x vector s.t. b=A*x 7. solve using cholesky, l_solve, u_solve. L2 error should be roundoff. L=cholesky(A); y=l_solve(L,b); x0=u_solve(L',y); norm(x-x0)