Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 63 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 65 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 67 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 69 Notice: Undefined variable: base in /home/misc/mst/public_html/lecture.php on line 36 Obliczenia naukowe – 5. Hodowla zwierząt – MIM UW

Zagadnienia

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]:

  1. Wyznacz rozkład Cholesky'ego A=LL^{T}.

  2. 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}.
  3. 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 SVD13Wię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.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.