Boost library. Jacobian and system setup

Question

Let's assume I have a system of equations:

x^2 + y^2 = 1
x - y = 0

I want to solve this system by Newton's iterative method: J(x) * delta_x = -f(x) Where J(x) - the Jacobian of this system.

The question is:What is the best way to calculate the Jacobian and set the system? I have a working version, but it is awkward to set equations and a system in it:

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/math/differentiation/autodiff.hpp>

using namespace boost::numeric::ublas;
using namespace boost::math::differentiation;

template<typename X, typename Y>
promote<X, Y> f0(const X &x, const Y &y) {
    return x * x + y + y - 1;
}

template<typename X, typename Y>
promote<X, Y> f1(const X &x, const Y &y) {
    return x - y;
}

vector<double> f(const vector<double> &var) {
    vector<double> f(2, 0);
    f[0] = f0(var[0], var[1]);
    f[1] = f1(var[0], var[1]);
    return f;
}

matrix<double> J(const vector<double> &var) {
    matrix<double> J(2, 2, 0);
    
    auto const variables = make_ftuple<double, 1, 1>(var[0], var[1]);
    
    auto const &x = std::get<0>(variables);
    auto const &y = std::get<1>(variables);
    
    J(0, 0) = f0(x, y).derivative(1, 0);
    J(0, 1) = f0(x, y).derivative(0, 1);
    J(1, 0) = f1(x, y).derivative(1, 0);
    J(1, 1) = f1(x, y).derivative(0, 1);
    
    return J;
}

int main() {
    vector<double> x(2, 0);
    x[0] = 0.5;
    x[1] = 0.5;
    double eps = 1.e-6;
    int iter = 0;
    
    vector<double> fx, delta_x;
    matrix<double> Jx;
    
    while (iter < 1.e3) {
        Jx = J(x);
        fx = f(x);
        delta_x = -fx;
        // Solve the linear system J(x) * delta_x = -f(x)
        
        permutation_matrix<std::size_t> pm(Jx.size1());
        lu_factorize(Jx, pm);
        lu_substitute(Jx, pm, delta_x);
        
        x += delta_x;
        if (norm_2(delta_x) < eps) break;
        
        if (++iter % 10 == 0)
            std::cout << "[" << iter << "] " << "Solution: " << x << std::endl;
    }
    std::cout << "Final answer:" << std::endl;
    std::cout << "[" << iter << "] " << "Solution: " << x << std::endl;
    return 0;
}

Answer 1

Although you didn't really explain what is "awkward" about the code.

I looked at it for a bit and was able to come up with some optimizations and simplifications (e.g. using C++17 language features).

Note, for clarity (and safety) I use namespace aliases:

namespace U = boost::numeric::ublas;
namespace D = boost::math::differentiation;

The optimizations is to avoid unbounded array storage for vectors/matrices because they are fixed size:

template <size_t N, typename T = double> using FixedVector = U::vector<T, std::array<T, N>>;
using Vec = FixedVector<2>;
using Mat = U::matrix<double, U::row_major, std::array<double, 4>>;

I opted for std::array because it's trivial and constexpr friendly, and (contrary to U::bounded_array<>) allows initializer-list construction. This will allow us more natural assignments:

Vec f(Vec const& var) {
    auto& [x, y] = var.data();
    return Vec(2, {f0(x, y), f1(x, y)});
}

Mat J(Vec const& var) {
    auto const& [x, y] = D::make_ftuple<double, 1, 1>(var[0], var[1]);
    return Mat{2, 2,
        {
            f0(x, y).derivative(1, 0),
            f0(x, y).derivative(0, 1),
            f1(x, y).derivative(1, 0),
            f1(x, y).derivative(0, 1),
        }};
}

The other part is to use structed bindings instead of the get<n> or [n] dances.

Looking at this from an efficiency standpoint, I'd verify that the compiler eliminates the common subexpressions in the matrix initializers.

For the f0/f1 functions I've let type deduction do the work:

constexpr auto f0(auto const& x, auto const& y) { return x * x + y + y - 1; }
constexpr auto f1(auto const& x, auto const& y) { return x - y; }

The main driver hasn't changed much, though I reshaped the loop so it's more conventional:

Live On Coliru

#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/operations.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/math/differentiation/autodiff.hpp>

constexpr auto f0(auto const& x, auto const& y) { return x * x + y + y - 1; }
constexpr auto f1(auto const& x, auto const& y) { return x - y; }

namespace U = boost::numeric::ublas;
namespace D = boost::math::differentiation;

template <size_t N, typename T = double> using FixedVector = U::vector<T, std::array<T, N>>;
using Vec = FixedVector<2>;
using Mat = U::matrix<double, U::row_major, std::array<double, 4>>;

Vec f(Vec const& var) {
    auto& [x, y] = var.data();
    return Vec(2, {f0(x, y), f1(x, y)});
}

Mat J(Vec const& var) {
    auto const& [x, y] = D::make_ftuple<double, 1, 1>(var[0], var[1]);
    return Mat{2, 2,
        {
            f0(x, y).derivative(1, 0),
            f0(x, y).derivative(0, 1),
            f1(x, y).derivative(1, 0),
            f1(x, y).derivative(0, 1),
        }};
}

int main() {
    constexpr double eps  = 1.e-6;
    Vec x(2, {0.5, 0.5});

    U::vector<double> fx, delta_x;
    U::matrix<double> Jx;

    for (auto iteration = 0; iteration < 1'000; ++iteration) {
        Jx      = J(x);
        fx      = f(x);
        delta_x = -fx;
        // Solve the linear system J(x) * delta_x = -f(x)

        U::permutation_matrix<std::size_t> pm(Jx.size1());
        lu_factorize(Jx, pm);
        lu_substitute(Jx, pm, delta_x);
        
        x += delta_x;
        std::cout << "[" << iteration << "] Solution: " << x << std::endl;

        if (norm_2(delta_x) < eps)
            break;
    }
    std::cout << "Solution: " << x << std::endl;
}

Printing

[0] Solution: [2](0.416667,0.416667)
[1] Solution: [2](0.414216,0.414216)
[2] Solution: [2](0.414214,0.414214)
[3] Solution: [2](0.414214,0.414214)
Solution: [2](0.414214,0.414214)

Summary

I expect the changes to make the code marginally faster (it'd be interesting to see a micro-benchmark without the I/O).

I'm not sure whether any of this addressed some of the "awkwardness" that you were referring to in the question.