A Bayesian Factorised Covariance Model for Image Analysis

1 Introduction

Statistical learning is about inducing trend from data. Given observations of a set of attributes, trend discovery means to learn the relationship between those interested attributes. General relationship is an elusive notion, which may involve attribute groups of any size and can manifest themselves in arbitrary patterns in the observed data. Nevertheless, one can learn rich information about the population of the data from the knowledge about how attribute pairs vary in a coordinated manner. That is, given instances of a set of variables, a useful analysis is to estimate the covariance matrix of the variables.

Learning the covariance is theoretically straightforward but practically difficult. The sample covariance matrix asymptotically approximates the population covariance. But when the number of variables, p, is large compared with the sample size N, the estimation is unreliable [Marcenko and Pastur, 1967].

∗This work was supported by the Australian Research Council under Discovery Project ARC DP-120103730.

This is hardly surprising. A non-degenerated population of p variables has a covariance matrix lying on a Stiefel manifold characterised by p(p-1)/2 parameters. In the context of modern high-dimensional data analysis, even for a relatively moderate p, the N required for the naive covariance estimator to work properly could be prohibitively large. Therefore, given limited observations, it is necessary to exploit side information about the data population to obtain a reliable estimate of covariance.

In this paper, we are concerned with a special setting for the covariance estimating problem, where the variables are organised as a matrix. In particular, we propose a Bayesian model, where the covariance between two variables factorises into two parts: (i) the covariance between the two corresponding rows of the matrix and (ii) that of the two columns. The row- and column-wise covariance matrices are considered as random variables and learned by Bayesian inference. Therefore, we employ the prior knowledge about the matrix structure of the variables as the side information regularising the learning, and use the name of Bayesian factorised covariance model (BFCM). The inference in BFCM is more manageable than directly learning the original covariance matrix, which is of much higher dimension than either of the two factor covariance matrices. The Bayesian treatment, i.e. prior distribution of the factor covariance, provides additional regularisation.

BFCM is particularly suited to images. Generally, direct estimate of pixel-to-pixel covariance is difficult, because practical images often contain a large number of pixels. The hypothesis of BFCM states that the correlation between two pixels is induced by their respective positions in the image, i.e.

The hypothesis is sensible, because a category of images is often characterised by the relations between the rows and those between the columns in those images.

More generally, the proposed model benefits covariance estimation for matrix data in two aspects. First, the factorisation structure significantly reduces the degree of freedom in the covariance estimator and improves its scalability. Second, the factorisation structure makes the covariance estimator more robust. For example, if the information source is faulty, e.g. due to a damaged sensor, some variables can be missing in all observations, it would be infeasible to estimate element-wise covariance, but relations between rows and columns can still be inferred.

2 Bayesian factorised covariance model

This section specifies the definition and inference of BFCM. Frequently used denotations are listed as follows. The symbol ⊗ represents the Kronecker product of two matrices, operator vec (·) transforms a matrix to a column vector by stacking the columns of the matrix, vec (A) = [a11 . . . am1, . . . , a1n . . . amn]T , and function etr[A] := exp{tr(A)}, where tr(·) is matrix trace. The formulations of standard probability density functions follow those in [Gelman et al., 2003].

2.1 Model Definition

BFCM focuses on the relations between the entries in a randommatrix. Fig. 1 shows a graphical illustration of the probabilistic structure of BFCM, and the Bayesian model is defined as follows. The data are N samples of random matrices. We consider an observed matrix Y(n) ∈ Rnr×nc to be consisting of an underlying pattern F(n) ∈ Rnr×nc and independent Gaussian noises

for i = 1 . . . nr and j = 1 . . . nc, where σ2 represents the variance of the noise. Note that the sample index (n) will be omitted when there is no ambiguity. Then the covariance between the elements of F is of interest. According to BFCM, the covariance has a factorised form

According to the discussion in Sec. 1, Ωri and Ωcj denote the indexing sets of the rows and the columns. Note that these indexes only represent conceptual objects as place holders in cov [·, ·] function to give rise to their respective covariance matrices. The objects themselves are not of interest, and will have no significance in the following development of BFCM. Let the row- and column-covariance matrices be R and C, respectively:

Then the factorisation (2) states that the element-wise covariance of F, cov [vec (F) , vec (F)], is the Kronecker product of R and C. Given R and C, the prior of F can be specified by the following Gaussian distribution

An equivalent form of (3) is a matrix variate normal distribution [Gupta and Nagar, 2002],

where det(·) denotes matrix determinant. BFCM explicitly treats R and C as latent random variables, and defines priors and performs inferences on them. Technically, the model is formulated in terms of R-1 and C-1 . The prior on R-1 is defined as

where W(.) and L(.) represent the Wishart and Laplace density functions and ‖·‖1 is the 1-norm of a matrix. With respect to R-1 , the prior has the conjugate form to (4). Furthermore, the Laplace term encourages R-1 to be sparse. The prior belief of sparsity on R-1 is sensible and helping inference. An entry in R-1 corresponds to direct connections between two rows, which is a strong relation: if [R-1 ]i1,i2 ̸= 0, then knowing the i1-th row provides unique information about the i2-th row, which is not provided by knowning all the remaining rows, and vice versa. The sparse prior encourages the model to focus on significant relations, which reduces the risk of over-fitting and improves efficiency of evaluating functions like (4). The distribution of C is defined in the same form as that of R.

We employ the conjugate prior for the noise variance σ-2 , which is a gamma distribution

where G(.) stands for the gamma density function. Therefore the Bayesian model structure in Fig. 1 is completed by the joint density

where we use concise denotations for the priors: p(R), p(C) are formulated by (5) and p(σ2 ) by (6).

2.2 Inference

In this subsection, we are concerned with the posterior distributions of the latent variables given the data. The posterior is determined by the joint probability. Substituting (4), (5), (6)

and (1) in 7, and taking the logarithmic form, we have

The joint posterior w.r.t. all the latent variables is not computationally tractable. Thus the inference is performed in an alternating manner, where we compute the posterior of each variable in turn by conditioning (8) on the remaining variables.

Posterior distributions

Posterior of F(n) In (8), F (omitting the sample index) appears in the product of two Gaussian density functions p(F|R,C)p(Y|F), which is specified by (1) and (3), respectively. Thus the posterior of F given R, C and Y is a Gaussian

Posterior of R and C Since the prior (5) is derived from the conjugate family of (4), the posterior of R has the same functional form as the prior. Considering the relevant terms in (8) gives the posterior as

where ˜S = τ0I + ∑n F(n)C-1F(n)T . The posterior of C can be derived similarly.

Posterior of σ-2 The gamma prior of σ-2 is conjugate to the likelihood (1), thus the posterior is also a gamma distribution:

where ‖ · ‖F denotes the Frobenius norm.

Computation procedure

We compute the maximum a posteriori (MAP) solution to the posterior distributions of the latent variables. We choose MAP not only because there are off-the-shelf optimisation algorithms providing convenient and fast computations, but also because the posterior modes of R-1 and C-1 are sparse, which improves model robustness and interpretability. We implement MAP as an alternating optimisation of (8) in terms of F, σ-2 , R-1 and C-1 . The computation of F and that of R-1 and C-1 need special attention.

Compute posterior mode of F The posterior of F is a standard Gaussian distribution given all the other variables. The corresponding mode is the mean parameter of the density function (9). However, F is a multivariate distribution of nr ×nc variables, and the naive computation of the posterior mode following (10) and (11) involves an expensive matrix inversion. Thus we compute efficiently [Yan et al., 2011] by using SVD of the factor covariance R and C. Let

Then ΣF can be obtained by reformulating (11)

where U = UC ⊗ UR, S = SC ⊗ SR, and the diagonal matrix Λ is given by

Comparing (10), we can compute µF efficiently as

By the relation of Kronecker product and matrix vectorisation, we have

and σ-2Λ × vec (UTRYUC ) can be arranged as a elementwise product as

