“W4 方法”的实现（newton-raphson 扩展）

Question

我目前正在实施一种 newton-raphson 求根方法，该方法可确保在多维环境中收敛（不是作业！）。目前它找到了 x 的根，但没有找到 y。我还观察到一个奇怪的行为，f1 和 f2 被均衡为相同的数字。例如，在 2000 次迭代之后，两者都是 ≈ 560.0。我认为f1 和f2 都需要接近0。至少，这是使用经典的newton-raphson 方法的工作原理。

谁能看出是什么原因造成的？我需要第二双眼睛。

论文：https://arxiv.org/pdf/1809.04495.pdf 和附录：https://arxiv.org/pdf/1809.04358.pdf（D.2 节 -> 包括附加的数学）

注意：U、L 是雅可比矩阵的上下三角矩阵（偏导数矩阵）。

我当前的实现如下所示（使用了 Eigen，但很清楚它的作用）。目前有些奇怪

#include "../../Eigen/Eigen/Core"
#include "../../Eigen/Eigen/LU"
#include <iostream>

int main(){
    double eps = 1e-4;
    Eigen::Vector2d p(0.0, 0.0);

    double x = 0.1;
    double y = 1.0;
    double f1 = 1e9;
    double f2 = 1e9;

    unsigned int count = 0;

    while (count < 2000 && f1 > eps){
        std::cout << "count : " << count << std::endl;

        f1 = x*x - 10*x + y*y - 10*y + 34;
        f2 = x*x - 22*x + y*y - 10*y + 130;

        std::cout << "f1: " << f1 << ", f2: " << f2 << std::endl;

        double A = 2*x - 10;
        double B = 2*y - 10;
        double C = 2*x - 22;
        double D = 2*y - 10;

        Eigen::Matrix2d J;
        J << A, B, C, D;

        Eigen::Matrix2d J_U_inv;
        J_U_inv << J(0,0), J(0,1), 0.0, J(1,1);
        J_U_inv = J_U_inv.inverse();

        Eigen::Matrix2d J_L_inv;
        J_L_inv << J(0,0), 0.0, J(1,0), J(1,1);
        J_L_inv = J_L_inv.inverse();

        Eigen::Vector2d f3(f1, f2);
        Eigen::Vector2d T(x, y);

        if (count == 0){
            p = -0.5 * J_U_inv * f3;
        }

        Eigen::Vector2d E = T + 0.5 * J_L_inv * p;

        p = -0.5 * J_U_inv * f3;

        x = E(0);
        y = E(1);

        std::cout << "x, y: " << x << ", " << y << std::endl;

        ++count;

    }
}

Answer 1

我似乎不知道进行矩阵分解的正确方法。

以下是二维系统的 W4 方法的工作示例。

#include "../../Eigen/Eigen/Core"
#include "../../Eigen/Eigen/LU"
#include <iostream>

int main(){
    double eps = 1e-4;
    Eigen::Vector2d p(0.0, 0.0);

    double x = 0.1;
    double y = 1.0;
    double f1 = 1e9;
    double f2 = 1e9;

    unsigned int count = 0;

    while (std::abs(f1) > eps && std::abs(f2) > eps){
        std::cout << "count : " << count << std::endl;

        f1 = x*x - 10*x + y*y - 10*y + 34;
        f2 = x*x - 22*x + y*y - 10*y + 130;

        std::cout << "f1: " << f1 << ", f2: " << f2 << std::endl;

        double A = 2*x - 10;
        double B = 2*y - 10;
        double C = 2*x - 22;
        double D = 2*y - 10;

        Eigen::Matrix2d J;
        J << A, B, C, D;

        Eigen::Matrix2d J_U_inv;
        J_U_inv << J(0,0) -J(0,1)*J(1,0)/J(1,1),   J(0,1),
                     0.0,                          J(1,1);
        J_U_inv = J_U_inv.inverse().eval();

        Eigen::Matrix2d J_L_inv;
        J_L_inv << 1.0,                             0.0, 
                   J(1,0)/J(1,1),                   1.0;
        J_L_inv = J_L_inv.inverse().eval();


        Eigen::Vector2d f3(f1, f2);
        Eigen::Vector2d T(x, y);

        if (count == 0){
            p = -0.5 * J_U_inv * f3;
        }

        Eigen::Vector2d E = T + 0.5 * J_L_inv * p;

        p = -0.5 * J_U_inv * f3;

        x = E(0);
        y = E(1);

        std::cout << "x, y: " << x << ", " << y << std::endl;

        ++count;

    }
}