LRP-Zerlegung

4. Beseitigung des Schönheitsfehlers: LRP-Zerlegung

Wir kommen nun auf das Problem der Division durch Null (vom Ende des 1. Abschnitts) zurück. Bei der folgenden Matrix versagt unser Algorithmus, obwohl sie invertierbar ist:

~~A = \begin{pmatrix} 0 & 5 & 2 \cr 2 & 6 & 4 \cr 2 & 1 & 1 \end{pmatrix}~~

Wir können die Division durch Null hier verhindern, indem wir die ersten beiden Zeilen der Matrix vertauschen:

~~A' = \begin{pmatrix} 2 & 6 & 4 \cr 0 & 5 & 2 \cr 2 & 1 & 1 \end{pmatrix}~~

Die Vertauschung der der ersten beiden Zeilen lässt sich erreichen, indem man A von links mit der entsprechenden Permutationsmatrix multipliziert. Diese entsteht aus der Einheitsmatrix durch Vertauschung der der entsprechenden Zeilen:

A' = \begin{pmatrix} 2 & 6 & 4 \cr 0 & 5 & 2 \cr 2 & 1 & 1 \end{pmatrix} =\begin{pmatrix} 0 & 1 & 0 \cr 1 & 0 & 0 \cr 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 0 & 5 & 2 \cr 2 & 6 & 4 \cr 2 & 1 & 1 \end{pmatrix}

In einer invertierbaren Matrix lassen sich die Zeilen immer so permutieren, dass unser Algorithmus anschließend erfolgreich ist. Es gilt nämlich:

Satz

Zu jeder invertierbaren Matrix A gibt es

eine Permutationsmatrix P,
eine linke Dreiecksmatrix L (nur Einsen auf der Hauptdiagonalen),
eine rechte Dreiecksmatrix R (keine Null auf der Hauptdiagonalen),

so dass PA = LR.

Wir können das lineare Gleichungssystem auf beiden Seiten mit P multiplizieren, und dann mit der Zerlegung der der permutierten Matrix arbeiten:

~~A\textbf x = \textbf b \quad \Rightarrow\quad \green{PA}\textbf x = P\textbf b \quad \Rightarrow\quad \green{LR}\textbf x = P\textbf b~~

Dieses Gleichungssystem lässt sich in zwei Schritten lösen:

1. Schritt:

Mit dem Ansatz ~~\textbf x = \textbf y~~ wird die Gleichung zu ~~L\textbf y = P\textbf b~~. Diese Gleichung nach ~~\textbf y~~ auflösen!

2. Schritt:

Mit dem nun bekannten ~~\textbf y~~ kann der Ansatz ~~R\textbf x = \textbf y~~ nach x aufgelöst werden.

In der nun folgenden Implementierung vertauschen wir Zeilen nicht nur, um Divisionen durch Null zu verhindern, sondern bringen immer das dem Betrage nach größte Element an die Divisionsposition. Damit erhöht sich die numerische Stabilität beim Rechnen mit Gleitkommazahlen.

function lup(A) {
	let C = A.copy() // Kopie, damit übergebene Matrix unverändert bleibt
	const n = C.height
	let P = Matrix.identity(n) // Ausgangsbasis für Permutationsmatrix
	for (k of fromTo(0,n-2)) {
		// Suche das dem Betrag nach größte Element
		// der k-ten Spalte ab der k-ten Zeile:
		const [absMax, kMax] = fromTo(k,n-1)
			.map(i => [abs(C[i][k]), i])
			.reduce((a,b) => greater(a[0], b[0]) ? a : b)
		assert(greater(absMax, 0))

// Vertausche die Zeile des größten Elements mit der k-ten Zeile...
		P.swap(k, kMax)	// ... in der Permutation
		C.swap(k, kMax)	// ... und in der Matrix

for (const i of fromTo(k+1,n-1)) {
			// Dividiere die restlichen Elemente
			// der k-ten Spalte durch das „Kopfelement“:
			C[i][k] = div(C[i][k], C[k][k])
			// Subtrahiere in der Restmatrix überall das Produkt
			// aus den „Randelementen“ links und oberhalb:
			for (j of fromTo(k+1,n-1))
				C[i][j] = sub(C[i][j], mul(C[i][k], C[k][j]))
		}
	}
	// C in linke und rechte Dreiecksmatrix zerlegen
	const L = Matrix.fromFunction(n, n, (i,k) => (k<i) ? C[i][k] : (k==i) ? 1 : 0)
	const R = Matrix.fromFunction(n, n, (i,k) => (k>=i) ? C[i][k] : 0)
	return [L, R, P]
}

Die Funktion lup ist in mathGUIde bereits eingebaut, allerdings nicht als globale Funktion sondern als Matrix-Methode. Sie wird demnach mit A.lup() statt lup(A) aufgerufen.