<hr />
<p><strong>The Answer That was Here Originally</strong></p>
<p>The &quot;surface&quot; of the SDF is given as the set <span class="math-container">$S = f^{-1}[0] =\{p\in\mathbb{R}^n\,:\, f(p) =0 \}$</span>. Thus you can generate a number of points <span class="math-container">$p_1,\ldots,p_m\in\mathbb{R}^n$</span> and set them as initial guesses to a Newton iteration:
<span class="math-container">\begin{align}
p^{k+1}_i &amp;= p^k_i - \theta^k_i (J_f(p^k_i))^+f(p^k_i) \\
&amp;= 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)\\
&amp;= 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}\\
&amp;= p^k_i +\alpha^k_i d^k_i. 
\end{align}</span></p>
<p>For an actual sdf <span class="math-container">$\|\nabla f(p^k_i)\|_2 = 1$</span>, and if you set the damping factor <span class="math-container">$\theta^k_i=1$</span>, then you get <span class="math-container">$$p^{k+1}_i = p^k_i - f(p^k_i)\nabla f(p^k_i).$$</span></p>
<p>Note that <span class="math-container">$f(p^k_i)$</span> is the signed distance to the surface. So you are doing sphere tracing but with the direction from the negated gradient.</p>
<p>Another option, if you have a point <span class="math-container">$p$</span> alread on the surface, is to construct an orthonormal basis for the space orthogonal to <span class="math-container">$\nabla f(p)$</span>. 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 <span class="math-container">$\theta_k$</span> larger and sometimes smaller in order to boost convergence speed.</p>
<p>One could equivalently look at this problem as <span class="math-container">$$\min_{p}\frac{1}{2}\|f(p)\|^2_2.$$</span>
The gradient of the above is precisely <span class="math-container">$f(p)\nabla f(p)$</span>, so you could just as well interpret this as a gradient descent for the above energy.</p>
<p>If you don't have an sdf but a function that underestimates the distance then I would use <span class="math-container">$$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}.$$</span> The convergence would be slower depending on how much <span class="math-container">$g$</span> 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.</p>
<hr />
<p><strong>Summary for Implementation Purposes</strong></p>
<ul>
<li>For a density function <span class="math-container">$g$</span>, generate points from a normal distribution <span class="math-container">$\mathcal{N}(\mu,\sigma^2)$</span>. Do a subsequent pass to filter out the boundary points: <span class="math-container">$x,y$</span> are (fuzzy) boundary points if <span class="math-container">$\|x-y\|&lt;2R$</span> and <span class="math-container">$g(x)&lt;0$</span>, <span class="math-container">$g(y)&gt; 0$</span> or <span class="math-container">$g(x)&gt;0$</span>, <span class="math-container">$g(y)&lt;0$</span>. You can refine the boundary by studying <span class="math-container">$z = (1-t)x + ty$</span> and performing a bisection search for <span class="math-container">$t\in[0,1]$</span> with the goal to tighten the bound. If <span class="math-container">$f$</span> is an SDF (or not too far from an exact SDF) then <span class="math-container">$|g(x)|&lt;R$</span> can be used as a criterion that does not refer to other points, then you can reject all other points and you don't need two passes.
<ul>
<li>Pros: Makes no additional assumptions.</li>
<li>Cons: Can happen to be extremely slow.</li>
<li>Modifications: Larger <span class="math-container">$R$</span> allows working with fewer points but makes the boundary fuzzier. Low-discrepancy sequences may help if the sets <span class="math-container">$\{x\,|\,g(x)&lt;0\}$</span> and <span class="math-container">$\{x\,|\, g(x)&gt;0\}$</span> are regular enough. A good choice of <span class="math-container">$\mu,\sigma^2$</span> can improve the convergence speed. If you can estimate a bounding sphere/box that contains one of the above sets you can use a uniform distribution and generate points only there. In the latter case you could also use a regular grid or an adaptive structure such as an octree.</li>
</ul>
</li>
<li>If your function is totally differentiable you can use Newton-Raphson/gradient descent as <span class="math-container">$x^{k+1} = x^k - \frac{\theta^k g(x^k)}{\|\nabla g(x^k)\|^2_2}\nabla g(x^k)$</span>. For an SDF <span class="math-container">$\|\nabla g(x^k)\|_2=1$</span> so the denominator disappears. If <span class="math-container">$g$</span> is an SDBF then using <span class="math-container">$x^{k+1} = x^k - \frac{\theta^k g(x^k)}{\|\nabla g(x^k)\|^2}\nabla g(x^k)$</span> may be reasonable. If the iterates oscillate too much you can damp by picking <span class="math-container">$\theta^k&lt;1$</span>.
<ul>
<li>Pros: It should be quite performant. You could initialize with the random approach from above, and then update the points by using Newton-Raphson.</li>
<li>Cons: Assumes total differentiability and requires evaluating the gradient. Newton/gradient descent does not necessarily converge if the step size is too large - you can damp through <span class="math-container">$\theta^k$</span>. Even if you damp the step size the method may still get stuck in a local minimum for nasty densities or SDBFs.</li>
<li>Modifications: You can look into conjugate gradients (for nasty densities that may have features like the Rosencbrock function) and stochastic gradient descent. By the latter I don't mean it in the ML sense, I mean it in adding randomness to potentially avoid getting stuck at minima, but add the noise such that in expectation it is the true gradient. This is like optimizing a function that has been convolved with the probability density of the randomness. You decrease the randomness in later iterations.</li>
</ul>
</li>
<li>If your function is totally differentiable and you have a point on the boundary then you could try to construct new points on the surface starting from it. This is similar to <a href="https://en.wikipedia.org/wiki/Numerical_continuation" rel="nofollow noreferrer">numerical continuation</a>. You generate some points in the tangent space (the tangent space at <span class="math-container">$x$</span> is <span class="math-container">$T_xg^{-1}[0] =\{y\,|\,\nabla g(x)^T(y-x) = 0\}$</span>) and then project them down to the surface (e.g. using Newton/gradient descent again). And you repeat that process recursively.
<ul>
<li>Pros: More efficient than just Newton applied to random points.</li>
<li>Cons: This is a rather sequential approach. But you can run several realizations simultaneously starting from different points.</li>
<li>Modifications: How far the points are generated in different directions on the tangent space should ideally depend on the curvature, but you need the second derivatives for that. You can also just generate points in a ball around <span class="math-container">$x$</span> instead of in the tangent space if that is too tedious.</li>
</ul>
</li>
<li>If your function is not differentiable but it is locally Lipschitz continuous (a standard setting for SDBFs) then you can use the Clarke subgradient instead of the gradient. You can find an illustration of this as well as automatic differentiation here <a href="https://www.shadertoy.com/view/cs3XRB" rel="nofollow noreferrer">https://www.shadertoy.com/view/cs3XRB</a> and here <a href="https://www.shadertoy.com/view/4dVGzw" rel="nofollow noreferrer">https://www.shadertoy.com/view/4dVGzw</a>
<ul>
<li>Pros: Works even if your function is not totally differentiable (e.g. if you use <span class="math-container">$\min$</span>, <span class="math-container">$\max$</span>, <span class="math-container">$|\cdot|$</span>). For more details see: <a href="https://arxiv.org/abs/2006.14901" rel="nofollow noreferrer">&quot;Understanding Notions of Stationarity in Non-Smooth Optimization&quot;</a> by Li et al.</li>
<li>Cons: Convergence is actually not guaranteed unless you have some regularity assumptions. Probably nothing much to worry about in your setting though.</li>
</ul>
</li>
<li>You can use finite differences to approximate the gradient and use that for the &quot;gradient&quot; descent. Even if the function is not totally differentiable you could also just generate points in the neighbourhood of the current one and check at which the function decreases, and use that. The latter is like using directional derivatives for the descent.
<ul>
<li>Pros: You don't need to know anything except the function <span class="math-container">$g$</span>, so no automatic differentiation or analytical derivations of the gradient required.</li>
<li>Cons: If you don't pick your <span class="math-container">$\epsilon$</span> in the finite differences appropriately it can have a strong negative effect. You require multiple function evaluations for estimating the gradient.</li>
</ul>
</li>
</ul>
<hr />
<p><strong>0-th Order Oracle</strong></p>
<p>Suppose you have a function <span class="math-container">$g:\mathbb{R}^n\to\mathbb{R}$</span> that you can only evaluate the value of at arbitrary points (i.e. a zero-th order oracle) and you want to construct a point cloud for the set where <span class="math-container">$g$</span> is negative:
<span class="math-container">\begin{align}
\Omega &amp;= \{x\in\mathbb{R}^n\,|\, g(x) &lt; 0\}.
\end{align}</span>
Note - here I don't even assume that <span class="math-container">$g^{-1}[0] = \partial\Omega$</span>.</p>
<p>With only access to <span class="math-container">$g$</span> and without any additional assumptions you can perform rejection sampling. For instance you could sample points from the normal distribution <span class="math-container">$x_i \sim \mathcal{N}(\mu, \sigma^2 I)$</span>. If <span class="math-container">$g(x_i)\leq 0$</span> you accept the point <span class="math-container">$x_i$</span>, otherwise you reject it. If <span class="math-container">$\Omega$</span> has a non-zero <span class="math-container">$n$</span>-volume, then in infinite time you will be able to (non-uniformly) fill <span class="math-container">$\Omega$</span> with points up to sets of measure zero, since the normal distribution has infinite support (so it doesn't matter where your set is situated). Of course, in practice you don't have infinite time and the convergence speed depends on the set <span class="math-container">$\Omega$</span> and how you choose <span class="math-container">$\mu$</span> and <span class="math-container">$\sigma^2$</span>. To maximize <span class="math-container">$P(\mu,\sigma^2;\Omega)$</span> you can take the derivatives w.r.t. <span class="math-container">$\mu$</span> and <span class="math-container">$\theta$</span> and set the resulting expressions to zero, then you get:</p>
<p><span class="math-container">$$\mu = \frac{\int_{\Omega} x p(\mu,\sigma^2;x)\,dx}{\int_{\Omega} p(\mu,\sigma^2;x)\,dx}, \quad \sigma^2 = \frac{\int_{\Omega} \|x-\mu\|^2 p(\mu,\sigma^2;x)\,dx}{\int_{\Omega} p(\mu,\sigma^2;x)\,dx}.$$</span></p>
<p>In general if your probability density depends smoothly on some parameter <span class="math-container">$\theta$</span> and <span class="math-container">$\partial_{\theta} p(\theta;x) = C(\theta)(g(\theta;x)-\theta)p(\theta;x)$</span> then the <span class="math-container">$\theta$</span> that maximizes <span class="math-container">$P(\theta; \Omega)$</span> satisfies:</p>
<p><span class="math-container">$$\theta = \frac{\int_{\Omega}g(\theta;x)p(\theta;x)\,dx}{\int_{\Omega}p(\theta;x)\,dx}.$$</span></p>
<p>Note that even if the above gives you the optimal <span class="math-container">$\mu,\sigma^2$</span> (or in the general case <span class="math-container">$\theta$</span>), they are expressed in terms of integral equations that depend on <span class="math-container">$\Omega$</span>. If you had a uniform distribution however, then <span class="math-container">$p$</span> is constant, and the expression for the optimal mean reduces to the barycenter of the set:</p>
<p><span class="math-container">$$\mu = \frac{\int_{\Omega} x p\,dx}{\int_{\Omega} p\,dx} = \frac{\int_{\Omega} x\,dx}{|\Omega|}.$$</span></p>
<p>So if you have further details on the approximate location of the barycenter and some approximate radius you can generate points uniformly in the ball with center <span class="math-container">$\mu$</span> and with the assumed radius. To generate uniform points in a <span class="math-container">$n$</span>-D ball you can first generate points from the uniform distribution, then normalize to get uniform points on the sphere, and then rescale them randomly in [0,R] with radial density proportional to <span class="math-container">$r^{n-1}$</span>.</p>
<p>The runtime would depend on how expensive <span class="math-container">$g$</span> is to evaluate at a point and the probability that a point ends up in the set <span class="math-container">$\Omega$</span>.
If <span class="math-container">$\Omega$</span> is &quot;regular enough&quot; you could get faster convergence by also replacing the PRNG uniform points in <span class="math-container">$[0,1]^n$</span>, that you warp using your probability density, with quasi-random points (e.g. Halton, Sobol, or other low-discrepancy sequences).</p>
<p>If the set is regular enough and you know a bounding box, you could also use a regular or adaptive grid, or an octree. Marching cubes and friends take such an approach. The looser your bounding box estimate, the better adaptive structures would perform compared to a regular grid.</p>
<p>The &quot;fuzzy&quot; boundary can then be estimated as the set of points than have about twice fewer neighbours than the average in some radius <span class="math-container">$R$</span> (the larger this radius is the fewer points you need to have introduced, but the fuzzier your boundary). Note that this approach works regardless of whether <span class="math-container">$g$</span> is a signed distance function (SDF), a signed distance bound function (SDBF), or just a density that specifies the interior of the set through <span class="math-container">$g(x)&lt;0$</span>.</p>
<hr />
<p><strong>Differentiable Functions: First Order Oracle</strong></p>
<p>Suppose that <span class="math-container">$g:\mathbb{R}^n\to\mathbb{R}$</span> is differentiable and that the boundary that you are looking for is given as its zero set: <span class="math-container">$$\partial\Omega = g^{-1}[0] = \{x\in\mathbb{R}^n\,:\, g(x) = 0\}.$$</span></p>
<p>Note that for general functions <span class="math-container">$\partial\Omega \ne g^{-1}[0]$</span> since <span class="math-container">$\partial\Omega$</span> is the boundary for <span class="math-container">$\Omega = \{x\in\mathbb{R}^n\,|\,g(x)&lt;0\}$</span>, while for arbitrary functions the set <span class="math-container">$g^{-1}[0]$</span> does not have to be of measure zero, nor do points from it require neighbours that are negative (e.g. you could have a purely positive function with zeroes somewhere). So already assuming that <span class="math-container">$g^{-1}[0]$</span> gives your surface is pretty strong, and assuming that <span class="math-container">$g$</span> is totally differentiable is also pretty strong.</p>
<p>A standard approach to find the roots of differentiable functions is Newton's method. You start with some initial guess <span class="math-container">$x^0$</span> and proceed as follows:</p>
<p><span class="math-container">$$x^{k+1} = x^k - \theta^k (J_g(x^k))^+g(x^k).$$</span></p>
<p>In your case the function is scalar so the Moore-Penrose inverse of the Jacobian is the gradient divided by its squared magnitude:</p>
<p><span class="math-container">$$x^{k+1} = x^k - \frac{\theta^k g(x^k)}{\|\nabla g(x^k)\|^2}\nabla g(x^k).$$</span></p>
<p>Note that Newton can get stuck near local minima of your function similar to gradient descent. Even if the function monotonically decreasing (in 1D) if Newton is started too far away from where the function has its roots then <a href="https://math.stackexchange.com/questions/3279513/convergence-of-newton-method-for-monotonic-polynomials">it may still not converge</a>. In the latter case one may damp the iteration by setting <span class="math-container">$\theta^k$</span> to something smaller than <span class="math-container">$1$</span>, which would guarantee convergence for monotone functions. It will not guarantee convergence if you're stuck in a local minimum however (taking smaller steps won't help there, taking larger ones will). You can equivalently interpret the Newton iteration as gradient descent for the following minimization problem:</p>
<p><span class="math-container">$$\min_{x\in\mathbb{R}^n}\frac{1}{2}\|g(x)\|^2_2 \implies \nabla \frac{1}{2}\|g(x)\|^2_2 = g(x)\nabla g(x).$$</span></p>
<p>Thus the Newton iteration is equivalent to gradient descent with step size <span class="math-container">$\alpha^k = \frac{\theta^k}{\|\nabla g(x^k)\|^2_2}.$</span></p>
<p>If the computation of the Jacobian is very expensive or hard to compute one typically considers finite difference approximations of it (the issue is choosing the step size for the finite difference). In your case this is just the finite difference approximation of the gradient. If the updates are also done iteratively from the previous approximations then you get quasi-Newton methods such as Broyden's method.</p>
<hr />
<p><strong>Non-Differentiability</strong></p>
<p>If your sdf is non-differentiable then maybe subgradient descent may be useful.</p>
<p>Here are some examples of automatic differentiation routines in shadertoy:</p>
<ul>
<li><a href="https://www.shadertoy.com/view/cs3XRB" rel="nofollow noreferrer">https://www.shadertoy.com/view/cs3XRB</a></li>
<li><a href="https://www.shadertoy.com/view/4dVGzw" rel="nofollow noreferrer">https://www.shadertoy.com/view/4dVGzw</a></li>
<li><a href="https://www.shadertoy.com/view/Mdl3Ws" rel="nofollow noreferrer">https://www.shadertoy.com/view/Mdl3Ws</a></li>
</ul>
<p>You will notice that they don't care too much for the non-differentiable point in max and min, and in the absolute value they just pick the subgradient <span class="math-container">$0$</span>. Even if that may not be the most theoretically ideal thing, you can see that it works just fine in practice. I assume this is due to the nature of the problem, where the functions are continuously differentiable on a dense subset.</p>
<p>Just for completeness they use the standard calculus rules, with something extra for non-differentiable functions like <span class="math-container">$\min$</span>, <span class="math-container">$\max$</span>, and <span class="math-container">$|\cdot|$</span>:</p>
<ul>
<li><span class="math-container">$\nabla (f+g)(x) = \nabla f(x) + \nabla g(x)$</span></li>
<li><span class="math-container">$\nabla (fg)(x) = g(x)\nabla f(x) + f(x)\nabla g(x)$</span></li>
<li><span class="math-container">$\nabla (f/g)(x) = (\nabla f(x))/g(x) - f(x) (\nabla g(x))/g^2(x)$</span></li>
<li><span class="math-container">$\nabla (f\circ g)(x) = J_g(x)^T\nabla f(g(x))$</span></li>
<li><span class="math-container">$\nabla \min(f,g)(x) := \begin{cases} \nabla f(x), &amp;f(x)&lt;g(x), \\ \nabla g(x), &amp;f(x)\geq g(x)\end{cases}$</span></li>
<li><span class="math-container">$\nabla \max(f,g)(x) := \begin{cases} \nabla f(x), &amp;f(x)&gt;g(x), \\ \nabla g(x), &amp;f(x)\leq g(x)\end{cases}$</span></li>
<li><span class="math-container">$\nabla |f|(x) := \text{sgn}(f(x)) \nabla f(x)$</span></li>
</ul>
<p>So they pretend for the most part that the functions are differentiable, except for <span class="math-container">$|x|$</span> where the derivative at zero is set to zero. The choices above for <span class="math-container">$\min$</span>, <span class="math-container">$\max$</span>, and <span class="math-container">$|\cdot|$</span> at the non-differentiable poinrt are consistent with a specific choice of a subgradient at that location. There are other options, but I guess the above are the most straightforward. You do not have rock-solid guarantees with the above, but it seems to work just fine in the slinked shadertoys.</p>
<hr />
<p><strong>A Bit Of Theory</strong></p>
<hr />
<p><strong>Lipschitz Continuous Functions and Signed Distance Bound Functions</strong></p>
<p>Suppose <span class="math-container">$S$</span> is a closed surface, i.e., <span class="math-container">$S = \partial\Omega$</span> for some set <span class="math-container">$\Omega$</span> (this is the &quot;inside&quot; of the SDF). Define <span class="math-container">$d(x,S) = \inf_{y\in S}d(x,y)$</span> for some metric <span class="math-container">$d:\mathbb{R}^n\times\mathbb{R}^n\to[0,\infty)$</span>. The function <span class="math-container">$f:\mathbb{R}^n\to\mathbb{R}$</span> is a signed distance function:
<span class="math-container">$$f(x) = 
\begin{cases} -d(x,S), &amp;x\in\Omega,\\ 
d(x,S), &amp;x\not\in\Omega.\end{cases}$$</span></p>
<p>A signed distance bound function (SDBF) <span class="math-container">$g:\mathbb{R}^n\to\mathbb{R}$</span> for <span class="math-container">$\Omega$</span> is a function such that <span class="math-container">$g^{-1}[0] = S$</span> but unlike <span class="math-container">$f$</span> it underestimates the true distance:</p>
<p><span class="math-container">$$\begin{cases} 0\leq g(x) \leq f(x) = d(x,S), &amp;x\in\Omega, \\
0\leq -g(x) \leq -f(x) = d(x,S), &amp;x\not\in\Omega.\end{cases}$$</span></p>
<p>Note that just with the above definition the SDBF can be quite pathological. E.g. consider the sphere in 1D <span class="math-container">$\mathbb{S}^0 = \{-1,+1\}$</span>: <span class="math-container">$f(x) = |x|-1$</span>. Now consider an SDBF that satisfies the above definition:</p>
<p><span class="math-container">$$g(x) = 
\begin{cases}
|x|-1, &amp; x&lt;3, \\
2+\sin(\pi+2\pi (x-3)), &amp; x\geq 3. 
\end{cases}$$</span></p>
<p>You can even check that this function is Lipschitz continuous. However it is a terrible SDBF after <span class="math-container">$x=3$</span>, and would make approaches such as Newton get stuck. In practice you typically don't have such nasty settings though.</p>
<p>The metric I consider from now on is given through the Euclidean norm: <span class="math-container">$d(x,y) = \|x-y\|_2$</span>.</p>
<p>A natural candidate to construct signed distance bound functions from are Lipschitz continuous functions.
For a Lipschitz continuous function <span class="math-container">$h:\mathbb{R}^n\to\mathbb{R}$</span> there exists a constant <span class="math-container">$L(h)\geq 0$</span> such that:
<span class="math-container">$$|h(x)-h(y)|\leq L(h)\|x-y\|.$$</span>
Note that <span class="math-container">$h\equiv 0$</span> is Lipschitz continuous, with <span class="math-container">$L(h) = 0$</span>, but it is not very useful. On the other hand it does not satisfy <span class="math-container">$S=h^{-1}[0]$</span>.</p>
<p>Suppose <span class="math-container">$h$</span> satisfies <span class="math-container">$h^{-1}[0]=S$</span> and that it is Lipschitz continuous,. and it is positive on <span class="math-container">$\Omega$</span> and negative on <span class="math-container">$(\mathbb{R}^n\setminus S)\setminus \Omega$</span>. Then <span class="math-container">$g(x) = h(x) / \lambda$</span> where <span class="math-container">$\lambda \geq L(h)$</span> is a signed distance bound function for <span class="math-container">$\Omega$</span>.</p>
<p>The proof is simple and can also be found in <a href="https://link.springer.com/article/10.1007/s003710050084" rel="nofollow noreferrer">Hart's paper</a>.
Let <span class="math-container">$y\in S$</span> be such that <span class="math-container">$y\in \arg\min_{z\in S}\|x-z\|$</span>, i.e., <span class="math-container">$\|x-y\| = d(x,S)$</span>. We can show:
<span class="math-container">$$|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).$$</span></p>
<hr />
<p><strong>Standard Functions Used in Constructing an SDBF</strong></p>
<p>Typical functions used in the constructions of signed distance bound functions are usually totally differentiable, e.g., <span class="math-container">$\cos$</span>, <span class="math-container">$\sin$</span>, <span class="math-container">$\exp$</span>, <span class="math-container">$\ln$</span>, linear combinations (this includes affine transforms), polynomials, powers. Inigo Quilez even has some examples involving computation of gradients: <a href="https://iquilezles.org/articles/distgradfunctions2d/" rel="nofollow noreferrer">see this for 2D SDFs\ gradients</a>, and <a href="https://iquilezles.org/articles/morenoise/" rel="nofollow noreferrer">this for noise gradients</a>.</p>
<p>The issue would arise with functions such as absolute value, min, max, <code>clamp(x,a,b)=min(max(x,a),b)</code>, <code>mod(x,y)</code>, and worst of all functions involving a <code>floor</code> (e.g. <code>mod</code> and <code>fract</code>), <code>ceil</code>, or <code>round</code>, or potentially an <code>if</code> 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 <code>floor</code> they are not even continuous (so definitely not locally Lipschitz continuous, however their composition with other functions may be continuous and even differentiable).</p>
<p>From those <code>abs</code> and <code>max</code> are the nicest since both are subdifferentially regular, so you could use subgradient calculus with them without much issue. The <code>min</code> 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 href="https://iquilezles.org/articles/smin/" rel="nofollow noreferrer">a smooth minimum</a>. The worst offender is the <code>floor</code> 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 <code>floor</code>, 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. <code>sin(2*PI*fract(x))</code> is virtually the same as <code>sin(2*PI*x)</code> 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 <code>floor</code> is discontinuous.</p>
<p>Next I discuss some theory on nonsmooth optimization, although I feel like this is likely overkill.</p>
<hr />
<p><strong>Lipschitz Continuous Functions and Differentiation</strong></p>
<p>Lipschitz continuous functions do not have to be differentiable. On the other hand they are differentiable almost everywhere (<a href="https://en.wikipedia.org/wiki/Rademacher%27s_theorem" rel="nofollow noreferrer">see Rademacher's theorem</a>). 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 href="https://cs.nyu.edu/%7Eoverton/papers/pdffiles/gradsamp.pdf" rel="nofollow noreferrer">&quot;A ROBUST GRADIENT SAMPLING ALGORITHM FOR NONSMOOTH, NONCONVEX OPTIMIZATION&quot;</a> by Burke et al.</p>
<p>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 <a href="https://arxiv.org/abs/2006.14901" rel="nofollow noreferrer">&quot;Understanding Notions of Stationarity in Non-Smooth Optimization&quot;</a> by Li et al.</p>
<p>The Bouligand subgradient (this is often termed subdifferential, but that terminology is misleading) is defined as follows:</p>
<p><span class="math-container">$$\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\}.$$</span></p>
<p>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 <span class="math-container">$f(x) = |x|$</span>, i.e. <span class="math-container">$\partial_B f(0) = \{-1,1\}$</span>, there is the undesirable property that <span class="math-container">$0\not\in \partial_B f(0)$</span>, despite of the fact that <span class="math-container">$0$</span> is the global minimum of the above function. Remember that in basic calculus you have <span class="math-container">$\nabla f(x) = 0$</span> at extrema. To rectify this, one may consider the Clarke subgradient (see <a href="https://www.ams.org/journals/tran/1975-205-00/S0002-9947-1975-0367131-6/S0002-9947-1975-0367131-6.pdf" rel="nofollow noreferrer">&quot;GENERALIZED GRADIENTS AND APPLICATIONS&quot;</a> for the original work, but reading Li's exposition feels much nicer):</p>
<p><span class="math-container">$$\partial^{\circ} f(x) = \text{conv}(\partial_B f(x)),$$</span>
which takes the convex hull of the Bouligand subgradient. Then it follows that <span class="math-container">$\partial^{\circ} f(0) = [-1,1]$</span> for <span class="math-container">$f(x) = |x|$</span> and <span class="math-container">$0\in\partial^{\circ} f(0)$</span>. The Clarke subgradient also happens to be a generalization of the convex subgradient.</p>
<hr />
<p><strong>The Clarke (Generalized) Directional Derivative and Subgradient</strong></p>
<p>You can find the definition of the Clarke directional derivative and gradient <a href="https://en.wikipedia.org/wiki/Clarke_generalized_derivative#Definitions" rel="nofollow noreferrer">on wikipedia</a>, or in his book <a href="https://epubs.siam.org/doi/book/10.1137/1.9781611971309" rel="nofollow noreferrer">&quot;Optimization and Nonsmooth Analysis&quot;</a> (page 10):</p>
<p><span class="math-container">\begin{align}f^{\circ}(x; v) &amp;:= \lim_{y\to x}\sup_{t\downarrow 0}\frac{f(y+t v)-f(y)}{t} = \inf_{\substack{\epsilon&gt;0,\\ \delta&gt;0}}\sup_{\substack{\|y-x\|\leq\epsilon,\\ 0&lt;t&lt;\delta}}\frac{f(y+t v)-f(y)}{t},\tag{Clarke derivative} \\
\partial^{\circ} f(x) &amp;:= \{u\in\mathbb{R}^n\,:\, u^Tv\leq f^{\circ}(x; v), \,\, \forall v\in\mathbb{R}^n\} \tag{Clarke subgradient}.\end{align}</span></p>
<p>You have to actually be careful with the definitions, since for example in his other book (<a href="https://link.springer.com/book/10.1007/b97650" rel="nofollow noreferrer">&quot;Nonsmooth Analysis and Control Theory&quot;</a> 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:
<span class="math-container">$$
\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}$$</span>
The relation between the two is <span class="math-container">$\zeta(v) = u^Tv$</span> (i.e. like between the differential and the gradient: <span class="math-container">$df(x)(v) = \nabla f(x)^Tv$</span>). The subdifferential is more general and can be defined even without an inner product, but in <span class="math-container">$\mathbb{R}^n$</span> it doesn't really matter. In either case, I'll stick to the convention that <span class="math-container">$\partial^{\circ} f$</span> refers to the subgradient.</p>
<p>The Clarke (directional) derivative is a generalization of the one-sided directional derivative <span class="math-container">$\partial_v f(x)$</span> and the Clarke subgradient is a generalization of the gradient. Remember <span class="math-container">$\nabla f(x)^T v = \partial_v f(x)$</span> for totally differentiable functions, while here <span class="math-container">$u^Tv \leq f^{\circ}(x;v)$</span> for <span class="math-container">$u\in\partial^{\circ} f(x)$</span>. The Clarke derivative is positive homogeneous and sublinear. That is, if <span class="math-container">$v,w\in\mathbb{R}^n$</span> and <span class="math-container">$\alpha\geq 0$</span>  then
<span class="math-container">\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}</span></p>
<hr />
<p><strong>Clarke Calculus</strong></p>
<p>I list some properties of the Clarke gradient that hold for locally Lipschitz functions (you can find them in Chapter 2.2 of <a href="https://link.springer.com/book/10.1007/b97650" rel="nofollow noreferrer">Clarke's book</a>, or chapter 2.3 of <a href="https://epubs.siam.org/doi/book/10.1137/1.9781611971309" rel="nofollow noreferrer">his other book</a>, or some in <a href="https://arxiv.org/abs/2006.14901" rel="nofollow noreferrer">Li et al.'s exposition</a>):</p>
<ul>
<li>If <span class="math-container">$x^*$</span> is a local minimum or maximum for <span class="math-container">$f$</span> then <span class="math-container">$0\in\partial^{\circ} f(x^*)$</span> (note that the converse is not necessarily true, but at least for SDBFs you have an easy way to check whether you are at an actual minimum).</li>
<li><span class="math-container">$f\in C^1 \implies \partial^{\circ} f(x) = \{\nabla f(x)\}$</span></li>
<li><span class="math-container">$\partial^{\circ}(\lambda f)(x) = \lambda \partial^{\circ} f(x)$</span></li>
<li><span class="math-container">$\partial^{\circ}(f+g)(x) \subseteq \partial^{\circ}f(x) + \partial^{\circ}g(x)$</span></li>
<li><span class="math-container">$\partial^{\circ} (fg)(x) \subseteq f(x)\partial^{\circ} g(x) + g(x)\partial^{\circ} f(x)$</span></li>
<li><span class="math-container">$\partial^{\circ} (f/g)(x) \subseteq \frac{g(x)\partial^{\circ} f(x) - f(x)\partial^{\circ} g(x)}{g^2(x)}$</span></li>
<li><span class="math-container">$F\in C^1$</span> and <span class="math-container">$g$</span> Lipschitz, then <span class="math-container">$\partial^{\circ} (g\circ F)(x) \subseteq J_F(x)^T\partial^{\circ} g(F(x))$</span>.</li>
<li><span class="math-container">$\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))\}.$</span></li>
<li><span class="math-container">$\partial^{\circ}\max_jf_j(x) \subseteq \text{conv}\{\partial^{\circ}f_i(x)\,:\,i\in I(x)\}$</span> where <span class="math-container">$I(x)$</span> is the set of indices of functions such that <span class="math-container">$f_i(x) = \max_j f_j(x).$</span></li>
</ul>
<p>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 <span class="math-container">$f(x),g(x)\geq 0$</span> for the product rule; and <span class="math-container">$f(x)\geq 0$</span>, <span class="math-container">$g(x)&gt;0$</span>, and that <span class="math-container">$f$</span> and <span class="math-container">$-g$</span> are regular at <span class="math-container">$x$</span> for the quotient rule).</p>
<p>A function <span class="math-container">$f$</span> is regular at <span class="math-container">$x$</span> if for all <span class="math-container">$v$</span> the directional derivatives <span class="math-container">$\partial_v f(x)$</span> exist and agree with <span class="math-container">$f^{\circ}(x; v)$</span>. 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.</p>
<p>Finally here are some subgradients you may find useful:</p>
<p><span class="math-container">\begin{align}
\partial^{\circ}|x| &amp;= \begin{cases} \{1\}, &amp;x&gt;0, \\ [-1,1], &amp;x=0, \\ \{-1\}, &amp;x&lt;0.\end{cases} \\
\partial^{\circ}\max\{x,c\} &amp;= \begin{cases} \{0\}, &amp;x&lt;c, \\ [0,1], &amp;f(x)=c, \\ \{1\}, &amp;c&lt;x.\end{cases} \\
\partial^{\circ}\min\{x,c\} &amp;= \begin{cases} \{1\}, &amp;x&lt;c, \\ [0,1], &amp;x=c, \\ \{0\}, &amp;c&lt;x.\end{cases}
\end{align}</span></p>
<p>For the last two note that <span class="math-container">$\partial^{\circ} (\min{x,c} + \max{x,c})(c) = \{1\}$</span>  but <span class="math-container">$\partial^{\circ} \max\{x,c\}(c) + \partial^{\circ} \min\{x,c\}(c) = [0,2]$</span> which goes to illustrate the idea that
<span class="math-container">$$\partial^{\circ} (f+g)(x) \subseteq \partial^{\circ} f(x) + \partial^{\circ} g(x).$$</span>
Notably, here <span class="math-container">$\min$</span> is not a regular function which spoils the equality. In either case I highly recommend reading <a href="https://arxiv.org/abs/2006.14901" rel="nofollow noreferrer">&quot;Understanding Notions of Stationarity in Non-Smooth Optimization&quot;</a> 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.</p>
<p>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 <a href="https://arxiv.org/abs/1910.13604" rel="nofollow noreferrer">&quot;Pathological subgradient dynamics&quot;</a> by Daniilidis et al. If you want to see a more applied discussion for automatic differentiation see the paper <a href="https://arxiv.org/abs/1809.08530" rel="nofollow noreferrer">&quot;Provably Correct Automatic Subdifferentiation for Qualified Programs&quot;</a> and the references within.</p>

