Zagadnienia

5.1 Dyskusja problemu
5.2 Implementacja
5.3 Wykorzystanie specyfiki zadania
5.4 Przypadek dużego $n$

5. Hodowla zwierząt

Aby ocenić wartość hodowlaną $n$ zwierząt na podstawie $m\leq n$ pomiarów wartości pewnej cechy, stosuje się pewien model liniowy, prowadzący do następującego zagadnienia, opisanego w [29, rozdział 5]:

$\displaystyle\begin{pmatrix}X^{T}X&X^{T}Z\\ Z^{T}X&Z^{T}Z+kA^{{-1}}\end{pmatrix}\begin{pmatrix}b\\ a\end{pmatrix}=\begin{pmatrix}X^{T}y\\ Z^{T}y\end{pmatrix}$

(5.1)

Szukane wartości hodowlane zwierząt oznaczone są jako wektor $a\in\mathbb{R}^{n}$ (zatem $i$ -te zwierzę ma wartość $a_{i}$ ). Niewiadomymi pomocniczymi są parametry wpływu płciowości $b\in\mathbb{R}^{2}$ .

Pozostałe parametry występujące w (5.1) — $A$ , $X$ , $Z$ , $y$ , $k$ — są zadane jako parametry modelu. Jak to często się zdarza w przypadku obliczeń naukowych, będziemy dysponować jedynie częściową informacją o danych modelu, a naszym zadaniem będzie po prostu wskazanie sensownej metody numerycznego rozwiązywania układu równań liniowych (5.1). Co zatem na początek wiemy o parametrach modelu? Nasz ,,zleceniodawca” z pewnością zwróci nam uwagę na to, że macierze $X$ oraz $Z$ są macierzami zerojedynkowymi. Macierz $X$ rozmiaru $m\times 2$ określa płeć badanego zwierzęcia, a macierz $Z$ , rozmiaru $m\times n$ , odpowiada za ,,wkład” danego zwierzęcia do badanej cechy hodowlanej.

Wreszcie, o macierzy $A$ rozmiaru $n\times n$ — tzw. macierzy addytywnych pokrewieństw — wiadomo, że ma elementy nieujemne, jest symetryczna i dodatnio określona.

W praktyce [19, Table 1] spotyka się zadania dla $n$ z zakresu od $10^{1}$ do $10^{6}$ . Macierz addytywnych pokrewieństw $A$ w niektórych modelach może być pominięta (co odpowiada wartości parametru skalującego $k=0$ ), a w ogólności ze względu na swoją naturę powinna ona być dosyć rozrzedzona (hodowcy dążą do tego, by ograniczyć pokrewieństwa pomiędzy osobnikami).

Poniżej zacytujemy konkretne zadanie modelowe opisane w [29].

Przykład 5.1 (Przykład 62 z [29])

W oparciu o metodę BLUP wykonać ocenę wartości hodowlanej zwierząt dla cechy masa cieląt przy odsadzeniu w oparciu o następujące informacje:

Cielę	Płeć	Ojciec	Matka	Waga przy odsadzeniu (kg)
4	buhajek	1	-	4,5
5	cieliczka	3	2	2,9
6	cieliczka	1	2	3,9
7	buhajek	4	5	3,5
8	buhajek	3	6	5,0

Zebrane informacje fenotypowe prowadzą do zadania (5.1), w którym

$X=\begin{pmatrix}1&\\ &1\\ &1\\ 1&\\ 1&\\ \end{pmatrix}_{{5\times 2}},\quad Z=\begin{pmatrix}&&&1&&&&\\ &&&&1&&&\\ &&&&&1&&\\ &&&&&&1&\\ &&&&&&&1\\ \end{pmatrix}_{{5\times 8}},\quad y=\begin{pmatrix}4.5\\ 2.9\\ 3.9\\ 3.5\\ 5.0\end{pmatrix},$

