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ń

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

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

— 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 XZ. 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,

Biorąc więc S=k1/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:

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 ST⁢S=k⁢A-1. Biorąc więc (łatwo obliczalny) rozkład Cholesky'ego–Banachiewicza macierzy A,

dostajemy S=k1/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,

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),

Oznaczając g=L-1⁢a, mamy równoważnie

— 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=L⁢LT.

2. Wyznacz rozwiązanie b,g∈R2×Rn zadania najmniejszych kwadratów

   y0-B⁢bg2→min!

   (5.4)

   gdzie B jest prostokątną macierzą rozmiaru m+n×2+n postaci

   B=XZ⁢L0k1/2⁢I.

3. Oblicz a=L⁢g.

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 UT⁢U=A, czyli innymi słowy, U=LT. 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

co spowoduje, że iloczyn Z⁢L będzie wyznaczalny zerowym kosztem obliczeniowym.

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=L⁢LT, lub, unikając kosztownych pierwiastków, A=L⁢D⁢LT)

• 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⁢n2 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:

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ą,

Widzimy więc, że macierz A-1 nie dość, że jest banalna do wyznaczenia, to od razu jest dana w postaci rozkładu k⁢A-1=ST⁢S, w którym co prawda S=k1/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 2⁢n 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)];

Niektóre badania mogą dotyczyć n=106 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.