CC BY-SA 4.0
2024-11-16T16:03:09.203Z
Your answers are always mind blowing.
S. M.
CC BY-SA 4.0
2024-11-16T17:22:44.277Z
@S.M. If you take a course on numerical methods or read a book on it, it would not feel so. Newton and gradient descent are standard and very accessible concepts that you can read about even on wikipedia. Implicit surfaces as the set of zeroes of some function are also standard, especially if one already knows what an sdf is. The only not immediately obvious thing from the above is that $\|\nabla f\|_2=1$ for sdfs, but you can read about that on wikipedia also. Regarding the tangent space - any introductory material on differential geometry should do. There's nothing original in my answer.
lightxbulb
CC BY-SA 4.0
2024-11-17T06:50:16.263Z
Note that $\| \nabla f \|_2 = 1$ is not true in general for SDF's.

If you think of SDFs as solutions to the eikonal queation, then the right hand side doesn't even need to be a constant, and in practice many SDFs are not exact, e.g. the SDF for an ellipsiod.
Makogan
CC BY-SA 4.0
2024-11-17T10:36:21.037Z
@Makogan If $\|\nabla f\|\ne 1$
 it is not a signed distance function. By definition a signed distance function gives you the signed distance to the zero set. When you have a function that underestimates the distance it is typically termed a signed distance bound function, for the latter the identity does not have to hold, which is why in my answer I dicussed the case with $g$.
