NumPy和SciPy都提供了线性代数函数库linalg，SciPy的线性代数库比NumPy更加全面。

一、解线性方程组

numpy.linalg.solve(A, b)和scipy.linalg.solve(A, b)可以用来解线性方程组Ax=b，也就是计算x=A-1b。这里，A是m×n的方形矩阵，x和b是长为m的向量。有时候A是固定的，需要对多组b进行求解，因此第二个参数也可以是m×n的矩阵B。这样计算出来的X也是m×n的矩阵，相当于计算A-1B。

在一些矩阵公式中经常会出现类似于A-1B的运算，它们都可以用solve(A,B)计算，这要比直接计算逆矩阵然后做矩阵乘法更快捷一些，下面的程序比较solve( )和逆矩阵的运算速度：

import numpy as np
from scipy import linalg
－^1－﹣﹣
m, n = 500, 50
A = np.random.rand(m, m)
B = np.random.rand(m, n)
X1 = linalg.solve(A, B)
X2 = np.dot(linalg.inv(A), B)
print np.allclose(X1, X2)
%timeit linalg.solve(A, B)
%timeit np.dot(linalg.inv(A), B)
True
100 loops, best of 3: 10.1 ms per loop
10 loops, best of 3: 20 ms per loop

若需要对多组B进行求解，但是又不好将它们合并成一个矩阵，例如某些矩阵公式中可能会有A-1B、A-1C、A-1D等乘法，而B、C、D是通过某种方式逐次计算的。这时可以采用lu_factor( )和lu_solve( )。先调用lu_factor(A)对矩阵A进行LU分解，得到一个元组：（LU矩阵，排序数组）。将这个元组传递给lu_solve( )，即可对不同的B进行求解。由于已经对A进行了LU分解，lu_solve( )能够很快得出结果。

luf = linalg.lu_factor(A)
X3 = linalg.lu_solve(luf, B)
np.allclose(X1, X3)
True

除了使用lu_factor( )和lu_solve( )之外，可以先通过inv( )计算逆矩阵，然后通过dot( )计算矩阵乘积。下面比较二者的速度，可以看出lu_factor( )比inv( )要快很多，而lu_solve( )和dot( )的运算速度几乎相同：

M, N = 1000, 100
np.random.seed(0)
A = np.random.rand(M, M)
B = np.random.rand(M, N)
Ai = linalg.inv(A)
luf = linalg.lu_factor(A)   
%timeit linalg.inv(A)
%timeit np.dot(Ai, B)
%timeit linalg.lu_factor(A)    
%timeit linalg.lu_solve(luf, B)
10 loops, best of 3: 131 ms per loop
100 loops, best of 3: 9.65 ms per loop
10 loops, best of 3: 52.6 ms per loop
100 loops, best of 3: 13.8 ms per loop

二、最小二乘解

lstsq( )比solve( )更一般化，它不要求矩阵A是正方形的，也就是说方程的个数可以少于、等于或多于未知数的个数。它找到一组解x，使得‖b-Ax‖最小。我们称得到的结果为最小二乘解，即它使得所有等式的误差的平方和最小。下面以求解离散卷积的逆运算为例，介绍lstsq( )的用法。

首先简单介绍一下离散卷积的相关知识和计算方法。对于离散的线性时不变系统h，如果它的输入是x，那么其输出y可以用x和h的卷积表示：y=x*h。

离散卷积的计算公式如下：

y[n]=∑h[m]x[n-m]

假设h的长度为n，x的长度为m，则卷积计算所得到的y的长度将为n+m-1。它的每个值都是按照下面的公式计算得到的：

y[0] = h[0]*x[0]

y[1] = h[0]*x[1] + h[1]*x[0]

y[2] = h[0]*x[2] + h[1]*x[1] + h[2]*x[0]

y[3] = h[0]*x[3] + h[1]*x[2] + h[2]*x[1]

...

y[n+m-1] = h[n-1]*x[m-1]

所谓卷积的逆运算就是指：假设已知x和y，需要求解h。由于h的长度为n，于是有n个未知数，而由于y的长度为n+m-1，因此这n个未知数需要满足n+m-1个线性方程。由于方程数比未知数多，卷积的逆运算不一定有精确解，因此问题就变成了找到一组h，使得x*h与y之间的误差最小，显然它就是最小二乘解。下面的程序演示了如何使用lstsq( )计算卷积的逆运算：

首先make_data( )创建所需的数据，它使用随机数函数standard_normal( )初始化数组x和h。在实际的系统中h通常是未知的，并且值会逐渐衰减。make_data( )返回系统的输入信号x以及添加了随机噪声的输出信号yn。为了和最小二乘法的结果相比较，我们同时也输出了系统的系数h。

