CODEWITHSUNDEEP
c++boostmatrix
BartoszKP
BartoszKP

Reputation: 35901

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?

Upvotes: 6

Views: 3128

Answers (2)

Cloud Cho
Cloud Cho

Reputation: 1774

The boost doesn't have document of LU factorization (a lower triangular matrix L and upper triangular matrix U), but the source code shared with the public.

If the code is hard to follow, please check the webpage by Nick Higham. It had an detailed explanation. Here are an example from the link:

Let's say we need to solve Ax = b.
  (1) Make LU from input matrix, A
[3 -1 1  1]
[-1  3 1 -1] ->
[-1 -1 3  1]
[1  1 1  3]

Low
[1     0    0    0]
[-1/3   1   0   0]
[-1/3 -1/2 1 0]
[1/3    1/2  0 1]

Upper
[3    -1   1   1]
[0 8/3 4/3 -2/3]
[0   0   4    1]
[0   0   0    3]

   This example looks straight forward to human but algorithm wise could be numerous steps. This is why LU Factorization came. Methodically, Relation with Gaussian Elimination, Schur Complements, and Block Implementations are some.

  (2) Solve the triangular systems Ly = b and Ux = y, since then b = L(Ux).

Upvotes: 0

BartoszKP
BartoszKP

Reputation: 35901

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.

Upvotes: 8

Related Questions