# A Linearly Convergent Algorithm for Rotationally Invariant $\ell_{1}$-Norm
Principal Component Analysis

Taoli Zheng, Peng Wang, and Anthony Man-Cho So, Senior Member, IEEE
T. Zheng and A. M.-C. So are with the Department of Systems Engineering and
Engineering Management, The Chinese University of Hong Kong, Hong Kong SAR,
China (e-mails: {tlzheng,manchoso}@se.cuhk.edu.hk). P. Wang is with the
Department of Electrical Engineering and Computer Science, University of
Michigan, Ann Arbor (e-mail: pengwa@umich.edu).

###### Abstract

To do dimensionality reduction on the datasets with outliers, the
$\ell_{1}$-norm principal component analysis (L1-PCA) as a typical robust
alternative of the conventional PCA has enjoyed great popularity over the past
years. In this work, we consider a rotationally invariant L1-PCA, which is
hardly studied in the literature. To tackle it, we propose a proximal
alternating linearized minimization method with a nonlinear extrapolation for
solving its two-block reformulation. Moreover, we show that the proposed
method converges at least linearly to a limiting critical point of the
reformulated problem. Such a point is proved to be a critical point of the
original problem under a condition imposed on the step size. Finally, we
conduct numerical experiments on both synthetic and real datasets to support
our theoretical developments and demonstrate the efficacy of our approach.

###### Index Terms:

Rotational invariance, $\ell_{1}$-norm principal component analysis, proximal
alternating linearized minimization, non-linear extrapolation.

##  I Introduction

With the widespread availability of high-dimensional data, principal component
analysis (PCA) [1] as a typical dimensionality reduction technique plays an
increasingly important role in data analysis. It has found extensive
applications in diverse fields such as computer vision [2, 3, 4], image
processing [5, 6], and network analysis [7, 8], to name a few. Due to the fact
that the conventional PCA (also known as L2-PCA) is sensitive to corruptions
(generically referred to as outliers [9]) in the datasets, L1-PCA as a robust
alternative of L2-PCA has recently received significant attention; see, e.g.,
[10, 11]. In this work, we consider a particular form of L1-PCA as follows:

$\max\left\\{\|\bm{Q}\bm{Q}^{T}\bm{X}\|_{1}:\bm{Q}\in{\rm St}(d,K)\right\\}.$
(1)

Here, $\bm{X}=\left[\bm{x}_{1},\dots,\bm{x}_{n}\right]\in\mathbb{R}^{d\times
n}$ is the data matrix with the sample mean being zero, where $n$ and $d$
respectively denote the number of samples and dimension of the data points,
$\|\bm{A}\|_{1}=\sum_{i,j}|a_{ij}|$ denotes the $\ell_{1}$-norm of the matrix
$\bm{A}$, $K$ is the dimension of subspace with $K\leq\min\\{n,d\\}$, and
${\rm St}(d,K)=\left\\{\bm{Q}\in\mathbb{R}^{d\times
K}:\bm{Q}^{T}\bm{Q}=\bm{I}_{K}\right\\}$ denotes the Stiefel manifold.
Generally, this problem can be used to find a low-dimensional subspace
underlying a corrupted dataset by minimizing the variation of the projections
of data points onto the subspace measured by the $\ell_{1}$-norm. Different
from the formulation of Problem (1), there are actually many other
formulations for L1-PCA that have been well studied in the literature; see,
e.g., [12, 13, 14, 15, 16, 17]. However, as pointed out in [11], no work
focuses on Problem (1). To fill this gap, this work is devoted to develop an
efficient algorithm for solving Problem (1) with convergence analysis.

Now, we review different formulations for L1-PCA in the literature. It is
known that PCA admits two common interpretations for dimension reduction. One
is that PCA aims to find a low-dimensional subspace for which the projections
of data points onto it preserve the most variance. Specifically, the
corresponding formulation known as L2-PCA is given by

$\max\left\\{\|\bm{Q}\bm{Q}^{T}\bm{X}\|_{F}=\|\bm{Q}^{T}\bm{X}\|_{F}:\bm{Q}\in{\rm
St}(d,K)\right\\},$ (2)

The other sheds light on another perspective of PCA that aims to find a low-
dimensional subspace for which the projections of data points onto it best
approximate these data points. This leads to a minimization formulation of the
form

$\min\left\\{\|\bm{X}-\bm{Q}\bm{Q}^{T}\bm{X}\|_{F}:\bm{Q}\in{\rm
St}(d,K)\right\\}.$ (3)

One can verify that this problem is equivalent to Problem (2). In particular,
both of them admit a closed-form solution that can be solved efficiently by
computing the singular value decomposition (SVD) of the data matrix $\bm{X}$.
Moreover, the subspace spanned by the columns of the obtained solution
possesses many nice properties [1]. Nevertheless, the above two formulations
have an essential defect that they are sensitive to outliers in the data
matrix due to the fact that the $\ell_{2}$-norm is not a robust scale
function; see, e.g., [10]. The outliers could significantly disturb the
singular values of the data matrix and mislead the learning of the desired
subspace. To tackle this problem, a recent line of research on robust subspace
recovery becomes more and more popular [11] and a great deal of work adopting
different norms as scale functions have been done in the literature; see,
e.g., [18, 17, 16, 13, 15]. Notably, taking the $\ell_{1}$ norm as the scale
function gives L1-PCA. On one hand, along the line of maximizing the variance
of projected data as in Problem (2), we can naturally obtain Problem (1) and

$\max\left\\{\|\bm{Q}^{T}\bm{X}\|_{1}:\bm{Q}\in{\rm St}(d,K)\right\\}.$ (4)

Compared to L2-PCA, the $\ell_{1}$-norm endows PCA with robustness but
introduces non-smoothness into the objective function simultaneously.
Consequently, Problems (1) and (4) are both non-convex and non-smooth, and
have no closed-form solutions. Unlike the equivalence of L2-PCA in Problem
(2), Problems (1) and (4) are generally not equivalent. Indeed, the objective
of the latter problem is the pointwise maximum of linear functions, while the
former is the pointwise maximum of quadratic functions. In particular, Problem
(4) is widely studied in the literature. For example, Kwak [19] proposed a
fixed-point method for solving Problem (4) in a greedy manner and later Nie et
al. [20] improved it to a non-greedy manner. Recently, Markopoulos et al. [21]
developed an algorithm based on bit-flipping iteration. More recently, Wang et
al. [22, 14] proposed a proximal alternating minimization method with
extrapolation for solving it. We refer the reader to [23] for more algorithms
for solving Problem (4). Besides, it is known that L2-PCA is rotationally
invariant in the sense that the objective value remains unchanged under
arbitrary subspace rotations. Specifically, the objective values of Problems
(2) and (3) are unchanged if we replace $\bm{Q}$ by $\bm{Q}\bm{U}$ for any
orthogonal matrix $\bm{U}\in\mathbb{R}^{K\times K}$. However, this property
generally does not hold for the L1-PCA formulation (4). Fortunately, the
formulation (1) keeps this property. This also motivates us to study Problem
(1). In the following, we refer to it as rotationally invariant L1-PCA. On the
other hand, along the line of minimizing the approximation error of data
points as in Problem (3), we can obtain another L1-PCA formulation of the form

$\min\left\\{\|\bm{X}-\bm{Q}\bm{Q}^{T}\bm{X}\|_{1}:\bm{Q}\in{\rm
St}(d,K)\right\\}.$ (5)

This problem has been considered in the literature. For example, Baccini et
al. [12] proposed a heuristic approach based on a canonical correlation
analysis. Recently, Park et al. [16] proposed an iteratively reweighted least
squares (IRLS) methods and provided convergence analysis. More related
algorithms can be found in [24, 25].

Then, we review some other robust alternatives of L2-PCA closely related to
L1-PCA. A popular one is the low-rank matrix factorization formulation of the
form

$\min\left\\{\|\bm{X}-\bm{U}\bm{V}\|_{1}:\bm{U}\in\mathbb{R}^{d\times K},\
\bm{V}\in\mathbb{R}^{K\times n}\right\\}.$ (6)

Ke and Kanade [17] proposed an alternating convex minimization method for
solving it. Later, Yu et al. [26] considered an augmented lagrange multiplier
method and Eriksson and Hengel proposed a generalization of the Wiberg
algorithm [27] for tackling it. Up to now, the above formulations all utilize
the $\ell_{1}$-norm as the scale function to achieve robustness. Actually,
many other different norms have been employed in the literature. A notable one
is the $\ell_{2,1}$-norm, which is defined as
$\|\bm{A}\|_{2,1}=\sum_{i}(\sum_{j}a_{ij}^{2})^{1/2}$. This is also known as
the $R_{1}$-norm in [13]. This norm, together with the second interpretation
of PCA, motivates the R1-PCA formulation as follows:

$\min\left\\{\|\bm{X}-\bm{Q}\bm{Q}^{T}\bm{X}\|_{2,1}:\bm{Q}\in{\rm
St}(d,K)\right\\}.$ (7)

It is noteworthy that this problem not only satisfies the mentioned subspace
rotational invariance but also satisfy another rotationally invariant property
in the sense that it is invariant to the change of coordinates in
$\mathbb{R}^{d}$. Specifically, under data transformations
$\bm{X}\rightarrow\bm{U}\bm{X}$ for any orthogonal matrix
$\bm{U}\in\mathbb{R}^{d\times d}$, the objective value is unchanged when
$\bm{Q}\rightarrow\bm{U}\bm{Q}$. That is to say that when the sample space is
rotated, the feature subspace will be rotated in the same way. For the sake of
clarity, in our paper, we refer to rotationally invariance as invariance with
respect to _subspace rotations_. To solve this problem, Ding et al. [13]
proposed a subspace iteration algorithm. Later, Wang et al. [28] generalized
this robust formulation via replacing $\ell_{2,1}$-norm with $\ell_{2,p}$-norm
for any $p\in(0,2)$. Recently, Nie et al. [15] considered a robust PCA model
using the $\ell_{2,1}$-norm based on the first interpretation of PCA, i.e.,

$\max\left\\{\|\bm{Q}^{T}\bm{X}\|_{2,1}:\bm{Q}\in{\rm St}(d,K)\right\\}.$ (8)

They developed an efficient non-greedy method to solve this problem. We refer
the reader to [11] for more formulations of robust PCA.

###  I-A Contribution

In this work, our first contribution is to propose a _proximal alternating
linearized minimization_ (PALM) algorithm for solving a two-block
reformulation (see Problem (10)) of Problem (1). Motivated by the encouraging
performance of extrapolation techniques for accelerating non-convex
optimization problems (see, e.g., [29, 30, 31, 14]), we incorporate a
quadratic extrapolation into the update of a block variable to possibly
accelerate the PALM method. It is worth noting that compared to the widely
used linear extrapolation in the literature (see, e.g., [29, 30, 31, 14]), our
proposed quadratic extrapolation (see (13)) seems to be new. We refer to the
resulting method as proximal alternating linearized minimization with
extrapolation (PALMe). PALMe is computationally efficient since the update of
each block of PALMe admits a closed-form solution and per-iteration cost is
$\mathcal{O}(ndK+dK^{2})$. Our second contribution is to show that the
proposed method converges at least linearly to a limiting critical point (see
Definition 1) of the reformulation problem and this point is proved to be a
critical point of Problem (1) (see Definition 2) under a verifiable condition
imposed on the step size. To this end, we show that the Kurdyka-Łojasiewicz
(KŁ) exponent (see Definition 3) of Problems (10) is $1/2$. With this
characterization, we employ the convergence analysis framework in [32] to show
the linear convergence of PALMe. Our third contribution is to conduct
numerical experiments on synthetic and real data sets to support our
theoretical results. Specifically, the experimental results demonstrate that
our approach is competitive, in terms of both numerical efficiency and
clustering accuracy, with other PALM-type methods for solving Problem (1). We
also compare our approach with various robust PCA approaches on image
reconstruction, which illustrates the robustness to outliers of our approach.

Lastly, let us highlight the differences between this work and a closely
related one [14]. First, the proximal alternating minimization method with
extrapolation (PAMe) in [14] is developed for solving Problem (4), in which
the objective is a pointwise maximum of linear functions. By contrast, the
objective of our considered Problem (10) is a pointwise maximum of quadratic
functions and thus PAMe is not applicable to ours. Second, we apply a
quadratic extrapolation step to the update of the block variable $\bm{P}$,
which highly relies on the quadratic form of the objective function, while the
authors in [14] utilize the standard linear extrapolation step. Moreover, this
extrapolation scheme also gives rise to the third difference between our
convergence analysis and that in [14]. Specifically, the quadratic
extrapolation step makes it more complicated to verify the decreasing
property.

###  I-B Notation and Definitions

Let $\mathbb{R}^{n}$ be the $n$-dimensional Euclidean space. We write the
matrices in bold capital letter like $\bm{Q}$, vectors in bold low-case
letters like $\bm{q}$, and scalars in plain letters like $q$. Given a matrix
$\bm{X}\in\mathbb{R}^{d\times n}$, we use $\|\bm{X}\|_{F}$ to denote its
Frobenius norm, $\|\bm{X}\|$ its spectral norm, and $x_{ij}$ its $(i,j)$-th
element. For any $x\in\mathbb{R}$, let

