.. _LUdecomp: LU Decomposition ================ .. index:: lu A variant of Gaussian elimination called *LU decomposition* (for Lower–Upper) factors a matrix into products of two matrices. The first is a lower triangular matrix, and the second is an upper triangular matrix. The LU decomposition algorithm can be implemented in software to achieve the best possible accuracy from an elimination-based algorithm. [1]_ It is used internally by solver functions to solve square matrix systems of equations and to calculate determinants. The ``lu`` from the ``scipy.linalg`` module function finds a matrix’s LU decomposition. Then triangular solvers are used to find the solution to a system of equations, The ``lu_factor`` embeds the triangular factors in its returned data and then the ``lu_solve`` function uses those results to solve a system of equations. Recall that to solve a system of equations, we can use Gaussian elimination with an augmented matrix so that changes made to matrix :math:`\bf{A}` are also applied to vector :math:`\bm{b}`. .. math:: \begin{array}{lcl} \mathbf{A}\,\bm{x} = \bm{b} & \mapsto & \mathbf{U}\,\bm{x} = \bm{c} \end{array} In the modified equation, :math:`\bf{U}` is an upper triangular matrix for which simple back substitution may be used to solve for the unknown vector :math:`\bm{x}`. The LU decomposition method differs in that it operates only on the :math:`\bf{A}` matrix to find a matrix factoring, :math:`\mathbf{A} = \bf{L\,U}`, where :math:`\bf{L}` is lower triangular and :math:`\bf{U}` is upper triangular. Thus the substitution phase for solving :math:`\mathbf{A}\,\bm{x} = \bm{b}` has two steps—solving :math:`\mathbf{L}\,\bm{y} = \bm{b}` and then solving :math:`\mathbf{U}\,\bm{x} = \bm{y}`. LU decomposition changes the equations used to solve systems of equations as follows. In addition to :math:`\bf{L}` and :math:`\bf{U}`, we have an elimination matrix, :math:`\bf{E}`. Its inverse is :math:`\bf{L}`. We also introduce the :math:`\bf{P}` matrix, the permutation matrix that orders the rows to match the row exchanges. [2]_ .. math:: \begin{array}{lcl} \mathbf{A} & \mapsto & \mathbf{E\,P}^T \mathbf{A} = \mathbf{U} \\ \mathbf{P}^T\mathbf{A} = \mathbf{E}^{-1}\mathbf{U} & \mapsto & \mathbf{A} = \mathbf{P\,L\,U} \\ \mathbf{A}\,\bm{x} = \bm{b} & \mapsto & \mathbf{P\,L\,U}\,\bm{x} = \bm{b} \\ \bm{x} = \mathbf{A}^{-1}\bm{b} & \mapsto & \bm{x} =\mathbf{U}^{-1}\mathbf{L}^{-1}\mathbf{P}^T \bm{b} \end{array} .. math:: :label: eq-xlu \boxed{\bm{x} =\mathbf{U}^{-1}\mathbf{L}^{-1}\mathbf{P}^T \bm{b}} .. topic:: Which algorithm is faster? There are some applications where there is a need to solve several systems of equations with the same :math:`\mathbf{A}` matrix, but with different :math:`\bm{b}` vectors. It is faster to find the LU factors once and then use forward and backward substitution with each :math:`\bm{b}` vector than to do full elimination each time. The combined forward and backward substitution for LU decomposition has a run time of :math:`\mathcal{O}(n^2)`. Elimination with an augmented matrix and the LU decomposition algorithm have the same order of magnitude run time performance, which, depending on the exact implementation, is approximately :math:`\mathcal{O}(n^3)`. But more precisely, the cost of computation for LU decomposition is :math:`(2/3)n^3 + n^2` floating point operations (flops) [HIGHAM11]_, which explains why experimentation shows that LU decomposition is faster than the ``rref`` function, especially for larger matrices. Here is a quick experiment using MATLAB on a personal computer. :: >> A = rand(100); b = rand(100,1); >> tic; [L,U,P] = lu(A); x = U\(L\(P*b)); toc; Elapsed time is 0.000605 seconds. >> tic; x = rref([A b]); toc; Elapsed time is 0.069932 seconds. .. index:: LU decomposition .. _LU-example: LU Example ---------- Before demonstrating how to use the ``lu`` function, we will work on part of a problem and then turn to Python to help with the computation. We will use the same example solved in :ref:`elimination`. For purposes of understanding how the :math:`\bf{U}` and :math:`\bf{L}` matrices are factors of the :math:`\bf{A}` matrix, we will follow a traditional elimination strategy to find :math:`\bf{U}`. We will form :math:`\bf{L}` from a record of the row operations needed to convert :math:`\bf{A}` into :math:`\bf{U}`. We will not make row exchanges in this example, but will use row exchanges in :ref:`LUdecomp`. .. math:: \left\{\begin{aligned} 2x - 3y {\quad} {\quad} &= 3 \\ 4x - 5y + z {\:} &= 9 \\ 2x - y - 3z {\,} &= -1 \end{aligned}\right. The needed row operations are the same as noted in :ref:`elimination`. Each row operation is represented by a matrix :math:`\mathbf{E}_{i,j}` that can multiply with :math:`\bf{A}` to affect the change. These matrices start with an identity matrix and then change one element to represent each operation. We call these matrices *Elementary Matrices*. Then the product of the elementary matrices is called the *elimination matrix*. .. note:: The next few examples demonstrate how LU decompositions rather than how to use LU decomposition with Python, so the examples come from my previous MATLAB notes. #. Add :math:`-1` of *row 1* to *row 3*, called :math:`\mathbf{E}_{1}`. #. Add :math:`-2` of *row 1* to *row 2*, called :math:`\mathbf{E}_{2}`. #. Add :math:`-2` of *row 2* to *row 3*, called :math:`\mathbf{E}_{3}`. :: >> E1 E1 = 1 0 0 0 1 0 -1 0 1 >> E2 E2 = 1 0 0 -2 1 0 0 0 1 >> E3 E3 = 1 0 0 0 1 0 0 -2 1 The combined row operations can be found by multiplying them together. The order of matrices must be such that the first operation applied is next to :math:`\bf{A}` in the equation :math:`\mathbf{U} = \bf{E\,A}`. :: >> E = E3*E2*E1 E = 1 0 0 -2 1 0 -1 -2 1 >> U = E*A U = 2 -3 0 0 1 1 0 0 -5 We can also use an augmented matrix to find E and U. Start with :math:`\left[\,\mathbf{A}\,|\, \mathbf{I}\, \right]` and perform row operations until the matrix in the position of :math:`\bf{A}` is upper triangular. The resulting augmented matrix is then :math:`\left[\,\mathbf{U}\,|\, \mathbf{E}\, \right]` [KUTZ13]_. .. math:: \left[\begin{array}{ccc|ccc} 2 & -3 & 0 & 1 & 0 & 0 \\ 4 & -5 & 1 & 0 & 1 & 0 \\ 2 & -1 & -3 & 0 & 0 & 1 \end{array}\right] % \end{equation*}% % \begin{equation*} \; \mapsto \; \left[\begin{array}{ccc|ccc} 2 & -3 & 0 & 1 & 0 & 0 \\ 0 & 1 & 1 & -2 & 1 & 0 \\ 0 & 0 & -5 & -1 & -2 & 1 \end{array}\right] The :math:`\bf{L}` matrix is the inverse of :math:`\bf{E}`. Recall that to take the inverse of a product of matrices, we reverse the order of the matrices—:math:`(\mathbf{A\,B})^{-1} = \mathbf{B}^{-1}\mathbf{A}^{-1}`. Note that each elementary matrix only differs from the identity matrix at one element. The inverse of each elementary matrix is found by changing the sign of the additional nonzero value in the identity matrix. :: >> inv(E1) ans = 1 0 0 0 1 0 1 0 1 >> inv(E2) ans = 1 0 0 2 1 0 0 0 1 >> inv(E3) ans = 1 0 0 0 1 0 0 2 1 >> L = inv(E1)*inv(E2)*inv(E3) L = 1 0 0 2 1 0 1 2 1 >> inv(E) ans = 1.0000 0 0 2.0000 1.0000 0 1.0000 2.0000 1.0000 Also note that the :math:`\bf{L}` matrix does the reverse of :math:`\bf{E}`, so we could find :math:`\bf{L}` directly from the list of row operations. Now, with :math:`\bf{L}` and :math:`\bf{U}` found, we can find the solution to :math:`\mathbf{L\,U}\,\bm{x} = \bm{b}`. Since these are lower and upper triangular matrices, the left-divide operator will quickly solve them. If several solutions are needed, the one-time calculation of :math:`\bf{L}` and :math:`\bf{U}` will substantially expedite the calculations. :: In [17]: import numpy as np In [18]: from scipy import linalg In [19]: U = np.array([[2, -3, 0], [0, 1, 1], [0, 0, -5]], dtype=float) In [20]: L = np.array([[1, 0, 0], [2, 1, 0], [1, 2, 1]], dtype=float) In [21]: b = np.array([[3, 9, -1]], dtype=float).T In [22]: x = linalg.solve_triangular(U, linalg.solve_triangular(L, b, lower=True)); x Out[22]: array([[3.], [1.], [2.]]) .. note:: The substitution phase has two steps: solving :math:`\mathbf{L}\,\bm{y} = \bm{b}` and :math:`\mathbf{U}\,\bm{x} = \bm{y}`. Either another vector, :math:`\bm{y}`, is needed, or the triangular solving functions may be nested such that the returned vector :math:`\bm{y}` is passed as an argument to the solver finding :math:`\bm{x}`. .. index:: elimination matrix .. index:: elementary matrices Example with Row Exchanges -------------------------- This example shows the LU decomposition with partial pivoting (row exchanges). We will show the row exchanges with a permutation matrix. A permutation matrix has the rows of an identity matrix, but the order of the rows shows what row exchanges were performed. .. math:: \mathbf{A} = \begin{bmatrix} 0 & 12 & -3 \\ 8 & -4 & -6 \\ -4 & -2 & 12 \end{bmatrix} Looking at the first column, we will use the largest value, 8, as our first pivot, so we exchange rows one and two and define a permutation matrix. The :math:`\mathbf{U}` matrix will reflect the row exchanges. The :math:`\mathbf{L}` matrix begins as an identity matrix, and the ones on the diagonal will remain. We will document our elimination row operations directly in the :math:`\mathbf{L}` matrix as the negative values of what we would add to the elementary matrices. .. math:: \mathbf{P} = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix} .. math:: \begin{array}{ll} \mathbf{L} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix} &\mathbf{U} = \begin{bmatrix} 8 & -4 & -6 \\ 0 & 12 & -3 \\ -4 & -2 & 12 \end{bmatrix} \end{array} Add :math:`1/2` of *row 1* to *row 3*. .. math:: \begin{array}{ll} \mathbf{L} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -1/2 & 0 & 1 \end{bmatrix} &\mathbf{U} = \begin{bmatrix} 8 & -4 & -6 \\ 0 & 12 & -3 \\ 0 & -4 & 9 \end{bmatrix} \end{array} We will keep the pivot of :math:`12` in *row 2*, so no change to the permutation matrix is needed. Add :math:`1/3` of *row 2* to *row 3*. .. math:: \begin{array}{ll} \mathbf{L} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -1/2 & -1/3 & 1 \end{bmatrix} &\mathbf{U} = \begin{bmatrix} 8 & -4 & -6 \\ 0 & 12 & -3 \\ 0 & 0 & 8 \end{bmatrix} \end{array} Python Examples of LU --------------------- We will accept a permutation matrix as an output matrix. The relationship of the matrices is :math:`\mathbf{A} = \mathbf{P\,L\,U}`. Note that the permutation matrix merely swaps the rows of whatever vector or matrix it is multiplied by. Like the identity matrix, permutation matrices are orthogonal, thus their inverse is just their transpose. So the solution comes from equation :eq:`eq-xlu`. In Python, it is necessary to use two triangular solvers. :: In [1]: import numpy as np In [2]: from scipy import linalg In [3]: A = np.random.randint(-10, 15, (4, 4)).astype(float) In [4]: P, L, U = linalg.lu(A) In [5]: P Out[5]: array([[1., 0., 0., 0.], [0., 1., 0., 0.], [0., 0., 0., 1.], [0., 0., 1., 0.]]) In [6]: L Out[6]: array([[ 1. , 0. , 0. , 0. ], [ 0.69230769, 1. , 0. , 0. ], [-0.30769231, -0.40540541, 1. , 0. ], [ 0.76923077, 0.48648649, 0.60566449, 1. ]]) In [7]: U Out[7]: array([[ 13. , -6. , -2. , 12. ], [ 0. , -2.84615385, 12.38461538, 3.69230769], [ 0. , 0. , 12.40540541, 7.18918919], [ 0. , 0. , 0. , -17.38126362]]) In [8]: x = linalg.solve_triangular(U, linalg.solve_triangular(L, (P.T @ b), lower=True)); x Out[8]: array([[2.61017799], [4.70569065], [0.4728002 ], [0.2706192 ]]) # Verify solution In [8]: y = linalg.solve(A, b); y Out[8]: array([[2.61017799], [4.70569065], [0.4728002 ], [0.2706192 ]]) Another option is to use ``lu_factor`` and ``lu_solve``. As we see from the example below. The ``lu_factor`` function returns two arrays in a tuple. The first array is the values of ``L`` and ``U`` combined. The second array holds the pivot indices representing the permutation matrix :math:`\bf{P}`. :: In [16]: lu_and_pivot = linalg.lu_factor(A) In [17]: lu_and_pivot Out[17]: (array([[ 13. , -6. , -2. , 12. ], [ 0.69230769, -2.84615385, 12.38461538, 3.69230769], [ -0.30769231, -0.40540541, 12.40540541, 7.18918919], [ 0.76923077, 0.48648649, 0.60566449, -17.38126362]]), array([0, 1, 3, 3], dtype=int32)) In [18]: x = linalg.lu_solve(lu_and_pivot, b); x Out[18]: array([[2.61017799], [4.70569065], [0.4728002 ], [0.2706192 ]]) .. _LUdet: Determinant Shortcut -------------------- Another application of LU decomposition is a faster algorithm for computing the determinant of a matrix. The product of the determinants of a matrix factorization is equal to the determinant of the original matrix. Thus, .. math:: |\mathbf{A}| = |\mathbf{P}|\,|\mathbf{L}|\,|\mathbf{U}|. Because of the zeros in all three of these matrices, their determinants are quick to compute. The determinant of :math:`\bf{P}` is either 1 or -1. Since each row has only a single one, the determinant is quick to calculate. The determinant of :math:`\bf{L}` is always 1, which is the product of the ones along the diagonal. The determinant of :math:`\bf{U}` is the product of the values on its diagonal. .. index:: det Refer to the normal method for calculating the determinant in :ref:`determinant` and take a moment to see why the determinant of an upper or lower triangular matrix is the product of the diagonal values. The determinant of a matrix is computed as follows. :: In [33]: P, L, U = linalg.lu(A) In [34]: determinant = linalg.det(P) * np.prod(np.diag(U)) .. index:: determinant .. rubric:: Footnotes .. [1] Matrices in a special form, such as diagonal, triangular, or positive definite, have more efficient solvers. .. [2] Equation :eq:`eq-xlu` is expressed with a pure mathematics notation using the inverse of matrices. The expression that is preferred with MATLAB is an expression that directly uses back substitution on triangular matrices (``x = U\(L\(P*b))``). Note that the parentheses are important to ensure that substitution is applied in the correct order.