lightxbulb
CC BY-SA 4.0
2024-11-16T15:14:55.513Z
2024-11-20T23:21:26.577Z
2024-11-20T23:21:26.577Z
lightxbulb
<p>Note that I am well aware of dual contouring, marching cubes and friends.</p>
<p>What I am interested about is if there is an algorithm for placing point samples on the surface of the SDF, i.e. just to generate a point cloud.</p>
<p>Although lightbulb's answer is very nice mathematically, it has some issues from a practical perspective:</p>
<p>&quot;Thus you can generate a number of points&quot;</p>
<p>Generating a number of points is in <span class="math-container">$\mathbb{R}^n$</span> and of itself a challenge. The most reasonable approach would be to have a bounding box for the SDF, but I would prefer not assuming an input BB as one of the problems I need to solve is creating the BB from a blackbox function.</p>
<p>Using a newton iteration will work for many analytic or close to analytic SDFs, but SDFs in the wild are very discontinuous and have unreliable gradients.  For example in this data set:</p>
<p><a href="https://tovacinni.github.io/sdf-explorer/" rel="nofollow noreferrer">https://tovacinni.github.io/sdf-explorer/</a></p>
<p>Look at any of the fractal SDFs or even the dreads on the girl model. The gradients near those regions don't behave nicely.</p>
<p>The suggestion on subgradient descent seems promising. I will need to read further on that.</p>

