# Randomized Rank-Revealing UZV Decomposition for Low-Rank Approximation of
Matrices

Maboud F. Kaloorazi, and Rodrigo C. de Lamare

Pontifical Catholic University of Rio de Janeiro, Brazil
E-mail<EMAIL_ADDRESS>Centre for Telecommunications
Studies (CETUC)

###### Abstract

Low-rank matrix approximation plays an increasingly important role in signal
and image processing applications. This paper presents a new rank-revealing
decomposition method called randomized rank-revealing UZV decomposition (RRR-
UZVD). RRR-UZVD is powered by randomization to approximate a low-rank input
matrix. Given a large and dense matrix ${\bf A}\in\mathbb{R}^{m\times n}$
whose numerical rank is $k$, where $k$ is much smaller than $m$ and $n$, RRR-
UZVD constructs an approximation $\hat{\bf A}$ such as $\hat{\bf A}={\bf
UZV}^{T}$, where ${\bf U}$ and ${\bf V}$ have orthonormal columns, the
leading-diagonal block of ${\bf Z}$ reveals the rank of $\bf A$, and its off-
diagonal blocks have small $\ell_{2}$-norms. RRR-UZVD is simple, accurate, and
only requires a few passes through $\bf A$ with an arithmetic cost of $O(mnk)$
floating-point operations. To demonstrate the effectiveness of the proposed
method, we conduct experiments on synthetic data, as well as real data in
applications of image reconstruction and robust principal component analysis.

###### Index Terms:

Low-rank approximation, randomized algorithms, rank-revealing decompositions,
matrix computations, image reconstruction, robust PCA, dimension reduction.

##  I Introduction

Low rank matrix approximation is approximating an input matrix by one of lower
rank. The goal is to compactly represent the matrix with limited loss of
information. Such a representation can provide a significant reduction in
memory requirements as well as computational costs. Given a real large and
dense $m\times n$ matrix ${\bf A}$ with numerical rank $k$ and $m\geq n$, its
singular value decomposition (SVD) [1] is written as follows:

${\bf A}={\bf U}_{A}{\bf\Sigma}_{A}{\bf V}_{A}^{T}=\begin{bmatrix}{{\bf
U}_{k}\quad{\bf U}_{0}}\end{bmatrix}\begin{bmatrix}{\bf\Sigma}_{k}&0\\\
0&{\bf\Sigma}_{0}\end{bmatrix}\begin{bmatrix}{{\bf V}_{k}\quad{\bf
V}_{0}}\end{bmatrix}^{T},$ (1)

where the matrices ${\bf U}_{k}\in\mathbb{R}^{m\times k}$, ${\bf
U}_{0}\in\mathbb{R}^{m\times n-k}$, ${\bf V}_{k}\in\mathbb{R}^{n\times k}$ and
${\bf V}_{0}\in\mathbb{R}^{n\times n-k}$ are orthonormal, and
${\bf\Sigma}_{k}\in\mathbb{R}^{k\times k}$ and
${\bf\Sigma}_{0}\in\mathbb{R}^{n-k\times n-k}$ are diagonal with entries
$\rho_{i}$s as singular values. The SVD constructs the best rank-$k$
approximation ${\bf B}$ to ${\bf A}$ [1], i.e.,

$\displaystyle\ \|{\bf A}-{\bf B}\|_{2}=\rho_{k+1},$ (2) $\displaystyle\
{\|{\bf A}-{\bf
B}\|_{F}}=\Big{(}\rho_{k+1}^{2}+...+\rho_{n}^{2}\Big{)}^{1/2},$

where ${\bf B}={\bf U}_{k}{\bf\Sigma}_{k}{\bf V}_{k}^{T}$, ${\|{\cdot}\|_{2}}$
and ${\|{\cdot}\|_{F}}$ denote the spectral norm and the Frobenius norm,
respectively. In this paper we focus on the matrix $\bf A$ defined above.
Another method to construct low-rank approximations of matrices is column-
pivoted QR decomposition, or rank-revealing QR (RRQR) [2]. The RRQR [2] is a
special QR with column pivoting (QRCP) which reveals the rank, i.e., the gap
in the singular value spectrum, of the input matrix. The RRQR factors the
matrix $\bf A$ as the product:

${\bf A}={\bf Q}{\bf R}{\bf\Pi}^{T}={\bf Q}\begin{bmatrix}{\bf R}_{11}&{\bf
R}_{12}\\\ {\bf 0}&{\bf R}_{22}\end{bmatrix}{\bf\Pi}^{T},$ (3)

where ${\bf Q}\in\mathbb{R}^{m\times n}$ is an orthonormal matrix, ${\bf
R}\in\mathbb{R}^{n\times n}$ is an upper triangular matrix where ${\bf
R}_{11}\in\mathbb{R}^{k\times k}$ is well-conditioned with
$\rho_{\text{min}}({\bf R}_{11})=O(\rho_{k})$, and ${\bf
R}_{22}\in\mathbb{R}^{n-k\times n-k}$ has a sufficiently small
$\ell_{2}$-norm. The matrix ${\bf\Pi}\in\mathbb{R}^{n\times n}$ is a
permutation matrix.

Low-rank matrices arise in many applications such as ranking and collaborative
filtering [3], background modeling [4, 5, 6, 7, 8], image reconstruction [9,
10, 11], system identification [12], Internet Protocol network anomaly
detection [13, 14], sensor and multichannel signal processing [15], and
biometrics [16]. These two widely used methods, however, are computationally
prohibitive for large matrices. In addition, their computations using standard
schemes are challenging to be parallelized on advanced computational
architectures [17, 18, 19]. Recently, low-rank approximation methods based on
randomization have been developed [20, 21, 22, 23, 24, 17, 25, 18, 6].
Randomized methods first transform an input matrix into a lower dimensional
space by means of random matrices, and next apply traditional methods to
further process the reduced-size matrix. These methods have been demonstrated
to be remarkably efficient, accurate, and robust. The performances of
randomized methods are known to be superior to those of the classical methods
in many practical circumstances. The advantage of randomized methods over
their classical counterparts is twofold: i) they are computationally more
efficient, and ii) they readily lend themselves to a parallel implementation
in order to exploit parallel architectures.

  * Structure

The remainder of this paper is organized as follows. In Section II we discuss
prior works, the problem this work is concerned with, and our contributions.
In Section III we describe our proposed method in detail. In Section IV we
present the mathematical analysis of RRR-UZVD. In Section V we present and
discuss our experimental results, and concluding remarks are given in Section
VI.

##  II Prior Works and Problem Statement

Economical variants of SVD and RRQR are partial SVD and truncated RRQR [1].
Both methods can compute an approximation of a matrix $\bf A$ with $O(mnk)$
floating-point operations (flops). Partial SVD is computed by invoking Krylov
subspace methods, such as the Lanczos and Arnoldi algorithms. These methods,
however, have two disadvantages: i) inherently, they are numerically unstable
[1], and ii) they are difficult to parallelize [17, 18], making them
unsuitable to apply on advanced computational platforms. On the other hand,
the major shortcoming in computing RRQR is pivoting strategy in which the
columns of the large input matrix need to be permuted, i.e., RRQR is not pass-
efficient [17]. To address this concern, recently several RRQR algorithms
based on randomization have been proposed in order to carry out the
factorization with minimum communication costs [26, 10, 27]. However, as we
will show, RRQR does not provide highly accurate approximations.

After emerging the works in [20, 21], many randomized algorithms have been
proposed for computing low-rank approximations of matrices. The algorithms in
[22, 28, 29], first sample columns of an input matrix according to a
probability proportional to either their $\ell_{2}$-norm or leverage scores,
leading to a compact represention of the matrix. The submatrix is then used
for further computation using deterministic algorithms such as the SVD and
pivoted QR [1] to construct the final low-rank approximation. Rokhlin et al.
[24] proposed to first apply a random Gaussian embedding matrix to compress
the input matrix in order to obtain an orthnormal basis for the range of the
matrix. Next, the matrix is transformed into a lower-dimensional subspace by
means of the basis. The low-rank approximation is then obtained through
computations by means of an SVD on the reduced-size matrix. Sarlós [23]
proposed a method based on the Johnson-Lindenstrauss (JL) lemma. He showed
that random linear combinations of rows render a good approximation to a low-
rank matrix. The work in [30], built on Sarlós’s idea, computed a low-rank
approximation using subspace embedding. Halko et al. [17] developed an
algorithm based on randomized sampling techniques to compute an approximate
SVD of a matrix. Their method, randomized SVD, for the matrix ${\bf A}$ and
integer $0<k\leq\ell<n$ is computed as described in Alg. 1.

Algorithm 1 Randomized SVD (R-SVD) [17]

0: Matrix $\ {\bf A}\in\mathbb{R}^{m\times n}$, integers $k>0$, $k\leq\ell<n$.

0: A rank-$\ell$ approximation.

1: Draw a random matrix ${\bf\Gamma}\in\mathbb{R}^{n\times\ell}$;

2: Compute the matrix product ${\bf Y}={\bf A}{\bf\Gamma}$;

3: Compute a QR factorization ${\bf Y}={\bf Q}{\bf R}$;

4: Compute the matrix product ${\bf B}={\bf Q}^{T}\bf A$;

5: Compute an SVD ${\bf B}=\bar{\bf U}\bar{\bf\Sigma}\bar{\bf V}^{T}$;

6: Form the low-rank approximation $\hat{\bf A}_{\text{RSVD}}=({\bf Q}\bar{\bf
U})\bar{\bf\Sigma}\bar{\bf V}^{T}$.

The R-SVD computes an SVD of an $n\times\ell$ matrix. For large matrices,
specifically when $m\approx n$, this computation can be burdensome in terms of
the communication cost, i.e., data movement either between different levels of
a memory hierarchy or between processors [31]. This makes the R-SVD unsuitable
on modern computational environments. To address this issue, in this paper we
develop a randomized algorithm that only utilizes the QR factorization to
construct the approximation. The proposed algorithm, due to recently developed
Communication-Avoiding QR (CAQR) algorithms [31], can be organized to exploit
modern architectures, thereby being computed efficiently.

Gu [18] slightly modified the R-SVD algorithm and applied it to improve
subspace iteration methods. The work in [6] proposed an algorithm termed
subspace-orbit randomized SVD (SOR-SVD). SOR-SVD alternately projects the
input matrix onto its column and row space via a random Gaussian matrix. The
matrix is then transformed into a lower dimensional subspace, and a truncated
SVD follows to construct an approximation. The work in [9] proposed a rank-
revealing algorithm based on randomized techniques termed compressed
randomized UTV (CoR-UTV) decomposition. CoR-UTV first compresses the matrix
through approximate orthonormal bases. Next, a QRCP is applied on the
compressed matrix to give the final low-rank approximation.

  * Our Contributions

Driven by recent developments, this paper introduces a rank-revealing
algorithm based on randomized sampling techniques termed randomized rank-
revealing UZV decomposition (RRR-UZVD). RRR-UZVD constructs an approximation
$\hat{\bf A}$ to the matrix ${\bf A}$ such as:

$\hat{\bf A}={\bf UZV}^{T},$ (4)

where ${\bf U}$ and ${\bf V}$ are orthonormal matrices, the leading-diagonal
block of ${\bf Z}$ reveals the rank of $\bf A$, and its off-diagonal blocks
are sufficiently small in magnitude. RRR-UZVD requires a few passes through
$\bf A$, and runs in $O(mnk)$ flops. The main operations of RRR-UZVD consist
of matrix-matrix multiplications and QR factorizations. Recently,
Communication-Avoiding QR algorithms [31] have been developed. These
algorithms can carry out a QR computation with optimal communication costs.
Thus, RRR-UZVD can be optimized for peak machine performance on high
performance computing devices. We establish a theoretical analysis for RRR-
UZVD in which the rank-revealing property of the algorithm is proved (Theorem
15). Moreover, we apply RRR-UZVD to reconstruct a low-rank image. We further
apply RRR-UZVD, as a surrogate to the expensive SVD, to solve the robust
principal component analysis (robust PCA) problem [4, 32, 33] in applications
of background/foreground separation in video sequences as well as removing
shadows and specularities from face images.

##  III Randomized Rank-Revealing UZV Decomposition

A rank-revealing decomposition is any decomposition in which the rank of a
matrix is revealed [2, 34, 35]. The crème de la crème of rank-revealing
decompositions is the SVD [1]. Other deterministic rank-revealing methods
include rank-revealing QR decomposition [2], URV decomposition (URVD) [34],
and ULV decomposition (ULVD) [35]. The drawback of these deterministic
methods, however, includes computational costs, both arithmetic and, more
importantly, communication [31, 36, 10].

The randomized rank-revealing UZV decomposition (RRR-UZVD) furnishes
information on singular values and singular vectors of an input matrix through
randomized schemes. Given ${\bf A}$, and integer $0<k\leq\ell<n$, RRR-UZVD
computes an approximation $\hat{\bf A}$ to ${\bf A}$ such that:

$\hat{\bf A}={\bf U}{\bf Z}{\bf V}^{T}={\bf U}\begin{bmatrix}{\bf Z}_{k}&\bf
G\\\ \bf H&\bf E\end{bmatrix}{\bf V}^{T},$ (5)

where ${\bf U}$ and ${\bf V}$ are orthonormal matrices of size $m\times\ell$
and $n\times\ell$, respectively. The matrix ${\bf Z}_{k}$ is of order $k$ and
nonsingular. The diagonals of ${\bf Z}_{k}$ are estimations to the first $k$
singular values of $\bf A$. The matrices ${\bf G}$, ${\bf H}$, and ${\bf E}$
are of size ${k\times\ell-k}$, ${\ell-k\times k}$, and $\ell-k\times\ell-k$,
respectively all having small $\ell_{2}$-norms. We term the diagonal elements
of ${\bf Z}$, $\bf Z$-$values$ of the matrix ${\bf A}$. The RRR-UZVD is rank-
revealer in that the submatrix ${\bf Z}_{k}$ reveals the numerical rank $k$ of
$\bf A$, and the other submatrices of ${\bf Z}$ have sufficiently small
$\ell_{2}$-norms; see Theorem 15. The definition given for RRR-UZVD is
analogous to those of rank-revealing algorithms in the literature [2, 34, 35,
37, 9]. To compute RRR-UZVD, the input matrix is first transformed into a
lower dimensional space by means of random sampling techniques. Next, the
entries of reduced-size matrix is manipulated. Lastly, the matrix is projected
to the original space. For the matrix $\bf A$, RRR-UZVD is constructed by
taking the following six steps:

  1. 1.

Form a real $n\times\ell$ random test matrix ${\bf\Theta}$.

  2. 2.

Compute the $m\times\ell$ matrix product:

${\bf F}={\bf A}{\bf\Theta}.$ (6)

The matrix $\bf F$ is a projection onto a subspace spanned by $\bf A$’s
columns.

  3. 3.

Compute the $n\times\ell$ matrix product:

${\bf T}={\bf A}^{T}{\bf F}.$ (7)

The matrix $\bf T$ is a projection onto a subspace spanned by $\bf A$’s rows.

  4. 4.

Form QR factorizations of ${\bf F}$ and ${\bf T}$:

${\bf F}={\bf U}{\bf R},\quad\text{and}\quad{\bf T}={\bf V}{\bf S},$ (8)

