Non linear solvers#

XLiFE++ provides some tools to solve non linear equations:

\[f(x)=0 \text{ with } f:\mathbb{K}^n\rightarrow \mathbb{K}^m,\quad \mathbb{K}=\mathbb{R} \text{ or } \mathbb{C}.\]

The following solvers are available:

  • Newton method that requires derivatives of \(f\):

    \[x_{k+1}=x_k-J_f(x_k)^{-1}f(x_k).\]

    When \(m\ne n\), \(J_f(x_k)^{-1}\) is replaced by its Moore-Penrose pseudoinverse \(J_f(x_k)^{\dagger}=(J_f(x_k)^*J_f(x_k))^{-1}J_f(x_k)^*\).

  • quasi-Newton (Broyden) method that does not require derivatives of \(f\) but approximates them using finite differences (\(B_k\) an approximation of the Jacobian \(J_f(x_k)\)):

    \[x_{k+1}=x_k-B_k^{-1}f(x_k)\text{ with } B_{k+1}=B_k+\frac{y_k-B_ks_k}{\|xs_k\|^2},\quad y_k=f(x_{k+1})-f(x_k),\ s_k=x_{k+1}-x_k.\]

  • descent method that minimizes \(\frac{1}{2}\|f(x)\|^2\) :

    \[x_{k+1}=x_k-\rho_k d_k \text{ with } d_k=-J_f(x_k)^*f(x_k)\]

    and \(\rho_k\) a step size given either as a constant or by a line search method:

    • Armijo method: if \(f(x_k-\rho_k*d_k)>f(x_k)-\alpha*\rho_k|d_k|^2\) reduce step size \(\rho_k\) by a factor \(\beta\) until the condition is satisfied. In practice, \(\alpha=10^{-4}\) and \(\beta=0.5\) are good values.

    • Barzilai-Borwein \(\rho_k=|s_k|^2/(s_k,y_k)\) where \(s_k=x_{k}-x_{k-1}\), \(y_k=J_f(x_k)^*f(x_k)-J_f(x_{k-1})^*f(x_{k-1})\)

  • Levenberg-Marquardt method that combines the Newton method and the descent method:

    \[x_{k+1}=x_k-(J_f(x_k)^*J_f(x_k)+\lambda I)^{-1}J_f(x_k)^*f(x_k)\]

    where \(\lambda\) is a positive parameter that is reduced when the iteration converges and increased otherwise.

    Newton and quasi-Newton methods are more efficient when the initial guess is close to the solution but may diverge otherwise. Descent and Levenberg-Marquardt methods are more robust but may require more iterations to converge.

    Non linear solvers are implemented in template classes that can be used with any type of scalar (real or complex) and any type of vector classes. The only requirement is that the vector classes must support basic operations such as resize, addition, subtraction, scalar multiplication and dot product. The hierarchy of classes is as follows:

nonlinearsolver_classes

The NonLinearSolverT class is the base class for all non linear solvers. It defines the interface for the solvers and manages main parametrers of the solvers:

template <typename T> class NonLinearSolverT
{
  protected:
   T&(*f_)(const T&, T&, DiffOpType d);    // f(x) : E -> F  (E,F either R, C, R^n, C^n, ...)
   number_t dimIn=1, dimOut=1;             // dim of x, dim of f(x)
  public:
   T x0, x;                      // initial guess, current iterate and solution
   number_t max_iter= 1000;      // maximum of iteration in the itetative process
   number_t k=0;                 // iteration number
   real_t tol=1.E-6;             // tolerance used to stop the itetative process
   bool storeXk = false;         // if true store intermediate xk
   std::ostream* out=nullptr;    // if not nullptr, display intermediate states on *out
   std::list<T> xs;              // list of intermediate xk
   ComputationInfo status=_noConvergence; // final computation status (_success, _noConvergence, ...)
   bool scalar=true;              // true if scalar equation
   string_t name="";              // name of the solver (for display purposes)

In particular, the function \(f\) is a member function pointer that takes a value of type T as input and returns a value of type T. The DiffOpType parameter is an enum that indicates whether the function should return the value of f(x) (when d==_id) or the Jacobian of f at x (when d==_grad).

Inheriting classes manage additional parameters, propose basic constructors, and implement the specific algorithms for each solver.

NewtonSolverT class:

template <typename T> class NewtonSolverT : public NonLinearSolverT<T>
{public:
  NewtonSolverT(T&(*f)(const T&, T&, DiffOpType d), T& x_, number_t nmax, real_t eps, bool stoxk=false);
  T& solve(number_t itmax=0, real_t eps=0);

QuasiNewtonSolverT class:

template <typename T> class QuasiNewtonSolverT : public NonLinearSolverT<T>
{public:
  bool damping = false;  // optionnal damping
  QuasiNewtonSolverT(T&(*f)(const T&, T&, DiffOpType d), T& x_, number_t nmax, real_t eps, bool stoxk=false);
  T& solve(number_t itmax=0, real_t eps=0);

GradientSolverT class:

template <typename T> class GradientSolverT : public NonLinearSolverT<T>
{public:
  real_t alpha=0.;    // constant gradient step
  real_t rho=0.5;     // shrink factor in Armijo rule
  real_t c=1E-4;      // coef in Armijo rule
  StepRule stepRule = _BarzilaiBorweinStepRule;  // step size rule (constant, Armijo, Barzilai-Borwein)
  GradientSolverT(T&(*f)(const T&, T&, DiffOpType d), T& x, number_t nmax, real_t eps, bool stoxk=false);
  T& solve(number_t itmax=0, real_t eps=0);

LevenbergMarquardtSolverT class:

template <typename T> class LevenbergMarquardtSolverT : public NonLinearSolverT<T>
{ public:
  real_t lambda=0.;               // Tikonov regularisation in pseudo inverse of Jacobian
  LevenbergMarquardtSolverT(T&(*f)(const T&, T&, DiffOpType d), T& x, number_t nmax, real_t eps, bool stoxk=false);
  T& solve(number_t itmax=0, real_t eps=0);

Each solver class has a solve method as input the maximum number of iterations and the tolerance for convergence. If process converges,’solve’ returns the solution x that satisfies f(x)=0 within the specified tolerance.

For real scalar case, template solvers classes are aliases with NewtonSolver, QuasiNewtonSolver, GradientSolver and LevenbergMarquardtSolver.