%% Blatt 4 - Programmieraufgaben:
%% Pivotisierte Cholesky-Zerlegung (4 Punkte)
% Vorname Nachname Mat.-Nr.
%% Pivotisierte Cholesky-Zerlegung
% Die pivotisierte Cholesky-Zerlegung liefert für symmetrisch positiv semi-definite Matrizen im j-ten Eliminationsschritt eine gute Approximation bezüglich der Spurnorm, \Vert \mathbf{P}\mathbf{A}\mathbf{P}^\top - \mathbf{L}^{(j)} (\mathbf{L}^{(j)})^\top \Vert_{\text{tr}} < \varepsilon. Ist für ein vorgegebenes \varepsilon > 0 hierbei j \ll n, so spricht man auch von einer Niedrigrangapproximation an \mathbf{A}. Um eine solche direkt berechnen zu können, wandelt man den Algorithmus weiter ab, sodass beim Erreichen der gewünschten Fehlertoleranz \varepsilon die Berechnung abgebrochen wird. Dabei lässt sich der Algorithmus gestalten, dass man die Matrix \mathbf{L} separat speichert und die Einträge permutiert, damit direkt \mathbf{A} \approx \mathbf{L} \mathbf{L}^\top gilt. Andererseits muss die Matrix \mathbf{A} dann nicht mehr explizit gespeichert werden, sondern es reicht den Zugriff auf die aktuell benötigten Matrixeinträge (z.B. per Berechnung) zu gewährleisten. Das heisst, statt einer Matrix \mathbf{A} übergeben wir ein function handle \mathtt{a}, welches durch den Aufruf von \mathtt{a}(i, j) den Eintrag a_{i, j} berechnet; speziell sollen dabei auch Blöcke auswertbar sein, also
% \mathtt{a}(i:l, j:k) = \mathbf{A}_{i:l, j:k}.
% Dies führt auf den folgenden Algorithmus:
% Input: function handle \mathtt{a}, mit \mathtt{a}(i, j) = a_{i, j} einer symmetrisch positiv semi-definiten Matrix \mathbf{A} = [a_{i, j}]_{i, j=1}^n \in \mathbb{R}^{n \times n}, Dimension n und Toleranz \varepsilon
% Output: pivotisierte Dreiecksmatrix \mathbf{L} = [\ell_{i, j}]_{i, j=1}^{n, k} \in \mathbb{R}^{n \times k} der Niedrigrang-Approximation \mathbf{L} \mathbf{L}^\top \approx \mathbf{A}
% initialisiere \mathbf{L} := []
% initialisiere \mathbf{p} := [1, 2, \ldots, n]^\top
% initialisiere \mathbf{d} := [\mathtt{a}(1,1), \mathtt{a}(2,2), \ldots, \mathtt{a}(n,n)]^\top
% setze \delta := \sum_{r=1}^{n} d_r
% setze j := 1
% while j \leq n & \delta \geq \varepsilon do
% \qquadsetze k := {\arg\max}_{r \in \{j, j+1, \ldots, n\}} d_{p_r}
% \qquadvertausche die Einträge p_k und p_j
% \qquadsetze \ell_{p_j, j} := \Big(d_{p_j} \Big)^{1/2}
% \qquadsetze \mathbf{L}_{\mathbf{p}_{j+1:n}, j} := \mathtt{a}(\mathbf{p}_{j+1:n}, p_j) / \ell_{p_j, j}
% \qquadupdate \mathbf{L}_{\mathbf{p}_{j+1:n},j} := \mathbf{L}_{\mathbf{p}_{j+1:n},j}-\mathbf{L}_{\mathbf{p}_{j+1:n},1:j-1} \mathbf{L}_{p_j,1:j-1}^\top / \ell_{p_j,j}\)
% \qquadupdate \mathbf{d}_{\mathbf{p}_{j+1:n}} := \mathbf{d}_{\mathbf{p}_{j+1:n}} - [\ell_{p_{j+1},j}^2, \ell_{p_{j+2},j}^2, \ldots, \ell_{p_n,j}^2]^\top
% \qquadupdate \delta := \sum_{r=j+1}^{n} d_{p_r}
% \qquadsetze j := j+1
% end for
% Schreiben Sie am Ende dieses Live Scripts eine MATLAB-Funktion
% [L, p] = choleskyLowrank(a, n, epsilon) ,
% welche diesen Algorithmus implementiert.
%
% Damit lässt sich das folgende Testproblem rechnen. Hierbei fällt einerseits auf, dass diese Matrizen (meist) keine echten Niedrigrangapproximationen besitzen, andernseits aber, dass der Algorithmus stabil ist:

R = rand(2^10, 2^8);
A = R * R.';
A = A ./ max(abs(A), [], "all");

n = size(A, 1);
epsilon = 1e-1;

tic;
[L, p] = choleskyLowrank(A, n, epsilon);
toc
max(abs(A - L*L.'), [], 'all')
rank(L)

% Ein interessanteres Beispiel ergibt die Betrachtung des Exponentialkerns
% k(x, y) = e^{-(x-y)^2} auf [0, 1]^2,
% diskretisiert an (n+1) \times (n+1) äquidistanten Gitterpunkten (x_i, y_j) mit x_i = (i-1)/n und y_j = (j-1)/n. Die Matrix \mathbf{A} = [k(x_i, y_j)]_{i,j=1}^{n+1} besitzt dann eine hoch effiziente Niedrigrangapproximation. Wir berechnen diejenige mit n=2^{12} und für \varepsilon = 10^{-d} mit d = 3, 5, \ldots, 13:

n = 2^12;

a = @(is, js) exp(-((is.'-1)/n - (js-1)/n).^2);
A = a(1:n+1, 1:n+1);

for d=3:2:13
    epsilon = 10^(-d);
    disp(['For Precision: ' num2str(epsilon)]);

    tic;
    [L, p] = choleskyLowrank(a, n+1, epsilon);
    Ap = L * L.';
    t = toc;
    disp(['  Time: ' num2str(t) 's']);
    disp(['  Rank: ' num2str(size(L, 2))]);
    disp(['  Error: ' num2str(max(abs(A - Ap), [], 'all'))]);
end

%
%
% Funktionen:

% TODO: implement choleskyLowrank