Talk:LU decomposition

2004

Computer Science Student 13:52, 17 February 2013 (UTC)Algorithm 2 does not appear to me to be an algorithm, but rather a description of the goal of an algorithm. Too Old 17:17, 16 Jul 2004 (UTC)

I resumed work on this article. The following sections were ripped nearly verbatim from [[1]] so I deleted them.

Here is a direct method for decomposing an $N\times N$ matrix $\mathbf {A}$ into the product of a lower triangular matrix $\mathbf {L}$ and an upper triangular matrix $\mathbf {U}$ ,

\mathbf {L} \,\mathbf {U} =\mathbf {A}

For a $3\times 3$ matrix, this becomes:

{\begin{bmatrix}l_{11}&0&0\\l_{12}&l_{22}&0\\l_{13}&l_{23}&l_{33}\\\end{bmatrix}}\cdot {\begin{bmatrix}u_{11}&u_{12}&u_{13}\\0&u_{22}&u_{23}\\0&0&u_{33}\\\end{bmatrix}}={\begin{bmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\\\end{bmatrix}}

This gives three types of equations

i<j\qquad l_{i1}u_{1j}+l_{i2}u_{2j}+\dots +l_{ii}u_{ij}=a_{ij}

i=j\qquad l_{i1}u_{1j}+l_{i2}u_{2j}+\dots +l_{ii}u_{jj}=a_{ij}

i>j\qquad l_{i1}u_{1j}+l_{i2}u_{2j}+\dots +l_{ii}u_{jj}=a_{ij}

This gives $N^{2}$ equations for $N^{2}+N$ unknowns (the decomposition is therefore not unique); the system can be solved using Crout's method.

Given a matrix equation and the LU decomposition of A

\mathbf {A} \mathbf {x} =(\mathbf {L} \mathbf {U} )\mathbf {x} =\mathbf {L} (\mathbf {U} \mathbf {x} )=\mathbf {b} ,

first solve $\mathbf {L} \mathbf {y} =\mathbf {b}$ for $\mathbf {y}$ . This can be done by forward substitution

y_{1}={\frac {b_{1}}{l_{11}}}

y_{i}={\frac {1}{l_{ii}}}\left(b_{i}-\sum _{j=1}^{i-1}l_{ij}y_{j}\right)

for i=2,...,N. Then solve $\mathbf {U} \mathbf {x} =\mathbf {y}$ for $\mathbf {x}$ . This can be done by back substitution

x_{N}={\frac {y_{N}}{u_{NN}}}

x_{i}={\frac {1}{u_{ii}}}\left(y_{i}-\sum _{j=i+1}^{N}u_{ij}x_{j}\right)

for i=N-1,...,1

This algorithm for solving matrix equations takes O(n²) operations (excluding the LU decomposition) and is therefore much faster than using Gauss algorithm. MathMartin 16:01, 18 Jul 2004 (UTC)

I suggest the present LU decomposition article should be divided into an article which gives a general explanation of what is an LU decomposition and another article on the doolittle method for LU decomposition. The general LU decomposition article would work as an "anchor" for all the LU decomposition algorythms' articles, including the doolittle.

--Maciel 14:32, 14 Aug 2004 (UTC)

I did some work on the article recently but I am not really satisfied with the result. If you think you can improve the article go ahead and do it. But I would only split the article if it got too long, because it is easier to keep the material consistent if kept on one page. If the page grows too big you could always do a split later. MathMartin 14:53, 14 Aug 2004 (UTC)

On the LU decomposition page, what I don't like is that a matrix A is denoted A (italic font) when not in a math formula, and \mathbf{A} (bold font) when in a math formula. This makes things inconsistent. How to fix this, one way or another? I would suggest making everything italic, why on earth would one want a bold matrix? --Olegalexandrov 21:26, 5 Dec 2004 (UTC)

R in Existence and Uniqueness

In the Existence and Uniqueness section, what is the matrix R? It looks like it's referring to something defined earlier, that is no longer there. --Ben Kovitz 20:16, 4 September 2005 (UTC)[reply]

R should be U. In some languages, including Dutch and German, LU is sometimes called LR (R for right); this probably confused whoever wrote that sentence (which may well have been me). Now fixed, thanks for catching this. -- Jitse Niesen (talk) 21:29, 4 September 2005 (UTC)[reply]

Pivoting Explained? Order of the calculation?

There is mention of PLU and PLUQ decompositions, but no explanation of their use. Should someone mention that the purpose is to have P'AQ'=LU giving P'AQ' with pivots in Gaussian elimination that lead to the least error? (Those ' are transposes, of course.) 149.43.x.x 11:10, 29 September 2006 (UTC)[reply]

Oops, I replaced PLU by LUP, its transpose, and deleted PLUQ. PLUQ does not have any purpose in exact computation. If someone wants to restore PLUQ, then he/she should mention its purpose for approximate computations. Someone@somewhere, 13 June 2007.

Computational complexity

Intro. to Algorithms by Cormen et al. says LU decomposition is Θ(n³) for (I think) an n×n matrix. Did I get that right. If so, it should appear prominently in this article. —Ben FrantzDale 17:18, 1 February 2007 (UTC)[reply]

That's right. More precisely, Golub & van Loan (reference as in Givens rotation) say that you need 2n³ / 3 operations, neglecting lower-order terms. This remains the case if you're doing partial pivoting, but full pivoting is more expensive. Please do put it in. -- Jitse Niesen (talk) 01:44, 2 February 2007 (UTC)[reply]

I added a bit along these lines. Perhaps not prominently, as you proposed; feel free to edit (as always). -- Jitse Niesen (talk) 07:24, 11 March 2007 (UTC)[reply]

There is room for improvement, to

n^{2.376}

, using matrix multiplication with this complexity. You perform the LUP-algorithm (upper) half by (lower) half instead of row by row. I guess you will find it in Modern computer algebra by Joachim von zur Gathen and Jürgen Gerhard--131.174.22.39 (talk) 07:31, 31 March 2009 (UTC)[reply]

Explicit Form

I added a brief paragraph about the explicit form of the LDU factors. Although it is not practical for computation (like the explicit formula of the inverse), I think it can potentially have some great theoretical value (it does in my own research). Householder is the only reference I've found that makes mention of this. Seems that the existence of the explicit form is not so well known today. I don't have the Householder book with me at the moment. If somebody wants to look it up and write out the actual formulas for the factors, knock yourself out. —Preceding unsigned comment added by 76.120.34.175 (talk) 22:49, 12 April 2008 (UTC)[reply]

Inverting triangular matrices

So I checked this article, triangular matrix, and matrix inverse, but I don't see anything specific describing how to efficiently invert a triangular matrix. I assume there must be a method more efficient than Gaussian elimination. How's it work? Dcoetzee 04:02, 23 January 2009 (UTC)[reply]

The LU decomposition is sometimes called Gaussian elimination. It is not matrix inversion, it is rapidly solvable linear systems, that is it. Try this, write out a triangular matrix times a vector equals another vector and try to solve it--then you will get it. PDBailey (talk) 04:17, 23 January 2009 (UTC)[reply]

I know about forward and backward substitution, that's obvious (and described here). I'm asking if there's a more efficient method for inverting a triangular matrix. This article suggests as much when it says that the LU decomposition is used to help compute inverses. Dcoetzee 21:43, 23 January 2009 (UTC)[reply]

Ah, think about AA^-1=I which you could also write as AX^-1=I and solve for X (that is via n individual vector solves). Trefethen & Bau (pp161) give an operation count of 8/3 m³ and say that, "Gaussian elimination can be used to compute the inverse [...] though it is rarely really necessary to do so." PDBailey (talk) 01:47, 24 January 2009 (UTC)[reply]

Yeah, that method of inversion seems kind of pointless when Gaussian elimination uses only 2/3 m³ operations (ref Gaussian_elimination#Analysis). I guess that raises the question, why would you ever compute A^-1 as U^-1L^-1? For numerical stability? If there is no practical reason to do so it might make sense to remove it from Applications. Dcoetzee 02:36, 24 January 2009 (UTC)[reply]

(de-indent) The 2/3 m³ operations is the cost for solving linear systems of m equations with m unknowns (Ax = b). It is the cost for Gaussian elimination on the augmented matrix [A b] of size (m+1) × m. To compute the inverse, you need to do Gaussian elimination on the matrix [A I] which has size (2m) × m. I don't know the cost for this, but I assume it is the same 8/3 m³ as PDBailey quotes. Trefethen and Bau say on the next page that a slightly smarter method, making use of the sparsity of I, costs only 2m³ operations.

I think forward and backward substitution is the fastest method for computing the inverse of a triangular matrix. I think the cost is 1/2 m³ (I find counting operations surprisingly hard so please be suitably sceptical about this). So, the cost of computing A^-1 as U^-1L^-1 is 2/3 m³ (for computing L and U) + 1/2 m³ (for inverting L) + 1/2 m³ (for inverting U) + 1/2 m³ (for multiplying U^-1 with L^-1) = 13/6 m³, so slightly more. I also don't think this method is more stable; in fact, I'd expect the contrary. So indeed, I see no reason to compute the matrix inverse in this way.

The article on Gauss-Jordan elimination says that "it is well suited for calculating the matrix inverse". Embarassingly. when checking the history I found out that I actually added that myself, but I have no idea now why I did that. I don't even remember whether I wanted to say that Gaussian elimination and Gauss-Jordan are equally well suited for calculating the matrix inverse or whether Gauss-Jordan is better than Gaussian elimination. -- Jitse Niesen (talk) 15:29, 24 January 2009 (UTC)[reply]

I just found "Stability of Methods for Matrix Inversion" (PDF), doi:10.1093/imanum/12.1.1, which is probably a good source for information. I haven't read it yet. -- Jitse Niesen (talk) 15:44, 24 January 2009 (UTC)[reply]

In response to why would you do this, It is typical to calculate Ax=b using OLS and then want to know about the variance covariance matrix, which is the inverse of A^TA. So you might first solve for b and then later (having done the LU) invert A^TA. PDBailey (talk) 15:49, 24 January 2009 (UTC)[reply]

LU factorization? Huh?

The article does not define what LU factorization is. —Preceding unsigned comment added by 84.50.107.182 (talk) 09:39, 20 October 2009 (UTC)[reply]

I was wondering the same thing. Is factorization == decomposition? 68.174.98.250 (talk) 18:48, 1 March 2010 (UTC)[reply]

10 years later: Yes. --Svennik (talk) 19:59, 10 December 2020 (UTC)[reply]

Shortcut method for LU decomposition of a square matrix

Here's a link which i found very helpful for decomposition of a square matrix in LU form.

LU Decomposition Video Shortcut method for LU decomposition. — Preceding unsigned comment added by J.Revati (talk • contribs) 13:47, 17 February 2013 (UTC)[reply]

Every invertible matrix A—square or not—?

That shows up in the "Existence and uniqueness" section. Aren't all invertible matrices square? If so, what does this mean? Maybe it should be written in a clearer way. Sinan Taifour (talk) 08:43, 13 November 2010 (UTC)[reply]

Coefficient domain?

This article says nothing about where the coefficients of the matrix are supposed to live. A little reflection tells me there is little hope of finding a LU decomposition for genral matrices with entries in a commutative ring (say polynomials) without extending the coefficient domain. So I think this article is assuming throughout that coefficients are in a field, but it should clearly say so. This also impacts efficiency considerations: if you're dealing with a matrix with integer entries, the computational cost rises notably if you need to compute over the rational numbers (quite likely with bignum numerator/denominator). Even more so for the case of polynomials/rational functions. Marc van Leeuwen (talk) 08:04, 16 May 2011 (UTC)[reply]

fortran example

i wrote a short fortran program to implement LU with pivoting for 2x2 matrices, built from the example in the article. i am happy to submit it here if someone wants to insert it into the main article. 128.83.67.213 (talk) 23:44, 21 December 2011 (UTC)[reply]

The above got a positive response (I think), so here is the code. I don't know how to format code in Wikipedia, or if putting this on a new page of algorithms is better than adding it to the "definition" page. Feel free to use, modify, whatever to your heart's content. The code standard is Fortran 2011 because of the "intent"s and "parameter"s. Take them out and it's Fortran 77. (But don't quote me on that.) 128.83.67.99 (talk) 08:29, 12 July 2015 (UTC)[reply]

       function solve2x2system(K,b) result(x)
       ! for fun, LU decomposition with pivoting
       implicit none

       real*8, intent(in) :: K(2,2)
       real*8, intent(in) :: b(2)
       real*8             :: x(2)

       real*8, parameter :: l11=1d0  ! L as a unit one matrix
       real*8            :: l21      ! 
       real*8, parameter :: l22=1d0  ! 
       real*8 :: u11,u12,u22         ! U
       real*8 :: y1,y2

       ! for pivoting
       real*8, parameter :: small = 1d-14 ! this is what we consider small
       integer p(2)                       ! pivot vector

       ! BEGIN code
       p(1) = 1
       p(2) = 2

       ! check to see if we need to pivot
       if (K(1,1).lt.small) then
         if (abs(K(2,1)).gt.abs(K(1,1))) then
           !write(*,*) "we are pivoting"
           p(1) = 2
           p(2) = 1
         else
           ! don't pivot; K21 is smaller!
         endif
       endif

       ! compute LU decomposition of K
       ! (assuming L's main diagonal is all 1)
       u11 =  K(p(1),1)
       u12 =  K(p(1),2)
       l21 =  K(p(2),1)/u11
       u22 = (K(p(2),2) - l21*u12)

       !write(*,*) "K(p(1),:) = ", K(p(1),1), K(p(1),2)
       !write(*,*) "K(p(2),:) = ", K(p(2),1), K(p(2),2)
       !write(*,*) "p(:)      = ", p
       !write(*,*) "l21 = ", l21
       !write(*,*) "u11 = ", u11
       !write(*,*) "u12 = ", u12
       !write(*,*) "u22 = ", u22

       ! solve Ly=b
       y1 = b(p(1))
       y2 = b(p(2)) - y1*l21
       !write(*,*) "y1,2 = ", y1, y2

       ! solve Ux=y
       x(2) = y2/u22
       x(1) = (y1 - u12*x(2))/u11
       !write(*,*) "x1,2 = ", x(1), x(2)

       return
       end function solve2x2system

Missing word?

A, but these decompositions can all be generalized to rectangular matrices as well.

This sounds strange to me. Is there something missing after "A,"? --MartinThoma (talk) 09:08, 20 May 2013 (UTC)[reply]

You're right, some text was inadvertently deleted a few months ago, and it seems nobody who noticed bothered to fix it. Thanks very much for your comment. It should be better now. -- Jitse Niesen (talk) 00:43, 22 May 2013 (UTC)[reply]

ratio of the principal submatrix?

The article says:

$D_{i}$ is the ratio of the $i^{th}$ principal submatrix to the $(i-1)^{th}$ principal submatrix.

What's the ratio of 2 matrices? Shouldn't this be the ratio of the principal minors?

Bogus explanation

Resolved

– I've replaced this algorithm by the actual one by Cormen et al. Felixytliu (talk) 03:59, 2 June 2020 (UTC)[reply]

Current, Cormen's supposed recursive algorithm for LUP is described as

   Crout and LUP algorithms
   
   The LUP decomposition algorithm by Cormen et al. generalizes Crout matrix decomposition. It can be described as follows.
   
   If A {\displaystyle A} A has a nonzero entry in its first row, then take a permutation matrix P 1 {\displaystyle P_{1}} P_{1} such that A P 1 {\displaystyle AP_{1}} A P_1 has a nonzero entry in its upper left corner. Otherwise, take for P 1 {\displaystyle P_{1}} P_{1} the identity matrix. Let A 1 = A P 1 {\displaystyle A_{1}=AP_{1}} A_1 = A P_1.
   Let A 2 {\displaystyle A_{2}} A_{2} be the matrix that one gets from A 1 {\displaystyle A_{1}} A_{1} by deleting both the first row and the first column. Decompose A 2 = L 2 U 2 P 2 {\displaystyle A_{2}=L_{2}U_{2}P_{2}} A_2 = L_2 U_2 P_2 recursively. Make L {\displaystyle L} L from L 2 {\displaystyle L_{2}} L_{2} by first adding a zero row above and then adding the first column of A 1 {\displaystyle A_{1}} A_{1} at the left.
   Make U 3 {\displaystyle U_{3}} U_3 from U 2 {\displaystyle U_{2}} U_{2} by first adding a zero row above and a zero column at the left and then replacing the upper left entry (which is 0 at this point) by 1. Make P 3 {\displaystyle P_{3}} P_{3} from P 2 {\displaystyle P_{2}} P_{2} in a similar manner and define A 3 = A 1 / P 3 = A P 1 / P 3 {\displaystyle A_{3}=A_{1}/P_{3}=AP_{1}/P_{3}} A_3 = A_1 / P_3 = A P_1 / P_3. Let P {\displaystyle P} P be the inverse of P 1 / P 3 {\displaystyle P_{1}/P_{3}} P_1 / P_3.
   At this point, A 3 {\displaystyle A_{3}} A_{3} is the same as L U 3 {\displaystyle LU_{3}} L U_3, except (possibly) at the first row. If the first row of A {\displaystyle A} A is zero, then A 3 = L U 3 {\displaystyle A_{3}=LU_{3}} A_3 = L U_3, since both have first row zero, and A = L U 3 P {\displaystyle A=LU_{3}P} A = L U_3 P follows, as desired. Otherwise, A 3 {\displaystyle A_{3}} A_{3} and L U 3 {\displaystyle LU_{3}} L U_3 have the same nonzero entry in the upper left corner, and A 3 = L U 3 U 1 {\displaystyle A_{3}=LU_{3}U_{1}} A_3 = L U_3 U_1 for some upper triangular square matrix U 1 {\displaystyle U_{1}} U_{1} with ones on the diagonal ( U 1 {\displaystyle U_{1}} U_{1} clears entries of L U 3 {\displaystyle LU_{3}} L U_3 and adds entries of A 3 {\displaystyle A_{3}} A_{3} by way of the upper left corner). Now A = L U 3 U 1 P {\displaystyle A=LU_{3}U_{1}P} A = L U_3 U_1 P is a decomposition of the desired form.

It is easy to demonstrate that this doesn't work. There is no explanation of how the matrix U1 is to be calculated, and it fact it is easy to demonstrate with basic counter-examples that this supposed U1 matrix cannot be determined. Also, this explanation doesn't seem to be the one in Cormen's well-known textbook. Is there are source for this dodgy hand-waving explanation of the recursive algorithm ? Lathamibird (talk) 10:51, 30 September 2016 (UTC)[reply]

Apparently this algorithm was the first contribution of a netizen way back in 2007, including a reference to CRLS's Intro to Algorithms [2]. It was never revised since then. As you said, this algorithm isn't the one that appears in Cormen's textbook; and the hand-wave in the last step is clearly wrong, since we can't produce a U1 matrix such as described. I'm still not sure how to correct this (I've just stumbled into this situation), but anyone should feel free to intervene. Felixytliu (talk) 04:50, 1 June 2020 (UTC)[reply]

Fundamental omission

This article discusses two things: 1) The LU matrix factorization, and 2) algorithms for calculating the LU factorization.

