Implementation of "The W4 Method" (newton-raphson extension)

Question

I am currently implementing a newton-raphson root-finding method that guarantees convergence in a multidimensional setting (not homework!). Currently it finds the root for x, but not for y. I also observed a strange behaviour where f1 and f2 get equalised to the same number. For example, after 2000 iterations both are ≈ 560.0. I think f1 and f2 both need to approach 0. At least, this is how it works using the classical newton-raphson method.

Can anyone see what could cause this? I need a second pair of eyes.

paper: https://arxiv.org/pdf/1809.04495.pdf and addendum: https://arxiv.org/pdf/1809.04358.pdf (section D.2 -> includes the attached math)

Note: U, L are the upper and lower triangular matrices of the Jacobian (matrix of partial derivatives).

My current implementation looks like the following (Eigen is used, but it's pretty clear what it does). Currently something strange

#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;

    }
}

zeno · Accepted Answer · 2020-03-14T14:12:42.157

It seems I was not aware of the proper way to do the matrix decomposition.

Below is a working example of the W4 method for a 2dimensional system.

#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;

    }
}

Implementation of "The W4 Method" (newton-raphson extension)

1 Answers1