obtained by adding just prefixes:

$\displaystyle|\mathcal{S}_{k}(S,\Gamma)|$
$\displaystyle\geq\left|\left\\{\psi(P\oplus\phi(S)),\forall
P\in\\{P^{\prime}\in\mathcal{S}(\Gamma\cup\\{\xi\\}),~{}|P^{\prime}|\leq
k\\}\right\\}\right|$ $\displaystyle=\left|\left\\{\psi(P),\forall
P\in\\{P^{\prime}\in\mathcal{S}(\Gamma\cup\\{\xi\\}),~{}|P|\leq
k\\}\right\\}\right|$
$\displaystyle=\left|\bigcup_{i=0}^{k}\\{P^{\prime}\in\mathcal{S}(\Gamma),~{}|P^{\prime}|=i\\}\right|$
[Disjoint sets]
$\displaystyle=\sum_{i=0}^{k}\left|\\{P^{\prime}\in\mathcal{S}(\Gamma),~{}|P^{\prime}|=i\\}\right|$
$\displaystyle=\sum_{i=0}^{k}|\Gamma|^{i}$ [Geometric series]
$\displaystyle=\left\\{\begin{matrix}\frac{1-|\Gamma|^{k+1}}{1-|\Gamma|}&\text{if
}|\Gamma|>1\\\ k+1&\text{if }|\Gamma|=1\\\ \end{matrix}\right.$

∎

##  Appendix D Alternative attack designs

In this section, we cover alternative algorithmic designs called PGA-Charmer
for solving Eq. 1 Specifically, we study relaxing the binary constraints in
Eq. BP in order to perform a Projected Gradient Ascent (PGA) procedure.

###  D.1 PGA-Charmer

Let $\bm{E}^{(i)}=\textbf{Token}(S^{(i)})\bm{T}\in\mathbb{R}^{l_{i}\times d}$
be the embeddings of tokens for any sentence
$S^{(i)}\in\mathcal{S}^{\prime}\subseteq\mathcal{S}(\Gamma)$ with
$i\in\left[\left|\mathcal{S}^{\prime}\right|\right]$. The zero-padded
embeddings for the sentences in the set $\mathcal{S}^{\prime}$ become:

$\hat{\bm{E}}^{(i)}=\bm{E}^{(i)}\oplus\bm{0}_{(\overline{l}-l_{i})\times
d}\in\mathbb{R}^{\overline{l}\times d},\quad\forall
i\in\left[\left|\mathcal{S}^{\prime}\right|\right]\,,$

where
$\overline{l}=\max\\{l_{i}:i\in\left[\left|\mathcal{S}^{\prime}\right|\right]\\}$
and $\oplus$ is the concatenation operator along the first dimension.

######  Remark S1 (Model output after zero padding).

Given a function $f$, the output before and after zero padding is unchanged,
i.e.,
$\bm{f}(\hat{\bm{E}}^{(i)})=\bm{f}(\bm{E}^{(i)})~{}~{}\forall\bm{E}^{(i)}\in\mathcal{S}^{\prime}$.

We can reformulate the problem in Eq. 1 as:

$\begin{matrix}\underset{\bm{u}\in\mathbb{R}^{|\mathcal{S}_{k}(S,\Gamma)|}}{\max}&\mathcal{L}\left(\bm{f}\left(\sum_{i=1}^{|\mathcal{S}_{k}(S,\Gamma)|}u_{i}\cdot\hat{\bm{E}}^{(i)}\right),y\right)\\\
\text{s.t.}&u_{i}\in\\{0,1\\}~{}~{}\forall
i\in[|\mathcal{S}_{k}(S,\Gamma)|],\quad\left|\left|\bm{u}\right|\right|_{1}=1\end{matrix}\,,$
(BP)

which is a constrained binary optimization problem. Note that given
$\bm{u}^{\text{BP}}$ a maximizer of Eq. BP with
$i^{\text{BP}}:=\operatorname*{arg\max}_{i\in[|\mathcal{S}_{k}(S,\Gamma)|]}u_{i}^{\text{BP}}$,
we have that the sentence $S^{(i^{\text{BP}})}\in\mathcal{S}_{k}(S,\Gamma)$ is
a maximizer of Eq. 1.

However, solving Eq. BP is as hard as solving Eq. 1 because of the exponential
size of $\mathcal{S}_{k}(S,\Gamma)$, see Corollary 3.9. Alternatively, we can
relax the binary constraints from the $\bm{u}$ vector and solve:

$\begin{matrix}\max_{\bm{u}\in\mathbb{R}^{|\mathcal{S}_{k}(S,\Gamma)|}}&\mathcal{L}\left(\bm{f}\left(\sum_{i=1}^{|\mathcal{S}_{k}|}u_{i}\cdot\hat{\bm{E}}^{(i)}\right),y\right)\\\
\text{s.t.}&u_{i}\in[0,1]~{}~{}\forall
i\in[|\mathcal{S}_{k}(S,\Gamma)|],\quad\left|\left|\bm{u}\right|\right|_{1}=1.\end{matrix}$
(SP)

In this case, given $\bm{u}^{\text{SP}}$ a maximizer of Eq. SP, we know:

$\leavevmode\resizebox{208.13574pt}{}{$\mathcal{L}\left(\bm{f}\left(\sum_{i=1}^{|\mathcal{S}_{k}(S,\Gamma)|}u^{\text{SP}}_{i}\cdot\hat{\bm{E}}^{(i)}\right),y\right)\geq\mathcal{L}\left(\bm{f}\left(\sum_{i=1}^{|\mathcal{S}_{k}(S,\Gamma)|}u^{\text{BP}}_{i}\cdot\hat{\bm{E}}^{(i)}\right),y\right)$}\,.$

Note that the embeddings
$\sum_{i=1}^{|\mathcal{S}_{k}(S,\Gamma)|}u^{\text{SP}}_{i}\cdot\hat{\bm{E}}^{(i)}$
have no correspondence to any sentence in $\mathcal{S}_{k}(S,\Gamma)$.
However, we can still take
$i^{\text{SP}}=\operatorname*{arg\max}_{i\in[|\mathcal{S}_{k}(S,\Gamma)|]}u_{i}^{\text{SP}}$
and hopefully $S^{(i^{\text{SP}})}\in\mathcal{S}_{k}(S,\Gamma)$ is an
adversarial example. To solve Eq. SP, we employ projection gradient ascent
with step-size $\eta$ as follows:

$\begin{split}&\bm{u}^{t+1}=\Pi_{\Delta}(\bm{u}^{t}+\eta\nabla{\mathcal{L}_{\bm{u}}(\bm{u}^{t})})\,,\end{split}$

where $\Pi_{\Delta}(\cdot)$ is the projection function. Let us denote by
$\hat{\bm{u}}:=\bm{u}^{t}+\eta\nabla{\mathcal{L}_{\bm{u}}(\bm{u}^{t})}$ for
notation simplification, then the projection step essentially aims to solve
the following quadratic programming problem:

$\begin{split}&\bm{u}^{t+1}=\operatorname*{arg\min}_{\bm{u}}\frac{1}{2}||\bm{u}-\hat{\bm{u}}||_{2}^{2},\\\
&\text{subject to}:\sum u_{i}=1,\quad u_{i}\geq 0.\end{split}$ (QP)

The Lagrangian associated with Eq. QP is as follows:

$\mathcal{L}(\bm{u},\lambda,\bm{v}):=\frac{1}{2}||\bm{u}-\hat{\bm{u}}||_{2}^{2}+\lambda(\bm{u}^{\top}\bm{1}-1)-\bm{v}^{\top}\hat{\bm{u}}\,,$

where $\lambda\in\mathbb{R},\bm{v}\in\mathbb{R}^{|S_{\text{all}}|}$ are the
Lagrange multipliers. The Karush-Kuhn-Tucker optimality conditions are
necessary and sufficient for solving Eq. QP, that is:

$\displaystyle\nabla_{\bm{u}}\mathcal{L}(\bm{u},\lambda,\bm{v})=\bm{u}-\hat{\bm{u}}+\lambda\bm{1}-\bm{v}=\bm{0},$
(3) $\displaystyle u_{i}\geq 0\,,$ (4) $\displaystyle\sum u_{i}-1=0\,,$ (5)
$\displaystyle v_{i}\geq 0\,,$ (6) $\displaystyle v_{i}u_{i}=0\,.$ (7)

Clearly, given any $\lambda$, if we set
$u_{i}=\max(\hat{u}_{i}-\lambda,0),\quad v_{i}=\max(\lambda-\hat{u}_{i},0),$
then Eqs. 3, 4, 6 and 7 can be satisfied. Therefore, the remaining problem
reduces to find a $\lambda$ that satisfies Eq. 5, i.e.,

$\sum u_{i}-1=\sum\max(\hat{u}_{i}-\lambda,0)-1=0\,.$

We employ the algorithm proposed in Held et al. (1974) to solve it, as
presented in Algorithm 3. Lastly, we select the
$\operatorname*{arg\max}_{j}(u^{\star}_{j})$ element in $S_{\text{all}}$ as
the attack sentence $S^{\prime}$.

Algorithm 3 Projection into simplex (Held et al., 1974)

Input:
$\bm{\hat{u}}:=\bm{u}^{t}+\eta\nabla{\mathcal{L}_{\bm{u}}(\bm{u}^{t})}\in\mathbb{R}^{|S_{\text{all}}|}$.
Sort $\bm{\hat{u}}$ such that
$\hat{u}_{1}\leq\hat{u}_{2}\leq\cdots\leq\hat{u}_{|S_{\text{all}}|}$.
Set
$J_{0}:=\max({J:\frac{-1+\sum_{i=J+1}^{|S_{\text{all}}|}\hat{u}_{i}}{|S_{\text{all}}|-J}>\hat{u}_{J}})$.
Calculate
$\lambda=\frac{-1+\sum_{i=J_{0}+1}^{|S_{\text{all}}|}\hat{u}_{i}}{|S_{\text{all}}|-J}.$
Set $u^{t+1}_{i}=\max(\hat{u_{i}}-\lambda,0)$.
Output: $\bm{u}^{t+1}$

###  D.2 Comparison between PGA-Charmer and query-based Charmer

In this section, we experimentally validate the efficiency of PGA-Charmer and
compare it against query-based Charmer, which is proposed in the main body.
The result in Table S25 shows that PGA-Charmer can efficiently reduce the
runtime as it does not require the forward pass over a mini-batch of sentences
after position selection. However, PGA-Charmer is worse than Charmer on other
metrics, e.g., ASR, Levenshtine distance and similarity are degraded. We
believe that combining the efficiency of PGA-Charmer and high ASR in -Charmer
holds promise for future research endeavors.

Table S25: Comparison between our PGA-Charmer and query-based Charmer (proposed in the main body) in BERT: The best method is highlighted in bold. While the PGA-Charmer strategy can noticeably improve the runtime, the ASR, Levenshtine distance and similarity are considerably degraded. | Method | ASR (%) $\bm{\uparrow}$ | $d_{\text{lev}}(S,S^{\prime})~{}\bm{\downarrow}$ | $\text{Sim}(S,S^{\prime})~{}\bm{\uparrow}$ | $\text{Time (s)}~{}\bm{\downarrow}$
---|---|---|---|---|---
AG-News | PGA-Charmer | $86.94$ | $7.33_{\pm(5.01)}$ | $0.87_{\pm(0.11)}$ | $\mathbf{8.15}_{\pm(7.04)}$
Charmer | $\mathbf{98.51}$ | $\mathbf{3.68}_{\pm(3.08)}$ | $\mathbf{0.95}_{\pm(0.06)}$ | $8.74_{\pm(11.10)}$
MNLI-m | PGA-Charmer | $99.05$ | $2.11_{\pm(1.53)}$ | $0.79_{\pm(0.19)}$ | $\mathbf{0.85}_{\pm(0.72)}$
Charmer | $\mathbf{100.00}$ | $\mathbf{1.14}_{\pm(0.42)}$ | $\mathbf{0.85}_{\pm(0.13)}$ | $1.45_{\pm(0.81)}$
QNLI | PGA-Charmer | $81.19$ | $3.46_{\pm(2.32)}$ | $0.89_{\pm(0.12)}$ | $\mathbf{5.15}_{\pm(4.60)}$
Charmer | $\mathbf{97.68}$ | $\mathbf{1.94}_{\pm(1.48)}$ | $\mathbf{0.94}_{\pm(0.07)}$ | $9.19_{\pm(9.60)}$
RTE | PGA-Charmer | $72.64$ | $\mathbf{1.53}_{\pm(1.45)}$ | $\mathbf{0.86}_{\pm(0.11)}$ | $\mathbf{0.65}_{\pm(0.71)}$
Charmer | $\mathbf{97.01}$ | $1.55_{\pm(1.42)}$ | $\mathbf{0.86}_{\pm(0.13)}$ | $2.50_{\pm(2.33)}$
SST-2 | PGA-Charmer | $97.52$ | $2.68_{\pm(1.82)}$ | $0.84_{\pm(0.16)}$ | $\mathbf{1.09}_{\pm(0.89)}$
Charmer | $\mathbf{100.00}$ | $\mathbf{1.47}_{\pm(0.74)}$ | $\mathbf{0.90}_{\pm(0.11)}$ | $1.27_{\pm(0.84)}$


# Piece-wise quadratic approximations of arbitrary error functions for fast
and robust machine learning

A.N. Gorban<EMAIL_ADDRESS>E.M. Mirkes<EMAIL_ADDRESS>A. Zinovyev
<EMAIL_ADDRESS>Department of Mathematics, University of Leicester,
Leicester, LE1 7RH, UK  Institut Curie, PSL Research University, Mines Paris
Tech, Inserm, U900, F-75005, Paris, France

###### Abstract

Most of machine learning approaches have stemmed from the application of
minimizing the mean squared distance principle, based on the computationally
efficient quadratic optimization methods. However, when faced with high-
dimensional and noisy data, the quadratic error functionals demonstrated many
weaknesses including high sensitivity to contaminating factors and
dimensionality curse. Therefore, a lot of recent applications in machine
learning exploited properties of non-quadratic error functionals based on
$L_{1}$ norm or even sub-linear potentials corresponding to quasinorms $L_{p}$
($0<p<1$). The back side of these approaches is increase in computational cost
for optimization. Till so far, no approaches have been suggested to deal with
arbitrary error functionals, in a flexible and computationally efficient
framework. In this paper, we develop a theory and basic universal data
approximation algorithms ($k$-means, principal components, principal manifolds
and graphs, regularized and sparse regression), based on piece-wise quadratic
error potentials of subquadratic growth (PQSQ potentials). We develop a new
and universal framework to minimize arbitrary sub-quadratic error potentials
using an algorithm with guaranteed fast convergence to the local or global
error minimum. The theory of PQSQ potentials is based on the notion of the
cone of minorant functions, and represents a natural approximation formalism
based on the application of min-plus algebra. The approach can be applied in
most of existing machine learning methods, including methods of data
approximation and regularized and sparse regression, leading to the
improvement in the computational cost/accuracy trade-off. We demonstrate that
on synthetic and real-life datasets PQSQ-based machine learning methods
achieve orders of magnitude faster computational performance than the
corresponding state-of-the-art methods, having similar or better approximation
accuracy.

###### keywords:

data approximation , nonquadratic potential , principal components ,
clustering , regularized regression , sparse regression

##  1 Introduction

Modern machine learning and artificial intelligence methods are
revolutionizing many fields of science today, such as medicine, biology,
engineering, high-energy physics and sociology, where large amounts of data
have been collected due to the emergence of new high-throughput computerized
technologies. Historically and methodologically speaking, many machine
learning algorithms have been based on minimizing the mean squared error
potential, which can be explained by tractable properties of normal
distribution and existence of computationally efficient methods for quadratic
optimization. However, most of the real-life datasets are characterized by
strong noise, long-tailed distributions, presence of contaminating factors,
large dimensions. Using quadratic potentials can be drastically compromised by
all these circumstances: therefore, a lot of practical and theoretical efforts
have been made in order to exploit the properties of non-quadratic error
potentials which can be more appropriate in certain contexts. For example,
methods of regularized and sparse regression such as lasso and elastic net
based on the properties of $L_{1}$ metrics [1, 2] found numerous applications
in bioinformatics [3], and $L_{1}$ norm-based methods of dimension reduction
are of great use in automated image analysis [4]. Not surprisingly, these
approaches come with drastically increased computational cost, for example,
connected with applying linear programming optimization techniques which are
substantially more expensive compared to mean squared error-based methods.

In practical applications of machine learning, it would be very attractive to
be able to deal with arbitrary error potentials, including those based on
$L_{1}$ or fractional quasinorms $L_{p}$ ($0<p<1$), in a computationally
efficient and scalable way. There is a need in developing methods allowing to
tune the computational cost/accuracy of optimization trade-off accordingly to
various contexts.

In this paper, we suggest such a universal framework able to deal with a large
family of error potentials. We exploit the fact that finding a minimum of a
piece-wise quadratic function, or, in other words, a function which is the
minorant of a set of quadratic functionals, can be almost as computationally
efficient as optimizing the standard quadratic potential. Therefore, if a
given arbitrary potential (such as $L_{1}$-based or fractional quasinorm-
based) can be approximated by a piece-wise quadratic function, this should
lead to relatively efficient and simple optimization algorithms. It appears
that only potentials of quadratic or subquadratic growth are possible in this
approach: however, these are the most usefull ones in data analysis. We
introduce a rich family of piece-wise quadratic potentials of subquadratic
growth (PQSQ-potentials), suggest general approach for their optimization and
prove convergence of a simple iterative algorithm in the most general case. We
focus on the most used methods of data dimension reduction and regularized
regression: however, potential applications of the approach can be much wider.

Data dimension reduction by constructing explicit low-dimensional
approximators of a finite set of vectors is one of the most fundamental
approach in data analysis. Starting from the classical data approximators such
as $k$-means [5] and linear principal components (PCA) [6], multiple
generalizations have been suggested in the last decades (self-organizing maps,
principal curves, principal manifolds, principal graphs, principal trees,
etc.)[7, 8] in order to make the data approximators more flexible and suitable
for complex data structures.

We solve the problem of approximating a finite set of vectors
${\vec{x}_{i}}\in R^{m},\,i=1,\ldots,N$ (data set) by a simpler object $L$
embedded into the data space, such that for each point $\vec{x}_{i}$ an
approximation error $err(\vec{x}_{i},L)$ function can be defined. We assume
this function in the form

$err(\vec{x}_{i},L)=\min_{y\in L}\sum_{k}u(x_{i}^{k}-y^{k}),$ (1)

where the upper $k=1,\ldots,m$ stands for the coordinate index, and $u(x)$ is
a monotonously growing symmetrical function, which we will be calling the
error potential. By data approximation we mean that the embedment of $L$ in
the data space minimizes the error

$\sum_{i}err(\vec{x}_{i},L)\rightarrow\min.$

Note that our definition of error function is coordinate-wise (it is a sum of
error potential over all coordinates).

The simplest form of the error potential is quadratic $u(x)=x^{2}$, which
leads to the most known data approximators: mean point ($L$ is a point),
principal points ($L$ is a set of points) [9], principal components ($L$ is a
line or a hyperplane) [6]. In more advanced cases, $L$ can posses some regular
properties leading to principal curves ($L$ is a smooth line or spline) [10],
principal manifolds ($L$ is a smooth low-dimensional surface) and principal
graphs (eg., $L$ is a pluri-harmonic graph embedment) [11, 7].

There exist multiple advantages of using quadratic potential $u(x)$, because
it leads to the most computationally efficient algorithms usually based on the
splitting schema, a variant of expectation-minimization approach [7]. For
example, $k$-means algorithm solves the problem of finding the set of
principal points and the standard iterative Singular Value Decomposition finds
principal components. However, quadratic potential is known to be sensitive to
outliers in the data set. Also, purely quadratic potentials can suffer from
the curse of dimensionality, not being able to robustly discriminate ‘close’
and ‘distant’ point neighbours in a high-dimensional space [12].

There exist several widely used ideas for increasing approximator’s robustness
in presence of strong noise in data such as: (1) using medians instead of mean
values, (2) substituting quadratic norm by $L_{1}$ norm (e.g. [13, 14]), (3)
outliers exclusion or fixed weighting or iterative reweighting during
optimizing the data approximators (e.g. [15, 16, 17]), and (4) regularizing
the PCA vectors by $L_{1}$ norm [18, 19, 20]. In some works, it was suggested
to utilize ‘trimming’ averages, e.g. in the context of the $k$-means
clustering or some generalizations of PCA [21, 14]). In the context of
regression, iterative reweighting is exploited to mimic the properties of
$L_{1}$ norm [22]. Several algorithms for constructing PCA with $L_{1}$ norm
have been suggested [23, 24, 25] and systematically benchmarked [26, 27]. Some
authors go even beyond linear metrics and suggests that fractional quasinorms
$L_{p}$ ($0<p<1$) can be more appropriate in high-dimensional data
approximation [12].

