Moore-Penrose continuation

• Moore-Penrose continuation
• • Predictor
  • Corrector
  • Linear Algebra
  • • Norm
    • Linear problem
  • Algorithm for solving (MS)
  • • Direct case
    • Iterative case
  • Step size control
  • References

This is one of the various continuation methods implemented in BifurcationKit.jl. It is set by the option alg = MoorePenrose() in continuation. See also MoorePenrose for more information.

For solving

\[\mathbb R^n\ni F(x,p) = 0 \quad\tag{E}\]

using a Newton algorithm, we miss an equation. Hence, we proceed as follows ^[Meijer]. Starting from a predictor $(x_1,p_1)$, we look for the solution to (E) that is closest to $(x_1,p_1)$. Hence, we optimise

\[\min_{(x,p)} \{ \|(x,p)-(x_1,p_1)\| \text{ such that } F(x,p)=0\} \tag{MS}\]

It can be interpreted as a PALC in which the hyperplane is adapted at every step.

Predictor

The possible predictors alg::AbstractContinuationAlgorithm are listed in Predictors - Correctors. They can be used to create a Moore-Penrose algorithm like MoorePenrose(predictor = PALC()), the default being predictor = PALC() whose own tangent (e.g. Secant(), Bordered(), ...) is used to predict the next point.

Corrector

The corrector is the Gauss Newton algorithm applied to (MS).

Linear Algebra

Norm

The option normC of continuation specifies the norm used to evaluate the distance in (MS). The dot product (resp. norm) used in the (iterative) linear solvers is LinearAlgebra.dot (resp. LinearAlgebra.norm). It can be changed by importing these functions and redefining it. Note that by default, the $L^2$ norm is used.

Linear problem

The linear solver for the linear problem associated to (MS) is set by the option linear_algo in continuation or, more generally, in the constructor MoorePenrose(; predictor, method, ls). It is one of Bordered linear solvers (BLS).

Algorithm for solving (MS)

Let us write $y\equiv(x,p)\in\mathbb R^{N+1}$. In order to solve for the argmin, we apply the newton algorithm with jacobian belonging to $\mathbb R^{N\times (N+1)}$:

\[y^{k+1} = y^k -d_yF(y^k)^+F(y^k)\]

where the superscript $^+$ indicates the Moore-Penrose pseudoinverse of rank $N$.

Direct case

This is the default, triggered by the option MoorePenrose(method = BifurcationKit.direct): the Newton step is obtained by solving the underdetermined linear system $d_yF(y^k)\delta = F(y^k)$ in the least-norm sense, i.e. with the backslash operator \.

For the option MoorePenrose(method = BifurcationKit.pInv), the pseudoinverse is computed with pinv. These two methods are suited for small to moderate dimensional problems.

Iterative case

In this case, triggered by the option MoorePenrose(method = BifurcationKit.iterative), the pseudoinverse is computed with an iterative method described in ^[Meijer]:

\[\left\{\begin{array}{l} y_{1}^{j+1}=y_{1}^{j}-\left(\begin{array}{c} F_{y}\left(y_{1}^{j}\right) \\ \left(\phi_{1}^{j}\right)^{\top} \end{array}\right)^{-1}\left(\begin{array}{c} F\left(y_{1}^{j}\right) \\ 0 \end{array}\right) \\ \phi_{1}^{j+1}=\left(\begin{array}{c} F_{y}\left(y_{1}^{j+1}\right) \\ \left(\phi_{1}^{j}\right)^{\top} \end{array}\right)^{-1}\left(\begin{array}{l} 0 \\ 1 \end{array}\right), \quad j=0,1,2, \ldots \end{array}\right.\]

We initialise $\phi_1^0$ with the tangent. The linear systems above are solved with the bordered linear solver ls passed to the constructor (default MatrixBLS() when not specified).

Step size control

Each time the corrector fails, the step size $ds$ is halved. This has the disadvantage of having lost Newton iterations (which costs time) and imposing small steps (which can be slow as well). To prevent this, the step size is controlled internally with the idea of having a constant number of Newton iterations per point. This is in part controlled by the aggressiveness factor a in ContinuationPar.

References