I have a smooth function $F: \mathbb{R}^n \to \mathbb{R}^{n-1}$ and points $x_0, y_0$ with $F(x_0) = y_0$. For theoretical reasons, I know that $y_0$ is a regular value of $F$, which means that the Jacobian at every point that maps to $y_0$ has full rank $n{-}1$. Since $y_0$ is a regular value, the preimage $F^{-1}(y_0)$ is a 1-dimension manifold in $\mathbb{R}^n$. For theoretical reasons, I know that this preimage is a smooth loop.

I wonder whether there is a numerical routine that can take these as input:

  • a function that computes $F(x)$
  • a function that computes the Jacobian $DF(x)$
  • $x_0$ and $y_0$
  • an integer $M$

and then returns $M$ points, in order, on the smooth loop $F^{-1}(y_0)$. The returned set should begin with the given point $x_0$, and be uniformly distributed on the loop in some natural way.

Thank you so much

Added later: @lightxbulb. Following your approach, I have written a routine that is working well enough for me.

  1. This step computes the "extra singular vector" $u$ of $DF(x_k)$, and interprets it as a unit tangent vector to the loop at $x_k$. There are 2 such vectors, so to pick one for consistency, it does a sign test. It adds $u$, as a new bottom row, to the $(n{-}1) \times n$ Jacobian matrix $DF(x_k)$. and then computes the determinant of the new $n \times n$ matrix. If this is negative, $u$ is reversed. Set $z = x_k + \tau u$. I think it may be possible to devise a different sign test that uses an $(n{-}1) \times (n{-}1)$ determinant, but I haven't tried that yet.
  2. For the projection step, instead of using the pseudo-inverse, it uses an n'th constraint: the projected point must lie on the hyperplane through $z$ with normal $u$. It then uses a "canned" Newton-Raphson n-dimensional rootfinder to find $x_{k+1}$.

In summary, to go from $x_k$ to $x_{k+1}$, first there is a move in the direction $u$, and then a move orthogonal to $u$. The intermediate point $z$ is discarded.

I am not using the integer $M$ at this time.


1 Answer 1


I haven't read any literature on the subject, so there are likely much better approaches, but here's a sketch of what you can do. You have a set defined as $S = \{x \in\mathbb{R}^n \,:\, F(x)-y_0 = 0\}$. The kernels of $DF$ at points of $S$ give you the tangent spaces to $S$. Using this knowledge you can iterate the following process:

  1. Compute $DF(x_k)$, and do a step $z_{k+1} = x_k + \tau n$ where $n$ is a unit vector in the kernel of $DF(x_k)$ (e.g. if you compute the SVD that is the extra singular vector that corresponds to no singular value). $(x_k, x_k+s n)$ for $s\in\mathbb{R}$ gives you the tangent space at $x_k$.
  2. Project $z_{k+1}$ onto $S$ by iteratively solving $F(z) = y_0$ with an initial guess $z = z_{k+1}$. Set $x_{k+1}=z$.

You can try to use Newton-Rhapson for finding the projection in step 2: $z = z - (DF(z))^{+}F(z)$, with initial guess $z=z_{k+1}$. Here $(DF(z))^{+}$ is to be understood as the Moore-Penrose pseudoinverse, e.g. you could use CGNR to compute its application to $F(z)$.

Step one is essentially something like explicit Euler, where you say $\frac{dx}{dt} = n(x)$ and you discretize this as $z_{k+1} = x_k + \tau n(x_k)$. If you don't do step two and you set $x_{k+1} = z_{k+1}$ then assuming your steps are small enough you should get something close to the curve, but explicit Euler has stability issues. Step 2 is implementing something like implicit Euler with a guess given by explicit Euler. You can of course generalize that to an arbitrary time stepping scheme. Just pick you favourite discretisation of $\frac{dx}{dt} = n(x)$ where $n(x)$ is the singular vector corresponding to the kernel for $DF(x)$. After you are done you can pick out the $M$ points you desire, potentially by sampling some interpolant of the intermediate points you got from the time-stepping.

Generally $n(x_k)$ can be chosen in two ways when $DF(x_k)$ is of rank $n-1$, and you would generally choose the one consistent with the previous step. If $DF(x_k)$ is of rank less than $n-1$ then we can choose $n(x_k)$ in infinitely many ways. You could try to make it smoothly varying, e.g. by projecting the $n(x_{k-1})$ from the previous step onto the kernel of $DF(x_k)$, and find the new $n(x_k)$ in such a way. In the worst case scenario the projection may be zero, but then I would argue that you took a step that is too large. There may be some pathological scenarios where this doesn't work even when the step is small enough (e.g. if $S$ is not differentiable or is even discontinuous) but then $DF$ would not really be a representation of the total derivative either.

Also since you know that's a loop, at the initial point you can start two evolutions $\frac{dx}{dt} = n$ and $\frac{dx}{dt} = -n$ where initially you have $n(x_0)$ of course. Evolving both may allow you to connect them at the other side of the loop and may be more accurate, sine you know that $x_0 \in S$.

  • $\begingroup$ It certainly looks like this approach will work. But I was hoping for an existing routine in a public package. If one could always "choose the one [the unit vector] with the previous step", then it appears that we have defined a smooth unit vector field. And we are solving a special ODE - an ODE with n-1 "conservation laws". Does that make sense ? Maybe the place to look for packages is ODE solvers that support many conservation laws ? $\endgroup$ Jan 27 at 14:23
  • $\begingroup$ Oops ! I should have written ""choose the one [the unit vector] consistent with the previous step". $\endgroup$ Jan 27 at 14:34
  • $\begingroup$ @GlennDavis I think your interpretation makes sense but I don't know of any packages. The way I described the ODE problem it is certainly an evolution along a smooth vector field. Of course maybe there is another way to constrain the choice of $n$ but I felt that the smoothness assumption is the most natural, hence why I suggested projecting $n(x_k)$ onto the tangent space at $x_{k+1}$ in order to find $n(x_{k+1})$. You can achieve the latter by using again CGNR by the way (if you solve $Aw = n(x_k)$ with CGNR with an initial guess $x=0$ it gives you $w = A^+n(x_k)$, then $n_{k+1} = AA^+n_k$). $\endgroup$
    – lightxbulb
    Jan 27 at 14:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.