%% Chapter 6 Exercise Solutions
%
%% Exercise 6.1

A = sym([2, 3, -1 ; 4, 9, -1 ; 0, 3, 2]);
b = sym([4 ; 18 ; 11]);

L = eye(3);
U = zeros(3);

% j = 1
U(1,1) = A(1,1);
L(2,1) = A(2,1) / U(1,1);
L(3,1) = A(3,1) / U(1,1);

% j = 2
U(1,2) = A(1,2);
U(2,2) = A(2,2) - L(2,1) * U(1,2);
L(3,2) = 1 / U(2,2) * (A(3,2) - L(3,1) * U(1,2));

% j = 3
U(1,3) = A(1,3);
U(2,3) = A(2,3) - L(2,1) * U(1,3);
U(3,3) = A(3,3) - L(3,1) * U(1,3) - L(3,2) * U(2,3);

L
U

% Solving L\mathbf{y} = \mathbf{b}

y = L \ b

% Solving U \mathbf{x} = \mathbf{y}

x = U \ y

%
%% Exercise 6.2

A = sym([1, 3, 2 ; 2, 1, 0 ; 3, 1 -1])
b = [16 ; 0 ; -5];
P = sym(eye(3));

% Apply partial pivoting to calculate P

A([1,3],:) = A([3,1],:)
A([2,3],:) = A([3,2],:)
P([1,3],:) = P([3,1],:);
P([2,3],:) = P([3,2],:)

% Calculate LUP decomposition

L = eye(3);
U = zeros(3);


% j = 1
U(1,1) = A(1,1);
L(2,1) = A(2,1) / U(1,1);
L(3,1) = A(3,1) / U(1,1);

% j = 2
U(1,2) = A(1,2);
U(2,2) = A(2,2) - L(2,1) * U(1,2);
L(3,2) = 1 / U(2,2) * (A(3,2) - L(3,1) * U(1,2));

% j = 3
U(1,3) = A(1,3);
U(2,3) = A(2,3) - L(2,1) * U(1,3);
U(3,3) = A(3,3) - L(3,1) * U(1,3) - L(3,2) * U(2,3);

L
U

% Solve L \mathbf{y} = P \mathbf{b}

y = L \ (P * b)

% Solve U\mathbf{x} = \mathbf{y}

x = U \ y

%
%% Exercise 6.3

clear
A = sym([16, 16, 4 ; 16, 25, 10 ; 4, 10, 6])
b = [-8 ; -47 ; -30];

% Check whether A is positive definte.

A'
det(A(1,1))
det(A(1:2,1:2))
det(A)

% Calculate Cholesky decomplosition

L(1,1) = sqrt(A(1,1));
L(2,1) = 1 / L(1,1) * (A(2,1));
L(3,1) = 1 / L(1,1) * (A(3,1));
L(2,2) = sqrt(A(2,2) - L(2,1)^2);
L(3,2) = 1 / L(2,2) * (A(3,2) - L(3,1) * L(2,1));
L(3,3) = sqrt(A(3,3) - L(3,1)^2 - L(3,2)^2)

% Solve L \mathbf{y} = \mathbf{b}

y = sym(zeros(3,1));
y(1) = b(1) / L(1,1);
y(2) = (b(2) - L(2,1) * y(1)) / L(2,2);
y(3) = (b(3) - L(3,1) * y(1) - L(3,2) * y(2)) / L(3,3)

% Solve L^T \mathbf{x} = \mathbf{y}

L = L';
x = sym(zeros(3,1));
x(3) = y(3) / L(3,3);
x(2) = (y(2) - L(2,3) * x(3)) / L(2,2);
x(1) = (y(3) - L(1,2) * x(2) - L(1,3) * x(3)) / L(1,1)

%
%% Exercise 6.4
% Exercise 6.1 solution

function [L, U] = lu(A)

n = size(A, 1);
L = eye(n);
U = zeros(n);

