# Correlated Wishart Matrices Classification via an Expectation-Maximization
Composite Likelihood-Based Algorithm

Zhou Lan<EMAIL_ADDRESS>
Center for Outcome Research and Evaluation
Yale University
New Haven, CT 06510, USA

###### Abstract

Positive-definite matrix-variate data is becoming popular in computer vision.
The computer vision data descriptors in the form of Region Covariance
Descriptors (RCD) are positive definite matrices, which extract the key
features of the images. The RCDs are extensively used in image set
classification. Some classification methods treating RCDs as Wishart
distributed random matrices are being proposed. However, the majority of the
current methods preclude the potential correlation among the RCDs caused by
the so-called non-voxel information (e.g., subjects’ ages and nose widths,
etc). Modeling correlated Wishart matrices is difficult since the joint
density function of correlated Wishart matrices is difficult to be obtained.
In this paper, we propose an Expectation-Maximization composite likelihood-
based algorithm of Wishart matrices to tackle this issue. Given the numerical
studies based on the synthetic data and the real data (Chicago face data-set),
our proposed algorithm performs better than the alternative methods which do
not consider the correlation caused by the so-called non-voxel information.
All these above demonstrate our algorithm’s compelling potential in image set
classification in the coming future.

Keywords: Correlated Wishart Matrices, Composite Likelihood, Computer Vision,
Expectation–maximization algorithm, Image Set Classification, Region
Covariance Descriptor

##  1 Introduction

Positive definite matrix-variate data is a type of data in which each variable
is a positive definite matrix other than a scalar or a vector. The type of
data is involved in several applications, e.g., computer vision (see Tuzel et
al., 2006). The computer vision data descriptor in the form of Region
Covariance Descriptors (RCD) are positive definite matrices (Tuzel et al.,
2006). Generally, the specification of the RCD can be given as follows. Let
$I$ be an RGB color image. $I$ is usually represented as a $W\times H\times 3$
array, where $W$ and $H$ are the width and height of an image, respectively,
and $3$ is the number of color channels, i.e., red (R), green (G), and blue
(B). Let $F$ be the $W\times H\times D$ dimensional feature image extracted
from $I$. That is, the (x,y,d) element of the array $F$ is
$F(x,y,d)=\phi_{d}(x,y,{I})$. The function $\phi_{d}()$ is a mapping function
to extract intensity, color, gradient, filter response, etc. at the voxel
$(x,y)$ of the image $I$. Let $\bm{F}(x,y)$ be
$\bm{F}(x,y)=[F(x,y,1),...,F(x,y,D)]^{T}$. The RCD of the image $I$ is111In
several studies, the covariance matrix is calculated within certain areas of
an image.

$\bm{C}_{I}=\frac{1}{WH-1}\sum_{x=1}^{W}\sum_{y=1}^{H}(\bm{F}(x,y)-\bm{\mu})(\bm{F}(x,y)-\bm{\mu})^{T},$

where $\bm{\mu}$ is the sample mean of all $\bm{F}(x,y)$. $\bm{C}_{I}$, a
$D\times D$ covariance matrix, is a positive definite matrix-based descriptor
of the image $I$.

Compared to the traditional vector-based descriptors, many computer vision
studies have found that the RCD is a very useful quantity to describe the
distinguishing features of an image (e.g., Cherian et al., 2016), especially
for the applications to image set classification. This includes identifying
objects (e.g., Chen et al., 2020), identifying textures (e.g., Díaz-Pernas et
al., 2009), and facial recognition (e.g., Pang et al., 2008). Several
classification methods are proposed by using the RCDs of images as inputs.
Although the current works enjoy good classification, they may preclude the
available non-voxel information. The so-called non-voxel information is
defined as the information which is image-specific but not voxel-specific.
Taking the Chicago face data-set (Ma et al., 2015) as an example, besides the
RGB images of the subjects’ headshots, we also have rich image-level
information such as the subjects’ ages, nose widths, etc. We find that the
RCDs are correlated dependent on these pieces of non-voxel information. For
example, in Figure 1, we observe the RCDs’ largest eigenvalues varies
dependent on age and nose width. Therefore, the RCDs are correlated across the
non-voxel information. Although the correlation may or may not be explained
scientifically, we conjecture that modeling the variations of the RCDs across
the non-voxel information may increase the model fitness. In light of the
scientific objective, the better-fitted model may also produce more accurate
classification results.

Figure 1: The 3-D scatter plot shows the RCDs’ largest eigenvalues varies
dependent on age and nose width. The Z-axis is for the largest eigenvalues;
The X-axis and Y-axis are for the normalized ages and normalized nose widths,
respectively. The points represent the observations. The fitted surface is
fitted by a polynomial regression model.

The most commonly used probability distribution of positive definite matrix-
variate data is the Wishart distribution (Dryden et al., 2009; Cherian et al.,
2016; Lee and Schwartzman, 2017; Lan et al., 2020+). Modeling correlated
Wishart matrices is difficult. The most challenging component is, unlike
Gaussian distributed variables which you can feasibly use a multivariate
normal distribution, the joint density function of correlated Wishart matrices
is challenging to be obtained (Viraswami, 1991; Blumenson and Miller, 1963;
Smith and Garth, 2007). A desirable but challenging means is the so-called
(spatial) Wishart process (Gelfand et al., 2004; Lan et al., 2021+). The model
has a very elegant construction, that is to induce dependence of Wishart
matrices via the latent Gaussian processes. Unfortunately, this density
function of the (spatial) Wishart process is not analytically available under
general conditions (Viraswami, 1991; Blumenson and Miller, 1963; Smith and
Garth, 2007). Fortunately, standing on the shoulders of the giants who
provided the prodigious results (James et al., 1964), we can obtain the
analytic expression of the bivariate density function.

To mitigate this inference issue caused by the analytic form of the joint
density function, we propose a computationally-feasible composite likelihood-
based inference equipped with an expectation-maximization (EM) algorithm. We
first propose a hierarchical model based on the (spatial) Wishart process.
Given our proposed hierarchical model, we use the EM algorithm for model
fitting. However, due to the analytic form of the joint density function is
not available, both the E-step and M-step cannot be executed in a standard
way. Therefore, we replace the full likelihood with the composite likelihood
in the algorithm. The composite likelihood is an inference function derived by
multiplying a collection of component likelihood (Lindsay, 1988), which is
primarily to resolve the issue that the full likelihood is unavailable.
Because the analytic density function of two correlated Wishart matrices can
be obtained (James et al., 1964), the composite likelihood is in a pairwise
setting, which is then derived by multiplying a collection of possible pairs.
The estimators based on composite likelihood is asymptotically unbiased (Varin
et al., 2011), but the efficiency of the estimators is determined by the
likelihood-weights (Padoan et al., 2010; Bevilacqua et al., 2012; Bai et al.,
2012, 2014). Considering a trade-off between parameter estimation and
computational cost, we further proposed a novel likelihood-weight function,
enjoying accurate parameter estimation and fast computational speed.

In light of our goal of classification, we use both the simulated data and the
real data to demonstrate that our proposed method produces a more accurate
classification of positive definite matrices than the benchmark methods. The
simulated data are generated from the proposed model. The real data is the
Chicago face data (Ma et al., 2015). The benchmark methods include a wide
range of methods, including Wishart Bayesian nonparametrics (Wishart BNP)
(Cherian et al., 2016), Log-Euclidean metric method (Huang et al., 2015),
K-means, and Gaussian mixture model. The numerical results show that our
proposed method produces a more accurate classification of positive definite
matrices, demonstrating its compelling potential in image set classification
in the coming future.

In the rest of the paper, we first give our hierarchical model in Section 2.
The estimation method is provided in Section 3. We provide simulation studies
in Section 4, showing that our algorithm performs better than the other
benchmark methods. A real data application using the Chicago face database (Ma
et al., 2015) with the objective of race classification is provided in Section
5. In Section 6, we conclude with a discussion. The essential derivations and
the codes written in MATLAB for implementing our method are given in the
appendices.

##  2 Model

In this section, we illustrate our model construction step-by-step. First, we
give the basic model setup in Section 2.1. The basic model describes how the
correlation of positive definite matrices across the non-voxel information is
constructed. In Section 2.2, we further introduce the latent variables which
can be used for the purpose of classification. Finally in Section 2.3, we give
a model summary and introduce its relationship to the RCD-based image set
classification.

Our model is based on the the parameterized Wishart distribution. Note that
the parameterized Wishart distribution is slightly different from the Wishart
distribution introduced in the classic textbooks. The parameterized Wishart
distribution is constructed as follows. $\bm{Z}_{j}$ independently follows a
mean-zero normal distribution with the covariance matrix $\bm{\Sigma}$,
denoted as $\bm{Z}_{j}\stackrel{{\scriptstyle
ind}}{{\sim}}\mathcal{N}(\bm{0},\bm{\Sigma})$.
$\bm{W}=\sum_{j=1}^{M}\bm{Z}_{j}\bm{Z}_{j}^{T}/M$ follows a parameterized
Wishart distribution with the mean matrix $\bm{\Sigma}$ and the degrees of
freedom $M$, denoted as $\bm{W}\sim\mathcal{W}_{p}(\bm{\Sigma},M)$. The
density function of $\bm{W}$ is given in Appendix A.

###  2.1 Basic Model Setup

Let $\bm{A}_{t}$ be a $p\times p$ positive definite matrix-random variable
following a parameterized Wishart distribution with the mean matrix
$\bm{\Sigma}_{t}$ and the degrees of freedom $M$, denoted as

$\bm{A}_{t}\sim\mathcal{W}_{p}(\bm{\Sigma}_{t},M).$ (1)

$\bm{A}_{t}$ can be treated as the RCD of the image $t$. Our goal is to build
a model which induces the correlation between $\bm{A}_{t}$ and $\bm{A}_{s}$
($t\not=s$). To achieve this goal, We decompose $\bm{A}_{t}$ as
$\bm{A}_{t}=\bm{L}_{t}\bm{U}_{t}\bm{L}_{t}^{T}$, where $\bm{L}_{t}$ is the
lower Cholesky factor of $\bm{\Sigma}_{t}$ such as
$\bm{\Sigma}_{t}=\bm{L}_{t}\bm{L}_{t}^{T}$. Since $\bm{\Sigma}_{t}$ determines
the mean of $\bm{A}_{t}$, we treat $\bm{U}_{t}$ as the terms inducing the
correlation. A more approachable interpretation of this decomposition is that
$\bm{\Sigma}_{t}$ is the mean term driving the mean and $\bm{U}_{t}$ is the
residual term explaining the correlation among the matrices (Lan et al.,
2021+). Next, we give the details of constructing these two terms.

A means to construct the correlation among $\bm{U}_{t}$ has been known as
(spatial) Wishart process (Gelfand et al., 2004; Lan et al., 2021+). Let
$\bm{Z}_{t}$ be a mean-zero Gaussian process with the covariance as
$\text{cov}(\bm{Z}_{t},\bm{Z}_{s})=\rho_{ts}\times\bm{\Sigma}$, where
$\rho_{ts}\in(0,1)$ is the correlation, and $\bm{\Sigma}$ is the cross
covariance matrix of the Gaussian process. In the rest of the paper, we simply
give $\bm{\Sigma}=\bm{I}$, where $\bm{I}$ is an identity matrix. If we have
independent and identically distributed realizations of this process, denoted
as $\bm{Z}_{jt}$, where $j$ is the index for realizations and $t$ is the index
for the Gaussian process, then
$\bm{U}_{t}=\sum_{j=1}^{M}\bm{Z}_{jt}\bm{Z}_{jt}^{T}/M$ follows a so-called
Wishart process, denoted as
$\bm{U}_{t}\sim\mathcal{WP}(M,\rho_{ts}\otimes\bm{I})$. The correlation
$\rho_{ts}$ can be expressed as a correlation function, i.e.,
$\rho_{ts}=\mathcal{K}(||\bm{X}_{t}-\bm{X}_{s}||;\bm{\phi})$, where
$\mathcal{K}(d;\bm{\phi})$ is a correlation function (e.g., Matern,
exponential) of $d$, and $\bm{\phi}$ is a vector of the function parameters.
Let $\bm{X}_{t}$ be a $d\times 1$ vector representing the non-voxel
information of the image $t$ (e.g., age, nose width, etc). The distance
$||\bm{X}_{t}-\bm{X}_{s}||$ is the Euclidean distance of the two images’ non-
voxel information.

We emphasize that this model construction has two important properties: (1)
$\bm{A}_{t}$ preserves to follow $\mathcal{W}_{p}(\bm{\Sigma}_{t},M)$
marginally; (2) $\bm{A}_{t}$ and $\bm{A}_{s}$ are correlated ($t\not=s$) due
to the residual terms. The property (1) can be validated by giving the lemma
below.

###### Lemma 1

Let $\bm{\Sigma}$ be a $p\times p$ positive definite matrix. $\bm{L}$ is the
lower Cholesky factor of $\bm{\Sigma}$ such as $\bm{\Sigma}=\bm{L}\bm{L}^{T}$.
If $\bm{U}\sim\mathcal{W}_{p}(\bm{I},M)$, then
$\bm{W}=\bm{L}\bm{U}\bm{L}^{T}\sim\mathcal{W}_{p}(\bm{\Sigma},M)$.

This lemma is borrowed from Proposition 8.1 of Eaton (1983, Chapter 8: The
Wishart Distribution) with a proof, and it can also be found in many other
multivariate textbooks. To validate the property (2), we can quantify the
residual dependence of the Wishart process by using the expected squared
Frobenius norm (Lan et al., 2021+), expressed as