where ${\bf U}\in\mathbb{R}^{m\times\ell}$ and ${\bf
V}\in\mathbb{R}^{n\times\ell}$ are orthonormal matrices, which approximate the
bases for $\mathcal{R}({\bf A})$ and $\mathcal{R}({\bf A}^{T})$, respectively.
The notation $\mathcal{R}(\cdot)$ refers to the range of a matrix. The
matrices ${\bf R}$, ${\bf S}\in\mathbb{R}^{\ell\times\ell}$ are upper
triangular.

  5. 5.

Compute the ${\ell\times\ell}$ matrix product:

${\bf Z}={\bf U}^{T}{\bf A}{\bf V}.$ (9)

The matrix ${\bf Z}$ is constructed by right and left multiplications of $\bf
A$ through the approximate bases.

  6. 6.

Form the low-rank approximation of $\bf A$ by projecting the reduced matrix
back to the original space:

$\hat{\bf A}={\bf U}{\bf Z}{\bf V}^{T}.$ (10)

The RRR-UZVD needs three passes over the data. But, it is possible to modify
it in order for the algorithm to revisit the matrix ${\bf A}$ only once. For
this purpose, the matrix $\bf Z$ (9) is approximated through the available
matrices: both sides of ${\bf Z}$ are right-multiplied by ${\bf
V}^{T}{\bf\Theta}$, thus ${\bf Z}{\bf V}^{T}{\bf\Theta}={\bf U}^{T}{\bf A}{\bf
V}{\bf V}^{T}{\bf\Theta}$. Considering ${\bf A}\approx{\bf A}{\bf V}{\bf
V}^{T}$, and ${\bf F}={\bf A}{\bf\Theta}$, ${\bf Z}$ is approximated by

${\bf Z}_{\text{approx}}={\bf U}^{T}{\bf F}({\bf
V}^{T}{\bf\Theta})^{\dagger},$ (11)

where $\dagger$ refers to the Moore–Penrose inverse. There exist two concerns
associated with the basic version of RRR-UZVD: i) RRR-UZVD may provide poor
approximate (leading) singular values and singular vectors, compared to those
of the SVD, and ii) the columns of $\bf U$ and $\bf V$ may not be in a
contributing order. Accordingly, the $\bf Z$-$values$ may not be sorted in a
decreasing order. We propose two methods to cope with these issues:

  * 1.

Power iterations. We incoroprate a few steps of the power method [24, 17],
which can substantially improve the performance of RRR-UZVD.

  * 2.

Column permutation. We implement a column pivoting technique as follows:
first, sort the $\bf Z$-$values$ in a non-increasing order, delivering a
permutation matrix ${\bf P}$. Second, right-multiply the matrices $\bf U$ and
$\bf V$ by ${\bf P}$ such as ${\bf U}_{s}={\bf UP}$, and ${\bf V}_{s}=\bf VP$.

The modified RRR-UZVD is presented in Alg. 2. Note that in Alg. 2 (when the
power iteration technique is utilized), a non-updated $\bf T$ must be used to
form the approximation ${\bf Z}_{\text{approx}}$.

We also make use of column permutation technique, a shortcoming of QRCP when
implemented on parallel architectures. However, there are fundamental
differences between QRCP and RRR-UZVD in this regard. In RRR-UZVD:

  * •

The column permutation strategy is implemented on matrices that are much
smaller than the input matrix since $\ell\ll\text{min}\\{m,n\\}$.

  * •

The column permutation technique does not need access columns of the input
matrix itself, i.e., it does not need passes over the data.

Algorithm 2 The RRR-UZVD Algorithm

0: Data matrix $\ {\bf A}\in\mathbb{R}^{m\times n}$, integers $k>0$,
$k\leq\ell<n$, and power iteration factor $q$.

0: A low-rank approximation.

1: Draw a random test matrix ${\bf T}\in\mathbb{R}^{n\times\ell}$;

2: for $j=$ 1 to $q+1$ do

3: Compute the matrix product ${\bf F}={\bf A}{\bf T}$;

4: Compute the matrix product ${\bf T}={\bf A}^{T}{\bf F}$;

5: end for

6: Form QR factorizations ${\bf F}={\bf U}{\bf R}$ and ${\bf T}={\bf V}{\bf
S}$;

7: Compute ${\bf Z}_{\text{approx}}={\bf U}^{T}{\bf F}({\bf V}^{T}{\bf}{\bf
T})^{\dagger}$, (or ${\bf Z}={\bf U}^{T}{\bf A}{\bf V}$);

8: Carry out the column pivoting technique, delivering ${\bf U}_{s}$, ${\bf
V}_{s}$, ${\bf Z}_{\text{approx}}={\bf U}_{s}^{T}{\bf Z}_{1}({\bf
V}_{s}^{T}{\bf}{\bf Z}_{2})^{\dagger}$;

9: Construct the low-rank approximation $\hat{\bf A}={\bf U}_{s}{\bf
Z}_{\text{approx}}{\bf V}_{s}^{T}$.

The key difference between RRR-UZVD with SOR-SVD [6] and CoR-UTV [9] is that
both SOR-SVD and CoR-UTV apply a deterministic algorithm, the former applies
the SVD and latter applies column-pivoted QR, on the reduced matrix in their
procedures to factor an input matrix. However, RRR-UZVD only uses column
permutation and matrix-matrix multiplications on the small matrix. The
deterministic methods employed can be challenging to parallelize when
implemented on high performance computing devices, while matrix-matrix
multiplications can readily be implemented in parallel, as will be explained
later. In addition, SOR-SVD produces a rank-$k$ while RRR-UZVD constructs a
rank-$\ell$ approximation.

##  IV Analysis of RRR-UZV Decomposition

In this section, we show that RRR-UZVD has the rank-revealing property, and
discuss its computational cost.

###  IV-A Rank-Revealing Property

The RRR-UZVD, with partitioned matrices, takes the form:

$\hat{\bf A}={\bf U}{\bf Z}{\bf V}^{T}=\begin{bmatrix}{{\bf U}_{1}\quad{\bf
U}_{2}}\end{bmatrix}\begin{bmatrix}{\bf Z}_{k}&{\bf G}\\\ {\bf H}&{\bf
E}\end{bmatrix}\begin{bmatrix}{{\bf V}_{1}\quad{\bf V}_{2}}\end{bmatrix}^{T},$
(12)

where ${\bf U}_{1}$ and ${\bf U}_{2}$ are orthonormal matrices of size
$m\times k$ and $m\times\ell-k$, respectively. ${\bf Z}_{k}$ is of order $k$
and invertible whose diagonals being approximations of $\ell$ leading singular
values of $\bf A$. The matrices ${\bf V}_{1}$ and ${\bf V}_{2}$ are
orthonormal of size $n\times k$ and ${n\times\ell-k}$, respectively. We
demonstrate that ${\bf Z}_{k}$ reveals the numerical rank of $\bf A$, and
submatrices ${\bf G}\in\mathbb{R}^{k\times\ell-k}$, ${\bf
H}\in\mathbb{R}^{\ell-k\times k}$, ${\bf E}\in\mathbb{R}^{\ell-k\times\ell-k}$
are small in magnitude. We show that if approximate bases for
$\mathcal{R}({\bf A})$ and $\mathcal{R}({\bf A}^{T})$ are obtained through the
procedure described, the middle matrix $\bf Z$ has the rank-revealing
property. The theorem below sets forth the rank-revealing property of the RRR-
UZVD algorithm. This result is new.

###### Theorem 1

Suppose that ${\bf A}$ is a real $m\times n$ matrix, where $m\geq n$, and its
numerical rank is $k$. Further, suppose that the SVD of $\bf A$ is defined as
in (1), and its RRR-UZVD is defined as in (12). Then,

$\rho_{\text{min}}({\bf Z}_{k})=O(\rho_{k}),$ (13) $\|[{\bf H}\quad{\bf
E}]\|_{2}=O(\rho_{k+1}),$ (14) $\|[{\bf G}^{T}\quad{\bf
E}^{T}]^{T}\|_{2}=O(\rho_{k+1}).$ (15)

Proof. The proof is given in the Appendix.

###  IV-B Computational Cost

The basic version of RRR-UZVD in order to construct an approximation to $\bf
A$ has the following costs:

  * •

Forming a random matrix $\bf\Phi$, e.g., standard Gaussian, in Step 1 costs
$O(n\ell)$.

  * •

Computing the matrix product in Step 2 costs $O(mn\ell)$.

  * •

Computing the matrix product in Step 3 costs $O(mn\ell)$.

  * •

Forming the QR factorizations in Step 4 costs $O(m\ell^{2}+n\ell^{2})$.

  * •

Computing the matrix product in Step 5 costs $O(mn\ell+m\ell^{2})$ (if in this
step the matrix $\bf Z$ is estimated by ${\bf Z}_{\text{approx}}$ (11), the
cost would be $O(m\ell^{2}+n\ell^{2}+\ell^{3})$).

The dominant cost of Step 1 through Step 6 occurs when $\bf A$ and ${\bf
A}^{T}$ are multiplied by the corresponding matrices. Hence

$C_{\text{UZV}}=O(mn\ell).$ (16)

The cost of the algorithm increases when the modifications (power iteration
and column permutation techniques) are applied. The column reordering
technique costs $O(m\ell)$. RRR-UZVD needs either three or two passes
(considering $\bf Z$ being approximated by ${\bf Z}_{\text{approx}}$) over
data to decompose $\bf A$. With the power method used (Alg. 2), RRR-UZVD needs
either $(2q+3)$ passes over data with flop count of $(2q+3)C_{\text{UZV}}$, or
$(2q+2)$ passes over data (considering $\bf Z$ being approximated by ${\bf
Z}_{\text{approx}}$) with flop count of $(2q+2)C_{\text{UZV}}$. The sample
size parameter $\ell$ is usually close to the exact numerical rank $k$.

To compute any algorithm one needs to consider two costs: i) arithmetic, that
is, the number of flops, and ii) communication, that is, data flow either
between different levels of a memory hierarchy or between processors [31]. The
communication cost becomes considerably more expensive compared to the
arithmetic when carrying out the computations for a large data matrix stored
externally on modern computational platforms [31, 38]. The RRR-UZVD carries
out several matrix-matrix multiplications which can be readily implemented in
parallel. RRR-UZVD further carries out QR factorizations on two $m\times\ell$
and $n\times\ell$ matrices. The R-SVD [17], however, carries out one QR
factorization on a matrix of size $m\times\ell$ and one SVD on a matrix of
size $n\times\ell$. Recently, Communication-Avoiding QR (CAQR) algorithms [31]
have been developed, which carry out the orthogonal triangularization with
optimal communication costs. However, computing an SVD using standard schemes
are difficult to parallelize [17, 18, 19]. Hence, the operations of RRR-UZVD
can be arranged to produce a low-rank approximation with optimal communication
costs. This is an advantage of RRR-UZVD over R-SVD. Compared to CoR-UTV [9],
RRR-UZVD only makes use of the matrix-vector multiplication to manipulate the
reduced matrix $\bf Z$. While, CoR-UTV employs a deterministic QRCP for
processing the compressed matrix, which may impose considerable communication
cost on the algorithm. This is an advantage of RRR-UZVD over CoR-UTV.

##  V Simulations

In this section, we assess the empirical performance of the RRR-UZVD
algorithm. We first, through numerical examples, show that RRR-UZVD i) is a
rank-revealing algorithm, and ii) furnishes estimates of singular values,
which with remarkable fidelity track the exact singular values of the input
matrix. To make a fair judgment of the behaviour of RRR-UZVD, we include QR
with column pivoting (QRCP), the optimal SVD, R-SVD [17], and CoR-UTV [9] in
our comparison. Next, we treat an image reconstruction problem; we reconstruct
a low-rank image of a differential gear through RRR-UZVD. Lastly, we devise a
robust PCA algorithm by using RRR-UZVD. We experimentally investigate the
efficacy and efficiency of the proposed algorithm on synthetic and real-time
data. The experiments were conducted in MATLAB on a PC with a 3 GHz intel Core
i5-4430 processor and 8 GB of RAM. To compute RRR-UZVD and R-SVD, we generate
a random matrix whose entries are independent, identically distributed
Gaussian random variables of zero mean and unit variance.

###  V-A Rank-Revealing Property and Singular Value Approximation

We first assess the performance of RRR-UZVD on synthetically generated data.
For the first test, we construct a rank-$k$ matrix ${\bf A}$ of order 1000
perturbed with Gaussian noise such that $\bf A={\bf A}_{1}+{\bf A}_{2}$ [39].
${\bf A}_{1}={\bf U}_{{\bf A}_{1}}{\bf\Sigma}_{{\bf A}_{1}}{\bf V}_{{\bf
A}_{1}}^{T}$, where ${\bf U}_{{\bf A}_{1}}$ and ${\bf V}_{{\bf A}_{1}}$ are
random orthonormal matrices, ${\bf\Sigma}_{{\bf A}_{1}}$ only has non-zero
elements $\rho_{i}$s on the diagonal as singular values which decrease
linearly from $1$ to $10^{-9}$, and $\rho_{k+1}=...=\rho_{1000}=0$. The matrix
${\bf A}_{2}$ is Gaussian normalized to have $\ell_{2}$-norm
$\text{gap}\times\rho_{k}$. We set $k=20$, and gap = 0.15.

In the second test, we consider a challenging matrix called the devil’s stairs
[39]; a matrix ${\bf A}$ of order 1000 whose singular value spectrum has
multiple gaps. The singular values are arranged analogous to a descending
staircase. Each step comprises $d=10$ equal singular values.

We factor the test matrices through RRR-UZVD (Alg. 2), the SVD, QRCP, R-SVD,
and CoR-UTV. We compare the results to assess the quality of singular values.
For the randomized methods, we (arbitrarily) set $\ell=2k$, and the power
method factor $q=1$. The results are shown in Fig. 1. We make two
observations: i) for the first test matrix, RRR-UZVD strongly reveals the gap
in $\rho_{20}$ and $\rho_{21}$, and provides estimations to singular values
that show no loss of accuracy compared with the optimal SVD, while QRCP fails
to reveal the rank. ii) For the second test matrix, RRR-UZVD reveals the gaps
between the singular values and, moreover, excellently tracks them. RRR-UZVD
furnishes highly accurate approximations to the singular values of the devil’s
stairs, whereas QRCP fails to reveal the gaps, provides poor estimations to
the singular values, and is not able to track them.

$10$$20$$30$$0$$0.2$$0.4$$0.6$$0.8$$1$Approximation rankMagnitudeNoisy low-
rank$5$$10$$15$$20$$25$$30$$0$$0.2$$0.4$$0.6$$0.8$$1$Approximation rankThe
devil’s stairsSVDQRCPR-SVDCoR-UTVRRR-UZVD Figure 1: Comparison of singular
values of matrices considered computed by the SVD, QRCP, R-SVD [17], CoR-UTV
[9], and proposed RRR-UZVD (Alg. 2).

###  V-B Low-Rank Image Reconstruction

In this experiment, we reconstruct a gray-scale image of a differential gear
with dimension $1280\times 804$ [10] in order to evaluate the quality of
approximation constructed by RRR-UZVD. In our comparison, we consider
truncated QRCP, R-SVD, and the truncated SVD of PROPACK package [40]. The
PROPACK function efficiently computes a specified number of leading singular
pairs, suitable for approximating large matrices.

(a) $Rank=25$ reconstruction

