text:
\section{Introduction}

Symmetric positive definite (SPD) matrices have been used in many different contexts. In diffusion tensor imaging for instance, a diffusion tensor is a 3-dimensional SPD matrix \cite{lenglet_statistics_2006,pennec_riemannian_2006,Fletcher07}; in brain-computer interfaces (BCI) \cite{Barachant13}, in functional MRI \cite{Deligianni11} or in computer vision \cite{Cheng13}, an SPD matrix can represent a covariance matrix of a feature vector, for example a spatial covariance of electrodes or a temporal covariance of signals in BCI. In order to make statistical operations on SPD matrices like interpolations, computing the mean
or performing a principal component analysis, it has been proposed to consider the set of SPD matrices as a manifold and to provide it with some geometric structures like a \textit{Riemannian metric}, a transitive \textit{group action} or some \textit{symmetries}. These structures can be more or less natural depending on the context of the applications, and they can provide closed-form formulas and consistent algorithms \cite{pennec_riemannian_2006,Dryden09}.

Many Riemannian structures have been introduced over the manifold of SPD matrices \cite{Dryden09}: Euclidean, log-Euclidean, affine-invariant, Cholesky, square root, power-Euclidean, Procrustes... Each of them has different mathematical properties that can fit the data in some problems but can be inappropriate in some other contexts: for example the curvature can be null, positive, negative, constant, not constant, covariantly constant... These properties on the curvature have some important consequences on the way we interpolate two points, on the consistence of algorithms, and more generally on every statistical operation one could want to do with SPD matrices. Therefore, a natural question one can ask is: given the practical context of an application, how should one choose the metric on SPD matrices? Are there some relations between the mathematical properties of the geometric structure and the intrinsic properties of the data?

In this context, the affine-invariant metric \cite{pennec_riemannian_2006,Lenglet06,Fletcher07} was introduced to give an invariant computing framework under affine transformations of the variables. This metric endows the manifold of SPD matrices with a structure of a Riemannian symmetric space. Such spaces have a covariantly constant curvature, thus they share some convenient properties with constant curvature spaces but with less constraints. It was actually shown that there exists not only one but a one-parameter family that is invariant under these affine transformations \cite{Pennec08}. More recently, \cite{Su12,Su14,Zhang18} introduced another Riemannian symmetric structure that does not belong to the previous one-parameter family: the polar-affine metric.

In this work, we unify these two frameworks by showing that the polar-affine metric is a square deformation of the affine-invariant metric (Section 2). We generalize in Section 3.1 this construction to a family of power-affine metrics that comprises the two previous metrics, and in Section 3.2 to the wider family of deformed-affine metrics. Finally, we propose in Section 4 a theoretical approach in the choice of subfamilies of the deformed-affine metrics with relevant properties.

\enlargethispage{8mm}

\section{Affine-invariant versus polar-affine}

The affine-invariant metric \cite{pennec_riemannian_2006,Lenglet06,Fletcher07} and the polar-affine metric \cite{Zhang18} are different but they both provide a Riemannian symmetric structure to the manifold of SPD matrices. Moreover, both claim to be very naturally introduced. The former uses only the action of the real general linear group ${GL}_n$ on covariance matrices. The latter uses the canonical left action of ${GL}_n$ on the left coset space ${GL}_n/O_n$ and the polar decomposition ${GL}_n\simeq{SPD}_n\times O_n$, where $O_n$ is the orthogonal group. Furthermore, the affine-invariant framework is exhaustive in the sense that it provides \textit{all} the metrics invariant under the chosen action \cite{Pennec08} whereas the polar-affine framework only provides \textit{one} invariant metric.

In this work, we show that the two frameworks coincide on the same quotient manifold ${GL}_n/O_n$ but differ because of the choice of the diffeomorphism between this quotient and the manifold of SPD matrices. In particular, we show that there exists a one-parameter family of polar-affine metrics and that any polar-affine metric is a square deformation of an affine-invariant metric.

In 2.1 and 2.2, we build the affine-invariant metrics $g^1$ and the polar-affine metric $g^2$ in a unified way, using indexes $i\in\{1,2\}$ to differentiate them. First, we give explicitly the action $\eta^i:{GL}_n\times{SPD}_n\longrightarrow{SPD}_n$ and the quotient diffeomorphism $\tau^i:{GL}_n/O_n\longrightarrow{SPD}_n$; then, we explain the construction of the orthogonal-invariant scalar product $g^i_{I_n}$ that characterizes the metric $g^i$; finally, we give the expression of the metrics $g^1$ and $g^2$. In 2.3, we summarize the results and we focus on the Riemannian symmetric structures of ${SPD}_n$.

\subsection{The one-parameter family of affine-invariant metrics}

\subsubsection{Affine action and quotient diffeomorphism}
In many applications, one would like the analysis of covariance matrices to be invariant under affine transformations $X\longmapsto AX+B$ of the random vector $X\in\mathbb{R}^n$, where $A\in{GL}_n$ and $B\in\mathbb{R}^n$. Then the covariance matrix $\Sigma=\mathrm{Cov}(X)$, is modified under the transformation $\Sigma\longmapsto A\Sigma A^\top$. This transformation can be thought as a transitive Lie group action $\eta^1$ of the general linear group on the manifold of SPD matrices:
\begin{equation}
    \eta^1:\fun{{GL}_n\times{SPD}_n}{{SPD}_n}{(A,\Sigma)}{\eta^1_A(\Sigma)=A\Sigma A^\top}.
\end{equation}
This transitive action induces a diffeomorphism between the manifold ${SPD}_n$ and the quotient of the acting group ${GL}_n$ by the stabilizing group $\mathrm{Stab}^1(\Sigma)=\{A\in{GL}_n,\eta^1(A,\Sigma)=\Sigma\}$ at any point $\Sigma$. It reduces to the orthogonal group at $\Sigma=I_n$ so we get the quotient diffeomorphism $\tau^1$:
\begin{equation}
    \tau^1:\fun{{GL}_n/O_n}{{SPD}_n}{[A]=A.O_n}{\eta^1(A,I_n)=AA^\top}.
\end{equation}

\subsubsection{Orthogonal-invariant scalar product}
We want to endow the manifold $\mathcal{M}={SPD}_n$ with a metric $g^1$ invariant under the affine action $\eta^1$, i.e. an affine-invariant metric. As the action is transitive, the metric at any point $\Sigma$ is characterized by the metric at one given point $I_n$. As the metric is affine-invariant, this scalar product $g_{I_n}$ has to be invariant under the stabilizing group of $I_n$. As a consequence, the metric $g^1$ is characterized by a scalar product $g^1_{I_n}$ on the tangent space $T_{I_n}\mathcal{M}$ that is invariant under the action of the orthogonal group.

The tangent space $T_{I_n}\mathcal{M}$ is canonically identified with the vector space $\mathrm{Sym}_n$ of symmetric matrices by the differential of the canonical embedding $\mathcal{M}\hookrightarrow\mathrm{Sym}_n$. Thus we are now looking for all the scalar products on symmetric matrices that are invariant under the orthogonal group. Such scalar products are given by the following formula \cite{Pennec08}, where $\alpha>0$ and $\beta>-\frac{\alpha}{n}$: for all tangent vectors $V_1,W_1\in T_{I_n}\mathcal{M}$, $g^1_{I_n}(V_1,W_1)=\alpha\,\mathrm{tr}(V_1W_1)+\beta\,\mathrm{tr}(V_1)\mathrm{tr}(W_1)$.

\subsubsection{Affine-invariant metrics}
To give the expression of the metric, we need a linear isomorphism between the tangent space $T_\Sigma\mathcal{M}$ at any point $\Sigma$ and the tangent space $T_{I_n}\mathcal{M}$. Since the action $\eta^1_{\Sigma^{-1/2}}$ sends $\Sigma$ to $I_n$, its differential given by $T_\Sigma\eta^1_{\Sigma^{-1/2}}:V\in T_\Sigma\mathcal{M}\longmapsto V_1=\Sigma^{-1/2}V\Sigma^{-1/2}\in T_{I_n}\mathcal{M}$ is such a linear isomorphism.
Combining this transformation with the expression of the metric at $I_n$ and reordering the terms in the trace, we get the general expression of the affine-invariant metric: for all tangent vectors $V,W\in T_\Sigma\mathcal{M}$,
\begin{equation}
    g^1_\Sigma(V,W)=\alpha\,\mathrm{tr}(\Sigma^{-1}V\Sigma^{-1}W)+\beta\,\mathrm{tr}(\Sigma^{-1}V)\mathrm{tr}(\Sigma^{-1}W).
