A tutorial on task-parameterized movement learning and retrieval

Abstract

Task-parameterized models of movements aim at automatically adapting movements to new situations encountered by a robot. The task parameters can, for example, take the form of positions of objects in the environment or landmark points that the robot should pass through. This tutorial aims at reviewing existing approaches for task-adaptive motion encoding. It then narrows down the scope to the special case of task parameters that take the form of frames of reference, coordinate systems or basis functions, which are most commonly encountered in service robotics. Each section of the paper is accompanied by source codes designed as simple didactic examples implemented in Matlab with a full compatibility with GNU Octave, closely following the notation and equations of the article. It also presents ongoing work and further challenges that remain to be addressed, with examples provided in simulation and on a real robot (transfer of manipulation behaviors to the Baxter bimanual robot). The repository for the accompanying source codes is available at http://www.idiap.ch/software/pbdlib/.

Appendices

Appendix 1: Expectation-maximization for TP-GMM parameters estimation

In order to estimate the parameters of a TP-GMM, the following two steps are repeated until convergence:

E-step:

$$\begin{aligned} h_{t,i} = \frac{\pi _i \prod \nolimits _{j=1}^P \mathcal {N}\Big ({\varvec{X}}^{(j)}_t \Big |\; {\varvec{\mu }}^{(j)}_i,{\varvec{{\varSigma }}}^{(j)}_i \Big )}{\sum _{k=1}^K \pi _k \prod \nolimits _{j=1}^P \mathcal {N}\Big ({\varvec{X}}^{(j)}_t \Big |\; {\varvec{\mu }}^{(j)}_k,{\varvec{{\varSigma }}}^{(j)}_k \Big )}. \end{aligned}$$

(49)

M-step:

$$\begin{aligned}&\pi _i \leftarrow \frac{\sum _{t=1}^N h_{t,i}}{N}, \end{aligned}$$

(50)

$$\begin{aligned}&{\varvec{\mu }}^{(j)}_i \leftarrow \frac{ \sum _{t=1}^N h_{t,i}\; {\varvec{X}}^{(j)}_t }{\sum _{t=1}^N h_{t,i}},\end{aligned}$$

(51)

$$\begin{aligned}&{\varvec{{\varSigma }}}^{(j)}_i \leftarrow \frac{\sum _{t=1}^N h_{t,i}\; \Big ({\varvec{X}}^{(j)}_t-{\varvec{\mu }}^{(j)}_i\Big ) \Big ({\varvec{X}}^{(j)}_t-{\varvec{\mu }}^{(j)}_i \Big )^{\scriptscriptstyle \top }}{\sum _{t=1}^N h_{t,i}}. \end{aligned}$$

(52)

In practice, it is recommended to start EM from a coarse estimate of the parameters. For example, based on an equal split in time of motion segments, based on a geometric segmentation with k-means [59], based on moments or spectral approaches with circular covariances [42, 53, 84] or based on an iterative clustering algorithm [83].

Model selection (i.e., determining the number of Gaussians in the GMM) is compatible with the techniques employed in standard GMM, such as the use of a Bayesian information criterion [82], Dirichlet process [22, 50, 65, 74], iterative pairwise replacement [83], spectral clustering [53, 69, 84] or based on segmentation points [56]. Model selection in mixture modeling shares a similar core challenge as that of data-driven sparse kernel regression techniques, which requires to find the right bandwidth parameters to select a subset of existing/new datapoints that are the most representatives of the dataset.

Appendix 2: Expectation-maximization for TP-MFA and TP-MPPCA parameters estimation

In TP-MFA, the generative model for the jth frame and ith mixture component assumes that a D-dimension random vector ${\varvec{X}}^{(j)}$ is modeled using a d-dimension vector of latent (unobserved) factors ${\varvec{z}}^{(j)}$

$$\begin{aligned} {\varvec{X}}^{(j)} = {\varvec{{\varLambda }}}^{(j)}_i {\varvec{z}}^{(j)} + {\varvec{\mu }}^{(j)}_i + {\varvec{\epsilon }}^{(j)}_i, \end{aligned}$$

(53)

where ${\varvec{\mu }}^{(j)}_i\in \mathbb {R}^D$ is the mean vector of the ith factor analyzer, ${\varvec{z}}^{(j)}\sim \mathcal {N}({\varvec{0}},{\varvec{I}})$ (the factors are assumed to be distributed according to a zero-mean normal with unit variance), and ${\varvec{\epsilon }}^{(j)}_i\sim \mathcal {N}({\varvec{0}},{\varvec{{\varPsi }}}^{(j)}_i)$ is a centered normal noise with diagonal covariance ${\varvec{{\varPsi }}}^{(j)}_i$.

