%% Harmadik gyakorlat
%% Gauss elimináció, LU-felbontás
%
% 1. feladat: Írjunk programot, mely kiszámítja az alábbi egyenletrendszer megoldását a Cramer-szabály, illetve a Gauss-módszer segítségével!
% \left[\begin{array}{cccc}
% 8 & -6 & -7 & 5 \\
% 2 & -1 & -2 & 1 \\
% 6 & 1 & -9 & 2 \\
% -2 & -5 & 4 & 0 \\
% \end{array}\right]
% \left[ \begin{array}{c}
% x_1 \\ x_2 \\ x_3 \\ x_4
% \end{array} \right] = \left[ \begin{array}{c}
% -28 \\ -7 \\ -23 \\ 8
% \end{array} \right]
% Végezetül oldjuk meg az egyenletrendszert a rref parancs segítségével is.

% Maga a mátrix:
A = [8 -6 -7  5; 
     2 -1 -2  1; 
     6  1 -9  2; 
    -2 -5  4  0]

% A jobb oldal:
b = [-28 -7 -23 8]'

cramer_szabaly(A, b)

%

format short
% Gauss-elimináció
M = [A b]
% U meghatározása
for i=1:4
    % i. sor normalizálása
    M(i, i:end) = M(i, i:end) / M(i, i);
    
    % i. sor megfelelő együtthatókkal súlyozott változatainak
    % kivonása az alatta lévő sorokból
    M((i+1):end, :) = M((i+1):end, :) - M((i+1):end, i) * M(i, :);
end

disp('M immár felsőháromszög alakú, leszámítva az utolsó oszlopot:')
M

% Visszahelyettesítés
for j=4:-1:1
    % j. oszlopbeli elemek kiütése;
    % ugyanezen múvelet elvégzése a hozzácsatolt, utolsó oszlopon
    M(1:(j-1), end) = M(1:(j-1), end) - M(j, end)*M(1:(j-1), j);
    % ismerve a kimenetelt, nem muszáj elvégeznünk a műveletet
    M(1:(j-1), j) = 0;
end

disp('A Gauss(-Jordan) elimináció (ideértve a visszahelyettesítést) eredménye:')
M

%

disp('A beépített parancsokkal kapott megoldások:')
% Beépített Matlab parancs
x = A\b
% reduced row echelon form
rref([A b])

%
% 2. feladat: Adjuk meg az előző feladat mátrixának LU-felbontását!
% Megoldás: A Gauss-elimináció adja meg az U mátrixot, míg az adott lépésekben a szorzótényezők az L elemeit.
% Tétel: Egy adott A mátrixnak pontosan akkor létezik LU-felbontása, ha egyik főminorja se zérus.
% Ellenőrizzük ezt a tulajdonságot!

nincs_kozel_nulla_fominor(A)

ones(2)
nincs_kozel_nulla_fominor(ones(2))

% Maga a felbontás:

% Beépített paranccsal:
help lu
[L_, U_, P] = lu(A)
P'*L_*U_

% U és normált L meghatározása
L=eye(4);
U=A;
for j=1:3
    for i=j+1:4
        m = U(i,j)/U(j,j);
        L(i,j) = m;
        U(i,:) = U(i,:) - m*U(j,:);
    end 
end

L
U
L*U

% Vektorizalt megoldas:

[L, U] = lu_vec(A)

% Ezt is szokás LU-felbontásnak nevezni, az előadáson ugyanakkor az U volt normálva, nem az L. Ezt a problémát úgy tudjuk megoldani, ha bevezetjük az alábbi felbontást:
% A = L D U,
% ahol D egy diagonális mátrix, főátlójában az U főátlójával, az U pedig így normált lesz (az L-el egyetemben). Ezután ha összeszorozzuk a D mátrixot az L mátrixszal, megkapjuk az előadáson tanult formulát.

% Normált U és L meghatározása az L*D*U felírásával

