Broyden 的方法在 Python 中的实现

Question

我正在尝试用 Python 编写 Broyden 方法的实现。我在Matlab中有一个例子，试图重写它但没有成功。 +A 是雅可比矩阵 F-4雄辩矩阵：

function F=f(X)
 F(1)=X(1)+2*X(2)+X(3)+4*X(4)-20.7;
 F(2)=X(1)^2+2*X(1)*X(2)+X(4)^3-15.88;
 F(3)=X(1)^3+X(3)^2+X(4)-21.218;
 F(4)=3*X(2)+X(3)*X(4)-7.9;
 F=F(:);

for iii=1:itmax 
    deltax=-A\f(x);

    x=x+deltax; fi1=f(x); 
    %A=A+(fi1-fi-A*deltax)*deltax'/(deltax'*deltax);
    a1 = A*deltax;
    a2 = fi1-fi;
    a3 = a2-a1;
    a4 = a3*deltax';
    a5 = deltax'*deltax;
    a6 = a4/a5;
    A=A+a6;

    tikslumas=norm(deltax)/(norm(x)+norm(deltax));

    if tikslumas < eps
        fprintf(1,'\n x ='); fprintf(1,'  %g',x);
        fprintf(1,'\n  f ='); fprintf(1,'  %g',f(x));
        break
    elseif iii == itmax
        fprintf(1,'\n ****failed: x ='); fprintf(1,'  %g',x);
        fprintf(1,'\n  f ='); fprintf(1,'  %g',f(x));
        break
fi=fi1;
end

我正在尝试在 Python 中做同样的事情：

for i in range(1, itmax):
    deltax,resid,rank,s = np.linalg.lstsq(A, f(x))
    for j in range(len(deltax)):
        deltax[j] = -deltax[j]
    for j in range(len(deltax)):
        x[j] = x[j] + deltax[j]        
    fi1=f(x);
    deltaxT = np.matrix(deltax).T
    a1 = np.multiply(A, deltax)
    a1 = a1.sum(axis=1)
    a2 = np.subtract(fi1,fi)
    a3 = np.subtract(a2,a1)
    a4 = np.multiply(a3,deltaxT)
    a4 = np.matrix(a4).T
    a5 = (np.multiply(deltax,deltax)).sum(axis=0) #???? gal???
    a6 = np.divide(a4,a5)
    A=A+a6;

但是3-4个for循环后的结果不一样，最终结果完全不同

Answer 1

主要错误

您的主要错误是混淆了逐点乘法和矩阵乘法。您编写的 Broyden 算法的一个步骤可以如下所示：

def broyden_step(f, fi, A, x):
    # A = A +(fi1-fi-A @ deltax) @ deltax'/(deltax, deltax);
    deltax = -np.linalg.solve(A, fi)
    x = x + deltax
    fi1 = f(x)
    a1 = np.dot(A, deltax)
    a2 = fi1 - fi
    a3 = a2 - a1
    a4 = np.outer(a3, deltax)
    a5 = np.dot(deltax, deltax)
    a6 = a4/a5
    A = A + a6
    return fi1, A, x, deltax

算法生成器

实现这种算法的Pythonic解决方案是使用生成器。这是一个例子：

def broyden(f, A, x0, itmax=100):
    x = x0
    fi = f(x)
    for i in range(itmax):
        fi, A, x, deltax = broyden_step(f, fi, A, x)
        crit = np.linalg.norm(deltax)/(np.linalg.norm(x) + np.linalg.norm(deltax))
        yield crit, x

您可以使用

data = list(broyden(f, A, x0, itmax=200))

运行它。 您应该将其包装在一个循环中，该循环根据

crit

决定何时停止算法。

额外：具有自动微分功能的雅可比行列式

为此我将使用 JAX。

import jax
import jax.numpy as np

这是示例中的非线性函数：

def f(x):
    zero = np.zeros(4)
    res = (
        zero
        .at[0].set(x[0] + 2*x[1] + x[2] + 4*x[3])
        .at[1].set(x[0]**2 + 2*x[0]*x[1] + x[3]**3)
        .at[2].set(x[0]**3 + x[2]**2 + x[3])
        .at[3].set(3*x[1] + x[2]*x[3])
    )
    shift = np.array([20.7, 15.88, 21.218, 7.9])
    return res - shift

现在雅可比矩阵可以用

计算

J = jax.jacfwd(f)

您可以在算法中使用它，例如

A = J(x0)

x0 = np.ones(4)
broyden(f, J(x0), x0, itmax=200)