Essays/Triangular Matrix Inverse

From J Wiki
Jump to: navigation, search

A matrix is triangular if it satisfies one of the following propositions:

upperQ=: 0&= >: >/&i.@$
lowerQ=: 0&= >: </&i.@$

A recursive method for inverting a square triangular matrix reveals itself by considering it as a 2-by-2 matrix of matrices and substituting matrix operations for scalar ones in the inversion of a 2-by-2 triangular matrix of scalars.

uinv=: 3 : 0
 if. 1>:#y do. %y return. end.
 m =. >.-:#y
 n =. <.-:#y
 ai=. uinv (m,m){.y
 di=. uinv (m,m)}.y
 b =. (m,-n){.y
 x =. - ai +/ .* b +/ .* di
 (ai,.x) , (-m+n){."1 di
)

linv=: 3 : 0
 if. 1>:#y do. %y return. end.
 m =. >.-:#y
 n =. <.-:#y
 ai=. linv (m,m){.y
 di=. linv (m,m)}.y
 c =. ((-n),m){.y
 y =. - di +/ .* c +/ .* ai
 ((m+n){."1 ai) , (y,.di)
)

Since uinv -: linv&.|: and linv -: uinv&.|: , one of the verbs can be defined simply in terms of the other.

Some examples using Pascal's triangle and Stirling numbers:

   ] A=: !/~ i.7x
1 1 1 1 1  1  1
0 1 2 3 4  5  6
0 0 1 3 6 10 15
0 0 0 1 4 10 20
0 0 0 0 1  5 15
0 0 0 0 0  1  6
0 0 0 0 0  0  1
   uinv A
1 _1  1 _1  1  _1   1
0  1 _2  3 _4   5  _6
0  0  1 _3  6 _10  15
0  0  0  1 _4  10 _20
0  0  0  0  1  _5  15
0  0  0  0  0   1  _6
0  0  0  0  0   0   1
   A +/ .* uinv A
1 0 0 0 0 0 0
0 1 0 0 0 0 0
0 0 1 0 0 0 0
0 0 0 1 0 0 0
0 0 0 0 1 0 0
0 0 0 0 0 1 0
0 0 0 0 0 0 1

   s1=: 1: ` ([     ((0:,]) + <:@[ * ],0:) $:@<:) @. *
   s2=: 1: ` (i.@>: ((0:,]) + [    * ],0:) $:@<:) @. *

   s1"0 i.7x
1   0   0   0  0  0 0
0   1   0   0  0  0 0
0   1   1   0  0  0 0
0   2   3   1  0  0 0
0   6  11   6  1  0 0
0  24  50  35 10  1 0
0 120 274 225 85 15 1
   s2"0 i.7x
1 0  0  0  0  0 0
0 1  0  0  0  0 0
0 1  1  0  0  0 0
0 1  3  1  0  0 0
0 1  7  6  1  0 0
0 1 15 25 10  1 0
0 1 31 90 65 15 1
   | linv s1"0 i.7x
1 0  0  0  0  0 0
0 1  0  0  0  0 0
0 1  1  0  0  0 0
0 1  3  1  0  0 0
0 1  7  6  1  0 0
0 1 15 25 10  1 0
0 1 31 90 65 15 1



See also


Contributed by Roger Hui.