I am trying to get the LDLt factorization of a given symmetric matrix with SciPy's Python bindings to LAPACK using the dsysv routine which actually solves linear systems using this matrix factorization.
I have tried the following:
import numpy as np
from scipy.linalg.lapack import dsysv
A = np.random.randint(1, 1000, size=(5, 5))
A = (A + A.T)
b = np.random.randn(5)
lult, piv, x, _ = dsysv(A, b, lower=1)
Where x would be the solution for the above linear system and lult and piv contain information about the factorization.
How can I reconstruct LDLt from it? Sometimes negative values are contained in piv and from the docs I was not able to understand their meaning.
LAPACK's sytrf actually computes this factorization (without solving any linear system) but it does not seem available via SciPy.
There is an example here with the output I am interested in (see eq. 3-23).
All the required information is found in the documentation of systrf. But admittedly, it is a somewhat verbose.
So just give me the code:
import numpy as np
from scipy.linalg.lapack import dsysv
def swapped(i, k, n):
"""identity matrix where ith row and column are swappend with kth row and column"""
P = np.eye(n)
P[i, i] = 0
P[k, k] = 0
P[i, k] = 1
P[k, i] = 1
return P
# example
n = 5
A = np.random.rand(n, n)
A = (A + A.T)
b = np.random.randn(n)
lult, piv, x, _ = dsysv(A, b, lower=1)
# reconstruct L and D
D = np.zeros_like(A, dtype=float)
L = np.eye(n)
k = 0
while k < n:
i = piv[k]
if i < 0:
s = 2
else:
s = 1
if s == 1:
i = i - 1
D[k, k] = lult[k, k] # D(k) overwrites A(k,k)
Pk = swapped(k, i, n)
v = lult[k+1:n, k] # v overwrites A(k+1:n,k)
Lk = np.eye(n)
Lk[k+1:n, k] = v
else:
m = -i - 1
D[k:k+2, k:k+2] = lult[k:k+2, k:k+2] # the lower triangle of D(k) overwrites A(k,k), A(k+1,k), and A(k+1,k+1)
D[k, k+1] = D[k+1, k] # D is symmeric
Pk = swapped(k+1, m, n)
v = lult[k+2:n, k:k+2] # v overwrites A(k+2:n,k:k+1)
Lk = np.eye(n)
Lk[k+2:n, k:k+2] = v
L = L.dot(Pk).dot(Lk)
if s == 1:
k += 1
else:
k += 2
print(np.max(np.abs(A - L.dot(D).dot(L.T)))) # should be close to 0
The snipped above reconstructs L and D from the decomposition (it would need to be adapted to reconstruct U from an UDUt decomposition). I will try to explain below. First a quote from the documentation:
... additional row interchanges are required to recover U or L explicitly (which is seldom necessary).
Reconstructing L (or U) requires a number of iterations with row exchanging operations and matrix multiplication. This is not very efficient (less so when done in Python) but luckily this reconstruction is seldom necessary. So make sure you really have to do this!
We reconstruct L from L = P(1)*L(1)* ... *P(k)*L(k)*...,. (Fortran indices are 1-based). So we need to iterate k from 0 to n, obtain K and L in each step and multiply them.
P is a permutation matrix, defined by piv. A positive value of piv is straight-forward (i = piv[k]). It means that the ith and kth row/column were swapped in A before performing the operation. In this case the kth diagonal element of lult corresponds to the kth diagonal element of D. L(k) contains the kth column of the lower diagonal matrix - after the swapping.
A negative value of piv means that the corresponding element of D is a 2x2 block instead of just one element, and L(k) corresponds to two columns of the lower diagonal matrix.
Now for each step in k we obtain L(k), apply the swapping operation P(k), and combine it with the existing L. We also obtain the 1x1 or 2x2 block of D and correspondingly increase k by 1 or 2 for the next step.
I won't blame anyone for not comprehending my explanation. I simply wrote it down as I figured it out... Hopefully, the combination of the code snippet, the description, and the original documentation prove useful :)