Show code cell source Hide code cell source

import numpy as np;
import scipy as sc;
import matplotlib.pyplot as plt;
from scipy.linalg import lu, solve_triangular

%matplotlib inline

2.1. Gaussove eliminacije i LU faktorizacija

Gaussove eliminacije

Opis problema. Problem rješavanja sustava linearnih jednadžbi terminima linearne algebre može se opisati kao problem u kojem za zadanu matricu \(A \in \R^{n \times n}\) realnih brojeva i vektor desne strane \(b\in \mathbb R^n\) treba pronaći vektor \(x \in \mathbb R^n\) takav da vrijedi \(A x = b\). Iako teorija linearne algebre pokriva i općenitije slučajeve (mnogi od njih mogu se pokriti i algoritmima numeričke matematike), za potrebe ovog kolegija fokusiramo se na probleme kod kojih pretpostavljamo da je A regularna kvadratna matrica. U tom slučaju, sustav \(Ax = b\) ima jedinstveno rješenje za svaku desnu stranu \(b\).

Kako rješavamo sustav linearnih jednadžbi, i kako bismo ga implementirali u računalu?

• Jedan od načina je odrediti inverz matrice \(A^{-1}\), jer se tada rješenje dobiva po formuli \(x=A^{-1}b\) (na kraju treba pomnožiti matricu vektorom, što je trivijalno). To je teži problem od onog s kojim smo krenuli budući da je \(i\)-ti stupac matrice \(A^{-1}\) rješenje \(x_{e_i}\) sustava \(Ax_{e_i} = e_{i}\), gdje je \(e_i\) \(i\)-ti vektor kanonske baze.

• Drugi način je koristeći Cramerovo pravilo: \(i\)-ta komponenta (označimo ju s \(x_i\)) rješenja sustava \(x\) računa se po formuli \(x_i = \dfrac{\det A_i}{\det A}\). Ovdje je matrica \(A_i\) dobivena tako da je umjesto \(i\)-tog stupca matrice \(A\) zapisan vektor \(b\). Ponovno imamo eksplicitnu formulu za rješenje, samo treba izračunati determinantu od ukupno \(n+1\) matrica. Pomoću Laplaceovog razvoja možemo izračunati determinantu proizvoljne matrice, no složenost takvog računanja je eksponencijalna: problem određivanja determinante matrice reda \(n\) svodi se na problem određivanja determinanti \(n\) matrica reda \(n-1\), dakle \(\bigO(n!)\) operacija, što je ekstremno složeno. Zato odbacujemo i ovaj način rješavanja.

• Treći način je onaj kojim smo rješavali sustava na kolegijma Linearna algebra 1 i 2: Gaussove eliminacije.

Primjer 2.1

Gaussovim eliminacijama riješimo sustav \(Ax = b\), gdje je

\[\begin{split}A = \begin{bmatrix} 1 & 1 & 1 \\ 2 & 4 & 2 \\ -1 & 5 & -4 \end{bmatrix}, \quad b = \begin{bmatrix} 6 \\ 16 \\ -3 \end{bmatrix}.\end{split}\]

Tom sustavu odgovara proširena matrica \(\hat{A}^{(1)} = \left[ A \mid b \right] \) sustava, i ona glasi

\[\begin{split}\hat{A}^{(1)} = \begin{bmatrix} \begin{array}{ccc:c} 1 & 1 & 1 & 6 \\ 2 & 4 & 2 & 16 \\ -1 & 5 & -4 & -3 \end{array} \end{bmatrix}.\end{split}\]

U \(i\)-tom koraku Gaussovih eliminacija promatrat ćemo \(i\)-ti stupac matrice \(\hat{A}^{(i)}\) te pomoću elementa na dijagonali \(a_{i, i}^{(i)}\) poništiti (postaviti na nulu) sve elemente ispod dijagonale u \(i\)-tom stupcu.

U prvom koraku, elementom \(a_{1,1}^{(1)} = 1\) poništimo preostale elemente u prvom stupcu koji se nalaze ispod tog elementa. Prvi redak ne mijenjamo, drugi redak ćemo umanjiti za prvi redak pomnožen s \(m_{2}^{(1)}:= \dfrac{a_{2,1}^{(1)}}{a_{1,1}^{(1)}} = 2\), a treći redak ćemo umanjiti za prvi redak pomnožen s \(m_{3}^{(1)}:= \dfrac{a_{3,1}^{(1)}}{a_{1,1}^{(1)}} = -1\). Izvedene transformacije (dodavanje/oduzimanje jednadžbi pomnoženih skalarom različitim od nule) ne mijenjaju skup rješenja sustava. Te transformacije moramo raditi na proširenoj matrici.

Novodobivena proširena matrica glasi

\[\begin{split}\hat{A}^{(2)} = \begin{bmatrix} \begin{array}{ccc:c} 1 & 1 & 1 & 6 \\ 0 & 2 & 0 & 4 \\ 0 & 6 & -3 & 3 \end{array} \end{bmatrix}.\end{split}\]

Sada prelazimo na drugi korak Gaussovih eliminacija. Prvi redak sada više ne diramo. Elementom \(a^{(2)}_{2,2} = 2\) poništavamo preostale retke u drugom stupcu koji se nalaze ispod tog elementa. To znači da ćemo samo treći redak umanjiti za drugi redak pomnožen s \(m_{3}^{(2)}:=\dfrac{a^{(2)}_{3,2}}{a^{(2)}_{2,2}} = 3\). Nova proširena matrica sustava glasi

\[\begin{split}\hat{A}^{(3)} = \begin{bmatrix} \begin{array}{ccc:c} 1 & 1 & 1 & 6 \\ 0 & 2 & 0 & 4 \\ 0 & 0 & -3 & -9 \end{array} \end{bmatrix}.\end{split}\]

Algoritam nam garantira da matrica sustava koju dobivamo na kraju bude gornjetrokutasta matrica:

\[\begin{split}U := \begin{bmatrix} 1 & 1 & 1 \\ 0 & 2 & 0 \\ 0 & 0 & -3 \end{bmatrix}.\end{split}\]

Sustave s gornjetrokutastim matricama je lako riješiti supstitucijom unatrag. U zadnjem retku odgovarajuća je jednadžba s jednom nepoznanicom: \(-3x_3 = -9\), odakle je \(x_3=3\). S tim rješenjem ulazimo u predzadnji redak (tj., ušli bismo da element \(u_{2,3}\) nije jednak nuli). Odgovarajuća je jednadžba \(2x_2 = 4\), odakle je \(x_2=2\). Konačno, treba odrediti preostalu nepoznanicu iz prve jednadžbe \(x_1+x_2+x_3 = 6\), tj.

\[x_1 = 6 - x_2 - x_3 = 1.\]

Konačno rješenje sustava je \(x = (1,2,3)^{T}\).

Općenito, algoritam za rješavanje sustava Gaussovim eliminacijama glasi ovako:

# Rješavanje sustava Ax=b Gaussovim eliminacijama.
def GaussoveEliminacije(A, b):
    n = A.shape[0]

    # Gaussove eliminacije.
    for i in range(n):
        # Promatramo i-ti stupac matrice A.
        for j in range(i+1, n):
            # Multiplikator: poništavamo element A[j, i].
            m = A[j, i] / A[i, i]

            # Od svakog elementa u j-tom retku oduzmemo m*pripadni element u i-tom retku.
            A[j, :] -= m * A[i, :]
            b[j] -= m * b[i]

    # Rješavanje trokutastog sustava: supstitucija unatrag.
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = b[i]
        for j in range(i+1, n):
            x[i] -= A[i, j]*x[j]
        x[i] /= A[i, i]

    return x

Kod možemo provjeriti na gornjem primjeru:

A = np.array([[1., 1., 1.], [2., 4., 2.], [-1., 5., -4.]])
b = np.array([6., 16., -3.])

x = GaussoveEliminacije(A, b)
print( f'Rješenje danog sustava je {x}.' );

Rješenje danog sustava je [1. 2. 3.].

Gaussove eliminacije s parcijalnim pivotiranjem

Postavlja se pitanje možemo li ovaj algoritam primijeniti svaki put prilikom rješavanja sustava \(Ax=b\), gdje je \(A\) proizvoljna regularna kvadratna matrica. Odgovor je negativan, budući da se može dogoditi da je element \(a^{(i)}_{i,i}\) u nekoj matrici \(\hat{A}^{(i)}\) jednak nuli. Dapače, to se može dogoditi odmah na početku. Sustav

(2.1)\[\begin{split}\begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \end{bmatrix}\end{split}\]

očito ima jedinstveno rješenje \(x_1=x_2 = 1\), ali ga ne možemo riješiti Gaussovim eliminacijama kao u prošlom primjeru. Ono što možemo napraviti je mijenjati poredak jednadžbi kako bismo na element na dijagonali kojim poništavamo sve ostale elemente ispod njega doveli element koji nije jednak nuli.

Dok gornji primjer pokazuje slučaj kada u egzaktnoj aritmetici imamo probleme kada ne mijenjamo retke matrice sustava, primjer koji slijedi pokazuje da treba biti oprezan i s nekim primjerima koji se u egzaktnoj aritmetici mogu riješiti, jer mogu prouzročiti velike greške u računalu.

Primjer 2.2

Gaussovim eliminacijama riješimo sustav \(Ax = b\), gdje je

\[\begin{split}A = \begin{bmatrix} 10^{-6} & 1 \\ 1 & 1 \end{bmatrix}, \quad b = \begin{bmatrix} 1 \\ 2 \end{bmatrix}.\end{split}\]

Rješavamo u aritmetici računala s \(3\) dekadske znamenke mantise (iza decimalne točke). Pretpostavimo da jedina pogreška koju naše računalo radi prilikom svih operacija je samo u činjenici da koeficijente jednadžbi i rješenja pamti kao njemu prikazive brojeve.

U takvoj aritmetici prikazivi realni brojevi su oblika

\[\pm b_0.b_1b_2b_3 \cdot 10^e,\]

gdje su \(b_0,b_1,b_2,b_3 \in \{ 0,1,\dots, 9 \}\). Svi brojevi u početnom sustavu su prikazivi, pa zasad nema pogrešaka. Egzaktno rješenje sustava je