solve_h( )使用最小二乘法计算系统的参数h，因为通常我们不知道未知系统的系数的长度，因此这里用N表示所求系数的长度。

观察前面的卷积方程组可知，在n+m-1个方程中，中间的n-m+1个方程使用了h的所有系数。为了程序计算方便，我们对这m-n+1个方程进行最小二乘运算。根据h的长度，需要将一维数组x变换成一个形状为(m-n+1, n)的二维数组X，它的每行相对于上一行都左移了一个元素。这个二维数组可以很容易地使用第2章中介绍过的as_strided( )得到。我们取出输出数组y中与数组X每行对应的部分，然后调用lstsq( )对这m-n+1个方程进行最小二乘运算。lstsq( )返回一个元组，它的第0个元素是最小二乘解，注意得到的结果顺序是颠倒的，因此还需要对其进行翻转。

from numpy.lib.stride_tricks import as_strided

def make_data(m, n, noise_scale): 
np.random.seed(42)
x = np.random.standard_normal(m)
h = np.random.standard_normal(n)
y = np.convolve(x, h)
yn = y + np.random.standard_normal(len(y)) * noise_scale * np.max(y)
return x, yn, h

def solve_h(x, y, n):      
X = as_strided(x, shape=(len(x)-n+1, n), strides=(x.itemsize, x.itemsize)) 
Y = y[n-1:len(x)]      
h = linalg.lstsq(X, Y) 
return h[0][::-1]      

接下来对长度为100的未知系统系数h，分别计算长度为80和120的最小二乘解。由于我们对系统的输出添加了一些噪声信号，因此二者并不完全吻合。图1比较了这两个解与真实系数。

x, yn, h = make_data(1000, 100, 0.4)   
H1 = solve_h(x, yn, 120)
H2 = solve_h(x, yn, 80)

print "Average error of H1:", np.mean(np.abs(H[:100] - h))
print "Average error of H2:", np.mean(np.abs(h[:80] - H2))
Average error of H1: 0.301548854044
Average error of H2: 0.295842215834

图1　实际的系统参数与最小二乘解的比较

三、特征值和特征向量

n×n的矩阵A可以看作n维空间中的线性变换。若x为n维空间中的一个向量，那么A与x的矩阵乘积就是对x进行线性变换之后的向量。如果x是线性变换的特征向量，那么经过这个线性变换之后，得到的新向量仍然与原来的x保持在同一方向上，但其长度也许会改变。特征向量的长度在该线性变换下缩放的比例称为其特征值。即特征向量x满足如下等式，λ的值可以是一个任意复数：

Ax=λx

下面以二维平面上的线性变换矩阵为例，演示特征值和特征变量的几何含义。通过linalg.eig(A)计算矩阵A的两个特征值evalues和特征向量evectors，在evectors中，每一列是一个特征向量。

A = np.array([[1, -0.3], [-0.1, 0.9]])
evalues, evectors = linalg.eig(A)
evalues                          evectors          
----------------------------------  ----------------------------
[ 1.13027756+0.j,  0.76972244+0.j]  [[ 0.91724574,  0.79325185],
[-0.3983218 ,  0.60889368]]

图2显示了变换前后的向量。在图中，蓝色箭头为变换之前的向量，红色箭头为变换之后的向量。粗箭头为变换前后的特征向量。可以看出特征向量变换前后方向没有改变，只是长度发生了变化。长度的变化倍数由特征值决定：一个变为原来的1.13倍长，一个变为原来的0.77倍长。

图2　线性变换将蓝色箭头变换为红色箭头

numpy.linalg模块中也有eig( )函数，与之不同的是，scipy.linalg模块中的eig( )函数支持计算广义特征值和广义特征向量，它们满足如下等式，其中B是一个n×n的矩阵：

Ax=λBx

广义特征向量可以用于椭圆拟合，椭圆拟合的公式与原理可以参考下面的论文：

http://research.microsoft.com/pubs/67845/ellipse-pami.pdf

用广义特征向量计算椭圆拟合。

椭圆上的点满足如下方程，其中a,b,c,d,e,f为椭圆的参数，(x,y)为平面上的坐标点：

f(x,y)=ax2+bxy+cy2+dx+ey+f=0

所谓椭圆拟合，就是指给出一组平面上的点(xi,yi)，找到一组椭圆参数，使得∑f(xi,yi)2最小。显然这是一个最小化问题，可以使用上节介绍的优化算法optimize.leastsq( )求解。为了避免参数全为0的平凡解，需要一点小技巧，请读者自行演练一下。下面给出论文中用广义特征向量计算椭圆拟合的方法：