\end{equation}

As the geometry of the manifold is not much affected by a scalar multiplication of the metric, we often drop the parameter $\alpha$, as if it were equal to 1, and we consider that this is a one-parameter family indexed by $\beta>-\frac{1}{n}$.

\subsection{The polar-affine metric}

\subsubsection{Quotient diffeomorphism and affine action}
In \cite{Zhang18}, instead of defining a metric directly on the manifold of SPD matrices, a metric is defined on the left coset space ${GL}_n/O_n=\{[A]=A.O_n,A\in{GL}_n\}$, on which the general linear group ${GL}_n$ naturally acts by the left action $\eta^0:(A,[A'])\longmapsto[AA']$. Then this metric is pushed forward on the manifold ${SPD}_n$ into the polar-affine metric $g^2$ thanks to the polar decomposition $\mathrm{pol}:A\in{GL}_n\longmapsto(\sqrt{AA^\top},{\sqrt{AA^\top}}^{-1}A)\in{SPD}_n\times O_n$ or more precisely by the quotient diffeomorphism $\tau^2$:
\begin{equation}
    \tau^2:\fun{{GL}_n/O_n}{{SPD}_n}{A.O_n}{\sqrt{AA^\top}}.
\end{equation}
This quotient diffeomorphism induces an action of the general linear group ${GL}_n$ on the manifold ${SPD}_n$, under which the polar-affine metric will be invariant:
\begin{equation}
    \eta^2:\fun{{GL}_n\times{SPD}_n}{{SPD}_n}{(A,\Sigma)}{\eta^2_A(\Sigma)=\sqrt{A\Sigma^2A^\top}}.
\end{equation}
It is characterized by $\eta^2(A,\tau^2(A'.O_n))=\tau^2(\eta^0(A,A'.O_n))$ for $A,A'\in{GL}_n$.

\subsubsection{Orthogonal-invariant scalar product}
The polar-affine metric $g^2$ is characterized by the scalar product $g^2_{I_n}$ on the tangent space $T_{I_n}\mathcal{M}$. This scalar product is obtained by pushforward of a scalar product $g^0_{[I_n]}$ on the tangent space $T_{[I_n]}({GL}_n/O_n)$. It is itself induced by the Frobenius scalar product on $\mathfrak{gl}_n=T_{I_n}{GL}_n$, defined by $\dotprod{v}{w}_\mathrm{Frob}=\mathrm{tr}(vw^\top)$, which is orthogonal-invariant. This is summarized on the following diagram.
$$
\begin{array}{ccccc}
{GL}_n & \overset{s}{\longrightarrow} & {GL}_n/O_n & \overset{\tau^2}{\longrightarrow} & \mathcal{M}={SPD}_n\\ \relax
A & \longmapsto & A.O_n & \longmapsto & \sqrt{AA^\top}\\ \relax
\dotprod{\cdot}{\cdot}_\mathrm{Frob} & & g^0_{[I_n]} & & g^2_{I_n}
\end{array}
$$
Finally, we get the scalar product $g^2_{I_n}(V_2,W_2)=\mathrm{tr}(V_2W_2)$ for $V_2,W_2\in T_{I_n}\mathcal{M}$.

\subsubsection{Polar-affine metric}
Since the action $\eta^2_{\Sigma^{-1}}$ sends $\Sigma$ to $I_n$, a linear isomorphism between tangent spaces is given by the differential of the action
$T_\Sigma\eta^2_{\Sigma^{-1}}:V\in T_\Sigma M\longrightarrow V_2=\Sigma^{-1}T_\Sigma\mathrm{pow}_2(V)\Sigma^{-1}\in T_{I_n}\mathcal{M}$. Combined with the above expression of the scalar product at $I_n$, we get the following expression for the polar affine metric: for all tangent vectors $V,W\in T_\Sigma\mathcal{M}$,
\begin{equation}
    g^2_\Sigma(V,W)=\mathrm{tr}(\Sigma^{-2}\,T_\Sigma\mathrm{pow}_2(V)\,\Sigma^{-2}\,T_\Sigma\mathrm{pow}_2(W)).
\end{equation}

\subsection{The underlying Riemannian symmetric manifold}

In the affine-invariant framework, we started from defining the affine action $\eta^1$ (on covariance matrices) and we inferred the quotient diffeomorphism $\tau^1:({GL}_n/O_n,\eta^0)\longrightarrow({SPD}_n,\eta^1)$. In the polar-affine framework, we started from defining the quotient diffeomorphism $\tau^2:{GL}_n/O_n\longrightarrow{SPD}_n$ (corresponding to the polar decomposition) and we inferred the affine action $\eta^2$. The two actually correspond to the same underlying affine action $\eta^0$ on the quotient ${GL}_n/O_n$. Then there is also a one-parameter family of affine-invariant metrics on the quotient ${GL}_n/O_n$ and a one-parameter family of polar-affine metrics on the manifold ${SPD}_n$. This is stated in the following theorems.

\begin{theorem}[Polar-affine is a square deformation of affine-invariant]
\begin{enumerate}
    \item There exists a one-parameter family of affine-invariant metrics on the quotient ${GL}_n/O_n$
    \item This family is in bijection with the one-parameter family of affine-invariant metrics on the manifold of SPD matrices thanks to the diffeomorphism $\tau^1:A.O_n\longmapsto AA^\top$. The corresponding action is $\eta^1:(A,\Sigma)\longmapsto A\Sigma A^\top$.
    \item This family is also in bijection with a one-parameter family of polar-affine metrics on the manifold of SPD matrices thanks to the diffeomorphism $\tau^2:A.O_n\longmapsto\sqrt{AA^\top}$. The corresponding action is $\eta^2:(A,\Sigma)\longmapsto\sqrt{A\Sigma^2A^\top}$.

    \item The diffeomorphism $\mathrm{pow}_2:\fun{({SPD}_n,4g^2)}{({SPD}_n,g^1)}{\Sigma}{\Sigma^2}$ is an isometry between polar-affine metrics $g^2$ and affine-invariant metrics $g^1$.
\end{enumerate}
\end{theorem}

In other words, performing statistical analyses (e.g. a principal component analysis) with the polar-affine metric on covariance matrices is equivalent to performing these statistical analyses with the classical affine-invariant metric on the \textit{square} of our covariance matrix dataset.\\

All the metrics mentioned in Theorem 1 endow their respective space with a structure of a Riemannian symmetric manifold. We recall the definition of that geometric structure and we give the formal statement.

\begin{definition}[Symmetric manifold, Riemannian symmetric manifold]
A manifold $\mathcal{M}$ is symmetric if it is endowed with a family of involutions $(s_x)_{x\in\mathcal{M}}$ called symmetries such that $s_x\circ s_y\circ s_x=s_{s_x(y)}$ and $x$ is an isolated fixed point of $s_x$. It implies that $T_xs_x=-\mathrm{Id}_{T_x\mathcal{M}}$. A Riemannian manifold $(\mathcal{M},g)$ is symmetric if it is endowed with a family of symmetries that are isometries of $\mathcal{M}$, i.e. that preserve the metric: $g_{s_x(y)}(T_ys_x(v),T_ys_x(w))=g_y(v,w)$ for $v,w\in T_y\mathcal{M}$.
\end{definition}

\begin{theorem}[Riemannian symmetric structure on ${SPD}_n$]
The Riemannian manifold $({SPD}_n,g^1)$, where $g^1$ is an affine-invariant metric, is a Riemannian symmetric space with symmetry $s_\Sigma:\Lambda\longmapsto\Sigma\Lambda^{-1}\Sigma$.
The Riemannian manifold $({SPD}_n,g^2)$, where $g^2$ is a polar-affine metric, is also a Riemannian symmetric space whose symmetry is $s_\Sigma:\Lambda\longmapsto\sqrt{\Sigma^2\Lambda^{-2}\Sigma^2}$.
\end{theorem}

This square deformation of affine-invariant metrics can be generalized into a power deformation to build a family of affine-invariant metrics that we call power-affine metrics. It can even be generalized into any diffeomorphic deformation of SPD matrices. We now develop these families of affine-invariant metrics.

\section{Families of affine-invariant metrics}

There is a theoretical interest in building families comprising some of the known metrics on SPD matrices to understand how one can be deformed into another. For example, power-Euclidean metrics \cite{Dryden10} comprise the Euclidean metric and tends to the log-Euclidean metric \cite{Fillard07} when the power tends to 0. We recall that the log-Euclidean metric is the pullback of the Euclidean metric on symmetric matrices by the symmetric matrix logarithm $\log:{SPD}_n\longrightarrow\mathrm{Sym}_n$. There is also a practical interest in defining families of metrics: for example, it is possible to optimize the power to better fit the data with a certain distribution \cite{Dryden10}.

First, we generalize the square deformation by deforming the affine-invariant metrics with a power function $\mathrm{pow}_\theta:\Sigma\in{SPD}_n\longmapsto\Sigma^\theta=\exp(\theta\log\Sigma)$ to define the power-affine metrics. Then we deform the affine-invariant metrics by any diffeomorphism $f:{SPD}_n\longrightarrow{SPD}_n$ to define the deformed-affine metrics.

\subsection{The two-parameter family of power-affine metrics}

We recall that $\mathcal{M}={SPD}_n$ is the manifold of SPD matrices. For a power $\theta\ne 0$, we define the $\theta$-power-affine metric $g^\theta$ as the pullback by the diffeomorphism $\mathrm{pow}_\theta:\Sigma\longmapsto\Sigma^\theta$ of the affine-invariant metric, scaled by a factor ${1}/{\theta^2}$.

Equivalently, the $\theta$-power-affine metric is the metric invariant under the $\theta$-affine action $\eta^\theta:(A,\Sigma)\longmapsto(A\Sigma^\theta A^\top)^{1/\theta}$ whose scalar product at $I_n$ coincides with the scalar product $g^1_{I_n}:(V,W)\longmapsto\alpha\,\mathrm{tr}(VW)+\beta\,\mathrm{tr}(V)\mathrm{tr}(W)$. The $\theta$-affine action induces an isomorphism $V\in T_\Sigma M\longmapsto V_\theta=\frac{1}{\theta}\Sigma^{-\theta/2}\,\partial_V\mathrm{pow}_\theta(\Sigma)\,\Sigma^{-\theta/2}\in T_{I_n}\mathcal{M}$ between tangent spaces. The $\theta$-power-affine metric is given by:
\begin{equation}
    g^\theta_\Sigma(V,W)=\alpha\,\mathrm{tr}(V_\theta W_\theta)+\beta\,\mathrm{tr}(V_\theta)\mathrm{tr}(W_\theta).
\end{equation}

Because a scaling factor is of low importance, we can set $\alpha=1$ and consider that this family is a two-parameter family indexed by $\beta>-{1}/{n}$ and $\theta\ne 0$.

We have chosen to define the metric $g^\theta$ so that the power function $\mathrm{pow}_\theta:(\mathcal{M},\theta^2g^\theta)\longrightarrow(\mathcal{M},g^1)$ is an isometry. Why this factor $\theta^2$? The first reason is for consistence with previous works: the analogous power-Euclidean metrics have been defined with that scaling \cite{Dryden10}. The second reason is for continuity: when the power tends to 0, the power-affine metric tends to the log-Euclidean metric.

\begin{theorem}[Power-affine tends to log-Euclidean for $\theta\rightarrow 0$]
Let $\Sigma\in\mathcal{M}$ and $V,W\in T_\Sigma\mathcal{M}$.
Then $\lim_{\theta\rightarrow 0}{g^\theta_\Sigma(V,W)}=g^{LE}_\Sigma(V,W)$
where the log-Euclidean metric is $g^{LE}_\Sigma(V,W)=\alpha\, \mathrm{tr}(\partial_V\!\log(\Sigma)\,\partial_W\!\log(\Sigma))+\beta\,\mathrm{tr}(\partial_V\!\log(\Sigma))\mathrm{tr}(\partial_W\!\log(\Sigma))$.
\end{theorem}

\subsection{The continuum of deformed-affine metrics}
In the following, we call  a diffeomorphism $f:{SPD}_n\longrightarrow{SPD}_n$ a deformation. We define the $f$-deformed-affine metric $g^f$ as the pullback by the diffeomorphism $f$ of the affine-invariant metric, so that $f:(\mathcal{M},g^f)\longrightarrow(\mathcal{M},g^1)$ is an isometry. (Regarding the discussion before the Theorem 3, $g^{\mathrm{pow}_\theta}=\theta^2g^\theta$.)

The $f$-deformed-affine metric is invariant under the $f$-affine action $\eta^f:(A,\Sigma)\longmapsto f^{-1}(Af(\Sigma)A^\top)$. It is given by $g^f_\Sigma(V,W)=\alpha\mathrm{tr}(V_fW_f)+\beta\mathrm{tr}(V_f)\mathrm{tr}(W_f)$
where $V_f=f(\Sigma)^{-1/2}\, \partial_Vf(\Sigma) \, f(\Sigma)^{-1/2}$. The basic Riemannian operations are obtained by pulling back the affine-invariant operations.

\begin{theorem}[Basic Riemannian operations]
For SPD matrices $\Sigma,\Lambda\in\mathcal{M}$ and a tangent vector $V\in T_\Sigma\mathcal{M}$, we have at all time $t\in\mathbb{R}$:
\begin{center}
$\begin{array}{|c|l|}
\hline
\mathrm{Geodesics} & \gamma^f_{(\Sigma,V)}(t)=f^{-1}(f(\Sigma)^{1/2}\exp(tf(\Sigma)^{-1/2}T_\Sigma f(V)f(\Sigma)^{-1/2})f(\Sigma)^{1/2}) \\
\hline
\mathrm{Logarithm} & \mathrm{Log}^f_\Sigma(\Lambda)=(T_\Sigma f)^{-1}(f(\Sigma)^{1/2}\log(f(\Sigma)^{-1/2}f(\Lambda)f(\Sigma)^{-1/2})f(\Sigma)^{1/2})\\
\hline
\mathrm{Distance} & d_f(\Sigma,\Lambda)=d_1(f(\Sigma),f(\Lambda))=\sum_{k=1}^n{(\log\lambda_k)^2}\\
\hline
\end{array}$
\end{center}
where $\lambda_1,...,\lambda_n$ are the eigenvalues of the symmetric matrix $f(\Sigma)^{-1/2}f(\Lambda)f(\Sigma)^{-1/2}$.
\end{theorem}

All tensors are modified thanks to the pushforward $f_*$ and pullback $f^*$ operators, e.g. the Riemann tensor of the $f$-deformed metric is $R^f(X,Y)Z=f^*(R(f_*X,f_*Y)(f_*Z))$. As a consequence, the deformation $f$ does not affect the values taken by the sectional curvature and these metrics are negatively curved.

From a computational point of view, it is very interesting to notice that  the identification $\mathcal{L}'_\Sigma:V\in T_\Sigma\mathcal{M}\longmapsto V'=T_\Sigma f(V)\in T_{f(\Sigma)}\mathcal{M}$ simplifies the above expressions by removing the differential $T_\Sigma f$. This change of basis can prevent from numerical approximations of the differential but one must keep in mind that $V\ne V'$ in general. This identification was already used for the polar-affine metric ($f=\mathrm{pow}_2$) in \cite{Zhang18} without explicitly mentioning.

\section{Interesting subfamilies of deformed-affine metrics}

Some deformations have already been used in applications. For example, the family $A_r:\mathrm{diag}(\lambda_1,\lambda_2,\lambda_3)\longmapsto\mathrm{diag}(a_1(r)\lambda_1,a_2(r)\lambda_2,a_3(r)\lambda_3)$ where $\lambda_1\geqslant\lambda_2\geqslant\lambda_3>0$ was proposed to map the anisotropy of water measured by diffusion tensors to the one of the diffusion of tumor cells in tumor growth modeling \cite{Jbabdi05}. The inverse function $\mathrm{inv}=\mathrm{pow}_{-1}:\Sigma\longmapsto\Sigma^{-1}$ or the adjugate function $\mathrm{adj}:\Sigma\longmapsto\det(\Sigma)\Sigma^{-1}$ were also proposed in the context of DTI \cite{Lenglet04,Fuster16}. Let us find some properties satisfied by some of these examples. We define the following subsets of the set $\mathcal{F}=\mathrm{Diff}({SPD}_n)$ of diffeomorphisms of ${SPD}_n$.

(Spectral) $\mathcal{S}=\{f\in\mathcal{F}|\forall U\in O_n,\forall D\in\mathrm{Diag}^{++}_n,f(UDU^\top)=Uf(D)U^\top\}$.
Spectral deformations are characterized by their values on sorted diagonal matrices so the deformations described above are spectral: $A_r,\mathrm{adj},\mathrm{pow}_\theta\in\mathcal{S}$.\\
For a spectral deformation $f\in\mathcal{S}$, $f(\mathbb{R}_+^*I_n)=\mathbb{R}_+^*I_n$ so we can unically define a smooth diffeomorphism $f_0:\mathbb{R}_+^*\longrightarrow\mathbb{R}_+^*$ by $f(\lambda I_n)=f_0(\lambda)I_n$.

(Univariate) $\mathcal{U}=\{f\in\mathcal{S}|f(\mathrm{diag}(\lambda_1,...,\lambda_n))=\mathrm{diag}(f_0(\lambda_1),...,f_0(\lambda_n))\}$.
The power functions are univariate. Any polynomial $P=\lambda X\prod_{k=1}^p(X-a_i)$ null at 0, with non-positive roots $a_i\leqslant 0$ and positive coefficient $\lambda>0$, also gives rise to a univariate deformation.

(Diagonally-stable) $\mathcal{D}=\{f\in\mathcal{F}|f(\mathrm{Diag}^{++}_n)\subset\mathrm{Diag}^{++}_n\}$.
The deformations described above $A_r,\mathrm{adj},\mathrm{pow}_\theta$ and the univariate deformations are clearly diagonally-stable: $A_r,\mathrm{adj},\mathrm{pow}_\theta\in\mathcal{D}$ and $\mathcal{U}\subset\mathcal{D}\cap\mathcal{S}$.

(Log-linear) $\mathcal{L}=\{f\in\mathcal{F}|\log_*f=\log\circ\, f\circ\exp\,\,\text{is linear}\}$.
The adjugate function and the power functions are log-linear deformations. More generally, the functions $f_{\lambda,\mu}:\Sigma\longmapsto(\det\Sigma)^{\frac{\lambda-\mu}{n}}\Sigma^\mu$ for $\lambda,\mu\ne0$, are log-linear deformations. We can notice that the $f_{\lambda,\mu}$-deformed-affine metric belongs to the one-parameter family of $\mu$-power-affine metrics with $\beta=\frac{\lambda^2-\mu^2}{n\mu^2}>-\frac{1}{n}$.

The deformations $f_{\lambda,\mu}$ just introduced are also spectral and the following result states that they are the only spectral log-linear deformations.

\begin{theorem}[Characterization of the power-affine metrics]
If $f\in\mathcal{S}\cap\mathcal{L}$ is a spectral log-linear diffeomorphism, then there exist real numbers $\lambda,\mu\in\mathbb{R}^*$ such that $f=f_{\lambda,\mu}$ and the $f$-deformed-affine metric is a $\mu$-power-affine metric.
\end{theorem}

The interest of this theorem comes from the fact that the group of spectral deformations and the vector space of log-linear deformations have large dimensions while their intersection is reduced to a two-parameter family. This strong result is a consequence of the theory of Lie group representations because the combination of the spectral property and the linearity makes $\log_*f$ a homomorphism of $O_n$-modules (see the sketch of proof below).\\

\textit{Sketch of the proof.}
Thanks to Lie group representation theory, the linear map $F=\log_*f:{Sym}_n\longrightarrow{Sym}_n$ appears as a homomorphism of $O_n$-modules for the representation $\rho:P\in O_n\longmapsto(V\longmapsto PVP^\top)\in GL({Sym}_n)$. Once shown that ${Sym}_n=\mathrm{span}(I_n)\oplus\ker{\mathrm{tr}}$ is a $\rho$-irreducible decomposition of ${Sym}_n$ and that each one is stable by $F$, then according to Schur's lemma, $F$ is homothetic on each subspace, i.e. there exist $\lambda,\mu\in\mathbb{R}^*$ such that for $V\in\mathrm{Sym}_n$, $F(V)=\lambda\frac{\mathrm{tr}(V)}{n}I_n+\mu\left(V-\frac{\mathrm{tr}(V)}{n}I_n\right)=\log_*f_{\lambda,\mu}(V)$, so $f=f_{\lambda,\mu}$.\\

\section{Conclusion}
We have shown that the polar-affine metric is a square deformation of the affine-invariant metric and this process can be generalized to any power function or any diffeomorphism on SPD matrices. It results that the invariance principle of symmetry is not sufficient to distinguish all these metrics, so we should find other principles to limit the scope of acceptable metrics in statistical computing. We have proposed a few characteristics (spectral, diagonally-stable, univariate, log-linear) that include some functions on tensors previously introduced.
Future work will focus on studying the effect of such deformations on real data and on extending this family of metrics to positive semi-definite matrices. Finding families that comprise two non-cousin metrics could also help understand the differences between them and bring principles to make choices in applications.

\paragraph{Acknowledgements.} This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant G-Statistics agreement No 786854). This work has been supported by the French government, through the UCAJEDI Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01.