L
D = diag(diag(U))
U = D\U

L*D*U

%
% 3. feladat: Írjuk fel az alábbi mátrix LU- és Cholesky-felbontását!
% A =  \left[\begin{array}{ccc}
% 4 & 2 & 1\\
% 2 & 17 & 13/2 \\
% 1 & 13/2 & 23/2
% \end{array}\right]
% Elméleti összefoglaló:
% Tétel: Tegyük fel, hogy az A
%  mátrix szimmetrikus, pozitív definit. Ekkor létezik pontosan egy olyan pozitív diagonálisú G alsó háromszögmátrix, amelyre A = G G^T. Ezt a felbontást Cholesky-felbontásnak nevezzük.
% A Cholesky-felbontást meghatározhatjuk az LU-felbontásból is, vagy kiszámíthatjuk direkten is.
% Az előző feladat segítségével az LU-felbontás alakja:

A = [4 2       1;
     2 17   13/2; 
     1 13/2 23/2]

[L, U, P] = lu(A)
P
L*U

% Ellenőrizzük, hogy a mátrixnak valóban lesz Cholesky-felbontása! Látható, hogy a mátrix szimmetrikus, a sajátértékeket pedig matlab segítségével kaphatjuk meg:

eig(A)

% Egy másik lehetőség, ha a főminorjait vizsgáljuk meg: ha mindegyik pozitív, akkor rendben vagyunk.

for i=1:3
    det(A(1:i,1:i))
end

fominorokra_igaz(A, @(x)(x > eps))

% A Cholesky-t megkaphatjuk az előzőkben megismert LDU felbontás segítségével úgy, hogy a D minden eleméből gyököt vonok, majd az egyik gyökös mátrixot L-el, míg a másikat U-val szorzom.

sqrtD = diag(sqrt(diag(U)))

G=L*sqrtD

G*G'

% A Cholesky-t direkten is meghatározhatjuk, ha felírjuk az alábbi összefüggést:
% A =  \left[\begin{array}{ccc}
% 4 & 2 & 1\\
% 2 & 17 & 13/2 \\
% 1 & 13/2 & 23/2
% \end{array}\right] =  \left[\begin{array}{ccc}
% a & 0 & 0\\
% b & c & 0 \\
% d & e & f
% \end{array}\right] \cdot  \left[\begin{array}{ccc}
% a & b & c\\
% 0 & d & e \\
% 0 & 0 & f
% \end{array}\right]
% Látható, hogy a bal felső elem miatt a^2=4, azaz a=2
% :
% A =  \left[\begin{array}{ccc}
% 4 & 2 & 1\\
% 2 & 17 & 13/2 \\
% 1 & 13/2 & 23/2
% \end{array}\right] =  \left[\begin{array}{ccc}
% 2 & 0 & 0\\
% b & c & 0 \\
% d & e & f
% \end{array}\right] \cdot  \left[\begin{array}{ccc}
% 2 & b & c\\
% 0 & d & e \\
% 0 & 0 & f
% \end{array}\right]
% Az első sor második elemére 2b=2 azaz b=1 adódik:
% A =  \left[\begin{array}{ccc}
% 4 & 2 & 1\\
% 2 & 17 & 13/2 \\
% 1 & 13/2 & 23/2
% \end{array}\right] =  \left[\begin{array}{ccc}
% 2 & 0 & 0\\
% 1 & c & 0 \\
% d & e & f
% \end{array}\right] \cdot  \left[\begin{array}{ccc}
% 2 & 1 & c\\
% 0 & d & e \\
% 0 & 0 & f
% \end{array}\right]
% Az első sor harmadik eleme miatt pedig 2c=1, azaz c=1/2. Így tovább haladva kapjuk az előzővel megegyező eredményt.
%
% Szorgalmi feladat: Tekintsük az alábbi mátrixot, mely egy mátrix LU-felbontását tartalmazza olyan módon, hogy a főátló "alatt" találhatók a normált L főátló alatti elemei, míg a többi helyeken U értékei szerepelnek:
% A =  \left[\begin{array}{ccc}
% 5 & 4 & 1 \\
% 0.8 & 6.8 & 5.2 \\
% 0.2 & 13/17 & 65/17
% \end{array}\right]
% Létezik-e az eredeti A mátrixnak Cholesky felbontása?
% Ha igen, adjuk meg az ehhez tartozó G mátrixot a gyakorlaton megismert mindkét módszerrel is!
% Adjuk meg az eredeti mátrixot is!