首先定义。D是一个n×6的矩阵，其中n为点的个数，D中的每一行与一个坐标点相对应。a为拟合椭圆的系数：a=[a,b,c,d,e,f]T。则a满足如下方程：

DTDa=λCa

其中C是一个6×6的矩阵：

显然上式符合广义特征向量的等式，因此可以用linalg.eig( )求解。下面首先使用椭圆的参数方程计算某个椭圆上随机的60个点，并引入一些随机噪声：

np.random.seed(42)
t = np.random.uniform(0, 2*np.pi, 60)

alpha = 0.4
a = 0.5
b = 1.0
x = 1.0 + a*np.cos(t)*np.cos(alpha) - b*np.sin(t)*np.sin(alpha)
y = 1.0 + a*np.cos(t)*np.sin(alpha) - b*np.sin(t)*np.cos(alpha)
x += np.random.normal(0, 0.05, size=len(x))
y += np.random.normal(0, 0.05, size=len(y))

当传递第二个参数时，eig( )计算广义特征值和向量。evectors中共有6个特征向量，将这6个特征向量代入椭圆方程中，计算平均误差，并挑选误差最小的特征向量作为椭圆的参数p。图3显示了参数p所表示的椭圆以及数据点。

D = np.c_[x**2, x*y, y**2, x, y, np.ones_like(x)]
A = np.dot(D.T, D)
C = np.zeros((6, 6))
C[[0, 1, 2], [2, 1, 0]] = 2, -1, 2
evalues, evectors = linalg.eig(A, C)     
evectors = np.real(evectors)
err = np.mean(np.dot(D, evectors)**2, 0) 
p = evectors[:, np.argmin(err) ]         
print p
[-0.55214278  0.5580915  -0.23809922  0.54584559 -0.08350449 -0.14852803]

图3　用广义特征向量计算的拟合椭圆

四、奇异值分解-SVD

奇异值分解是线性代数中一种重要的矩阵分解，在信号处理、统计学等领域都有重要应用。假设M是一个m×n阶矩阵，存在一个分解使得：M=UΣV*。其中U是m×m阶酉矩阵；Σ是半正定m×n阶对角矩阵；而V*，即V的共轭转置，是n×n阶酉矩阵。这样的分解就称作M的奇异值分解。Σ对角线上的元素为M的奇异值。通常奇异值按照从大到小的顺序排列。

奇异值的数学描述读起来有些难懂，让我们通过一个实例说明奇异值分解的用途。下面的例子对一幅灰度图像进行奇异值分解，然后从三个分解矩阵中提取奇异值较大的部分数据还原原始图像。首先读入一幅图像，并通过其红绿蓝三个通道计算灰度图像img，图像的宽为375个像素、高为505个像素。

r, g, b = np.rollaxis(pl.imread("vinci_target.png"), 2).astype(float)
img = 0.2989 * r + 0.5870 * g + 0.1140 * b
img.shape
(505, 375)

调用scipy.linalg.svd( )对图像矩阵进行奇异值分解，得到三个分解部分：

·U：对应于公式中的U。

·s：对应于公式中的Σ，由于它是一个对角矩阵，只有对角线上的元素非零，因此s是一个一维数组，保存对角线上的非零元素。

·Vh：对应于公式中的V*。

下面的程序查看这三个数组的形状：

U, s, Vh = linalg.svd(img)
U.shape       s.shape      Vh.shape
----------     -------     ----------
(505, 505)     (375,)      (375, 375)

s中的每个值与Vh的行向量以及U中的列向量对应，默认按照从大到小的顺序排列，它表示与其对应的向量的重要性。由图4可知s中的奇异值大小差别很大，注意Y轴是对数坐标系。

pl.semilogy(s, lw=3)

图4　按从大到小排列的奇异值

下面的composite( )选择U和Vh中的前n个向量重新合成矩阵，当用上所有向量时，重新合成的矩阵和原始矩阵相同：

def composite(U, s, Vh, n):
return np.dot(U[:, :n], s[:n, np.newaxis] * Vh[:n, :])

print np.allclose(img, composite(U, s, Vh, len(s)))
True

下面演示选择前10个、20个以及50个向量合成时的效果，如图5所示，可以看到使用的向量越多，结果就越接近原始图像：

img10 = composite(U, s, Vh, 10)
img20 = composite(U, s, Vh, 20)
img50 = composite(U, s, Vh, 50)

%array_image img; img10; img20; img50

图5　原始图像以及使用10个、20个、50个向量合成的图像（从左到右）