Implicitly Defined Manifolds
An implicitly defined manifold $M$ is the set of points (the Solution Set) satisfying a set of equations $${\bf F}({\bf u})=0,\qquad\qquad{\bf F}:R^n\rightarrow R^{n-k}$$ with ${\bf F}$ a smooth mapping.
Continuation Methods were originally created for one dimensional manifolds and are sometimes called path following algorithms. They were formulated as a way of finding solutions of $$G({\bf u},\lambda)=0,\qquad F:R^n\times R\rightarrow R^n.$$ This can obviously be rewritten in the form above.
Implicitly defined manifolds include a large range of surfaces. The simple sphere has an implicit representation $${\bf u}^T{\bf u} = R^2$$. Wikipedia has a list of other examples, and there are other lists here, here, and others. They are popular in computer graphics as isosurfaces because of their flexibility and occur naturally in scientific visualization and medical imaging.
However, they are many more examples in higher dimensions. Many laws in physics are conservation laws and those define implicit manifolds, as points in infinite dimensional spaces with parameters. In dynamics special trajectories can be posed as two point boundary value problems. Fixed points, periodic orbits, homo and heteroclinic orbits, but also problems in control theory, deformation of elastic beams, robotics and orbital mechanics, all can be written as implicit manifolds. You might say that most problems in applied math are solved by trying to find a regular system of equations whose solution is the object you're trying to find. We'd prefer an explicit answer $x=something$, but instead we have to be satisfied with a list of properties the object has $f(x)=something$ and be satisfied that the only object with those properties is the one we're looking for (i.e. the system is regular). And while I usually fall back on Newton's method, there is a large catalogue of iterative methods for solving these systems.
If the Jacobian ${\bf F}_{\bf u}$ at a point is full rank ($n-k$) the Implicit Function Theorem says that the solution set near that point is a smooth $k$ dimensional manifold. If we parameterize using a mapping from ${\bf s}\in R^k$ to a point ${\bf u}$ on $M$, we can develop a Taylor expansion for the parameterization ${\bf u}({\bf s})$ with ${\bf F}({\bf u}(0))=0$ and a parameterization $$\begin{equation} \begin{cases} {\bf F}({\bf u}({\bf s}))=0 \\ \\ {\left( {\bf u}_{\bf s}(0)\right)}^T \left( {\bf u}({\bf s}) - {\bf u}(0) \right) = {\bf s}. \end{cases} \end{equation}$$ We can keep differentiating the composition of ${\bf F}$ and ${\bf u}({\bf s})$ to get linear equations at ${\bf s}=0$ for the derivatives of ${\bf u}$ $$\begin{align} &\begin{cases} {\bf F}_{\bf u}{\bf u}_{\bf s}=0\\ \\ {\bf u}_{\bf s}^T {\bf u}_{\bf s} = I \end{cases}\\ \\ &\begin{cases} {\bf F}_{\bf u}{\bf u}_{\bf ss}+{\bf F}_{\bf uu}{\bf u}_{\bf s}{\bf u}_{\bf s}=0\\ \\ {\bf u}_{\bf s}^T {\bf u}_{\bf ss} = 0 \end{cases}\\ \\ &\begin{cases} {\bf F}_{\bf u} {\bf u}_{\bf sss} + 3 {\bf F}_{\bf uu} {\bf u}_{\bf s} {\bf u}_{\bf ss} + {\bf F}_{\bf uuu} {\bf u}_{\bf s} {\bf u}_{\bf s} {\bf u}_{\bf s}=0\\ \\ {\bf u}_{\bf s}^T {\bf u}_{\bf sss} = 0 \end{cases} \end{align}$$ Here ${\bf F}$ and its derivatives are evaluated at ${\bf u}(0)$ and on ${\bf u}$ and its derivatives ated at ${\bf s}=0$. If ${\bf F}_{\bf u}$ is full rank and ${\bf u}_{\bf s}$ is an orthonormal basis for its null space these are regular systems for the derivatives of ${\bf u}$.
Rearranging to collect the highest derivatives of ${\bf u}$ on the left \begin{align} &\begin{cases} {\bf F}_{\bf u}{\bf u}_{\bf s}=0\\ \\ {\bf u}_{\bf s}^T {\bf u}_{\bf s} = I \end{cases}\\ \\ &\begin{cases} {\bf F}_{\bf u}{\bf u}_{\bf ss} = - {\bf F}_{\bf uu}{\bf u}_{\bf s}{\bf u}_{\bf s}\\ \\ {\bf u}_{\bf s}^T {\bf u}_{\bf ss} = 0 \end{cases}\\ \\ &\begin{cases} {\bf F}_{\bf u} {\bf u}_{\bf sss} = - 3 {\bf F}_{\bf uu} {\bf u}_{\bf s} {\bf u}_{\bf ss} + {\bf F}_{\bf uuu} {\bf u}_{\bf s} {\bf u}_{\bf s} {\bf u}_{\bf s}\\ \\ {\bf u}_{\bf s}^T {\bf u}_{\bf sss} = 0 \end{cases} \end{align} It gets a little messy, but systems for the higher order derivatives can also be written. The coefficients of the derivatives of ${\bf F}$ come from the derivative of a composition (``Faá di Bruno's Formula''). A little care is needed, as these are tensors. For example \begin{align*} 3 {\bf F}_{\bf uu} {\bf u}_{\bf s} {\bf u}_{\bf ss} & = \\ & \sum_p^{p\lt n}\sum_q^{q\lt n} \left( \frac{\partial F^i}{\partial u^p\partial u^q}\frac{\partial u^p}{\partial s^j}\frac{\partial^2 u^q}{\partial s^k \partial s^l} +\frac{\partial F^i}{\partial u^p\partial u^q}\frac{\partial u^p}{\partial s^k}\frac{\partial^2 u^q}{\partial s^j \partial s^l} +\frac{\partial F^i}{\partial u^p\partial u^q}\frac{\partial u^p}{\partial s^l}\frac{\partial^2 u^q}{\partial s^j \partial s^k} \right) \end{align*}
To third order $$ {\bf u}({\bf s})\sim {\bf u}(0) + {\bf u}_{\bf s}(0) {\bf s} + \frac{1}{2} {\bf u}_{\bf ss}(0) {\bf s}{\bf s} + \frac{1}{3!} {\bf u}_{\bf sss}(0) {\bf s}{\bf s}{\bf s}.$$ As in other polynomial approximations, higher order approximations tend to sacrifice convergence radius for order. So although a third order approximation might be more accurate near the origin than a first order approximation, it might only be better nearer ${\bf s}=0$. The question of how to locally approximate ${\bf u}({\bf s})$ is regular perturbation theory, or asymptotic analysis, and it can get complicated. The group down the hall when I worked at IBM in the 1980's worked on how to evaluate things like square roots and trig functions, and they had many tricks.
This local approximation is enough to define a chart centered at a point ${\bf u}_0\in M$. An orthonormal basis ${\bf u}_{\bf s}(0)$ for the tangent space is an orthonormal basis for the null space of the Jacobian ${\bf F}_{\bf u}$. A point ${\bf s}\in R^k$, $|{\bf s}|\leq R$ to be projected to a point ${\bf u}({\bf s})\in M$ with ${\bf u}({\bf s})={\bf u}$ by solving the nonlinear system $$\begin{align} &{\bf F} ({\bf u}({\bf s}))=0 \\ &{({\bf u}_{\bf s}(0)}^T ( {\bf u}({\bf s}) - {\bf u}(0) ) = {\bf s} \end{align}$$ My favorite is to use Newton's method with an initial guess of ${\bf u}({\bf s})\sim {\bf u}(0) + \Phi {\bf s}$ where $\Phi$ is the matrix whose columns are the tangent vectors (an $n\times k$ matrix). Newton converges if initial guess is close to the solution (the Kantorovich condition). Keeping the radius $R$ small enough that $$ \frac{1}{2} |{\bf u}_{\bf ss}| R^2 \lt \epsilon $$ is a good idea. This can be interpreted geometrically as the distance from the tangent space to $M$. Another approach is to say the radius if more than three iterations are required and too small if fewer than two.
This is when things get interesting. If the left and right null spaces of the Jacobian are finite dimensional, with orthonormal bases $\Phi$ and $\Psi$ respectively, the Lyapunov-Schmidt reduction can be used to separate the problem into a linear part and a smaller nonlinear part. The chart mapping is no longer guaranteed to be a mapping from $R^k\rightarrow M$, but may be one to many.
The Lyapunov-Schmidt decomposition splits the system ${\bf F}({\bf u})=0$ into the equivalent system $$\begin{align} & (I-\Psi\Psi^T){\bf F}({\bf u}(0) + \Phi {\bf s} + {\bf a}({\bf s}))=0\\ \\ & \Psi\Psi^T{\bf F}({\bf u}(0) + \Phi {\bf s} + {\bf a}({\bf s}))=0\\ \\ & {\bf a}(0) = 0 \\ \\ & \Phi^T {\bf a}({\bf s}) =0 \end{align}$$ A big improvement, no? Well, yes as a matter of fact. The system $$\begin{align} & (I-\Psi\Psi^T){\bf F}({\bf u}(0) + \Phi {\bf s} + {\bf a})=0\\ \\ & \Phi^T {\bf a} = 0 \end{align}$$ for ${\bf a}$ has a regular solution ${\bf a}=0$ if ${\bf F}({\bf u}(0))=0$. The inverse of the Jacobian is the pseudo-inverse of the Jacobian of ${\bf F}$. Another way to view it is via the fundemental theorem of linear algebra. The restriction of the linear operator to a map from the complement of its null space (spanned by $(I-\Phi\Phi^T)$) to its range (spanned by $(I-\Psi\Psi^T)$) is full rank.
So the Implicit Function Theorem says that in some neighborhood of ${\bf s}=0$ there is a unique mapping ${\bf a}({\bf s})$ with ${\bf a}(0) = 0$ which satisfies $(I-\Psi\Psi^T){\bf F}({\bf u}(0) + \Phi {\bf s} + {\bf a}({\bf s}))=0$. We can churn out the coefficients in a Taylor series. At ${\bf s}=0$ They satisfy $$\begin{align} & {\bf a} = 0 \\ \\ &\begin{cases} (I-\Psi\Psi^T) {\bf F}_{\bf u} \left( \Phi + {\bf a}_{\bf s} \right) = 0 \\ \\ \Phi^T {\bf a}_{\bf s} = 0 \end{cases}\\ \\ & \begin{cases} (I-\Psi\Psi^T) {\bf F}_{\bf u} {\bf a}_{\bf ss} + (I-\Psi\Psi^T){\bf F}_{\bf uu}\left( \Phi + {\bf a}_{\bf s} \right)\left( \Phi + {\bf a}_{\bf s} \right)= 0 \\ \\ \Phi^T {\bf a}_{\bf ss} = 0 \end{cases}\\ \\ &\begin{cases} (I-\Psi\Psi^T) {\bf F}_{\bf u} {\bf a}_{\bf sss} + 3 (I-\Psi\Psi^T){\bf F}_{\bf uu}\left( \Phi + {\bf a}_{\bf s} \right) {\bf a}_{\bf ss} + (I-\Psi\Psi^T){\bf F}_{\bf uuu}\left( \Phi + {\bf a}_{\bf s} \right)\left( \Phi + {\bf a}_{\bf s} \right)\left( \Phi + {\bf a}_{\bf s} \right)= 0 \\ \\ \Phi^T {\bf a}_{\bf sss} = 0 \end{cases} \end{align}$$ Collecting the highest derivatives of ${\bf a}$ on the left and noting that ${\bf F}_{\bf u} \Phi = 0$: $$\begin{align} & {\bf a} = 0 \\ \\ & \begin{cases} (I-\Psi\Psi^T) {\bf F}_{\bf u} {\bf a}_{\bf s} = 0 \\ \\ \Phi^T {\bf a}_{\bf s} = 0 \end{cases} \\ \\ &\begin{cases} (I-\Psi\Psi^T) {\bf F}_{\bf u} {\bf a}_{\bf ss} = - (I-\Psi\Psi^T){\bf F}_{\bf uu} \Phi\Phi \\ \\ \Phi^T {\bf a}_{\bf ss} = 0 \end{cases}\\ \\ &\begin{cases} (I-\Psi\Psi^T) {\bf F}_{\bf u} {\bf a}_{\bf sss} = - 3 (I-\Psi\Psi^T){\bf F}_{\bf uu} \Phi {\bf a}_{\bf ss} - (I-\Psi\Psi^T){\bf F}_{\bf uuu}\Phi\Phi\Phi \\ \\ \Phi^T {\bf a}_{\bf sss} = 0 \end{cases} \end{align}$$ The first pair means that ${\bf a}_{\bf s}=0$. I've used that in the third and fourth pairs. The linear systems for the derivatives of ${\bf a}$ are identical, with different right hand sides.
That leaves us with the Bifurcation Equations $$\Psi^T{\bf F}({\bf u} + \Phi {\bf s} + {\bf a}({\bf s}))=0$$ with the known function ${\bf a}({\bf s})$. The range of $\Psi^T{\bf F}$ is $dim(\Psi)$ dimensional, and the domain is $dim(\Phi)$ dimensional. The Jacobian w.r.t. ${\bf s}$ is $\Psi^T{\bf F}_{\bf u}\left(\Phi+{\bf a}_{\bf s}\right)$, which is identically zero. The Bifucration Equations are usually solved by reducing them to an equivalent algebraic system called the algebraic bifurcation equations. This can be done by introducing a parameter $\epsilon$ and scaling ${\bf s}=\epsilon\tilde{\bf s}$ and ${\bf F}$ by $\epsilon^{-m} $, then considering the Jacobian w.r.t $\epsilon$ of $$\begin{equation} \epsilon^{-m} \Psi^T{\bf F}({\bf u}(0) + \epsilon \Phi \tilde {\bf s} + {\bf a}(\epsilon\tilde{\bf s}))=0 \end{equation}$$ If $m$ is the leading order of the system (i.e. quadratic $m=2$), this brings the higher order derivatives down to order zero and you choose $\tilde {\bf s}$ so that term is zero (the algebraic system) and get conditions on the new linear term being full rank. For example, with $m=2$ \begin{align} \frac{1}{\epsilon^2} \Psi^T{\bf F}({\bf u}(0) + \epsilon \Phi \tilde {\bf s} + {\bf a}(\epsilon\tilde{\bf s})) \sim & \frac{1}{2} \Psi^T{\bf F}_{\bf uu} \left( \Phi \tilde s \right) \left( \Phi \tilde s \right)\\ \\ + & \frac{1}{6}\epsilon\left( \Psi^T{\bf F}_{\bf uuu} \left( \Phi \tilde s \right) \left( \Phi \tilde s \right) \left( \Phi \tilde s \right) +2\Psi^T{\bf F}_{\bf uu} \left( \Phi \tilde s \right) \left( {\bf a}_{\bf ss} \tilde {\bf s} \tilde {\bf s} \right) \right) + O(\epsilon^2). \end{align} We choose $\tilde{\bf s}$ to make the leading term zero -- the Algebraic Bifurcation Equations, in this case a system of quadratics. If the term that is linear in $\epsilon$ is non-zero the IFT can be used to give solution curves $\tilde{\bf s}(\epsilon)$, one for each root of the Algebraic Bifurcation Equations. (This implies that the Bifurcation Equations are contact equivalent to the Algebraic Bifurcation Equations.)
I should be able to write equations for a Taylor expansion of the functions $\tilde{\bf s}(\epsilon)$. Can I just say that there's no room in the margins?
Michael E. Henderson michael.e.henderson.10590 at gmail.com