c++boostmatrix

What does lu_factorize return?


boost::number::ublas contains the M::size_type lu_factorize(M& m) function. Its name suggests that it performs the LU decomposition of a given matrix m, i.e. should produce two matrices that m = L*U. There seems to be no documentation provided for this function.

It is easy to deduce that it returns 0 to indicate successful decomposition, and a non-zero value when the matrix is singular. However, it is completely unclear where is the result. Taking the matrix by reference suggests that it works in-place, however it should produce two matrices (L and U) not one. So what does it do?


Solution

  • There is no documentation in boost, but looking at the documentation of SciPy's lu_factor one can see, that it's not uncommon to return one result for the LU decomposition.

    This is enough, because in a typical approach to LU decomposition, L's diagonal consists of ones only, as presented in this answer from Mathematics, for example.

    So, it is possible to fit both L and U into one matrix, putting L in result's lower part, omitting the diagonal (which is assumed to contain only ones), and U in the upper part. For example, for a 3x3 problem the result is:

        u11 u12 u13
    m = l21 u22 u23
        l31 l32 u33
    

    which implies:

         1    0   0
    L =  l21  1   0
         l31  l32 1
    

    and

        u11 u12 u13
    U = 0   u22 u23
        0   0   u33
    

    Inspecting boost's void lu_substitute(const M& m, vector_expression<E>& e) function, from the same namespace seems to confirm this. It solves the equation LUx = e, where both L and U are contained in its m argument in two steps.

    First solve Lz = e for z, where z = Ux, using lower part of m:

    inplace_solve(m, e, unit_lower_tag ());
    

    then, having computed z = Ux (with e modified in place), Ux = e can be solved, using upper part of m:

    inplace_solve(m, e, upper_tag ());
    

    inplace_solve is mentioned in the documentation, and it:

    Solves a system of linear equations with triangular form, i.e. A is triangular.

    So everything seems to make sense.