# On Generalization and Computation of Tukey’s Depth: Part II

Yiyuan She, Shao Tang, and Jingze Liu
Department of Statistics, Florida State University

###### Abstract

This paper studies how to generalize Tukey’s depth to problems defined in a
restricted space that may be curved or have boundaries, and to problems with a
nondifferentiable objective. First, using a manifold approach, we propose a
broad class of Riemannian depth for smooth problems defined on a Riemannian
manifold, and showcase its applications in spherical data analysis, principal
component analysis, and multivariate orthogonal regression. Moreover, for
nonsmooth problems, we introduce additional slack variables and inequality
constraints to define a novel slacked data depth, which can perform center-
outward rankings of estimators arising from sparse learning and reduced rank
regression. Real data examples illustrate the usefulness of some proposed data
depths.

_Keywords_ : Riemannian depth, principal component analysis, slacked data
depth, reduced rank regression, sparsity-promoting regularizers.

##  1 Introduction

Tukey’s half-space depth (Tukey,, 1975) can be generalized to a polished
subspace depth, as shown in our companion paper (She et al., 2021a, ). The
basic Tukeyfication process there assumes a simple problem structure in the
sense that one can directly write down some sample-additive estimating
equations. Modern statistical applications however pose new challenges.

First, the parameter space $\Omega$ can be curved or have boundaries, so that
evaluating the gradient in the ambient Euclidean space may not directly
deliver reasonable influences. Second, with a regularizer in use, the
objective function is typically non-differentiable. Some specific examples are
given as follows.

######  Example 1.

(Watson depth) Assume that all data points lie on an $m$-dimensional sphere
$\boldsymbol{z}_{i}\in\mathbb{S}^{m-1}$, and $\pm\boldsymbol{z}_{i}$ are
deemed equivalent. This kind of data are typically referred to as axially
symmetric data. They have recently received attention in clustering and
directional statistics (Dhillon et al.,, 2003; Bijral et al.,, 2007; Sra and
Karp,, 2013). To characterize the distribution of such data, a commonly used
one is the Watson distribution with density (Watson,, 1965; Mardia and Jupp,,
1999)

$p(\boldsymbol{z};\boldsymbol{\mu},\kappa)\propto
e^{\kappa(\boldsymbol{\mu}^{T}\boldsymbol{z})^{2}}.$

Here, $\boldsymbol{\mu}\in\mathbb{S}^{m-1}$ gives the mean direction, $\kappa$
is the so-called concentration parameter, and the normalizing constant does
not depend on $\boldsymbol{\mu}$. We require $\kappa\neq 0$ (otherwise
$\boldsymbol{\mu}$ is not an effective parameter to introduce depth). When
$\kappa>0$, the data points concentrate around $\boldsymbol{\mu}$, and when
$\kappa<0$, the data spread around the great circle orthogonal to
$\boldsymbol{\mu}$. How to “Tukeyfy” more complex distributions defined on a
sphere (such as the Fisher-Bingham distribution) is nontrivial, but could give
rise to more useful spherical data depths. We will see that depth-enhanced
principal component analysis to be introduced in Section 2.2 poses a similar
manifold challenge.

######  Example 2.

(Nonnegative regression depth) As an extension of the celebrated regression
depth (Rousseeuw and Hubert,, 1999), let us consider data depth in a setting
where all coefficients are nonnegative. This corresponds to the nonnegative
least squares problem:

$\min_{\boldsymbol{\beta}}\lVert\boldsymbol{y}-\boldsymbol{X}\boldsymbol{\beta}\rVert_{2}^{2}\;\textrm{s.t.}\;\beta_{j}\geq
0,1\leq j\leq p,$

where we can denote the constraints by
$\boldsymbol{\beta}\in\mathbb{R}_{+}^{p}$ with $\mathbb{R}_{+}=[0,\infty)$.
Clearly, the closed parameter space $\mathbb{R}_{+}^{p}$ has boundary points.
Regression depth can be simply applied if $\boldsymbol{\beta}^{\circ}$ is an
interior point, but if $\boldsymbol{\beta}^{\circ}$ lies on the boundary,
i.e., $\beta_{j}^{\circ}=0$ for some $j$, which is of practical interest in
significance tests, regular depth does not apply, and the normal-equation
based influences must be corrected—but how?

######  Example 3.

(Sparsity depth) Consider a sparse learning problem

$\min_{\boldsymbol{\beta}}\sum_{i=1}^{n}\allowbreak
l_{0}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta})+\sum_{j=1}^{p}P(|\beta_{j}|;\lambda),$

where $l_{0}$ is a loss function defined on the systematic component
$\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}$ and $P$ is a penalty function to
promote sparsity in $\boldsymbol{\beta}$. Examples of $P$ include $\ell_{1}$,
$\ell_{0}$, SCAD (Fan and Li,, 2001), and MCP (Zhang,, 2010), among many
others that are popularly used in high dimensional statistics for building a
parsimonious model. We assume that the regularization parameter $\lambda$ is
given, either by theory—see, e.g., Cai et al., (2009), or by tuning, like
cross-validation (She and Tran,, 2019), so that the criterion is fully
specified. A new class of depths like $\ell_{1}$-depth or $\ell_{0}$-depth
would be helpful for high-dimensional robust inference, but the
nondifferentiability and nonconvexity of $P$ make it difficult to obtain
sample-additive estimation equations.

To tackle the challenges, we propose two approaches based on manifolds and
slack variables, respectively, to extend Tukey’s depth to Riemannian depth and
slacked depth. The rest of the paper is organized as follows. Section 2
studies how to handle a smooth problem defined on a Riemannian manifold. The
resulting Riemannian depth finds applications in spherical data analysis,
principal component analysis, and multivariate orthogonal regression. Section
3 uses slack variables to cope with parameter spaces with boundaries and
nondifferentiable objectives. A novel class of slacked data depth can perform
center-outward rankings of estimators arising from sparse learning and reduced
rank regression. Section 4 performs computer experiments on some real data
examples. We conclude the whole work in Section 5.

#### Notation.

We use bold symbols to denote vectors and matrices. A matrix
$\boldsymbol{X}\in\mathbb{R}^{n\times p}$ is frequently partitioned into rows
$\boldsymbol{X}=[\boldsymbol{x}_{1}\ldots\boldsymbol{x}_{n}]^{T}$ with
$\boldsymbol{x}_{i}\in\mathbb{R}^{p}$. The vectorization of $\boldsymbol{X}$
is denoted by $\textrm{vec}(\boldsymbol{X})\in\mathbb{R}^{np}$. Let
$\mathbb{R}_{+}=[0,+\infty]$. We use $\boldsymbol{X}[\mathcal{I},\mathcal{J}]$
to denote a submatrix of $\boldsymbol{X}$ with rows and columns indexed by
$\mathcal{I}$ and $\mathcal{J}$, respectively, and occasionally abbreviate
$\boldsymbol{X}[,\mathcal{J}]$ to $\boldsymbol{X}_{\mathcal{J}}$ by selecting
the corresponding columns. Given $\boldsymbol{X}\in\mathbb{R}^{n\times p}$,
$\lVert\boldsymbol{X}\rVert_{F}$ and $\lVert\boldsymbol{X}\rVert_{2}$ denote
its Frobenius norm and spectral norm, respectively,
$\|\boldsymbol{X}\|_{\max}\triangleq\max_{1\leq i\leq n,1\leq j\leq
p}|x_{ij}|$, and $\mbox{rank}(\boldsymbol{X})$ denotes its rank. The Moore-
Penrose inverse of $\boldsymbol{X}$ is denoted by $\boldsymbol{X}^{+}$. The
inner product of two matrices $\boldsymbol{X}$ and $\boldsymbol{Y}$ (of the
same size) is defined as
$\langle\boldsymbol{X},\boldsymbol{Y}\rangle=\textrm{Tr}(\boldsymbol{X}^{T}\boldsymbol{Y})$
and their element-wise product (Hadamard product) is
$\boldsymbol{X}\circ\boldsymbol{Y}$. The Kronecker product is denoted by
$\boldsymbol{X}\otimes\boldsymbol{Y}$ (where $\boldsymbol{X}$ and
$\boldsymbol{Y}$ need not have the same dimensions). Given a set
$\mathcal{A}\subset\mathbb{R}^{p\times m}$ and a matrix
$\boldsymbol{T}\in\mathbb{R}^{n\times p}$,
$\boldsymbol{T}\circ\mathcal{A}=\\{\boldsymbol{T}\boldsymbol{A}:\boldsymbol{A}\in\mathcal{A}\\}$.
We use $\mathbb{O}^{m\times r}$ to represent the set of all $m\times r$
matrices $\boldsymbol{V}$ satisfying the orthogonality constraint
$\boldsymbol{V}^{T}\boldsymbol{V}=\boldsymbol{I}$. For a vector
$\boldsymbol{a}=[a_{1},\ldots,a_{n}]^{T}$ $\in\mathbb{R}^{n}$,
$\mbox{diag}\\{\boldsymbol{a}\\}$ is defined as an $n\times n$ diagonal matrix
with diagonal entries given by $a_{1},\ldots,a_{n}$, and for a square matrix
$\boldsymbol{A}=[a_{ij}]_{n\times n}$,
$\mbox{diag}(\boldsymbol{A}):=\mbox{diag}\\{a_{11},\ldots,a_{nn}\\}$. The
indicator function $1_{\mathcal{A}}(t)$ means $1_{\mathcal{A}}(t)=1$ if
$t\in\mathcal{A}$ and $0$ otherwise. Given $f:\mathbb{R}^{n\times
p}\rightarrow\mathbb{R}$, $f\in\mathcal{C}^{1}$ means that its Euclidean
gradient $\nabla f(\boldsymbol{X})$, an $n\times p$ matrix with the $(i,j)$
element ${\partial f}/{\partial x_{ij}}$, exists and is continuous for any
$\boldsymbol{X}\in\mathbb{R}^{n\times p}$. Given two vectors
$\boldsymbol{\alpha},\boldsymbol{\beta}\in\mathbb{R}^{p}$,
$\boldsymbol{\alpha}\succeq\boldsymbol{\beta}$ means
$\alpha_{j}\geq\beta_{j},1\leq j\leq p$ and
$\boldsymbol{\alpha}\succ\boldsymbol{\beta}$ means $\alpha_{j}>\beta_{j},1\leq
j\leq p$. Finally, $a\wedge b=\min\\{a,b\\}$.

##  2 A Manifold Approach

When the problem is defined on a Riemannian manifold (without boundaries), we
can introduce Riemannian influences, along with defining a proper influence
space to complete the definition of Riemannian depth. In contrast, the
commonly used methods to deal with constraints (such as the elimination
approach for (13) below) may be infeasible in higher dimensions. We will see
the important role of the influence space $\mathcal{G}$ introduced in Section
LABEL:sec:Tukeyfication of She et al., 2021a , since a Riemannian gradient
always lies in a tangent space.

###  2.1 Riemannian depth

We begin with Example 1 to motivate the main idea. Starting from such an
example, the Watson depth will be introduced. Useful for the analysis of axial
data, it is defined on a Riemmanian manifold and it will be a special case of
a more generic Riemannian depth.

For the MLE problem

$\min_{\boldsymbol{\mu}}-\kappa\sum_{i}\langle\boldsymbol{\mu},\boldsymbol{z}_{i}\rangle^{2}+c(\kappa;m)\;\textrm{s.t.}\;\lVert\boldsymbol{\mu}\rVert_{2}^{2}=1,$
(1)

Lagrange multiplier or eigenvalue decomposition can be used to solve for
$\boldsymbol{\mu}$, but they do not yield a simple set of estimating equations
like (LABEL:eq:sample_addtivity) in She et al., 2021a  to be conveniently used
for the purpose of data depth.

Instead, we view (1) as an unconstrained problem on the sphere
$\mathbb{S}^{m-1}$, which is a Riemannian manifold. Then the Riemannian
gradient with respect to $\boldsymbol{\mu}$ can be calculated to define the
desired (Riemannian) influence function. We show some detailed derivation to
give the reader more intuition.

Concretely, adopting the canonical metric induced by the inner product
$G_{\boldsymbol{\mu}}(\boldsymbol{u}_{1},\boldsymbol{u}_{2})=\boldsymbol{u}_{1}^{T}(\boldsymbol{I}-\boldsymbol{\mu}\boldsymbol{\mu}^{T}/2)\boldsymbol{u}_{2}$
(Edelman et al.,, 1998), the Riemannian gradient of
$l_{i}:=-\kappa\langle\boldsymbol{\mu},\boldsymbol{z}_{i}\rangle^{2}$ with
respect to $\boldsymbol{\mu}$, denoted by
$\boldsymbol{g}_{i}(\boldsymbol{\mu})$, is defined as the unique element in
the tangent space

$\mathcal{T}_{\boldsymbol{\mu}}(\mathbb{S}^{m-1})\triangleq\\{\boldsymbol{u}\in\mathbb{R}^{m}:\boldsymbol{u}^{T}\boldsymbol{\mu}=0\\}$

satisfying