CC BY-SA 4.0
2024-11-17T10:21:07.510Z
If your function is Lipschitz continuous it is differentiable almost everywhere by Rademacher's theorem. Suppose you are at a point $p$ where $f$ is not differentiable, then move by a very small random offset. The probability that you end up at a differentiable point is 100%. Then continue your gradient descent. Since you would work up to some precision $\epsilon$ the offset should not matter. As for the generation of points - if you just have a black box function I think the best you can do is trial and error, a grid, or the approach where you start on the surface.
lightxbulb
CC BY-SA 4.0
2024-11-17T10:41:03.157Z
Note that generating random points in $\mathbb{R}^n$ is trivial. For example you can generate them from a uniform distribution, or from a normal distribution. The worst that could happen is that the boundary turns out to be far away. But you could start with 1 point, let it reach the boundary, and proceed with the second approach. Or just generate normally distributed points in the tangent space around it, and again move them to the surface with gradient descent.
lightxbulb
CC BY-SA 4.0
2024-11-17T16:16:34.987Z
Also I don't get exactly what you mean by "unreliable gradients". Is the issue that the function is not differentiable? Then my above suggestion would probably still work because it should be differentiable almost everywhere if you consider Lipschitz continuous functions. Is the issue that you do not have an analytic expression for the gradient? Maybe automatic differentiation would help in this case or a finite difference approximation. Is the issue that the field is pathological like a Rosenbrock function? The non-linear conjugate gradient may help.
lightxbulb
CC BY-SA 4.0
2024-11-17T19:00:51.520Z
@lightxbulb What I mean by the gradient being problematic is that many SDF's are either fractals or exhibit fractal like behaviour. Look at the towaki data set I linked in my edit.
Makogan
CC BY-SA 4.0
2024-11-17T19:13:00.073Z
@lightxbulb Generating the points is trivial if one knows the bounding box of the shape a pirori. Which is an assumption I am trying to avoid.
Makogan
CC BY-SA 4.0
2024-11-17T19:15:44.750Z
The points are not required to be in the bounding box of the shape - it's nice if they are because then you can expect a somewhat uniform distribution and faster convergence, but there is no requirement that they need to be in a bounding box. Clearly normally distributed points are not in any box. As far as fractals are concerned - if it's a recursive iterate of a known shape you can compute the gradient. If it's something very fancy you can still use automatic differentiation or finite differences. Finite differences are not ideal, but it should work if you pick the step well.
lightxbulb
CC BY-SA 4.0
2024-11-19T00:46:08.843Z
I updated my answer with some additional information. However if you just have the sdbf as a black box function $f$ then the best you can do is estimate the gradient numerically, or construct an adaptive grid or octree such that you make use of $f$ in order to approximate the boundary better and better. Another option is to just uniformly sample points, and potentially use some metropolis algorithm to try and find the boundary. But considering that most sdbfs have an analytical gradient almost everywhere it is a pity to not use it.
lightxbulb
<rendering><computational-geometry><signed-distance-field>
Algorithm to sample an SDF?
Makogan
# Algorithm to sample an SDF?