(b) $Rank=85$ reconstruction

Figure 2: Low-rank image reconstruction through the truncated SVD [40],
truncated QRCP, R-SVD [17], and RRR-UZVD with (a) $rank=25$, and (b)
$rank=85$.

The results are displayed in Figs. 2 and 3. Fig. 2 shows the differential gear
images reconstructed using $rank=25$ and $rank=85$, respectively via the
methods considered. From Fig. 2(a), RRR-UZVD and R-SVD with no power
iterations ($q=0$) demonstrate the poorest qualities. A better reconstruction
is shown by the truncated QRCP. However, RRR-UZVD and R-SVD using one step of
power iterations ($q=1$) construct approximations as accurate as the truncated
SVD. In this case, these methods outperform the truncated QRCP. The rank-$85$
approximations in Fig. 2(b) show tiny artifacts in images reconstructed
through RRR-UZVD and R-SVD with $q=0$ as well as truncated QRCP. If one step
of power iterations is incorporated, images reconstructed by truncated SVD,
RRR-UZVD and R-SVD are visually identical.

Fig. 3(a) shows the approximation error in terms of the Frobenius norm against
the approximation rank. The error is calculated as follows:

$\zeta=\|{\bf A}-\hat{\bf A}_{\text{approx}}\|_{F},$ (17)

where $\hat{\bf A}_{\text{approx}}$ is a low-rank approximation produced by
each algorithm. In Fig. 3(b) the execution times of the considered algorithms
are compared, except for truncated QRCP; we have excluded this algorithm
because there does not exist an optimized LAPACK function to compute QRCP with
an assigned rank. The figure shows how the runtime of truncated SVD
considerably grows as the approximation rank increases. While, one step of a
power iteration together with column pivoting technique barely adds to the
runtime of RRR-UZVD. This demonstrates the RRR-UZVD ability to provide
comparable approximations with truncated SVD with much lower computational
cost. R-SVD with $q=0$ appears to be slightly faster than RRR-UZVD with $q=0$
in our experiment. However, as RRR-UZVD’s operations can be arranged to be
carried out with minimum communication costs (subsection IV-B), we expect that
on current and future advanced computational devices the RRR-UZVD algorithm to
be faster than both truncated SVD and R-SVD, where the cost of communication
is a serious bottleneck on performance of any algorithm.

$40$$120$ $200$$280$$10^{1}$$10^{2}$Approximation rank$\|{\bf A}-\hat{\bf
A}_{\text{approx}}\|_{F}$Tr SVDTr QRCPR-SVD q=0R-SVD q=1UZVD q=0UZVD q=1 (a)
Approximation error

$40$$120$ $200$$280$$10^{-2}$$10^{-1}$$10^{0}$Approximation rankTime (seconds)
(b) Computational time

Figure 3: (a) Approximation errors by the methods studied in reconstructing
the low-rank differential gear image. (b) Runtime in seconds for different
methods.

###  V-C Robust PCA with RRR-UZVD

Robust PCA [4, 32, 33] represents an $m\times n$ low-rank matrix ${\bf A}$
whose fraction of entries being grossly perturbed as a linear combination of a
clean low-rank matrix ${\bf B}$ and a sparse matrix of errors ${\bf C}$ such
as ${\bf A=B+C}$. Robust PCA solves the convex program:

$\displaystyle{\text{minimize}_{\bf(B,C)}}\ {\|{\bf B}\|_{*}+\gamma\|{\bf
C}\|_{1}}$ (18) $\displaystyle{\text{subject to}}\ {\bf A}={\bf B}+{\bf C},$

where ${\|\mbox{\bf M}\|_{*}}\triangleq\sum_{i}\rho_{i}(\mbox{\bf M})$ is
defined as the nuclear norm of any matrix M, ${\|\mbox{\bf
M}\|_{1}}\triangleq\sum_{ij}|\mbox{\bf M}_{ij}|$ is defined as the
$\ell_{1}$-norm of M, and $\gamma>0$ is a weighting parameter. The method of
augmented Lagrange multipliers (ALM) [41] is an efficient method to solve
(18), which yields the optimal solution. However, it has a major bottleneck
which is computing an SVD at each iteration for approximating the low-rank
component $\bf B$ of $\bf A$. The work in [42] proposes several techniques
such as predicting the dimension of principal singular space, a continuation
method [43], and a truncated SVD by making use of PROPACK package [40] to
speed up the convergence of the ALM method. The bottleneck of the proposed
algorithm [42], however, is that the employed truncated SVD [40] uses the
lanczos method that is i) numerically unstable, and ii) has poor performance
on modern computing devices because of limited data reuse in its operations
[1, 17, 18]. To address this issue, by retaining the original objective
function [4, 32, 33, 42], we apply RRR-UZVD to solve the robust PCA problem.
We further adopt the continuation technique [43, 42]. We call our proposed
method ALM-UZVD given in Alg. 3. In Alg. 3, for any matrix $\bf M$ with an
RRR-UZVD described in Section III, $\mathcal{Z}_{\nu}({\bf M})$ denotes a UZV
thresholding operator defined as:

$\mathcal{Z}_{\nu}({\bf M})={\bf U}(:,1:s){\bf Z}(1:s,:){\bf V}^{T},$ (19)

where $s$ is the number of diagonal elements of $\bf Z$ larger than $\nu$,
$\mathcal{S}_{\nu}(x)={\text{sgn}(x)\text{max}}(|x|-\nu,0)$ is a shrinkage
operator, and $\gamma$, $\eta_{0}$, ${\bar{\eta}}$, $\tau$, ${\bf Y}_{0}$, and
${\bf C}_{0}$ are initial values.

Algorithm 3 Robust PCA solved via ALM-UZVD

0: Matrix ${\bf A},\gamma,\eta_{0},{\bar{\eta}},\tau,{\bf Y}_{0},{\bf
C}_{0},i=0$;

0: Low-rank plus sparse matrix

1: while the algorithm does not converge do

2: Compute ${\bf B}_{i+1}=\mathcal{Z}_{\eta_{i}^{-1}}({\bf A}-{\bf
C}_{i}+\eta_{i}^{-1}{\bf Y}_{i})$;

3: Compute ${\bf C}_{i+1}=\mathcal{S}_{\gamma\eta_{i}^{-1}}({\bf A}-{\bf
B}_{i+1}+\eta_{i}^{-1}{\bf Y})$;

4: Compute${\bf Y}_{i+1}={\bf Y}_{i}+\eta_{i}({\bf A}-{\bf B}_{i+1}-{\bf
C}_{i+1})$;

5: Update $\eta_{i+1}=\text{max}(\tau\eta_{i},{\bar{\eta}})$;

6: end while

7: return $\bf B^{*}$ and $\bf C^{*}$

####  V-C1 Synthetic Low-Rank and Sparse Matrix Recovery

We construct a matrix $\bf A=B+C$ as a sum of a low-rank matrix ${\bf
B}\in\mathbb{R}^{n\times n}$ and a sparse matrix of errors ${\bf
C}\in\mathbb{R}^{n\times n}$. The matrix ${\bf B}$ is generated as a product
of two standard Gaussian matrices ${\bf W}$, ${\bf Q}\in\mathbb{R}^{n\times
k}$ such as ${\bf B}={\bf W}{\bf Q}^{T}$. The matrix ${\bf C}$ has $c$ non-
zero entries drawn independently from the set $\\{$-100, 100$\\}$. We consider
the rank $k=\text{rank}({\bf B})=0.05\times n$, and $c=\|{\bf
C}\|_{0}=0.05\times n^{2}$, where $\|\cdot\|_{0}$ refers to the
$\ell_{0}$-norm.

We apply the ALM-UZVD and efficiently implemented robust PCA algorithm in
[42], hereafter InexactALM, to $\bf A$ to recover ${\bf B}^{*}$ and ${\bf
C}^{*}$. The numerical results are summarized in Table I. We adopt the initial
values proposed in [42] in our experiments. The algorithms are terminated when
${\|{\bf A}-{\bf B}^{sol}-{\bf C}^{sol}\|_{F}}<10^{-4}{{\|{\bf A}\|_{F}}}$ is
satisfied, where $({\bf B}^{sol},{\bf C}^{sol})$ is the solution pair of
either algorithm. In Table I, $Time$ denotes the computational time in
seconds, $Iter.$ denotes the number of iterations, and $\xi={\|{\bf A}-{\bf
B}^{sol}-{\bf C}^{sol}\|_{F}}/{{\|{\bf A}\|_{F}}}$ is the relative error. RRR-
UZVD needs a preassigned rank $\ell$ to carry out the decomposition. We
arbitrarily set $\ell=2k$, and also $q=2$.

Judging from the results presented in Table I, we make the following
observations: i) in all cases, ALM-UZVD successfully detects the exact rank
$k$ of the matrices , ii) ALM-UZVD provides the optimal solution, while it
needs one more iteration in comparison with InexactALM, and iii) ALM-UZVD
outperforms InexactALM in runtime, with speedups of up to 5 times.

TABLE I: Numerical results of InexactALM [42] and proposed ALM-UZVD (Alg. 3) for synthetic matrix recovery. |  |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---
$n$ |  $k$ |  $c$ |  Algorithms |  $\hat{k}$ |  $\hat{c}$ |  Time |  Iter. |  $\xi$
|  |  |  |  |  |  |  |
1e3 |  50 |  5e4 |  ALM-UZVD 50 5e4 0.6 10 4.2e-5 InexactALM 50 5e4 2.5 9 3.1e-5 |  |  |  |  |
2e3 |  100 |  2e5 |  ALM-UZVD 100 2e5 4.4 10 4.1e-5 InexactALM 100 2e5 17.6 9 4.9e-5 |  |  |  |  |
3e3 |  150 |  45e4 |  ALM-UZVD 150 45e4 10.5 10 5.3e-5 InexactALM 150 45e4 52.9 9 5.2e-5 |  |  |  |  |
|  |  |  |  |  |  |  |

####  V-C2 Background Subtraction in Surveillance Video

In this experiment, we employ ALM-UZVD in order to separate the background and
foreground of a video sequence. We use one surveillance video from [44]. The
video has 200 grayscale frames each with dimension ${256\times 320}$, taken in
a shopping mall. We form an $81920\times 200$ matrix ${\bf A}$ through
stacking individual frames as its columns.

To approximate the low-rank component (background), the RRR-UZVD incorporated
in ALM-UZVD requires a preassigned rank $\ell$. We use the following bound [1]
that relates the rank $k$ of any matrix $\bf M$ with the Frobenius and nuclear
norms in order to determine $\ell$:

$\sqrt{k}\geq\dfrac{\|{\bf M}\|_{*}}{\|{\bf M}\|_{F}}.$ (20)

We assign $\ell=k+p$, where $k$ is the minimum value satisfying (20), and
$p=2$ is an oversampling factor. We again set $q=2$ for RRR-UZVD. Some video
frames with recovered backgrounds and foregrounds are shown in Fig. 4(a). It
is observed that ALM-UZVD recovers the low-rank and sparse components of the
video successfully. Table II presents the results, showing that ALM-UZVD
outperforms InexactALM in terms of runtime.

(a) Shopping mall

(b) Yale B02 face

Figure 4: (a) Original video frames are shown in column 1. Recovered
backgrounds ${\bf B}^{*}$ and foregrounds ${\bf C}^{*}$, by ALM-UZVD are shown
in columns 2 and 3, respectively. (b) Cropped images of a face under different
illuminations are shown in column 1. Clean images and errors associated with
shadows and specularities recovered by ALM-UZVD are shown in columns 2 and 3,
respectively.

####  V-C3 Removing Shadows and Specularities from Face Images

In this experiment, we apply ALM-UZVD in order to remove shadows and
specularities of face Images. We use a face image of dimenstion ${192\times
168}$ with a total of 64 illuminations from the Yale B face database [45]. The
individual images are stacked as the columns to form a $32256\times 64$ matrix
${\bf A}$. We apply ALM-UZVD to $\bf A$. The results are shown in Fig. 4(b),
and Table II. It is observed that the sparse components contain the shadows
and specularities of the face images that have been effectively extracted by
ALM-UZVD. The ALM-UZVD method, we conclude, successfully recovers the face
images taken under distant illumination from the dataset considered nearly two
times faster than InexactALM.

TABLE II: Numerical results of InexactALM [42] and ALM-UZVD (Alg. 3) for real data recovery. |  |
---|---|---
Dataset |  ALM-UZVD Time Iter. $\xi$ |  InexactALM Time Iter. $\xi$
|  |
Shopping mall |  15.8 23 8.1e-5 47.8 23 7.1e-5 |
Yale B02 |  2.0 21 9.1e-5 3.5 21 7.5e-5 |
|  |

##  VI Conclusion

In this paper we proposed RRR-UZVD. The RRR-UZVD method makes use of
randomized sampling to construct an approximation to a low-rank matrix. We
provided theoretical analysis for RRR-UZVD, showing the algorithm is rank-
revealer. Through numerical experiments we showed that RRR-UZVD considerably
outperforms QRCP in revealing the numerical rank of a low-rank matrix, and
provides results with no loss of accuracy compared with those of the optimal
SVD. The performance of RRR-UZVD exceeds that of QRCP in computing low-rank
approximation when the power iteration and column permutation schemes are
employed. In this case, RRR-UZVD furnishes results as accurate as those of the
SVD. Compared to the SVD and QRCP, RRR-UZVD is more efficient in terms of
arithmetic cost. In addition, RRR-UZVD can take advantage of modern computer
architectures better than SVD, QRCP, as well as competing R-SVD and CoR-UTV.
We further applied RRR-UZVD in applications of robust PCA. Our results show
that ALM-UZVD renders the optimal solution, and is considerably faster than
the efficiently implemented InexactALM method.

##  VII Appendix

Proof of Theorem 15.

To prove the first relation in (13), consider ${\bf B}\in\mathbb{R}^{m\times
m}$ and ${\bf C}\in\mathbb{R}^{n\times n}$ to be orthonormal matrices (i.e.,
${\bf B}^{T}{\bf B}={\bf I}_{m}$ and ${\bf C}^{T}{\bf C}={\bf I}_{n}$, where
${\bf I}_{s}$ denotes an $s\times s$ identity matrix). The matrices ${\bf A}$
and $\bf BAC$ have exactly the same singular values [1]. Let a matrix ${\bf
D}$ be a submatrix of $\bf BAC$ such that:

${\bf D}={\bf B}_{1}{\bf A}{\bf C}_{1},$ (21)

where ${\bf B}_{1}\in\mathbb{R}^{k\times m}$ consists of the first $k$ rows of
$\bf B$, and ${\bf C}_{1}\in\mathbb{R}^{n\times k}$ consists of the first $k$
columns of $\bf C$. Suppose that ${\bf D}$ has singular values such as
$\delta_{1}\geq\delta_{2}\geq...\geq\delta_{k}$. We form the following
product:

${\bf DD}^{T}={\bf B}_{1}{\bf MM}^{T}{\bf B}_{1}^{T},$ (22)

where the matrix ${\bf M}={\bf A}{\bf C}_{1}$ is of size $m\times k$. Thus,
${\bf DD}^{T}$ is a $k\times k$ principal submatrix of the $m\times m$
symmetric ${\bf B}{\bf MM}^{T}{\bf B}^{T}$. Suppose that ${\bf M}$ has
singular values such as $\mu_{1}\geq\mu_{2}\geq...\geq\mu_{k}$. Thus, the
eigenvalues of ${\bf MM}^{T}$ satisfy

