14

I see cholesky decomposition in numpy.linalg.cholesky, but could not find a LDU decompositon. Can anyone suggest a function to use?

7
  • Take a look here: docs.scipy.org/doc/numpy/reference/generated/… Commented Aug 2, 2017 at 1:51
  • 7
    This is not an off topic request, there is a function in scipy which does this. Whoever voted to close - you don't seem to know that, you probably shouldn't be viewing this tag.
    – cs95
    Commented Aug 2, 2017 at 1:57
  • @cᴏʟᴅsᴘᴇᴇᴅ It does fall under "..request for API.."; although that's sort of nebulous. Commented Aug 2, 2017 at 1:59
  • @user2357112 I assumed pivot and normal were the same since the permutation matrix was indeed being set to the identity. Thanks for correcting me.
    – cs95
    Commented Aug 2, 2017 at 2:04
  • 3
    @user2357112 Just get a username already, or I'll start calling you Twickler.
    – cs95
    Commented Aug 2, 2017 at 2:05

2 Answers 2

20

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)
>>> P
array([[ 0.,  1.,  0.],
       [ 0.,  0.,  1.],
       [ 1.,  0.,  0.]])
>>> L
array([[ 1.        ,  0.        ,  0.        ],
       [ 0.5       ,  1.        ,  0.        ],
       [ 0.25      ,  0.83333333,  1.        ]])
>>> U
array([[ 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)))    # Check
array([[ 2.,  4.,  5.],
       [ 1.,  3.,  2.],
       [ 4.,  2.,  1.]])
5
  • It'd be useful to demonstrate how to perform the normalization. Commented Aug 2, 2017 at 2:41
  • @user2357112 Done.
    – Praveen
    Commented Aug 2, 2017 at 3:01
  • 4
    This looks like the best available built-in, but it's disappointing that it gives a non-identity permutation matrix for an input that looks like it could be LU factorized without one. Commented Aug 2, 2017 at 3:20
  • This only works for invertible matrices. When a is singular in your example U will have 0 entries on the diagonal, and the normalization procedure does not work. Commented Feb 26, 2021 at 9:56
  • I might be missing something, but P @ L just destroys the low-triangle structure?
    – C. Li
    Commented Aug 29, 2022 at 23:07
-1

Try this:

import numpy as np

A = 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)
1
  • 3
    This isn't an LDU decomposition: it's just pulling out the upper and lower triangles and the diagonal of A. For an actual LDU decomposition, when you multiply L, D and U, you should recover A, i.e., A = L @ D @ U
    – Praveen
    Commented Apr 8, 2022 at 15:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.