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