This diagonality is a key assumption in factor analysis. Namely, the observed variables are independent given the factors, and the goal of TP-MFA is to best model the covariance structure of ${\varvec{X}}^{(j)}$. It follows from this model that the marginal distribution of ${\varvec{X}}^{(j)}$ for the ith component is

$$\begin{aligned} {\varvec{X}}^{(j)} \;\sim \; \mathcal {N}\left( {\varvec{\mu }}^{(j)}_i,\; {\varvec{{\varLambda }}}^{(j)}_i {{\varvec{{\varLambda }}}^{(j)}_i}^{\scriptscriptstyle \top }+ {\varvec{{\varPsi }}}^{(j)}_i \right) , \end{aligned}$$

(54)

and the joint distribution of ${\varvec{X}}^{(j)}$ and ${\varvec{z}}^{(j)}$ is

$$\begin{aligned} \left[ \begin{matrix} {\varvec{X}}^{(j)} \\ {\varvec{z}}^{(j)} \end{matrix}\right] \;\sim \; \mathcal {N}\left( \left[ \begin{matrix} {\varvec{\mu }}^{(j)}_i \\ {\varvec{0}} \end{matrix}\right] ,\; \left[ \begin{matrix} {\varvec{{\varLambda }}}^{(j)}_i {{\varvec{{\varLambda }}}^{(j)}_i}^{\scriptscriptstyle \top }+ {\varvec{{\varPsi }}}^{(j)}_i &{} \;{\varvec{{\varLambda }}}_i \\ {{\varvec{{\varLambda }}}^{(j)}_i}^{\scriptscriptstyle \top }&{} \;{\varvec{I}} \end{matrix}\right] \right) . \end{aligned}$$

(55)

The above can be used to show that the d factors are informative projections of the data, which can be computed by Gaussian conditioning, corresponding to the affine projection

$$\begin{aligned} {\varvec{z}}^{(j)} | {\varvec{X}}^{(j)} \;\sim \; \overbrace{{{\varvec{{\varLambda }}}^{(j)}_i}^{\scriptscriptstyle \top }{\Big ({\varvec{{\varLambda }}}^{(j)}_i {{\varvec{{\varLambda }}}^{(j)}_i}^{\scriptscriptstyle \top }+ {\varvec{{\varPsi }}}^{(j)}_i\Big )}^{-1}}^{{\varvec{B}}^{(j)}_i}\; \Big ({\varvec{\mu }}^{(j)}_i-{\varvec{X}}^{(j)}\Big ). \end{aligned}$$

(56)

As highlighted by [32], the same process can be used to estimate the second moment of the factors $\mathbb {E}\Big ({\varvec{z}}^{(j)}{{\varvec{z}}^{(j)}}^{\scriptscriptstyle \top }|{\varvec{X}}^{(j)}\Big )$, which provides a measure of uncertainty in the factors that has no analogue in PCA. This relation can be exploited to derive an EM algorithm (see for example [32] or [62]) to train a TP-MFA model of K components with parameters $\big \{\pi _i,\{{\varvec{\mu }}^{(j)}_i,{\varvec{{\varLambda }}}^{(j)}_i,{\varvec{{\varPsi }}}^{(j)}_i\}_{j=1}^P\big \}_{i=1}^K$, yielding an EM parameters estimation strategy.

The following two steps are repeated until convergence:

E-step:

$$\begin{aligned} h_{t,i} = \frac{\pi _i \prod \nolimits _{j=1}^P \mathcal {N}\Big ({\varvec{X}}^{(j)}_t \Big | \;{\varvec{\mu }}^{(j)}_i,\; {\varvec{{\varLambda }}}^{(j)}_i {{\varvec{{\varLambda }}}^{(j)}_i}^{\scriptscriptstyle \top }+ {\varvec{{\varPsi }}}^{(j)}_i \Big )}{\sum _{k=1}^K \pi _k \prod \nolimits _{j=1}^P \mathcal {N}\Big ({\varvec{X}}^{(j)}_t \Big | \;{\varvec{\mu }}^{(j)}_k,\; {\varvec{{\varLambda }}}^{(j)}_k {{\varvec{{\varLambda }}}^{(j)}_k}^{\scriptscriptstyle \top }+ {\varvec{{\varPsi }}}^{(j)}_k \Big )}. \end{aligned}$$