\[\begin{split}\begin{align*} x_1 &= 1+\dfrac{1}{10^6-1} = 1.\overline{000001}, \\ x_2 &= 1-\dfrac{1}{10^6-1} = 0.\overline{999998}. \end{align*}\end{split}\]

U „našem” računalu ti brojevi nisu prikazivi, pa se jedino možemo nadati da ćemo postupkom rješavanja sustava na računalu dobiti rješenje \(fl(x_1) = 1\) i \(fl(x_2)=1\).

Počnimo s rješavanjem. Drugoj jednadžbi oduzmimo prvu pomnoženu s \(10^6\). Druga jednadžba tada postaje

\[\begin{split}\begin{align*} (1-10^6\cdot 10^{-6})x_1 + (1-10^6\cdot 1)x_2 &= 2-10^6\cdot 1, \quad \text{odnosno} \\ (1-10^6)x_2 &= 2-10^6. \end{align*}\end{split}\]

Koeficijenti zadnje jednadžbe nalaze se između dva uzastopna prikaziva broja:

\[\begin{split}\begin{alignat*}{2} -1.000\cdot 10^6 &= \, \fl(-10^6) = -10^6 && \quad \text{prikaziv}\\ -9.99999 \cdot 10^5 &= 1-10^6 && \quad \text{nije prikaziv} \\ -9.99998 \cdot 10^5 &= 2-10^6 && \quad \text{nije prikaziv} \\ -9.999 \cdot 10^5 &= \, \fl(100-10^6) = 100-10^6 && \quad \text{prikaziv}. \end{alignat*}\end{split}\]

Računalo koeficijente te jednadžbe pohranjuje zaokruživanjem na najbliži prikazivi broj, a to je \(\fl(1-10^6)=\fl(2-10^6)=-10^6\). Zato je dobivena jednadžba u našem računalu zapravo \(-10^6 x_2 = -10^6\), odakle dobivamo rješenje \(x_2=1\).

S tim rješenjem ulazimo u prvu jednadžbu:

\[\begin{split}\begin{align*} 10^{-6} \cdot x_1 + x_2 &= 1, \quad \text{odnosno,} \\ 10^{-6} \cdot x_1 &= 0, \end{align*}\end{split}\]

odakle je \(x_1=0\). Budući da je pravo rješenje \(x_1\) približno jednako \(1\), dobili smo ogromnu grešku (relativna greška je \(1\)).

Ponukani prošlim primjerom, pokušajmo sada u istom računalu riješiti sustav u kojem su jednadžbe u drugačijem poretku:

\[\begin{split}A = \begin{bmatrix} 1 & 1 \\10^{-6} & 1 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}= \begin{bmatrix} 2 \\ 1 \end{bmatrix}.\end{split}\]

Drugoj jednadžbi oduzmimo prvu pomnoženu s \(10^{-6}\). Druga jednadžba tada postaje

\[\begin{split}\begin{align*} (10^{-6}-10^{-6}\cdot 1)x_1 + (1-10^{-6}\cdot 1)x_2 &= 1-10^{-6}\cdot 2, \quad \text{odnosno} \\ (1-10^{-6})x_2 &= 1-2\cdot 10^{-6}. \end{align*}\end{split}\]

Koeficijenti zadnje jednadžbe nalaze se između dva uzastopna prikaziva broja:

\[\begin{split}\begin{alignat*}{3} 9.999 &\cdot 10^{-1} &&= \, \fl(1-100 \cdot 10^{-6}) = 1-100\cdot 10^{-6} && \quad \text{prikaziv} \\ 9.99998 &\cdot 10^{-1} &&= 1-2\cdot 10^{-6} && \quad \text{nije prikaziv} \\ 9.99999 &\cdot 10^{-1} &&= 1-10^{-6} && \quad \text{nije prikaziv} \\ 1.000 &\cdot 10^0 &&= \, \fl(1) = 1 && \quad \text{prikaziv} \end{alignat*}\end{split}\]

Računalo koeficijente te jednadžbe pohranjuje zaokruživanjem na najbliži prikazivi broj, a to je \(\fl(1-10^{-6})=\fl(1-2\cdot10^{-6})=1\). Ponovno dobivamo rješenje \(x_2 = 1\), i ponovno to nije pogrešno (s obzirom na aritmetiku računala).

S tim rješenjem ulazimo u prvu jednadžbu:

\[\begin{split}\begin{align*} x_1 + x_2 &= 2, \quad \text{odnosno,} \\ x_1 &= 2-1 = 1. \end{align*}\end{split}\]

Sada nema pogreške: uzimajući u obzir prikazivost realnih brojeva u računalu, dobili smo najtočnije moguće rješenje sustava.

Ovaj primjer uči nas nekoliko stvari:

• Algoritam proveden u egzaktnoj aritmetici zakazuje kada je prvi dijagonalni element u matrici \(A\) jednak nula, a isti taj algoritam proveden u konačnoj aritmetici računala zakazuje kada je taj element blizu nule. Uzimajući neku intuiciju da se ideje u matematici ponašaju neprekidno, ovo nas ne bi smjelo iznenaditi.

• Ovaj primjer sam po sebi pokazuje da postoji primjer sustava u kojem Gaussove eliminacije loše funkcioniraju kada ne dopuštamo promjenu poretka jednadžbi u aritmetici računala. Možemo se pitati je li ovakav slučaj rijedak, odnosno je li on razlog da ozbiljno razmislimo o algoritmu rješavanja linearnih sustava. Naime, jedan od koeficijenata u sustavu je manji od jedinične greške zaokruživanja za to računalo – koliko je često da ćemo se susretati s takvim sustavima? Odgovor leži u tome što se u praksi ne pojavljuju nužno tako veliki omjeri među koeficijentima u jednadžbi, ali sustavi su mnogo veće dimenzije od \(n=2\). Zato se greška koju radimo u eliminirajući svaki redak jednadžbe može gomilati eksponencijalno.

• Do greške je došlo zato što smo drugoj jednadžbi oduzeli prvu pommnoženu velikim faktorom (\(10^6\)). Zbog toga je ta druga jednadžba postala „nebitna”, njeni koeficijenti izgubili su značaj. To pokazuje da, ako možemo birati, bolje je da su multiplikatori \(m_{i}^{(j)}\) što manji, a ne što veći. To ćemo postići tako da element na dijagonali kojim eliminiramo sve koeficijente ispod njega bude što veći (budući da se on nalazi u nazivniku izraza za \(m_{i}^{(j)}\)).

Ti zaključci dovode nas do sljedećeg algoritma.

# Rješavanje sustava Ax=b Gaussovim eliminacijama s parcijalnim pivotiranjem.
def GaussoveEliminacijeSPivotiranjem(A, b):
    n = A.shape[0]

    # Gaussove eliminacije.
    for i in range(n):
        # Promatramo i-ti stupac matrice A.
        # Pronalaženje maksimalnog elementa u i-tom stupcu ispod dijagonale.
        max_abs_val = abs(A[i, i])
        max_row = i
        for j in range(i + 1, n):
            if abs(A[j, i]) > max_abs_val:
                max_abs_val = abs(A[j, i])
                max_row = j

        # Zamjena i-tog redaka s retkom max_row.
        if max_row != i:
            A[[i, max_row], :] = A[[max_row, i], :]
            b[[i, max_row]] = b[[max_row, i]]

        # Eliminacija na isti način kao ranije.
        for j in range(i + 1, n):
            m = A[j, i] / A[i, i]
            A[j, :] -= m * A[i, :]
            b[j] -= m * b[i]

    # Rješavanje trokutastog sustava: supstitucija unatrag kao ranije.
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = b[i]
        for j in range (i+1, n):
            x[i] -= A[i, j]*x[j]
        x[i] /= A[i, i]

    return x

Razlika u odnosu na prošli algoritam: prije eliminacije pronalazimo najveći element u ostatku stupca po apsolutnoj vrijednosti i vršimo zamjenu tog retka s trenutnim.

Općenito, zamjenu redaka (ili stupaca) u matrici sustava nazivamo pivotiranje. Element koji dovodimo na dijagonalu i kojim poništavamo sve elemente ispod njega nazivamo pivotnim elementom. Iako u egzaktnoj aritmetici možemo odabrati bilo koji pivotni element različit od nule, u aritmetici računala biramo najveći po apsolutnoj vrijednosti kako bismo postigli da su multiplikatori kojim poništavamo preostale retke što manji, a time i račun što točniji.

Primjer 2.3

Gaussovim eliminacijama s parcijalnim pivotiranjem riješimo sustav \(Ax = b\), gdje je

\[\begin{split}A = \begin{bmatrix} 1 & 1 & 1 \\ 2 & 4 & 2 \\ -1 & 5 & -4 \end{bmatrix}, \quad b = \begin{bmatrix} 6 \\ 16 \\ -3 \end{bmatrix}.\end{split}\]

Tom sustavu odgovara proširena matrica \(\hat{A}^{(1)} = \left[ A \mid b \right] \). sustava glasi

\[\begin{split}\hat{A}^{(1)} = \begin{bmatrix} \begin{array}{ccc:c} 1 & 1 & 1 & 6 \\ 2 & 4 & 2 & 16 \\ -1 & 5 & -4 & -3 \end{array} \end{bmatrix}.\end{split}\]

Prvi korak zamjena je redaka tako da na mjesto \((1,1)\) dovedemo broj iz prvog stupca koji je najveći po apsolutnoj vrijednosti. To je broj \(2\), pa vršimo zamjenu prvog i drugog retka:

\[\begin{split}\tilde{A}^{(1)} = \begin{bmatrix} \begin{array}{ccc:c} 2 & 4 & 2 & 16 \\ 1 & 1 & 1 & 6 \\ -1 & 5 & -4 & -3 \end{array} \end{bmatrix}.\end{split}\]

Pivotnim elementom \(\tilde{a}_{1,1}^{(1)} = 2\) poništimo preostale elemente u prvom stupcu koji se nalaze ispod tog elementa. Prvi redak ne mijenjamo, drugi redak ćemo umanjiti za prvi redak pomnožen s \(m_{2}^{(1)}:= \dfrac{1}{2}\), a treći redak ćemo umanjiti za prvi redak pomnožen s \(m_{3}^{(1)}:= \dfrac{-1}{2}\). Novodobivena proširena matrica glasi

