%% Hatodik gyakorlat
%% Sajátértékek numerikus közelítése
% Feladat: Meg szeretnénk határozni (közelítőleg) egy valós mátrix sajátértékeit.
%
% Emlékeztető:
% Def.: A és B négyzetes mátrixok hasonlók, ha létezik egy olyan S invertálható mátrix, hogy
% B = S^{-1} A S.
% Def.: Egy mátrix diagonizálható, ha hasonló egy diagonális mátrixhoz.
% Állítás: Ha az A valós mátrix szimmetrikus, akkor diagonizálható ortogonális mátrixszal, és
% A = V D V^{-1} = VDV^T,
% ahol D egy diagonális mátrix az A sajátértékeivel, V oszlopai pedig az A mátrix sajátvektorai.
% Állítás: Ha egy n \times n-es valós mátrixnak n  darab különböző sajátértéke van, akkor az diagonizálható.
% Miért?: A teljes n-dimenziós teret kifeszíti az n darab (lineárisan független) sajátvektor.
% Állítás: Ha egy n \times n-es valós mátrix diagonizálható, abból nem következik, hogy n  darab különböző sajátértéke van (pl. az identitásmátrix).
% Megjegyzés: Ha egy B mátrix ortogonális, akkor \mathrm{cond}_2(B)=1.
% Biz.: Ha B ortogonális, akkor B B^T=B^T B = I, tehát
% \mathrm{cond}_2(B)=\Vert B \Vert_2 \Vert B^{-1}\Vert_2 =\sqrt{\lambda_{max}(B^T B)} \sqrt{\lambda_{max}(B B^T)} = 1.
%
%% Bauer-Fike tétel
%
% Tétel (Bauer-Fike): Legyen A diagonizálható négyzetes mátrix, azaz amely felírható A = V D V^{-1} alakban (ahol legyen V a normált sajátvektorokból álló mátrix), és jelöljük A
%  sajátértékeit \lambda_i-vel. Legyen továbbá \delta A egy tetszőleges másik mátrix, és legyen \mu az A + \delta A mátrix egy sajátértéke. Ekkor
% \min_{i} |\mu - \lambda_i| \leq \mathrm{cond}_p(V) \Vert \delta A\Vert_p.
% Picit átfogalmazva: HaAdiagonalizálható -azaz van olyan V sajátvektoraiból álló mátrix, mellyel a V^{-1}AV mátrix diagonális- akkor egy B mátrix bármely sajátértékéhez van A-nak legfeljebb \text{cond}_p(V) \|B-A\|_p távolságra lévő sajátértéke.
% Állítás: Ha a fenti tételben A normális (tehát V^{-1} = V^T), akkor p=2 választással a távolság \|B-A\|_2-ra egyszerűsödik -hiszen V ortogonális, így \text{cond}_2(V)=1.
% Gondoljuk meg, hogy ha A = 0 ésB diagonális, akkor B-nek természetesen van sajátértéke a nullától legfeljebb\|B\|_2  = |\lambda_{max}(B)|távolságra, de ez az állítás általában nem túl erős - persze ha |\lambda_\min(B)| = |\lambda_\max(B)|, akkor ennél erősebb nem is lehetne.
%
% 1. feladat: Tekintsük az alábbi mátrixot:
% A = \left(\matrix{2 & 1 \cr 1 & 2}\right).
% a) A Bauer-Fike tétel segítségével és p=2 választással becsüljük meg, hogy mennyivel változhatnak ezen mátrix sajátértékei, ha a mátrix bal alsó eleme megnő 0.1-gyel!
% b) Nézzük meg MATLAB segítségével, hogy az új sajátértékek mennyiben térnek el a régiektől.
% Megoldás: Ez szimmetrikus mátrix, így diagonizálható ortogonális mátrix segítségével.
%     Ezek alapján tehát \mathrm{cond}_2 (V) = 1, és így a becslés:
% \min_{i} |\mu - \lambda_i| \leq 1 \cdot \left\Vert \left( \matrix{0 & 0 \cr 0.1 & 0 }\right)\right\Vert_2.
%     A kettes normához szükségünk lesz a (\delta A)^T \delta A  sajátértékére:
% (\delta A)^T \delta A = \left( \matrix{0 & 0.1 \cr 0 & 0 }\right) \left( \matrix{0 & 0 \cr 0.1 & 0 }\right) =\left( \matrix{0.01 & 0 \cr 0 & 0 }\right).
%     Ez egy diagonális mátrix, így a két sajátértéke 0 és 0.01, azaz a spektrálsugara 0.01. Ezek alapján tehát:
% \left\Vert \left( \matrix{0 & 0 \cr 0.1 & 0 }\right)\right\Vert_2 = \sqrt{0.01} = 0.1.
%     Így a mátrix sajátértékei maximum 0.1-gyel változhatnak - a változhatnak szó fenti tételbeli értelmében.
%     Kód segítségével ellenőrizzük le a kapott eredményt!