where Ψ is a nr × nc matrix of the diagonal of Λ and λr and λc are column vectors of the diagonal of SR and SC, respectively. The operator ⊙ represents element-wise product. Employing the same trick as (14) completes the computation of the posterior mean µF

Compute posterior model of R and C The posterior of R is given by (12). We consider the following change of variables in the equation, letting Q = R-1 , ˜S′ = ˜S/(ν0 +

Nnc) and ρ′ = 2ρ/(ν0 +Nnc). Then maximising (12) w.r.t. R-1 is equivalent to

The maximisation can be solved efficiently by employing the graph Lasso algorithm [Friedman et al., 2008].

In theory, the time complexity for each iteration is O(ˆn3 ), where ˆn = max(nr, nc,N). In practice, the computation is dominated by the computations of R-1 and C-1 , for which the worst scenario time complexities are O(n3r ) and O(n3c ) [Friedman et al., 2008].

written digits.

solution to the optimisation problem has been developed in

the factorised covariance structure

There are several efforts addressing the small sample size (SSS) problem related to covariance estimation for highdimensional data. For images, [Shashua and Levin, 2001] pointed that for compression, a model should account for not only redundancies between images, but also the spatial redundancies (row/column-correlations) within individual images, and they represented a set of images as a 3-D array. The model relieves the SSS problem, since it treats rows and columns, rather than entire images as units of analysis - it is possible to compute a valid model from a single image. In [Shashua and Levin, 2001], the computation is based on addition of rank-1 arrays, while in [Vasilescu and Terzopouexperiment, BFCM is used to analyse images of handwritten digits from the USPS dataset [Hull, 1994]. A BFCM is learned for each digit by using 10 randomly chosen images. Fig. 2 shows the row- and column-wise covariance matrices learned by the models. A brief investigation of the entries of the resultant covariance suggests that they are meaningful and representing the characteristic structural information of the corresponding image classes. For example, Fig. 3 illustrates the interpretation of two entries of the discovered covariance: one from the row-covariance for digit "1" and the other from the column-covariance for digit "8". The resultant rowcovariance of "1" suggests that many rows of a "1"-image should be highly correlated, which is consistent with our truth values according to the figure visible artefacts, because the obse insufficient compared to the numb 40 = 2, 000) for learning the proba To address the small sample siz sent the images as a multi-dimensi and Levin, 2001] by using high-or [Vasilescu and Terzopoulos, 2002]. missing values is an iterative proce performing HOSVD on complete d ing the data by using the prediction cedure starts with randomly initia The third row shows the predicted which supresses over-fitting and imp they are meaningful and representing the characteristic structural information of the corresponding image classes. For example, Fig. 3 illustrates the interpretation of two entries of the discovered covariance: one from the row-covariance for digit "1" and the other from the column-covariance for digit "8". The resultant rowcovariance of "1" suggests that many rows of a "1"-image should be highly correlated, which is consistent with our knowledge of how "1" is written and is shown by the example images on the top row in Fig. 3. On the other a [ m p in c T w A ces to be fixed in learning, in contrast, BFCM explicitly employs a Bayesian model for the factor covariance matrices. We impose priors and perform inference on (R,R-1 ) and (C,C-1 ) directly. Probabilistic treatment is also employed by [Yu et al., 2006], where an inverse Wishart density is imposed on Gaussian process kernels. However, the kernels in that model operate separately on two sets of factors generating the relational matrix, which differs from the direct covariance learning in BFCM. 4 Experiments

3 RelatedWork Factorised covariance model is related to several existing lines of research in machine learning and statistics. The covariance matrix plays an essential role in latent factor models [Jolliffe, 2002; Tao et al., 2009], which analyse the covariance to design a small number of new variables (latent factors), so that varying the latent factors reproduces the joint heterogeneity of the observed variables. BFCM can be specialised as a covariance estimation step that precedes a traditional latent factor analysis. In particular, we can stack each matrix into a long column vector so all the data are transformed into a single matrix YD of [nv × N], where nv = nr×nc. BFCM is applied to YD: we solve for R and fix C = I so that the N samples are independent. Then the learned R summarises the correlations between all observed variables without considering the matrix structure, and a compact parameterisation of R leads to the traditional latent factors, e.g. a reduced rank SVD of R represents a Bayesian implementation of probabilistic PCA [Tipping and Bishop, 1999; Li and Tao, 2013]. Figure 3: Some revealed relations in the images 4 Experiments In this section, we demonstrate the usage and effectiveness of BFCM in three image analysis tasks. First, we use BFCM to learn the geometric structures of handwritten digits. Second, we predict pixels that are missing at different conditions. The third task is to learn new representation from a single image. The second and the third task show the advantage of learning factorised covariance in challenging problem settings. Discover image structure We first examine whether the model is effective in discovering knowledge about the population of an image class from a relatively small number of observations. In the Predict missing pixels In the second example, we learn from facial images to predict missi uses 50 facial images from the AR prediction is to relate the unkno One strategy is to exploit the con samples. For example, the problem sions in facial images has been cast regression task in [Wright et al., 2 formulation is effective when a full database is available. If all observ tain missing values, a probabilistic scale the analysis of connections fro to the constituent variables. For e ploy the probabilistic PCA model [ 1999] and learn the correlations bet els. Then the missing pixels can observed ones. Examples of applying this scheme (I) in Fig. 4. In the top row, the faces containing 10% randomly mis by dark dots). The second row sho by a learned probabilistic PCA m missing pixels are roughly consiste Figure 3: Some revealed relations in the images 4 Experiments In this section, we demonstrate the usage and effectiveness of BFCM in three image analysis tasks. First, we use BFCM to learn the geometric structures of handwritten digits. Second, we predict pixels that are missing at different conditions. The third task is to learn new representation from a single image. The second and the third task show the advantage of learning factorised covariance in challenging problem settings. Discover image structure We first examine whether the model is effective in discovering knowledge about the population of an image class from a relatively small number of observations. In the experiment, BFCM is used to analyse images of handwritten digits from the USPS dataset [Hull, 1994]. A BFCM is learned for each digit by using 10 randomly chosen images. Fig. 2 shows the row- and column-wise covariance matrices learned by the models. A brief investigation of the entries of the resultant covariance suggests that u p O s si re fo d t sc t pl 1 el o (I fa b b m tr vi in 4 se [Friedman et al., 2008], which is employed for the inference of BFCM, as well as many other regularised models [Guan et al., 2012]. In addition to sparse inverse, BFCM exploits the matrix form of the data, and assumes the convariance matrix is the Kronecker product of two covariance matrices: one of the rows and the other of the columns. Factorisation-models have also been employed in relational data analysis, where the observed array consists of records of the interactions among multiple groups of entities [Nickel et al., 2011]. In [Yu et al., 2006], a relational data model has been developed by defining multiple Gaussian processes on the participating entity groups. As a result, the joint covariance is given by the product of the covariance of those constituent Gaussian processes. The model can be alternatively motivated by considering the observations as the result of applying a (possibly infinite) multilinear transformation to a random core matrix of independent unit Gaussian entries [Yan et al., 2011]. Extensions to tensor data has also been considered in [Chu et al., 2009] and [Xu et al., 2012]. Most relational data models treat R and C as kernel matri-

los, 2002; Tao et al., 2007b; 2007a], it has been shown that

high-order SVD [Tucker, 1966] provides a natural and effec-

tive computational tool for multi-dimensional arrays of image

data.

Alternatively, reliable covariance estimators have been proposed by regularising the problem based on sparsity assumptions. For example, banded structure and thresholds have been proposed in [Bickel and Levina, 2008b] and [Bickel and Levina, 2008a], respectively. Regularisation of pairwise correlations has been proposed in our previous work [Li and Tao, 2012]. Learning a sparse inverse covariance has also been studied [Rothman et al., 2008], which corresponds to trimmed connections between the variables from the viewpoint of graphical models. Based on the efficient gradient

knowledge of how "1" is written and is shown by the example images on the top row in Fig. 3. On the other 1 http://www2.ece.ohio-state.ed ARdatabase.html Discover image structure We first show the model is effective in discovering knowledge about an image class from a relatively small number of samples. BFCM is used to analyse hand-written digits from USPS dataset [Hull, 1994]. A BFCM is learned for each digit by using 10 randomly chosen images.

Fig. 2 shows the learned row- and column-wise covariance matrices. A brief investigation of the covariance matrices suggests that they are meaningful and representing structural information of the corresponding image classes. For example, Fig. 3 illustrates the interpretation of two entries of the covariance: one from the row-covariance for digit "1" and the other from the column-covariance for digit "8". The rowcovariance of "1" suggests that most rows of a "1"-image are highly correlated, which is consistent with how "1" is written and is shown by the example images in Fig. 3. On the other hand, the marked entry of the column-covariance of "8" indicates that the corresponding pair of columns are negatively correlated. This is validated by the (transposed) example images , where the pen stokes of "8" cover different parts on the two columns. It is also noticeable that the symmetrical structures of the images manifest themselves in the anti-diagonal entries of the covariance matrices. Moreover, the images of "5/6" has similar row-covariance, which contrasts to that of "7/9".

In this experiment, the pixels are normalised so the average variance is 1. Diffuse prior parameters are used: α0 = 1, β0 = 0.01, τ0 = 1, ν0 = 1, and we let ρ = 0 to disable sparsity, because the demonstration focuses on the factorised covariance structure for images.

Predict missing pixels Secondly, we learn probabilistic models from facial images to predict missing pixels. The test uses 50 facial images from the AR dataset 1 . Essentially, prediction is to relate the unknowns to observations. One strategy is to exploit the connections between the samples. For example, the problem of overcoming occlusions in facial images has been cast as an ℓ1 regularised regression task in [Wright et al., 2009]. The regression formulation is effective when a fully observed reference database is available. If all observed samples may contain missing values, a probabilistic model is necessary to scale the analysis of connections from the samples down to the constituent variables. E.g., probabilistic PCA (PPCA) [Tipping and Bishop, 1999] can learn the correlations between individual pixels, so that missing pixels can be inferred from observed ones.

Examples of applying this scheme are shown in Group (I) in Fig. 4. In the top row, the five images contain 10% randomly missing pixels. The second row shows the images given by PPCA. The estimated missing pixels are roughly consistent with their ground-truth values according to the figure. However, there are visible artefacts, because the observed samples (50) is insufficient compared to the number of variables (50 × 40 = 2, 000) for learning PPCA reliably.

1 www2.ece.ohio-state.edu/∼aleix/ARdatabase.html

To address the small sample size problem, we represent the images as a multi-dimensional array by using high-order SVD (HOSVD) [Vasilescu and Terzopoulos, 2002]. The prediction of missing values is an iterative procedure: (i) performing HOSVD on complete data and (ii) completing the data by using the prediction of HOSVD. The procedure starts with randomly initialised missing pixels. The third row shows the results by HOSVD, which supresses over-fitting and improves the prediction.

We also test the matrix Gaussian process (MGP) [Yu et al., 2006; Yan et al., 2011] and BFCM, which also produce better result than PCA. Similar to HOSVD, MGP and BFCM avoid over-fitting by accounting for the array-structure of images and dealing with a smaller learning problem. On the other hand, MGP and BFCM consider the correlations in both the rows and the columns, and thus give more precise predictions. In Table 1, the left column compares the MSE of the missing value predictions, where the average variance of the pixels is normalised to 1.

Besides the over-fitting issue, the strategy of focusing relations between rows or columns instead of those between pixels helps the learning by expanding the applicability. A noticeable prerequisite of learning the pixel-wise covariance is that all pixels must be observed at least once. This is appropriate if the missing is at random, e.g. when data are lost due a poor communication channel. However, the assumption fails for an interesting scenario of image processing: when the images are acquired by a compromised sensor, which results missing pixels at fixed positions in all images.

Group (II) in Fig. 4 displays examples of predicting missing pixels at fixed positions. As above, the second row shows predicted images by PPCA. The model stops working, because it is infeasible to explicitly learn the covariance between pixel pairs, if one of them has never been observed. On the other hand, HOSVD, MGP and BFCM are not affected by fixing missing positions. The second column of Table 1 lists the MSEs. The predictions by PCA are invalid, and those by HOSVD, MGP and BFCM are similar to the previous test.

We also test the model by letting the missing happen in a fixed contiguous rectangle. Although the missing rectangle of 20 × 12 contains a similar number of missing pixels as in the previous tests, the task is more challenging. An image pixel is most likely to be related to its neighbours, and if a contiguous area has never been observed, it is difficult to obtain useful knowledge about the inner pixels.

Group (III) shows the example images and the prediction results. HOSVD gives unreliable predictions in this test. In principle, HOSVD considers the joint correlations between columns and rows. However, in implementation, when computing the correlations in one mode, e.g. those of rows, all rows in the images are treated as independent samples. The computation becomes unreliable when significant number of rows contain missing values at the same positions. MGP and BFCM suffer less from the contiguous missing pattern. In particular, BFCM explicitly learns the covariance and is better at discovering and exploiting the symmetrical structure between the columns in facial images. The visual assessment of the results is supported by the MSEs in the third column of Table 1.

In the experiment, we determine the model parameters by the performance in the first test (missing at random), including the rank of PCA and HOSVD, the kernel width and the latent dimension in MGP, and τ0 in BFCM (ν0, α0 and β0 are set as previously).

Learn representation by a single example Thirdly, we test model-based representation for the facial images. In particular, we challenge the model by using only a single image for training. This setting of sample scarcity seems to be somewhat extreme. However, it shows the strength of a factorised structure. One image contains many columns and rows. Thus it is possible to learn both the row and column covariance from a single image, which is an ill-defined problem for learning pixel-wise covariance. Once a BFC model is learned, a projection can be derived by following the rule of maximum variance, which resembles standard PCA. Specifically, the projection is given by performing SVD on the estimated covariance. After R and C are estimated, the SVD and the projection is obtained by (14). A similar projection can be derived for MGP by using the estimated kernels of the row and column Gaussian processes. We also include PCA 2 and HOSVD in this experiment. For classification, the simple nearest neighbour rule (NN) is used. The test is done on three face datasets (AR, YALE 3 and ORL 4 ) and the results are listed in Table 2. As a baseline, we also include results of applying NN on the raw data. The parameters of the algorithms in this experiment is determined by cross-validation.

The results shows that BFCM learns effective covariance from a single observation and leads to new representation. For discrimination, the new representation is superior to the raw data for the AR and YALE datasets, which contain relatively large appearance variations. On the other hand, the other models do not provide representations that are significantly more discriminative than the raw features. Hence direct regularisation of the covariance enables BFCM to learn a useful facial image model from a single example.

5 Conclusion

A Bayesian covariance estimator for matrix data has been developed. A factorised structure is adopted. The structure is motivated by meaningful relations in the row- and columnindexes of the matrix, and it helps reduce the covariance learning problem to a manageable scale for practical matrix data sizes.

Different from existing relational data models using factorised covariance, this work explicitly deals with the factor covariance matrices, treats them as random variables in a Bayesian model, imposes meaningful and sparsity-inducing priors and takes advantage of a fast optimisation algorithm for inference. In particular, the priors are defined w.r.t. the inverse of the factor covariance, consisting of a Wishart and a Laplace term. The components in the prior jointly encourage independency, which further regularises the learning problem and makes the model applicable for scarce observations.

Experiments demonstrate that the model is capable of learning useful relations for image understanding. The knowledge can be applied to infer missing data and to extract informative representation for subsequent classification tasks. Noticeably, in the former task, the model recovered a pixel even if it is missing in all images, and its relation to the remaining pixels is not directly available; in the latter task, the model learned effective representations from a single image.

2 PCA does not work on a single training image, and is computed using up to 100 images from 5 subjects excluded from test set.

3 cvc.yale.edu/projects/yalefaces/yalefaces.html

4 www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html