\[\begin{split}\hat{A}^{(2)} = \begin{bmatrix} \begin{array}{ccc:c} 2 & 4 & 2 & 6 \\ 0 & -1 & 0 & -2 \\ 0 & 7 & -3 & 5 \end{array} \end{bmatrix}.\end{split}\]

Ako je potrebno (a jest), ponovno na pivotno mjesto dovodimo najveći element po aspolutnoj vrijednosti u drugom stupcu (prvi redak ne promatramo). To je broj \(7\), pa treba zamijeniti drugi i treći redak:

\[\begin{split}\tilde{A}^{(2)} = \begin{bmatrix} \begin{array}{ccc:c} 2 & 4 & 2 & 6 \\ 0 & 7 & -3 & 5\\ 0 & -1 & 0 & -2 \end{array} \end{bmatrix}.\end{split}\]

Elementom \(\tilde{a}^{(2)}_{2,2} = 7\) poništavamo zadnji redak tako da ćemo ga umanjiti za drugi redak pomnožen s \(m_{3}^{(2)}:=\dfrac{1}{7}\). Nova proširena matrica sustava glasi

\[\begin{split}\hat{A}^{(3)} = \begin{bmatrix} \begin{array}{ccc:c} 2 & 4 & 2 & 16 \\ 0 & 7 & -3 & 5 \\ 0 & 0 & -\dfrac{3}{7} & -\dfrac{9}{7} \end{array} \end{bmatrix}.\end{split}\]

Ponovno je dobivena matrica sustava gornjetrokutasta:

\[\begin{split}U := \begin{bmatrix} 2 & 4 & 2 \\ 0 & 7 & -3 \\ 0 & 0 & -\dfrac{3}{7} \end{bmatrix},\end{split}\]

te sustav rješavamo povratnom supstitucijom. U zadnjem retku odgovarajuća je jednadžba s jednom nepoznanicom: \(-\dfrac 37 x_3 = -\dfrac 97\), odakle je \(x_3=3\). S tim rješenjem ulazimo u predzadnji redak. Odgovarajuća je jednadžba

\[\begin{split}\begin{align*} 7x_2 -3 x_3 &= 5, \quad \text{odnosno}\\ 7x_2 &= 5 + 3 \cdot 3 = 14, \end{align*}\end{split}\]

odakle je \(x_2=2\). Konačno, treba odrediti preostalu nepoznanicu iz prve jednadžbe \(2x_1+4x_2+2x_3 = 16\), tj.

\[x_1 = 16 - 4x_2 - 2x_3 = 1.\]

Konačno rješenje sustava je \(x = (1,2,3)^{T}\).

Rezultat možemo provjeriti u Pythonu našom funkcijom.

A = np.array([[1., 1., 1.], [2., 4., 2.], [-1., 5., -4.]])
b = np.array([6., 16., -3.])

x = GaussoveEliminacijeSPivotiranjem(A, b)
print( f'Rješenje danog sustava je {x}.' );

Rješenje danog sustava je [1. 2. 3.].

Primijetimo da su svi multiplikatori \(m_{i}^{(j)}\) po apsolutnoj vrijednosti manji ili jednaki \(1\). To je svojstvo parcijalnog pivotiranja, budući da se oni računaju kao omjer u kojem je u nazivniku najveći broj tog stupca po apsolutnoj vrijednosti.

Prokomentirajmo složenost algoritama:

• Za svaki \(i\) pronalazak najvećeg elementa u stupcu je složenosti \(\mathcal O (n)\) - to nam sugerira petlja po \(j\).

• Zamjena redaka je također složenosti \(\mathcal O(n)\) za svaki \(i\), jer do promjene dolazi u svakom stupcu tih dvaju redaka.

• Provođenje eliminacije je složenosti \(\mathcal O (n^2)\) za svaki \(i\): imamo petlju po \(j\), no unutar nje naredba \(\texttt{A[j, :] -= m * A[i, :]}\) označava ažuriranje cijelog retka - ta naredba sama je složenosti \(\mathcal O (n)\).

• Uzimajući u obzir da sve opisane dijelove algoritma treba promatrati u kontekstu da su unutar petlje po \(i\), ukupna složenost za provođenje samih Gaussovih eliminacija je \(\mathcal O (n^3)\).

• Rješavanje samog trokutastog sustava je \(\mathcal O (n^2)\), jer imamo dvije petlje po \(i\) i \(j\).

Dakle, ukupna složenost rješavanja sustava Gaussovim eliminacijama s parcijalnim pivotiranjem iznosi \(\mathcal O (n^3)\).

LU faktorizacija

Vratimo se ponovno na Primjer 2.1. Svaku od transformacija u kojoj jednom retku oduzimamo drugi pomnožen skalarom moguće je opisati elementarnom matricom oblika

\[\begin{split}L_{(i,j,m)} = \begin{bmatrix} 1 & & & & & & \\ & \ddots & & & & & \\ & & 1 & & & & \\ & & & \ddots & & & \\ & & m & & 1 & & \\ & & & & & \ddots & \\ & & & & & & 1 \end{bmatrix}\end{split}\]

(matrica je jednaka identiteti, osim na poziciji \((i,j)\) na kojoj se nalazi vrijednost \(m\)) gdje su \(i\) i \(j\) indeksi takvi da je \(1 \leq j<i\leq n\), a \(m\) proizvoljan realni broj. Zaista, redom smo imali

\[\begin{split}\begin{align*} \begin{bmatrix} 1 & 1 & 1 \\ 0 & 2 & 0 \\ 0 & 6 & -3 \end{bmatrix} &= {A}^{(2)} = L_{(3,1,1)} L_{(2,1,-2)} A^{(1)} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 1 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ -2 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 1 & 1 \\ 2 & 4 & 2 \\ -1 & 5 & -4 \end{bmatrix}, \\ \begin{bmatrix} 1 & 1 & 1 \\ 0 & 2 & 0 \\ 0 & 0 & -3 \end{bmatrix} &= U = {A}^{(3)} = L_{(3,2,-3)} A^{(2)} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -3 & 1 \end{bmatrix} \begin{bmatrix} 1 & 1 & 1 \\ 0 & 2 & 0 \\ 0 & 6 & -3 \end{bmatrix} \end{align*}\end{split}\]

Dakle, matrica \(U\) dobivena je kao umnožak nekoliko matrica \(L_{(i,j,m)}\) i matrice \(A=A^{(1)}\):

\[U = L_{(3,2,-3)} L_{(3,1,1)} L_{(2,1,-2)} A.\]

Matrice \(L_{(i,j,m)}\) su regularne i njihov inverz je \(L_{(i,j,-m)}\) (provjerite sami za vježbu). Zato možemo pisati

(2.2)\[A = L_{(2,1,2)} L_{(3,1,-1)} L_{(3,2,3)} U.\]

Sve matrice oblika \(L_{(i,j,m)}\) su donjetrokutaste s jedinicama na dijagonali. Njihov umnožak ponovno je donjetrokutasta matrica s jedinicama na dijagonali (provjerite i to sami). Označimo li taj umnožak s \(L\), matricu \(A\) zapisali smo u obliku \(A = LU\).

Definicija 2.4

Neka je \(A\) kvadratna matrica reda \(n\). Kažemo da je prikaz \(A=LU\) LU faktorizacija matrice \(A\) ako je \(L\) donjetrokutasta matrica s jedinicama na dijagonali, a \(U\) gornjetrokutasta matrica, pri čemu su obje matrice \(L\) i \(U\) kvadratne reda \(n\).

Vratimo se na jednakost (2.2). U tom primjeru je

\[\begin{split}L = L_{(2,1,2)} L_{(3,1,-1)} L_{(3,2,3)} = \begin{bmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & 3 & 1 \end{bmatrix} .\end{split}\]

Dakle, na svakom mjestu \((i,j)\), \(i>j\), u matrici \(L\) je zapisan upravo multiplikator \(m_{i}^{(j)}\) iz elementarne matrice \(L_{(i,j,m_{i}^{(j)})}\). To nije slučajno, to je još jedno svojstvo elementarnih matrica koje možete pokazati sami i za proizvoljan \(n\) (ukoliko se množenje vrši u pravom redoslijedu).

Primjer (2.1) pokazuje da na nekim matricama nije moguće izvršiti Gaussove eliminacije. Za istu matricu može se pokazati da postoje matrice za koje ne postoji LU faktorizacija. Naime, kada bi postojala, vrijedilo bi

\[\begin{split}\begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ l_{2,1} & 1 \end{bmatrix} \begin{bmatrix} u_{1,1} & u_{2,1} \\ 0 & u_{2,2} \end{bmatrix}.\end{split}\]

Jednakost za element \((1,1)\) nalaže da je \(0=u_{1,1}\), no tada bi zadnja matrica bila singularna (prvi stupac joj je jednak nula), čime dolazimo do kontradikcije jer matrica na lijevoj strani nije singularna.

Sljedeći teorem daje nužan i dovoljan uvjet da bi regularna matrica imala LU faktorizaciju.

Teorem 2.5

Neka je \(A\) kvadratna regularna matrica reda \(n\). Tada postoji i jedinstvena je LU faktorizacija matrice \(A\) ako i samo ako su podmatrice \(A(1:k,1:k)\) regularne za sve \(k = 1,\dots,n-1\).

U gornjem teoremu, \(A(1:k,1:k)\) označava matricu reda \(k\) dobivenu iz \(A\) brisanjem svih redaka osim prvih \(k\) te brisanjem svih stupaca osim prvih \(k\).

Oznaku \(A(i_1:i_2, j_1:j_2)\) za podmatricu od \(A\) koja se sastoji od svih njezinih elemenata iz redaka \(i_1, i_1+1, \ldots, i_2\) i stupaca \(j_1, j_1+1, \ldots, j_2\) nazivamo Matlab notacija.

Dokaz.

Dokaz provodimo matematičkom indukcijom. Za \(n=1\) svaka matrica \(A=[a]\) ima jedinstvenu \(LU\) faktorizaciju \(LU = [1] \cdot [a]\).

Pretpostavimo da tvrdnja vrijedi za svaku regularnu matricu reda \(n-1\) i neka je \(A\) neka regularna matrica reda \(n\). Za nju trebamo dokazati tvrdnju iz teorema u oba smjera. Umjesto klasičnog označavanja matrice po elementima

\[\begin{split}A = \begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\ a_{2,1} & a_{2,2} & \dots & a_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n,1} & a_{n,2} & \dots & a_{n,n} \end{bmatrix},\end{split}\]

