# MML Probabilistic Principal Component Analysis

Enes Makalic
Melbourne School of Population and Global Health
University of Melbourne
Carlton, VIC 3053
<EMAIL_ADDRESS>
& Daniel F. Schmidt
Faculty of Information Technology
Monash University
Clayton, VIC 3168
<EMAIL_ADDRESS>

###### Abstract

Principal component analysis (PCA) is perhaps the most widely method for data
dimensionality reduction. A key question in PCA decomposition of data is
deciding how many factors to retain. This manuscript describes a new approach
to automatically selecting the number of principal compoents based on the
Bayesian minimum message length method of inductive inference. We also derive
a new estimate of the isotropic residual variance and demonstrate, via
numerical experiments, that it improves on the usual maximum likelihood
approach.

_Keywords_ minimum message length $\cdot$ factor analysis $\cdot$ model
selection

##  1 Introduction

The principal component analysis (PCA) model postulates that $N$ independent
realisations of $K$-dimensional data ${\bf x}_{i}\in\mathbb{R}^{K}$
($i=1,\ldots,N$) are described as

${\bf x}_{i}=v_{i1}{\bf a}_{1}+\cdots+v_{iJ}{\bf
a}_{J}+\bm{\epsilon}_{i}=\left(\sum_{j=1}^{J}v_{ij}{\bf
a}_{j}\right)+\bm{\epsilon}_{i},\quad\bm{\epsilon}\sim{\rm N}_{K}({\bf
0}_{K},\sigma^{2}{\bf I}_{K}),$ (1)

where $\\{{\bf a}_{1},\ldots,{\bf a}_{J}\\}$ are the $J(<K)$ latent
(unobserved) factor loadings with each factor loading ${\bf
a}_{j}\in\mathbb{R}^{K}$, and ${\bf v}_{j}\sim{\rm N}({\bf 0}_{N},{\bf
I}_{N})$ are the factor scores distributed as per the standard normal
distribution. It is assumed that the residuals follow an isotropic zero mean
normal distribution with the variance-covariance matrix $\sigma^{2}{\bf
I}_{K}$. Without loss of generality, we assume that the data has been centered
to the mean. We can write this PCA model in matrix notation as follows

${\bf x}_{i}={\bf A}{\bf v}_{i}+\bm{\varepsilon},\quad{\bf
A}\in\mathbb{R}^{K\times J},\quad{\bf
v}_{i}\in\mathbb{R}^{J},\quad\bm{\varepsilon}\sim N_{K}({\bf 0},\sigma^{2}{\bf
I}_{K}),$ (2)

where $i=(1,\ldots,N)$, ${\bf A}=({\bf a}_{1},\ldots,{\bf a}_{J})$ and ${\bf
V}=({\bf v}_{1},\ldots,{\bf v}_{N})\in\mathbb{R}^{J\times N}$. Integrating out
the factor scores yields the multivariate Gaussian marginal distribution of
the data

${\bf x}_{i}\sim N_{K}({\bf 0}_{K},\bm{\Sigma}_{0}),\quad\bm{\Sigma}={\bf
A}{\bf A}^{\prime}+\sigma^{2}{\bf I}_{K}.$ (3)

The probabilistic principal component model as described suffers from
identifiability constraints (Anderson and Rubin, 1956; Lawley and Maxwell,
1971). A key reason for this is that the latent factors affect the likelihood
function only through their outer product ${\bf A}{\bf A}^{\prime}$, which
implies that an estimate of the factors can only be determined up to a
rotation. To ensure the model is not over parametarised, the maximum number of
latent factors to be estimated cannot exceed

$J_{\text{MAX}}\leq K+\frac{1}{2}\left(1-\sqrt{8k+1}\right),$ (4)

see (Beal, 2003) (pp. 108) for details. For example, when $K=4,5,6$ we have
$J_{\text{MAX}}=1,2,3$, respectively. Tipping and Bishop (1999) showed how to
interpret standard PCA model in a probabilistic framework and obtained maximum
likelihood estimates of the latent factors and residual variance.

This manuscript examines estimation of the probabilistic PCA model under the
Bayesian minimum message length (MML) inductive inference framework. Although
single and multiple factor analysis has been examined within the MML framework
by Wallace and Freeman (1992) and Wallace (1998) respectively, this manuscript
departs from the earlier work in the following:

  * •

We consider the marginal distribution of the data (3) rather than the model
(1) analysed by Wallace (1998).

  * •

Using polar decomposition, we write the factor load matrix ${\bf A}$ as a
product of an orthogonal matrix and a diagonal matrix representing the
direction and length of the loadings, respectively. Unlike the earlier MML
approaches, we parameterize the orthogonal matrix via Givens rotations to
explicitly capture orthogonality constraints.

  * •

We use matrix polar decomposition to develop a prior distribution for the
latent factors ${\bf A}$ that is a product of a matrix variate Cauchy
distribution and a uniform distribution over the corresponding Stiefel
manifold.

  * •

We obtain analytic MML estimates of the parameters and find a polynomial whose
roots yield the MML estimate of the residual variance.

##  2 Maximum likelihood estimation

This section summarises the results of Tipping and Bishop (1999). The negative
log-likelihood of the data under the probabilistic PCA model is

$\ell(\bm{\theta})=\frac{NK}{2}\log(2\pi)+\frac{N}{2}\log|\bm{\Sigma}|+\frac{N}{2}{\rm
tr}\left(\bm{\Sigma}^{-1}{\bf S}_{x}\right)$ (5)

where ${\bf S}_{x}=\frac{1}{N}\sum_{i}{\bf x}_{i}{\bf x}_{i}^{\prime}$ is the
sample variance-covariance matrix. We have the observed data ${\bf X}$ and
wish to estimate the number of latent factors $J$ and all parameters
$\bm{\theta}=\left\\{{\bf A},\sigma^{2}\right\\}$. Differentiating the
negative log-likelihood with respect to the factor loads

$\displaystyle\partial\ell(\bm{\theta})$ $\displaystyle=N{\rm tr}\,{\bf
A}^{\prime}\bm{\Sigma}^{-1}(\partial{\bf A})-N{\rm tr}\left({\bf
A}^{\prime}\bm{\Sigma}^{-1}{\bf S}_{x}\bm{\Sigma}^{-1}(\partial{\bf
A})\right)$

and setting the derivatives to zero we get