However, most of the suggested approaches exploiting properties of non-
quadratic metrics either represent useful but still arbitrary heuristics or
are not sufficiently scalable. The standard approach for minimizing
$L_{1}$-based norm consists in solving a linear programming task. Despite
existence of many efficient linear programming optimizer implementations, by
their nature these computations are much slower than the iterative methods
used in the standard SVD algorithm or $k$-means.

In this paper, we provide implementations of the standard data approximators
(mean point, $k$-means, principal components) using a PQSQ potential. As an
other application of PQSQ-based framework in machine learning, we develop
PQSQ-based regularized and sparse regression (imitating the properties of
lasso and elastic net).

##  2 Piecewise quadratic potential of subquadratic growth (PQSQ)

###  2.1 Definition of the PQSQ potential

Let us split all non-negative numbers $x\in R_{\geq 0}$ into $p+1$ non-
intersecting intervals
$R_{0}=[0;r_{1}),R_{1}=[r_{1};r_{2}),\ldots,R_{k}=[r_{k};r_{k+1}),\ldots,R_{p}=[r_{p};\infty)$,
for a set of thresholds $r_{1}<r_{2}<\ldots<r_{p}$. For convenience, let us
denote $r_{0}=0,r_{p+1}=\infty$. Piecewise quadratic potential is a continuous
monotonously growing function $u(x)$ constructed from pieces of centered at
zero parabolas $y=b_{k}+a_{k}x^{2}$, defined on intervals
$x\in[r_{k},r_{k+1})$, satisfying $y(r_{i})=f(r_{i})$ (see Figure 1):

