S2Fun.optimalSample

optimal discrete sampling points of a spherical density function

Description

optimalSample behaves similarly as discreteSample with the difference, that the sampling points are optimized to reproduce the input density function as perfect as possible. The price you pay is time. optimalSample computes the sampling points by solving a minimization problem, which becomes harder the more points you want to generate.

The density function is represented by the discrete measure

\[ \mu = \lambda \, \sum_{j=1}^M c_j \, \delta_{v_j}, \qquad \lambda = \int_{S^2} f(v) \,dv, \]

where \(f\) may be any nonnegative function - it does not have to be normalized, since scaling \(f\) multiplies the functional below by \(\lambda^2\) and leaves both the optimal directions and the optimal weights exactly where they are.

The directions \(v_j\) are moved such that the kernel discrepancy

\[ J(v,c) = \| \mu - f \|_{\psi}^2 = \sum_{n=1}^N \frac{4\pi A_n}{2n+1} \sum_{k=-n}^{n} | \hat{\mu}_n^{k} - \hat{f}_n^{k} |^2 \]

becomes small. Here \(A_n\) are the Legendre coefficients of the restricted distance kernel and \(N\) is the bandwidth, i.e. only the harmonic degrees up to \(N\) are taken into account - choose it to match the intended use of the points. The minimization is done by a limited memory BFGS iteration along the geodesics of the sphere together with an Armijo line search.

Why L-BFGS and not gradient descent. \(J\) is badly conditioned and the convergence rate of gradient descent degrades with exactly that condition number, while collecting the curvature of the steps already taken costs nothing beyond a few vectors of memory. In contrast to the rotation group the sphere is no group, so the tangent plane turns with the node and the stored steps have to be parallel transported along every step before they may be used again.

Both descent methods are available and are selected by method,

v = optimalSample(sF,n,'method','lbfgs')            % the default
v = optimalSample(sF,n,'method','steepestDescent')

where 'steepestDescent' is the plain gradient descent described above.

Note that the sum above starts at degree 1. The restricted distance kernel is only conditionally positive definite, its Legendre coefficient of degree 0 being negative, \(A_0 = -\frac43 < 0\). Since \(Y_0^0\) is constant, the degree 0 term of \(\mu - f\) equals \(\lambda (\sum_j c_j - 1)/\sqrt{4\pi}\) and hence vanishes for weights that sum up to 1 - which is why it is simply dropped.

Optimizing the weights. By default the weights are fixed, \(c_j = 1/M\), and the directions are the only unknowns. If the weights are asked for as a second output

[v,c] = optimalSample(sF,n)

they are optimized alongside the directions. Since a fixed set of directions gains an additional \(M\) degrees of freedom this way, far fewer points are needed to represent \(f\) up to a given accuracy. The resulting weights are volume fractions and may be passed on directly, e.g. by calcDensity(v,'weights',c).

The weights have to stay on the probability simplex, i.e.

\[ c_j \ge 0, \qquad \sum_{j=1}^M c_j = 1 . \]

This is not merely a cosmetic normalization. Keeping the weights there removes the only degree with a negative kernel coefficient, see above - without the constraint \(\sum_j c_j = 1\) the functional would be unbounded from below and the iteration would simply inflate the weights.

Both constraints are met by construction, since the weights enter through a softmax, \(c_j = \exp(z_j) / \sum_k \exp(z_k)\). The unknowns are then \((v,z)\) without any constraint left, and directions and weights are moved by one L-BFGS iteration - parallel transport concerns the directions, the softmax block is euclidean. The gradient with respect to the weights,

\[ \frac{\partial J}{\partial c_j} = 2 \lambda \, (\psi * (\mu - f))(v_j), \]

is the value of the very convolution whose gradient moves the directions and costs one additional transform per iteration.

The first warmUp iterations move the points only. Started on a grid that ignores the density, e.g. equispacedS2Grid, the weights would otherwise concentrate on the few points that happen to lie well before the points ever had a chance to move. Letting the points settle first is measurably better: on abs(S2Fun.smiley) with 100 points and bandwidth 32 the warm up buys a factor 2.5 in \(J\).

The weights help most where the points are far from optimal. On abs(S2Fun.smiley) with 200 points and bandwidth 64 the weights alone reduce \(J\) by a factor of 16 on the raw discreteSample nodes, while on nodes that are already driven to convergence they only add a few percent - the points have then essentially solved the problem on their own. The benefit of the weights is therefore that a given accuracy is reached with fewer points and fewer iterations, not that the converged optimum is much lower.

Note that the weights buy their accuracy up to bandwidth partly at the expense of the higher degrees, since concentrating the mass on fewer points makes the sample less uniform.

For more details, see

Gräf, Manuel; Potts, Daniel; Steidl, Gabriele (2012). Quadrature Errors, Discrepancies, and Their Relations to Halftoning on the Torus and the Sphere. SIAM Journal on Scientific Computing, 34(5), A2760–A2791. doi:10.1137/100814731

Knezevic, Marko; Landry, Nicholas W. (2015). Procedures for reducing large datasets of crystal orientations using generalized spherical harmonics. Mechanics of Materials, 88, 73–86. doi:10.1016/j.mechmat.2015.04.014

Syntax

v = optimalSample(sF,n)
v = optimalSample(sF,v)
v = optimalSample(sF,n,'bandwidth',128)
v = optimalSample(sF,n,'method','steepestDescent')
v = optimalSample(sF,n,'maxIter',1000,'tol',0.005*degree)

[v,c] = optimalSample(sF,n)                   % optimize the weights, too
[v,c] = optimalSample(sF,n,'minWeight',1e-4)

Input

Output

Options

See also

S2Fun.discreteSample S2RestrictedDistanceKernel