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
)

Essays/LU Decomposition

See also