$\displaystyle{\bf S}_{x}{\bf\Sigma}^{-1}{\bf A}$ $\displaystyle={\bf A}$

Consider the singular value decomposition ${\bf A}={\bf U}{\bf L}{\bf
V}^{\prime}$, where ${\bf U}\in\mathbb{R}^{K\times J}$, ${\bf L}={\rm
diag}(l_{1},\ldots,l_{j})$ and ${\bf V}\in\mathbb{R}^{J\times J}$ is an
orthogonal matrix. Noting that $\bm{\Sigma}^{-1}{\bf A}={\bf U}{\bf L}({\bf
L}^{2}+\sigma^{2}{\bf I}_{J})^{-1}{\bf V}^{\prime}$, we have

$\displaystyle{\bf S}_{x}{\bf U}$ $\displaystyle={\bf U}({\bf
L}^{2}+\sigma^{2}{\bf I}_{J})$ $\displaystyle{\bf S}_{x}{\bf u}_{j}$
$\displaystyle=(l_{j}^{2}+\sigma^{2}){\bf u}_{j},\quad(j=1,\ldots,J),$

which is an example of the eigenvalue problem. That is, ${\bf U}$ is a
$(K\times J)$ matrix whose columns are the top $J$ eigenvectors of the sample
covariance matrix ${\bf S}_{x}$ corresponding to the $J$ largest eigenvalues

$\displaystyle\delta_{j}=\hat{l}^{2}_{j}+\sigma^{2},\quad j=1,\ldots,J,$ (6)

where $\hat{l}_{j}=(\delta_{j}-\sigma^{2})^{\frac{1}{2}}$ is the $j$-th
largest singular value of ${\bf A}$. Without loss of generality we assume that
$\delta_{1}>\delta_{2}>\ldots>\delta_{K}>0$ throughout the manuscript. This
implies that the maximum likelihood estimate is

$\displaystyle\hat{\bf A}_{\text{ML}}={\bf U}(\bm{\Delta}-\sigma^{2}{\bf
I}_{J})^{\frac{1}{2}}{\bf O},\quad\bm{\Delta}={\rm
diag}(\delta_{1},\ldots,\delta_{J})$ (7)

where ${\bf O}$ is an arbitrary (orthogonal) rotation matrix and $\bm{\Delta}$
is a diagonal matrix with the $J$-th largest eigenvalues of ${\bf S}_{x}$.
Substituting the maximum likelihood estimate of the factor loads into the
negative log-likelihood we have

$\displaystyle\ell(\sigma,\hat{\bf
A}_{\text{ML}})=\frac{NK}{2}\log(2\pi)+\frac{N}{2}\sum_{i=1}^{J}\log\delta_{j}+\frac{N(K-J)}{2}\log\sigma^{2}+\frac{NJ}{2}+\frac{N}{2\sigma^{2}}\sum_{j=J+1}^{K}\delta_{j}.$
(8)

The concentrated negative log-likelihood is minimised by

$\displaystyle\hat{\sigma}^{2}=\frac{1}{K-J}\sum_{j=J+1}^{K}\delta_{j}$ (9)

which is the sum of the (K-J) smallest eigenvalues of the sample variance-
covariance matrix. Tipping and Bishop (1999) show that these estimates
minimise the negative log-likelihood and discuss other saddle points of the
log-likelihood function.

##  3 Minimum message length analysis of the PCA model

Under suitable regularity conditions (Wallace, 2005)) (pp. 226), the MML87
codelength approximation for data ${\bf x}$ is

$\mathcal{I}_{87}({\bf
x},\bm{\theta})=\underbrace{-\log\pi(\bm{\theta})+\frac{1}{2}\log\lvert{\bf
J}_{\bm{\theta}}(\bm{\theta})\rvert+\frac{P}{2}\log\kappa_{P}}_{\rm
assertion}+\underbrace{\frac{P}{2}-\log p({\bf x}|\bm{\theta})}_{\rm detail}$
(10)

where $\pi_{\bm{\theta}}(\bm{\theta})$ is the prior distribution for the
parameters $\bm{\theta}$, $\lvert{\bf J}_{\bm{\theta}}(\bm{\theta})\rvert$ is
the determinant of the expected Fisher information matrix, $p({\bf
x}|\bm{\theta})$ is the likelihood function of the model and $\kappa_{P}$ is a
quantization constant (Conway and Sloane, 1998; Agrell and Eriksson, 1998);
for small $P$ we have

$\kappa_{1}=\frac{1}{12},\quad\kappa_{2}=\frac{5}{36\sqrt{3}},\quad\kappa_{3}=\frac{19}{192\times
2^{1/3}},$ (11)

while, for large $P$, $\kappa_{P}$ is well-approximated by (Wallace, 2005):

$\frac{p}{2}(\log\kappa_{p}+1)\approx-\frac{p}{2}\log 2\pi+\frac{1}{2}\log
p\pi-\gamma,$ (12)

where $\gamma\approx 0.5772$ is the Euler–Mascheroni constant.

###  3.1 Orthogonality constraints

As discussed in Section 1, it is well-known that the PCA model is not
identifiable given the data. A key reason for this is that the latent vectors
affect the likelihood only through their outer product ${\bf A}{\bf
A}^{\prime}=\sum_{j=1}^{J}{\bf a}_{j}{\bf a}^{\prime}_{j}.$ However, there are
infinitely many sets of vectors that could generate the same matrix. To
resolve this ambiguity, it is a convention to estimate the factor load vectors
to be mutually orthogonal; that is,

${\bf A}^{\prime}{\bf A}=\bm{\alpha}={\rm
diag}(\alpha^{2}_{1},\ldots,\alpha^{2}_{J}),\quad\alpha_{j}=({\bf
a}_{j}^{\prime}{\bf a}_{j})^{\frac{1}{2}},\quad(j=1,\ldots,J),$ (13)

where $\alpha_{j}$ denote the length of the $j$-th load vector. We enforce
orthogonality constraints by parameterizing the matrix ${\bf A}$ in terms of
Givens rotations Pourzanjani et al. (2021). Specifically, we can write ${\bf
A}$ as