(Asked by: username_0 on 2024-11-16)

Note that I am well aware of dual contouring, marching cubes and friends.

What I am interested about is if there is an algorithm for placing point samples on the surface of the SDF, i.e. just to generate a point cloud.

Although lightbulb's answer is very nice mathematically, it has some issues from a practical perspective:

"Thus you can generate a number of points"

Generating a number of points is in $\mathbb{R}^n$ and of itself a challenge. The most reasonable approach would be to have a bounding box for the SDF, but I would prefer not assuming an input BB as one of the problems I need to solve is creating the BB from a blackbox function.

Using a newton iteration will work for many analytic or close to analytic SDFs, but SDFs in the wild are very discontinuous and have unreliable gradients. For example in this data set:

<https://tovacinni.github.io/sdf-explorer/>

Look at any of the fractal SDFs or even the dreads on the girl model. The gradients near those regions don't behave nicely.

The suggestion on subgradient descent seems promising. I will need to read further on that.

**Tags:** rendering, computational-geometry, signed-distance-field

> **Comments:**
> **username_1**: If your function is Lipschitz continuous it is differentiable almost everywhere by Rademacher's theorem. Suppose you are at a point $p$ where $f$ is not differentiable, then move by a very small random offset. The probability that you end up at a differentiable point is 100%. Then continue your gradient descent. Since you would work up to some precision $\epsilon$ the offset should not matter. As for the generation of points - if you just have a black box function I think the best you can do is trial and error, a grid, or the approach where you start on the surface.  
> **username_1**: Note that generating random points in $\mathbb{R}^n$ is trivial. For example you can generate them from a uniform distribution, or from a normal distribution. The worst that could happen is that the boundary turns out to be far away. But you could start with 1 point, let it reach the boundary, and proceed with the second approach. Or just generate normally distributed points in the tangent space around it, and again move them to the surface with gradient descent.  
> **username_1**: Also I don't get exactly what you mean by "unreliable gradients". Is the issue that the function is not differentiable? Then my above suggestion would probably still work because it should be differentiable almost everywhere if you consider Lipschitz continuous functions. Is the issue that you do not have an analytic expression for the gradient? Maybe automatic differentiation would help in this case or a finite difference approximation. Is the issue that the field is pathological like a Rosenbrock function? The non-linear conjugate gradient may help.  
> **username_0**: @username_1 What I mean by the gradient being problematic is that many SDF's are either fractals or exhibit fractal like behaviour. Look at the towaki data set I linked in my edit.  
> **username_0**: @username_1 Generating the points is trivial if one knows the bounding box of the shape a pirori. Which is an assumption I am trying to avoid.  
> **username_1**: The points are not required to be in the bounding box of the shape - it's nice if they are because then you can expect a somewhat uniform distribution and faster convergence, but there is no requirement that they need to be in a bounding box. Clearly normally distributed points are not in any box. As far as fractals are concerned - if it's a recursive iterate of a known shape you can compute the gradient. If it's something very fancy you can still use automatic differentiation or finite differences. Finite differences are not ideal, but it should work if you pick the step well.  
> **username_1**: I updated my answer with some additional information. However if you just have the sdbf as a black box function $f$ then the best you can do is estimate the gradient numerically, or construct an adaptive grid or octree such that you make use of $f$ in order to approximate the boundary better and better. Another option is to just uniformly sample points, and potentially use some metropolis algorithm to try and find the boundary. But considering that most sdbfs have an analytical gradient almost everywhere it is a pity to not use it.  

