Häufig ist in der numerischen Mathematik das Gleichungssystem
A
x
=
b
{\displaystyle Ax=b}
A
∈
R
n
×
n
{\displaystyle A\in \mathbb {R} ^{n\times n}}
Bearbeiten
In einem solchen Fall ist es ineffizient , für jede neue rechte Seiten wieder den gesamten Gauß-Algorithmus durchzuführen, da er bei jedem Durchlauf in Bezug auf
A
{\displaystyle A}
A
{\displaystyle A}
Dieses kann in Form einer Zerlegung von
A
{\displaystyle A}
P
A
=
L
R
{\displaystyle PA=LR}
L
∈
R
n
×
n
{\displaystyle L\in \mathbb {R} ^{n\times n}}
R
∈
R
n
×
n
{\displaystyle R\in \mathbb {R} ^{n\times n}}
P
∈
R
n
×
n
{\displaystyle P\in \mathbb {R} ^{n\times n}}
Permutationsmatrix ist.
Eine solche Zerlegung bzw. Faktorisierung von
A
{\displaystyle A}
Ein Permutationsmatrix ist invertierbar und
P
−
1
{\displaystyle P^{-1}}
P
⋅
A
=
L
⋅
R
{\displaystyle P\cdot A=L\cdot R}
A
{\displaystyle A}
A
=
P
−
1
⋅
L
⋅
R
{\displaystyle A=P^{-1}\cdot L\cdot R}
Liegt eine solche Faktorisierung vor, so kann man das Gleichungssystem
A
x
=
b
{\displaystyle Ax=b}
L
R
x
=
P
b
{\displaystyle LRx=Pb}
b
∈
R
n
{\displaystyle b\in \mathbb {R} ^{n}}
L
y
=
P
b
,
R
x
=
y
{\displaystyle Ly=Pb,\quad Rx=y}
löst, wobei die Lösung
y
{\displaystyle y}
Bearbeiten
Hat man das System
A
x
=
b
{\displaystyle Ax=b}
b
{\displaystyle b}
P
A
=
L
R
{\displaystyle PA=LR}
x
{\displaystyle x}
L
y
=
P
b
,
R
x
=
y
{\displaystyle Ly=Pb,\quad Rx=y}
In der obigen Beschreibung haben wir die Notwendigkeit identifiziert, dass man für das Lösungsverfahren Zeilen vertauschen müssen. Eine einzelne Zeilvertauschung kann durch ein Matrix
Z
(
i
,
j
)
{\displaystyle Z_{(i,j)}}
Z
(
1
,
4
)
=
(
0
0
0
1
0
0
1
0
0
0
0
0
1
0
0
1
0
0
0
0
0
0
0
0
1
)
∈
R
5
×
5
{\displaystyle Z_{(1,4)}={\begin{pmatrix}0&0&0&1&0\\0&1&0&0&0\\0&0&1&0&0\\1&0&0&0&0\\0&0&0&0&1\end{pmatrix}}\in \mathbb {R} ^{5\times 5}}
Geben Sie die
Z
(
1
,
4
)
{\displaystyle Z_{(1,4)}}
Z
{\displaystyle Z}
Maxima CAS ein
A
=
(
1
2
3
4
5
0
6
0
0
0
0
0
7
0
0
0
0
0
8
0
0
0
0
0
9
)
∈
R
5
×
5
{\displaystyle A={\begin{pmatrix}1&2&3&4&5\\0&6&0&0&0\\0&0&7&0&0\\0&0&0&8&0\\0&0&0&0&9\end{pmatrix}}\in \mathbb {R} ^{5\times 5}}
und berechnen Sie das Matrizenprodukt über Z.A mit "." als Operatorsymbol für die Matrixmultiplikation.
Jede Matrix
P
∈
R
n
×
n
{\displaystyle P\in \mathbb {R} ^{n\times n}}
Permutationsmatrix .
Die folgende Matrix stellt eine Permutationsmatrix dar:
P
:=
(
0
1
0
0
0
0
1
0
1
0
0
0
0
0
0
1
)
∈
R
4
×
4
.
{\displaystyle P:={\begin{pmatrix}0&1&0&0\\0&0&1&0\\1&0&0&0\\0&0&0&1\end{pmatrix}}\in \mathbb {R} ^{4\times 4}.}
Jede bijektive Abbildung
π
:
{
1
,
.
.
.
,
n
}
→
{
1
,
.
.
.
,
n
}
{\displaystyle \pi :\{1,...,n\}\to \{1,...,n\}}
Permutation .
Zusammenhang Permutation - Permutationsmatrix
Bearbeiten
Offensichtlich ist
P
∈
R
n
×
n
{\displaystyle P\in \mathbb {R} ^{n\times n}}
π
:
{
1
,
.
.
.
,
n
}
→
{
1
,
.
.
.
,
n
}
{\displaystyle \pi :\{1,...,n\}\to \{1,...,n\}}
P
=
(
e
π
(
1
)
.
.
.
e
π
(
n
)
)
{\displaystyle P={\begin{pmatrix}e^{\pi (1)}&...&e^{\pi (n)}\end{pmatrix}}}
gilt, wobei
e
j
∈
R
n
{\displaystyle e^{j}\in \mathbb {R} ^{n}}
j
{\displaystyle j}
Beispiel - Permutation - Permutationmatrix
Bearbeiten
Die im obigen Beispiel genannte Permutationmatrix ist die durch folgende Permutation definiert
π
(
1
)
:=
3
,
π
(
2
)
:=
1
,
π
(
3
)
:=
2
,
π
(
4
)
:=
4
,
{\displaystyle \pi (1):=3,\quad \pi (2):=1,\quad \pi (3):=2,\quad \pi (4):=4,}
wobei
π
:
{
1
,
.
.
.
,
n
}
→
{
1
,
.
.
.
,
n
}
{\displaystyle \pi :\{1,...,n\}\to \{1,...,n\}}
P
=
(
e
3
e
1
e
2
e
4
)
=
(
e
π
(
1
)
e
π
(
2
)
e
π
(
3
)
e
π
(
4
)
)
.
{\displaystyle P={\begin{pmatrix}e^{3}&e^{1}&e^{2}&e^{4}\end{pmatrix}}={\begin{pmatrix}e^{\pi (1)}&e^{\pi (2)}&e^{\pi (3)}&e^{\pi (4)}\end{pmatrix}}.}
Mittels Zusammenhang von Permutation und Permutationsmatrix kann man nun die Aussagen des folgenden Satzes erschließen.
Geben Sie in Maxima die Permutationsmatrix an, die die obige Permutation durchführt.
Sei
P
∈
R
n
×
n
{\displaystyle P\in \mathbb {R} ^{n\times n}}
π
:
{
1
,
.
.
.
,
n
}
→
{
1
,
.
.
.
,
n
}
{\displaystyle \pi :\{1,...,n\}\to \{1,...,n\}}
P
{\displaystyle P}
(i)
P
{\displaystyle P}
P
−
1
=
P
T
{\displaystyle P^{-1}=P^{T}}
(ii) Es gilt
P
−
1
=
P
T
=
(
e
π
−
1
(
1
)
.
.
.
e
π
−
1
(
n
)
)
{\displaystyle P^{-1}=P^{T}={\begin{pmatrix}e^{\pi ^{-1}(1)}&...&e^{\pi ^{-1}(n)}\end{pmatrix}}}
(iii) Für
α
:=
(
α
i
)
∈
R
n
{\displaystyle \alpha :=(\alpha _{i})\in \mathbb {R} ^{n}}
P
α
=
(
α
π
−
1
(
1
)
⋮
α
π
−
1
(
n
)
)
,
α
T
P
=
(
α
π
(
1
)
.
.
.
α
π
(
n
)
)
.
{\displaystyle P\alpha ={\begin{pmatrix}\alpha _{\pi ^{-1}(1)}\\\vdots \\\alpha _{\pi ^{-1}(n)}\end{pmatrix}},\quad \alpha ^{T}P={\begin{pmatrix}\alpha _{\pi (1)}&...&\alpha _{\pi (n)}\end{pmatrix}}.}
(iv) Für jede Matrix
A
∈
R
n
×
n
{\displaystyle A\in \mathbb {R} ^{n\times n}}
(
z
i
)
T
{\displaystyle (z^{i})^{T}}
a
i
{\displaystyle a^{i}}
(
i
=
1
,
.
.
.
,
n
)
{\displaystyle (i=1,...,n)}
(3.6)
P
A
=
(
(
z
π
−
1
(
1
)
)
T
⋮
(
z
π
−
1
(
n
)
)
T
)
,
A
P
=
(
a
π
(
1
)
.
.
.
a
π
(
n
)
)
.
{\displaystyle PA={\begin{pmatrix}\left(z^{\pi ^{-1}(1)}\right)^{T}\\\vdots \\\left(z^{\pi ^{-1}(n)}\right)^{T}\end{pmatrix}},\quad AP={\begin{pmatrix}a^{\pi (1)}&...&a^{\pi (n)}\end{pmatrix}}.}
(i) folgt mit (3.5) aus
(
P
T
P
)
i
j
=
[
(
(
e
π
(
1
)
)
T
⋮
(
e
π
(
n
)
)
T
)
(
e
π
(
1
)
.
.
.
e
π
(
n
)
)
]
i
j
=
[
(
e
π
(
i
)
)
T
e
π
(
j
)
]
i
j
=
(
δ
π
(
i
)
,
π
(
j
)
)
i
j
=
(
δ
i
j
)
i
j
=
I
.
{\displaystyle \left(P^{T}P\right)_{ij}=\left[{\begin{pmatrix}\left(e^{\pi (1)}\right)^{T}\\\vdots \\\left(e^{\pi (n)}\right)^{T}\end{pmatrix}}{\begin{pmatrix}e^{\pi (1)}&...&e^{\pi (n)}\end{pmatrix}}\right]_{ij}=\left[\left(e^{\pi (i)}\right)^{T}e^{\pi (j)}\right]_{ij}=(\delta _{\pi (i),\pi (j)})_{ij}=(\delta _{ij})_{ij}=I.}
Die Behauptung in (ii) ergibt sich aus
[
(
e
π
−
1
(
1
)
.
.
.
e
π
−
1
(
n
)
)
(
e
π
(
1
)
.
.
.
e
π
(
n
)
)
]
i
j
=
∑
ℓ
=
1
n
(
e
π
−
1
(
ℓ
)
)
i
(
e
π
(
j
)
)
ℓ
=
∑
ℓ
=
1
n
δ
i
,
π
−
1
(
ℓ
)
δ
π
(
j
)
,
ℓ
{\displaystyle \left[{\begin{pmatrix}e^{\pi ^{-1}(1)}&...&e^{\pi ^{-1}(n)}\end{pmatrix}}{\begin{pmatrix}e^{\pi (1)}&...&e^{\pi (n)}\end{pmatrix}}\right]_{ij}=\sum _{\ell =1}^{n}\left(e^{\pi ^{-1}(\ell )}\right)_{i}\left(e^{\pi (j)}\right)_{\ell }=\sum _{\ell =1}^{n}\delta _{i,\pi ^{-1}(\ell )}\delta _{\pi (j),\ell }}
=
∑
k
=
1
n
δ
i
,
k
δ
π
(
j
)
,
π
(
k
)
=
δ
π
(
j
)
,
π
(
i
)
=
δ
j
i
=
δ
i
j
.
{\displaystyle =\sum _{k=1}^{n}\delta _{i,k}\delta _{\pi (j),\pi (k)}=\delta _{\pi (j),\pi (i)}=\delta _{ji}=\delta _{ij}.}
Die erste Behauptung in (iii) folgt aus
P
α
=
∑
j
=
1
n
α
j
e
π
(
j
)
=
∑
ℓ
=
1
n
α
π
−
1
(
ℓ
)
e
ℓ
=
(
α
π
−
1
(
1
)
,
.
.
.
,
α
π
−
1
(
n
)
)
T
,
{\displaystyle P\alpha =\sum _{j=1}^{n}\alpha _{j}e^{\pi (j)}=\sum _{\ell =1}^{n}\alpha _{\pi ^{-1}(\ell )}e^{\ell }=(\alpha _{\pi ^{-1}(1)},...,\alpha _{\pi ^{-1}(n)})^{T},}
die zweite folgt mit (ii) und letzterer Identität aus
α
T
P
=
(
P
T
α
)
T
=
(
∑
j
=
1
n
α
j
e
π
−
1
(
j
)
)
T
=
(
α
π
(
1
)
,
.
.
.
,
α
π
(
n
)
)
.
{\displaystyle \alpha ^{T}P=\left(P^{T}\alpha \right)^{T}=\left(\sum _{j=1}^{n}\alpha _{j}e^{\pi ^{-1}(j)}\right)^{T}=(\alpha _{\pi (1)},...,\alpha _{\pi (n)}).}
Die entsprechenden Matrixversionen in (iv) folgen analog aus
P
A
=
∑
j
=
1
n
e
π
(
j
)
(
z
j
)
T
=
∑
ℓ
=
1
n
e
ℓ
(
z
π
−
1
(
ℓ
)
=
(
(
z
π
−
1
(
1
)
)
T
⋮
(
z
π
−
1
(
n
)
)
T
)
{\displaystyle PA=\sum _{j=1}^{n}e^{\pi (j)}(z^{j})^{T}=\sum _{\ell =1}^{n}e^{\ell }(z^{\pi ^{-1}(\ell )}={\begin{pmatrix}\left(z^{\pi ^{-1}(1)}\right)^{T}\\\vdots \\\left(z^{\pi ^{-1}(n)}\right)^{T}\end{pmatrix}}}
sowie unter Verwendung dieser Beziehung und (ii) aus
A
P
=
(
P
T
A
T
)
T
=
[
P
T
(
(
a
1
)
T
⋮
(
a
n
)
T
)
]
T
=
(
(
a
π
(
1
)
)
T
⋮
(
a
π
(
n
)
)
T
)
T
.
{\displaystyle AP=\left(P^{T}A^{T}\right)^{T}=\left[P^{T}{\begin{pmatrix}(a^{1})^{T}\\\vdots \\(a^{n})^{T}\end{pmatrix}}\right]^{T}={\begin{pmatrix}(a^{\pi (1)})^{T}\\\vdots \\(a^{\pi (n)})^{T}\end{pmatrix}}^{T}.}
q.e.d.
Man beachte, dass die Indizes
π
−
1
(
1
)
,
.
.
.
,
π
−
1
(
n
)
{\displaystyle \pi ^{-1}(1),...,\pi ^{-1}(n)}
P
T
{\displaystyle P^{T}}
P
{\displaystyle P}
A
{\displaystyle A}
A
{\displaystyle A}
P
{\displaystyle P}
Seien
P
:=
(
0
0
1
1
0
0
0
1
0
)
,
A
:=
(
1
2
3
4
5
6
7
8
9
)
.
{\displaystyle P:={\begin{pmatrix}0&0&1\\1&0&0\\0&1&0\end{pmatrix}},\quad A:={\begin{pmatrix}1&2&3\\4&5&6\\7&8&9\end{pmatrix}}.}
Dann folgt
P
A
=
(
0
0
1
1
0
0
0
1
0
)
(
1
2
3
4
5
6
7
8
9
)
=
(
7
8
9
1
2
3
4
5
6
)
,
{\displaystyle PA={\begin{pmatrix}0&0&1\\1&0&0\\0&1&0\end{pmatrix}}{\begin{pmatrix}1&2&3\\4&5&6\\7&8&9\end{pmatrix}}={\begin{pmatrix}7&8&9\\1&2&3\\4&5&6\end{pmatrix}},}
A
P
=
(
1
2
3
4
5
6
7
8
9
)
(
0
0
1
1
0
0
0
1
0
)
=
(
2
3
1
5
6
4
8
9
7
)
.
{\displaystyle AP={\begin{pmatrix}1&2&3\\4&5&6\\7&8&9\end{pmatrix}}{\begin{pmatrix}0&0&1\\1&0&0\\0&1&0\end{pmatrix}}={\begin{pmatrix}2&3&1\\5&6&4\\8&9&7\end{pmatrix}}.}
In numerischen Implementierungen erfolgt die Abspeicherung einer Permutationsmatrix mit der zugehörigen Permutation
π
{\displaystyle \pi }
(
π
−
1
(
1
)
,
.
.
.
,
π
−
1
(
n
)
)
T
∈
R
n
{\displaystyle (\pi ^{-1}(1),...,\pi ^{-1}(n))^{T}\in \mathbb {R} ^{n}}
(
π
(
1
)
,
.
.
.
,
π
(
n
)
)
T
∈
R
n
{\displaystyle (\pi (1),...,\pi (n))^{T}\in \mathbb {R} ^{n}}
Eine besondere Rolle spielen Elementarpermutationen
π
:
{
1
,
.
.
.
,
n
}
→
{
1
,
.
.
.
,
n
}
{\displaystyle \pi :\{1,...,n\}\to \{1,...,n\}}
i
,
r
∈
{
1
,
.
.
.
,
n
}
{\displaystyle i,r\in \{1,...,n\}}
(3.7)
π
(
i
)
=
r
,
π
(
r
)
=
i
,
π
(
j
)
=
j
(
j
∉
{
i
,
r
}
)
.
{\displaystyle \pi (i)=r,\quad \pi (r)=i,\quad \pi (j)=j\quad (j\notin \{i,r\}).}
Hier gilt wegen
π
−
1
(
r
)
=
i
=
π
(
r
)
,
π
−
1
(
i
)
=
r
=
π
(
i
)
,
π
−
1
(
j
)
=
j
=
π
(
j
)
(
j
∉
{
i
,
r
}
)
{\displaystyle \pi ^{-1}(r)=i=\pi (r),\quad \pi ^{-1}(i)=r=\pi (i),\quad \pi ^{-1}(j)=j=\pi (j)\quad (j\notin \{i,r\})}
die Identität
π
−
1
=
π
{\displaystyle \pi ^{-1}=\pi }
π
{\displaystyle \pi }
P
∈
R
n
×
n
{\displaystyle P\in \mathbb {R} ^{n\times n}}
(3.8)
P
−
1
=
P
T
=
P
.
{\displaystyle P^{-1}=P^{T}=P.}
Wir betrachten eine weitere wichtige Klasse von Matrizen.
Sei
k
∈
{
1
,
.
.
.
,
n
}
{\displaystyle k\in \{1,...,n\}}
(3.9)
(
1
⋱
1
−
l
k
+
1
,
k
1
⋮
⋱
−
l
n
,
k
1
)
∈
R
n
×
n
{\displaystyle {\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&-l_{k+1,k}&1&&\\&&\vdots &&\ddots &\\&&-l_{n,k}&&&1\end{pmatrix}}\in \mathbb {R} ^{n\times n}}
heißt Frobenius-Matrix vom Index
k
{\displaystyle k}
.
Eine Frobenius-Matrix vom Index
k
{\displaystyle k}
k
{\displaystyle k}
k
{\displaystyle k}
k
{\displaystyle k}
z
j
∈
R
n
{\displaystyle z^{j}\in \mathbb {R} ^{n}}
(
j
=
1
,
.
.
.
,
n
)
{\displaystyle (j=1,...,n)}
(3.10)
(
1
⋱
1
−
l
k
+
1
,
k
1
⋮
⋱
−
l
n
,
k
1
)
(
(
z
1
)
T
⋮
(
z
n
)
T
)
=
(
(
z
1
)
T
⋮
(
z
k
)
T
(
z
k
+
1
)
T
−
l
k
+
1
,
k
(
z
k
)
T
⋮
(
z
n
)
T
−
l
n
,
k
(
z
k
)
T
)
{\displaystyle {\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&-l_{k+1,k}&1&&\\&&\vdots &&\ddots &\\&&-l_{n,k}&&&1\end{pmatrix}}{\begin{pmatrix}\left(z^{1}\right)^{T}\\\vdots \\\left(z^{n}\right)^{T}\end{pmatrix}}={\begin{pmatrix}\left(z^{1}\right)^{T}\\\vdots \\\left(z^{k}\right)^{T}\\\left(z^{k+1}\right)^{T}-l_{k+1,k}\left(z^{k}\right)^{T}\\\vdots \\\left(z^{n}\right)^{T}-l_{n,k}\left(z^{k}\right)^{T}\end{pmatrix}}}
Offenbar lässt sich die Frobenius-Matrix in (3.9) mit
(3.11)
f
k
:=
(
0
,
.
.
.
,
0
,
l
k
+
1
,
k
,
.
.
.
,
l
n
,
k
)
T
∈
R
n
{\displaystyle f^{k}:=(0,...,0,l_{k+1,k},...,l_{n,k})^{T}\in \mathbb {R} ^{n}}
wegen
f
k
(
e
k
)
T
=
(
0
⋮
0
0
l
k
+
1
,
k
⋮
l
n
,
k
)
(
0
,
.
.
.
,
0
,
1
⏞
k
-te Stelle
,
0
,
.
.
.
,
0
)
=
(
0
.
.
.
0
0
0
.
.
.
0
⋮
⋱
⋮
⋮
⋮
⋱
⋮
0
.
.
.
0
0
0
.
.
.
0
0
.
.
.
0
0
0
.
.
.
0
0
.
.
.
0
l
k
+
1
,
k
0
.
.
.
0
⋮
⋱
⋮
⋮
⋮
⋱
⋮
0
.
.
.
0
l
n
,
k
0
.
.
.
0
)
{\displaystyle f^{k}(e^{k})^{T}={\begin{pmatrix}0\\\vdots \\0\\0\\l_{k+1,k}\\\vdots \\l_{n,k}\end{pmatrix}}(0,...,0,\overbrace {1} ^{k{\text{-te Stelle}}},0,...,0)={\begin{pmatrix}0&...&0&0&0&...&0\\\vdots &\ddots &\vdots &\vdots &\vdots &\ddots &\vdots \\0&...&0&0&0&...&0\\0&...&0&0&0&...&0\\0&...&0&l_{k+1,k}&0&...&0\\\vdots &\ddots &\vdots &\vdots &\vdots &\ddots &\vdots \\0&...&0&l_{n,k}&0&...&0\end{pmatrix}}}
in der Form
(3.12)
F
k
=
I
−
f
k
(
e
k
)
T
{\displaystyle F_{k}=I-f^{k}(e^{k})^{T}}
darstellen, wobei
I
∈
R
n
×
n
{\displaystyle I\in \mathbb {R} ^{n\times n}}
e
k
∈
R
n
{\displaystyle e^{k}\in \mathbb {R} ^{n}}
k
{\displaystyle k}
I
{\displaystyle I}
Für
k
=
1
,
.
.
.
,
n
−
1
{\displaystyle k=1,...,n-1}
F
k
∈
R
n
×
n
{\displaystyle F_{k}\in \mathbb {R} ^{n\times n}}
k
{\displaystyle k}
k
=
1
,
.
.
.
,
n
−
1
{\displaystyle k=1,...,n-1}
(3.13)
F
k
−
1
=
(
1
⋱
1
l
k
+
1
,
k
1
⋮
⋱
l
n
,
k
1
)
=
I
+
f
k
(
e
k
)
T
{\displaystyle F_{k}^{-1}={\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&l_{k+1,k}&1&&\\&&\vdots &&\ddots &\\&&l_{n,k}&&&1\end{pmatrix}}=I+f^{k}(e^{k})^{T}}
sowie
(3.14)
F
1
−
1
⋯
F
n
−
1
−
1
=
(
1
l
21
1
⋮
l
32
1
⋮
⋮
⋱
⋱
l
n
1
l
n
2
.
.
.
l
n
,
n
−
1
1
)
=
I
+
∑
j
=
1
n
−
1
f
j
(
e
j
)
T
.
{\displaystyle F_{1}^{-1}\cdots F_{n-1}^{-1}={\begin{pmatrix}1&&&&\\l_{21}&1&&&\\\vdots &l_{32}&1&&\\\vdots &\vdots &\ddots &\ddots &\\l_{n1}&l_{n2}&...&l_{n,n-1}&1\end{pmatrix}}=I+\sum _{j=1}^{n-1}f^{j}(e^{j})^{T}.}
Für
F
k
{\displaystyle F_{k}}
[
I
+
f
k
(
e
k
)
T
]
[
I
−
f
k
(
e
k
)
T
]
⏟
=
F
k
=
I
+
f
k
(
e
k
)
T
−
f
k
(
e
k
)
T
−
f
k
[
(
e
k
)
T
f
k
]
⏟
=
0
(
e
k
)
T
=
I
{\displaystyle \left[I+f^{k}(e^{k})^{T}\right]\underbrace {\left[I-f^{k}(e^{k})^{T}\right]} _{=F_{k}}=I+f^{k}(e^{k})^{T}-f^{k}(e^{k})^{T}-f^{k}\underbrace {\left[(e^{k})^{T}f^{k}\right]} _{=0}(e^{k})^{T}=I}
folgt
det
(
I
+
f
k
(
e
k
)
T
)
det
(
I
+
f
k
(
e
k
)
T
)
=
1
,
{\displaystyle \det \left(I+f^{k}(e^{k})^{T}\right)\det \left(I+f^{k}(e^{k})^{T}\right)=1,}
also die Regularität von
F
k
{\displaystyle F_{k}}
F
k
−
1
{\displaystyle F_{k}^{-1}}
(3.15)
F
1
−
1
⋯
F
k
−
1
=
I
+
∑
j
=
1
k
f
j
(
e
j
)
T
,
k
=
1
,
.
.
.
,
n
−
1
{\displaystyle F_{1}^{-1}\cdots F_{k}^{-1}=I+\sum _{j=1}^{k}f^{j}(e^{j})^{T},\quad k=1,...,n-1}
nachgewiesen werden, welche im Fall
k
:=
n
−
1
{\displaystyle k:=n-1}
Die Darstellung in (3.15) ist sicher richtig für
k
:=
1
{\displaystyle k:=1}
k
∈
{
1
,
.
.
.
,
n
−
2
}
{\displaystyle k\in \{1,...,n-2\}}
k
+
1
{\displaystyle k+1}
F
1
−
1
⋯
F
k
−
1
F
k
+
1
−
1
=
[
I
+
∑
j
=
1
k
f
j
(
e
j
)
T
]
[
I
+
f
k
+
1
(
e
k
+
1
)
T
]
{\displaystyle F_{1}^{-1}\cdots F_{k}^{-1}F_{k+1}^{-1}=\left[I+\sum _{j=1}^{k}f^{j}(e^{j})^{T}\right]\left[I+f^{k+1}(e^{k+1})^{T}\right]}
=
I
+
f
k
+
1
(
e
k
+
1
)
T
+
∑
j
=
1
k
f
j
(
e
j
)
T
+
∑
j
=
1
k
f
j
[
(
e
j
)
T
f
k
+
1
]
⏟
=
0
(
e
k
+
1
)
T
=
I
+
∑
j
=
1
k
+
1
f
j
(
e
j
)
T
.
{\displaystyle =I+f^{k+1}(e^{k+1})^{T}+\sum _{j=1}^{k}f^{j}(e^{j})^{T}+\sum _{j=1}^{k}f^{j}\underbrace {\left[(e^{j})^{T}f^{k+1}\right]} _{=0}(e^{k+1})^{T}=I+\sum _{j=1}^{k+1}f^{j}(e^{j})^{T}.}
q.e.d.
Sei
F
k
{\displaystyle F_{k}}
k
{\displaystyle k}
P
{\displaystyle P}
π
{\displaystyle \pi }
i
,
r
∈
{
k
+
1
,
.
.
.
,
n
}
{\displaystyle i,r\in \{k+1,...,n\}}
P
F
k
P
{\displaystyle PF_{k}P}
F
k
{\displaystyle F_{k}}
i
{\displaystyle i}
r
{\displaystyle r}
k
{\displaystyle k}
P
F
k
P
=
I
−
(
P
f
k
)
(
e
k
)
T
.
{\displaystyle PF_{k}P=I-\left(Pf^{k}\right)(e^{k})^{T}.}
Die Aussage ergibt sich unmittelbar aus
P
F
k
P
=
P
[
I
−
f
k
(
e
k
)
T
]
P
=
P
2
⏟
=
I
−
[
P
f
k
]
[
(
e
k
)
T
P
⏟
=
(
e
k
)
T
]
,
{\displaystyle PF_{k}P=P\left[I-f^{k}(e^{k})^{T}\right]P=\underbrace {P^{2}} _{=I}-[Pf^{k}][\underbrace {(e^{k})^{T}P} _{=(e^{k})^{T}}],}
wobei (3.8) und Satz 3.10 (iii) sowie die Forderungen für
i
{\displaystyle i}
r
{\displaystyle r}
q.e.d.
LR -Zerlegung mittels Gauß-Algorithmus
Bearbeiten
Im Folgenden wird die allgemeine Vorgehensweise beim Gauß-Algorithmus zur sukzessiven Erzeugung von Matrizen
A
(
k
)
∈
R
n
×
n
{\displaystyle A^{(k)}\in \mathbb {R} ^{n\times n}}
(3.16)
A
(
k
)
=
(
a
11
(
k
)
a
12
(
k
)
.
.
.
.
.
.
.
.
.
a
1
n
(
k
)
a
22
(
k
)
.
.
.
.
.
.
.
.
.
a
1
n
(
k
)
⋱
⋮
a
k
k
(
k
)
.
.
.
a
k
n
(
k
)
⋮
⋱
⋮
a
n
k
(
k
)
.
.
.
a
n
n
(
k
)
)
∈
R
n
×
n
{\displaystyle A^{(k)}={\begin{pmatrix}a_{11}^{(k)}&a_{12}^{(k)}&...&...&...&a_{1n}^{(k)}\\&a_{22}^{(k)}&...&...&...&a_{1n}^{(k)}\\&&\ddots &&&\vdots \\&&&a_{kk}^{(k)}&...&a_{kn}^{(k)}\\&&&\vdots &\ddots &\vdots \\&&&a_{nk}^{(k)}&...&a_{nn}^{(k)}\end{pmatrix}}\in \mathbb {R} ^{n\times n}}
mit Spaltenpivotsuche als Folge spezieller Matrix-Operationen beschrieben. Und zwar wird im
k
{\displaystyle k}
A
(
k
)
→
P
k
A
(
k
)
{\displaystyle A^{(k)}\to P_{k}A^{(k)}}
P
k
∈
R
n
×
n
{\displaystyle P_{k}\in \mathbb {R} ^{n\times n}}
k
{\displaystyle k}
p
k
≥
k
{\displaystyle p_{k}\geq k}
A
(
k
)
{\displaystyle A^{(k)}}
p
k
=
k
{\displaystyle p_{k}=k}
P
k
=
I
{\displaystyle P_{k}=I}
(3.17)
A
(
k
+
1
)
=
F
k
P
k
A
(
k
)
{\displaystyle A^{(k+1)}=F_{k}P_{k}A^{(k)}}
mit
(3.18)
F
k
:=
(
1
⋱
1
−
l
k
+
1
,
k
1
⋮
⋱
−
l
n
,
k
1
)
{\displaystyle F_{k}:={\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&-l_{k+1,k}&1&&\\&&\vdots &&\ddots &\\&&-l_{n,k}&&&1\end{pmatrix}}}
für
(3.19)
l
p
k
,
k
:=
a
k
k
(
k
)
a
p
k
,
k
(
k
)
,
l
i
k
=
a
i
k
(
k
)
a
p
k
,
k
(
k
)
(
i
=
k
+
1
,
.
.
.
,
n
,
i
≠
p
k
)
{\displaystyle l_{p_{k},k}:={\frac {a_{kk}^{(k)}}{a_{p_{k},k}^{(k)}}},\quad l_{ik}={\frac {a_{ik}^{(k)}}{a_{p_{k},k}^{(k)}}}\quad (i=k+1,...,n,i\neq p_{k})}
sowie
(3.20)
P
k
:=
(
1
⋱
1
0
1
1
⋱
1
1
0
1
⋱
1
)
{\displaystyle P_{k}:={\begin{pmatrix}1&&&&&&&&&&\\&\ddots &&&&&&&&&\\&&1&&&&&&&&\\&&&0&&&&1&&&\\&&&&1&&&&&&\\&&&&&\ddots &&&&&\\&&&&&&1&&&&\\&&&1&&&&0&&&\\&&&&&&&&1&&\\&&&&&&&&&\ddots &\\&&&&&&&&&&1\end{pmatrix}}}
←
Zeile
k
←
Zeile
p
k
{\displaystyle {\begin{array}{ll}\leftarrow &{\text{Zeile }}k\\&\\&\\&\\\leftarrow &{\text{Zeile }}p_{k}\end{array}}}
Der Index
p
k
{\displaystyle p_{k}}
A
(
k
)
{\displaystyle A^{(k)}}
Mit den Definitionen (3.16) - (3.20) gilt die Identität
P
A
=
L
R
{\displaystyle PA=LR}
P
:=
P
n
−
1
⋯
P
1
,
R
:=
A
(
n
)
{\displaystyle P:=P_{n-1}\cdots P_{1},\quad R:=A^{(n)}}
und
(3.21)
L
:=
(
1
l
^
21
1
⋮
l
^
32
1
⋮
⋮
⋱
⋱
l
^
n
1
l
^
n
2
.
.
.
l
^
n
,
n
−
1
1
)
{\displaystyle L:={\begin{pmatrix}1&&&&\\{\hat {l}}_{21}&1&&&\\\vdots &{\hat {l}}_{32}&1&&\\\vdots &\vdots &\ddots &\ddots &\\{\hat {l}}_{n1}&{\hat {l}}_{n2}&...&{\hat {l}}_{n,n-1}&1\end{pmatrix}}}
mit
(3.22)
(
0
⋮
0
1
l
^
k
+
1
,
k
⋮
l
^
n
,
k
)
:=
P
n
−
1
⋯
P
k
+
1
(
0
⋮
0
1
l
k
+
1
,
k
⋮
l
n
,
k
)
{\displaystyle {\begin{pmatrix}0\\\vdots \\0\\1\\{\hat {l}}_{k+1,k}\\\vdots \\{\hat {l}}_{n,k}\end{pmatrix}}:=P_{n-1}\cdots P_{k+1}{\begin{pmatrix}0\\\vdots \\0\\1\\l_{k+1,k}\\\vdots \\l_{n,k}\end{pmatrix}}}
Für
k
=
1
,
2
,
.
.
.
{\displaystyle k=1,2,...}
A
(
2
)
=
F
1
P
1
A
=
F
1
(
P
1
A
)
,
{\displaystyle A^{(2)}=F_{1}P_{1}A=F_{1}(P_{1}A),}
A
(
3
)
=
F
2
P
2
A
(
2
)
=
F
2
(
P
2
F
1
P
2
)
(
P
2
⏟
=
I
P
1
A
)
,
{\displaystyle A^{(3)}=F_{2}P_{2}A^{(2)}=F_{2}(P_{2}F_{1}\underbrace {P_{2})(P_{2}} _{=I}P_{1}A),}
A
(
4
)
=
F
3
P
3
A
(
3
)
=
F
3
(
P
3
F
2
P
3
)
(
P
3
⏟
=
I
P
2
F
1
P
2
P
3
)
(
P
3
⏟
=
I
P
2
P
1
A
)
{\displaystyle A^{(4)}=F_{3}P_{3}A^{(3)}=F_{3}(P_{3}F_{2}\underbrace {P_{3})(P_{3}} _{=I}P_{2}F_{1}P_{2}\underbrace {P_{3})(P_{3}} _{=I}P_{2}P_{1}A)}
und so weiter, was schließlich
(3.23)
R
=
A
(
n
)
=
F
^
n
−
1
⋯
F
^
1
P
A
{\displaystyle R=A^{(n)}={\hat {F}}_{n-1}\cdots {\hat {F}}_{1}PA}
ergibt mit
P
:=
P
n
−
1
⋯
P
1
{\displaystyle P:=P_{n-1}\cdots P_{1}}
und den Frobenius-Matrizen
F
^
n
−
1
=
F
n
−
1
{\displaystyle {\hat {F}}_{n-1}=F_{n-1}}
F
^
k
:=
P
n
−
1
⋯
P
k
+
1
F
k
P
k
+
1
⋯
P
n
−
1
=
(
1
⋱
1
−
l
^
k
+
1
,
k
⋱
⋮
⋱
−
l
^
n
k
1
)
{\displaystyle {\hat {F}}_{k}:=P_{n-1}\cdots P_{k+1}F_{k}P_{k+1}\cdots P_{n-1}={\begin{pmatrix}1&&&&&\\&\ddots &&&&\\&&1&&&\\&&-{\hat {l}}_{k+1,k}&\ddots &&\\&&\vdots &&\ddots &\\&&-{\hat {l}}_{nk}&&&1\end{pmatrix}}}
für
k
=
1
,
.
.
.
,
n
−
2
{\displaystyle k=1,...,n-2}
P
A
=
(
F
^
1
−
1
.
.
.
F
^
n
−
1
−
1
)
R
=
L
R
,
{\displaystyle PA={\begin{pmatrix}{\hat {F}}_{1}^{-1}&...&{\hat {F}}_{n-1}^{-1}\end{pmatrix}}R=LR,}
wobei die letzte Gleichheit mit Lemma 3.13 folgt. Damit ist alles bewiesen.
q.e.d.
Man beachte, dass die Matrix
L
{\displaystyle L}
l
^
i
k
{\displaystyle {\hat {l}}_{ik}}
(
i
>
k
)
{\displaystyle (i>k)}
A
{\displaystyle A}
L
{\displaystyle L}
A
{\displaystyle A}
R
{\displaystyle R}
P
{\displaystyle P}
(
e
r
i
)
T
,
i
=
1
,
.
.
.
,
n
{\displaystyle (e^{r_{i}})^{T},i=1,...,n}
r
∈
N
n
{\displaystyle r\in \mathbb {N} ^{n}}
(
r
1
⋮
r
n
)
=
P
(
1
⋮
n
)
{\displaystyle {\begin{pmatrix}r_{1}\\\vdots \\r_{n}\end{pmatrix}}=P{\begin{pmatrix}1\\\vdots \\n\end{pmatrix}}}
Wir wollen die Vorgehensweise an einem Beispiel vorführen.
Gegeben sei die Matrix
A
:=
(
0
0
1
1
2
2
2
2
1
2
2
2
1
2
3
6
)
.
{\displaystyle A:={\begin{pmatrix}0&0&1&1\\2&2&2&2\\1&2&2&2\\1&2&3&6\end{pmatrix}}.}
Nach Anhängen des für die Speicherung der Zeilenpermutationen zuständigen Buchhaltervektors liefert der Algorithmus folgendes (unterhalb der Treppe ergeben sich sukzessive die Einträge von
L
{\displaystyle L}
(
0
0
1
1
2
2
2
2
1
2
2
2
1
2
3
6
)
,
(
1
2
3
4
)
Zeilentausch
⟶
(
2
2
2
2
0
0
1
1
1
2
2
2
1
2
3
6
)
,
(
2
1
3
4
)
{\displaystyle {\begin{pmatrix}0&0&1&1\\2&2&2&2\\1&2&2&2\\1&2&3&6\end{pmatrix}},{\begin{pmatrix}1\\2\\3\\4\end{pmatrix}}{\underset {\longrightarrow }{\text{ Zeilentausch }}}{\begin{pmatrix}2&2&2&2\\0&0&1&1\\1&2&2&2\\1&2&3&6\end{pmatrix}},{\begin{pmatrix}2\\1\\3\\4\end{pmatrix}}}
Elimination
⟶
(
2
2
2
2
0
0
1
1
1
2
1
_
1
1
1
2
1
2
5
)
,
(
2
1
3
4
)
Zeilentausch
⟶
(
2
2
2
2
1
2
1
1
1
0
0
1
1
1
2
1
2
5
)
,
(
2
3
1
4
)
{\displaystyle {\underset {\longrightarrow }{\text{Elimination }}}{\begin{pmatrix}2&2&2&2\\0&0&1&1\\{\frac {1}{2}}&{\underline {1}}&1&1\\{\frac {1}{2}}&1&2&5\end{pmatrix}},{\begin{pmatrix}2\\1\\3\\4\end{pmatrix}}{\underset {\longrightarrow }{\text{ Zeilentausch }}}{\begin{pmatrix}2&2&2&2\\{\frac {1}{2}}&1&1&1\\0&0&1&1\\{\frac {1}{2}}&1&2&5\end{pmatrix}},{\begin{pmatrix}2\\3\\1\\4\end{pmatrix}}}
Elimination
⟶
(
2
2
2
2
1
2
1
1
1
0
0
1
_
1
1
2
1
1
4
)
,
(
2
3
1
4
)
Elimination
⟶
(
2
2
2
2
1
2
1
1
1
0
0
1
1
1
2
1
1
3
)
,
(
2
3
1
4
)
{\displaystyle {\underset {\longrightarrow }{\text{Elimination }}}{\begin{pmatrix}2&2&2&2\\{\frac {1}{2}}&1&1&1\\0&0&{\underline {1}}&1\\{\frac {1}{2}}&1&1&4\end{pmatrix}},{\begin{pmatrix}2\\3\\1\\4\end{pmatrix}}{\underset {\longrightarrow }{\text{ Elimination }}}{\begin{pmatrix}2&2&2&2\\{\frac {1}{2}}&1&1&1\\0&0&1&1\\{\frac {1}{2}}&1&1&3\end{pmatrix}},{\begin{pmatrix}2\\3\\1\\4\end{pmatrix}}}
Dabei ist das jeweils gewählte Pivotelement unterstrichen. Der letzte Permutationsvektor
(
2
,
3
,
1
,
4
)
T
{\displaystyle (2,3,1,4)^{T}}
π
(
1
)
=
3
,
π
(
2
)
=
1
,
π
(
3
)
=
2
,
π
(
4
)
=
4
{\displaystyle \pi (1)=3,\quad \pi (2)=1,\quad \pi (3)=2,\quad \pi (4)=4}
für die zu
P
{\displaystyle P}
π
{\displaystyle \pi }
P
A
{\displaystyle PA}
A
{\displaystyle A}
A
{\displaystyle A}
P
A
=
(
2
2
2
2
1
2
2
2
0
0
1
1
1
2
3
6
)
=
(
1
1
2
1
0
0
1
1
2
1
1
1
)
(
2
2
2
2
1
1
1
1
1
3
)
=
L
R
.
{\displaystyle PA={\begin{pmatrix}2&2&2&2\\1&2&2&2\\0&0&1&1\\1&2&3&6\end{pmatrix}}={\begin{pmatrix}1&&&\\{\frac {1}{2}}&1&&\\0&0&1&\\{\frac {1}{2}}&1&1&1\end{pmatrix}}{\begin{pmatrix}2&2&2&2\\&1&1&1\\&&1&1\\&&&3\end{pmatrix}}=LR.}
Man beachte, dass die Zerlegung
P
A
=
L
R
{\displaystyle PA=LR}
L
{\displaystyle L}
α
≠
0
{\displaystyle \alpha \neq 0}
R
{\displaystyle R}
1
/
α
{\displaystyle 1/\alpha }
Mit dem Gauß-Algorithmus kann man bekanntlich auch die Determinante von
A
{\displaystyle A}
P
A
=
L
R
{\displaystyle PA=LR}
det
(
P
)
=
(
−
1
)
σ
{\displaystyle \det(P)=(-1)^{\sigma }}
und
det
(
P
A
)
=
det
(
P
)
det
(
A
)
=
det
(
L
)
det
(
R
)
=
∏
i
=
1
n
r
i
i
,
{\displaystyle \det(PA)=\det(P)\det(A)=\det(L)\det(R)=\prod _{i=1}^{n}r_{ii},}
wobei
σ
{\displaystyle \sigma }
I
{\displaystyle I}
P
{\displaystyle P}
r
i
i
{\displaystyle r_{ii}}
R
{\displaystyle R}
det
(
A
)
=
(
−
1
)
σ
∏
i
=
1
n
r
i
i
.
{\displaystyle \det(A)=(-1)^{\sigma }\prod _{i=1}^{n}r_{ii}.}
Aufgrund von Rundungsfehlern errechnet man in der Praxis nicht eine Zerlegung
L
R
{\displaystyle LR}
L
~
R
~
{\displaystyle {\tilde {L}}{\tilde {R}}}
P
A
{\displaystyle PA}
P
A
≈
L
~
R
~
{\displaystyle PA\approx {\tilde {L}}{\tilde {R}}}
gilt. Statt der (hier als eindeutig vorausgesetzten Lösung)
x
{\displaystyle x}
A
x
=
b
{\displaystyle Ax=b}
P
A
x
=
P
b
{\displaystyle PAx=Pb}
L
~
{\displaystyle {\tilde {L}}}
R
~
{\displaystyle {\tilde {R}}}
x
0
{\displaystyle x^{0}}
d
0
:=
P
b
−
P
A
x
0
≠
0.
{\displaystyle d^{0}:=Pb-PAx^{0}\neq 0.}
Daher ist die folgende Nachiteration sinnvoll: wiederum unter Verwendung der vorliegenden Zerlegung
L
~
R
~
{\displaystyle {\tilde {L}}{\tilde {R}}}
P
A
{\displaystyle PA}
δ
x
0
{\displaystyle \delta x^{0}}
(3.24)
P
A
δ
x
=
d
0
{\displaystyle PA\delta x=d^{0}}
und setzt man anschließend
x
1
:=
x
0
+
δ
x
0
.
{\displaystyle x^{1}:=x^{0}+\delta x^{0}.}
Bei exakter Lösung des Systems (3.24) (mit
P
A
{\displaystyle PA}
L
~
R
~
{\displaystyle {\tilde {L}}{\tilde {R}}}
P
A
x
1
=
P
A
x
0
+
P
A
δ
x
0
=
P
b
−
d
0
+
d
0
=
P
b
.
{\displaystyle PAx^{1}=PAx^{0}+PA\delta x^{0}=Pb-d^{0}+d^{0}=Pb.}
Da man i. a. jedoch nicht exakt rechnet, könnte man diesen Prozess wiederholen. Normalerweise genügt es jedoch,
d
0
{\displaystyle d^{0}}
δ
x
0
{\displaystyle \delta x^{0}}
Wir betrachten das Gleichungssystem
A
x
=
b
{\displaystyle Ax=b}
A
:=
(
1.05
1.02
1.04
1.02
)
,
b
:=
(
1
2
)
.
{\displaystyle A:={\begin{pmatrix}1.05&1.02\\1.04&1.02\end{pmatrix}},\quad b:={\begin{pmatrix}1\\2\end{pmatrix}}.}
Die Lösung
x
{\displaystyle x}
x
≈
(
−
100
,
103.921
569
)
T
.
{\displaystyle x\approx (-100,103.921\ 569)^{T}.}
Gauß-Elimination ohne Zeilenvertauschung liefert bei 3-stelliger Rechnung
L
~
:=
(
1
0
.990
1
)
,
R
~
:=
(
1.05
1.02
0
0.01
)
{\displaystyle {\tilde {L}}:={\begin{pmatrix}1&0\\.990&1\end{pmatrix}},\quad {\tilde {R}}:={\begin{pmatrix}1.05&1.02\\0&0.01\end{pmatrix}}}
sowie
L
~
R
~
−
A
=
(
0
0
5
10
−
4
2
10
−
4
)
.
{\displaystyle {\tilde {L}}{\tilde {R}}-A={\begin{pmatrix}0&0\\5_{10}-4&2_{10}-4\end{pmatrix}}.}
Man errechnet
x
0
:=
(
−
97.1
,
101
)
T
{\displaystyle x^{0}:=(-97.1,101)^{T}}
mit dem Defekt
d
0
:=
b
−
A
x
0
=
{
(
0
,
0
)
T
(
3
-stellig
)
,
(
0.065
,
0.035
)
T
(
6
-stellig
)
.
{\displaystyle d^{0}:=b-Ax^{0}={\begin{cases}(0,0)^{T}&(3{\text{-stellig}}),\\(0.065,0.035)^{T}&(6{\text{-stellig}}).\end{cases}}}
Nachiteration mit 6-stelliger Rechnung ergibt
x
1
:=
x
0
+
δ
x
0
=
(
−
99.9
,
104
)
T
.
{\displaystyle x^{1}:=x^{0}+\delta x^{0}=(-99.9,104)^{T}.}
In gewissen Situationen ist es möglich und zwecks Ausnutzung vorhandener Strukturen der Matrix
A
:=
(
a
i
j
)
∈
R
n
×
n
{\displaystyle A:=(a_{ij})\in \mathbb {R} ^{n\times n}}
L
:=
(
l
i
j
)
∈
R
n
×
n
{\displaystyle L:=(l_{ij})\in \mathbb {R} ^{n\times n}}
R
:=
(
r
i
j
)
∈
R
n
×
n
{\displaystyle R:=(r_{ij})\in \mathbb {R} ^{n\times n}}
LR -Zerlegung der Form
(3.25)
(
a
11
a
12
.
.
.
a
1
n
a
21
a
22
.
.
.
a
2
n
⋮
⋮
⋱
⋮
a
n
1
a
n
2
.
.
.
a
n
n
)
=
(
1
l
21
1
⋮
⋱
⋱
l
n
1
.
.
.
l
n
,
n
−
1
1
)
(
r
11
r
12
.
.
.
r
1
n
r
22
.
.
.
r
2
n
⋱
⋮
r
n
n
)
{\displaystyle {\begin{pmatrix}a_{11}&a_{12}&...&a_{1n}\\a_{21}&a_{22}&...&a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&...&a_{nn}\end{pmatrix}}={\begin{pmatrix}1&&&\\l_{21}&1&&\\\vdots &\ddots &\ddots &\\l_{n1}&...&l_{n,n-1}&1\end{pmatrix}}{\begin{pmatrix}r_{11}&r_{12}&...&r_{1n}\\&r_{22}&...&r_{2n}\\&&\ddots &\vdots \\&&&r_{nn}\end{pmatrix}}}
auf direkte Weise zu bestimmen. Eine solche Zerlegung einer regulären Matrix
A
{\displaystyle A}
det
(
H
k
)
≠
0
{\displaystyle \det(H_{k})\neq 0}
(
k
=
1
,
.
.
.
,
n
)
{\displaystyle (k=1,...,n)}
Hauptuntermatrizen
(3.26)
H
k
:=
(
a
11
.
.
.
a
1
k
⋮
⋱
⋮
a
k
1
.
.
.
a
k
k
)
∈
R
k
×
k
,
k
=
1
,
.
.
.
,
n
{\displaystyle H_{k}:={\begin{pmatrix}a_{11}&...&a_{1k}\\\vdots &\ddots &\vdots \\a_{k1}&...&a_{kk}\end{pmatrix}}\in \mathbb {R} ^{k\times k},\quad k=1,...,n}
von
A
{\displaystyle A}
det
(
A
)
≠
0
{\displaystyle \det(A)\neq 0}
det
(
R
)
≠
0
{\displaystyle \det(R)\neq 0}
r
i
i
≠
0
{\displaystyle r_{ii}\neq 0}
(
i
=
1
,
.
.
.
,
n
)
{\displaystyle (i=1,...,n)}
Existiert eine LR -Zerlegung von
A
{\displaystyle A}
n
2
{\displaystyle n^{2}}
r
i
k
{\displaystyle r_{ik}}
(
i
≤
k
)
{\displaystyle (i\leq k)}
l
i
k
{\displaystyle l_{ik}}
(
i
>
k
)
{\displaystyle (i>k)}
n
2
{\displaystyle n^{2}}
a
i
k
=
∑
j
=
1
n
l
i
j
r
j
k
,
i
,
k
=
1
,
.
.
.
,
n
{\displaystyle a_{ik}=\sum _{j=1}^{n}l_{ij}r_{jk},\quad i,k=1,...,n}
zu erhalten, welche wegen
l
i
j
=
0
{\displaystyle l_{ij}=0}
i
<
j
{\displaystyle i<j}
r
j
k
=
0
{\displaystyle r_{jk}=0}
j
>
k
{\displaystyle j>k}
(3.27)
a
i
k
=
∑
j
=
1
min
(
i
,
k
)
l
i
j
r
j
k
,
i
,
k
=
1
,
.
.
.
,
n
{\displaystyle a_{ik}=\sum _{j=1}^{\min(i,k)}l_{ij}r_{jk},\quad i,k=1,...,n}
identisch sind. Dabei gibt es verschiedene Möglichkeiten, wie man aus den Gleichungen (3.27) die Einträge von
L
{\displaystyle L}
R
{\displaystyle R}
R
{\displaystyle R}
L
{\displaystyle L}
(0) Gib
A
:=
(
a
i
j
)
∈
R
n
×
n
{\displaystyle A:=(a_{ij})\in \mathbb {R} ^{n\times n}}
det
(
A
)
≠
0
{\displaystyle \det(A)\neq 0}
i
:=
1
{\displaystyle i:=1}
(1) Berechne
r
i
k
:=
a
i
k
−
∑
j
=
1
i
−
1
l
i
j
r
j
k
(
k
=
i
,
.
.
.
,
n
)
,
{\displaystyle r_{ik}:=a_{ik}-\sum _{j=1}^{i-1}l_{ij}r_{jk}\quad (k=i,...,n),}
l
k
i
:=
(
a
k
i
−
∑
j
=
1
i
−
1
l
k
j
r
j
i
)
/
r
i
i
(
k
=
i
+
1
,
.
.
.
,
n
)
.
{\displaystyle l_{ki}:=\left(a_{ki}-\sum _{j=1}^{i-1}l_{kj}r_{ji}\right)/r_{ii}\quad (k=i+1,...,n).}
(2) Falls
i
=
n
{\displaystyle i=n}
(3) Setze
i
:=
i
+
1
{\displaystyle i:=i+1}
Für eine solche direkte LR -Zerlegung sind insgesamt
(
n
3
/
3
)
(
1
+
0
(
1
/
n
)
)
{\displaystyle (n^{3}/3)(1+0(1/n))}
Cholesky-Zerlegung positiv definiter Matrizen
Bearbeiten
In diesem Abschnitt zeigen wir, dass die eben vorgestellte LR -Zerlegungen für positiv definite Matrizen besonders attraktiv ist.
Eine Matrix
A
∈
R
n
×
n
{\displaystyle A\in \mathbb {R} ^{n\times n}}
positiv definit , falls
A
{\displaystyle A}
A
=
A
T
{\displaystyle A=A^{T}}
⟨
x
,
A
x
⟩
=
x
T
A
x
>
0
,
x
∈
R
n
∖
{
0
}
.
{\displaystyle \langle x,Ax\rangle =x^{T}Ax>0,\quad x\in \mathbb {R} ^{n}\setminus \{0\}.}
Mit dieser Definition selbst lässt sich, wie es ähnlich häufig in der Mathematik der Fall ist, nur schwer arbeiten. Deshalb sind äquivalente Bedingungen von Interesse.
Hesse-Matrix und ein lokales Minimum einer Fehlerfunktion
Bearbeiten
Sei
f
:
U
⊂
R
n
→
R
{\displaystyle f\colon U\subset \mathbb {R} ^{n}\to \mathbb {R} }
zweimal stetig differenzierbare Funktion . Dann ist die Hesse-Matrix von
f
{\displaystyle f}
x
=
(
x
1
,
…
,
x
n
)
∈
U
{\displaystyle x=(x_{1},\ldots ,x_{n})\in U}
H
f
(
x
)
:=
(
∂
2
f
∂
x
i
∂
x
j
(
x
)
)
i
,
j
=
1
,
…
,
n
=
(
∂
2
f
∂
x
1
∂
x
1
(
x
)
∂
2
f
∂
x
1
∂
x
2
(
x
)
⋯
∂
2
f
∂
x
1
∂
x
n
(
x
)
∂
2
f
∂
x
2
∂
x
1
(
x
)
∂
2
f
∂
x
2
∂
x
2
(
x
)
⋯
∂
2
f
∂
x
2
∂
x
n
(
x
)
⋮
⋮
⋱
⋮
∂
2
f
∂
x
n
∂
x
1
(
x
)
∂
2
f
∂
x
n
∂
x
2
(
x
)
⋯
∂
2
f
∂
x
n
∂
x
n
(
x
)
)
.
{\displaystyle \operatorname {H} _{f}(x):=\left({\frac {\partial ^{2}f}{\partial x_{i}\partial x_{j}}}(x)\right)_{i,j=1,\dots ,n}={\begin{pmatrix}{\frac {\partial ^{2}f}{\partial x_{1}\partial x_{1}}}(x)&{\frac {\partial ^{2}f}{\partial x_{1}\partial x_{2}}}(x)&\cdots &{\frac {\partial ^{2}f}{\partial x_{1}\partial x_{n}}}(x)\\[0.5em]{\frac {\partial ^{2}f}{\partial x_{2}\partial x_{1}}}(x)&{\frac {\partial ^{2}f}{\partial x_{2}\partial x_{2}}}(x)&\cdots &{\frac {\partial ^{2}f}{\partial x_{2}\partial x_{n}}}(x)\\\vdots &\vdots &\ddots &\vdots \\{\frac {\partial ^{2}f}{\partial x_{n}\partial x_{1}}}(x)&{\frac {\partial ^{2}f}{\partial x_{n}\partial x_{2}}}(x)&\cdots &{\frac {\partial ^{2}f}{\partial x_{n}\partial x_{n}}}(x)\end{pmatrix}}.}
Ist die Hesse-Matrix an einer Stelle
x
=
(
x
1
,
…
,
x
n
)
∈
U
{\displaystyle x=(x_{1},\ldots ,x_{n})\in U}
f
{\displaystyle f}
Bearbeiten
Folgende Aussagen sind für eine Matrix
A
:=
(
a
i
j
)
∈
R
n
×
n
{\displaystyle A:=(a_{ij})\in \mathbb {R} ^{n\times n}}
A
=
A
T
{\displaystyle A=A^{T}}
(i)
A
{\displaystyle A}
(ii) Alle Eigenwerte
λ
i
{\displaystyle \lambda _{i}}
(
i
=
1
,
.
.
.
,
n
)
{\displaystyle (i=1,...,n)}
A
{\displaystyle A}
(iii) Die Determinanten der
n
{\displaystyle n}
H
k
{\displaystyle H_{k}}
(
k
=
1
,
.
.
.
,
n
)
{\displaystyle (k=1,...,n)}
A
{\displaystyle A}
Bearbeiten
Den Beweis des Lemmas findet man in Büchern der Linearen Algebra. Im Folgenden beweisen wir Eigenschaften positiv definite Matrizen
A
{\displaystyle A}
Lemma - positiv Definitheit von Untermatrizen
Bearbeiten
Die Matrix
A
:=
(
a
i
j
)
∈
R
n
×
n
{\displaystyle A:=(a_{ij})\in \mathbb {R} ^{n\times n}}
(i)
A
{\displaystyle A}
(ii) alle Untermatrizen von
A
{\displaystyle A}
(3.29)
(
a
k
k
.
.
.
a
k
,
k
+
ℓ
⋮
⋱
⋮
a
k
+
ℓ
,
k
.
.
.
a
k
+
ℓ
,
k
+
ℓ
)
∈
R
(
ℓ
+
1
)
×
(
ℓ
+
1
)
{\displaystyle {\begin{pmatrix}a_{kk}&...&a_{k,k+\ell }\\\vdots &\ddots &\vdots \\a_{k+\ell ,k}&...&a_{k+\ell ,k+\ell }\end{pmatrix}}\in \mathbb {R} ^{(\ell +1)\times (\ell +1)}}
sind positiv definit, (iii)
det
(
A
)
>
0
{\displaystyle \det(A)>0}
Der Beweis erfolgt in der Reihenfolge der drei Eigenschaften von positiv definiten Matrizen aus dem Lemma.
(i) Wäre
A
{\displaystyle A}
x
∈
R
n
∖
{
0
}
{\displaystyle x\in \mathbb {R} ^{n}\setminus \{0\}}
A
x
=
0
{\displaystyle Ax=0}
x
T
A
x
=
0
{\displaystyle x^{T}Ax=0}
A
{\displaystyle A}
(ii) Sei nun
B
∈
R
(
ℓ
+
1
)
×
(
ℓ
+
1
)
{\displaystyle B\in \mathbb {R} ^{(\ell +1)\times (\ell +1)}}
x
:=
(
x
i
)
∈
R
ℓ
+
1
∖
{
0
}
{\displaystyle x:=(x_{i})\in \mathbb {R} ^{\ell +1}\setminus \{0\}}
A
{\displaystyle A}
B
{\displaystyle B}
z
:=
(
z
i
)
∈
R
n
{\displaystyle z:=(z_{i})\in \mathbb {R} ^{n}}
z
i
:=
{
x
i
+
1
−
k
,
k
≤
i
≤
k
+
ℓ
0
,
sonst
{\displaystyle z_{i}:={\begin{cases}x_{i+1-k},&k\leq i\leq k+\ell \\0,&{\text{sonst}}\end{cases}}}
gilt dann
z
≠
0
{\displaystyle z\neq 0}
x
T
B
x
=
∑
i
,
j
=
k
k
+
ℓ
a
i
j
x
i
+
1
−
k
x
j
+
1
−
k
=
∑
i
,
j
=
k
k
+
ℓ
a
i
j
z
i
z
j
=
∑
i
,
j
=
1
n
a
i
j
z
i
z
j
=
z
T
A
z
>
0.
{\displaystyle x^{T}Bx=\sum _{i,j=k}^{k+\ell }a_{ij}x_{i+1-k}x_{j+1-k}=\sum _{i,j=k}^{k+\ell }a_{ij}z_{i}z_{j}=\sum _{i,j=1}^{n}a_{ij}z_{i}z_{j}=z^{T}Az>0.}
(iii) Die Eigenwerte
λ
1
,
.
.
.
,
λ
n
{\displaystyle \lambda _{1},...,\lambda _{n}}
A
{\displaystyle A}
x
(
k
)
∈
R
n
∖
{
0
}
{\displaystyle x^{(k)}\in \mathbb {R} ^{n}\setminus \{0\}}
A
x
(
k
)
=
λ
k
x
(
k
)
{\displaystyle Ax^{(k)}=\lambda _{k}x^{(k)}}
0
<
(
x
(
k
)
)
T
A
x
(
k
)
=
λ
k
[
(
x
(
k
)
)
T
x
(
k
)
]
,
k
=
1
,
.
.
.
,
n
{\displaystyle 0<(x^{(k)})^{T}Ax^{(k)}=\lambda _{k}\left[(x^{(k)})^{T}x{(k)}\right],\quad k=1,...,n}
(siehe auch Lemma 3.19). Weiter ist die Matrix
A
{\displaystyle A}
S
∈
R
n
×
n
{\displaystyle S\in \mathbb {R} ^{n\times n}}
S
A
S
−
1
=
D
,
{\displaystyle SAS^{-1}=D,}
wobei
D
{\displaystyle D}
D
:=
diag
(
λ
1
,
.
.
.
,
λ
n
)
∈
R
n
×
n
{\displaystyle D:=\operatorname {diag} (\lambda _{1},...,\lambda _{n})\in \mathbb {R} ^{n\times n}}
ist. Somit folgt
det
(
A
)
=
det
(
S
−
1
)
det
(
D
)
det
(
S
)
=
det
(
D
)
=
∏
k
=
1
n
λ
k
>
0.
{\displaystyle \det(A)=\det(S^{-1})\det(D)\det(S)=\det(D)=\prod _{k=1}^{n}\lambda _{k}>0.}
q.e.d.
Satz - Produktzerlegung mit unterer Dreiecksmatrix
Bearbeiten
Die Matrix
A
∈
R
n
×
n
{\displaystyle A\in \mathbb {R} ^{n\times n}}
L
:=
(
l
k
j
)
∈
R
n
×
n
{\displaystyle L:=(l_{kj})\in \mathbb {R} ^{n\times n}}
l
k
k
>
0
{\displaystyle l_{kk}>0}
(
k
=
1
,
.
.
.
,
n
)
{\displaystyle (k=1,...,n)}
(3.30)
A
=
L
L
T
.
{\displaystyle A=LL^{T}.}
Der Beweis wird mit vollständiger Induktion geführt. Für
n
:=
1
{\displaystyle n:=1}
A
=
(
α
)
∈
R
1
×
1
{\displaystyle A=(\alpha )\in \mathbb {R} ^{1\times 1}}
α
>
0
{\displaystyle \alpha >0}
α
=
l
⋅
l
,
l
=
α
{\displaystyle \alpha =l\cdot l,\quad l={\sqrt {\alpha }}}
geschrieben werden. Wir nehmen nun an, dass die Behauptung für positiv definite Matrizen bis zur Dimension
n
−
1
{\displaystyle n-1}
A
∈
R
n
×
n
{\displaystyle A\in \mathbb {R} ^{n\times n}}
A
n
−
1
∈
R
(
n
−
1
)
×
(
n
−
1
)
{\displaystyle A_{n-1}\in \mathbb {R} ^{(n-1)\times (n-1)}}
s
∈
R
n
{\displaystyle s\in \mathbb {R} ^{n}}
A
=
(
A
n
−
1
s
s
T
a
n
n
)
{\displaystyle A={\begin{pmatrix}A_{n-1}&s\\s^{T}&a_{nn}\end{pmatrix}}}
partitionieren, wobei
A
n
−
1
{\displaystyle A_{n-1}}
L
n
−
1
:=
(
l
k
j
)
∈
R
(
n
−
1
)
×
(
n
−
1
)
{\displaystyle L_{n-1}:=(l_{kj})\in \mathbb {R} ^{(n-1)\times (n-1)}}
l
k
k
>
0
{\displaystyle l_{kk}>0}
(
k
=
1
,
.
.
.
,
n
−
1
)
{\displaystyle (k=1,...,n-1)}
(3.31)
A
n
−
1
=
L
n
−
1
L
n
−
1
T
.
{\displaystyle A_{n-1}=L_{n-1}L_{n-1}^{T}.}
Für die gesuchte Matrix
L
∈
R
n
×
n
{\displaystyle L\in \mathbb {R} ^{n\times n}}
L
:=
(
L
n
−
1
0
c
T
α
)
{\displaystyle L:={\begin{pmatrix}L_{n-1}&0\\c^{T}&\alpha \end{pmatrix}}}
und versuchen wir
c
∈
R
n
−
1
{\displaystyle c\in \mathbb {R} ^{n-1}}
α
∈
R
{\displaystyle \alpha \in \mathbb {R} }
(3.32)
A
=
(
A
n
−
1
s
s
T
a
n
n
)
=
(
L
n
−
1
0
c
T
α
)
(
L
n
−
1
T
c
0
α
)
{\displaystyle A={\begin{pmatrix}A_{n-1}&s\\s^{T}a_{nn}\end{pmatrix}}={\begin{pmatrix}L_{n-1}&0\\c^{T}&\alpha \end{pmatrix}}{\begin{pmatrix}L_{n-1}^{T}&c\\0&\alpha \end{pmatrix}}}
gilt. Gleichheit in (3.32) hat man nun wegen (3.31) genau dann, wenn
(3.33)
L
n
−
1
c
=
s
,
{\displaystyle L_{n-1}c=s,}
(3.34)
c
T
c
+
α
2
=
a
n
n
{\displaystyle c^{T}c+\alpha ^{2}=a_{nn}}
gilt. Die Gleichung (3.33) besitzt sicher eine eindeutige Lösung
c
=
L
n
−
1
−
1
s
{\displaystyle c=L_{n-1}^{-1}s}
L
n
−
1
∈
R
(
n
−
1
)
×
(
n
−
1
)
{\displaystyle L_{n-1}\in \mathbb {R} ^{(n-1)\times (n-1)}}
α
∈
C
{\displaystyle \alpha \in \mathbb {C} }
det
(
A
)
=
det
(
L
n
−
1
0
c
T
α
)
det
(
L
n
−
1
T
c
0
α
)
=
α
2
[
det
(
L
n
−
1
)
]
2
,
{\displaystyle \det(A)=\det {\begin{pmatrix}L_{n-1}&0\\c^{T}&\alpha \end{pmatrix}}\det {\begin{pmatrix}L_{n-1}^{T}&c\\0&\alpha \end{pmatrix}}=\alpha ^{2}[\det(L_{n-1})]^{2},}
so dass wegen
det
(
A
)
>
0
{\displaystyle \det(A)>0}
det
(
L
n
−
1
)
=
∏
k
=
1
n
l
k
k
>
0
{\displaystyle \det(L_{n-1})=\prod _{k=1}^{n}l_{kk}>0}
α
2
>
0
{\displaystyle \alpha ^{2}>0}
α
>
0
{\displaystyle \alpha >0}
q.e.d.
Die Zerlegung
A
=
L
L
T
{\displaystyle A=LL^{T}}
A
:=
(
a
i
j
)
∈
R
n
×
n
{\displaystyle A:=(a_{ij})\in \mathbb {R} ^{n\times n}}
Cholesky-Zerlegung von
A
{\displaystyle A}
A
=
L
L
T
{\displaystyle A=LL^{T}}
(
a
11
a
12
.
.
.
a
1
n
a
21
a
22
.
.
.
a
2
n
⋮
⋮
⋱
⋮
a
n
1
a
n
2
.
.
.
a
n
n
)
=
(
l
11
l
21
l
22
⋮
⋱
⋱
l
n
1
.
.
.
l
n
,
n
−
1
l
n
n
)
(
l
11
l
21
.
.
.
l
n
1
l
22
.
.
.
l
n
2
⋱
⋮
l
n
n
)
{\displaystyle {\begin{pmatrix}a_{11}&a_{12}&...&a_{1n}\\a_{21}&a_{22}&...&a_{2n}\\\vdots &\vdots &\ddots &\vdots \\a_{n1}&a_{n2}&...&a_{nn}\end{pmatrix}}={\begin{pmatrix}l_{11}&&&\\l_{21}&l_{22}&&\\\vdots &\ddots &\ddots &\\l_{n1}&...&l_{n,n-1}&l_{nn}\end{pmatrix}}{\begin{pmatrix}l_{11}&l_{21}&...&l_{n1}\\&l_{22}&...&l_{n2}\\&&\ddots &\vdots \\&&&l_{nn}\end{pmatrix}}}
als
n
(
n
+
1
)
/
2
{\displaystyle n(n+1)/2}
n
(
n
+
1
)
/
2
{\displaystyle n(n+1)/2}
l
i
k
{\displaystyle l_{ik}}
(
i
≥
k
)
{\displaystyle (i\geq k)}
a
i
k
=
∑
j
=
1
k
l
i
j
l
k
j
,
1
≤
k
≤
i
≤
n
(
i
=
1
,
.
.
.
,
n
)
.
{\displaystyle a_{ik}=\sum _{j=1}^{k}l_{ij}l_{kj},\quad 1\leq k\leq i\leq n\quad (i=1,...,n).}
Spaltenweise Berechnung der Einträge der unteren Dreiecksmatrix
L
∈
R
n
×
n
{\displaystyle L\in \mathbb {R} ^{n\times n}}
(0) Gib eine positiv definite Matrix
A
:=
(
a
i
j
)
∈
R
n
×
n
{\displaystyle A:=(a_{ij})\in \mathbb {R} ^{n\times n}}
k
:=
1
{\displaystyle k:=1}
(1) Berechne
l
k
k
:=
(
a
k
k
−
∑
j
=
1
k
−
1
l
k
j
2
)
1
/
2
,
{\displaystyle l_{kk}:=\left(a_{kk}-\sum _{j=1}^{k-1}l_{kj}^{2}\right)^{1/2},}
l
i
k
:=
(
a
i
k
−
∑
j
=
1
k
−
1
l
i
j
l
k
j
)
/
l
k
k
(
i
=
k
+
1
,
.
.
.
,
n
)
.
{\displaystyle l_{ik}:=\left(a_{ik}-\sum _{j=1}^{k-1}l_{ij}l_{kj}\right)/l_{kk}\quad (i=k+1,...,n).}
(2) Falls
k
=
n
{\displaystyle k=n}
(3) Setze
k
:=
k
+
1
{\displaystyle k:=k+1}
Gegeben sei die positiv definite Matrix
A
:=
(
a
i
j
)
=
(
1
2
1
2
5
2
1
2
10
)
.
{\displaystyle A:=(a_{ij})={\begin{pmatrix}1&2&1\\2&5&2\\1&2&10\end{pmatrix}}.}
Dann errechnet man für den Spaltenindex
k
:=
1
{\displaystyle k:=1}
l
11
:=
a
11
=
1
=
1
,
{\displaystyle l_{11}:={\sqrt {a_{11}}}={\sqrt {1}}=1,}
l
21
:=
a
21
/
l
11
=
2
/
1
=
2
,
{\displaystyle l_{21}:=a_{21}/l_{11}=2/1=2,}
l
31
:=
a
31
/
l
11
=
1
/
1
=
1
,
{\displaystyle l_{31}:=a_{31}/l_{11}=1/1=1,}
für den Spaltenindex
k
:=
2
{\displaystyle k:=2}
l
22
:=
a
22
−
l
21
2
=
5
−
4
=
1
,
{\displaystyle l_{22}:={\sqrt {a_{22}-l_{21}^{2}}}={\sqrt {5-4}}=1,}
l
32
:=
(
a
32
−
l
31
l
21
)
/
l
22
=
(
2
−
1
⋅
2
)
/
1
=
0
{\displaystyle l_{32}:=(a_{32}-l_{31}l_{21})/l_{22}=(2-1\cdot 2)/1=0}
und schließlich für den Spaltenindex
k
:=
3
{\displaystyle k:=3}
l
33
:=
a
33
−
l
31
2
−
l
32
2
=
10
−
1
2
−
0
2
=
3.
{\displaystyle l_{33}:={\sqrt {a_{33}-l_{31}^{2}-l_{32}^{2}}}={\sqrt {10-1^{2}-0^{2}}}=3.}
Somit erhält man für
A
{\displaystyle A}
A
=
(
1
0
0
2
1
0
1
0
3
)
(
1
2
1
0
1
0
0
0
3
)
.
{\displaystyle A={\begin{pmatrix}1&0&0\\2&1&0\\1&0&3\end{pmatrix}}{\begin{pmatrix}1&2&1\\0&1&0\\0&0&3\end{pmatrix}}.}
Eine Cholesky-Zerlegung erfordert insgesamt die folgende Anzahl von Multiplikationen, Divisionen und Berechnungen von Quadratwurzeln:
∑
k
=
1
n
[
k
+
∑
i
=
k
+
1
n
k
]
=
∑
k
=
1
n
[
k
+
(
n
−
k
)
k
]
=
∑
k
=
1
n
k
+
n
∑
k
=
1
n
k
−
∑
k
=
1
n
k
2
{\displaystyle \sum _{k=1}^{n}\left[k+\sum _{i=k+1}^{n}k\right]=\sum _{k=1}^{n}[k+(n-k)k]=\sum _{k=1}^{n}k+n\sum _{k=1}^{n}k-\sum _{k=1}^{n}k^{2}}
=
n
(
n
+
1
)
2
+
n
n
(
n
+
1
)
2
−
n
(
n
+
1
)
(
2
n
+
1
)
6
=
n
3
6
(
1
+
0
(
1
n
)
)
.
{\displaystyle ={\frac {n(n+1)}{2}}+n{\frac {n(n+1)}{2}}-{\frac {n(n+1)(2n+1)}{6}}={\frac {n^{3}}{6}}\left(1+0\left({\frac {1}{n}}\right)\right).}
Dies sind etwa halb so viele wesentliche Rechenoperationen, wie sie der Gauß-Algorithmus bzw. eine direkte LR -Zerlegung für eine beliebige reguläre Matrix erfordern.
Bei der Diskretisierung von gewöhnlichen und partiellen Differentialgleichungen oder auch der Berechnung der Momente kubischer Splines (vgl. Abschnitt 7) ergeben sich lineare Gleichungssysteme
A
x
=
b
{\displaystyle Ax=b}
A
=
(
a
i
j
)
∈
R
n
×
n
{\displaystyle A=(a_{ij})\in \mathbb {R} ^{n\times n}}
det
(
A
)
≠
0
{\displaystyle \det(A)\neq 0}
Bandmatrix der Bandbreite
p
+
q
+
1
{\displaystyle p+q+1}
A
{\displaystyle A}
(3.35)
A
:=
(
a
11
.
.
.
a
1
,
q
+
1
⋮
⋱
⋱
a
p
+
1
,
1
⋱
⋱
⋱
⋱
a
n
−
q
,
n
⋱
⋱
⋮
a
n
,
n
−
p
.
.
.
a
n
n
)
{\displaystyle A:={\begin{pmatrix}a_{11}&...&a_{1,q+1}&&&\\\vdots &\ddots &&\ddots &&\\a_{p+1,1}&&\ddots &&\ddots &\\&\ddots &&\ddots &&a_{n-q,n}\\&&\ddots &&\ddots &\vdots \\&&&a_{n,n-p}&...&a_{nn}\end{pmatrix}}}
hat und somit
a
i
k
:=
0
,
i
=
1
,
.
.
.
,
n
,
1
≤
k
<
i
−
p
≤
n
,
1
≤
i
+
q
<
k
≤
n
{\displaystyle a_{ik}:=0,\quad i=1,...,n,\quad 1\leq k<i-p\leq n,\quad 1\leq i+q<k\leq n}
mit gewissen
p
,
q
∈
{
1
,
.
.
.
,
n
−
1
}
{\displaystyle p,q\in \{1,...,n-1\}}
p
=
q
=
1
{\displaystyle p=q=1}
A
:=
(
a
11
a
12
a
21
a
22
a
23
a
32
a
33
⋱
⋱
⋱
a
n
−
2
,
n
−
1
a
n
−
1
,
n
−
2
a
n
−
1
,
n
−
1
a
n
−
1
,
n
a
n
,
n
−
1
a
n
n
)
{\displaystyle A:={\begin{pmatrix}a_{11}&a_{12}&&&&\\a_{21}&a_{22}&a_{23}&&&\\&a_{32}&a_{33}&\ddots &&\\&&\ddots &\ddots &a_{n-2,n-1}&\\&&&a_{n-1,n-2}&a_{n-1,n-1}&a_{n-1,n}\\&&&&a_{n,n-1}&a_{nn}\end{pmatrix}}}
von einer Tridiagonalmatrix .
Bei Gleichungssystemen mit Bandmatrizen lässt sich der zu betreibende Aufwand bei allen in diesem Kapitel angesprochenen Methoden verringern, außer bei denen mit Pivotstrategien, da diese die Bandstruktur im Allgemeinen zerstören. Exemplarisch soll das Vorgehen für Bandmatrizen am Beispiel der direkten LR -Zerlegung demonstriert werden. Wenn eine LR -Zerlegung für
A
{\displaystyle A}
det
(
A
)
≠
0
{\displaystyle \det(A)\neq 0}
L
{\displaystyle L}
A
=
(
1
l
21
⋱
⋮
⋱
⋱
l
p
+
1
,
1
⋱
⋱
⋱
⋱
⋱
l
n
,
n
−
p
.
.
.
l
n
,
n
−
1
1
)
(
r
11
.
.
.
r
1
,
q
+
1
⋱
⋱
⋱
r
n
−
q
,
n
⋱
⋮
r
n
n
)
{\displaystyle A={\begin{pmatrix}1&&&&&\\l_{21}&\ddots &&&&\\\vdots &\ddots &\ddots &&&\\l_{p+1,1}&&\ddots &\ddots &&\\&\ddots &&\ddots &\ddots &\\&&l_{n,n-p}&...&l_{n,n-1}&1\end{pmatrix}}{\begin{pmatrix}r_{11}&...&r_{1,q+1}&&\\&\ddots &&\ddots &\\&&\ddots &&r_{n-q,n}\\&&&\ddots &\vdots \\&&&&r_{nn}\end{pmatrix}}}
(siehe z. B. Satz 2.29 bei Kanzow). Komponentenschreibweise geschrieben heißt dies
a
i
k
=
∑
j
=
j
0
min
{
i
,
k
}
l
i
j
r
j
k
,
i
=
1
,
.
.
.
,
n
,
k
=
max
{
1
,
i
−
p
}
,
.
.
.
,
min
{
n
,
i
+
q
}
,
j
0
:=
max
{
1
,
i
−
p
,
k
−
q
}
,
{\displaystyle a_{ik}=\sum _{j=j_{0}}^{\min\{i,k\}}l_{ij}r_{jk},\quad i=1,...,n,\quad k=\max\{1,i-p\},...,\min\{n,i+q\},\quad j_{0}:=\max\{1,i-p,k-q\},}
was bei einer Parkettierung wie in (3.28) auf den folgenden Algorithmus zur Bestimmung der LR -Zerlegung der Bandmatrix
A
{\displaystyle A}
LR -Zerlegung für Bandmatrizen)
Bearbeiten
(0) Gib eine Matrix
A
:=
(
a
i
j
)
∈
R
n
×
n
{\displaystyle A:=(a_{ij})\in \mathbb {R} ^{n\times n}}
a
i
k
:=
0
,
i
=
1
,
.
.
.
,
n
,
1
≤
k
<
i
−
p
≤
n
,
1
≤
i
+
q
<
k
≤
n
{\displaystyle a_{ik}:=0,\quad i=1,...,n,\quad 1\leq k<i-p\leq n,\quad 1\leq i+q<k\leq n}
für gegebene
p
,
q
∈
{
1
,
.
.
.
,
n
−
1
}
{\displaystyle p,q\in \{1,...,n-1\}}
i
:=
1
{\displaystyle i:=1}
(1) Für
k
=
i
,
.
.
.
,
min
{
i
+
q
,
n
}
{\displaystyle k=i,...,\min\{i+q,n\}}
j
0
=
max
{
1
,
i
−
p
,
k
−
q
}
{\displaystyle j_{0}=\max\{1,i-p,k-q\}}
r
i
k
:=
a
i
k
−
∑
j
=
j
0
i
−
1
l
i
j
r
j
k
.
{\displaystyle r_{ik}:=a_{ik}-\sum _{j=j_{0}}^{i-1}l_{ij}r_{jk}.}
(2) Für
k
=
i
+
1
,
.
.
.
,
min
{
i
+
p
,
n
}
{\displaystyle k=i+1,...,\min\{i+p,n\}}
j
0
=
max
{
1
,
i
−
p
,
k
−
q
}
{\displaystyle j_{0}=\max\{1,i-p,k-q\}}
l
k
i
:=
(
a
k
i
−
∑
j
=
j
0
i
−
1
l
k
j
r
j
i
)
/
r
i
i
.
{\displaystyle l_{ki}:=\left(a_{ki}-\sum _{j=j_{0}}^{i-1}l_{kj}r_{ji}\right)/r_{ii}.}
(3) Falls
i
=
n
{\displaystyle i=n}
(4) Setze
i
:=
i
+
1
{\displaystyle i:=i+1}
Ist
A
{\displaystyle A}
(3.36)
(
α
1
β
2
γ
2
α
2
⋱
⋱
⋱
β
n
γ
n
α
n
)
=
(
1
l
2
1
⋱
⋱
l
n
1
)
(
d
1
r
2
⋱
⋱
⋱
r
n
d
n
)
,
{\displaystyle {\begin{pmatrix}\alpha _{1}&\beta _{2}&&\\\gamma _{2}&\alpha _{2}&\ddots &\\&\ddots &\ddots &\beta _{n}\\&&\gamma _{n}&\alpha _{n}\end{pmatrix}}={\begin{pmatrix}1&&&\\l_{2}&1&&\\&\ddots &\ddots &\\&&l_{n}&1\end{pmatrix}}{\begin{pmatrix}d_{1}&r_{2}&&\\&\ddots &\ddots &\\&&\ddots &r_{n}\\&&&d_{n}\end{pmatrix}},}
so vereinfacht sich Algorithmus 4 zu
LR -Zerlegung Tridiagonalmatrizen)
Bearbeiten
(0) Gib eine Tridiagonalmatrix
A
∈
R
n
×
n
{\displaystyle A\in \mathbb {R} ^{n\times n}}
A
,
L
{\displaystyle A,L}
R
{\displaystyle R}
d
1
:=
α
1
,
r
2
:=
β
2
{\displaystyle d_{1}:=\alpha _{1},\quad r_{2}:=\beta _{2}}
und
i
:=
2
{\displaystyle i:=2}
(1) Berechne
l
i
:=
γ
i
/
d
i
−
1
,
d
i
:=
α
i
−
l
i
r
i
,
r
i
+
1
:=
β
i
+
1
.
{\displaystyle l_{i}:=\gamma _{i}/d_{i-1},\quad d_{i}:=\alpha _{i}-l_{i}r_{i},\quad r_{i+1}:=\beta _{i+1}.}
(2) Falls
i
=
n
−
1
{\displaystyle i=n-1}
l
n
:=
γ
n
/
d
n
−
1
,
d
n
:=
α
n
−
l
n
r
n
{\displaystyle l_{n}:=\gamma _{n}/d_{n-1},\quad d_{n}:=\alpha _{n}-l_{n}r_{n}}
und stoppe!
(3) Setze
i
:=
i
+
1
{\displaystyle i:=i+1}
Man kann zeigen, dass im Fall einer Tridiagonalmatrix eine LR -Zerlegung wie in (3.36) möglich ist, wenn gilt (Lemma 2.28 bei Kanzow):
|
α
1
|
>
|
β
2
|
,
|
α
i
|
≥
|
γ
i
|
+
|
β
i
+
1
|
,
i
=
2
,
.
.
.
,
n
−
1
,
{\displaystyle |\alpha _{1}|>|\beta _{2}|,\quad |\alpha _{i}|\geq |\gamma _{i}|+|\beta _{i+1}|,\quad i=2,...,n-1,}
|
α
n
|
≥
|
γ
n
|
,
γ
i
≠
0
,
i
=
2
,
.
.
.
,
n
.
{\displaystyle |\alpha _{n}|\geq |\gamma _{n}|,\quad \gamma _{i}\neq 0,\quad i=2,...,n.}
Diese Bedingungen besagen offenbar, dass man für die erste Zeile strikte und für die anderen Zeilen nur normale Diagonaldominanz fordern muss. Die Forderung
γ
i
≠
0
{\displaystyle \gamma _{i}\neq 0}
(
i
=
2
,
.
.
.
,
n
)
{\displaystyle (i=2,...,n)}
LR -Zerlegung mit
L
:=
I
{\displaystyle L:=I}
R
:=
A
{\displaystyle R:=A}
LR -Zerlegung einer Tridiagonalmatrix sind offenbar nur
(
n
−
2
)
⋅
2
+
2
=
2
n
−
2
{\displaystyle (n-2)\cdot 2+2=2n-2}
wesentliche Rechenoperationen erforderlich.
![{\displaystyle \left(P^{T}P\right)_{ij}=\left[{\begin{pmatrix}\left(e^{\pi (1)}\right)^{T}\\\vdots \\\left(e^{\pi (n)}\right)^{T}\end{pmatrix}}{\begin{pmatrix}e^{\pi (1)}&...&e^{\pi (n)}\end{pmatrix}}\right]_{ij}=\left[\left(e^{\pi (i)}\right)^{T}e^{\pi (j)}\right]_{ij}=(\delta _{\pi (i),\pi (j)})_{ij}=(\delta _{ij})_{ij}=I.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2f8414a5c451049975f3fa469fe0e7c28a3a4a71)
![{\displaystyle \left[{\begin{pmatrix}e^{\pi ^{-1}(1)}&...&e^{\pi ^{-1}(n)}\end{pmatrix}}{\begin{pmatrix}e^{\pi (1)}&...&e^{\pi (n)}\end{pmatrix}}\right]_{ij}=\sum _{\ell =1}^{n}\left(e^{\pi ^{-1}(\ell )}\right)_{i}\left(e^{\pi (j)}\right)_{\ell }=\sum _{\ell =1}^{n}\delta _{i,\pi ^{-1}(\ell )}\delta _{\pi (j),\ell }}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6a82bb20bf68e30395a5f608f97eb7ae49de5cae)
![{\displaystyle AP=\left(P^{T}A^{T}\right)^{T}=\left[P^{T}{\begin{pmatrix}(a^{1})^{T}\\\vdots \\(a^{n})^{T}\end{pmatrix}}\right]^{T}={\begin{pmatrix}(a^{\pi (1)})^{T}\\\vdots \\(a^{\pi (n)})^{T}\end{pmatrix}}^{T}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/af85693684d30d3ef64bfa65333595c28df915b2)
![{\displaystyle \left[I+f^{k}(e^{k})^{T}\right]\underbrace {\left[I-f^{k}(e^{k})^{T}\right]} _{=F_{k}}=I+f^{k}(e^{k})^{T}-f^{k}(e^{k})^{T}-f^{k}\underbrace {\left[(e^{k})^{T}f^{k}\right]} _{=0}(e^{k})^{T}=I}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e870e7ac1f00f358d93fcc7a4245c6a30dd7ba15)
![{\displaystyle F_{1}^{-1}\cdots F_{k}^{-1}F_{k+1}^{-1}=\left[I+\sum _{j=1}^{k}f^{j}(e^{j})^{T}\right]\left[I+f^{k+1}(e^{k+1})^{T}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/afd246c821197f1e30fb6e011ed727cc42011212)
![{\displaystyle =I+f^{k+1}(e^{k+1})^{T}+\sum _{j=1}^{k}f^{j}(e^{j})^{T}+\sum _{j=1}^{k}f^{j}\underbrace {\left[(e^{j})^{T}f^{k+1}\right]} _{=0}(e^{k+1})^{T}=I+\sum _{j=1}^{k+1}f^{j}(e^{j})^{T}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/95f6231d92b264c225b270ae5e4b899aa4ede67a)
![{\displaystyle PF_{k}P=P\left[I-f^{k}(e^{k})^{T}\right]P=\underbrace {P^{2}} _{=I}-[Pf^{k}][\underbrace {(e^{k})^{T}P} _{=(e^{k})^{T}}],}](https://wikimedia.org/api/rest_v1/media/math/render/svg/5cd1b573f617436a4fc459d196b1bc495b6ddca6)
![{\displaystyle \operatorname {H} _{f}(x):=\left({\frac {\partial ^{2}f}{\partial x_{i}\partial x_{j}}}(x)\right)_{i,j=1,\dots ,n}={\begin{pmatrix}{\frac {\partial ^{2}f}{\partial x_{1}\partial x_{1}}}(x)&{\frac {\partial ^{2}f}{\partial x_{1}\partial x_{2}}}(x)&\cdots &{\frac {\partial ^{2}f}{\partial x_{1}\partial x_{n}}}(x)\\[0.5em]{\frac {\partial ^{2}f}{\partial x_{2}\partial x_{1}}}(x)&{\frac {\partial ^{2}f}{\partial x_{2}\partial x_{2}}}(x)&\cdots &{\frac {\partial ^{2}f}{\partial x_{2}\partial x_{n}}}(x)\\\vdots &\vdots &\ddots &\vdots \\{\frac {\partial ^{2}f}{\partial x_{n}\partial x_{1}}}(x)&{\frac {\partial ^{2}f}{\partial x_{n}\partial x_{2}}}(x)&\cdots &{\frac {\partial ^{2}f}{\partial x_{n}\partial x_{n}}}(x)\end{pmatrix}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/04d705e48acb7356d6f3cf9dd39c135f342ed55a)
![{\displaystyle 0<(x^{(k)})^{T}Ax^{(k)}=\lambda _{k}\left[(x^{(k)})^{T}x{(k)}\right],\quad k=1,...,n}](https://wikimedia.org/api/rest_v1/media/math/render/svg/88b7d6990937c77661f541ab33cd3fbb99f5bff5)
![{\displaystyle \det(A)=\det {\begin{pmatrix}L_{n-1}&0\\c^{T}&\alpha \end{pmatrix}}\det {\begin{pmatrix}L_{n-1}^{T}&c\\0&\alpha \end{pmatrix}}=\alpha ^{2}[\det(L_{n-1})]^{2},}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fc01a06be683cc8ccb1b57d2f6440a6abd10283f)
![{\displaystyle \sum _{k=1}^{n}\left[k+\sum _{i=k+1}^{n}k\right]=\sum _{k=1}^{n}[k+(n-k)k]=\sum _{k=1}^{n}k+n\sum _{k=1}^{n}k-\sum _{k=1}^{n}k^{2}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/24b9edaba99ab30f2104df91489f6376c44e9bd2)