for j = 1 : n
    for i = 1 : j
        sum_ = 0;
        for k = 1 : i - 1
            sum_ = sum_ + L(i,k) * U(k,j);
        end
        U(i,j) = A(i,j) - sum_;
    end

    for i = j + 1 : n
        sum_ = 0;
        for k = 1 : i - 1
            sum_ = sum_ + L(i,k) * U(k,j);
        end
        L(i,j) = (A(i,j) - sum_) / U(j,j);
    end
end

end

function x = forward_substitution(L, b)

n = size(L, 1);
x = zeros(n, 1);
for i = 1 : n
    sum_ = 0;
    for j = 1 : i - 1
        sum_ = sum_ + L(i,j) * x(j);
    end
    x(i) = (b(i) - sum_) / L(i,i);
end

end

function x = back_substitution(U, b)

n = size(U, 1);
x = zeros(n, 1);
x(n) = b(n) / U(n,n);
for i = n - 1 : -1 : 1
    sum_ = 0;
    for j = i + 1 : n
        sum_ = sum_ + U(i,j) * x(j);
    end
    x(i) = (b(i) - sum_) / U(i,i);
end

end

A = [ 2, 3, -1 ; 4, 9, -1, ; 0, 3, 2 ];
b = [ 4 ; 18 ; 11 ];

[L, U] = lu(A)
y = forward_substitution(L, b)
x = back_substitution(U, y)

% Exercise 6.2 solution

function P = partial_pivot(A)

n = size(A, 1);
P = eye(n);
for j = 1 : n
    maxpivot = abs(A(j,j));
    maxpivotrrow = j;
    for i = j + 1 : n
        if abs(A(i,j)) > maxpivot
            maxpivot = abs(A(i,j));
            maxpivotrow = i;
        end
    end
    A([j,maxpivotrow],:) = A([maxpivotrow,j],:);
    P([j,maxpivotrow],:) = P([maxpivotrow,j],:);
end

end

A = [1, 3, 2 ; 2, 1, 0 ; 3, 1, -1];
b = [16, ; 0 ; -5];

P = partial_pivot(A)
[L,U] = lu(P * A)
y = forward_substitution(L, P * b)
x = back_substitution(U, y)

% Exercise 6.3 solution

function L = cholesky(A)

n = size(A, 1);
L = zeros(n);

for j = 1 : n
    sum_ = 0;
    for k = 1 : j - 1
        sum_ = sum_ + L(j,k) ^ 2;
    end
    L(j,j) = sqrt(A(j,j) - sum_);

    for i = j + 1 : n
        sum_ = 0;
        for k = 1 : j - 1
            sum_ = sum_ + L(i,k) * L(j,k);
        end
        L(i,j) = (A(i,j) - sum_) / L(j,j);
    end
end

end


A = [16, 16, 4; 16, 25, 10; 4, 10, 6];
b = [-8 ; -47 ; -30];