$\displaystyle{\bf A}$ $\displaystyle=\left[R_{12}(\phi_{1,2})\cdots
R_{1,K}(\phi_{1,K})R_{2,3}(\phi_{2,3})\cdots R_{2,K}(\phi_{2,K})\cdots
R_{J,J+1}(\phi_{J,J+1})\cdots R_{J,K}(\phi_{J,K}){\bf
I}_{K,J}\right]\bm{\alpha}$ (14) $\displaystyle={\bf R}\,\bm{\alpha}$ (15)

where ${\bf I}_{K,J}$ is the first $J$ columns of a $K\times K$ identity
matrix and $R_{i,j}(\phi_{i,j})$ is a $(K\times K)$ rotation matrix that is
equal to the identity matrix except for the $(i,i)$ and $(j,j)$ positions
which are replaced by $\cos(\phi_{i,j})$, and the $(i,j)$ and $(j,i)$
positions which are replaced by $-\sin(\phi_{i,j})$ and $\sin(\phi_{i,j})$
respectively. Thus ${\bf R}\in\mathbb{R}^{K\times J}$ and
$\bm{\alpha}\in\mathbb{R}^{J\times J}$ denote the orientations and lengths of
the factor load vectors, respectively. For example, when $K=J=2$, we have

${\bf
A}=\left(\begin{array}[]{cc}\cos\left(\phi_{1,2}\right)&-\sin\left(\phi_{1,2}\right)\\\
\sin\left(\phi_{1,2}\right)&\cos\left(\phi_{1,2}\right)\\\
\end{array}\right)\left(\begin{array}[]{cc}1&0\\\ 0&1\\\
\end{array}\right)\left(\begin{array}[]{c}\alpha_{1}\\\ \alpha_{2}\\\
\end{array}\right).$ (16)

This parametarization explicitly takes into account that the estimated factor
loads are pairwise orthogonal. The model parameters are now

  * •

the lengths of the latent factors
$\bm{\alpha}=(\alpha_{1},\ldots,\alpha_{J})\in\mathbb{R}_{+}^{J}$,

  * •

the orientation of the factor load vectors as captured by the $D=JK-J(J+1)/2$
angles

$\bm{\phi}=(\phi_{1,2},\ldots,\phi_{1,K},\phi_{2,3},\ldots,\phi_{2,K},\ldots,\phi_{J,J+1},\ldots,\phi_{J,K})$

  * •

and the residual variance $\sigma^{2}>0$.

###  3.2 Fisher information

Following lengthy and tedious algebra, the expected Fisher information matrix
is seen to be block diagonal with determinant

$\displaystyle|{\bf J}(\bm{\alpha},\sigma,\bm{\phi})|$
$\displaystyle=N^{P}|{\bf J}(\bm{\alpha},\sigma)|\,|{\bf J}(\bm{\phi})|$ (17)
$\displaystyle|{\bf J}(\bm{\alpha},\sigma)|$
$\displaystyle=\frac{2^{J+1}(K-J)}{{\sigma^{2}}}\prod_{j=1}^{J}\frac{\alpha_{j}^{2}}{\left(\alpha_{j}^{2}+\sigma^{2}\right)^{2}}$
(18) $\displaystyle|{\bf J}(\bm{\phi})|$ $\displaystyle=|J_{{\bf
A}\to\bm{\phi}}|^{2}\left(\prod_{i=1}^{J}\left(\frac{\alpha_{i}^{4}}{\sigma^{2}}\right){}^{K-J}\frac{1}{\left(\alpha_{i}^{2}+\sigma^{2}\right){}^{K-1}}\right)\prod_{j<k}\left(\alpha_{j}^{2}-\alpha_{k}^{2}\right)^{2}$
(19)

where $|J_{{\bf A}\to\bm{\phi}}|$ is the transformation of measure under the
Givens representation (Pourzanjani et al., 2021)

$|J_{{\bf
A}\to\bm{\phi}}|=\prod_{i=1}^{J}\prod_{j=i+1}^{K}(\cos\phi_{i,j})^{j-i-1},$

and $P=(D+J+1)$ is the total number of free parameters. Combining all the
terms we have

$\displaystyle|{\bf J}(\bm{\alpha},\sigma,\bm{\phi})|$
$\displaystyle=N^{P}\frac{\left(2^{J+1}(K-J)\right)}{{\sigma^{2}}}|J_{{\bf
A}\to\bm{\phi}}|^{2}\left(\prod_{i=1}^{J}\frac{\alpha_{i}^{2}}{\left(\alpha_{i}^{2}+\sigma^{2}\right)^{2}}\left(\frac{\alpha_{i}^{4}}{\sigma^{2}}\right){}^{K-J}\frac{1}{\left(\alpha_{i}^{2}+\sigma^{2}\right){}^{K-1}}\right)\prod_{j<k}\left(\alpha_{j}^{2}-\alpha_{k}^{2}\right)^{2}$
$\displaystyle=\frac{N^{P}2^{J+1}(K-J)|J_{{\bf
A}\to\bm{\phi}}|^{2}}{{\sigma^{2(J(K-J)+1)}}}\prod_{i=1}^{J}\frac{\alpha_{i}^{4(K-J)+2}}{\left(\alpha_{i}^{2}+\sigma^{2}\right)^{K+1}}\prod_{j<k}^{J}\left(\alpha_{j}^{2}-\alpha_{k}^{2}\right)^{2}$
(20)

###  3.3 Prior information

The prior distribution for the residual standard deviation
$\bm{\Sigma}^{\frac{1}{2}}$ is chosen to be the scale-invariant density

$\pi_{\sigma}(\sigma)\propto\sigma^{-1},$ (21)

defined over some suitable range.

The prior distribution for the matrix of factor loads ${\bf
A}\in\mathbb{R}^{K\times J}$ is not immediately obvious as the estimates of
the factor loads are enforced to be mutually orthogonal. Ideally, we would
like a prior distribution that is uniform over the direction of the $J$
factors, while the distribution of the lengths of these vectors should be
heavy tailed to allow for a wide range of lengths. We follow a similar
approach to Wallace (1998) and assume a prior distribution over the unknown
latent vectors that is then transformed to account for the estimated factors
being mutally orthogonal. Further, as in Wallace (1998), we shall consider a
prior distribution for the scaled factors

$\displaystyle{\bf b}_{j}=\left(\frac{{\bf
a}_{j}}{\sigma}\right),\quad\beta_{j}=({\bf b}_{j}^{\prime}{\bf
b}_{j})^{\frac{1}{2}},\quad(j=1,\ldots,J)$