lams1 = eig([2 1;
             1 2])

lams2 = eig([2   1;
             1.1 2])

abs(sort(lams1) - sort(lams2))

%
% 2. feladat: Tekintsük az alábbi mátrixot!
% A = \left(\matrix{2 & 9 \cr 1 & 2}\right).
% a) A Bauer-Fike tétel segítségével és p=\infty választással becsüljük meg, hogy mennyivel változhatnak ezen mátrix sajátértékei, ha a mátrix bal alsó eleme megnő 0.1-gyel!
% b) Nézzük meg MATLAB segítségével, hogy az új sajátértékek mennyiben térnek el a régiektől.
% Megoldás: Ez a mátrix nem szimmetrikus, de ettől még lehet diagonizálható, például ha a sajátértékei különböző (valós) értékek.
% A sajátértékeket és a sajátvektorokból álló V mátrixot meghatározhatjuk kézzel, vagy kóddal is:

A = [2 9;
     1 2];

[V, lam] = eig(A)
V ./ (V(2, :))

% Tehát s_1 = \left( \matrix{3 \cr 1} \right), \lambda_1 = 5 és s_2 = \left( \matrix{-3 \cr 1} \right), \lambda_1 = -1. Ezek alapján V = \left( \matrix{-3 & 3 \cr 1 & 1} \right).
%  A becsléshez szükségünk lesz még ezen mátrix végtelen normában számított kondíciószámára is. Emlékeztető: a végtelen norma a sorok abszolútértékösszegeinek maximuma. Mivel V^{-1}= \frac{1}{-6} \left( \matrix{1 & -3 \cr 1 & -3} \right), így
% \mathrm{cond}_{\infty} (V) = \Vert V \Vert_{\infty}  \Vert V^{-1} \Vert_{\infty} = 6 \cdot \frac{4}{6}=4.
% A becsléshez szükségünk van még a változás mátrixának végtelen normájára, de ez könnyen láthatóan \left\Vert \left( \matrix{0 & 0 \cr 0.1 & 0 }\right)\right\Vert_{\infty} = 0.1. Ezek alapján a becslés:
% \min_{i} |\mu - \lambda_i| \leq 4 \cdot 0.1 = 0.4.
% Az eredményt kóddal is leellenőrizhetjük:

abs(sort(eig([2 9; 1 2])) - sort(eig([2 9; 1.1 2])))

