LU-faktorisering og Gauss-eliminasjon

LU-faktorisering og Gauss-eliminasjon#

TMA4400 Matematikk 1: Kalkulus og lineær algebra

Dato: 01. september 2025

Eksempel#

Vi kan bruke standard Gauss-eliminasjon på følgende utvidede matrise

\[\begin{split}(A|\mathbf{b}) = \left[\begin{array}{ccc|c} 1 & 1 & 1 & 1 \\ 4 & 3 & -1 & 1 \\ 3 & 5 & 3 & 1 \end{array}\right]\end{split}\]

for å få

\[\begin{split}(A|\mathbf{b}) = \left[\begin{array}{ccc|c} 1 & 1 & 1 & 1 \\ 0 & -1 & -5 & -3 \\ 0 & 0 & -10 & -8 \end{array}\right].\end{split}\]

Om vi gjør det samme på

\[\begin{split}(A|\mathbf{b}) = \left[\begin{array}{ccc|c} 1 & 1 & 1 & 1 \\ 4 & 3 & -1 & 2 \\ 3 & 5 & 3 & 3 \end{array}\right]\end{split}\]

får vi

\[\begin{split}(A|\mathbf{b}) = \left[\begin{array}{ccc|c} 1 & 1 & 1 & 1 \\ 0 & -1 & -5 & -2 \\ 0 & 0 & -10 & -4 \end{array}\right].\end{split}\]

Matrisa \(A\) får altså den samme radreduserte formen, men vi må gjøre Gauss-elminiasjonen på nytt for hver av de to utvidede matrisene.

Finnes det en måte å lagre Gauss-eliminasjonen av \(A\) slik at vi ikke trenger å gjøre Gauss-elminiasjonen på nytt for hele systemet?

Først kan vi legge merke til at den radreduserte matrisa som svarer til \(A\),

\[\begin{split}\left[\begin{array}{ccc} 1 & 1 & 1 \\ 0 & -1 & -5 \\ 0 & 0 & -10 \end{array}\right],\end{split}\]

er ei øvre triangulær matrise (elementene under diagonalen er \(0\)). Dette gjør at ligningssystemet til den radreduserte \((A|\mathbf{b})\) er veldig enkelt å løse siden vi raskt leser ut at

\[\begin{align*} x_3&=\frac{-4}{-10}=\frac{2}{5},\\ x_2&=2-5\cdot x_3=0,\\ x_1&=1-x_2-x_3=1-0-\frac{2}{5}=\frac{3}{5}. \end{align*}\]

På samme måte er ligningssystemer i form av nedre triangulære matriser (elementene over diagonalen er \(0\)) også lette å løse.

LU-faktorisering#

La \(A\) være ei generell kvadratisk \(n\times n\)-matrise, \(\mathbf{x}\) en vektor i \(\mathbb{R}^n\) og \(\mathbf{b}\) en vektor i \(\mathbb{R}^n\). Vi skal nå løse følgende problem:

\[ \text{Finn $\mathbf{x}$ slik at $A\mathbf{x}=\mathbf{b}$.} \]

Anta at vi kan skrive \(A=LU\), der \(n\times n\)-matrisa \(L\) er nedre triangulær og \(n\times n\)-matrisa \(U\) er øvre triangulær. Da kan vi omformulere problemet på følgende måte:

\[\begin{align*} \text{Steg $1.$}&\quad \text{Finn $\mathbf{y}$ slik at $L\mathbf{y}=\mathbf{b}$.}\\ \text{Steg $2.$}&\quad \text{Finn $\mathbf{x}$ slik at $U\mathbf{x}=\mathbf{y}$.}\\ \end{align*}\]

Vi har altså byttet ut det å radredusere den utvidede matrisa \((A|\mathbf{b})\) med å løse to enkle triangulære ligningssystem \(L\mathbf{y}=\mathbf{b}\) og \(U\mathbf{x}=\mathbf{y}\). Fordelen er at vi nå ikke trenger å radredusere \((A|\mathbf{b})\) hver gang vi får gitt en ny \(\mathbf{b}\).

Men hvordan finner vi \(L\) og \(U\)? For å finne \(U\) velger vi rett og slett den radreduserte matrisa til \(A\).

import numpy as np

# Gauss-eliminasjon (antar at pivotelementene er forskjellig fra null i alle rader!)
def GE(B):
    A = np.copy(B) # vi lager en kopi av den gitte matrisa B, for å unngå å endre denne
    n,n = np.shape(B) # n=rader, n=kolonner; kvadratisk matrise
    for k in range(0,n-1):
        for j in range(k+1,n):
            l = A[j,k]/A[k,k]
            for p in range(k,n):
                A[j,p] = A[j,p]-l*A[k,p]
    U = A
    return U