$\displaystyle\mathbb{E}||\bm{U}_{t}-\bm{U}_{s}||_{F}^{2}=\gamma(M)(1-\rho_{ts}^{2})$
(2)
$\displaystyle\gamma(M)=\frac{2}{M}\Big{(}Tr(\bm{I}\bm{I})+Tr(\bm{I})Tr(\bm{I})\Big{)},$

where $||\bm{U}_{t}-\bm{U}_{s}||_{F}^{2}$ is the squared Frobenius norm of the
term $\bm{U}_{t}-\bm{U}_{s}$. The norm can also be treated as a variogram
(Cressie, 1992). Given the variogram, it is more transparent to claim that a
larger correlation of the underlying Gaussian process leads to a larger
residual dependence, and a larger value of $M$ leads to a smaller local
variation.

###  2.2 Latent Variable Setup

In light of our scientific objective to classify the RCDs, the mean matrix
$\bm{\Sigma}_{t}$ is specified as a term dependent on a latent group variable.
It is specified as $\bm{\Sigma}_{t}|[Z_{t}=k]\equiv\mathbf{S}_{k}$, where
$Z_{t}\in\\{1,2,...,K\\}$ is the latent group index. The latent group variable
$Z_{t}$ independently follows a discrete distribution such that
$Pr(Z_{t}=k)=\omega_{k}$, denoted as

$Z_{t}\stackrel{{\scriptstyle
ind}}{{\sim}}\text{Discrete}([1,2,..,K],[\omega_{1},\omega_{2},...,\omega_{K}]).$
(3)

The hyperparameters $\\{\mathbf{S}_{k}:k=1,2,...,K\\}$ are the known and given
values describing the features of each group ($k\in\\{1,2,...,K\\}$) in terms
of positive definite matrices. The hyperparameters
$\\{\mathbf{S}_{k}:k=1,2,...,K\\}$ can be obtained from the training data
and/or the experienced values. The number of components $K$ is also given and
fixed. Based on this hierarchical model setup, $Z_{t}$ can be used as a group
indicator to classify the image $t$ to a group $k\in\\{1,2,...,K\\}$.

###  2.3 Model Summary and the RCD-Based Image Set Classification

To this end, our proposed hierarchical model framework is summarized as

$\displaystyle\textbf{Basic Level}:\
\bm{A}_{t}=\bm{L}_{t}\bm{U}_{t}\bm{L}_{t}^{T}\quad\bm{\Sigma}_{t}=\bm{L}_{t}\bm{L}_{t}^{T}\quad\bm{U}_{t}\sim\mathcal{WP}(M,\rho_{ts}\otimes\bm{I})\quad\rho_{ts}=\mathcal{K}(||\bm{X}_{t}-\bm{X}_{s}||;\bm{\phi})$
(4) $\displaystyle\textbf{Latent Level:}\
\bm{\Sigma}_{t}|[Z_{t}=k]\equiv\mathbf{S}_{k}\quad
Z_{t}\stackrel{{\scriptstyle
ind}}{{\sim}}\text{Discrete}([1,2,..,K],[\omega_{1},\omega_{2},...,\omega_{K}]).$

The above illustration gives a classification model for the correlated
positive definite matrices, which may be still abstract to the practitioners.
Next, we elaborate our proposed model by giving its relationship to the RCD-
based image set classification. The image set classification is an important
research field in computer vision. The objective is to classify the pictures
into different topics (e.g., Cherian et al., 2016). We use the Chicago face
database (Ma et al., 2015) as an illustration example. The data-set provides
high-resolution, standardized headshots of male and female faces of different
races between the ages of 17-65, as well as the subjects’ non-voxel
information (e.g., age, nose width, etc). We have a scientific goal to
classify the pictures by their races, i.e., Black, Latino, and White, which is
one of the common scientific goals in computer vision studies (e.g., Ou et
al., 2005; Roomi et al., 2011; Fu et al., 2014).

Let $\bm{A}_{t}$ be the RCD of the face picture $t$. First, the residual terms
in the basic level capture the variations caused by the non-voxel information.
Because the non-voxel information (e.g., age, nose width, etc) is usually
irrelevant to the races, our model specification indicates that the residual
dependence is separated from the effect of races which are driven by the mean
terms.

Given that the potential correlation among the random positive definite
matrices caused by the non-voxel information has been captured by our proposed
model, we can classify the data with fewer noises. We use $k=\\{1,2,3\\}$ to
represent Black, Latino, and White, respectively. The hyperparameters
$\\{\mathbf{S}_{k}:k=1,2,...,K\\}$ are plug-in values that can be obtained
from the training data and/or the experienced value, describing the typical
features of a race in terms of RCD. The maximum a posteriori estimator of
$Z_{t}$ can be used for the classification of images, resolving the scientific
objective.

##  3 Estimation

In this section, we describe our approach to obtain the parameter estimates of
the unknowns. Our estimation approach is essentially a hybrid method composed
of the composite likelihood method and the EM algorithm.

###  3.1 Standard EM Algorithm

Given Model (4), we first provide the EM Algorithm under the standard
situations. Rigorously, we use $[\bm{a}_{1},...,\bm{a}_{t}]$ to represent the
observed positive definite matrices of the random matrices
$[\bm{A}_{1},...,\bm{A}_{t}]$. We give that
$\bm{a}=[\bm{a}_{1},...,\bm{a}_{t}]$ and $\bm{A}=[\bm{A}_{1},...,\bm{A}_{t}]$.
We also give $\bm{Z}=[Z_{1},...,Z_{T}]$ to denote all latent group indexes.
The parameters to be estimated are
$\bm{\omega}=[\omega_{1},...,\omega_{k},...,\omega_{K}]$ and
$\bm{\Phi}=[\bm{\phi}^{T},M]^{T}$, where
$[\omega_{1},...,\omega_{k},...,\omega_{K}]$ are the group-weights,
$\bm{\phi}$ is a vector of the parameters in the correlation function
$\mathcal{K}(d;\bm{\phi})$, and $M$ is the degrees of freedom. We give the
notation $\bm{\theta}=[M,\bm{\omega}^{T},\bm{\Phi}^{T}]^{T}$ which contains
all the parameters needed to be estimated. The full log-likelihood is

$\displaystyle\ell_{Full}(\bm{\theta}|\bm{a},\bm{Z})$
$\displaystyle=\log\Bigg{[}\prod_{j_{1}=1}^{K}...\prod_{j_{t}=1}^{K}...\prod_{j_{T}=1}^{K}\bigg{(}f_{\bm{A}}(\bm{a}|\bm{\Sigma}_{1}=\bm{S}_{j_{1}},...,\bm{\Sigma}_{t}=\bm{S}_{j_{t}},...,\bm{\Sigma}_{T}=\bm{S}_{j_{T}},\bm{\Phi})$
(5)
$\displaystyle\times\prod_{t=1}^{T}\omega_{j_{t}}\bigg{)}^{\mathbb{I}(Z_{1}=j,,...,Z_{t}=j_{t},...,Z_{T}=j_{T})}\Bigg{]},$

where
$f_{\bm{A}}(\bm{a}|\bm{\Sigma}_{1}=\bm{S}_{j_{1}},...,\bm{\Sigma}_{t}=\bm{S}_{j_{t}},...,\bm{\Sigma}_{T}=\bm{S}_{j_{T}},\bm{\theta})$
is the joint density function of $\bm{A}$ conditional on that
$[Z_{1}=j_{1},...,Z_{t}=j_{t},...,Z_{T}=j_{T}]$.

Our goal is to obtain the maximum likelihood estimates for $\bm{\theta}$.
Therefore, we have to maximize the likelihood whose the latent variables
$\bm{Z}$ are integrated out. The log-likelihood $\ell(\bm{\theta}|\bm{a})$
whose latent variables are integrated out is expressed as

$\displaystyle\ell(\bm{\theta}|\bm{a})$
$\displaystyle=\int...\int...\int\log\Bigg{[}\prod_{j_{1}=1}^{K}...\prod_{j_{t}=1}^{K}...\prod_{j_{T}=1}^{K}\bigg{(}f_{\bm{A}}(\bm{a}|\bm{\Sigma}_{1}=\bm{S}_{j_{1}},...,\bm{\Sigma}_{t}=\bm{S}_{j_{t}},...,\bm{\Sigma}_{T}=\bm{S}_{j_{T}},\bm{\Phi})$
(6)
$\displaystyle\times\prod_{t=1}^{T}\omega_{j_{t}}\bigg{)}^{\mathbb{I}(Z_{1}=j,...,Z_{t}=j_{t},...,Z_{T}=j_{T})}\Bigg{]}dZ_{1}...dZ_{t}...dZ_{T},$

The hurdle caused by the unobserved latent variables $\bm{Z}$ can be easily
overcome through the famous EM algorithm. The iterative scheme of the
algorithm is given as follows. If we have the current parameter
$\bm{\theta}^{(r)}$, the likelihood is
$\ell_{Full}(\bm{\theta}^{(r)}|\bm{a},\bm{Z})$. We first compute the expected
value of the likelihood function with respect to the current conditional
distribution of $\bm{Z}$ given the data $\bm{a}$ and the current parameter
estimates $\bm{\theta}^{(r)}$, denoted as
$Q(\bm{\theta}|\bm{\theta}^{(r)})=\mathbb{E}_{\bm{Z}|\bm{A},\bm{\theta}^{(r)}}\ell_{Full}(\bm{\theta}|\bm{a},\bm{Z})$.
We call this step as E-step since we take the expectation. Next we maximize
the term $Q(\bm{\theta}|\bm{\theta}^{(r)})$ with respect to $\bm{\theta}$ to
obtain the next estimates $\bm{\theta}^{(r+1)}$. This step is called M-step
because we maximize the function. We iteratively repeat the two steps until
convergence. Finally, $\hat{\bm{\theta}}$ is obtained as the maximum
likelihood estimates for $\bm{\theta}$.

###  3.2 Composite Likelihood-Based Inference

The standard EM algorithm has a bottleneck which is caused by the joint
density function
$f_{\bm{A}}(\bm{a}|\bm{\Sigma}_{1}=\bm{S}_{j_{1}},...,\bm{\Sigma}_{t}=\bm{S}_{j_{t}},...,\bm{\Sigma}_{T}=\bm{S}_{j_{T}},\bm{\Phi})$.
Given several reports (Viraswami, 1991; Blumenson and Miller, 1963; Smith and
Garth, 2007), an analytic expression of the joint density function is not
available under general conditions. However, standing on the shoulders of the
giants who provided the prodigious results (James et al., 1964), we are able
to obtain the bivariate density, i.e., the conditional density function of
$[\bm{A}_{t},\bm{A}_{s}]$ given $[Z_{t},Z_{s}]$, denoted as
$f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi})$.
The details in obtaining the density function are summarized in Appendix B.
The expression of
$f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi})$
is written as

$\displaystyle
f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi})$
(7) $\displaystyle=\
_{0}F_{1}\left(\frac{1}{2}M;\frac{1}{4}\left(\frac{M\rho_{ts}}{1-\rho_{ts}^{2}}\right)^{2}\bm{Q}_{j_{s}}^{-1}\bm{a}_{s}(\bm{Q}_{j_{s}}^{T})^{-1}\bm{Q}_{j_{t}}^{-1}\bm{a}_{t}(\bm{Q}_{j_{t}}^{T})^{-1}\right)$
$\displaystyle\times
2^{-Mp}\times\left[\Gamma_{p}(\frac{1}{2}M)\right]^{-2}\times|\bm{S}_{j_{t}}^{-1}\bm{a}_{t}|^{\frac{1}{2}(M-p-1)}\times|\bm{S}_{j_{s}}^{-1}\bm{a}_{s}|^{\frac{1}{2}(M-p-1)}$
$\displaystyle\times\left(\frac{1-\rho_{ts}^{2}}{M^{2}}\right)^{-\frac{pM}{2}}\times\exp\left(-\frac{1}{2}tr\left[M\frac{1}{1-\rho_{ts}^{2}}\left(\bm{S}_{j_{s}}^{-1}\bm{a}_{s}+\bm{S}_{j_{t}}^{-1}\bm{a}_{t}\right)\right]\right)$
$\displaystyle\times|\bm{S}_{j_{t}}^{-1}|^{\frac{p+1}{2}}\times|\bm{S}_{j_{s}}^{-1}|^{\frac{p+1}{2}},$

where ${}_{0}F_{1}(;)$ is a hypergeometric function of matrix argument (See
James et al., 1964, Section 6) and $\bm{Q}_{j_{t}}$ is the lower Cholesky
decomposition factor of $\bm{S}_{j_{t}}$ such as
$\bm{S}_{j_{t}}=\bm{Q}_{j_{t}}\bm{Q}_{j_{t}}^{T}$. The value of this function
${}_{0}F_{1}(;)$ can be numerically evaluated (Koev and Edelman, 2006).

Here, we want to highlight two important properties about this density
function: (1) when $\rho_{ts}=0$,
$f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi})$
is the product of the probability density functions of
$\mathcal{W}_{p}(\bm{S}_{j_{t}},M)$ and $\mathcal{W}_{p}(\bm{S}_{j_{s}},M)$
(see the density function of parameterized Wishart distribution in Appendix
A); (2) conversely, $\rho_{ts}=1$ leads to that
$f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi})$
is an improper density function. The property (1) is consistent with our basic
model setup because $\rho_{ts}=0$ implies that there is no residual
dependence. The property (2) may be analogised to the similar issue
encountered in Gaussian process modeling (Williams and Rasmussen, 2006), that
$\rho_{ts}=1$ leads to the improper multivariate Gaussian density function.

