Matrizen mit dem LU Verfahren und partieller Pivotisierung zerlegen

Die Faktorisierung

Die $\mathrm{LU}$ -Zerlegung speichert die Gauß-Elimination in zwei Dreiecksmatrizen. Mit partieller Pivotisierung wird eine invertierbare quadratische Matrix $A \in \mathbb{R}^{n\times n}$ faktorisiert als

PA=LU

Die drei Faktoren haben folgende Aufgaben:

$P$ eine Permutationsmatrix ist, welche die Zeilenvertauschungen erfasst
$L$ eine untere Dreiecksmatrix mit Einsen auf der Diagonalen ist
$U$ eine obere Dreiecksmatrix ist

Im Eliminationsschritt $k$ wählt die partielle Pivotisierung in Spalte $k$ den betragsmäßig größten Eintrag auf oder unterhalb der Diagonalen. Dadurch wird $eine Division durch einen Nullpivot vermieden$ und das Wachstum numerischer Fehler gewöhnlich begrenzt.

Für eine Permutationsmatrix gilt $P^{-1}=P^{\mathsf T}$ . Die Gleichheit $P^{-1}=P$ gilt nur für spezielle Permutationen, etwa eine einzelne Zeilenvertauschung. Bei einer zusammengesetzten Permutation kann sie falsch sein.

Wie die Elimination die Faktoren erzeugt

Nach der Pivotwahl in Spalte $k$ wird jeder Eintrag darunter mit dem Multiplikator

\ell_{ik}=\frac{u_{ik}}{u_{kk}},\qquad i=k+1,\ldots,n

Anschließend wenden wir die folgende Zeilenaktualisierung an.

\begin{aligned} &u_{ij}\leftarrow u_{ij}-\ell_{ik}u_{kj}, \\ &j=k,\ldots,n \end{aligned}

Diese Aktualisierung eliminiert den Eintrag unterhalb des Pivots. Die Multiplikatoren werden unterhalb der Diagonalen von $L$ gespeichert. Werden nach früheren Eliminationsschritten zwei Zeilen vertauscht, müssen auch die $bereits gefüllten Spalten$ von $L$ vertauscht werden. Andernfalls gilt $PA=LU$ nicht mehr.

Eine Matrix mit partieller Pivotisierung faktorisieren

Wir verfolgen $beide Zeilenvertauschungen$ , damit die endgültige Permutation zu den Faktoren passt.

A=\begin{pmatrix}0&2&1\\1&1&0\\2&1&1\end{pmatrix}

Erster Pivot

Der größte Betrag in der ersten Spalte ist $2$ . $Vertausche die erste und dritte Zeile$ und eliminiere dann den Eintrag der zweiten Zeile mit $\ell_{21}=1/2$ . Die dritte Zeile beginnt bereits mit null.

\begin{pmatrix}2&1&1\\1&1&0\\0&2&1\end{pmatrix}\longrightarrow\begin{pmatrix}2&1&1\\0&1/2&-1/2\\0&2&1\end{pmatrix}

Zweiter Pivot

Der nächste Pivot hat den Betrag $2$ und steht in Zeile drei. $Vertausche deshalb die zweite und dritte Zeile$ , einschließlich ihrer bereits gespeicherten Multiplikatoren in der ersten Spalte von $L$ . Für die letzte Zeile gilt danach $\ell_{32}=(1/2)/2=1/4$ .

u_{33}=-\frac12-\frac14\cdot1=-\frac34

Damit ergeben sich die Faktoren

P=\begin{pmatrix}0&0&1\\1&0&0\\0&1&0\end{pmatrix}

Die Multiplikation der Faktoren bestätigt die Zerlegung:

LU=\begin{pmatrix}2&1&1\\0&2&1\\1&1&0\end{pmatrix}=PA

Ein lineares System lösen

Aus $PA=LU$ wird die Gleichung $Ax=b$ zu

LUx=Pb

Wir lösen dieses System in $zwei Schritten$ mit Dreiecksmatrizen.

Vorwärtssubstitution

Zuerst lösen wir $Ly=Pb$ . Hier gilt

Pb=\begin{pmatrix}7\\7\\3\end{pmatrix}

Für eine allgemeine untere Dreiecksmatrix lautet die Rekursion

\begin{aligned} &y_i=\frac{1}{\ell_{ii}}\left(c_i-\sum_{j=1}^{i-1}\ell_{ij}y_j\right), \\ &c=Pb \end{aligned}

Die $Diagonaleinträge$ dieses $L$ sind alle eins.

Rückwärtssubstitution

Anschließend lösen wir $Ux=y$ von der $letzten Zeile nach oben$ :

x_i=\frac{1}{u_{ii}}\left(y_i-\sum_{j=i+1}^{n}u_{ij}x_j\right)

Für das Beispiel erhalten wir

x=\begin{pmatrix}1\\2\\3\end{pmatrix}

Direktes Einsetzen bestätigt $Ax=b$ .

Aufwand und Wiederverwendung

Eine dichte $\mathrm{LU}$ -Faktorisierung benötigt ungefähr

\frac{2}{3}n^3+O(n^2)