(57)

M-step:

$$\begin{aligned} \pi _i&\leftarrow \frac{\sum _{t=1}^N h_{t,i} }{N},\end{aligned}$$

(58)

$$\begin{aligned} {\varvec{\mu }}^{(j)}_i&\leftarrow \frac{\sum _{t=1}^N h_{t,i} {\varvec{X}}^{(j)}_t}{\sum _{t=1}^N h_{t,i} },\end{aligned}$$

(59)

$$\begin{aligned} {\varvec{{\varLambda }}}^{(j)}_i&\leftarrow {\varvec{S}}^{(j)}_i {{\varvec{B}}^{(j)}_i}^{\scriptscriptstyle \top }{\Big ({\varvec{I}}-{\varvec{B}}^{(j)}_i {\varvec{{\varLambda }}}^{(j)}_i + {\varvec{B}}^{(j)}_i {\varvec{S}}^{(j)}_i {{\varvec{B}}^{(j)}_i}^{\scriptscriptstyle \top }\Big )}^{-1}, \end{aligned}$$

(60)

$$\begin{aligned} {\varvec{{\varPsi }}}^{(j)}_i&\leftarrow {{\mathrm {diag}}}\Big ({{\mathrm {diag}}}\big ( {\varvec{S}}^{(j)}_i - {\varvec{{\varLambda }}}^{(j)}_i {\varvec{B}}^{(j)}_i {\varvec{S}}^{(j)}_i \big )\Big ), \end{aligned}$$

(61)

computed with the help of the intermediary variables

$$\begin{aligned} {\varvec{S}}^{(j)}_i&= \frac{\sum _{t=1}^N h_{t,i} \Big ({\varvec{X}}^{(j)}_t-{\varvec{\mu }}^{(j)}_i\Big ) {\Big ({\varvec{X}}^{(j)}_t-{\varvec{\mu }}^{(j)}_i\Big )}^{\scriptscriptstyle \top }}{\sum _{t=1}^N h_{t,i} }, \end{aligned}$$

(62)

$$\begin{aligned} {\varvec{B}}^{(j)}_i&= {{\varvec{{\varLambda }}}^{(j)}_i}^{\scriptscriptstyle \top }{\Big ({\varvec{{\varLambda }}}^{(j)}_i {{\varvec{{\varLambda }}}^{(j)}_i}^{\scriptscriptstyle \top }+ {\varvec{{\varPsi }}}^{(j)}_i\Big )}^{-1}. \end{aligned}$$

(63)

Alternatively, an update step simultaneously computing ${\varvec{\mu }}^{(j)}_i$ and ${\varvec{{\varLambda }}}^{(j)}_i$ can be derived, see [32] for details.

Similarly, the M-step in TP-MPPCA is given by

$$\begin{aligned} {\varvec{\tilde{\varLambda }}}^{(j)}_i&\leftarrow {\varvec{S}}^{(j)}_i {\varvec{{\varLambda }}}^{(j)}_i {\Big ({\varvec{I}}{\sigma ^{(j)}_i}^2 + {{\varvec{M}}^{(j)}_i}^{-1} {{\varvec{{\varLambda }}}^{(j)}_i}^{\scriptscriptstyle \top }{\varvec{S}}^{(j)}_i {\varvec{{\varLambda }}}^{(j)}_i\Big )}^{-1}, \end{aligned}$$

(64)

$$\begin{aligned} {\varvec{{\varPsi }}}^{(j)}_i&\leftarrow {\varvec{I}} {\sigma ^{(j)}_i}^2, \end{aligned}$$

(65)

computed with the help of the intermediary variables

$$\begin{aligned} {\varvec{S}}^{(j)}_i&= \frac{\sum _{t=1}^N h_{t,i} \Big ({\varvec{\xi }}^{(j)}_t - {\varvec{\mu }}^{(j)}_i\Big ) {\Big ({\varvec{\xi }}^{(j)}_t - {\varvec{\mu }}^{(j)}_i\Big )}^{\scriptscriptstyle \top }}{\sum _{t=1}^N h_{t,i} }, \end{aligned}$$

(66)

