3.3.2 最小二乘解

lstsq( )比solve( )更一般化，它不要求矩阵A是正方形的，也就是说方程的个数可以少于、等于或多于未知数的个数。它找到一组解x，使得‖b-A x‖最小。我们称得到的结果为最小二乘解，即它使得所有等式的误差的平方和最小。下面以求解离散卷积的逆运算为例，介绍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的最小二乘解。由于我们对系统的输出添加了一些噪声信号，因此二者并不完全吻合。图3-6比较了这两个解与真实系数。

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

    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