LU-decompositie, Cholesky en QR-decompositie

Op de site http://wiki.sagemath.org/quickref kan men een beknopte lijst vinden van de beschikbare bewerkingen voor lineaire algebra: quickref-linalg.pdf. We zullen nader kennismaken met LU(), cholesky() en QR(). Verder wordt verwacht dat je experimenteert met een aantal andere procedures.

def hilbert_matrix(n): return matrix([[1/i for i in [j..j+(n-1)]] for j in [1..n]]) H = hilbert_matrix(6); show(H)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} \\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} \\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} \\ \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} \\ \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} \\ \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} & \frac{1}{10} & \frac{1}{11} \end{array}\right)$

Help voor LU() vind je door de volgende opdracht uit te voeren. Je vindt er informatie over het aantal argumenten, mogelijke opties, etc.

H.LU?

We voeren de LU decompositie uit op $H$:

H.LU()

(
[1 0 0 0 0 0]  [    1     0     0     0     0     0]
[0 1 0 0 0 0]  [  1/2     1     0     0     0     0]
[0 0 0 1 0 0]  [  1/6   5/7     1     0     0     0]
[0 0 0 0 1 0]  [  1/3     1 14/25     1     0     0]
[0 0 0 0 0 1]  [  1/4  9/10 21/25   6/7     1     0]
[0 0 1 0 0 0], [  1/5   4/5 24/25   3/7   8/9     1],

[         1        1/2        1/3        1/4        1/5        1/6]
[         0       1/12       1/12       3/40       1/15       5/84]
[         0          0      5/504       1/63      2/105   100/4851]
[         0          0          0    -1/1800     -1/875     -1/616]
[         0          0          0          0   -1/49000   -1/19404]
[         0          0          0          0          0 -1/1746360]
)

(
[1 0 0 0 0 0]  [    1     0     0     0     0     0]
[0 1 0 0 0 0]  [  1/2     1     0     0     0     0]
[0 0 0 1 0 0]  [  1/6   5/7     1     0     0     0]
[0 0 0 0 1 0]  [  1/3     1 14/25     1     0     0]
[0 0 0 0 0 1]  [  1/4  9/10 21/25   6/7     1     0]
[0 0 1 0 0 0], [  1/5   4/5 24/25   3/7   8/9     1],

[         1        1/2        1/3        1/4        1/5        1/6]
[         0       1/12       1/12       3/40       1/15       5/84]
[         0          0      5/504       1/63      2/105   100/4851]
[         0          0          0    -1/1800     -1/875     -1/616]
[         0          0          0          0   -1/49000   -1/19404]
[         0          0          0          0          0 -1/1746360]
)

Na uitvoering verschijnen er drie matrices, $p$, $l$ en $u$, met $H = p \cdot l \cdot u$. Hierbij is $l$ de benedentriangulaire matrix (met 1 op de diagonaal) en $u$ de boventriangulaire matrix; $p$ is de permutatiematrix die nodig is bij het pivoteren. Ter verificatie:

[Hp, Hl, Hu] = H.LU() show(Hp*Hl*Hu) Hp*Hl*Hu == H

True

True

Als rekenhulp bij oefening E2.1 (zie pakket extra oefeningen) kan je de LU-decompositie van de volgende matrix A beschouwen:

var('a'); A = matrix([[1,a,a],[a,1,a],[a,a,1]]); show(A)

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1 & a & a \\ a & 1 & a \\ a & a & 1 \end{array}\right)$

A.LU()

Traceback (click to the left of this block for traceback)
...
TypeError: base ring of the matrix must be exact, not Symbolic Ring

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_7.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("QS5MVSgp"),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/var/sage/tmpCFFI9D/___code___.py", line 2, in <module>
    exec compile(u'A.LU()
  File "", line 1, in <module>
    
  File "sage/matrix/matrix2.pyx", line 11551, in sage.matrix.matrix2.Matrix.LU (/opt/sage/sage-6.8/src/build/cythonized/sage/matrix/matrix2.c:75899)
TypeError: base ring of the matrix must be exact, not Symbolic Ring

De error vertelt ons dat dit niet is toegelaten. De LU implementatie van Sage laat (nog) geen parameters toe. Dit is wel toegelaten in Maple, we tonen de resultaten van de decompositie uit Maple:

Al = matrix([[1,0,0],[a,1,0],[a,a/(a+1),1]]); show(Al); Au = matrix([[1,a,a],[0,1-a^2, a - a^2],[0,0,(-a-1+2*a^2)/-(a+1)]]); show(Au); show((Al*Au).simplify_rational())

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1 & 0 & 0 \\ a & 1 & 0 \\ a & \frac{a}{a + 1} & 1 \end{array}\right)$

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1 & a & a \\ 0 & -a^{2} + 1 & -a^{2} + a \\ 0 & 0 & -\frac{2 \, a^{2} - a - 1}{a + 1} \end{array}\right)$

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1 & a & a \\ a & 1 & a \\ a & a & 1 \end{array}\right)$

We bekijken een voorbeeld waarbij de matrix pivotelementen heeft die nul worden:

B = matrix([[0,0,1],[0,1,1],[1,1,1]]) show(B) B.LU()

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 0 & 0 & 1 \\ 0 & 1 & 1 \\ 1 & 1 & 1 \end{array}\right)$