where the residual variance is used as a default scale. Let $\tilde{\bf
B}\in\mathbb{R}^{K\times J}$ denote the matrix containing the $J$ true
(unknown) scaled factors. We assume $\tilde{\bf B}$ to follow a matrix variate
Cauchy distribution (Bandekar and Nagar, 2003) with probability density
function

$\pi_{\tilde{A}}(\tilde{\bf
B})=\frac{\Gamma_{K}((K+J)/2)}{\pi^{KJ/2}\Gamma_{K}(K/2)}{\rm det}({\bf
I}_{K}+\tilde{\bf B}\tilde{\bf B}^{\prime})^{-(K+J)/2}.$ (22)

This is a reasonable choice as the matrix variate Cauchy is spherically
symmetric and has appropriately heavy tails. Further, our choice of the prior
distribution implies that $\tilde{\bf B}^{\prime}\in\mathbb{R}^{J\times K}$
follows a matrix variate Cauchy distribution with density

$\pi_{\tilde{B}^{\prime}}(\tilde{\bf
B}^{\prime})=\frac{\Gamma_{J}((K+J)/2)}{\pi^{KJ/2}\Gamma_{J}(J/2)}{\rm
det}({\bf I}_{J}+\tilde{\bf B}^{\prime}\tilde{\bf B})^{-(K+J)/2}.$ (23)

Consider the unique matrix polar decomposition

$\tilde{\bf B}^{\prime}={\bf W}_{B}^{\frac{1}{2}}\,{\bf H}_{B},\quad{\bf
W}_{B}=\tilde{\bf B}^{\prime}\tilde{\bf B},\quad{\bf H}_{B}=(\tilde{\bf
B}^{\prime}\tilde{\bf B})^{-\frac{1}{2}}\tilde{\bf B}^{\prime}$ (24)

where ${\bf H}_{B}$ is defined over the Stiefel manifold
$\mathcal{V}_{J}(\mathbb{R}^{K})$ and ${\bf W}_{B}$ is a symmetric positive
definite matrix. We may think of the matrix ${\bf H}_{B}$ as the orientation
matrix, while the matrix ${\bf W}_{B}$ determines the squared lengths of the
true scaled latent vectors. If $\tilde{\bf B}^{\prime}$ follows a matrix
variate Cauchy distribution, it is known that ${\bf H}_{B}$ is distributed
uniformly over the Stiefel manifold with density function (Bandekar and Nagar,
2003):

$\pi_{H}({\bf H}_{B})=\frac{1}{{\rm
Vol}(\mathcal{V}_{J}(\mathbb{R}^{K}))},\quad{\rm
Vol}(\mathcal{V}_{J}(\mathbb{R}^{K}))=\frac{2^{J}\pi^{KJ/2}}{\Gamma_{J}(K/2)},$
(25)

where $\Gamma_{p}(y)$ is the multivariate Gamma function

$\displaystyle\Gamma_{J}(y)=\pi^{J(J-1)/4}\prod_{j=1}^{J}\Gamma(y+(1-j)/2).$

Further, the random variable ${\bf W}_{B}$ representing the squared lengths of
the true scaled factors is independent of ${\bf H}_{B}$ with probability
density function Bandekar and Nagar (2003)

$\displaystyle\pi_{W}({\bf W}_{B})$ $\displaystyle\propto$ $\displaystyle{\rm
det}({\bf W}_{B})^{(K-J-1)/2}{\rm det}({\bf I}_{K}+{\bf W}_{B})^{-(K+J)/2}$
(26)

which is a matrix variate beta type II distribution ${\bf W}_{B}\sim
B_{J}^{II}(K/2,J/2)$ with parameters $(K/2,J/2)$ (see Gupta and Nagar (1999),
pp. 166, for further details); this is also known as the matrix variate $F$
distribution. Recall that the estimated (scaled) factor load vectors obey

${\bf S}_{B}=\tilde{\bf B}\tilde{\bf
B}^{\prime}=\sum_{j=1}^{J}\tilde{\bm{\beta}}_{j}\tilde{\bm{\beta}}_{j}^{\prime}=\sum_{j=1}^{J}\bm{\beta}_{j}\bm{\beta}_{j}^{\prime}={\bf
B}{\bf B}^{\prime},\quad\bm{\beta}_{j}^{\prime}\bm{\beta}_{k\neq j}=0.$ (27)

where ${\bf S}_{B}$ is a $(K\times K)$ symmetric matrix of rank $J$. The
distribution of the squared scaled lengths $\beta_{j}^{2}$ of the estimated
latent vectors can then be taken as the joint distribution of the $J$
eigenvalues of ${\bf S}_{B}$ which is (see Appendix F)

$\displaystyle\pi_{\bm{\beta}^{2}}(\beta_{1}^{2},\ldots,\alpha_{J}^{2})$
$\displaystyle=\frac{\pi^{J^{2}/2}}{\Gamma_{J}(J/2)\mathcal{B}_{J}(K/2,J/2)}\prod_{j=1}^{J}\beta_{j}^{(K-J-1)}(1+\beta_{j}^{2})^{-(K+J)/2}\prod_{j<k}^{J}|\beta_{j}^{2}-\beta_{k}^{2}|,$

where $\mathcal{B}_{p}(a,b)$ denote the multivariate beta function

$\displaystyle\mathcal{B}_{J}(a,b)=\frac{\Gamma_{J}(a)\Gamma_{J}(b)}{\Gamma_{J}(a+b)}.$

This implies that the prior distribution of the lengths of the scaled latent
factors is

$\displaystyle\pi_{\bm{\beta}}(\beta_{1},\ldots,\beta_{J})$
$\displaystyle=\frac{2^{J}\pi^{J^{2}/2}}{\Gamma_{J}(J/2)\mathcal{B}_{J}(K/2,J/2)}\prod_{j=1}^{J}\beta_{j}^{(K-J)}(1+\beta_{j}^{2})^{-(K+J)/2}\prod_{j<k}^{J}|\beta_{j}^{2}-\beta_{k}^{2}|.$
(28)

Finally, the prior distribution for the lengths of the (unscaled) latent
factors is

