提升库。雅可比和系统设置-解网

问：

假设我有一个方程组：

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

我想用牛顿的迭代方法求解这个系统：其中 J（x） - 这个系统的雅可比。J(x) * delta_x = -f(x)

问题是：计算雅可比和设置系统的最佳方法是什么？我有一个工作版本，但是在其中设置方程式和系统很尴尬：

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

C++ 数学提升

在指定系统的方程时，您已经必须做很多与问题相关的工作。我很欣赏你可以有很多变量，但如果你能用代码写出各个方程，那么你就可以对雅可比式做同样的事情。一个更大的问题可能是牛顿-拉夫森可能会失败，以及设定一个好的初始猜测的挑战。例如，您的系统示例有两个解决方案（1/sqrt（2），1/sqrt（2））和（-1/sqrt（2），-1/sqrt（2）），但是如果您从初始猜测（0,0）开始，那么它将在第一步失败。

0赞 Robert Dodier 11/10/2023

有趣的问题;你能澄清一下你的任务是什么吗？我的意思是，目标是求解方程，还是目标是使用 Boost 求解方程？如果您不必使用 Boost，还有其他（更简单的）方法可以实现 Newton 方法，但也许您需要使用 Boost。

0赞 sehe 11/10/2023

@RobertDodier 我阅读它的方式，可能不是关于解决特定示例，而是如何以一种不那么笨拙/乏味的方式使用 Boost 接口。

答：

3赞 sehe 11/9/2023 #1

虽然你没有真正解释代码的“尴尬”之处。

我看了一会儿，能够提出一些优化和简化（例如，使用 C++17 语言功能）。

请注意，为了清楚起见（和安全起见），我使用命名空间别名：
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>>;

我之所以选择它，是因为它简单且对 constexpr 友好，并且（与相反）允许初始值设定项列表构造。这将允许我们进行更自然的任务：std::arrayU::bounded_array<>

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),
        }};
}

另一部分是使用结构化绑定而不是 or dances。get<n>[n]

从效率的角度来看，我会验证编译器是否消除了矩阵初始值设定项中的常见子表达式。

对于 / 函数，我让类型推导来完成工作：f0f1

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

主驱动程序没有太大变化，尽管我重新塑造了循环，使其更加传统：

在 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;
}

印刷

