[Scipy课后习题]

Scipy课后习题

Exercise 10.1: Least squares

Generate matrix A ∈ Rm×n with m > n. Also generate some vector b ∈ Rm. Now find x = argminxkAx−bk2. Print the norm of the residual.

from scipy import optimize as soe
import numpy as np
import numpy.matlib as npm
import numpy.linalg as npl
import math
m = 15
n = 10
A = np.matlib.rand((m, n))
b = np.matlib.rand((m, 1))
p0 = np.matlib.rand((1, n))
def func(param):
    p = np.matlib.mat(param).T
    return (b - A*p).T.tolist()[0]
r = soe.leastsq(func, p0)
p_com = np.matlib.mat(r[0]).T
print('By call leastsq, norm = ' + str(npl.norm(A*p_com-b)))
x = npl.solve(A.T*A, A.T*b)
print('By slove equation : A.T*A*x = A.T*b, norm = ' 
    + str(npl.norm(A*x-b)))

实验结果：

利用数值方法，解A.T*A*x = A.T*b得到的结果x，计算残差的二范数与之前的残差的二范数进行比较验证，可以看出结果非常接近，结果正确。

Exercise 10.2: Optimization

Find the maximum of the function

f(x) = sin2(x−2)e−x2

方法：通过basinhopping找出-f(x)的最小值，取反即为f(x)最大值

from scipy import optimize as soe
import math

def func(x):
    """
    -f(x)=-sin^2(x-2)e^(-x^2)
    """
    return -math.pow(math.sin(x[0]-2), 2)*math.exp(-math.pow(x[0], 2))

x0 = [0]
minimizer_kwargs = {"method": "BFGS"}
ret = soe.basinhopping(func, x0, minimizer_kwargs=minimizer_kwargs,
    niter=200)
print("global maximum: x = %.4f, f(x0) = %.4f" % (ret.x, -ret.fun))

结果：

Exercise 10.3: Pairwise distances

Let X be a matrix with n rows and m columns. How can you compute the pairwise distances between every two rows?

As an example application, consider n cities, and we are given their coordinates in two columns. Now we want a nice table that tells us for each two cities, how far they are apart.

Again, make sure you make use of Scipy’s functionality instead of writing your own routine.

方法：利用函数pdist

import scipy.spatial.distance as ssdis
import numpy.matlib as npm

X = npm.rand((6, 6))
print('matrix X:')
print(X)
print("matrix X's Pairwise distances:")
print(ssdis.pdist(X, 'euclidean'))

结果：