A New Bayesian Approach to Global Optimization on Parametrized Surfaces in $\mathbb {R}^{3}$

Abstract

This work introduces a new Riemannian optimization method for registering open parameterized surfaces with a constrained global optimization approach. The proposed formulation leads to a rigorous theoretic foundation and guarantees the existence and the uniqueness of a global solution. We also propose a new Bayesian clustering approach where local distributions of surfaces are modeled with spherical Gaussian processes. The maximization of the posterior density is performed with Hamiltonian dynamics which provide a natural and computationally efficient spherical Hamiltonian Monte Carlo sampling. Experimental results demonstrate the efficiency of the proposed method.

Appendices

Appendix

More results

See Figs. 8 and  9

Fig. 9

The Markov chain values obtained by HMC a, the normalized histogram b, and the nonparametric density estimation c of $a^1_{1,1}$ (top) and $b^1_{1,1}$ (bottom)

Proofs

Proof of Theorem 3.1

Let $\gamma \in \varGamma$ be defined as follows

$$\begin{aligned} \gamma : [0,1]^{2}\rightarrow & {} [0,1]^{2}\\ \xi\rightarrow & {} \left( \gamma _{1}(\xi ), \gamma _{2}(\xi ) \right) , \end{aligned}$$

where $\gamma _{1}(\xi )=\xi _{1}$ and $\gamma _{2}(\lambda , \xi _{2})=\alpha ^{\lambda }(\xi _{2})$ for $\lambda \in [0,1]$. Hence, the determinant of the Jacobi of $\gamma$ is

$$\begin{aligned} \det \left( \begin{bmatrix} \frac{d\gamma _{1}}{d\xi _{1}} &{} \frac{d\gamma _{1}}{d\xi _{2}} \\ \frac{d\gamma _{2}}{d\xi _{1}} &{} \frac{d\gamma _{2}}{d\xi _{2}} \end{bmatrix} \right)= & {} \det \left( \begin{bmatrix} 1 &{} 0 \\ \frac{d\gamma _{2}}{d\xi _{2}} &{} \frac{d\gamma _{2}}{d\xi _{2}} \end{bmatrix} \right) \\ {}= & {} \frac{d \gamma _{2}}{d \xi _{2}}. \end{aligned}$$

Now let $q^{\lambda }(\xi _{2})$ and $p^{\lambda }(\xi _{2})$ denote $Q_{1}(\lambda , \xi _{2})$ and $Q_{2}(\lambda , \xi _{2})$ respectively. We have

$$\begin{aligned} h^{\lambda }(\alpha ^{\lambda })= & {} \int _{0}^{1} \Vert Q_{1}(\lambda , \xi _{2}) - (Q_{2},\gamma )(\lambda ,\xi _{2})\Vert ^{2}_{2} d\xi _{2}, \\= & {} \int _{0}^{1} \Vert Q_{1}(\lambda , \xi _{2}) - \sqrt{|J_{\gamma }|}(Q_{2} \circ \gamma )(\lambda , \xi _{2}) \Vert _{2}^{2} d\xi _{2}, \\= & {} \int _{0}^{1} \Vert q^{\lambda }(\xi _{2}) - \sqrt{\dot{\alpha }^{\lambda }} p^{\lambda } \circ \alpha ^{\lambda }(\xi _{2}) \Vert _{2}^{2} d\xi _{2}, \end{aligned}$$

which completes the proof for $h^{\lambda }$. Likewise we obtain a proof for $k^{\lambda }$. $\square$

Proof of Theorem 3.2

Substituting $\gamma$ with its truncated version $\gamma _{m}$ obtained from the K–L expansion of $\psi _j$, we can rewrite (4) as

$$\begin{aligned} h^{\lambda }(\alpha ^{\lambda }_{2,m})= & {} {\small \int _{0}^{1} } || Q_1(\lambda ,\xi _2)-\sqrt{\int _{0}^{\lambda } {\psi _{2,m}(t_1,\xi _2)}^2 dt_1} Q_2 \circ \gamma _m(\lambda ,\xi _2) ||_2^2 d \xi _2, \nonumber \\= & {} \int _{0}^{1} || q^{\lambda }(\xi _2)-\sqrt{ \dot{\alpha }^{\lambda }_{2,m} (\xi _2) } p^{\lambda } \circ \alpha ^{\lambda }_{2,m}(\xi _2) ||_2^2 d \xi _2. \end{aligned}$$

(A1)