parametr $k=2$ , a symetryczna macierz addytywnych pokrewieństw jest zadana przez

$A=\begin{pmatrix}1&&&1/2&&1/2&1/4&1/4\\ &1&&&1/2&1/2&1/4&1/4\\ &&1&&1/2&&1/4&1/2\\ &&&1&&1/4&1/2&1/8\\ &sym&&&1&1/4&1/2&3/8\\ &&&&&1&1/4&1/2\\ &&&&&&1&1/4\\ &&&&&&&1\end{pmatrix}_{{8\times 8}}.$

5.1. Dyskusja problemu

Zadanie wydaje się łatwe do rozwiązania: dane są macierze i parametry, wskazany jest układ $n+2$ równań do rozwiązania (5.1) — więc wystarczy zbudować jego macierz i go rozwiązać. Gdy zadanie nie jest zbyt wielkiego rozmiaru ( $n$ rzędu tysiąca?) — możemy je rozwiązać wprost w Octave. Rzeczywiście, [29, rozdział 10] podaje gotowy skrypt:

G = [X'*X, X'*Z ; Z'*X, Z'*Z + k*inv(A)];
r = [X'*y; Z'*y];
s = inv(G)*r;
b = s(1:2); a = s(3:end);

Przedyskutujmy wady i zalety powyższego rozwiązania. Tym, co od razu kłuje na w oczy, jest używanie funkcji inv do wyznaczania macierzy odwrotnej do $G$ i do $A$ . O ile macierz odwrotna do $A$ występuje w samym sformułowaniu zadania, o tyle wyznaczenie rozwiązania układu równań

$Gs=r$

metodą s = inv(G)*r powinno zjeżyć nam włos na głowie. Oczywiście, choć matematycznie jest to poprawne, w realizacji numerycznej nie powinniśmy wyznaczać wprost macierzy odwrotnej! Znacznie lepiej dokonać rozkładu LDL macierzy $G$ (wszak jest symetryczna!) i następnie rozwiązać dwa układy równań z macierzami trójkątnymi i jedną diagonalną. Ten algorytm realizuje w Octave operator ,,dzielenia” macierzowego:

s = G\r;

Zapatrzeni w odwracanie macierzy, możemy przeoczyć inną niepokojącą cechę układu (5.1): jeśli bowiem $k=0$ , to macierz naszego układu przyjmuje postać:

$\begin{pmatrix}X^{T}X&X^{T}Z\\ Z^{T}X&Z^{T}Z\end{pmatrix}.$

To jest przecież nic innego, jak macierz równań normalnych dla zadania najmniejszych kwadratów

$\| y-\begin{pmatrix}X&Z\end{pmatrix}\begin{pmatrix}b\\ a\end{pmatrix}\| _{2}\rightarrow\min!$

— a zadanie najmniejszych kwadratów, jak wiemy, bezpieczniej rozwiązywać metodami innymi niż poprzez układ równań normalnych: na przykład, opartymi na rozkładzie QR macierzy $\begin{pmatrix}X&Z\end{pmatrix}$ . Najpierw jednak musimy zadać sobie pytanie, czy również oryginalny układ (5.1) odpowiada jakiemuś zadaniu najmniejszych kwadratów? Możemy domyślać się, że tak (wiedząc o tym, jaki jest jego rodowód). I rzeczywiście, przecież dla dowolnej macierzy symetrycznej $S$ ,

$\begin{pmatrix}X&Z\\ 0&S\end{pmatrix}^{T}\cdot\begin{pmatrix}X&Z\\ 0&S\end{pmatrix}=\begin{pmatrix}X^{T}X&X^{T}Z\\ Z^{T}X&Z^{T}Z+S^{2}\end{pmatrix}.$

Biorąc więc $S=k^{{1/2}}A^{{-1/2}}$ (istnieje, bo $A$ jest symetryczna i dodatnio określona) dostajemy, że (5.1) jest układem równań normalnych dla zadania najmniejszych kwadratów:

