Merge branch 'master' of tea.maggioni.xyz:maggicl/ICS

hw3/assignment3.m

@@ -11,13 +11,24 @@
 % To generate a pdf for submission in your current directory, use the following
 % three lines of code at the command window:
 %
-%  >> options.format = 'pdf'; options.outputDir = pwd; publish('assignment2.m', options)
+%  >> options.format = 'pdf'; options.outputDir = pwd; publish('assignment3.m', options)
 %
 
 %% Problem 3
-format rational
+format long
 A = [4 3 2 1; 8 8 5 2; 16 12 10 5; 32 24 20 11];
-[L,U,P] = pivotedOuterProductLU(A)
+[L,U,P] = pivotedOuterProductLU(A);
+display(L);
+display(U);
+display(P);
+
+%% Problem 6
+[X,Y] = meshgrid((0:2000)/2000); 
+A = exp(-sqrt((X-Y).^2));
+L = outerProductCholesky(A);
+disp(norm(L*L'-A, 'fro'));
+
+%% Problem 3 (continued)
 
 function [L,U,P] = pivotedOuterProductLU(A)
     dimensions = size(A);
@@ -33,6 +44,14 @@ function [L,U,P] = pivotedOuterProductLU(A)
         [~, p_k] = max(values);
         k = find(p == p_k);
         
+        if k == -Inf
+            disp("Matrix is singular");
+            L = [];
+            U = [];
+            P = [];
+            return
+        end
+        
         p(k) = p(i);
         p(i) = p_k;
         
@@ -47,6 +66,23 @@ function [L,U,P] = pivotedOuterProductLU(A)
         P(:,i) = I(:,p(i));
     end
     P = transpose(P);
-    L
     L = P * L;   
-end
+end
+
+%% Problem 6 (continued)
+
+function [L] = outerProductCholesky(A)
+    if ~all(eig(A) >= 0)
+        disp("matrix is not positive definite");
+        L = [];
+        return;
+    end
+    dimensions = size(A);
+    n = dimensions(1);
+    L = zeros(n);
+    
+    for i = 1:n
+        L(:, i) = A(:, i) / sqrt(A(i,i));
+        A = A - L(:, i) * transpose(L(:, i));
+    end
+end