$u(x)=b_{k}+a_{k}x^{2},\mbox{if }r_{k}\leq|x|<r_{k+1},\,k=0,\ldots,p,$ (2)

$a_{k}=\frac{f(r_{k})-f(r_{k+1})}{r_{k}^{2}-r_{k+1}^{2}},$ (3)

$b_{k}=\frac{f(r_{k+1})r_{k}^{2}-f(r_{k})r_{k+1}^{2}}{r_{k}^{2}-r_{k+1}^{2}},$
(4)

where $f(x)$ is a majorating function, which is to be approximated (imitated)
by $u(x)$. For example, in the simplest case $f(x)$ can be a linear function:
$f(x)=x$, in this case, $\sum_{k}u(x^{k})$ will approximate the $L_{1}$-based
error function.

Note that accordingly to (3,4), $b_{0}=0,a_{p}=0,b_{p}=f(r_{p})$. Therefore,
the choice of $r_{p}$ can naturally create a ‘trimmed’ version of error
potential $u(x)$ such that some data points (outliers) do not have any
contribution to the gradient of $u(x)$, hence, will not affect the
optimization procedure. However, this set of points can change during
minimizaton of the potential.

The condition of subquadratic growth consists in the requirement $a_{k+1}\leq
a_{k}$ and $b_{k+1}\geq b_{k}$. To guarantee this, the following simple
condition on $f(x)$ should be satisfied:

$f^{\prime}>0,\>\>\>f^{\prime\prime}x\leq f^{\prime}.$ (5)

Therefore, $f(x)$ is a monotonic concave function of $q=x^{2}$:

$\frac{d^{2}f(\sqrt{q})}{dq^{2}}=\frac{1}{4x^{2}}f^{\prime\prime}(x)-\frac{1}{4x^{3}}f^{\prime}(x)\leq
0.$

In particular, $f(x)$ should grow not faster than any parabola
$ax^{2}+c,\,c>0$, which is tangent to $f(x)$.

Figure 1: Trimmed piecewise quadratic potential of subquadratic growth $u(x)$
(solid blue line) defined for the majorating function $f(x)$ (red dashed line)
and several thresholds $r_{k}$. Dotted lines show the parabolas which
fragments are used to construct $u(x)$. The last parabola is flat ($a_{p}=0$)
which corresponds to trimmed potential.

###  2.2 Basic approach for optimization

In order to use the PQSQ potential in an algorithm, a set of $p$ interval
thresholds $r_{s}^{k},\,s=1,\ldots,p$ for each coordinate $k=1,\ldots,m$
should be provided. Matrices of $a$ and $b$ coefficients defined by (3,4)
based on interval definitions: $a_{s}^{k},b_{s}^{k}$, $s=0,\ldots,p$,
$k=1,\ldots,m$ are computed separately for each coordinate $k$.

Minimization of PQSQ-based functional consists in several basic steps which
can be combined in an algorithm:

1) For each coordinate $k$, split all data point indices into non-overlapping
sets $\mathcal{R}_{s}^{k}$:

$\mathcal{R}_{s}^{k}=\\{i:r_{s}^{k}\leq|x_{i}^{k}-\beta^{k}_{i}|<r_{s+1}^{k}\\},\,s=0,\ldots,p,$
(6)

where $\bf{\beta}$ is a matrix which depends on the nature of the algorithm.

2) Minimize PQSQ-based functional where each set of points
$\\{x_{i\in\mathcal{R}_{s}^{k}}\\}$ contributes to the functional
quadratically with coefficient $a_{s}^{k}$. This is a quadratic optimization
task.

3) Repeat (1)-(2) till convergence.

##  3 General theory of the piece-wise convex potentials as the cone of
minorant functions

In order to deal in most general terms with the data approximation algorithms
based on PQSQ potentials, let us consider a general case where a potential can
be constructed from a set of functions $\\{q_{i}(x)\\}$ with only two
requirements: 1) that each $q_{i}(x)$ has a (local) minimum; 2) that the whole
set of all possible $q_{i}(x)$s forms a cone. In this case, instead of the
operational definition (2) it is convenient to define the potential $u(x)$ as
the minorant function for a set of functions as follows. For convenience, in
this section, $x$ will notify a vector $\vec{x}\in R^{m}$.

Let us consider a generating cone of functions $Q$. We remind that the
definition of a cone implies that for any $q(x)\in Q,p(x)\in Q$, we have
$\alpha q(x)+\beta p(x)\in Q$, where $\alpha\geq 0,\beta\geq 0$.

For any finite set of functions

${q_{1}(x)\in Q,q_{2}(x)\in Q,\ldots,q_{s}(x)\in Q},$

we define the minorant function (Figure 2):

$u_{q_{1},q_{2},\ldots,q_{s}}(x)=\min(q_{1}(x),q_{2}(x),\ldots,q_{s}(x)).$ (7)

It is convinient to introduce a multiindex

$I_{q_{1},q_{2},\ldots,q_{s}}(x)$

indicating which particular function(s) $q_{i}$ corresponds to the value of
$u(x)$, i.e.

$I_{q_{1},q_{2},\ldots,q_{s}}(x)=\\{i|u_{q_{1},q_{2},\ldots,q_{s}}(x)=q_{i}(x)\\}.$
(8)

For a cone $Q$ let us define a set of all possible minorant functions
$\mathbb{M}(Q)$

$\begin{split}\mathbb{M}(Q)=\\{u_{q_{i_{1}},q_{i_{2}},\ldots,q_{i_{n}}}|q_{i_{1}}\in
Q,q_{i_{2}}\in Q,\\\ q_{i_{n}}\in Q,\,n=1,2,3,\dots\\}.\end{split}$ (9)

###### Proposition 1

$\mathbb{M}(Q)$ is a cone.

  * Proof

For any two minorant functions

$u_{q_{i_{1}},q_{i_{2}},\ldots,q_{i_{k}}},u_{q_{j_{1}},q_{j_{2}},\ldots,q_{j_{s}}}\in\mathbb{M}(Q)$

we have

$\begin{split}\alpha u_{q_{i_{1}},q_{i_{2}},\ldots,q_{i_{k}}}+\beta
u_{q_{j_{1}},q_{j_{2}},\ldots,q_{j_{s}}}=\\\ u_{\\{\alpha q_{i_{p}}+\beta
q_{j_{r}}\\}}\in\mathbb{M}(Q),\\\ p=1,\dots,k,\,r=1,\dots,s,\end{split}$ (10)

where ${\\{\alpha q_{i_{p}}+\beta q_{j_{r}}\\}}$ is a set of all possible
linear combinations of functions from
$\\{q_{i_{1}},q_{i_{2}},\ldots,q_{i_{k}}\\}$ and
$\\{q_{j_{1}},q_{j_{2}},\ldots,q_{j_{s}}\\}$.

###### Proposition 2

Any restriction of $\mathbb{M}(Q)$ onto a linear manifold $L$ is a cone.

  * Proof

Let us denote $q(x)|_{L}$ a restriction of $q(x)$ function onto $L$, i.e.
$q(x)|_{L}=\\{q(x)|x\in L\\}$. $q(x)|_{L}$ is a part of $Q$. Set of all
$q(x)|_{L}$ forms a restriction $Q|_{L}$ of $Q$ onto $L$. $Q|_{L}$ is a cone,
hence, $\mathbb{M}(Q)|_{L}=\mathbb{M}(Q|_{L})$ is a cone (Proposition 1).

  * Definition