If $B \in \mathcal {S}^{m^2-1}$ then $\alpha ^{\lambda }_{2,m}$ is a 1D reparametrization which gives the optimal minimizer $\hat{B}$, see [15] for further details about solution of optimization problem (A1). Similar considerations apply to (5) to obtain the optimal minimizer $\hat{A}$.

Proof of Theorem 4.1

From (9), the complete log-likelihood term is

$$\begin{aligned}{} & {} \log p (\textbf{D}|A^1,\dots ,A^K,B^1,\dots ,B^K,\pi _1,\dots ,\pi _K,\tilde{Q}^{1,m}(\varvec{\xi }),\dots , \tilde{Q}^{K,m}(\varvec{\xi }), \sigma ^2 ) \nonumber \\{} & {} \propto \sum _{i=1}^{N} \log \Big ( \sum _{k=1}^K \pi _k \exp \big (- \frac{1}{2 \sigma ^2} ||Q^{*,m}_i(\varvec{\xi }) - \tilde{Q}^{k,m}(\varvec{\xi })||_2^2 \big ) \Big ). \end{aligned}$$

(A2)

From (10), we can write the constrained log-prior as

$$\begin{aligned} \log p(A^1,\dots ,A^K,B^1,\dots ,B^K) \propto - \frac{1}{2} \sum _{k=1}^{K} \sum _{l_1,l_2=1}^{m} \frac{(a_{l_1,l_2}^{k})^2+(b_{l_1,l_2}^{k})^2}{\lambda _{l_1,l_2}}, \end{aligned}$$

(A3)

under the constraint that $A^1,\dots , A^K, B^1,\dots , B^K$ belong to $\mathcal {S}^{m^2-1}$. The desired result in (11) yields by plugging (A2) and (A3) into the log-posterior probability term. $\square$

Optimal Reparametrization Between Two Distinct Curves

Let I be a certain univariate interval of parametrization and let ${{\,\mathrm{\mathbb {L}^2}\,}}(I, \mathbb {R}^{n})$ be the set of square integrable functions from I to $\mathbb {R}^{n}$. Let $\beta : I \rightarrow \mathbb {R}^{n}$ denote a parametrized curve in ${{\,\mathrm{\mathbb {L}^2}\,}}(\mathcal {D}, \mathbb {R}^{n})$, where $\mathcal {D}$ is $I=[0,1]$ for open curves. We are going to restrict to those $\beta$ satisfying: (i) $\beta$ is absolutely continuous, (ii) $\dot{\beta } \in {{\,\mathrm{\mathbb {L}^2}\,}}(\mathcal {D}, \mathbb {R}^{n})$. Note that absolute continuity is equivalent to $\dot{\beta }$ exists for almost $t \in \mathcal {D}$, that $\dot{\beta }$ is summable and that $\beta (t) =\int _{0}^{t} \dot{\beta }(s) ds$. In [41], Srivastava et al. proposed the Square Root Velocity Function (SRVF) for comparing absolutely continuous curves in $\mathbb {R}^{n}$. The SRVF representation $q: \rightarrow \mathbb {R}^{n}$ of $\beta$ is defined as follows