L = cholesky(A)
y = forward_substitution(L, b)
x = back_substitution(L', y)

%
%% Exercise 6.5
% (a)

clear

A = sym([1, 1 ; -1, 0]);

% j = 1
u1 = A(:,1);
R(1,1) = norm(u1);
Q(:,1) = u1 / R(1,1);

% j = 2
R(1,2) = dot(Q(:,1), A(:,2));
u2 = A(:,2) - R(1,2) * Q(:,1);
R(2,2) = norm(u2);
Q(:,2) = u2 / R(2,2);
Q
R

% (b) A = \pmatrix{ 6 & 6 & 1 \cr 3 & 6 & 1 \cr 2 & 1 & 1 }

clear

A = sym([6, 6, 1 ; 3, 6, 1 ; 2, 1, 1]);

% j = 1
u1 = A(:,1);
R(1,1) = norm(u1);
Q(:,1) = u1 / R(1,1);

% j = 2
R(1,2) = dot(Q(:,1), A(:,2));
u2 = A(:,2) - R(1,2) * Q(:,1);
R(2,2) = norm(u2);
Q(:,2) = u2 / R(2,2);

% j = 3
R(1,3) = dot(Q(:,1), A(:,3));
R(2,3) = dot(Q(:,2), A(:,3));
u3 = A(:,3) - R(1,3) * Q(:,1) - R(2,3) * Q(:,2);
R(3,3) = norm(u3);
Q(:,3) = u3 / R(3,3);
Q
R

% (c) A = \pmatrix{ 1 & 2 & 1 \cr 1 & 4 & 3 \cr 1 & -4 & 6 \cr 1 & 2 & 1 }

clear

A = sym([1, 2, 1 ; 1, 4, 3 ; 1, -4, 6 ; 1, 2, 1]);

% j = 1
u1 = A(:,1);
R(1,1) = norm(u1);
Q(:,1) = u1 / R(1,1);

% j = 2
R(1,2) = dot(Q(:,1), A(:,2));
u2 = A(:,2) - R(1,2) * Q(:,1);
R(2,2) = norm(u2);
Q(:,2) = u2 / R(2,2);

% j = 3
R(1,3) = dot(Q(:,1), A(:,3));
R(2,3) = dot(Q(:,2), A(:,3));
u3 = A(:,3) - R(1,3) * Q(:,1) - R(2,3) * Q(:,2);
R(3,3) = norm(u3);
Q(:,3) = u3 / R(3,3);
Q
R

%
%% Exercise 6.6
% (a)

clear

A = sym([3, 0 ; 4, 5]);

Q = eye(2);
R = A;

% j = 1
x = R(:,1);
v = x - norm(x) * [1 ; 0];
H = eye(2) - 2 * v * v' / (v' * v)
R = H * R
Q = Q * H

% Second main diagonal element of R is negative, so multiply row 2 of R and column 2 of Q by -1.

R(2,:) = -R(2,:)
Q(:,2) = -Q(:,2)

% (b)

clear

A = sym([1, 2, 4 ; 0, 0, 5 ; 0, 3, 6]);
Q = eye(3);
R = A

% First column of R is already parallel to \mathbf{e}_1 so move to the second column.

% j = 2
x = R(2:3,2);
v = x - norm(x) * [1 ; 0]
Hdash = eye(2) - 2 * v * v' / (v' * v);
H = eye(3);
H(2:3,2:3) = Hdash
R = H * R
Q = Q * H

%
%% Exercise 6.7
% Exercise 6.5 solution

function [Q, R] = qr_gramschmidt(A)

[~, n] = size(A);
Q = A;
R = zeros(n, n);

for j = 1 : n
    for i = 1 : j - 1
        R(i,j) = dot(Q(:,i), A(:,j));
        Q(:,j) = Q(:,j) - R(i,j) * Q(:,i);
    end

    R(j,j) = norm(Q(:,j));
    Q(:,j) = Q(:,j) / R(j,j);
end

for i = 1 : n
    if R(i,i) < 0
        R(i,:) = -R(i,:);
        Q(:,i) = -Q(:,i);
    end
end

end


A = [1, 1 ; -1, 0];
[Q, R] = qr_gramschmidt(A)

% (b)

A = [6, 6, 1 ; 3, 6, 1 ; 2, 1 ,1];
[Q, R] = qr_gramschmidt(A)

% (c)

A = [1, 2, 1 ; 1, 4, 3 ; 1, -4, 6 ; 1, 2, 1];
[Q, R] = qr_gramschmidt(A)

% Exercise 6.6 solution
% (a)

function [Q, R] = qr_householder(A)

[m, n] = size(A);
Q = eye(m);
R = A;

for i = 1 : m - 1
    x = R(i:end,i);
    if sign(x(1)) > 0
        v = x + norm(x) * [1 ; zeros(m-i,1)];
    else
        v = x - norm(x) * [1 ; zeros(m-i,1)];
    end
    H = eye(m);
    H(i:end,i:end) = H(i:end,i:end) - 2 * v * v' / dot(v, v);
    R = H * R;
    Q = Q * H;
end

for i = 1 : n
    if R(i,i) < 0
        R(i,:) = -R(i,:);
        Q(:,i) = -Q(:,i);
    end
end

end


A = [3, 0 ; 4, 5];
[Q, R] = qr_householder(A)

% (b)

A = [1, 2, 4 ; 0, 0, 5 ; 0, 3, 6];
[Q, R] = qr_householder(A)

%