How to understand the pivot matrix of scipy.linalg.lu_factor?

Question

How can I manually reconstruct a matrix A that was factorized by lu_factor? (A = PLU)

My current attempts all failed due to the setup of matrix P. Here is what I have so far:

A = np.random.rand(3,3)
lu, piv = lu_factor(A)

U = np.triu(lu)
L = np.tril(lu, -1)
L[np.diag_indices_from(L)] = 1.0

I am looking for the matrix P that makes this line print True:

print np.allclose(A, np.dot(P, np.dot(L, U)))

Any hint/link/suggestion is appreciated!

Answer 1

The permutation vector needs to be interpreted in sequence. If piv=[1,2,2] then the following needs to be done in sequence (with zero-based indexing):

1. Row 0 changes with Row 1
2. The new Row 1 changes with Row 2 and
3. The new Row 2 stays the same.

In code this would do the trick:

P = np.eye(3)
for i, p in enumerate(piv):
    Q = np.eye(3,3)
    q = Q[i,:].copy()
    Q[i,:] = Q[p,:]
    Q[p,:] = q
    P = np.dot(P, Q)

For piv=[1,2,2] P is

[[ 0.  0.  1.]
 [ 1.  0.  0.]
 [ 0.  1.  0.]]

This is probably not a very fast way of computing P but it does the trick and answers the question.