$\mu_{1}^{2}\geq\mu_{2}^{2}\geq...\geq\mu_{k}^{2}\geq\mu_{k+1}^{2}=...=\mu_{m}^{2}=0.$
(23)

By invoking the formulas that relate the eigenvalues of a principal submatrix
with those of a symmetric matrix [1], we obtain

$\mu_{1}^{2}\geq\delta_{1}^{2}\geq...\geq\mu_{k}^{2}\geq\delta_{k}^{2}.$ (24)

We now form the following product:

${\bf M}^{T}{\bf M}={\bf C}_{1}^{T}{\bf A}^{T}{\bf A}{\bf C}_{1}.$ (25)

The nonzero eigenvalues of ${\bf M}^{T}{\bf M}$ coincide with those of ${\bf
MM}^{T}$. By (25) we observe that ${\bf M}^{T}\bf M$ is a $k\times k$
principal submatrix of the $n\times n$ symmetric ${\bf C}^{T}{\bf A}^{T}{\bf
A}{\bf C}$. Hence

$\rho_{1}^{2}\geq\mu_{1}^{2}\geq...\geq\rho_{k}^{2}\geq\mu_{k}^{2},$ (26)

and accordingly,

$\rho_{1}^{2}\geq\mu_{1}^{2}\geq\delta_{1}^{2}\geq...\geq\rho_{k}^{2}\geq\mu_{k}^{2}\geq\delta_{k}^{2}.$
(27)

Substituting ${\bf B}_{1}$ and ${\bf C}_{1}$ in (21) with ${\bf U}_{1}^{T}$
and ${\bf V}_{1}$ in (12), respectively, results in

$\rho_{k}({\bf Z})\leq\rho_{k}.$ (28)

Note that for any decomposition to be rank-revealing, there must be a well-
defined gap in the singular value spectrum of the input matrix [2, 37].

To prove (14), with an SVD of $\bf A$ defined in (1), let the matrix $\hat{\bf
A}$ formed by RRR-UZVD have an SVD such that $\hat{\bf A}=\hat{\bf
U}\hat{\bf\Sigma}\hat{\bf V}^{T}$. We write $\hat{\bf A}=\hat{\bf
U}_{k}\hat{\bf\Sigma}_{k}\hat{\bf V}_{k}^{T}+\hat{\bf
U}_{0}\hat{\bf\Sigma}_{0}\hat{\bf V}_{0}^{T}$. By the argument given in the
proof of the first bound (13), the following relation must hold

$\rho_{k+1}=\|{\bf A}^{T}{\bf U}_{0}\|_{2}\geq\hat{\rho}_{k+1}=\|\hat{\bf
A}^{T}\hat{\bf U}_{0}\|_{2}.$ (29)

We further have

$\hat{\rho}_{k+1}=\|\hat{\bf A}^{T}\hat{\bf U}_{0}\|_{2}\leq\|\hat{\bf
A}^{T}{\bf U}_{2}\|_{2}.$ (30)

The relation in (30) holds since ${\bf U}_{2}$ provides an approximation to
$\hat{\bf U}_{0}$. However, in practice, combining the power iteration and
column pivoting techniques leads to a ${\bf U}_{2}$ which spans the null space
of $\hat{\bf A}^{T}$. Thus

$\hat{\rho}_{k+1}\approx\|\hat{\bf A}^{T}{\bf U}_{2}\|_{2}.$ (31)

By substituting $\hat{\bf A}$ (12) into (31), it follows that

$\displaystyle\hat{\rho}_{k+1}$ $\displaystyle{\approx\|\begin{bmatrix}{{\bf
V}_{1}\quad{\bf V}_{2}}\end{bmatrix}\begin{bmatrix}{\bf Z}_{k}^{T}&{\bf
H}^{T}\\\ {\bf G}^{T}&{\bf E}^{T}\end{bmatrix}\begin{bmatrix}{\bf
U}_{1}^{T}\\\ {\bf U}_{2}^{T}\end{bmatrix}}{\bf U}_{2}\|_{2}$ (32)
$\displaystyle{\approx\|[{\bf H}\quad{\bf E}]\|_{2}}.$

To prove (15), with an analogous argument as in the proof of (14), we have

$\displaystyle\rho_{k+1}=\|{\bf A}{\bf
V}_{0}\|_{2}\geq\hat{\rho}_{k+1}=\|\hat{\bf A}\hat{\bf V}_{0}\|_{2}.$ (33)

and consequently,

$\hat{\rho}_{k+1}=\|\hat{\bf A}\hat{\bf V}_{0}\|_{2}\leq\|\hat{\bf A}{\bf
V}_{2}\|_{2}.$ (34)

The relation in (34) holds since ${\bf V}_{2}$ is an approximation to
$\hat{\bf V}_{0}$. In practice, combining the power iteration and column
permutation techniques leads to a ${\bf V}_{2}$ which spans the null space of
$\hat{\bf A}$. Thus

$\displaystyle\hat{\rho}_{k+1}\approx\|\hat{\bf A}{\bf V}_{2}\|_{2}.$ (35)

By substituting $\hat{\bf A}$ (12) into (35), it follows that

$\displaystyle\hat{\rho}_{k+1}$ $\displaystyle{\approx\|\begin{bmatrix}{{\bf
U}_{1}\quad{\bf U}_{2}}\end{bmatrix}\begin{bmatrix}{\bf Z}_{k}&{\bf G}\\\ {\bf
H}&{\bf E}\end{bmatrix}\begin{bmatrix}{\bf V}_{1}^{T}\\\ {\bf
V}_{2}^{T}\end{bmatrix}}{\bf V}_{2}\|_{2}$ (36) $\displaystyle{\approx\|[{\bf
G}^{T}\quad{\bf E}^{T}]^{T}\|_{2}}.$

This completes the proof.

The authors would like to thank…

## References

  * [1] G. H. Golub and C. F. van Loan, _Matrix Computations_. 3nd ed.: Johns Hopkins Univ. Press, Baltimore, MD, 1996\.
  * [2] T. F. Chan, “Rank revealing QR factorizations,” _Linear Algebra and its Applications_ , vol. 88–89, pp. 67–82, Apr 1987.
  * [3] N. Srebro, N. Alon, and T. Jaakkola, “Generalization error bounds for collaborative prediction with low-rank matrices,” in _NIPS’04_ , 2004, pp. 5–27.
  * [4] J. Wright, Y. Peng, Y. Ma, A. Ganesh, and S. Rao, “Robust Principal Component Analysis: Exact Recovery of Corrupted Low-Rank Matrices,” _Advances in Neural Information Processing Systems (NIPS)_ , pp. 2080–2088, 2009.
  * [5] T. Bouwmans, F. Porikli, B. Höferlin, and A. Vacavant, _Background Modeling and Foreground Detection for Video Surveillance_. CRC Press, 2014.
  * [6] M. F. Kaloorazi and R. C. de Lamare, “Subspace-Orbit Randomized Decomposition for Low-Rank Matrix Approximations,” _IEEE Trans. Signal Process._ , vol. 66, no. 16, pp. 4409–4424, Aug 2018.
  * [7] T. Bouwmans, S. Javed, H. Zhang, H. Lin, and R. Otazo, “On the Applications of Robust PCA in Image and Video Processing,” _Proceedings of the IEEE_ , vol. 106, no. 8, pp. 1427–1457, Aug 2018.
  * [8] S. Wang, Y. Wang, Y. Chen, P. Pan, Z. Sun, and G. He, “Robust PCA Using Matrix Factorization for Background/Foreground Separation,” _IEEE Access_ , vol. 6, pp. 18 945–18 953, Mar 2018.
  * [9] M. F. Kaloorazi and R. C. de Lamare, “Compressed Randomized UTV Decompositions for Low-Rank Matrix Approximations,” _IEEE Journal of Selected Topics in Signal Process._ , vol. 12, no. 6, pp. 1–15, Dec 2018.
  * [10] J. A. Duersch and M. Gu, “Randomized QR with Column Pivoting,” _SIAM J. Sci. Comput._ , vol. 39, no. 4, pp. C263–C291, 2017.
  * [11] H. Liu, P. Sun, Q. Du, Z. Wu, and Z. Wei, “Hyperspectral image restoration based on low-rank recovery with a local neighborhood weighted spectral- spatial total variation model,” _IEEE Transactions on Geoscience and Remote Sensing_ , pp. 1–14, 2018.
  * [12] M. Fazel, T. K. Pong, D. Sun, and P. Tseng, “Hankel Matrix Rank Minimization with Applications to System Identification and Realization,” _SIAM. J. Matrix Anal. & Appl._, vol. 34, no. 3, pp. 946––977, Apr 2013.
  * [13] M. F. Kaloorazi and R. C. de Lamare, “Anomaly Detection in IP Networks Based on Randomized Subspace Methods,” in _ICASSP 2017_ , Mar 2017.
  * [14] A. M. J. Niyaz Hussain and M. Priscilla, “A Survey on Various Kinds of Anomalies Detection Techniques in the Mobile Adhoc Network Environment,” _Int. J. S. Res. CSE & IT._, vol. 3, no. 3, pp. 1538–1541, 2018.
  * [15] R. de Lamare and R. Sampaio-Neto, “Adaptive Reduced-Rank Processing Based on Joint and Iterative Interpolation, Decimation, and Filtering,” _IEEE Transactions on Signal Processing_ , vol. 57, no. 7, pp. 2503–2514, Jul 2009.
  * [16] B. Victor, K. Bowyer, and S. Sarkar, “An evaluation of face and ear biometrics,” in _ICPR 2002_ , 2002, pp. 429–432.
  * [17] N. Halko, P.-G. Martinsson, and J. Tropp, “Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions,” _SIAM Review_ , vol. 53, no. 2, pp. 217–288, Jun 2011\.
  * [18] M. Gu, “Subspace Iteration Randomization and Singular Value Problems,” _SIAM Journal on Scientific Computing_ , vol. 37, no. 3, pp. A1139–A1173, 2015.
  * [19] H. Anzt, J. Dongarra, M. Gates, K. J, P. Luszczek, S. Tomov, and I. Yamazaki, _Bringing High Performance Computing to Big Data Algorithms_. Handbook of Big Data Technologies. Springer, 2017\.
  * [20] A. Frieze, R. Kannan, and S. Vempala, “Fast monte-carlo algorithms for finding low-rank approximations,” _J. ACM_ , vol. 51, no. 6, pp. 1025–1041, Nov. 2004.
  * [21] S. Goreinov, E. Tyrtyshnikov, and N. Zamarashkin, “A theory of pseudoskeleton approximations,” _Linear Algebra and its Applications_ , vol. 261, pp. 1–21, Aug 1997.
  * [22] P. Drineas, R. Kannan, and M. Mahoney, “Fast Monte Carlo Algorithms for Matrices II: Computing a Low-Rank Approximation to a Matrix,” _SIAM Journal on Computing_ , vol. 36, no. 1, pp. 158––183, Jul 2006.
  * [23] T. Sarlós, “Improved approximation algorithms for large matrices via random projections,” in _FOCS ’06._ , vol. 1, Oct. 2006.
  * [24] V. Rokhlin, A. Szlam, and M. Tygert, “A randomized algorithm for principal component analysis,” _SIAM J. on Matrix Anal. & Applications_, vol. 31, no. 3, pp. 1100–1124, 2009.
  * [25] C. Boutsidis and A. Gittens, “Improved Matrix Algorithms via the Subsampled Randomized Hadamard Transform,” _SIAM J. Matrix Anal. & Appl._, vol. 34, no. 3, pp. 1301–1340, 2013.
  * [26] J. Demmel, L. Grigori, M. Gu, and H. Xiang, “Communication Avoiding Rank Revealing QR Factorization with Column Pivoting,” _SIAM J. Matrix Anal. & Appl._, vol. 36, no. 1, pp. 55––89, 2015.
  * [27] P.-G. Martinsson, G. Quintana Ortí, N. Heavner, and R. van de Geijn, “Householder QR Factorization With Randomization for Column Pivoting (HQRRP),” _SIAM J. Sci. Comput._ , vol. 39, no. 2, pp. C96–C115, 2017.
  * [28] A. Deshpande and S. Vempala, “Adaptive Sampling and Fast Low-Rank Matrix Approximation,” _Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques_ , vol. 4110, pp. 292–303, 2006.
  * [29] M. Rudelson and R. Vershynin, “Sampling from large matrices: An approach through geometric functional analysis,” _J. ACM_ , vol. 54, no. 4, Jul. 2007\.
  * [30] K. Clarkson and D. Woodruff, “Low-rank approximation and regression in input sparsity time,” _J. ACM_ , vol. 63, no. 6, pp. 54:1–54:45, Jan. 2017.
  * [31] J. Demmel, L. Grigori, M. Hoemmen, and J. Langou, “Communication-optimal Parallel and Sequential QR and LU Factorizations,” _SIAM Journal on Scientific Computing_ , vol. 34, no. 1, pp. A206––A239, 2012.
  * [32] V. Chandrasekaran, S. Sanghavi, P. Parrilo, and A. Willsky, “Rank-Sparsity Incoherence for Matrix Decomposition,” _SIAM J Opt_ , vol. 21, no. 2, pp. 572–596, 2009.
  * [33] E. J. Candès, X. Li, Y. Ma, and J. Wright, “Robust principal component analysis?” _J. ACM_ , vol. 58, no. 3, pp. 1–37, May 2011.
  * [34] G. Stewart, “An updating algorithm for subspace tracking,” _IEEE Trans. Signal Process._ , vol. 40, no. 6, pp. 1535–1541, Jun 1992.
  * [35] G. W. Stewart, “Updating a Rank-Revealing ULV Decomposition,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 14, no. 2, pp. 494––499, 1993.
  * [36] B. C. Gunter and R. A. Van De Geijn, “Parallel out-of-core computation and updating of the qr factorization,” _ACM Trans. Math. Softw._ , vol. 31, no. 1, pp. 60–78, Mar. 2005.
  * [37] S. Chandrasekaran and I. Ipsen, “On rank-revealing QR factorizations,” _SIAM. J. Matrix Anal. & Appl._, vol. 15, no. 2, pp. 946––977, 1994.
  * [38] J. Dongarra, S. Tomov, P. Luszczek, J. Kurzak, M. Gates, I. Yamazaki, H. Anzt, A. Haidar, and A. Abdelfattah, “With Extreme Computing, the Rules Have Changed,” _Computing in Science & Engineering_, vol. 19, no. 3, pp. 52–62, May 2017.
  * [39] G. W. Stewart, “The QLP Approximation to the Singular Value Decomposition,” _SIAM Journal on Scientific Computing_ , vol. 20, no. 4, pp. 1336–1348, 1999\.
  * [40] R. M. Larsen, _Efficient algorithms for helioseismic inversion_ , PhD Thesis, University of Aarhus, Denmark (1998).
  * [41] X. Yuan and J. Yang, “Sparse and low-rank matrix decomposition via alternating direction methods,” _Pacific Journal of Opt._ , vol. 9, no. 1, pp. 167–180, 2009.
  * [42] Z. Lin, R. Liu, and Z. Su, “Linearized Alternating Direction Method with Adaptive Penalty for Low-Rank Representation,” in _NIPS_ , no. 1, 2011, pp. 1–9.
  * [43] K. C. Toh and S. Yun, “An accelerated proximal gradient algorithm for nuclear norm regularized linear least squares problems,” _Pacific J. of Opt._ , vol. 6, no. 3, pp. 615–640, 2010.
  * [44] L. Li, W. Huang, I.-H. Gu, and Q. Tian, “Statistical Modeling of Complex Backgrounds for Foreground Object Detection,” _IEEE Trans. Image Process._ , vol. 13, no. 11, pp. 1459–1472, Nov 2004.
  * [45] A. S. Georghiades, P. N. Belhumeur, and D. J. Kriegman, “From few to many: illumination cone models for face recognition under variable lighting and pose,” _IEEE Trans. on Pattern Anal Mach Intell_ , vol. 23, no. 6, pp. 643–660, 2001.


