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
)
```