Obtaining the expression of the bivariate joint density is a milestone which
allows us to proceed the next composite likelihood-based inference. Composite
likelihood is an inference function derived by multiplying a collection of
component likelihood. For example, we denote
${\mathcal{A}_{1},...,\mathcal{A}_{C}}$ as a set of marginal or conditional
events with associated likelihoods $\mathcal{L}_{c}(\bm{\theta};\bm{y})\propto
f(\bm{y}\in\mathcal{A}_{c};\bm{\theta})$, where $f()$ is a density function,
$\bm{y}$ is the data, and $\bm{\theta}$ is the parameters. Following Lindsay
(1988), a composite likelihood is a weighted product such as
$\prod_{c=1}^{C}\mathcal{L}_{c}(\bm{\theta};\bm{y})^{p_{c}}$, where
$\\{p_{c}:k=1,2,...,C\\}$ are non-negative weights. Under the usual regularity
conditions, the composite likelihood-based inference can provide
asymptotically unbiased parameter estimates when standard likelihood
estimators are not available (Varin et al., 2011). For example, Padoan et al.
(2010) proposed a pairwise setting for max-stable processes by using composite
likelihood-based inference, since only the analytically tractable form can be
obtained in their case.

We define the events $\mathcal{A}_{c}=\\{\bm{A}_{t},\bm{A}_{s}\\}$ as the sets
of bivariate positive definite matrices taken over all $T(T-1)/2$ distinct
pairs $t$ and $s$. Next we replace the full log-likelihood
$\ell_{Full}(\bm{\theta}|\bm{a},\bm{Z})$ with the composite log-likelihood
$\ell_{Pair}(\bm{\theta}|\bm{a},\bm{Z})$. The log composite likelihood
$\ell_{Pair}(\bm{\theta}|\bm{a},\bm{Z})$ in the pairwise setting is written as

$\displaystyle\ell_{Pair}(\bm{\theta}|\bm{a},\bm{Z})=\sum_{t=1}^{T-1}\sum_{s=t+1}^{T}p_{ts}\ell_{ts}(\bm{\theta}|\bm{a},\bm{Z}),$
(8)

where

$\ell_{ts}(\bm{\theta}|\bm{a},\bm{Z})=\log\Bigg{[}\prod_{j_{t}=1}^{K}\prod_{j_{s}=1}^{K}\bigg{(}f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi})\times\omega_{j_{t}}\omega_{j_{s}}\bigg{)}^{\mathbb{I}(Z_{t}=j_{t},Z_{s}=j_{s})}\Bigg{]},$

$p_{ts}>0$ is its likelihood-weight. We usually scale the likelihood-weights
so that $\sum_{t<s}p_{ts}=1$.

The estimators obtained via composite likelihood-based inference are
asymptotically unbiased but not asymptotically efficient (Varin et al., 2011;
Padoan et al., 2010). A key to resolve this issue is the choice of the
likelihood-weights $\\{p_{ts}\\}$ that minimize the total amount of variation
of the asymptotic variance (Heyde, 2008, Section 2.3). Padoan et al. (2010)
proposed the $0/1$ likelihood-weights where a threshold $\delta>0$ is
introduced. If $||\bm{X}_{t}-\bm{X}_{s}||>\delta$, then $p_{ts}=0$, otherwise,
$p_{ts}=1$ (not scaled). This method is computationally-feasible because many
pairs are excluded from computation. Bai et al. (2012) and Bevilacqua et al.
(2012) further propose adaptive weights which rely on the analytic expression
of the Fisher information of estimators.

To combine the two ideas above, we propose a novel smooth likelihood-weight
function which is simple and computationally-fast for the practitioners but
produces relatively efficient estimators. The proposed likelihood-weights are
$p_{ts}\propto\exp(-||\bm{X}_{t}-\bm{X}_{s}||/\lambda)\times u_{ts}$. The
variable $u_{ts}\in\\{0,1\\}$ is pre-specified: among all $u_{ts}$, we
randomly select $\left\lfloor u\times T(T-1)/2\right\rfloor$ pairs which have
$u_{ts}=1$, where $u\in(0,1)$; the rest pairs have $u_{ts}=0$. This
likelihood-weight function has two tuning parameters. $\lambda>0$ controls the
decay rate of the likelihood-weight. That is, a smaller $\lambda$ leads to
more homogeneous non-zero likelihood-weights (See Figure 2). $u\in(0,1)$
controls the computational cost and reduces the pairs to be evaluated from
$T(T-1)/2$ to $u\times T(T-1)/2$. In the numerical studies, we will
investigate the roles of the two tuning parameters on parameter estimation and
classification accuracy.

###  3.3 Composite Likelihood-Based EM Algorithm

Combining all the results above, we give the final EM composite likelihood-
based algorithm, which is a hybrid method composed of the EM algorithm and
composite likelihood-based inference. The key idea of this hybrid method is to
replace $\ell_{Full}(\bm{\theta}|\bm{a},\bm{Z})$ by
$\ell_{Pair}(\bm{\theta}|\bm{a},\bm{Z})$ in the algorithm. Thus, the step of
conditional expectation is with the posterior based on composite likelihood,
that is
$p(\bm{Z}|\bm{A},\bm{\theta}^{(r)})\propto\ell_{Pair}(\bm{\theta}^{(r)}|\bm{a},\bm{Z})p(\bm{Z})$.
The corresponding E-step and M-step are described as follows:

#### E-Step:

In this step, we aim to obtain
$\mathbb{E}_{\bm{Z}|\bm{A},\bm{\theta}^{(r)}}\ell_{Pair}(\bm{\theta}|\bm{a},\bm{Z})$.
It is expressed as

$\displaystyle
Q_{Pair}(\bm{\theta}|\bm{\theta}^{(r)})=\mathbb{E}_{\bm{Z}|\bm{A},\bm{\theta}^{(r)}}\ell_{Pair}(\bm{\theta}|\bm{a},\bm{Z})$
(9)
$\displaystyle=\mathbb{E}_{\bm{Z}|\bm{A},\bm{\theta}^{(r)}}\Bigg{(}\sum_{t=1}^{T-1}\sum_{s=t+1}^{T}p_{ts}\log\Bigg{[}\prod_{j_{t}=1}^{K}\prod_{j_{s}=1}^{K}\bigg{(}f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi})\times\omega_{j_{t}}\omega_{j_{s}}\bigg{)}^{\mathbb{I}(Z_{t}=j_{t},Z_{s}=j_{s})}\Bigg{]}\Bigg{)}$
$\displaystyle=\mathbb{E}_{\bm{Z}|\bm{A},\bm{\theta}^{(r)}}\Bigg{(}\sum_{t=1}^{T-1}\sum_{s=t+1}^{T}p_{ts}\sum_{j_{t}=1}^{K}\sum_{j_{s}=1}^{K}\log\bigg{(}f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi})\times\omega_{j_{t}}\omega_{j_{s}}\bigg{)}{\mathbb{I}(Z_{t}=j_{t},Z_{s}=j_{s})}\Bigg{)}$
$\displaystyle=\sum_{t=1}^{T-1}\sum_{s=t+1}^{T}p_{ts}\sum_{j_{t}=1}^{K}\sum_{j_{s}=1}^{K}\log\bigg{(}f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi})\times\omega_{j_{t}}\omega_{j_{s}}\bigg{)}\mathbb{E}_{\bm{Z}|\bm{A},\bm{\theta}^{(r)}}{\mathbb{I}(Z_{t}=j_{t},Z_{s}=j_{s})}.$

This means that we only need to calculate
$T_{\binom{j_{t},j_{s}}{t,s}}^{(r)}=\mathbb{E}_{\bm{Z}|\bm{A},\bm{\theta}^{(r)}}\mathbb{I}(Z_{t}=j_{t},Z_{s}=j_{s})$
and plug them in. The term $T_{\binom{j_{t},j_{s}}{t,s}}^{(r)}$ is explicitly
expressed as

$\displaystyle T_{\binom{j_{t},j_{s}}{t,s}}^{(r)}$
$\displaystyle=\mathbb{E}_{\bm{Z}|\bm{A},\bm{\theta}^{(r)}}\mathbb{I}(Z_{t}=j_{t},Z_{s}=j_{s})$
(10)
$\displaystyle=\frac{\Big{(}f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi}^{(r)})\times\omega_{j_{t}}^{(r)}\omega_{j_{s}}^{(r)}\Big{)}^{p_{ts}}}{\sum_{j_{t}=1}^{K}\sum_{j_{s}=1}^{K}\bigg{(}f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi}^{(r)})\times\omega_{j_{t}}^{(r)}\omega_{j_{s}}^{(r)}\bigg{)}^{p_{ts}}}.$

By plugging in $T_{\binom{j_{t},j_{s}}{t,s}}^{(r)}$, we have the updated
$Q_{Pair}(\bm{\theta}|\bm{\theta}^{(r)})$ expressed as

$\displaystyle
Q_{Pair}(\bm{\theta}|\bm{\theta}^{(r)})=\mathbb{E}_{\bm{Z}|\bm{A},\bm{\theta}^{(r)}}\ell_{Pair}(\bm{\theta}|\bm{a},\bm{Z})$
(11)
$\displaystyle=\sum_{t=1}^{T-1}\sum_{s=t+1}^{T}p_{ts}\log\Bigg{[}\prod_{j_{t}=1}^{K}\prod_{j_{s}=1}^{K}\bigg{(}f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi})\times\omega_{j_{t}}\omega_{j_{s}}\bigg{)}^{T_{\binom{j_{t},j_{s}}{t,s}}^{(r)}}\Bigg{]}$

#### M-Step:

We obtain $\bm{\theta}^{(r+1)}$ via maximizing
$Q_{Pair}(\bm{\theta}|\bm{\theta}^{(r)})$ with respect to $\bm{\theta}$. To be
specific, $\bm{\omega}$ and $\bm{\Phi}$ can be parallelly estimated. The
estimate of the weight $\bm{\omega}$ has an analytic expression, expressed as

$\omega_{k}^{(r+1)}=\frac{\sum_{t=1}^{T-1}\sum_{s=t+1}^{T}\sum_{j_{t}=1}^{K}T_{\binom{j_{t},k}{t,s}}^{(r)}}{\sum_{t=1}^{T-1}\sum_{s=t+1}^{T}\sum_{j_{t}=1}^{K}\sum_{j_{s}=1}^{K}T_{\binom{j_{t},j_{s}}{t,s}}^{(r)}}.$
(12)

$\bm{\Phi}^{(r+1)}$ is obtained via maximizing the function

$Q_{Pair}(\bm{\Phi}|\bm{\Phi}^{(r)})=\sum_{t=1}^{T-1}\sum_{s=t+1}^{T}p_{ts}\log\Bigg{[}\prod_{j_{t}=1}^{K}\prod_{j_{s}=1}^{K}\bigg{(}f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi})\bigg{)}^{T_{\binom{j_{t},j_{s}}{t,s}}^{(r)}}\Bigg{]},$
(13)

with respect to $\bm{\Phi}$. The function
$Q_{Pair}(\bm{\Phi}|\bm{\Phi}^{(r)})$ can be feasibly maximized using the
Quasi-Newton method (Dennis and Moré, 1977).

We obtain the maximum likelihood estimate $\hat{\bm{\theta}}$ until
convergence. The conditional expectation of $\mathbb{I}(Z_{t}=k)$, denoted as
$G_{t,k}$, can be used as the classifier to show the probability that the
image $t$ belongs to the group $k$, expressed as

$G_{t,k}=\mathbb{E}_{\bm{Z}|\bm{A},\hat{\bm{\theta}}}\mathbb{I}(Z_{t}=k)=\frac{\sum_{s=1}^{T}\sum_{j_{s}=1}^{K}\bigg{(}f(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{k},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\hat{\bm{\Phi}})\times\hat{\omega}_{j_{t}}\hat{\omega}_{k}\bigg{)}^{p_{ts}}}{\sum_{s=1}^{T}\sum_{j_{t}=1}^{K}\sum_{j_{s}=1}^{K}\bigg{(}f(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{s}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\hat{\bm{\Phi}})\times\hat{\omega}_{j_{t}}\hat{\omega}_{j_{s}}\bigg{)}^{p_{ts}}}.$
(14)

We classify the observation $t$ to the group
$k_{t}=\max_{k=1,2,...,K}G_{t,k}$.

The composite likelihood-based EM algorithm borrows the ideas from the two
popular statistical methods. We are concerned about if our algorithm still
enjoys the properties of the original methods. Fortunately, there have been
many relevant works (Choi and Briollais, 2011; Gao and Song, 2011; Vasdekis et
al., 2012; Chen et al., 2016) endorsing our proposed algorithm. The major
concern is that whether the hybrid method still enjoys the three key
properties of standard EM algorithm. As shown in detail by Gao and Song
(2011), it is easy to prove the three key properties which are (i) the ascent
property that is
$\ell_{Pair}(\bm{\theta}^{{(r+1)}}|\bm{a})>\ell_{Pair}(\bm{\theta}^{(r)}|\bm{a})$,
where
$\ell_{Pair}(\bm{\theta}|\bm{a})=\int\ell_{Pair}(\bm{\theta}^{{(r+1)}}|\bm{a},\bm{Z})d\bm{Z}$,
(ii) convergence to a stationary point of the objective function, and (iii)
convergence rate depending on the curvature of the objective function. The
MATLAB codes implementing this hybrid algorithm are attached in the Appendix
C.