---

## 1 answers

### Answer by username_1 on 2024-11-16. Score: 4 | Accepted answer

---

**The Answer That was Here Originally**

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.

---

**Summary for Implementation Purposes**

* For a density function $g$, generate points from a normal distribution $\mathcal{N}(\mu,\sigma^2)$. Do a subsequent pass to filter out the boundary points: $x,y$ are (fuzzy) boundary points if $\|x-y\|<2R$ and $g(x)<0$, $g(y)> 0$ or $g(x)>0$, $g(y)<0$. You can refine the boundary by studying $z = (1-t)x + ty$ and performing a bisection search for $t\in[0,1]$ with the goal to tighten the bound. If $f$ is an SDF (or not too far from an exact SDF) then $|g(x)|<R$ can be used as a criterion that does not refer to other points, then you can reject all other points and you don't need two passes.
  + Pros: Makes no additional assumptions.
  + Cons: Can happen to be extremely slow.
  + Modifications: Larger $R$ allows working with fewer points but makes the boundary fuzzier. Low-discrepancy sequences may help if the sets $\{x\,|\,g(x)<0\}$ and $\{x\,|\, g(x)>0\}$ are regular enough. A good choice of $\mu,\sigma^2$ can improve the convergence speed. If you can estimate a bounding sphere/box that contains one of the above sets you can use a uniform distribution and generate points only there. In the latter case you could also use a regular grid or an adaptive structure such as an octree.