A0 = np.array([[1.,1.,1.],[4.,3.,-1.],[3.,5.,3.]])
U0 = GE(A0)
print('\nVi begynner med \n{0}\nog bruker Gauss-eliminasjon som gitt over: \n{1}'.format(np.matrix(A0),np.matrix(U0)))

Vi begynner med 
[[ 1.  1.  1.]
 [ 4.  3. -1.]
 [ 3.  5.  3.]]
og bruker Gauss-eliminasjon som gitt over: 
[[  1.   1.   1.]
 [  0.  -1.  -5.]
 [  0.   0. -10.]]

For å få en entydig bestemt måte å skrive \(A=LU\) velger vi at \(L\) skal ha en diagonal med bare \(1\)-ere. Det gjenstår altså å bestemme elementene som befinner seg under diagonalen.

Hvis vi tar utgangspunkt i eksempelet over, skal vi altså finne

\[\begin{split}L=\left[\begin{array}{ccc} 1 & 0 & 0 \\ l_{21} & 1 & 0 \\ l_{31} & l_{32} & 1 \end{array}\right].\end{split}\]

Legg merke til at \(L\) er slik at \(A=LU\), og \(U\) er den radreduserte matrisa. Da må \(L\) sine gjenværende elementer ha en verdi som kommer fra radreduseringen. Vi modifiserer koden på følgende måte:

import numpy as np

# Gauss-eliminasjon (antar at pivotelementene er forskjellig fra null i alle rader!)
def GE(B):
    A = np.copy(B)  # vi lager en kopi av den gitte matrisa B, for å unngå å endre denne
    n, n = np.shape(B) # n=rader, n=kolonner; kvadratisk matrise
    L = np.eye(n)   # L begynner som identitetsmatrisa, og så fyller vi inn elementene under diagonalen
    for k in range(n-1):
        for j in range(k+1,n):
            L[j,k] = A[j,k]/A[k,k]  # vi lagrer multipliseringen i L
            for p in range(k,n):
                A[j,p] = A[j,p]-L[j,k]*A[k,p]
    U = A
    return L, U

A1 = np.array([[1.,1.,1.],[4.,3.,-1.],[3.,5.,3.]])
L1, U1 = GE(A1)
print('\nVi begynner med \n{0}\nog ender opp med følgende matriser: \n{1}\n{2}'.format(np.matrix(A1),np.matrix(L1),np.matrix(U1)))

Vi begynner med 
[[ 1.  1.  1.]
 [ 4.  3. -1.]
 [ 3.  5.  3.]]
og ender opp med følgende matriser: 
[[ 1.  0.  0.]
 [ 4.  1.  0.]
 [ 3. -2.  1.]]
[[  1.   1.   1.]
 [  0.  -1.  -5.]
 [  0.   0. -10.]]

Vi har altså at

\[\begin{split}L=\left[\begin{array}{ccc} 1 & 0 & 0 \\ 4 & 1 & 0 \\ 3 & -2 & 1 \end{array}\right].\end{split}\]

La oss teste dette:

L2 = np.array([[1.,0.,0.],[4.,1.,0.],[3.,-2.,1.]])
U2 = np.array([[1.,1.,1.], [0.,-1.,-5.],[0.,0.,-10.]])
A2 = np.dot(L2,U2) # notasjonen for matrisemultiplikasjon i python
print(A2)

[[ 1.  1.  1.]
 [ 4.  3. -1.]
 [ 3.  5.  3.]]

Numerisk ustabilitet#

La oss gjøre \(LU\)-faktorisering på følgende matrise:

\[\begin{split}A = \left[\begin{array}{ccc} 10^{-12} & 4 & 1 \\ 2 & -1 & -2 \\ 1 & 3 & 2 \end{array}\right].\end{split}\]
A3 = np.array([[1.e-12,4.,1.],[2.,-1.,-2.],[1.,3.,2.]])
L3, U3 = GE(A3)
print('\nVi begynner med \n{0}\nog ender opp med følgende matriser: \n{1}\n{2}'.format(np.matrix(A3),np.matrix(L3),np.matrix(U3)))

Vi begynner med 
[[ 1.e-12  4.e+00  1.e+00]
 [ 2.e+00 -1.e+00 -2.e+00]
 [ 1.e+00  3.e+00  2.e+00]]
og ender opp med følgende matriser: 
[[1.e+00 0.e+00 0.e+00]
 [2.e+12 1.e+00 0.e+00]
 [1.e+12 5.e-01 1.e+00]]