##  4 Simulation Studies

In this section, we use synthetic data to study the performance of our
proposed algorithm. The synthetic data are generated based on the full model
(Model 4), where the correlation function is the exponential function, i.e.,
$\mathcal{K}(d;\phi)=\exp(-d/\phi)$. The range parameter $\phi>0$ controls the
correlation, and a larger $\phi$ leads to a larger correlation. In each
replication, $T=50$ positive definite matrices are generated for
classification. The residual terms $U_{t}$ are generated by setting the
degrees of freedom as $M=5$ and setting range parameter as $\phi=1$. The
latent variable $Z_{t}$ are sampled equally from $\\{1,2,3\\}$, indicating
$K=3$ and $\omega_{1}=\omega_{2}=\omega_{3}=\frac{1}{3}$. In each replication,
the mean $\bm{S}_{k}$ is independently sampled from
$\mathcal{W}_{3}(\bm{I},3)$. The covariate is
$\bm{X}_{t}=[X_{t1},...,X_{td},...,X_{t,10}]$, and $X_{td}$ is sampled from a
uniform distribution ranging from 0 to 1. Subsequently, all the entries in the
distance matrix measuring all the pairs of the subjects are also scaled to
$[0,1]$. For each comparison, $100$ replications are generated.

To have a fair comparison, we do not use the true means
$\\{\bm{S}_{k}:k=1,2,...,K\\}$ as the input of our algorithm. Instead, we give
the trained mean as followed. For each given $\bm{S}_{k}$, we generate $10$
positive definite matrices according to our model, and then we calculate their
sample mean. The sample means, denoted as
$\\{\tilde{\bm{S}}_{k}:k=1,2,...,K\\}$, are the algorithm inputs. Using
$\\{\tilde{\bm{S}}_{k}:k=1,2,...,K\\}$ is also better mimicking the real
world, because the practitioners can never obtain the true mean
$\\{\bm{S}_{k}:k=1,2,...,K\\}$ from the real data but from the training data.

Huang et al. (2015) uses the Log-Euclidean metric for image set
classification. The Log-Euclidean metric has already been used as a metric
measuring the distance between positive definite matrices in many other
studies (e.g., Arsigny et al., 2006). In our study, we implement this method
in a simplified way, that is to give the group of the image $t$ via finding a
$k$ which minimizes the function
$||\log\bm{A}_{t}-\log\tilde{\bm{S}}_{k}||_{F}$. Since both our algorithm and
the method of Log-Euclidean metric utilize the information from the training
data, i.e., $\\{\tilde{\bm{S}}_{k}:k=1,2,...,K\\}$, the method of Log-
Euclidean metric can be a benchmark method to tell the improvement caused by
modeling the residual correlation caused by the non-voxel information.

A compelling benchmark method is the Wishart BNP proposed by Cherian et al.
(2016), which is a Dirichlet process of Wishart matrices. Similar to our
algorithm, the Wishart BNP was also proposed with the application to the RCD-
based image set classification. The Wishart BNP does not require training data
and produces classification with an undetermined number of groups.

The methods we give above utilize the whole matrix. Alternatively, we may only
use the eigenvalues of an RCD. By fixing the number of groups as $K=3$, we can
classify the images by the K-means and the Gaussian mixture model, treating
the eigenvalues as the multivariate responses. The properties of our proposed
method (EM Hybrid) and the benchmark methods are summarized in Table 1,
helping us learn the key differences among the methods.

Table 1: The proposed method (EM Hybrid) and the benchmark methods are summarized. The properties of these methods are number of groups (fixed v.s. flexible), training data request (yes v.s. no), model based (yes v.s. no), and information used (whole matrix v.s. eigenvalues).  Method |  | Number
---
of Groups
| Training Data
---
Request
| Model
---
Based
| Information
---
Used
EM Hybrid | Fixed | Yes | Yes | Whole Matrix
Wishart BNP | Flexible | No | Yes | Whole Matrix
Log-Euclidean Metric | Fixed | Yes | No | Whole Matrix
K-means | Fixed | No | No | Eigenvalues
Gaussian Mixture Model | Fixed | No | Yes | Eigenvalues

We inspect the performance of our algorithm in two aspects: (a) the
classification accuracy and (b) the role of the tuning parameters $\lambda$
and $u$ in the likelihood-weight function. To explore the roles of the tuning
parameters, we give the combinations of $\lambda\in\\{5\times 2^{-5},5\times
2^{-4},5\times 2^{-3},5\times 2^{-2},5\times 2^{-1},5\times 2^{0}\\}$, and
$u\in\\{0.2,0.4,0.6,0.8\\}$. The exponential component in the weight function
is displayed in Figure 2.

Figure 2: The exponential component in the weight function is displayed.

The methods to be compared include the supervised methods (i.e., EM Hybrid and
Log-Euclidean Metric), the semi-supervised methods (i.e., Gaussian mixture
Model and K-means), and the unsupervised method (i.e., Wishart BNP), To make
the classification accuracy measured under the same system, we use the
commonly-used Rand index as metric to describe the classification accuracy
(Hubert and Arabie, 1985). The Rand index (Rand, 1971) is a value ranging from
0 to 1, and a larger value indicates that the classification is more accurate.

First, we compare our algorithm to the other benchmark methods, in terms of
classification accuracy. The side-by-side boxplots of the Rand indices are
displayed in Figure 3. From the results, our algorithm with different settings
of tuning parameters produces a more accurate classification result than the
others. We also note that the method using only eigenvalues both produce worse
results, compared to the ones which utilize the whole matrix. The performance
of the Log-Euclidean metric method shows that the inclusion of the variations
caused by the non-voxel information is important. As we stated before, the
Wishart BNP is a compelling method. In our simulation study, the performance
of our algorithm is better than the Wishart BNP.

Figure 3: The side-by-side boxplots of the Rand indices are displayed. The
title of each subplot gives the value of $u$. Under a certain $u$, the x-axis
is for the values of $\lambda$ and the y-axis is for the Rand index. The Rand
indices produced by a method of all replications are visualized by a boxplot.
For each method, the medians of the boxplots are connected by a dashed line.
The legends on the right side give the colors for each method. To clarify, the
benchmark methods do not rely on $\lambda$ and $u$ but are only visualized at
the locations to compare with the corresponding Rand indices produced by our
algorithm.

A nice feature of our algorithm is that our algorithm is robust to the tuning
parameters, in terms of classification accuracy. However, we still recommend
to give a larger $u$ if the large computational cost allows in a certain
problem, because we observe that a larger $u$ leads to larger Rand indices in
Figure 4 where the Rand indices are aggregated over all $\lambda\in\\{5\times
2^{-5},5\times 2^{-4},5\times 2^{-3},5\times 2^{-2},5\times 2^{-1},5\times
2^{0}\\}$. Given the demo script provided in Appendix C, the typical
computational time of $u=0.1$ and $T=50$ is 1 minute.

Figure 4: The side-by-side boxplots of the Rand indices are displayed. The
x-axis is for the values of $u$ and the y-axis is for the Rand index. Given a
$u$, the rand indices are aggregated over all $\lambda\in\\{5\times
2^{-5},5\times 2^{-4},5\times 2^{-3},5\times 2^{-2},5\times 2^{-1},5\times
2^{0}\\}$. The aggregated Rand indices are visualized by a boxplot. The
medians of the boxplots are connected by a dashed line.

Furthermore, we want to investigate the role of the tuning parameters on
parameter estimation. The parameter estimations are evaluated by the mean
squared error (MSE). The MSEs are defined as follows:

  * •

Group-Weights: $\frac{1}{K}\sum_{k=1}^{K}(\hat{\omega}_{k}-\omega_{k})^{2}$

  * •

Range Parameter: $(\hat{\phi}-\phi)^{2}$

  * •

Degrees of Freedom: $(\hat{M}-M)^{2}$,

where the terms without hats are the true values of the parameters and the
terms with hats are the estimates. We use dot plots to show the relationship
between the MSE and the log-likelihood (Figure 5). Generally, we observe an
negative association between the MSE and the log-likelihood. In light of the
negative association, we recommend reporting the parameter estimation whose
log-likelihood value is largest.

Figure 5: We use scatter plot to show the relationship between the MSE and the
log-likelihood. The title of each subplot gives the value $u$. The x-axis is
for the values of MSE and the y-axis is for the values of log-likelihood. The
legends on the right side give the colors for each parameter.

##  5 Application to the Chicago Face Database

We use the Chicago face database (Ma et al., 2015) as the motivating data-set.
The data-set is open to the public. It can be requested from the website
https://chicagofaces.org. The data-set has headshots of the subjects of
different races, together with the subjects’ non-voxel information (e.g., age,
nose width, etc). Our goal is to classify the images by their races (i.e.,
Black, Latino, and White), through using their RCDs. Race classification is
one of the common scientific goals in computer vision studies (e.g., Ou et
al., 2005; Roomi et al., 2011; Fu et al., 2014). The true races of the
subjects are known in our motivating data-set, thus we can evaluate the
classification accuracy by comparing the inferred race classification and the
true race classification.

The RCDs are calculated as follows. First we only use the voxels which cover
the skin to calculate the RCD. This step can be simply implemented by the
MATLAB codes provided by Jain (2020). The feature vector is then defined as
$\bm{F}(v)=[F(v,R),F(v,G),F(v,B)]^{T}$, where $F(v,R),F(v,G),F(v,B)$ are the
magnitudes of red, green, and blue at the voxel $v$. We further normalized the
magnitudes by dividing the image’s sample standard error of the magnitudes at
the blue channel, denoted as
$\tilde{\bm{F}}(v)=\bm{F}(v)/\sigma_{B}=[F(v,R),F(v,G),F(v,B)]^{T}/\sigma_{B}$.
Finally, the RCD is computed by taking the sample covariance of all
$\tilde{\bm{F}}(v)$.

Note that the original RCD paper (Tuzel et al., 2006) only introduced a big
picture, and the scheme to generate the RCD varies under different scenarios.
In short, the scheme to generate the RCDs relies on the relevant domain
knowledge. Our paper primarily focuses on introducing the proposed
methodology. Thus we simply use the scheme described above, and the scheme is
considerably useful in our classification problem. However, other appropriate
schemes (e.g., Roomi et al., 2011; Cherian et al., 2016) to generate the RCDs
can also be used, if the features extracted are useful in the specific
classification problem. More importantly, the other schemes also have no
hurdle to be used in our algorithm.

After obtaining the RCDs of the subjects, we select the training data to
obtain the trained means $\\{\tilde{\bm{S}}_{k}:k=1,2,...,K\\}$. For each race
(i.e., Black, Latino, and White), we randomly select one-fifth of the subjects
of each race and calculate their sample mean treated as the trained mean
$\\{\tilde{\bm{S}}_{k}:k=1,2,...,K\\}$. We use all the available numerical
non-voxel information to construct $\bm{X}_{t}$. Each covariate is normalized,
and the entries in the distance matrix measuring all the pairs of the subjects
are also scaled to $[0,1]$.

We continue to use the methods listed in Table 1. We also continue to give the
combinations of $\lambda\in\\{5\times 2^{-5},5\times 2^{-4},5\times
2^{-3},5\times 2^{-2},5\times 2^{-1},5\times 2^{0}\\}$, and
$u\in\\{0.2,0.4,0.6,0.8\\}$. The Rand index is used to measure the
classification accuracy too. The classification accuracy in terms of the Rand
index is visualized in Figure 6. The real data result is consistent with what
we have had in the simulation study (Section 4). First, our algorithm still
holds the first place among all the methods. Also, the classification result
is robust to the tuning parameters in the likelihood-weight function. In the
real data analysis, the model-based methods (EM Hyrbid, Wishart BNP, and
Gaussian Mixture Model) are better than the non-model-based methods (Log-
Euclidean Metric and K-means). This may imply the importance of including
statistical uncertainties in modeling real data. The Gaussian mixture model
performs relatively well. This may be caused by that the eigenvalues provide
distinguishing information in this classification problem.

Figure 6: The side-by-side dot-plots of the Rand indices are displayed. The
title of each subplot gives the value of $u$. Under a certain $u$, the x-axis
is for the values of $\lambda$ and the y-axis is for the Rand index. The Rand
index produced by a method is visualized by a dot. For each method, the dots
are connected by a dashed line. The legends on the right side give the colors
for each method. To clarify, the benchmark methods do not rely on $\lambda$
and $u$ but are only visualized at the locations to compare with the
corresponding Rand index produced by our algorithm.

In light of the suggestion provided in Section 4, the parameter estimation are
given as follows. The estimation of the weights are $\hat{\omega}_{1}=0.3333$,
$\hat{\omega}_{2}=0.3334$, and $\hat{\omega}_{3}=0.3333$; the parameter
estimation for the range parameter is $\hat{\phi}=0.5647$; the parameter
estimation for the degrees of freedom is $\hat{M}=39.1730$. The associated
tuning parameters are $u=0.4$ and $\lambda=5\times 2^{-5}$, whose log-
likelihood which is the largest among all settings.

##  6 Conclusion and Future Directions

In this paper, we propose a novel EM composite likelihood-based algorithm of
Wishart matrices to classify the Wishart matrices which are correlated caused
by the non-voxel information. The methodological details are explicitly given
in the paper. The algorithm also performs better than the other methods in
terms of image set classification. In summary, our paper has significant
contributions in both the methodological development of matrix-variate
analysis and the application of computer vision.

