【矩阵分解】Python下基于Numpy的四种矩阵基本分解的实现

0x00 需求

完成课堂上讲的关于矩阵分解的

· LU、

· QR（Gram-Schmidt）

· Orthogonal Reduction

Householder reduction

Givens reduction

程序实现，要求一个综合程序，根据选择参数的不同，实现不同的矩阵分解。

反正也是要写，就顺手做成了实现类，可以import调用的那种，为了写作业方便，也设置了输出中间过程，方便拿过程分。

0x01 作者

陈点 学号 201618013229031

个人邮箱 chendian@baidu.com

个人主页 blog.csdn.net/okcd00

Github github.com/okcd00

0x02 文件树

E:\UCAS\矩阵分析与应用\CD_MatrixDecomp 的目录

2016/11/28  23:11                27 A.txt       % 用于测试各种不同种类输入的文件
2016/11/28  23:11                28 B.txt       % 用于测试非方阵进行LU分解时的文件
2016/11/23  15:51                40 LU.txt      % 用于测试LU分解的输入示例
2016/11/28  23:17            10,560 MatrixDecomp.py % 主程序，python MatrixDecomp.py即可
2016/11/23  15:50                27 QR.txt      % 用于测试QR分解的输入示例
2016/11/27  16:01             2,606 Readme.txt      % 说明文档

0x03 程序输入

封装为类，返回值实际为一个dict，在LU分解时返回的dict包含三个键值对：P,A,T，类似的QR分解时范围的是Q与R，value均为np.ndarray格式矩阵

在import该类进行矩阵计算的时候可以直接调用各类函数，将该类为主进程打开时会通过简单的main函数的提示信息引导展示输入输出

E:\UCAS\矩阵分析与应用\BigHomework>python MatrixDecomp.py
> Current Selection is: <Default>
> Please show me the Matrix for Decomposition
> It can be a list or path to a Matrix_File
> Example: [[1,0],[0,1]] or "A.txt"
The Matrix is:

可直接输入形似[[1,0,0],[0,1,0],[0,0,1]]的list，也可以输入文件地址，需要注意地址应该加上引号，如”A.txt”

目录下提供了课件上的LU分解的实例与QR分解时使用的矩阵

> Current Selection is: <Mode>
> Please Select Decomposition Type
> Example: LU GS HH or GV
The Matrix's Size:

此处可选择LU GS HH or GV，分别代表LU分解，Gram-Schmidt、Householder、与Givens的QR分解

0x04 参数

mdp.Show_Process = False

（强烈推荐调成True跑一下看看效果，也不枉coding那么久QvQ）

该参数控制是否输出形似如下形式的中间计算过程：

Calculation[0]:

Q = xxx

R = xxx

Calculation[1]:

Q = xxx

R = xxx

Calculation[x]:

mdp.setMatA(mdp.getInput(‘Default’))

目前已经编码的合法参数为default, random, mode，其中random为随机生成矩阵，根据提示输入矩阵的行列数，值为0-9之间随机选择

0x05 程序源码

#coding=utf8
# ========================================================
#   Copyright (C) 2016 All rights reserved.
#
#   filename : MatrixDecomp.py
#   author   : okcd00 / chendian@baidu.com
#   date     : 2016-11-23
#   desc     : Matrix Decomposition
#   homepage : blog.csdn.net/okcd00
# ========================================================

# Basic package
import random
import os,sys,time

# Numpy for Matrix Calculation
import numpy as np
from numpy import *
from numpy import linalg as La

class MatrixDecomp:

Time = 0
Mode = "NULL" # LU QR(Gram-Schmidt\HouseHolder\Given)
MatA = "No Input" # Matrix_A for Calculation
bak_MatA = "temp" # Matrix_A Back_Up
Show_Process = False

def __init__(self):
self.Time = time.localtime(time.time())

def setMatA(self, inp):
# Judge the type of inp, then achieve the Matrix
if isinstance(inp, np.ndarray):
self.MatA = inp
elif isinstance(inp, list):
self.MatA = np.array(inp)
elif isinstance(inp, str):
if os.path.exists(inp):
self.MatA = np.array(self.readFile(inp))
elif os.path.exists(inp + '.txt'):
self.MatA = np.array(self.readFile(inp + '.txt'))
else:
print "Invalid Input"
self.bak_MatA = self.MatA

def MatDecomp(self, inp):
self.Mode = inp
try:
if inp.upper() == "LU": return self.LU_Decomp(self.MatA)
if inp.upper() == "GS": return self.GS_Decomp(self.MatA)
if inp.upper() == "HH": return self.HH_Decomp(self.MatA)
if inp.upper() == "GV": return self.GV_Decomp(self.MatA)
return "Invalid Decomp Type. (LU/GS/HH/GV)"
except Exception,e:
return "Decomposition Error for %s" % str(e)