Splitting algorithm minimizing

$u_{q_{1},q_{2},\ldots,q_{n}}(x)$

is defined as Algorithm 1.

Algorithm 1 Finding local minimum of a minorant function
$u_{q_{1},q_{2},\ldots,q_{n}}(x)$

1:procedure Minimizing minorant function

2: $\textit{initialize }x\leftarrow x_{0}$

3:  _repeat until stopping criterion has been met_ :

4: $\textit{compute multiindex }I_{q_{1},q_{2},\ldots,q_{s}}(x)$

5: for all  $i\in I_{q_{1},q_{2},\ldots,q_{s}}(x)$

6: $x_{i}=\operatorname*{arg\,min}\,q_{i}(x)$

7: end for

8: $\textit{select optimal }x_{i}:$

9: $x_{opt}\leftarrow\operatorname*{arg\,min}_{x_{i}}u(x_{i})$

10: $x\leftarrow x_{opt}$

11: stopping criterion: check if the multiindex
$I_{q_{1},q_{2},\ldots,q_{s}}(x)$ does not change compared to the previous
iterationf

###### Theorem 3.1

Splitting algorithm (Algorithm 1) for minimizing
$u_{q_{1},q_{2},\ldots,q_{n}}(x)$ converges in a finite number of steps.

  * Proof

Since the set of functions $\\{q_{1},q_{2},\ldots,q_{n}\\}$ is finite then we
only have to show that at each step the value of the function
$u_{q_{1},q_{2},\ldots,q_{n}}(x)$ can not increase. For any $x$ and the value
$x^{\prime}=\operatorname*{arg\,min}q_{i}(x)$ for $i\in
I_{q_{1},q_{2},\ldots,q_{s}}(x)$ we can have only two cases:

(1) Either
$I_{q_{1},q_{2},\ldots,q_{s}}(x)=I_{q_{1},q_{2},\ldots,q_{s}}(x^{\prime})$
(convergence, and in this case $q_{i^{\prime}}(x^{\prime})=q_{i}(x^{\prime})$
for any $i^{\prime}\in I_{q_{1},q_{2},\ldots,q_{s}}(x^{\prime})$);

(2) Or
$u_{q_{1},q_{2},\ldots,q_{n}}(x^{\prime})<u_{q_{1},q_{2},\ldots,q_{n}}(x)$
since, accordingly to the definition (7),
$q_{i^{\prime}}(x^{\prime})<q_{i}(x)$, for any $i^{\prime}\in
I_{q_{1},q_{2},\ldots,q_{s}}(x^{\prime}),i\in I_{q_{1},q_{2},\ldots,q_{s}}(x)$
(see Figure 2).

Note that in Algorithm 1 we do not specify exactly the way to find the local
minimum of $q_{i}(x)$. To be practical, the cone $Q$ should contain only
functions for which finding a local minimum is fast and explicit. Evident
candidates for this role are positively defined quadratic functionals
$q(x)=q_{0}+(\vec{q_{1}},x)+(x,\mathbb{Q}_{2}x)$, where $\mathbb{Q}_{2}$ is a
positively defined symmetric matrix. Any minorant function (7) constructed
from positively defined quadratic functions will automatically provide
subquadratic growth, since the minorant can not grow faster than any of the
quadratic forms by which it is defined.

Operational definition of PQSQ given above (2), corresponds to a particular
form of the quadratic functional, with $\mathbb{Q}_{2}$ being diagonal matrix.
This choice corresponds to the coordinate-wise definition of data
approximation error function (1) which is particularly simple to minimize.
This circumstance is used in Algorithms 2,3.

Figure 2: Optimization of a one-dimensional minorant function $u(x)$, defined
by three functions $q_{1}(x),q_{2}(x),q_{3}(x)$ each of which has a local
minimum. Each optimization step consists in determining which
$q_{I(x)}(x)=u(x)$ and making a step into the local minimum of $q_{I(x)}$.

##  4 Commonly used data approximators with PQSQ potential

###  4.1 Mean value and $k$-means clustering in PQSQ approximation measure

Mean vector $\bar{X}_{L}$ for a set of vectors $X=\\{x_{i}^{k}\\}$,
$i=1,\ldots,N,\,k=1,\ldots,m$ and an approximation error defined by potential
$f(x)$ can be defined as a point minimizing the mean error potential for all
points in $X$:

$\sum_{i}\sum_{k}f(x_{i}^{k}-\bar{X}^{k})\rightarrow\min.$ (11)

For Euclidean metrics $L_{2}$ ($f(x)=x^{2}$) it is the usual arithmetric mean.

For $L_{1}$ metrics ($f(x)=|x|$), (11) leads to the implicit equation
$\\#(x_{i}^{k}>\bar{X}^{k})=\\#(x_{i}^{k}<\bar{X}^{k})$, where $\\#$ stands
for the number of points, which corresponds to the definition of median. This
equation can have a non-unique solution in case of even number of points or
when some data point coordinates coincide: therefore, definition of median is
usually accompanied by heuristics used for breaking ties, i.e. to deal with
non-uniquely defined rankings. This situation reflects the general situation
of existence of multiple local minimuma and possible non-uniqueness of global
minimum of (11) (Figure 3).

For PQSQ approximation measure (2) it is difficult to write down an explicit
formula for computing the mean value corresponding to the global minimum of
(11). In order to find a point $\bar{X}_{PQSQ}$ minimizing mean $PQSQ$
potential, a simple iterative algorithm can be used (Algorithm 2). The
suggested algorithm converges to the local minimum which depends on the
initial point approximation.

Algorithm 2 Computing PQSQ mean value

1:procedure PQSQ Mean Value

2: $\textit{define intervals }r_{s}^{k},s=0,\ldots,p,k=1,\ldots,m$

3: $\textit{compute coefficients }a_{s}^{k}$

4: $\textit{initialize }\bar{X}_{PQSQ}\linebreak\textit{eg., by arithmetic
mean}$

5:  _repeat till convergence of $\bar{X}_{PQSQ}$_:

6: for each  coordinate  $k$

7: define sets of indices
$\begin{split}\mathcal{R}_{s}^{k}=\\{i:r_{s}^{k}\leq|x_{i}^{k}-\bar{X}_{PQSQ}^{k}|<r_{s+1}^{k}\\},\\\
s=0,\dots,p\end{split}$

8: compute new approximation for  $\bar{X}_{PQSQ}$:

9:
$\bar{X}_{PQSQ}^{k}\leftarrow\frac{\sum_{s=1,\ldots,p}a_{s}^{k}\sum_{i\in\mathcal{R}_{s}^{k}}x_{i}^{k}}{\sum_{s=1,\ldots,p}a_{s}^{k}|\mathcal{R}_{s}^{k}|}$

10: end for

11: goto _repeat till convergence_

Based on the PQSQ approximation measure and the algorithm for computing the
PQSQ mean value (Algorithm 2), one can construct the PQSQ-based $k$-means
clustering procedure in the usual way, splitting estimation of cluster
centroids given partitioning of the data points into $k$ disjoint groups, and
then re-calculating the partitioning using the PQSQ-based proximity measure.

Figure 3: Minimizing the error to a point (finding the mean value) for a set
of 4 points (shown by black circles). Solid red line corresponds to
$L_{1}$-based error. Thing blue lines correspond to PQSQ error potential
imitating the $L_{1}$-based error. Several choices of PQSQ potential for
different numbers of intervals (indicated by a number put on top of the line)
is illustrated. On the right panel a zoom of a particular region of the left
plot is shown. Neither function ($L_{1}$-based or PQSQ-based) possesses a
unique local minimum. Moreover, $L_{1}$-based error function has infinite
number of points corresponding to the global minimum (any number between 3 and
4), while PQSQ error function has several local minimuma in [3;4] interval
which exact positions are sensitive to the concrete choice of PQSQ parameters
(interval definitions).

###  4.2 Principal Component Analysis (PCA) in PQSQ metrics

Accordingly to the classical definition of the first principal component, it
is a line best fit to the data set $X$ [6]. Let us define a line in the
parametric form $\vec{y}=\vec{V}\nu+\vec{\delta}$, where $\nu\in R^{1}$ is the
parameter. Then the first principal component will be defined by vectors
$\vec{V},\vec{\delta}$ satisfying

$\sum_{i}\sum_{k}u(x_{i}^{k}-V^{k}\nu_{i}-\delta^{k})\rightarrow\min,$ (12)

where

$\nu_{i}=\arg\min_{s}\sum_{k}u(x_{i}^{k}-V^{k}s-\delta^{k}).$ (13)

The standard first principal component (PC1) corresponds to $u(x)=x^{2}$ when
the vectors $\vec{V},\vec{\delta}$ can be found by a simple iterative
splitting algorithm for Singular Value Decomposition (SVD). If $X$ does not
contain missing values then $\vec{\delta}$ is the vector of arithmetic mean
values. By contrast, computing $L_{1}$-based principal components ($u(x)=|x|$)
represents a much more challenging optimization problem [25]. Several
approximative algorithms for computing $L_{1}$-norm PCA have been recently
suggested and benchmarked [23, 24, 25, 26, 27]. To our knowledge, there have
not been a general efficient algorithm suggested for computing PCA in case of
arbitrary approximation measure for some monotonous function $u(x)$.

Computing PCA based on PQSQ approximation error is only slightly more
complicated than computing the standard $L_{2}$ PCA by SVD. Here we provide a
pseudo-code (Algorithm 3) of a simple iterative algorithm (similar to
Algorithm 2) with guaranteed convergence (see Section 3).

Algorithm 3 Computing PQSQ PCA

1:procedure PQSQ First Principal Component

2: $\textit{define intervals }r_{s}^{k},s=0,\ldots,p,\,k=1,\ldots,m$

3: $\textit{compute coefficients }a_{s}^{k}$

4: $\vec{\delta}\leftarrow\bar{X}_{PQSQ}$

5: $\textit{initialize }\vec{V}\textit{ : eg., by $L_{2}$-based PC1}$

