# Is there a software package that can compute the 1-dimensional preimage of a point?

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.

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$$.

• 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 ? Jan 27 at 14:23
• Oops ! I should have written ""choose the one [the unit vector] consistent with the previous step". Jan 27 at 14:34
• @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$). Jan 27 at 14:45