def Row_Swap(self, mat, ra, rb):
ret = mat
if mat.ndim == 1:
ret[ra], ret[rb] = mat[rb], mat[ra]
if mat.ndim == 2:
ret[[ra, rb],:] = mat[[rb, ra],:]
return ret

def Col_Swap(self, mat, ca, cb):
ret = mat
if mat.ndim == 1:
ret[ca], ret[cb] = mat[cb], mat[ca]
if mat.ndim == 2:
ret[:,[ca, cb]] = mat[:,[cb, ca]]
return ret

def MaxLine(self, colomn, row):
ret = row
for idx in range(row, colomn.__len__()):
if abs(colomn[idx]) > abs(colomn[ret]):
ret = idx
return ret

def LU_Operation(self, A, cur):
(rSize, cSize) = A.shape
for r in range(cur+1, rSize):
A[r][cur] = A[r][cur] / A[cur][cur]
for c in range(cur+1, cSize):
A[r][c] = A[r][c] - A[r][cur] * A[cur][c]
return A

def LU_GetAns(self, P1D, A):
(rSize, cSize) = A.shape
# Need to transform P from 1D to 2D
P = np.zeros([rSize, rSize])
for idx in range(rSize):
P[idx][P1D[idx]-1] = 1
L = np.eye(rSize, cSize)
U = np.zeros([rSize, cSize])
# Split MatrixA into L-Lower & U-Upper
for r in range(rSize):
for c in range(cSize):
if r <= c : U[r][c] = A[r][c]
else : L[r][c] = A[r][c]
return {'P':P, 'L':L, 'U':U}

def LU_Decomp(self, A): # PA = LU
(rSize, cSize) = A.shape
if rSize!=cSize :
print "> LU_Decomp needs a Nonsingular Square Matrix."
print "> Extend Matrix into a Square Matrix filled by zero."
Size = max(rSize, cSize)
Zero = np.zeros([Size,Size])
Zero[:rSize,:cSize] = np.copy(A)
A, (rSize, cSize) = np.copy(Zero), (Size, Size)
print "> Current Matrix_A = \n", A
P = np.arange(rSize) + 1
for r in range(rSize):
# Swap MaxLine(current_Colomn's Max abs_Value) to the top
idxML = self.MaxLine(A[:,r], r)
A = self.Row_Swap(A, idxML, r)
P = self.Col_Swap(P, idxML, r)
A = self.LU_Operation(A, r)
if self.Show_Process:
print 'Calculation[%d]:\nP = ' % r, P, '^T\nA = \n', A
return self.LU_GetAns(P,A)

def GS_Decomp(self, A): # A = QR
(rSize, cSize) = A.shape
Q, R = np.copy(A), np.zeros([rSize, cSize])
for c in range(cSize):
for r in range(c):
if r < c:
# R_rc = qr^T * Ac
R[r][c] = np.dot(np.transpose(Q[:,r]), A[:,c])
Q[:,c] = Q[:,c] - R[r][c] * Q[:,r]
R[c][c] = La.norm(Q[:,c])
Q[:,c] = Q[:,c] / R[c][c]
if self.Show_Process:
print 'Calculation[%d]:\nQ = \n' % c, Q, '\nR = \n', R
return {'Q':Q, 'R':R}

def HH_Decomp(self, A): # A = QR / PA = T
(rSize, cSize) = A.shape
P = np.eye(rSize, cSize)
for c in range(cSize):
MatA, MatU = np.copy(A[c:,c:]), np.copy(A[c:,c]) # MatU = MatA[:,0]
MatU[0] = MatU[0]+La.norm(MatU) if MatU[0]<0 else MatU[0]-La.norm(MatU)
MatU.shape = (1, MatU.shape[0])
MatU = np.transpose(MatU)
MatR = np.eye(MatU.shape[0])
# Transposing a 1-D array returns an unchanged view of the original array / Note for np.transpose
UTU = np.dot(np.transpose(MatU), MatU)
MatR = MatR - 2.0 * (
(np.dot(MatU, np.transpose(MatU)) / UTU) if UTU!=0 else 0
)
MatA = np.dot(MatR, MatA)
R = np.eye(rSize, cSize)
R[c:,c:] = np.copy(MatR)
P = np.dot(R, P)
A[c:,c:] = np.copy(MatA)
if self.Show_Process:
print 'Calculation[%d]:\nR%d = \n' % (c+1,c+1), MatR, '\nR%dA%d = \n' % (c+1,c+1), MatA, '\nCurrent P = \n', P
return {'Q':np.transpose(P), 'R':A }#, 'T':A}

