I see cholesky decomposition in numpy.linalg.cholesky, but could not find a LDU decompositon. Can anyone suggest a function to use?
-
Take a look here: docs.scipy.org/doc/numpy/reference/generated/…– Dalton CézaneCommented Aug 2, 2017 at 1:51
-
7This 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.– cs95Commented Aug 2, 2017 at 1:57
-
@cᴏʟᴅsᴘᴇᴇᴅ It does fall under "..request for API.."; although that's sort of nebulous.– user2864740Commented 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.– cs95Commented Aug 2, 2017 at 2:04
-
3@user2357112 Just get a username already, or I'll start calling you Twickler.– cs95Commented Aug 2, 2017 at 2:05
|
Show 2 more comments
2 Answers
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.]])
-
It'd be useful to demonstrate how to perform the normalization. Commented Aug 2, 2017 at 2:41
-
-
4This 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. LiCommented Aug 29, 2022 at 23:07
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)
-
3This 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
– PraveenCommented Apr 8, 2022 at 15:15