$\displaystyle
G_{\boldsymbol{\mu}}(\boldsymbol{g}_{i}(\boldsymbol{\mu}),\boldsymbol{u})=\boldsymbol{u}^{T}\nabla
l_{i},\
\forall\boldsymbol{u}\in\mathcal{T}_{\boldsymbol{\mu}}(\mathbb{S}^{m-1})$ (2)

where $\nabla l_{i}$ is the Euclidean gradient. It follows that

$\boldsymbol{g}_{i}(\boldsymbol{\mu})=[\nabla
l_{i}\boldsymbol{\mu}^{T}-\boldsymbol{\mu}(\nabla
l_{i})^{T}]\boldsymbol{\mu}=-2\kappa\langle\boldsymbol{z}_{i},\boldsymbol{\mu}\rangle(\boldsymbol{z}_{i}-\langle\boldsymbol{z}_{i},\boldsymbol{\mu}\rangle\boldsymbol{\mu}).$

From Boothby, (1986) (and $\kappa\neq 0$), the optimal ${\boldsymbol{\mu}}$
satisfies

$\sum_{i}\langle\boldsymbol{z}_{i},\boldsymbol{\mu}\rangle(\boldsymbol{z}_{i}-\langle\boldsymbol{z}_{i},\boldsymbol{\mu}\rangle\boldsymbol{\mu})=\boldsymbol{0}.$
(3)

Given ${\boldsymbol{\mu}}^{\circ}\in\mathbb{S}^{m-1}$, the Riemannian
influence
$\langle\boldsymbol{z}_{i},\boldsymbol{\mu}^{\circ}\rangle(\boldsymbol{z}_{i}-\langle\boldsymbol{z}_{i},\boldsymbol{\mu}^{\circ}\rangle\boldsymbol{\mu}^{\circ})$,
denoted by $\boldsymbol{T}^{\mbox{\tiny
R}}(\boldsymbol{\mu}^{\circ};\boldsymbol{z}_{i})$, is no longer
$\boldsymbol{z}_{i}-\boldsymbol{\mu}^{\circ}$ as in location depth. Notably,
$\boldsymbol{T}^{\mbox{\tiny R}}(\boldsymbol{\mu}^{\circ};\boldsymbol{z}_{i})$
vanishes when $\theta_{i}=j\pi/2$ $(j=0,1,2,3)$ with
$\cos\theta_{i}=\langle\boldsymbol{z}_{i},{\boldsymbol{\mu}}^{\circ}\rangle$,
corresponding to various circumstances with $\kappa>0$ and $\kappa<0$.

Not only does the manifold perspective provide the desirable estimation
equations, but it defines an important influence space
$\mathcal{G}=\mathcal{T}_{\boldsymbol{\mu}^{\circ}}(\mathbb{S}^{m-1})$ to
restrict $\boldsymbol{v}$. Accordingly, our Watson depth considers all one-
dimensional projections tangentially passing through
$\boldsymbol{\mu}^{\circ}$:

$d_{01}^{\mbox{\tiny
W}}(\boldsymbol{\mu}^{\circ})=\min_{\boldsymbol{v}}\sum_{i}1_{\geq
0}(G_{\boldsymbol{\mu}^{\circ}}(\boldsymbol{v},(\boldsymbol{z}_{i}^{T}\boldsymbol{\mu}^{\circ})[\boldsymbol{z}_{i}-(\boldsymbol{z}_{i}^{T}\boldsymbol{\mu}^{\circ})\boldsymbol{\mu}^{\circ}]))\;\textrm{s.t.}\;\boldsymbol{v}^{T}\boldsymbol{\mu}^{\circ}=0,\boldsymbol{v}^{T}\boldsymbol{v}=1$
(4)

or equivalently

$d_{01}^{\mbox{\tiny
W}}(\boldsymbol{\mu}^{\circ})=\min_{\boldsymbol{v}}\sum_{i}1_{\geq
0}(\langle\boldsymbol{v},(\boldsymbol{z}_{i}^{T}\boldsymbol{\mu}^{\circ})\boldsymbol{z}_{i}\rangle)\;\textrm{s.t.}\;\boldsymbol{v}^{T}\boldsymbol{\mu}^{\circ}=0,\boldsymbol{v}^{T}\boldsymbol{v}=1,$
(5)

regardless of the Riemannian metric, as an outcome of (2). The factor
$\boldsymbol{z}_{i}^{T}\boldsymbol{\mu}^{\circ}$ in (5), possibly negative,
amounts to replacing $\boldsymbol{z}_{i}$ by
$\operatorname{sgn}\langle\boldsymbol{z}_{i}^{T},\boldsymbol{\mu}^{\circ}\rangle\cdot\boldsymbol{z}_{i}$.
This is in accordance with the Watson distribution for axially symmetric
spherical data. The algorithms in Section LABEL:sec:Data-depth-computation of
our companion paper can be applied, after a simple reparametrization of
$\boldsymbol{v}$ in the orthogonal complement space of
${\boldsymbol{\mu}}^{\circ}{\boldsymbol{\mu}}^{\circ T}$.

The above derivation is standard and can be generalized to introduce a
Riemannian depth for the Tukeyfication of a differentiable loss $l$ on a
Riemannian manifold $\mathcal{M}$ of an Euclidean space:
$\min_{\boldsymbol{B}}\sum_{i}l(\boldsymbol{B};\boldsymbol{x}_{i},\boldsymbol{y}_{i})$
s.t. $\boldsymbol{B}\in\mathcal{M}$. Given a point
$\boldsymbol{B}^{\circ}\in\mathcal{M}$ of interest, letting
$\boldsymbol{T}_{i}^{\circ}=\boldsymbol{T}(\boldsymbol{B}^{\circ};\boldsymbol{x}_{i},\boldsymbol{y}_{i})=\nabla_{\boldsymbol{B}}l(\boldsymbol{B}^{\circ};\boldsymbol{x}_{i},\boldsymbol{y}_{i})$
as before and considering all directional derivatives of $l$ in the directions
of $\boldsymbol{V}\in\mathcal{T}_{\boldsymbol{B}^{\circ}}(\mathcal{M})$, we
define

$\displaystyle\mbox{{Riemannian depth:}}\ \ d_{01}^{\mbox{\tiny
R}}(\boldsymbol{B}^{\circ})=\ $ $\displaystyle\
\min_{\boldsymbol{V}}\sum_{i}1_{\geq
0}(\langle\boldsymbol{V},\boldsymbol{T}_{i}^{\circ}\rangle)$ $\displaystyle\
\textrm{ s.t.
}\boldsymbol{V}\in\mathcal{T}_{\boldsymbol{B}^{\circ}}(\mathcal{M}),\lVert\boldsymbol{V}\rVert_{F}=1.$
(6)

Eqn. (6) performs location depth of Riemannian influences in the tangent space
$\mathcal{T}_{\boldsymbol{B}^{\circ}}(\mathcal{M})$. Because
$\mathcal{T}_{\boldsymbol{B}^{\circ}}(\mathcal{M})$ is linear, the restricted
Procrustes rotation in Section LABEL:sec:Data-depth-computation of She et al.,
2021a  applies with no difficulty in optimization.

When $\mathcal{M}$ is compact and/or $l$ is nonconvex, it becomes necessary to
exclude locally maximal solutions in the estimating equations. We give an
“order-2 Tukeyfication” as follows. Given
$\boldsymbol{B}^{\circ}\in\mathcal{M}$ and
$\boldsymbol{V}\in\mathcal{T}_{\boldsymbol{B}^{\circ}}(\mathcal{M})$, let
$\gamma$ be the geodesic satisfying $\gamma(0)=\boldsymbol{B}^{\circ}$ and
$\gamma^{\prime}(0)=\boldsymbol{V}$. The first step is to restrict $l$ to the
geodesic and define

$g_{i}=\frac{\,\mathrm{d}}{\,\mathrm{d}t}l(\gamma(t);\boldsymbol{x}_{i},\boldsymbol{y}_{i})\big{\rvert}_{t=0,},\quad
h_{i}=\frac{\,\mathrm{d}^{2}}{\,\mathrm{d}t^{2}}l(\gamma(t);\boldsymbol{x}_{i},\boldsymbol{y}_{i})\big{\rvert}_{t=0,}$

where $g_{i}$ simplifies to
$\langle\boldsymbol{V},\boldsymbol{T}_{i}^{\circ}\rangle$ and $h_{i}$ can be
calculated via Riemannian Hessian. (Our companion paper mostly considers an
Euclidean $\mathcal{M}$, where a line restriction
$l(\boldsymbol{B}^{\circ}+t\boldsymbol{V})$ with
$\boldsymbol{V}\neq\boldsymbol{0}$ is used, and $g_{i}$ and $h_{i}$ only
involve the ordinary gradient and Hessian of $l$.) The second step robustly
measures how well the following two optimality condition are obeyed:

$\sum_{i=1}^{n}g_{i}=0,\quad\sum_{i=1}^{n}h_{i}\geq 0.$

Concretely, changing the one-dimensional averages to medians motivates us to
adopt $\sum\left(1_{=0}+2(1_{<0}\wedge 1_{>0})\right)(g_{i})$ and $\sum
1_{\geq 0}(h_{i})$ to quantify to what extent the two condition are satisfied,
respectively, in the possible occurrence of extreme outliers. Finally,
combining the two measures leads to

$\displaystyle\mbox{{Riemannian depth (order 2):}}\ d_{01}^{\mbox{\tiny
R2}}(\boldsymbol{B}^{\circ})=\ $ $\displaystyle\
\min_{\boldsymbol{V}}\sum_{i}1_{\gtrsim 0}(g_{i})\sum_{i}1_{\geq 0}(h_{i})$
$\displaystyle\textrm{s.t.
}\boldsymbol{V}\in\mathcal{T}_{\boldsymbol{B}^{\circ}}(\mathcal{M}),\lVert\boldsymbol{V}\rVert_{F}=1,$

where $1_{\gtrsim 0}:=0.5\cdot 1_{=0}+1_{>0}$ replaces $1_{=0}+2(1_{<0}\wedge
1_{>0})$ in the optimization because $g_{i}$ is linear in $\boldsymbol{V}$,
the Riemannian Hessian is a bilinear map, and
$T_{\boldsymbol{B}^{\circ}}(\mathcal{M})$ is a linear space. A more aggressive
proposal is to use $\sum_{i}1_{\gtrsim 0}(g_{i})1_{\geq 0}(h_{i})$ as the
criterion. (Notice the mild difference between $1_{\gtrsim 0}$ and $1_{\geq
0}$; the first seems to be more appropriate to deal with equality-type
optimality conditions in defining a $d_{01}$-type data depth.) When $l$ is
(geodesically) convex, $h_{i}\geq 0$ and thus the associated factor with
proper scaling, $\sum_{i}1_{\geq 0}(h_{i})/n$, will not affect the depth.

######  Remark 1.

If we Tukeyfy the basic _von Mises-Fisher distribution_ (vMF) (Mardia and
Jupp,, 1999), with the density given by
$p(\boldsymbol{z};\boldsymbol{\mu},\kappa)\propto
e^{\kappa\boldsymbol{\mu}^{T}\boldsymbol{z}}$, where
${\boldsymbol{\mu}}:\|{\boldsymbol{\mu}}\|_{2}=1$ is the mean direction,
$\kappa>0$ and the normalizing constant does not depend on
${\boldsymbol{\mu}}$, (6) yields $d_{01}^{\mbox{\tiny
R}}(\boldsymbol{\mu}^{\circ})=\min_{\boldsymbol{v}}\sum_{i}1_{\geq
0}(\langle\boldsymbol{v},\boldsymbol{z}_{i}\rangle)\;\textrm{s.t.}\;\boldsymbol{v}^{T}\boldsymbol{\mu}^{\circ}=0,\boldsymbol{v}^{T}\boldsymbol{v}=1.$
This is closely related to but different from the angular Tukey’s depth that
can be defined as $\min_{\boldsymbol{v}}\sum_{i}1_{\geq
0}(\langle\boldsymbol{v},\boldsymbol{z}_{i}\rangle)\;\textrm{s.t.}\;\boldsymbol{v}^{T}\boldsymbol{\mu}^{\circ}\geq
0,\boldsymbol{v}^{T}\boldsymbol{v}=1$ for $m$-dimensional spherical data (see
Liu and Singh, (1992) for some theoretical studies when $m=2,3$). The order-2
depths involve
$h_{i}=\langle\boldsymbol{\mu}^{\circ},\boldsymbol{z}_{i}\rangle$ which are
independent of $\boldsymbol{v}$ in this case.

More interesting notions of spherical data depth can be induced by some more
flexible distributions through our manifold framework, such as the _Kent
distribution_ and the more general _Fisher-Bingham distribution_ whose
quadratic exponential form is more powerful than vMF for statistical modeling
in bioinformatics, meteorology, and computer vision.

###  2.2 Depth-enhanced principal component analysis

This part uses the Riemannian depth introduced in the last subsection to
Tukeyfy the well-known principal component analysis (PCA). Let
$\boldsymbol{Z}=[\boldsymbol{z}_{i},\ldots,\allowbreak\boldsymbol{z}_{n}]^{T}\in\mathbb{R}^{n\times
m}$ be a data matrix. The PCA model can be stated as