The proposed hierarchical model is set up in a relatively simplified way. More
complicated features can be added given specific computer vision problems. For
example, the mean of the positive definite matrices may also be determined by
the subject-level covariates, but not necessarily be plug-in values. However,
the added complicated features shall make the statistical inference more
challenging. Statistical inference tools such as variational Bayesian methods
and Markov chain Monte Carlo methods are encouraged to be implemented in the
future.

Acknowledgments

We would like to acknowledge support for this project from research funds from
Center of Outcome Research and Evaluation, Yale University. We thank Professor
Brian J Reich at North Carolina State University for providing valuable
suggestions. We finally thank the Center for Decision Research, The University
of Chicago for providing the motivating data, allowing methodological
development in computer vision.

##  A Density Function of Parameterized Wishart Distribution

Let $\bm{W}\sim\mathcal{W}_{p}(\bm{\Sigma},M)$. The probability density
function of $\bm{W}$ is

${\displaystyle
f_{\mathbf{W}}(\mathbf{w})={\frac{|\mathbf{w}|^{(M-p-1)/2}e^{-\operatorname{tr}(M\mathbf{\Sigma}^{-1}\mathbf{w})/2}}{2^{\frac{Mp}{2}}|{\mathbf{\Sigma}/M}|^{M/2}\Gamma_{p}({\frac{M}{2}})}}}.$
(15)

##  B Derivations of the Bivariate Density

This scheme to obtain the joint density borrows the idea from (Smith and
Garth, 2007) in which the joint density for complex Wishart matrices is
derived. First, we give the joint density function of
$[\bm{U}_{t},\bm{U}_{s}]$ specified in our model (Model 4), denoted as
$f_{\bm{U}_{t},\bm{U}_{s}}[\bm{u}_{t},\bm{u}_{s}]$. The function
$f_{\bm{U}_{t},\bm{U}_{s}}(\bm{u}_{t},\bm{u}_{s})$ can be decomposed as

$f_{\bm{U}_{t},\bm{U}_{s}}(\bm{u}_{t},\bm{u}_{s})=f_{\bm{U}_{t}|\bm{U}_{s}}(\bm{u}_{t}|\bm{u}_{s})\times
f_{\bm{U}_{s}}(\bm{u}_{s}),$ (16)

where $f_{\bm{U}_{t}|\bm{U}_{s}}(\bm{u}_{t}|\bm{u}_{s})$ is the conditional
density function of $[\bm{U}_{t}|\bm{U}_{s}]$ and $f_{\bm{U}_{s}}(\bm{u}_{s})$
is the marginal density of $\bm{U}_{s}$. Given James et al. (1964, Equation
67), we can give the explicit expression of
$f_{\bm{U}_{t}|\bm{U}_{s}}(\bm{u}_{t}|\bm{u}_{s})$, which is expressed as

$\displaystyle
f_{\bm{U}_{t}|\bm{U}_{s}}(\bm{u}_{t}|\bm{u}_{s})=\exp\left[-\frac{1}{2}tr\left(\frac{M\rho_{ts}}{1-\rho_{ts}^{2}}\bm{u}_{s}\right)\right]\times\
_{0}F_{1}\left(\frac{1}{2}M;\frac{1}{4}\left(\frac{M\rho_{ts}}{1-\rho_{ts}^{2}}\right)^{2}\bm{u}_{s}\bm{u}_{t}\right)$
(17)
$\displaystyle\times\frac{1}{2^{Mp/2}\times\Gamma_{p}(\frac{1}{2}M)\times|\frac{1-\rho_{ts}^{2}}{M}\bm{I}|^{M/2}}\times\exp\left[-\frac{1}{2}tr\left(\frac{M}{1-\rho_{ts}^{2}}\bm{u}_{s}\right)\right]\times|\bm{u}_{t}|^{\frac{1}{2}(M-p-1)},$

where is ${}_{0}F_{1}(;)$ a hypergeometric function of matrix argument. Then
the joint density function $f_{\bm{U}_{t},\bm{U}_{s}}(\bm{u}_{t},\bm{u}_{s})$
is

$\displaystyle
f_{\bm{U}_{t},\bm{U}_{s}}(\bm{u}_{t},\bm{u}_{s})=f_{\bm{U}_{t}|\bm{U}_{s}}(\bm{u}_{t}|\bm{u}_{s})\times
f_{\bm{U}_{s}}(\bm{u}_{s})$ (18)
$\displaystyle=\exp\left[-\frac{1}{2}tr\left(\frac{M\rho_{ts}^{2}}{1-\rho_{ts}^{2}}\bm{u}_{s}\right)\right]\times\
_{0}F_{1}\left(\frac{1}{2}M;\frac{1}{4}\left(\frac{M\rho_{ts}}{1-\rho_{ts}^{2}}\right)^{2}\bm{u}_{s}\bm{u}_{t}\right)$
$\displaystyle\times\frac{1}{2^{Mp/2}\times\Gamma_{p}(\frac{1}{2}M)\times|\frac{1-\rho_{ts}^{2}}{M}\bm{I}|^{M/2}}\times\exp\left[-\frac{1}{2}tr\left(\frac{M}{1-\rho_{ts}^{2}}\bm{u}_{t}\right)\right]\times|\bm{u}_{t}|^{\frac{1}{2}(M-p-1)}$
$\displaystyle\times|\bm{u}_{s}|^{\frac{1}{2}(M-p-1)}\times\exp\left(-\frac{1}{2}tr\left(\bm{u}_{s}M\right)\right)\times\left(2^{Mp/2}\times|\bm{I}/M|^{M/2}\times\Gamma_{p}\left(\frac{1}{2}M\right)\right)^{-1}$
$\displaystyle=\
_{0}F_{1}\left(\frac{1}{2}M;\frac{1}{4}\left(\frac{M\rho_{ts}}{1-\rho_{ts}^{2}}\right)^{2}\bm{u}_{s}\bm{u}_{t}\right)\times
2^{-Mp}\times\left[\Gamma_{p}(\frac{1}{2}M)\right]^{-2}\times|\bm{u}_{t}|^{\frac{1}{2}(M-p-1)}\times|\bm{u}_{s}|^{\frac{1}{2}(M-p-1)}$
$\displaystyle\times\left(\frac{1-\rho_{ts}^{2}}{M^{2}}\right)^{-pM/2}\times\exp\left(-\frac{1}{2}tr\left[M\frac{1}{1-\rho_{ts}^{2}}(\bm{u}_{s}+\bm{u}_{t})\right]\right).$

Next, we want to derive the expression of the joint density function
$f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi})$.
This can be done via variable change. We have the relationship that
$\bm{U}_{t}=\bm{Q}_{j_{t}}^{-1}\bm{A}_{t}(\bm{Q}_{j_{t}}^{T})^{-1}$. Be very
careful, the Jacobian is
$|\bm{J}|=|\bm{\Sigma}_{t}^{-1}|^{(p+1)/2}\times|\bm{\Sigma}_{s}^{-1}|^{(p+1)/2}$,
given Proposition 5.11 of Eaton (1983, Chapter 5: Matrix Factorizations and
Jacobians)! Therefore, the expression of
$f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi})$
is

$\displaystyle
f_{\bm{A}_{t},\bm{A}_{s}}(\bm{a}_{t},\bm{a}_{s}|\bm{\Sigma}_{t}=\bm{S}_{j_{t}},\bm{\Sigma}_{s}=\bm{S}_{j_{s}},\bm{\Phi})$
(19) $\displaystyle=\
_{0}F_{1}\left(\frac{1}{2}M;\frac{1}{4}\left(\frac{M\rho_{ts}}{1-\rho_{ts}^{2}}\right)^{2}\bm{Q}_{j_{s}}^{-1}\bm{a}_{s}(\bm{Q}_{j_{s}}^{T})^{-1}\bm{Q}_{j_{t}}^{-1}\bm{a}_{t}(\bm{Q}_{j_{t}}^{T})^{-1}\right)$
$\displaystyle\times
2^{-Mp}\times\left[\Gamma_{p}(\frac{1}{2}M)\right]^{-2}\times|\bm{S}_{j_{t}}^{-1}\bm{a}_{t}|^{\frac{1}{2}(M-p-1)}\times|\bm{S}_{j_{s}}^{-1}\bm{a}_{s}|^{\frac{1}{2}(M-p-1)}$
$\displaystyle\times\left(\frac{1-\rho_{ts}^{2}}{M^{2}}\right)^{-\frac{pM}{2}}\times\exp\left(-\frac{1}{2}tr\left[M\frac{1}{1-\rho_{ts}^{2}}\left(\bm{S}_{j_{s}}^{-1}\bm{a}_{s}+\bm{S}_{j_{t}}^{-1}\bm{a}_{t}\right)\right]\right)$
$\displaystyle\times|\bm{S}_{j_{t}}^{-1}|^{\frac{p+1}{2}}\times|\bm{S}_{j_{s}}^{-1}|^{\frac{p+1}{2}},$

##  C Codes

The MATLAB codes are in EM$\\_$Hybrid$\\_$Codes.zip. The instructions of
implementing the codes are provided readme.txt. Some example scripts are also
given.

