daily pastebin goal
83%
SHARE
TWEET

Untitled

a guest Nov 14th, 2017 57 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
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top