Most of the article completely fails to distinguish between these two subjects. (For example, the phrase "the LU factorization with pivots": Is this a kind of factorization or just a kind of algorithm for calculating the LU factorization? I think I know the answer, but these things need to be made much clearer".)108.245.209.39 (talk) 20:55, 7 June 2018 (UTC)[reply]

Code examples don't belong in a section about applications

The code examples should be placed in a separate section. --Svennik (talk) 02:14, 9 December 2020 (UTC)[reply]

Done. --Svennik (talk) 19:55, 10 December 2020 (UTC)[reply]

Added Fortran example above to the article [since removed]

Added the Fortran example by 128.83.67.99. Can someone check if it's right? I don't know Fortran. --Svennik (talk) 17:04, 11 December 2020 (UTC)[reply]

From the looks of it, it's only for 2x2 systems. Should it be removed? --Svennik (talk) 09:14, 12 December 2020 (UTC)[reply]

Done. --Svennik (talk) 10:05, 12 December 2020 (UTC)[reply]

Forward v. Backward Looking + Error Analysis

I don't quite know how to fit it in, but would it be worthwhile to distinguish between forward and backward looking LU decomposition? Also, to talk about BLAS-3 v. BLAS-2 v. BLAS-1 algorithms, and error analysis? Fephisto (talk) 05:25, 7 November 2022 (UTC)[reply]

missing subsection?

The page "doolitle algorithm" redirects here but the named subsection "#doolittle..." no longer exists. I'm open to suggestions as to how to properly supply info on said algorithm . Cellocgw (talk) 22:12, 10 December 2023 (UTC)[reply]