Moore-Penrose continuation

WIP

This is work in progress. The interface will be improved inn the future.

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 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 tangent::AbstractTangentPredictor are listed in Predictors - Correctors. They can be used to create a Moore-Penrose algorithm like MoorePenrose(tangent = PALC())

Corrector

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

Linear Algebra

Let us discuss more about the norm and dot product. First, the option normC 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.

The linear solver for the linear problem associated to (MS) is set by the option linearAlgo in continuation: 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

In this case, triggered by the option MoorePenrose(directLS = true), the pseudoinverse is computed with pinv.

Iterative case

In this case, triggered by the option MoorePenrose(directLS = false), 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.

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

  • Meijer

    Meijer, Dercole, and Oldeman, “Numerical Bifurcation Analysis.”