$\displaystyle\|\begin{pmatrix}y\\ 0\end{pmatrix}-\begin{pmatrix}X&Z\\ 0&S\end{pmatrix}\begin{pmatrix}b\\ a\end{pmatrix}\| _{2}\rightarrow\min!$

(5.2)

Problem z ostatnim sformułowaniem problemu polega na tym, że aby go postawić, musimy wyznaczyć $A^{{-1/2}}$ , czyli — w rzeczywistości — rozwiązać pełne zagadnienie własne dla macierzy $A$ (znaleźć wszystkie wektory i wartości własne).

Jednak po chwili namysłu możemy zobaczyć, że $S$ nie musi być symetryczna, wystarczy tylko, żeby $S^{T}S=kA^{{-1}}$ . Biorąc więc (łatwo obliczalny) rozkład Cholesky'ego–Banachiewicza macierzy $A$ ,

$A=LL^{T},$

dostajemy $S=k^{{1/2}}L^{{-1}}$ . Ponieważ $L$ jest macierzą dolną trójkątną, macierz $L^{{-1}}$ można wyznaczyć przyzwoitym kosztem.

Ostatecznie więc, zamiast zadania (5.1) należy rozwiązać równoważne zadanie najmniejszych kwadratów,

$\displaystyle\|\begin{pmatrix}y\\ 0\end{pmatrix}-\begin{pmatrix}X&Z\\ 0&\sqrt{k}L^{{-1}}\end{pmatrix}\begin{pmatrix}b\\ a\end{pmatrix}\| _{2}\rightarrow\min!$

(5.3)

Jeśli obecność macierzy odwrotnej w powyższym sformułowaniu wciąż nas niepokoi, możemy drążyć dalej. Rozpisując, dostajemy inną postać minimalizowanego wyrażenia (5.3),

$\| y-Xb-Za\| _{2}^{2}+k\| L^{{-1}}a\| _{2}^{2}\rightarrow\min!$

Oznaczając $g=L^{{-1}}a$ , mamy równoważnie

$\| y-Xb-ZLg\| _{2}^{2}+k\| g\| _{2}^{2}\rightarrow\min!$

