Solvers¶

In Ikarus, there are essentially two types of solvers: linear and non-linear.

Linear solver¶

It solves for the vector \( \boldsymbol{x} \) in

$$ \boldsymbol{A} \boldsymbol{x} = \boldsymbol{b} $$ where \(\boldsymbol{A} \) and \(\boldsymbol{b}\) are any matrix or vector, respectively. These solvers can be direct or iterative. Furthermore, they depend on the underlying structure of the matrix \(\boldsymbol{A} \), i.e., whether it is stored in a dense or sparse format.

Currently, only the linear solvers provided by the Eigen library are supported.

Linear solvers can be constructed by calling the constructor:

ILinearSolver(const SolverTypeTag& solverTypeTag)

There exists an enum type SolverTypeTag with the following values:

  enum class SolverTypeTag {
    si_ConjugateGradient,
    si_LeastSquaresConjugateGradient,
    si_BiCGSTAB,
    sd_SimplicialLLT,
    sd_SimplicialLDLT,
    sd_SparseLU,
    sd_SparseQR,
    sd_CholmodSupernodalLLT,
    sd_UmfPackLU,
    sd_SuperLU,
    d_PartialPivLU,
    d_FullPivLU,
    d_HouseholderQR,
    d_ColPivHouseholderQR,
    d_FullPivHouseholderQR,
    d_CompleteOrthogonalDecomposition,
    d_LLT,
    d_LDLT,
  };

The prefixes s_ and d_ indicate whether the linear solver can be used for sparse or dense matrices, respectively. Furthermore, there is also a second prefix for sparse solvers: d and i for direct solvers and for iterative solvers, respectively. Thus, using si_ConjugateGradient means that this solver is for sparse matrices and is an iterative solver.

The naming of the solvers is the same as in the Eigen library. For more details on the solvers, we refer to Eigen's documentation for dense decompositions and to Eigen's documentation for sparse decompositions.

Interface¶

Similar to Eigen's interface, the following functions are provided:

void analyzePattern(const MatrixType& A);  // (1)!
void factorize(const MatrixType& A); // (2)!
ILinearSolver& compute(const MatrixType& A); // (3)!
void solve(Eigen::VectorX<ScalarType>&x, const Eigen::VectorX<ScalarType>& b); // (4)!

1. If the matrix is sparse, Eigen can collect information on the sparsity pattern of the matrix for a faster solve step. This pattern does not change if the non-zero entries are modified.
2. This method applies a decomposition technique for direct solvers, e.g., LU decomposition. For iterative solvers, the method is a non-linear operator.
3. It calls both the functions analyzePattern and factorize.
4. Solves the problem and stores the result in x.

Nonlinear solver¶

Non-linear solvers are usually used to solve a kind of optimization problem, e.g., root-finding or minimization problems like: \begin{align} \boldsymbol{R}(\boldsymbol{x}) \stackrel{!}{=} \boldsymbol{0} \quad \text{or} \quad \min_{\boldsymbol{x} \in \mathcal{M}} f(\boldsymbol{x} ) \end{align}

It has the following interface:

void setup(const NewtonRaphsonSettings& p_settings); // (1)!
SolverInformation solve(const SolutionType& dx_predictor = NoPredictor{}); // (2)!
auto& nonLinearOperator(); // (3)!

1. With this function, several properties of the nonlinear solver can be set. E.g., residual tolerance or maximum number of iterations
2. Solves the non-linear problem. An initial guess to the function can be passed, otherwise, a zero vector is assumed. It returns the SolverInformation which contains information like the success of the solution step and more.
3. Just returns the underlying nonLinearOperator, see link.

Implementations¶

To see the Newton-Raphson implementation, we refer to the tests inside nonLinearOperatorTest.cpp, and for the trust-region method, trustRegionTest.cpp.