6: $\textit{initialize }\\{\nu_{i}\\}\textit{ : eg., by
}\linebreak$$\nu_{i}=\frac{\sum_{k}V^{k}(x_{i}^{k}-\delta^{k})}{\sum_{k}(V^{k})^{2}}$

7:  _repeat till convergence of $\vec{V}$_:

8: $\textit{normalize }\vec{V}\textit{ :
}\vec{V}\leftarrow\frac{\vec{V}}{\|\vec{V}\|}$

9: for each  coordinate  $k$

10: define sets of indices
$\begin{split}\mathcal{R}_{s}^{k}=\\{i:r_{s}^{k}\leq|x_{i}^{k}-V^{k}\nu_{i}-\delta^{k}|<r_{s+1}^{k}\\},\\\
s=0,\ldots,p\end{split}$

11: end for

12: for each  data point  $i$ and coordinate  $k$

13: find all $s_{i,k}$ such that $i\in\mathcal{R}_{s_{i,k}}^{k}$

14: if all $a^{k}_{s_{i,k}}=0$ then $\nu^{\prime}_{i}\leftarrow 0$ else

15:
$\nu^{\prime}_{i}\leftarrow\frac{\sum_{k}a_{s_{i,k}}^{k}V^{k}(x_{i}^{k}-\delta^{k})}{\sum_{k}a_{s_{i,k}}^{k}(V^{k})^{2}}$

16: end for

17: for each  coordinate  $k$
$V^{k}\leftarrow\frac{\sum_{s}a_{s}^{k}\sum_{i\in\mathcal{R}_{s}^{k}}(x_{i}^{k}-\delta^{k})\nu_{i}}{\sum_{s}a_{s}^{k}\sum_{i\in\mathcal{R}_{s}^{k}}(\nu_{i})^{2}}$

18: end for

19: for each  $i$ :

20: $\nu_{i}\leftarrow\nu^{\prime}_{i}$

21: end for

22: goto _repeat till convergence_

Computation of second and further principal components follows the standard
deflation approach: projections of data points onto the previously computed
component are subtracted from the data set, and the algorithm is applied to
the residues. However, as it is the case in any non-quadratic metrics, the
resulting components can be non-orthogonal or even not invariant with respect
to the dataset rotation. Moreover, unlike $L_{2}$-based principal components,
the Algorithm 3 does not always converge to a unique global minimum; the
computed components can depend on the initial estimate of $\vec{V}$. The
situation is somewhat similar to the standard $k$-means algorithm. Therefore,
in order to achieve the least possible approximation error to the linear
subspace, $\vec{V}$ can be initialized randomly or by data vectors
$\vec{x}_{i}$ many times and the deepest in PQSQ approximation error (1)
minimum should be selected.

How does the Algorithm 1 serve a more abstract version of the Algorithms 2,3?
For example, the ‘variance’ function
$m(\vec{x})=\frac{1}{N}\sum_{j}u(\vec{x}_{j}-\vec{x})$ to be minimized in
Algorithm 2 uses the generating functions in the form
$Q=\\{b_{ji}^{k}+\sum_{k}a_{ji}^{k}(x^{k}-x_{j}^{k})^{2}\\}$, where $i$ is the
index of the interval in (2). Hence, $m(x)$ is a minorant function, belonging
to the cone $\mathbb{M}(Q)$, and must converge (to a local minimum) in a
finite number of steps accordingly to Theorem 3.1.

###  4.3 Nonlinear methods: PQSQ-based Principal
Graphs and Manifolds

In a series of works, the authors of this article introduced a family of
methods for constructing principal objects based on graph approximations
(e.g., principal curves, principal manifolds, principal trees), which allows
constructing explicit non-linear data approximators (and, more generally,
approximators with non-trivial topologies, suitable for approximating, e.g.,
datasets with branching or circular topology) [28, 29, 30, 31, 11, 8, 7, 32].
The methodology is based on optimizing a piece-wise quadratic elastic energy
functional (see short description below). A convenient graphical user
interface was developed with implementation of some of these methods [33].

Let $G$ be a simple undirected graph with set of vertices $Y$ and set of edges
$E$. For $k\geq 2$ a $k$-star in $G$ is a subgraph with $k+1$ vertices
$y_{0,1,\ldots,k}\in Y$ and $k$ edges $\\{(y_{0},y_{i})\ |\
i=1,\ldots,k\\}\subset E$. Suppose for each $k\geq 2$, a family $S_{k}$ of
$k$-stars in $G$ has been selected. We call a graph $G$ with selected families
of $k$-stars $S_{k}$ an elastic graph if, for all $E^{(i)}\in E$ and
$S^{(j)}_{k}\in S_{k}$, the correspondent elasticity moduli $\lambda_{i}>0$
and $\mu_{kj}>0$ are defined. Let $E^{(i)}(0),E^{(i)}(1)$ be vertices of an
edge $E^{(i)}$ and $S^{(j)}_{k}(0),\ldots,S^{(j)}_{k}(k)$ be vertices of a
$k$-star $S^{(j)}_{k}$ (among them, $S^{(j)}_{k}(0)$ is the central vertex).

For any map $\phi:Y\to R^{m}$ the energy of the graph is defined as

$\begin{split}U^{\phi}{(G)}:=\sum_{E^{(i)}}\lambda_{i}\left\|\phi(E^{(i)}(0))-\phi(E^{(i)}(1))\right\|^{2}+\\\
+\sum_{S^{(j)}_{k}}\mu_{kj}\left\|\sum_{i=1}^{k}\phi(S^{(j)}_{k}(i))-k\phi(S^{(j)}_{k}(0))\right\|^{2}.\end{split}$

For a given map $\phi:Y\to R^{m}$ we divide the dataset $D$ into node
neighborhoods $K^{y},\,y\in Y$. The set $K^{y}$ contains the data points for
which the node $\phi(y)$ is the closest one in $\phi(y)$. The energy of
approximation is:

$U^{\phi}_{A}(G,D)=\sum_{y\in Y}\sum_{x\in K^{y}}w(x)\|x-\phi(y)\|^{2},$ (14)

where $w(x)\geq 0$ are the point weights. Simple and fast algorithm for
minimization of the energy

$U^{\phi}=U^{\phi}_{A}(G,D)+U^{\phi}{(G)}$ (15)

is the splitting algorithm, in the spirit of the classical $k$-means
clustering: for a given system of sets $\\{K^{y}\ |\ y\in Y\\}$ we minimize
$U^{\phi}$ (optimization step, it is the minimization of a positive quadratic
functional), then for a given $\phi$ we find new $\\{K^{y}\\}$ (re-
partitioning), and so on; stop when no change.

Application of PQSQ-based potential is straightforward in this approach. It
consists in replacing (14) with

$U^{\phi}_{A}(G,D)=\sum_{y\in Y}\sum_{x\in
K^{y}}w(x)\sum_{k}u(x^{k}-\phi(y^{k})),$

where $u$ is a chosen PQSQ-based error potential. Partitioning of the dataset
into $\\{K^{y}\\}$ can be also based on calculating the minimum PQSQ-based
error to $y$, or can continue enjoying nice properties of $L_{2}$-based
distance calculation.

##  5 PQSQ-based regularized regression

###  5.1 Regularizing linear regression with PQSQ potential

One of the major application of non-Euclidean norm properties in machine
learning is using non-quadratic terms for penalizing large absolute values of
regression coefficients [1, 2]. Depending on the chosen penalization term, it
is possible to achieve various effects such as sparsity or grouping
coefficients for redundant variables. In a general form, regularized
regression solves the following optimization problem

$\frac{1}{N}\sum_{i=1}^{N}\left(y_{i}-\sum_{k=1}^{m}\beta^{k}x_{i}^{k}\right)^{2}+\lambda
f(\vec{\beta})\rightarrow\min,$ (16)

where $N$ is the number of observations, $m$ is the number of independent
variables in the matrix $\\{x_{i}^{k}\\}$, $\\{y_{i}\\}$ are dependent
variables (to be predicted), $\lambda$ is an internal parameter controlling
the strength of regularization (penalty on the amplitude of regression
coefficients $\beta$), and $f(\vec{z})$ is the regularizer function, which is
$f(\vec{z})=\|\vec{z}\|_{L2}^{2}$ for ridge regression,
$f(\vec{z})=\|\vec{z}\|_{L1}$ for lasso and
$f(\vec{z})=\frac{1-\alpha}{\alpha}\|\vec{z}\|_{L2}^{2}+\alpha\|\vec{z}\|_{L1}$
for elastic net methods correspondingly.

Here we suggest to imitate $f(x)$ with a $PQSQ$ potential function, i.e.
instead of (16) solving the problem

$\frac{1}{N}\sum_{i=1}^{N}\left(y_{i}-\sum_{k=1}^{m}\beta^{k}x_{i}^{k}\right)^{2}+\lambda\sum_{k=1}^{m}u(\beta^{k})\rightarrow\min,$
(17)

where $u(\beta)$ is a PQSQ potential imitating arbitrary subquadratic
regression regularization penalty.

Solving (17) is equivalent to iteratively solving a system of linear equations

$\begin{split}&\frac{1}{N}\sum_{k=1}^{m}\beta^{k}\sum_{i=1}^{N}x_{i}^{k}x_{i}^{j}+\lambda
a_{I(\beta^{j})}\beta^{j}\\\
&\;\;\;\;\;=\sum_{i=1}^{N}y_{i}x_{i}^{j},\,j=1,\dots,m,\end{split}$ (18)

where $a_{I(\beta^{j})}$ constant (where $I$ index is defined from
$r_{I}\leq\beta^{j}<r_{I+1}$) is computed accordingly to the definition of
$u(x)$ function (see (3)), given the estimation of $\beta^{k}$ regression
coefficients at the current iteration. In practice, iterating (18) converges
in a few iterations, therefore, the algorithm can work very fast and
outperform the widely used least angle regression algorithm for solving (16)
in case of $L_{1}$ penalties.

###  5.2 Introducing sparsity by ‘black hole’ trick