def GV_Rotate(self, A, i, j):
(rSize, cSize) = A.shape
ret = np.eye(rSize, cSize)
upValue = sum(item**2 for item in A[j:i,j])
c = upValue**0.5 / (upValue + A[i][j]**2)**0.5
s = A[i][j] / (upValue + A[i][j]**2)**0.5
ret[i][i], ret[j][j] = c, c
ret[i][j], ret[j][i] = -s, s
return ret

def GV_Decomp(self, A): # A = QR
(rSize, cSize) = A.shape
U = np.eye(rSize, cSize)
for c in range(cSize):
for r in range(c+1, rSize):
if A[r,c] != 0:
rot = self.GV_Rotate(A,r,c)
U = np.dot(rot, U)
A = np.dot(rot, A)
if self.Show_Process:
print 'Calculation[%d,%d]:\nU%d%d = \n' % (r+1,c+1,r+1,c+1,), rot, '\nCurrent U = \n', U, '\nCurrent A = \n', A
return {'Q':np.transpose(U), 'R':A }

def readFile(self, filename):
# The Matrix File needs to be splited
with open(filename,'r') as f:
ret = []
lines = [ line for line in f.readlines() ]
for each in lines :
line = [ float(num) for num in each.split() ]
ret.append(line)
return ret

def getInput(self, inp='Default'):
print "> Current Selection is: <%s>" % inp
# Get varName or filePath from console for Original Matrix
if inp.upper() == 'DEFAULT':
# print "MatrixDecomp 1.0.7 (v1.0.7, Nov 27 2016, 22:54:40) Type \"help\" for more information."
print "> Please show me the Matrix for Decomposition"
print "> It can be a list or path to a Matrix_File"
print "> Example: [[1,0],[0,1]] or \"A.txt\", \"LU\" etc."
ret = input("The Matrix is: ")
elif inp.upper() == 'RANDOM':
# Define constants
print "> Please show me the Matrix's Size, split by \',\' "
print "> Example: 5,3 or 7,7"
sz = raw_input("The Matrix's Size: ").split(',')
r,c = int(sz[0]), int(sz[1])
# Generate a target Matrix
ret = random.randint(0,9, size=(r,c))
elif inp.upper() == 'MODE':
# Define constants
print "> Please Select Decomposition Type"
print "> Example: LU GS HH or GV"
ret = raw_input("Type of the Matrix Decomposition is: ")
elif inp.upper() == 'HELP':
print """
> Help v1.0.0 Authured by Chendian / okcd00

> mdp.Show_Process
> 该参数控制是否输出中间计算过程, 默认为False, 可在Main函数中改为True

> mdp.setMatA(mdp.getInput('xxx'))
> 目前已经编码的合法参数为default, random, mode, help
"""
mdp.setMatA(mdp.getInput('Default'))
# Can also be used for passing vars,
# i.e. list or numpyArray
else: ret = inp
return ret

if __name__ == "__main__":
np.set_printoptions(suppress=True)
mdp = MatrixDecomp()
mdp.Show_Process = False # True, show mid-calculation
mdp.setMatA(mdp.getInput('Default'))
# mdp.setMatA(mdp.getInput('RANDOM'))
print mdp.MatA
Ans = mdp.MatDecomp(mdp.getInput('Mode'))
try:
print '==========Answer Sheet=========='
for (k,v) in Ans.items():
print '> Matrix', k, '=\n', v
except Exception,e:
print e, '\n', Ans

"""
E:\UCAS\矩阵分析与应用\BigHomework>python MatrixDecomp.py
> Current Selection is: <Default>
> Please show me the Matrix for Decomposition
> It can be a list or path to a Matrix_File
> Example: [[1,0],[0,1]] or "A.txt", "LU" etc.
The Matrix is: "A.txt"
[[  0. -20. -14.]
[  3.  27.  -4.]
[  4.  11.  -2.]]
> Current Selection is: <Mode>
> Please Select Decomposition Type
> Example: LU GS HH or GV
Type of the Matrix Decomposition is: GS
==========Answer Sheet==========
> Matrix Q =
[[ 0.   -0.8  -0.6 ]
[ 0.6   0.48 -0.64]
[ 0.8  -0.36  0.48]]
> Matrix R =
[[  5.  25.  -4.]
[  0.  25.  10.]
[  0.   0.  10.]]
"""

0xFF 测试数据

LU.txt

1 2 -3 4

4 8 12 -8

2 3 2 1

-3 -1 1 -4

QR.txt

0 -20 -14

3 27 -4

4 11 -2

LU.txt for non-square

1 2 -3 4

4 8 12 -8

2 3 2 1