$\displaystyle\boldsymbol{Z}=\boldsymbol{1}{\boldsymbol{\mu}}^{*T}+\boldsymbol{A}^{*}\boldsymbol{U}^{*T}+\boldsymbol{E},$
(7)

with ${\boldsymbol{\mu}}^{*}\in\mathbb{R}^{m}$,
$\boldsymbol{A}^{*}\in\mathbb{R}^{n\times r}$, and
$\boldsymbol{U}^{*}\in\mathbb{O}^{n\times r}$ all unknown. Eqn. (7) means that
the $n$ data points, after some proper translation, all approximately
concentrate in an $r$-dimensional subspace, and $r$ is typically much lower
than $m$ and $n$. The columns of $\boldsymbol{U}^{*}$ are often called the
principal component (PC) loading directions. Assuming that the entries of
$\boldsymbol{E}$ are i.i.d. Gaussian, we can estimate the intercept vector and
the low-dimensional subspace by

$\displaystyle\min_{(\boldsymbol{U},\,\boldsymbol{\mu})}\lVert(\boldsymbol{Z}-\boldsymbol{1}\boldsymbol{\mu}^{T})(\boldsymbol{I}-\boldsymbol{U}\boldsymbol{U}^{T})\rVert_{F}^{2}\quad\textrm{s.t.}\quad\boldsymbol{U}^{T}\boldsymbol{U}=\boldsymbol{I}_{r\times
r}.$ (8)

The solution is given by standard PCA, which is however sensitive to outliers.

One may want to estimate $\boldsymbol{U},{\boldsymbol{\mu}}$ more robustly
through a depth enhancement. Here, the orthogonality constraint
$\boldsymbol{U}^{T}\boldsymbol{U}=\boldsymbol{I}$ may appear more complex than
that in the spherical problem (1), but (8) is a smooth problem on a Stiefel
manifold. Therefore, we can define a Riemannian depth for any
$({\boldsymbol{\mu}}^{\circ},\boldsymbol{U}^{\circ})\in\mathbb{R}^{m}\times\mathbb{O}^{m\times
r}$ based on Section 2.1, which we call the principal component (PC) depth, as
follows

$\displaystyle\mbox{{PC-depth}:}\ \ \min_{(\boldsymbol{v},\boldsymbol{V})}$
$\displaystyle\sum_{i}1_{\geq
0}(\langle\boldsymbol{v},(\boldsymbol{I}-\boldsymbol{U}^{\circ}\boldsymbol{U}^{\circ
T})({\boldsymbol{\mu}}^{\circ}-\boldsymbol{z}_{i})\rangle-\langle\boldsymbol{V},({\boldsymbol{\mu}}^{\circ}-\boldsymbol{z}_{i})({\boldsymbol{\mu}}^{\circ}-\boldsymbol{z}_{i})^{T}\boldsymbol{U}^{\circ}\rangle)$
$\displaystyle\textrm{ s.t.
}\boldsymbol{V}^{T}\boldsymbol{U}^{\circ}+\boldsymbol{U}^{\circ
T}\boldsymbol{V}=\boldsymbol{0},\
\|\boldsymbol{v}\|_{2}^{2}+\|\boldsymbol{V}\|_{F}^{2}=1.$ (9)

All matrix differentiation details are omitted. (9) may need an order-2
modification though, which will be clearly revealed by comparing it to (14)
later.

PCA is also helpful when ranking observations in ultra-high dimensions. It is
well known that the curse of dimensionality may make every observation look
like a corner point, thus harmful to describing data depth. Fortunately, under
(7), the true signals concentrate in the PC subspace determined by
$\boldsymbol{U}^{*}$; so to check a given point’s centrality or extremity, it
is helpful to project it onto the orthogonal complement (OC) subspace to
reveal its outlyingness. See She et al., (2016) for more discussions.
Specifically, letting $\bar{\boldsymbol{U}}^{*}\in\mathbb{O}^{m\times\bar{r}}$
($\bar{r}\leq m-r$) that is orthogonal to $\boldsymbol{U}^{*}$, we can obtain
from (7)

$\displaystyle\boldsymbol{Z}\bar{\boldsymbol{U}}^{*}=\boldsymbol{1}\bar{\boldsymbol{\mu}}^{*T}+\bar{\boldsymbol{E}},$
(10)

where
$\bar{\boldsymbol{\mu}}^{*}=\bar{\boldsymbol{U}}^{*T}{\boldsymbol{\mu}}^{*}$,
$\bar{\boldsymbol{E}}=\boldsymbol{E}\bar{\boldsymbol{U}}^{*}$. Eqn. (10) is in
the typical location estimation setting except that $\bar{\boldsymbol{U}}^{*}$
is unknown, which motivates us to consider

$\displaystyle\min_{(\bar{\boldsymbol{U}},\,\bar{\boldsymbol{\mu}})}\lVert\boldsymbol{Z}\bar{\boldsymbol{U}}-\boldsymbol{1}\bar{\boldsymbol{\mu}}^{T}\rVert_{F}^{2}\quad\textrm{s.t.}\quad\bar{\boldsymbol{U}}^{T}\bar{\boldsymbol{U}}=\boldsymbol{I}_{\bar{r}\times\bar{r}}.$
(11)

Interestingly, (11) can also be viewed as a _multivariate_ extension, of rank
$\bar{r}$, of the orthogonal regression due to Mizera, (2002):

$\displaystyle\min_{\mu\in\mathbb{R},\boldsymbol{u}\in\mathbb{R}^{p+1}}\left\|\,[\boldsymbol{X}\;\boldsymbol{y}]\,\boldsymbol{u}-\boldsymbol{1}\,\mu\,\right\|_{2}^{2}\mbox{
s.t. }\|\boldsymbol{u}\|_{2}^{2}=1.$ (12)

Moreover, when $\boldsymbol{Z}=[\boldsymbol{X}\;\boldsymbol{Y}]$, setting
$\bar{\boldsymbol{U}}=[\boldsymbol{B}^{T}\;\boldsymbol{\Gamma}^{T}]^{T}$ gives
a model
$\boldsymbol{Y}\boldsymbol{\Gamma}+\boldsymbol{X}\boldsymbol{B}-\boldsymbol{1}\bar{\boldsymbol{\mu}}^{T}=\boldsymbol{E}$
for canonical correlation analysis.

How to introduce an operational depth for (11) is a meaningful problem.
Indeed, with a deep $\bar{\boldsymbol{U}}$ provided, one would be able to rank
high-dimensional samples in a lower dimensional subspace.

Restricting to a naive case of (12) with a single predictor and a single
response:

$\displaystyle\min_{(\alpha,\beta,\mu)}\lVert\beta\boldsymbol{x}+\alpha\boldsymbol{y}-\boldsymbol{1}{\mu}\rVert_{F}^{2}\;\textrm{s.t.}\;\alpha^{2}+\beta^{2}=1,$
(13)

one can eliminate the constraint by, say, $\alpha=-\sin t,\beta=\cos t$ with a
free parameter $t$, and then take the Euclidean gradient with respect to
$(t,\mu)$ to define a tangent depth (Mizera,, 2002). Nevertheless, the
elimination method encounters difficulties when considering multiple
predictors, let alone a general $\bar{r}$. As far as we know, there exists no
commonly acknowledged multivariate orthogonal regression depth in the
literature.

Our manifold approach provides a systematic treatment of (11) for all $p$,
$m$, and $\bar{r}$. We call the resulting Riemannian depth the orthogonal
complement (OC) depth. It pursues an $\bar{r}$-dimensional subspace in the
original input space to rank the observations effectively. The influence space
here is
$\mathbb{R}^{\bar{r}}\times\mathcal{T}_{\bar{\boldsymbol{U}}}(\mathbb{O}^{m\times\bar{r}})\
\mbox{ with }\
\mathcal{T}_{\bar{\boldsymbol{U}}}(\mathbb{O}^{m\times\bar{r}})=\\{\boldsymbol{V}:\bar{\boldsymbol{U}}^{T}\boldsymbol{V}+\boldsymbol{V}^{T}\bar{\boldsymbol{U}}=\boldsymbol{0}\\},$
and the OC depth for any given
$(\bar{\boldsymbol{\mu}}^{\circ},\bar{\boldsymbol{U}}^{\circ})\in\mathbb{R}^{\bar{r}}\times\mathbb{O}^{m\times\bar{r}}$
is

$\displaystyle\begin{split}\mbox{{OC-depth}:}\ \
\min_{(\boldsymbol{v},\boldsymbol{V})}&\sum_{i}1_{\geq
0}(\langle\boldsymbol{v},\bar{\boldsymbol{\mu}}^{\circ}-\bar{\boldsymbol{U}}^{\circ
T}\boldsymbol{z}_{i}\rangle+\langle\boldsymbol{V},\boldsymbol{z}_{i}\boldsymbol{z}_{i}^{T}\bar{\boldsymbol{U}}^{\circ}-\boldsymbol{z}_{i}\bar{\boldsymbol{\mu}}^{\circ
T}\rangle)\\\ &\textrm{ s.t.
}{\boldsymbol{V}}^{T}\bar{\boldsymbol{U}}^{\circ}+\bar{\boldsymbol{U}}^{\circ
T}\boldsymbol{V}=\boldsymbol{0},\
\|\boldsymbol{v}\|_{2}^{2}+\|\boldsymbol{V}\|_{F}^{2}=1.\end{split}$ (14)

The derivations are similar to the PC depth and are omitted. Note that the
influence space constraint has a multivariate form but is linear. Eqn. (14)
also gives a multivariate orthogonal regression depth.

On the other hand, with $r=\bar{r}$ and no intercepts, the (order-1) PC-depth
and OC-depth coincide, since in this case the two losses in (8) and (11) only
differ by a minus sign and the influence space is linear in $\boldsymbol{V}$.
This means that the most and least informative subspaces will have the same
depth. As aforementioned, an order-2 Riemannian depth would be able to
distinguish between minimization and maximization problems, which deserves
further investigation.

##  3 The Slack Variable Approach

The nondifferentiability issue in the other two examples in Section 1 is much
trickier to cope with. In more detail, Example 2 has a closed parameter space
$\mathbb{R}_{+}^{p}$ with boundaries, which makes gradient-based influences
improper at any boundary point; Example 3 has a nonsmooth regularizer commonly
seen in high dimensional statistics, and sometimes regularization can be
imposed in a constrained manner.

Following Rousseeuw and Hubert, (1999), the first step to define a data depth
is to characterize a reasonable “fit”, or a class of reasonable estimators,
under a given model or method. It turns out that for such nonsmooth problems,
we can derive local optimality conditions in form of inequalities or obtain
some nonlinear fixed-point equations by use of a surrogate function, neither
of which however results in sample-additive estimating equations directly. The
good news is that we can then utilize some “slack variables” subject to proper
(convex) inequality and equality constraints to offer a universal solution,
which leads to a novel class of slacked data depth.

###  3.1 Slacked data depth and sparse learning

To begin with, let us consider $\min
f(\boldsymbol{\beta})\triangleq\sum_{i}l(\boldsymbol{\beta};\boldsymbol{x}_{i},y_{i})\mbox{
s.t. }\boldsymbol{\beta}\succeq\boldsymbol{0}$ or
$\boldsymbol{\beta}\in\mathbb{R}_{+}^{p}$, where $l$ is differentiable in the
augmented parameter space $\mathbb{R}^{p}$ but not necessarily convex. Because
$f$ is directionally differentiable in $\mathbb{R}_{+}^{p}$, any optimal
solution $\hat{\boldsymbol{\beta}}$ must obey

$\displaystyle\mathrm{D}_{\boldsymbol{u}}f(\hat{\boldsymbol{\beta}})\geq
0\mbox{ for all feasible }\boldsymbol{u}$

where $\mathrm{D}_{\boldsymbol{u}}f(\boldsymbol{\beta})$ denotes the one-sided
directional derivative of $f$ at $\boldsymbol{\beta}$ with increment
$\boldsymbol{u}$, namely,
$\mathrm{D}_{\boldsymbol{u}}f(\boldsymbol{\beta})=\lim_{\epsilon\rightarrow
0+}[f(\boldsymbol{\beta}+\epsilon\boldsymbol{u})-f(\boldsymbol{\beta})]/\epsilon$.
Nevertheless, unlike equalities that are maintained after projection (i.e.,
$\sum{\boldsymbol{T}}_{i}(\boldsymbol{B})=\boldsymbol{0}\Rightarrow\sum\langle\boldsymbol{V},{\boldsymbol{T}}_{i}(\boldsymbol{B})\rangle=\langle\boldsymbol{V},\sum{\boldsymbol{T}}_{i}(\boldsymbol{B})\rangle=0,\forall\boldsymbol{V}\in\mathcal{G}$),
applying the same operation on inequalities may destroy their meanings totally
during the process of Tukeyfication.

Our proposal is to associate each inequality with an additional slack
variable, and append a nonnegative constraint when performing projection and
error measurement. Let $\boldsymbol{e}_{j}$ be a vector with the $j$th
component $1$ and the remaining 0. In Example 2, taking
$\boldsymbol{u}=\pm\boldsymbol{e}_{j}$ for $j\in\mathcal{J}=\\{j:\beta_{j}\neq
0\\}$ and $\boldsymbol{u}=\boldsymbol{e}_{j}$ for $j\in{\mathcal{J}}^{c}$
leads to the following slacked estimating equation:

$\sum_{i}(\nabla
l(\boldsymbol{\beta};\boldsymbol{x}_{i},y_{i})-\boldsymbol{s}/n)=\boldsymbol{0},$

where $\boldsymbol{s}_{\mathcal{J}^{c}}\succeq\boldsymbol{0}$ and
$\boldsymbol{s}_{\mathcal{J}}=\boldsymbol{0}$. The ordinary Tukeyfication now
goes through, and we obtain a depth optimization problem for any
$\boldsymbol{\beta}^{\circ}\succeq\boldsymbol{0}$:

$\min_{(\boldsymbol{v},\boldsymbol{s})\in\mathbb{R}^{p}\times\mathbb{R}^{p}}\sum_{i}1_{\geq
0}(\langle\boldsymbol{v},\allowbreak\nabla
l(\boldsymbol{\beta}^{\circ};\boldsymbol{x}_{i},y_{i})-\boldsymbol{s}/n\rangle)\;\textrm{s.t.}\;\lVert\boldsymbol{v}\rVert_{2}=1,\;\boldsymbol{s}\circ\boldsymbol{\beta}^{\circ}=\boldsymbol{0},\;\boldsymbol{s}\succeq\boldsymbol{0}.$

When
$l(\boldsymbol{\beta};\boldsymbol{x}_{i},y_{i})=(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}-y_{i})^{2}/2$,
we get the nonnegative regression depth

$\begin{gathered}\min_{(\boldsymbol{v},\boldsymbol{s})\in\mathbb{R}^{p}\times\mathbb{R}^{p}}\sum_{i}1_{\geq
0}(\langle\boldsymbol{v},\boldsymbol{x}_{i}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}^{\circ}-y_{i})-\boldsymbol{s}/n\rangle)\textrm{
s.t.
}\lVert\boldsymbol{v}\rVert_{2}=1,\boldsymbol{s}\circ\boldsymbol{\beta}^{\circ}=\boldsymbol{0},\boldsymbol{s}\succeq\boldsymbol{0}.\end{gathered}$
(15)

Recall that $\circ$ denotes the elementwise product. When
$\boldsymbol{\beta}^{\circ}\succ\boldsymbol{0}$,
$\boldsymbol{s}=\boldsymbol{0}$, and (15) becomes the regression depth. In
general, the inclusion of $\boldsymbol{s}$ in the minimization, as an outcome
of the nonnegativity restriction, often results in a lower depth value.

The slack-variable technique can introduce useful depth notions for sparse
learning that is at the core of high dimensional statistics:

$\min_{\boldsymbol{\beta}}f(\boldsymbol{\beta})\triangleq\bar{l}(\boldsymbol{X}\boldsymbol{\beta};\boldsymbol{y})+\sum_{j=1}^{p}P(\lvert\beta_{j}\rvert;\lambda),$
(16)

where
$\bar{l}(\boldsymbol{X}\boldsymbol{\beta};\boldsymbol{y})=\sum_{i}l_{0}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta};\boldsymbol{y}_{i})$
with $l_{0}$ differentiable. Here, we assume that $P$ is sparsity-promoting in
the sense that it is induced by a thresholding rule $\Theta(\cdot;\lambda)$
with $\lambda$ as the threshold (see She, (2012) for the rigorous definition
and more details): $P(t;\lambda)=P_{\Theta}(t;\lambda)+q(t;\lambda)$, where

$P_{\Theta}(t;\lambda)=\int_{0}^{\lvert
t\rvert}(\Theta^{-1}(u;\lambda)-u)\,\mathrm{d}u\mbox{ with
}\Theta^{-1}(u;\lambda)=\sup\\{t:\Theta(t;\lambda)\leq u\\}$

and $q$ is an arbitrary nonnegative function satisfying $q(t;\lambda)=0$ if
$t=\Theta(s;\lambda)$ for some $s\in\mathbb{R}$. Hence if
$\Theta(\cdot;\lambda)$ is a continuous function, $q$ must be identical to
zero, but if $\Theta$ has discontinuities, the mapping from $P$ to $\Theta$ is
many-to-one. The universal $\Theta$-$P$ framework covers many practically used
penalties such as $\ell_{r}$ ($0\leq r\leq 1$), SCAD, MCP, which can be
nonconvex. For centered response and predictors, (16) suffices; when centering
the response is inappropriate, an intercept $\alpha$ subject to no
regularization should often be added in the systematic component. For clarity,
we assume $\alpha=0$ in the following derivation, but the extension to
$\boldsymbol{X}\boldsymbol{\beta}+\alpha\boldsymbol{1}$ is straightforward.

For penalties with $q\equiv 0$ (continuous $\Theta$), like $\ell_{1}$ and
SCAD, we can use the directional derivatives along $\pm\boldsymbol{e}_{j}$ to
show that any locally optimal $\hat{\boldsymbol{\beta}}$ satisfies the
thresholding equation (She,, 2016)

${\boldsymbol{\beta}}={\Theta}({\boldsymbol{\beta}}-\boldsymbol{X}^{T}\nabla{\bar{l}}(\boldsymbol{X}{\boldsymbol{\beta}});\lambda),$
(17)

under the mild assumption that $\Theta(\cdot;\lambda)$ is continuous at
$\hat{\boldsymbol{\beta}}-\boldsymbol{X}^{T}\nabla{\bar{l}}(\boldsymbol{X}\hat{\boldsymbol{\beta}})$.
But nontrivial $q$’s and discontinuous $\Theta$’s constitute an important
class of nonsmooth penalties, including, in particular, the discontinuous
$\ell_{0}$ penalty

$\frac{\lambda^{2}}{2}\|\boldsymbol{\beta}\|_{0},$

for which $\Theta$ is the hard-thresholding
$\Theta_{H}(t;\lambda)=t1_{|t|>\lambda}$, and
$q(t;\lambda)=(1/2)(\lambda-|t|)^{2}1_{0<|t|<\lambda}$. In such scenarios, if
$\nabla\bar{l}$ is $L$-Lipschitz continuous, the solutions can be
characterized by the _fixed points_ of an iterative optimization algorithm
based on a surrogate function $g$:

$\beta\in\arg\min
g(\cdot,\boldsymbol{\beta}^{-})|_{\boldsymbol{\beta}^{-}=\boldsymbol{\beta}}$

where
$g(\boldsymbol{\beta},\boldsymbol{\beta}^{-})=\bar{l}(\boldsymbol{X}\boldsymbol{\beta}^{-})+\langle\nabla\bar{l}(\boldsymbol{X}\boldsymbol{\beta}^{-}),\boldsymbol{X}(\boldsymbol{\beta}-\boldsymbol{\beta}^{-})\rangle+\sum
P(\lvert\beta_{j}\rvert;\lambda)+\rho\|\boldsymbol{\beta}-\boldsymbol{\beta}^{-}\|_{2}^{2}/2$
with $\rho\geq L\|\boldsymbol{X}\|_{2}^{2}$, and we can prove that they all
enjoy nearly minimax error rate under a proper choice of $\lambda$ and some
regularity conditions (She et al., 2021b, ). It can be verified that the
fixed-point estimators satisfy (17) as well, under the same continuity
assumption, as long as $\boldsymbol{X}$ has been properly scaled:
$\|\boldsymbol{X}\|_{2}\leq 1/\sqrt{L}$. Unfortunately, the thresholding
equation does not belong to the estimating equation framework examined in our
companion paper—specifically, the nonlinear thresholding effect desired in
sparse learning means that (17) is not sample additive.

Below we introduce $p$ additional slack variables to find a proper substitute
for (17) so that one can define data depth for (16) given an arbitrary
thresholding $\Theta$. Let $\boldsymbol{\beta}$ be a locally optimal solution
to the problem as $P=P_{\Theta}$, or a fixed-point solution as
$P=P_{\Theta}+q$. Define $\mathcal{J}=\\{j:\beta_{j}\neq 0\\}$ and
$\mathcal{J}^{c}=\\{j:\beta_{j}=0\\}$, and denote by $\boldsymbol{X}[,j]$ the
$j$th column vector of $\boldsymbol{X}$. Using the directional derivatives of
$f$ when $P=P_{\Theta}$ (see, e.g., the proof of Theorem 1 in She, (2016)), or
the directional derivatives of $g-q$ under the continuity assumption when
$P=P_{\Theta}+q$, we get

$\displaystyle{\Theta}^{-1}(|{\beta}_{j}|;\lambda)\operatorname{sgn}({\beta}_{j})={\beta}_{j}-\boldsymbol{X}[,j]^{T}\nabla{\bar{l}}(\boldsymbol{X}\boldsymbol{\beta}),\
\forall j\in\mathcal{J},$ (18)

which holds even if $\Theta$ is not strictly increasing in a neighborhood of
$|\beta_{j}|$ ($j\in\mathcal{J}$), while for $j\in\mathcal{J}^{c}$,
$\beta_{j}=0$, and so

$\displaystyle-\lambda\leq\boldsymbol{X}[,j]^{T}\nabla{\bar{l}}(\boldsymbol{X}\boldsymbol{\beta})\leq\lambda,\forall
j\in\mathcal{J}^{c}.$ (19)

Next, define $\boldsymbol{\gamma}(\boldsymbol{\beta})=[\gamma_{j}]$ with

$\displaystyle\gamma_{j}=\begin{cases}\Theta^{-1}(|\beta_{j}|;\lambda)\operatorname{sgn}(\beta_{j})-\beta_{j}&\mbox{
if }j\in\mathcal{J}\\\ 0&\mbox{ if }j\in\mathcal{J}^{c}.\end{cases}$ (20)

It follows from (18), (19) that
$\boldsymbol{X}^{T}\nabla{\bar{l}}(\boldsymbol{X}\boldsymbol{\beta})+\boldsymbol{\gamma}(\boldsymbol{\beta})+\boldsymbol{s}=\boldsymbol{0}$
for some $\boldsymbol{s}\in\mathbb{R}^{p}$,
$\boldsymbol{s}_{\mathcal{J}}=\boldsymbol{0}$, and $|s_{j}|\leq\lambda$,
$j\in\mathcal{J}^{c}$.

Now, given a penalty induced by a thresholding rule $\Theta(\cdot;\lambda)$
and a point of interest $\boldsymbol{\beta}^{\circ}\in\mathbb{R}^{p}$, the
slacked data depth resulting from (16), which we call “$\Theta$-depth”, can be
cast as a joint optimization problem with respect to direction
$\boldsymbol{v}\in\mathbb{R}^{p}$ and slack variables
$\boldsymbol{s}=[{s}_{j}]\in\mathbb{R}^{p}$:

$\begin{split}\mbox{{$\Theta$-depth:}}\quad
d_{01}^{\Theta}(\boldsymbol{\beta}^{\circ})=\min_{(\boldsymbol{v},\boldsymbol{s})}&\sum_{i}1_{\geq
0}(\langle\boldsymbol{v},\boldsymbol{x}_{i}l_{0}^{\prime}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}^{\circ};y_{i})+(\boldsymbol{\gamma}^{\circ}+\boldsymbol{s})/n\rangle)\\\
&\textrm{ s.t.
}\lVert\boldsymbol{v}\rVert_{2}=1,\;\boldsymbol{s}\circ\boldsymbol{\beta}^{\circ}=\boldsymbol{0},\;\lVert\boldsymbol{s}\rVert_{\infty}\leq\lambda,\end{split}$
(21)

where
$\boldsymbol{\gamma}^{\circ}=\boldsymbol{\gamma}(\boldsymbol{\beta}^{\circ})$
and $\|\boldsymbol{s}\|_{\infty}=\max|s_{j}|$. When $\Theta$ is the hard
thresholding $\Theta_{H}$ (corresponding to the class of $\ell_{0}$
penalties), $\boldsymbol{\gamma}^{\circ}=\boldsymbol{0}$. The user should
specify a reasonably small $\lambda$ (otherwise extremely low depth values are
to be expected): a theoretical choice in sparse regression is
$\lambda=\sigma\sqrt{cn\log p}$ (with say $c=2$) where $\sigma$ is the Orlicz
$\psi_{2}$-norm of the noise, and a less conservative one can often be
obtained via cross-validation. A fascinating fact is that the slack variable
approach requires no convexity of either the loss or the penalty.

An important alternative to penalized sparse learning is to directly limit the
sparsity level: $\|\boldsymbol{\beta}\|_{0}\leq q$, instead of specifying a
penalty parameter $\lambda$. Due to the lack of nonsmoothness of

$\displaystyle\min_{\|\boldsymbol{\beta}\|_{0}\leq
q}\bar{l}(\boldsymbol{X}\boldsymbol{\beta};\boldsymbol{y}),$ (22)

