【矩阵分解】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 文件树

1
2
3
4
5
6
7
8
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函数的提示信息引导展示输入输出

1
2
3
4
5
6
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分解时使用的矩阵

1
2
3
4
> 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 程序源码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
#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