$\displaystyle\pi_{\bm{\alpha}}(\alpha_{1},\ldots,\alpha_{J})$
$\displaystyle=\frac{2^{J}\pi^{J^{2}/2}\sigma^{J^{2}}}{\Gamma_{J}(J/2)\mathcal{B}_{J}(K/2,J/2)}\prod_{j=1}^{J}\alpha_{j}^{(K-J)}(\sigma^{2}+\alpha_{j}^{2})^{-(K+J)/2}\prod_{j<k}^{J}|\alpha_{j}^{2}-\alpha_{k}^{2}|.$
(29)

The complete prior distribution over all model parameters is

$\displaystyle\pi(\bm{\alpha},\sigma,\bm{\phi})=\pi_{\sigma}(\bm{\Sigma}^{\frac{1}{2}})\,\pi_{\bm{\alpha}}(\alpha_{1},\ldots,\alpha_{J})|J_{{\bf
A}\to\bm{\phi}}|\,J!,$ (30)

where the term $J!$ is included because the labelling of the latent factors is
arbitrary and $|J_{{\bf A}\to\bm{\phi}}|$ is the transformation of measure
from the matrix parametrization ${\bf A}$ to the orthogonality-preserving
parameterization based on Givens rotations.

###  3.4 Codelength

Omitting constants, the Wallace and Freeman (1987) codelength for the
probabilistic PCA model is

$\displaystyle\mathcal{I}$
$\displaystyle\propto\frac{N}{2}\log|\bm{\Sigma}|+\frac{N}{2}{\rm
tr}\left(\bm{\Sigma}^{-1}{\bf
S}_{x}\right)-KJ\log(\sigma)+\frac{1}{2}\sum_{j=1}^{J}\log\left[\alpha_{j}^{2(K-J+1)}\left(\alpha_{j}^{2}+\sigma^{2}\right)^{(J-1)}\right]$

where ${\bf S}_{x}=\frac{1}{N}\sum_{i}{\bf x}_{i}{\bf x}_{i}^{\prime}$ is the
sample variance-covariance matrix. To obtain MML estimates, we start with the
Langrangian of the factor orientations

$\displaystyle\psi({\bf R})$ $\displaystyle=\log|\bm{\Sigma}|+{\rm
tr}\left(\bm{\Sigma}^{-1}{\bf S}_{x}\right)-{\rm tr}{\bf L}({\bf
R}^{\prime}{\bf R}-I),$

where ${\bf L}$ is a $J\times J$ symmetric matrix of Lagrange multipliers.
Clearly, minimising $\psi({\bf R})$ is equivalent to minimising the codelength
with respect to ${\bf R}$. The first differential of the Lagrangian is

$\displaystyle\partial\psi({\bf R})$
$\displaystyle=2\text{tr}\left[\bm{\alpha}{\bf
A}^{\prime}\left(\bm{\Sigma}^{-1}-\bm{\Sigma}^{-1}{\bf
S}_{x}\bm{\Sigma}^{-1}\right)(d{\bf R})\right]-2\text{tr}\left({\bf L}{\bf
R}^{\prime}(d{\bf R})\right),$

which implies the following first order conditions

$\displaystyle\bm{\alpha}{\bf
A}^{\prime}\left(\bm{\Sigma}^{-1}-\bm{\Sigma}^{-1}{\bf
S}_{x}\bm{\Sigma}^{-1}\right)$ $\displaystyle={\bf 0}$ (31) $\displaystyle{\bf
L}{\bf R}^{\prime}$ $\displaystyle={\bf 0}$ (32) $\displaystyle{\bf
R}^{\prime}{\bf R}$ $\displaystyle={\bf I}_{J}$ (33)

From (32) we have that ${\bf L}={\bf 0}$ and from (31)

$\displaystyle{\bf S}_{x}{\bf R}$ $\displaystyle={\bf
R}\,\text{diag}\left(\sigma^{2}+\alpha_{1}^{2},\ldots,\sigma^{2}+\alpha_{J}^{2}\right)$
$\displaystyle{\bf S}_{x}{\bf r}_{j}$ $\displaystyle={\bf
r}_{j}(\sigma^{2}+\alpha_{j}^{2}),\quad(j=1,\ldots,J).$

We see that, at the codelength minimum, the MML estimate of the factor
orientations is the matrix ${\bf R}$ whose columns are the top $J$
eigenvectors of the variance–covariance matrix ${\bf S}_{x}$ with eigenvalues
$\delta_{j}=(\sigma^{2}+\alpha_{j}^{2})$, for $j=1,\ldots J$. This is
identical to the corresponding maximum likelihood estimate. The concentrated
codelength, as a function of $\sigma^{2}$ is

$\displaystyle\mathcal{I}(\sigma)$
$\displaystyle\propto\frac{N}{2}\log\left((\sigma^{2})^{K-J}\prod_{j=1}^{J}(\alpha_{j}^{2}+\sigma^{2})\right)+\frac{N}{2\sigma^{2}}\left(\sum_{j=1}^{K}\delta_{j}\right)-\frac{N}{2\sigma^{2}}\sum_{j=1}^{J}\alpha_{j}^{2}$
$\displaystyle\quad-
KJ\log(\sigma)+\frac{1}{2}\sum_{j=1}^{J}\log\left[\alpha_{j}^{2(K-J+1)}\left(\alpha_{j}^{2}+\sigma^{2}\right)^{(J-1)}\right]$
$\displaystyle=\frac{N(K-J)-KJ}{2}\log\left(\sigma^{2}\right)+\frac{N}{2\sigma^{2}}\left(\sum_{j=1}^{K}\delta_{j}\right)-\frac{N}{2\sigma^{2}}\sum_{j=1}^{J}(\delta_{j}-\sigma^{2})+\frac{(K-J+1)}{2}\sum_{j=1}^{J}\log\left(\delta_{j}-\sigma^{2}\right)$
(34)

The following theorem shows how to obtain the MML estimate of the residual
variance from the concentrated message length.

######  Theorem 1.

Let $\tau=\sigma^{2}$. The concentrated codelength (34) has $(J+1)$ stationary
points whose location are the roots of the $n=(J+1)$-degree polynomial

$a_{n}\tau^{n}+a_{n-1}\tau^{n-1}+\cdots+a_{1}\tau+a_{0},\quad(0<\tau<\delta_{J})$

with coefficients