%
%% Hatványmódszer
% Legyen A valós, diagonizálható mátrix, és tegyük fel, hogy van egy egyszeresen domináns sajátértéke. Tekintsük az alábbi algoritmust:
% for k=1:kmax
%     x(k) = A*y(k-1);
%     y(k) = x(k)/norm(x(k),2);
%     lambda(k)=(y(k))^T*A*y(k);
% end
%
% % Mivel általában nem érdekelnek minket a múltbeli iterációs lépések, így az alábbi alak is használható:
%
% % x_0 adott
% y = x_0 / norm(y_0)
% for k=1:kmax
%     y = A*x;
%     y = y / norm(y);
% end
%
% % a sajátértéket pedig ráérünk a végén megközelíteni, kivéve például akkor,
% % ha leállási feltételként a közelítésének a lépésenkénti megváltozásást
% % szeretnénk használni
% lambda = y' * A * y;
% Állítás: Ha a kezdeti vektor nem merőleges a domináns sajátértékhez tartozó sajátvektorhoz, akkor a fenti algoritmusban szereplő lambda a domináns sajátértékhez, az y pedig a megfelelő (normált) sajátvektorhoz tart. Itt érdemes lehet megjegyezni, hogy a normált sajátvektor még egyszeres algebrai multiplicitás esetén is csak egységhosszú szorzó erejéig egyértelmű.
% 3. feladat: Tekintsük az alábbi mátrixot!
% A = \left( \matrix{2 & 1 & 0 \cr 1 & 2 & 1 \cr 0 & 1 & 1} \right)
% Vizsgáljuk meg, hogy alkalmazható-e a hatványmódszer, és ha igen, akkor határozzuk meg a domináns sajátértéket legfeljebb \varepsilon = 10^{-4} hibával abban az értelemben, hogy az algoritmust addig futtassuk, amíg a sajátérték két egymásutáni közelítésének távolsága ezen \epsilon érték alá nem csökken!
% Legyen a kezdővektorunk \left( \matrix{6 \cr 7 \cr 3} \right). Mennyit kellett lépnünk? És ha az  \left( \matrix{1 \cr 1 \cr 1} \right) vektorból indulunk?
% Megoldás: A mátrix szimmetrikus, így diagonizálható, azaz alkalmazható a hatványmódszer.
% A kódot úgy írjuk meg, hogy addig fusson, míg az egymást követő lépésekben kapott sajátértékek különbsége nagyobb, mint 10^{-4}.

format long

y = [6 7 3]';
%y = [1 1 1]';
A = [2 1 0; 
     1 2 1; 
     0 1 1];

yregi=[0 0 0]';

lambda = 0;
lambdaregi= 3;

k=0; % Lépésszám

while abs(lambda - lambdaregi)>10^(-4)
    x = A*y;
    y = x/norm(x,2);
    lambdaregi=lambda;
    lambda = y'*A*y
    k = k + 1; 
end
lambda
y
k

% Ellenőrzés:

format long

A = [2 1 0; 
     1 2 1; 
     0 1 1];
eig(A)

% Használhatjuk a függvénytemetőbeli kódot is:

[lambda, y, k] = hatvanymodszer(A, [6 7 3]', 1e-4)

%
%% Inverz iteráció
% Legyen A valós, szimmetrikus mátrix, és tegyük fel, hogy van egy egyszeres sajátértéke (legyen ez \lambda), mely legközelebb van egy adott \mu értékhez.
% Ekkor A - \mu I sajátértékei \lambda_i - \mu alakúak, ahol \lambda_i A egy sajátértéke. Ezek között \lambda - \mu lesz a legkisebb abszolútértékű. A sajátvektorok változatlanok.
% Következésképpen (A-\mu I)^{-1} sajátértékei között az \frac1{\lambda - \mu} dominál abszolútértékben. A sajátvektorok változatlanok.
% Tehát az (A-\mu I)^{-1} mátrixra alkalmazható a hatványmódszer. Tkp. ezt nevezzük inverz iterációnak.
% Ez az eljárás tehát (közelítőleg) meghatározza a (A-\mu I)^{-1} domináns sajátértékéhez tartozó sajátvektort, mely persze ugyanaz mint a \lambda értékhez tartozó sajátvektor, amivel tehát a \lambda értéket is megközelíthetjük.
% Tekintsük az alábbi algoritmust, ahol A a fenti mátrix, x és y pedig vektorokból álló sorozatok, Id pedig az identitásmátrix:
% for k=1:kmax
%     x(k)=(A-mu*Id)^(-1)*y(k-1);
%     y(k)=x(k)/norm(x(k),2);
%     lambda(k)=(y(k))^T*A*y(k);
% end
%
% % illetve a hatványmódszer pszeudókódjánál szereplő megjegyzések itt is érvényesek
% Állítás: A fenti algoritmusban szereplő lambda tart az adott \mu számhoz legközelebbi sajátértékhez.
%
% 4. feladat: Tekintsük újfent a 3. feladat mátrixát:
% A = \left( \matrix{2 & 1 & 0 \cr 1 & 2 & 1 \cr 0 & 1 & 1} \right)
% Határozzuk meg a mátrix \mu=1 értékhez legközelebbi sajátértékét az inverz iteráció segítségével legfeljebb \varepsilon = 10^{-4} hibával, abban az értelemben, hogy az algoritmust addig futtassuk, amíg a sajátérték két egymásutáni közelítésének távolsága ezen \epsilon érték alá nem csökken!
%  Legyen a kezdővektorunk y_0 = \left( \matrix{3 \cr -3 \cr 5} \right). Mennyit kellett lépnünk?
% Megoldás: Az inverz iterációhoz szükségünk lesz az A - \mu Id = \left( \matrix{1 & 1 & 0 \cr 1 & 1 & 1 \cr 0 & 1 & 0}\right) mátrix inverzére, de ezt a MATLAB is ki tudja számítani.

format long

y = [3 -3 5]';
A = [2 1 0; 1 2 1; 0 1 1];
Ainv=inv(A-eye(3));

yregi=[0 0 0]';

lambda = 0;
lambdaregi= 3;

k=0; % Lépésszám

while abs(lambda - lambdaregi)>10^(-4)
    x = Ainv*y;
    y = x/norm(x,2);
    lambdaregi=lambda;
    lambda = y'*A*y
    k = k + 1; 
end
lambda
y
k

% Ellenőrzés:

format long

A = [2 1 0; 1 2 1; 0 1 1];
eig(A)

% Használhatjuk a függvénytemetőbeli kódot is:

[lam, v, num_steps] = inverz_iteracio(A, [3, -3, 5]', 1, 1e-4)

%
% Szorgalmi feladat: Tekintsük egy olyan n \times n-es valós mátrixot, melynek főátlójában n, míg a többi helyén 1 szerepel! Keressük ezen mátrix 1-hez legközelebbi sajátértékét a fentiekben bemutatott módon, majd olyan módon is, hogy az inverzmátrix direkt kiszámítása helyett először (a MATLAB beépített parancsának segítségével) meghatározzuk a mátrix LU-felbontását, s minden lépésben az inverzmátrix helyett U^{-1} L^{-1}-el szorzunk! Az iterációban szereplő kezdővektor minden esetben legyen a csupa 1 vektor, és akkor állítsuk le az iterációt, ha a hiba legfeljebb \varepsilon = 10^{-4}.
% Írjunk egy olyan programot, melynek bemenete az n értéke, kimenete pedig kettő érték: egyrészt a korábbiakban bemutatott módon, inverzmátrixszal való szorzást használó iteráció futási ideje, valamint az LU-felbontást tartalmazó iteráció futási ideje! Az idő méréséhez használjuk a tic és toc parancsokat!

%

% Árbárzoljuk a futásidőket n függvényében n=3-tól n=400-ig. Mi olvasható le az ábráról?
%

%

% Függvénytemető:

A = diag([1,2,3]);
[v, num_steps] = pow_it(A, ones(3,1), @(x, y)(abs(x'*A*x - y'*A*y) < 1e-4));
% diag([1,2,3]) domináns sajátértéke (és a hozzátartozó s.vektor)
[lam, v, num_steps] = hatvanymodszer(diag([1,2,3]), ones(3,1), 1e-4)

% 1.6-hoz legközelebbi s. értéke diag([1,2,3])-nak
[lam, v, num_steps] = inverz_iteracio(diag([1,2,3]), ones(3,1), 1.6, 1e-4)

%

function [y, num_steps] = pow_it(A, y0, leallhat)
    y = y0/norm(y0);
    num_steps = 0;
    
    while true
        y_prev = y;
        y = A * y;
        y = y / norm(y);
        num_steps = num_steps + 1;

        if leallhat(y, y_prev)
            break;
        end
    end
end

function [lambda, v, num_steps] = hatvanymodszer(A, y0, atol)
    [v, num_steps] = pow_it(A, y0,  @(x, y)(abs(x'*A*x - y'*A*y) < atol));
    lambda = v' * A * v;
end

function [lambda, v, num_steps] = inverz_iteracio(A, y0, mu, atol)
    [n, m] = size(A);
    assert(n==m)
    B = inv(A-mu*eye(n));
    % B dominans sajatertekehez tartozo sajatvektorat kozelitjuk meg,
    % mely A-nak is sajátvektora, de leállási feltételként azt figyeljük,
    % hogy A-nak a megfelelő sajátértékének a közelítése mennyivel változik
    % meg egyik iterációs lépésről a másikra.
    [v, num_steps] = pow_it(B, y0,  @(x, y)(abs(x'*A*x - y'*A*y) < atol));
    lambda = v'*A*v;
end