Gleitkommaoperationen, wenn Multiplikationen und Subtraktionen getrennt gezählt werden. Sobald die Faktoren vorliegen, benötigt jede neue rechte Seite nur $O(n^2)$ Arbeit für die beiden Substitutionen. Die $Wiederverwendung derselben$ $\mathrm{LU}$ -Faktoren vermeidet daher eine erneute Faktorisierung mit $O(n^3)$ Aufwand für jede rechte Seite.

Auch die Determinante folgt aus den Faktoren:

\det(A)=\det(P)\prod_{i=1}^{n}u_{ii}

Implementierung mit NumPy und Pivotprüfung

plu_solve.py

import numpy as np


def plu_factor(matrix, tolerance=1e-12):
  """Return P, L, U such that P @ matrix == L @ U."""
  A = np.asarray(matrix)
  if np.iscomplexobj(A):
      raise ValueError("this example requires real entries")
  A = A.astype(float)
  if A.ndim != 2 or A.shape[0] != A.shape[1] or A.shape[0] == 0:
      raise ValueError("PLU factorization requires a nonempty square matrix")
  if not np.all(np.isfinite(A)):
      raise ValueError("matrix entries must be finite")
  if not np.isfinite(tolerance) or tolerance < 0:
      raise ValueError("relative tolerance must be finite and nonnegative")
  threshold = tolerance * np.linalg.norm(A, ord=np.inf)

  n = A.shape[0]
  U = A.copy()
  L = np.eye(n)
  P = np.eye(n)

  for k in range(n - 1):
      pivot = k + np.argmax(np.abs(U[k:, k]))
      if abs(U[pivot, k]) <= threshold:
          raise np.linalg.LinAlgError("matrix is singular to working precision")

      if pivot != k:
          U[[k, pivot]] = U[[pivot, k]]
          P[[k, pivot]] = P[[pivot, k]]
          L[[k, pivot], :k] = L[[pivot, k], :k]

      for i in range(k + 1, n):
          L[i, k] = U[i, k] / U[k, k]
          U[i, k:] -= L[i, k] * U[k, k:]

  if abs(U[-1, -1]) <= threshold:
      raise np.linalg.LinAlgError("matrix is singular to working precision")

  return P, L, U


def forward_substitution(L, right_hand_side):
  b = np.asarray(right_hand_side, dtype=float)
  y = np.empty_like(b)
  for i in range(len(b)):
      y[i] = (b[i] - L[i, :i] @ y[:i]) / L[i, i]
  return y


def backward_substitution(U, right_hand_side, tolerance=1e-12):
  y = np.asarray(right_hand_side, dtype=float)
  x = np.empty_like(y)
  threshold = tolerance * np.linalg.norm(U, ord=np.inf)
  for i in range(len(y) - 1, -1, -1):
      if abs(U[i, i]) <= threshold:
          raise np.linalg.LinAlgError("singular triangular system")
      x[i] = (y[i] - U[i, i + 1:] @ x[i + 1:]) / U[i, i]
  return x


def plu_solve(matrix, right_hand_side):
  P, L, U = plu_factor(matrix)
  b = np.asarray(right_hand_side)
  if np.iscomplexobj(b):
      raise ValueError("this example requires a real right-hand side")
  b = b.astype(float)
  if b.shape != (len(P),) or not np.all(np.isfinite(b)):
      raise ValueError("right-hand side must be a finite vector of matching size")
  y = forward_substitution(L, P @ b)
  x = backward_substitution(U, y)
  return x, P, L, U


A = np.array([[0, 2, 1], [1, 1, 0], [2, 1, 1]], dtype=float)
b = np.array([7, 3, 7], dtype=float)
x, P, L, U = plu_solve(A, b)

assert np.allclose(P @ A, L @ U)
assert np.allclose(A @ x, b)
assert np.allclose(x, np.array([1, 2, 3], dtype=float))

Die Prüfungen gehören zum Beispiel, weil eine falsche $Faktorordnung$ und fehlerhaftes Verbuchen von Zeilenvertauschungen häufige Implementierungsfehler sind.

Die Pivottoleranz bezieht sich auf die Unendlichnorm der Matrix. Eine gleichmäßige Skalierung eines gut konditionierten Systems löst deshalb nicht allein wegen kleiner Einträge einen Fehler aus. Der Schwellenwert ist eine numerische Entscheidung und kein Beweis exakter Singularität.

Gültigkeitsbereich der Methode

Diese Implementierung löst nur Gleichungssysteme mit $invertierbaren quadratischen Matrizen$ und behandelt einen numerischen Nullpivot als Fehler. Für rechteckige Ausgleichsprobleme eignet sich eine $\mathrm{QR}$ -Zerlegung. Bei rangdefizienten Problemen sind eine $\mathrm{QR}$ -Zerlegung mit Spaltenpivotisierung oder eine $\mathrm{SVD}$ erforderlich. Auch für eine singuläre quadratische Matrix kann eine $\mathrm{PLU}$ -Faktorisierung existieren, doch dieses Lösungsverfahren kann damit keine eindeutige Lösung bestimmen.