## References

  * Arsigny et al. (2006) Vincent Arsigny, Pierre Fillard, Xavier Pennec, and Nicholas Ayache.  Log-euclidean metrics for fast and simple calculus on diffusion tensors.  _Magnetic Resonance in Medicine: An Official Journal of the International Society for Magnetic Resonance in Medicine_ , 56(2):411–421, 2006.
  * Bai et al. (2012) Yun Bai, Peter X-K Song, and TE Raghunathan.  Joint composite estimating functions in spatiotemporal models.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 74(5):799–824, 2012.
  * Bai et al. (2014) Yun Bai, Jian Kang, and Peter X-K Song.  Efficient pairwise composite likelihood estimation for spatial-clustered data.  _Biometrics_ , 70(3):661–670, 2014.
  * Bevilacqua et al. (2012) Moreno Bevilacqua, Carlo Gaetan, Jorge Mateu, and Emilio Porcu.  Estimating space and space-time covariance functions for large data sets: a weighted composite likelihood approach.  _Journal of the American Statistical Association_ , 107(497):268–280, 2012.
  * Blumenson and Miller (1963) LE Blumenson and KS Miller.  Properties of generalized rayleigh distributions.  _The Annals of Mathematical Statistics_ , 34(3):903–910, 1963.
  * Chen et al. (2016) Jiahua Chen, Yi Huang, and Peiming Wang.  Composite likelihood under hidden markov model.  _Statistica Sinica_ , pages 1569–1586, 2016.
  * Chen et al. (2020) Kai-Xuan Chen, Jie-Yi Ren, Xiao-Jun Wu, and Josef Kittler.  Covariance descriptors on a gaussian manifold and their application to image set classification.  _Pattern Recognition_ , page 107463, 2020.
  * Cherian et al. (2016) Anoop Cherian, Vassilios Morellas, and Nikos Papanikolopoulos.  Bayesian nonparametric clustering for positive definite matrices.  _IEEE Transactions on Pattern Analysis & Machine Intelligence_, (1):1–1, 2016.
  * Choi and Briollais (2011) Y Choi and L Briollais.  An em composite likelihood approach for multistage sampling of family data.  _Statistica Sinica_ , pages 231–253, 2011.
  * Cressie (1992) Noel Cressie.  Statistics for spatial data.  _Terra Nova_ , 4(5):613–617, 1992.
  * Dennis and Moré (1977) John E Dennis, Jr and Jorge J Moré.  Quasi-newton methods, motivation and theory.  _SIAM review_ , 19(1):46–89, 1977.
  * Díaz-Pernas et al. (2009) FJ Díaz-Pernas, Miriam Antón-Rodríguez, JF Díez-Higuera, Mario Martínez-Zarzuela, D González-Ortega, and Daniel Boto-Giralda.  Texture classification of the entire brodatz database through an orientational-invariant neural architecture.  In _International Work-Conference on the Interplay Between Natural and Artificial Computation_ , pages 294–303. Springer, 2009.
  * Dryden et al. (2009) Ian L Dryden, Alexey Koloydenko, and Diwei Zhou.  Non-Euclidean statistics for covariance matrices, with applications to diffusion tensor imaging.  _The Annals of Applied Statistics_ , pages 1102–1123, 2009.
  * Eaton (1983) Morris L Eaton.  Multivariate statistics: a vector space approach.  _JOHN WILEY & SONS, INC., 605 THIRD AVE., NEW YORK, NY 10158, USA, 1983, 512_, 1983.
  * Fu et al. (2014) Siyao Fu, Haibo He, and Zeng-Guang Hou.  Learning race from face: A survey.  _IEEE transactions on pattern analysis and machine intelligence_ , 36(12):2483–2509, 2014.
  * Gao and Song (2011) Xin Gao and Peter X-K Song.  Composite likelihood em algorithm with applications to multivariate hidden markov model.  _Statistica Sinica_ , pages 165–185, 2011.
  * Gelfand et al. (2004) Alan E Gelfand, Alexandra M Schmidt, Sudipto Banerjee, and CF Sirmans.  Nonstationary multivariate process modeling through spatially varying coregionalization.  _Test_ , 13(2):263–312, 2004.
  * Heyde (2008) Christopher C Heyde.  _Quasi-likelihood and its application: a general approach to optimal parameter estimation_.  Springer Science & Business Media, 2008.
  * Huang et al. (2015) Zhiwu Huang, Ruiping Wang, Shiguang Shan, Xianqiu Li, and Xilin Chen.  Log-euclidean metric learning on symmetric positive definite manifold with application to image set classification.  In _International conference on machine learning_ , pages 720–729, 2015.
  * Hubert and Arabie (1985) Lawrence Hubert and Phipps Arabie.  Comparing partitions.  _Journal of classification_ , 2(1):193–218, 1985\.
  * Jain (2020) Gaurav Jain, 2020.  URL https://www.mathworks.com/matlabcentral/fileexchange/28565-skin-detection.  Accessed Dec. 10, 2020.
  * James et al. (1964) Alan T James et al.  Distributions of matrix variates and latent roots derived from normal samples.  _The Annals of Mathematical Statistics_ , 35(2):475–501, 1964.
  * Koev and Edelman (2006) Plamen Koev and Alan Edelman.  The efficient evaluation of the hypergeometric function of a matrix argument.  _Mathematics of Computation_ , 75(254):833–846, 2006.
  * Lan et al. (2020+) Zhou Lan, Brian J Reich, and Dipankar Bandyopadhyay.  A spatial bayesian semiparametric mixture model for positive definite matrices with applications to diffusion tensor imaging.  _Canadian Journal of Statistics_ , Accepted, 2020+.
  * Lan et al. (2021+) Zhou Lan, Brian J Reich, Joseph Guinness, Dipankar Bandyopadhyay, Liangsuo Ma, and F Gerard Moeller.  The spatial Wishart process and its application to diffusion tensor imaging.  _Biometrics_ , Minor Revision, 2021+.
  * Lee and Schwartzman (2017) Han Na Lee and Armin Schwartzman.  Inference for eigenvalues and eigenvectors in exponential families of random symmetric matrices.  _Journal of Multivariate Analysis_ , 162:152–171, 2017\.
  * Lindsay (1988) Bruce G Lindsay.  Composite likelihood methods.  _Contemporary mathematics_ , 80(1):221–239, 1988\.
  * Ma et al. (2015) Debbie S Ma, Joshua Correll, and Bernd Wittenbrink.  The chicago face database: A free stimulus set of faces and norming data.  _Behavior research methods_ , 47(4):1122–1135, 2015.
  * Ou et al. (2005) Yongsheng Ou, Xinyu Wu, Huihuan Qian, and Yangsheng Xu.  A real time race classification system.  In _2005 IEEE International Conference on Information Acquisition_ , pages 6–pp. IEEE, 2005.
  * Padoan et al. (2010) Simone A Padoan, Mathieu Ribatet, and Scott A Sisson.  Likelihood-based inference for max-stable processes.  _Journal of the American Statistical Association_ , 105(489):263–277, 2010.
  * Pang et al. (2008) Yanwei Pang, Yuan Yuan, and Xuelong Li.  Gabor-based region covariance matrices for face recognition.  _IEEE Transactions on circuits and systems for video technology_ , 18(7):989–993, 2008.
  * Rand (1971) William M Rand.  Objective criteria for the evaluation of clustering methods.  _Journal of the American Statistical Association_ , 66(336):846–850, 1971.
  * Roomi et al. (2011) S Md Mansoor Roomi, SL Virasundarii, S Selvamegala, S Jeevanandham, and D Hariharasudhan.  Race classification based on facial features.  In _2011 third national conference on computer vision, pattern recognition, image processing and graphics_ , pages 54–57. IEEE, 2011.
  * Smith and Garth (2007) Peter J Smith and Lee M Garth.  Distribution and characteristic functions for correlated complex Wishart matrices.  _Journal of Multivariate Analysis_ , 98(4):661–677, 2007.
  * Tuzel et al. (2006) Oncel Tuzel, Fatih Porikli, and Peter Meer.  Region covariance: A fast descriptor for detection and classification.  In _European conference on computer vision_ , pages 589–600. Springer, 2006.
  * Varin et al. (2011) Cristiano Varin, Nancy Reid, and David Firth.  An overview of composite likelihood methods.  _Statistica Sinica_ , pages 5–42, 2011.
  * Vasdekis et al. (2012) Vassilis GS Vasdekis, Silvia Cagnone, and Irini Moustaki.  A composite likelihood inference in latent variable models for ordinal longitudinal responses.  _Psychometrika_ , 77(3):425–441, 2012.
  * Viraswami (1991) Kalyanee Viraswami.  On multivariate gamma distributions.  Master’s thesis, McGill University, 1991.
  * Williams and Rasmussen (2006) Christopher KI Williams and Carl Edward Rasmussen.  _Gaussian processes for machine learning_ , volume 2.  MIT press Cambridge, MA, 2006.


# Fine-grained anomaly detection via multi-task self-supervision

Loïc Jézéquel 1,2 Ngoc-Son Vu 1 Jean Beaudet 2 Aymeric Histace 1

