【Hello Python World】Week 13：SciPy：Fast and scientific

今天的任务是关于scipy的练习，上周的练习是matplotlib的练习，题目有不少是承接这一节的练习的，上一周的难题解决了，这一周的练习相对比较轻松。

最小二乘法

分析

关于公式$\boldsymbol{\hat{x}}=\arg\min_{\boldsymbol{x}}\left \| A\boldsymbol{x}-\boldsymbol{b} \right \|_2$，相信不少人会和我一开始想的一样，想要用Python中的argmin方法来解决。但是这个地方并不适合使用这些统计的方法，因为没有一个可供筛选的数据集啊！什么意思呢？想想一下使用argmin或者min的时候，有一个参数就是要有可迭代的数据集(比如列表)，答案已经在这个数据集里了，这样解更像是一道选择题，而在这里我们是不知道答案是什么的，甚至连范围都不知道(你也许会说：x不是已经知道了吗？可是那是拿来验证的啊，不能当做已知条件！)，解决起来更像是一道填空题，答案需要自己推导(想用穷举法？星辰大海，从何开始？)。

求解的方法是最小二乘法，这是多元线性回归里最常用的方法之一，定义可见wiki，我们直接使用结论：x估计解为$\hat {\boldsymbol{x}}=(A^{T}A)^{-1}A^{T}{\boldsymbol {b}}$,，具体的推导过程可见这篇博客，推导过程详尽完整。

当然，Python中也有提供相应的方法可以直接求解最小二乘问题，方法是scipy.optimize.leastsq，方法原型

scipy.optimize.leastsq(func, x0, args=(), Dfun=None, full_output=0, 
                        col_deriv=0, ftol=1.49012e-08, xtol=1.49012e-08,
                        gtol=0.0, maxfev=0, epsfcn=None, factor=100, 
                        diag=None)

一般我们只要指定前三个参数就可以了，在本问题中使用的形式是leastsq(error, x0, (A,b))，error是自定义的估计误差的函数，第二个元素是寻找x的迭代起始点b0，返回值是一个二元组，第一个元素是即是x，第二个元素用不到。完整的用法可见官方文档。

代码

# coding: utf-8
from numpy import *
from scipy.optimize import *
from scipy.linalg import *
#from matplotlib.pyplot import *

A = random.randint(-5, 5, size=(20,10))
x = random.randint(-5, 5, size=10)
b = random.random(size=20)
y = dot(A,x) + b

#方法1：公式法
e_x = dot(dot(inv(dot(A.T, A)), A.T), b)

#方法2：使用现成方法
def error(x_e, A, b):
    #return norm(dot(A, x_e)-b, 2) 错误！只需要做差就行了，不需要把完整的范数写出来
    return dot(A, x_e)-b

x0 = random.random(10) #需要指定一个迭代的初始点
ret = leastsq(error, x0, (A,b))
e_x = ret[0] #取二元组的第一个参数就是解了

print('Result:\nx = ', end = '')
print(e_x)
print('The norm of the residual = ', end = '')
print(norm(error(e_x, A, b), 2))

结果
Result:
x = [ 0.07502006 -0.08488566 -0.00700185  0.00694222  0.00683927  0.01354523
  0.05560228 -0.00382892 -0.0730227  -0.03361303]
The norm of the residual = 1.2709049389807172

最优化问题

分析

这道题的图形化解决已经在上周的练习中展示出来了(点此了解)，本以为这一周的任务轻松一点，只需要求出极值点的数值即可，但实际上这个过程不是很轻松。

对于使用的库，答案是明确的：scipy.optimize，至于方法和思路，我前后尝试了好几种：
1. fmim(了解详情)和fminbound(了解详情)：只能求局部，且需要像在plot时一样声明离散的点，适合画图而不适合这里的精确计算
2. minimize(了解详情)：不能区分全局最小值和局部最小值，初始值选得不好，很容易就落到局部最小值里出不来
3. 在查找方法的时候，还有一种最小值的方法：basinhopping，官方教程在此，这个函数返回的东西有一大串，我们需要提取的是最小值点的列表，这里的用法如下：

#参数是单参函数func，迭代起始点x0，返回列表x的第一个元素就是“解”了
answer = basinhopping(func, x0).x[0] 

可是发现这个函数球的还是局部最小值，在面对不用的起始点的时候结果还是不一样的，因此也摈弃了。

在用符号化求导diff和求解方程solve结合的过程中发现Python的solve挺不实用的，什么东西只要复杂一点就没办法求数值解的零点了，遂放弃。
数学上分析，$e^{x^2}$在$x\rightarrow \pm \infty$的时候都趋近于$0$，而$sin^2(x-2)$是有界的，所以最大值在有穷的区间里，我们可以把区域限制在$[-20,20]$的区间里，用fminbound来求解即可。
以上函数都是求最小值的，那我们要求的是最大值，那么我们如何将求解最小值的问题转化成求解最大值的问题呢？

答案很简单，只需要加个负号就行了。

代码

# coding: utf-8
from numpy import *
from scipy.optimize import *

def fun(x):
    return -((sin(x-2))**2)*exp(-x*x)

x_max = fminbound(fun, -50, 50) #这个函数，撑死不会超过1，且最大值肯定离原点不远

print("The maximum of the function is f(%f) = %f"
    %(x_max, -fun(x_max)))

结果如下
The maximum of the function is f(0.216241) = 0.911685
结果和画图出来的结果一致

点对距离

分析

这道题考察的是关于scipy.spatial里的distance子库，这个库定义的是各种关于几何空间里距离的计算，本题只用到了欧氏距离的计算，对应的范数是2-范数。具体到方法是pdist()方法，这个函数的简单用法是(模式默认是欧氏距离，当然还有其他实用的距离，诸如曼哈顿距离、闵可夫斯基距离等,了解更多)

pdist(X) 

X是一个n行m列的矩阵，一行表示一个m维空间里的点坐标，由n个点的坐标构成了这个矩阵，方法返回的是关于这n个点两两之间的欧式距离计算。
为了使结果更整齐，可以对返回结果的格式进行调整，用squareform()可以将结果转化为n*n距离矩阵。

代码

from numpy import *
from scipy.spatial.distance import *

pair = random.randint(0,50,size = (5,2)) 
print("The coordinates of the 5 points are:")

city_no_ = 1
for i in range(0,5):
    print("City %d : (%d, %d)" % (city_no_, pair[i][0], pair[i][1]))
    city_no_ += 1

print("\nThe distance between any two cities are listed as follow:")
dist = squareform(pdist(pair))
#题目要求对结果进行信息提取
for i in range(0,5):
    for j in range(0,i):
        print("City %d <--> %d : %f" % (i+1,j+1,dist[i][j]))

结果如下
The coordinates of the 5 points are:
City 1 : (3, 33)
City 2 : (44, 47)
City 3 : (9, 30)
City 4 : (20, 38)
City 5 : (32, 45)

The distance between any two cities are listed as follow:
City 2 <--> 1 : 43.324358
City 3 <--> 1 : 6.708204
City 3 <--> 2 : 38.910153
City 4 <--> 1 : 17.720045
City 4 <--> 2 : 25.632011
City 4 <--> 3 : 13.601471
City 5 <--> 1 : 31.384710
City 5 <--> 2 : 12.165525
City 5 <--> 3 : 27.459060
City 5 <--> 4 : 13.892444