— a tu już nie występuje macierz odwrotna do $L$ ! Ostatecznie, dostajemy następujące sformułowanie zadania wyjściowego (por. [15, zadanie 11.3.4]:

Wyznacz rozkład Cholesky'ego $A=LL^{T}$ .
Wyznacz rozwiązanie $(b,g)\in\mathbb{R}^{2}\times R^{n}$ zadania najmniejszych kwadratów

$\displaystyle\|\begin{pmatrix}y\\ 0\end{pmatrix}-B\begin{pmatrix}b\\ g\end{pmatrix}\| _{2}\rightarrow\min!$ (5.4)

gdzie $B$ jest prostokątną macierzą rozmiaru $(m+n)\times(2+n)$ postaci

$B=\begin{pmatrix}X&ZL\\ 0&k^{{1/2}}I\end{pmatrix}.$
Oblicz $a=Lg$ .

5.2. Implementacja

Dzięki odpowiedniemu potraktowaniu problemu, całkowicie uniknęliśmy wyznaczania macierzy odwrotnych oraz niepotrzebnego przekształcania zadania najmniejszych kwadratów do postaci normalnej.

Rozkład Cholesky'ego macierzy $A$ wyznaczymy korzystając z funkcji Octave U = chol(A). Jednak musimy pamiętać, że wynikiem działania chol(A) jest macierz trójkątna górna taka, że $U^{T}U=A$ , czyli innymi słowy, $U=L^{T}$ . Dlatego musimy odpowiednio dostosować macierz zadania (5.4):

B = [X Z*U'; zeros(n,2), sqrt(k)*eye(size(U))];

Dalej, wystarczy rozwiązać zadanie najmniejszych kwadratów z macierzą $B$ :

f = [y; zeros(n,1)];
q = B \ f;
b = q(1:2); a = U'*q(3:end);

Przypomnijmy, że w Octave zadanie najmniejszych kwadratów rozwiązuje się, korzystając z tego samego operatora ,,dzielenia”, \, który służy do rozwiązywania układów równań.

Jednak nasze zadanie jest bardzo szczególnej postaci blokowej: blok $(2,2)$ macierzy $B$ jest macierzą diagonalną, co może nam pozwolić osiągnąć dalszą redukcję kosztów obliczeń [15, rozdział 11.3]. Ponadto, można od razu tak ponumerować równania, by macierz $Z$ była postaci

$Z=\begin{pmatrix}0&I\end{pmatrix},$

co spowoduje, że iloczyn $ZL$ będzie wyznaczalny zerowym kosztem obliczeniowym.

5.3. Wykorzystanie specyfiki zadania

Dotychczas atakowaliśmy postawione zadanie przy minimalnej wiedzy o jego naturze. Staraliśmy się przyjąć neutralny punkt widzenia numeryka, pozwalający nam dostrzec w zadaniu pewne typowe cechy samego zadania obliczeniowego. Jednak tym, co naprawdę jest piękne w obliczeniach naukowych jest to, że zadania — choć na swój sposób typowe — mają swoje niuanse, które powodują, że czasem warto zmienić swój punkt widzenia i dopasować używane metody do tego, co więcej wiemy o charakterze zadania!

Jak dotąd, braliśmy pod uwagę jedynie fakt, że macierz $A$ jest z góry zadana, dosyć rzadka, symetryczna i dodatnio określona. Nie chcieliśmy wyznaczać $A^{{-1}}$ wiedząc, jak niezręcznie jest numerycznie korzystać z takiej macierzy.

Przypomnijmy powody:

Wyznaczenie macierzy odwrotnej $A^{{-1}}$ jest procesem bardziej kosztownym niż wyznaczenie współczynników jej rozkładu (np. Cholesky'ego, $A=LL^{T}$ , lub, unikając kosztownych pierwiastków, $A=LDL^{T}$ )
Aby dla danego $y$ wyznaczyć $A^{{-1}}y$ wystarczy rozwiązać dwa układy równań z macierzą $L$ , każdy kosztem co najwyżej $O(n^{2})$ działań (i ewentualnie $D$ , kosztem liniowym). Tak wyznaczone rozwiązanie numeryczne jest rozwiązaniem pewnego zadania sąsiedniego, o ile tylko użylismy dobrego algorytmu wyznaczania rozkładu macierzy.
Nawet jeśli $A$ jest macierzą rzadką, to $A^{{-1}}$ zazwyczaj jest macierzą gęstą.

Tymczasem okazuje się, że spotykane w praktyce modele (5.1) korzystają z macierzy $A$ , która jest tak zwaną macierzą addytywnych pokrewieństw. Mając zadaną taką macierz, możemy wyznaczyć macierz do niej odwrotną, ale na pierwszy rzut oka nie widać w niej jakiejś wyrazistej regularności:

A = [0 0 0 1/2 0 1/2 1/4 1/4;
 0 0 0 0 1/2 1/2 1/4 1/4;
 0 0 0 0 1/2 0   1/4 1/2;
 0 0 0 0 0   1/4 1/2 1/8;
 0 0 0 0 0   1/4 1/2 3/8;
 0 0 0 0 0   0   1/4 1/2;
 0 0 0 0 0   0   0   1/4;
 0 0 0 0 0   0   0   0;];
n = size(A,1);
A = A + A' + eye(n);
inv(A)

Jednak, jeśli sprawdzić w fachowej literaturze (np. w [19]) definicję macierzy $A$ to okaże się, że jej elementy wyznacza się z prostego rekurencyjnego wzoru, który bazuje na znanym z danych hodowlanych rodowodzie każdego zwierzęcia. Spojrzenie na wzór określający $A$ w terminach operacji macierzowych, może słabszych psychicznie zwalić z nóg:

$A=(I-P)^{{-1}}D(I-P)^{{-T}},$

gdzie $P$ jest pewną macierzą o co najwyżej dwóch niezerowych elementach w wierszu! Co więcej, macierz $P$ bardzo łatwo wyznaczyć z danych rodowodowych, natomiast $D$ jest zadaną macierzą diagonalną [19]. Stąd oczywiście dostajemy natychmiast bardzo łatwo wyliczalną macierz odwrotną,

$A^{{-1}}=(I-P)^{{T}}D^{{-1}}(I-P).$

Widzimy więc, że macierz $A^{{-1}}$ nie dość, że jest banalna do wyznaczenia, to od razu jest dana w postaci rozkładu $kA^{{-1}}=S^{T}S$ , w którym co prawda $S=k^{{1/2}}D^{{-1/2}}(I-P)$ nie musi być trójkątna górna, ale za to na pewno jest bardzo rozrzedzona (ma tylko około $2n$ niezerowych elementów.

To odkrycie zmienia nasz punkt widzenia! Jeśli bowiem mamy do dyspozycji rozrzedzoną macierz $P$ i diagonalę $d$ macierzy $D$ , konstrukcja macierzy zadania najmniejszych kwadratów (5.2) upraszcza się do utworzenia bloku $S$ , co możemy zaimplementować w Octave na przykład w poniższy sposób:

B = [X Z; zeros(n,2), spdiag(sqrt(k*d))*(speye(n)-P)];

5.4. Przypadek dużego $n$

Niektóre badania mogą dotyczyć $n=10^{6}$ zwierząt ( $m$ , czyli zbiór danych pomiarowych) jest wtedy zwykle dużo mniejsze), więc każdy sposób na to, by obniżyć koszt pamięciowy i obliczeniowy wyznaczenia rozwiązania jest wart uwagi.

Macierze rzadkie

Przede wszystkim, należy wykorzystać fakt, że macierze $X,Z,A^{{-1}}$ są w praktyce macierzami mocno rozrzedzonymi.

Jeśli więc utworzymy $X,Z,P$ jako macierze rzadkie (poleceniem sparse), to macierz $B$ też będzie rzadka. W przypadku, gdy operator \ zostanie przyłożony do prostokątnej macierzy rzadkiej, spowoduje wywołanie specjalizowanej funkcji bibliotecznej z pakietu CXSPARSE, wykonującej rzadki rozkład QR.

Zmniejszenie rozmiaru zadania przez powrót do układu równań normalnych

Równań normalnych nie musimy się obawiać, gdy uwarunkowanie macierzy $B^{T}B$ jest nieduże. W niektórych przypadkach tak rzeczywiście będzie nawet dla bardzo dużych $n$ . Należy jednak pamiętać, że jeśli $m\ll n$ , to zastąpienie macierzy prostokątnej $B$ rozmiaru $(m+n)\times(2+n)$ macierzą kwadratową $B^{T}B$ rozmiaru $(2+n)\times(2+n)$ nie musi dawać znaczących zysków, zwłaszcza, że $B^{T}B$ będzie mniej rozrzedzona niż $B$ .

Użycie metody iteracyjnej

Gdy $n$ jest na tyle duże, że układ równań normalnych $B^{T}B$ stanowiłby poważne wyzwanie dla metody bezpośredniej, można byłoby wówczas skorzystać z metody iteracyjnej rozwiązywania układu $B^{T}B$ — na przykład, moglibyśmy tu zastosować metodę PCG, jednak wyłącznie wówczas, gdy wyjściowe zadanie jest bardzo dobrze uwarunkowane.

Stosując metodę iteracyjną, można przy okazji uniknąć składania całej wielkiej macierzy $A^{{-1}}$ : wystarczy, zgodnie z powyższym przepisem, przyłożyć macierz do wektora, co możemy zapisać stosunkowo prostą funkcją.

Alternatywą dla rozwiązania układu równań normalnych metodą PCG możne być rozwiązanie (nieco lepiej uwarunkowanego) zadania z (pozornie!) wielką macierzą kwadratową $M$ rozmiaru $2n+m+2$ , o bardzo specjalnej, prostej strukturze:

$M\begin{pmatrix}p\\ q\end{pmatrix}\equiv\begin{pmatrix}\alpha I&B\\ B^{T}&\end{pmatrix}\begin{pmatrix}p\\ q\end{pmatrix}=\begin{pmatrix}y\\ 0\end{pmatrix},$

gdzie $\alpha>0$ jest zadanym parametrem dobranym do zadania.

Rzeczywiście, $(p,q)$ jest rozwiązaniem powyższego równania wtedy i tylko wtedy, gdy $q$ jest rozwiązaniem zadania najmniejszych kwadratów $||Bq-y||_{2}=\min{}$ , a $p=(y-Bq)/\alpha$ jest przeskalowanym wektorem residuum.

Ponieważ macierz $M$ jest symetryczna, ale nie jest dodatnio określona, należy zastosować do niej metodę PCR, dostępną m.in. w Octave. Prostoduszna implementacja

M = [alpha*speye(size(B,1)), B; B', spalloc(size(B,2),size(B,2))];
[X,info,relres] = pcr(M,[y;zeros(size(B,2),1)]);
info
q = X(size(B,1)+1:end);

może efektywnością ustąpić miejsca bardziej wyrafinowanej, korzystającej z operatorowej definicji $M$ :

Mmult = @(x) [alpha*x(1:size(B,1)) + B*x(size(B,1)+1:end); B'*x(1:size(B,1))];
[X,info,relres] = pcr(Mmult,[y;zeros(size(B,2),1)]);
info
q = X(size(B,1)+1:end);

Koszt mnożenia wektora przez $M$ możemy jeszcze bardziej zredukować, korzystając z wiedzy o strukturze $B$ : jest ona macierzą blokową trójkątną górną, a i bloki mają specyficzną strukturę, upraszczającą mnożenie przez wektor.

Obniżenie kosztu mnożenia przez $A^{{-1}}$

Możemy też próbować mnożenie przez $A^{{-1}}$ uczynić tańszym, zwłaszcza, gdy implementujemy nasz program w języku C lub podobnym. Z pomocą przychodzi tu doświadczenie… programowania metody elementu skończonego (technika stosowana w numerycznych obliczeniach inżynierskich!). Ponieważ w wierszu $P$ odpowiadającym $i$ -temu zwierzęciu znajdują się co najwyżej dwie niezerowe wartości:

$p_{{i,s(i)}}=p_{{i,d(i)}}=1/2,\qquad d_{i}=1/2,$

gdy ojcem $i$ jest zwierzę o numerze $s(i)$ , a matką — zwierzę o numerze $d(i)$ , a w przypadku, gdy znane jest tylko jedno z rodziców, $r(i)$ , tylko jeden element jest niezerowy,

$p_{{i,r(i)}}=1/2,\qquad d_{i}=3/4,$

a gdy żadne z rodziców zwierzęcia $i$ nie jest znane, cały $i$ -ty wiersz $P$ jest zerowy, natomiast $d_{i}=1$ .

Stąd wynika, że $A^{{-1}}=\sum _{{i}}A_{i}^{{-1}}$ , gdzie $A_{i}^{{-1}}$ jest wkładem $i$ -tego zwierzęcia. Na przykład, gdy znani są oboje rodzice zwierzęcia,

$A_{i}^{{-1}}=\begin{pmatrix}1\\ -1/2\\ -1/2\end{pmatrix}d_{i}^{{-1}}\begin{pmatrix}1&-1/2&-1/2\end{pmatrix},$

przy czym elementy tej macierzy należy rozrzucić na trójkę $i,s(i),d(i)$ współrzędnych $A^{{-1}}$ .

Ćwiczenie 5.1

Niech będzie dana tabela $R$ , rozmiaru $n\times 2$ , określająca bezpośrednie relacje pomiędzy zwierzętami:

$r_{{i,1}}=\begin{cases}s(i),\text{ gdy ojcem $i$-tego zwierzęcia był osobnik o numerze $s(i)$},\\ 0.\end{cases}$

Podobnie określamy $r_{{i,2}}$ , identyfikator matki zwierzęcia $i$ . Niech $A$ rozmiaru $n\times n$ będzie macierzą addytywnych pokrewieństw pomiędzy zwierzętami wyznaczoną przez $R$ . Napisz możliwie szybko działającą procedurę, obliczającą $A^{{-1}}y$ na zadanym wektorze $y$ przy minimalnym zapotrzebowaniu na pamięć roboczą.

Jeśli chcesz zbliżyć się do realnych warunków pracy, przyjmij założenie, że dane z tablicy $R$ są zapisane w arkuszu kalkulacyjnym Excela (lub w jakimś bardziej egzotycznym formacie): konwersja formatów danych to jeden z pomijanych, a bardzo uciążliwych programistycznie, aspektów obliczeń naukowych.

Ćwiczenie 5.2

W bardzo licznych zastosowaniach statystycznych należy wskazać zestaw tych wartości własnych (i, przy okazji, wektorów własnych) macierzy kowariancji $X^{T}X$ , które są powyżej zadanego progu $\eta^{2}$ . To zadanie nazywa się analizą głównych składowych (principal component analysis, PCA), a celem może być zmniejszenie rozmiaru zbioru danych statystycznych poprzez eliminację mało istotnych parametrów.

W sformułowaniu matematycznym, mając zadaną macierz $X\in\mathbb{R}^{{n\times m}}$ — przy czym $n\geq m$ — musimy wyznaczyć rozkład spektralny macierzy (symetrycznej i nieujemnie określonej) $X^{T}X$ :

$X^{T}X=V\Lambda V^{T},$

gdzie $\Lambda$ jest macierzą diagonalną, zawierającą wartości własne, a $V$ jest macierzą ortogonalną, $V^{T}V=I$ , złożoną z wektorów własnych macierzy $X^{T}X$ , a następnie zwrócić te wektory własne $v_{i}$ , które spełniają warunek $\lambda _{i}\geq\eta^{2}$ .

Napisz funkcję, która wykona to zadanie.

Rozwiązanie:

Funkcja Octave, realizująca nasze zadanie mogłaby mieć postać:

function [V, L] = pca1(X, eta)
[V, L] = eig(X'*X);
[L, I] = sort(diag(L), 'descend'); V = V(:,I);
I = find(L>eta);
V = V(:,I); L = L(I);
end

Dodatkowo, nasza funkcja sortuje wektory i wartości własne w kolejności od największej, do najmniejszej.

Jednak ma ona tę wadę, że formując macierz $X^{T}X$ tracimy informację zawartą w $X$ : macierz $X^{T}X$ jest wymiaru tylko $m\times m$ . Dlatego bezpieczniej skorzystać z rozkładu SVD¹³Więcej o rozkładzie według wartości szczególnych (SVD) dowiesz się z wykładu z Matematyki Obliczeniowej II. macierzy $X$ :

$X=U\Sigma V^{T}$

i faktu, że $X^{T}X=V\Sigma^{2}V^{T}$ . Ponieważ funkcja svd zwraca macierze $U,\Sigma,V$ , lepsza numerycznie wersja naszej funkcji miałaby następującą postać:

function [V, L] = pca(X, eta)
[U L V] = svd(X,0); % economy-version
[L, I]  = sort(diag(L).^2, 'descend'); V = V(:,I);
I = find(L>eta);
V = V(:,I); L = L(I);
end

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Obliczenia naukowe