Pairing the $Q_{\Delta}$ functions the other way around, we get the same bound
but with the norm
$\|(z_{it_{2}},z_{jt_{2}})-(z^{\prime\prime}_{it_{1}},z^{\prime\prime}_{jt_{1}})\|$,

$\displaystyle
D^{t_{2}}_{(z^{\prime\prime}_{it_{2}},z^{\prime\prime}_{jt_{2}})}D^{t_{1}}_{(z^{\prime}_{it_{2}},z^{\prime}_{jt_{2}})}B_{\Delta}({\bf
z}_{i},{\bf z}_{j})$
$\displaystyle\leq\frac{\|(z_{it_{2}},z_{jt_{2}})-(z^{\prime\prime}_{it_{1}},z^{\prime\prime}_{jt_{1}})\|}{\Delta}\left(\frac{2}{m^{2}}+\frac{6n}{m^{2}}+\frac{4n(n-1)}{m^{2}}\right).$

Averaging the above two bounds and using that

$\|(z_{it_{1}},z_{jt_{1}})-(z^{\prime}_{it_{1}},z^{\prime}_{jt_{1}})\|+\|(z_{it_{2}},z_{jt_{2}})-(z^{\prime\prime}_{it_{1}},z^{\prime\prime}_{jt_{1}})\|\leq\sqrt{2}\|{\bf
z}-{\bf z}^{\prime\prime}\|,$

by the definition of $J(B_{\Delta})$ we get

$J(B_{\Delta})~{}\leq~{}\frac{4n^{2}+2n+2}{\sqrt{2}\Delta
m}\leq\frac{(2n+1)^{2}}{\sqrt{2}\Delta m}.$ (47)

#### Bounding
$\mathbb{E}_{\tilde{S}_{i},\tilde{S}_{j}}[R(\mathcal{F}^{\alpha}(\tilde{S}_{i}\cup\tilde{S}_{j})))]$.

Similarly to (31),

$\mathbb{E}_{\tilde{S}_{i},\tilde{S}_{j}}[R(\mathcal{F}^{\alpha}(\tilde{S}_{i}\cup\tilde{S}_{j})))]~{}\leq~{}2\alpha
B\sqrt{m(q\log(2)+\log(p))}$ (48)

as $f(\tilde{S}_{i}\cup\tilde{S}_{j})$ is a collection of $2m$ $p$-dimensional
vectors with norm at most $B$.

Substituting the bounds (44), (45), (47) and (48) into (38), we conclude that
with probability at least $1-\delta^{\prime}$ over the selection of
$\tilde{S}_{i}$ and $\tilde{S}_{j}$, for any $f\in\mathcal{F}^{\prime}$ with
$\sup_{f\in\mathcal{F}^{\prime}}\mathcal{C}(f)\leq\alpha$,

$\displaystyle\ell_{\Delta}(f;\tilde{P}_{i},\tilde{P}_{j}){}$
$\displaystyle\leq~{}4\left(1+\frac{n}{m-n}\right)\cdot\ell_{\Delta}(f;\tilde{S}_{i},\tilde{S}_{j})$
(49) $\displaystyle\qquad+16\alpha
B\frac{\sqrt{m}(2n+4)^{2}}{(m-n)\Delta}\sqrt{\pi(q\log(2)+\log(p))\log(2m)}$
$\displaystyle\qquad+32\alpha
B\frac{\sqrt{m}}{(m-n)\Delta}\sqrt{2\log({\textstyle\frac{1}{\delta^{\prime}}})}$

To finish the proof, we extend the above bound to all $\Delta>0$ and all
$f\in\mathcal{F}$ with $\alpha$ replaced with the complexity $\mathcal{C}(f)$.
We proceed exactly as in the proof of Lemma 4: We apply (49) to the function
class
$\mathcal{F}^{\prime}=\mathcal{F}^{\alpha}:=\\{f\in\mathcal{F}:\mathcal{C}(f)\leq\alpha\\}$
with positive integers $\alpha$ and for an exponential grid for $\Delta$
defined as $\Delta_{t}=2^{-t}B$ for any integer $t$ and
$\delta^{\prime}=\delta_{t,\alpha}$ where
$\delta_{t,\alpha}:={\textstyle\frac{\delta}{3\alpha(\alpha+1)|t|(|t|+1)}}$
for all $t\neq 0$ and
$\delta_{0,\alpha}:={\textstyle\frac{\delta}{3\alpha(\alpha+1)}}$. Note that
by definition,
$\sum_{t=-\infty}^{\infty}\sum_{\alpha=1}^{\infty}\delta_{t,\alpha}=\delta$.
Therefore, by the union bound, with probability at least $1-\delta$ over the
selection of $\tilde{S}_{i}$ and $\tilde{S}_{j}$, for any
$f\in\mathcal{F}^{\alpha}$, (49) holds for $\Delta=\Delta_{t}$ and
$\delta^{\prime}=\delta_{t,\alpha}$ simultaneously for all integers $t$ and
positive integers $\alpha$.

The final bound of the lemma follows for any $\Delta>0$ and $f\in\mathcal{F}$
by instantiating this bound with $\alpha=\lceil\mathcal{C}(f)\rceil$ and $t$
such that $\Delta\leq\Delta_{t}\leq 2\Delta$ (which always exists by the
definition of the grid $\\{\Delta_{t}\\}$), and then applying the bound
$\ell_{\Delta}(f;\tilde{P}_{i},\tilde{P}_{j})\leq\ell_{\Delta_{t}}(f;\tilde{P}_{i},\tilde{P}_{j})\leq\ell_{2\Delta}(f;\tilde{P}_{i},\tilde{P}_{j})$
which follows from the monotonicity of $\ell_{\Delta}$, bounding
$1/\Delta_{t}$ by $1/\Delta$, and bounding $\delta_{t,\alpha}$ using that
$|t|\leq|\log({\textstyle\frac{\Delta}{B}})|+1$ (which again holds by the
definition of $\Delta_{t}$).

##  C Proofs of Lemmas 6-8

See 6

Proof  For $c\in\\{i,j\\}$, let $\hat{S}_{c}\sim Q^{n}_{c}$ and $x_{c}\sim
Q_{c}$, $u_{c}=f(x_{c})$, $\mu_{c}=\mu_{f}(Q_{c})$,
$\hat{\mu}_{c}=\mu_{f}(\hat{S}_{c})$ and
$\sigma^{2}_{c}=\mathrm{Var}(\hat{u}_{c})$ and let $s$ be a variable that we
can choose to be $p$ if $\\{f\circ Q_{i},f\circ Q_{j}\\}$ are spherically
symmetric distributions and $s=1$ otherwise. Let
$\Lambda_{ij}:=\|\mu_{i}-\mu_{j}\|$ denote the considered centers’ true
distance, which sets the overall scale, and
$\hat{\Lambda}_{ij}:=\|\hat{\mu}_{i}-\hat{\mu}_{j}\|$ its estimate. Let
$\gamma_{1},\gamma_{2},\gamma_{3},\gamma_{4}\in(0,1)$ be constants that
satisfy $\gamma_{4}-2\gamma_{3}-4\gamma_{1}-2\gamma^{2}_{2}/s>0$ (to be
selected by the end of the proof). Moreover, let
$\alpha_{ij}:=\Lambda_{ij}^{2}\gamma_{1}/\Delta$, where
$\Delta\leq\Lambda_{ij}\gamma_{2}/\sqrt{s}$. Since
$\mathbb{I}[r<\Delta]\leq\ell_{\Delta}(-r)$, we have

$\displaystyle\ell_{\Delta}(f;Q_{i},Q_{j})~{}$
$\displaystyle=~{}\mathbb{E}_{\hat{S}_{i},\hat{S}_{j}}\mathbb{E}_{x_{i}}[\ell_{\Delta}(\|u_{i}-\hat{\mu}_{i}\|\leq\|u_{i}-\hat{\mu}_{j}\|)]$
$\displaystyle\leq~{}\Pr\left[\|u_{i}-\hat{\mu}_{j}\|\leq\|u_{i}-\hat{\mu}_{i}\|+\Delta\right].$

$\hat{e}_{ij}$$\perp\\!\\!\hat{e}_{ij}$$B$$\alpha_{ij}$ErrU${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}B}\\!\cap\\!{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\textnormal{ErrU}}$$H$$\hat{\mu}_{i}$$\mu_{i}$$\mu_{j}$$\hat{\mu}_{j}$$u_{i}$$\hat{\Lambda}_{ij}$$\hat{\beta}_{ij}$
Figure 6: A geometric illustration of the various components in the proof of
Lemma 6. The blue section denotes the ball $B$ of radius $\alpha_{ij}$ and the
red section denotes the set ErrU.

The proof idea is graphically illustrated in Figure 6: Consider the error
event
$\text{Err}:=[\|u_{i}-\hat{\mu}_{j}\|\leq\|u_{i}-\hat{\mu}_{i}\|+\Delta]$. For
$\Delta=0$, Err is a $p$-dimensional half-space orthogonal to the unit vector
$\hat{e}_{ij}:=(\hat{\mu}_{j}-\hat{\mu}_{i})/\hat{\Lambda}_{ij}$, and to
determine its probability, we can simply project all quantities onto the line
$\mathbb{R}\cdot\hat{e}_{ij}$. This makes the problem one-dimensional and
improves the bound by a factor of $p^{-1}$ in the spherically symmetric case.

However, one complication is that $\hat{e}_{ij}$ itself is random, so in order
to achieve a scale of $p^{-1}$, one has to carefully avoid using triangle
inequalities involving $u_{i}$, which complicates the analysis, and also gives
an extra $1/n$ term.

Another complicated arises from assuming $\Delta>0$: Consider the error set
$\text{ErrU}(\hat{\mu}_{i},\hat{\mu}_{i}):=\\{u_{i}:\|u_{i}-\hat{\mu}_{j}\|\leq\|u_{i}-\hat{\mu}_{i}\|+\Delta\\}$.
Its boundary $\partial\text{ErrU}$ is now (one sheet of an) elliptic
hyperboloid, and direct projection onto $\hat{e}_{ij}$ does not work anymore.
If we confine $u_{i}$ to the hyperball
$B:=\\{u_{i}:\|u_{i}-\hat{\mu}_{i}\|\leq\alpha_{ij}\\}$, then $B$ cuts ErrU
into a finite moon-shaped set, and the projection method works again, while
the probability of $u_{i}\not\in B$ scales with $1/p$.

We derive a sequence of bounds (A)-(H) which we chain all together at the end.

#### (A)

As a first step we determine the halfspace (given $\hat{\mu}_{i}$ and
$\hat{\mu}_{j}$) that covers the moon $B\cap\text{ErrU}$. Assume that
$\|u_{i}-\hat{\mu}_{j}\|\leq\|u_{i}-\hat{\mu}_{i}\|+\Delta$ and
$\|u_{i}-\hat{\mu}_{i}\|\leq\alpha_{ij}$. Then,

$\displaystyle\|u_{i}-\hat{\mu}_{i}\|^{2}+2\alpha_{ij}\Delta+\Delta^{2}~{}$
$\displaystyle\geq~{}\|u_{i}-\hat{\mu}_{i}\|^{2}+2\|u_{i}-\hat{\mu}_{i}\|\cdot\Delta+\Delta^{2}$
$\displaystyle=~{}(\|u_{i}-\hat{\mu}_{i}\|+\Delta)^{2}$
$\displaystyle\geq~{}\|u_{i}-\hat{\mu}_{j}\|^{2}$
$\displaystyle=~{}\|(u_{i}-\hat{\mu}_{i})-(\hat{\mu}_{j}-\hat{\mu}_{i})\|^{2}$
$\displaystyle=~{}\|u_{i}-\hat{\mu}_{i}\|^{2}-2(u_{i}-\hat{\mu}_{i})^{\top}(\hat{\mu}_{j}-\hat{\mu}_{i})+\hat{\Lambda}_{ij}^{2}.$

Therefore, we have
$2(u_{i}-\hat{\mu}_{i})^{\top}(\hat{\mu}_{j}-\hat{\mu}_{i})~{}\geq~{}\hat{\Lambda}_{ij}^{2}-2\alpha_{ij}\Delta-\Delta^{2}$,
which can also be written as follows

$\displaystyle(u_{i}-\hat{\mu}_{i})^{\top}\hat{e}_{ij}~{}\geq~{}\hat{\beta}_{ij}~{}:=~{}{\textstyle\frac{1}{2}}\hat{\Lambda}_{ij}-\frac{2\alpha_{ij}\Delta+\Delta^{2}}{2\hat{\Lambda}_{ij}}.$

This is the halfspace $H$ that covers $B\cap\text{ErrU}$ most tightly.

#### (B)

We can use this to bound the probability of the event Err in terms of
$\alpha_{ij}$ and $\Delta$ by conditioning:

$\displaystyle\Pr[\text{Err}]~{}$
$\displaystyle=~{}\Pr[\|u_{i}-\hat{\mu}_{j}\|\leq\|u_{i}-\hat{\mu}_{i}\|+\Delta]$
$\displaystyle=~{}\Pr[\|u_{i}-\hat{\mu}_{j}\|\leq\|u_{i}-\hat{\mu}_{i}\|+\Delta~{}\land~{}\|u_{i}-\hat{\mu}_{i}\|\leq\alpha_{ij}]$
$\displaystyle~{}~{}~{}+\Pr[\|u_{i}-\hat{\mu}_{j}\|\leq\|u_{i}-\hat{\mu}_{i}\|+\Delta~{}\land~{}\|u_{i}-\hat{\mu}_{i}\|>\alpha_{ij}]$
$\displaystyle\smash{\stackrel{{\scriptstyle(A)}}{{\leq}}}~{}\Pr[(u_{i}-\hat{\mu}_{i})^{\top}\hat{e}_{ij}\geq\hat{\beta}_{ij}]~{}+~{}\Pr[\|u_{i}-\hat{\mu}_{i}\|>\alpha_{ij}]~{}=:~{}(C)~{}+~{}(F).$

#### (C)

Next we upper bound the term
$\Pr[(u_{i}-\hat{\mu}_{i})^{\top}\hat{e}_{ij}\geq\hat{\beta}_{ij}]$. Namely,