Any PQSQ potential $u(x)$ is characterized by zero derivative for $x=0$ by
construction: $u^{\prime}(x)|_{x=0}=0$, which means that the solution of (17)
does not have to be sparse for any $\lambda$. Unlike pure $L_{1}$-based
penalty, the coefficients of regression diminish with increase of $\lambda$
but there is nothing to shrink them to exact zero values, similar to the ridge
regression. However, it is relatively straightforward to modify the algorithm,
to achieve sparsity of the regression solution. The ’black hole’ trick
consists in eliminating from regression training after each iteration (18) all
regression coefficients $\beta^{k}$ smaller by absolute value than a given
parameter $\epsilon$ (’black hole radius’). Those regression coefficients
which have passed the ’black hole radius’ are put to zero and do not have any
chance to change their values in the next iterations.

Figure 4: ’Black hole trick’ for introducing sparsity into the PQSQ-based
regularized regression. Here PQSQ function imitates $L_{1}$ norm (for
illustration only three intervals are used to define PQSQ function). Black
hole trick consists in introducing an $\epsilon$ zone (hatched territory on
the plot) of the potential in the vicinity of zero such that any coefficient
of regression falling down into this zone is set to zero and eliminated from
further learning. It is convenient to define $\epsilon$ as the midst of the
smallest interval as it is shown in this plot.

The optimal choice of $\epsilon$ value requires a separate study. From general
considerations, it is preferable that the derivative
$u^{\prime}(x)|_{x=\epsilon}$ would not be very close to zero. As a pragmatic
choice for the numerical examples in this article, we define $\epsilon$ as the
midst of the smallest interval in the definition of PQSQ potential (see Figure
4), i.e. $\epsilon=r_{1}/2$, which guarantees far from zero
$u^{\prime}(x)|_{x=\epsilon}$. It might happen that this value of $\epsilon$
would collapse all $\beta^{k}$ to zero even without regularization (i.e., with
$\lambda=0$). In this case, the ’black hole radius’ is divided by half
$\epsilon\leftarrow\epsilon/2$ and it is checked that for $\lambda=0$ the
iterations would leave at list half of the regression coefficients. If it is
not the case then the process of diminishing the ’black hole radius’ repeated
recursively till meeting the criterion of preserving the majority of
regression coefficients. In practice, it requires only few (very fast)
additional iterations of the algorithm.

As in the lasso methodology, the problem (17) is solved for a range of
$\lambda$ values, calibrated such that the minimal $\lambda$ would select the
maximum number of regression variables, while the maximum $\lambda$ value
would lead to the most sparse regression (selecting only one single non-zero
regression coefficient).

##  6 Numerical examples

###  6.1 Practical choices of parameters

The main parameters of PQSQ are (a) majorating function $f(x)$ and (b)
decomposition of each coordinate range into $p+1$ non-overlapping intervals.
Depending on these parameters, various approximation error properties can be
exploited, including those providing robustness to outlier data points.

When defining the intervals $r_{j},\,j=1,\dots,p$, it is desirable to achieve
a small difference between $f(\Delta x)$ and $u(\Delta x)$ for expected
argument values $\Delta x$ (differences between an estimator and the data
point), and choose the suitable value of the potential trimming threshold
$r_{p}$ in order to achieve the desired robustness properties. If no trimming
is needed, then $r_{p}$ should be made larger than the maximum expected
difference between coordinate values (maximum $\Delta x$).

In our numerical experiments we used the following definition of intervals.
For any data coordinate $k$, we define a characteristic difference $D^{k}$,
for example

$D^{k}=\alpha_{scale}(max_{i}(x_{i}^{k})-min_{i}(x_{i}^{k})),$ (19)

where $\alpha_{scale}$ is a scaling parameter, which can be put at 1 (in this
case, the approximating potential will not be trimmed). In case of existence
of outliers, for defining $D^{k}$, instead of amplitude one can use other
measures such as the median absolute deviation (MAD):

$D^{k}=\alpha_{scale}median_{i}(|x_{i}^{k}-median(\\{x_{i}^{k}\\})|);$ (20)

in this case, the scaling parameter should be made larger, i.e.
$\alpha_{scale}=10$, if no trimming is needed.

After defining $D^{k}$ we use the following definition of intervals:

$r_{j}^{k}=D^{k}\frac{j^{2}}{p^{2}}\,,\,j=0,\ldots,p.$ (21)

More sophisticated approaches are also possible to apply such as, given the
number of intervals $p$ and the majorating function $f(x)$, choose
$r_{j},\,j=1,\dots,p$ in order to minimize the maximal difference

$d=\max_{x}|f(x)-u(x)|\rightarrow\min.$

The calculation of intervals is straightforward for a given value of $d$ and
many smooth concave functions $f(x)$ like $f(x)=x^{p}$ ($0<p\leq 1$) or
$f(x)=\ln(1+x)$.

###  6.2 Implementation

We provide Matlab implementation of PQSQ approximators (in particular, PCA)
together with the Matlab and R code used to generate the example figures in
this article at ‘PQSQ-DataApproximators’ GitHub
repository111https://github.com/auranic/PQSQ-DataApproximators and Matlab
implementation of PQSQ-based regularized regression with build-in imitations
of $L_{1}$ (lasso-like) and $L_{1}\&L_{2}$ mixture (elastic net-like)
penalties at ‘PQSQ-regularized-regression’ GitHub
repository222https://github.com/Mirkes/PQSQ-regularized-regression/. The Java
code implementing elastic graph-based non-linear approximator implementations
is available from the authors on request.

###  6.3 Motivating example: dense two-cluster distribution contaminated by
sparse noise

We demonstrate the value of PQSQ-based computation of $L_{1}$-based PCA by
constructing a simple example of data distribution consisting of a dense two-
cluster component superimposed with a sparse contaminating component with
relatively large variance whose co-variance does not coincide with the dense
signal (Figure 5). We study the ability of PCA to withstand certain level of
sparse contamination and compare it with the standard $L_{2}$-based PCA. In
this example, without noise the first principal component coincides with the
vector connecting the two cluster centers: hence, it perfectly separates them
in the projected distribution. Noise interferes with the ability of the first
principal component to separate the clusters to the degree when the first
principal component starts to match the principal variance direction of the
contaminating distribution (Figure 5A,B). In higher dimensions, not only the
first but also the first two principal components are not able to distinguish
two clusters, which can hide an important data structure when applying the
standard data visualization tools.

In the first test we study a switch of the first principal component from
following the variance of the dense informative distribution (abscissa) to the
sparse noise distribution (ordinate) as a function of the number of
contaminating points, in $R^{2}$ (Figure 5A-C). We modeled two clusters as two
100-point samples from normal distribution centered in points $[-1;0]$ and
$[1;0]$ with isotropic variance with the standard deviation 0.1. The sparse
noise distribution was modeled as a $k$-point sample from the product of two
Laplace distributions of zero means with the standard deviations 2 along
abscissa and 4 along ordinate. The intervals for computing the PQSQ functional
were defined by thresholds $R=\\{0;0.01;0.1;0.5;1\\}$ for each coordinate.
Increasing the number of points in the contaminating distribution diminishes
the average value of the abscissa coordinate of PC1, because the PC1 starts to
be attracted by the contaminating distribution (Figure 5C). However, it is
clear that on average PQSQ $L_{1}$-based PCA is able to withstand much larger
amplitude of the contaminating signal (very robust up to 20-30 points, i.e.
10-20% of strong noise contamination) compared to the standard $L_{2}$-based
PCA (which is robust to 2-3% of contamination).

In the second test we study the ability of the first two principal components
to separate two clusters, in $R^{100}$ (Figure 5D-F). As in the first test, we
modeled two clusters as two 100-point samples from normal distribution
centered in points $[-1;0;\dots;0]$ and $[1;0;\dots;0]$ with isotropic
variance with the standard deviation 0.1 in all 100 dimensions. The sparse
noise distribution is modeled as a $k$-point sample from the product of 100
Laplace distributions of zero means with the standard deviations 1 along each
coordinate besides the third coordinate (standard deviation of noise is 2) and
the forth coordinate (standard deviation of noise is 4). Therefore, the first
two principal component in the absence of noise are attracted by the
dimensions 1 and 2, while in the presence of strong noise they are be
attracted by dimensions 3 and 4, hiding the cluster structure of the dense
part of the distribution. The definitions of the intervals were taken as in
the first test. We measured the ability of the first two principal components
to separate clusters by computing the $t$-test between the two clusters
projected in the 2D-space spanned by the first principal components of the
global distribution (Figure 5D-F). As one can see, the average ability of the
first principal components to separate clusters is significantly stronger in
the case of PQSQ $L_{1}$-based PCA which is able to separate robustly the
clusters even in the presence of strong noise contamination (up to 80 noise
points, i.e. 40% contamination).

Figure 5: Comparing $L_{2}$\- and PQSQ $L_{1}$-based PCA using example of two-
cluster distribution (100 black circles and 100 squares) contaminated by
sparse noise (red crosses). A) Dense two cluster distribution contaminated by
sparse distribution (20 points) of large variance. In the presence of noise,
the abscissa coordinate $x$ of PC1 vector is slightly less than 1. B) Same as
A) but in the case of strong contamination (60 points). The value of $x$ is
much smaller in this case. C) Average absolute value of the abscissa
coordinate of PC1 $|x|$ (thick lines) shown with standard interval (thin
lines) for 100 samples of $k$ contaminating points. D) Projection of the data
distribution on the first two principal components computed with the standard
$L_{2}$ PCA algorithm. The number of contaminating points is 40. The cluster
structure of the dense part of the distribution is completely hidden as shown
in the zoom window. E) Same as in D) but computed with PQSQ $L_{1}$-based
algorithm. The cluster structure is perfectly separable. F) The value of
$t$-test computed based on the known cluster labels of the dense part of the
distribution, in the projections onto the first two principal components of
the global distribution. As in C), the mean values of 100 contamination
samples together with confidence intervals are shown.

###  6.4 Performance/stability trade-off benchmarking of $L_{1}$-based PCA

