%% Ötödik gyakorlat
%% Gradiens módszerek
%
%% Klasszikus gradiens módszer
% Meg szeretnénk oldani az
% Ax = b
% egyenletet, ahol A szimmetrikus, pozitív definit.
% Állítás: Ha x megoldása a fentinek, akkor ez ekvivalens azzal, hogy x a minimuma a \Phi(x)=\frac12 x^T A x - x^T b függvénynek.
% Alapgondolat: Keressük a fenti függvény minimumát úgy, hogy egy adott p irányba elindulva megnézzük, azon az egyenesen hol lesz ez a minimum, aztán innen haladunk tovább.
% Állítás: Egy adott x pontból kiinduló p irányú egyenesen a fenti \Phi függvény minimuma az alábbi pontban található:
% x + \frac{r^T r}{p^T A p} p
% ahol r=b-Ax (a "hiba", "maradék").
% Megjegyzés: Egy lépés után mindig r irányába megyünk tovább, amely mindig merőleges az előző p irányra.
% Ezek alapján maga a gradiens módszer pszeudokódja:
% x = [0, 0, ... 0];
% r = b;
% p = r;
% % k a lépések száma
% for i=1:k
%    alfa=(r^T*r)/(p^T*A*p);
%    x = x + alfa*p;
%    r = b - A*x;
%    p = r;
% end
% 1. feladat:  Tekintsük az alábbi egyenletrendszert:
% \begin{cases}
% 4x + 2y = 7,\\
% 2x+3y=10.
% \end{cases}
% Keressük meg az egyenlet megoldását gradiens módszer segítségével az origóból indulva! Végezzünk el egy lépést papíron, majd programozzuk le a műveletet és így lépjünk többet! Hány lépés után lesz a végtelennormában vett hiba kisebb, mint 10^{-6} (a pontos megoldás x=0.125, y=3.25)?
% Megoldás: Most A = \left[\matrix{4 & 2 \cr 2 & 3}\right] és b = \left[ \matrix{7 \cr 10}\right]. Mivel kezdetben x = \left[ \matrix{0 \cr 0}\right], így r=b és legyen p=r. Az első alfa érték:
% \alpha = \frac{\left[ \matrix{7 &10}\right] \cdot \left[ \matrix{7 \cr 10}\right]}{\left[\matrix{7 &10}\right] \left[\matrix{4 & 2 \cr 2 & 3}\right] \left[ \matrix{7 \cr 10}\right]}
% =
% \frac{49+100}{\left[ \matrix{7 &10}\right] \cdot \left[ \matrix{48 \cr 44}\right]}
% =
% \frac{149}{776}=0,192
% Ezek alapján a lépés a p irányba:
% x_{új}= \left[ \matrix{0 \cr 0}\right] + 0,192 \left[ \matrix{7 \cr 10}\right] =\left[ \matrix{1,344 \cr 1,92}\right]
% A következő irányt kijelölő maradékvektor pedig:
% r_{új} = \left[ \matrix{7 \cr 10}\right] - \left[\matrix{4 & 2 \cr 2 & 3}\right] \left[ \matrix{1,344 \cr 1,92}\right] = \left[ \matrix{-2,217 \cr 1,55}\right]
% A további lépéseket a MATLAB segítségével végezzük el:

A = [4 2; 2 3];
b = [7; 10];

x = [0; 0];
r=b;
p=r;

% Lépésszám
k=0; 
pontos = [0.125; 3.25];
while norm(pontos-x,'Inf')>10^(-6)
    alfa=(r'*r)/(p'*A*p);
    x = x + alfa*p;
    r = b - A*x;
    p=r;
    k=k+1;
end
x
k

%
%% Konjugált gradiens módszer
% Kérdés: hogyan lehetne az új keresési irányt jobban választani?
% Állítás: Egy jobb keresési irány az lesz, ha az új irány r - \beta p, ahol
% \beta = \frac{p^T A r}{p^T A p}.
% Tehát ekkor az új, ún. konjugált gradiens algoritmus:
% x = [0, 0, ... 0];
% r = b;
% p = r;
% % k a lépések száma
% for i=1:k
%    alfa=(r^T*r)/(p^T*A*p);
%    x = x + alfa*p;
%    r = b - A*x;
%    beta = (p^T*A*r)/(p^T*A*p);
%    p = r - beta*p;
% end
% Állítás: Egy n-változós egyenlet esetén a konjugált gradiens módszer véget ér n lépésben.
% 2. feladat: Tekintsük az 1. feladat egyenletrendszerét, és oldjuk meg konjugált gradiens módszerrel! Végezzünk el két lépést kézzel! Mit tapasztalunk? MATLAB kód segítségével ellenőrizzük is le az eredményt!
% Megoldás: Az első feladathoz hasonlóan most A = \left[\matrix{4 & 2 \cr 2 & 3}\right], b = \left[ \matrix{7 \cr 10}\right], x = \left[ \matrix{0 \cr 0}\right], r=b és  p=r, továbbá \alpha = 0.192, azaz az első lépés ugyanaz, és amit végül kapunk:  x_{új}=\left[ \matrix{1,344 \cr 1,92}\right] és r_{új} = \left[ \matrix{-2,217 \cr 1,55}\right].
% A különbség az új keresési irány számításában van:
% \beta = \frac{\left[ \matrix{7 &10}\right] \left[\matrix{4 & 2 \cr 2 & 3}\right] \left[ \matrix{-2,217 \cr 1,55}\right]}{\left[\matrix{7 &10}\right] \left[\matrix{4 & 2 \cr 2 & 3}\right] \left[ \matrix{7 \cr 10}\right]}
% =
% \frac{\left[ \matrix{7 &10}\right] \left[ \matrix{-5,768 \cr 0,216}\right]}{776}=-0.049
% Azaz így az új keresési irány:
% p = \left[ \matrix{-2,217 \cr 1,55}\right] - (-0.049)\left[ \matrix{7 \cr 10}\right]
% =
% \left[ \matrix{-1,874 \cr 2,04}\right] .
% A következő lépés:
% \alpha = \frac{\left[ \matrix{-2,217 &1,55}\right] \cdot \left[ \matrix{-2,217 \cr 1,55}\right]}{\left[\matrix{-1,874 &2,04}\right] \left[\matrix{4 & 2 \cr 2 & 3}\right] \left[ \matrix{1,874 \cr 2,04}\right]}
% =0.651
% x_{új}= \left[ \matrix{1,344 \cr 1,92}\right] + 0,651 \left[ \matrix{-1,874 \cr 2,04}\right] =\left[ \matrix{1,25 \cr 3,25}\right]
% A maradékvektor:
% r_{új}= \left[ \matrix{7 \cr 10}\right] - \left[\matrix{4 & 2 \cr 2 & 3}\right] \left[ \matrix{1,25 \cr 3,25}\right] = \left[ \matrix{0 \cr 0}\right],
% ami azt jelenti, hogy megkaptuk a pontos megoldást (ami szinkronban van a korábban említett állítással).
% Ugyanez MATLAB segítségével:

A = [4 2; 2 3];
b = [7; 10];

x = [0; 0];
r=b;
p=r;

for i=1:2
    alfa = (r'*r)/(p'*A*p);
    x = x + alfa*p;
    r = b - A*x;
    beta = (p'*A*r)/(p'*A*p);
    p=r - beta*p;
end
x
r

% A MATLAB-ban nem pont a nullavektort kaptuk meg, hiszen vannak gépi számítási hibák, de nagyon kicsi a maradékvektor, ezért feltételezhetjük, hogy ez már a valódi megoldás.
%
% Szorgalmi feladat: Tekintsük egy n\times n-es négyzetes, tridiagonális A mátrixot, melynek főátlójában 2-esek, főátlója alatti és feletti mellékátlóiban pedig -1-esek állnak. Tegyük fel, hogy meg szeretnénk oldani egy olyan A x = b egyenletet, ahol b a csupa egyesekből álló vektor.
% (a) Lássuk be, hogy alkalmazható ezen rendszerre a gradiens módszer! Azaz, ellenőrizzük a függvény konvexitását, azaz A pozitív szemidefinitását!
%       Eljárhatnánk a következőképpen, ha ez a megoldás nem lett volna már véletlenül kiadva:
% x^T A x = 2 \sum_{i=1}^n x_i^2 - 2 \sum_{i \neq j} x_i x_j = x_1^2 + \sum_{i=1}^{n-1} (x_{i+1} - x_i)^2 + x_n^2 \geq 0.
% Ehelyett csináljuk a következőt:
% legyen A_n a fenti n\times n-es mátrix, legyen A_0 az üres mátrix, melynek determinánsa legyen (def. szerint) 1. Mutassuk meg, hogy
% az A_n bal felső k\times k-as részmátrixai pont A_k-k; majd
% például sor szerinti kifejtéssel fejezzük ki A_n determinánsát nála kisebb indexű A_j mátrixok determinánsával; majd
% néhány konkrét n-re számoljuk ki \det(A_n)-et, fogalmazzunk meg ezekről egy hipotézist, végül (például indukcióval) bizonyítsuk is azt; majd
% a fentiek segítségével indokoljuk meg, hogy miért alkalmazható a gradiens-módszer.

%

% (b) Írjunk programot, mely adott n értékre visszaadja azt a lépésszámot, amely esetén a módszer hibája kisebb, mint 10^{-6} (a pontos megoldás kiszámításához használjuk az x=b/A parancsot (ahol b az egyenlet jobb oldala))! A lépésszám értékeit ábrázoljuk n függvényében (n \in \{ 3, \dots 20 \})!

%

% Szorgalmi feladat:
% Tegyük fel, hogy az
% Ap = b
% egyenletrendszert szeretnénk iteratív módon megoldani, ahol az alkalmazás által diktált feltétel, hogy a p vektor koordinátái nemnegatívak és egy összegűek.
% a) Ábrázoljuk a háromdimenziós tér azon legszűkebb részét, melynek a fentip eleme. Egydimenzióban egy pontról van szó, kétdimenzióban egy egyenesről. Felhasználhatjuk a linspace, meshgrid és surf parancsokat.

%

% b) Módosítsuk a gradiens-módszert úgy, hogy az
% x_{n+1} = x_n - \alpha_n r_n
% sorra alkalmazzuk koordinátánként az exponenciális függvényt, majd az ezt követő sorban normáljuk le úgy az x_{n+1} vektort, hogy a kapott vektor koordinátáinak összege már 1 legyen. Legyen az így kapott vektor az új közelítő megoldás, x_{n+1}.
% Kódoljuk le ezt az algoritmust.

%

% c) Alkalmazzuk az algoritmust akkor, ha A és b

A = [3, 1, 2;
     1, 4, 1;
     2, 1, 3]/6;

b = [1,1,1]'/3;

% továbbá x_0 kezdővektor

x0 = [1/8, 1/4, 1/2]';

% és addig lépjünk, amíg két szomszédos x és x' iteráltra nem teljesül, hogy
%  \sum\limits_{i=1}^3 x_i \left( \log(x_i) - \log(x_i')\right) \leq 10^{-3},
% ahol x_i most az x iterált i-edik koordinátája.
% A fenti képletben ha x_i = 0, akkor az összeg i-edik tagja 0, azt hagyjuk el. Ha x_i' = 0, akkor a kifejezés értéke \infty.

%