we take the surrogate route. Statistically accurate estimates can be obtained
from the resulting iterative quantile-thresholding algorithm (She et al.,
2021c, ), which all satisfy
${\boldsymbol{\beta}}={\Theta}^{\\#}({\boldsymbol{\beta}}-(1/\rho)\boldsymbol{X}^{T}\nabla{\bar{l}}(\boldsymbol{X}{\boldsymbol{\beta}});q),$
assuming no ties occur and $\rho$ is large enough (e.g.,
$L\|\boldsymbol{X}\|_{2}^{2}$). Here, the quantile thresholding
$\Theta^{\\#}(\boldsymbol{\alpha};q)$ for any
$\boldsymbol{\alpha}\in\mathbb{R}^{p}$ is a vector $\boldsymbol{\zeta}$ with
$\zeta_{(j)}=\alpha_{(j)}\mbox{ if }1\leq j\leq q\mbox{, and }0\mbox{
otherwise,}$ where $\alpha_{(1)},\ldots,\alpha_{(p)}$ are the order statistics
of $\alpha_{1},\ldots,\alpha_{p}$ satisfying
$|\alpha_{(1)}|\geq\cdots\geq|\alpha_{(p)}|$.
$\Theta^{\\#}(\boldsymbol{\alpha};q)$ can be viewed as a variant of
$\Theta_{H}(\boldsymbol{\alpha};\lambda)$ (by setting
$\lambda=|\alpha_{(q+1)}|$, say), but it uses an adaptive threshold. Again, we
suppose that the regularization parameter $q$ is already given.

By use of slack variables to rewrite the $\Theta^{\\#}$-equation (details
omitted), we can define the $q$-sparse constrained $\ell_{0}$-depth (which we
call “$\Theta^{\\#}$-depth”) for any
$\boldsymbol{\beta}^{\circ}:\|\boldsymbol{\beta}^{\circ}\|_{0}=q$ as

$\displaystyle\begin{split}\mbox{{$\Theta^{\\#}$-depth:}}&\
\min_{(\boldsymbol{v},\boldsymbol{s})\in\mathbb{R}^{p}\times\mathbb{R}^{p}}\sum_{i}1_{\geq
0}(\langle\boldsymbol{v},\boldsymbol{x}_{i}l_{0}^{\prime}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}^{\circ};y_{i})+\boldsymbol{s}/n\rangle)\\\
&\mbox{ s.t.
}\lVert\boldsymbol{v}\rVert_{2}=1,\;\boldsymbol{s}\circ\boldsymbol{\beta}^{\circ}=\boldsymbol{0},\;\lVert\boldsymbol{s}\rVert_{\infty}\leq\|\boldsymbol{X}_{({\mathcal{J}}^{\circ})^{c}}^{T}\nabla{\bar{l}}(\boldsymbol{X}{\boldsymbol{\beta}}^{\circ})\|_{\infty},\end{split}$
(23)

where ${\mathcal{J}}^{\circ}=\\{j:\beta_{j}^{\circ}=0\\}$,
$\boldsymbol{X}_{({\mathcal{J}}^{\circ})^{c}}$ is a submatrix of
$\boldsymbol{X}$ by selecting the columns corresponding to the complement of
${\mathcal{J}}^{\circ}$, and both $\boldsymbol{\gamma}^{\circ}$ and $\rho$ (as
long as $\rho>0$) disappear in the $\ell_{0}$-constrained depth, just like in
the $\ell_{0}$-penalized case. Of all the constraints on $\boldsymbol{s}$, the
equality ones are affine, and the inequality ones are convex. The deepest
$q$-sparse estimate is defined as the saddle point that maximizes (23) over
all $\boldsymbol{\beta}^{\circ}:\|\boldsymbol{\beta}^{\circ}\|_{0}=q$ (cf.
(LABEL:eqcompositedep) in She et al., 2021a ).

Clearly, in the special case of $q=p$, all slack variables are removed, but as
$q<p$, the constrained problem (22) results in more stringent estimating
equations that are easier to violate, compared with the plain (non-
regularized) problem. This is reflected by the inclusion of $\boldsymbol{s}$
during the minimization, thereby lower depth values. The same conclusion holds
for the $\Theta$-depth (21) due to the existence of additional slack
variables. On the other hand, sparsity depths may be very low for large $p$
and $\lambda$ in (21) or large $p$ and small $q$ in (23). To alleviate the
issue, it is beneficial to change the crude “0-1 loss” to some more elegant
$\varphi$, as discussed in our companion paper (She et al., 2021a, ). For
example, (23) could be replaced by
$\min_{(\boldsymbol{v},\boldsymbol{s})}\sum_{i}\varphi(\langle\boldsymbol{v},\boldsymbol{x}_{i}l_{0}^{\prime}(\boldsymbol{x}_{i}^{T}\boldsymbol{\beta}^{\circ};y_{i})/\rho+\boldsymbol{s}/n\rangle)\mbox{
s.t.
}\lVert\boldsymbol{v}\rVert_{2}=1,\;\boldsymbol{s}\circ\boldsymbol{\beta}^{\circ}=\boldsymbol{0},\;\lVert\boldsymbol{s}\rVert_{\infty}\leq\|\boldsymbol{X}_{({\mathcal{J}}^{\circ})^{c}}^{T}\nabla{\bar{l}}(\boldsymbol{X}{\boldsymbol{\beta}}^{\circ})/\rho\|_{\infty}$,
with $\varphi(\cdot)$ nonzero for mild or moderate negative inputs, which
warrants further investigation in the future. A tight upper bound of sparsity
depths is also worth studying in theory.

Slacked data depth can be introduced for groupwise variable selection and low-
rank matrix estimation (She,, 2012, 2013) as well; see, e.g., Section 3.2.

######  Remark 2 (Computation of slacked depth).

A simple alternating optimization or _block coordinate descent_ (BCD)
algorithm can be used to to compute slacked data depth. Take the
$\varphi$-form of (21) as an example. Given $\boldsymbol{s}$, the optimization
problem for $\boldsymbol{v}$,

$\min_{\boldsymbol{v}\in\mathbb{R}^{p},\lVert\boldsymbol{v}\rVert_{2}=1}\allowbreak\sum_{i}\varphi(\langle\boldsymbol{v},r_{i}\boldsymbol{x}_{i}+(\boldsymbol{\gamma}^{\circ}+\boldsymbol{s})/n\rangle)$

where
$\boldsymbol{r}=\nabla_{\boldsymbol{\Theta}}\bar{l}\lvert_{\boldsymbol{\Theta}=\boldsymbol{X}\boldsymbol{\beta}^{\circ}}$,
has been investigated in She et al., 2021a . Fixing $\boldsymbol{v}$, we can
rewrite the $\boldsymbol{s}$-problem as

$\min_{\boldsymbol{s}\in\mathbb{R}^{p}}\sum_{i}\varphi(\langle\boldsymbol{v},r_{i}\boldsymbol{x}_{i}+\boldsymbol{\gamma}^{\circ}/n+\boldsymbol{s}/n\rangle)\;\textrm{
s.t.
}\;\boldsymbol{s}_{\mathcal{J}}=\boldsymbol{0},\lVert\boldsymbol{s}\rVert_{\infty}\leq\lambda.$

The problem has a differentiable criterion in $\boldsymbol{s}$ and some simple
box constraints, and conventional numerical methods apply, including L-BFGS-B,
interior point, and proximal gradient descent algorithms (Byrd et al.,, 1995;
Boyd and Vandenberghe,, 2004; Parikh and Boyd,, 2014).

###  3.2 Reduced-rank regression depth

Applying ordinary least squares on multiple responses may easily result in a
large number of unknowns. Researchers often prefer adding a low-rank
constraint in estimating the coefficient matrix, leading to the celebrated
reduced-rank regression (RRR) (Anderson,, 1951)

$\displaystyle\min_{\boldsymbol{B}\in\mathbb{R}^{p\times
m}}f(\boldsymbol{B};\boldsymbol{X},\boldsymbol{Y})\triangleq\frac{1}{2}\|\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{B}\|_{F}^{2}\mbox{
s.t. }\mbox{rank}(\boldsymbol{B})\leq r,$ (24)

where
$\boldsymbol{Y}=[\boldsymbol{y}_{1}\ldots\boldsymbol{y}_{n}]^{T}\in\mathbb{R}^{n\times
m}$ and
$\boldsymbol{X}=[\boldsymbol{x}_{1}\ldots\boldsymbol{x}_{n}]^{T}\in\mathbb{R}^{n\times
p}$ are the (centered) response and predictor matrices. A weighted criterion
to account for the dependency between the responses can be given, but the
problem can be converted to (24) with a simple reparametrization. If the
variables are not centered, an intercept term
$\boldsymbol{1}\boldsymbol{\alpha}^{T}$ should be added in the loss, but the
depth derivation below carries over (cf. Section 2.2). We assume that
$(\boldsymbol{x}_{i},\boldsymbol{y}_{i})$ are i.i.d. (or in an approximate
sense), and so the data depth in this subsection does not apply to PCA where
$\boldsymbol{X}=\boldsymbol{I}$, thus distinct from the PC-depth and OC-depth
introduced earlier; see some related discussions in Section LABEL:subsec:ex of
our companion paper.

RRR provides a low-dimensional projection space to view and analyze supervised
multivariate data, and finds widespread applications in machine learning and
econometrics (Reinsel and Velu,, 1998; Izenman,, 2008). In fact, once an
estimate $\boldsymbol{B}$ of rank $r$ is obtained, we can write
$\boldsymbol{B}=\boldsymbol{B}_{1}\boldsymbol{B}_{2}^{T}$ for
$\boldsymbol{B}_{1}\in\mathbb{R}^{p\times r}$,
$\boldsymbol{B}_{2}\in\mathbb{R}^{m\times r}$. This suggests that $r$ factors
can be constructed by $\boldsymbol{X}\boldsymbol{B}_{1}$ from $p$ predictors
to explain all response variables. The number of factors required in real
applications is often much smaller than the number of input $x$-variables.

Limiting the rank of the matrix estimators at $r$, how to perform a “center-
outward” ranking in high dimensions, or more generally, test

$H_{0}:\boldsymbol{B}\in\Omega_{0}\cap\\{\mbox{rank}(\boldsymbol{B})=r\\}\mbox{
vs.
}H_{a}:\boldsymbol{B}\in\Omega_{0}^{c}\cap\\{\mbox{rank}(\boldsymbol{B})=r\\},$

where the set or event $\Omega_{0}$ is not necessarily a singleton (cf. Remark
LABEL:rem:depth_median_computation in She et al., 2021a ), is an intriguing
and open question.

In the following, we extend multivariate regression depth (Rousseeuw and
Hubert,, 1999; Bern and Eppstein,, 2002) to the reduced-rank regression depth
(27), using the techniques developed in the last subsection. Toward this, we
first give a fixed-point formulation of all RRR estimators. Define a matrix
version of the $\Theta^{\\#}$ introduced in the last subsection

$\Theta^{\sigma\\#}(\boldsymbol{B};r)\triangleq\boldsymbol{U}{\rm{diag}}\\{\Theta^{\\#}([\sigma_{i}{(\boldsymbol{B})}];r)\\}\boldsymbol{V}^{T},\
\ \ \forall\boldsymbol{B}\in\mathbb{R}^{p\times m}$ (25)

where $\boldsymbol{U}$, $\boldsymbol{V}$, and
diag$\\{\sigma{(\boldsymbol{B})}_{i}\\}$ are from the SVD of
$\boldsymbol{B}=\boldsymbol{U}\text{diag}\\{\sigma_{i}{(\boldsymbol{B})}\\}\boldsymbol{V}$,
and $\Theta^{\\#}$ is applied to the vector $[\sigma_{i}{(\boldsymbol{B})}]$,
with $\sigma_{i}{(\boldsymbol{B})}$ denoting the $i$th largest singular value
of $\boldsymbol{B}$.

Construct a surrogate function

$g(\boldsymbol{B},\boldsymbol{B}^{-})=f(\boldsymbol{B}^{-})+\langle\nabla
f(\boldsymbol{B}^{-}),\boldsymbol{B}-\boldsymbol{B}^{-}\rangle+\rho\|\boldsymbol{B}-\boldsymbol{B}^{-}\|_{F}^{2}/2,$

where $\rho$ is larger than $\|\boldsymbol{X}\|_{2}^{2}$. Let
$\hat{\boldsymbol{B}}_{\mbox{\tiny rrr}}$ be an RRR estimator that solves
(24). Then for $\tilde{\boldsymbol{B}}\in\arg\min_{\boldsymbol{B}:\mbox{\tiny
rank}(\boldsymbol{B})\leq r}g(\boldsymbol{B},\hat{\boldsymbol{B}}_{\mbox{\tiny
rrr}})$, it follows from the chain inequalities
$f(\hat{\boldsymbol{B}}_{\mbox{\tiny
rrr}})-f(\tilde{\boldsymbol{B}})=g(\hat{\boldsymbol{B}}_{\mbox{\tiny
rrr}},\hat{\boldsymbol{B}}_{\mbox{\tiny rrr}})-f(\tilde{\boldsymbol{B}})\geq
g(\tilde{\boldsymbol{B}},\hat{\boldsymbol{B}}_{\mbox{\tiny
rrr}})-f(\tilde{\boldsymbol{B}})\geq(\rho-\|\boldsymbol{X}\|_{2}^{2})\|\tilde{\boldsymbol{B}}-\hat{\boldsymbol{B}}_{\mbox{\tiny
rrr}}\|_{F}^{2}/2$ that $\hat{\boldsymbol{B}}_{\mbox{\tiny
rrr}}=\tilde{\boldsymbol{B}}$. On the other hand, it is easy to show that
$\tilde{\boldsymbol{B}}=\Theta^{\sigma\\#}(\hat{\boldsymbol{B}}_{\mbox{\tiny
rrr}}-\frac{1}{\rho}\boldsymbol{X}^{T}(\boldsymbol{X}\hat{\boldsymbol{B}}_{\mbox{\tiny
rrr}}-\boldsymbol{Y});r)$ (She,, 2013), and so
$\hat{\boldsymbol{B}}_{\mbox{\tiny rrr}}$ satisfies the matrix thresholding
equation

$\boldsymbol{B}=\Theta^{\sigma\\#}({\boldsymbol{B}}-\frac{1}{\rho}\boldsymbol{X}^{T}(\boldsymbol{X}\boldsymbol{B}-\boldsymbol{Y});r).$
(26)

(In fact, under the mild condition that
$\boldsymbol{Y}^{T}\boldsymbol{X}(\boldsymbol{X}^{T}\boldsymbol{X})^{+}\boldsymbol{X}^{T}\boldsymbol{Y}$
has distinct eigenvalues, the RRR estimator is unique (Reinsel and Velu,,
1998), and $\rho$ can be way smaller than $\|\boldsymbol{X}\|_{2}^{2}$.)
Perform a compact SVD of $\boldsymbol{B}$:
$\boldsymbol{B}=\boldsymbol{P}\boldsymbol{D}\boldsymbol{Q}^{T}$ with
$\boldsymbol{P}\in\mathbb{O}^{p\times r}$ and
$\boldsymbol{Q}\in\mathbb{O}^{m\times r}$, and denote by
$\boldsymbol{P}_{\perp}\in\mathbb{O}^{p\times(p-r)}$ and
$\boldsymbol{Q}_{\perp}\in\mathbb{O}^{m\times(m-r)}$ their orthogonal
complements (which can be obtained from the full SVD of $\boldsymbol{B}$).
Like in the $\Theta^{\\#}$-case, based on (26), we work on
$({1}/{\rho})\boldsymbol{X}^{T}(\boldsymbol{X}\boldsymbol{B}-\boldsymbol{Y})+\boldsymbol{S}=\boldsymbol{0}$
for a slack matrix $\boldsymbol{S}$ satisfying

$\boldsymbol{P}^{T}\boldsymbol{S}=\boldsymbol{0},\boldsymbol{S}\boldsymbol{Q}=\boldsymbol{0},\|\boldsymbol{S}\|_{2}\leq\|\frac{1}{\rho}\boldsymbol{P}_{\perp}\boldsymbol{P}_{\perp}^{T}\boldsymbol{X}^{T}(\boldsymbol{X}\boldsymbol{B}-\boldsymbol{Y})\boldsymbol{Q}_{\perp}\boldsymbol{Q}_{\perp}^{T}\|_{2}.$

Now, given a regularization parameter $r:1\leq r\leq p\wedge m$ and a matrix
of interest $\boldsymbol{B}^{\circ}\in\mathbb{R}^{p\times
m}:\mbox{rank}(\boldsymbol{B}^{\circ})=r$, obtain the associated
$\boldsymbol{P}_{\perp}^{\circ}\in\mathbb{O}^{p\times(p-r)}$,
$\boldsymbol{Q}_{\perp}^{\circ}\in\mathbb{O}^{m\times(m-r)}$ as above; the
rank-$r$ RRR depth of $\boldsymbol{B}^{\circ}$ is defined by

$\displaystyle
d_{01}^{\mbox{{RRR}}}(\boldsymbol{B}^{\circ})=\min_{(\boldsymbol{V},\boldsymbol{S})}\sum_{i}1_{\geq
0}(\langle\boldsymbol{V},\frac{1}{\rho}\boldsymbol{x}_{i}(\boldsymbol{x}^{T}_{i}\boldsymbol{B}^{\circ}-\boldsymbol{y}_{i}^{T})+\frac{1}{n}\boldsymbol{S}\rangle)$
$\displaystyle\ \text{ s.t. }\|\boldsymbol{V}\|_{F}=1,\boldsymbol{P}^{\circ
T}\boldsymbol{S}=\boldsymbol{0},\boldsymbol{S}\boldsymbol{Q}^{\circ}=\boldsymbol{0},\|\boldsymbol{S}\|_{2}\leq\|\frac{1}{\rho}\boldsymbol{P}_{\perp}^{\circ
T}\boldsymbol{X}^{T}(\boldsymbol{X}\boldsymbol{B}^{\circ}-\boldsymbol{Y})\boldsymbol{Q}_{\perp}^{\circ}\|_{2},$

or equivalently,

RRR-depth: $\displaystyle
d_{01}^{\mbox{{RRR}}}(\boldsymbol{B}^{\circ})=\min_{(\boldsymbol{V},\boldsymbol{L})}\sum_{i}1_{\geq
0}(\langle\boldsymbol{V},\boldsymbol{x}_{i}(\boldsymbol{x}^{T}_{i}\boldsymbol{B}^{\circ}-\boldsymbol{y}_{i}^{T})+\frac{1}{n}\boldsymbol{P}_{\perp}^{\circ}\boldsymbol{L}\boldsymbol{Q}^{\circ
T}_{\perp}\rangle)$ $\displaystyle\ \text{ s.t.
}\|\boldsymbol{V}\|_{F}=1,\|\boldsymbol{L}\|_{2}\leq\|\boldsymbol{P}_{\perp}^{\circ
T}\boldsymbol{X}^{T}(\boldsymbol{X}\boldsymbol{B}^{\circ}-\boldsymbol{Y})\boldsymbol{Q}_{\perp}^{\circ}\|_{2},$
(27)

where $\rho$ vanishes due to the scale invariance of $1_{\geq 0}$, regardless
of how small or large $\rho$ is. Clearly, in the full rank case $r=m\wedge p$,
either $\boldsymbol{P}_{\perp}^{\circ}$ or $\boldsymbol{Q}_{\perp}^{\circ}$
must vanish, and so $\boldsymbol{L}=\boldsymbol{0}$, meaning that (27) reduces
to the multivariate regression depth (Bern and Eppstein,, 2002).

######  Remark 3 (Combined treatment).

The manifold approach and slack variable approach can be combined together to
define data depth for some challenging problems. Consider a _sparse_ RRR (one
of the variants in She, (2017)) that constructs $r$ predictive factors from a
subset of predictors

$\displaystyle\min_{\boldsymbol{S}\in\mathbb{R}^{p\times
r},\boldsymbol{U}\in\mathbb{R}^{m\times
r}}\frac{1}{2}\|\boldsymbol{Y}-\boldsymbol{X}\boldsymbol{A}\boldsymbol{U}^{T}\|_{F}^{2}\mbox{
s.t. }\|\operatorname{\mbox{vec}\,}(\boldsymbol{A})\|_{0}\leq q,\
\boldsymbol{U}^{T}\boldsymbol{U}=\boldsymbol{I}_{r\times r}.$ (28)

The overall coefficient matrix
$\boldsymbol{B}=\boldsymbol{A}\boldsymbol{U}^{T}$ has rank at most $r$ as in
RRR, but sparsity is imposed on the loading matrix $\boldsymbol{A}$. By use of
a slack matrix $\boldsymbol{S}$ for $\boldsymbol{A}$, and a Riemannian tangent
space for $\boldsymbol{U}$, the depth for
$(\boldsymbol{A}^{\circ},\boldsymbol{U}^{\circ})$:
$\boldsymbol{A}^{\circ}\in\mathbb{R}^{p\times
r},\|\operatorname{\mbox{vec}\,}(\boldsymbol{A}^{\circ})\|_{0}=q$,
$\boldsymbol{U}^{\circ}\in\mathbb{O}^{m\times r}$ is given by
$\min_{\boldsymbol{W}\in\mathbb{R}^{m\times
r},\boldsymbol{V}\in\mathbb{R}^{p\times
r},\boldsymbol{S}\in\mathbb{R}^{p\times r}}\allowbreak\sum_{i=1}^{n}1_{\geq
0}(-\langle\boldsymbol{W},\boldsymbol{y}_{i}\boldsymbol{x}_{i}^{T}\boldsymbol{A}^{\circ}\rangle+\langle\boldsymbol{V},\boldsymbol{x}_{i}(\boldsymbol{x}_{i}^{T}\boldsymbol{A}^{\circ}-\boldsymbol{y}_{i}^{T}\boldsymbol{U}^{\circ})+\boldsymbol{S}/n\rangle)$
s.t.
$\|\boldsymbol{W}\|_{F}^{2}+\|\boldsymbol{V}\|_{F}^{2}=1,\boldsymbol{V}^{T}\boldsymbol{U}^{\circ}+\boldsymbol{U}^{\circ
T}\boldsymbol{V}=\boldsymbol{0},\operatorname{\mbox{vec}\,}(\boldsymbol{A}^{\circ})\circ\operatorname{\mbox{vec}\,}(\boldsymbol{S})=\boldsymbol{0},\|\boldsymbol{S}\|_{\max}\leq\lambda^{\circ}$,
with
$\lambda^{\circ}=\|\operatorname{\mbox{vec}\,}(\boldsymbol{X}^{T}(\boldsymbol{X}\boldsymbol{A}^{\circ}-\boldsymbol{Y}\boldsymbol{U}^{\circ}))[{({\mathcal{J}}^{\circ})^{c}}]\|_{\infty}=\|(\boldsymbol{I}\otimes\boldsymbol{X}^{T})[{({\mathcal{J}}^{\circ})^{c}},]\operatorname{\mbox{vec}\,}(\boldsymbol{X}\boldsymbol{A}^{\circ}-\boldsymbol{Y}\boldsymbol{U}^{\circ})\|_{\infty}$
and
${\mathcal{J}}^{\circ}=\\{j:\operatorname{\mbox{vec}\,}(\boldsymbol{A}^{\circ})[j]\neq
0,1\leq j\leq pr\\}$.

##  4 Experiments

This section performs real data experiments to illustrate the usefulness of
some new notions of depth.

###  4.1 Reduced-rank depth in time series

We consider the 52 weekly stock log-return data for nine of the ten largest
American corporations in 2004 (Rothman et al.,, 2010), with
$\boldsymbol{y}_{t}\in\mathbb{R}^{9}$ ($t=1,\ldots,T$) and $T=52$.

For the purpose of constructing market factors that drive general stock
movements, a reduced-rank vector autoregressive (VAR) model can be used, i.e.,
$\boldsymbol{y}_{t+1}=\boldsymbol{B}^{*T}\boldsymbol{y}_{t}+\boldsymbol{e}_{t}$,
with $\boldsymbol{B}^{*}$ of low rank. By conditioning on the initial state
$\boldsymbol{y}_{0}$ and assuming the normality of $\boldsymbol{e}_{t}$, the
conditional likelihood function leads to a least squares criterion, so the
estimation of $\boldsymbol{B}^{*}$ can be formulated as a reduced-rank
regression problem; see Lütkepohl, (2007) for more details. We fit the
reduced-rank VAR with $r=6$. The optimization algorithm for (27) (implemented
based on Remark 2) however shows that the objective function can reach zero
for some feasible $(\boldsymbol{V},\boldsymbol{L})$. Hence, although the
standard RRR approach is widely used in multivariate times series and
econometrics, our analysis revealed a perhaps surprisingly low data depth on
this financial series dataset.

We then considered the Cauchy-based reduced rank regression (Zhao and
Palomar,, 2017; Yang and Zhao,, 2020) (denoted by C-RRR) and robust reduced
rank regression (She and Chen,, 2017) (denoted by R4, with 5% of data treated
as outliers), as well as a deeper estimate obtained by random sampling
(denoted by D-RRR). The rank-6 depth values of these estimates are 0.02, 0.08
and 0.12, respectively, suggesting more reliable fitted models than the plain
RRR from the perspective of data depth.

Figure 1: The fitted models $y_{1,t+1}\sim x_{4,t}$ using different methods in
the low-rank VAR(1), to demonstrate how the log-return of Walmart is related
to that of Ford in the previous week. Notice the right-most point that has a
high leverage.

To further illustrate the differences between the estimates, we plot the
fitted models of Ford ($x_{4,t}$) in response to Walmart ($y_{1,t+1}$) in
Figure 1. Notably, the right-most point has high leverage, and the RRR model
passes close to that particular observation. In contrast, D-RRR seems to fit
better the majority of the sample.

A careful examination of the series shows the point corresponds to the log-
return of Ford at week 17, a real major market disturbance attributed to the
auto industry. Several other stock returns experienced dramatic short-term
changes as well, and we occasionally observe that the slopes obtained from RRR
and its robust counterparts can have opposite signs. Financial time series
often contain anomalies or demonstrate heavier tails than those of a normal
distribution due to extreme market movements. The issue may jeopardize the
recovery of common market behaviors and asset return forecasting: the
autoregressive structure can make any outlier in the time series also a
leverage point in the covariates. Although an elaborate robustification of the
low-rank VAR merits further investigation, our depth-based analysis seems to
offer an effective fix in this regard.

###  4.2 Sparsity depth for performance evaluation

Data depth provides a nonparametric means of performance evaluation. In this
experiment, we use the sparsity depth defined in (23) to conduct a comparison
between some commonly used sparse learning methods on the Boston housing
dataset (Harrison and Rubinfeld,, 1978). The dataset was collected by the U.S.
Census Service and consists of 13 predictors regarding socioeconomic and
environmental conditions for 506 neighborhoods in the Boston area. The
response is the median value of owner-occupied homes in the area.

We compare Lasso (Tibshirani,, 1996), SCAD (Fan and Li,, 2001), sparse LTS
(S-LTS) (Alfons et al.,, 2013), quantile-SCAD (Q-SCAD) (Sherwood and Maidman,,
2019) and PIQ (She et al., 2021c, ), in terms of data centrality defined in
(23) for the same given support size $q$. More concretely, assuming that the
observations are i.i.d., we split the dataset in halves, fit the methods on
the first half, and then evaluate their performance via sparsity depth on the
rest half. The whole procedure is repeated 20 times.

(a) $q=7$ (b) $q=8$ (c) $q=9$ (d) $q=10$

Figure 2: Sparsity depth comparison between Lasso, SCAD, sparse LTS (S-LTS),
quantile SCAD (Q-SCAD) and PIQ estimates with respect to the support size $q$.

Figure 2 shows a series of radar plots for the median $\Theta^{\\#}$-depths in
respect to the number of selected variables. The depth values are small but
present useful ranking information. The estimates are quite different seen
from the data depth comparison: Lasso and SCAD exhibit lower depth in most
cases, Q-SCAD and PIQ often give deeper estimates, and S-LTS is unstable (and
costly) in our experiments. The last three methods all use a more robust loss,
as well as a nonconvex regularizer, while Lasso solves a convex optimization
problem with the ordinary $\ell_{2}$-loss and $\ell_{1}$-penalty. The depth
differences between these sparse learning methods indicate that the data must
deviate from Gaussianity and may contain anomalies, and incorporating the
desired type of regularization into data depth can provide a helpful tool for
robust performance evaluation.

##  5 Summary

Our work investigated Tukey’s notion of depth for robustifying a given
optimization criterion, an estimating equation, or an algorithm in statistical
inference and estimation. In Part I, we introduced a polished subspace depth
framework, where the elements like the influence space constraint, rectified
redescending discrepancy measures, and subspace projection are new to the best
of our knowledge. In Part II, we proposed two novel approaches based on
manifolds and slack variables to extend the concept to problems defined in
some restricted parameter spaces or with a nonsmooth regularizer. Our matrix
formulation of the problems, together with state-of-the-art optimization
techniques (particularly momentum-based acceleration), gave rise to a new
class of efficient algorithms that has guaranteed convergence and scales up
with problem dimensions. The efficient computation of the deepest point or
composite depth (cf. Remark LABEL:rem:depth_median_computation of Part I) is
yet more difficult, and recent advances in nonconvex min-max optimization
(Razaviyayn et al.,, 2020) may shed new light on the topic.

The proposed computational inference tool caters to machine learning
applications beyond the standard likelihood setup. For example, given a
feedforward neural network, it can be used to evaluate the reliability of a
given estimate, or an event concerned with some properties of the unknowns,
which only requires the gradient information that can be obtained from back
propagation with ease. Moreover, the influence-driven deepest estimation
provides a universal means of accommodating distortions and anomalies given
any criterion or estimation equations. We hope that the work is helpful to
advance the practice of data depth in sophisticated setups and in higher
dimensions.

## References

  * Alfons et al., (2013) Alfons, A., Croux, C., and Gelper, S. (2013).  Sparse least trimmed squares regression for analyzing high-dimensional large data sets.  The Annals of Applied Statistics, 7(1):226–248.
  * Anderson, (1951) Anderson, T. W. (1951).  Estimating linear restrictions on regression coefficients for multivariate normal distributions.  Annals of Mathematical Statistics, 22:327–351.
  * Bern and Eppstein, (2002) Bern and Eppstein (2002).  Multivariate regression depth.  Discrete & Computational Geometry, 28(1):1–17.
  * Bijral et al., (2007) Bijral, A. S., Breitenbach, M., and Grudic, G. (2007).  Mixture of Watson distributions: A generative model for hyperspherical embeddings.  In Proceedings of the Eleventh International Conference on Artificial Intelligence and Statistics, pages 35–42, San Juan, Puerto Rico.
  * Boothby, (1986) Boothby, W. M. (1986).  An introduction to differentiable manifolds and Riemannian geometry.  Academic Press, Orlando, FL.
  * Boyd and Vandenberghe, (2004) Boyd, S. and Vandenberghe, L. (2004).  Convex optimization.  Cambridge University Press, New York, NY.
  * Byrd et al., (1995) Byrd, R., Lu, P., Nocedal, J., and Zhu, C. (1995).  A limited memory algorithm for bound constrained optimization.  SIAM Journal on Scientific Computing, 16(5):1190–1208.
  * Cai et al., (2009) Cai, T. T., Xu, G., and Zhang, J. (2009).  On recovery of sparse signals via $\ell_{1}$ minimization.  IEEE Transactions on Information Theory, 55(7):3388–3397.
  * Dhillon et al., (2003) Dhillon, I. S., Marcotte, E. M., and Roshan, U. (2003).  Diametrical clustering for identifying anti-correlated gene clusters.  Bioinformatics, 19(13):1612–1619.
  * Edelman et al., (1998) Edelman, A., Arias, T. A., and Smith, S. T. (1998).  The geometry of algorithms with orthogonality constraints.  SIAM Journal on Matrix Analysis and Applications, 20(2):303–353.
  * Fan and Li, (2001) Fan, J. and Li, R. (2001).  Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American Statistical Association, 96(456):1348–1360.
  * Harrison and Rubinfeld, (1978) Harrison, D. and Rubinfeld, D. (1978).  Hedonic housing prices and the demand for clean air.  Journal of Environmental Economics and Management, 5:81–102.
  * Izenman, (2008) Izenman, A. (2008).  Modern Multivariate Statistical Techniques: Regression, Classification and Manifold Learning.  Springer, New York.
  * Liu and Singh, (1992) Liu, R. Y. and Singh, K. (1992).  Ordering directional data: concepts of data depth on circles and spheres.  The Annals of Statistics, 20(3):1468–1484.
  * Lütkepohl, (2007) Lütkepohl, H. (2007).  New Introduction to Multiple Time Series Analysis.  Springer-Verlag Berlin Heidelberg.
  * Mardia and Jupp, (1999) Mardia, K. V. and Jupp, P. E. (1999).  Directional statistics.  John Wiley & Sons, Hoboken, NJ.
  * Mizera, (2002) Mizera, I. (2002).  On depth and deep points: a calculus.  Ann. Statist., 30(6):1681–1736.
  * Parikh and Boyd, (2014) Parikh, N. and Boyd, S. (2014).  Proximal algorithms.  Foundations and Trends in Optimization, 1(3):127–239.
  * Razaviyayn et al., (2020) Razaviyayn, M., Huang, T., Lu, S., Nouiehed, M., Sanjabi, M., and Hong, M. (2020).  Nonconvex min-max optimization: Applications, challenges, and recent theoretical advances.  IEEE Signal Processing Magazine, 37(5):55–66.
  * Reinsel and Velu, (1998) Reinsel, G. and Velu, R. (1998).  Multivariate Reduced-Rank Regression: Theory and Applications.  Springer, New York.
  * Rothman et al., (2010) Rothman, A. J., Levina, E., and Zhu, J. (2010).  Sparse multivariate regression with covariance estimation.  Journal of Computational and Graphical Statistics, 19(4):947–962.
  * Rousseeuw and Hubert, (1999) Rousseeuw, P. J. and Hubert, M. (1999).  Regression depth.  Journal of the American Statistical Association, 94(446):388–402.
  * She, (2012) She, Y. (2012).  An iterative algorithm for fitting nonconvex penalized generalized linear models with grouped predictors.  Computational Statistics & Data Analysis, 56(10):2976–2990.
  * She, (2013) She, Y. (2013).  Reduced rank vector generalized linear models for feature extraction.  Statistics and Its Interface, 6:197–209.
  * She, (2016) She, Y. (2016).  On the finite-sample analysis of $\Theta$-estimators.  Electron. J. Statist., 10(2):1874–1895.
  * She, (2017) She, Y. (2017).  Selective factor extraction in high dimensions.  Biometrika, 104(1):97–110.
  * She and Chen, (2017) She, Y. and Chen, K. (2017).  Robust reduced-rank regression.  Biometrika, 104(3):633–647.
  * She et al., (2016) She, Y., Li, S., and Wu, D. (2016).  Robust orthogonal complement principal component analysis.  Journal of the American Statistical Association, 111(514):763–771.
  * (29) She, Y., Tang, S., and Liu, L. (2021a).  On Generalization and Computation of Tukey’s Depth: Part I.  Journal of Data Science, Statistics, and Visualisation.  To appear.
  * She and Tran, (2019) She, Y. and Tran, H. (2019).  On cross-validation for sparse reduced rank regression.  Journal of the Royal Statistical Society: Series B, 81:145–161.
  * (31) She, Y., Wang, Z., and Jin, J. (2021b).  Analysis of Generalized Bregman Surrogate Algorithms for Nonsmooth Nonconvex Statistical Learning.  The Annals of Statistics.  To appear.
  * (32) She, Y., Wang, Z., and Shen, J. (2021c).  Gaining Outlier Resistance with Progressive Quantiles: Fast Algorithms and Theoretical Studies.  Journal of the American Statistical Association.  To appear.
  * Sherwood and Maidman, (2019) Sherwood, B. and Maidman, A. (2019).  rqPen: Penalized Quantile Regression.  R package version 2.1.
  * Sra and Karp, (2013) Sra, S. and Karp, D. (2013).  The multivariate watson distribution: Maximum-likelihood estimation and other aspects.  Journal of Multivariate Analysis, 114:256 – 269.
  * Tibshirani, (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological), 58(1):267–288.
  * Tukey, (1975) Tukey, J. W. (1975).  Mathematics and the picturing of data.  In Proceedings of the international congress of mathematicians, volume 2.
  * Watson, (1965) Watson, G. S. (1965).  Equatorial distributions on a sphere.  Biometrika, 52(1/2):193–201.
  * Yang and Zhao, (2020) Yang, Y. and Zhao, Z. (2020).  RRRR: Online Robust Reduced-Rank Regression Estimation.
  * Zhang, (2010) Zhang, C.-H. (2010).  Nearly unbiased variable selection under minimax concave penalty.  The Annals of Statistics, 38(2):894–942.
  * Zhao and Palomar, (2017) Zhao, Z. and Palomar, D. P. (2017).  Robust maximum likelihood estimation of sparse vector error correction model.  In 2017 IEEE Global Conference on Signal and Information Processing (GlobalSIP), pages 913–917. IEEE.


# TensorLy-Quantum: Quantum Machine Learning with Tensor Methods

Taylor L. Patti1,2<EMAIL_ADDRESS>
Jean Kossaifi1<EMAIL_ADDRESS>
Susanne F. Yelin2<EMAIL_ADDRESS>
Anima Anandkumar1,3<EMAIL_ADDRESS>
1NVIDIA
2Harvard University
3 California Institute of Technology

###### Abstract

Simulation is essential for developing quantum hardware and algorithms.
However, simulating quantum circuits on classical hardware is challenging due
to the exponential scaling of quantum state space. While factorized tensors
can greatly reduce this overhead, tensor network-based simulators are
relatively few and often lack crucial functionalities. To address this
deficiency, we created TensorLy-Quantum, a Python library for quantum circuit
simulation that adopts the PyTorch API. Our library leverages the optimized
tensor methods of the existing TensorLy ecosystem to represent, simulate, and
manipulate large-scale quantum circuits. Through compact tensor
representations and efficient operations, TensorLy-Quantum can scale to
hundreds of qubits on a single GPU and thousands of qubits on multiple GPUs.
TensorLy-Quantum is open-source and accessible at
https://github.com/tensorly/quantum

##  1 Introduction

The size and intricacies of quantum systems have long motivated quantum
computing (QC) research (Montanaro (2016)). Recently, quantum algorithms have
incorporated variational techniques (Peruzzo et al. (2014); Farhi et al.
(2014)), creating the subfield of quantum machine learning (QML) Mehta et al.
(2019), which includes applications like optimization and quantum chemistry.
Much like its classical counterpart, QML leverages gradient-based methods to
seek optimal solutions, but with the added benefit of an exponentially large
solution kernel (Schuld (2021)). For instance, the most common form of QML
uses quantum expectation values to define objective functions and classical
gradient descent to optimize quantum networks. Both QC and QML ultimately seek
to demonstrate a tangible improvement over classical computing for useful
tasks, a concept known as quantum advantage.

As QC and QML are still nascent fields, simulation remains essential to both
the characterization of quantum systems at scale (Carleo et al. (2019)), as
well as the design heuristic algorithms (Cerezo et al. (2021)). However,
simulating quantum algorithms on informative scales is a challenging task, as
quantum state space overwhelms classical resources with even relatively few
quantum bits (qubits). This limitation is partially overcome with networks of
factorized tensors (Bridgeman and Chubb (2017)), which are typically more
efficient than the traditional vector representation of quantum states,
providing up to an exponential reduction in overhead with respect to the
number of qubits $n$. In tensor network-based simulation, both quantum states,
gates, and operators are expressed in a factorized form, which compress the
quantum circuit by including only relevant ranks. Operations can be
efficiently done directly in the factorized form without having to reconstruct
any dense tensor. The individual factors are contracted, or summed over at
corresponding indices, along an optimized contraction path. These paths are
typically chosen to minimize either memory or runtime (Smith and Gray (2018)).
Quantum tensor contractions consist of highly parallelizable and
arithmetically intensive operations, making them prime candidates for GPU
acceleration (NVIDIA (2019)). A similar approach has proven useful in deep
learning, for instance by expressing the weights of linear (Novikov et al.
(2015)) and convolutional (Kossaifi et al. (2020)) layers of deep nets in
factorized form. However, this approach remains under-explored in quantum
simulation.

Figure 1: TensorLy-Quantum sits atop the hierarchy of TensorLy libraries,
inheriting its extensive tensor methods and providing full Autograd support
for the features (circuits, operations, density operators, and algorithms)
required by quantum simulation.

While a variety of quantum simulators exist (Nguyen et al. (2021); Anis et al.
(2021); pyQuil (2019)), relatively few utilize tensor methods. Among those
that do, integration with essential machine learning tools, such as automatic
differentiation (i.e., PyTorch Autograd), is rare (Gray (2018)). Moreover,
advanced functionalities, like tensor regression, while used in machine
learning (Panagakis et al. (2021)), are so far unused in quantum simulation.

To address such deficits, we have created TensorLy-Quantum, a PyTorch API for
efficient tensor-based simulation of QC and QML protocols. It is a member of
the TensorLy family of libraries and makes use of its extensive features and
optimized implementations (Kossaifi et al. (2019)). As a result, it is the
first quantum library with direct support for tensor decomposition,
regression, and algebra, which have proven fruitful in a myriad of fields and
stand to enrich QC and QML research. Uniquely, TensorLy-Quantum provides
built-in support for Multi-Basis Encoding for MaxCut problems (Patti et al.
(2021a)) and was used to develop Markov Chain Monte Carlo-based QML (Patti et
al. (2021b)). Through these features, TensorLy-Quantum aims to provide
extensive tensor-based quantum simulation capabilities, providing a simple and
flexible API for quantum algorithms research. Moreover, its optimized
functionality and lightweight API would make it an efficient QML backend to
supplement other quantum simulators.

##  2 TensorLy-Quantum

TensorLy-Quantum offers functionalities essential to quantum computation (Fig.
1). The library provides a high-level interface for quantum circuit
simulation, with an API that follows the PyTorch Module structure and offers
end-to-end support for automatic differentiation via PyTorch Autograd. Users
can seamlessly design quantum circuits by combining either pre-defined or
customizable quantum gates, operators, and states. The library also provides
extensive circuit operations, ranging from pre-contraction techniques that
simplify contraction path search, to dynamic partial traces that compactly
evaluate quantum circuit outputs. In addition to specializing in factorized
representations, such as Matrix-Product State (known as tensor-train in the
Machine-Learning literature Oseledets (2011)), TensorLy-Quantum also supports
efficient analysis on quantum density operators, both pure and reduced,
including partial traces and information metrics. TensorLy-Quantum is designed
to bridge the gap between practitioners of quantum and classical machine
learning, providing an intuitive and Pythonic interface that is supplemented
with extensive documentation and 97% unit-test coverage.

TensorLy-Quantum leverages the deep tensorized network capabilities of
TensorLy and TensorLy-Torch, using these to mitigate the significant
computational overhead posed by quantum state space. It is suitable for both
CPUs and GPUs, and that it acquires the GPU acceleration of these parent
libraries. Likewise, while TensorLy-Quantum is PyTorch-based, TensorLy’s
flexible backend structure enables dynamic transitions between many of the
most utilized Python libraries for machine learning and numerical analysis.
Due to its strategic location atop the TensorLy ecosystem (Fig. 1), TensorLy-
Quantum is exceptionally positioned to accelerate and innovate quantum
simulation. In what follows, we illustrate the scalability and speed of our
library, particularly on GPU.

##  3 Performance

Due to its efficiency and scalability, TensorLy-Quantum holds world records
for the number of qubits simulated in a successful quantum optimization
algorithm. These records include both the largest single-qubit implementation
of MaxCut (Patti et al. (2021a)) and a multi-GPU scaleup that used cuQuantum
as a backend for tensor network contraction (NVIDIA (2021a, b)).

(a)

(b)

Figure 2: Logarithmic-scale runtime comparison for (left) full-rank partial
trace from $n$ to $m$ qubits and (right) contraction operations on large-scale
networks of factorized tensors.

We highlight this performance with numerical experiments, showcasing both
operations on the full-rank density operator (a reshaped, matrix-like quantum
state representation) and tensor contraction functionalities of TensorLy-
Quantum. In the density operator experiments, TensorLy-Quantum outperforms the
leading software, QuTip (Johansson et al. (2012)), by two orders of magnitude
on CPU for partial traces of $n=20$ to $m=14,10$ qubit systems and four orders
of magnitude on GPU. TensorLy-Quantum can also complete full-rank calculations
on larger quantum systems than its predecessors through the use of compact
tensor algebra.

Experiments on networks of factorized tensors include a forward pass of the
expectation value $\langle H\rangle$, and the full gradient calculation
$\nabla\langle H\rangle$, where $H$ is a transverse-field Ising model
Hamiltonian of $n=500$ qubits. The circuit ansatze contain $5000$ gates and
the gradient calculation constitutes a full backpropagation with PyTorch
Autograd. GPU acceleration provides a $20$x speedup over the CPU
implementation. Moreover, we emphasize that without factorized tensor methods,
simulations of this size are impossible, as matrix-encoded Hamiltonians
require $\sim 10^{6}$ GB of memory for systems of merely $n\sim 25$ qubits.
Contraction was accomplished with the Opt-Einsum library (Smith and Gray
(2018)) and all GPU experiments were done on an NVIDIA A100 GPU.

##  4 Conclusions

TensorLy-Quantum is an open-source library designed to streamline the workflow
of QC and QML researchers. It provides highly optimized operations for large-
scale and compute heavy quantum circuit simulations, buttressed by extensive
documentation and high-coverage unit-testing. Its API seamlessly integrates
with PyTorch and provides an interface that is amenable to diverse scientific
backgrounds, accessibly packaging state-of-the-art tensor network operations
alongside quantum protocols. Moreover, its lightweight API and optimized
functions would make it an excellent backend for existing quantum simulators.
As TensorLy-Quantum uses the TensorLy libraries as a backend, it offers direct
access to tools unavailable in other quantum APIs, like tensor regression and
decomposition, as well as convenient conversion to the native data structures
of numerous backends. In future releases, we will expand both the efficiency
and scope of TensorLy-Quantum, adding features such as causality-simplified
contraction, quantum state compression, and novel quantum protocols, as well
as support for more factorized representations, using TensorLy’s existing
tensor decomposition.

## References

  * Anis et al. (2021) M.D. Sajid Anis, Héctor Abraham, Adu Offei, and et al.  Qiskit: An open-source framework for quantum computing, 2021.
  * Bridgeman and Chubb (2017) Jacob C Bridgeman and Christopher T Chubb.  Hand-waving and interpretive dance: an introductory course on tensor networks.  _Journal of Physics A: Mathematical and Theoretical_ , 50:223001, 2017.
  * Carleo et al. (2019) Giuseppe Carleo, Ignacio Cirac, Kyle Cranmer, Laurent Daudet, Maria Schuld, Naftali Tishby, Leslie Vogt-Maranto, and Lenka Zdeborová.  Machine learning and the physical sciences.  _Rev. Mod. Phys._ , 91:045002, Dec 2019.  doi: 10.1103/RevModPhys.91.045002.  URL https://link.aps.org/doi/10.1103/RevModPhys.91.045002.
  * Cerezo et al. (2021) Marco Cerezo, Andrew Arrasmith, Ryan Babbush, Simon C Benjamin, Suguru Endo, Keisuke Fujii, Jarrod R McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio, et al.  Variational quantum algorithms.  _Nature Reviews Physics_ , pages 1–20, 2021.
  * Farhi et al. (2014) Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm.  _arXiv preprint arXiv:1411.4028_ , 2014.
  * Gray (2018) Johnnie Gray.  quimb: A python package for quantum information and many-body calculations.  _Journal of Open Source Software_ , 3(29):819, 2018.  doi: 10.21105/joss.00819.  URL https://doi.org/10.21105/joss.00819.
  * Johansson et al. (2012) J.R. Johansson, P.D. Nation, and Franco Nori.  Qutip: An open-source python framework for the dynamics of open quantum systems.  _Computer Physics Communications_ , 183(8):1760–1772, 2012.  ISSN 0010-4655.  doi: https://doi.org/10.1016/j.cpc.2012.02.021.  URL https://www.sciencedirect.com/science/article/pii/S0010465512000835.
  * Kossaifi et al. (2019) Jean Kossaifi, Yannis Panagakis, Anima Anandkumar, and Maja Pantic.  Tensorly: Tensor learning in python.  _Journal of Machine Learning Research_ , 20(26):1–6, 2019.  URL http://jmlr.org/papers/v20/18-277.html.
  * Kossaifi et al. (2020) Jean Kossaifi, Antoine Toisoul, Adrian Bulat, Yannis Panagakis, Timothy M Hospedales, and Maja Pantic.  Factorized higher-order cnns with an application to spatio-temporal emotion estimation.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 6060–6069, 2020.
  * Mehta et al. (2019) Pankaj Mehta, Marin Bukov, Ching-Hao Wang, Alexandre G.R. Day, Clint Richardson, Charles K. Fisher, and David J. Schwab.  A high-bias, low-variance introduction to machine learning for physicists.  _Physics Reports_ , 810:1–124, 2019.  ISSN 0370-1573.  doi: https://doi.org/10.1016/j.physrep.2019.03.001.  URL https://www.sciencedirect.com/science/article/pii/S0370157319300766.  A high-bias, low-variance introduction to Machine Learning for physicists.
  * Montanaro (2016) Ashley Montanaro.  Quantum algorithms: an overview.  _npj Quantum Information_ , page 15023, 2016.
  * Nguyen et al. (2021) Thien Nguyen, Dmitry Lyakh, Eugene Dumitrescu, David Clark, Jeff Larkin, and Alexander McCaskey.  Tensor network quantum virtual machine for simulating quantum circuits at exascale, 2021.
  * Novikov et al. (2015) Alexander Novikov, Dmitry Podoprikhin, Anton Osokin, and Dmitry Vetrov.  Tensorizing neural networks, 2015.
  * NVIDIA (2019) NVIDIA.  Tensor linear algebra on nvidia gpus.  https://developer.nvidia.com/cutensor, 2019.
  * NVIDIA (2021a) NVIDIA.  Nvidia sets world record for quantum computing simulation with cuquantum running on dgx superpod.  https://blogs.nvidia.com/blog/2021/11/09/cuquantum-world-record/, 2021a.
  * NVIDIA (2021b) NVIDIA.  cuquantum: Accelerate quantum information science.  https://developer.nvidia.com/cuquantum-sdk, 2021b.
  * Oseledets (2011) Ivan V Oseledets.  Tensor-train decomposition.  _SIAM Journal on Scientific Computing_ , 33(5):2295–2317, 2011.
  * Panagakis et al. (2021) Yannis Panagakis, Jean Kossaifi, Grigorios G Chrysos, James Oldfield, Mihalis A Nicolaou, Anima Anandkumar, and Stefanos Zafeiriou.  Tensor methods in computer vision and deep learning.  _Proceedings of the IEEE_ , 109(5):863–890, 2021\.
  * Patti et al. (2021a) Taylor L. Patti, Jean Kossaifi, Anima Anandkumar, and Susanne F. Yelin.  Variational quantum optimization with multi-basis encodings, 2021a.
  * Patti et al. (2021b) Taylor L. Patti, Omar Shehab, Khadijeh Najafi, and Susanne F. Yelin.  Markov chain monte-carlo enhanced variational quantum algorithms, 2021b.
  * Peruzzo et al. (2014) Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J Love, Alán Aspuru-Guzik, and Jeremy L O’brien.  A variational eigenvalue solver on a photonic quantum processor.  _Nature communications_ , 5:4213, 2014.
  * pyQuil (2019) pyQuil.  https://pyquil-docs.rigetti.com/en/v2.7.0/qvm.html, 2019.
  * Schuld (2021) Maria Schuld.  Supervised quantum machine learning models are kernel methods, 2021.
  * Smith and Gray (2018) Daniel G. A. Smith and Johnnie Gray.  opt_einsum - a python package for optimizing contraction order for einsum-like expressions.  _Journal of Open Source Software_ , 3(26):753, 2018.  doi: 10.21105/joss.00753.  URL https://doi.org/10.21105/joss.00753.