In order to compare the computation time and the robustness of PQSQ-based PCA
algorithm for the case $u(x)=|x|$ with existing R-based implementations of
$L_{1}$-based PCA methods (pcaL1 package), we followed the benchmark described
in [26]. We compared performance of PQSQ-based PCA based on Algorithm 3 with
several $L_{1}$-based PCA algorithms: L1-PCA* [25], L1-PCA [23], PCA-PP [34],
PCA-L1 [24]. As a reference point, we used the standard PCA algorithm based on
quadratic norm and computed using the standard SVD iterations.

The idea of benchmarking is to generate a series of datasets of the same size
($N=1000$ objects in $m=10$ dimensions) such that the first 5 dimensions would
be sampled from a uniform distribution $U(-10,10)$. Therefore, the first 5
dimensions represent ‘true manifold’ sampled by points.

The values in the last 5 dimensions represent ‘noise+outlier’ signal. The
background noise is represented by Laplacian distribution of zero mean and 0.1
variance. The outlier signal is characterized by mean value $\mu$, dimension
$p$ and frequency $\phi$. Then, for each data point with a probability $\phi$,
in the first $p$ outlier dimensions a value is drawn from $Laplace(\mu,0.1)$.
The rest of the values is drawn from background noise distribution.

As in [26], we’ve generated 1300 test sets corresponding to $\phi=0.1$, with
100 examples for each combination of $\mu=1,5,10,25$ and $p=1,2,3$.

For each test set 5 first principal components $\vec{V}_{1},\dots,\vec{V}_{5}$
of unit length were computed, with corresponding point projection
distributions $U^{1},\dots,U^{5}$ and the mean vector $\vec{C}$. Therefore,
for each initial data point $\vec{x}_{i}$, we have the ‘restored’ data point

$P(\vec{x}_{i})=\sum_{k=1,\dots,5}U^{k}_{i}\vec{V}_{k}+\vec{C}.$

For computing the PQSQ-based PCA we used 5 intervals without trimming.
Changing the number of intervals did not significantly changed the
benchmarking results.

Two characteristics were measured: (1) computation time measured as a ratio to
the computation of 5 first principal components using the standard
$L_{2}$-based PCA and (2) the sum of absolute values of the restored point
coordinates in the ‘outlier’ dimensions normalized on the number of points:

Figure 6: Benchmarking several algorithms for constructing $L_{1}$-based PCA,
using synthetic datasets representing ‘true’ five-dimensional linear manifold
contaminated by noise and outliers (located in other five dimensions). The
abscissa is the error of detecti tnghe ‘true’ manifold by a particular method
and the ordinate is the computational time relative to the standard SVD
($L_{2}$-based PCA) computation, in logarithmic scale. The computational cost
of application of linear programming methods instead of simpler iterative
methods is approximately shown by an arrow.

$\sigma=\frac{1}{N}\sum_{i=1,\dots,N}\sum_{k=6,\dots,10}|P^{k}(\vec{x}_{i})|.$
(22)

Formally speaking, $\sigma$ is $L_{1}$-based distance from the point
projection onto the first five principal components to the ‘true’ subspace. In
simplistic terms, larger values of $\sigma$ correspond to the situation when
the first five principal components do not represent well the first ‘true’
dimensions having significant loadings into the ‘outlier dimensions’.
$\sigma=0$ if and only if the first five components do not have any non-zero
loadings in the dimensions $6,\dots,10$.

The results of the comparison averaged over all 1300 test sets are shown in
Figure 6. The PQSQ-based computation of PCA outperforms by accuracy the
existing heuristics such as PCA-L1 but remains computationally efficient. It
is 100 times faster than L1-PCA giving almost the same accuracy. It is almost
500 times faster than the L1-PCA* algorithm, which is, however, better in
accuracy (due to being robust with respect to strong outliers). One can see
from Figure 6 that PQSQ-based approach is the best in accuracy among fast
iterative methods. The detailed tables of comparison for all combinations of
parameters are available on GitHub333http://goo.gl/sXBvqh. The scripts used to
generate the datasets and compare the results can also be found
there444https://github.com/auranic/PQSQ-DataApproximators.

###  6.5 Comparing performances of PQSQ-based regularized regression and lasso
algorithms

We compared performance of PQSQ-based regularized regression method imitating
$L_{1}$ penalty with lasso implementation in Matlab, using 8 datasets from UCI
Machine Learning Repository [35], Regression Task section. In the selection of
datasets we chose very different numbers of objects and variables for
regression construction (Table 1). All table rows containing missing values
were eliminated for the tests.

Table 1: Comparing time performance (in seconds, on ordinary laptop) of lasso vs PQSQ-based regularized regression imitating $L_{1}$ penalty. Average acceleration of PQSQ-based method vs lasso in these 8 examples is 120 fold with comparable accuracy.  Dataset | Objects | Variables | lasso | PQSQ | Ratio
---|---|---|---|---|---
Breast cancer | 47 | 31 | 10.50 | 0.05 | 233.33
Prostate cancer | 97 | 8 | 0.07 | 0.02 | 4.19
ENB2012 | 768 | 8 | 0.53 | 0.03 | 19.63
Parkinson | 5875 | 26 | 20.30 | 0.04 | 548.65
Crime | 1994 | 100 | 10.50 | 0.19 | 56.24
Crime reduced | 200 | 100 | 17.50 | 0.17 | 102.94
Forest fires | 517 | 8 | 0.05 | 0.02 | 3.06
Random regression (1000$\times$250) | 1000 | 250 | 2.82 | 0.58 | 4.86

We observed up to two orders of magnitude acceleration of PQSQ-based method
compared to the lasso method implemented in Matlab (Table 1), with similar
sparsity properties and approximation power as lasso (Figure 7).

While comparing time performances of two methods, we’ve noticed that lasso (as
it is implemented in Matlab) showed worse results when the number of objects
in the dataset approaches the number of predictive variables (see Table 1). In
order to test this observation explicitly, we took a particular dataset
(’Crime’) containing 1994 observations and 100 variables and compared the
performance of lasso in the case of complete table and a reduced table (’Crime
reduced’) containing only each 10th observation. Paradoxically, lasso
converges almost two times slower in the case of the smaller dataset, while
the PQSQ-based algorithm worked slightly faster in this case.

It is necessary to stress that here we compare the basic algorithms without
many latest technical improvements which can be applied both to $L_{1}$
penalty and its PQSQ approximation (such as fitting the whole lasso path).
Detailed comparison of all the existent modifications if far beyond the scope
of this work.

For comparing approximation power of the PQSQ-based regularized regression and
lasso, we used two versions of PQSQ potential for regression coefficients:
with and without trimming. In order to represent the results, we used the
‘Number of non-zero parameters vs Fraction of Variance Unexplained (FVU)’
plots (see two representative examples at Figure 7). We suggest that this type
of plot is more informative in practical applications than the ’lasso plot’
used to calibrate the strength of regularization, since it is a more explicit
representation for optimizing the accuracy vs complexity ratio of the
resulting regression.

From our testing, we can conclude that PQSQ-based regularized regression has
similar properties of sparsity and approximation accuracy compared to lasso.
It tends to slightly outperform lasso (to give smaller FVU) in case of
$N\approx P$. Introducing trimming in most cases does not change the best FVU
for a given number of selected variables, but tends to decrease its variance
(has a stabilization effect). In some cases, introducing trimming is the most
advantageous method (Figure 7B).

The GitHub ‘PQSQ-regularized-regression’ repository contains exact dataset
references and more complete report on comparing approximation ability of
PQSQ-based regularized regression with lasso555https://github.com/Mirkes/PQSQ-
regularized-regression/wiki.

Figure 7: Number of non-zero regression coefficients vs FVU plot for two
example real-life datasets (A - Breast cancer Wisconsin dataset from UC Irvine
Machine Learning Repository, B - original prostate cancer example from the
seminal lasso paper[1]). Each cross shows a particular solution of the
regularized regression problem. Solid lines show the best (minimal) obtained
FVU within the same number of selected variables.

##  7 Conclusion

In this paper we develop a new machine learning framework (theory and
application) allowing one to deal with arbitrary error potentials of not-
faster than quadratic growth, imitated by piece-wise quadratic function of
subquadratic growth (PQSQ error potential).

We develop methods for constructing the standard data approximators (mean
value, $k$-means clustering, principal components, principal graphs) for
arbitrary non-quadratic approximation error with subquadratic growth and
regularized linear regression with arbitrary subquadratic penalty by using a
piecewise-quadratic error functional (PQSQ potential). These problems can be
solved by applying quasi-quadratic optimization procedures, which are
organized as solutions of sequences of linear problems by standard and
computationally efficient algorithms.

Table 2: List of methods which can use PQSQ-based error potentials  Data
approximation/Clustering/Manifold learning
---
Principal Component Analysis |  Includes robust trimmed version of PCA, $L_{1}$-based PCA, regularized PCA, and many other PCA modifications
Principal curves and manifolds |  Provides possibility to use non-quadratic data approximation terms and trimmed robust version
Self-Organizing maps |  Same as above
Principal graphs/trees |  Same as above
K-means |  Can include adaptive error potentials based on estimated error distributions inside clusters
High-dimensional data mining
Use of fractional quasinorms $L_{p}$ ($0<p<1$) |  Introducing fractional quasinorms in existing data-mining techniques can potentially deal with the curse of dimensionality, helping to better distinguish close from distant data points [12]
Regularized and sparse regression
Lasso |  Application of PQSQ-based potentials leads to speeding up in case of large and $N\approx P$ datasets
Elastic net |  Same as above

The suggested methodology have several advantages over existing ones:

(a) Scalability: the algorithms are computationally efficient and can be
applied to large data sets containing millions of numerical values.

(b) Flexibility: the algorithms can be adapted to any type of data metrics
with subquadratic growth, even if the metrics can not be expressed in explicit
form. For example, the error potential can be chosen as adaptive metrics [36,
37].

(c) Built-in (trimmed) robustness: choice of intervals in PQSQ can be done in
the way to achieve a trimmed version of the standard data approximators, when
points distant from the approximator do not affect the error minimization
during the current optimization step.