* If your function is totally differentiable you can use Newton-Raphson/gradient descent as $x^{k+1} = x^k - \frac{\theta^k g(x^k)}{\|\nabla g(x^k)\|^2\_2}\nabla g(x^k)$. For an SDF $\|\nabla g(x^k)\|\_2=1$ so the denominator disappears. If $g$ is an SDBF then using $x^{k+1} = x^k - \frac{\theta^k g(x^k)}{\|\nabla g(x^k)\|^2}\nabla g(x^k)$ may be reasonable. If the iterates oscillate too much you can damp by picking $\theta^k<1$.
  + Pros: It should be quite performant. You could initialize with the random approach from above, and then update the points by using Newton-Raphson.
  + Cons: Assumes total differentiability and requires evaluating the gradient. Newton/gradient descent does not necessarily converge if the step size is too large - you can damp through $\theta^k$. Even if you damp the step size the method may still get stuck in a local minimum for nasty densities or SDBFs.
  + Modifications: You can look into conjugate gradients (for nasty densities that may have features like the Rosencbrock function) and stochastic gradient descent. By the latter I don't mean it in the ML sense, I mean it in adding randomness to potentially avoid getting stuck at minima, but add the noise such that in expectation it is the true gradient. This is like optimizing a function that has been convolved with the probability density of the randomness. You decrease the randomness in later iterations.