1 ETIS UMR 8051 (CY Cergy Paris Université, ENSEA, CNRS) F-95000
2 Idemia Identity & Security, 95520 Osny France
{loic.jezequel, son.vu<EMAIL_ADDRESS>
###### Abstract

Detecting anomalies using deep learning has become a major challenge over the
last years, and is becoming increasingly promising in several fields. The
introduction of self-supervised learning has greatly helped many methods
including anomaly detection where simple geometric transformation recognition
tasks are used. However these methods do not perform well on fine-grained
problems since they lack finer features. By combining both high-scale shape
features and low-scale fine features in a multi-task framework, our method
greatly improves fine-grained anomaly detection. It outperforms state-of-the-
art with up to 31% relative error reduction measured with AUROC on various
anomaly detection problems including one-vs-all, out-of-distribution detection
and face presentation attack detection.

978-1-6654-3396-9/21/$31.00 ©2021 IEEE

Detecting anomalies straying apart from a well-defined normal situation has
always been a major challenge in many fields such as video surveillance [33,
43], intrusion detection [15], fraud detection [42], medical imaging [20] and
more recently adversarial attack detection [25]. Deep visual anomaly detection
has been introduced to tackle this problem and has proven to be more robust
and reliable than classical binary classification. Rather than directly try to
discriminate anomalies from normal samples, we only learn the normal class
boundary and deem as anomalous any observation outside.

Recently, the introduction of self-supervised learning has greatly improved
many one-class anomaly detection learning methods. It enables to discriminate
anomalies from normal samples by learning to solve simple tasks such as
geometric transformation classification [10]. However, even if this approach
has greatly improved anomaly detection performance, it still suffers from
limitations on more challenging problems with local and fine-grained
differences between anomalies and normal samples.

In this given context, our main contributions in this paper are the following:

  * •

We improve the detection of fine-grained anomalies by independently solving in
a multi-task self-supervised fashion high-scale geometric task and low-scale
jigsaw puzzle task.

  * •

We validate the efficiency of the proposed method using an exhaustive protocol
for anomaly detection on one-vs-all, out-of-distribution detection and anti-
spoofing problems.

  * •

The proposed method obtains better overall results with up to 31% AUROC
relative improvement from state of the art methods.

##  1 Related work

###  1.1 Anomaly detection

The main goal in anomaly detection is to classify a sample as normal or
anomalous. Formally, we predict $P({\mathbf{x}}\in\mathcal{X}_{\text{norm}})$
for an observation ${\mathbf{x}}$ and a normal (or positive) class
$\mathcal{X}_{\text{norm}}$. In practice, a proxy anomaly score function
$s_{a}({\mathbf{x}})$ is usually estimated instead. The anomalous (or
negative) class is then defined implicitly as the complementary of the normal
class in image space. We can generally categorize anomalies into three
families:

  1. 1.

Object anomaly: any object which is not included in the positive class, e.g.,
a cat is an object anomaly in regards to dogs.

  2. 2.

Style anomaly: observations representing the same object as the positive class
but with a different style or support, e.g., a realistic mask or a printed
face represent faces but with a visible different style.

  3. 3.

Local anomaly: observations representing and sharing the same style as the
positive class, however a localized part of the image is different. Most of
the time, these anomalies are the superposition of two generative processes,
e.g., a fake nose on a real face is a local anomaly.

Usually, we assume in anomaly detection that only normal samples are available
during training, meaning that most methods are part of one-class learning
scheme. The first introduced methods simply used a pre-trained neural network
to extract features, on which a classical algorithm such as One-Class SVM [31]
(OCSVM) or Isolation Forest [17] (IF) were trained.

There have been also semi-supervised anomaly detection methods such as DeepSAD
[27] or deviation networks [23] where we assume some of the anomalies
representing a few modes are available. These methods can achieve better
accuracy on borderline cases given enough diverse anomalies, which is often
less manageable in practice. In particular, these two methods directly learn
representations by minimizing the distance of normal sample features to an
hypersphere center, while maximizing the distance to the anomalies. It follows
the compactness principle, where we minimize the normal class representations
variance and maximize the inter-class representations variance.

###  1.2 Self-supervised learning

Self Supervised Learning (SSL) is a part of representation learning, where we
want to learn useful and general representations from an unlabeled dataset
$\mathcal{X}=\\{{\mathbf{x}}_{i}\\}_{1}^{N}$. We can then use the learned
features for a different task such as classification.

We learn representations by solving from the data an auxiliary task
$\mathcal{T}$, which is often unrelated to the final one. Therefore SSL
consists of two steps:

  1. 1.

Generating a labeled dataset $\mathcal{X_{\mathcal{T}}}$ aligned with
$\mathcal{T}$, which for classification is usually done by applying $c$
transformations $T_{j}$ to our unlabeled samples

$\mathcal{X}_{\mathcal{T}}=\\{(T_{j}({\mathbf{x}}_{i}),j)\\}_{i,j}$ (1)

  2. 2.

Training a classification or regression network on this generated labeled set.

One of the final layers $\phi_{\mathcal{T}}$ can thus be used as a feature
extractor. Some commonly used tasks are: 90° rotation prediction [9], jigsaw
puzzle [22], distortions [7], colorization [41], image inpainting [24] or
relative patches prediction [6].

###  1.3 SSL anomaly detection

Very recently, SSL has been adapted to the one-class anomaly detection
framework. First we learn to solve an auxiliary task $\mathcal{T}$ in a SSL
fashion to obtain a pre-trained network $\phi_{\mathcal{T}}$. Then, to
classify at inference time an observation ${\mathbf{x}}$ as anomalous or
normal, we evaluate how well the network can solve the task. Indeed, the main
assumption is that the network will perform relatively well on normal samples
but will fail on anomalies. A task-independent metric $L$ is computed on the
generated labeled samples to compute the anomaly score function:

$s_{a}({\mathbf{x}})=\left\\{L(\phi_{\mathcal{T}}(T_{i}({\mathbf{x}})),i)|i\in\llbracket
1,c\rrbracket\right\\}$ (2)

Unlike SSL, we are not directly interested in the intermediate features, but
rather the final task outputs.

In GeoTrans [10], the auxiliary task is to classify which geometrical
transformation has been applied to the input. A set of 72 transformations
including identity is randomly sampled over all possible compositions of
translations, rotations and symmetries. At the end of training, 72 Dirichlet
distributions respectively parameterized by $\tilde{\boldsymbol{\alpha}}_{i}$
are fitted over the normal class softmax outputs
$\mathbf{y}\left(T_{i}({\mathbf{x}})\right)$ for each transformation. The log-
likelihood can then be used during inference as the task-independent metric
$L$:

$s_{a}({\mathbf{x}})=\sum_{i=1}^{72}\left(\tilde{\boldsymbol{\alpha}}_{i}-1\right)\cdot\log\mathbf{y}\left(T_{i}({\mathbf{x}})\right)$
(3)

In MHRot [12], the task is to simultaneously classify three types of
transformations, each modeled by a softmax head: vertical translation,
horizontal translations and 90° rotations. Accordingly, we are trying to
predict the three following variables: vertical translations
(0,$-t_{y}$,$+t_{y}$), horizontal translations (0,$-t_{x}$,$+t_{x}$) and 90°
rotations (0°,90°,180°,270°).

During inference, we sum the three softmax of the known transformations for
each transformation combination:

$s_{a}(x)=\quad\sum\
\sum_{{\tiny\mathclap{\begin{subarray}{c}r\in\\{0,90,180,270\\}\\\
s\in\\{0,-t_{x},+t_{x}\\}\\\ t\in\\{0,-t_{y},+t_{y}\\}\end{subarray}}}}\
\sum\quad\mathbf{y}(T_{\scriptscriptstyle r,s,t}(x))_{r,s,t}$ (4)

##  2 Method overview

###  2.1 Anomaly detection pretext task

We present here a general rule of thumb regarding the choice of tasks for SSL
anomaly detection. It is generally more restrictive than for simple
representation learning [2].

Let $\mathcal{T}$ be a task along its training loss $L_{\mathcal{T}}$. On the
one hand, if the task is too hard on normal samples, meaning that the accuracy
of our network remains close to random predictor throughout training (or that
$\|\nabla_{\phi}L_{\mathcal{T}}\|$ is always small and that the minimum of $L$
is high), then no meaningful representation will be reached at convergence.
This will also result in poor accuracy on anomalies (Fig.1.c) and yield
unpredictable results during anomaly detection. On the other hand, if the task
is too easy on normal samples, meaning that our model will converge to a
perfect predictor in the first epochs, then the task loss will be minimized by
many representations including trivial ones. Thus the network is more likely
to learn such representations which will be unspecific to the normal class and
encode very generic visual features. Since many anomalies will share these
features, the task accuracy will be high on anomalies as well (Fig.1.a).

To observe these effects, we train a network on several isolated tasks as
described in Section 1.3. By monitoring its task classification accuracy on
evaluation normal data and anomalous data during the first epochs, we
empirically measure how well-suited a pretext task is for anomaly detection on
a given dataset. We show that even though 90° rotation is more adapted than
translations on coarse anomaly detection, it ultimately fails on fine-grained
anomaly detection (Fig.1.b). This confirms that basic geometric transformation
recognition tasks, such as 90° rotations or translations, are only suited to
simple object anomaly detection. Indeed, since these tasks are solvable
accurately by learning high scale and shape features, it is unlikely the
network will use finer characteristics that allow discriminating normal
samples from more subtle anomalies.

Figure 1: Tasks accuracy during training on coarse object AD (CIFAR-10 in
plain line) and fine-grained AD (CaltechBirds in dotted line): (a) vertical
translations, (b) 90° rotation, (c) unsolvable task

###  2.2 Method overview

Finding a single task satisfying all the previous conditions is difficult, and
must be highly dataset dependent. Therefore we resort to ensemble methods [5]
by allowing the network to learn $N$ tasks and merge their decision at
inference. We learn richer features via multi-task learning [3], by sharing a
common representation across all tasks. Our model is accordingly composed of a
main feature extractor network $\phi$ and $N$ dense layers
$f_{\mathcal{T}_{1}},\cdots,f_{\mathcal{T}_{N}}$, where
$f_{\mathcal{T}_{i}}(\phi({\mathbf{x}}))$ is the output for the ith task.

Figure 2: Multi-task self-supervised anomaly detection. In dotted line are
additional steps during inference

During inference, we aggregate the anomaly scores of all tasks into the final
anomaly detection score. The whole training and inference scheme is summarized
in Figure 2. For each classification task $\mathcal{T}_{i}$, the task-
independent metric $L_{\mathcal{T}_{i}}$ chosen is the softmax score
corresponding to the true known class and we sum up these scores using the
mean:

$s_{a}({\mathbf{x}})=\frac{1}{N}\sum_{i=1}^{N}s_{a}^{(\mathcal{T}_{i})}({\mathbf{x}})$
(5)

where $s_{a}^{(\mathcal{T}_{i})}$ is the anomaly score of the ith task:

$s_{a}^{(\mathcal{T}_{i})}({\mathbf{x}})=\sum_{j}{\operatorname{softmax}}(\phi\circ
f_{\mathcal{T}_{i}}(T_{j}^{(i)}({\mathbf{x}})))_{j}$ (6)

We note that there is a caveat using the mean as anomaly score: adding new
tasks can have a negative impact on the model performance. In practice if the
task is not well suited to the normal class, it will add significant noise to
the anomaly score and ultimately harm the anomaly detection accuracy.

To prevent our multi-task from being too easy on fine-grained problems, we
introduce more challenging tasks. We choose here a simplified version of the
jigsaw puzzle task. The jigsaw puzzle task consists in splitting an image into
a grid of $n_{h}\times n_{w}$ patches, then randomly shuffling the different
patches. The task is then to predict the original order of each patch. This
task has proven in representation learning to provide a great challenge for
extracting more local and finer features [22]. To avoid trivial solutions and
force our model to understand pieces neighborhood, we are careful to add a
margin between each patch with a random small offset.

Permutation 1

Permutation 2

Permutation 3

Figure 3: Example of simplified jigsaw puzzle task for $k=3$

Since we chose the softmax truth as the task-independent metric, we need to
re-frame it into a classification problem by considering each permutation as a
single class. This would greatly increase our model complexity, effectively
adding $(n_{w}\cdot n_{h})!$ classes. Therefore, we only consider
$k<(n_{w}\cdot n_{h})!$ randomly chosen permutations including the identity
permutation. This quantity $k$ becomes an additional parameter controlling the
task difficulty.

The complete training loss for a single sample ${\mathbf{x}}$ becomes

$\displaystyle\textstyle
L({\mathbf{x}})=\textstyle\sum_{i=1}^{3}L_{CE}(\phi\circ
f_{\text{v}}(T_{i}^{(\text{v})}({\mathbf{x}})),i)\,+$
$\displaystyle\textstyle\sum_{i=1}^{3}L_{CE}(\phi\circ
f_{\text{h}}(T_{i}^{(\text{h})}({\mathbf{x}})),i)\,+$
$\displaystyle\textstyle\sum_{i=1}^{4}L_{CE}(\phi\circ
f_{\text{rot}}(T_{i}^{(\text{rot})}({\mathbf{x}})),i)\,+$
$\displaystyle\textstyle\sum_{i=1}^{k}L_{CE}(\phi\circ
f_{\text{puzz}}(T_{i}^{(\text{puzz})}({\mathbf{x}})),i)$ (7)

where $L_{CE}$ is the cross-entropy and
$f_{\text{v}},f_{\text{h}},f_{\text{rot}},f_{\text{puzz}}$ are respectively
the dense layers for the vertical translations, horizontal translations,
rotations and puzzle tasks. Compared to the MHRot anomaly score in Equation 4
we evaluate each task independently, which allows us to greatly reduce the
required amount of network forward pass to compute the anomaly score during
inference. As a consequence, our method inference step is roughly 10 times
faster.

By combining the base geometrical transformation recognition task with the
jigsaw task, we allow the model to learn high-scale shape features more suited
toward object anomaly detection as well as low-scale fine features more suited
toward style anomaly and local anomaly detection.

##  3 Implementation details

The geometrical transformation task is composed as in [12] of horizontal
translations, vertical translations and 90° rotations. As for the jigsaw
puzzle task, we found best results with $n_{w}=n_{h}=3$ and $k=3$.

Regarding network architecture, we use a 16-4 WideResNet [40] ($\approx 10M$
parameters with a depth of 16) for the feature extractor network $\phi$, along
with two dense softmax layers respectively of size $10$ for the geometrical
transformation task and size $3$ for the jigsaw puzzle task. Each of these
dense layers have a dropout rate of 0.3 [32]. Training is performed under SGD
optimizer with nesterov momentum [34], using a batch size of 32.

##  4 Results

###  4.1 Evaluation protocol

Until now, most of the anomaly detection literature have adopted the one-vs-
all protocol to evaluate their method. In the one-vs-all protocol, we consider
one class of a multi-class dataset, originally created for object recognition,
as the normal class. All the other classes are then considered as anomalous,
and we can in a leave-one-out cross-validation fashion evaluate the model on
each possible normal class. The final reported result is the mean of each run.

Even though such datasets are easier to acquire and result in a highly multi-
modal anomaly class, these might not be enough to fully evaluate anomaly
detection methods. Indeed, these only cover coarse object anomalies which are
now becoming too easy for state-of-the-art methods, and do not reflect
realistic anomaly detection challenges.

| Dataset | Anomaly type
---|---|---
| Object | Style | Local
{31.7cm[  Obj.classif ] | MNIST | ✓ | ✗ | ✗
F-MNIST | ✓ | ✗ | ✗
CIFAR-100 | ✓ | ✗ | ✗
{21.7cm[  Fine-grained ] | Caltech-Birds | ✓ | ✓ | ✗
FounderType | ✗ | ✓ | ✗
{11.7cm[  Anti-spoofing ] | SiW-M | ✓ | ✓ | ✓
Table 1: Summary of evaluation datasets. | Model | CIFAR-100 | MNIST | F-MNIST | Caltech-Birds 200 | Fonts | SiW-M
---|---|---|---|---|---|---|---
{31.9cm[  Semi-Supervised ] | Deep-SAD (75%) [27] | 88.7 | 99.9 | 98.1 | 73.6 | 99.8 | 85.4
Deep-SAD (25%) | 87.9 | 98.5 | 95.4 | 70.9 | 99.4 | 76.0
Deep-SAD (10%) | 89.1 | 96.5 | 88.2 | 66.1 | 98.0 | 80.6
{91.9cm[  One-class ] | ADGAN [4] | 54.7 | 94.7 | 88.4 | - | - | -
GANomaly [1] | 56.5 | 92.8 | 80.9 | - | - | -
ARNet [8] | 78.8 | 98.3 | 93.9 | - | - | -
2-8 | OCSVM [31] | - | 84.7 | 74.2 | 76.3 | - | -
| IF [17] | - | 87.1 | 84.0 | 74.2 | - | -
| PIAD [36] | 78.8 | 98.1 | 94.3 | 63.5 | - | 81.2
| GeoTrans [10] | 84.7 | 96.9 | 92.6 | 66.6 | 92.3 | 81.1
| MHRot [12] | 83.6 | 95.2 | 92.5 | 77.6 | 96.7 | 83.1
2-8 | Ours | 85.8 | 96.0 | 92.8 | 83.2 | 96.9 | 88.4
Table 2: Comparison with the state-of-the-art AUROC over several datasets,
underline indicates best result, bold indicates best one-class learning
result. We re-implemented all the methods except the three one-class methods
in the first block (results are from original papers [4, 8]).
---
Model | Airplane | Automobile | Bird | Cat | Deer | Dog | Frog | Horse | Ship | Truck | Avg
VAE [13] | 70.0 | 38.6 | 67.9 | 53.5 | 74.8 | 52.3 | 68.7 | 49.3 | 69.6 | 38.6 | 58.3
OCSVM [31] | 63.0 | 44.0 | 64.9 | 48.7 | 73.5 | 50.0 | 72.5 | 53.3 | 64.9 | 50.8 | 58.5
AnoGAN [30] | 67.1 | 54.7 | 52.9 | 54.5 | 65.1 | 60.3 | 58.5 | 62.5 | 75.8 | 66.5 | 61.8
PixelCNN [37] | 53.1 | 99.5 | 47.6 | 51.7 | 73.9 | 54.2 | 59.2 | 78.9 | 34.0 | 66.2 | 61.8
Deep-SVDD [28] | 61.7 | 65.9 | 50.8 | 59.1 | 60.9 | 65.7 | 67.7 | 67.3 | 75.9 | 73.1 | 64.8
OCGAN [26] | 75.7 | 53.1 | 64.0 | 62.0 | 72.3 | 62.0 | 72.3 | 57.5 | 82.0 | 55.4 | 65.6
Puzzle-AE [29] | 78.9 | 78.0 | 69.9 | 54.8 | 75.4 | 66.0 | 74.7 | 73.3 | 83.3 | 69.9 | 72.4
DROCC [11] | 81.7 | 76.7 | 66.7 | 67.1 | 73.6 | 74.4 | 74.4 | 71.4 | 80.0 | 76.2 | 74.2
GeoTrans [10] | 74.7 | 95.7 | 78.1 | 72.4 | 87.8 | 87.8 | 83.4 | 95.5 | 93.3 | 91.3 | 86.0
1-12 Ours | 75.1 | 96.3 | 84.8 | 74.2 | 91.1 | 89.9 | 88.7 | 95.5 | 94.7 | 91.9 | 88.2
Table 3: Detailed comparison with one-class state-of-the-art AUROC on CIFAR-10
dataset.

Thus we propose to use fine-grained classification datasets in the same one-
vs-all protocol. Since discrimination between theses classes is mostly based
on local and fine patterns, we can have a good coverage of style anomalies and
local anomalies. Also we note that because of the increased shift in object
recognition toward fined-grain classification, such datasets have become
readily available. For one-vs-all datasets, we used MNIST [16], Fashion MNIST
[39], CIFAR-100 [14]. For the fine-grained dataset, we chose the Caltech-Birds
200 database [38].

We also put forward datasets from real anomaly detection problems over
different fields. First, we use font recognition challenges as they provide
shape-focused style anomaly detection. Indeed two different fonts represent
the same characters albeit with a distinctive style. Even though these images
lie on a low dimensional manifold compared to natural images, they still
provide insight into how well the model can capture small shape hints. In
particular, we use FounderType-200 [18] introduced for novelty detection and
containing 6700 images per font. Furthermore, we choose a dataset from face
anti-spoofing, where the goal is to discriminate real faces from fake
representations of someone’s face. Due to the richness and high variability of
such frauds, this problem effectively encompasses all three types of
anomalies. We use here the Spoof in the Wild Multiple (SiW-M) [19] database
which contains more than 1600 short videos of real faces and presentation
attacks. There are 493 real identities along with several types of attacks:
paper print, screen replay, masks and partial attacks where only a localized
area of the face is fake. The masks are composed of half-masks, paper masks,
silicone mask and transparent masks. All evaluation datasets are summarized in
Table 1.

We additionally evaluate our anomaly detection model on out-of-distribution
(OOD) protocol. OOD detection, which is broader than anomaly detection, aims
at discriminating the training dataset from other data distributions. The
”normal” distribution in OOD is therefore usually more diverse and highly
multi-modal. We also have a greater overlap in term of class between the in-
distribution samples and out-of-distribution samples compared to anomaly
detection. Nevertheless it gives us great insight into the multi-modality
limits of our model. The most common evaluation setup is to discriminate one
training multi-class dataset from other datasets. Here we choose to learn on
CIFAR-10 and discriminate CIFAR-100 and the easier Street View House Numbers
(SVHN) dataset [21].

For all of the evaluations, the metric used is the area under the ROC curve
(AUROC), averaged over all possible normal classes in the case of one-vs-all
datasets.

###  4.2 Ablation study

We evaluate in Table 4 how combining the two tasks of geometric transformation
recognition and jigsaw puzzle improves the anomaly detection. We drastically
improve performances with a relative error reduction regarding AUROC of 13% on
CIFAR-100, 25% on Caltech-Birds 200 and 31% on SIW-M. This validates our
statement in Section 2: the finer the differences between anomaly and normal
class, the greater the improvement is by adding the jigsaw task.

---
Auxiliary Task | CIFAR-100 | Caltech-Birds | SiW-M
Geometric (G) | 83.6 | 77.6 | 83.1
Jigsaw (J) | 80.1 | 78.5 | 76.3
Ours (G+J) | 85.8 | 83.2 | 88.4
Table 4: AUROC for different tasks, best result is in bold.

###  4.3 Comparison to the state-of-the-art

We compare our method with different one-class learning state-of-the-art
approaches to anomaly detection: reconstruction error generative models with
the PIAD model [36], self-supervised methods with GeoTrans [10] and MHRot
[12]. As an addition, we include a semi-supervised learning anomaly detection
method DeepSAD [27], which has access to a portion of the anomalies during
training. As such, we train it with the same normal samples but three
different ratio of the anomaly subclasses: 10%, 25% and 75%.

For the sake of fair comparison in the same conditions, we take the existing
implementations or re-implement each method and evaluate each, except for the
ADGAN, GANomaly and ARNet which we reference results from their original
papers [4, 1, 8].

The results are gathered in Table 2 and 3. First of all, we can see our method
generally maintains among the best accuracies on simple object anomaly
detection, and even improves it on more challenging datasets such as
CIFAR-100. Moreover, it greatly improves fine-grained anomaly detection and
outperforms state-of-the-art methods which could not be realistically be used
for this problem. We also show that our method, without further tuning,
improves anti-spoofing detection performances on SiW-M. Finally, we notice our
one-class learning model generally reduces the gap with semi-supervised
method, and even outperforms these on Caltech-Birds 200 and SiW-M, even though
these take advantage of a significant amount of additional anomalous data.

---
Metrics | AUROC | EER |  | APCER
---
(5%BPCER)
MHRot [12] | 83.0 | 21.6 | 77.5
Ours | 88.4 | 18.7 | 39.1
Table 5: AUROC, EER and APCER at 5% BPCER of MHRot and our model with jigsaw
task on SiW-M dataset, best result is in bold.

We compare in Table 5 our method with the second best self-supervised method
MHRot on SiW-M. We use metrics more adapted to face presentation attack
detection with equal error rate (EER) and the false acceptance rate for the
rate of false reject fixed at 5% (APCER@5%BPCER). Our comparison does not
include other face anti-spoofing methods since we only use real faces training
images while these all use a set of presentation attacks during training.
Using our method, the APCER@5%BPCER drops from 77.5 to 39.1 thus also showing
promising usage of anomaly detection methods in fraud detection.

---
OOD | SVHN | CIFAR-100 | Avg
VAE [13] | 2.4 | 52.8 | 27.6
Deep-SVDD [28] | 14.5 | 52.1 | 33.3
PixelCNN [37] | 15.8 | 52.4 | 34.1
RotNet [9] | 97.9 | 81.2 | 89.5
Ours | 98.8 | 83.4 | 91.1
1-4 CSI [35] | 99.8 | 89.2 | 94.5
Table 6: Comparison with state-of-the-art on the Out-Of-Distribution detection
protocol with CIFAR-10 as in-distribution, best result is underlined, best
pretext task driven method is in bold.

Lastly, we compare our model with state-of-the-art on OOD detection in Table
6. Although not designed specifically for such complex normal class, we obtain
better detection rates than other self-supervised anomaly detection methods
with pretext tasks.

##  5 Conclusion and Future Work

In this paper, we investigate the power of multi-task self supervision for
anomaly detection and show the limits of simple geometric tasks. In more
details, we combine two complementary tasks of jigsaw puzzle and geometric
transformation recognition. Through an ablation study, we show that this
enables it to learn much complex and finer features and therefore better
detect anomalies. Finally, we provide a more comprehensive evaluation protocol
than previously used datasets in the anomaly detection literature. It presents
more challenging datasets and covers object, style and local anomalies. Our
method outperforms state-of-the-art, including a semi-supervised method, on
most of the fine-grained datasets.

For future work we could explore the combination of more tasks, including
generative tasks (in contrast to discriminative tasks used here). Such tasks
could range from re-colorization to image in-painting.

## References

  * [1] S. Akcay, A. Atapour-Abarghouei, and T. P. Breckon.  Ganomaly: Semi-supervised anomaly detection via adversarial training.  In Asian Conference on Computer Vision, pages 622–637. Springer, 2018.
  * [2] Y. M. Asano, C. Rupprecht, and A. Vedaldi.  A critical analysis of self-supervision, or what we can learn from a single image.  In International Conference on Learning Representations, 2020.
  * [3] R. Caruana.  Multitask learning.  Machine Learning, 28(1):41–75, Jul 1997.
  * [4] L. Deecke, R. Vandermeulen, L. Ruff, S. Mandt, and M. Kloft.  Image Anomaly Detection with Generative Adversarial Networks.  In M. Berlingerio, F. Bonchi, T. Gärtner, N. Hurley, and G. Ifrim, editors, Machine Learning and Knowledge Discovery in Databases, Lecture Notes in Computer Science, pages 3–17. Springer International Publishing, 2019.
  * [5] T. G. Dietterich.  Ensemble methods in machine learning.  In Multiple Classifier Systems, pages 1–15, Berlin, Heidelberg, 2000. Springer Berlin Heidelberg.
  * [6] C. Doersch, H. Mulam, and A. Efros.  Unsupervised visual representation learning by context prediction.  IEEE International Conference on Computer Vision (ICCV), pages 1422–1430, 05 2015.
  * [7] A. Dosovitskiy, J. Springenberg, M. Riedmiller, and T. Brox.  Discriminative unsupervised feature learning with exemplar convolutional neural networks.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 1, 06 2014.
  * [8] Y. Fei, C. Huang, C. Jinkun, M. Li, Y. Zhang, and C. Lu.  Attribute Restoration Framework for Anomaly Detection.  IEEE Transactions on Multimedia, pages 1–1, 2020.
  * [9] S. Gidaris, P. Singh, and N. Komodakis.  Unsupervised representation learning by predicting image rotations.  In 6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings, 2018.
  * [10] I. Golan and R. El-Yaniv.  Deep anomaly detection using geometric transformations.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems 31, pages 9758–9769. Curran Associates, Inc., 2018.
  * [11] S. Goyal, A. Raghunathan, M. Jain, H. V. Simhadri, and P. Jain.  Drocc: Deep robust one-class classification.  In International Conference on Machine Learning, 2020.
  * [12] D. Hendrycks, M. Mazeika, S. Kadavath, and D. Song.  Using self-supervised learning can improve model robustness and uncertainty.  Advances in Neural Information Processing Systems (NeurIPS), 2019\.
  * [13] D. Kingma and M. Welling.  Auto-encoding variational bayes.  In Conference proceedings: papers accepted to the International Conference on Learning Representations (ICLR), 2014.
  * [14] A. Krizhevsky.  Learning multiple layers of features from tiny images.  2009\.
  * [15] D. Kwon, H. Kim, J. Kim, S. Suh, I. Kim, and K. Kim.  A survey of deep learning-based network anomaly detection.  Cluster Computing, 22, 01 2019.
  * [16] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [17] F. T. Liu, K. Ting, and Z.-H. Zhou.  Isolation forest.  In 2008 Eighth IEEE International Conference on Data Mining, pages 413 – 422, 01 2009.
  * [18] J. Liu, Z. Lian, Y. Wang, and J. Xiao.  Incremental kernel null space discriminant analysis for novelty detection.  In 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 4123–4131, 07 2017.
  * [19] Y. Liu, J. Stehouwer, A. Jourabloo, and X. Liu.  Deep tree learning for zero-shot face anti-spoofing.  In In Proceeding of IEEE Computer Vision and Pattern Recognition, Long Beach, CA, June 2019.
  * [20] A. S. Lundervold and A. Lundervold.  An overview of deep learning in medical imaging focusing on mri.  Zeitschrift für Medizinische Physik, 29(2):102 – 127, 2019.  Special Issue: Deep Learning in Medical Physics.
  * [21] Y. Netzer, T. Wang, A. Coates, A. Bissacco, B. Wu, and A. Ng.  Reading digits in natural images with unsupervised feature learning.  In NIPS Workshop on Deep Learning and Unsupervised Feature Learning, 2011.
  * [22] M. Noroozi and P. Favaro.  Unsupervised learning of visual representations by solving jigsaw puzzles.  In B. Leibe, J. Matas, N. Sebe, and M. Welling, editors, Computer Vision – ECCV 2016, pages 69–84, Cham, 2016. Springer International Publishing.
  * [23] G. Pang, C. Shen, and A. V. D. Hengel.  Deep anomaly detection with deviation networks.  Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, 2019.
  * [24] D. Pathak, P. Krähenbühl, J. Donahue, T. Darrell, and A. A. Efros.  Context encoders: Feature learning by inpainting.  2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 2536–2544, 2016.
  * [25] A. Paudice, L. Muñoz-González, A. György, and E. C. Lupu.  Detection of adversarial training examples in poisoning attacks through anomaly detection.  CoRR, abs/1802.03041, 2018.
  * [26] P. Perera, R. Nallapati, and B. Xiang.  Ocgan: One-class novelty detection using gans with constrained latent representations.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 2898–2906, 2019.
  * [27] L. Ruff, R. A. Vandermeulen, N. Görnitz, A. Binder, E. Müller, K.-R. Müller, and M. Kloft.  Deep semi-supervised anomaly detection.  In International Conference on Learning Representations, 2020.
  * [28] L. Ruff, R. A. Vandermeulen, N. Görnitz, L. Deecke, S. A. Siddiqui, A. Binder, E. Müller, and M. Kloft.  Deep one-class classification.  In Proceedings of the 35th International Conference on Machine Learning, volume 80, pages 4393–4402, 2018.
  * [29] M. Salehi, A. Eftekhar, N. Sadjadi, M. H. Rohban, and H. R. Rabiee.  Puzzle-ae: Novelty detection in images through solving puzzles.  abs/2008.12959.
  * [30] T. Schlegl, P. Seeböck, S. M. Waldstein, U. Schmidt-Erfurth, and G. Langs.  Unsupervised anomaly detection with generative adversarial networks to guide marker discovery.  In International conference on information processing in medical imaging, pages 146–157. Springer, 2017.
  * [31] B. Schölkopf, R. Williamson, A. Smola, J. Shawe-Taylor, and J. Platt.  Support vector method for novelty detection.  In Proceedings of the 12th International Conference on Neural Information Processing Systems, NIPS’99, page 582–588, Cambridge, MA, USA, 1999\. MIT Press.
  * [32] N. Srivastava, G. E. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov.  Dropout: a simple way to prevent neural networks from overfitting.  J. Mach. Learn. Res., 15:1929–1958, 2014.
  * [33] W. Sultani, C. Chen, and M. Shah.  Real-World Anomaly Detection in Surveillance Videos.  In 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 6479–6488, 06 2018.  ISSN: 2575-7075.
  * [34] I. Sutskever, J. Martens, G. Dahl, and G. Hinton.  On the importance of initialization and momentum in deep learning.  In S. Dasgupta and D. McAllester, editors, Proceedings of the 30th International Conference on Machine Learning, volume 28 of Proceedings of Machine Learning Research, pages 1139–1147, Atlanta, Georgia, USA, 17–19 Jun 2013. PMLR.
  * [35] J. Tack, S. Mo, J. Jeong, and J. Shin.  CSI: novelty detection via contrastive learning on distributionally shifted instances.  Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, virtual, 2020.
  * [36] N. Tuluptceva, B. Bakker, I. Fedulova, and A. Konushin.  Perceptual Image Anomaly Detection.  In S. Palaiahnakote, G. Sanniti di Baja, L. Wang, and W. Q. Yan, editors, Pattern Recognition, Lecture Notes in Computer Science, pages 164–178. Springer International Publishing, 2020.
  * [37] A. Van den Oord, N. Kalchbrenner, L. Espeholt, O. Vinyals, A. Graves, et al.  Conditional image generation with pixelcnn decoders.  In Advances in neural information processing systems, pages 4790–4798, 2016.
  * [38] P. Welinder, S. Branson, T. Mita, C. Wah, F. Schroff, S. Belongie, and P. Perona.  Caltech-UCSD Birds 200.  Technical Report CNS-TR-2010-001, California Institute of Technology, 2010\.
  * [39] H. Xiao, K. Rasul, and R. Vollgraf.  Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms.  CoRR, abs/1708.07747, 2017.
  * [40] S. Zagoruyko and N. Komodakis.  Wide residual networks.  In R. C. Wilson, E. R. Hancock, and W. A. P. Smith, editors, Proceedings of the British Machine Vision Conference 2016, BMVC 2016, York, UK, September 19-22, 2016. BMVA Press, 2016.
  * [41] R. Zhang, P. Isola, and A. A. Efros.  Colorful image colorization.  In B. Leibe, J. Matas, N. Sebe, and M. Welling, editors, Computer Vision - ECCV 2016 - 14th European Conference, Amsterdam, The Netherlands, October 11-14, 2016, Proceedings, Part III, volume 9907 of Lecture Notes in Computer Science, pages 649–666. Springer, 2016.
  * [42] Z. Zhang, X. Zhou, X. Zhang, L. Wang, and P. Wang.  A model based on convolutional neural network for online transaction fraud detection.  Security and Communication Networks, 2018:1–9, 08 2018.
  * [43] S. Zhu, C. Chen, and W. Sultani.  Video anomaly detection for smart surveillance.  CoRR, abs/2004.00222, 2020.