$\displaystyle\Pr[(u_{i}-\hat{\mu}_{i})^{\top}\hat{e}_{ij}\geq\hat{\beta}_{ij}]$
$\displaystyle=~{}\Pr[(u_{i}-\mu_{i})^{\top}\hat{e}_{ij}+(\mu_{i}-\hat{\mu}_{i})^{\top}\hat{e}_{ij}\geq\hat{\beta}_{ij}]$
$\displaystyle\leq~{}\Pr[(u_{i}-\mu_{i})^{\top}\hat{e}_{ij}+(\mu_{i}-\hat{\mu}_{i})^{\top}\hat{e}_{ij}\geq\gamma_{3}\Lambda_{ij}~{}\lor~{}\hat{\beta}_{ij}\leq\gamma_{3}\Lambda_{ij}]$
$\displaystyle\leq~{}\Pr[(u_{i}-\mu_{i})^{\top}\hat{e}_{ij}+(\mu_{i}-\hat{\mu}_{i})^{\top}\hat{e}_{ij}\geq\gamma_{3}\Lambda_{ij}]~{}+~{}\Pr[\hat{\beta}_{ij}\leq\gamma_{3}\Lambda_{ij}]~{}=:~{}(D)~{}+~{}(H),$

for $\gamma_{3}\in(0,1)$.

#### (D)

By applying the union bound and use that a projection onto $\hat{e}_{ij}$
decreases the length of the vector $\mu_{i}-\hat{\mu}_{i}$, we get

$\displaystyle(D)~{}$
$\displaystyle=~{}\Pr\left[(u_{i}-\mu_{i})^{\top}\hat{e}_{ij}+(\mu_{i}-\hat{\mu}_{i})^{\top}\hat{e}_{ij}\geq\gamma_{3}\Lambda_{ij}\right]$
$\displaystyle\leq~{}\Pr\left[(u_{i}-\mu_{i})^{\top}\hat{e}_{ij}\geq\gamma_{3}\Lambda_{ij}/2~{}\lor~{}(\mu_{i}-\hat{\mu}_{i})^{\top}\hat{e}_{ij}\geq\gamma_{3}\Lambda_{ij}/2\right]$
$\displaystyle\leq~{}\Pr\left[(u_{i}-\mu_{i})^{\top}\hat{e}_{ij}\geq\gamma_{3}\Lambda_{ij}/2\right]+\Pr\left[\|\mu_{i}-\hat{\mu}_{i}\|\geq\gamma_{3}\Lambda_{ij}/2\right]~{}=:~{}(E)~{}+~{}(G).$

#### (E)

Next we analyze each term separately. For a fixed pair $\hat{\mu}_{i}$ and
$\hat{\mu}_{j}$, the variance of the random variable
$(u_{i}-\mu_{i})^{\top}\hat{e}_{ij}$ is at most $\mathrm{Var}(u_{i})/s$, where
we have exploited the projection in the spherically symmetric case. Therefore,
by Markov’s inequality,

$\displaystyle\Pr\left[(u_{i}-\mu_{i})^{\top}\hat{e}_{ij}\geq\Lambda_{ij}\gamma_{3}/2\right]~{}$
$\displaystyle\leq~{}\mathbb{E}_{\hat{\mu}_{i},\hat{\mu}_{j}}\left[\textstyle{\Pr_{u_{i}}}\left[(u_{i}-\mu_{i})^{\top}\hat{e}_{ij}\geq\Lambda_{ij}\gamma_{3}/2\right]\right]$
$\displaystyle\leq~{}\mathbb{E}_{\hat{\mu}_{i},\hat{\mu}_{j}}\left[\frac{4\gamma_{3}^{-2}\mathrm{Var}(u_{i})}{s\cdot\Lambda_{ij}^{2}}\right]~{}=~{}\frac{4\mathrm{Var}(u_{i})}{s\cdot\gamma^{2}_{3}\Lambda_{ij}^{2}}~{}.$

#### (F)

Again, by Markov’s inequality, and using
$\Delta\leq\Lambda_{ij}\gamma_{2}/\sqrt{s}$, hence
$\alpha_{ij}=\Lambda_{ij}^{2}\gamma_{1}/\Delta\geq\Lambda_{ij}\sqrt{s}(\gamma_{1}/\gamma_{2})$,
we have

$\displaystyle\Pr\left[\|u_{i}-\hat{\mu}_{i}\|\geq\alpha_{ij}\right]~{}$
$\displaystyle\leq~{}\Pr\left[\|\hat{\mu}_{i}-\mu_{i}\|\geq\alpha_{ij}/2\right]+\Pr\left[\|u_{i}-\mu_{i}\|\geq\alpha_{ij}/2\right]$
$\displaystyle\leq~{}\frac{4\mathrm{Var}(\hat{\mu}_{i})}{\alpha^{2}_{ij}}+\frac{4\mathrm{Var}(u_{i})}{\alpha^{2}_{ij}}$
$\displaystyle\leq~{}\frac{4\gamma^{2}_{2}\mathrm{Var}(\hat{\mu}_{i})}{\gamma^{2}_{1}\cdot
s\cdot\Lambda_{ij}^{2}}+\frac{4\gamma^{2}_{2}\mathrm{Var}(u_{i})}{\gamma^{2}_{1}\cdot
s\cdot\Lambda_{ij}^{2}}~{},$

#### (G)

and

$\Pr\left[\|\hat{\mu}_{i}-\mu_{i}\|\geq\gamma_{3}\Lambda_{ij}/2\right]\leq\frac{4\mathrm{Var}(\hat{\mu}_{i})}{\gamma^{2}_{3}\Lambda_{ij}^{2}}=\frac{4\mathrm{Var}(u_{i})}{n\gamma^{2}_{3}\Lambda_{ij}^{2}}~{}.$

#### (H)

Next, by conditioning,

$\displaystyle\Pr[\hat{\beta}_{ij}\leq\gamma_{3}\Lambda_{ij}]~{}$
$\displaystyle=~{}\Pr[\hat{\beta}_{ij}\leq\gamma_{3}\Lambda_{ij}~{}\land~{}\hat{\Lambda}_{ij}\leq\gamma_{4}\Lambda_{ij}]+\Pr[\hat{\beta}_{ij}\leq\gamma_{3}\Lambda_{ij}~{}\land~{}\hat{\Lambda}_{ij}>\gamma_{4}\Lambda_{ij}]$
$\displaystyle\leq~{}\Pr[\hat{\Lambda}_{ij}\leq\gamma_{4}\Lambda_{ij}]+\Pr[\hat{\beta}_{ij}\leq\gamma_{3}\Lambda_{ij}~{}\land~{}\hat{\Lambda}_{ij}>\gamma_{4}\Lambda_{ij}]$
$\displaystyle\leq~{}\Pr[\|\mu_{j}-\hat{\mu}_{j}\|+\|\mu_{i}-\hat{\mu}_{i}\|\geq(1-\gamma_{4})\Lambda_{ij}]+\Pr[\hat{\beta}_{ij}\leq\gamma_{3}\Lambda_{ij}~{}\land~{}\hat{\Lambda}_{ij}>\gamma_{4}\Lambda_{ij}],$

where the last inequality follows from solving the triangle inequality

$\displaystyle\Lambda_{ij}\leq\|\mu_{i}-\hat{\mu}_{i}\|+\hat{\Lambda}_{ij}+\|\hat{\mu}_{j}-\mu_{j}\|.$

We can bound the first term by

$\displaystyle\Pr[\|\mu_{j}-\hat{\mu}_{j}\|+\|\mu_{i}-\hat{\mu}_{i}\|\geq(1-\gamma_{4})\Lambda_{ij}]$
$\displaystyle\leq~{}\Pr[\|\mu_{j}-\hat{\mu}_{j}\|\geq\gamma_{4}\Lambda_{ij}/2]+\Pr[\|\mu_{i}-\hat{\mu}_{i}\|\geq(1-\gamma_{4})\Lambda_{ij}/2]$
$\displaystyle\leq~{}\frac{4(\mathrm{Var}(\hat{\mu}_{j})+\mathrm{Var}(\hat{\mu}_{i}))}{(1-\gamma_{4})^{2}\Lambda_{ij}^{2}}.$

We notice that if $\hat{\beta}_{ij}\leq\gamma_{3}\Lambda_{ij}$ and
$\hat{\Lambda}_{ij}>\gamma_{4}\Lambda_{ij}$, then

$\displaystyle\gamma_{3}\Lambda_{ij}~{}\geq~{}\hat{\beta}_{ij}~{}$
$\displaystyle>~{}\frac{\gamma_{4}\Lambda_{ij}}{2}-\frac{2\alpha_{ij}\Delta+\Delta^{2}}{\Lambda_{ij}}$
$\displaystyle\geq~{}\frac{\gamma_{4}\Lambda_{ij}}{2}-\frac{2\gamma_{1}\Lambda_{ij}^{2}+\gamma^{2}_{2}\Lambda_{ij}^{2}/p}{\Lambda_{ij}}$
$\displaystyle\geq~{}\Lambda_{ij}(\gamma_{4}/2-2\gamma_{1}-\gamma^{2}_{2}/p),$

which is a contradiction to
$\gamma_{4}-2\gamma_{3}-4\gamma_{1}-2\gamma^{2}_{2}/s>0$ for appropriate
choices of $\gamma_{1},\dots,\gamma_{4}\in(0,1)$. Hence

$\displaystyle\Pr[\hat{\beta}_{ij}\leq\gamma_{3}\Lambda_{ij}]~{}$
$\displaystyle\leq~{}\frac{4(\mathrm{Var}(\hat{\mu}_{j})+\mathrm{Var}(\hat{\mu}_{i}))}{(1-\gamma_{4})^{2}\Lambda_{ij}^{2}}$

#### Combining everything together

$(B)\leq(C)+(F)\leq(D)+(H)+(F)\leq(G)+(E)+(F)+(H)$ and using
$\mathrm{Var}(\hat{\mu}_{i})=\mathrm{Var}(u_{i})/n$ we finally obtain

$\displaystyle\Pr[\text{Err}]~{}$
$\displaystyle\leq~{}\frac{4\mathrm{Var}(u_{i})}{n\gamma^{2}_{3}\Lambda_{ij}^{2}}+\frac{4\mathrm{Var}(u_{i})}{s\cdot\gamma^{2}_{3}\Lambda_{ij}^{2}}+\frac{4\gamma^{2}_{2}\mathrm{Var}(u_{i})}{n\gamma^{2}_{1}\cdot
s\cdot\Lambda_{ij}^{2}}+\frac{4\gamma^{2}_{2}\mathrm{Var}(u_{i})}{\gamma^{2}_{1}\cdot
s\cdot\Lambda_{ij}^{2}}+\frac{4(\mathrm{Var}(u_{j})+\mathrm{Var}(u_{i}))}{n(1-\gamma_{4})^{2}\Lambda_{ij}^{2}}.$

Finally, by choosing $\gamma_{1}=\gamma_{2}=0.1$, $\gamma_{3}=0.18$ and
$\gamma_{4}=0.8$ (satisfying
$\gamma_{4}-2\gamma_{3}-4\gamma_{1}-2\gamma^{2}_{2}/s>0$), we obtain

$\displaystyle\Pr[\text{Err}]~{}$
$\displaystyle\leq~{}\frac{124\mathrm{Var}(u_{i})}{n\Lambda_{ij}^{2}}+\frac{124\mathrm{Var}(u_{i})}{s\Lambda_{ij}^{2}}+\frac{4\mathrm{Var}(u_{i})}{n\Lambda_{ij}^{2}}+\frac{4\mathrm{Var}(u_{i})}{s\Lambda_{ij}^{2}}+\frac{100(\mathrm{Var}(u_{j})+\mathrm{Var}(u_{i}))}{n\Lambda_{ij}^{2}}$
$\displaystyle=~{}\frac{128\mathrm{Var}(u_{i})}{n\Lambda_{ij}^{2}}+\frac{128\mathrm{Var}(u_{i})}{s\Lambda_{ij}^{2}}+\frac{100(\mathrm{Var}(u_{j})+\mathrm{Var}(u_{i}))}{n\Lambda_{ij}^{2}}.$

See 7

Proof  For $c\in\\{i,j\\}$, let $\hat{S}_{c}\sim Q^{n}_{c}$ and $x_{c}\sim
Q_{c}$, $u_{c}=f(x_{c})$, $\mu_{c}=\mu_{f}(Q_{c})$,
$\hat{\mu}_{c}=\mu_{f}(\hat{S}_{c})$ and $\Lambda_{ij}=\|\mu_{i}-\mu_{j}\|$.
Since $\mathbb{I}[r<\Delta]\leq\ell_{\Delta}(-r)$, we have

$\displaystyle\ell_{\Delta}(f;Q_{i},Q_{j})=\mathbb{E}_{\hat{S}_{i},\hat{S}_{j}}\mathbb{E}_{x_{i}}[\ell_{\Delta}(\|u_{i}-\hat{\mu}_{i}\|-\|u_{i}-\hat{\mu}_{j}\|)]\leq\Pr\left[\|u_{i}-\hat{\mu}_{j}\|\leq\|u_{i}-\hat{\mu}_{i}\|+\Delta\right].$

Next, by the triangle inequality and the union bound,

$\displaystyle\Pr\left[\|u_{i}-\hat{\mu}_{j}\|\leq\|u_{i}-\hat{\mu}_{i}\|+\Delta\right]$
$\displaystyle\leq~{}\Pr\left[\|u_{i}-\mu_{j}\|\leq\|u_{i}-\mu_{i}\|+\|\hat{\mu}_{i}-\mu_{i}\|+\|\hat{\mu}_{j}-\mu_{j}\|+\Delta\right]$
$\displaystyle\leq~{}\Pr\left[\|u_{i}-\mu_{j}\|\leq\|u_{i}-\mu_{i}\|+2\Delta~{}\lor~{}\|\hat{\mu}_{i}-\mu_{i}\|\geq
0.5\Delta~{}\lor~{}\|\hat{\mu}_{j}-\mu_{j}\|\geq 0.5\Delta\right]$
$\displaystyle\leq~{}\Pr\left[\|u_{i}-\mu_{j}\|\leq\|u_{i}-\mu_{i}\|+2\Delta\right]+\Pr\left[\|\hat{\mu}_{i}-\mu_{i}\|\geq
0.5\Delta\right]+\Pr\left[\|\hat{\mu}_{j}-\mu_{j}\|\geq 0.5\Delta\right]$
$\displaystyle=~{}E_{2\Delta}(f;Q_{i},Q_{j})+\Pr\left[\|\hat{\mu}_{i}-\mu_{i}\|\geq
0.5\Delta\right]+\Pr\left[\|\hat{\mu}_{j}-\mu_{j}\|\geq 0.5\Delta\right].$

By Markov’s inequality,

$\Pr\left[\|\hat{\mu}_{i}-\mu_{i}\|\geq
0.5\Delta\right]\leq\frac{4\mathrm{Var}(\hat{\mu}_{i})}{\Delta^{2}}=\frac{4\mathrm{Var}(u_{i})}{n^{2}\Delta^{2}}=\frac{\|\mu_{i}-\mu_{j}\|^{2}}{\Delta^{2}}\cdot\frac{4\mathrm{Var}(u_{i})}{n^{2}\|\mu_{i}-\mu_{j}\|^{2}}~{},$

and similarly

$\displaystyle\Pr\left[\|\hat{\mu}_{j}-\mu_{j}\|\geq
0.5\Delta\right]\leq\frac{\|\mu_{i}-\mu_{j}\|^{2}}{\Delta^{2}}\cdot\frac{4\mathrm{Var}(u_{j})}{n^{2}\|\mu_{i}-\mu_{j}\|^{2}}~{}.$

In particular, we obtain that