$\displaystyle a_{0}$ $\displaystyle=-\hat{\tau}_{{\rm ML}}e_{J},$
$\displaystyle a_{j}$ $\displaystyle=(-1)^{j+1}\left[\hat{\tau}_{{\rm
ML}}\,e_{J-j}+\left(1-\frac{(j-1)(J-1)}{N(K-J)}-\frac{K(J-j+1)}{N(K-J)}\right)e_{J-j+1}\right],\quad(1\leq
j\leq J)$ $\displaystyle a_{n}$
$\displaystyle=(-1)^{J}\left[1-\frac{J(J-1)}{N(K-J)}\right]$

where the terms $e_{t}$ refer to elementary symmetric polynomials
$e_{t}(\delta_{1},\ldots,\delta_{J})$ in $J$ variables
$(\delta_{1},\ldots,\delta_{J})$. For example, for $J=3$, we have the
following four elementary symmetric polynomials

$\displaystyle e_{0}(\delta_{1},\delta_{2},\delta_{3})$ $\displaystyle=1$
$\displaystyle e_{1}(\delta_{1},\delta_{2},\delta_{3})$
$\displaystyle=\delta_{1}+\delta_{2}+\delta_{3}$ $\displaystyle
e_{2}(\delta_{1},\delta_{2},\delta_{3})$
$\displaystyle=\delta_{1}\delta_{2}+\delta_{1}\delta_{3}+\delta_{2}\delta_{3}$
$\displaystyle e_{3}(\delta_{1},\delta_{2},\delta_{3})$
$\displaystyle=\delta_{1}\delta_{2}\delta_{3}.$

MML estimate of the residual variance is the stationary point in the domain
$0<\tau<\delta_{J}$ that yields the shortest codelength. MML estimates of the
factor lengths can be obtained from
$\hat{\alpha}_{j}=(\delta_{j}-\hat{\sigma}^{2})^{\frac{1}{2}}$ for all
$j=1,\ldots,J$.

Importantly, all $(J+1)$ elementary symmetric polynomials can be efficiently
computed in $O(J\log^{2}J)$ time by Fast Fourier Transform (FFT) polynomial
multiplication, following an algorithm widely attributed to Ben-Or.

Example 1: For a single latent factor ($J=1$), the stationary points of the
concentrated codelength are the roots of the quadratic polynomial in $\tau$:

$\displaystyle-\delta_{1}\hat{\tau}_{\text{ML}}+\left(\hat{\tau}_{\text{ML}}+c\delta_{1}\right)\tau-\tau^{2}=0,\quad
c=1-\frac{K}{N(K-1)},$ (35)

given by

$\displaystyle\frac{1}{2}\left(\hat{\tau}_{\text{ML}}+c\,\delta_{1}\pm\Delta^{\frac{1}{2}}\right),\quad\Delta=c^{2}\delta_{1}^{2}+2(c-2)\delta_{1}\hat{\tau}_{\text{ML}}+\hat{\tau}_{\text{ML}}^{2}.$

The quadratic polynomial has no real roots if:

$\displaystyle-\frac{c+2\sqrt{1-c}-2}{c^{2}}<\frac{\delta_{1}}{\hat{\tau}_{\text{ML}}}<\frac{-c+2\sqrt{1-c}+2}{c^{2}}$
(36)

in which case the MML solution is the uncorrelated multivariate Gaussian
model. For example, when $N=25$ and $K=4$, the quadratic will have no real
roots if

$\displaystyle
0.219<\frac{\delta_{1}}{\delta_{2}+\delta_{3}+\delta_{4}}<0.564.$

In the limit as $N\to\infty$, we have

$\displaystyle\lim_{N\to\infty}\pm\frac{c+2\sqrt{1-c}-2}{c^{2}}=1.$

so that both the lower and upper bound in (36) approach 1.

Example 2: For a PCA model with two latent factors ($J=2$), the stationary
points of the concentrated codelength are the roots of the cubic polynomial in
$\tau$:

$\displaystyle-\delta_{1}\delta_{2}\hat{\tau}_{\text{ML}}+\left((\delta_{1}+\delta_{2})\hat{\tau}_{\text{ML}}+c_{1}\delta_{1}\delta_{2}\right)\tau-\left(\hat{\tau}_{\text{ML}}+\left(c_{0}+c_{1}\right)(\delta_{1}+\delta_{2})\right)\tau^{2}+(2c_{0}+c_{1})\tau^{3}$

where the constants

$\displaystyle c_{0}=\frac{K-1}{N(K-2)},\quad c_{1}=1-\frac{2K}{N(K-2)},$

depend only on $N$ and $K$.

##  4 Experiments

###  4.1 Parameter estimation

$N$ | $K$ | $J$ | $S_{1}$ |  | $S_{2}$ |  | KL Divergence
---|---|---|---|---|---|---|---
|  |  | MLE | MML87 |  | MLE | MML87 |  | MLE | MML87
25 | 5 | 1 | -0.072 | -0.023 |  | 0.011 | 0.007 |  | 0.167 | 0.130
2 | -0.093 | 0.001 |  | 0.021 | 0.009 |  | 0.279 | 0.179
8 | 1 | -0.069 | -0.030 |  | 0.008 | 0.004 |  | 0.244 | 0.194
2 | -0.120 | -0.021 |  | 0.018 | 0.005 |  | 0.463 | 0.285
4 | -0.103 | 0.038 |  | 0.020 | 0.007 |  | 0.647 | 0.353
16 | 1 | -0.066 | -0.034 |  | 0.006 | 0.003 |  | 0.407 | 0.332
2 | -0.111 | -0.040 |  | 0.014 | 0.003 |  | 0.818 | 0.557
4 | -0.207 | -0.012 |  | 0.045 | 0.003 |  | 1.809 | 0.764
50 | 5 | 1 | -0.035 | -0.009 |  | 0.004 | 0.003 |  | 0.074 | 0.065
2 | -0.056 | 0.006 |  | 0.008 | 0.004 |  | 0.127 | 0.096
8 | 1 | -0.034 | -0.012 |  | 0.003 | 0.002 |  | 0.113 | 0.100
2 | -0.060 | -0.008 |  | 0.005 | 0.002 |  | 0.208 | 0.159
4 | -0.069 | 0.034 |  | 0.009 | 0.004 |  | 0.324 | 0.209
16 | 1 | -0.033 | -0.015 |  | 0.002 | 0.001 |  | 0.209 | 0.186
2 | -0.056 | -0.018 |  | 0.004 | 0.001 |  | 0.401 | 0.323
4 | -0.107 | -0.010 |  | 0.012 | 0.001 |  | 0.785 | 0.489
100 | 5 | 1 | -0.017 | -0.004 |  | 0.002 | 0.001 |  | 0.033 | 0.031
2 | -0.031 | 0.004 |  | 0.003 | 0.002 |  | 0.058 | 0.049
8 | 1 | -0.016 | -0.005 |  | 0.001 | 0.001 |  | 0.051 | 0.048
2 | -0.029 | -0.002 |  | 0.002 | 0.001 |  | 0.095 | 0.082
4 | -0.049 | 0.021 |  | 0.004 | 0.002 |  | 0.159 | 0.118
16 | 1 | -0.016 | -0.006 |  | 0.001 | 0.000 |  | 0.100 | 0.094
2 | -0.027 | -0.007 |  | 0.001 | 0.000 |  | 0.193 | 0.171
4 | -0.054 | -0.004 |  | 0.003 | 0.000 |  | 0.367 | 0.283
Table 1: Performance metrics for maximum likelihood (MLE) and MML87 estimates
of residual variance $\sigma^{2}$ computed over $10^{5}$ simulations.