$\operatorname{sgn}(x)\in\left\\{\begin{array}[]{c@{\,\,\,}l}\\{x/|x|\\},&x\not=0,\\\
\vskip 3.0pt plus 1.0pt minus 1.0pt\cr\\{-1,1\\},&x=0\end{array}\right.$

denote its sign function.

Next, we introduce some standard concepts in non-smooth analysis for our
development. The details can be found in, e.g., [33]. For a non-empty closed
set $\mathcal{S}\subseteq\mathbb{R}^{p}$, the _indicator function_
$\delta_{\mathcal{S}}:\mathbb{R}^{p}\rightarrow\\{0,+\infty\\}$ associated
with $\mathcal{S}$ is defined as

$\delta_{\mathcal{S}}(\bm{x})=\left\\{\begin{array}[]{c@{\,\,\,}l}0,&\bm{x}\in\mathcal{S},\\\
+\infty,&\mbox{otherwise}.\end{array}\right.$

Let $f:\mathbb{R}^{p}\rightarrow(-\infty,+\infty]$ be a given function with
$\operatorname{dom}(f)=\\{\bm{x}\in\mathbb{R}^{p}:f(\bm{x})<+\infty\\}$. The
function $f$ is said to be _proper_ if $\operatorname{dom}(f)\not=\emptyset$.
A vector $\bm{v}\in\mathbb{R}^{p}$ is said to be a _Fréchet subgradient_ of
$f$ at $\bm{x}\in\operatorname{dom}(f)$ if

$\liminf_{\bm{y}\rightarrow\bm{x},\atop\bm{y}\not=\bm{x}}\frac{f(\bm{y})-f(\bm{x})-\langle\bm{v},\bm{y}-\bm{x}\rangle}{\|\bm{y}-\bm{x}\|_{F}}\geq
0.$ (9)

The set of vectors $\bm{v}\in\mathbb{R}^{p}$ satisfying eq. 9 is called the
_Fréchet subdifferential_ of $f$ at $\bm{x}\in\operatorname{dom}(f)$ and
denoted by $\widehat{\partial}f(\bm{x})$. The _limiting subdifferential_ , or
simply the _subdifferential_ , of $f$ at $\bm{x}\in\operatorname{dom}(f)$ is
defined as

$\partial f(\bm{x})=\left\\{\begin{split}&\
\bm{v}\in\mathbb{R}^{p}:\exists\bm{x}^{k}\rightarrow\bm{x},\,\bm{v}^{k}\rightarrow\bm{v}\,\\\
&\ \mbox{ with }\,f(\bm{x}^{k})\rightarrow
f(\bm{x}),\,\bm{v}^{k}\in\widehat{\partial}f(\bm{x}^{k})\end{split}\right\\}.$

Then, the limiting critical point of a function can be defined as follows.

######  Definition 1 (limiting critical point).

Suppose that the function $f:\mathbb{R}^{p}\rightarrow(-\infty,+\infty]$ is
proper and lower semicontinuous. We say that a point $\bm{x}\in\mathbb{R}^{p}$
is a limiting critical point of $f$ if $\bm{0}\in\partial f(\bm{x})$.

Remark that by the generalized Fermat rule (see, e.g., [34, Theorem 10.1]), a
local minimizer of $f$ is a limiting critical point of $f$. By convention, if
$\bm{x}\not\in\operatorname{dom}(f)$, then $\partial f(\bm{x})=\emptyset$. The
_domain_ of $\partial f$ is defined as $\operatorname{dom}(\partial
f)=\\{\bm{x}\in\mathbb{R}^{p}:\partial f(\bm{x})\not=\emptyset\\}$. The
limiting subdifferential of indicator function
$\delta_{\mathcal{S}}:\mathbb{R}^{p}\rightarrow\\{0,+\infty\\}$ is given as
follows.

$\displaystyle\
\widehat{\partial}\delta_{\mathcal{S}}(\bm{x})=\left\\{\bm{v}\in\mathbb{R}^{p}:\limsup_{\bm{y}\rightarrow\bm{x},\,\bm{y}\in\mathcal{S}\atop\bm{y}\not=\bm{x}}\frac{\langle\bm{v},\bm{y}-\bm{x}\rangle}{\|\bm{y}-\bm{x}\|_{F}}\leq
0\right\\}$ $\displaystyle\
\quad\mbox{and}\quad\partial\delta_{\mathcal{S}}(\bm{x})=\mathcal{N}_{\mathcal{S}}(\bm{x}),\forall\bm{x}\in\mathcal{S},$

where $\mathcal{N}_{\mathcal{S}}(\bm{x})$ is the _normal cone_ to
$\mathcal{S}$ at $\bm{x}$.

###  I-C Organization

The rest of this paper is organized as follows. In Section II, we introduce
the proposed PALMe method for solving Problem (1). Then, we prove the
sufficient decrease property and relative error property in Section III-A,
estimate the KŁ exponent of Problem (1) in Section III-B, and prove the main
theorem of this work in Section III-C. In Section IV, we conduct numerical
experiments to validate our theoretical results and compare the proposed
approach to the existing approaches. We end with some concluding remarks in
Section V.

##  II Algorithm Design

We begin by reformulating Problem (1) as a two-block form. Noting that
$|x|=\max\\{x,-x\\}$ for any $x\in\mathbb{R}$, we can reformulate Problem (1)
as

$\displaystyle\min\ $ $\displaystyle
H(\bm{P},\bm{Q}):=-\langle\bm{P},\bm{X}^{T}\bm{Q}\bm{Q}^{T}\rangle$ (10)
$\displaystyle\mathrm{s.t.}\ $ $\displaystyle\bm{P}\in{\mathcal{B}}(n,d),\
\bm{Q}\in{\rm St}(d,K),$ (11)

where $\langle\bm{A},\bm{B}\rangle={\rm tr}(\bm{A}^{T}\bm{B})$ denotes the
Euclidean inner product of two matrices $\bm{A},\bm{B}$ of the same size and
${\mathcal{B}}(n,d)=\left\\{\bm{P}\in\mathbb{R}^{n\times d}:p_{ij}\in\\{\pm
1\\},\ i=1,\ldots,n,\ j=1,\ldots,d\right\\}$ denotes the set of all $n\times
d$ matrices with elements $\pm 1$. Observing that this problem has two
separate blocks of variables $\bm{P}$ and $\bm{Q}$, one can apply the PALM
method (see, e.g., [35, 32, 36]) for solving it. Specifically, given the
current iterate
$\left(\bm{P}^{k},\bm{Q}^{k}\right)\in\mathcal{B}(n,d)\times{\rm St}(d,K)$,
the method generates the next iterate
$\left(\bm{P}^{k+1},\bm{Q}^{k+1}\right)\in{\mathcal{B}}(n,d)\times{\rm
St}(d,K)$ via

$\bm{P}^{k+1}\in\operatorname*{argmin}\left\\{\begin{split}&\
-\langle\bm{P},\bm{X}^{T}\bm{Q}^{k}\bm{Q}^{k^{T}}\rangle\\\ &\
+\frac{\alpha_{k}}{2}\|\bm{P}-\bm{P}^{k}\|^{2}_{F}:\
\bm{P}\in\mathcal{B}(n,d)\end{split}\right\\},$ (12a)
$\bm{Q}^{k+1}\in\operatorname*{argmin}\left\\{\begin{split}&\
-\langle\bm{Q},(\bm{X}\bm{P}^{k+1}+\bm{P}^{{k+1}^{T}}\bm{X}^{T})\bm{Q}^{k}\rangle\\\
&\ +\frac{\beta_{k}}{2}\|\bm{Q}-\bm{Q}^{k}\|^{2}_{F}:\ \bm{Q}\in{\rm
St}(d,K)\end{split}\right\\},$ (12b)

where $\alpha_{k},\beta_{k}>0$ are step-size parameters. Recently, many
different extrapolation techniques have been successfully applied to
accelerate proximal algorithms for convex and non-convex optimization
problems; see, e.g., [37, 38, 29, 30, 31, 14]. This motivates us to
incorporate an extrapolation step into the update of $\bm{P}$ to achieve
possible acceleration of the PALM iterations. Specifically, we replace (12a)
with

$\displaystyle\
\bm{E}^{k}=\bm{Q}^{k}\bm{Q}^{k^{T}}+\gamma_{k}\left(\bm{Q}^{k}\bm{Q}^{k^{T}}-\bm{Q}^{k-1}\bm{Q}^{{k-1}^{T}}\right),$
(13) $\displaystyle\
\bm{P}^{k+1}\in\operatorname*{argmin}\left\\{\begin{split}&\
-\langle\bm{P},\bm{X}^{T}\bm{E}^{k}\rangle\\\ &\
+\frac{\alpha_{k}}{2}\|\bm{P}-\bm{P}^{k}\|^{2}_{F}:\
\bm{P}\in{\mathcal{B}}(n,d)\end{split}\right\\},$

where $\bm{E}^{k}\in\mathbb{R}^{d\times K}$ is the point extrapolated from
$\bm{Q}^{k}\bm{Q}^{k^{T}}$ and $\bm{Q}^{k-1}\bm{Q}^{{k-1}^{T}}$ and
$\gamma_{k}\in[0,1)$ is the parameter controlling extrapolation step-size. It
is worth noting that the above iterates (12b) and (13) admit closed-form
solutions. Specifically, the update (12a) in its most simplified form reads as

$\bm{P}^{k+1}\in\operatorname{sgn}\left(\bm{P}^{k}+\bm{X}^{T}\bm{E}^{k}/\alpha_{k}\right).$
(14)

On the other hand, the update (12b) is essentially an instance of the
orthogonal Procrustes problem [39], whose solution is given by

$\bm{Q}^{k+1}=\bm{U}^{k+1}\bm{V}^{{k+1}^{T}},$

where $\bm{U}^{k+1}\in\mathrm{St}(d,K)$ and $\bm{V}^{k+1}\in\mathrm{St}(K,K)$
are obtained by a thin SVD
$\bm{U}^{k+1}\bm{\Sigma}^{k+1}\bm{V}^{{k+1}^{T}}=\bm{Q}^{k}+(\bm{X}\bm{P}^{k+1}\bm{Q}^{k}+{\bm{P}^{k+1}}^{T}\bm{X}^{T}\bm{Q}^{k})/\beta_{k}$.
Now, we summarize the proposed method in Algorithm 1. One can verify that the
per-iteration cost of the proposed method is $\mathcal{O}(ndK+dK^{2})$, which
is cheap when $K\ll\min\left\\{n,d\right\\}$.

Input: $\bm{X}\in\mathbb{R}^{d\times n},\ \bm{P}^{0}\in{\mathcal{B}}(n,d),\
\bm{Q}^{-1}=\bm{Q}^{0}\in{\rm St}(d,K)$.

1 for  _$k=0,1,2,\cdots$_ do

2 choose step sizes $\alpha_{k},\beta_{k}>0$ and extraploation parameter
$\gamma_{k}\in[0,1]$;

3 set
$\bm{E}^{k}=\bm{Q}^{k}\bm{Q}^{k^{T}}+\gamma_{k}\left(\bm{Q}^{k}\bm{Q}^{k^{T}}-\bm{Q}^{k-1}\bm{Q}^{{k-1}^{T}}\right)$;

4 pick
$\bm{P}^{k+1}\in\operatorname{sgn}(\bm{P}^{k}+\bm{X}^{T}\bm{E}^{k}/\alpha_{k})$;

5 compute a thin SVD
$\bm{U}^{k+1}\bm{\Sigma}^{K+1}\bm{V}^{{k+1}^{T}}=\bm{Q}^{k}+\left(\bm{X}\bm{P}^{k+1}\bm{Q}^{k}+\bm{P}^{{k+1}^{T}}\bm{X}^{T}\bm{Q}^{k}\right)/\beta_{k}$;

6 Set $\bm{Q}^{k+1}=\bm{U}^{k+1}\bm{V}^{{k+1}^{T}}$;

7 Terminate if stopping criteria are met;

8 end for

Algorithm 1 PALM with a quadratic extrapolation (PALMe) for L1-PCA

We should emphasize that our algorithm has some essential differences with
other block coordinate descent-type methods. First, most of the methods
perform an extrapolation step on each block (see, e.g., [40, 41, 42, 43]),
while PALMe only takes the extrapolation on one block $\bm{Q}$. Experimental
results in Section IV validate the effectiveness of the latter approach for
Problem (1). Second, although the update (13) seems similar to that in [14],
our quadratic extrapolation step is totally new compared to the linear
extrapolation step in [14]. As a result, the update of $\bm{Q}$ is completely
different. To the best of our knowledge, our work is the first one to design
and analyze the quadratic extrapolation step for non-smooth and non-convex
problems.

##  III Convergence Analysis

Our goal in this section is to establish the convergence result of Algorithm
1. Generally, it is not easy to analyze the convergence behavior of algorithms
for solving a non-smooth and non-convex problem. Fortunately, Attouch et al.
[35] developed a unified framework that sheds light on analyzing the
convergence behavior of proximal algorithms for solving non-smooth and non-
convex problems. This framework has been widely used in the literature; see,
e.g., [14, 44, 31, 45]. Specifically, one needs to verify that the sequence of
iterates generated by the considered algorithm satisfies the _sufficient
decrease property_ and the _relative error property_ in terms of a proper
potential function. Moreover, one needs to verify that the considered
potential function satisfies the _KŁ property with the associated exponent_.
Note that the first two properties are algorithm-dependent, and we can
directly verify them by analyzing the updates of Algorithm 1. By contrast, the
last one depends on the function itself rather than the algorithm according to
the definition of the KŁ property. To simplify our development, we define

$\displaystyle\ell(\bm{Q})=-\|\bm{Q}\bm{Q}^{T}\bm{X}\|_{1}+\delta_{{\rm
St}(d,K)}(\bm{Q}),$ (15)

and

$\displaystyle
h(\bm{P},\bm{Q})=H(\bm{P},\bm{Q})+\delta_{\mathcal{B}(n,K)}(\bm{P})+\delta_{{\rm
St}(d,K)}(\bm{Q}).$ (16)

Moreover, we can define the critical points of $\ell$ as follows.

######  Definition 2 (Critical point of $\ell$).

We say that $\bm{Q}\in{\rm St}(d,K)$ is a critical point of $\ell$ if

$\displaystyle\mathbf{0}\in$
$\displaystyle-\bm{X}\operatorname{sgn}(\bm{X}^{T}\bm{Q}\bm{Q}^{T})\bm{Q}-\operatorname{sgn}(\bm{Q}\bm{Q}^{T}\bm{X})\bm{X}^{T}\bm{Q}$
$\displaystyle+\mathcal{N}_{{\rm St}(d,K)}(\bm{Q}).$ (17)

Invoking the subdifferential calculus rules in [34, Chapter 10B], we remark
that this is a necessary condition for local optimality. Besides, it should be
noted that every limiting critical point of $\ell$ is a critical point of
$\ell$, but the converse is not known to hold.

###  III-A Sufficient Decrease and Relative Error Properties

To prove the sufficient decrease and relative error properties, we need to
choose a proper potential function. One immediate choice of the potential
function is the function $h$ defined in (16). However, due to the
extrapolation step in Algorithm 1, it is not clear how to verify the mentioned
properties in terms of $h$. Motivated by the potential functions constructed
in [29, 41, 14], we consider the potential function
$\Phi_{\beta}(\bm{P},\bm{Q},\bm{Q}^{\prime}):\mathbb{R}^{n\times
d}\times\mathbb{R}^{d\times K}\times\mathbb{R}^{d\times
K}\rightarrow(-\infty,+\infty]$ defined by

$\Phi_{\beta}(\bm{P},\bm{Q},\bm{Q}^{\prime})=h(\bm{P},\bm{Q})+\frac{\beta}{2}\|\bm{Q}-\bm{Q}^{\prime}\|^{2}_{F}.$
(18)

Then, when the step-size and the extrapolation parameters in Algorithm 1 are
suitably chosen, we can verify the desired properties in terms of
$\Phi_{\beta}$ as follows.

######  Lemma 1.

Let $\bm{C}^{k}=\left(\bm{P}^{k},\bm{Q}^{k},\bm{Q}^{k-1}\right)$ for all
$k\geq 0$. Suppose that the step-size and extrapolation parameters in
Algorithm 1 satisfy for all $k\geq 0$,

$\displaystyle\ \alpha_{*}\leq\alpha_{k}\leq\alpha^{*},\
\frac{3}{2}\beta_{*}+2\|\bm{X}\bm{P}^{k}\|\leq\beta_{k}\leq\beta^{*},$ (19)
$\displaystyle\
0\leq\gamma_{k}<\gamma^{*}=\min\left\\{1,\frac{\alpha_{*}\beta_{*}}{8\|\bm{X}\|^{2}}\right\\}$
(20)

for some $\alpha_{*},\alpha^{*},\beta_{*},\beta^{*}\in(0,\infty)$. Then, the
following statements hold:
(i) The sequence $\\{\bm{C}^{k}\\}_{k\geq 0}$ is bounded.
(ii) There exists a constant $\kappa_{1}>0$ such that for all $k\geq 0$,

$\displaystyle\Phi_{\beta_{*}}(\bm{C}^{k+1})-\Phi_{\beta_{*}}(\bm{C}^{k})\leq-\kappa_{1}\|{\bm{C}}^{k+1}-{\bm{C}}^{k}\|_{F}^{2}.$
(21)

(iii) There exists a constant $\kappa_{2}>0$ such that for all $k\geq 0$,

$\displaystyle\operatorname{dist}\left(\bm{0},\partial\Phi_{\beta_{*}}(\bm{C}^{k+1})\right)\leq\kappa_{2}\|{\bm{C}}^{k+1}-{\bm{C}}^{k}\|_{F}.$
(22)

We defer the detailed proof to Section -A of the appendix. Thanks to the
boundness of ${\mathcal{B}}(n,d)$ and ${\rm St}(d,K)$, one can easily verify
the boundness of the sequence $\\{\bm{C}^{k}\\}_{k\geq 0}$ as in (i). We refer
to (ii) as _sufficient decrease property_ , which follows from the updates in
(12b) and (13), and the Lipschitz smoothness of $H(\bm{P},\bm{Q})$ on a
bounded set. It ensures that the potential function
$\Phi_{\beta_{*}}(\bm{C}^{k})$ is monotonically decreasing. This, together
with (i), implies that the sequence $\\{\Phi_{\beta_{*}}(\bm{C}^{k})\\}_{k\geq
0}$ converges and $\|{\bm{C}}^{k+1}-{\bm{C}}^{k}\|_{F}$ goes to zero as $k$
increases. We refer to (iii) as _relative error property_ , which is also
known as safeguard in [46, 47]. It is proved by the optimality condition of
the updates in (12b) and (13). Armed with this lemma, we can conclude that
every accumulation point of $\\{{\bm{C}}^{k}\\}_{k\geq 0}$ is a critical point
of $\Phi_{\beta_{*}}$, which is essentially _subsequence convergence_. Remark
that compared to the proof in [14, Proposition 6], this proof is different and
more involved. Indeed, the objective function of Problem (10) is quadratic in
terms of $\bm{Q}$, while that is linear in [14]. Moreover, our extrapolation
step in (13) is non-linear, while that is linear in [14].

###  III-B KŁ Property

In this subsection, we verify the KŁ property of the constructed potential
function $\Phi_{\beta}$ and derive the corresponding KŁ exponent. For the
definition of the KŁ property, we refer the reader to [35, Definition 2.5] for
the details. The KŁ property is widely used for studying convergence behavior
of various first-order methods; see, e.g., [48, 35, 32]. It is known that the
KŁ property holds for many extended-valued lower-semicontinuous functions,
such as semi-algebraic functions. According to [48, Section 4.3] and [32,
Theorem 3], one can easily verify that $\Phi_{\beta}$ satisfies the KŁ
property. Here, we are more interested in the KŁ exponent of $\Phi_{\beta}$,
since it determines the convergence rate of the proposed method.

######  Definition 3 (KŁ exponent).

Suppose that $f:\mathbb{R}^{d}\rightarrow(-\infty,\infty]$ is proper and lower
semicontinuous. We say that the function $f$ has a _KŁ exponent_ of
$\theta\in[0,1)$ at the point $\bar{\bm{x}}\in\operatorname{dom}(\partial f)$
if there exist constants $\epsilon,\eta>0$, $\nu\in(0,+\infty]$ such that

$\operatorname{dist}(\bm{0},\partial
f(\bm{x}))\geq\eta(f(\bm{x})-f(\bar{\bm{x}}))^{\theta}$

whenever $\|\bm{x}-\bar{\bm{x}}\|\leq\epsilon$ and
$f(\bar{\bm{x}})<f(\bm{x})<f(\bar{\bm{x}})+\nu$.

Now, we are devoted to estimating the KŁ exponent of $\Phi_{\beta}$. Let
$\bm{P}_{1},\ldots,\bm{P}_{2^{nd}}$ be an enumeration of the elements in
$\mathcal{B}(n,d)$. By definition of the $\ell_{1}$-norm, we can express the
objective function $\ell$ of Problem (15) as the pointwise minimum of finitely
many proper and lower semicontinuous functions:

$\ell(\bm{Q})=\min_{i\in\\{1,\ldots,2^{nd}\\}}\big{\\{}\underbrace{\langle\bm{X}\bm{P}_{i},\bm{Q}\bm{Q}^{\top}\rangle+\delta_{{\rm
St}(d,K)}(\bm{Q})}_{\ell_{i}(\bm{Q})}\big{\\}}.$

Using the same argument in [14, Section 3.1], we can elucidate that estimating
the KŁ exponent of $\ell$ boils down determining the KŁ exponent of
$\ell_{1},\dots,\ell_{2^{nd}}$. Based on [44, Theorem 1] and [49, Lemma 2.1],
we can characterize the KŁ exponent of the quadratic optimization problem with
orthogonality constraint (QP-OC).

######  Fact 1.

For the QP-OC problem

$\displaystyle\min\\{\langle\bm{A},\bm{Q}\bm{B}\bm{Q}^{T}\rangle:\
\bm{Q}\in{\rm St}(d,K)\\},$

where $\bm{A}\in\mathbb{R}^{d\times d}$ and $\bm{B}\in\mathbb{R}^{K\times K}$
are symmetric, the KŁ exponent is $1/2$.

Using [49, Theorem 3.6], the KŁ exponent of $h$, and relations between $h$ and
$\Phi_{\beta_{*}}$, we can obtain the following corollary. The detailed proof
is provided in Section -B of the appendix.

######  Corollary 1.

The KŁ exponent of $\Phi_{\beta_{*}}$ is $1/2$.

###  III-C Linear Convergence of Algorithm 1 and Properties of Limit Points

Equipped with Lemma 1, Corollary 1, and [35, Theorem 2.9], we are ready to
prove that the sequence $\\{\bm{C}^{k}\\}_{k\geq 0}$ converges at least
linearly to a limiting critical point of $\Phi_{\beta_{*}}$. However, our
goals are to prove the convergence rate of Algorithm 1 and figure out under
what condition the limit point of the sequence $\\{\bm{Q}^{k}\\}$ is a
critical point of Problem (1). To tackle the issues, we need to characterize
the relationships among $\ell$, $h$, and $\Phi_{\beta}$.

######  Lemma 2.

Let $\beta>0$ be given and $\ell$, $h$, and $\Phi_{\beta}$ be respectively
defined in (15), (16), and (18). Then, the following statements hold:
(i) If $(\bm{P},\bm{Q})$ with
$\bm{P}\in\operatorname{sgn}(\bm{X}^{T}\bm{Q}\bm{Q}^{T})$ is a limiting
critical point of $h$, then $\bm{Q}\in{\rm St}(d,K)$ is a critical point of
$\ell$.
(ii) Suppose that
$(\bm{P},\bm{Q},\bm{Q}^{\prime})\in\mathcal{B}(n,d)\times{\rm
St}(d,K)\times{\rm St}(d,K)$ is a limiting critical point of $\Phi_{\beta}$.
Then, we have $\bm{Q}=\bm{Q}^{\prime}$. Moreover, $(\bm{P},\bm{Q},\bm{Q})$ is
a limiting critical point of $\Phi_{\beta}$ if and only if $(\bm{P},\bm{Q})$
is a limiting critical point of $h$.

We defer the detailed proof to Section -C of the appendix. Armed with this
lemma, 1, Corollary 1, and [35, Theorem 2.9], we can establish the main
theorem of this work, which provides convergence analysis of the proposed
method.

######  Theorem 1.

Let $\\{(\bm{P}^{k},\bm{Q}^{k})\\}_{k\geq 0}$ be the sequence of iterates
generated by Algorithm (1) under the conditions (i)
$\alpha_{*}\leq\alpha_{k}\leq\alpha^{*}$ for some
$\alpha^{*},\alpha_{*}\in(0,+\infty)$, (ii)
$3\beta_{*}/2+2\|\bm{X}\bm{P}^{k}\|\leq\beta_{k}\leq\beta^{*}$ for some
$\beta^{*},\beta_{*}\in(0,+\infty)$, and (iii) $\gamma_{k}$ satisfy
$0\leq\gamma_{k}<\gamma^{*}=\min\left\\{1,\alpha_{*}\beta_{*}/({8\|\bm{X}\|^{2}})\right\\}$.
Then, the sequence $\\{(\bm{P}^{k},\bm{Q}^{k})\\}_{k\geq 0}$ converges at
least linearly to a limiting critical point $(\bm{P}^{*},\bm{Q}^{*})$ of
Problem (10). Moreover, if $\alpha_{k}=\alpha_{*}$ for all $k\geq 0$ such that

$0<\alpha_{*}<\min\left\\{\begin{split}&\
|(\bm{X}^{T}\bm{Q}^{*}\bm{Q}^{*\top})_{ij}|:(\bm{X}^{T}\bm{Q}^{*}\bm{Q}^{*T})_{ij}\not=0,\\\
&\ \quad i=1,\dots,n,\,j=1,\dots,d\end{split}\right\\},$ (23)

then $\bm{Q}^{*}$ is a critical point of Problem (1).

###### Proof.

Let $\bm{C}^{k}=\left(\bm{P}^{k},\bm{Q}^{k},\bm{Q}^{k-1}\right)$ for all
$k\geq 0$. It follows from Lemma 1, Corollary 1, and [32, Theorem 1] that the
sequence $\\{\bm{C}^{k}\\}_{k\geq 0}$ converges at least linearly to a
limiting critical point $\left(\bm{P}^{*},\bm{Q}^{*},\bm{Q}^{*}\right)$ of
$\Phi_{\beta_{*}}$. This, together with (ii) of Lemma 2, gives that the
sequence $\\{(\bm{P}^{k},\bm{Q}^{k})\\}_{k\geq 0}$ converges at least linearly
to the limiting critical point $(\bm{P}^{*},\bm{Q}^{*})$ of $h$.

Suppose that $\alpha_{k}=\alpha_{*}$ for all $k\geq 0$. According to (13) and
(14), we have
$\bm{E}^{k}=\bm{Q}^{k}\bm{Q}^{k^{T}}+\gamma_{k}(\bm{Q}^{k}\bm{Q}^{k^{T}}-\bm{Q}^{k-1}\bm{Q}^{{k-1}^{T}})$
and
$\bm{P}^{k+1}\in\operatorname{sgn}(\bm{P}^{k}+\bm{X}^{T}\bm{E}^{k}/\alpha_{k})$.
This, together with
$(\bm{P}^{k},\bm{Q}^{k})\rightarrow(\bm{P}^{*},\bm{Q}^{*})$ and
$\alpha_{k}=\alpha_{*}$, yields that

$\displaystyle\bm{P}^{*}\in\operatorname{sgn}(\bm{P}^{*}+\bm{X}^{T}\bm{Q}^{*}\bm{Q}^{*^{T}}/\alpha_{*}).$

Combining this with the result that $(\bm{P}^{*},\bm{Q}^{*})$ is a limiting
critical point of $h$ gives

$\begin{split}\mathbf{0}\in&\
-\bm{X}\operatorname{sgn}(\bm{P}^{*}+\bm{X}^{T}\bm{Q}^{*}\bm{Q}^{*^{T}}/\alpha_{*})\bm{Q}^{*}\\\
&\
-\operatorname{sgn}(\bm{P}^{*^{T}}+\bm{Q}^{*}\bm{Q}^{*^{T}}\bm{X}^{T}/\alpha_{*})\bm{X}\bm{Q}^{*}+\mathcal{N}_{{\rm
St}(d,K)}(\bm{Q}^{*}).\end{split}$

In particular, noting that $\bm{P}^{*}\in\mathcal{B}(n,d)$ and $\alpha_{*}$
satisfies (23), we have
$\operatorname{sgn}(\bm{P}^{*}+\bm{X}^{T}\bm{Q}^{*}\bm{Q}^{*^{T}}/\alpha_{*})\subseteq\operatorname{sgn}(\bm{X}^{T}\bm{Q}^{*}\bm{Q}^{*^{T}})$.
These further implies

$\begin{split}\mathbf{0}\in&\
-\bm{X}\operatorname{sgn}(\bm{X}^{T}\bm{Q}^{*}\bm{Q}^{*^{T}})\bm{Q}^{*}-\operatorname{sgn}(\bm{Q}^{*}\bm{Q}^{*^{T}}\bm{X}^{T})\bm{X}^{T}\bm{Q}^{*}\\\
&\ +\mathcal{N}_{{\rm St}(d,K)}(\bm{Q}^{*}).\end{split}$

Then, we complete the proof.  $\sqcup\hbox to0.0pt{\hss$\sqcap$}$

The first part of this theorem demonstrates that under some conditions imposed
on step-size and extrapolation parameters, the iterates generated by Algorithm
1 linearly converges to a limiting critical point of the reformulated L1-PCA
problem (10). The conditions imposed on $\alpha_{k},\beta_{k}$ and
$\gamma_{k}$ are introduced by proving the sufficient decreasing of objective
function value. Interestingly, we found that in experiments, these parameters
can be chosen flexibly even if they may violate aforementioned conditions. We
leave this theory-practice gap for further research. The second part of
Theorem 1 shows that the condition (23) is sufficient to guarantee that the
point returned by Algorithm 1 is not only a limiting critical point of Problem
(10), but also a critical point of Problem (1). It is worth noting that we can
efficiently verify whether this condition holds after obtaining a limit point
$\bm{Q}^{*}$.

##  IV Experiment results

In this section, we conduct numerical experiments on synthetic and real
datasets to demonstrate the efficacy of our studied approach. We first compare
the proposed method with other existing ones for solving Problem (1) in terms
of convergence performance and solution quality in Section IV-A and clustering
accuracy in Section IV-B. In Section IV-C, we then compare the image
reconstruction performance of our studied approach with the surveyed other
L1-PCA approaches in Section I. All the experiments are performed on a PC
running Windows 10 with an Intel® Core™ i5-8600 3.10GHz CPU and 16GB memory.
Our codes are implemented in MATLAB R2021a and can be found at
https://github.com/TaoliZheng/RIL1PCA.

TABLE I: Step-size parameters of the tested methods $(\alpha,\beta)$ | PALMe | PALM | iPALM | GiPALM | pDCAe
---|---|---|---|---|---
synthetic $(n,d)=(5000,1000)$ | $(10^{-7},100)$ | $(10^{-7},10)$ | $(10^{-7},10)$ | $(10^{-7},100)$ | $(-,0.05)$
synthetic $(n,d)=(1000,5000)$ | $(10^{-6},1)$ | $(10^{-6},10)$ | $(10^{-6},0.1)$ | $(10^{-6},1)$ | $(-,0.1)$
_colon-cancer_ $(n,d)=(62,2000)$ | $(10^{-10},100)$ | $(10^{-10},500)$ | $(10^{-10},500)$ | $(10^{-10},100)$ | $(-,100)$
TABLE II: TEV of the tested methods $(\alpha,\beta)$ | PALMe | PALM | iPALM | GiPALM | pDCAe
---|---|---|---|---|---
synthetic $(n,d)=(5000,1000)$ | 0.978176 | 0.973894 | 0.973955 | 0.974290 | 0.973938
synthetic $(n,d)=(1000,5000)$ | 0.955969 | 0.940970 | 0.940904 | 0.943625 | 0.941027
_colon-cancer_ $(n,d)=(62,2000)$ | 0.925389 | 0.921079 | 0.921001 | 0.928077 | 0.924373

(a) synthetic: $(n,d)=(1000,5000)$

(b) synthetic: $(n,d)=(5000,1000)$

(c) real _colon-cancer_ : $(n,d)=(62,2000)$

Figure 1: Convergence performance of function values: The $x$-axis is number
of iterations, and the $y$-axis is function value gap
$h(\bm{P}^{k},\bm{Q}^{k})-h(\bm{P}^{*},\bm{Q}^{*})$, where
$(\bm{P}^{*},\bm{Q}^{*})$ is the last iterate of the tested method.

###  IV-A Convergence Performance and Solution Quality

We first examine the convergence performance and solution quality of the
proposed method on synthetic and real datasets. We also compare it with some
existing methods that can be applied to solve Problem (1) or its reformulation
(10), which includes the standard PALM in [32], the inertial proximal
alternating linearized minimization (iPALM) method in [41], the Gauss-Seidel-
type iPALM (GiPALM) method in [42], and the proximal difference-of-convex with
extrapolation (pDCAe) method in [29]. As demonstrated in [50, 14], we can use
the following measure named _total explained variation_ (TEV) to compare the
quality of solutions returned by the tested methods:

$\displaystyle{\rm
TEV}=\frac{\|\bm{X}^{T}\bm{Q}\|_{F}^{2}}{\|\bm{X}^{T}\bar{\bm{Q}}\|_{F}^{2}},$

where $\bm{Q}\in\mathbb{R}^{d\times K}$ is the solution returned by the tested
method and $\bar{\bm{Q}}\in\mathbb{R}^{d\times K}$ is the matrix formed by the
eigenvectors associated with the leading $K$ eigenvalues of
$\bm{X}\bm{X}^{T}$. In general, it is assumed that the larger TEV leads to the
better solution.

In the tests, we employ the fixed effect model in [12] to generate the
synthetic data. We use the same approach as that in [14] to generate a
$K$-dimensional subspace and a Laplacian noise with mean zero and variance
$\sigma^{2}$. As a result, the generated data points lie near the subspace
perturbed by the noise. We set $\sigma=0.5,K=50$ and generate two synthetic
datasets with dimensions $\left(n,d\right)=\left(1000,5000\right)$ and
$\left(n,d\right)=\left(5000,1000\right)$, respectively. As for the real
dataset, we use the dataset _colon-cancer_ downloaded from LIBSVM
[51]111https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/ with
dimensions $\left(n,d\right)=\left(62,2000\right)$ and set $K=20$.

(a) synthetic: $(n,d)=(1000,5000)$

(b) synthetic: $(n,d)=(5000,1000)$

(c) _colon-cancer_ : $(n,d)=(62,2000)$

Figure 2: Convergence performance of iterates: The $x$-axis is number of
iterations, and the $y$-axis is iterate gap $\|\bm{Q}^{k}-\bm{Q}^{*}\|$, where
$\bm{Q}^{*}$ is the last iterate of the tested method.

We then present the setting of parameters for different algorithms in the
tests. To be fair, all the parameters are tuned based on numerous runs to
balance the solution quality and running time. We first specify the step-size
parameters of the tested methods. We set
$(\alpha_{k},\beta_{k})=(\alpha,\beta)$ for all $k\geq 0$, whose values are
provided in Table I. Note that it is not required to choose $\alpha_{k}$ for
pDCAe since it only has one block of variables. We next specify the
extrapolation parameters of the tested methods. For PALMe, we set it as $1$.
We observe that it works surprisingly well, even though such a choice could
violate the condition in Theorem 1. For pDCAe, we set the extrapolation
parameter as 0.2. For iPALM, we set the extrapolation parameters for updating
the block variables $\bm{P}$ and $\bm{Q}$ both as 0.2. For GiPALM, these two
parameters are respectively set as $1/2$ and $1/4$ for all $k\geq 0$. In each
test, we employ the same starting point for the tested algorithms. For each
algorithm, we terminate it when the Frobenious norm of the difference between
two consecutive iterates is less than $10^{-6}$.

To compare the convergence performance of the tested algorithms, we plot the
the function value gap $h(\bm{P}^{k},\bm{Q}^{k})-h(\bm{P}^{*},\bm{Q}^{*})$ and
iterates gap $\|\bm{Q}^{k}-\bm{Q}^{*}\|_{F}$ against the iteration number for
all the tested algorithms in Figs. 1 and 2, respectively. For each tested
algorithm, we choose its last iterate as $(\bm{P}^{*},\bm{Q}^{*})$. According
to these figures, we can observe that the sequences of function value gaps and
iterate gaps converge linearly, which supports our theoretical result in
Theorem 1. Moreover, it can be observed that PALMe enjoys a substantially
faster linear convergence compared to the other methods. To compare the
solution quality, we report the value of TEV of the tested methods in Table II
averaged over 10 runs. We can observe that the solution quality of PALMe is
comparable to those of the other methods.

###  IV-B Clustering Accuracy

In this subsection, we compare the clustering performance of our proposed
method with PALM, iPALM, GiPALM, and pDCAe, which are introduced in Section
IV-A. A good way of evaluating the clustering performance of a algorithm for
solving Problem (1) is to study its clustering accuracy when applied to
clustering on a subspace; see, e.g., [13, 14]. Specifically, we obtain a label
vector by applying the _k-means_ clustering to the projections of data points
onto the subspace returned by a algorithm for solving Problem (1) and compute
the clustering accuracy by comparing the returned label and the true label. In
our experiments, we use the real-world datasets _a6a_ , a9a, gisette, and
ijcnn1 downloaded from LIBSVM [51], whose dimensions can be found in Table
III. The dimension $K$ of the subspace used by L1-PCA is chosen such that
$\sum_{k=1}^{K}\sigma^{2}_{k}\geq 0.8\sum_{k=1}^{p}\sigma^{2}_{k}$, where
$p=\min\\{n,d\\}$ and $\sigma_{1}\geq\dots\geq\sigma_{p}\geq 0$ are singular
values of $\bm{X}$. We list the value of $K$ for each dataset in Table III.
When the difference between two consecutive iterates is less than $10^{-6}$,
we stop the iterations.

In the tests, we set the step-size and extrapolation parameters of the tested
algorithms as in Table III. The fixed step-size is still used in experiments.
To better evaluate the clustering performance, each algorithm is tested for 10
times. Then, we plot the clustering accuracy and the running time for the
tested algorithms in Fig. 3 averaged over 10 runs. It is worth noting that the
point in the top-left corner of the axes converges faster and achieves higher
clustering accuracy than that in the bottom-right corner. Then, we can observe
from Fig. 3 that PALMe can generally achieve a comparable clustering accuracy
to those of the other methods with less running time, which demonstrates its
efficacy.

(a) a6a

(b) a9a

(c) gisette

(d) ijcnn

Figure 3: Clustering accuracy and running time on LIBSVM datasets TABLE III:
Dimension of datasets, dimension of subspace, step-size and extrapolation
parameters of the tested methods

$(\alpha,\beta,\gamma_{P},\gamma_{Q})$ | PALMe | PALM | iPALM | GiPALM | pDCAe
---|---|---|---|---|---
|  _a6a_ dataset
---
$(n,d,K)=(11220,122,6)$
$(10^{-7},10^{5},0,1)$ | $(10^{-8},10^{5},0,0)$ | $(10^{-7},10^{5},0.2,0.3)$ | $(10^{-6},5*10^{4},0.5,0.25)$ | $(-,5*10^{-6},-,1)$
|  _a9a_ dataset
---
$(n,d,K)=(32561,123,6)$
$(10^{-8},10^{6},0,0.5)$ | $(10^{-7},10^{6},0,0)$ | $(10^{-8},10^{6},0.2,0.2)$ | $(10^{-8},10^{6},0.3,0.4)$ | $(-,10^{-6},-,0.3)$
|  _gisette_ dataset
---
$(n,d,K)=(6000,5000,1)$
$(10^{-5},5000,0,1)$ | $(10^{-6},1000,0,0)$ | $(10^{-6},10^{4},0.05,0.05)$ | $(10^{-5},10^{4},0.05,0.05)$ | $(-,100,-,0.2)$
|  _ijcnn1_ dataset
---
$(n,d,K)=(49990,22,8)$
$(10^{-7},10^{4},0,0.1)$ | $(10^{-7},10^{4},0,0)$ | $(10^{-6},1.5*10^{4},0.05,0.05)$ | $(10^{-6},10^{4},0.2,0.2)$ | $(-,10^{-4},-,0.2)$

###  IV-C Image Reconstruction

As suggested in [52, 15], the robustness to outliers can be visualized by
image reconstruction. In this section, we conduct image reconstruction
experiments on 6 gray-scale human face images (see the first row of Fig. 4)
downloaded from _AT &T_ Database of
Faces222https://www.kaggle.com/datasets/kasikrit/att-database-of-
faces?resource=download to test robustness of different models. We assume that
the clean image $\bm{A}$ is not available and instead we have 9 corrupted
versions $\bm{A}_{1},\dots,\bm{A}_{9}$. Here, these corrupted instances are
generated by adding outliers to a block as follows. We firstly partition the
original image $\bm{A}$ into 9 equal-sized blocks, say
$\bm{B}_{1},\dots,\bm{B}_{9}$. Next, we select one block $\bm{B}_{i}$, divide
it into 4 smaller blocks of equal size, add outliers drawn from the discrete
uniform distribution in range $[1,200]$ to its diagonal blocks, and normalize
it into range $[0,255]$ to generate $\bm{A}_{i}$ for $i=1,\dots,9$. One
instance of the corrupted image is shown in the second row of Fig. 4. Then, we
can vectorize these 9 corrupted images by stacking them to form the data
matrix $\bm{X}$, i.e.,

$\displaystyle\bm{X}=\begin{bmatrix}\mathrm{vec}(\bm{A}_{1})&\dots&\mathrm{vec}(\bm{A}_{9}).\end{bmatrix}$

Now, we apply various PCA models to do image reconstruction on $\bm{X}$,
including our proposed approach, L1-PCA-based approach in [14], R1-PCA-based
approach in [13], matrix factorization based L1-PCA (MF-L1PCA) approach in
[26], and L12-PCA-based approach in [15]. We simply denote our approach by
RI-L1PCA. The projection of $\bm{X}$ onto the learned subspace is referred as
the reconstructed image, i.e, $\hat{\bm{X}}=\bm{Q}\bm{Q}^{T}\bm{X}$, where
$\bm{Q}$ is the returned subspace of the tested approach.

In the experiments, the maximal iteration is set as 1000 and $K=2$. For each
algorithm, we stop it when the Frobenious norm of the difference between two
consecutive iterates is less then $10^{-3}$. Then, we can obtain the
reconstructed images learned by different approaches in Fig. 4. We can observe
that the quality of the reconstructed images returned by our approach has
comparable performance with others.

Figure 4: Image reconstruction of different models

##  V Conclusions

In this paper, we proposed a proximal alternating linearzed minimization
method with a non-linear extrapolation to solve the rotationally invariant
L1-PCA problem. To the best of our knowledge, this is the first work that
presents a fast iterative method to tackle this non-smooth non-convex problem.
We proved that our method converges at least linearly to a critical point of
the L1-PCA problem under a mild condition. To demonstrate the efficacy of our
method, extensive experimental results on both synthetic and real datasets are
reported. In particular, we compare our approach with some existing robust PCA
approaches to elucidate its robustness to outliers. Motivated by our proposed
non-linear extrapolation, one future direction is to extend such extrapolation
scheme to more general function class. Another possible direction is to design
other optimization algorithms to tackle the rotationally invariant L1-PCA
problem and analyze its convergence behavior.

## References

  * [1] I. T. Jolliffe and J. Cadima, “Principal component analysis: a review and recent developments,” _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , vol. 374, no. 2065, p. 20150202, 2016.
  * [2] J. Yang, D. Zhang, A. F. Frangi, and J.-y. Yang, “Two-dimensional pca: a new approach to appearance-based face representation and recognition,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 26, no. 1, pp. 131–137, 2004.
  * [3] H. Moon and P. J. Phillips, “Computational and performance aspects of pca-based face-recognition algorithms,” _Perception_ , vol. 30, no. 3, pp. 303–321, 2001.
  * [4] M. O. Faruqe and M. A. M. Hasan, “Face recognition using pca and svm,” in _2009 3rd International Conference on Anti-counterfeiting, Security, and Identification in Communication_. IEEE, 2009, pp. 97–101.
  * [5] Y. Ke and R. Sukthankar, “Pca-sift: A more distinctive representation for local image descriptors,” in _Proceedings of the 2004 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2004. CVPR 2004._ , vol. 2. IEEE, 2004, pp. II–II.
  * [6] T. Bouwmans, S. Javed, H. Zhang, Z. Lin, and R. Otazo, “On the applications of robust pca in image and video processing,” _Proceedings of the IEEE_ , vol. 106, no. 8, pp. 1427–1457, 2018.
  * [7] A. Lakhina, M. Crovella, and C. Diot, “Diagnosing network-wide traffic anomalies,” _ACM SIGCOMM computer communication review_ , vol. 34, no. 4, pp. 219–230, 2004.
  * [8] J. Lei and A. Rinaldo, “Consistency of spectral clustering in stochastic block models,” _The Annals of Statistics_ , vol. 43, no. 1, pp. 215–237, 2015\.
  * [9] V. Barnett and T. Lewis, “Outliers in statistical data,” _Wiley Series in Probability and Mathematical Statistics. Applied Probability and Statistics_ , 1984.
  * [10] S. J. Devlin, R. Gnanadesikan, and J. R. Kettenring, “Robust estimation of dispersion matrices and principal components,” _Journal of the American Statistical Association_ , vol. 76, no. 374, pp. 354–362, 1981.
  * [11] G. Lerman and T. Maunu, “An overview of robust subspace recovery,” _Proceedings of the IEEE_ , vol. 106, no. 8, pp. 1380–1410, 2018.
  * [12] A. Baccini, P. Besse, and A. de Falguerolles, “A $L_{1}$-norm PCA and a heuristic approach,” in _Ordinal and Symbolic Data Analysis_ , E. Diday, Y. Lechevallier, and O. Opitz, Eds. Berlin: Springer-Verlag, 1996, pp. 359–368.
  * [13] C. Ding, D. Zhou, X. He, and H. Zha, “R1-PCA: rotational invariant L1-norm principal component analysis for robust subspace factorization,” in _Proceedings of the 23rd international conference on Machine learning_ , 2006, pp. 281–288.
  * [14] P. Wang, H. Liu, and A. M.-C. So, “Linear convergence of a proximal alternating minimization method with extrapolation for $ell\\_1$-norm principal component analysis,” _arXiv preprint arXiv:2107.07107_ , 2021.
  * [15] F. Nie, L. Tian, H. Huang, and C. Ding, “Non-greedy l21-norm maximization for principal component analysis,” _IEEE Transactions on Image Processing_ , vol. 30, pp. 5277–5286, 2021.
  * [16] Y. W. Park and D. Klabjan, “Iteratively reweighted least squares algorithms for l1-norm principal component analysis,” in _2016 IEEE 16th International Conference on Data Mining (ICDM)_. IEEE, 2016, pp. 430–438.
  * [17] Q. Ke and T. Kanade, “Robust $L_{1}$ norm factorization in the presence of outliers and missing data by alternative convex programming,” in _Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on_ , vol. 1. IEEE, 2005, pp. 739–746.
  * [18] J. P. Brooks, J. H. Dulá, and E. L. Boone, “A pure l1-norm principal component analysis,” _Computational statistics & data analysis_, vol. 61, pp. 83–98, 2013.
  * [19] N. Kwak, “Principal component analysis based on l1-norm maximization,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 30, no. 9, pp. 1672–1680, 2008.
  * [20] F. Nie, H. Huang, C. Ding, D. Luo, and H. Wang, “Robust principal component analysis with non-greedy $\ell_{1}$-norm maximization,” in _Proceedings of the 22nd International Joint Conference on Artificial Intelligence (IJCAI 2011)_ , vol. 2, 2011, pp. 1433–1438.
  * [21] P. P. Markopoulos, S. Kundu, S. Chamadia, and D. A. Pados, “Efficient L1-norm principal-component analysis via bit flipping,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 16, pp. 4252–4264, 2017.
  * [22] P. Wang, H. Liu, and A. M.-C. So, “Globally convergent accelerated proximal alternating maximization method for l1-principal component analysis,” in _ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2019, pp. 8147–8151.
  * [23] P. P. Markopoulos, S. Kundu, S. Chamadia, N. Tsagkarakis, and D. A. Pados, “Outlier-resistant data processing with L1-norm principal component analysis,” in _Advances in Principal Component Analysis_. Springer, 2018, pp. 121–135.
  * [24] B. Minnehan and A. Savakis, “Grassmann manifold optimization for fast $l\\_1$-norm principal component analysis,” _IEEE Signal Processing Letters_ , vol. 26, no. 2, pp. 242–246, 2018.
  * [25] K. Liu and Y. Cao, “Robust principal component analysis: A construction error minimization perspective,” _arXiv preprint arXiv:2111.12132_ , 2021.
  * [26] L. Yu, M. Zhang, and C. Ding, “An efficient algorithm for L1-norm principal component analysis,” in _2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2012, pp. 1377–1380.
  * [27] A. Eriksson and A. Van Den Hengel, “Efficient computation of robust weighted low-rank matrix approximations using the l_1 norm,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 34, no. 9, pp. 1681–1690, 2012.
  * [28] Q. Wang, Q. Gao, X. Gao, and F. Nie, “$\ell_{2,p}$-norm based PCA for image recognition,” _IEEE Transactions on Image Processing_ , vol. 27, no. 3, pp. 1336–1346, 2017.
  * [29] B. Wen, X. Chen, and T. K. Pong, “A proximal difference-of-convex algorithm with extrapolation,” _Computational optimization and applications_ , vol. 69, no. 2, pp. 297–324, 2018.
  * [30] H. Li and Z. Lin, “Accelerated proximal gradient methods for nonconvex programming,” in _Advances in neural information processing systems_ , 2015, pp. 379–387.
  * [31] Z. Lu, Z. Zhou, and Z. Sun, “Enhanced proximal DC algorithms with extrapolation for a class of structured nonsmooth DC minimization,” _Mathematical Programming_ , vol. 176, no. 1-2, pp. 369–401, 2019.
  * [32] J. Bolte, S. Sabach, and M. Teboulle, “Proximal alternating linearized minimization for nonconvex and nonsmooth problems,” _Mathematical Programming_ , vol. 146, no. 1, pp. 459–494, 2014.
  * [33] R. T. Rockafellar and R. J.-B. Wets, _Variational analysis_. Springer Science & Business Media, 2009, vol. 317\.
  * [34] ——, _Variational Analysis_ , 2nd ed., ser. Grundlehren der mathematischen Wissenschaften. Berlin Heidelberg: Springer–Verlag, 2004, vol. 317.
  * [35] H. Attouch, J. Bolte, and B. F. Svaiter, “Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward–backward splitting, and regularized gauss–seidel methods,” _Mathematical Programming_ , vol. 137, no. 1, pp. 91–129, 2013.
  * [36] P. Tseng, “Convergence of a block coordinate descent method for nondifferentiable minimization,” _Journal of optimization theory and applications_ , vol. 109, no. 3, pp. 475–494, 2001.
  * [37] Y. Nesterov, “A method for unconstrained convex minimization problem with the rate of convergence o (1/k^ 2),” in _Doklady an ussr_ , vol. 269, 1983, pp. 543–547.
  * [38] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” _SIAM journal on imaging sciences_ , vol. 2, no. 1, pp. 183–202, 2009.
  * [39] P. H. Schönemann, “A generalized solution of the orthogonal procrustes problem,” _Psychometrika_ , vol. 31, no. 1, pp. 1–10, 1966.
  * [40] H. Le, N. Gillis, and P. Patrinos, “Inertial block proximal methods for non-convex non-smooth optimization,” in _International Conference on Machine Learning_. PMLR, 2020, pp. 5671–5681.
  * [41] T. Pock and S. Sabach, “Inertial proximal alternating linearized minimization (ipalm) for nonconvex and nonsmooth problems,” _SIAM Journal on Imaging Sciences_ , vol. 9, no. 4, pp. 1756–1787, 2016.
  * [42] X. Gao, X. Cai, and D. Han, “A gauss–seidel type inertial proximal alternating linearized minimization for a class of nonconvex optimization problems,” _Journal of Global Optimization_ , vol. 76, no. 4, pp. 863–887, 2020.
  * [43] Y. Xu and W. Yin, “A globally convergent algorithm for nonconvex optimization based on block coordinate update,” _Journal of Scientific Computing_ , vol. 72, no. 2, pp. 700–734, 2017.
  * [44] H. Liu, A. M.-C. So, and W. Wu, “Quadratic optimization with orthogonality constraint: explicit Łojasiewicz exponent and linear convergence of retraction-based line-search and stochastic variance-reduced gradient methods,” _Mathematical Programming_ , pp. 1–48, 2019.
  * [45] L. Zeng, P. Yu, and T. K. Pong, “Analysis and algorithms for some compressed sensing models based on l1/l2 minimization,” _SIAM Journal on Optimization_ , vol. 31, no. 2, pp. 1576–1603, 2021.
  * [46] Z. Zhou and A. M.-C. So, “A unified approach to error bounds for structured convex optimization problems,” _Mathematical Programming_ , vol. 165, no. 2, pp. 689–728, 2017.
  * [47] H. Liu, M.-C. Yue, and A. M.-C. So, “On the estimation performance and convergence rate of the generalized power method for phase synchronization,” _SIAM Journal on Optimization_ , vol. 27, no. 4, pp. 2426–2446, 2017.
  * [48] H. Attouch, J. Bolte, P. Redont, and A. Soubeyran, “Proximal alternating minimization and projection methods for nonconvex problems: An approach based on the Kurdyka-Łojasiewicz inequality,” _Mathematics of Operations Research_ , vol. 35, no. 2, pp. 438–457, 2010.
  * [49] G. Li and T. K. Pong, “Calculus of the exponent of kurdyka–łojasiewicz inequality and its applications to linear convergence of first-order methods,” _Foundations of computational mathematics_ , vol. 18, no. 5, pp. 1199–1232, 2018.
  * [50] C. Kim and D. Klabjan, “A simple and fast algorithm for L1-norm kernel PCA,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 42, no. 8, pp. 1842–1855, 2020.
  * [51] C.-C. Chang and C.-J. Lin, “Libsvm: a library for support vector machines,” _ACM transactions on intelligent systems and technology (TIST)_ , vol. 2, no. 3, pp. 1–27, 2011.
  * [52] P. P. Markopoulos, G. N. Karystinos, and D. A. Pados, “Optimal algorithms for $L_{1}$-subspace signal processing,” _IEEE Transactions on Signal Processing_ , vol. 62, no. 19, pp. 5046–5058, 2014.

###  -A Proof of Lemma 1

###### Proof.

We first prove (i). Due to $\bm{C}^{k}\in{\mathcal{B}}(n,d)\times{\rm
St}(d,K)\times{\rm St}(d,K)$ for all $k\geq 0$ and the boundness of
${\mathcal{B}}(n,d)$ and ${\rm St}(d,K)$, the sequence
$\\{\bm{C}^{k}\\}_{k\geq 0}$ is automatically bounded.

We next prove (ii). To simplify the notation, let

$\begin{split}&\ \bm{\Delta}_{\bm{P}}^{k+1}=\bm{P}^{k+1}-\bm{P}^{k},\
\bm{\Delta}_{\bm{Q}}^{k+1}=\bm{Q}^{k+1}-\bm{Q}^{k},\\\ &\
L_{k}=\|\bm{X}\bm{P}^{k}+\bm{P}^{k^{T}}\bm{X}^{T}\|,\ \forall\ k\geq
0.\end{split}$

According to the updates (13), we have

$\langle\bm{X}(\bm{P}^{k+1}-\bm{P}^{k}),\bm{E}^{k}\rangle\geq\frac{\alpha_{k}}{2}\|\bm{P}^{k+1}-\bm{P}^{k}\|^{2}_{F}.$
(24)

This, together with
$H(\bm{P},\bm{Q})=-\langle\bm{P},\bm{X}^{T}\bm{Q}\bm{Q}^{T}\rangle$ and
$\bm{E}^{k}=\bm{Q}^{k}\bm{Q}^{k^{T}}+\gamma_{k}(\bm{Q}^{k}\bm{Q}^{k^{T}}-\bm{Q}^{k-1}\bm{Q}^{{k-1}^{T}})$,
implies

$\begin{split}&\
H(\bm{P}^{k+1},\bm{Q}^{k})-H(\bm{P}^{k},\bm{Q}^{k})\leq-\frac{\alpha_{k}}{2}\|\bm{\Delta}_{\bm{P}}^{k+1}\|_{F}^{2}\\\
&\
+\gamma_{k}\langle\bm{X}(\bm{P}^{k+1}-\bm{P}^{k}),\bm{Q}^{k}\bm{Q}^{k^{T}}-\bm{Q}^{k-1}\bm{Q}^{{k-1}^{T}}\rangle\\\
&\
\leq-\frac{\alpha_{k}}{2}\left(1-\gamma_{k}\right)\|\bm{\Delta}_{\bm{P}}^{k+1}\|_{F}^{2}+\frac{2\gamma_{k}\|\bm{X}\|^{2}}{\alpha_{k}}\|\bm{\Delta}_{\bm{Q}}^{k}\|^{2}_{F},\end{split}$

where the second inequality is due to the following result:

$\begin{split}&\
\langle\bm{X}(\bm{P}^{k+1}-\bm{P}^{k}),\bm{Q}^{k}\bm{Q}^{k^{T}}-\bm{Q}^{k-1}\bm{Q}^{{k-1}^{T}}\rangle\\\
=&\
\langle\bm{P}^{k+1}-\bm{P}^{k},\bm{X}^{T}\bm{Q}^{k-1}(\bm{Q}^{k^{T}}-\bm{Q}^{{k-1}^{T}})\\\
+&\
\bm{X}^{T}(\bm{Q}^{k}-\bm{Q}^{k-1})\bm{Q}^{k^{T}}\rangle\leq\frac{\alpha_{k}}{4}\|\bm{\Delta}_{\bm{P}}^{k+1}\|_{F}^{2}\\\
+&\
\frac{\|\bm{X}\bm{Q}^{k-1}\|^{2}}{\alpha_{k}}\|\bm{\Delta}_{\bm{Q}}^{k}\|^{2}_{F}+\frac{\alpha_{k}}{4}\|\bm{\Delta}_{\bm{P}}^{k+1}\|_{F}^{2}+\frac{\|\bm{X}\bm{Q}^{k}\|^{2}}{\alpha_{k}}\|\bm{\Delta}_{\bm{Q}}^{k}\|^{2}_{F}\\\
\leq&\
\frac{\alpha_{k}}{2}\|\bm{\Delta}_{\bm{P}}^{k+1}\|_{F}^{2}+\frac{2\|\bm{X}\|^{2}}{\alpha_{k}}\|\bm{\Delta}_{\bm{Q}}^{k}\|^{2}_{F},\end{split}$

where the first inequality uses the fact that
$2\langle\bm{A},\bm{B}\rangle\leq\rho\|\bm{A}\|_{F}^{2}+\|\bm{B}\|_{F}^{2}/\rho$
for any $\rho>0$ and $\|\bm{A}\bm{B}\|_{F}\leq\|\bm{A}\|_{F}\|\bm{B}\|_{F}$.
The last inequality is due to $\|\bm{X}\bm{Q}\|_{F}\leq\|\bm{X}\|_{F}$ for any
$\bm{Q}\in{\rm St}(d,K)$. Moreover, according to the update (12b) and [32,
Lemma 2], we have

$H(\bm{P}^{k+1},\bm{Q}^{k+1})-H(\bm{P}^{k+1},\bm{Q}^{k})\leq-\frac{1}{2}\left(\beta_{k}-L_{k}\right)\|\bm{\Delta}_{\bm{Q}}^{k+1}\|_{F}^{2}.$
(25)

Then, we have

$\begin{split}&\
\Phi_{\beta_{*}}(\bm{C}^{k+1})-\Phi_{\beta_{*}}(\bm{C}^{k})\\\ =&\
H(\bm{P}^{k+1},\bm{Q}^{k+1})-H(\bm{P}^{k},\bm{Q}^{k})+\frac{\beta_{*}}{2}\|\bm{\Delta}_{\bm{Q}}^{k+1}\|^{2}_{F}-\frac{\beta_{*}}{2}\|\bm{\Delta}_{\bm{Q}}^{k}\|^{2}_{F}\\\
\leq&\
-\frac{\alpha_{k}}{2}\left(1-\gamma_{k}\right)\|\bm{\Delta}_{\bm{P}}^{k+1}\|_{F}^{2}-\frac{1}{2}\left(\beta_{k}-L_{k}-\beta_{*}\right)\|\bm{\Delta}_{\bm{Q}}^{k+1}\|_{F}^{2}\\\
&\
-\frac{1}{2}\left(\beta_{*}-\frac{4\gamma_{k}\|\bm{X}\|^{2}}{\alpha_{k}}\right)\|\bm{\Delta}_{\bm{Q}}^{k}\|^{2}_{F}\leq-\kappa_{1}\|\bm{C}^{k+1}-\bm{C}^{k}\|_{F}^{2},\end{split}$

where the first inequality follow from (-A) and (25) and the second inequality
is due to (19) and
$\kappa_{1}=\min\\{\alpha_{*}(1-\gamma^{*})/2,\beta^{*}/4\\}$.

Finally, we prove (iii). According to the updates (13), (12b), we have

$\mathbf{0}\in-\bm{X}^{T}\bm{E}^{k}+\alpha_{k}(\bm{P}^{k+1}-\bm{P}^{k})+\mathcal{N}_{\mathcal{B}(n,d)}(\bm{P}^{k+1})$
(26) $\displaystyle\mathbf{0}\in$ $\displaystyle\
-(\bm{X}\bm{P}^{k+1}+\bm{P}^{{k+1}^{T}}\bm{X}^{T})\bm{Q}^{k}\bm{Q}^{k^{T}}$
(27) $\displaystyle\ +\beta_{k}(\bm{Q}^{k+1}-\bm{Q}^{k})+\mathcal{N}_{{\rm
St}(d,K)}(\bm{Q}^{k+1}).$

This, together with [48, Proposition 2.1], yields

$\displaystyle\
\partial\Phi_{\beta_{*}}\left(\bm{P}^{k+1},\bm{Q}^{k+1},\bm{Q}^{k}\right)$
(28) $\displaystyle=$ $\displaystyle\
\left\\{-\bm{X}^{T}\bm{Q}^{k+1}\bm{Q}^{{k+1}^{T}}+\mathcal{N}_{\mathcal{B}(n,d)}(\bm{P}^{k+1})\right\\}$
$\displaystyle\times$ $\displaystyle\ \left\\{\begin{split}&\
-(\bm{X}\bm{P}^{k+1}+\bm{P}^{{k+1}^{T}}\bm{X}^{T})\bm{Q}^{k+1}\bm{Q}^{{k+1}^{T}}\\\
&\ +\beta_{*}(\bm{Q}^{k+1}-\bm{Q}^{k})+\mathcal{N}_{{\rm
St}(d,K)}(\bm{Q}^{k+1})\end{split}\right\\}$ $\displaystyle\times$
$\displaystyle\left\\{\beta_{*}(\bm{Q}^{k}-\bm{Q}^{k+1})\right\\}.$

Combining this with (26) and (27) yields

$\begin{split}&\
\operatorname{dist}^{2}\left(\mathbf{0},\partial\Phi_{\beta_{*}}\left(\bm{P}^{k+1},\bm{Q}^{k+1},\bm{Q}^{k}\right)\right)\\\
\leq&\
\|\bm{X}^{T}(\bm{E}^{k}-\bm{Q}^{k+1}\bm{Q}^{{k+1}^{T}})+\alpha_{k}\bm{\Delta}_{\bm{P}}^{k+1}\|_{F}^{2}\\\
+&\
\|(\bm{X}\bm{P}^{k+1}+\bm{P}^{{k+1}^{T}}\bm{X}^{T})(\bm{Q}^{k}\bm{Q}^{{k}^{T}}-\bm{Q}^{k+1}\bm{Q}^{{k+1}^{T}})+\\\
&\
(\beta_{*}-\beta_{k})\bm{\Delta}_{\bm{Q}}^{k+1}\|_{F}^{2}+\|\beta_{*}\bm{\Delta}_{\bm{Q}}^{k+1}\|_{F}^{2}\\\
\leq&\
4\|\bm{X}\|^{2}(\|\bm{Q}^{k+1}\bm{Q}^{{k+1}^{T}}-\bm{Q}^{k}\bm{Q}^{{k}^{T}}\|_{F}^{2}\\\
+&\
\gamma_{k}^{2}\|\bm{Q}^{k}\bm{Q}^{{k}^{T}}-\bm{Q}^{k-1}\bm{Q}^{{k-1}^{T}}\|_{F}^{2})+2\alpha_{k}^{2}\|\bm{\Delta}_{\bm{P}}^{k+1}\|_{F}^{2}\\\
+&\
2\|\bm{X}\bm{P}^{k+1}+\bm{P}^{{k+1}^{T}}\bm{X}^{T}\|^{2}\|\bm{Q}^{k}\bm{Q}^{{k}^{T}}-\bm{Q}^{k-1}\bm{Q}^{{k-1}^{T}}\|_{F}^{2}\\\
+&\
2(\beta_{k}-\beta_{*})^{2}\|\bm{\Delta}_{\bm{Q}}^{k+1}\|_{F}^{2}+\beta_{*}^{2}\|\bm{\Delta}_{\bm{Q}}^{k+1}\|_{F}^{2}\\\
\leq&\
16\|\bm{X}\|^{2}(\|\bm{\Delta}^{k+1}_{\bm{Q}}\|_{F}^{2}+\gamma^{2}_{k}\|\bm{\Delta}^{k}_{\bm{Q}}\|_{F}^{2})+2\alpha_{k}^{2}\|\bm{\Delta}_{\bm{P}}^{k+1}\|_{F}^{2}\\\
+&\
8nd\|\bm{X}\|^{2}\|\bm{\Delta}^{k}_{\bm{Q}}\|_{F}^{2}+(2(\beta_{k}-\beta_{*})^{2}+\beta_{*}^{2})\|\bm{\Delta}_{\bm{Q}}^{k+1}\|_{F}^{2}\\\
\leq&\
2\alpha_{k}^{2}\|\bm{\Delta}_{\bm{P}}^{k+1}\|_{F}^{2}+8(2\gamma^{2}_{k}+nd)\|\bm{X}\|^{2}\|\bm{\Delta}^{k}_{\bm{Q}}\|_{F}^{2}\\\
+&\
(16\|\bm{X}\|^{2}+2(\beta_{k}-\beta_{*})^{2}+\beta_{*}^{2})\|\bm{\Delta}_{\bm{Q}}^{k+1}\|_{F}^{2}.\end{split}$

By taking
$\kappa_{2}=\max\\{\sqrt{2}\alpha^{*},2\sqrt{4\gamma^{*^{2}}+2nd}\|\bm{X}\|,\\\
\sqrt{16\|\bm{X}\|^{2}+2(\beta_{k}-\beta_{*})^{2}+\beta_{*}^{2}}\\}$, we
obtain the desired result (22).  $\sqcup\hbox to0.0pt{\hss$\sqcap$}$

###  -B Proof of Corollary 1

###### Proof.

Using Fact 1 with setting $\bm{B}=\bm{I}$ and
$\bm{A}=(\bm{X}\bm{P}_{i}+\bm{P}_{i}^{T}\bm{X}^{T})/{2}$, we have that the KŁ
exponent of $\ell_{i}\left(Q\right)$ is $1/2$ for all $i=1,\dots,2^{nd}$.
Using the same argument in the proof in [14, Section 3.3] yields that the KŁ
exponent of $h$ is $1/2$. This, together with (18) and [49, Theorem 3.6],
implies the KŁ exponent of $\Phi_{\beta}$ is $1/2$.  $\sqcup\hbox
to0.0pt{\hss$\sqcap$}$

###  -C Proof of Lemma 2

###### Proof.

We first prove (i). According to [48, Proposition 2.1], we have

$\begin{split}\partial h(\bm{P},\bm{Q})=&\
\left\\{-\bm{X}^{T}\bm{Q}\bm{Q}^{T}+\mathcal{N}_{\mathcal{B}(n,d)}(\bm{P})\right\\}\\\
\times&\ \\{-(\bm{X}\bm{P}+\bm{P}^{T}\bm{X}^{T})\bm{Q}+\mathcal{N}_{{\rm
St}(d,K)}(\bm{Q})\\}.\end{split}$ (29)

This, together with $(\mathbf{0},\mathbf{0})\in\partial h(\bm{P},\bm{Q})$,
$\bm{P}\in\operatorname{sgn}(\bm{X}^{T}\bm{Q}\bm{Q}^{T})$, and (2), implies
the desired result. Next, we prove (ii). According to
$\mathbf{0}\in\partial\Phi_{\beta}(\bm{P},\bm{Q},\bm{Q}^{\prime})$ and (28),
we obtain $\bm{Q}^{\prime}=\bm{Q}$. This, together with (29), implies that
$\mathbf{0}\in\partial\Phi_{\beta}(\bm{P},\bm{Q},\bm{Q})$ if and only if
$\mathbf{0}\in\partial h(\bm{P},\bm{Q})$.  $\sqcup\hbox to0.0pt{\hss$\sqcap$}$


# Explicitly Encouraging Low Fractional Dimensional Trajectories Via
Reinforcement Learning

Sean Gillen
Department of Electrical Engineering
University of California Santa Barbara United States
<EMAIL_ADDRESS>
&Katie Byl
Department of Electrical Engineering
University of California Santa Barbara United States
<EMAIL_ADDRESS>

###### Abstract

A key limitation in using various modern methods of machine learning in
developing feedback control policies is the lack of appropriate methodologies
to analyze their long-term dynamics, in terms of making any sort of guarantees
(even statistically) about robustness. The central reasons for this are
largely due to the so-called curse of dimensionality, combined with the black-
box nature of the resulting control policies themselves. This paper aims at
the first of these issues. Although the full state space of a system may be
quite large in dimensionality, it is a common feature of most model-based
control methods that the resulting closed-loop systems demonstrate dominant
dynamics that are rapidly driven to some lower-dimensional sub-space within.
In this work we argue that the dimensionality of this subspace is captured by
tools from fractal geometry, namely various notions of a fractional dimension.
We then show that the dimensionality of trajectories induced by model free
reinforcement learning agents can be influenced adding a post processing
function to the agents reward signal. We verify that the dimensionality
reduction is robust to noise being added to the system and show that that the
modified agents are more actually more robust to noise and push disturbances
in general for the systems we examined.

> Keywords: Locomotion, Reinforcement Learning

##  1 Introduction

The availability of computation as a resource has been growing exponentially
since at least the 1970s, and there is every indication that this resource
will continue to become cheaper and more available well into the conceivable
future. Researchers have been able to leverage the large amounts compute
available to better control robotic systems, and advances in computational
capacity and algorithmic development continue to open up new domains. One
promising manifestation of this is model-free reinforcement learning, a branch
of machine learning which allows an agent to interact with its environment and
autonomously learn how to maximize some measure of reward. The promise here is
to allow researchers to solve problems for systems that are hard to model,
and/or that the user doesn’t know how to solve themselves. Recent examples in
the context of robotics include controlling a 47 DOF humanoid to navigate a
variety of obstacles [heess_emergence_2017], dexterously manipulating objects
with a 24 DOF robotic hand [openai_learning_2018], and allowing a physical
quadruped robot to run [hwangbo_learning_2019], and recover from falls
[lee_robust_2019].

In this paper we study legged locomotion. This class of problems is
notoriously difficult, and as a result reinforcement learning is a popular
tool to throw at it. We would argue that hand-designed, model based control
still represents the state of the art (a la Boston Dynamics), but RL has been
a fruitful approach. There are examples of learned policies outperforming hand
designed ones [hwangbo_learning_2019], and there is good reason to believe
these learning methods will continue their current trajectory of increasing
performance gains and ease of use. But these algorithms have a serious draw
back in that they are mostly black boxes. It is an open challenge to figure
out what exactly it is that your RL agent has learned. If all you know is that
one of your agents achieved very high reward, it is not clear how to verify
that this system is safe and sensible in all the regions of state space it
will visit during its life. Nor can we necessarily say anything about the
stability or robustness properties of the system. Recent work [Taleledeep] has
used so-called mesh-based tools to examine precisely these questions.

However, utility of any mesh-based tool to accurately discretize a state-space
is limited, due to the curse of dimensionality. In practice, these methods are
only able to work on relatively high dimensional systems if the reachable
space grows at a rate that is much smaller than the exponential growth of the
full state space the system within which it is embedded. To expand these
methods to higher dimensional systems. We will need to find ways to keep the
volume of visited states from expanding commensurately. One way to quantify
this rate of growth is by using one of the several notions of ”fractional
dimensions” from fractal geometry.

In this work, we discuss an efficient meshing algorithm, which we call box
meshing. We show that this approach makes calculating the so called mesh
dimension feasible in the context of reinforcement learning. We also propose
using other notions of fractional dimension from the literature as a proxy for
the property we care about. We then show that reinforcement learning agents
can be trained to shrink these measures by post processing their reward
function. We present the results of this training, and finally present some
brief analysis of the resulting structure for select policies.

##  2 Meshing & Fractional Dimensions

(a) Scaling in different dimensions

(b) A non uniform mesh

(c) Calculating the mesh dimension

Figure 1: Image credit for sub-figure a: [BrendanRyan/Publicdomain2020], image
credit for sub-figures b and c: [Talelepush]

Let’s say we have a continuous set S that we want to approximate by selecting
a discrete set M composed of regions in S. We will call this set M a mesh of
our space. Figure 1(a) shows some examples of this: a line is broken into
segments, a square into grid spaces, and so on. The question is: as we
increase the resolution of these regions, how many more regions N do we need?
Again, figure 1(a) shows us some very simple examples. For a D dimensional
system, if we go from regions of size d to d/k, then we would expect the
number of mesh points to scale as $N\propto k^{D}$. But not all systems will
scale like this, as 1(b) and 1(c) illustrate. Figure 1(b) is an example of a
curve embedded in a two dimensional space. The question of how many mesh
points are required must be answered empirically. Going backwards, we can use
this relationship to assign a notion of ”dimension” to the curve.

$D_{f}=-\lim_{k\rightarrow 0}\frac{\log N(k)}{\log k}$ (1)

What we are talking about is called the Minkowski–Bouligand dimension, also
known as the box counting dimension. This dimension need not be an integer,
hence the name ”fractional dimension”. As a practical matter, we use the slope
of the log-log plot of mesh sizes over d to calculate this, rather than taking
a limit. This is one of many measures of ”fractional dimension” that that
emerged from the study of fractal geometry. Although these measures were
invented to study fractals, they can still be usefully applied to non-fractal
sets.

In [OguzSaglam2015], Saglam and Byl introduced a technique that is able to
simultaneously build a non-uniform mesh of a reachable state space while
developing robust policies for a bipedal walker on rough terrain. Having a
discrete mesh allows for value iteration over several candidate controllers,
which found a robust control policy. In addition this mesh allows for the
construction of a state transition matrix, which was used to calculate the
mean first passage time [Byl2009], a metric that quantifies the expected
number of steps a meta-stable system can take before falling.

Since its introduction, meshing in this fashion has been used for designing
walking controllers robust to push disturbances [Talelepush], to design agile
gaits for a quadruped [Byl2017], and to analyze hybrid zero dynamics (HZD)
controllers [Saglam2015]. There has also been recent work to use these tools
to analyze policies trained by deep reinforcement learning [Taleledeep]. A
long term goal and motivation for this work is to take a high performance
controller obtained via reinforcement learning, and extract from it a mesh-
based policy that is both explainable and amenable to analysis.

###  2.1 Box Meshing

Our primary improvement to the prior work on meshing is to introduce something
we call box meshing. Prior, a new mesh point could take any value in the state
space. To determine if a new state is already in the mesh, we would compute a
distance metric to every point in the mesh, and check if the minimum was below
our threshold. Thus, building the mesh was an $O(n)^{2}$ algorithm. By
contrast, in box meshing we apriori divide the space uniformly into boxes with
side length $d$. We identify any state s with a key obtained by:
$\text{key}=\text{round}(\frac{s}{d})d$, where round performs an element-wise
rounding to the nearest integer. We can then use these keys to store mesh
points in a hash table. Using this data structure, we can still store the mesh
compactly, only keeping the points we come across. However, insertion and
search are now $O(1)$, and so building the mesh is $O(n)$. This is very
similar to non-hierarchical bucket methods, which are well studied spatial
data structure [Samet1990], although we are using them for data compression
here. In the prior meshing work, this sort of speedup would be minor, the run-
time is dominated by the simulator or robot. However, this speedup does open
some new possibilities: most poignantly, it makes calculating the mesh
dimension during reinforcement learning plausible.

###  2.2 Algorithmic Box Mesh Dimension

The ”mesh dimension” is the quantity extracted from the slope of the log log
plot of mesh sizes vs d values. For this paper, it is assumed that the mesh
algorithm being used for this calculation is the box mesh. Automatically
computing the mesh dimension of a data set generated from learning agent with
speed, accuracy, and robustness is very challenging. A single trajectory
provides only a small amount of data, which adds significant noise to the mesh
sizes. Agents might do things like fall over and generate extremely short
trajectories, or learn a trajectory that ”stands in place”, which can lead to
numerical errors. Finally, every decision is a trade-off between accuracy and
speed. Model-free RL is predicated on having a huge number of rollouts to
learn from, and we would like for any mesh-dimension quantification algorithm
to be fast enough so as to not dominate the total learning time. With these
factors in mind, we introduce two box mesh dimensions. The lower mesh
dimension does the linear fit, but intentionally errs on the side of including
flat parts of the graph, and therefore tends to underestimate the true mesh
dimension. We then have the upper mesh dimension, which takes the largest
slope in the log log relationship, thus tending to overestimate the true mesh
dimension. Neither of these measures are correct, but taken together they can
bound the mesh dimension, and as we will see they can be useful on their own.

Input: State set $S$, box size d.

Output: Mesh size m.

Initialize: Empty hash table M.

for  _s $\in$ S_ do

$\bar{s}=\text{Normalize(s)}$

key = round($\bar{s}$ / d)d

if  _s $\in$ M_ then

M[key]++

else

M[key]=1

end if

end for

Return: M

Algorithm 1 Create Box Mesh, see section 2.1

Input: State set $S$.

Output: Mesh M.

Hyperparameters: scaling factor f, initial box size $d_{0}$.

Initialize: Empty list of mesh sizes H, empty list of d values D.

m = Size(CreateBoxMesh(S, $d_{0}$))

d = $d_{0}$

Append m to H, append d to D.

while  _m $<$ size(S)_ do

d = d/f

m = Size(CreateBoxMesh(S, d))

Prepend m to H, prepend d to D.

end while

while  _m $\neq$ 1_ do

d = d*f

m = Size(CreateBoxMesh(S,d))

Append m to H, append d to D.

end while

X = $\log$ d

Y = -$\log$ m

Lower Mesh Dim: fit Y = gX + b, Return: g

Upper Mesh Dim: w = greatest slope in Y over X Return: w

Algorithm 2 Compute Box Mesh Dimension, see section 2.2

##  3 Reinforcement Learning

The goal of reinforcement learning is to train an agent acting in an
environment to maximize some reward function. At every timestep
$t\in\mathbb{Z}$, the agent receives the current state $s_{t}\in R^{n}$, uses
that to compute an action $a_{t}\in\mathbb{R}^{b}$, and the receives the next
state $s_{t+1}$, which is used to calculate a reward
$r:\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{n}\rightarrow\mathbb{R}$.
The objective is to find a policy
$\pi_{\theta}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$

$\operatorname*{arg\,max}_{\theta}\mathop{\mathbb{E}}_{\eta}\left[\sum_{t=0}^{T}r(s_{t},a_{t},s_{t+1})\right]$
(2)

Where $\theta\in\mathbb{R}^{d}$ is a set that parameterizes the policy, and
$\eta$ is a parameter representing the randomness in the environment. This
includes the random initial conditions for episodes.

###  3.1 Post Processing Rewards

In order to influence the dimensionality of the resulting policies, we
introduce various postprocessors, which act on the reward signals before
passing them to the agent. These obviously modify the problem: in some sense
the postprocessed environment is a completely different problem from the
original. However our meta-goal is to train agents that achieve reasonable
rewards in the base environment, while simultaneously exhibiting reduced
dimensionality we are looking for. These postprocessors take the form:

$R_{*}(\mathbf{s},\mathbf{a})=\frac{1}{D_{*}(\mathbf{s})}\sum_{t=0}^{T}r(s_{t},a_{t},s_{t+1})$
(3)

Where $\mathbf{s},\mathbf{a}$ are understood to be an entire trajectory of
state action pairs, and $D_{*}$ is some measure of fractional dimension. Some
measures of dimensionality can be inserted here directly (See 7.1). However
the mesh dimensions computed by algorithm 2 require a little more care. We
must first define a clipped dimension:

$D^{c}_{*}=\text{clip}[D_{*}(\mathbf{s}_{t>Tr}),1,D_{t}/2)]$ (4)

where $D_{t}$ is the topological dimension, equal to the number of states in
the system. $Tr$ is a fixed timestep chosen to exclude the initial transients
resulting from a system moving from rest to into a quasi-cyclical “gait”. In
this paper we set $Tr$ = 200 for all experiments. For comparison, the nominal
episode length is 1000. The clipping is to ensure that the pathological
trajectories that and RL agent sometimes generates don’t interfere with the
training. It will also clip trajectories that terminate early, to prevent
agents learning to fall over immediately to “game the system”. Half of the
topological dimension proved to be a decent upper bound for the worst case
dimensionality of each system in practive. The mesh dimension postprocessors
use the clipped dimension. Finally, when $D_{*}=1$ is used, we call the result
is the identity post processor, since in this case the total reward is
completely unchanged.

###  3.2 Environments

We examine a subset of the popular OpenAI Mujoco locomotion environments
introduced in [1606.01540]. In particular, we evaluate our work on
HalfCheetah-v2, Hopper-v2, and Walker2d-v2. These environments were chosen
because they have a relatively high dimensionality (11-17 DOF), yet we believe
can be made feasible for meshing based approaches. The state space consists of
all joint / base positions and velocities, with the x (the ”forward”) position
being held out, because we want a policy that is invariant along that
dimension.

###  3.3 Augmented Random Search

In [Mania2018] Mania et al introduce Augmented Random Search (ARS) which
proved to be efficient and effective on the locomotion tasks. Rather than a
neural network, ARS used static linear policies, and compared to most modern
reinforcement learning, the algorithm is very straightforward. The algorithm
operates directly on the policy weights, each epoch the agent perturbs it’s
current policy N times, and collects 2N rollouts using the modified policies.
The rewards from these rollouts are used to update the current policy weights,
repeat until completion. The algorithm is known to have high variance; not all
seeds obtain high rewards, but to our knowledge their work in many ways
represents the state of the art on these benchmarks. Mania et al introduce
several small modifications of the algorithm in their paper, our
implementation corresponds to the version they call ARS-V2t.

###  3.4 Training

To keep things simple, we wanted to find one set of hyper parameters for all
environments and post processors. These parameters were chosen by hand, with
the parameters reported in Table 9 from [Mania2018] as a starting point. We
tuned until our unprocessed learning achieved satisfactory results across all
tasks. Again, ARS is known to have high variance between random seeds, and
some seeds never learn to gather a large reward. The parameters we found are
able to consistently solve the cheetah and walker; for the hopper, the
algorithm learns a policy with high reward around half the time. This seems
consistent with the performance reported in [Mania2018]. We train each
postprocessor on 10 random seeds, the evaluation metrics are averages over 5
rollouts from each seed, and for the dimension metrics we use extended
episodes of length 10,000 to get a more accurate measurement. The reported
returns, and the training, both use the normal 1,000 step episodes. We found
that the mesh postprocessors were getting very poor performance when trained
from a random policy. However, we found that we saw good results when these
trials were initialized with a working policy. Therefore we trained agents for
750 epochs without post processing, and used that to initialize the mesh
dimension policies. The mesh policies were then trained for an additional 250
epochs and the results reported.

##  4 Results

###  4.1 Mesh Dimension Postprocessors

(a) HalfCheetah

(b) Hopper

(c) Walker

Figure 2: Reward curves for mesh dimension postproccesor runs.

Environment | Postprocessor | Lower Mesh Dim. | Upper Mesh Dim. | Return
---|---|---|---|---
HalfCheetah-v2 | Identity | 2.31 $\pm$ 0.71 | 7.34 $\pm$ 1.56 | 5469 $\pm$ 823
Lower Mesh Dim | 0.66 $\pm$ 0.51 | 2.55 $\pm$ 1.52 | 4962 $\pm$ 598
Upper Mesh Dim. | 1.06 $\pm$ 1.13 | 2.83 $\pm$ 1.27 | 4432 $\pm$ 539
Hopper-v2 | Madogram∗ | 1.62 $\pm$ .27 | 4.68 $\pm$ 0.82 | 3461 $\pm$ 119
Lower Mesh Dim. | 1.13 $\pm$ .02 | 3.54 $\pm$ 0.96 | 2941 $\pm$ 538
Upper Mesh Dim. | 1.27 $\pm$ .50 | 2.98 $\pm$ 1.48 | 3020 $\pm$ 337
Walker2d-v2  (walking seeds)∗∗ | Identity | 2.13 $\pm$ 0.31 | 4.62 $\pm$ 1.03 | 3758 $\pm$ 1037
Lower Mesh Dim. | 1.21 $\pm$ 0.06 | 4.09 $\pm$ 1.03 | 3339 $\pm$ 887
Upper Mesh Dim. | 1.89 $\pm$ 0.42 | 3.10 $\pm$ 0.93 | 3359 $\pm$ 903
Walker2d-v2  (all seeds)∗∗ | Identity | 2.13 $\pm$ 0.31 | 4.62 $\pm$ 1.03 | 3758 $\pm$ 1037
Lower Mesh Dim. | 1.04 $\pm$ 0.53 | 4.45 $\pm$ 1.19 | 3034 $\pm$ 1086
Upper Mesh Dim. | 1.48 $\pm$ 0.67 | 2.27 $\pm$ 0.95 | 2556 $\pm$ 1378

Table 1:  Mesh dimensions and returns for trajectories after training. See 3.4
for details * Because ARS with our chosen hyper parameters does not
consistently produce 10 seeds that perform well on the hopper, we instead use
madodiv (see the 7.1) for the seed policies. * See 4.1

For all environments the mesh post processors had a significant impact on the
mesh dimensions. It’s important to remember here, the dimensions reported
represent lower and upper bounds for the actual mesh dimensions. There was
also a corresponding and significant decrease in the unprocessed rewards.
However with our meta goal of training agents that have acceptable reward but
which are more amenable to meshing, this is a more than acceptable trade. In
the case of walker, several seeds (4 for the upper dim., 3 for the lower mesh
dim.), ”forget how to walk”, and learn a policy that stands in place. This
certainly has a low dimensionality, but is not very useful, to be complete we
include statistics from the seeds that learned a gait, and for all 10 seeds,
including the standing policies.

##  5 Analysis

We now examine the learned behavior for one of the more notable policies. By
far the most dramatic effect from the tables above was the mesh dimension
postprocessors on the cheetah. Both measures of dimension shrunk by 2-4 times.
Figure 7 presents data for this case.

Figure 3: Top: mesh sizes vs log of the box size for the cheetah environment.
Lower Left: Every five frames overlaid for the an identity policy on the
cheetah. Lower Right: Every five frames of cheetah after the lower mesh
dimension training.

Toward more intuitively understanding this data, a few comments are worth
making, first. We have discussed the mesh dimension rather abstractly so far.
In visualizing what this really means, imagine two different gait cycles. In
one case, there is a general pattern to the motion, but it wanders in a noisy-
looking way, like a “signature” that does not quite match up, cycle after
cycle. As motions become closer to being exact limit cycles, there is a more
clear pattern of repetition, exactly analogous to re-tracing the same path,
again and again, within the state space. Such a more tightly-structured limit
cycle nature in turn results in a significantly lower-dimensional set of
states being visited.

We can see from the mesh size curves in Figure 7 that there is an overwhelming
difference in the mesh sizes between the lower mesh dimension post processor
and the other two. To put this in perspective, before the extra 250 epochs of
training, if given a box size of .01, the agent would need a unique mesh point
for every single state in the 10,000 state trajectory. After the additional
training, the agent can represent all 10,000 points with just 5 mesh points!
In this case it appears both agents learned a quasi periodic gait with period
5. In figure 3 we present an overlay of the agents rendered every 5 steps. The
results show us that the mesh agent has learned an extremely tight limit
cycle. It’s a bit of a strange limit cycle, being only 5 timesteps long, but
nonetheless we think this is interesting and surprising behavior.

The behavior displayed in figure 3 is clearly something that can only happen
in a noiseless simulation, so we also measured the mesh dimensions of our
policies when subjected to noise during rollouts using the noise values from
the robustness experiments (see table LABEL:tab:robust)to inform parameter
values. The difference in fractal dimension is less pronounced than the no
noise case but is still a clear improvement. Furthermore we expect that if we
were to add noise at training time, that the learning may be able to find ways
to lower the mesh dimension that is more robust to noise.

MeshDim: f = 1.5, $d_{0}$ = 1e-2

###  7.1 Variation Estimators

As discussed, computing the mesh dimension automatically is fraught with
peril, in many practical scenarios. But there are also many other, different
metrics one might consider, to give various approximations to the fractional
dimension we seek to estimate. Gneiting et al. [Gneiting2012] compare a number
of these estimators, and submit that the variation estimator [Emery2005]
offers a very good trade off between speed and robustness. To obtain this
estimator, first define the power variation of order p as:

$P_{p}(X,l)=\frac{1}{(2n-l)}\sum_{i=l}^{n}|X_{i}-X_{i-l}|^{p}$ (5)

Then, we define the variation estimator of order p as:

$Dv_{p}(X)=2-\frac{\log P_{p}(X,2)-\log P_{p}(X,1)}{p\log 2}$ (6)

The madogram estimator is the special case of (6) where p = 1, and the
variogram is where p = 2.

###  7.2 Variational Postprocessors

(a) HalfCheetah

(b) Hopper

(c) Walker

Figure 4: Reward curves for the variation postprocessors Environment | Postprocessor | Variogram | Madogram | Lower Mesh Dim. | Return
---|---|---|---|---|---
HalfCheetah-v2 | Identity | 1.71 $\pm$ .03 | 1.42 $\pm$ .05 | 2.36 $\pm$ .61 | 5545 $\pm$ 593
Variogram | 1.68 $\pm$ .01 | 1.36 $\pm$ .02 | 2.06 $\pm$ .60 | 5136 $\pm$ 851
Madogram | 1.65 $\pm$ .02 | 1.31 $\pm$ .04 | 2.09 $\pm$ .64 | 5234 $\pm$ 950
Hopper-v2 | Identity∗ | 1.61 $\pm$ .14 | 1.22 $\pm$ .28 | 1.03∗ $\pm$ .71 | 2063 $\pm$ 1052
Variogram | 1.51 $\pm$ .02 | 1.03 $\pm$ .04 | 1.58 $\pm$ .54 | 3299 $\pm$ 711
Madogram | 1.51 $\pm$ .002 | 1.02 $\pm$ .004 | 1.57 $\pm$ .36 | 3449 $\pm$ 146
Walker2d-v2 | Identity | 1.68 $\pm$ .35 | 1.36 $\pm$ .71 | 2.14 $\pm$ .29 | 3742 $\pm$ 1038
Variogram | 1.54 $\pm$ .07 | 1.07 $\pm$ .01 | 1.85 $\pm$ .54 | 3779 $\pm$ 894
Madogram | 1.53 $\pm$ .01 | 1.06 $\pm$ .02 | 1.99 $\pm$ .53 | 3414 $\pm$ 1025
Table 4:  Mesh dimensions and returns for trajectories after training. See 3.4
for details * This includes policies which learned to ”stand still”, which
lowers the average mesh dimension considerably see discussion

It seems that the variational postprocessors had a modest effect the
variational dimension, but that does not seem to correlate to a smaller mesh
dimension, despite what our preliminary tests had led us to believe. The
hopper and walker did have remarkable consistency in the variation dimensions
they found; possibly this could be used to lower the variance in ARS. The fact
that the variogram and madogram also got higher performance on the hopper task
could support this claim. However without running many more trials and hyper
parameter sweeps, that’s not a claim that can be substantiated. These
experiments show that 1) measures for fractional dimension can be influenced
without adversely effecting the reward, and 2) that it is possible for an
agent to shrink it’s variogram and madogram dimensions without a large impact
on its mesh dimension.

