SO3Fun.optimalSample

optimal discrete sampling points of an orientation 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_{R_j}, \qquad \lambda = \int_{SO(3)} f® \,dR, \]

and the orientations \(R_j\) are moved such that the kernel discrepancy

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

becomes small. Here \(A_n\) are the Chebyshev 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 orientations. The minimization is done by a limited memory BFGS iteration with Armijo line search.

Why L-BFGS and not gradient descent. \(J\) is badly conditioned - in a minimizer of a 92 point sample of the SO3Fun.dubna ODF the Hessian has a condition number of roughly 2000 - and the convergence rate of gradient descent degrades with exactly that number. Collecting curvature from the steps already taken costs nothing beyond a few vectors of memory and buys about an order of magnitude in time to a given discrepancy. A true Newton method does not pay on top of that: the Hessian couples every pair of orientations, so it is a dense \(3M \times 3M\) matrix, and in a measurement its iteration count was no better than the one of L-BFGS.

Both descent methods are available and are selected by method,

ori = optimalSample(f,n,'method','lbfgs')            % the default
ori = optimalSample(f,n,'method','steepestDescent')

where 'steepestDescent' is the plain gradient descent described above. It is kept because it is the method the cited literature uses and because it is the reference the L-BFGS iteration is compared against; there is no reason to prefer it in production.

Note that the sum above starts at degree 1. The restricted distance kernel is only conditionally positive definite, its Chebyshev coefficient of degree 0 being negative, \(A_0 = -\frac{16\sqrt2}{3\pi} < 0\). Since \(D_0\) is constant, the degree 0 term of \(\mu - f\) equals \(\lambda (\sum_j c_j - 1)/(\sqrt{8}\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 orientations are the only unknowns. If the weights are asked for as a second output

[ori,c] = optimalSample(f,n)

they are optimized alongside the orientations. Since a fixed set of orientations gains an additional \(M\) degrees of freedom this way, far fewer orientations 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(ori,'weights',c).

The weights are then restricted to 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 on the simplex 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.

The problem is solved by alternating minimization. For fixed orientations \(J\) is a convex least squares functional in \(c\) and is decreased with mlsq, which preserves both \(\sum_j c_j\) and \(c_j \ge 0\). For fixed weights the orientations are moved as described above.

Choose the bandwidth to match the intended use. The weights buy their accuracy up to bandwidth partly at the expense of the higher degrees, since concentrating the mass on fewer orientations makes the sample less uniform. On the SO3Fun.dubna ODF with 200-300 orientations the weighted sample beats the equally weighted one in the L1 error of the reconstructed density by 4-8 percent at every kernel halfwidth between 5 and 12 degree - but only if bandwidth is chosen large enough (64 in that test). With the default of 32 the advantage is present for halfwidth 10 degree and above and turns into a disadvantage below. Note also that the improvement is in the L1 error; measured in the L2 error, which punishes the deviation of the peaks much harder, the equally weighted sample can remain the better choice at halfwidths well below the optimized bandwidth.

Starting from a grid that ignores the density - e.g. equispacedSO3Grid instead of the default discreteSample - the first weight step is drastic and kills nodes before they ever had a chance to move. Since the update of mlsq is multiplicative, a weight that reaches 0 stays 0, and since the gradient with respect to \(R_j\) carries the factor \(c_j\), such a node is frozen in place as well and is lost for good. Use warmUp to move the orientations only for the first iterations in that case.

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

ori = optimalSample(f,n)
ori = optimalSample(f,ori)
ori = optimalSample(f,n,'bandwidth',32)
ori = optimalSample(f,n,'method','steepestDescent')
ori = optimalSample(f,n,'maxIter',1000,'tol',0.05*degree)

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

Input

Output

Options

See also

SO3Fun.discreteSample mlsq SO3RestrictedDistanceKernel