This section compares the newly derived MML parameter estimates for the
probabilistic PCA model to the standard approach based on the maximum
likelihood estimator. Since MML and maximum likelihood estimates of the the
factor lengths (for a given $\sigma^{2}$) and factor orientations are
identical, the key difference between to two approaches is in the estimation
of the residual variance. Our simulation experiments are loosely based on
Section 6 in (Wallace, 1998). We conducted $10^{5}$ simulations for each
combination of the sample size $N\in\\{25,50,100\\}$, the dimensionality of
the data $K\in\\{5,8,16\\}$ and the number of estimated latent factors
$J\in\\{1,2,4\\}$. As both maximum likelihood and MML are scale invariant, the
residual variance was set to $\sigma^{2}=1$ and the factors lengths were
$\alpha_{j}=1$ ($j=1,\ldots,J$) for each simulation run, without loss of
generality. The factor directions were randomly sampled from a unit
$K$-sphere.

We used the three performance metrics discussed in (Wallace, 1998) to evaluate
the estimators:

$\displaystyle S_{1}=\log\hat{\sigma}_{i},\quad
S_{2}=(\log\hat{\sigma}_{i})^{2},$

and the Kullback–Leibler (KL) divergence (Kullback and Leibler, 1951) between
two multivariate Gaussian distributions

$\text{KL}(\bm{\Sigma}_{0},\bm{\Sigma}_{1})=\frac{1}{2}\left(\text{tr}\left(\bm{\Sigma}^{-1}_{1}\bm{\Sigma}_{0}\right)+\log\left(\frac{\bm{|\Sigma}_{1}|}{|\bm{\Sigma}_{0}|}\right)-K\right).$

The first metric $S_{1}$ is a measure of bias, while $S_{2}$ measures error in
any direction. Both $S_{1}$ and $S_{2}$ are zero for exact estimates as the
true residual variance was $\sigma^{2}=1$ in all experiments. The error
measures were specifically chosen as they do not depend on the number of
estimated latent vectors $J$. Simulation results averaged over $10^{5}$
iterations are shown in Table 1.

The MML estimate of the residual variance was found to be superior to the
usual maximum likelihood estimate for all tested combinations of sample sizes,
data dimensionality and the number of latent vectors. Maximum likelihood
appeared to underestimate the residual variance more strongly compared to the
minimum message length estimate. The differences in the performances of the
two estimates were most pronounced when the sample size and data
dimensionality was small ($N\leq 50$, $K=5$).

###  4.2 Model selection

We have also compared the performance of MML model selection against the
highly popular Bayesian information criterion (BIC) and Laplace’s method for
approximating the marginal distribution of the data (Minka, 2000), referred to
as ‘Bayes’ henceforth. Using numerical experiments, Minka (2000) demonstrated
that approximating Bayesian evidence is superior to methods like cross
validation. The simulation setup was identical to Section 4.1 except the
sample size was $N\in\\{50,100\\}$, the dimensionality of the data $K=10$ and
the number of estimated latent factors $J\in\\{1,2,4\\}$. Simulation results,
averaged over $10^{5}$ iterations, are shown in Table 2. As expected, both MML
and the Bayes method have similar performance and both improve significantly
over the popular BIC criterion.

$N$ | $J$ | Method | KL Divergence |  | Model Selection (%)
---|---|---|---|---|---
|  |  |  |  | $<J$ | $=J$ | $>J$
50 | 1 | MML | 0.060 |  | – | 99.55 | 0.45
BIC | 0.063 |  | – | 100.00 | 0.00
Bayes | 0.066 |  | – | 97.19 | 2.81
2 | MML | 0.126 |  | 70.40 | 28.32 | 1.27
BIC | 0.145 |  | 96.40 | 3.60 | 0.00
Bayes | 0.129 |  | 52.52 | 45.17 | 2.31
4 | MML | 0.198 |  | 81.00 | 5.78 | 13.22
BIC | 0.257 |  | 99.99 | 0.01 | 0.00
Bayes | 0.209 |  | 91.09 | 7.92 | 0.98
100 | 1 | MML | 0.028 |  | – | 99.74 | 0.26
BIC | 0.029 |  | – | 100.00 | 0.00
Bayes | 0.030 |  | – | 97.86 | 2.14
2 | MML | 0.062 |  | 30.47 | 68.49 | 1.05
BIC | 0.094 |  | 74.75 | 25.25 | 0.00
Bayes | 0.060 |  | 17.96 | 79.84 | 2.20
4 | MML | 0.096 |  | 66.32 | 22.70 | 10.98
BIC | 0.166 |  | 99.48 | 0.52 | 0.00
Bayes | 0.103 |  | 72.71 | 25.93 | 1.36
Table 2: Model selection simulation results for minimum message length (MML),
Bayesian information criterion (BIC) and Laplace’s method for estimating
Bayesian evidence averaged over $10^{5}$ simulations. In all experiments, data
dimensionality was $K=10$.

## Appendix A. Joint eigenvalue distribution for the central $F$ matrix

A $(p\times p)$ random symmetric positive definite matrix has a matrix variate
beta type II distribution with parameters $(a,b)$ if it has the probability
density function

