• API
• FAQ
• Tools
• Archive
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.

Top