Source code for connectivipy.mvar.comp
# -*- coding: utf-8 -*-
#!/usr/bin/env python
from __future__ import absolute_import
from __future__ import print_function
import numpy as np
from six.moves import range
[docs]def ldl(A):
"""
LDL decomposition (implementation from *en.wikipedia.org/wiki/Cholesky_decomposition*)
Args:
*A* : numpy.array
matrix kXk
Returns:
*L*, *D*, *LT* : np.array
*L* is a lower unit triangular matrix, *D* is a diagonal matrix
and *LT* is a transpose of *L*.
"""
n = A.shape[1]
L = np.eye(n)
D = np.zeros(n)
for j in range(n):
D[j] = A[j, j] - np.dot(L[j, :j]**2, D[:j])
for i in range(j+1, n):
L[i, j] = (A[i, j] - np.dot(L[i, :j]*L[j, :j], D[:j]))/D[j]
D = np.diag(D)
return L, D, L.T
if __name__ == '__main__':
# test of wikipedia data
A = np.array([[4, 12, -16], [12, 37, -43], [-16, -43, 98]], dtype=float)
l, d, lt = ldl(A)
print(l)
print(d)
print(lt)