线性代数

线性代数

numpy 和 scipy 中，负责进行线性代数部分计算的模块叫做 linalg。

In [1]:

import numpy as np
import numpy.linalg
import scipy as sp
import scipy.linalg
import matplotlib.pyplot as plt
from scipy import linalg
 
%matplotlib inline

numpy.linalg VS scipy.linalg

一方面scipy.linalg 包含 numpy.linalg 中的所有函数，同时还包含了很多 numpy.linalg 中没有的函数。

另一方面，scipy.linalg 能够保证这些函数使用 BLAS/LAPACK 加速，而 numpy.linalg 中这些加速是可选的。

因此，在使用时，我们一般使用 scipy.linalg 而不是 numpy.linalg。

我们可以简单看看两个模块的差异：

In [2]:

print "number of items in numpy.linalg:", len(dir(numpy.linalg))
print "number of items in scipy.linalg:", len(dir(scipy.linalg))

number of items in numpy.linalg: 36
number of items in scipy.linalg: 115

numpy.matrix VS 2D numpy.ndarray

线性代数的基本操作对象是矩阵，而矩阵的表示方法主要有两种：numpy.matrix 和 2D numpy.ndarray。

numpy.matrix

numpy.matrix 是一个矩阵类，提供了一些方便的矩阵操作：

支持类似 MATLAB 创建矩阵的语法
矩阵乘法默认用 * 号
.I 表示逆，.T 表示转置可以用 mat 或者 matrix 来产生矩阵：

In [3]:

A = np.mat("[1, 2; 3, 4]")
print repr(A)
 
A = np.matrix("[1, 2; 3, 4]")
print repr(A)

matrix([[1, 2],
        [3, 4]])
matrix([[1, 2],
        [3, 4]])

转置和逆：

In [4]:

print repr(A.I)
print repr(A.T)

matrix([[-2. ,  1. ],
        [ 1.5, -0.5]])
matrix([[1, 3],
        [2, 4]])

矩阵乘法：

In [5]:

b = np.mat('[5; 6]')
print repr(A * b)

matrix([[17],
        [39]])

2 维 numpy.ndarray

虽然 numpy.matrix 有着上面的好处，但是一般不建议使用，而是用 2 维 numpy.ndarray 对象替代，这样可以避免一些不必要的困惑。

我们可以使用 array 复现上面的操作：

In [6]:

A = np.array([[1,2], [3,4]])
print repr(A)

array([[1, 2],
       [3, 4]])

逆和转置：

In [7]:

print repr(linalg.inv(A))
print repr(A.T)

array([[-2. ,  1. ],
       [ 1.5, -0.5]])
array([[1, 3],
       [2, 4]])

矩阵乘法：

In [8]:

b = np.array([5, 6])
 
print repr(A.dot(b))

array([17, 39])

普通乘法：

In [9]:

print repr(A * b)

array([[ 5, 12],
       [15, 24]])

scipy.linalg 的操作可以作用到两种类型的对象上，没有区别。

基本操作

求逆

矩阵 $\mathbf{A}$ 的逆 $\mathbf{B}$ 满足：$\mathbf{BA}=\mathbf{AB}=I$，记作 $\mathbf{B} = \mathbf{A}^{-1}$。

事实上，我们已经见过求逆的操作，linalg.inv 可以求一个可逆矩阵的逆：

In [10]:

A = np.array([[1,2],[3,4]])
 
print linalg.inv(A)
 
print A.dot(scipy.linalg.inv(A))

[[-2.   1. ]
 [ 1.5 -0.5]]
[[  1.00000000e+00   0.00000000e+00]
 [  8.88178420e-16   1.00000000e+00]]

求解线性方程组

