tr_solver.hpp File Reference

Trust region algorithms for determining step size. More...

#include <optional>
#include <utility>
#include <Eigen/Cholesky>
#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
#include "smooth/version.hpp"

Include dependency graph for tr_solver.hpp:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
template<typename D2 , typename D3 >
SMOOTH_BEGIN_NAMESPACE auto	solve_linear_ldlt (const auto &J, const Eigen::MatrixBase< D2 > &d, const Eigen::MatrixBase< D3 > &r, const double lambda, std::optional< std::reference_wrapper< double > > dphi={}) -> Eigen::Vector< typename std::decay_t< decltype(J)>::Scalar, std::decay_t< decltype(J)>::ColsAtCompileTime >
	Solve least-squares problem using LDL' decomposition.
template<typename D2 , typename D3 >
auto	solve_trust_region (const auto &J, const Eigen::MatrixBase< D2 > &d, const Eigen::MatrixBase< D3 > &r, const double Delta) -> std::pair< Eigen::Vector< typename std::decay_t< decltype(J)>::Scalar, std::decay_t< decltype(J)>::ColsAtCompileTime >, double >
	Approximately solve trust-region step determination problem.

Detailed Description

Trust region algorithms for determining step size.

Definition in file tr_solver.hpp.

Function Documentation

solve_linear_ldlt()

template<typename D2 , typename D3 >

SMOOTH_BEGIN_NAMESPACE auto solve_linear_ldlt	(	const auto &	J,
		const Eigen::MatrixBase< D2 > &	d,
		const Eigen::MatrixBase< D3 > &	r,
		const double	lambda,
		std::optional< std::reference_wrapper< double > >	dphi = {}
	)		-> Eigen::Vector<typename std::decay_t<decltype(J)>::Scalar, std::decay_t<decltype(J)>::ColsAtCompileTime>

Solve least-squares problem using LDL' decomposition.

The problem is

[    J      ] dx  = [ -r ]                            (1)
[ sqrt(λ) D ]       [  0 ]

Parameters

[in]	J	matrix of size M x N
[in]	d	positive vector of size N representing diagonal of D
[in]	r	vector of size M
[in]	lambda	non-zero number
[out]	dphi	optionally calculate derivative ϕ'(λ)

ϕ is the norm of the scaled solution as a function of λ:

ϕ(λ) = || D dx || = || D (J' J + λ D'D)^{-1} J' r ||

The system (1) is solved via LDLt factorization of the left-hand side of the normal equations

(J' J + λ D' D) dx = -J' r.

Definition at line 48 of file tr_solver.hpp.

solve_trust_region()

template<typename D2 , typename D3 >

auto solve_trust_region	(	const auto &	J,
		const Eigen::MatrixBase< D2 > &	d,
		const Eigen::MatrixBase< D3 > &	r,
		const double	Delta
	)		-> std:: pair<Eigen::Vector<typename std::decay_t<decltype(J)>::Scalar, std::decay_t<decltype(J)>::ColsAtCompileTime>, double>

Approximately solve trust-region step determination problem.

The step determination problem is to find the minimizing dx of

min 0.5 || J dx  + r ||^2  s.t. || D dx || ≤ Δ          (1)

Parameters

J	matrix of size M x N
d	vector of size N representing diagonal of D
r	vector of size M
Delta	trust region size

Returns

{dx, λ} where dx is the (approximate) minimizer of (1) and λ is the corresponding Lagrange multiplier for the inequality constraints.

Note

The current implementation sets the Lagrange multiplier to λ = 1 / Δ.

Theory

It can be shown that (1) is equivalent to

min 0.5 ( || J dx + r ||^2 + λ || D dx ||^2 )           (2)

for some λ that satisfies the complimentarity condition λ (|| D dx || - Δ) = 0.

in turn, (2) is equivalent to

min || [    J      ] dx  + [ r ] ||^2                   (3)
    || [ sqrt(λ) D ]       [ 0 ] ||

which is the least-squares solution to

[    J      ] dx  = [ -r ]                              (4)
[ sqrt(λ) D ]       [  0 ]

with normal equations

(J' J + λ D' D) dx = - J' r                             (5)

Definition at line 132 of file tr_solver.hpp.