Man zerlege die Koeffizientenmatrix A in zwei Dreiecksmatrizen L und U, wobei L eine untere und U eine obere Dreiecksmatrix darstellt, so dass gilt: A = L⋅U und L auf der Hauptdiagonalen nur Einsen stehen hat. Dann kann in zwei Schritten das Gleichungssystem Ax=b gelöst werden. zuerst bestimmt man die Lösung zum Gleichungssystem L⋅z=b und dann U⋅x=z. Da die Matrizen L und U in Dreiecksform vorliegen, lässt sich die jeweilige Rückwertsersetzeung recht schnell durchführen.
Eine nähere Beschreibung dieser Faktorisierung findet man auf Wikipedia

• A = input(' Matrix A eingeben: ');
• b = input(' Vektor b eingeben: ');
• b2=b;
• n=length(b);
•  
• %Permutationsvektor
• p=[1:n];
•  
• % In Dreiecksform umwandeln
• for k=1:n-1
• % Zeilenpivotisierung
• m=k;
• for l=k:n
• if(abs(A(l,k))>abs(A(m,k)))
• m=l;
• end
• end
• if(m~=p(k))
• t = p(k);
• p(k) = p(m);
• p(m) = t;
• end
•  
• % Elemination
• for i=k+1:n
• if(A(p(k),k) == 0)
• stop;
• end
• A(p(i),k) = A(p(i),k) / A(p(k),k);
• b(p(i)) = b(p(i)) - b(p(k))*A(p(i),k);
• for j=k+1:n
• A(p(i),j) = A(p(i),j) - A(p(k),j)*A(p(i),k);
• end
• end
• end
•  
• A = A(p(1:n),:);
• b = b2(p(1:n));
•  
• % Rueckwertssubstitution in zwei Schritten
• z=zeros(n,1);
• x=zeros(n,1);
• z(1) = b(1);
• for i=2:n
• z(i) = (b(i)-A(i,1:i-1)*z(1:i-1));
• end
•  
• for i=n:-1:1
• x(i) = (z(i)-A(i,i+1:n)*x(i+1:n)) / A(i,i);
• end
• x

Da die Matrizen L und U direkt aus der durch den Gauß-Algorithmus entstehenden Dreiecksmatrix (wenn man die Eliminationsfaktoren in die Matrix schreibt) direkt ablesbar sind, werden die beiden Matrizen nicht extra gespeichert. Dies vermindert den Speicherverbrauch.

Der Aufwand des Verfahrens ist nicht kleiner als der des Gauß-Algorithmus. Allerdings hat es den Vorteil, dass wenn man mehrere Gleichungssystem mit verschiedenen rechten Seiten b lösen muss, die Faktorisierung nicht nochmals durchgeführt werden braucht, sondern nur die Rückwärtssubstitution. Diese Situation tritt in vielen numerischen Algorithmen auf, z.B. beim Lösen nichtlinearer Gleichungssystem mit einem vereinfachten Newtonverfahren, oder beim iterativen berechnen von Eigenwerten einer Matrix. Beim vereinfachten Newtonverfahren führt man eine Iteration nach der folgenden Vorschrift aus, solange bis die gewünschte Genauigkeit erreicht ist:

xk+1 = xk - [f(x0)']-1⋅f(xk)

Wobei man in dem Verfahren nicht die Inverse der Jacobimatrix berechnen sollte, sondern das Gleichungssystem f(x0)'⋅Δx = -f(xk) löst und dann xk+1=xk+Δx setzt. Da man in jedem Iterationsschritt die selbe Koeffizientenmatrix f(x0)' verwendet, lohnt es sich eine Faktorisierung der Matrix zu bestimmen und dann in jedem Schritt nur noch 2 Rücksubstitutionen zu berechnen sind mit einem Aufwand von O(n2).

Es gibt noch andere Faktorisierungen der Matrix A, die auch mit einem kleineren Aufwand bestimmt werden können, oder die andere/besondere Eigenschaften haben. Siehe dazu z.B. Cholesky-Zerlegung oder QR-Faktorisierung