The "surface" of the SDF is given as the set $S = f^{-1}[0] =\{p\in\mathbb{R}^n\,:\, f(p) =0 \}$. Thus you can generate a number of points $p_1,\ldots,p_m\in\mathbb{R}^n$ and set them as initial guesses to a Newton iteration:
\begin{align}
p^{k+1}_i &= p^k_i - \theta^k_i (J_f(p^k_i))^+f(p^k_i) \\
&= p_i^k -\theta^k_i(J_f(p_i^k)^TJ_f(p_i^k))^{-1}J_f(p^k_i)^Tf(p^k_i)\\
&= p^k_i + \underbrace{\frac{f(p_i^k)\theta^k_i}{\|\nabla f(p^k_i)\|_2}}_{\alpha^k_i}\underbrace{\left(-\frac{\nabla f(p^k_i)}{\|\nabla f(p^k_i)\|_2}\right)}_{d^k_i}\\
&= p^k_i +\alpha^k_i d^k_i.
\end{align}
For an actual sdf $\|\nabla f(p^k_i)\|_2 = 1$, and if you set the damping factor $\theta^k_i=1$, then you get $$p^{k+1}_i = p^k_i - f(p^k_i)\nabla f(p^k_i).$$
Note that $f(p^k_i)$ is the signed distance to the surface. So you are doing sphere tracing but with the direction from the negated gradient.
Another option, if you have a point $p$ alread on the surface, is to construct an orthonormal basis for the space orthogonal to $\nabla f(p)$. The latter would be the tangent space. And then take steps along each axis of this space and subsequently project onto the surface. The size of the steps could use information from the Hessian in order to figure out the curvature in each direction. This is more in line with the idea of tracing out the surface starting from a point, while the former is based on arbitrary points traveling down to the surface along the gradient field induced by the sdf. On that note, sometimes it may be beneficial to take $\theta_k$ larger and sometimes smaller in order to boost convergence speed.
One could equivalently look at this problem as $$\min_{p}\frac{1}{2}\|f(p)\|^2_2.$$
The gradient of the above is precisely $f(p)\nabla f(p)$, so you could just as well interpret this as a gradient descent for the above energy.
If you don't have an sdf but a function that underestimates the distance then I would use $$p^{k+1}_i = p^k_i - \theta^k_i g(p^k_i) \frac{\nabla g(p^k_i)}{\|\nabla g(p^k_i)\|_2}.$$ The convergence would be slower depending on how much $g$ underestimates the real distance. To accelerate both methods you can use something like a heavy-ball method, or also look at (nonlinear) conjugate gradients to handle more unpleasant sdfs. If your sdf is non-differentiable then maybe subgradient descent may be useful.
Edit: I provided some more details in the non-differentiable setting below.
---
**Lipschitz Continuous Functions and Signed Distance Bound Functions**
Suppose $S$ is a closed surface, i.e., $S = \partial\Omega$ for some set $\Omega$ (this is the "outside" of the SDF). Define $d(x,S) = \inf_{y\in S}d(x,y)$ for some metric $d:\mathbb{R}^n\times\mathbb{R}^n\to[0,\infty)$. The function $f:\mathbb{R}^n\to\mathbb{R}$ is a signed distance function:
$$f(x) =
\begin{cases} d(x,S), &x\in\Omega,\\
-d(x,S), &x\not\in\Omega.\end{cases}$$
A signed distance bound function (SDBF) $g:\mathbb{R}^n\to\mathbb{R}$ for $\Omega$ is a function such that $g^{-1}[0] = S$ but unlike $f$ it underestimates the true distance:
$$\begin{cases} 0\leq g(x) \leq f(x) = d(x,S), &x\in\Omega, \\
0\leq -g(x) \leq -f(x) = d(x,S), &x\not\in\Omega.\end{cases}$$
The metric I consider from now on is given through the Euclidean norm: $d(x,y) = \|x-y\|_2$.
A natural candidate to construct signed distance bound functions from are Lipschitz continuous functions.
For a Lipschitz continuous function $h:\mathbb{R}^n\to\mathbb{R}$ there exists a constant $L(h)\geq 0$ such that:
$$|h(x)-h(y)|\leq L(h)\|x-y\|.$$
Note that $h\equiv 0$ is Lipschitz continuous, with $L(h) = 0$, but it is not very useful. On the other hand it does not satisfy $S=h^{-1}[0]$.
Suppose $h$ satisfies $h^{-1}[0]=S$ and that it is Lipschitz continuous,. and it is positive on $\Omega$ and negative on $(\mathbb{R}^n\setminus S)\setminus \Omega$. Then $g(x) = h(x) / \lambda$ where $\lambda \geq L(h)$ is a signed distance bound function for $\Omega$.
The proof is simple and can also be found in [Hart's paper][1].
Let $y\in S$ be such that $y\in \arg\min_{z\in S}\|x-z\|$, i.e., $\|x-y\| = d(x,S)$. We can show:
$$|g(x)|\stackrel{y\in S}{=}|g(x)-g(y)| \leq \frac{L(h)}{\lambda} \|x-y\|\stackrel{L(f)/\lambda \leq 1}{\leq} \|x-y\| = d(x,S).$$
---
**Standard Functions Used in Constructing an SDBF**
Typical functions used in the constructions of signed distance bound functions are usually totally differentiable, e.g., $\cos$, $\sin$, $\exp$, $\ln$, linear combinations (this includes affine transforms), polynomials, powers. Inigo Quilez even has some examples involving computation of gradients: [see this for 2D SDFs\ gradients][2], and [this for noise gradients][3].
The issue would arise with functions such as absolute value, min, max, `clamp(x,a,b)=min(max(x,a),b)`, `mod(x,y)`, and worst of all functions involving a `floor` (e.g. `mod` and `fract`), `ceil`, or `round`, or potentially an `if` statement, where you can't really argue about what happens on the boundary cases without some manual work. Note that the presence of these functions doesn't necessarily mean you cannot compute a reasonable gradient (see the examples from Inigo Quilez above), but those functions are not differentiable, and in the case of `floor` they are not even continuous (so definitely not locally Lipschitz continuous, however their composition with other functions may be continuous and even differentiable).
From those `abs` and `max` are the nicest since both are subdifferentially regular, so you could use subgradient calculus with them without much issue. The `min` function is not subdifferentially regular, so it is not very nice to work with, although its Clarke subgradient exists. On the other hand the subgradient may be unnecessary here as you can use [a smooth minimum][4]. The worst offender is the `floor` function which is discontinuous at integers and I don't think we can do anything with subdifferential calculus to it. This doesn't mean that the resulting function from a composition with a floor function is discontinuous. You can even have an infinitely smooth periodic function as a result of `floor`, but the moment it appears it disallows you from exploiting that analytically (if you just have a black box view of the whole thing) - e.g. `sin(2*PI*fract(x))` is virtually the same as `sin(2*PI*x)` but an automatic differentiation approach wouldn't necessarily know that. Note that a numerically computed gradient has no issue with that whatsoever - it cares not that `floor` is discontinuous.
Next I discuss some theory on nonsmooth optimization, although I feel like this is likely overkill.
---
**Lipschitz Continuous Functions and Differentiation**
Lipschitz continuous functions do not have to be differentiable. On the other hand they are differentiable almost everywhere ([see Rademacher's theorem][5]). As discussed, one option would be then to apply random small shifts when you are not at a differentiable point and take the gradient that you sampled. There are much more involved elaboration similar to this idea, e.g. see ["A ROBUST GRADIENT SAMPLING ALGORITHM FOR NONSMOOTH, NONCONVEX OPTIMIZATION"][6] by Burke et al.
Another option is to try to generalize the gradient at the non-differentiable points. What follows below is essentially a discussion similar to the one in ["Understanding Notions of Stationarity in Non-Smooth Optimization"][7] by Li et al.
The Bouligand subgradient (this is often termed subdifferential, but that terminology is misleading) is defined as follows:
$$\partial_B f(x) = \{s\in\mathbb{R}^n\,|\, \exists x^k\to x, \,:\, \exists \nabla f(x_k), \, \nabla f(x_k)\to s\}.$$
It takes the limits of nearby gradients (at differentiable points) as they approach the non-differentiable point and makes a set of those.
However, when you compute the gradient of $f(x) = |x|$, i.e. $\partial_B f(0) = \{-1,1\}$, there is the undesirable property that $0\not\in \partial_B f(0)$, despite of the fact that $0$ is the global minimum of the above function. Remember that in basic calculus you have $\nabla f(x) = 0$ at extrema. To rectify this, one may consider the Clarke subgradient (see ["GENERALIZED GRADIENTS AND APPLICATIONS"][8] for the original work, but reading Li's exposition feels much nicer):
$$\partial^{\circ} f(x) = \text{conv}(\partial_B f(x)),$$
which takes the convex hull of the Bouligand subgradient. Then it follows that $\partial^{\circ} f(0) = [-1,1]$ for $f(x) = |x|$ and $0\in\partial^{\circ} f(0)$. The Clarke subgradient also happens to be a generalization of the convex subgradient.
---
**The Clarke (Generalized) Directional Derivative and Subgradient**
You can find the definition of the Clarke directional derivative and gradient [on wikipedia][9], or in his book ["Optimization and Nonsmooth Analysis"][10] (page 10):
\begin{align}f^{\circ}(x; v) &:= \lim_{y\to x}\sup_{t\downarrow 0}\frac{f(y+t v)-f(y)}{t} = \inf_{\substack{\epsilon>0,\\ \delta>0}}\sup_{\substack{\|y-x\|\leq\epsilon,\\ 0<t<\delta}}\frac{f(y+t v)-f(y)}{t},\tag{Clarke derivative} \\
\partial^{\circ} f(x) &:= \{u\in\mathbb{R}^n\,:\, u^Tv\leq f^{\circ}(x; v), \,\, \forall v\in\mathbb{R}^n\} \tag{Clarke subgradient}.\end{align}
You have to actually be careful with the definitions, since for example in his other book (["Nonsmooth Analysis and Control Theory"][11] page 70 and 72) Clarke uses the exact same symbol and terminology (i.e. generalized gradient) to refer to the subdifferential instead of the subgradient:
$$
\partial^{\circ} f(x) := \{\zeta\in (\mathbb{R}^n)^*\,:\, \zeta(v)\leq f^{\circ}(x; v), \,\, \forall v\in\mathbb{R}^n\}. \tag{Clarke subdifferential}$$
The relation between the two is $\zeta(v) = u^Tv$ (i.e. like between the differential and the gradient: $df(x)(v) = \nabla f(x)^Tv$). The subdifferential is more general and can be defined even without an inner product, but in $\mathbb{R}^n$ it doesn't really matter. In either case, I'll stick to the convention that $\partial^{\circ} f$ refers to the subgradient.
The Clarke (directional) derivative is a generalization of the one-sided directional derivative $\partial_v f(x)$ and the Clarke subgradient is a generalization of the gradient. Remember $\nabla f(x)^T v = \partial_v f(x)$ for totally differentiable functions, while here $u^Tv \leq f^{\circ}(x;v)$ for $u\in\partial^{\circ} f(x)$. The Clarke derivative is positive homogeneous and sublinear. That is, if $v,w\in\mathbb{R}^n$ and $\alpha\geq 0$ then
\begin{align}
f^{\circ}(x,\alpha v) = \alpha f^{\circ}(x,v), \tag{positive homogeneity}\\
f^{\circ}(x,v+w) \leq f^{\circ}(x,v) + f^{\circ}(x,w). \tag{sublinearity}
\end{align}
---
**Clarke Calculus**
I list some properties of the Clarke gradient that hold for locally Lipschitz functions (you can find them in Chapter 2.2 of [Clarke's book][11], or chapter 2.3 of [his other book][10], or some in [Li et al.'s exposition][7]):
- If $x^*$ is a local minimum or maximum for $f$ then $0\in\partial^{\circ} f(x^*)$.
- $f\in C^1 \implies \partial^{\circ} f(x) = \{\nabla f(x)\}$
- $\partial^{\circ}(\lambda f)(x) = \lambda \partial^{\circ} f(x)$
- $\partial^{\circ}(f+g)(x) \subseteq \partial^{\circ}f(x) + \partial^{\circ}g(x)$
- $\partial^{\circ} (fg)(x) \subseteq f(x)\partial^{\circ} g(x) + g(x)\partial^{\circ} f(x)$
- $\partial^{\circ} (f/g)(x) \subseteq \frac{g(x)\partial^{\circ} f(x) - f(x)\partial^{\circ} g(x)}{g^2(x)}$
- $F\in C^1$ and $g$ Lipschitz, then $\partial^{\circ} (g\circ F)(x) \subseteq J_F(x)^T\partial^{\circ} g(F(x))$.
- $\partial^{\circ} (g\circ (h_1,\ldots,h_m))(x) \subseteq \text{conv}\{\sum_{j=1}^m \alpha_j u_j\,:\, u_j\in \partial^{\circ}h_j(x), \,\, \alpha \in \partial^{\circ} g(h(x))\}.$
- $\partial^{\circ}\max_jf_j(x) \subseteq \text{conv}\{\partial^{\circ}f_i(x)\,:\,i\in I(x)\}$ where $I(x)$ is the set of indices of functions such that $f_i(x) = \max_j f_j(x).$
Already here you can see an issue with the calculus, namely that you have set inclusion but not necessarily equality. This is not nice in practice, since you may pick a vector from the larger set and it may turn out to not be a vector from the smaller set. However, if your functions are saubdifferentially regular, then all of the above turn into equalities (you additionally need that $f(x),g(x)\geq 0$ for the product rule; and $f(x)\geq 0$, $g(x)>0$, and that $f$ and $-g$ are regular at $x$ for the quotient rule).
A function $f$ is regular at $x$ if for all $v$ the directional derivatives $\partial_v f(x)$ exist and agree with $f^{\circ}(x; v)$. For example convex Lipschiz continuous functions are regular. A linear combination (with non-negative weights) of regular functions is regular. The maximum of regular functions is regular.
Finally here are some subgradients you may find useful:
\begin{align}
\partial^{\circ}|x| &= \begin{cases} \{1\}, &x>0, \\ [-1,1], &x=0, \\ \{-1\}, &x<0.\end{cases} \\
\partial^{\circ}\max\{x,c\} &= \begin{cases} \{0\}, &x<c, \\ [0,1], &f(x)=c, \\ \{1\}, &c<x.\end{cases} \\
\partial^{\circ}\min\{x,c\} &= \begin{cases} \{1\}, &x<c, \\ [0,1], &x=c, \\ \{0\}, &c<x.\end{cases}
\end{align}
For the last two note that $\partial^{\circ} (\min{x,c} + \max{x,c})(c) = \{1\}$ but $\partial^{\circ} \max\{x,c\} + \partial^{\circ} \min\{x,c\} = [0,1]$ which goes to illustrate the idea that
$$\partial^{\circ} (f+g)(x) \subseteq \partial^{\circ} f(x) + \partial^{\circ} g(x).$$
Notably, here $\min$ is not a regular function which spoils the equality. In either case I highly recommend reading ["Understanding Notions of Stationarity in Non-Smooth Optimization"][7] by Li et al. if you want more examples and applicable ideas. Though they don't have all of Clarke's calculus rules in their paper, so you will have to look at his books if you need that. But I believe I extracted the main properties you would use.
Maybe as a final note - you can also find pathological settings for subgradient descent, but I don't think those apply to your case. I am referring to ["Pathological subgradient dynamics"][12] by Daniilidis et al.
[1]: https://link.springer.com/article/10.1007/s003710050084
[2]: https://iquilezles.org/articles/distgradfunctions2d/
[3]: https://iquilezles.org/articles/morenoise/
[4]: https://iquilezles.org/articles/smin/
[5]: https://en.wikipedia.org/wiki/Rademacher%27s_theorem
[6]: https://cs.nyu.edu/~overton/papers/pdffiles/gradsamp.pdf
[7]: https://arxiv.org/abs/2006.14901
[8]: https://www.ams.org/journals/tran/1975-205-00/S0002-9947-1975-0367131-6/S0002-9947-1975-0367131-6.pdf
[9]: https://en.wikipedia.org/wiki/Clarke_generalized_derivative#Definitions
[10]: https://epubs.siam.org/doi/book/10.1137/1.9781611971309
[11]: https://link.springer.com/book/10.1007/b97650
[12]: https://arxiv.org/abs/1910.13604