$$\begin{aligned} {\varvec{M}}^{(j)}_i&= {{\varvec{{\varLambda }}}^{(j)}_i}^{\scriptscriptstyle \top }{\varvec{{\varLambda }}}^{(j)}_i + {\varvec{I}} {\sigma ^{(j)}_i}^2, \end{aligned}$$

(67)

$$\begin{aligned} {\sigma ^{(j)}_i}^2&= \frac{1}{D} {{\mathrm {tr}}}\Big ({\varvec{S}}^{(j)}_i - {\varvec{S}}^{(j)}_i {\varvec{{\varLambda }}}^{(j)}_i {{\varvec{M}}^{(j)}_i}^{-1} {{\varvec{\tilde{\varLambda }}}^{(j)}_i}^{\scriptscriptstyle \top }\Big ), \end{aligned}$$

(68)

where ${\varvec{{\varLambda }}}^{(j)}_i$ is replaced by ${\varvec{\tilde{\varLambda }}}^{(j)}_i$ at each iteration, see [93] for details.

Appendix 3: Gaussian mixture regression approximated by a single normal distribution

Let us consider a datapoint ${\varvec{\xi }}_t$ distributed as in Eq. (6), with $\mathcal {P}({\varvec{\xi }}_t) = \mathcal {P}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t,{\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t)$ being the joint distribution describing the data. The conditional probability of an output given an input is

$$\begin{aligned} \mathcal {P}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t | {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t) = \frac{\mathcal {P}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t,{\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t)}{\mathcal {P}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t)} = \frac{\sum _{i=1}^K \mathcal {P}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t,{\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t | z_i) \mathcal {P}(z_i)}{\mathcal {P}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t)}, \end{aligned}$$

(69)

where $z_i$ represents the ith component of the GMM. Namely,