例如，下列方程组\begin{eqnarray} x + 3y + 5z & = & 10 \2x + 5y + z & = & 8 \2x + 3y + 8z & = & 3\end{eqnarray}的解为：\begin{split}\left[\begin{array}{c} x\ y\ z\end{array}\right]=\left[\begin{array}{ccc} 1 & 3 & 5\ 2 & 5 & 1\ 2 & 3 & 8\end{array}\right]^{-1}\left[\begin{array}{c} 10\ 8\ 3\end{array}\right]=\frac{1}{25}\left[\begin{array}{c} -232\ 129\ 19\end{array}\right]=\left[\begin{array}{c} -9.28\ 5.16\ 0.76\end{array}\right].\end{split}

我们可以使用 linalg.solve 求解方程组，也可以先求逆再相乘，两者中 solve 比较快。

In [11]:

import time
 
A = np.array([[1, 3, 5],
              [2, 5, 1],
              [2, 3, 8]])
b = np.array([10, 8, 3])
 
tic = time.time()
 
for i in xrange(1000):
    x = linalg.inv(A).dot(b)
 
print x
print A.dot(x)-b
print "inv and dot: {} s".format(time.time() - tic)
 
tic = time.time()
 
for i in xrange(1000):
    x = linalg.solve(A, b)
 
print x
print A.dot(x)-b
print "solve: {} s".format(time.time() - tic)

[-9.28  5.16  0.76]
[  0.00000000e+00  -1.77635684e-15  -8.88178420e-16]
inv and dot: 0.0353579521179 s
[-9.28  5.16  0.76]
[  0.00000000e+00  -1.77635684e-15  -1.77635684e-15]
solve: 0.0284671783447 s

计算行列式

方阵的行列式为\left|\mathbf{A}\right|=\sum{j}\left(-1\right)^{i+j}a{ij}M_{ij}.

其中 $a{ij}$ 表示 $\mathbf{A}$ 的第 $i$ 行第 $j$ 列的元素，$M{ij}$ 表示矩阵 $\mathbf{A}$ 去掉第 $i$ 行第 $j$ 列的新矩阵的行列式。

例如，矩阵\begin{split}\mathbf{A=}\left[\begin{array}{ccc} 1 & 3 & 5\ 2 & 5 & 1\ 2 & 3 & 8\end{array}\right]\end{split}的行列式是：\begin{eqnarray} \left|\mathbf{A}\right| & = & 1\left|\begin{array}{cc} 5 & 1\ 3 & 8\end{array}\right|-3\left|\begin{array}{cc} 2 & 1\ 2 & 8\end{array}\right|+5\left|\begin{array}{cc} 2 & 5\ 2 & 3\end{array}\right|\ & = & 1\left(5\cdot8-3\cdot1\right)-3\left(2\cdot8-2\cdot1\right)+5\left(2\cdot3-2\cdot5\right)=-25.\end{eqnarray}

可以用 linalg.det 计算行列式：

In [12]:

A = np.array([[1, 3, 5],
              [2, 5, 1],
              [2, 3, 8]])
 
print linalg.det(A)

-25.0

计算矩阵或向量的模

