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. 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" 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. 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:
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.
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, and this for noise gradients.
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. 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). 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" 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" 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" 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, or in his book "Optimization and Nonsmooth Analysis" (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" 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, or chapter 2.3 of his other book, or some in Li et al.'s exposition):
- 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" 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" 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" and the references within.