[[ 1.000e-12  4.000e+00  1.000e+00]
 [ 0.000e+00 -8.000e+12 -2.000e+12]
 [ 0.000e+00  0.000e+00  2.125e+00]]

Vi tester så om vi får riktig produkt:

L4 = np.array([[1.,0.,0.],[2.e+12,1.,0.],[1.e+12,0.5,1.]])
U4 = np.array([[1.e-12,4.,1.], [0.,-8.e+12,-2.e+12],[0.,0.,2.125]])
A4 = np.dot(L4,U4) # notasjonen for matrisemultiplikasjon i python
print(A4)

[[1.000e-12 4.000e+00 1.000e+00]
 [2.000e+00 0.000e+00 0.000e+00]
 [1.000e+00 0.000e+00 2.125e+00]]

Vi ser at matrisa vi ender opp med,

\[\begin{split}A = \left[\begin{array}{ccc} 10^{-12} & 4 & 1 \\ 2 & 0 & 0 \\ 1 & 0 & 17/8 \end{array}\right],\end{split}\]

ikke er den samme som vi begynte med. Vi trenger altså en forbedret algoritme.

LU-faktorisering med delvis pivotering#

Vi tar utgangspunkt i det vi gjorde med Gauss-eliminasjon, nemlig delvis pivotering. Vi gjenbruker koden og gjør de nødvendige endringene.

# Gauss-eliminasjon med delvis pivotering
def GEPP(B):
    A = np.copy(B).astype(float) # vi lager en kopi av den gitte matrisa B, for å unngå å endre denne. i tillegg konverterer vi til desimaltall
#   A = np.copy(B) # prøv gjerne uten for å se hva som skjer 
    n,n = np.shape(A) # n=rader, n=kolonner; kvadratisk matrise
    L = np.eye(n)   # L begynner som identitetsmatrisa, og så fyller vi inn elementene under diagonalen
    for k in range(0,n-1):
        a = np.abs(A[k,k])
        pk = k
        for j in range(k+1,n):
            if a < np.abs(A[j,k]):
                a = np.abs(A[j,k])
                pk = j
        if pk != k:
            for p in range(k,n): # vi må bytte radene i A
                r = A[k,p] 
                A[k,p] = A[pk,p]
                A[pk,p] = r
            for p in range(0,k): # vi må bytte de aktuelle radelementene i L
                r = L[k,p] 
                L[k,p] = L[pk,p]
                L[pk,p] = r
        for j in range(k+1,n):
            L[j,k] = A[j,k]/A[k,k]
            for p in range(k,n):
                A[j,p] = A[j,p]-L[j,k]*A[k,p]
    U = A
    return L, U

Igjen prøver vi oss på \(LU\)-faktorisering av følgende matrise:

\[\begin{split}A = \left[\begin{array}{ccc} 10^{-12} & 4 & 1 \\ 2 & -1 & -2 \\ 1 & 3 & 2 \end{array}\right].\end{split}\]
A5 = np.array([[1.e-12,4.,1.],[2.,-1.,-2.],[1.,3.,2.]])
L5, U5 = GEPP(A5)
print('\nVi begynner med \n{0}\nog ender opp med følgende matriser: \n{1}\n{2}'.format(np.matrix(A5),np.matrix(L5),np.matrix(U5)))

Vi begynner med 
[[ 1.e-12  4.e+00  1.e+00]
 [ 2.e+00 -1.e+00 -2.e+00]
 [ 1.e+00  3.e+00  2.e+00]]
og ender opp med følgende matriser: 
[[1.00e+00 0.00e+00 0.00e+00]
 [5.00e-13 1.00e+00 0.00e+00]
 [5.00e-01 8.75e-01 1.00e+00]]
[[ 2.    -1.    -2.   ]
 [ 0.     4.     1.   ]
 [ 0.     0.     2.125]]

Vi tester igjen om vi får riktig produkt:

L6 = np.array([[1.,0.,0.],[5.e-13,1.,0.],[1/2.,0.875,1.]])
U6 = np.array([[2.,-1.,-2.], [0.,4.,1.],[0.,0.,17/8.]])
A6 = np.dot(L6,U6) # notasjonen for matrisemultiplikasjon i python
print(A6)

[[ 2.e+00 -1.e+00 -2.e+00]
 [ 1.e-12  4.e+00  1.e+00]
 [ 1.e+00  3.e+00  2.e+00]]

Vi får nå matrisa vi begynte med, men med rad 1 og 2 byttet om. Hvis vi hadde lagret permutasjonene som matrisa

\[\begin{split}P = \left[\begin{array}{ccc} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{array}\right]\end{split}\]

ville vi fått at \(PA=LU\). Dette går vi ikke videre inn på her.

Contents