$$\begin{aligned} q(t) = {\left\{ \begin{array}{ll} \frac{\dot{\beta }(t)}{\sqrt{||\dot{\beta }(t)||_2}}, &{} \quad \text {if} \quad ||\dot{\beta }(t)||_2 \ne 0. \\ 0, &{}\quad \text {otherwise}. \\ \end{array}\right. } \end{aligned}$$

where $||. ||_{2}$ denotes the 2-norm in $\mathbb {R}^{n}$. Conversely, given $q \in {{\,\mathrm{\mathbb {L}^2}\,}}(I, \mathbb {R}^n)$, the original curve $\beta$ can be recovered using $\beta (t) = \beta _{0} + \int _{0}^{t} q(s)||q(s)||_2 ds$, where $\beta _{0}\in \mathbb {R}^{n}$ is the starting point of $\beta$. Let $\bar{M}_I$ denote the set of all SRVFs. It is easily seen that in this case of open curves $\bar{M}_I= {{\,\mathrm{\mathbb {L}^2}\,}}(I,\mathbb {R}^n)$. The reparameterization group for open curves in $\mathbb {R}^n$ is

$$\begin{aligned} \varGamma _{I}= \lbrace \gamma :I \rightarrow I | \gamma \; \text {is absolutely continuous}, \gamma (0)=0, \\ \gamma (1)=1,\; \text {and} \quad \dot{\gamma }(t)>0 \; \text {almost everywhere} \rbrace , \end{aligned}$$

and its action on the curve $(\beta , \gamma )= \beta \circ \gamma$ is equivalent to the action $(q,\gamma )= \sqrt{\dot{\gamma }}q \circ \gamma$. The shape space for open curves, given by the quotient: $\mathcal {M}_{I}= \lbrace [q]| q \in \bar{M}_I \rbrace$, where $[q]= \lbrace (q,\gamma )| \gamma \in \varGamma _{I} \rbrace$, is an orbit. It has been shown in [41], that the group action on q is isometric with respect to the ${{\,\mathrm{\mathbb {L}^2}\,}}$ metric. The isometry property is that if any two curves are reparameterized by the same function $\gamma$, then the resulting distance between them under the elastic metric does not change. Hence, $\mathcal {M}_{I}$ equipped with the ${{\,\mathrm{\mathbb {L}^2}\,}}$ metric is a well defined metric space. Thus, the cost function to find the optimal reparametrization between two distinct curves is,

$$\begin{aligned} C: \varGamma _{I} \rightarrow \mathbb {R}_+: \gamma\mapsto & {} || q_{1} - (q_{2},\gamma )||_{{{\,\mathrm{\mathbb {L}^2}\,}}}^2 = \int _{I} || q_{1} - \sqrt{\dot{\gamma (t)}}q_{2}\circ \gamma (t)||_{2}^{2} dt. \nonumber \\ \end{aligned}$$

(A4)

The Hilbert Sphere

We list some analytical expressions that are useful for analyzing elements of the Hilbert sphere $\mathcal {H}$.

• Geodesic path. Given $\psi \in \mathcal {H}$ and a vector $g \in T_{\psi }(\mathcal {H})$, the geodesic path with initial condition $\psi$ and velocity g at any time instant t can be parameterized as

  $$\begin{aligned} \psi (t)=\cos \big (t|| g||_{{{\,\mathrm{\mathbb {L}^2}\,}}} \big ) \psi + \sin \big (t|| g||_{{{\,\mathrm{\mathbb {L}^2}\,}}}\big ) \frac{g}{|| g||_{{{\,\mathrm{\mathbb {L}^2}\,}}}}. \end{aligned}$$

• Geodesic distance. The arc length of the geodesic path in $\mathcal {H}$ between two functions $\psi _1$ and $\psi _2$, called geodesic distance, is given by

  $$\begin{aligned} \text {dist}\big (\psi _1,\psi _2\big )_{\mathcal {H}}=\arccos \big (\big <\psi _1,\psi _2\big >_{{{\,\mathrm{\mathbb {L}^2}\,}}}\big ). \end{aligned}$$

• Exponential map. Let $\psi$ be any element of $\mathcal {H}$ and $g \in T_{\psi }(\mathcal {H})$. We define the exponential map as an isometry from $T_{\psi }(\mathcal {H})$ to $\mathcal {H}$, satisfying

  $$\begin{aligned} \exp _{\psi }(g)=\cos \big (|| g||_{{{\,\mathrm{\mathbb {L}^2}\,}}}\big ) \psi + \sin \big (|| g||_{{{\,\mathrm{\mathbb {L}^2}\,}}}\big ) \frac{g}{|| g||}_{{{\,\mathrm{\mathbb {L}^2}\,}}}. \end{aligned}$$

  The exponential map is a bijection between the tangent space and the unit sphere if we restrict $|| g||_{{{\,\mathrm{\mathbb {L}^2}\,}}}$ so that $|| g||_{{{\,\mathrm{\mathbb {L}^2}\,}}} \in [0, \pi )$.

• Log map. For $\psi _1,\psi _2 \in \mathcal {H}$, we define $\zeta \in T_{\psi _1}(\mathcal {H})$ to be the inverse exponential (log) map of $\psi _2$ if $\exp _{\psi _1}(\zeta )=\psi _2$. We use the notation $\zeta =\log _{\psi _1}(\psi _2)$ where

  $$\begin{aligned} \zeta =\frac{\alpha }{||\alpha ||}_{{{\,\mathrm{\mathbb {L}^2}\,}}} \text {dist}\big (\psi _1,\psi _2\big )_{\mathcal {H}} \; \text {and} \; \alpha =\psi _2-\psi _1\big <\psi _2,\psi _1\big >_{{{\,\mathrm{\mathbb {L}^2}\,}}} \psi _1. \end{aligned}$$