Strictly Triangular LU

Strictly Triangular L.U formulation of Nested Factorization

by John Appleyard B.A., Ph D, Polyhedron Software Ltd, May 2004

 

The nested factorization approximation, P, for a 2D matrix is given by

 

P = (T+m)T^-1(T+v)

where   T = (γ+ l)γ^-1(γ+u)

 

so         P = ((γ+ l)γ^-1(γ+u)+ m) .(γ+ u) ^-1.γ. (γ+l) ^-1.((γ+ l)γ^-1(γ+u)+v)

P = (γ + l + m .(I+ U) ^-1).γ^-1.(γ + u + (I+L) ^-1v)

 

where I is the unit matrix, L= lγ^-1  and  U= γ^-1 u

 

P= (γ + l + m - m.U + m.U² -  m.U³ +....).γ^-1(γ + u + v - L.v ₊ L².v- L³.v+....)

 

So P can be expressed as the product of strictly upper and lower triangular factors.  However, when compared to a zero-fill ILU method, the factors have extra fill in bands of the form  m Uⁿ and Lⁿ v which are exactly like those inserted in more accurate ILU decompositions.  The difference is that nested factorization produces a complete subset of the required fill-ins with little extra computation, and no extra storage.

 

Sparsity Pattern of L factor

 

For 3D matrices, a similar analysis gives

 

B = [γ+l+ m(I+U)^-1+n(I+U+(I+L)^-1V)^-1] γ^-1 [γ+u+(I+L)^-1v +(I+L+M(I+U)^-1)^-1w]

 

where M= mγ^-1  and  V= γ^-1 v.  The sparsity pattern of the L factor in this strictly triangular formulation is shown above.

 

An Alternative Implementation of Nested Factorization

 

It is possible to implement nested factorization using strictly triangular factors, without explicitly computing all the fill-in elements.  For example, in the 2D factorization, the block of m elements which connects to the previous line is replaced by

 

 

 

When evaluating the product (k) of this with a vector x, we can use a recurrence formula to include the extra bands.  Instead of 

 

k_i  = m_nx+i . x_i

 

we work up from the bottom of the using

 

temp  =  -  U_i (x_i+1  +  temp ) 

 k_i  = m_nx+i . (x_{i +}  temp)

 

where   temp is initially 0  

 

For the corresponding v elements, which connect each line to the next line, we use

 

 

 

and compute k by working from the top using the recurrence formula

 

k_i   = v_i x_i+nx – L_i.k­_i-1 

 

where   k₀ is 0