SHARE
TWEET

Untitled

a guest Nov 14th, 2017 56 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. %% Q3b
  2. % Name: William Qi
  3. % Date: Nov. 13, 2017
  4. % Assignment 4
  5.  
  6. % Set up A matrices with varying epsilon values in the lower rows
  7. A1 = [1 1 1; 2*sqrt(eps) 0 0; 0 2*sqrt(eps) 0; 0 0 2*sqrt(eps)];
  8. A2 = [1 1 1; sqrt(eps) 0 0; 0 sqrt(eps) 0; 0 0 sqrt(eps)];
  9. A3 = [1 1 1; 2*eps 0 0; 0 2*eps 0; 0 0 2*eps];
  10. A4 = [1 1 1; eps 0 0; 0 eps 0; 0 0 eps];
  11. b = [1; 0; 0; 0];
  12.  
  13. % Set up solutions based off of exact solution from 3A
  14. x1 = 1/(3 + (2*sqrt(eps))^2)
  15. x2 = 1/(3 + (sqrt(eps))^2)
  16. x3 = 1/(3 + (2*eps)^2)
  17. x4 = 1/(3 + (eps)^2)
  18.  
  19. % i
  20. x1_i_err = norm(normalSolve(A1, b) - x1);
  21. x2_i_err = norm(normalSolve(A2, b) - x2);
  22. x3_i_err = norm(normalSolve(A3, b) - x3);
  23. x4_i_err = norm(normalSolve(A4, b) - x4);
  24.  
  25. % ii
  26. x1_ii_err = norm(classicalGSSolve(A1, b) - x1);
  27. x2_ii_err = norm(classicalGSSolve(A2, b) - x2);
  28. x3_ii_err = norm(classicalGSSolve(A3, b) - x3);
  29. x4_ii_err = norm(classicalGSSolve(A4, b) - x4);
  30.  
  31. % iii
  32. x1_iii_err = norm(modifiedGSSolve(A1, b) - x1);
  33. x2_iii_err = norm(modifiedGSSolve(A2, b) - x2);
  34. x3_iii_err = norm(modifiedGSSolve(A3, b) - x3);
  35. x4_iii_err = norm(modifiedGSSolve(A4, b) - x4);
  36.  
  37. % iv
  38. x1_iv_err = norm(householderSolve(A1, b) - x1);
  39. x2_iv_err = norm(householderSolve(A2, b) - x2);
  40. x3_iv_err = norm(householderSolve(A3, b) - x3);
  41. x4_iv_err = norm(householderSolve(A4, b) - x4);
  42.  
  43. % v
  44. x1_v_err = norm(givensSolve(A1, b) - x1);
  45. x2_v_err = norm(givensSolve(A2, b) - x2);
  46. x3_v_err = norm(givensSolve(A3, b) - x3);
  47. x4_v_err = norm(givensSolve(A4, b) - x4);
  48.  
  49. % vi
  50. x1_vi_err = norm(qrSolve(A1, b) - x1);
  51. x2_vi_err = norm(qrSolve(A2, b) - x2);
  52. x3_vi_err = norm(qrSolve(A3, b) - x3);
  53. x4_vi_err = norm(qrSolve(A4, b) - x4);
  54.  
  55. % vii
  56. x1_vii_err = norm(backslashSolve(A1, b) - x1);
  57. x2_vii_err = norm(backslashSolve(A2, b) - x2);
  58. x3_vii_err = norm(backslashSolve(A3, b) - x3);
  59. x4_vii_err = norm(backslashSolve(A4, b) - x4);
  60.  
  61. % viii
  62. x1_viii_err = norm(svdSolve(A1, b) - x1);
  63. x2_viii_err = norm(svdSolve(A2, b) - x2);
  64. x3_viii_err = norm(svdSolve(A3, b) - x3);
  65. x4_viii_err = norm(svdSolve(A4, b) - x4);
  66.  
  67. %ix
  68. x1_ix_err = norm(tsvdSolve(A1, b) - x1);
  69. x2_ix_err = norm(tsvdSolve(A2, b) - x2);
  70. x3_ix_err = norm(tsvdSolve(A3, b) - x3);
  71. x4_ix_err = norm(tsvdSolve(A4, b) - x4);
  72.  
  73. %% Function Definitions
  74. % i
  75. function x = normalSolve(A, b)
  76.     % Set up system of normal equations and use backslash to solve
  77.     x = (A.' * A) \ (A.' * b);
  78. end
  79.  
  80. % ii
  81. function x = classicalGSSolve(A, b)
  82.     % Compute QR decomposition using classical Gram-Schmidt orthgonalization
  83.     [~, n] = size(A);
  84.     R = zeros(n);
  85.  
  86.     % r_11 = ||q_1||
  87.     R(1,1) = norm(A(:,1));
  88.     % q_1 = q_1/r_11
  89.     Q(:,1) = A(:,1)/R(1,1);
  90.     for j = 2:n
  91.         R(1:j-1, j) = Q(:, 1:j-1)' * A(:, j);
  92.         Q(:, j) = A(:, j) - Q(:, 1:j-1) * R(1:j-1, j);
  93.         % r_jj = ||q_j||
  94.         R(j, j)= norm(Q(:, j));
  95.         % q_j = q_j/r_jj
  96.         Q(:, j) = Q(:, j) / R(j, j);
  97.     end
  98.    
  99.     % Solve system Rx=Q^T*b
  100.     y = Q.' * b;
  101.     x = R \ y;
  102. end
  103.  
  104. % iii
  105. function x = modifiedGSSolve(A, b)
  106.     % Compute QR decomposition using modified Gram-Schmidt orthgonalization
  107.     [m, n] = size(A);
  108.     Q = zeros(m, n);
  109.     R = zeros(n, n);
  110.  
  111.     for j = 1:n
  112.         % q_j = a_j
  113.         Q(:, j) = A(:, j);
  114.         for i = 1:j-1
  115.             % r_ij = <q_j, q_i>
  116.             R(i, j) = Q(:,i)' * Q(:, j);
  117.             % q_j = q_j = r_ij * q_i
  118.             Q(:, j) = Q(:, j) - R(i, j) * Q(:, i);
  119.         end
  120.         % r_jj = ||q_j||
  121.         R(j, j) = norm(Q(:, j));
  122.         % q_j = q_j/r_jj
  123.         Q(:, j) = Q(:, j) / R(j, j);
  124.     end
  125.    
  126.     % Solve system Rx=Q^T*b
  127.     y = Q.' * b;
  128.     x = R \ y;
  129. end
  130.  
  131. % iv
  132. function x = householderSolve(A, b)
  133.     [m, n] = size(A);
  134.     p = zeros(1, n);
  135.    
  136.     % Compute QR decomposition using Householder reflections
  137.     for k = 1:n
  138.         % Define u of length = m-k+1
  139.         z = A(k:m, k);
  140.         e1 = [1; zeros(m-k, 1)];
  141.         u = z + sign(z(1))*norm(z)*e1;
  142.         u = u / norm(u);
  143.         % Update nonzero part of A by I-2uu^T
  144.         A(k:m, k:n) = A(k:m, k:n) - 2*u*(u'*A(k:m, k:n));
  145.         % Store u
  146.         p(k) = u(1);
  147.         A(k+1:m, k) = u(2:m-k+1);
  148.     end
  149.    
  150.     % Use QR decomposition to solve least squares
  151.     y = b(:);
  152.     % Transform b
  153.     for k=1:n
  154.         u = [p(k);A(k+1:m, k)];
  155.         y(k:m) = y(k:m) - 2*u*(u'*y(k:m));
  156.     end
  157.     % Form upper triangular R and solve
  158.     R = triu(A(1:n, :));
  159.     x = R \ y(1:n);
  160. end
  161.  
  162. % v
  163. function x = givensSolve(A, b)
  164.     [m,n] = size(A);
  165.     Q = eye(m);
  166.     R = A;
  167.  
  168.     % Use Givens rotation to zero out a single entry at a time
  169.     for j = 1:n
  170.         for i = m:-1:(j+1)
  171.             G = eye(m);
  172.             % Set up G_ij by computing c, s
  173.             [c,s] = givensrotation(R(i-1, j), R(i, j));
  174.             G([i-1, i],[i-1, i]) = [c -s; s c];
  175.            
  176.             % Apply G_ij on A
  177.             R = G' * R;
  178.             Q = Q * G;
  179.         end
  180.     end
  181.  
  182.     % Solve system Rx=Q^T*b
  183.     y = Q' * b;
  184.     x = R \ y;
  185. end
  186.  
  187. function [c,s] = givensrotation(a,b)
  188.   % Check if a_ij is 0
  189.   if b == 0
  190.     c = 1;
  191.     s = 0;
  192.   else
  193.     % r = sqrt(a^2 + b^2)
  194.     % c = a/r
  195.     % s = b/r
  196.     if abs(b) > abs(a)
  197.       r = a / b;
  198.       s = 1 / sqrt(1 + r^2);
  199.       c = s*r;
  200.     else
  201.       r = b / a;
  202.       c = 1 / sqrt(1 + r^2);
  203.       s = c*r;
  204.     end
  205.   end
  206. end
  207.  
  208. % vi
  209. function x = qrSolve(A, b)
  210.     % Obtain QR decomposition using matlab builtin function
  211.     [Q,R] = qr(A);
  212.    
  213.     % Solve system Rx=Q^T*b
  214.     y = Q' * b;
  215.     x = R \ y;
  216. end
  217.  
  218. % vii
  219. function x = backslashSolve(A, b)
  220.     % Solve least squares using backslash
  221.     x = A \ b;
  222. end
  223.  
  224. % viii
  225. function x = svdSolve(A, b)
  226.     % Obtain SVD decomposition using matlab builtin function
  227.     [U, S, V] = svd(A);
  228.    
  229.     % Format y and z with the correct dimensions
  230.     y = diag(S).^-1;
  231.     z = U' * b;
  232.     z = z(1:numel(y));
  233.    
  234.     % x = V * S^-1 * U^T * b
  235.     x = V * (diag(y) * z);
  236. end
  237.  
  238. % ix
  239. function x = tsvdSolve(A, b)
  240.     % Obtain SVD decomposition using matlab builtin function
  241.     [U, S, V] = svd(A);
  242.     % Set singluar value truncation threshold to 10^-6
  243.     threshold = 10^-6;
  244.    
  245.     % Find the index of the last element in S that is above the threshold
  246.     r = find(diag(S) < threshold, 1) - 1;
  247.    
  248.     % Format y, z, and V with the correct truncated dimensions
  249.     y = diag(S).^-1;
  250.     z = U' * b;
  251.     y = y(1:r);
  252.     z = z(1:r);
  253.     V = V(:, 1:r);
  254.    
  255.     % x = V * S^-1 * U^T * b with truncated dimensions
  256.     x = V * (diag(y) * z);
  257. end
RAW Paste Data
Top