uvedimo drugačije oznake. Neka su \(A_1 \in \R^{(n-1) \times (n-1)}\), \(a,b \in \R^{n-1}\) te \(\alpha \in \R\) takvi da je

\[\begin{split}A_1 = \begin{bmatrix} a_{1,1} & a_{1,2} & \cdots & a_{1,n-1} \\ a_{2,1} & a_{2,2} & \dots & a_{2,n-1} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n-1,1} & a_{n-1,2} & \dots & a_{n-1,n-1} \end{bmatrix}, \quad a = \begin{bmatrix} a_{1,n} \\ a_{2,n} \\ \vdots \\ a_{n-1,n} \end{bmatrix}, \quad b = \begin{bmatrix} a_{n,1} \\ a_{n,2} \\ \vdots \\ a_{n,n-1} \end{bmatrix}, \quad \alpha = a_{n,n},\end{split}\]

odnosno da vrijedi

(2.3)\[\begin{split}A = \begin{bmatrix} A_1 & a \\ b^T & \alpha \end{bmatrix}.\end{split}\]

Proizvoljnu donjetrokutastu matricu \(L \in \R^{n \times n}\) s jedinicama na dijagonali i gornjetrokutastu matricu \(U \in \R^{n \times n}\) možemo rastaviti na analogan način: za svaku takvu \(L\) i \(U\) postoje donjetrokutasta matrica s jedinicama na dijagonali \(L_1 \in \R^{(n-1) \times (n-1)}\), gornjetrokutasta matrica \(U_1 \in \R^{(n-1) \times (n-1)}\), vektori \(c,d \in \mathbb R^{n-1}\) i skalar \(\beta \in \mathbb R\) takvi da je

\[\begin{split}L = \begin{bmatrix} L_1 & 0 \\ c^T & 1 \end{bmatrix}, \quad U = \begin{bmatrix} U_1 & d \\ 0 & \beta \end{bmatrix}.\end{split}\]

Matrica \(A\) ima LU faktorizaciju ako i samo ako postoje matrice \(L\) i \(U\) iz definicije LU faktorizacije takve da je \(A=LU\). Uzimajući gore prikazani rastav, možemo reći da matrica \(A\) prikazana kao u (2.3) ima LU faktorizaciju ako i samo ako je

(2.4)\[\begin{split}A = \begin{bmatrix} A_1 & a \\ b^T & \alpha \end{bmatrix} = \begin{bmatrix} L_1 & 0 \\ c^T & 1 \end{bmatrix} \begin{bmatrix} U_1 & d \\ 0 & \beta \end{bmatrix} = LU\end{split}\]

Pravila množenja matrica po komponentama ovakvih blok-matrica su analogna množenju \(2\times 2\) matrica. Konačno zaključujemo: matrica \(A\) iz (2.3) ima LU faktorizaciju ako i samo ako postoje donjetrokutasta matrica s jedinicama na dijagonali \(L_1 \in \R^{(n-1) \times (n-1)}\), gornjetrokutasta matrica \(U_1 \in \R^{(n-1) \times (n-1)}\), vektori \(c,d \in \mathbb R^{n-1}\) i skalar \(\beta \in \mathbb R\) takvi da je

(2.5)\[\begin{split}\begin{matrix} A_1 = L_1 U_1, & a = L_1 d, \\ b^T = c^T U_1, & \alpha = c^Td + \beta. \end{matrix}\end{split}\]

Prijeđimo napokon na dokaz svake od implikacija u iskazu teorema. Pretpostavimo prvo da regularna matrica \(A\) ima \(LU\) faktorizaciju, i da je ona oblika (2.4). Dokažimo da podmatrice \(A(1:k, 1:k)\) moraju biti regularne za sve \(k=1, 2, \ldots, n-1\). Provodeći Laplaceov razvoj po zadnjem retku matrice \(U\) slijedi da je \(\det U = \beta \det U_1\). Slično, \(\det L = \det L_1\). Primjenom Binet-Cauchyjevog teorema na \(A=LU\) slijedi

\[\det A = \det L \cdot \det U = \beta \det L_1 \cdot \det U_1.\]

Kako je matrica \(A\) regularna, \(\det A \not=0\), a to znači i da nijedan od faktora na desnoj strani nije jednak nula, dakle \(L_1\) i \(U_1\) su regularne matrice. Iz \(A_1 = L_1 U_1\) slijedi da je i matrica \(A_1\) regularna, dakle \(A(1:n-1,1:n-1)\) je regularna matrica. Konačno, jednakost \(A_1 = L_1 U_1\) dokazuje da matrica \(A_1\) ima LU faktorizaciju, što korištenjem pretpostavke indukcije znači da su sve matrice \(A_1(1:k,1:k) = A(1:k,1:k)\), za \(k=1,\dots,n-2\), regularne. Time je dokazana jedna implikacija iz teorema u koraku indukcije.

Pretpostavimo sada da su sve matrice \(A(1:k,1:k)\) regularne za \(k = 1,\dots,n-1\) i dokažimo da \(A\) dopušta LU fakorizaciju, tj. da postoje objekti \(L_1,U_1,b,c,\beta\) takvi da vrijede (2.5). Iz pretpostavke posebno slijedi da je \(A_1 = A(1:n-1,1:n-1)\) je regularna matrica i sve podmatrice \(A_1(1:k,1:k)\) su regularne za \(k = 1,\dots,n-2\). Po pretpostavci indukcije, \(A_1\) ima jedinstvenu LU faktorizaciju, pa postoje jedinstvena donjetrokutasta matrica s jedinicama na dijagonali \(L_1 \in \R^{(n-1) \times (n-1)}\) i gornjetrokutasta matrica \(U_1 \in \R^{(n-1) \times (n-1)}\) takve da je \(A_1 = L_1 U_1\). Po Binet-Cauchyjevom teoremu, \(\det A_1 = \det L_1 \cdot \det U_1\). Kako je \(\det A_1 \not=0\), i \(L_1\) i \(U_1\) su regularne matrice.

Jednadžba \(a = L_1 d\) sustav je linearnih jednadžbi s poznatim vektorom \(a\) i nepoznatim vektorom \(d\). Kako je \(L_1\) regularna matrica, postoji jedinstveno rješenje \(d\). Jednadžba \(b^T = c^T U_1\) ekvivalentna je s \(b = U_1^T c\), pa iz analognog razloga definira jedinstveni vektor \(c\). Konačno, iz jednadžbe \(\alpha = c^Td + \beta\) dobivamo jedinstvenu vrijednost za \(\beta\). Svi elementi na desnoj strani jednadžbe (2.4) su definirani, pa to dokazuje da \(A\) ima jedinstvenu LU faktorizaciju, što je i trebalo dokazati.

Za neku matricu, njezina \(LU\) faktorizacija, osim što nudi lakši zapis i analizu svojstava procesa rješavanja sustava Gaussovim eliminacijama, nudi i alternativan način rješavanja sustava. Naime, za dani sustav \(Ax = b\), ukoliko znamo LU faktorizaciju \(A = LU \), možemo pisati

\[LUx = Ax = b,\]

te uz supstituciju \(y=Ux\) dobiti sustav \(Ly = b\). U stvari, imamo dva trokutasta sustava. U tom slučaju, algoritam za rješavanje sustava LU faktorizacijom glasi

1. Provesti LU faktorizaciju matrice \(A\), odakle dobivamo matrice \(L\) i \(U\) \(\; \ldots \; \bigO(n^3)\);

2. Odrediti vektor \(y\) koji zadovoljava sustav \(Ly=b\) \(\; \ldots \; \bigO(n^2)\);

3. Odrediti vektor \(x\) koji zadovoljava sustav \(Ux=y\) \(\; \ldots \; \bigO(n^2)\).

Matrica sustava \(Ly=b\) je donjetrokutasta. Kao što sustave s gornjetrokutastom matricom rješavamo supstitucijom unatrag, tako sustave s donjetrokutastom matricom rješavamo supstitucijom unaprijed: prvo određujemo prvu komponentu vektora \(y\), pa drugu itd.

Ovaj način rješavanja nudi još jednu prednost u odnosu na Gaussove eliminacije. Zamislimo da trebamo riješiti nekoliko sustava \(Ax=b^{(k)}\) za razne desne strane \(b^{(k)}\). Budući da proces rješavanja Gaussovim eliminacijama počiva na transformacijama koje provodimo na proširenoj matrici sustava, sve korake moramo ponovno provoditi za svaki \(k\). S druge strane, kada jednom napravimo \(LU\) faktorizaciju (koja ovisi samo o matrici \(A\), ne i o vektoru desne strane), za svaku desnu stranu \(b^{(k)}\) dovoljno je riješiti dva trokutasta sustava. To je dosta brže, jer složenost jedne \(LU\) faktorizacije je \(\bigO(n^3)\), dok je složenost rješavanja svakog od trokutastog sustava \(\bigO(n^2)\). Dakle, nakon što jednom napravimo (skupu) LU-faktorizaciju matrice \(A\) složenosti \(\bigO(n^3)\), rješavanje svakog sustava s matricom \(A\) i bilo kojom desnom stranom ima vremensku složenost od samo \(\bigO(n^2)\).

Rezultate o stabilnosti rješavanja sustava linearnih jednadžbi ćemo također moći iskazati koristeći faktore \(L\) i \(U\).

LU faktorizacija s parcijalnim pivotiranjem

Kao što Gaussove eliminacije bez pivotiranja vode na LU faktorizaciju, slično vrijedi i za Gaussove eliminacije s parcijalnim pivotiranjem. Promjenu poretka općenito možemo opisati matricom \(P_{i,j}\) oblika

\[\begin{split}\begin{bmatrix} 1 & & & & & & \\ & \ddots & & & & & \\ & & 0 & & 1 & & \\ & & & \ddots & & & \\ & & 1 & & 0 & & \\ & & & & & \ddots & \\ & & & & & & 1 \end{bmatrix}\end{split}\]

