Gauss-Newton (via normal eqs)#

Header: <pressio/solvers_nonlinear_gaussnewton.hpp>

API#

namespace pressio{

template<class SystemType, class LinearSolverType>
#ifdef PRESSIO_ENABLE_CXX20
  requires nonlinearsolvers::RealValuedNonlinearSystemFusingResidualAndJacobian<SystemType>
  && nonlinearsolvers::valid_state_for_least_squares<typename SystemType::state_type>::value
  && (Traits<typename SystemType::state_type>::rank    == 1)
  && (Traits<typename SystemType::residual_type>::rank == 1)
  && (Traits<typename SystemType::jacobian_type>::rank == 2)
  && requires(
	    typename SystemType::state_type    & x,
      const typename SystemType::jacobian_type & J,
      const typename SystemType::residual_type & r,
      nonlinearsolvers::normal_eqs_default_hessian_t<typename SystemType::state_type>  & H,
      nonlinearsolvers::normal_eqs_default_gradient_t<typename SystemType::state_type> & g,
      LinearSolverType && linSolver)
  {
    { ::pressio::ops::norm2(r) } -> std::same_as< nonlinearsolvers::scalar_of_t<SystemType> >;
    { ::pressio::ops::product(transpose(), nontranspose(), 1, J, 0, H) };
    { ::pressio::ops::product(transpose(), 1, J, r, 0, g) };
    { linSolver.solve(std::as_const(H), std::as_const(g), x) };
  }
#endif
auto create_gauss_newton_solver(const SystemType & system,           /*(1)*/
				LinearSolverType && linSolver)
  where: delta = x_k+1 - x_k, H = J^T_r* W * J_r, and g = J^T_r * W * r
  where: delta = x_k+1 - x_k, H = J^T_r* W * J_r, and g = J^T_r * W * r
*/

template<class SystemType, class LinearSolverType, class WeightingOpType>
#ifdef PRESSIO_ENABLE_CXX20
  requires nonlinearsolvers::RealValuedNonlinearSystemFusingResidualAndJacobian<SystemType>
  && nonlinearsolvers::valid_state_for_least_squares<typename SystemType::state_type>::value
  && (Traits<typename SystemType::state_type>::rank    == 1)
  && (Traits<typename SystemType::residual_type>::rank == 1)
  && (Traits<typename SystemType::jacobian_type>::rank == 2)
  && requires(
	     typename SystemType::state_type    & x,
       const typename SystemType::residual_type & r,
       const typename SystemType::jacobian_type & J,
       typename SystemType::residual_type & Wr,
       typename SystemType::jacobian_type & WJ,
       nonlinearsolvers::normal_eqs_default_hessian_t<typename SystemType::state_type>  & H,
       nonlinearsolvers::normal_eqs_default_gradient_t<typename SystemType::state_type> & g,
       WeightingOpType && W,
       LinearSolverType && linSolver)
  {
    { ::pressio::ops::norm2(std::as_const(r)) }
      -> std::same_as< nonlinearsolvers::scalar_of_t<SystemType> >;
    { ::pressio::ops::dot(std::as_const(r), std::as_const(Wr)) }
      -> std::same_as< nonlinearsolvers::scalar_of_t<SystemType> >;

    { W(r, Wr) };
    { W(J, WJ) };
    { ::pressio::ops::product(transpose(), nontranspose(), 1, J, std::as_const(WJ), 0, H) };
    { ::pressio::ops::product(transpose(), 1, J, std::as_const(Wr), 0, g) };

    { linSolver.solve(std::as_const(H), std::as_const(g), x) };
  }
#endif
auto create_gauss_newton_solver(const SystemType & system,           /*(2)*/
  using reg_t    = nonlinearsolvers::impl::RegistryGaussNewtonQr<SystemType, QRSolverType>;
  using scalar_t = nonlinearsolvers::scalar_of_t<SystemType>;

Parameters#


`system`	your problem instance
`linSolver`	linear solver to solve the normal equations at each nonlinear iteration
`weigher`	the weighting operator for doing weighted least-squares

Constraints#

Concepts are documented here. Note: constraints are enforced via proper C++20 concepts when PRESSIO_ENABLE_CXX20 is enabled, otherwise via SFINAE and static asserts.

Examples#

Demos

full demo