ICS/hw3/assignment3.m

%% Assignment 2
% Name: Claudio Maggioni
%
% Date: 19/3/2019
%
% This is a template file for the first assignment to get started with running
% and publishing code in Matlab.  Each problem has its own section (delineated
% by |%%|) and can be run in isolation by clicking into the particular section
% and pressing |Ctrl| + |Enter| (evaluate current section).
%
% 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('assignment3.m', options)
%

%% Problem 3
format long
A = [4 3 2 1; 8 8 5 2; 16 12 10 5; 32 24 20 11];
[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);
    n = dimensions(1);
    
    p = 1:n;
    L = zeros(n);
    U = zeros(n);

    for i = 1:n
        values = A(:,i);
        values(values == 0) = -Inf;
        [~, 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;
        
        L(:,i) = A(:,i) / A(p(i),i);
        U(i,:) = A(p(i),:);
        A = A - L(:,i) * U(i,:);
    end
    
    I = eye(n);
    P = zeros(n);
    for i = 1:n
        P(:,i) = I(:,p(i));
    end
    P = transpose(P);
    L = P * L;   
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