(matrica je jednaka identiteti, osim na mjestima \((i,i), (i,j), (j,i), (j,j)\)), pri čemu su \(i\) i \(j\) različiti indeksi. Svaka matrica \(P_{i,j}\) poseban je slučaj matrice permutacije.

Definicija 2.6

Kažemo da je kvadratna matrica \(P\) matrica permutacije ako u svakom retku ima točno jednu jedinicu, a sve su ostalo nule, te ako u svakom stupcu ima točno jednu jedinicu, a sve su ostalo nule.

Primijetite da je svaka matrica permutacija \(P\) ortogonalna matrica, pa posebno vrijedi

• regularna je,

• \(\det P \in \{ \pm 1 \}\),

• \(P^{-1} = P^T\).

Promotrimo sada Primjer 2.3. Matricu \(\tilde{A}^{(1)}\) dobili smo permutiranjem redaka matrice \(A\), a zatim smo na nju djelovali matricama \(L_{(i,j,m)}\) da dobijemo \(\hat{A}^{(2)}\). Dakle, imamo

\[\begin{split}\begin{align*} \begin{bmatrix} 2 & 4 & 2 \\ 0 & -1 & 0 \\ 0 & 7 & -3 \end{bmatrix} &= \hat{A}^{(2)} = L_{\left(3,1,\frac{1}{2}\right)}L_{\left(2,1,-\frac{1}{2}\right)}P_{1,2} A \\ &= \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ \dfrac{1}{2} & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ -\dfrac{1}{2} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 1 & 1 \\ 2 & 4 & 2 \\ -1 & 5 & -4 \end{bmatrix}. \end{align*}\end{split}\]

Nakon toga smo zamijenili drugi i treći redak, pa izvršili transformaciju matricom \(L_{(3,2,\frac{1}{7})}\):

\[\begin{split}\begin{align*} \begin{bmatrix} 2 & 4 & 2 \\ 0 & 7 & -3 \\ 0 & 0 & -\dfrac{3}{7} \end{bmatrix} &= U = \hat{A}^{(3)} = L_{\left(3,2,\frac{1}{7}\right)}P_{2,3} \hat{A}^{(2)}\\ &= \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & \dfrac{1}{7} & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} 2 & 4 & 2 \\ 0 & -1 & 0 \\ 0 & 7 & -3 \end{bmatrix}. \end{align*}\end{split}\]

Dakle, vrijedi

\[U = L_{\left(3,2,\frac{1}{7}\right)}P_{2,3} L_{\left(3,1,\frac{1}{2}\right)}L_{\left(2,1,-\frac{1}{2}\right)}P_{1,2} A.\]

Iz ove jednakosti htjeli bismo dobiti jednakost što sličniju \(A=LU\). Jasno je da matrice permutacije tu „smetaju”, ali i da trebaju ostati na kraju u jednakosti. To ćemo postići tako da matricu \(P_{2,3}\) dovedemo desno uz matricu \(P_{2,1}\). Naime, vrijedi

\[\begin{split}\begin{align*} &P_{2,3} L_{\left(3,1,\frac{1}{2}\right)}L_{\left(2,1,-\frac{1}{2}\right)} \\&= \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ \dfrac{1}{2} & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ -\dfrac{1}{2} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \\ &= \begin{bmatrix} 1 & 0 & 0 \\ \dfrac{1}{2} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ -\dfrac{1}{2} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \\ &= \begin{bmatrix} 1 & 0 & 0 \\ \dfrac{1}{2} & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -\dfrac{1}{2} & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \end{bmatrix} \\ &= L_{\left(2,1,\frac{1}{2}\right)}L_{\left(3,1,-\frac{1}{2}\right)}P_{2,3}. \end{align*}\end{split}\]

Dakle, vrijedi

\[U = L_{\left(3,2,\frac{1}{7}\right)}L_{\left(2,1,\frac{1}{2}\right)}L_{\left(3,1,-\frac{1}{2}\right)}P_{2,3} P_{1,2} A.\]

Prebacimo li matrice \(L_{(i,j,m)}\) na lijevu stranu (uzimajući im inverz) i taj umnožak označimo s \(L\), a umnožak permutacijskih matrica s \(P\), dobili smo \(LU=PA\).

Definicija 2.7

Neka je \(A\) kvadratna matrica reda \(n\). Kažemo da je prikaz \(PA=LU\) LU faktorizacija s parcijalnim pivotiranjem matrice \(A\) ako je \(P\) matrica permutacije, \(L\) donjetrokutasta matrica s jedinicama na dijagonali, a \(U\) gornjetrokutasta matrica, pri čemu su sve matrice \(P\), \(L\) i \(U\) kvadratne reda \(n\).

U gornjem primjeru pokazali smo kako konkretno matricu \(P_{2,3}\) možemo dovesti na desnu stranu umnoška elementarnih matrica. To se može ponoviti i u općenitom slučaju za proizvoljnu dimenziju matrice \(n\). Dakle, čini se da ovu faktorizaciju možemo uvijek provesti. Zaista, imamo sljedeći teorem.

Teorem 2.8

Neka je \(A\) kvadratna matrica reda \(n\). Tada postoje gornjetrokutasta matrica \(U\), donjetrokutasta matrica s jedinicama na dijagonali \(L\) i permutacijska matrica \(P\), sve reda \(n\), takve da je \(PA=LU\).

Dokaz.

Tvrdnju ćemo dokazati matematičkom indukcijom. Za \(n=1\), tvrdnja je trivijalna: za \(A:=[a]\) imamo odabir \(P=[1]\), \(L=[1]\), \(U=[a]\).

Pretpostavimo da tvrdnja vrijedi za svaku matricu reda \(n-1\), i neka je \(A\) proizvoljna kvadratna matrica reda \(n\). Dokažimo da ta matrica dopušta LU fakorizaciju s parcijalnim pivotiranjem.

Za tu matricu provedimo jedan korak LU faktorizacije s parcijalnim pivotiranjem. Prvo na mjesto \((1,1)\) dovedimo najveći po apsolutnoj vrijednosti element u prvom stupcu. Ako se nalazi u retku \(r\), tada ćemo matricu \(A\) pomnožiti s matricom \(P_{1,r}\) slijeva i dobivamo \(\tilde A := P_{1,r} A\). Matricu \(\tilde A\) rastavimo slično kao u dokazu Teorema 2.5. Postoje matrica \(B \in \R^{(n-1) \times (n-1)}\), vektori \(a,b \in \mathbb R^{n-1}\) i skalar \(\alpha \in \R\) takvi da je

\[\begin{split}\tilde A = \begin{bmatrix} \alpha & b^T \\ a & B \end{bmatrix}.\end{split}\]

gdje su \(\alpha \in \R\), \(a,b \in \R^{n-1}\) i \(B \in \R^{(n-1) \times (n-1)}\). Broj \(\alpha\) upravo je najveći po apsolutnoj vrijednosti u prvom stupcu.

Pretpostavimo prvo da je \(\alpha \not = 0\). Poništavajući sve elemente ispod pivotnog elementa u prvom stupcu, dobivamo multiplikatore koje smještamo u vektor \(c \in \mathbb R^{n-1}\) koji se nalaze u prvom stupcu donjetrokutaste matrice s jedinicama na dijagonali. U novoj matrici prvi redak matrice ostaje isti, dok su elementi prvog stupca ispod prvog elementa jednaki nula. Preciznije, dobili smo neki vektor \(c\in \R^{n-1}\) i matricu \(C \in \R^{(n-1) \times (n-1)}\) takve da je

(2.6)\[\begin{split}\begin{bmatrix} \alpha & b^T \\ a & B \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ c & I \end{bmatrix} \begin{bmatrix} \alpha & b^T \\ 0 & C \end{bmatrix}.\end{split}\]

Ovdje je \(I\) identiteta reda \(n-1\). Uspoređujući po komponentama, \(c\) i \(C\) zadovoljavaju jednakosti

(2.7)\[\begin{matrix} a = \alpha c, & B = ab^T + C. \end{matrix}\]

Ako \(\alpha\) jest jednak nuli, kako je to najveći element po apsolutnoj vrijednosti u prvom stupcu matrice \(\tilde A\), zaključujemo da je i \(a=0\). Ako definiramo \(c=0\) i \(C=B\) jednakost (2.6), odnosno obje jednakosti (2.7), i dalje su zadovoljene.

Za matricu \(C = B - ab^T\) primijenimo pretpostavku indukcije: postoje \(P_1\), \(L_1\) i \(U_1\) (permutacijska, donjetrokutasta s jedinicama na dijagonali i gornjetrokutasta redom, sve reda \(n-1\)) takve da je \(P_1 C = L_1 U_1\), odnosno \(C =P_1^{T} L_1 U_1\). Direktnom provjerom po komponentama kao gore, možemo se uvjeriti u niz jednakosti

\[\begin{split}\begin{align*} \tilde A &= \begin{bmatrix} 1 & 0 \\ c & I \end{bmatrix} \begin{bmatrix} \alpha & b^T \\ 0 & C \end{bmatrix} \\ &= \begin{bmatrix} 1 & 0 \\ c & I \end{bmatrix} \begin{bmatrix} \alpha & b^T \\ 0 & P_1^{T}L_1U_1 \end{bmatrix} \\ &= \begin{bmatrix} 1 & 0 \\ c & I \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & P_1^T \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & L_1 \end{bmatrix} \begin{bmatrix} \alpha & b^T \\ 0 & U_1 \end{bmatrix} \\ &= \begin{bmatrix} 1 & 0 \\ 0 & P_1^T \end{bmatrix} \begin{bmatrix} 1 & 0 \\ P_1c & I \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & L_1 \end{bmatrix} \begin{bmatrix} \alpha & b^T \\ 0 & U_1 \end{bmatrix} \end{align*}\end{split}\]

Kako je matrica

\[\begin{split}U:=\begin{bmatrix} \alpha & b^T \\ 0 & U_1 \end{bmatrix}\end{split}\]

gornjetrokutasta, umnožak matrica

\[\begin{split}L:=\begin{bmatrix} 1 & 0 \\ P_1c & I \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & L_1 \end{bmatrix}\end{split}\]

