LU-Zerlegung mit partieller Pivotisierung

Die 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

$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, und nicht für jede zusammengesetzte Permutation.

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

und der Zeilenaktualisierung

u_{ij}\leftarrow u_{ij}-\ell_{ik}u_{kj},\qquad j=k,\ldots,n

eliminiert. 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.

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

In der ersten Spalte liegt der größte Betrag $2$ in der dritten Zeile. Wir vertauschen die erste und dritte Zeile. Der erste Eliminationsmultiplikator ist $\ell_{31}=1/2$ ; der andere Eintrag unterhalb des Pivots ist bereits null.

Der größte verbleibende Eintrag in der zweiten Spalte ist $2$ , daher ist keine zweite Vertauschung nötig. Der Multiplikator für die letzte Zeile ist $\ell_{32}=1/4$ .

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

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

Wir lösen sie in zwei Dreiecksschritten.

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

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

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

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}

und direktes Einsetzen bestätigt $Ax=b$ .

Eine dichte 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. Deshalb ist LU besonders nützlich, wenn viele Systeme dieselbe Matrix verwenden.

Auch die Determinante folgt aus den Faktoren:

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

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, dtype=float)
  if A.ndim != 2 or A.shape[0] != A.shape[1]:
      raise ValueError("PLU factorization requires a square matrix")

  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]) <= tolerance:
          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]) <= tolerance:
      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)
  for i in range(len(y) - 1, -1, -1):
      if abs(U[i, i]) <= tolerance:
          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)
  y = forward_substitution(L, P @ right_hand_side)
  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.

Diese Implementierung akzeptiert bewusst nur invertierbare quadratische Matrizen. Ein numerischer Nullpivot wird gemeldet und nicht stillschweigend übersprungen. Für rechteckige Ausgleichsprobleme oder rangdefiziente Matrizen verwendet man eine rangenthüllende QR-Zerlegung oder eine SVD. Diese Probleme besitzen andere Ausgabeverträge und sollten nicht in einem quadratischen PLU-Löser verborgen werden.

LU-Zerlegung mit partieller Pivotisierung

Die Faktorisierung

LU-Zerlegung mit partieller Pivotisierung

Die Faktorisierung

Wie die Elimination die Faktoren erzeugt

Vollständig durchgerechnetes Beispiel

Erster Pivot

Zweiter Pivot

Ein lineares System lösen

Vorwärtssubstitution

Rückwärtssubstitution

Aufwand und Wiederverwendung

Abgesicherte NumPy-Implementierung

Gültigkeitsbereich der Methode