Python——Matrix Decomposition

Python的Numpy库提供了很多矩阵分解的算法，比如QR分解，Cholesky分解，LU分解。我们可以调用Numpy/Scipy的linalg library

import pprint
import scipy
import scipy.linalg   # SciPy Linear Algebra Library

A = scipy.array([[6, 3, 4, 8], [3, 6, 5, 1], [4, 5, 10, 7], [8, 1, 7, 25]])
L = scipy.linalg.cholesky(A, lower=True)
U = scipy.linalg.cholesky(A, lower=False)

print "A:"
pprint.pprint(A)

print "L:"
pprint.pprint(L)

print "U:"
pprint.pprint(U)

import pprint
import scipy
import scipy.linalg   # SciPy Linear Algebra Library

A = scipy.array([ [7, 3, -1, 2], [3, 8, 1, -4], [-1, 1, 4, -1], [2, -4, -1, 6] ])
P, L, U = scipy.linalg.lu(A)  

print "A:"
pprint.pprint(A)

print "P:"
pprint.pprint(P)

print "L:"
pprint.pprint(L)

print "U:"
pprint.pprint(U)

PA=LU，LU分解的含义

In order to solve for the lower triangular matrix, we will make use of the Cholesky-Banachiewicz Algorithm. First, we calculate the values for L on the 

main diagonal. Subsequently, we calculate the off-diagonals for the elements below the diagonal:

lkklik==akk−∑j=1k−1l2kj−−−−−−−−−−⎷1lkk⎛⎝aik−∑j=1k−1lijlkj⎞⎠,i>k

from math import sqrt
from pprint import pprint
 
def cholesky(A):
    """Performs a Cholesky decomposition of A, which must 
    be a symmetric and positive definite matrix. The function
    returns the lower variant triangular matrix, L."""
    n = len(A)

    # Create zero matrix for L
    L = [[0.0] * n for i in xrange(n)]

    # Perform the Cholesky decomposition
    for i in xrange(n):
        for k in xrange(i+1):
            tmp_sum = sum(L[i][j] * L[k][j] for j in xrange(k))
            
            if (i == k): # Diagonal elements
                # LaTeX: l_{kk} = \sqrt{ a_{kk} - \sum^{k-1}_{j=1} l^2_{kj}}
                L[i][k] = sqrt(A[i][i] - tmp_sum)
            else:
                # LaTeX: l_{ik} = \frac{1}{l_{kk}} \left( a_{ik} - \sum^{k-1}_{j=1} l_{ij} l_{kj} \right)
                L[i][k] = (1.0 / L[k][k] * (A[i][k] - tmp_sum))
    return L
 
A = [[6, 3, 4, 8], [3, 6, 5, 1], [4, 5, 10, 7], [8, 1, 7, 25]]
L = cholesky(A)

print "A:"
pprint(A)

print "L:"
pprint(L)

这里要注意几个细节，子循环语句中是xrange(i+1)，表示从0,1,..,i，这样就可以遍历矩阵的下三角，如果是xrange(i)，那么取不到对角线元素。且j一直到k-1也是类似的，for j in xrange(k)