LU decomposition with 5 output arguments

Dear all,

I am trying to solve a huge linear system Ax = b using LU decomposition with five output arguments:

[L, U, p, q, D] = lu(A, 'vector');

My matrix A is sparse with nnz(A) =1e+9 and size(A) = (1e+7, 1e+7)

However, I cannot find any official help to use these 5 arguments to solve the system.

Does anyone have any idea how to ultilize these arguments?

Best,

Doan

Answer 1

[Christine Tobler]

As Matt says, the simplest is likely to just use decomposition object.

In terms of how to solve a linear system, the doc page gives the formula

[L,U,P,Q,D] = lu(S) satisfies P*(D\S)*Q = L*U

from this we can get a formula for the matrix S by just reformulating:

S = D*P'*L*U*Q'

So you can write the initial formula to solve

S*x = b

and insert that formula

(D*P'*L*U*Q') * x = b

now applying each factor's inverse from the left results in

x = Q*(U\(L\(P*(D\x))));

To use the 'vector' syntax instead, each of P*... and Q*... will become an indexing operation instead.

[Duy Hai]

Thank a lot. That is the answer I am looking for.

Small typo, but I am sure @Christine Tobler really means

x = Q*(U\(L\(P*(D\b))));

[Torsten]

@Christine Tobler

And how to do it when the "vector" option is used ? Do the vectors P and Q have to be rewritten as matrices, or can the solve be done directly in the vector format ?

[Duy Hai]

When the "vector" option is used, the solving process should be

c = D\b;
x(q) = U\(L\c(p));

[Torsten]

Thanks. Except that x is a row vector instead of a column vector, this seems to work.

Answer 2

[Matt J]

dA=decomposition(A,'ldl');

x=dA\b;

[Duy Hai]

Thanks, but my matrix is not Hermitian so that LDL doesn't apply.

decomposition function is quite handy but I just want to save memory space using vector permutations p and q instead of matrices. I have no idea on which option that decomposition uses: 'vector' format or the default format 'matrix'.

[Torsten]

[L,U,P,Q,D] = lu(S) also returns a diagonal scaling matrix D such that P*(D\S)*Q = L*U. Typically, the row-scaling leads to a sparser and more stable factorization.

outputForm — Shape of permutation outputs

'matrix' (default) | 'vector'

Shape of permutation outputs, specified as 'matrix' or 'vector'. This flag controls whether lu returns the row permutations P and column permutations Q as permutation matrices or permutation vectors.

As matrices, the outputs P and Q satisfy these identities:

• Three outputs — P satisfies P*A = L*U.
• Four outputs — P and Q satisfy P*S*Q = L*U.
• Five outputs — P, Q, and D satisfy P*(D\S)*Q = L*U.

As vectors, the outputs P and Q satisfy these identities:

• Three outputs — P satisfies A(P,:) = L*U
• Four outputs — P and Q satisfy S(P,Q) = L*U
• Five outputs — P, Q, and D satisfy D(:,P)\S(:,Q) = L*U.

Example: [L,U,P] = lu(A,'vector')

[Matt J]

Thanks, but my matrix is not Hermitian so that LDL doesn't apply.

Well, you can just do,

dA=decomposition(A)

and let it choose what is best for thte matrix

I have no idea on which option that decomposition uses: 'vector' format or the default format 'matrix'.

I don't think you need to know. You can just check how much memory the dA object consumes.