donjetrokutasta matrica s jedinicama na dijagonali, te umnožak

\[\begin{split}P:= \begin{bmatrix} 1 & 0 \\ 0 & P_1 \end{bmatrix} P_{1,r}\end{split}\]

permutacijska matrica (sve po pretpostavci indukcije), slijedi da je \(PA=LU\) LU faktorizacija matrice \(A\) s parcijalnim pivotiranjem, čime je dokaz indukcijom završen.

Ovakvu faktorizaciju također možemo direktno iskoristiti za rješavanje sustava \(Ax=b\). Iz \(Pb = PAx = LUx\) slijedi algoritam:

1. Odrediti LU faktorizaciju s parcijalnim pivotiranjem matrice \(A\), odakle dobivamo matrice \(L\), \(U\) i \(P\);

2. Odrediti vektor \(y\) koji zadovoljava sustav \(Ly=Pb\);

3. Odrediti vektor \(x\) koji zadovoljava sustav \(Ux=y\).

Pogledajmo primjer u Pythonu:

Show code cell source Hide code cell source

A = np.array([[1., 1., 1.], [2., 4., 2.], [-1., 5., -4.]])
b = np.array([6., 16., -3.])

# LU faktorizacija s parcijalnim pivotiranjem.
# Funkcija lu(A) iz biblioteke scipy.linalg vraća matrice takve da je A = Pt * L * U.
(Pt, L, U) = lu(A)
# Za potrebe naše notacije nužno je transponirati matricu P.
P = Pt.T
# Sada vrijedi P*A = L*U.

print( f'Matrica L:\n {L}\n' )
print( f'Matrica U:\n {U}\n' )
print( f'Matrica P:\n {P}\n' )

# Rješavanje donjetrokutastog sustava Ly = P * b.
# Koristimo gotovu funkciju solve_triangular iz scipy.linalg.
y = solve_triangular( L, P @ b, lower=True )
print( f'Međurješenje y:\n {y}\n' )

# Rješavanje gornjetrokutastog sustava Ux = y.
x = solve_triangular( U, y, lower=False )
print( f'Rješenje sustava x:\n {x}\n' )

Matrica L:
 [[ 1.          0.          0.        ]
 [-0.5         1.          0.        ]
 [ 0.5        -0.14285714  1.        ]]

Matrica U:
 [[ 2.          4.          2.        ]
 [ 0.          7.         -3.        ]
 [ 0.          0.         -0.42857143]]

Matrica P:
 [[0. 1. 0.]
 [0. 0. 1.]
 [1. 0. 0.]]

Međurješenje y:
 [16.          5.         -1.28571429]

Rješenje sustava x:
 [1. 2. 3.]

Kompaktan zapis LU faktorizacije

Zbog prirode algoritma, kod LU faktorizacije matrice \(A\) nije potrebno zauzeti novih \(2n^2\) mjesta u memoriji za čuvanje matrica \(L\) i \(U\). Kao prvo, matrica \(U\) je gornjetrokutasta, a matrica \(L\) je donjetrokutasta s jedinicama na dijagonali. Svi elementi tih dviju matrica mogu se čuvati u jednoj matrici (ne moramo pamtiti dijagonalu matrice \(L\)). Kao drugo, u \(k\)-tom koraku provođenja LU faktorizacije elementi iznad dijagonale u prvih \(k-1\) redaka jednaki su elementima matrice \(U\), elementi u presjeku zadnjih \(n-k\) stupaca i \(n-k\) redaka čekaju ostatak algoritma, a elementi ispod dijagonale u prvih \(k-1\) stupaca su jednaki nulama. Kako znamo da su to nule, ne moramo ih pamtiti eksplicitno, nego ta mjesta možemo iskoristiti da memoriramo multiplikatore koji će završiti u matrici \(L\). Ako radimo tako, na kraju algoritma u memoriji gdje je nekoć pisala matrica \(A\) pišu sve informacije o matricama \(L\) i \(U\), i pri tome nismo alocirali novu memoriju.

Za matricu permutacije \(P\) također ne moramo alocirati \(n^2\) mjesta u memoriji budući da su svi osim \(n\) elemenata te matrice jednaki nulama. Najčešće se za nju pamti vektor \(p\) takav da je \(p_i = j \iff P_{i,j}=1\).

Više primjera obradit će se na vježbama.

Potpuno pivotiranje

Za kraj ovog dijela, spomenimo da uz parcijalno pivotiranje postoji i potpuno pivotiranje: u svakom koraku na dijagonalno mjesto \((i,i)\) dovodimo po apsolutnoj vrijednosti najveći element u matrici \(\tilde{A}^{(i)}(i:n,i:n)\) (umjesto onog koji je najveći u preostalom \(i\)-tom stupcu). U praksi se pokazuje ovakvo pivotiranje usporava trošak LU faktorizacije (ukupni trošak pivotiranja također je \(\bigO(n^3)\)), a parcijalno pivotiranje je već samo po sebi daje dovoljno točne rezultate, pa se potpuno pivotiranje rijetko koristi.

Posebne familije matrica

Promotrimo neke matrice \(A\) posebnog oblika.

U praksi se pojavljuju rijetko popunjene matrice velikih dimenzija (\(n \times n\) matrice kod kojih je samo \(\bigO(n)\) elemenata različito od nule). Za takve matrice prilagođen je način memoriranja matrice (kako se ne bi zauzelo \(\bigO(n^2)\) mjesta u memoriji za \(\bigO(n)\) elemenata) pa su i algoritmi zbrajanja i množenja prilagođeni. Tako treba prilagoditi i algoritme za rješavanje sustava. Ukoliko ne postoji struktura kako je ta matrica popunjena, LU faktorizacija može biti loš izbor algoritma za rješavanje linearnog sustava budući da ne čuva tu rijetku popunjenost (neke nule mogu postati ne-nule i drastično povećati utrošak memorije i složenost).

S druge strane, može postojati neka posebna struktura rijetko popunjenih matrica. Najpoznatije su vrpčaste matrice. Za neki \(k \in \mathbb N\), to su matrice kod kojih su svi nenul elementi za najviše \(k\) mjesta udaljeni od dijagonale. Za \(k=1\), svi nenul elementi su na dijagonali, neposredno ispod ili neposredno ispod. Takve matrice zovu se tridijagonalnim matricama. LU faktorizacija bez pivotiranja čuva strukturu vrpčastih matrica, odnosno i matrice \(L\) i \(U\) su također vrpčaste s jednakim \(k\). To smanjuje složenost izvođenja \(LU\) faktorizacije.

Ako vrpčasta matrica ima LU faktorizaciju, onda su i faktori također vrpčaste matrice.

Na slici dolje, širina vrpce je \(k=2\) (tj. \(A\) je 5-dijagonalna matrica). Lako se vidi (pokažite!) da tada:

• \(L\) može imati samo 2 pod-dijagonale s ne-nul elementima;

• \(U\) može imati samo 2 nad-dijagonale s ne-nul elementima.

\[\begin{split} A = \smb \xx & \xx & \xx & 0 & 0 & 0 & 0 & 0 \\ \xx & \xx & \xx & \xx & 0 & 0 & 0 & 0 \\ \xx & \xx & \xx & \xx & \xx & 0 & 0 & 0 \\ 0 & \xx & \xx & \xx & \xx & \xx & 0 & 0 \\ 0 & 0 & \xx & \xx & \xx & \xx & \xx & 0 \\ 0 & 0 & 0 & \xx & \xx & \xx & \xx & \xx \\ 0 & 0 & 0 & 0 & \xx & \xx & \xx & \xx \\ 0 & 0 & 0 & 0 & 0 & \xx & \xx & \xx \sme = \underbrace{\smb 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \xx & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ \xx & \xx & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & \xx & \xx & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & \xx & \xx & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \xx & \xx & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & \xx & \xx & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & \xx & \xx & 1 \sme}_L \cdot \underbrace{\smb \xx & \xx & \xx & 0 & 0 & 0 & 0 & 0 \\ 0 & \xx & \xx & \xx & 0 & 0 & 0 & 0 \\ 0 & 0 & \xx & \xx & \xx & 0 & 0 & 0 \\ 0 & 0 & 0 & \xx & \xx & \xx & 0 & 0 \\ 0 & 0 & 0 & 0 & \xx & \xx & \xx & 0 \\ 0 & 0 & 0 & 0 & 0 & \xx & \xx & \xx \\ 0 & 0 & 0 & 0 & 0 & 0 & \xx & \xx \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \xx \sme}_U. \end{split}\]

Problem može jedino nastati samo ako vrpčaste matrice treba pivotirati kako bi algoritam zadržao točnost. Promotrimo dvije familije matrica kod kojih pivotiranje nije potrebno.

Definicija 2.9

Kažemo da je kvadratna matrica \(A\) reda \(n\) dijagonalno dominantna po stupcima ako za svaki \(j \in \{ 1, \dots, n \}\) vrijedi

\[\begin{split}|a_{j,j}| \geq \sum_{\substack{i=1 \\ i\not=j}}^{n} | a_{i,j}|.\end{split}\]

Kažemo da je kvadratna matrica \(A\) reda \(n\) dijagonalno dominantna po retcima ako za svaki \(i \in \{ 1, \dots, n \}\) vrijedi

\[\begin{split}|a_{i,i}| \geq \sum_{\substack{j=1 \\ i\not=j}}^{n} | a_{i,j}|.\end{split}\]

Reći ćemo da je matrica strogo dijagonalno dominantna po retcima/stupcima ako vrijedi stroga nejednakost u svih \(n\) nejednakosti.

Teorem 2.10

• Neka je \(A\) regularna i dijagonalno dominantna matrica po stupcima. Tada je prilikom LU faktorizacije s parcijalnim pivotiranjem za matricu \(A\) matrica permutacije jednaka identiteti, odnosno, \(A=LU\) ujedno je i LU faktorizacija s parcijalnim pivotiranjem. Drugim riječima, u tom slučaju matricu \(A\) ne treba pivotirati.

• Neka je \(A\) regularna i dijagonalno dominantna matrica po retcima. Tada matrica \(A\) dopušta \(LU\) faktorizaciju bez pivotiranja.