(
[0 0 1]  [1 0 0]  [1 1 1]
[0 1 0]  [0 1 0]  [0 1 1]
[1 0 0], [0 0 1], [0 0 1]
)

(
[0 0 1]  [1 0 0]  [1 1 1]
[0 1 0]  [0 1 0]  [0 1 1]
[1 0 0], [0 0 1], [0 0 1]
)

Wat gebeurt er als de matrix singulier is?

C=matrix([[0,1,0,1],[0,1,0,1],[0,0,0,1],[0,0,0,1]]); show(C); C.LU()

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} 0 & 1 & 0 & 1 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 \end{array}\right)$

(
[1 0 0 0]  [1 0 0 0]  [0 1 0 1]
[0 1 0 0]  [0 1 0 0]  [0 1 0 1]
[0 0 1 0]  [0 0 1 0]  [0 0 0 1]
[0 0 0 1], [0 0 0 1], [0 0 0 1]
)

(
[1 0 0 0]  [1 0 0 0]  [0 1 0 1]
[0 1 0 0]  [0 1 0 0]  [0 1 0 1]
[0 0 1 0]  [0 0 1 0]  [0 0 0 1]
[0 0 0 1], [0 0 0 1], [0 0 0 1]
)

Vervolgens beschouwen we de procedure cholesky. Pas ze opnieuw toe op de Hilbert-matrix H en verifieer het resultaat.

G = H.cholesky(); show(G);

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 1 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{2} & 0.2886751345948129? & 0 & 0 & 0 & 0 \\ \frac{1}{3} & 0.2886751345948129? & 0.0745355992499930? & 0 & 0 & 0 \\ \frac{1}{4} & 0.2598076211353316? & 0.11180339887498949? & 0.01889822365046137? & 0 & 0 \\ \frac{1}{5} & 0.2309401076758503? & 0.1277753129999880? & 0.03779644730092273? & \frac{1}{210} & 0 \\ \frac{1}{6} & 0.2061965247105807? & 0.1330992843749875? & 0.05249506569572600? & \frac{1}{84} & 0.001196473589594301? \end{array}\right)$

Sage is er niet in geslaagd om de berekende waarden zuiver symbolisch te houden (De cholesky methode laat dat ook niet toe indien men dat eist). De geëvalueerde waarden zijn echter geen vlottende-puntgetallen maar een ander type van getal:

parent(G[3,2])

Algebraic Field

Algebraic Field

#verificatie: show(G * G.transpose())

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrrrr} 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} \\ \frac{1}{2} & 0.3333333333333334? & 0.2500000000000000? & 0.2000000000000000? & 0.1666666666666667? & 0.1428571428571429? \\ \frac{1}{3} & 0.2500000000000000? & 0.2000000000000000? & 0.1666666666666667? & 0.1428571428571429? & 0.1250000000000000? \\ \frac{1}{4} & 0.2000000000000000? & 0.1666666666666667? & 0.1428571428571429? & 0.1250000000000000? & 0.11111111111111111? \\ \frac{1}{5} & 0.1666666666666667? & 0.1428571428571429? & 0.1250000000000000? & 0.11111111111111111? & 0.10000000000000000? \\ \frac{1}{6} & 0.1428571428571429? & 0.1250000000000000? & 0.11111111111111111? & 0.10000000000000000? & 0.0909090909090909? \end{array}\right)$

Er is een verlies aan precisie opgetreden bij het berekenen van de Cholesky decompositie, we berekenen het verschil met de oorspronkelijke matrix:

norm(G * G.transpose() - H)

1.993793885322399e-20

1.993793885322399e-20

De QR-decompositie heeft opnieuw een aantal optionele argumenten. We beschouwen opnieuw een voorbeeld en passen de procedure QR toe.

A = matrix([[39,-42,-3,16],[3,6,4,-18],[3,-4,4,-18],[31,-28,3,64]]); show(A); #informatie over de QR methode: A.QR?

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} 39 & -42 & -3 & 16 \\ 3 & 6 & 4 & -18 \\ 3 & -4 & 4 & -18 \\ 31 & -28 & 3 & 64 \end{array}\right)$

We verifiëren: A == Aq * Ar ?

A = matrix([[39,-42,-3,16],[3,6,4,-18],[3,-4,4,-18],[31,-28,3,64]]); [Aq , Ar] = A.QR(); show(Aq); show(Ar);

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} \frac{39}{50} & -\frac{3}{10} & -\frac{3}{10} & -\frac{23}{50} \\ \frac{3}{50} & \frac{9}{10} & -\frac{1}{10} & -\frac{21}{50} \\ \frac{3}{50} & -\frac{1}{10} & \frac{9}{10} & -\frac{21}{50} \\ \frac{31}{50} & \frac{3}{10} & \frac{3}{10} & \frac{33}{50} \end{array}\right)$

$\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} 50 & -50 & 0 & 50 \\ 0 & 10 & 5 & 0 \\ 0 & 0 & 5 & 0 \\ 0 & 0 & 0 & 50 \end{array}\right)$

We verifiëren of Aq orthogonaal is:

Aq * Aq.transpose()

[1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1]

[1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 1]

Experimenteer nu met een aantal andere procedures.

[NumAn] PC-Oef4

2274 days ago by Niels.Neirynck

LU-decompositie, Cholesky en QR-decompositie