Essays/LU Decomposition
Given a matrix A where <:/$A , the LU decomposition produces a permutation (vector) p , a unit lower triangular matrix L , and an upper triangular matrix U , such that
mp=: +/ . * NB. matrix product A -: p {"1 L mp U
The permutation p is necessary as otherwise some non-singular matrices such as 2 2$0 1 1 0 would not have a LU decomposition. The permutation can be converted into the more conventional permutation matrix P=:(i.#p)=/p , whence the equivalence becomes A -: L mp U mp P .
A recursive computations of the LU decomposition reveals itself by considering A as a 2-row matrix of matrices and substituting matrix operations for scalar and vector ones in the LU decomposition of a 2-row matrix of scalars. The algorithm is described in section 6.4 of Aho, Hopcroft, and Ullman, The Design and Analysis of Computer Algorithms, Addison-Wesley, 1974.
LU=: 3 : 0 'm n'=. $ A=. y if. 1=m do. p ; (=1) ; p{"1 A [ p=. C. (n-1);~.0,(0~:,A)i.1 else. m2=. >.m%2 'p1 L1 U1'=. LU m2{.A D=. (/:p1) {"1 m2}.A F=. m2 {."1 D E=. m2 {."1 U1 FE1=. F mp %. E G=. m2}."1 D - FE1 mp U1 'p2 L2 U2'=. LU G p3=. (i.m2),m2+p2 H=. (/:p3) {"1 U1 (p1{p3) ; (L1,FE1,.L2) ; H,(-n){."1 U2 end. )
The following table summarizes the local names in LU , with m2=.>.m%2
name type shape A matrix m,n p1 permutation n L1 L m2,m2 U1 U m2,n D matrix (m-m2),n F matrix (m-m2),m2 E U m2,m2 FE1 matrix (m-m2),m2 G matrix (m-m2),n-m2 p2 permutation n-m2 L2 L (m-m2),m-m2 U2 U (m-m2),n-m2 p3 permutation n H U m2,n
The verb is now applied to an example:
] A=: (6 8 ?.@$ 20x) * 0=6 8 ?.@$ 3 6 0 0 0 0 19 0 0 0 0 6 0 0 0 0 0 0 0 0 2 0 0 0 4 4 0 0 0 16 0 0 0 0 8 2 0 0 0 19 0 1 0 0 0 17 0 0 13 LU A ┌───────────────┬───────────────────┬───────────────────────┐ │0 4 1 2 3 5 6 7│ 1 0 0 0 0 0│6 0 0 0 0 19 0 0│ │ │ 0 1 0 0 0 0│0 6 0 0 0 0 0 0│ │ │ 0 0 1 0 0 0│0 0 2 0 0 0 0 4│ │ │2r3 0 0 1 0 0│0 0 0 16 0 _38r3 0 0│ │ │ 0 1r3 0 0 1 0│0 0 0 0 8 0 19 0│ │ │1r6 0 0 17r16 0 1│0 0 0 0 0 247r24 0 13│ └───────────────┴───────────────────┴───────────────────────┘ 'p L U'=: LU A A -: p {"1 L mp U 1 ] P=: (i.#p)=/p 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 A -: L mp U mp P 1
LUcheck has the same argument and result as LU , but incorporates checks:
LUcheck=: 3 : 0 A=. y assert. <:/$A 'p L U'=.z=. LU A 'm n'=. $A assert. (($p) -: ,n) *. (i.n) e. p assert. (($L) -: m,m) *. ((0~:L)<:(i.m)>:/i.m) *. 1=(<0 1)|:L assert. (($U) -: m,n) *. (0~:U)<:(i.m)<:/i.n assert. A -: p {"1 L mp U z )
See also
- Cholesky Decomposition
- LU Decomposition
- QR Decomposition
- Matrix Inverse
- Triangular Matrix Inverse
- Determinant
- Minors
- Hilbert Matrix
- Block Matrix Inverse
- Kronecker Product
Contributed by Roger Hui.