Dokaz.

Dokažimo tvrdnju samo za matricu koja je dijagonalno dominantno po stupcima, na kraju ćemo komentirati kako se dokazuje tvrdnja za matrice koja je dijagonalna po retcima.

Primijetimo da je element na mjestu \((1,1)\) sigurno različit od nule, u suprotnom ne bi mogao biti po apsolutnoj vrijednosti veći od sume apsolutnih vrijednosti elemenata u prvom stupcu osim ako je cijeli stupac nulstupac, što se u kontradikciji s regularnošću matrice \(A\). Zato je prvi korak LU faktorizacije sigurno moguće provesti. Dapače, kako je taj element po apsolutnoj vrijednosti veći od zbroja svih elemenata u prvom stupcu (po dijagonalnoj dominantnosti u prvom stupcu), posebno je veći od svakog drugog elementa po apsolutnoj vrijednosti u prvom stupcu, pa je to ujedno pivotni element.

Neka je \(\hat{A}\) matrica dobivena jednim korakom LU faktorizacije (svi elementi u prvom stupcu ispod dijagonale Gaussovim su eliminacijama poništeni). Dokažimo da je matrica \(\hat{A}(2:n,2:n)\) ponovno dijagonalno dominantna po stupcima. To će induktivno dovršiti dokaz prve tvrdnje. Naime, za matricu \(\hat{A}(2:n,2:n)\) ponovno vrijedi da je element na mjestu \((1,1)\) te matrice različit od nule najveći po apsolutnoj vrijednosti u tom stupcu, pa može biti pivotni element kojim nastavljamo poništavati elemente ispod njega, a preostala matrica je ponovno dijagonalno dominantna, itd.

Odredimo vrijednosti \(\hat{a}_{i,j}\) matrice \(\hat{A}\) u ovisnosti o elementima \(a_{i,j}\) matrice \(A\). Prvi redak matrice \(\hat{A}\) je isti kao prvi redak matrice \(A\), dok su svi elementi prvog stupca matrice \(\hat{A}\) osim prvog jednaki nula. Element \(\hat{a}_{i,j}\) dobivamo tako da elementu \(a_{i,j}\) oduzmemo element \(a_{1,j}\) iz prvog retka pomnožen s \(\dfrac{a_{i,1}}{a_{1,1}}\):

\[\hat{a}_{i,j} = a_{i,j} - \frac{a_{i,1}}{a_{1,1}} a_{1,j}.\]

Potrebno je dokazati dijagonalnu dominantnost matrice \(\hat{A}(2:n,2:n)\) koja je ekvivalentna uvjetu

\[\begin{split}\left| \hat{a}_{j,j} \right| = \left| a_{j,j} - \frac{a_{j,1}}{a_{1,1}} a_{1,j} \right| \geq \sum_{\substack{i=2 \\ i\not= j}}^n \left| a_{i,j} - \frac{a_{i,1}}{a_{1,1}} a_{1,j} \right| = \sum_{\substack{i=2 \\ i\not= j}}^n \left| \hat{a}_{i,j} \right|\end{split}\]

za svaki \(j\geq 2\). Primijetimo da u gornjoj sumi nema sumanada za \(i=1\) budući da gledamo podmatricu matrice \(\hat{A}\).

Iz dijagonalne dominantnosti za matricu \(A\) za svaki \(j\geq 2\) slijedi

\[\begin{split}|a_{j,j}| \geq \sum_{\substack{i=1 \\ i\not=j}}^{n} | a_{i,j}|, \text{ odnosno } |a_{j,j}| - |a_{1,j}| \geq \sum_{\substack{i=2 \\ i\not= j}}^n | a_{i,j}|.\end{split}\]

Sada imamo

\[\begin{split}\begin{align*} \sum_{\substack{i=2 \\ i\not= j}}^n \left| \hat{a}_{i,j} \right| &= \sum_{\substack{i=2 \\ i\not=j}} \left| a_{i,j} - \frac{a_{i,1}}{a_{1,1}} a_{1,j} \right| \\ &\leq \sum_{\substack{i=2 \\ i\not=j}} \left| a_{i,j}\right| + \sum_{\substack{i=2 \\ i\not=j}} \left|\frac{a_{i,1}}{a_{1,1}} a_{1,j} \right| \\ &= \sum_{\substack{i=2 \\ i\not=j}} \left| a_{i,j}\right| + \left|\frac{a_{1,j}}{a_{1,1}}\right| \sum_{\substack{i=2 \\ i\not=j}} |a_{i,1}| \\ &\leq ( |a_{j,j}| - |a_{1,j}| ) + \frac{|a_{1,j}|}{|a_{1,1}|} ( |a_{1,1}| - |a_{j,1}| ) \\ & = |a_{j,j}| - |a_{1,j}| + |a_{1,j}| - \frac{|a_{1,j}|}{|a_{1,1}|} |a_{j,1}| \\ & \leq \left| a_{j,j} - \frac{a_{1,j}}{a_{1,1}} a_{j,1} \right| = \left| \hat{a}_{j,j} \right| \end{align*}\end{split}\]

(prva nejednakost i zadnja nejednakost su nejednakost trokuta, a druga nejednakost je dijagonalna dominatnost za matricu \(A\)). Dakle, preostala matrica je dijagonalno dominantna, čime induktivno zaključujemo tvrdnju.

Ako je matrica dijagonalno dominantna po retcima, ponovno lako pokažemo da je element na mjestu \((1,1)\) različit od nule jer bi zbog dijagonalne dominantnosti jedino bilo moguće da je cijeli prvi red jednak nuli, što se kosi s regularnošću matrice \(A\). Nakon provođenja prvog koraka LU fakorizacije na isti način se pokaže da je matrica \(\hat{A}(2:n,2:n)\) ponovno dijagonalno dominantna po retcima, što na isti način dovršava tvrdnju.

Posebno, vrpčaste dijagonalno dominantne matrice po stupcima ne treba pivotirati, pa onda LU faktorizacija sigurno čuva njihovu strukturu.

Sada promotrimo zadnju familiju matrica. Možemo se pitati mogu li faktori \(L\) i \(U\) u LU faktorizaciji biti jednaki (nakon transponiranja). Uvjet da su na dijagonali matrice \(L\) nužno jedinice zanemarujemo. Dakle, tražimo regularne kvadratne matrice \(A\) reda \(n\) za koje postoji gornjetrokutasta matrica \(R\) reda \(n\) takva da je \(A = R^TR\).

Prvo što možemo očito zaključiti je da je \(A\) simetrična. Drugo, što je možda malo teže zaključiti, je da je matrica \(A\) pozitivno definitna.

Definicija 2.11

Kvadratna simetrična matrica \(A \in \R^{n \times n} \) je pozitivno definitna ako za svaki vektor \(x \in \R^n\), \(x\not = 0\) vrijedi \(x^T A x > 0\).

Poznate su karakterizacije pozitivne definitnosti, navodimo ih bez dokaza.

Propozicija 2.12

• Simetrična matrica je pozitivno definitna ako i samo ako joj svaka vodeća podmatrica \(A(1:k,1:k)\) ima pozitivnu determinantu, pri čemu je \(k=1,2,\dots,n\).

• Simetrična matrica je pozitivno definitna ako i samo ako su joj sve svojstvene vrijednosti pozitivne.

Ono što možemo dokazati su neki nužni ili samo dovoljni uvjeti za pozitivnu definitnost matrice.

Propozicija 2.13

Simetrična strogo dijagonalno dominantna matrica s pozitivnim dijagonalnim elementima je pozitivno definitna.

Budući da je matrica simetrična, dijagonalna dominantnost po retcima ekvivalentna je dijagonalnoj dominantnosti po stupcima, pa zato to ne navodimo eksplicitno u iskazu gornjeg rezultata.

Dokaz.

Dokazat ćemo po definiciji pozitivne definitnosti. Neka je \(x\) proizvoljan vektor različit od nule, te \(A\in \R^{n \times n}\) matrica koja zadovoljava uvjete propozicije. Svaku nejednakost iz definicije stroge pozitivne definitnosti

\[\begin{split}a_{i,i} > \sum_{\substack{j=1 \\ j\not=i}}^{n} | a_{i,j}|\end{split}\]

množimo s \(x_i^2\) i zbrojimo. Budući da neka od koordinata \(x_j\) može biti jednaka nula, u nekoj od pomnoženih nejednakosti možemo izgubiti strogu nejednakost, no kako nisu svi \(x_i\) jednaki nula, u zbroju svih tih nejednakosti i dalje vrijedi stroga nejednakost. Dakle:

\[\begin{split}\sum_{i=1}^n a_{i,i} x_i^2 > \sum_{i=1}^n \sum_{\substack{j=1 \\ j\not=i}}^{n} | a_{i,j}|x_i^2.\end{split}\]

Sada raspišimo uvjet pozitivne definitnosti:

\[\begin{split}x^T A x &= \sum_{i=1}^n x_i (Ax)_{i} = \sum_{i=1}^n x_i \left( \sum_{j=1}^n a_{i,j} x_j \right)\\ &= \sum_{i=1}^n \sum_{j=1}^n a_{i,j}x_ix_j\\ &= \sum_{i=1}^n a_{i,i} x_i^2 + \sum_{i=1}^n \sum_{\substack{j=1 \\ j\not=i}}^{n} a_{i,j}x_ix_j\\ &> \sum_{i=1}^n \sum_{\substack{j=1 \\ j\not=i}}^{n} | a_{i,j}|x_i^2 + \sum_{i=1}^n \sum_{\substack{j=1 \\ j\not=i}}^{n} a_{i,j}x_ix_j\\ &= \sum_{1\leq i < j \leq n} |a_{i,j}| (x_i^2+x_j^2) + \sum_{1\leq i < j \leq n}2a_{i,j}x_ix_j \\ &= \sum_{1\leq i < j \leq n} |a_{i,j}| (x_i \pm x_j)^2 \geq 0\end{split}\]

(gdje je znak \(\pm\) za par \((i,j)\) jednak plusu ako je \(a_{i,j} \geq 0\), a minusu inače), što je i trebalo dokazati.

Propozicija 2.14