[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/O 的微基准测试会很有趣）。

我不确定这是否解决了您在问题中提到的一些“尴尬”。

我不喜欢函数 f0、f1 的单独声明。我想将案例扩展到大量变量，而代码中的更改却很少。例如，将函数 f 设为数组，并设置例如 20 个函数。接下来，雅可比将自动计算为 20 个具有 20 个变量的函数（2 个令人不快的地方是“提升”会很大）。但我真的很喜欢你的代码改进，你使用了我以前从未见过的构造。我非常感谢你向他们展示=）

2赞 lastchance 11/10/2023 #2

我想你可以避免提升，只滚动你自己的。

使用泛函设置的方程。雅各派通过有限的差异来制定普遍性，但可以进行改进。

#include <iostream>
#include <vector>
#include <functional>
#include <utility>
#include <cmath>
using namespace std;

using FUNC = function< double( vector<double> ) >;         // double-valued function with vector argument

bool solve( vector<FUNC> eqns, vector<double> &x );
vector<double> operator -( vector<double> a, const vector<double> &b );
vector<double> getFunction( vector<FUNC> eqns, const vector<double> &x );
vector< vector<double> > getJacobian( vector<FUNC> eqns, const vector<double> &x );
double getResiduals( const vector<double> &f );
bool GaussElimination( vector< vector<double> > A, vector<double> B , vector<double> &X );

//======================================================================

int main()
{
   vector<FUNC> eqns = {                                                  // equations to be solved
      []( vector<double> x ){ return x[0] * x[0] + x[1] * x[1] - 1; },
      []( vector<double> x ){ return x[0] - x[1];                   }
                       };
   vector<double> x = { 1.0, 1.0 };                                       // initial guess (gives positive root)
// vector<double> x = {-1.0,-1.0 };                                       // alternative (gives negative root)

   if ( solve( eqns, x ) )
   {
      for ( auto e : x ) cout << e << '\n';
   }
   else
   {
      cout << "Failed\n";
   }
}

//======================================================================

bool solve( vector<FUNC> eqns, vector<double> &x )
{
   const int MAXITER = 1000;
   const double TOL = 1.0e-6;
   double residuals = 1.0;
   for ( int iter = 1; iter <= MAXITER && residuals > TOL; iter++ )
   {
      auto f = getFunction( eqns, x );                       // compute f(x*)
      if ( residuals < TOL ) break;                          // if converged, finish

      auto J = getJacobian( eqns, x );                       // computed J(x*)
      vector<double> deltax( x.size() );
      if ( !GaussElimination( J, f, deltax ) ) break;        // update x -> x* - J^-1.f(x*)
      x = x - deltax;                                        //

      residuals = getResiduals( f );                         // test for convergence
   }
   return residuals <= TOL;
}

//======================================================================

vector<double> operator -( vector<double> a, const vector<double> &b )         // utility routine to subtract vectors
{
   for ( int i = 0; i < b.size(); i++ ) a[i] -= b[i];
   return a;
}

//======================================================================

vector<double> getFunction( vector<FUNC> eqns, const vector<double> &x )       // fill in current function values
{
   vector<double> f;
   for ( FUNC e : eqns ) f.push_back( e( x ) );
   return f;
}

//======================================================================

vector< vector<double> > getJacobian( vector<FUNC> eqns, const vector<double> &x )   // create Jacobian from FINITE DIFFERENCES
{
   const double SMALL = 1.0e-6;                                                      // needs improving for generality
   int N = x.size();
   vector< vector<double> > J( N, vector<double>( N ) );
   for ( int i = 0; i < N; i++ )
   {
      for ( int j = 0; j < N; j++ ) 
      {
         auto xplus  = x;   xplus [j] = x[j] + 0.5 * SMALL;
         auto xminus = x;   xminus[j] = x[j] - 0.5 * SMALL;
         J[i][j] = ( eqns[i]( xplus ) - eqns[i]( xminus ) ) / SMALL;                 // finite-difference approximation for df_i/dx_j
      }
   }
   return J;
}

//======================================================================

double getResiduals( const vector<double> &x )                                       // L2 norm of vector x
{
   double res = 0.0;
   for ( auto e : x ) res += e * e;
   return sqrt( res );
}

//======================================================================

bool GaussElimination( vector< vector<double> > A, vector<double> B , vector<double> &X )
//-------------------------------------------------------------------
// Solve AX = B by Gaussian elimination
// Note: copies made of A and B
//-------------------------------------------------------------------
{
   const double SMALL = 1.0e-20;
   int n = A.size();

   // Row operations for i = 0, ,,,, n - 2 (n-1 not needed)
   for ( int i = 0; i < n - 1; i++ )
   {
      // Pivot: find row r below with largest element in column i
      int r = i;
      double maxA = abs( A[i][i] );
      for ( int k = i + 1; k < n; k++ )
      {
         double val = abs( A[k][i] );
         if ( val > maxA )
         {
            r = k;
            maxA = val;
         }
      }
      if ( r != i )
      {
         for ( int j = i; j < n; j++ ) swap( A[i][j], A[r][j] );
         swap( B[i], B[r] );
      }

      // Row operations to make upper-triangular
      double pivot = A[i][i];
      if ( abs( pivot ) < SMALL ) return false;            // Singular matrix

      for ( int r = i + 1; r < n; r++ )                    // On lower rows
      {
         double multiple = A[r][i] / pivot;                // Multiple of row i to clear element in ith column
         for ( int j = i; j < n; j++ ) A[r][j] -= multiple * A[i][j];
         B[r] -= multiple * B[i];
      }
   }

   // Check last row
   if ( abs( A[n-1][n-1] ) < SMALL ) return false;         // Singular matrix

   // Back-substitute
   X = B;
   for ( int i = n - 1; i >= 0; i-- )
   {
      for ( int j = i + 1; j < n; j++ )  X[i] -= A[i][j] * X[j];
      X[i] /= A[i][i];
   }

   return true;
}

输出：

0.707107
0.707107

谢谢你，@sehe。正如您所说，矩阵求解器中的 SMALL 可能太小了（！），但无论如何都会使用库例程。雅可比例程中的 SMALL 更麻烦：如果变量范围是 O（1）就可以了，但如果变量范围是 O（1e-10）或 O（1e10），则不行。我从来没有找到一个可靠的答案（除了更喜欢使用无量纲变量，这些变量在构造上是 O（1））。

0赞 sehe 11/11/2023

据我所知，获得“域相对 epsilon”的最佳方法是从 std：：nextafter[f] 和朋友推断。事实上，使用某种任意/高精度的实数类型会非常难降低性能：添加了cpp_bin_float比较图表的代码 i.imgur.com/uRqSlAH.png

上一个：浮点数学坏了吗？

下一个：为什么我不能在 C/C++ 中用常量替换变量？

提升库。雅可比和系统设置

Boost library. Jacobian and system setup

评论

总结

评论

评论