$\frac{\Gamma_{p}(a+b)}{\Gamma_{p}(a)\Gamma_{p}(b)}{\rm det}({\bf
V})^{a-(p+1)/2}{\rm det}({\bf I}_{p}+{\bf V})^{-(a+b)},{\bf V}>0,$ (37)

where $a>(p-1)/2$ and $b>(p-1)/2$. We write ${\bf V}\sim B_{p}^{II}(a,b)$ to
denote this distribution, which is also known as the matrix variate $F$
distribution. Let $\mathcal{B}_{p}(a,b)$ denote the multivariate beta function

$\mathcal{B}_{p}(a,b)=\frac{\Gamma_{p}(a)\Gamma_{p}(b)}{\Gamma_{p}(a+b)}.$
(38)

Consider the random variable ${\bf V}\sim B_{p}^{II}(n_{1}/2,n_{2}/2)$ and the
transformation ${\bf V}={\bf H}\bm{\Lambda}{\bf H}^{\prime}$ from ${\bf V}$ to
its eigenvalues $\bm{\Lambda}={\rm diag}(\lambda_{1},\ldots,\lambda_{p})$ and
eigenvectors ${\bf H}$, where ${\bf H}\in O(p)$ is in the orthogonal group
with the $j$-th column being the normalized eigenvector of ${\bf V}$
corresponding to the eigenvalue $\lambda_{j}$. The joint distribution of the
$p$ eigenvalues $\bm{\Lambda}$ of ${\bf V}$ is ((Muirhead, 1982), Theorem
3.2.17, pp. 104)

$\displaystyle\pi_{\Lambda}(\lambda_{1},\ldots,\lambda_{p})$ $\displaystyle=$
$\displaystyle\frac{\pi^{p^{2}/2}}{\Gamma_{p}(p/2)}\prod_{i<j}|\lambda_{i}-\lambda_{j}|\int_{O(p)}f({\bf
H}\bm{\Lambda}{\bf H}^{\prime})(d{\bf H}).$ (39)

The integral can be evaluated as follows

$\displaystyle\int_{O(p)}f({\bf H}\bm{\lambda}{\bf H}^{\prime})(d{\bf H})$
$\displaystyle=$
$\displaystyle\frac{1}{\mathcal{B}_{p}(n_{1}/2,n_{2}/2)}\prod_{j=1}^{p}\lambda_{j}^{(n_{1}-p-1)/2}(1+\lambda_{j})^{-(n_{1}+n_{2})/2}\int_{O(p)}(d{\bf
H})$ (40) $\displaystyle=$
$\displaystyle\frac{1}{\mathcal{B}_{p}(n_{1}/2,n_{2}/2)}\prod_{j=1}^{p}\lambda_{j}^{(n_{1}-p-1)/2}(1+\lambda_{j})^{-(n_{1}+n_{2})/2}$

where (see Muirhead (1982), pp. 104)

$\int_{O(p)}(d{\bf H})=1.$ (41)

## References

  * Anderson and Rubin [1956] T. W. Anderson and Herman Rubin.  Statistical inference in factor analysis.  In _Berkeley Symposium on Mathematical Statistics and Probability_ , pages 111–150, 1 1956.
  * Lawley and Maxwell [1971] D. N. Lawley and A. E. Maxwell.  _Factor analysis as a statistical method_.  American Elsevier Publishing, 1971.
  * Beal [2003] M. J Beal.  _Variational Algorithms for Approximation Bayesian Inference_.  PhD thesis, 2003.
  * Tipping and Bishop [1999] M. E. Tipping and C. M. Bishop.  Probabilistic principal component analysis.  _Journal of the Royal Statistical Society (Series B)_ , 21(3):611–622, 1999.
  * Wallace and Freeman [1992] Chris S. Wallace and Peter R. Freeman.  Single-factor analysis by minimum message length estimation.  _Journal of the Royal Statistical Society (Series B)_ , 54(1):195–209, 1992.
  * Wallace [1998] C. S. Wallace.  Intrinsic classification of spatially correlated data.  _The Computer Journal_ , 41(8):602–611, 1998\.
  * Wallace [2005] Chris S. Wallace.  _Statistical and Inductive Inference by Minimum Message Length_.  Information Science and Statistics. Springer, first edition, 2005.
  * Conway and Sloane [1998] J. H. Conway and N. J. A. Sloane.  _Sphere Packing, Lattices and Groups_.  Springer-Verlag, third edition, December 1998.
  * Agrell and Eriksson [1998] E. Agrell and T. Eriksson.  Optimization of lattices for quantization.  _IEEE Transactions on Information Theory_ , 44(5):1814–1828, September 1998.
  * Pourzanjani et al. [2021] Arya A. Pourzanjani, Richard M. Jiang, Brian Mitchell, Paul J. Atzberger, and Linda R. Petzold.  Bayesian inference over the Stiefel manifold via the Givens representation.  _Bayesian Analysis_ , 16(2), 2021.  doi:10.1214/20-BA1202.
  * Bandekar and Nagar [2003] Rajesh R. Bandekar and Daya K. Nagar.  Matrix variate Cauchy distribution.  _Statistics: A Journal of Theoretical and Applied Statistics_ , 37:537–550, 2003.
  * Gupta and Nagar [1999] A. K. Gupta and D. K. Nagar.  _Matrix Variate Distributions_.  Chapman and Hall/CRC, 1999.
  * Wallace and Freeman [1987] Chris S. Wallace and Peter R. Freeman.  Estimation and inference by compact coding.  _Journal of the Royal Statistical Society (Series B)_ , 49(3):240–252, 1987.
  * Kullback and Leibler [1951] S. Kullback and R. A. Leibler.  On information and sufficiency.  _The Annals of Mathematical Statistics_ , 22(1):79–86, March 1951.
  * Minka [2000] Thomas Minka.  Automatic choice of dimensionality for PCA.  In T. Leen, T. Dietterich, and V. Tresp, editors, _Advances in Neural Information Processing Systems_ , volume 13. MIT Press, 2000.  URL https://proceedings.neurips.cc/paper/2000/file/7503cfacd12053d309b6bed5c89de212-Paper.pdf.
  * Muirhead [1982] Robb J. Muirhead.  _Aspects of Multivariate Statistical Theory_.  Wiley Series in Probability and Statistics. Wiley, 1982.