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;
}
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:
#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)
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.