Is there a built-in/easy LDU decomposition method in Numpy?

Scipy has an LU decomposition function: scipy.linalg.lu. Note that this also introduces a permutation matrix P into the mix. This answer gives a nice explanation of why this happens.

If you specifically need LDU, then you can just normalize the U matrix to pull out D.

Here's how you might do it:

>>> import numpy as np>>> import scipy.linalg as la>>> a = np.array([[2, 4, 5],                  [1, 3, 2],                  [4, 2, 1]])>>> (P, L, U) = la.lu(a)>>> Parray([[ 0.,  1.,  0.],       [ 0.,  0.,  1.],       [ 1.,  0.,  0.]])>>> Larray([[ 1.        ,  0.        ,  0.        ],       [ 0.5       ,  1.        ,  0.        ],       [ 0.25      ,  0.83333333,  1.        ]])>>> Uarray([[ 4. ,  2. ,  1. ],       [ 0. ,  3. ,  4.5],       [ 0. ,  0. , -2. ]])>>> D = np.diag(np.diag(U))   # D is just the diagonal of U>>> U /= np.diag(U)[:, None]  # Normalize rows of U>>> P.dot(L.dot(D.dot(U)))    # Checkarray([[ 2.,  4.,  5.],       [ 1.,  3.,  2.],       [ 4.,  2.,  1.]])

Try this:

import numpy as npA = np.array([[1,2,3,4],[1,2,3,4],[1,2,3,4],[1,2,3,4]])U = np.triu(A,1)L = np.tril(A,-1)D = np.tril(np.triu(A))print(A)print(L)print(D)print(U)