LUP decomposition







Numerical Methods

LUP decomposition

The stability of LU decomposition is improved if pivoting is used to maximize the absolute values of the diagonal elements of the upper triangular matrix $U$. The function LUP_decomp(A) performs LU-decomposition with partial pivoting. It is the same as the function LU_matrix(A) except a pivoting section has been added. It returns an object consisting of the LU matrix, the permutation matrix, and the number of row exchanges made during partial pivoting.

<script> function LUP_decomp(A) { //LU decomposition with partial pivoting //returns the LU matrix, the permutation matrix of matrix A, and the number of row swaps //This function overwrites A N = A.length; P = identity(N); exchanges = 0; //count the number of row swaps for (i=0; i<N; i++) { //start pivot section Umax = 0; for (r=i; r<N; r++) { Uii=A[r][i]; q = 0; while (q<i) { Uii -= A[r][q]A[q][r]; q++; } if (Math.abs(Uii)>Umax) { Umax = Math.abs(Uii); row = r; } } if (i!=row) {//swap rows exchanges++; for (q=0; q<N; q++) { tmp = P[i][q]; P[i][q]=P[row][q]; P[row][q]=tmp; tmp = A[i][q]; A[i][q]=A[row][q]; A[row][q]=tmp; } } //end pivot section j = i; while (j<N) { //Determine U across row i q = 0; while (q<i) { A[i][j] -= A[i][q]A[q][j]; q++; } j++; } j = i+1; while (j<N) { //Determine L down column i q = 0; while (q<i) { A[j][i] -= A[j][q]*A[q][i]; q++; } A[j][i] = A[j][i]/A[i][i]; j++; } } return {LU:A,P:P,exchanges:exchanges}; } function identity(N) { // returns an NxN identity matrix var I = []; for (j=0; j<N; j++) { I[j] = new Array(N); for (k=0; k<N; k++) { I[j][k] = 0; } I[j][j]=1; } return I; } function write_matrix(A) { rows = A.length; cols = A[0].length; for (i=0; i<rows; i++) { for (j=0; j<cols; j++) { document.write(A[i][j]+" "); } document.write("<\/br> "); } document.write("<\/br> "); } var A = [[2,7,6,2],[9,5,1,3],[4,3,8,4],[5,6,7,8]]; document.write("A = <br \/>"); write_matrix(A); LUP = LUP_decomp(A); document.write("LU = <br />"); write_matrix(LUP.LU); document.write("P = <br />"); write_matrix(LUP.P); </script>

It is clear that with the strategy of partial pivoting we avoid the problem caused by the presence of a zero element on the diagonal of $U$. The only exception occurs if during the calculation all the $U_{ij},\quad i=j,\ldots,n$ are simultaneously zero. In that case the matrix of coefficients is singular and the system of equations cannot be solved with this method.

Matrix inversion

The inverse of matrix A is defined by the equation,

\begin{equation} A \cdot A^{-1} = I, \end{equation}

where $I$ is the identity matrix. The column vectors of $A^{-1}$ can be found by solving $A\cdot\vec{x}=\vec{b}$ where $\vec{b}$ is the corresponding column vector of the identity matrix. We use LUP-decomposition to find the permutation matrix $P$ and the upper and lower triangular matricies $L$ and $U$ such that $P\cdot A= L\cdot U$. Then we can call LU_solve($LU$,$P\cdot I$) to find $A^{-1}$.

function inverse(A) { // returns the inverse of a square matrix
  LUP = LUP_decomp(A);
  Pb = dot(LUP.P,identity(A.length));
  Ainv = LU_solve(LUP.LU,Pb);
  return Ainv;
}

Calculation of the determinant of a matrix

The determinant of the product of two matrices is the product of their determinants. After LUP-decomposition we have, $\text{det}(P)\text{det}(A)=\text{det}(L)\text{det}(U)$. The determinant of the permutation matrix $P$ is $-1^n$ where $n$ is the number of row exchanges made during partial pivoting. The determinant of an upper or lower triangular matrix is the product of its diagonal elements: $\text{det}(L)=1$, $\text{det}(U)=\prod\limits_{i=1}^NU_{ii}$. The following code will calculate the determinant.

function det(A) { //returns the determinant of square matrix A
  LUP = LUP_decomp(A);
  D = 1;
  for (i=0; i < A.length; i++) {//product the diagonal elements of the LU matrix
    D = D*LUP.LU[i][i]; 
  }
  D = D*LUP.sgn;
  return D;
}

LUP-decomposition

This app performs LU decomposition of a square matrix with or without partial pivoting. It can solve a set of linear inhomogeneous equations, perform matrix multiplication, and find the determinant, transpose, or inverse of a matrix.