Dekomposisi LU dengan Pivoting Parsial

Faktorisasi

Dekomposisi $\mathrm{LU}$ menyimpan eliminasi Gauss dalam dua matriks segitiga. Dengan pivoting parsial, matriks persegi nonsingular $A \in \mathbb{R}^{n\times n}$ difaktorkan sebagai

PA=LU

Ketiga faktor memiliki peran berikut:

$P$ adalah matriks permutasi yang mencatat pertukaran baris
$L$ adalah matriks segitiga bawah dengan diagonal satu
$U$ adalah matriks segitiga atas

Pada langkah eliminasi $k$ , pivoting parsial memilih entri dengan nilai mutlak terbesar pada kolom $k$ yang berada pada atau di bawah diagonal. Langkah ini $mencegah pembagian dengan pivot nol$ dan biasanya membatasi pertumbuhan galat numerik.

Untuk matriks permutasi, berlaku $P^{-1}=P^{\mathsf T}$ . Kesamaan $P^{-1}=P$ hanya berlaku pada permutasi khusus, misalnya satu pertukaran baris, dan tidak benar untuk setiap permutasi gabungan.

Cara eliminasi membentuk faktor

Setelah memilih pivot pada kolom $k$ , hilangkan setiap entri di bawahnya dengan pengali

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

Kemudian terapkan pembaruan baris berikut.

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

Pengali disimpan di bawah diagonal $L$ . Jika dua baris ditukar setelah beberapa langkah eliminasi, $L$ juga harus ditukar. Jika tidak, kesamaan $PA=LU$ tidak lagi terpenuhi.

Memfaktorkan Matriks dengan Pivoting Parsial

Kita akan mencatat $kedua pertukaran baris$ agar permutasi akhir sesuai dengan faktor.

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

Pivot pertama

Entri terbesar pada kolom pertama memiliki magnitudo $2$ . $Tukar baris pertama dan ketiga$ , lalu hilangkan entri baris kedua dengan pengali $\ell_{21}=1/2$ . Baris ketiga sudah diawali nol.

\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}

Pivot kedua

Pivot berikutnya memiliki magnitudo $2$ di baris ketiga, sehingga $tukar baris kedua dan ketiga$ . Tukar juga pengali kolom pertama yang sudah disimpan pada $L$ . Pengali untuk baris terakhir kemudian menjadi $\ell_{32}=(1/2)/2=1/4$ .

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

Faktor yang dihasilkan adalah

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

Perkalian kedua faktor memverifikasi dekomposisi:

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

Menyelesaikan sistem linear

Dari $PA=LU$ , persamaan $Ax=b$ berubah menjadi

LUx=Pb

Kita menyelesaikan sistem ini dalam $dua tahap$ menggunakan matriks segitiga.

Substitusi maju

Pertama, selesaikan $Ly=Pb$ . Pada contoh ini,

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

Untuk matriks segitiga bawah umum,

\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}

$Entri diagonal$ $L$ pada dekomposisi ini semuanya satu.

Substitusi mundur

Kemudian selesaikan $Ux=y$ dari $baris terakhir ke atas$ :

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

Untuk contoh tersebut,

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

Substitusi langsung mengonfirmasi $Ax=b$ .

Biaya dan penggunaan ulang

Faktorisasi $\mathrm{LU}$ padat membutuhkan sekitar

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

Jumlah tersebut menyatakan operasi titik mengambang ketika perkalian dan pengurangan dihitung secara terpisah. Setelah faktor tersedia, setiap ruas kanan baru hanya membutuhkan $O(n^2)$ kerja untuk dua substitusi. $Pemakaian ulang$ faktor $\mathrm{LU}$ menghindari faktorisasi $O(n^3)$ untuk setiap ruas kanan.

Determinan juga dapat diperoleh dari faktor:

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

Implementasi NumPy dengan Pemeriksaan Pivot

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))

Pemeriksaan tersebut menjadi bagian dari contoh karena $urutan faktor$ dan pencatatan pertukaran baris merupakan sumber kesalahan implementasi yang umum.

Toleransi pivot dihitung relatif terhadap norma tak hingga matriks. Karena itu, penskalaan seragam pada sistem berkondisi baik tidak menyebabkan kegagalan hanya karena semua entri kecil. Ambang ini merupakan keputusan numerik, bukan bukti singularitas eksak.

Batas penggunaan metode

Implementasi ini hanya menyelesaikan sistem dengan $matriks persegi yang dapat dibalik$ . Program menganggap pivot numerik nol sebagai kegagalan. Gunakan $\mathrm{QR}$ untuk masalah kuadrat terkecil dengan matriks persegi panjang. Untuk matriks yang tidak berperingkat penuh, gunakan $\mathrm{QR}$ dengan pivot kolom atau $\mathrm{SVD}$ . Faktorisasi $\mathrm{PLU}$ masih dapat dibuat untuk matriks persegi singular, tetapi program ini tidak dapat memakainya untuk memperoleh satu solusi unik.