Neka je matrica \(A \in \R^{n \times n}\) simetrična pozitivno definitna. Tada su joj svi dijagonalni elementi \(a_{i,i}\) pozitivni.

Dokaz.

U definiciju pozitivne definitnosti uvrstimo kanonski vektor \(e_i\), za svaki \(i=1,\dots,n\). Dobivamo \(a_{i,i} = e_i^T Ae_i > 0 \), što je i trebalo dokazati.

Vratimo se na faktorizaciju \(A=R^TR\). Ako zaista vrijedi ta jednakost za neku regularnu matricu \(A\), onda za svaki vektor različit od nule imamo

\[x^TAx = x^TR^TRx = \| Rx \|_2^2 \geq 0.\]

Kada bi za neki vektor \(x \not= 0\) vrijedila jednakost \(\| Rx \| = 0\), to bi značilo da je matrica \(R\) singularna, a to bi onda značilo da je i matrica \(A\) singularna, čime dobivamo kontradikciju. Dakle, nužno je matrica pozitivno definitna.

Vrijedi i obrat koji ćemo dokazati u sljedećem teoremu.

Teorem 2.15

Neka je \(A \in \R^{n \times n}\) pozitivno definitna matrica. Tada postoje jedinstvena donjetrokutasta matrica s jedinicama na dijagonali \(L\), jedinstvena dijagonalna matrica \(D\) i jedinstvena gornjetrokutasta matrica \(R\) s pozitivnim vrijednostima na dijagonali (sve reda \(n\)) takve da vrijedi

\[A = LDL^T = R^TR.\]

Faktorizacija \(A=R^TR\) naziva se faktorizacija Choleskog.

Dokaz.

Matrica \(A\) je pozitivno definitna. Kako su joj po karakterizaciji pozitivne definitnosti determinante svih vodećih podmatrica \(A(1:k,1:k)\) pozitivne, te podmatrice su regularne, pa matrica \(A\) dopušta jedinstvenu \(LU\) faktorizaciju po Teoremu 2.5. Dakle, postoje matrice \(L\) i \(U\) takve da je \(A=LU\).

Prema Binet-Cauchyjevom teoremu vrijedi \(\det A = \det L \cdot \det U\), te su sve te determinante različite od nule. Kako za trokutaste matrice vrijedi da im je determinanta jednaka umnošku dijagonalnih elemenata, zaključujemo da su na dijagonalni matrice \(U\) elementi različiti od nula. Označimo s \(D\) (regularnu) dijagonalnu matricu dobivenu od dijagonalnih elemenata matrice \(U\) i definirajmo \(L_2:= U^T D^{-1}\). To je donjetrokutasta matrica (kao umnožak donjetrokutastih matrica) s jedinicama na dijagonali (množenjem \(D^{-1}\) zdesna dijelimo svaki stupac s odgovarajućim elementom matrice \(D\)).

Novom definicijom možemo pisati \(A = LU = LDL_2^T\). Kako je matrica \(A\) simetrična, iz \(A^T = L_2DL^T\) dobivamo

\[LDL_2^T = L_2DL^T.\]

Množenjem te jednakosti slijeva s \(L^{-1}\) i zdesna s \(L^{-T}\) (matrica \(L\) je regularna, pa ti inverzi postoje) dobvamo

\[DL_2^TL^{-T} = L^{-1}L_2D.\]

Na lijevoj strani nalazi se umnožak gornjetrokutastih matrica, što je gornjetrokutasta matrica. Na desnoj strani nalazi se umnožak donjetrokutastih matrica, što je donjetrokutasta matrica. To je jedino moguće samo ako su na obje strane dijagonalne matrice. Dodatno, kako se na dijagonalama svih matrica osim \(D\) nalaze jedinice, uspoređivanjem dijagonalnih elemenata dobivamo da su te matrice upravo jednake \(D\). Dakle, dobili smo

\[L^{-1}L_2D = D.\]

Množenjem slijeva s \(L\), i zdesna s \(D^{-1}\) slijedi \(L_2 = L\). Dakle, matricu \(A\) zaista možemo pisati kao \(A = LDL^T\).

Matrica \(A\) je pozitivno definitna, pa za svaki \(x\not= 0\) vrijedi

\[0< x^T A x = x^T LDL^T x.\]

Uvrstimo \(x = L^{-T} e_i\), gdje su \(e_i\) kanonski vektori za \(\R^n\). Dobivamo

\[0< e_i^T L^{-1} LDL^T L^{-T} e_i = e_i^T D e_i.\]

Zadnji broj upravo je \(i\)-ti element dijagonale \(D\). Ova nejednakost pokazuje da su elementi na dijagonali matrice \(D\) pozitivne, pa postoji dijagonalna matrica \(E\) kojoj su elementi jednaki drugim korijenima matrice \(D\). Posebno, vrijedi \(E^2 = D\). Tada definiramo gornjetrokutastu matricu \(R:=EL^T\), pa za nju vrijedi

\[A = LDL^T = LE^2L^T = (EL^T)^T (EL^T) = R^TR,\]

što je i trebalo dokazati. Jedinstvenost matrica \(L\), \(D\) i \(R\) slijedi iz jedinstvenosti \(LU\) faktorizacije za matricu \(A\).

Faktorizacija Choleskog, osim što se vrši na simetričnim matricama, ima dodatnu prednost u tome što se znaju eksplicitne formule za elemente matrice \(R\). Direktnom provjerom jednakosti \(A = R^T R\) na mjestu \((i,j)\), uz \(i\leq j \), dobivamo

\[a_{i,j} = e_i^T R^T R e_j = (Re_i)^T \cdot (R e_j) = \sum_{k=1}^n r_{k,i} r_{k,j} = \sum_{k=1}^i r_{k,i} r_{k,j},\]

gdje smo na kraju iskoristili da je matrica \(R\) gornje trokutasta pa su ispod dijagonale svi elementi jednaki nula. Prvo promotrimo ove jednakosti za one indekse \(i\) za koje sume na desnoj strani imaju što manje pribrojnika. Za \(i=1\) imamo

\[\begin{split}\begin{align*} a_{1,1} = r_{1,1}^2 &\implies r_{1,1} = \sqrt{ a_{1,1} } \\ \text{za } j>1: \ a_{1,j} = r_{1,1} r_{1,j} &\implies r_{1,j} = \frac{1}{r_{1,1}} a_{1,j}. \end{align*}\end{split}\]

Za \(i=2\) imamo

\[\begin{split}\begin{align*} a_{2,2} = r_{1,1}^2 + r_{2,2}^2 &\implies r_{2,2} = \sqrt{ a_{2,2} - r_{1,2}^2 } \\ \text{za } j>1: \ a_{2,j} = r_{1,2} r_{1,j} + r_{2,2} r_{2,j} &\implies r_{2,j} = \frac{1}{r_{2,2}} \left( a_{2,j} - r_{1,2} r_{1,j} \right). \end{align*}\end{split}\]

Za \(i=3\) imamo

\[\begin{split}\begin{align*} a_{3,3} = r_{1,1}^2 + r_{2,2}^2 + r_{3,3}^2 &\implies r_{3,3} = \sqrt{ a_{3,3} - r_{1,3}^2 -r_{2,3}^2 } \\ \text{za } j>1: \ a_{3,j} = r_{1,3} r_{1,j} + r_{2,3} r_{2,j} + r_{3,3} r_{3,j} &\implies r_{3,j} = \frac{1}{r_{3,3}} \left( a_{3,j} - r_{1,3} r_{1,j} - r_{2,3} r_{2,j} \right). \end{align*}\end{split}\]

Ovaj proces možemo nastaviti za sve veće \(i\). Dolazimo do algoritma koji pišemo u Pythonu:

# Faktorizacija Choleskog simetrične poz. def. matrice A 
def FaktorizacijaCholeskog(A):
    n = A.shape[0]
    R = np.zeros((n, n))

    for i in range(n):
        # Dijagonala od R.
        izraz = A[i, i] - np.sum(R[:i, i] ** 2)
        if izraz <= 0:
            return np.zeros((n,n))
        else:
            R[i, i] = np.sqrt(izraz)

        # Elementi iznad dijagonale od R.
        for j in range(i + 1, n):
            R[i,j] = (A[i,j] - np.sum(R[:i, j] * R[:i, i])) / R[i, i]

    return R

Testirano na konkretnom primjeru:

A = np.array([[1, -2, 1], [-2, 8, -8], [1, -8, 19]])
R = FaktorizacijaCholeskog(A)

print( f'Matrica A:\n {A}\n' )
print( f'Matrica R:\n {R}\n' )
print( f'Matrica R^T*R:\n {R.T @ R}\n' )

Matrica A:
 [[ 1 -2  1]
 [-2  8 -8]
 [ 1 -8 19]]

Matrica R:
 [[ 1. -2.  1.]
 [ 0.  2. -3.]
 [ 0.  0.  3.]]

Matrica R^T*R:
 [[ 1. -2.  1.]
 [-2.  8. -8.]
 [ 1. -8. 19.]]

U formulama na dva mjesta dobivamo mogući problem s računom: za dijagonalne elemente potrebno je izračunati određeni drugi korijen, dok za vandijagonalne elemente dijelimo s elementom \(r_{i,i}\). Teorem 2.15 nam garantira da pozitivno definitna matrica ima faktorizaciju Choleskog, pa nam to zapravo garantira da nećemo imati računskih problema, odnosno da ćemo gornji algoritam moći provesti. Vrijedi i obrat: ako matrica nije pozitivno definitna, znamo da nema faktorizaciju Choleskog. Kada bismo za takvu matricu mogli uspješno provesti gornji algoritam, to bi značilo da ona ima faktorizaciju Choleskog, time mora biti pozitivno definitna, čime dolazimo do neke kontradikcije. Poanta je: uspješno provođenje gore opisanog algoritma za račun matrice \(R\) iz faktorizacije Choleskog je nužan i dovoljan uvjet da şimetrična matrica bude pozitivno definitna.

Spomenimo još da je u Choleskyjevoj faktorizaciji ponekad potrebno pivotirati iz razloga sličnih pivotiranju u klasičnoj LU faktorizaciji. U tom slučaju pivotiranje se radi simetrično kako bi se održala simetrija matrice \(A\), no u ovom kolegiju to nećemo detaljnije analizirati.