(d) Guaranteed convergence: the suggested algorithms converge to local or
global minimum just as the corresponding predecessor algorithms based on
quadratic optimization and expectation/minimization-based splitting approach.

In theoretical perspective, using PQSQ-potentials in data mining is similar to
existing applications of min-plus (or, max-plus) algebras in non-linear
optimization theory, where complex non-linear functions are approximated by
infimum (or supremum) of finitely many ‘dictionary functions’ [38, 39]. We can
claim that just as using polynomes is a natural framework for approximating in
rings of functions, using min-plus algebra naturally leads to introduction of
PQSQ-based functions and the cones of minorants of quadratic dictionary
functions.

One of the application of the suggested methodology is approximating the
popular in data mining $L_{1}$ metrics. We show by numerical simulations that
PQSQ-based approximators perform as fast as the fast heuristical algorithms
for computing $L_{1}$-based PCA but achieve better accuracy in a previously
suggested benchmark test. PQSQ-based approximators can be less accurate than
the exact algorithms for optimizing $L_{1}$-based functions utilizing linear
programming: however, they are several orders of magnitude faster. PQSQ
potential can be applied in the task of regression, replacing the classical
Least Squares or $L_{1}$-based Least Absolute Deviation methods. At the same
time, PQSQ-based approximators can imitate a variety of subquadratic error
potentials (not limited to $L_{1}$ or variations), including fractional
quasinorms $L_{p}$ ($0<p<1$). We demonstrate that the PQSQ potential can be
easily adapted to the problems of sparse regularized regression with non-
quadratic penalty on regression coefficients (including imitations of lasso
and elastic net). On several real-life dataset examples we show that PQSQ-
based regularized regression can perform two orders of magnitude faster than
the lasso algorithm implemented in the same programming environment.

To conclude, in Table 2 we list possible applications of the PQSQ-based
framework in machine learning.

## Aknowledgement

This study was supported in part by Big Data Paris Science et Lettre Research
University project ‘PSL Institute for Data Science’.

## References

  * Tibshirani [1996] R. Tibshirani,  Regression shrinkage and selection via the lasso,  Journal of the Royal Statistical Society. Series B (Methodological) (1996) 267–288.
  * Zou and Hastie [2005] H. Zou, T. Hastie,  Regularization and variable selection via the elastic net,  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 67 (2005) 301–320.
  * Barillot et al. [2012] E. Barillot, L. Calzone, P. Hupe, J.-P. Vert, A. Zinovyev, Computational Systems Biology of Cancer, Chapman & Hall, CRC Mathemtical and Computational Biology, 2012.
  * Wright et al. [2010] J. Wright, Y. Ma, J. Mairal, G. Sapiro, T. S. Huang, S. Yan,  Sparse representation for computer vision and pattern recognition,  Proceedings of the IEEE 98 (2010) 1031–1044.
  * Lloyd [1957] S. Lloyd,  Last square quantization in pcm’s,  Bell Telephone Laboratories Paper (1957).
  * Pearson [1901] K. Pearson,  On lines and planes of closest fit to systems of points in space,  Philos. Mag. 2 (1901) 559–572.
  * Gorban and Zinovyev [2009] A. N. Gorban, A. Zinovyev,  Principal graphs and manifolds,  In Handbook of Research on Machine Learning Applications and Trends: Algorithms, Methods and Techniques, eds. Olivas E.S., Guererro J.D.M., Sober M.M., Benedito J.R.M., Lopes A.J.S. (2009).
  * Gorban et al. [2008] A. Gorban, B. Kegl, D. Wunsch, A. Zinovyev (Eds.), Principal Manifolds for Data Visualisation and Dimension Reduction, LNCSE 58, Springer, 2008.
  * Flury [1990] B. Flury,  Principal points,  Biometrika 77 (1990) 33–41.
  * Hastie [1984] T. Hastie,  Principal curves and surfaces,  PhD Thesis, Stanford University, California (1984).
  * Gorban et al. [2007] A. N. Gorban, N. R. Sumner, A. Y. Zinovyev,  Topological grammars for data approximation,  Applied Mathematics Letters 20 (2007) 382–386.
  * Aggarwal et al. [2001] C. C. Aggarwal, A. Hinneburg, D. A. Keim,  On the surprising behavior of distance metrics in high dimensional space,  in: Database Theory - ICDT 2001, 8th International Conference, London, UK, January 4-6, 2001, Proceedings, Springer, 2001, pp. 420–434.
  * Ding et al. [2006] C. Ding, D. Zhou, X. He, H. Zha,  R1-PCA: rotational invariant L1-norm principal component analysis for robust subspace factorization,  ICML (2006) 281–288.
  * Hauberg et al. [2014] S. Hauberg, A. Feragen, M. J. Black,  Grassmann averages for scalable robust pca,  in: Computer Vision and Pattern Recognition (CVPR), 2014 IEEE Conference on, IEEE, pp. 3810–3817.
  * Xu and Yuille [1995] L. Xu, A. L. Yuille,  Robust principal component analysis by self-organizing rules based on statistical physics approach,  Neural Networks, IEEE Transactions on 6 (1995) 131–143.
  * Allende et al. [2004] H. Allende, C. Rogel, S. Moreno, R. Salas,  Robust neural gas for the analysis of data with outliers,  in: Computer Science Society, 2004. SCCC 2004. 24th International Conference of the Chilean, IEEE, pp. 149–155.
  * Kohonen [2001] T. Kohonen, Self-organizing Maps, Springer Series in Information Sciences, Vol.30, Berlin, Springer, 2001.
  * Jolliffe et al. [2003] I. T. Jolliffe, N. T. Trendafilov, M. Uddin,  A modified principal component technique based on the lasso,  Journal of computational and Graphical Statistics 12 (2003) 531–547.
  * Candès et al. [2011] E. J. Candès, X. Li, Y. Ma, J. Wright,  Robust principal component analysis?,  Journal of the ACM (JACM) 58 (2011) 11.
  * Zou et al. [2006] H. Zou, T. Hastie, R. Tibshirani,  Sparse principal component analysis,  Journal of computational and graphical statistics 15 (2006) 265–286.
  * Cuesta-Albertos et al. [1997] J. Cuesta-Albertos, A. Gordaliza, C. Matrán, et al.,  Trimmed $k$-means: An attempt to robustify quantizers,  The Annals of Statistics 25 (1997) 553–576.
  * Lu et al. [2015] C. Lu, Z. Lin, S. Yan,  Smoothed low rank and sparse matrix recovery by iteratively reweighted least squares minimization.,  IEEE Trans Image Process 24 (2015) 646–654.
  * Ke and Kanade [2005] Q. Ke, 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, volume 1, IEEE, pp. 739–746.
  * Kwak [2008] N. Kwak,  Principal component analysis based on L1-norm maximization,  Pattern Analysis and Machine Intelligence, IEEE Transactions on 30 (2008) 1672–1680.
  * Brooks et al. [2013] J. Brooks, J. Dulá, E. Boone,  A pure L1-norm principal component analysis,  Comput Stat Data Anal 61 (2013) 83–98.
  * Brooks and Jot [2012] J. Brooks, S. Jot, PCAL1: An implementation in R of three methods for L1-norm principal component analysis, Optimization Online preprint, http://www.optimization-online.org/DB_HTML/2012/04/3436.html, 2012\.
  * Park and Klabjan [2014] Y. W. Park, D. Klabjan, Algorithms for L1-norm principal component analysis. Tutorial, http://dynresmanagement.com/uploads/3/3/2/9/3329212/algorithms_for_l1pca.pdf, 2014\.
  * Gorban’ and Rossiev [1999] A. Gorban’, A. Rossiev,  Neural network iterative method of principal curves for data with gaps,  Journal of Computer and Systems Sciences International c/c of Tekhnicheskaia Kibernetika 38 (1999) 825–830.
  * Gorban and Zinovyev [2001a] A. Gorban, A. Zinovyev,  Visualization of data by method of elastic maps and its applications in genomics, economics and sociology,  IHES Preprints (2001a).
  * Gorban and Zinovyev [2001b] A. Gorban, A. Y. Zinovyev,  Method of elastic maps and its applications in data visualization and data modeling,  International Journal of Computing Anticipatory Systems, CHAOS 12 (2001b) 353–369.
  * Gorban and Zinovyev [2005] A. Gorban, A. Zinovyev,  Elastic principal graphs and manifolds and their practical applications,  Computing 75 (2005) 359–379.
  * Gorban and Zinovyev [2010] A. N. Gorban, A. Zinovyev,  Principal manifolds and graphs in practice: from molecular biology to dynamical systems.,  Int J Neural Syst 20 (2010) 219–232.
  * Gorban et al. [2014] A. N. Gorban, A. Pitenko, A. Zinovyev,  ViDaExpert: user-friendly tool for nonlinear visualization and analysis of multidimensional vectorial data,  arXiv preprint arXiv:1406.5550 (2014).
  * Croux et al. [2007] C. Croux, P. Filzmoser, M. R. Oliveira,  Algorithms for projection–pursuit robust principal component analysis,  Chemometrics and Intelligent Laboratory Systems 87 (2007) 218–225.
  * Lichman [2013] M. Lichman, University of California, Irvine (UCI) Machine Learning Repository, http://archive.ics.uci.edu/ml, 2013\.
  * Yang and Jin [2006] L. Yang, R. Jin,  Distance metric learning: A comprehensive survey,  Michigan State Universiy 2 (2006).
  * Wu et al. [2009] L. Wu, R. Jin, S. C. Hoi, J. Zhu, N. Yu,  Learning Bregman distance functions and its application for semi-supervised clustering,  in: Advances in neural information processing systems, pp. 2089–2097.
  * Gaubert et al. [2011] S. Gaubert, W. McEneaney, Z. Qu,  Curse of dimensionality reduction in max-plus based approximation methods: theoretical estimates and improved pruning algorithms,  Arxiv preprint 1109.5241 (2011).
  * Magron et al. [2015] V. Magron, X. Allamigeon, S. Gaubert, B. Werner,  Formal proofs for nonlinear optimization,  Arxiv preprint 1404.7282 (2015).