LU/PLU分解

一、简介

在线性代数中， LU分解(LU Decomposition)是矩阵分解的一种，可以将一个矩阵分解为一个单位下三角矩阵和一个上三角矩阵的乘积（有时是它们和一个置换矩阵的乘积）。LU分解主要应用在数值分析中，用来解线性方程、求反矩阵或计算行列式。
详解

二、实现

# -*- coding: utf-8 -*-
"""
Created on Tue Dec 13 15:59:27 2016

    LU PLU矩阵分解 LU求解方程

@author: Administrator
"""

from numpy import *
import matplotlib.pyplot as plt

def lu_decomp(A):
    m,n = A.shape
    if m != n:
        return
    temp_A = A.copy()
    L = eye(n)
    U = zeros((m,n))
    for k in range(0,n-1):
        for i in range(k+1,n):
           if abs(temp_A[i,k]) > 1.0e-9:
               lam = temp_A[i,k] / temp_A[k,k]
               temp_A[i,k+1:n] = temp_A[i,k+1:n] - lam * temp_A[k,k+1:n]
               temp_A[i,k] = 0
               L[i,k] = lam
    U[:n,:n] = temp_A[:n,:n]
    return L,U

def plu_decomp(A):
    m,n = A.shape
    if m != n:
        return
    temp_A = A.copy()
    L = eye(n)
    U = zeros((m,n))
    index = [0 for i in range(n-1)]
    for k in range(0,n-1):
        index[k] = argmax(abs(temp_A[k:n,k])) + k
        if(index[k] != k):
            for j in range(n):
                temp_A[index[k],j],temp_A[k,j] = temp_A[k,j],temp_A[index[k],j]
        for i in range(k+1,n):
           if abs(temp_A[i,k]) > 1.0e-9:
               if temp_A[k,k] == 0.:
                   return
               lam = temp_A[i,k] / temp_A[k,k]
               temp_A[i,k+1:n] = temp_A[i,k+1:n] - lam * temp_A[k,k+1:n]
               temp_A[i,k] = 0
               L[i,k] = lam
    U[:n,:n] = temp_A[:n,:n]
    return U

#列主元求解方程 ex=true
def solve_ex(A,b,ex=True):
    m,n = A.shape
    if m != n:
        return
    temp_A = A.copy()
    x = zeros(n)
    index = [0 for i in range(n-1)]
    #消去
    for k in range(0,n-1):
        index[k] = argmax(abs(temp_A[k:n,k])) + k
        if(index[k] != k and ex):
            for j in range(n):
                temp_A[index[k],j],temp_A[k,j] = temp_A[k,j],temp_A[index[k],j]
            b[k],b[index[k]] = b[index[k]],b[k]
        for i in range(k+1,n):
           if abs(temp_A[i,k]) > 1.0e-9:
               lam = temp_A[i,k] / temp_A[k,k]
               temp_A[i,k+1:n] = temp_A[i,k+1:n] - lam * temp_A[k,k+1:n]
               temp_A[i,k] = 0
               b[i] = b[i] - lam * b[k]
    #回代
    x[n-1] = b[n-1] / temp_A[n-1,n-1]
    for k in range(n-2,-1,-1):
        temp = b[k]
        for j in range(k+1,n):
            temp = temp - temp_A[k,j] * x[j]
        x[k] = temp / temp_A[k,k]
    return x

if __name__ == '__main__':
    A = array([[1,-6,-1],[1,2,2],[2,-2,1]]) #square matrix
    b = array([-2,4,1],dtype=float)
    A = array([[2,1,1,0],[4,3,3,1],[8,7,9,5],[6,7,9,8]]) #square matrix
    A = matrix(A,dtype=float)
    m,n = A.shape
    L,U = lu_decomp(A)
    e = eye(n,dtype=float)
    A_ = []
    for i in range(0,n):
        A_.append(solve_ex(U,solve_ex(L,e[:,i]),ex=False))
    print matrix(A_,dtype=float)