$\ell_{\Delta}(f;Q_{i},Q_{j})\leq
E_{2\Delta}(f;Q_{i},Q_{j})+{\textstyle\frac{4\|\mu_{i}-\mu_{j}\|^{2}}{n^{2}\Delta^{2}}}\cdot(V_{f}(Q_{i},Q_{j})+V_{f}(Q_{j},Q_{i})).$

Finally, we want to bound $E_{2\Delta}(f;Q_{i},Q_{j})$:

#### (A)

Assume that $\|u_{i}-\mu_{j}\|\leq\|u_{i}-\mu_{i}\|+2\Delta$ and
$\|u_{i}-\mu_{i}\|\leq\alpha_{ij}$. Then,

$\displaystyle\|u_{i}-\mu_{i}\|^{2}+4\alpha_{ij}\Delta+4\Delta^{2}~{}$
$\displaystyle\geq~{}\|u_{i}-\mu_{i}\|^{2}+2\|u_{i}-\mu_{i}\|\cdot
2\Delta+4\Delta^{2}$ $\displaystyle=~{}(\|u_{i}-\mu_{i}\|+2\Delta)^{2}$
$\displaystyle\geq~{}\|u_{i}-\mu_{j}\|^{2}$
$\displaystyle=~{}\|(u_{i}-\mu_{i})-(\mu_{j}-\mu_{i})\|^{2}$
$\displaystyle=~{}\|u_{i}-\mu_{i}\|^{2}-2(u_{i}-\mu_{i})^{\top}(\mu_{j}-\mu_{i})+\Lambda_{ij}^{2}.$

Therefore, we have
$2(u_{i}-\mu_{i})^{\top}(\mu_{j}-\mu_{i})~{}\geq~{}\Lambda_{ij}^{2}-4\alpha_{ij}\Delta-4\Delta^{2}$,
which can also be written as follows

$\displaystyle(u_{i}-\mu_{i})^{\top}e_{ij}~{}\geq~{}\beta_{ij}~{}:=~{}{\textstyle\frac{1}{2}}\Lambda_{ij}-\frac{2(\alpha_{ij}\Delta+\Delta^{2})}{\Lambda_{ij}},$

where $e_{ij}:=(\mu_{j}-\mu_{i})/\|\mu_{j}-\mu_{i}\|$.

#### (B)

We can use this to bound the probability of the event Err in terms of
$\alpha_{ij}$ and $\Delta$ by conditioning:

$\displaystyle E_{2\Delta}(f;Q_{i},Q_{j})~{}$
$\displaystyle=~{}\Pr[\|u_{i}-\mu_{j}\|\leq\|u_{i}-\mu_{i}\|+2\Delta]$
$\displaystyle=~{}\Pr[\|u_{i}-\mu_{j}\|\leq\|u_{i}-\mu_{i}\|+2\Delta~{}\land~{}\|u_{i}-\mu_{i}\|\leq\alpha_{ij}]$
$\displaystyle\quad+\Pr[\|u_{i}-\mu_{j}\|\leq\|u_{i}-\mu_{i}\|+2\Delta~{}\land~{}\|u_{i}-\mu_{i}\|>\alpha_{ij}]$
$\displaystyle\smash{\stackrel{{\scriptstyle(A)}}{{\leq}}}~{}\Pr[(u_{i}-\mu_{i})^{\top}e_{ij}\geq\beta_{ij}]+\Pr[\|u_{i}-\mu_{i}\|>\alpha_{ij}]~{}=:~{}(C)~{}+~{}(D).$

#### (C,D)

Next we upper bound the terms
$\Pr[(u_{i}-\mu_{i})^{\top}e_{ij}\geq\beta_{ij}]$ and
$\Pr[\|u_{i}-\mu_{i}\|>\alpha_{ij}]$. Namely, by Markov’s inequality,

$\displaystyle\Pr[(u_{i}-\mu_{i})^{\top}e_{ij}\geq\beta_{ij}]\leq\frac{\mathrm{Var}(u_{i})}{s\cdot\beta^{2}_{ij}}~{},$

and

$\displaystyle\Pr[\|u_{i}-\mu_{i}\|>\alpha_{ij}]\leq\frac{\mathrm{Var}(u_{i})}{\alpha^{2}_{ij}}~{}.$

We choose $\alpha_{ij}={\textstyle\frac{1}{8}}\Lambda_{ij}^{2}/\Delta$, denote
$\gamma=\Delta/\Lambda_{ij}$ and assume $\Delta\leq\Lambda_{ij}/4$
($\gamma\leq 1/4$). In particular, $\alpha_{ij}=\Lambda_{ij}/(8\gamma)$ and
$\beta_{ij}\geq({\textstyle\frac{1}{4}}-2\gamma^{2})\Lambda_{ij}\geq\Lambda_{ij}/8$.
Hence,

$\displaystyle E_{2\Delta}(f;Q_{i},Q_{j}){}$
$\displaystyle\leq~{}\frac{\mathrm{Var}(u_{i})}{s\cdot({\textstyle\frac{1}{4}}-2\gamma^{2})^{2}\Lambda_{ij}^{2}}+\frac{64\gamma^{2}\mathrm{Var}(u_{i})}{\Lambda_{ij}^{2}}$
$\displaystyle\leq~{}\frac{64\mathrm{Var}(u_{i})}{s\Lambda_{ij}^{2}}+\frac{64\gamma^{2}\mathrm{Var}(u_{i})}{\Lambda_{ij}^{2}}$
$\displaystyle\leq~{}64({\textstyle\frac{1}{s}}+\gamma^{2})\cdot
V_{f}(Q_{i},Q_{j})$
$\displaystyle\leq~{}64({\textstyle\frac{1}{s}}+{\textstyle\frac{\Delta^{2}}{\Lambda_{ij}^{2}}})\cdot
V_{f}(Q_{i},Q_{j}).$

#### Combining everything together.

We obtain that

$\displaystyle\ell_{\Delta}(f;Q_{i},Q_{j}){}$
$\displaystyle\leq~{}E_{2\Delta}(f;Q_{i},Q_{j})+{\textstyle\frac{4\|\mu_{i}-\mu_{j}\|^{2}}{n^{2}\Delta^{2}}}\cdot(V_{f}(Q_{i},Q_{j})+V_{f}(Q_{j},Q_{i}))$
$\displaystyle\leq~{}64({\textstyle\frac{1}{s}}+{\textstyle\frac{\Delta^{2}}{\|\mu_{i}-\mu_{j}\|^{2}}})\cdot
V_{f}(Q_{i},Q_{j})+{\textstyle\frac{4\|\mu_{i}-\mu_{j}\|^{2}}{n^{2}\Delta^{2}}}\cdot(V_{f}(Q_{i},Q_{j})+V_{f}(Q_{j},Q_{i})).$

See 8

Proof  Let $\hat{S}_{c}\sim Q^{n}_{c}$ and $x_{c}\sim Q_{c}$,
$u_{c}=f(x_{c})$, $\mu_{c}=\mu_{f}(Q_{c})$,
$\hat{\mu}_{c}=\mu_{f}(\hat{S}_{c})$ and $\sigma^{2}_{c}=\mathrm{Var}(u_{c})$.
Let $\Lambda_{ij}:=\|\mu_{i}-\mu_{j}\|$ be the considered centers’ distance,
which sets the overall scale. By the triangle inequality,

$\displaystyle\Pr\left[\|u_{i}-\hat{\mu}_{j}\|\leq\|u_{i}-\hat{\mu}_{i}\|+\Delta\right]{}$
$\displaystyle\leq~{}\Pr\left[\Lambda_{ij}\leq
2\|u_{i}-\mu_{i}\|+\|\hat{\mu}_{i}-\hat{\mu}_{i}\|+\|\hat{\mu}_{j}-\mu_{j}\|+\Delta\right]$
(50)
$\displaystyle\leq~{}\Pr\left[(\Lambda_{ij}-\Delta)/4\leq\|u_{i}-\mu_{i}\|\right]$
$\displaystyle\quad+\Pr\left[(\Lambda_{ij}-\Delta)/4\leq\|\hat{\mu}_{i}-\hat{\mu}_{i}\|\right]$
$\displaystyle\quad+\Pr\left[(\Lambda_{ij}-\Delta)/4\leq\|\hat{\mu}_{j}-\mu_{j}\|\right].$

Since $u_{c}\sim\mathcal{N}(\mu_{c},\sigma_{c}^{2}\cdot\mathbb{I}/p)$ and
$\hat{\mu}_{c}\sim\mathcal{N}(\mu_{c},\sigma_{c}^{2}\cdot\mathbb{I}/(np))$,
the random variables $E_{i}=\frac{\|u_{i}-\mu_{i}\|^{2}p}{\sigma^{2}_{i}}$ and
$U_{i}=\frac{\|\hat{\mu}_{i}-\mu_{i}\|^{2}p}{n\sigma^{2}_{i}}$ are chi-squared
random variables with $p$ components. We can write

$\displaystyle\Pr\left[\|u_{i}-\mu_{i}\|\geq{\textstyle\frac{\Lambda_{ij}-\Delta}{4}}\right]=\Pr\left[E_{i}\geq{\textstyle\frac{(\Lambda_{ij}-\Delta)^{2}p}{16\sigma^{2}_{i}}}\right]\leq
1-F\left({\textstyle\frac{(\Lambda_{ij}-\Delta)^{2}p}{16\sigma^{2}_{i}}},p\right)=1-F\left(zp,p\right),$

where $F(x,p)$ is the CDF of the chi-squared distribution with $p$ components
and $z=\left(\frac{\Lambda_{ij}-\Delta}{4\sigma_{c}}\right)^{2}$. By
Chernoff’s bound, if $z>1$, we have

$1-F\left(zp,~{}p\right)\leq\left(z\exp(1-z)\right)^{p/2}\leq\frac{\exp\left(-p/(32V^{ij}_{f})\right)}{(\textnormal{e}\cdot
V^{ij}_{f})^{p/2}}$ (51)

Similarly, we have

$\Pr\left[\|\hat{\mu}_{i}-\mu_{i}\|\geq{\textstyle\frac{\Lambda_{ij}-\Delta}{4}}\right]\leq\frac{\exp\left(-np/(32V^{ij}_{f})\right)}{(\textnormal{e}\cdot
V^{ij}_{f})^{p/2}}\leq\frac{\exp\left(-p/(32V^{ij}_{f})\right)}{(\textnormal{e}\cdot
V^{ij}_{f})^{p/2}}.$ (52)

Finally, by combining (50)-(52) we obtain the desired inequality.

##  D List of Notation

$\mathcal{X}$, $\mathcal{Y}_{k}$ |  |  instances and labels spaces
---|---|---
$\tilde{P}$, $P$ |  |  source and target distributions
$\tilde{P}_{c}$, $P_{c}$ |  |  source and target class-conditional distributions
$\mathbb{I}$ |  |  indicator function
$L_{\tilde{P}}(h)$, $L_{P}(h)$ |  |  source and target test errors
$L_{\mathcal{D}}(f)$ |  |  transfer risk
$c\in[l]|[k]$ |  |  class index
$\mathcal{E}$ |  |  set of class-conditional distributions
$\mathcal{D}$ |  |  distribution over class-conditional distributions
$l\in\mathbb{N}$ |  |  number of classes in the train=source task
$k\in\mathbb{N}$ |  |  number of classes in the test=target task
$m,n$ |  |  number of samples per source/target class
$p\in\mathbb{N}$ |  |  features dimension
$d\in\mathbb{N}$ |  |  input dimension
$q\in\mathbb{N}$ |  |  depth of a neural network
$f\in\mathcal{F}$ |  |  feature map $\mathbb{R}^{d}\to\mathbb{R}^{p}$. Penultimate layer of a neural network
$\mathcal{C}(f)$ |  |  complexity of a neural network
$g\in\mathcal{G}$ |  |  mapping $\mathbb{R}^{p}\to\mathbb{R}^{k}$ from features to soft test labels
$\tilde{g}\in\tilde{\mathcal{G}}$ |  |  mapping $\mathbb{R}^{p}\to\mathbb{R}^{l}$ from features to soft train labels
$\lambda_{n}$ |  |  regularization parameter $\lambda_{n}=\alpha\sqrt{n}$
$\Delta>0$ |  |  margin
$B>0$ |  |  maximal norm of $x\in\mathcal{X}$
$\ell_{\Delta}$ |  |  the soft-margin loss function
$S,S_{c}$ |  |  target training data and its subset consisting of all samples from class $c$
$\tilde{S},\tilde{S}_{c}$ |  |  source training data and its subset consisting of all samples from class $c$
$\mu_{f}(P)$ |  |  mean of $f(x)$ for $x\sim P$
$\mathrm{Var}_{f}(P)$ |  |  variance of $f(x)$ for $x\sim P$
$V_{f}(Q_{1},Q_{2})$ |  |  class-distance normalized variance (CDNV)
$\Lambda>0$ |  |  minimal distance between the embeddings means of two different source classes
$E_{\Delta}(f;Q_{1},Q_{2})$ |  |  the expected (one sided) $\Delta$-margin NCC classification error of $f$
$\ell_{\Delta}(f;Q_{1},Q_{2})$ |  |  the expected (one sided) soft-margin NCC classification error of $f$
$R(A)$ |  |  Rademacher complexity of $A\subset\mathbb{R}^{d}$
$\operatorname{Avg}^{k}_{i=1}[a_{i}]$ |  |  average
$u(A)$ |  |  the projection of $u:B\to R$ over $A\subset B$
$\mathcal{U}(A)$ |  |  the set $\\{u(A):u\in\mathcal{U}\\}$

## Acknowledgements

The authors would like to thank Ilja Kuzborskij, Razvan Pascanu, Tomaso
Poggio, and Csaba Szepesvári for illuminating discussions during the
preparation of this manuscript, and Miruna Pîslar for her priceless technical
support. During parts of this work, Tomer Galanti was a Research Scientist
Intern at Deepmind, and then was supported by the Center for Minds, Brains and
Machines (CBMM), funded by NSF STC award CCF-1231216.