###  7.3 Mesh Dimension Examples

Figure 5 illustrates two examples of the curves used to compute the mesh
dimension. Recall that to compute the mesh dimension, we choose several values
for d, the box length, and for each d construct a mesh using that box size.
The x axis of these plots represents the log of the box length used, the y
axis represents the log of size of the mesh created. For each curve, we
display the lower bound and upper bound for the dimension as computed by
algorithm 2, as well as several hand fits of the data. We hope that figure 5a
makes clear what a close to ideal situation looks like, and provide intuition
as to why the upper and lower mesh dimension bound the quantity we are trying
to measure. Figure 5b serves to illustrate some of the problems with making an
algorithmic measure of the dimension. There is much less data to work with due
to performance constraints, which causes a large amount of noise on the
estimate of the mesh dimension. Indeed even fitting this data by hand becomes
a challenge and we provide two fits which can both be argued to be ”correct”.
Which of these two represents the quantity we care about depends on the exact
system being used and the purpose of the meshes we want to build with the
resulting policy.

(a) High resolution curve of a well behaved policy

(b) Run time resolution curve of a typical policy

Figure 5: Mesh curve and mesh dimension examples

###  7.4 Implementation Details

For performance reasons, the mesh dimension algorithm does not actually create
meshes until the mesh size equals the total data size, but rather until the
mesh size is 4/5 the total data size. Figure 8a shows a typical mesh curve,
and we can see the long tail of values with mesh sizes close to the maximum
value. Not much useful information is gained from this and it is wasting time,
so we stop early. We do not place the same limitation on the lower size of the
mesh, since typically the mesh size hits one much more rapidly, again figure
8a illustrates this. In addition implement a minimum size for d, set to 1e-9
in this work to avoid numerical errors.

The normalization done during box creation uses a running mean and standard
deviation of all states seen so far during training. These stats are saved and
used for evaluation as well, we found that the upper mesh dimension is very
sensitive to the normalization used, but that the other metrics where not.