%

% Szorgalmi feladat: Legyen A invertálható mátrix. Mutassuk meg, hogy megfelelő méretű u és v oszlopvektorok esetén teljesül, hogy
% (A + uv^T)^{-1} = A^{-1} - \frac{A^{-1} u v^T A^{-1}}{1 + v^T A^{-1} u.
% Követhetjük az alábbi lépéseket:
% Határozzuk meg az uv^T
% mátrix sajátértékeit és karakterisztikus polinomját.
% Mutassuk meg, hogy \det(I + uv^T) = (1+v^T u ).
% Mutassuk meg, hogy \det(A+ uv^T) = (1 + v^T A^{-1} u)\det A.
% Mutassuk meg, hogy ha A invertálható, akkor A+uv^T
%  is invertálható.
% Mutassuk meg (például bonyolultnak tűnő, de valójában nem horribilis számolással) hogy teljesül a feladat fő állítása.

%

% Szorgalmi feladat:
% Tekintsük az
% A_n = \left(\matrix{3 & 0 & 1 & 1 & \dots & 1\cr 0 & 3 & 0 & 1& \dots & 1 \cr 1 & 0 & 3 & 0 & \dots & 1 \cr \vdots & & \ddots & \ddots & \ddots& \vdots \cr 1 & \dots && 1& 0 & 3}\right)
% n\times n-es mátrix által megadott
% A_n x = b
% egyenletrendszert, ahol b minden eleme 1. Az előző szorgalmi segítségével készítsünk hatékony eljárást ezen egyenletrendszer megoldására, majd kódoljuk is ezt le, végezetül futtassuk n=10^4 esetén.
% Segítség: Milyen B_n mátrixra és u, v vektorokkal lesz A_n = B_n + uv^T, ahol a B_n mátrixszal megadott egyenletrendszer megoldása olcsón kiszámolható?

%

%
% Függvénytemető

function x = cramer_szabaly(A, b)
    [n, m] = size(A);
    [k, l] = size(b);
    assert(n == m && n == k && l == 1)
    
    detA = det(A);
    
    x = zeros(n, 1);
    for i = 1:n
        x(i) = det([A(:, 1:(i-1)), b, A(:, (i+1):n)]) / detA;
    end
end

%

function r = fominorokra_igaz(A, condition)
    [n, m] = size(A);
    assert(n == m)
    
    r = true;

    for i=1:n
        r = condition(det(A(1:i,1:i)));
        if ~r
            return
        end
    end
end

function r = nincs_kozel_nulla_fominor(A)
    r = fominorokra_igaz(A, @(x)(abs(x) > eps));
end

%

function [L, U] = lu_vec(A)
    [n, m] = size(A);
    
    assert(n==m)
    
    L = zeros(n);
    U = zeros(n);
    % A = L*U
    for j=1:n
        % A(j, :) = L(j, :) * U -bol kifejezve az U(j, j:end) reszt
        U(j, j:end) =  A(j, j:end) - L(j, 1:(j-1)) * U(1:(j-1), j:end);
        % A(:, j) = L * U(:, j) -bol kifejezve az L(j:end, j) reszt
        L(j:end, j) = (A(j:end, j) - L(j:end, 1:(j-1)) * U(1:(j-1), j))/U(j,j);
    end
end