Here the finite family \(h_k\)  of centroids in \(\mathbb {R}^{d}\)  together with the uniform
length-scale \(\sigma _\nu \)  fully determine \(\nu _h\) , thus greatly reducing the
complexity of sampling. We defer sampling from function space to future work.
Defining a density \(p_{sm} \propto e^{- S_{sm}^\nu }\)  over the space of triples
\((\nu ,\,\mathbf {q}_\nu ,\,\mathbf {p}_\nu )\)  leads to the preconditioned
Crank-Nicholson algorithm REF , see e.g. [1]}.
[h!]

MCMC for selective metamorphosis
mcmcSM\(N\) , \(K\) , \(\mathbf {q}_0\) , \(\mathbf {q}_1\) , \(\beta \in (0,1]\) 
\(j \leftarrow 1\) 
\(\nu ^j \leftarrow \text{initial guess in } \mathbb {R}^{d\times K}\) 
Solve (REF ) with \(\nu ^j\)  and \(\mathbf {q}_0,\,\mathbf {q}_1\) 
to obtain \(\omega ^j = (\mathbf {q}^j,\,\mathbf {p}^j,\, u^j)\) 
\(j < N\) 
Sample a random point \(\xi \in \mathcal {N}(0, \text{Id}_{\mathbb {R}^d})^K\) 
\(\nu \leftarrow \beta \xi + \sqrt{1-\beta ^2} \nu ^j\) 
Solve (REF ) with \(\nu \)  and \(\mathbf {q}_0,\,\mathbf {q}_1\) 
to obtain \(\omega = (\mathbf {q},\, \mathbf {p},\, u)\) 
randomUnit()\(\,< \min (1,\, e^{- S_{sm}^{\nu ^j}(\omega ^j) + S_{sm}^{\nu }(\omega )})\) 
\(\nu ^{j+1} \leftarrow \nu \) 
\(\omega ^{j+1} \leftarrow \omega \) 
\(\nu ^{j+1} \leftarrow \nu ^j\) 
\(\omega ^{j+1} \leftarrow \omega ^j\) 
\(j\leftarrow j+1\) 
\(\lbrace \nu ^j,\, \omega ^j\rbrace _{j=1}^N\)
