用Python从Matlab代码实现牛顿法反向推导

Question

我正在尝试把我在Matlab里写的牛顿法和回退步骤的代码转到Python上，但在Python的语法上遇到了一些麻烦。Matlab的代码大约需要5次迭代就能完成，但我的Python代码却一直循环到最大迭代次数1000次，而且出现了域错误，因为回退步骤的机制没有正常工作（它试图计算一个负数的对数）。我有一段时间没用Python了，所以我可能搞混了一些语法。

这是在Matlab中正确工作的代码：

x = 10;                                               %defines x
f = @(x) log(x);                                        %defines objective function
df = @(x) 1/x;                                          %defines first derivative
tol = .00001;                                           %defines our tolerance level
maxit = 1000;                                           %defines maximum iteration steps
maxsteps = 200;                                         %defines maximum backsteps
for i=1:maxit                                           %starts loop
    fval = f(x);                                        %value of function at f(x)
    fjac = df(x);                                       %value of jacobian at f(x)
    fnorm = norm(fval);                                 %calculates norm value at fval
    if fnorm<tol, return, end                           %if fnorm less than tol, end
    x
    d = -(fjac\fval);                                   %forms second part of iteration rule
    d
    fnormold = inf;                                     %sets arbitrary fnormold
    for backstep=1:maxsteps
        fvalnew = f(x+d);                               %calculates f(x+d)
        fnormnew = norm(fvalnew);                       %calculates norm of fvalnew
        if fnormnew<fnorm, break, end                   %implements 1st backstepping rule
        if fnormold<fnormnew, d=2*d; break, end         %implements 2nd backstepping rule
        fnormold = fnormnew;                            %updates fnormold
        d=d/2;
    end
    x=x+d;
end
disp(x)

这是我的Python代码：

from math import log

x = 10

def f(x):
    f = x* log(x)
    return f

def df(x):
    df = 1/x
    return df

tol = .00001
maxit = 1000
maxsteps = 200
maxsteps = 200

for i in range(1, maxit):
    fval = f(x)
    fjac = df(x)
    fnorm = abs(fval)
    if fnorm < tol:
        print x
    d = -(fjac/fval)
    fnormold = float('Inf')
    for backstep in range(1, maxsteps):
        fvalnew = f(x+d)
        fnormnew = abs(fvalnew)
        if fnormnew < fnorm:
            break
        if fnormold < fnormnew:
            d= 2*d
            break
        fnormold = fnormnew
        d = d/2
    x = x+d
print x

Answer 1

• 在数据框(df)中，1/x 可能会在大多数情况下变成 0，因为在 Python 2.x 中，整数相除会返回整数结果。
• 基于范围的循环少了一个索引。