矩阵的模定义如下：\begin{split}\left\Vert \mathbf{A}\right\Vert =\left{ \begin{array}{cc} \max{i}\sum{j}\left|a{ij}\right| & \textrm{ord}=\textrm{inf}\ \min{i}\sum{j}\left|a{ij}\right| & \textrm{ord}=-\textrm{inf}\ \max{j}\sum{i}\left|a{ij}\right| & \textrm{ord}=1\ \min{j}\sum{i}\left|a{ij}\right| & \textrm{ord}=-1\ \max\sigma{i} & \textrm{ord}=2\ \min\sigma{i} & \textrm{ord}=-2\ \sqrt{\textrm{trace}\left(\mathbf{A}^{H}\mathbf{A}\right)} & \textrm{ord}=\textrm{'fro'}\end{array}\right.\end{split}其中，$\sigma_i$ 是矩阵的奇异值。

向量的模定义如下：\begin{split}\left\Vert \mathbf{x}\right\Vert =\left{ \begin{array}{cc} \max\left|x{i}\right| & \textrm{ord}=\textrm{inf}\ \min\left|x{i}\right| & \textrm{ord}=-\textrm{inf}\ \left(\sum{i}\left|x{i}\right|^{\textrm{ord}}\right)^{1/\textrm{ord}} & \left|\textrm{ord}\right|<\infty.\end{array}\right.\end{split}

linalg.norm 可以计算向量或者矩阵的模：

In [13]:

A = np.array([[1, 2],
              [3, 4]])
 
print linalg.norm(A)
 
print linalg.norm(A,'fro') # frobenius norm 默认值
 
print linalg.norm(A,1) # L1 norm 最大列和
 
print linalg.norm(A,-1) # L -1 norm 最小列和
 
print linalg.norm(A,np.inf) # L inf norm 最大行和

5.47722557505
5.47722557505
6
4
7

最小二乘解和伪逆

问题描述

所谓最小二乘问题的定义如下：

假设 $y_i$ 与 $\mathbf{x_i}$ 的关系可以用一组系数 $c_j$ 和对应的模型函数 $f_j(\mathbf{x_i})$ 的模型表示：

y{i}=\sum{j}c{j}f{j}\left(\mathbf{x}{i}\right)+\epsilon{i}

其中 $\epsiloni$ 表示数据的不确定性。最小二乘就是要优化这样一个关于 $c_j$ 的问题：J\left(\mathbf{c}\right)=\sum{i}\left|y{i}-\sum{j}c{j}f{j}\left(x_{i}\right)\right|^{2}

其理论解满足：\frac{\partial J}{\partial c{n}^{*}}=0=\sum{i}\left(y{i}-\sum{j}c{j}f{j}\left(x{i}\right)\right)\left(-f{n}^{*}\left(x_{i}\right)\right)

改写为：\begin{eqnarray} \sum{j}c{j}\sum{i}f{j}\left(x{i}\right)f{n}^{}\left(x{i}\right) & = & \sum{i}y{i}f{n}^{}\left(x_{i}\right)\ \mathbf{A}^{H}\mathbf{Ac} & = & \mathbf{A}^{H}\mathbf{y}\end{eqnarray}

其中：\left{ \mathbf{A}\right} {ij}=f{j}\left(x_{i}\right).

当 $\mathbf{A^HA}$ 可逆时，我们有：\mathbf{c}=\left(\mathbf{A}^{H}\mathbf{A}\right)^{-1}\mathbf{A}^{H}\mathbf{y}=\mathbf{A}^{\dagger}\mathbf{y}

矩阵 $\mathbf{A}^{\dagger}$ 叫做 $\mathbf{A}$ 的伪逆。

问题求解

注意到，我们的模型可以写为：\mathbf{y}=\mathbf{Ac}+\boldsymbol{\epsilon}.

在给定 $\mathbf{y}$ 和 $\mathbf{A}$ 的情况下，我们可以使用 linalg.lstsq 求解 $\mathbf c$。

在给定 $\mathbf{A}$ 的情况下，我们可以使用 linalg.pinv 或者 linalg.pinv2 求解 $\mathbf{A}^{\dagger}$。

例子

假设我们的数据满足：\begin{align}y{i} & =c{1}e^{-x{i}}+c{2}x{i} \z{i} & = y_i + \epsilon_i\end{align}

其中 $x_i = \frac{i}{10},\ i = 1,\dots,10$，$c_1 = 5, c_2 = 2$，产生数据

In [14]:

c1, c2 = 5.0, 2.0
i = np.r_[1:11]
xi = 0.1*i
yi = c1*np.exp(-xi) + c2*xi
zi = yi + 0.05 * np.max(yi) * np.random.randn(len(yi))

构造矩阵 $\mathbf A$：

In [15]:

A = np.c_[np.exp(-xi)[:, np.newaxis], xi[:, np.newaxis]]
print A

[[ 0.90483742  0.1       ]
 [ 0.81873075  0.2       ]
 [ 0.74081822  0.3       ]
 [ 0.67032005  0.4       ]
 [ 0.60653066  0.5       ]
 [ 0.54881164  0.6       ]
 [ 0.4965853   0.7       ]
 [ 0.44932896  0.8       ]
 [ 0.40656966  0.9       ]
 [ 0.36787944  1.        ]]

求解最小二乘问题：

In [16]:

c, resid, rank, sigma = linalg.lstsq(A, zi)
 
print c

[ 4.87016856  2.19081311]

其中 c 的形状与 zi 一致，为最小二乘解，resid 为 zi - A c 每一列差值的二范数，rank 为矩阵 A 的秩，sigma 为矩阵 A 的奇异值。

查看拟合效果：

In [17]:

xi2 = np.r_[0.1:1.0:100j]
yi2 = c[0]*np.exp(-xi2) + c[1]*xi2
 
plt.plot(xi,zi,'x',xi2,yi2)
plt.axis([0,1.1,3.0,5.5])
plt.xlabel('$x_i$')
plt.title('Data fitting with linalg.lstsq')
plt.show()

广义逆

linalg.pinv 或 linalg.pinv2 可以用来求广义逆，其区别在于前者使用求最小二乘解的算法，后者使用求奇异值的算法求解。

矩阵分解

特征值和特征向量

问题描述

对于给定的 $N \times N$ 矩阵 $\mathbf A$，特征值和特征向量问题相当与寻找标量 $\lambda$ 和对应的向量 $\mathbf v$ 使得：\mathbf{Av} = \lambda \mathbf{v}

矩阵的 $N$ 个特征值（可能相同）可以通过计算特征方程的根得到：\left|\mathbf{A} - \lambda \mathbf{I}\right| = 0

然后利用这些特征值求（归一化的）特征向量。

问题求解

linalg.eig(A)
- 返回矩阵的特征值与特征向量
linalg.eigvals(A)
- 返回矩阵的特征值
linalg.eig(A, B)
- 求解 $\mathbf{Av} = \lambda\mathbf{Bv}$ 的问题

例子

矩阵为\begin{split}\mathbf{A}=\left[\begin{array}{ccc} 1 & 5 & 2\ 2 & 4 & 1\ 3 & 6 & 2\end{array}\right].\end{split}

特征多项式为：\begin{eqnarray} \left|\mathbf{A}-\lambda\mathbf{I}\right| & = & \left(1-\lambda\right)\left[\left(4-\lambda\right)\left(2-\lambda\right)-6\right]-\ & & 5\left[2\left(2-\lambda\right)-3\right]+2\left[12-3\left(4-\lambda\right)\right]\ & = & -\lambda^{3}+7\lambda^{2}+8\lambda-3.\end{eqnarray}

特征根为：\begin{eqnarray} \lambda{1} & = & 7.9579\ \lambda{2} & = & -1.2577\ \lambda_{3} & = & 0.2997.\end{eqnarray}

In [18]:

A = np.array([[1, 5, 2], 
              [2, 4, 1],
              [3, 6, 2]])
 
la, v = linalg.eig(A)
 
print la
 
 
# 验证是否归一化
print np.sum(abs(v**2),axis=0)
 
# 第一个特征值
l1 = la[0]
# 对应的特征向量
v1 = v[:, 0].T
 
# 验证是否为特征值和特征向量对
print linalg.norm(A.dot(v1)-l1*v1)

[ 7.95791620+0.j -1.25766471+0.j  0.29974850+0.j]
[ 1.  1.  1.]
3.23301824835e-15

奇异值分解

问题描述

$M \times N$ 矩阵 $\mathbf A$ 的奇异值分解为：\mathbf{A=U}\boldsymbol{\Sigma}\mathbf{V}^{H}

其中 $\boldsymbol{\Sigma}, (M \times N)$ 只有对角线上的元素不为 0，$\mathbf U, (M \times M)$ 和 $\mathbf V, (N \times N)$ 为正交矩阵。

其具体原理可以查看维基百科：https://en.wikipedia.org/wiki/Singular_value_decomposition

问题求解

U,s,Vh = linalg.svd(A)
- 返回 $U$ 矩阵，奇异值 $s$，$V^H$ 矩阵
Sig = linalg.diagsvd(s,M,N)
- 从奇异值恢复 $\boldsymbol{\Sigma}$ 矩阵

例子

奇异值分解：

In [19]:

A = np.array([[1,2,3],[4,5,6]])
 
U, s, Vh = linalg.svd(A)

$\boldsymbol{\Sigma}$ 矩阵：

In [20]:

M, N = A.shape
Sig = linalg.diagsvd(s,M,N)
 
print Sig

[[ 9.508032    0.          0.        ]
 [ 0.          0.77286964  0.        ]]

检查正确性：

In [21]:

print A
print U.dot(Sig.dot(Vh))

[[1 2 3]
 [4 5 6]]
[[ 1.  2.  3.]
 [ 4.  5.  6.]]

LU 分解

$M \times N$ 矩阵 $\mathbf A$ 的 LU 分解为：\mathbf{A}=\mathbf{P}\,\mathbf{L}\,\mathbf{U}

$\mathbf P$ 是 $M \times M$ 的单位矩阵的一个排列，$\mathbf L$ 是下三角阵，$\mathbf U$ 是上三角阵。

可以使用 linalg.lu 进行 LU 分解的求解：

具体原理可以查看维基百科：https://en.wikipedia.org/wiki/LU_decomposition

In [22]:

A = np.array([[1,2,3],[4,5,6]])
 
P, L, U = linalg.lu(A)
 
print P
print L
print U
 
print P.dot(L).dot(U)

[[ 0.  1.]
 [ 1.  0.]]
[[ 1.    0.  ]
 [ 0.25  1.  ]]
[[ 4.    5.    6.  ]
 [ 0.    0.75  1.5 ]]
[[ 1.  2.  3.]
 [ 4.  5.  6.]]

Cholesky 分解

Cholesky 分解是一种特殊的 LU 分解，此时要求 $\mathbf A$ 为 Hermitian 正定矩阵（$\mathbf A = \mathbf{A^H}$）。

此时有：\begin{eqnarray} \mathbf{A} & = & \mathbf{U}^{H}\mathbf{U}\ \mathbf{A} & = & \mathbf{L}\mathbf{L}^{H}\end{eqnarray}即\mathbf{L}=\mathbf{U}^{H}.

可以用 linalg.cholesky 求解。

QR 分解

$M×N$ 矩阵 $\mathbf A$ 的 QR 分解为：\mathbf{A=QR}

$\mathbf R$ 为上三角形矩阵，$\mathbf Q$ 是正交矩阵。

维基链接：https://en.wikipedia.org/wiki/QR_decomposition

可以用 linalg.qr 求解。

Schur 分解

对于 $N\times N$ 方阵 $\mathbf A$, Schur 分解要求找到满足下式的矩阵：\mathbf{A=ZTZ^H}

其中 $\mathbf Z$ 是正交矩阵，$\mathbf T$ 是一个上三角矩阵。

维基链接：https://en.wikipedia.org/wiki/Schur_decomposition

In [23]:

A = np.mat('[1 3 2; 1 4 5; 2 3 6]')
 
print A
 
T, Z = linalg.schur(A)
 
print T, Z
 
print Z.dot(T).dot(Z.T)

[[1 3 2]
 [1 4 5]
 [2 3 6]]
[[ 9.90012467  1.78947961 -0.65498528]
 [ 0.          0.54993766 -1.57754789]
 [ 0.          0.51260928  0.54993766]] [[ 0.36702395 -0.85002495 -0.37782404]
 [ 0.63681656 -0.06646488  0.76814522]
 [ 0.67805463  0.52253231 -0.51691576]]
[[ 1.  3.  2.]
 [ 1.  4.  5.]
 [ 2.  3.  6.]]

矩阵函数

考虑函数 $f(x)$ 的泰勒展开：f\left(x\right)=\sum_{k=0}^{\infty}\frac{f^{\left(k\right)}\left(0\right)}{k!}x^{k}

对于方阵，矩阵函数可以定义如下：f\left(\mathbf{A}\right)=\sum_{k=0}^{\infty}\frac{f^{\left(k\right)}\left(0\right)}{k!}\mathbf{A}^{k}

这也是计算矩阵函数的最好的方式。

指数和对数函数

指数

指数可以定义如下：e^{\mathbf{A}}=\sum_{k=0}^{\infty}\frac{1}{k!}\mathbf{A}^{k}

linalg.expm3 使用的是泰勒展开的方法计算结果：

In [24]:

A = np.array([[1, 2], [3, 4]])
 
print linalg.expm3(A)

[[  51.96890355   74.73648784]
 [ 112.10473176  164.07363531]]

另一种方法先计算 A 的特征值分解：\mathbf{A}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1}

然后有（正交矩阵和对角阵的性质）：e^{\mathbf{A}}=\mathbf{V}e^{\boldsymbol{\Lambda}}\mathbf{V}^{-1}

linalg.expm2 使用的就是这种方法：

In [25]:

print linalg.expm2(A)

[[  51.9689562    74.73656457]
 [ 112.10484685  164.07380305]]

最优的方法是用 Padé 近似实现，Padé 近似往往比截断的泰勒级数准确，而且当泰勒级数不收敛时，Padé 近似往往仍可行，所以多用于在计算机数学中。

linalg.expm 使用的就是这种方法：

In [26]:

print linalg.expm(A)

[[  51.9689562    74.73656457]
 [ 112.10484685  164.07380305]]

对数

指数的逆运算，可以用 linalg.logm 实现：

In [27]:

print A
print linalg.logm(linalg.expm(A))

[[1 2]
 [3 4]]
[[ 1.  2.]
 [ 3.  4.]]

三角函数

根据欧拉公式，其定义为：\begin{eqnarray} \sin\left(\mathbf{A}\right) & = & \frac{e^{j\mathbf{A}}-e^{-j\mathbf{A}}}{2j}\ \cos\left(\mathbf{A}\right) & = & \frac{e^{j\mathbf{A}}+e^{-j\mathbf{A}}}{2}.\end{eqnarray}

正切函数定义为：\tan\left(x\right)=\frac{\sin\left(x\right)}{\cos\left(x\right)}=\left[\cos\left(x\right)\right]^{-1}\sin\left(x\right)

因此矩阵的正切函数定义为：\left[\cos\left(\mathbf{A}\right)\right]^{-1}\sin\left(\mathbf{A}\right).

具体实现：

linalg.sinm
linalg.cosm
linalg.tanm

双曲三角函数

\begin{eqnarray} \sinh\left(\mathbf{A}\right) & = & \frac{e^{\mathbf{A}}-e^{-\mathbf{A}}}{2}\ \cosh\left(\mathbf{A}\right) & = & \frac{e^{\mathbf{A}}+e^{-\mathbf{A}}}{2}\ \tanh\left(\mathbf{A}\right) & = & \left[\cosh\left(\mathbf{A}\right)\right]^{-1}\sinh\left(\mathbf{A}\right).\end{eqnarray}

具体实现：

linalg.sinhm
linalg.coshm
linalg.tanhm

特殊矩阵

Scipy 提供了一些特殊矩阵的实现，具体可以参考：

http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html#special-matrices

原文: https://nbviewer.jupyter.org/github/lijin-THU/notes-python/blob/master/04-scipy/04.09-linear-algbra.ipynb

04.09 线性代数

线性代数

numpy.linalg VS scipy.linalg

numpy.matrix VS 2D numpy.ndarray

numpy.matrix

2 维 numpy.ndarray

基本操作

求逆

求解线性方程组

计算行列式

计算矩阵或向量的模

最小二乘解和伪逆

问题描述

问题求解

例子

广义逆

矩阵分解

特征值和特征向量

问题描述

问题求解

例子

奇异值分解

问题描述

问题求解

例子

LU 分解

Cholesky 分解

QR 分解

Schur 分解

矩阵函数

指数和对数函数

指数

对数

三角函数

双曲三角函数

特殊矩阵