# Robust Bilinear Probabilistic Principal Component Analysis
## 3 Robust BPPCA
### Estimation of the Parameters
The reason is that the calculation of \((\Sigma_{c}^{0})^{-1}\) and \((\Sigma_{r}^{0})^{-1}\) can be more efficiently performed by using \[\operatorname{etr}\Big{(}-\frac{1}{2}(\mu_{r,n}\Sigma_{r})^{-1}(Y ^{r}_{n})^{\text{T}}\mu_{c,n}^{-1}Y^{r}_{n}\Big{)}\Big{(}\frac{v}{2}\Big{)}^{ \frac{v}{2}}\mu_{n}^{\frac{v}{2}-1}\exp\Bigl{(}-\frac{\
u}{2}n_{n}\Big{)}\frac{ 1}{\Gamma(\frac{v}{2})}\Bigg{)}\]0 In its per-iteration of Algorithm 1, the most expensive computational cost is \(\mathcal{O}(Nd_{c}d_{r}(d_{c}+d_{r}))\) appearing on the formation of \(\rho_{n}^{0}\) with \(1\leq n\leq N\). Owing to introducing the new latent variable \(\mu_{n}\), RBPPCA is a little more time complex than BPPCA having a calculation cost of \(\mathcal{O}(Nd_{c}d_{r}(q_{c}+q_{r}))\). However, it will be shown in our numerical example that the RBPPCA algorithm presents less sensitivity to outliers. 2. Compared with the AECM algorithm of BPPCA in [12] which uses the centered data and estimates \(\{\sigma_{c},C\}\) and \(\{\sigma_{r},R\}\) based on the model \(X=CY^{\tau}+Y_{\varepsilon}^{\tau}\) and \(X=Y^{\tau}R^{\text{T}}+Y_{\varepsilon}^{\tau}\), respectively, two more parameters \(\
u\) and \(W\) are needed to be computed in the AECM iteration of RBPPCA. Notice that both the models (14) and (18) contain the parameters \(\
u\) and \(W\). Thus, we split the parameter set \(\theta\) into \(\theta_{c}=\{\
u,\sigma_{c},W,C\}\) and \(\theta_{r}=\{\
u,\sigma_{r},W,R\}\) which naturally leads to the parameters \(\
u\) and \(W\) being calculated twice in each loop of Algorithm 1. Though other partitions of the set \(\theta\), such as \(\{\
u,\sigma_{c},C\}\) and \(\{\sigma_{r},W,R\}\), are also available for the estimation of parameters, we prefer \(\theta_{c}\) and \(\theta_{r}\), because updating \(\
u\) and \(W\) one more time in each iteration can be obtained by adding a little more computational cost. 3. As stated in Section 3.3 of [24], any AECM sequence increases \(L(\theta)=\sum_{n=1}^{N}\ln(p(X_{n}|\theta))\) at each iteration, and converges to a stationary point of \(L(\theta)\). Notice that the convergence results of the AECM algorithm proved in Section 3.3 of [24] do not depend on the distributions of the data sets. Therefore, up to set a limit on the maximum number of steps \(iter_{max}\), we use the following relative change of log-likelihood as the stopping criterion, i.e., \[\operatorname{etr}\Big{(}-\frac{1}{2}(\mu_{r,n}\Sigma_{r})^{-1}(Y ^{r}_{n})^{\text{T}}\mu_{c,n}^{-1}Y^{r}_{n}\Big{)}\Big{(}\frac{v}{2}\Big{)}^{ \frac{v}{2}}\mu_{n}^{\frac{v}{2}-1}\exp\Bigl{(}-\frac{\
u}{2}n_{n}\Big{)}\frac{ 1}{\Gamma(\frac{v}{2})}\Bigg{)}\]1 (29) where \(\varepsilon\) is a specified tolerance used, which by default is set to \(10^{-5}\) in our numerical examples. 4. Based on the computed results of Algorithm 1, and similar to PPCA [8] and BPPCA [12], it is known that \[\operatorname{etr}\Big{(}-\frac{1}{2}(\mu_{r,n}\Sigma_{r})^{-1}(Y ^{r}_{n})^{\text{T}}\mu_{c,n}^{-1}Y^{r}_{n}\Big{)}\Big{(}\frac{v}{2}\Big{)}^{ \frac{v}{2}}\mu_{n}^{\frac{v}{2}-1}\exp\Bigl{(}-\frac{\
u}{2}n_{n}\Big{)}\frac{ 1}{\Gamma(\frac{v}{2})}\Bigg{)}\]2 can be considered as the compressed representation of \(X\). Hence, we can reconstruct \(X\) as \[\operatorname{etr}\Big{(}-\frac{1}{2}(\mu_{r,n}\Sigma_{r})^{-1}(Y ^{r}_{n})^{\text{T}}\mu_{c,n}^{-1}Y^{r}_{n}\Big{)}\Big{(}\frac{v}{2}\Big{)}^{ \frac{v}{2}}\mu_{n}^{\frac{v}{2}-1}\exp\Bigl{(}-\frac{\
u}{2}n_{n}\Big{)}\frac{ 1}{\Gamma(\frac{v}{2})}\Bigg{)}\]3 (30)  
```
0: Initialization \(\
u^{0}\), \(\sigma_{c}^{0}\), \(\sigma_{r}^{0}\), \(W^{0}\), \(C^{0}\), \(R^{0}\), and sample matrices \(\{X_{n}\}_{n=1}^{N}\). Compute \(\Phi_{c}^{0}\). 0: the converged \(\{\
u,\sigma_{c},\sigma_{r},W,C,R\}\) for\(i=0,1,2,\ldots\), until \(\zeta\) defined in (29) is smaller than a threshold do % Stage 1  Compute \(\Phi_{r}^{0}\), \(\rho_{n}^{0}\), \(E_{\mu_{n}}^{0}\) and \(E_{Y_{\varepsilon}^{0}}^{0}\) via (21). Compute \(\widehat{W}^{(0)}\), \(C^{(s-1)}\) and \(\sigma_{c}^{(s-1)}\) by (23) and (24).