$$\begin{aligned} {\mathcal {P}}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t | {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t)&= \sum _{i=1}^K {\mathcal {P}}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t|{\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t,z_i) \frac{{\mathcal {P}}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t | z_i) {\mathcal {P}}(z_i)}{{\mathcal {P}}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t)} \nonumber \\&= \sum _{i=1}^K h_i({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t)\; \mathcal {N}\left( {\varvec{\hat{\mu }}}^{\scriptscriptstyle {{\mathcal {O}}}}_i({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t),{\varvec{\hat{{\varSigma }}}}^{\scriptscriptstyle {{\mathcal {O}}}}_i\right) . \end{aligned}$$

(70)

The conditional mean can be computed as

$$\begin{aligned} {\varvec{\hat{\mu }}}{}^{\scriptscriptstyle {{\mathcal {O}}}}_t = \mathbb {E}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t | {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t)&= \int {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t \; \mathcal {P}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t | {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t) \; d{\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t \nonumber \\&= \int {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t \sum _{i=1}^K h_i({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t)\; \mathcal {N}\left( {\varvec{\hat{\mu }}}^{\scriptscriptstyle {{\mathcal {O}}}}_i({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t),{\varvec{\hat{{\varSigma }}}}^{\scriptscriptstyle {{\mathcal {O}}}}_i\right) \nonumber \\&= \sum _{i=1}^K h_i({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t) \; {\varvec{\hat{\mu }}}_i^{\scriptscriptstyle {{\mathcal {O}}}}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t). \end{aligned}$$

(71)

In order to evaluate the covariance, we calculate

$$\begin{aligned} {{\mathrm {cov}}}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t | {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t) = \mathbb {E}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t {{\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t}^{\scriptscriptstyle \top }| {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t ) - \mathbb {E}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t | {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t) \mathbb {E}({{\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t}^{\scriptscriptstyle \top }| {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t). \end{aligned}$$

(72)

We have that

$$\begin{aligned} \mathbb {E}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t {{\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t}^{\scriptscriptstyle \top }| {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t)&= \int {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t {{\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t}^{\scriptscriptstyle \top }\mathcal {P}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t | {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t) \; d{\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t \nonumber \\&= \int \sum _{i=1}^K h_i({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t) \; {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t {{\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t}^{\scriptscriptstyle \top }\mathcal {N}\left( {\varvec{\hat{\mu }}}^{\scriptscriptstyle {{\mathcal {O}}}}_i({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t),{\varvec{\hat{{\varSigma }}}}^{\scriptscriptstyle {{\mathcal {O}}}}_i\right) d{\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t \nonumber \\&= \sum _{i=1}^K h_i({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t) \int {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t {{\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t}^{\scriptscriptstyle \top }\mathcal {N}\left( {\varvec{\hat{\mu }}}^{\scriptscriptstyle {{\mathcal {O}}}}_i({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t),{\varvec{\hat{{\varSigma }}}}^{\scriptscriptstyle {{\mathcal {O}}}}_i\right) d{\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t. \end{aligned}$$

(73)

By using Eq. (72) with a Gaussian distribution, we obtain

$$\begin{aligned} \mathbb {E}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t {{\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t}^{\scriptscriptstyle \top }| {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t)= & {} \sum _{i=1}^K h_i({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t) {\varvec{\hat{{\varSigma }}}}_i^{\scriptscriptstyle {{\mathcal {O}}}}\nonumber \\&+ \sum _{i=1}^K h_i({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t) \; {\varvec{\hat{\mu }}}_i^{\scriptscriptstyle {{\mathcal {O}}}}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t) {\big ({\varvec{\hat{\mu }}}_i^{\scriptscriptstyle {{\mathcal {O}}}}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t)\big )}^{\scriptscriptstyle \top }. \end{aligned}$$

(74)

Combining (72) with (74) we finally have that (see also [90])

$$\begin{aligned} {\varvec{\hat{\varSigma }}}{}^{\scriptscriptstyle {{\mathcal {O}}}}_t= & {} {{\mathrm {cov}}}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {O}}}}_t | {\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t) = \sum _{i=1}^K h_i({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t) \Big ( {\varvec{\hat{{\varSigma }}}}_i^{\scriptscriptstyle {{\mathcal {O}}}} + {\varvec{\hat{\mu }}}_i^{\scriptscriptstyle {{\mathcal {O}}}}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t) \; {{\varvec{\hat{\mu }}}_i^{\scriptscriptstyle {{\mathcal {O}}}}({\varvec{\xi }}^{\scriptscriptstyle {{\mathcal {I}}}}_t)}^{\scriptscriptstyle \top }\Big )\nonumber \\&-\,{\varvec{\hat{\mu }}}^{\scriptscriptstyle {{\mathcal {O}}}}_t {{\varvec{\hat{\mu }}}^{\scriptscriptstyle {{\mathcal {O}}}}_t}^{\scriptscriptstyle \top }. \end{aligned}$$

(75)

Appendix 4: Expectation-maximization for parametric GMM parameters estimation

The following two steps are repeated until convergence, see [102] for details:

E-step:

$$\begin{aligned} h_{t,i} = \frac{\pi _i \mathcal {N}\big ({\varvec{\xi }}_t \big |\; {\varvec{\mu }}_{t,i},{\varvec{{\varSigma }}}_i \big )}{\sum _{k=1}^K \pi _k \mathcal {N}\big ({\varvec{\xi }}_t \big |\; {\varvec{\mu }}_{t,k},{\varvec{{\varSigma }}}_k \big )}. \end{aligned}$$

(76)

M-step:

$$\begin{aligned}&\pi _i \leftarrow \frac{\sum _{t=1}^N h_{t,i}}{N}, \end{aligned}$$

(77)

$$\begin{aligned}&{\varvec{Z}}_i \leftarrow \left( \sum _{t=1}^N h_{t,i} \; {\varvec{\xi }}_t \big [{\varvec{Q}}_t^{\scriptscriptstyle \top },\; 1\big ] \right) \left( \sum _{t=1}^N h_{t,i} {\big [{\varvec{Q}}_t^{\scriptscriptstyle \top },\; 1\big ]}^{\scriptscriptstyle \top }\big [{\varvec{Q}}_t^{\scriptscriptstyle \top },\; 1\big ] \right) ^{-1} ,\end{aligned}$$

(78)

$$\begin{aligned}&{\varvec{{\varSigma }}}_i \leftarrow \frac{\sum _{t=1}^N h_{t,i} \; ({\varvec{\xi }}_t-{\varvec{\mu }}_{t,i})({\varvec{\xi }}_t-{\varvec{\mu }}_{t,i} )^{\scriptscriptstyle \top }}{\sum _{t=1}^N h_{t,i}},\end{aligned}$$

(79)

$$\begin{aligned}&{{\mathrm {where}}}\quad {\varvec{\mu }}_{t,i}={\varvec{Z}}_i \; {\big [{\varvec{Q}}_t^{\scriptscriptstyle \top },\; 1\big ]}^{\scriptscriptstyle \top }. \end{aligned}$$

(80)