牛頓法python

Ⅰ 平面內一點到另兩點距離之和最小的求法

怎樣「求空間內一點到其它所有點的距離之和友塵嘩最小」？首先將這個問題形式化：
公式代碼：
\min f(x,y) = \min \sum_i \sqrt {(x - x_i)^2 + (y - y_i)^2}
這里是距離之和，而不是平方和。Kmeans聚類中用的評價標準是平方和，如果只有一個類中心，那麼可以直接求偏導得到使得平方和最小的點就是中心。這里問題與平方和的解是不是一樣的，比如三角形到三個頂點距離之和最短的點就是費馬點。
這里可以用最優化方法中的「搜索」來求解，這一系列方法包括了梯度下降法、牛頓法和共軛梯度法等。在這里用梯度下降法是最簡單的，通過這個例子我也明白了為什麼實際運用中梯度下降法是應用最廣的。相比梯度下降法，牛頓法需要求Hesse矩陣，還是相對麻 煩不少。梯度下降法搜索步驟就是每一步都向導數的逆方向將自變數前進一個步長（可變），在這里導數方向就是

公式代碼：
abla f(x,y) =
\left[
\begin{array} {lcr}
\displaystyle \sum_i \frac{x - x_i}{\sqrt{(x - x_i)^2 + (y - y_i)^2pan >}} \\
\displaystyle \sum_i \frac{y - y_i}{\sqrt{(x - x_i)^2 + (y - y_i)^2}}
\end{array}
\right]
梯度下兄裂降法也有它使用起來讓人比較為難的地方，那就是步長很難選取，課本上所給出的例子一般都是針對較簡單表達式提出的可變步長計算。在本問題的求解中為簡單起見，步長是取的定值。整個過程用Python3實現（起初想用R來做，但是R沒法調試……歸根結底還是功力不夠）實現，結合了scipy和matplotlib兩個好行包，結果看起來還是比較靠譜：

最後附上源代碼：
Python 3語言: 高亮代碼由發芽網提供
from scipy import *
import pylab
def f(p, pts):
return sum(sum((p - pts) ** 2, axis=1) ** 0.5)
def fd(p, pts):
dx = sum((p[0] - pts[:, 0]) / sum((p - pts) ** 2, axis=1) ** 0.5)
dy = sum((p[1] - pts[:, 1]) / sum((p - pts) ** 2, axis=1) ** 0.5)
s = (dx ** 2 + dy ** 2) ** 0.5
br> dx /= s
dy /= s
return array([dx, dy])

pts = rand(10, 2)
x = array([0, 0])
t = 0.1
xstep = x
for k in range(100):
y = f(x, pts)
xk = x - t * fd(x, pts)
yk = f(xk, pts)
if y - yk > 1e-8:
x = xk
y = yk
elif yk - y > 1e-8:
t *= 0.5
else:
break
xstep = vstack((xstep, x))
print(x, y)
pylab.plot(pts[:, 0], pts[:, 1], 'bo')
pylab.plot(xstep[:, 0], xstep[:, 1], 'ro')
pylab.plot(xstep[:, 0], xstep[:, 1], 'k-')
pylab.xlabel('iter = %d, Min = %.3f, p = (%.3f, %.3f), t = %f' % (k, y, x[0], x[1], t))
pylab.show()