维基百科示例中高斯-牛顿法的实现

Question

我对 Python 比较陌生，正在尝试实现 Gauss-Newton 方法，特别是 Wikipedia 页面上的示例（Gauss–Newton algorithm，3 示例）。以下是我到目前为止所做的：

import scipy
import numpy as np
import math
import scipy.misc

from matplotlib import pyplot as plt, cm, colors

S = [0.038,0.194,.425,.626,1.253,2.500,3.740]
rate = [0.050,0.127,0.094,0.2122,0.2729,0.2665,0.3317]
iterations = 5
rows = 7
cols = 2

B = np.matrix([[.9],[.2]]) # original guess for B

Jf = np.zeros((rows,cols)) # Jacobian matrix from r
r = np.zeros((rows,1)) #r equations


def model(Vmax, Km, Sval):
   return ((vmax * Sval) / (Km + Sval))

def partialDerB1(B2,xi):
   return round(-(xi/(B2+xi)),10)

def partialDerB2(B1,B2,xi):
   return round(((B1*xi)/((B2+xi)*(B2+xi))),10)

def residual(x,y,B1,B2):
   return (y - ((B1*x)/(B2+x)))


for i in range(0,iterations):

   sumOfResid=0
   #calculate Jr and r for this iteration.
   for j in range(0,rows):
      r[j,0] = residual(S[j],rate[j],B[0],B[1])
      sumOfResid = sumOfResid + (r[j,0] * r[j,0])
      Jf[j,0] = partialDerB1(B[1],S[j])
      Jf[j,1] = partialDerB2(B[0],B[1],S[j])

   Jft =  np.transpose(Jf)
   B = B + np.dot((np.dot(Jft,Jf)**-1),(np.dot(Jft,r)))

   print B

在每次迭代中，残差的平方和都会增加而不是趋向于 0，我得到的 B 向量也会增加。

我无法理解我的问题出在哪里，我们将不胜感激。

Answer 1

你在测试版更新的代码中出错了：应该是

B = B - np.dot(np.dot( inv(np.dot(Jft, Jf)), Jft), r)

而不是**-1在矩阵上计算逆矩阵

import scipy
import numpy as np
from numpy.linalg import inv
import math
import scipy.misc

#from matplotlib import pyplot as plt, cm, colors

S = [0.038,0.194,.425,.626,1.253,2.500,3.740]
rate = [0.050,0.127,0.094,0.2122,0.2729,0.2665,0.3317]
iterations = 5
rows = 7
cols = 2

B = np.matrix([[.9],[.2]]) # original guess for B
print(B)

Jf = np.zeros((rows,cols)) # Jacobian matrix from r
r = np.zeros((rows,1)) #r equations


def model(Vmax, Km, Sval):
   return ((Vmax * Sval) / (Km + Sval))

def partialDerB1(B2,xi):
   return round(-(xi/(B2+xi)),10)

def partialDerB2(B1,B2,xi):
   return round(((B1*xi)/((B2+xi)*(B2+xi))),10)

def residual(x,y,B1,B2):
   return (y - ((B1*x)/(B2+x)))

#
for _ in xrange(iterations):

   sumOfResid=0
   #calculate Jr and r for this iteration.
   for j in xrange(rows):
      r[j,0] = residual(S[j],rate[j],B[0],B[1])
      sumOfResid += (r[j,0] * r[j,0])
      Jf[j,0] = partialDerB1(B[1],S[j])
      Jf[j,1] = partialDerB2(B[0],B[1],S[j])

   Jft =  Jf.T
   B -= np.dot(np.dot( inv(np.dot(Jft,Jf)),Jft),r)

   print B