\bibliographystyle{unsrt}


text:
\section{Introduction}
\label{sec:1}
Facial landmark localization, which is localizing the facial key points (e.g., eye brows, eyes, nose, mouth and jaw), plays an important role in many computer vision tasks, such as face detection \cite{Yang_2015_ICCV}, face recognition \cite{7377089,Li2015Towards,7448432} and facial expression analysis \cite{6482635,Kamarol201725,Zhang20153191}. In recent years, facial landmark localization has been extensively studied and achieved remarkable performance on standard datasets and even on datasets collected in the wild \cite{6248014,6909614,Zhu_2015_CVPR,6751298,7780823,7911334}. However, it still has obstacles for faces with various variations in appearance including pose, expression, especially occlusions.
\begin{figure}[!tp]
	\centering
	\includegraphics[width=3.45in]{RCPRfail
	\centering
	\caption{Visual results of RCPR on COFW dataset (red: occluded, green: un-occluded). The initial shapes (the first row) and their localization results (the second row) of RCPR \cite{6751298}. Facial landmark localization usually fails when it begins with a bad initial shape.}
	\label{bad_initial}
\end{figure}
\par
Since Cascaded Pose Regression (CPR) was used to estimate facial shapes \cite{5540094}, the shape regression in a cascaded manner has emerged as one of the most popular approaches for facial landmark localization \cite{6248015,6751298,6909614,6909637,7084187,7298989,7332780,7762938,7299048,DBLP:journals/corr/SmithD16}. CPR and its variations typically begin with an initial shape, such as an average shape or a random shape of training samples, and then update the shape from coarse to fine through the trained regressors. Based on CPR, Burgos-Artizzu et al. proposed a scheme of Robust CPR (RCPR) \cite{6751298}, which is the first scheme explicitly detect occlusion state at the same time to estimate locations of landmarks. And they created a popular challenging dataset named Caltech Occluded Faces in the Wild (COFW) \cite{6751298}, where most faces in this dataset have occlusions. Researchers have used this dataset to study facial landmark localization under occlusions \cite{6751298,7084187,7762938,7780742,Yu2014,7346480}. Although these methods make some progress on facial landmark localization under partial occlusion, the occlusion problem is not essentially solved. The accuracy of occlusion prediction is still unsatisfactory. Since the occluded landmarks can hardly provide information for further analysis, it is significant to detect occlusion state of landmarks, furthermore, the occluded landmarks may reduce the accuracy of localizing the un-occluded landmarks.
\par
Since regression is initialization dependent where an improper initialization will significantly decrease the performance sharply. When the pose variation and occlusion appear simultaneously on a face, localization will fail if a bad initial shape is selected. As shown in Fig. \ref{bad_initial}, a bad initial shape usually leads to a failure of landmark localization and occlusion prediction. In this paper, we propose a Robust Initialization for CPR (RICPR)\footnote{The source code of the proposed scheme can be found at https://github.com/pervadepyy/robust-initialization-rcpr
} to avoid the bad initialization by examining texture and pose of testing face to get the texture correlated and the pose correlated initial shapes. Since texture is always related with occlusion, we select the initial shapes according to the texture correlation between the testing face and the training faces. We firstly compute Local Binary Patterns (LBP) histograms of all training faces and the testing face. Then we get correlation distance between histograms of the testing face and each training face. We choose the shapes of the most correlated training faces as the texture correlated initialization. On the other hand, the rough pose, which is represented by the rotation vector, of the testing face is used to obtain the pose correlated initial shapes. More specifically, we first estimate five fiducial landmarks including the pupils, the tip of the nose, and the corners of the mouth using Multi-Task Cascaded Convolutional Networks (MTCNN) \cite{7553523}. According to the five landmarks and a mean 3D face shape with 5 facial key points, the face pose values can be obtained. Then another 3D mean face shape, represented by 29 facial key points, can be projected to a set of corresponding 2D locations by the face pose, which can obtain the pose correlated initial shapes. Finally, the texture correlated initial shapes and the pose correlated initial shapes are taken together as the robust initialization for regression, as shown in Fig. \ref{robust_initial}, which is more relevant to the true shape of the testing face in location and occlusion. We evaluate RICPR on the challenging dataset of COFW. The experimental results show that the Normalized Mean Error (NME) is 6.64$\times10^{-2}$ and the accuracy of occlusion detection is 80/54.6\% precision/recall, which is better than that of the state-of-the-art schemes.
\par
Some of the ideas presented in this paper were initially reported in \cite{7961799}. In this paper, we report the full and new formulation and extensive experimental evaluation of our method. The initialization not only depends on texture correlation but also pose correlation and the accuracies of landmark localization and occlusion detection are further improved.
\par
The remainder of this paper are organized as follows. In Section \ref{sec:2}, the related works are briefly introduced. We review CPR and describe the proposed scheme in Section \ref{sec:3}. The experimental results on COFW dataset are given in Section \ref{sec:4}. Finally, conclusions are drawn in Section \ref{sec:5}.

\section{Related Work}
\label{sec:2}
The works to solve the problems of facial landmark localization and face alignment can be roughly divided into two groups: holistic based methods and local based methods. The holistic methods regard the shape as a whole, which usually align the face in an iterative or cascaded way. A typical holistic based method is the Active Appearance Model (AAM) \cite{Cootes1998,Matthews2004Active,Alabort-i-Medina_2014_CVPR,Antonakos2015Feature}. CPR \cite{5540094} is a similar method with a random fern regressor, which is a fast and accurate solution of computing the 2D shape of an object. Explicit Shape Regression (ESR) \cite{6248015} and RCPR \cite{6751298} extended the idea of CPR, which also use pixel difference features and fern regressor. A similar method called Supervised Descent Method (SDM) was proposed in \cite{Xiong_2013_CVPR}. This method used cascade regression with fast SIFT feature and solved localization using newton-type optimization on nonlinear least squares problem.
\par
Since occlusions are very common in the real applications of computer vision and the occluded landmarks usually cannot provide information, some works focus on facial landmark localization and occlusion detection jointly \cite{6751298,7084187,Ghiasi_2014_CVPR,7762938,Yu2014}. Burgos-Artizzu et al. first proposed to detect occlusions state at the same time of estimating landmarks in RCPR \cite{6751298}, where occlusion states are applied at each iteration to get visually different regressors. The outputs of regressors are merged with weights that depend on the occlusion prediction results. Considering occlusions often cover a region, instead of visibility annotation, Yang et, al. \cite{7084187} used the consistency of votes of the local regression forest in several over-segmented regions to get a confidence value of each pixel, which is called Regional Predictive Power (RPP). Compared with RCPR, RPP obtained a higher accuracy in landmark localization. Yu et al. \cite{Yu2014} proposed a Consensus of occlusion-robust Regression method (CoR) by forming a consensus from estimates arising from a set of occlusion-specific regressors. Each regressor is trained to estimate facial landmark locations under the precondition that a particular predefined region of the occluded face. CoR improved the accuracy of occlusion detection. Liu et, al. \cite{7762938} proposed Cascade Regression with Adaptive Shape Model (CRASM) for robust facial landmark localization. In each iteration, the shape-indexed appearance is used to estimate the occlusion level of each landmark, and each landmark is then weighted according to the occlusion level. Moreover, the occlusion level of the landmark acts as adaptive weights on the shape-indexed features to decrease the noise on the shape-indexed features, which improved the performance of facial landmark localization and occlusion detection. CRASM improved the performance of landmark localization and occlusion detection compared with other methods, which obtained a 80/48.45\% precision/recall for occlusion detection and NME is 6.68$\times10^{-2}$ for localization on COWF dataset.

\section{The Proposed Scheme}

\begin{figure*}[!tp]
	\centering
	\includegraphics[width=.99\linewidth]{RI
	\caption{The procedure of RICPR. The texture correlated initial shapes and the pose correlated initial shapes are calculated in parallel. The texture correlated initialization is based on the correlation of LBP histograms between the testing face and the training faces, while the pose correlated initialization is based on the evaluated rough face pose. These initial shapes are combined together as robust initializations for regression to get predictions. Finally, the reliability of each prediction is evaluated by variance to get the final output.}
	\label{robust_initial}
\end{figure*}

\label{sec:3}
In this section, we briefly review CPR and RCPR, and then describe in detail the proposed RICPR scheme for facial landmark localization under occlusions.

\subsection{Cascaded Shape Regression}

\begin{algorithm}[!htbp]
	\caption{Cascaded Pose Regression}
	\label{cpr}
	\begin{algorithmic}[1]
		\Require
		Image $I$, initial shape $S^0$, regressors $R^{1..T}$
		\For{$t=1$ to $T$}
		\State $f^t=h^t(I,S^{t-1})$ \qquad //Shape-indexed features
		\State $\Delta{S^t}=R^t(f^t)$ \qquad //Apply regressor $R^t$
		\State $S^t=S^{t-1}+\Delta{S^t}$\qquad //Update estimated shape
		\EndFor
		\Ensure
		Estimated shape $S^T$
	\end{algorithmic}
\end{algorithm}

The main steps of CPR \cite{5540094} can be described as Algorithm \ref{cpr}. It starts from a raw initial shape $S^0$, which is progressively refined in each iteration by applying a cascade of $T$ regressors $R^{t},t=1,..,T$ until the final shape $S^T$ is estimated. At each iteration, image features are calculated as $f^t=h^t(I,S^{t-1})$, where $I$ is the face image and $S^{t-1}$ is the previous iteration's shape. Based on the shape-indexed features and the regressor $R^t$, an update $\Delta{ S^t}$ is calculated. The update $\Delta{S^t}$ is combined with $S^{t-1}$ to form a new shape $S^t$. ESR \cite{6248015} proposed some improvements over CPR, which uses two-level cascaded regression to strengthen regressors. There are $K$ primitive fern regressors $R^t=(R^t_1,..,R^t_k,..,R^t_K)$ at each iteration, and the shape update $\Delta{S^t}$ is obtained by:
\begin{equation}
\label{eqn_1}
\Delta{S^t}= \sum_{k=1}^{K} \Delta{S_k^t}=\sum_{k=1}^{K} R^t_k(h^t(I,S^{t-1})).
\end{equation}
\par
Burgos-Artizzu et al. \cite{6751298} proposed a novel regression approach RCPR, to handle localization under occlusions, which divides the face image into 9 zones. At each iteration $t$, the occlusions presented in each one of the 9 zones can be estimated by projecting the current estimation $S^{t-1}$ in the image. Then, instead of training a single regressor, RCPR trains $\eta$ regressors in each primitive fern regressor $R^t_k$. Moreover, each regressor is allowed to draw features only from 1 of the 9 predefined zones. Finally, the updates of the regressors $\delta{S_{i},i=1,..,\eta}$ are combined through a weighted mean voting. The weight $w_i^k$ is inversely proportional to the occlusions presented in the zone. At the \emph{t}-th iteration, the \emph{k}-th update can be described as:
\begin{equation}
\label{eqn_2}
\Delta{S^t_k}= \sum_{i=1}^{\eta} w_i^k \delta{S_i^k},
\end{equation}
\subsection{Robust Initialization for Cascaded Pose Regression}
The procedure of the proposed RICPR is illustrated in Fig. \ref{robust_initial}. Firstly, we get the texture correlated initial shapes by calculating texture correlation between the testing face and the training faces, at the same time, the pose correlated initial shapes are obtained by examining rough face pose of the testing face. Then, these initial shapes are taken as the robust initialization for cascaded regression. We describe these two initialization methods in the subsections \ref{lbp-i} and \ref{pose-i}, respectively.

\subsection{ The Texture Correlated Initial Shapes }
\label{lbp-i}
Since occlusion and pose variation change the appearance of a face and texture descriptor captures the local appearance detail, we can select a texture correlated initial shape to consider the occlusion information of the testing face, rather than a random initial shape.
\par
We propose an initialization method based on texture correlation analysis between the testing face and the training faces. The shapes of the training faces which are most correlated with the testing face are chosen as the initialization, instead of a random one. The LBP operator was originally proposed for texture analysis which is widely used in computer vision \cite{OJALA199651}. It labels an image by thresholding the $3\times 3$-neighborhood of each pixel with the centre value, as shown in Fig. \ref{lbp} (a). The histogram of the labels can be used as a texture descriptor. To deal with the limitation of the basic LBP operator, the rotation-invariant LBP and uniform LBP were proposed \cite{1017623}, as shown in Fig. \ref{lbp} (b). In the proposed scheme, we choose the uniform LBP since it balances performance and speed. We use the notation $LBP^{u2}_{P,Q}$ to denote the LBP operator, where the subscript represents using the operator in a $(P,Q)$ neighborhood, while $(P,Q)$ means $P$ sampling points on a circle of radius $Q$. The superscript $u2$ means using only uniform patterns.
\begin{figure}[!htbp]
	\centering
	\includegraphics[width=3.45in]{lbp
	\centering
	\caption{(a) Basic LBP operator. (b) Examples of extended LBP operators: The circular (8,1), (16,2), and (8,2) neighborhoods.}
	\label{lbp}
\end{figure}

\begin{figure}[!htbp]
	\centering
	\includegraphics[width=3.45in]{lbp-histogram
	\centering
	\caption{The procedure of constructing histograms matrix.}
	\label{our_lbp}
\end{figure}
\par
Given an image, we divide the face into $m$ non-overlapping sub-blocks, as shown in Fig. \ref{our_lbp}. For each block $i$, we use $LBP^{u2}_{P,Q}$ to calculate LBP features, then a histogram of the labeled block $f_l(x,y)$ is computed as:
\begin{equation}
\label{eqn_3}
\begin{split}
H^i_j=\sum_{x,y} I\{f_l(x,y)=j\},I\{(x,y)\in Q_i\},&\\j=0,..,n-1,i=1,..,m,&
\end{split}
\end{equation}
where $n$ is the number of labels of each block produced by the LBP operator and $I\{ \cdot \}$ is defined as:
\begin{equation*}
\label{eqn_4}
I\{\emph{A}\}=
\begin{cases}
1,& \text{\emph{A} is true};\\
0,& \text{\emph{A} is false}.
\end{cases}
\end{equation*}
Finally, the $m$ histograms are combined yielding the histogram-matrix $H$ with a size of $m\times n$.
\par

\begin{figure*}[!tp]
	\centering
	\includegraphics[width=.90\linewidth]{compare-lbp-rcpr
	\centering
	\caption{We run RCPR based on 5 initial shapes selected randomly from the training set and 5 most correlated shapes from the training set by the proposed texture correlated initialization. The median of all predictions is taken as the final output. (a) The five initial shapes, where the images in the first row are the initial shapes using random initialization and the images in the second row are initial shapes using the texture correlated initialization. (b) The corresponding outputs of the two facial landmark localization methods.}
	\label{initial_compare}
\end{figure*}

During testing, histogram-matrices of the testing image and the training samples are computed by using the above scheme. It should be noticed that, to save the time cost of the testing, the histogram-matrices of the training samples can be computed offline before the testing. The best way to classify histogram-matrices is to use one of the histogram similarity measures, such as histogram intersection, log-likelihood or Chi-Square statics \cite{Ahonen2004}. Since our work aims to select proper initial shapes for regression and we hope to pick a few of the most relevant shapes with the testing face from training faces, we need a method to assess the correlation between the testing face and the training faces. In this paper, we choose the Pearson correlation coefficient \cite{10.2307/115794} to measure the correlation between the testing face and the training faces. The Pearson correlation coefficient between the testing face histogram-matrix $TM$ and each training face histogram-matrix $TRM^\gamma$ is calculated by:

\begin{algorithm}[!tp]
	\caption{Initialization based on texture correlation }
	\label{initial}
	\begin{algorithmic}[1]
		\Require
		testing face $TI$, training faces $\{TRI\}^{1..\Gamma}$, shapes of training faces $\{TRS\}^{1..\Gamma}$
		\Ensure
		Texture correlated initial shapes $\{TS_0\}^{1..l}$
		\State Compute LBP histograms for all faces
		\State Obtain histograms-matrices $TM$ and $\{TRM\}^{1..\Gamma}$ for $TI$ and $\{TRI\}^{1..\Gamma}$
		\For{$\gamma=1$ to $\Gamma$}
		\State Calculate correlation coefficient $\rho_\gamma$ between $TM$ and $TRM^\gamma$
		\State Get correlation distance $d_\gamma$
		\EndFor
		\State Search the $l$ smallest correlation distances
		\State Select corresponding $l$ shapes in $\{TRS\}^{1..\Gamma}$ as $\{TS_0\}^{1..l}$
	\end{algorithmic}
\end{algorithm}

\begin{equation}
\label{eqn_5}
\begin{aligned}
\rho_\gamma&=\frac{Cov(TM,TRM^\gamma)} {\sigma{(TM)}\sigma{(TRM^\gamma)}}\\&=\frac{E[(TM-E(TM))(TRM^\gamma-E(TRM^\gamma))]} {\sigma{(TM)}\sigma{(TRM^\gamma)}},\\&\gamma=1,..,\Gamma,
\end{aligned}
\end{equation}
where $Cov(\cdot,\cdot)$ is the covariance, $\sigma(\cdot)$ is the standard deviation and $\Gamma$ is the total number of training faces. Due to the size of each histogram-matrix is $m\times n$, $\rho_\gamma$ can be calculated as:
\begin{equation*}
\label{eqn_6}
\frac{\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{n} (TM_{ij}-\overline{TM})(TRM_{ij}^\gamma-\overline{TRM^\gamma})} {\sqrt{\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{n} (TM_{ij}-\overline{TM})^2} \sqrt{\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{n} (TRM_{ij}^\gamma-\overline{TRM^\gamma})^2},}
\end{equation*}
where $\overline{TM}$ and $\overline{TRM^\gamma}$ are mean values of matrices ${TM}$ and ${TRM^\gamma}$, respectively. Then the correlation coefficient $\rho^\gamma$ can be used to calculate correlation distance $d_\gamma$:
\begin{equation}
\label{eqn_7}
d_\gamma=1-\rho_\gamma, \gamma=1,..,\Gamma.
\end{equation}
\par
A smaller $d$ represents that the training face is more correlated with the testing face. We choose the $l$ most correlated faces from $\Gamma$ training faces and select their shapes as initial shapes for the testing face. The main procedure of initialization based on texture correlation is presented in Algorithm \ref{initial}. As shown in Fig. \ref{initial_compare}, a comparison between random initialization \cite{6751298} and texture correlation based initialization is illustrated. It can be found that the proposed texture correlation based initialization usually obtains more accurate initial shapes which improves the accuracy of landmark localization.

\subsection{The Pose Correlated Initial Shapes}
\label{pose-i}
In the above section, we describe how to select the texture correlated initial shapes considering the occlusion information but ignoring pose information of the testing face. Empirically, landmark distribution is highly correlated to head pose. To further make the initial shapes more robust to various poses, we choose some pose correlated initial shapes for regression.
\par
To obtain the pose correlated initial shapes, we estimate the rough pose of the testing face, which can be obtained by the five fiducial landmarks, i.e., the pupils, the tip of the nose, and the corners of the mouth. In this paper, we use MTCNN \cite{7553523} to detect the five fiducial landmarks, as shown in Fig. \ref{pose_initial}. Inspired by Perspective-n-Point (PnP) problem, which is the problem of estimating the pose of a calibrated camera given a set of 3D points and their corresponding 2D projections in the image \cite{Lepetit2009EPnP}. Given a 3D mean shape $\overline{S}$ with 5 facial key points and the five detected fiducial landmarks, a rough face pose can be estimated by:
\begin{equation}
\label{eqn_8}
\vec{\theta} = g( \overline{S} , p),
\end{equation}
where $\vec{\theta}$ is a rotation vector, which represents the face pose and $p$ represents the five fiducial landmarks detected by MTCNN, $g$ is the Efficient PnP \cite{Lepetit2009EPnP}.
\par

\begin{figure}[!tp]
	\centering
	\includegraphics[width=.99\linewidth]{poseInitial
	\caption{Illustration of generating the pose correlated shape. Given an image, we first detect five fiducial landmarks and estimate face pose. Then, according to the face pose, a 3D mean face shape with 29 facial key points, can be projected to a set of corresponding 2D locations, which has similar pose with testing image.}
	\label{pose_initial}
\end{figure}

\begin{figure*}[!tp]
	\centering
	\subfigure[]{\label{error_rank}
		\begin{minipage}[b]{0.4\textwidth}
			\includegraphics[width=1\textwidth]{lbp-i-error}
		\end{minipage}
	}
	\subfigure[]{\label{good_rank}
		\begin{minipage}[b]{0.4\textwidth}
			\includegraphics[width=1\textwidth]{lbp-i-good1}
		\end{minipage}
	}
	\caption{Comparisons between the texture correlated initialization based RCPR and the traditional random initialization based RCPR. (a) The NME on various initial shapes with different correlation distances in training and testing processes. (b) The number of ``good'' instances determined by variance after 10\% cascades of each prediction on various correlation distances. Correlation distances in (a) and (b) are ranked in ascending order.}
	\label{ini_result}
\end{figure*}

\begin{table*}[!tp]
	\renewcommand\arraystretch{1.5}
	\begin{center}
		\normalsize
		\caption{\label{feture_compare} Texture Correlated Initialization Using Different Features }
		\begin{tabular}{|p{2.4cm}<{\centering}|p{1.7cm}<{\centering}|p{1.7cm}<{\centering}| p{1.7cm}<{\centering}|p{1.7cm}<{\centering}|p{1.7cm}<{\centering}|p{1.7cm}<{\centering}|p{1.7cm}<{\centering}|}
			\hline
			\textbf{ Feature} &\textbf{LBP} &\textbf{LDP} &\textbf{Gabor} &\textbf{GMRF} &\textbf{GLDS} &\textbf{GLCM} &\textbf{Eigenface}\\
			\hline
			\textbf{NME($\times10^{-2}$)}& 7.35 & 7.75 & 7.87 & 8.28 & 8.19 & 8.06 &8.18\\
			\hline
			\textbf{Precision/recall}& 80/51.4\% & 80/48.7\% & 80/46.1\% & 80/45.6\% & 80/47.2\% & 80/46.5\% &80/47.6\%\\
			\hline
		\end{tabular}
	\end{center}
	\begin{tablenotes}
		\footnotesize
		\item Accuracy of facial landmark localization and occlusion detection based on texture correlated initialization using different features. The results indicate that the LBP performs better than the others.
	\end{tablenotes}
\end{table*}

Then, a 3D mean face shape, represented by 29 facial landmark locations, is projected to a set of corresponding 2D locations according to the testing face pose $\vec{\theta}$, as shown in Fig. \ref{pose_initial}. After that, the shape which has similar pose with the testing face is obtained. To get a reasonable initial shape for each image, we re-scale the corresponding 2D locations based on the face bounding box and the detected five fiducial landmarks $p$. The initial occlusion information of the pose correlated initial shape is distributed randomly as:
\begin{equation}
\label{eqn_9}
S_0 = G( \overline{S^*} , \vec{\theta}, b, p),
\end{equation}
where $b$ is the face bounding box, $\overline{S^*}$ is the 3D mean face shape with 29 points, and $S_0$ is the pose correlated initial shape. To achieve a better performance, we selected several frontal faces from the training set to augment the pose correlated initial shape. Referring to their true 2D shapes and the 3D mean face shape $\overline{S^*}$, we construct different 3D frontal face shapes that have little variation compared with the 3D mean face shape. Then, based on Eq. \ref{eqn_9}, different initial shapes can be generated by replacing $\overline{S^*}$ with the constructed 3D frontal face shapes.

\subsection{Variance Evaluation}
As stated in \cite{6751298}, due to the coarse to fine nature of CPR, even if a face image is initialized by several different shapes, the predictions should reach a similarity after iterations. Based on this principle, instead of taking the median of all predicted results as the final output, the variance is used to determine the reliability of two initialization methods' predictions.
\par
Firstly, after finishing the regression, the variance $v$ of all predictions is calculated. If the value of $v$ is below a certain threshold $\zeta$, it indicates that the predictions is a good solution. In this case, all predictions are considered as reliable, thus we take the median of all predictions as the final output. Otherwise, part of predictions belong to ``bad'' class, then the variances between predictions based on two initialization methods are computed and represented as $v_{lbp}$ and $v_{pose}$. If $v_{lbp}$ is less than $v_{pose}$, it indicates that the predictions based on the texture correlated initialization are more reliable than these based on the pose correlated initialization. Therefore, only considering the predictions by the texture correlated initial shapes, we abandon the predictions which make obvious variance variation and take the median of the rest of predictions as the final output. If $v_{lbp}$ is greater than $v_{pose}$, it indicates that the predictions based on the pose correlated initialization are more reliable than those based on the texture correlated initialization. Then, the median of the predictions based on the pose correlated initial shapes are taken as the final output.

\section{Experimental Results }
\label{sec:4}

\subsection{Dataset and Implementation}

\begin{table}[!tp]
	\begin{center}
		\begin{threeparttable}

			\renewcommand\arraystretch{1.4}
			\caption{\label{err/fail}
				Results on COFW Dataset.}
			\begin{tabular}{l | c | c }
			\hline
				\multirow{2}*{Methods} & {Landmark localization error} & {Occlusion prediction} \\
				\cline{2-3}
				 & {NME ($\times10^{-2}$)} & {Precision/Recall} \\
				\hline
				\textbf{RCPR} \cite{6751298} & 8.01 & 80/42\%\\
				\textbf{HPM} \cite{Ghiasi_2014_CVPR} & 7.46* & 80/37\%*\\
				\textbf{RPP} \cite{7084187} & 7.52* & 78/40\%*\\
				\textbf{SDM} \cite{Xiong_2013_CVPR} & 10.88 & -\\
				\textbf{TCDCN} \cite{7208848} & 8.05* & -\\
				\textbf{CRASM} \cite{7762938} & 6.68* & 80/48.45\%*\\
				\textbf{HORSD} \cite{7911334} & 6.8* & -\\
				\textbf{LBP-I-RCPR} & 7.35 & 80/51.4\%\\
				\textbf{RICPR} & 6.64 & 80/54.6\%\\
				\textbf{Human} \cite{6751298} & 5.6 & - \\
				\hline
			\end{tabular}

			\begin{tablenotes}
				\footnotesize
				\item Comparison of facial landmark localization and occlusion prediction on COFW dataset. The table lists the results of NME and occlusion detection. * indicates that the result is from the published paper.
			\end{tablenotes}
		\end{threeparttable}
	\end{center}
\end{table}

We evaluate the performance of the proposed scheme on the challenging dataset COFW \cite{6751298}, which is widely used to evaluate the robustness of facial landmark localization and occlusion detection. The face images in COFW have large variations in shape and occlusions due to differences in pose, expression, hairstyle, using of accessories such as sunglasses, hats and interactions with objects (e.g. food, hands, microphones, etc.). Each image is annotated with the location and occluded/un-occluded state of 29 facial landmarks. This dataset has 1852 face images in total, where 1345 and 507 images are used for training and testing respectively. The average occlusion rate of faces in COFW is over 23\%.
\par
To evaluate the performance of the proposed scheme, we implement the proposed scheme with two configurations. One is RCPR based on LBP histogram correlation initialization (LBP-I-RCPR), which represents the initialization only based on texture correlation, and the median of predictions is taken as the final output. The other is the full version of the proposed scheme RICPR, in which the initialization is based jointly on both texture correlation and pose correlation jointly. In texture correlation analysis, the face, whose location is provided by a face detector, is divided into $m=8\times 8=64$ non-overlapping sub-blocks. The uniform LBP in a $(P,Q)=(8,1)$ neighborhood is employed to obtain texture information of the face. Thus, the number of labels produced by the LBP operator is 59. In pose correlation analysis, we utilize MTCNN \cite{7553523} to predict five fiducial landmarks. A threshold $\zeta=0.08$ is used to determine whether the predictions lead a good result. Since the proposed scheme always has a good initialization without the need for the smart restarts and the number of initial shapes is set to 10 and $\eta$ is set to 4. We run RICPR and RCPR using the same configuration.
\par
We compare LBP-I-RCPR and RICPR with several state-of-the-art methods on COFW dataset using NME (Normalized Mean Error) defined by Eq. \ref{eqn_10}.
\begin{equation}
\label{eqn_10}
NME=\frac{1}{N}\sum\limits_{i}^{N} \frac{\frac{1}{M}\sum\limits_{j}^{M}|s^p_{i,j}-s^g_{i,j}|_2} {|gl_i-gr_i|_2},
\end{equation}
where $N$ is the number of images in the test set, $M$ is the number of landmarks in one image, $s^p_{i,j}$ is the predicted position of the $j^{th}$ landmark of the $i^{th}$ image, $s^g_{i,j}$ is the ground truth position of the $j^{th}$ landmark of the $i^{th}$ image, $gl_i$ and $gr_i$ are the ground truth positions of the left and right eye centres respectively.
\par
Based on NME, we can plot the Cumulative Error Distribution (CED) curves to further analyse the performance of the proposed scheme, which is calculated from the NME over each image. We also evaluate the speed of the proposed scheme on the COFW dataset. Speed is measured in Frames Per Second (FPS). All methods are implemented using Matlab R2015b and run on a PC with 3.60 GHz CPU and 64-bit Windows 7 operating system.

\begin{figure}[!tp]
	\centering
	\includegraphics[width=3.4in]{ced2
	\centering
	\caption{CED curves on the COFW dataset.}
	\label{ced}
\end{figure}

\subsection{Results}
\textbf{1) Analysis of initialization based on texture correlation:} Instead of randomly selecting shapes from the training set as the initialization in RCPR, we employ a texture correlated initialization (LBP-I-RCPR) by computing LBP histograms. To prove the effectiveness of the texture correlated initialization method, we compare the performance of LBP-I-RCPR with RCPR on the COFW dataset as shown in Fig. \ref{ini_result}.
\par
The NMEs on various initial shapes with different correlation distances are shown in Fig. \ref{error_rank}. The results show that the NME is reduced with decreasing correlation distance and LBP-I-RCPR can significantly reduces NME by at least 45\%. It indicates that the initial shapes which are selected from the training faces based on texture correlation is closer to the real shape of the testing face.
\par
Moreover, given different initial shapes for each image, the variance between their predictions is applied to determine whether the face belongs to a ``good'' class as stated in \cite{6751298}. As shown in Fig. \ref{good_rank}, it can be found that the number of ``good'' instances increases as correlation distance decreases and more images belong to ``good'' class among 507 testing images when using LBP-I-RCPR. The number of ``good'' instances dramatically increase by at least 45\%, and thus less bad initial shapes are selected. Furthermore, the number of ``good'' instances increases from 395 to 504 among the 507 images in RICPR scheme, which means fewer than 1\% instances are ``bad'', thus the initialization become more robust.
\par
We also initialize the shapes using other different features, including Local Derivative Pattern (LDP) \cite{Zhang2010Local}, Gabor, Gaussian Markov Random Field (GMRF), Gray-Level Difference Statistics (GLDS), Gray-Level Co-occurrence Matrix (GLCM), and Eigenface. We report the NME and occlusion detection of each feature respectively in Table \ref{feture_compare}. The results indicate that the initialization based on LBP histogram correlation performs better.

\begin{figure}[!tp]
	\centering
	\includegraphics[width=3.4in]{pre-recall
	\centering
	\caption{Occlusion detection result on the COFW dataset.}
	\label{re/pre}
\end{figure}

\par
\textbf{2) Facial landmark localization evaluation on COFW:} Many facial landmark localization methods perform not well on the COFW database due to the large variation in occlusion. To evaluate the proposed scheme, we compare the proposed scheme with several state-of-the-art methods including RCPR \cite{6751298}, RPP \cite{7084187}, SDM \cite{Xiong_2013_CVPR}, Tasks-Constrained Deep Convolutional Network (TCDCN) \cite{7208848}, Hierarchical Deformable Part Model (HPM) \cite{Ghiasi_2014_CVPR}, CRASM \cite{7762938} and Hierarchical Occlusion Stage-wise Relational Dictionary (HOSRD)\cite{7911334}. The comparisons of NME on COFW dataset are given in Table \ref{err/fail}.
\par
We can find that RICPR obtains the smallest NME. Compared to RCPR, the LBP-I-RCPR reduces the NME from 8.01$\times10^{-2}$ to 7.35$\times10^{-2}$ and the RICPR further reduces the NME to 6.64$\times10^{-2}$. The NME is reduced by 17.1\% in total. RICPR performs even better than the most recent CRASM method proposed in 2017. To get the pose correlated initial shapes, we use MTCNN to detect five fiducial landmarks. The accuracy of five fiducial landmarks plays a significant role on performance. If the ground-truth of the five fiducial landmarks is employed, the NME can reach 5.52$\times10^{-2}$, which demonstrates that the proposed scheme can obtain a admirable performance if the five fiducial landmarks are detected accurately.
\par
We also show the CED curves of the COFW dataset in Fig. \ref{ced}. As can be seen, more images perform better using the proposed scheme, it also demonstrates the superiority of the proposed scheme for facial landmark localization in face image with occlusions.

\begin{figure}[!tp]
	\centering
	\includegraphics[width=3.4in]{sdm-rcpr
	\centering
	\caption{Results of SDM and RCPR using the proposed initialization methods.}
	\label{sdm-rcpr}
\end{figure}

\par
\textbf{3) Occlusion detection on COFW:} Since the COFW dataset provides the ground truth of occlusion, we evaluate the occlusion detection on COFW and compare the proposed scheme with RCPR \cite{6751298}, HPM \cite{Ghiasi_2014_CVPR}, CoR \cite{Yu2014}, RPP \cite{7084187} and CRASM \cite{7762938}. The occlusion prediction results are shown in Table \ref{err/fail} and Fig. \ref{re/pre}. As can be seen, the proposed scheme also outperforms the state-of-the-art methods in occlusion detection.
\par
When we set the false alarm at 80\%, the proposed scheme achieves an accuracy of 54.6\%, which is higher than 42\% obtained by RCPR, 37\% obtained by HPM, 41.44\% obtained by CoR, 48.45\% obtained by CoR and 78/40\% precision/recall obtained by RPP. Even if only using LBP-I-RCPR scheme, the accuracy of detecting occlusion reaches 51.4\%. It demonstrates that the proposed scheme achieves a much higher accuracy of occlusion detection, which can provide significant benefits in real world application, such as image texture analysis, facial expression understanding and face recognition. Fig. \ref{sample_result} shows example images with the result obtained by the proposed RICPR.

\par
\textbf{4) Run time:} We record the speeds of RCPR, LBP-I-RCPR and RICPR on the COFW dataset. The speeds of these methods are 5.3 FPS, 4.1 FPS and 4.0 FPS, respectively. We can find that the proposed scheme takes some time on calculating the correlation. The speed can be improved by implementing it with C++ or using a powerful server. We will try to improve the efficiency of the proposed scheme in the future, for example, by reducing the number of face images used for texture correlation based initialization.

\subsection{Generalization of the Proposed Initialization Scheme}
The experimental results demonstrate that the proposed initialization scheme significantly improved the performance of RCPR in both localization and occlusion prediction. Since the initialization is usually independent to facial landmark localization, the proposed initialization scheme can be applied to other algorithms such as SDM. The results are shown in Fig.\ref{sdm-rcpr}, where baseline is the original SDM or RCPR, LBP-I+baseline is the texture correlated initialization scheme applied to SDM or RCPR, RI+baseline denotes the joint texture correlation and pose correlation initialization scheme applied to SDM or RCPR. Compared with the original SDM which is based on random initialization, the LBP-I-SDM and the RI-SDM reduce the NME by 14.6\% and 19\% respectively. The results indicate that the proposed initialization scheme can also improve the performance of SDM.

\begin{figure*}[!htbp]
	\centering
	\includegraphics[width=.95\linewidth]{example
	\caption{Example result of facial landmark localization and occlusion detection obtained by the proposed RICPR on the COFW dataset.}
	\label{sample_result}
\end{figure*}

\section{Conclusions}
\label{sec:5}
In this paper, we propose a robust initialization scheme to solve the initialization sensitive problem for the cascaded pose regression approach through jointly analyzing texture and pose of a testing face. By examining the correlation of local binary patterns histograms between the testing face and the training faces, the texture correlated shapes are selected instead of random shapes. At the same time, the pose correlated initialization is proposed to further improve the robustness of the initialization by estimating the face pose. Experimental results show that the proposed scheme obtains remarkably higher accuracies on both facial landmark localization and occlusion detection on facial images than the state-of-the-art benchmarks. Moreover, since the initialization is usually independent with facial landmark localization, the proposed initialization scheme has the potential to be extended and applied to other algorithms.

\section*{Acknowledgment}
This work is supported by the National Natural Science Foundation of China (Grant No. 61601337).

\bibliographystyle{IEEEtran}

redpajama_set_name: RedPajamaArXiv