* If your function is totally differentiable and you have a point on the boundary then you could try to construct new points on the surface starting from it. This is similar to [numerical continuation](https://en.wikipedia.org/wiki/Numerical_continuation). You generate some points in the tangent space (the tangent space at $x$ is $T\_xg^{-1}[0] =\{y\,|\,\nabla g(x)^T(y-x) = 0\}$) and then project them down to the surface (e.g. using Newton/gradient descent again). And you repeat that process recursively.
  + Pros: More efficient than just Newton applied to random points.
  + Cons: This is a rather sequential approach. But you can run several realizations simultaneously starting from different points.
  + Modifications: How far the points are generated in different directions on the tangent space should ideally depend on the curvature, but you need the second derivatives for that. You can also just generate points in a ball around $x$ instead of in the tangent space if that is too tedious.
* If your function is not differentiable but it is locally Lipschitz continuous (a standard setting for SDBFs) then you can use the Clarke subgradient instead of the gradient. You can find an illustration of this as well as automatic differentiation here <https://www.shadertoy.com/view/cs3XRB> and here <https://www.shadertoy.com/view/4dVGzw>
  + Pros: Works even if your function is not totally differentiable (e.g. if you use $\min$, $\max$, $|\cdot|$). For more details see: ["Understanding Notions of Stationarity in Non-Smooth Optimization"](https://arxiv.org/abs/2006.14901) by Li et al.
  + Cons: Convergence is actually not guaranteed unless you have some regularity assumptions. Probably nothing much to worry about in your setting though.
* You can use finite differences to approximate the gradient and use that for the "gradient" descent. Even if the function is not totally differentiable you could also just generate points in the neighbourhood of the current one and check at which the function decreases, and use that. The latter is like using directional derivatives for the descent.
  + Pros: You don't need to know anything except the function $g$, so no automatic differentiation or analytical derivations of the gradient required.
  + Cons: If you don't pick your $\epsilon$ in the finite differences appropriately it can have a strong negative effect. You require multiple function evaluations for estimating the gradient.

---

**0-th Order Oracle**

Suppose you have a function $g:\mathbb{R}^n\to\mathbb{R}$ that you can only evaluate the value of at arbitrary points (i.e. a zero-th order oracle) and you want to construct a point cloud for the set where $g$ is negative:
\begin{align}
\Omega &= \{x\in\mathbb{R}^n\,|\, g(x) < 0\}.
\end{align}
Note - here I don't even assume that $g^{-1}[0] = \partial\Omega$.

With only access to $g$ and without any additional assumptions you can perform rejection sampling. For instance you could sample points from the normal distribution $x\_i \sim \mathcal{N}(\mu, \sigma^2 I)$. If $g(x\_i)\leq 0$ you accept the point $x\_i$, otherwise you reject it. If $\Omega$ has a non-zero $n$-volume, then in infinite time you will be able to (non-uniformly) fill $\Omega$ with points up to sets of measure zero, since the normal distribution has infinite support (so it doesn't matter where your set is situated). Of course, in practice you don't have infinite time and the convergence speed depends on the set $\Omega$ and how you choose $\mu$ and $\sigma^2$. To maximize $P(\mu,\sigma^2;\Omega)$ you can take the derivatives w.r.t. $\mu$ and $\theta$ and set the resulting expressions to zero, then you get:

$$\mu = \frac{\int\_{\Omega} x p(\mu,\sigma^2;x)\,dx}{\int\_{\Omega} p(\mu,\sigma^2;x)\,dx}, \quad \sigma^2 = \frac{\int\_{\Omega} \|x-\mu\|^2 p(\mu,\sigma^2;x)\,dx}{\int\_{\Omega} p(\mu,\sigma^2;x)\,dx}.$$

In general if your probability density depends smoothly on some parameter $\theta$ and $\partial\_{\theta} p(\theta;x) = C(\theta)(g(\theta;x)-\theta)p(\theta;x)$ then the $\theta$ that maximizes $P(\theta; \Omega)$ satisfies:

$$\theta = \frac{\int\_{\Omega}g(\theta;x)p(\theta;x)\,dx}{\int\_{\Omega}p(\theta;x)\,dx}.$$

Note that even if the above gives you the optimal $\mu,\sigma^2$ (or in the general case $\theta$), they are expressed in terms of integral equations that depend on $\Omega$. If you had a uniform distribution however, then $p$ is constant, and the expression for the optimal mean reduces to the barycenter of the set:

$$\mu = \frac{\int\_{\Omega} x p\,dx}{\int\_{\Omega} p\,dx} = \frac{\int\_{\Omega} x\,dx}{|\Omega|}.$$

So if you have further details on the approximate location of the barycenter and some approximate radius you can generate points uniformly in the ball with center $\mu$ and with the assumed radius. To generate uniform points in a $n$-D ball you can first generate points from the uniform distribution, then normalize to get uniform points on the sphere, and then rescale them randomly in [0,R] with radial density proportional to $r^{n-1}$.

The runtime would depend on how expensive $g$ is to evaluate at a point and the probability that a point ends up in the set $\Omega$.
If $\Omega$ is "regular enough" you could get faster convergence by also replacing the PRNG uniform points in $[0,1]^n$, that you warp using your probability density, with quasi-random points (e.g. Halton, Sobol, or other low-discrepancy sequences).

If the set is regular enough and you know a bounding box, you could also use a regular or adaptive grid, or an octree. Marching cubes and friends take such an approach. The looser your bounding box estimate, the better adaptive structures would perform compared to a regular grid.

The "fuzzy" boundary can then be estimated as the set of points than have about twice fewer neighbours than the average in some radius $R$ (the larger this radius is the fewer points you need to have introduced, but the fuzzier your boundary). Note that this approach works regardless of whether $g$ is a signed distance function (SDF), a signed distance bound function (SDBF), or just a density that specifies the interior of the set through $g(x)<0$.

---

**Differentiable Functions: First Order Oracle**

Suppose that $g:\mathbb{R}^n\to\mathbb{R}$ is differentiable and that the boundary that you are looking for is given as its zero set: $$\partial\Omega = g^{-1}[0] = \{x\in\mathbb{R}^n\,:\, g(x) = 0\}.$$

Note that for general functions $\partial\Omega \ne g^{-1}[0]$ since $\partial\Omega$ is the boundary for $\Omega = \{x\in\mathbb{R}^n\,|\,g(x)<0\}$, while for arbitrary functions the set $g^{-1}[0]$ does not have to be of measure zero, nor do points from it require neighbours that are negative (e.g. you could have a purely positive function with zeroes somewhere). So already assuming that $g^{-1}[0]$ gives your surface is pretty strong, and assuming that $g$ is totally differentiable is also pretty strong.

A standard approach to find the roots of differentiable functions is Newton's method. You start with some initial guess $x^0$ and proceed as follows:

$$x^{k+1} = x^k - \theta^k (J\_g(x^k))^+g(x^k).$$

In your case the function is scalar so the Moore-Penrose inverse of the Jacobian is the gradient divided by its squared magnitude:

$$x^{k+1} = x^k - \frac{\theta^k g(x^k)}{\|\nabla g(x^k)\|^2}\nabla g(x^k).$$

Note that Newton can get stuck near local minima of your function similar to gradient descent. Even if the function monotonically decreasing (in 1D) if Newton is started too far away from where the function has its roots then [it may still not converge](https://math.stackexchange.com/questions/3279513/convergence-of-newton-method-for-monotonic-polynomials). In the latter case one may damp the iteration by setting $\theta^k$ to something smaller than $1$, which would guarantee convergence for monotone functions. It will not guarantee convergence if you're stuck in a local minimum however (taking smaller steps won't help there, taking larger ones will). You can equivalently interpret the Newton iteration as gradient descent for the following minimization problem:

$$\min\_{x\in\mathbb{R}^n}\frac{1}{2}\|g(x)\|^2\_2 \implies \nabla \frac{1}{2}\|g(x)\|^2\_2 = g(x)\nabla g(x).$$

Thus the Newton iteration is equivalent to gradient descent with step size $\alpha^k = \frac{\theta^k}{\|\nabla g(x^k)\|^2\_2}.$

If the computation of the Jacobian is very expensive or hard to compute one typically considers finite difference approximations of it (the issue is choosing the step size for the finite difference). In your case this is just the finite difference approximation of the gradient. If the updates are also done iteratively from the previous approximations then you get quasi-Newton methods such as Broyden's method.

---

**Non-Differentiability**

If your sdf is non-differentiable then maybe subgradient descent may be useful.

Here are some examples of automatic differentiation routines in shadertoy:

* <https://www.shadertoy.com/view/cs3XRB>
* <https://www.shadertoy.com/view/4dVGzw>
* <https://www.shadertoy.com/view/Mdl3Ws>

You will notice that they don't care too much for the non-differentiable point in max and min, and in the absolute value they just pick the subgradient $0$. Even if that may not be the most theoretically ideal thing, you can see that it works just fine in practice. I assume this is due to the nature of the problem, where the functions are continuously differentiable on a dense subset.

Just for completeness they use the standard calculus rules, with something extra for non-differentiable functions like $\min$, $\max$, and $|\cdot|$:

* $\nabla (f+g)(x) = \nabla f(x) + \nabla g(x)$
* $\nabla (fg)(x) = g(x)\nabla f(x) + f(x)\nabla g(x)$
* $\nabla (f/g)(x) = (\nabla f(x))/g(x) - f(x) (\nabla g(x))/g^2(x)$
* $\nabla (f\circ g)(x) = J\_g(x)^T\nabla f(g(x))$
* $\nabla \min(f,g)(x) := \begin{cases} \nabla f(x), &f(x)<g(x), \\ \nabla g(x), &f(x)\geq g(x)\end{cases}$
* $\nabla \max(f,g)(x) := \begin{cases} \nabla f(x), &f(x)>g(x), \\ \nabla g(x), &f(x)\leq g(x)\end{cases}$
* $\nabla |f|(x) := \text{sgn}(f(x)) \nabla f(x)$

So they pretend for the most part that the functions are differentiable, except for $|x|$ where the derivative at zero is set to zero. The choices above for $\min$, $\max$, and $|\cdot|$ at the non-differentiable poinrt are consistent with a specific choice of a subgradient at that location. There are other options, but I guess the above are the most straightforward. You do not have rock-solid guarantees with the above, but it seems to work just fine in the slinked shadertoys.

---

**A Bit Of Theory**

---

**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 "inside" 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}$$

Note that just with the above definition the SDBF can be quite pathological. E.g. consider the sphere in 1D $\mathbb{S}^0 = \{-1,+1\}$: $f(x) = |x|-1$. Now consider an SDBF that satisfies the above definition:

$$g(x) =
\begin{cases}
|x|-1, & x<3, \\
2+\sin(\pi+2\pi (x-3)), & x\geq 3.
\end{cases}$$

You can even check that this function is Lipschitz continuous. However it is a terrible SDBF after $x=3$, and would make approaches such as Newton get stuck. In practice you typically don't have such nasty settings though.

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](https://link.springer.com/article/10.1007/s003710050084).
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](https://iquilezles.org/articles/distgradfunctions2d/), and [this for noise gradients](https://iquilezles.org/articles/morenoise/).

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](https://iquilezles.org/articles/smin/). 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](https://en.wikipedia.org/wiki/Rademacher%27s_theorem)). 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"](https://cs.nyu.edu/%7Eoverton/papers/pdffiles/gradsamp.pdf) 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"](https://arxiv.org/abs/2006.14901) 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"](https://www.ams.org/journals/tran/1975-205-00/S0002-9947-1975-0367131-6/S0002-9947-1975-0367131-6.pdf) 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](https://en.wikipedia.org/wiki/Clarke_generalized_derivative#Definitions), or in his book ["Optimization and Nonsmooth Analysis"](https://epubs.siam.org/doi/book/10.1137/1.9781611971309) (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"](https://link.springer.com/book/10.1007/b97650) 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](https://link.springer.com/book/10.1007/b97650), or chapter 2.3 of [his other book](https://epubs.siam.org/doi/book/10.1137/1.9781611971309), or some in [Li et al.'s exposition](https://arxiv.org/abs/2006.14901)):

* If $x^\*$ is a local minimum or maximum for $f$ then $0\in\partial^{\circ} f(x^\*)$ (note that the converse is not necessarily true, but at least for SDBFs you have an easy way to check whether you are at an actual minimum).
* $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\}(c) + \partial^{\circ} \min\{x,c\}(c) = [0,2]$ 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"](https://arxiv.org/abs/2006.14901) 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"](https://arxiv.org/abs/1910.13604) by Daniilidis et al. If you want to see a more applied discussion for automatic differentiation see the paper ["Provably Correct Automatic Subdifferentiation for Qualified Programs"](https://arxiv.org/abs/1809.08530) and the references within.

> **Comments:**
> **username_2**: Your answers are always mind blowing.  
> **username_1**: @S.M. If you take a course on numerical methods or read a book on it, it would not feel so. Newton and gradient descent are standard and very accessible concepts that you can read about even on wikipedia. Implicit surfaces as the set of zeroes of some function are also standard, especially if one already knows what an sdf is. The only not immediately obvious thing from the above is that $\|\nabla f\|_2=1$ for sdfs, but you can read about that on wikipedia also. Regarding the tangent space - any introductory material on differential geometry should do. There's nothing original in my answer.  
> **username_0**: Note that $\| \nabla f \|_2 = 1$ is not true in general for SDF's.  If you think of SDFs as solutions to the eikonal queation, then the right hand side doesn't even need to be a constant, and in practice many SDFs are not exact, e.g. the SDF for an ellipsiod.  
> **username_1**: @username_0 If $\|\nabla f\|\ne 1$  it is not a signed distance function. By definition a signed distance function gives you the signed distance to the zero set. When you have a function that underestimates the distance it is typically termed a signed distance bound function, for the latter the identity does not have to hold, which is why in my answer I dicussed the case with $g$.  

---