## References

  * Bartlett et al. (2017) Peter L. Bartlett, Dylan J. Foster, and Matus Telgarsky.  Spectrally-normalized margin bounds for neural networks.  In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , page 6241–6250, 2017.
  * Baxter (2000) Jonathan Baxter.  A model of inductive bias learning.  _J. Artif. Int. Res._ , 12(1):149–198, March 2000.
  * Bengio (2012) Yoshua Bengio.  Deep learning of representations for unsupervised and transfer learning.  In _Proceedings of ICML Workshop on Unsupervised and Transfer Learning_ , volume 27 of _Proceedings of Machine Learning Research_ , pages 17–36, Bellevue, Washington, USA, 02 Jul 2012. PMLR.
  * Bertinetto et al. (2019) Luca Bertinetto, Joao F. Henriques, Philip Torr, and Andrea Vedaldi.  Meta-learning with differentiable closed-form solvers.  In _International Conference on Learning Representations_ , 2019.
  * Bommasani et al. (2021) Rishi Bommasani, Drew A. Hudson, Ehsan Adeli, Russ Altman, Simran Arora, Sydney von Arx, Michael S. Bernstein, Jeannette Bohg, Antoine Bosselut, Emma Brunskill, Erik Brynjolfsson, Shyamal Buch, Dallas Card, Rodrigo Castellon, Niladri Chatterji, Annie S. Chen, Kathleen Creel, Jared Quincy Davis, Dorottya Demszky, Chris Donahue, Moussa Doumbouya, Esin Durmus, Stefano Ermon, John Etchemendy, Kawin Ethayarajh, Li Fei-Fei, Chelsea Finn, Trevor Gale, Lauren Gillespie, Karan Goel, Noah D. Goodman, Shelby Grossman, Neel Guha, Tatsunori Hashimoto, Peter Henderson, John Hewitt, Daniel E. Ho, Jenny Hong, Kyle Hsu, Jing Huang, Thomas Icard, Saahil Jain, Dan Jurafsky, Pratyusha Kalluri, Siddharth Karamcheti, Geoff Keeling, Fereshte Khani, Omar Khattab, Pang Wei Koh, Mark S. Krass, Ranjay Krishna, Rohith Kuditipudi, and et al.  On the opportunities and risks of foundation models.  _CoRR_ , abs/2108.07258, 2021.
  * Brown et al. (2020) Tom Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared D Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, Sandhini Agarwal, Ariel Herbert-Voss, Gretchen Krueger, Tom Henighan, Rewon Child, Aditya Ramesh, Daniel Ziegler, Jeffrey Wu, Clemens Winter, Chris Hesse, Mark Chen, Eric Sigler, Mateusz Litwin, Scott Gray, Benjamin Chess, Jack Clark, Christopher Berner, Sam McCandlish, Alec Radford, Ilya Sutskever, and Dario Amodei.  Language models are few-shot learners.  In _Advances in Neural Information Processing Systems_ , pages 1877–1901, 2020.
  * Caruana (1995) Rich Caruana.  Learning many related tasks at the same time with backpropagation.  In _Advances in Neural Information Processing Systems_ , 1995.
  * Chen et al. (2019a) Chaofan Chen, Oscar Li, Daniel Tao, Alina Barnett, Cynthia Rudin, and Jonathan K Su.  This looks like that: Deep learning for interpretable image recognition.  In _Advances in Neural Information Processing Systems_ , 2019a.
  * Chen et al. (2018) Liang-Chieh Chen, Yukun Zhu, George Papandreou, Florian Schroff, and Hartwig Adam.  Encoder-decoder with atrous separable convolution for semantic image segmentation.  In _ECCV (7)_ , volume 11211 of _Lecture Notes in Computer Science_ , pages 833–851. Springer, 2018.
  * Chen et al. (2019b) Sihong Chen, Kai Ma, and Yefeng Zheng.  Med3d: Transfer learning for 3d medical image analysis.  _CoRR_ , abs/1904.00625, 2019b.
  * Cohen et al. (2017) Gregory Cohen, Saeed Afshar, Jonathan Tapson, and André van Schaik.  EMNIST: an extension of MNIST to handwritten letters.  _arXiv preprint arXiv:1702.05373_ , 2017.
  * Dhillon et al. (2020) Guneet Singh Dhillon, Pratik Chaudhari, Avinash Ravichandran, and Stefano Soatto.  A baseline for few-shot image classification.  In _International Conference on Learning Representations_ , 2020.
  * Du et al. (2021) Simon Shaolei Du, Wei Hu, Sham M. Kakade, Jason D. Lee, and Qi Lei.  Few-shot learning via learning the representation, provably.  In _International Conference on Learning Representations_ , 2021.
  * Finn et al. (2017) Chelsea Finn, Pieter Abbeel, and Sergey Levine.  Model-agnostic meta-learning for fast adaptation of deep networks.  In _Proceedings of the 34th International Conference on Machine Learning_ , pages 1126–1135, 2017.
  * Galanti et al. (2016) Tomer Galanti, Lior Wolf, and Tamir Hazan.  A theoretical framework for deep transfer learning.  _Information and Inference: A Journal of the IMA_ , 5:159–209, 2016.
  * Galanti et al. (2022a) Tomer Galanti, Liane Galanti, and Ido Ben-Shaul.  On the implicit bias towards minimal depth of deep neural networks, 2022a.  URL https://arxiv.org/abs/2202.09028.
  * Galanti et al. (2022b) Tomer Galanti, András György, and Marcus Hutter.  On the role of neural collapse in transfer learning.  In _International Conference on Learning Representations_ , 2022b.  URL https://openreview.net/forum?id=SwIp410B6aQ.
  * Gidaris and Komodakis (2018) Spyros Gidaris and Nikos Komodakis.  Dynamic few-shot visual learning without forgetting.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2018.
  * Goldblum et al. (2020) Micah Goldblum, Steven Reich, Liam Fowl, Renkun Ni, Valeriia Cherepanova, and Tom Goldstein.  Unraveling meta-learning: Understanding feature representations for few-shot tasks.  In _Proceedings of the 37th International Conference on Machine Learning_ , pages 3607–3616, 2020.
  * Golowich et al. (2017) Noah Golowich, Alexander Rakhlin, and Ohad Shamir.  Size-independent sample complexity of neural networks.  _Information and Inference: A Journal of the IMA_ , 9, 12 2017.  doi: 10.1093/imaiai/iaz007.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 770–778, 2016.  doi: 10.1109/CVPR.2016.90.
  * He et al. (2017) Kaiming He, Georgia Gkioxari, Piotr Dollar, and Ross Girshick.  Mask r-cnn.  In _Proceedings of the IEEE International Conference on Computer Vision (ICCV)_ , Oct 2017.
  * Howard et al. (2017) Andrew G. Howard, Menglong Zhu, Bo Chen, Dmitry Kalenichenko, Weijun Wang, Tobias Weyand, Marco Andreetto, and Hartwig Adam.  Mobilenets: Efficient convolutional neural networks for mobile vision applications, 2017.
  * Huang and Li (2020) Zixuan Huang and Yin Li.  Interpretable and accurate fine-grained recognition via region grouping.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 8662–8672, 2020.
  * Khodak et al. (2019) Mikhail Khodak, Maria-Florina F Balcan, and Ameet S Talwalkar.  Adaptive gradient-based meta-learning methods.  In _Advances in Neural Information Processing Systems_ , volume 32, 2019.
  * Krizhevsky (2012) Alex Krizhevsky.  Learning multiple layers of features from tiny images.  Technical report, University of Toronto, 05 2012.
  * Lee et al. (2019) Kwonjoon Lee, Subhransu Maji, Avinash Ravichandran, and Stefano Soatto.  Meta-learning with differentiable convex optimization.  In _CVPR_ , 2019.
  * Li et al. (2018) Xingguo Li, Junwei Lu, Zhaoran Wang, Jarvis Haupt, and Tuo Zhao.  On tighter generalization bound for deep neural networks: CNNs, ResNets, and beyond, 2018.  URL https://arxiv.org/abs/1806.05159.
  * Liu et al. (2019) Yanbin Liu, Juho Lee, Minseop Park, Saehoon Kim, Eunho Yang, Sungju Hwang, and Yi Yang.  Learning to propagate labels: Transductive propagation network for few-shot learning.  In _International Conference on Learning Representations_ , 2019.
  * Long et al. (2015) Jonathan Long, Evan Shelhamer, and Trevor Darrell.  Fully convolutional networks for semantic segmentation.  In _2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 3431–3440, 2015.  doi: 10.1109/CVPR.2015.7298965.
  * Maurer and Pontil (2019) Andreas Maurer and M. Pontil.  Uniform concentration and symmetrization for weak interactions.  In _COLT_ , 2019.
  * Maurer et al. (2016) Andreas Maurer, Massimiliano Pontil, and Bernardino Romera-Paredes.  The benefit of multitask representation learning.  _J. Mach. Learn. Res._ , 17(1):2853–2884, January 2016.
  * Mishra et al. (2018) Nikhil Mishra, Mostafa Rohaninejad, Xi Chen, and Pieter Abbeel.  A simple neural attentive meta-learner.  In _International Conference on Learning Representations_ , 2018.
  * Mixon et al. (2020) Dustin G. Mixon, Hans Parshall, and Jianzong Pi.  Neural collapse with unconstrained features.  _CoRR_ , abs/2011.11619, 2020.  URL https://arxiv.org/abs/2011.11619.
  * Mohri et al. (2018) Mehryar Mohri, Afshin Rostamizadeh, and Ameet Talwalkar.  _Foundations of Machine Learning_.  Adaptive Computation and Machine Learning. MIT Press, Cambridge, MA, 2nd edition, 2018.
  * Munkhdalai et al. (2018) Tsendsuren Munkhdalai, Xingdi Yuan, Soroush Mehri, and Adam Trischler.  Rapid adaptation with conditionally shifted neurons.  In _Proceedings of the 35th International Conference on Machine Learning_ , pages 3664–3673, 2018.
  * Neyshabur et al. (2015) Behnam Neyshabur, Ryota Tomioka, and Nathan Srebro.  Norm-based capacity control in neural networks.  In _Proceedings of The 28th Conference on Learning Theory_ , pages 1376–1401, 2015.
  * Oreshkin et al. (2018) Boris Oreshkin, Pau Rodríguez López, and Alexandre Lacoste.  Tadam: Task dependent adaptive metric for improved few-shot learning.  In _Advances in Neural Information Processing Systems_ , 2018.
  * Papyan et al. (2020) Vardan Papyan, X. Y. Han, and David L. Donoho.  Prevalence of neural collapse during the terminal phase of deep learning training.  _Proceedings of the National Academy of Sciences_ , 117(40):24652–24663, 2020.
  * Pentina and Lampert (2014) Anastasia Pentina and Christoph Lampert.  A pac-bayesian bound for lifelong learning.  In _Proceedings of the 31st International Conference on Machine Learning_ , pages 991–999, 2014.
  * Qiao et al. (2019) Limeng Qiao, Yemin Shi, Jia Li, Yaowei Wang, Tiejun Huang, and Yonghong Tian.  Transductive episodic-wise adaptive metric for few-shot learning.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV)_ , October 2019.
  * Qiao et al. (2018) Siyuan Qiao, Chenxi Liu, Wei Shen, and Alan L. Yuille.  Few-shot image recognition by predicting parameters from activations.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2018.
  * Ramesh et al. (2021) Aditya Ramesh, Mikhail Pavlov, Gabriel Goh, Scott Gray, Chelsea Voss, Alec Radford, Mark Chen, and Ilya Sutskever.  Zero-shot text-to-image generation.  In _Proceedings of the 38th International Conference on Machine Learning_ , pages 8821–8831, 2021.
  * Ravi and Larochelle (2017) Sachin Ravi and Hugo Larochelle.  Optimization as a model for few-shot learning.  In _International Conference on Learning Representations_ , 2017.
  * Ravichandran et al. (2019) Avinash Ravichandran, Rahul Bhotika, and Stefano Soatto.  Few-shot learning with embedded class models and shot-free meta training.  _2019 IEEE/CVF International Conference on Computer Vision (ICCV)_ , pages 331–339, 2019.
  * Redmon et al. (2016) Joseph Redmon, Santosh Divvala, Ross Girshick, and Ali Farhadi.  You only look once: Unified, real-time object detection.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2016.
  * Ren et al. (2015) Shaoqing Ren, Kaiming He, Ross Girshick, and Jian Sun.  Faster r-cnn: Towards real-time object detection with region proposal networks.  In _Advances in Neural Information Processing Systems_ , volume 28, 2015.
  * Russakovsky et al. (2015) Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, Alexander C. Berg, and Li Fei-Fei.  ImageNet Large Scale Visual Recognition Challenge.  _International Journal of Computer Vision (IJCV)_ , 115(3):211–252, 2015.
  * Rusu et al. (2019) Andrei A. Rusu, Dushyant Rao, Jakub Sygnowski, Oriol Vinyals, Razvan Pascanu, Simon Osindero, and Raia Hadsell.  Meta-learning with latent embedding optimization.  In _International Conference on Learning Representations_ , 2019.
  * Simonyan and Zisserman (2014) Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  _CoRR_ , abs/1409.1556, 2014.
  * Snell et al. (2017) Jake Snell, Kevin Swersky, and Richard Zemel.  Prototypical networks for few-shot learning.  In _Advances in Neural Information Processing Systems_ , 2017.
  * Sun et al. (2019) Qianru Sun, Yaoyao Liu, Tat-Seng Chua, and Bernt Schiele.  Meta-transfer learning for few-shot learning.  In _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019.
  * Sung et al. (2018) Flood Sung, Yongxin Yang, Li Zhang, Tao Xiang, Philip H.S. Torr, and Timothy M. Hospedales.  Learning to compare: Relation network for few-shot learning.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2018.
  * Sutskever et al. (2013) Ilya Sutskever, James Martens, George E. Dahl, and Geoffrey E. Hinton.  On the importance of initialization and momentum in deep learning.  In _Proceedings of the 30th International Conference on Machine Learning, ICML 2013, Atlanta, GA, USA, 16-21 June 2013_ , pages 1139–1147, 2013\.
  * Tian et al. (2020) Yonglong Tian, Yue Wang, Dilip Krishnan, Joshua B. Tenenbaum, and Phillip Isola.  Rethinking few-shot image classification: A good embedding is all you need?  In _Proceedings of the European Conference on Computer Vision (ECCV)_ , volume 12359 of _Lecture Notes in Computer Science_ , pages 266–282. Springer, 2020.
  * Tirer and Bruna (2022) Tom Tirer and Joan Bruna.  Extended unconstrained features model for exploring deep neural collapse.  In _Proceedings of the 39th International Conference on Machine Learning_ , pages 21478–21505, 2022.
  * Vinyals et al. (2016) Oriol Vinyals, Charles Blundell, Timothy Lillicrap, koray kavukcuoglu, and Daan Wierstra.  Matching networks for one shot learning.  In _Advances in Neural Information Processing Systems_ , 2016.
  * Wainwright (2019) Martin Wainwright.  _High-dimensional statistics : a non-asymptotic viewpoint._ Cambridge University Press, Cambridge, 2019.
  * Xu et al. (2021) Mengjia Xu, Akshay Rangamani, Andrzej Banburski, Qianli Liao, Tomer Galanti, and Tomaso Poggio.  Dynamics and neural collapse in deep classifiers trained with the square loss.  Technical report, Center for Brains, Minds and Machines (CBMM), 2021.
  * Yan et al. (2018) Yan Yan, Tianbao Yang, Zhe Li, Qihang Lin, and Yi Yang.  A unified analysis of stochastic momentum methods for deep learning.  In _Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence_ , pages 2955–2961, 2018.
  * Yang et al. (2018) Ze Yang, Tiange Luo, Dong Wang, Zhiqiang Hu, Jun Gao, and Liwei Wang.  Learning to navigate for fine-grained classification.  In _Proceedings of the European Conference on Computer Vision (ECCV)_ , September 2018.
  * Yosinski et al. (2014) Jason Yosinski, Jeff Clune, Yoshua Bengio, and Hod Lipson.  How transferable are features in deep neural networks?  In _Advances in Neural Information Processing Systems_ , 2014.
  * Zagoruyko and Komodakis (2016) Sergey Zagoruyko and N. Komodakis.  Wide residual networks.  _ArXiv_ , abs/1605.07146, 2016.
  * Zhu et al. (2021) Zhihui Zhu, Tianyu DING, Jinxin Zhou, Xiao Li, Chong You, Jeremias Sulam, and Qing Qu.  A geometric analysis of neural collapse with unconstrained features.  In _Advances in Neural Information Processing Systems_ , 2021.