5. SciPy线性方程组求解
5.1 LU分解
将系数矩阵A转变成等价两个矩阵L和U的乘积 ，其中L和U分别是单位下三角矩阵和上三角矩阵。当A的所有顺序主子式都不为0时，矩阵A可以分解为A=LU（所有顺序主子式不为0，矩阵不一定不可以进行LU分解）。其中L是下三角矩阵，U是上三角矩阵。
当需要求解 Ax=b 的时候，左边的矩阵 A 很多时候是不变的，而右边的 b 随着输入而变化。做 LU 分解时，只会用到矩阵 A ，所以可以预先准备好 L 与 U ，当有求解 b 的需求时，直接拿来用就好了。求解方差 $Ax=LUx=b$分两步(1)解方程 $Ly = b$ ，得到 y，(2)再通过 $Ux = y$ ，得到 x。
5.1.1 SciPy的LU分解
SciPy里的scipy.linalg.lu函数可以基本实现对Ax=b的LU分解，但scipy.linalg.lu函数的返回值有三个p'、l'、u',所以矩阵分解变为(P'L')U' = A。用$P'L'y = b$在用$U'x = y$，差不多。
from scipy.linalg import lu
import numpy as np
A = np.matrix([[2,3],[5,4]])
b = np.matrix([4,3])
p,l,u = lu(A)
print p, "#p"
print l, "#l"
print u, "#u"
print p.dot(l).dot(u),"#plu"
print A ,"#A"
程序执行结果：
[[ 0. 1.]
[ 1. 0.]] #p
[[ 1. 0. ]
[ 0.4 1. ]] #l
[[ 5. 4. ]
[ 0. 1.4]] #u
[[ 2. 3.]
[ 5. 4.]] #plu
[[2 3]
[5 4]] #A
5.1.1 SciPy的LU分解求方程组的解
还是上边的例子里的矩阵，方程如下： $$ 2x_1 + 3x_2 = 4 $$ $$ 5x_1 + 4x_2 = 3 $$ 系数矩阵A为
[[2 3]
[5 4]]
b为
[4, 3]
通过程序求p、l、u然后用pl和b求y，用u和y求x的值。
#coding:utf-8
from scipy.linalg import lu,solve
import numpy as np
A = np.array([[2,3],[5,4]])
b = np.array([4,3])
# 求的p l u
p,l,u = lu(A)
print p, "#p"
print l, "#l"
print u, "#u"
# 求ply = b的y
y = solve(p.dot(l), b)
print y,"#y"
# 求ux = y的x
x = solve(u, y)
print x,"#x"
程序执行结果：
[[ 0. 1.]
[ 1. 0.]] #p
[[ 1. 0. ]
[ 0.4 1. ]] #l
[[ 5. 4. ]
[ 0. 1.4]] #u
[ 3. 2.8] #y
[-1. 2.] #x
结果最后一行输出的是x的值，即$$x =( x_1,x_2) = ( -1,2)$$ $$ 2x_1 + 3x_2 = 2 \times (-1) + 3 \times 2 = -2 + 6 = 4 $$ $$ 5x_1 + 4x_2 = 5 \times (-1) + 4 \times 2 = -5 + 8 = 3 $$
5.2 Cholesky分解
要求解线性方程组$Ax = b$，其中为对称正定矩阵，又叫平方根法，是求解对称正定线性方程组最常用的方法之一，那么可通过下面步骤求解
（1）求的Cholesky分解，得到$A= LL^T$;
（2）求解$Ly = b$，得到y;
（3）求解$L^T x= y$，得到x。 下面使用scipy.linalg模块下的cholesky函数来对系数矩阵进行求cholesky分解。
from scipy.linalg import cholesky
import numpy as np
A = np.array([[1,2,3],[2,8,8],[3,8,35]])
b = np.array([1,8,20])
l = cholesky(A, lower=True)
print l, "# L"
print np.matmul(l,l.T), "# matmul"
print l.dot(l.T), "# dot"
y = solve(l, b)
print l.dot(y), b
x = solve(l.T, y)
print x
print l.T.dot(x),y
程序执行结果：
[[ 1. 0. 0.]
[ 2. 2. 0.]
[ 3. 1. 5.]] # L
[[ 1. 2. 3.]
[ 2. 8. 8.]
[ 3. 8. 35.]] # matmul
[[ 1. 2. 3.]
[ 2. 8. 8.]
[ 3. 8. 35.]] # dot
[ 1. 8. 20.] [ 1 8 20]
[-3.12 1.22 0.56]
[ 1. 3. 2.8] [ 1. 3. 2.8]
5.3 QR分解
QR分解法是三种将矩阵分解的方式之一。这种方式，把矩阵分解成一个正交矩阵与一个上三角矩阵的积。QR 分解经常用来解线性最小二乘法问题。 scipy.linalg模块下的qr函数可以对矩阵进行QR分解操作。
from scipy.linalg import qr
import numpy as np
aa = np.array([[0,3,1],[0,4,-2],[2,1,2]])
qq, rr = qr(aa)
print qq,"# Q"
print rr, "# R"
print aa,"# A"
print qq.dot(rr),"# QR"
程序执行结果：
[[ 0. -0.6 -0.8]
[-0. -0.8 0.6]
[-1. 0. 0. ]] # Q
[[-2. -1. -2.]
[ 0. -5. 1.]
[ 0. 0. -2.]] # R
[[ 0 3 1]
[ 0 4 -2]
[ 2 1 2]] # A
[[ 0. 3. 1.]
[ 0. 4. -2.]
[ 2. 1. 2.]] # QR
5.4 SVD奇异分解
svd是现在比较常见的算法之一，也是数据挖掘工程师、算法工程师必备的技能之一。 假设A是一个$M \times N$的矩阵，那么通过矩阵分解将会得到U，Σ，$V^T$（V的转置）三个矩阵，其中U是一个$M \times M$的方阵，被称为左奇异向量，方阵里面的向量是正交的；Σ是一个$M \times N$的对角矩阵，除了对角线的元素其他都是0，对角线上的值称为奇异值；$V^T$(V的转置)是一个$N \times N$的矩阵，被称为右奇异向量，方阵里面的向量也都是正交的。 $$ A_{m \times n} = U_{m \times m}\Sigma_{m \times n}V_{n \times n}^T $$
from scipy.linalg import qr,svd
import numpy as np
aa = np.array([[0,3,1],[0,4,-2],[2,1,2]])
u, e, v = svd(aa)
print u,"#u"
print e,"#e"
print v,"#v"
程序执行结果：
[[-0.54374742 0.34887107 -0.76330054]
[-0.82489643 -0.38965121 0.40953366]
[-0.15454653 0.85232676 0.49965434]] #u
[ 5.15671459 3.32363965 1.16692507] #e
[[-0.05993992 -0.98616559 0.15454653]
[ 0.51288759 0.10239833 0.85232676]
[ 0.85636063 -0.1303534 -0.49965434]] #v
感谢Klang(金浪)智能数据看板klang.org.cn鼎力支持！