# Analysis of the rSVDdpd Algorithm: A Robust Singular Value Decomposition
Method using Density Power Divergence

Subhrajyoty Roy
Interdisciplinary Statistical Research Unit
Indian Statistical Insititute, Kolkata
<EMAIL_ADDRESS>
&Abhik Ghosh
Interdisciplinary Statistical Research Unit
Indian Statistical Insititute, Kolkata
<EMAIL_ADDRESS>
&Ayanendranath Basu
Interdisciplinary Statistical Research Unit
Indian Statistical Insititute, Kolkata
<EMAIL_ADDRESS>

###### Abstract

The traditional method of computing singular value decomposition (SVD) of a
data matrix is based on a least squares principle, thus, is very sensitive to
the presence of outliers. Hence the resulting inferences across different
applications using the classical SVD are extremely degraded in the presence of
data contamination (e.g., video surveillance background modelling tasks,
etc.). A robust singular value decomposition method using the minimum density
power divergence estimator (rSVDdpd) has been found to provide a satisfactory
solution to this problem and works well in applications. For example, it
provides a neat solution to the background modelling problem of video
surveillance data in the presence of camera tampering. In this paper, we
investigate the theoretical properties of the rSVDdpd estimator such as
convergence, equivariance and consistency under reasonable assumptions. Since
the dimension of the parameters, i.e., the number of singular values and the
dimension of singular vectors can grow linearly with the size of the data, the
usual M-estimation theory has to be suitably modified with concentration
bounds to establish the asymptotic properties. We believe that we have been
able to accomplish this satisfactorily in the present work. We also
demonstrate the efficiency of rSVDdpd through extensive simulations.

_Keywords_ Singular value decomposition $\cdot$ Matrix factorization $\cdot$
Density power divergence

##  1 Introduction

Singular value decomposition (SVD) is a matrix factorization method that
breaks down a real or complex matrix into three parts: two orthogonal matrices
consisting of the singular vectors and one diagonal matrix with non-negative
diagonal entries made up of singular values. SVD is commonly viewed as a low-
rank approximation of a linear transformation, representing the transformation
as a sequence of rotations and dilations. The applications of SVD are diverse.
Its mathematical applications include computing pseudoinverses or Moore-
Penrose inverses of matrices, efficiently solving systems of homogeneous
linear equations, determining range space, null space, and rank of a linear
transformation, and finding ordinary least square solutions using the Golub
and Reinsch algorithm (Golub and Reinsch, 1970). SVD has also been widely
utilized in statistical and machine-learning methods for data analysis and
modelling in various real-world applications. One notable application is
principal component analysis (PCA), which employs SVD to decompose a data
matrix into a lower-dimensional representation that captures the maximum
variability in the original data using a reduced set of variables. Other
popular dimension reduction techniques such as correspondence analysis
(Greenacre, 2017), latent semantic indexing (Hofmann, 1999; Anandarajan et
al., 2019) and clustering techniques (Drineas et al., 2004; Cheng et al.,
2019) also rely on SVD as a fundamental component. SVD is extensively employed
in pattern recognition within signal, image, and video processing domains. Its
applications include image watermarking schemes (Dappuri et al., 2020), signal
denoising and feature enhancements (Zhao and Jia, 2017), audio watermarking
(Rezaei and Khalili, 2019), sound source localization (Grondin and Glass,
2019), and sound recovery techniques (Zhang et al., 2016). Moreover, SVD has
gained prominence in the field of bioinformatics, where it is used for
analyzing protein functional associations (Franceschini et al., 2015),
clustering gene expression data (Bustamam et al., 2018), and predicting
protein-coding regions (Das et al., 2017). In the realm of geographical
science, Kumar et al. (2011) employed SVD-based techniques, including its
robust variant, to generate accurate graphical representations of climate
data, mitigating the impact of extreme weather events such as thunderstorms
and heavy rainfall. Such a wide range of applications clearly underscore the
relevance of SVD as an extremely integral component of data analysis across a
multitude of disciplines.

However, as already indicated by several authors (Hawkins et al., 2001); (Liu
et al., 2003); (Kumar et al., 2011), the usual method of computing SVD is
highly susceptible to outliers present in the data matrix. In contrast, as the
data are becoming increasingly vast and complex in the recent era, it is also
susceptible to the inclusion of different forms of noises, corruptions and
contamination by outlying observations. This readily introduces the need for a
robust algorithm for the computation of SVD which remains unaffected (or is
minimally affected) due to the presence of outliers, and thus leads to more
stable and trustworthy solutions in different applications mentioned
previously.

###  1.1 Related Works and Our Contribution

Ammann (1993) was one of the early pioneers in developing a robust version of
the SVD. He treated it as a special case of the projection pursuit problem to
be solved using the transposed QR algorithm. Other researchers, such as
Hawkins et al. (2001); Liu et al. (2003) and Ke and Kanade (2005), approached
the computation of SVD as a least squares problem and proposed robust
extensions using alternating $L_{1}$ regression algorithms with the least
absolute deviation (LAD) loss function. However, a simple LAD approach is
sensitive to high leverage points, which led to the exploration of weighted
LAD approaches and the utilization of the Huber weight function (Jung, 2010).
In a different attempt, Rey (2007) introduced a robust method called “Total”
SVD, which employed Huber’s weight function and “Total” least squares
(Markovsky and Van Huffel, 2007). This approach accounted for errors in both
the data matrix and the singular vectors, in contrast to the usual least
squares where the only source of error is the response variable. Although this
resulted in a more robust SVD estimate, the method faced several convergence
issues as mentioned by Rey (2007) himself. Alternatively, Zhang et al. (2013)
incorporated the Huber weight function in the loss function and combined it
with a squared error-based penalty function for regularization, creating
another robust SVD estimator. Wang (2017) used an estimator derived from an
$\alpha$-stable distribution with a cost function $\rho(x)=\log(x^{2}+K^{2})$,
where the tuning parameter $K$ provides a balance between robustness and
efficiency. Nevertheless, finding the appropriate tuning parameter $K$ for the
estimator was challenging. Apart from a simple alternating $L_{1}$ regression
approach (Gabriel and Zamir, 1979), there is a lack of theoretical guarantees
or properties for the resulting SVD estimates in the literature. Convergence
and orthogonality of the singular vectors are not assured, and the
computational complexity poses a significant challenge, especially for large
data matrices.

Roy et al. (2021) introduced the rSVDdpd (Robust Singular Value Decomposition
using Density Power Divergence) estimator as an outlier-resistant matrix
factorization technique in the context of video surveillance background
modelling. They transformed the SVD computation into an alternating regression
problem and utilized the density power divergence (DPD) loss function (Basu et
al., 1998). The minimum DPD estimator (MDPDE) has exhibited strong robustness
and high efficiency in statistics and information theory (see, e.g., Basu et
al. (1998); a brief description is also provided in Appendix A). The use of
the DPD loss function in the proposed rSVDdpd estimator benefits from these
desirable statistical properties. This paper aims to provide theoretical
justifications for the rSVDdpd algorithm by establishing its convergence and
mathematical properties like equivariance and asymptotic consistency. The
primary challenge here is to extend the existing results on the MDPDE type
estimators (Ghosh and Basu, 2013) to the case where the dimension of the
parameter grows to infinity linearly in sample size, and thus, concentration
bounds are needed to ensure that the desirable asymptotic properties hold. We
also conduct extensive simulation studies to demonstrate its applicability as
a general-purpose robust matrix factorization technique and its performance
compared to the existing approaches by Zhang et al. (2013) and Hawkins et al.
(2001), where rSVDdpd outperforms them in most simulation setups.

##  2 Description of the rSVDdpd estimator

The problem of Singular Value Decomposition starts with a data matrix
$\boldsymbol{X}$ of dimension $n\times p$ ($n$ and $p$ may be different),
admitting an approximate low-rank representation of the form

$\boldsymbol{X}=\sum_{k=1}^{r}\lambda_{k}\boldsymbol{u}_{k}\boldsymbol{v}_{k}^{\intercal}+\boldsymbol{E},$

where $\\{\boldsymbol{u}_{k}\\}_{k=1}^{r}$ is a set of $r$ orthonormal vectors
of length $n$ and $\\{\boldsymbol{v}_{k}\\}_{k=1}^{r}$ is a set of $r$
orthonormal vectors of length $p$. The $n\times p$ dimensional matrix
$\boldsymbol{E}$ constitutes of the errors $e_{ij}$s, which are generally
expected to be smaller in magnitude than the corresponding entries of the data
matrix $\boldsymbol{X}$, when the data are not contaminated. The vectors
$\boldsymbol{u}_{k}$s and $\boldsymbol{v}_{k}$s are called the left and the
right singular vectors, respectively, whereas $\lambda_{k}$s are called the
singular values of the matrix $\boldsymbol{X}$. As in Rey (2007), the above
description of the SVD can also be reformalized into

$\boldsymbol{X}=\sum_{k=1}^{r}\boldsymbol{a}_{k}\boldsymbol{b}_{k}^{\intercal}+E,$
(2.1)

where $\boldsymbol{a}_{k}$s and $\boldsymbol{b}_{k}$s are still orthogonal
sets of vectors for $k=1,2,\dots r$, but not necessarily normalized. Once the
estimates of $\boldsymbol{a}_{k}$s and $\boldsymbol{b}_{k}$s are known, they
can be normalized to obtain the $\boldsymbol{u}_{k}$s and
$\boldsymbol{v}_{k}$s and the singular values are then given by
$\lambda_{k}=\|\boldsymbol{a}_{k}\|\|\boldsymbol{b}_{k}\|$ for each
$k=1,\ldots,r$, where $\|\cdot\|$ denotes the usual Euclidean ($L_{2}$) norm.

Roy et al. (2021) describe the rSVDdpd estimator for a rank one decomposition
of the form

$X_{ij}=a_{i}b_{j}+e_{ij},\qquad i=1,2,\dots n.$ (2.2)

Following Rey (2007), they view (2.2) as a linear regression equation instead
of a singular value decomposition problem. For a fixed index $j$, (2.2) can be
interpreted as a linear regression with $X_{ij}$s as an observed response,
$a_{i}$s as the covariate values, $b_{j}$s as the regression slope parameters
to be estimated and $e_{ij}$s as the random error component. As we vary the
column index $j=1,2,\dots p$, we are posed with $p$ such linear regression
problems, solving each of them will jointly yield an estimate of
$\boldsymbol{b}=(b_{1},\ldots,b_{p})$ given the values of $a_{i}$s. Now, one
can interchange the role of $a_{i}$s and $b_{j}$s and view (2.2) as regression
of $X_{ij}$ on $b_{j}$s for fixed $i$, and estimate
$\boldsymbol{a}=(a_{1},\ldots,a_{n})$. Then, one can alternatively estimate
$\boldsymbol{a}$ and $\boldsymbol{b}$ by solving these linear regression
problems, and the converged values of $(\boldsymbol{a},\boldsymbol{b})$ will
yield the required decomposition. Roy et al. (2021) aim to solve these
regression problems using the popular minimum density power divergence
estimator (MDPDE) introduced in Ghosh and Basu (2013). Modelling the error
components $e_{ij}$ as i.i.d. normal random variables with mean $0$ and
variance $\sigma^{2}$, the MDPDE aims to minimize the loss function

$H_{n,p}(\boldsymbol{a},\boldsymbol{b},\sigma^{2})=\dfrac{1}{np}\sum_{i=1}^{n}\sum_{j=1}^{p}V(x_{ij};a_{i},b_{j},\sigma^{2}),$
(2.3)

where

$V(y;c,d,\sigma^{2})=\frac{1}{(2\pi)^{\alpha/2}\sigma^{\alpha}\sqrt{1+\alpha}}-\frac{(1+\alpha)}{\alpha}\frac{1}{(2\pi)^{\alpha/2}\sigma^{\alpha}}\exp\left\\{-\alpha\frac{(y-cd)^{2}}{2\sigma^{2}}\right\\}$
(2.4)

for any $y,c,d\in\mathbb{R}$ and $\sigma>0$. So, instead of estimating the
elements of the scaled singular vectors $\boldsymbol{a}$ and $\boldsymbol{b}$
i.e. $a_{i}$s and $b_{j}$s in the alternating regression models, we can indeed
minimize the MDPDE objective function iteratively over either of the
parameters ($a_{i}$, $b_{j}$ or $\sigma^{2}$) given the most recent estimates
of the other parameters and normalize them later to obtain estimates of the
singular values. For any non-negative integer $t$, let us denote
$\boldsymbol{\theta}^{(t)}=\left(\lambda^{(t)},\\{a_{i}^{(t)}\\},\\{b_{j}^{(t)}\\},(\sigma^{2})^{(t)}\right)$
as the estimate obtained at $t$-th iteration of the alternating solution of
the regression problem through MDPDE. Clearly, these are related in the
following way

$\displaystyle\lambda^{(t+1/2)}a_{i}^{(t+1)}$
$\displaystyle=\operatorname*{arg\,min}_{a}\dfrac{1}{p}\sum_{j=1}^{p}V(x_{ij};a,b_{j}^{(t)},(\sigma^{2})^{(t)}),\qquad
i=1,\ldots n,\ t=0,1,\ldots,$ (2.5)
$\displaystyle\lambda^{(t+1)}b_{j}^{(t+1)}$
$\displaystyle=\operatorname*{arg\,min}_{b}\dfrac{1}{n}\sum_{i=1}^{n}V(x_{ij};a_{i}^{(t+1)},b,(\sigma^{2})^{(t)}),\qquad
j=1,\ldots p,\ t=0,1,\ldots,$ (2.6) $\displaystyle(\sigma^{2})^{(t+1)}$
$\displaystyle=\operatorname*{arg\,min}_{\sigma^{2}}\dfrac{1}{np}\sum_{i=1}^{n}\sum_{j=1}^{p}V(x_{ij};a_{i}^{(t+1)},b_{j}^{(t+1)},\sigma^{2}),\qquad
t=0,1,\ldots,$ (2.7)

where $x_{ij}$ is the $(i,j)$-th entry of the data matrix $\boldsymbol{X}$ and
$V$ being the DPD loss function given in (2.4). As will be shown later in
Theorem 3.7, these sequence of estimators $\boldsymbol{\theta}^{(t)}$
converges to $\boldsymbol{\theta}^{\ast}$ in the parameter space
$\boldsymbol{\Theta}$, which is a $(n+p)$ dimensional manifold in
$\mathbb{R}^{(n+p+2)}$, namely $[0,\infty)\times S_{n}\times
S_{p}\times[0,\infty)$, with $S_{n}$ and $S_{p}$ denoting the $n$ and
$p$-dimensional hyperspheres of unit radius centered at the origin
respectively. This estimator
$\boldsymbol{\theta}^{\ast}=\left(\lambda^{\ast},\\{a_{i}^{\ast}\\}_{i=1}^{n},\\{b_{j}^{\ast}\\}_{j=1}^{p},(\sigma^{2})^{\ast}\right)$
is the rSVDdpd estimator of the rank-one decomposition of the matrix
$\boldsymbol{X}$.

By standard differentiation, as in Ghosh and Basu (2013), the corresponding
estimating equations for these minimization problems in (2.5)-(2.6) turn out
to be

$\displaystyle\sum_{j=1}^{p}b_{j}(x_{ij}-a_{i}b_{j})e^{-\alpha\dfrac{(x_{ij}-a_{i}b_{j})^{2}}{2\sigma^{2}}}$
$\displaystyle=0,\qquad i=1,\ldots n,$ (2.8)
$\displaystyle\sum_{i=1}^{n}a_{i}(x_{ij}-a_{i}b_{j})e^{-\alpha\dfrac{(x_{ij}-a_{i}b_{j})^{2}}{2\sigma^{2}}}$
$\displaystyle=0,\qquad j=1,\ldots p,$ (2.9)
$\displaystyle\sum_{i=1}^{n}\sum_{j=1}^{p}\left[1-\dfrac{(x_{ij}-a_{i}b_{j})^{2}}{\sigma^{2}}\right]e^{-\alpha\dfrac{(x_{ij}-a_{i}b_{j})^{2}}{2\sigma^{2}}}$
$\displaystyle=\frac{\alpha}{(1+\alpha)^{3/2}}.$ (2.10)

######  Remark 2.1.

As the robustness parameter $\alpha$ tends towards zero, the quantity
$V(y;c,d,\sigma^{2})$ in (2.4) tends towards the normalized squared error
$(y-cd)^{2}/2\sigma^{2}$. Thus, as $\alpha\rightarrow 0$, the MDPDE objective
function in (2.3) will converge to a scaled version of the usual squared error
loss function between the elements $x_{ij}$ of the data matrix and the rank
one estimates $a_{i}b_{j}$, i.e.

$\lim_{\alpha\rightarrow
0}H_{n,p}(\boldsymbol{a},\boldsymbol{b},\sigma^{2})=\dfrac{1}{2np\sigma^{2}}\sum_{i=1}^{n}\sum_{j=1}^{p}(x_{ij}-a_{i}b_{j})^{2}$

The alternating estimating equations (2.8)-(2.10) then turn out to be
explicitly solvable, with

$a_{i}=\dfrac{\sum_{j}x_{ij}b_{j}}{\sum_{j}b_{j}^{2}},\ i=1,\ldots n;\quad
b_{j}=\dfrac{\sum_{i}x_{ij}a_{i}}{\sum_{i}a_{i}^{2}},\ j=1,\ldots
p;\quad\sigma^{2}=\dfrac{1}{np}\sum_{i=1}^{n}\sum_{j=1}^{n}(x_{ij}-a_{i}b_{j})^{2}.$

Hawkins et al. (2001) showed that such alternating estimating equations lead
to the usual least square SVD estimates. Thus, the proposal yields the usual
SVD estimators as the robustness parameter is decreased to zero.

Given that such rank-one decompositions as in (2.2) can be obtained, one can
stack the estimates of $a_{i}$s to get $\widehat{\boldsymbol{a}}_{1}$, and
combine estimates of $b_{j}$s to get $\widehat{\boldsymbol{b}}_{1}$, the
unnormalized first set of the singular vectors. Then these can be normalized
and one obtains an estimate of the first singular value
$\widehat{\lambda}_{1}=\|\widehat{\boldsymbol{a}}_{1}\|\|\widehat{\boldsymbol{b}}_{1}\|$.
For subsequent singular values, one can apply the same estimation algorithm on
the residual matrix
$\boldsymbol{X}-\widehat{\lambda}_{1}\widehat{\boldsymbol{a}}_{1}\widehat{\boldsymbol{b}}_{1}^{\intercal}$.
Such an iterative method of rank one approximation is quite common in the
matrix factorization literature (Hawkins et al., 2001; Cichocki et al., 2011).
However, this method does not guarantee that the subsequent singular vectors
remain orthonormal to the previous set of singular vectors. The orthogonality
property usually degrades as one estimates more singular values. Thus, Roy et
al. (2021) propose to use a Gram-Schmidt orthogonalization trick (Giraud et
al., 2005) in between the iterative equations. In particular, between
alternatively using (2.8)-(2.10), the estimates of the $k$-th singular vectors
$\boldsymbol{a}_{k}$ and $\boldsymbol{b}_{k}$ are updated as
$\boldsymbol{a}_{k}\leftarrow\boldsymbol{a}_{k}-\sum_{r=1}^{(k-1)}{\boldsymbol{a}_{k}}^{\intercal}\boldsymbol{a}_{r}$
and
$\boldsymbol{b}_{k}\leftarrow\boldsymbol{b}_{k}-\sum_{r=1}^{(k-1)}{\boldsymbol{b}_{k}}^{\intercal}\boldsymbol{b}_{r}$
for all but the first singular value. Further details about the estimation
algorithm can be found in Roy et al. (2021).

##  3 Mathematical Properties

Since the estimation process of the subsequent singular values and vectors
follows from the same estimation technique of rank one decomposition on the
residual matrix, we shall restrict our attention to the study of the
properties of the rSVDdpd estimator for the rank one decomposition only. In
order to make sure all the results developed in this section are true for the
subsequent singular values and vectors as well, the assumptions on the
distribution of the data matrix $\boldsymbol{X}$ must hold for the residual
matrix after subtracting the effect of previous singular values and vectors.
We assume that this is the case. This Assumption is indeed true for the
situations where the true distributions of the random variable $X_{ij}$
denoting the elements of the data matrix $\boldsymbol{X}$ belong to a location
family of distributions with location parameters
$\sum_{r=1}^{k}\lambda_{r}a_{ir}b_{jr}$, where the vectors
$\boldsymbol{a}_{1},\ldots\boldsymbol{a}_{r}$ and
$\boldsymbol{b}_{1},\ldots\boldsymbol{b}_{r}$ are of unit $L_{2}$ norm and
$\lambda_{r}$ denotes the true singular values of the matrix.

Before proceeding, we shall define a few notations used throughout this
section. Any scalar deterministic quantity is denoted by lowercase letters and
the random variables will be denoted by uppercase alphabets. The lowercase
bold symbols $\boldsymbol{v}$ will denote a column vector and the bold
uppercase letters (e.g., $\boldsymbol{M}$) denote matrices. The transpose of a
matrix $\boldsymbol{M}$ will be denoted by $\boldsymbol{M}^{\intercal}$. The
notation $\mathcal{N}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ will be used
to denote a $p$-variate normal distribution with mean vector
$\boldsymbol{\mu}$ and dispersion matrix $\boldsymbol{\Sigma}$. When there is
no ambiguity, we will drop the suffix $p$ for notational convenience. For any
$t=0,1,2,\dots$, we shall use
$\boldsymbol{\theta}^{(t)}=(\lambda^{(t)},\\{a_{i}^{(t)}\\}_{i=1}^{n},\\{b_{j}^{(t)}\\}_{j=1}^{p},(\sigma^{(t)})^{2})$
to denote the rSVDdpd estimator at $t$-th step of the iteration and
$\boldsymbol{\theta}^{\ast}=(\lambda^{\ast},\\{a_{i}^{\ast}\\}_{i=1}^{n},\\{b_{j}^{\ast}\\}_{j=1}^{p},(\sigma^{\ast})^{2})$
to denote the converged rSVDdpd estimator using the iteration rules
(2.5)-(2.7). The population counterpart
$\boldsymbol{\theta}^{g}=(\lambda^{g},\\{a_{i}^{g}\\}_{i=1}^{n},\\{b_{j}^{g}\\}_{j=1}^{p},(\sigma^{g})^{2})$
will denote the ultimate true value of the parameters that are being
estimated. As the minimization in (2.5)-(2.7) defines the rSVDdpd estimator,
the true value $\boldsymbol{\theta}^{g}$ is also expected to satisfy such
minimization criteria, in the sense of overall population based measures
rather than its empirical counterparts. With this in mind, we start by
defining these “best” fitting parameters for the particular setup of SVD under
consideration.

######  Definition 3.1.

Let, $\boldsymbol{X}$ be a data matrix of order $n\times p$ such that the
random variable $X_{ij}$, corresponding to its $(i,j)$-th entry $x_{ij}$,
follows a distribution with density function $g_{ij}$, for all $i=1,\ldots n$
and $j=1,\ldots p$ and $X_{ij}$s are independent to each other. Then,
$\boldsymbol{\theta}^{g}$ is called a “best” fitting parameter if the
following conditions hold

  1. 1.

$\sum_{i=1}^{n}(a_{i}^{g})^{2}=1\qquad\text{ and
}\qquad\sum_{j=1}^{p}(b_{j}^{g})^{2}=1.$ (3.1)

  2. 2.

For any $i=1,2,\dots n,\ j=1,2,\dots p$,

$\lambda^{g}a_{i}^{g}=\operatorname*{arg\,min}_{a}\int
V(x;a,b_{j}^{g},(\sigma^{g})^{2})g_{ij}(x)dx.$ (3.2)

  3. 3.

For any $i=1,2,\dots n,\ j=1,2,\dots p$,

$\lambda^{g}b_{j}^{g}=\operatorname*{arg\,min}_{b}\int
V(x;a_{i}^{g},b,(\sigma^{g})^{2})g_{ij}(x)dx.$ (3.3)

  4. 4.

For any $i=1,2,\dots n,\ j=1,2,\dots p$,

$(\sigma^{g})^{2}=\operatorname*{arg\,min}_{\sigma^{2}}\int
V(x;\lambda^{g}a_{i}^{g},b_{j}^{g},\sigma^{2})g_{ij}(x)dx=\operatorname*{arg\,min}_{\sigma^{2}}\int
V(x;a_{i}^{g},\lambda^{g}b_{j}^{g},\sigma^{2})g_{ij}(x)dx.$ (3.4)

where $V(\cdot)$ is defined as in (2.4).

Here, (3.2) shows that the minimizer of the quantity on the right-hand side of
the equation is always $\lambda^{g}a_{i}^{g}$, independent of the choice of
column index $j$. This Assumption holds if the true densities $g_{ij}$ are
actually densities of the normal distributions with the location parameters
being elements from the best rank one approximation of $\boldsymbol{X}$, i.e.
the entries of the data matrix $X_{ij}$s are normally distributed with mean
$\mu_{ij}$ and constant variance $\sigma^{2}$, and the matrix
$\boldsymbol{\mu}=\left(\mu_{ij}\right)_{i=1,j=1}^{n,p}$ is of rank $1$.

###  3.1 Uniqueness

Since the aim of the rSVDdpd algorithm is to robustly estimate the singular
values and the singular vectors of a given data matrix, it is required to show
that the “best” fitting parameters as introduced by Definition 3.1 resemble
the behaviour of the usual singular values and vectors. Regarding this, the
following theorem claims that if the elements of the data matrix
$\boldsymbol{X}$ are normally distributed, the “best” fitting parameter given
by Definition 3.1 matches with the usual singular values and vectors.

######  Theorem 3.2.

Let the data matrix $\boldsymbol{X}$ be such that
$X_{ij}\sim\mathcal{N}\left(\lambda^{\ast}a_{i}^{\ast}b_{j}^{\ast},(\sigma^{\ast})^{2}\right)$,
and the elements $X_{ij}$’s are independently distributed. Assume that
$\sum_{i}(a_{i}^{\ast})^{2}=\sum_{j}(b_{j}^{\ast})^{2}=1$. Then there exists a
best fitting parameter given by
$\boldsymbol{\theta}^{g}=(\lambda^{\ast},\\{a_{i}^{\ast}\\}_{i=1}^{n},\\{b_{j}^{\ast}\\}_{j=1}^{p},(\sigma^{\ast})^{2})$.
Moreover, this parameter is unique up to the magnitude of each of its
coordinates.

The uniqueness provided in Theorem 3.2 is only upto the magnitude of each
element of the parameter vector $\boldsymbol{\theta}$, hence if both
$\boldsymbol{\theta}_{1}=(\lambda_{1},\\{a_{i1}\\}_{i=1}^{n},\\{b_{j1}\\}_{j=1}^{p},\sigma_{1}^{2})$
and
$\boldsymbol{\theta}_{2}=(\lambda_{2},\\{a_{i2}\\}_{i=1}^{n},\\{b_{j2}\\}_{j=1}^{p},\sigma_{2}^{2})$
are best fitting parameters for the normal setup given in Theorem 3.2, then we
must have
$\left\lvert\lambda_{1}\right\rvert=\left\lvert\lambda_{2}\right\rvert,\left\lvert
a_{i1}\right\rvert=\left\lvert a_{i2}\right\rvert,\left\lvert
b_{j1}\right\rvert=\left\lvert b_{j2}\right\rvert$ and
$\sigma_{1}^{2}=\sigma_{2}^{2}$. Hence, we shall restrict our attention to
only the subset of the parameter space $\boldsymbol{\Theta}$ where each of
these coordinates is positive, so that the uniqueness of the “best” fitting
parameter can be asserted. Thus, we redefine our parameter space as
$\boldsymbol{\Theta}=[0,\infty)\times S_{n}^{+}\times
S_{p}^{+}\times[0,\infty)$, where $S_{n}^{+}$ and $S_{p}^{+}$ are respectively
the parts of the $n$ and $p$-dimensional unit hyperspheres centered at the
origin and lying in the positive orthant, i.e.

$S_{r}^{+}=\left\\{(x_{1},\ldots x_{r}):\sum_{i=1}^{r}x_{i}^{2}=1,x_{i}\geq
0,i=1,\ldots r\right\\}$ (3.5)

An immediate corollary of Theorem 3.2 for a deterministic data matrix
$\boldsymbol{X}$ follows from the observation that degenerate distribution is
one special case of the family of normal distributions with the variance
parameter being equal to $0$.

######  Corollary 1.

If the data matrix $\boldsymbol{X}$ is of rank $1$ such that
$X_{ij}=\lambda^{\ast}a_{i}^{\ast}b_{j}^{\ast}$ for all $i=1,2,\dots n,\
j=1,2,\dots p$, with
$\sum_{i}(a_{i}^{\ast})^{2}=\sum_{j}(b_{j}^{\ast})^{2}=1$, then there exists a
best fitting parameter given by
$\boldsymbol{\theta}^{g}=(\lambda^{\ast},\\{a_{i}^{\ast}\\}_{i=1}^{n},\\{b_{j}^{\ast}\\}_{j=1}^{p},0)$.
Moreover, this parameter is unique up to the magnitude of each of its
coordinates.

In the deterministic setup given in Corollary 1, since the true distribution
is a degenerate distribution at
$x_{ij}=\lambda^{\ast}a_{i}^{\ast}b_{j}^{\ast}$, it follows that $\int
V(x;c,d,\sigma^{2})g_{ij}(x)dx=V(\lambda^{\ast}a_{i}^{\ast}b_{j}^{\ast};c,d,\sigma^{2})$.
Therefore, the population version of the alternating estimating equations
given in (3.2)-(3.4) becomes the same as the iteration formulas given in
(2.5)-(2.7). In other words,
$\boldsymbol{\theta}^{g}=(\lambda^{\ast},\\{a_{i}^{\ast}\\}_{i=1}^{n},\\{b_{j}^{\ast}\\}_{j=1}^{p},0)$
becomes the unique fixed point of the alternating iteration rules, in the
restricted parameter space $\boldsymbol{\Theta}=[0,\infty)\times
S_{n}^{+}\times S_{p}^{+}\times[0,\infty)$. For fixed $n$ and $p$, if we
assume the convergence of the iteration rules (which can be assured later
using Theorem 3.7), then it follows that the converged estimator will be a
fixed point of the iteration rules (2.5)-(2.7). Since the fixed point is
unique, hence the converged estimator should also be unique in this particular
case.

###  3.2 Equivariance Properties

An equivariance property satisfied by the usual singular value decomposition
is that whenever the data matrix is multiplied by some scalar quantity, the
singular values are also multiplied by the same scalar. However, the singular
vectors remain unchanged. The following theorem presents this equivariance
property for a “best” fitting parameter.

######  Theorem 3.3.

If a “best” fitting parameter for matrix $\boldsymbol{X}$ is
$\boldsymbol{\theta}^{g}$, then a “best” fitting parameter for the matrix
$c\boldsymbol{X}$ is
$\widetilde{\boldsymbol{\theta}}^{g}=(c\lambda^{g},\\{a_{i}^{g}\\}_{i=1}^{n},\\{b_{j}^{g}\\}_{j=1}^{p},(c\sigma^{g})^{2})$
for any real constant $c$.

Another equivariance property of the usual singular value decomposition is
that under any row (or column) permutation of the data matrix, the entries of
the left (or right) singular vectors permutes accordingly, while the singular
values remain unaffected. Such a property also holds for a “best” fitting
parameter. Let $\pi_{R}$ and $\pi_{C}$ denote some such permutation on the row
and column indices of the matrix respectively.

######  Theorem 3.4.

If a “best” fitting parameter for the matrix $\boldsymbol{X}$ is
$\boldsymbol{\theta}^{g}$, then a “best” fitting parameter for matrix
$\boldsymbol{P}\boldsymbol{X}\boldsymbol{Q}^{\intercal}$ where
$\boldsymbol{P},\boldsymbol{Q}$ are permutation matrices corresponding to the
permutations $\pi_{R}$ and $\pi_{C}$, is
$\widetilde{\boldsymbol{\theta}}^{g}=(\lambda^{g},\\{a_{\pi_{R}(i)}^{g}\\}_{i=1}^{n},\\{b_{\pi_{C}(j)}^{g}\\}_{j=1}^{p},(\sigma^{g})^{2})$.

Theorem 3.3 and Theorem 3.4 provide an indication of similarity between the
usual singular value and a “best” fitting parameter defined in Definition 3.1.
In other words, a “best” fitting parameter can also be regarded as a
representative of singular value and singular vectors, as both share the same
equivariance properties. Also, in view of Theorem 3.2, the “best” fitting
parameter turns out to be exactly the same as the singular value for normal
distribution. In addition to this, the converged estimator, as obtained by
rSVDdpd algorithm, is also expected to obey such nice equivariance properties.

######  Theorem 3.5.

Let $\boldsymbol{\theta}^{\ast}$ be the converged rSVDdpd estimator for the
matrix $\boldsymbol{X}$, starting from an initial estimate
$\boldsymbol{\theta}^{(0)}$. Then, starting with the new initial estimate
$(c\lambda^{(0)},\\{a_{i}^{(0)}\\}_{i=1}^{n},\\{b_{j}^{(0)}\\}_{j=1}^{p},(c\sigma^{(0)})^{2})$,
the rSVDdpd estimator for the data matrix
$\boldsymbol{X}^{\prime}=c\boldsymbol{X}$ converges to
$(c\lambda^{\ast},\\{a_{i}^{\ast}\\}_{i=1}^{n},\\{b_{j}^{\ast}\\}_{j=1}^{p},(c\sigma^{\ast})^{2})$.

######  Theorem 3.6.

Let $\boldsymbol{\theta}^{\ast}$ be the converged rSVDdpd estimator for the
matrix $\boldsymbol{X}$, starting from an initial estimate
$\boldsymbol{\theta}^{(0)}$. Also, let $\boldsymbol{P}$ and $\boldsymbol{Q}$
be the permutation matrices corresponding to the permutations $\pi_{R}$ and
$\pi_{C}$ respectively. Then, starting with the new initial estimate
$(\lambda^{(0)},\\{a_{\pi_{R}(i)}^{(0)}\\}_{i=1}^{n},\\{b_{\pi_{C}(j)}^{(0)}\\}_{j=1}^{p},(\sigma^{(0)})^{2})$,
the rSVDdpd estimator for the data matrix $\boldsymbol{PXQ}^{\intercal}$
converges to
$(\lambda^{\ast},\\{a_{\pi_{R}(i)}^{\ast}\\}_{i=1}^{n},\\{b_{\pi_{C}(j)}^{\ast}\\}_{j=1}^{p},(\sigma^{\ast})^{2})$.

###  3.3 Convergence and Consistency of the rSVDdpd estimator

The iteration formulas given in (2.5)-(2.7) are essentially obtained from the
particular need to solve the empirical versions of the minimization problems
shown in (3.2), (3.3) and (3.4) subject to the most current estimate of those
parameters in an alternating fashion. In connection to these iteration rules
of the rSVDdpd method in (2.5)-(2.7), and a “best” fitting parameter
$\boldsymbol{\theta}^{g}$ as in Definition 3.1, it is necessary to know if the
iterations converge to some estimate, and whether that estimate is close to a
“best” fitting parameter.

The following theorem answers the first question about the convergence of the
iteration rules of the rSVDdpd algorithm.

######  Theorem 3.7.

For fixed $n$ and $p$, the sequence of estimates
$\\{\boldsymbol{\theta}^{(t)}\\}_{t=0}^{\infty}$ obtained through (2.5)-(2.7)
converges to a minimizer $\boldsymbol{\theta}^{\ast}$ of $\mathcal{H}_{n,p}$,
where

$\mathcal{H}_{n,p}(\boldsymbol{\theta})=\mathcal{H}_{n,p}(\lambda,\\{a_{i}\\}_{i=1}^{n},\\{b_{j}\\}_{j=1}^{p},\sigma^{2})=\frac{1}{np}\sum_{i,j}V(x_{ij};\lambda
a_{i},b_{j},\sigma^{2})=\frac{1}{np}\sum_{i,j}V(x_{ij};a_{i},\lambda
b_{j},\sigma^{2}).$

Answering the second question about the consistency of the converged rSVDdpd
estimate $\boldsymbol{\theta}^{\ast}$ (whose existence is given by Theorem
3.7), has two primary difficulties. Firstly, the dimension of the parameter
space $\boldsymbol{\Theta}$ is $(n+p+2)$, which grows linearly as the
dimension of the matrix ($n$ or $p$) increases to infinity. Secondly, the
parameter space $\boldsymbol{\Theta}$ is not necessarily convex due to the
presence of the coordinates related to singular vectors. The second problem
can be circumvented using an inverse stereographic projection which transforms
this non-convex parameter space $\boldsymbol{\Theta}$ into a convex parameter
space $\boldsymbol{\Xi}\subseteq\mathbb{R}^{(n+p)}$. We call this parameter
space $\boldsymbol{\Xi}$ as the natural parameter space in the given setup.
The one-one transformation $\mathcal{T}$ between these two parameter spaces
$\boldsymbol{\Theta}$ and $\boldsymbol{\Xi}$ are governed by the following two
equations

$\mathcal{T}(\lambda,\\{a_{i}\\}_{i=1}^{n},\\{b_{j}\\}_{j=1}^{p},\sigma^{2})=\left(\lambda,\left\\{\frac{a_{i}}{(1-a_{n})}\right\\}_{i=1}^{(n-1)},\left\\{\frac{b_{j}}{(1-b_{n})}\right\\}_{j=1}^{(p-1)},\sigma^{2}\right),$
(3.6)

and,

$\mathcal{T}^{(-1)}\left(\lambda,\\{\alpha_{i}\\}_{i=1}^{(n-1)},\\{\beta_{j}\\}_{j=1}^{(p-1)},\sigma^{2}\right)\\\
=\left(\lambda,\left\\{\frac{2\alpha_{i}}{U^{2}+1}\right\\}_{i=1}^{(n-1)},\frac{U^{2}-1}{U^{2}+1},\left\\{\frac{2\beta_{j}}{V^{2}+1}\right\\}_{j=1}^{(p-1)},\frac{V^{2}-1}{V^{2}+1},\sigma^{2}\right),$
(3.7)

where $U^{2}=\sum_{i=1}^{(n-1)}\alpha_{i}^{2}$ and
$V^{2}=\sum_{j=1}^{(p-1)}\beta_{j}^{2}$. Accordingly, we denote
$\boldsymbol{\eta}$ as an element of this natural parameter space
$\boldsymbol{\Xi}$, where the corresponding transformed parameter
$\boldsymbol{\theta}=\mathcal{T}^{(-1)}(\boldsymbol{\eta})$ denotes an element
of $\boldsymbol{\Theta}$.

Now that the necessary foundations are laid out, we can present the
consistency theorem, which claims that under some reasonable Assumptions, the
minimizer $\boldsymbol{\theta}^{\ast}$ of $\mathcal{H}_{n,p}$ as given in
Theorem 3.7, is a consistent estimator of the best fitting parameter
$\boldsymbol{\theta}^{g}$. The uniqueness of the best fitting parameter in
this setup is a consequence of Theorem 3.2 following from the Assumption of
normality as the true data generating distribution. However, note that the
description of a “best” fitting parameter indicated in Definition 3.1 is
applicable for fixed $n$ and $p$. In contrast, statistical consistency is an
asymptotic property, which requires the dimensions of the matrix $n$ and $p$
to tend towards infinity. To resolve this conflict in a unified setup, we
assume that a sequence of “best” fitting parameters exists for each fixed $n$
and $p$. We denote the $(i,j)$-th entry of the data matrix
$\boldsymbol{X}_{n,p}$ of order $n\times p$ by the random variable
$(X_{ij})_{n,p}$, and make several Assumptions about it as follows

  1. (A1)

There exists a sequence
$\boldsymbol{\theta}^{g}_{n,p}=(\lambda^{g},\\{a_{i,n}^{g}\\}_{i=1}^{n},\\{b_{j,p}^{g}\\}_{j=1}^{p},(\sigma^{g}_{n,p})^{2})$
of the best fitting parameters. In other words, for any fixed positive
integers $n$ and $p$, the random variables $(X_{ij})_{n,p}$ are independent
for $i=1,\dots n;j=1,\dots,p$, and
$X_{ij}\sim\mathcal{N}\left(\lambda^{g}a_{i,n}^{g}b_{j,p}^{g},(\sigma^{g}_{n,p})^{2}\right)$.

  2. (A2)

For each pair of positive integers $n$ and $p$, there exists an open rectangle
$S_{n,p}$ inside $\boldsymbol{\Xi}_{n,p}\subset\mathbb{R}^{(n+p)}$ containing
the natural parametrization of $\boldsymbol{\theta}^{g}_{n,p}$ i.e.
$\boldsymbol{\eta}^{g}_{n,p}$ such that $\mathcal{\tau}(S_{n,p})$ does not
contain any points with $\boldsymbol{a}=(0,\dots,0,1)$ or
$\boldsymbol{b}=(0,\dots,0,1)$.

  3. (A3)

The converged rSVDdpd estimate for the data matrix $\boldsymbol{X}_{n,p}$,
i.e, $\boldsymbol{\theta}_{n,p}^{\ast}$ (the minimizer of $\mathcal{H}_{n,p}$
as indicated by Theorem 3.7) satisfy
$(\boldsymbol{a}_{n}^{\ast})^{\intercal}\neq(0,\dots,0,1)$ and
$(\boldsymbol{b}_{p}^{\ast})^{\intercal}\neq(0,\dots,0,1)$ for all
sufficiently large $n$ and $p$.

  4. (A4)

The error variances in best fitting parameters satisfy
$(\sigma^{g}_{n,p})^{2}=\Theta\left((np)^{-2/(\alpha+2)}\right)$ for some
$0\leq\alpha\leq 1$. This means, there exists positive finite constants
$\gamma_{1},\gamma_{2}$ such that for sufficiently large $n$ and $p$,
$\gamma_{1}(np)^{-1/(\alpha+2)}<\sigma^{g}_{n,p}<\gamma_{2}(np)^{-1/(\alpha+2)}$.

  5. (A5)

There exists a positive finite constant $M$ independent of $n$ and $p$ such
that

$\mathbb{P}\left(\max_{1\leq i\leq n,1\leq j\leq
p}\left|\left(X_{ij}\right)_{n,p}\right|<M\right)\rightarrow 1$

as $n$ and $p$ both tend to infinity, i.e., the entries of the matrix
$\boldsymbol{X}$ remain uniformly bounded in probability.

The first Assumption (A1) is simply a description of the setup. This
Assumption allows one to perform differentiation under the integral sign and
allows the interchange of differentiation and expectation operators, which is
quite common. Assumptions (A2) and (A3) ensure that the natural
parametrization based on the stereographic projection described above remains
valid. Assumption (A4) is required to provide concentration bound on the tail
probabilities, which is not connected to the usual Assumptions of MDPD
estimation (Ghosh and Basu, 2013). It is well known from random matrix theory
that the Gaussian ensemble with each entry following a standard normal
distribution has the singular values asymptotically at the order of
$(\sqrt{n}+\sqrt{p})$ (Mehta, 2004; Tracy and Widom, 1993). However, since
Assumption (A1) indicates that the same $\lambda^{g}$ acts as the singular
value for any $n$ and $p$, the variance of the entries of the data matrix has
to go down asymptotically to ensure that the singular value does not grow with
increasing $n$ or $p$. Finally, Assumption (A5) can be enforced in the
algorithm itself by scaling the entries of the data matrix which is also
indicated by Roy et al. (2021) for the choice of the initial estimate
$\boldsymbol{\theta}^{(0)}$.

######  Theorem 3.8.

Under Assumptions (A1)-(A5), the rSVDdpd estimate
$\boldsymbol{\theta}_{n,p}^{\ast}$ corresponding to the same $\alpha$ from
Assumption (A4), is a consistent estimator of the sequence of the “best”
fitting parameter values $\boldsymbol{\theta}_{n,p}^{g}$, as both dimensions
of the data matrix $\boldsymbol{X}_{n,p}$ (i.e. $n$ and $p$) tend to infinity.
In other words,
$\left\|\boldsymbol{\theta}_{n,p}^{\ast}-\boldsymbol{\theta}_{n,p}^{g}\right\|\rightarrow
0$ in probability as both $n\rightarrow\infty$ and $p\rightarrow\infty$
subject to a constant ratio in limit, i.e.
$\lim_{\begin{subarray}{c}n\rightarrow\infty\\\
p\rightarrow\infty\end{subarray}}\frac{n}{p}=c$ for some $0<c<\infty$.

######  Remark 3.1.

Since Roy et al. (2021) derived rSVDdpd as an extension of MDPDE in linear
regression setup based on the work of Ghosh and Basu (2013), it is natural to
think that the result on the consistency of the MDPDE given in Theorem 3.1 of
the same paper can be imitated to deliver a proof of the consistency of the
proposed rSVDdpd estimators. However, there are several major complications
involved.

  1. 1.

The basic assumption required for consistency of the MDPDE in the INH
(independent and non-homogeneous) setup considered in Ghosh and Basu (2013) is
the existence of an open set in the parameter space. However, in this
particular setup, the parameter space $\boldsymbol{\Theta}$ by itself does not
contain any open neighbourhood. Therefore, all the necessary formulations are
required to be applied on the natural parametrization
$\boldsymbol{\eta}\in\boldsymbol{\Xi}$ instead, which converts the setup into
a nonlinear regression problem.

  2. 2.

In the case of SVD, the length of the singular vectors is not fixed and grows
with the dimension of the matrix. Thus, as $n,p\rightarrow\infty$, the
dimension of the parameter space also increases to infinity. This problem has
been of considerable interest to many authors under the M-estimation setup
(Huber, 1973; Portnoy, 1984). All of these results assume convexity of the
objective function (He and Shao, 2000), but that cannot be employed in our
case. Here, the objective function is convex in any of the parameters
individually when the other parameters are kept fixed, but becomes non-convex
if all the parameters are taken together.

  3. 3.

In each of the alternating iterations, the consistency ensures that the
resulting estimator based on the minimization of that particular iteration is
probabilistically close to the minimizer of the population version of the
iterative equation. However, the population version of the iterative equation
depends on the current estimates of the other parameters. Hence, in each of
the iterations, the empirical estimates are allowed to deviate from the
population estimates and these errors sum up as the number of iterations
increases. Hence, we require a bound on the tail of the distribution of such a
sum of errors to ensure the consistency of the rSVDdpd estimator.

Thus, a non-trivial modification of the existing proof technique is required.
Due to its length and complications, the proof of this theorem is deferred to
Appendix B.7.

######  Remark 3.2.

A particularly interesting observation in view of Assumption (A4) is that when
$\alpha=0$, the usual MLE of singular values are consistent if the error
variance decreases at the order of $\frac{1}{np}$. However, for higher values
of $\alpha$, the error variance may decrease at a lower rate, only at the
order of $(np)^{-2/(\alpha+2)}$. This provides an indication that with
increased $\alpha$, the rSVDdpd algorithm is able to tolerate a greater
variation in errors (i.e., a low signal-to-noise ratio).

An immediate corollary of Theorem 3.7 and Theorem 3.8 establishes the
correctness of the proposed rSVDdpd estimator in a probabilistic sense.

######  Corollary 2.

Under the Assumptions (A1)-(A5), when the true data generating distribution is
a normal distribution as mentioned in Theorem 3.8, the rSVDdpd estimator of
singular values and vectors converges as the number of iterations goes to
infinity, and the converged estimator is a consistent estimator of the best
fitting parameter as both the dimensions of the matrix increase to infinity
subject to the restriction $\lim_{n,p\rightarrow\infty}n/p=c$ for some
$0<c<\infty$.

######  Remark 3.3.

Theorem 3.5 and 3.6 are, respectively, the empirical counterparts of Theorem
3.3 and 3.4. In view of Theorem 3.5 and 3.6, the equivariance properties hold
given that the initial values of the iterations of rSVDdpd also satisfy such
equivariance properties. However, under Assumptions (A1)-(A5), the convergence
of the rSVDdpd algorithm is immediate and by the consistency of this converged
estimator, it follows that for large $n$ and $p$, the converged estimator can
be made arbitrarily close to the true “best” fitting parameters. Since these
“best” fitting parameters obey equivariance properties as assured in Theorem
3.3 and 3.4, it follows that for large $n$ and $p$, the converged estimator
will also approximately satisfy these equivariance properties, irrespective of
the equivariance of starting values.

##  4 Numerical Illustrations: Simulation Studies

In this section, we compare the performance of the rSVDdpd estimator with two
existing robust SVD estimators, namely the ones proposed by Hawkins et al.
(2001) and Zhang et al. (2013). Implementation of the robust SVD algorithm
proposed by Hawkins et al. (2001) (to be referred to here as pcaSVD) is
available as an R package pcaMethods (Stacklies et al., 2007), which outputs
all singular values and vectors of the input data matrix. The second algorithm
by Zhang et al. (2013) obtains the first pair of singular vectors based on the
minimization procedure

$(\widehat{\boldsymbol{u}},\widehat{\boldsymbol{v}})=\operatorname*{arg\,min}_{(\boldsymbol{u},\boldsymbol{v})}\left[\rho\left(\dfrac{\boldsymbol{X}-\boldsymbol{u}\boldsymbol{v}^{\intercal}}{\sigma}\right)+\mathcal{P}_{\lambda}(\boldsymbol{u},\boldsymbol{v})\right]$

where $\rho(\cdot)$ is a robust loss function (namely Huber’s loss function)
and $\mathcal{P}_{\lambda}$ is a regularization penalty term to motivate
smoothing in the entries of the singular vectors. For an extensive comparison,
we consider two variants of this algorithm. In one variation, we perform the
minimization with only Huber’s loss function without any penalty term, which
we shall refer to as RobSVD, while in the other variation, we follow the
recommended procedure of minimization with penalty term, which we shall call
RobRSVD. Implementation of both of these variants is available in the R
package RobRSVD (Zhang and Pan, 2013) which is programmed only to output the
first singular value and its corresponding singular vectors. Thus, in order to
compare the performances of the algorithms on an equal footing, we add a
wrapper outside this function to apply the same algorithm on the residual
matrix to output the subsequent singular values (see Hawkins et al. (2001) for
details). Along with these, we also consider an implementation of the usual
$L_{2}$ norm minimization based SVD procedure available in R base package
written using LAPACK Fortran library. For the rSVDdpd algorithm, we use the R
package “rsvddpd” as provided by Roy et al. (2021).

###  4.1 Simulation Setups

To compare the performance of the robust SVD approaches, we employ a Monte
Carlo method. We generate random errors and add them to the true matrix, then
apply each of the robust SVD algorithms and track the estimates obtained from
each sample. The mean squared error (MSE) and bias are computed for each
singular value based on $B=1000$ Monte Carlo samples, serving as the accuracy
measures. The sum of squared biases and the sum of MSE across all singular
values are calculated for each estimator.

For comparing the left and right singular vectors, we consider a dissimilarity
score between two normalized vectors, denoted as
$\text{Diss}(\boldsymbol{u},\boldsymbol{v})$, which is equal to
$1-|\langle\boldsymbol{u},\boldsymbol{v}\rangle|$. Here, $|x|$ represents the
absolute value of $x$, and $\langle\boldsymbol{u},\boldsymbol{v}\rangle$ is
the Euclidean inner product between two vectors $\boldsymbol{u}$ and
$\boldsymbol{v}$. The dissimilarity score is $0$ if
$\boldsymbol{u}=\boldsymbol{v}$ or $\boldsymbol{u}=(-\boldsymbol{v})$, and it
is equal to $1$ if $\boldsymbol{u}$ and $\boldsymbol{v}$ are orthogonal. The
average dissimilarity score between the estimated singular vectors and the
true singular vectors, computed over all Monte Carlo samples, is used as a
performance measure.

To construct the true data matrix $\boldsymbol{X}$ for a given singular value
decomposition, we use the coefficients of the first three orthogonal
polynomial contrasts of order $10$ and $4$ and arrange them as columns of
matrices $\boldsymbol{U}$ and $\boldsymbol{V}$ respectively. The resulting
data matrix $\boldsymbol{X}$ is then constructed as

$\boldsymbol{X}=\boldsymbol{U}\begin{bmatrix}10&0&0\\\ 0&5&0\\\ 0&0&3\\\
\end{bmatrix}\boldsymbol{V}^{\intercal},$ (4.1)

with the true singular values being $10,5$ and $3$. In each of the resamples,
errors following a pre-specified distribution are added to the entries of
$\boldsymbol{X}$. Based on the chosen distribution, we divide the simulation
scenarios broadly into $5$ categories denoted by (S1)-(S5), as follows.

  1. (S1)

The errors follow the standard normal distribution $\mathcal{N}(0,1)$ (no
outliers per se).

  2. (S2)

The errors are distributed according to a contaminated standard normal
distribution namely

$e_{ij}\sim(1-\epsilon)\mathcal{N}(0,1)+\epsilon\delta_{25},$

where $\epsilon$ is the amount of contamination and $\delta_{25}$ denotes the
degenerate distribution at $25$. Based on the amount of contamination, we
consider three sub-cases of this simulation setting.

    1. (S2a)

$\epsilon=0.05$, denotes only $5\%$ contamination, which corresponds to a
relatively light amount of outlying observations.

    2. (S2b)

$\epsilon=0.1$, denoting medium level contamination with the presence of
$10\%$ outlying values.

    3. (S2c)

$\epsilon=0.2$, denoting heavy contamination with approximately $20\%$ of the
entries being same as the outlying observation of $25$.

  3. (S3)

The errors are distributed according to a standard normal distribution with
block-based contamination as presented in Zhang et al. (2013). Here, in each
resample, the normally distributed errors are first added to each of the
entries of $\boldsymbol{X}$, and then a $2\times 2$ block of submatrix is
chosen randomly and all entries of that submatrix are substituted to $25$.

  4. (S4)

The errors are distributed according to a standard Cauchy distribution. This
setup helps us to study the effect of heavy tailed errors in the robust
estimation of SVD.

  5. (S5)

The errors are distributed according to a standard lognormal distribution,
which is used to study the effect of an asymmetric error distribution with
only positive support.

Table 1 summarizes the comparative results of the usual SVD method, the three
existing robust SVD algorithms (pcaSVD (Stacklies et al., 2007) and two
variants of RobRSVD (Zhang and Pan, 2013)) and the rSVDdpd method (Roy et al.,
2021) for different choices of robustness parameter $\alpha$, based on the
aforementioned performance measures, for different simulation setups
(S1)-(S5).

Table 1: Summary of performance measures (Total squared bias, Total MSE of
singular values and Total dissimilarity, denoted by Diss, of left and right
singular vectors) of different existing SVD and robust SVD algorithms

Simulation Setup | Measure | Existing methods for computing SVD | Choice of $\alpha$ in rSVDdpd
---|---|---|---
Usual SVD | pcaSVD | RobSVD | RobRSVD | $\alpha=0.1$ | $\alpha=0.3$ | $\alpha=0.5$ | $\alpha=0.7$ | $\alpha=1$
S1 | Sq. Bias | 7.957 | 15.066 | 8.808 | 1.684 | 7.952 | 7.94 | 7.925 | 7.894 | 7.764
Total MSE | 10.456 | 24.242 | 11.608 | 6.865 | 10.455 | 10.473 | 10.553 | 10.68 | 10.841
Diss (left) | 0.701 | 1.198 | 0.799 | 0.733 | 0.702 | 0.707 | 0.717 | 0.734 | 0.771
Diss (right) | 0.418 | 0.968 | 0.52 | 0.514 | 0.419 | 0.425 | 0.433 | 0.449 | 0.487
S2a | Sq. Bias | 519.35 | 350.402 | 523.859 | 17.672 | 294.047 | 18.138 | 10.877 | 9.433 | 8.779
Total MSE | 729.83 | 793.612 | 748.043 | 114.802 | 622.819 | 91.894 | 44.016 | 32.329 | 27.144
Diss (left) | 1.93 | 1.679 | 1.933 | 1.324 | 1.613 | 0.944 | 0.896 | 0.885 | 0.904
Diss (right) | 1.529 | 1.286 | 1.516 | 1.074 | 1.163 | 0.615 | 0.576 | 0.571 | 0.59
S2b | Sq. Bias | 1397.36 | 1086.341 | 1427.869 | 60.112 | 1039.977 | 130.472 | 44.154 | 30.039 | 23.262
Total MSE | 1661.82 | 1706.458 | 1722.445 | 278.585 | 1529.374 | 494.34 | 237.68 | 167.244 | 132.886
Diss (left) | 2.155 | 1.965 | 2.148 | 1.672 | 2.013 | 1.283 | 1.113 | 1.084 | 1.078
Diss (right) | 1.702 | 1.5 | 1.665 | 1.371 | 1.5 | 0.899 | 0.756 | 0.733 | 0.731
S2c | Sq. Bias | 2434.021 | 2110.99 | 2488.98 | 193.664 | 2094.11 | 662.604 | 298.008 | 196.514 | 141.064
Total MSE | 2712.69 | 2782.575 | 2803.595 | 617.094 | 2587.792 | 1453.404 | 938.705 | 739.152 | 612.155
Diss (left) | 2.206 | 2.11 | 2.203 | 1.932 | 2.17 | 1.666 | 1.444 | 1.363 | 1.313
Diss (right) | 1.73 | 1.633 | 1.688 | 1.598 | 1.651 | 1.227 | 1.042 | 0.965 | 0.93
S3 | Sq. Bias | 1677.949 | 1640.708 | 1679.881 | 1090.284 | 1355.176 | 114.714 | 28.128 | 21.594 | 18.048
Total MSE | 1686.361 | 1654.5 | 1688.708 | 1248.276 | 1634.409 | 458.468 | 156.986 | 113.33 | 96.675
Diss (left) | 2.052 | 2.003 | 1.941 | 2.162 | 2.007 | 1.208 | 1.024 | 1.012 | 1.002
Diss (right) | 1.924 | 1.836 | 1.832 | 2.175 | 1.844 | 0.895 | 0.692 | 0.678 | 0.667
S4 | Sq. Bias | 41825.788 | 14224.145 | 41779.81 | 540.799 | 469.522 | 265.135 | 198.809 | 169.082 | 140.645
Total MSE | 2171697.037 | 2163377.363 | 2171711.033 | 29149.731 | 1629.485 | 1197.881 | 859.039 | 807.869 | 842.77
Diss (left) | 2.089 | 1.989 | 2.095 | 1.707 | 1.97 | 1.877 | 1.838 | 1.809 | 1.786
Diss (right) | 1.603 | 1.489 | 1.602 | 1.303 | 1.496 | 1.403 | 1.367 | 1.355 | 1.324
S5 | Sq. Bias | 93.633 | 77.499 | 94.901 | 29.135 | 69.98 | 67.092 | 62.429 | 56.64 | 49.424
Total MSE | 146.187 | 108.514 | 149.13 | 49.046 | 85.102 | 83.272 | 78.83 | 72.46 | 65.682
Diss (left) | 2.02 | 2.034 | 1.969 | 1.967 | 1.968 | 1.954 | 1.943 | 1.939 | 1.934
Diss (right) | 1.785 | 1.81 | 1.734 | 1.824 | 1.754 | 1.744 | 1.734 | 1.731 | 1.726

As shown in Table 1, the usual SVD generally leads to a biased estimator of
the singular values for Gaussian errors which is also supported in well-
established theory (Rudelson and Vershynin, 2010). For the simulation setup
(S1), the pcaSVD algorithm is found to be the most biased, followed by RobSVD,
both of which have more bias and MSE than the usual SVD algorithm. Compared to
the usual SVD, rSVDdpd algorithm achieves lesser bias as the robustness
parameter $\alpha$ increases, but at the cost of higher variance and MSE.
RobRSVD achieves the minimum bias and MSE in this scenario, but it shows a
higher variance and a higher dissimilarity score in singular vectors than the
rSVDdpd algorithm.

Turning our attention to simulation setups (S2)(S2a), (S2)(S2b) and (S2)(S2c),
we see that the usual SVD and existing robust SVD algorithms pcaSVD and RobSVD
do not yield very reliable estimates of the singular values. Although RobRSVD
provides reasonable estimates, rSVDdpd achieves lower bias and MSE for some
choices of $\alpha$. In the presence of random outlying observations, as in
the case of simulation setups (S2)(S2a), (S2)(S2b) and (S2)(S2c), both the
bias and MSE for rSVDdpd show reductions as the robustness parameter $\alpha$
is increased from $0$ to $1$. The dissimilarities of singular vectors also
tend to decrease with an increase in $\alpha$.

For the block level contamination in simulation setup (S3), we find that
rSVDdpd has much better performances than the other robust SVD algorithms for
all performance metrics. With errors from a heavy-tailed distribution as
considered in the simulation setup (S4), the results remain very similar. The
rSVDdpd algorithm provides the least bias and MSE, and even with a small
robustness parameter $\alpha=0.1$, rSVDdpd outperforms the existing robust SVD
algorithms under consideration.

In simulation setup (S5) with lognormally distributed errors having positive
support, rSVDdpd outperforms the usual SVD, pcaSVD and RobSVD methods by
showing a reduction in both bias and MSE. However, as in the simulation setup
(S1), RobRSVD is again found to provide estimates with the least bias and MSE,
but at a cost of higher variance and dissimilarity scores than rSVDdpd.

Although RobRSVD outputs better singular value estimates than rSVDdpd under
normally and lognormally distributed errors, it does so at the cost of
extremely high computational complexity. This is precisely due to the matrix
inversion step to compute
$(\mathcal{V}^{\intercal}\mathcal{W}^{\ast}\mathcal{V}+2\Omega_{\boldsymbol{u}\mid\boldsymbol{v}})^{-1}$
(see Eq. (9) of Zhang et al. (2013)). Since the best known matrix inversion
algorithm, i.e., a variant of Coppersmith-Winograd algorithm (Alman and
Williams, 2021) achieves a computational complexity of
$\mathcal{O}(n^{2.3728596})$ for inverting an $n\times n$ matrix, it follows
that each iteration of the RobRSVD algorithm has time complexity
$\mathcal{O}(n^{2.3728596}+p^{2.3728596})$. On the other hand, each iteration
of rSVDdpd performs only a weighted average computation, which reduces its
computational budget to $\mathcal{O}(n^{2}+p^{2})$. To demonstrate this, we
consider $n\times p$ matrices with uniformly distributed entries for different
choices of $n$ and fixed $p=25$, and apply different methods of computing SVD
on them. Table 2 summarizes the time taken (in units of milliseconds) to
obtain the first singular value from different algorithms for different
choices of $n$, in a computer with Intel i5-8300H 2.30GHz processor with 8 GB
of RAM. As seen from Table 2, the computational budget of rSVDdpd is similar
to pcaSVD, which is lower by several orders of magnitude than RobSVD and
RobRSVD. This extremely high computational cost of RobRSVD can be circumvented
if the matrix
$(\mathcal{V}^{\intercal}\mathcal{W}^{\ast}\mathcal{V}+2\Omega_{\boldsymbol{u}\mid\boldsymbol{v}})$
becomes a diagonal matrix, which happens if the penalty parameter is taken as
zero and RobRSVD is reduced to its non-regularized variant RobSVD. However, as
Table 1 shows, the RobSVD algorithm without the regularization cannot provide
a reliable robust estimate of singular values, even using Huber’s robust loss
function.

Table 2: Summary of average time taken (in milliseconds) to obtain the first
singular value and vectors of an $n\times 25$ matrix with random entries from
$U(0,1)$ via different SVD algorithms

Number of rows ($n$) | Existing methods for computing SVD | Choice of $\alpha$ in rSVDdpd
---|---|---
Usual SVD | pcaSVD | RobSVD | RobRSVD | $\alpha=0.1$ | $\alpha=0.5$ | $\alpha=1$
5 | 0.014 | 2.209 | 3.098 | 73.666 | 0.545 | 0.841 | 1.417
10 | 0.011 | 4.388 | 3.966 | 99.077 | 0.975 | 1.508 | 2.100
25 | 0.008 | 4.965 | 3.989 | 81.045 | 3.861 | 6.416 | 11.394
50 | 0.013 | 8.519 | 7.824 | 149.4 | 4.396 | 7.104 | 12.526
100 | 0.017 | 15.696 | 34.649 | 826.44 | 5.136 | 8.683 | 14.362
250 | 0.026 | 32.066 | 402.125 | 7839.942 | 7.756 | 13.252 | 22.379
500 | 0.041 | 35.828 | 2948.001 | 54209.15 | 13.622 | 21.413 | 32.404
750 | 0.058 | 58.697 | 10363.441 | 210494.564 | 24.67 | 36.563 | 55.422
1000 | 0.072 | 69.76 | 26282.893 | 531362.110 | 27.727 | 40.309 | 62.234

##  5 Conclusion

As depicted in Section 1, a plethora of algorithms from an extensive range of
disciplines use singular value decomposition as a core component of the
methods. However, the increasing prevalence of big data has made it
challenging to ensure the accuracy and reliability of the data. The input data
for these algorithms are prone to contamination by noise and outliers, leading
to inaccurate results when using standard SVD. To address this issue, several
robust SVD methods have been proposed (see Section 1.1), but most of them are
not scalable to large matrices encountered in real-life applications. The lack
of theoretical guarantees of these algorithms has limited their widespread
adoption and hindered their application in critical domains. While Roy et al.
(2021) demonstrate an application of the “rSVDdpd” algorithm to solve a real-
life problem, in this paper, we provide the theoretical justification for its
reliability. The simulation results further validate the superiority of
rSVDdpd compared to the existing algorithms. However, more investigation is
needed to develop asymptotic distributions of the estimated singular values
and vectors, which can provide confidence interval estimates.

We believe that the “rSVDdpd” algorithm has potential applications beyond
video surveillance background modelling by using it as a replacement of the
standard SVD in various algorithms to handle data contamination. As an
example, we can use rSVDdpd for modelling genetic data, performing community
detection in networks, estimating latent semantic representation of text
documents from term-document matrices, etc. We hope to explore these
applications in future.

##  Appendix A A Brief Review of Minimum Density Power Divergence Estimator

Basu et al. (1998) introduced the density power divergence as a measure of
discrepancy between two probability density functions, which being an
M-estimator, as well as a minimum distance estimator, enjoys various
theoretical properties. The density power divergence between the densities $g$
and $f$ is defined as

$d_{\alpha}(g,f)=\begin{cases}\displaystyle\int\left\\{f^{1+\alpha}-\left(1+\dfrac{1}{\alpha}\right)f^{\alpha}g+\dfrac{1}{\alpha}g^{1+\alpha}\right\\}&\alpha>0\\\
\displaystyle\int g\ln\left(\dfrac{g}{f}\right)&\alpha=0\\\ \end{cases}.$

Here $\ln(\cdot)$ denotes the natural logarithm. The control parameter
$\alpha$ provides a smooth bridge between robustness and efficiency.

In case of independent and identically distributed observations,
$Y_{1},Y_{2},\dots Y_{n}$, with true distribution function $G$ and
corresponding density $g$, we model this unknown density by a parametric
family of densities
$\mathcal{F}_{\boldsymbol{\theta}}=\left\\{f_{\boldsymbol{\theta}}:\boldsymbol{\theta}\in\boldsymbol{\Theta}\right\\}$.
The estimator of $\theta$ is then obtained as

$\widehat{\boldsymbol{\theta}}=\arg\min_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}d_{\alpha}(dG_{n},f_{\boldsymbol{\theta}}),$

where $G_{n}$ is the empirical distribution function. This can be shown to be
equivalent to

$\widehat{\boldsymbol{\theta}}=\arg\min_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\left[\int
f_{\boldsymbol{\theta}}^{1+\alpha}-\left(1+\dfrac{1}{\alpha}\right)\dfrac{1}{n}\sum_{i=1}^{n}f_{\boldsymbol{\theta}}(Y_{i})^{\alpha}\right].$

Later, Ghosh and Basu (2013) extended this work by allowing independent but
not identically distributed data. In this case, the observed data $Y_{i}\sim
g_{i}$, where each $g_{i}$ is an unknown density. Each of the true density
$g_{i}$ is modeled by a corresponding parametric family of densities
$\mathcal{F}_{i,\theta}=\left\\{f_{i,\theta}:\theta\in\boldsymbol{\Theta}\right\\}$
for all $i=1,2,\dots n$. Finally, the proposed MDPD estimator is obtained as

$\widehat{\boldsymbol{\theta}}=\arg\min_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}\dfrac{1}{n}\sum_{i=1}^{n}\left[\int
f_{i,\boldsymbol{\theta}}^{1+\alpha}-\left(1+\dfrac{1}{\alpha}\right)f_{i,\boldsymbol{\theta}}(Y_{i})^{\alpha}\right].$

Various nice theoretical properties like consistency and asymptotic normality
of the above MDPD estimator have been proven by Ghosh and Basu (2013).

##  Appendix B Proofs of the Results

Before proving the individual results, we note that the expression
$V(x;c,d,\sigma^{2})$ in (2.4) can be expressed as

$V(y;c,d,\sigma^{2})=\int\phi\left(\dfrac{x-cd}{\sigma}\right)^{(1+\alpha)}dx-\left(1+\dfrac{1}{\alpha}\right)\phi\left(\dfrac{y-cd}{\sigma}\right)^{\alpha},$
(B.1)

where $\phi$ is the density function of the standard normal distribution. We
shall be using this expression of $V$ in the proof.

###  B.1 Proof of Theorem 3.2

Relation (3.1) is verified by the conditions given for $a_{i}^{\ast}$ and
$b_{j}^{\ast}$’s. To verify (3.2), note that with $b_{j}^{g}=b_{j}^{\ast}$ and
$\sigma^{g}=(\sigma^{\ast})$, the quantity in (3.2) is same as minimizing

$\int\phi\left(\dfrac{x-ab_{j}^{\ast}}{\sigma^{\ast}}\right)^{(1+\alpha)}dx-\left(1+\dfrac{1}{\alpha}\right)\int\phi\left(\dfrac{x-ab_{j}^{\ast}}{\sigma^{\ast}}\right)^{\alpha}g_{ij}(x)dx+\dfrac{1}{\alpha}\int
g_{ij}(x)^{(1+\alpha)}dx,$ (B.2)

as the last term is independent of the minimization over $a$. The quantity
given in (B.2) is the density power divergence between the normal density with
mean $ab_{j}^{\ast}$ and variance $\sigma^{\ast}$ and, the true density
$g_{ij}$ with mean $\lambda^{\ast}a_{i}^{\ast}b_{j}^{\ast}$ and variance
$\sigma^{\ast}$. From Theorem 2.1 of Basu et al. (1998), it follows that this
divergence is minimized if and only if two densities match, i.e.,
$a=\lambda^{\ast}a_{i}^{\ast}$. By exactly similar logic and interchanging the
roles of $a_{i}$ and $b_{j}$, (3.3) and (3.4) can also be verified. This
proves that
$\boldsymbol{\theta}^{\ast}=\left(\lambda^{\ast},\\{a_{i}^{\ast}\\}_{i=1}^{n},\\{b_{j}^{\ast}\\}_{j=1}^{p},(\sigma^{\ast})^{2}\right)$
is a “best” fitting parameter for the given setup.

In order to prove uniqueness, suppose
$\boldsymbol{\theta}=\left(\lambda,\\{a_{i}\\}_{i=1}^{n},\\{b_{j}\\}_{j=1}^{p},\sigma^{2}\right)$
be another “best” fitting parameter. Then again, the quantity

$\int\phi\left(\dfrac{x-ab_{j}}{\sigma}\right)^{(1+\alpha)}dx-\left(1+\dfrac{1}{\alpha}\right)\int\phi\left(\dfrac{x-ab_{j}}{\sigma}\right)^{\alpha}g_{ij}(x)dx+\dfrac{1}{\alpha}\int
g_{ij}(x)^{(1+\alpha)}dx,$ (B.3)

is minimized at $a=\lambda a_{i}$, independently of the choice of $j$.
However, (B.3) can be made equal to its minimum value $0$ if and only if
$\sigma^{2}=(\sigma^{\ast})^{2}$, and $\lambda
a_{i}b_{j}=\lambda^{\ast}a_{i}^{\ast}b_{j}^{\ast}$ (follows from Theorem 2.1
of Basu et al. (1998)). Since, both $\boldsymbol{\theta}$ and
$\boldsymbol{\theta}^{\ast}$ are “best” fitting parameters, they must satisfy
(3.1). Hence, $(\lambda a_{i})^{2}=\sum_{j}(\lambda
a_{i}b_{j})^{2}=\sum_{j}(\lambda^{\ast}a_{i}^{\ast}b_{j}^{\ast})^{2}=\lambda^{\ast}a_{i}^{\ast}$.
Taking sum over the row index $i$ now gives $|\lambda|=|\lambda^{\ast}|$.
Then, it easily follows that $|a_{i}|=|a_{i}^{\ast}|$ and
$|b_{j}|=|b_{j}^{\ast}|$.

###  B.2 Proof of Theorem 3.3

It is obvious that $\widetilde{\boldsymbol{\theta}}^{g}$ satisfies (3.1) as
$\boldsymbol{\theta}^{g}$ is given to be a “best” fitting parameter.
Considering the matrix $\boldsymbol{Y}=c\boldsymbol{X}$, let us denote the
true density of $Y_{ij}$ as $g_{ij}^{Y}(\cdot)$, as opposed to $g_{ij}(\cdot)$
denoting the true density of $X_{ij}$. A change of variable formula yields
that $g_{ij}^{Y}(y)=g_{ij}(y/c)$. Hence, from the substitution principle of
integration, it follows that

$\int V\left(y;a,b_{j}^{g},c^{2}(\sigma^{g})^{2}\right)g_{ij}^{Y}(y)dy=c\int
V\left(z;a/c,b_{j}^{g},(\sigma^{g})^{2}\right)g_{ij}(z)dz.$

Since the right-hand side is minimized at $a/c=\lambda^{g}a_{i}^{g}$, the
left-hand side is minimized at $a=c\lambda^{g}a_{i}^{g}$. This verifies (3.2).
Similar to this, (3.3) can also be established by interchanging the role of
$\boldsymbol{a}$ and $\boldsymbol{b}$ above. For the parameter $\sigma$, it
can be verified that a substitution principal again applies, and we obtain

$\int V(y;c\lambda^{g}a_{i}^{g},b_{j}^{g},\sigma^{2})g_{ij}^{Y}(y)dy=\int
V(y;a_{i}^{g},c\lambda^{g}b_{j}^{g},\sigma^{2})g_{ij}^{Y}(y)dy\\\ =c\int
V(z;\lambda a_{i}^{g},b_{j}^{g},\sigma^{2}/c^{2})g_{ij}(z)dz.$

Again by the hypothesis that $\boldsymbol{\theta}^{g}$ is a “best” fitting
parameter, the later is minimized when $\sigma/c=\sigma^{g}$, hence the former
is minimized at $\sigma^{2}=c^{2}(\sigma^{g})^{2}$. This verifies (3.4).

###  B.3 Proof of Theorem 3.4

Let, $\boldsymbol{Y}=\boldsymbol{PXQ}^{\intercal}$. Then, (3.1) is satisfied
for the new setup as
$\sum_{i=1}^{n}(a_{\pi_{R}(i)}^{g})^{2}=\sum_{i=1}^{n}(a_{i}^{g})^{2}=1$ and
similarly
$\sum_{j=1}^{p}(b_{\pi_{C}(j)}^{g})^{2}=\sum_{j=1}^{p}(b_{j}^{g})^{2}=1$. To
see that (3.2) hold for the new setup with $\widetilde{\theta}^{g}$, note that
for every $j=1,2,\dots p$, the minimizer of the integral in (3.2) is
$a_{i}^{g}$, independent of the choice of $j$. Now, considering (3.2) for
$\pi_{R}^{-1}(i)$ and $\pi_{C}^{-1}(j)$ instead of $i$ and $j$, we obtain

$\lambda^{g}a_{\pi_{R}^{-1}(i)}^{g}=\operatorname*{arg\,min}_{a}\int
V\left(x;a,b_{\pi_{C}^{-1}(j)}^{g},(\sigma^{g})^{2}\right)g_{\pi_{R}^{-1}(i),\pi_{C}^{-1}(j)}(x)dx.$
(B.4)

However, $a^{g}_{\pi_{R}^{-1}(i)}$ is the $i$-th entry of sequence
$\\{a_{\pi_{R}(i)}^{g}:i=1,2,\dots n\\}$, i.e. if we consider a vector
$\boldsymbol{a}^{g}$ with its entries $a_{i}^{g}$, then
$a^{g}_{\pi_{R}^{-1}(i)}$ is the $i$-th entry of $\boldsymbol{Pa}$. Similarly,
$b^{g}_{\pi_{C}^{-1}(j)}$ is the $j$-th entry of $\boldsymbol{Qb}$. And
finally, the elements of the new matrix are
$Y_{ij}=X_{\pi_{R}(i),\pi_{C}(j)}$, thus the density for the element $Y_{ij}$
is $g^{Y}_{ij}(y)=g_{\pi^{-1}_{R}(i),\pi^{-1}_{C}(j)}(y)$ which can be
verified by a change of variable formula. Combining these, (B.4) can be
reformulated as

$\lambda^{g}(\boldsymbol{Pa}^{g})_{i}=\operatorname*{arg\,min}_{a}\int
V\left(x;a,(\boldsymbol{Qb}^{g})_{j},(\sigma^{g})^{2}\right)g^{Y}_{ij}(x)dx.$

This shows that (3.2) holds for new matrix $\boldsymbol{PXQ}^{\intercal}$ with
the given best fitting parameter $\widetilde{\boldsymbol{\theta}}^{g}$. The
relation (3.3) holds by imitating the same proof, except interchanging the
role of $\boldsymbol{a}$ and $\boldsymbol{b}$.

Finally, relation (3.4) for the permuted matrix follows from noting that

$\int
V\left(x;\lambda^{g}a_{\pi_{R}^{-1}(i)}^{g},b_{\pi_{C}^{-1}(j)}^{g},\sigma^{2}\right)g_{\pi_{R}^{-1}(i),\pi_{C}^{-1}(j)}(x)dx.=\int
V\left(x;\lambda^{g}(\boldsymbol{Pa}^{g})_{i},(\boldsymbol{Qb}^{g})_{j},\sigma^{2}\right)g^{Y}_{ij}(x)dx.$

###  B.4 Proof of Theorem 3.5

Let us denote
$\boldsymbol{\theta}^{(t)}=\left(\lambda^{(t)},\\{a_{i}^{(t)}\\}_{i=1}^{n},\\{b_{j}^{(t)}\\}_{j=1}^{p},(\sigma^{(t)})^{2}\right)$
as the estimate at the $t$-th iteration when the data matrix is
$\boldsymbol{X}$. Also, let us denote
$\widetilde{\boldsymbol{\theta}}^{(t)}=\left(\widetilde{\lambda}^{(t)},\\{\widetilde{a}_{i}^{(t)}\\}_{i=1}^{n},\\{\widetilde{b}_{j}^{(t)}\\}_{j=1}^{p},(\widetilde{\sigma}^{(t)})^{2}\right)$
as the estimate in the $t$-th iteration when the data matrix is
$c\boldsymbol{X}$. Clearly, it is then enough to show that,

$\displaystyle\widetilde{\lambda}^{(t)}=c\lambda^{(t)},\
(\widetilde{\sigma}^{(t)})^{2}=c^{2}(\sigma^{(t)})^{2},\
\widetilde{a}_{i}^{(t)}=a_{i}^{(t)},\
\widetilde{b}_{j}^{(t)}=b_{j}^{(t)},\quad i=1,\dots n;\ j=1,\dots p;\
t=1,2,\ldots.$

We will show this by using the principle of mathematical induction. For $t=0$,
the claim is validated by the equivariance of the initial estimate. To show
the inductive step, the estimates at $(t+1)$-th iteration can be expressed as

$\displaystyle\widetilde{\lambda}^{(t+1/2)}\widetilde{a}_{i}^{(t+1)}$
$\displaystyle=\operatorname*{arg\,min}_{a}\frac{1}{p}\sum_{j=1}^{p}V(cx_{ij};a,\widetilde{b}_{j}^{(t)},\widetilde{\sigma}^{(t)})$
$\displaystyle=\operatorname*{arg\,min}_{a}\frac{1}{p}\sum_{j=1}^{p}V(cx_{ij};a,b_{j}^{(t)},c\sigma^{(t)}),\
\text{by induction hypothesis}$
$\displaystyle=\operatorname*{arg\,min}_{a}-\frac{1}{p}\sum_{j=1}^{p}\frac{(1+\alpha)}{\alpha(2\pi)^{\alpha/2}(c\sigma^{(t)})^{\alpha}}\exp\left\\{-\alpha\frac{(cx_{ij}-ab_{j}^{(t)})}{c\sigma^{(t)}}\right\\}\
\text{by~{}\eqref{eqn:v-theta-normal}}$
$\displaystyle=\operatorname*{arg\,min}_{a}\frac{1}{p}\sum_{j=1}^{p}V(x_{ij};a/c,b_{j}^{(t)},\sigma^{(t)})$
$\displaystyle=c\lambda^{(t+1/2)}a_{i}^{(t+1)},$

where the last line follows from the fact that the same iteration formula
yields the minimizer $(a/c)=a_{i}^{(t+1)}$ due to the very definition of
$a_{i}^{(t+1)}$. An identical chain of equality for (2.6) can be used to show
that $\lambda^{(t+1)}\tilde{b}_{j}^{(t+1)}=c\lambda^{(t+1)}b_{j}^{(t+1)}$,
while similar derivation from (2.7) establishes
$(\widetilde{\sigma}^{(t+1)})^{2}=(c\sigma^{(t+1)})^{2}$. Finally, since the
estimates of the singular vectors are normalized and restricted to be in the
parameter space $\boldsymbol{\Theta}=[0,\infty)\times S_{n}^{+}\times
S_{p}^{+}\times[0,\infty)$, the inductive step follows.

###  B.5 Proof of Theorem 3.6

This proof is very similar to the proof of Theorem 3.5. We shall again denote
$\boldsymbol{\theta}^{(t)}$ as the estimate at the $t$-th iteration for the
data matrix $\boldsymbol{X}$ and $\widetilde{\boldsymbol{\theta}}^{(t)}$ as
the estimate at the $t$-th iteration for the data matrix
$\boldsymbol{PXQ}^{\intercal}$. Again, it is enough to show that

$\displaystyle\widetilde{\lambda}^{(t)}=c\lambda^{(t)},\
(\widetilde{\sigma}^{(t)})^{2}=c^{2}(\sigma^{(t)})^{2},\
\widetilde{a}_{i}^{(t)}=a_{i}^{(t)},\
\widetilde{b}_{j}^{(t)}=b_{j}^{(t)},\quad i=1,\dots n;\ j=1,\dots p;\
t=1,2,\ldots,$

which we shall show using the principle of mathematical induction. The initial
case $t=0$ follows from the equivariance of the initial estimate. Starting
with the iteration formula (2.5), we have

$\displaystyle\widetilde{\lambda}^{(t+1/2)}\widetilde{a}_{i}^{(t+1)}$
$\displaystyle=\operatorname*{arg\,min}_{a}\dfrac{1}{p}\sum_{j=1}^{p}V(x_{\pi_{R}(i),\pi_{C}(j)};a,\widetilde{b}_{j}^{(t)},\widetilde{\sigma}^{(t)})$
$\displaystyle=\operatorname*{arg\,min}_{a}\dfrac{1}{p}\sum_{j=1}^{p}V(x_{\pi_{R}(i),\pi_{C}(j)};a,b_{\pi_{C}(j)}^{(t)},\sigma^{(t)}),\
\text{by induction hypothesis}$
$\displaystyle=\operatorname*{arg\,min}_{a}\dfrac{1}{p}\sum_{j=1}^{p}V(x_{\pi_{R}(i),j};a,b_{j}^{(t)},\sigma^{(t)}),\
\text{since, }\pi_{C}\text{ is a permutation}$
$\displaystyle=\lambda^{(t+1/2)}a_{\pi_{R}(i)}^{(t+1)},$

where the last equality follows from the fact that the minimization criteria
satisfied by $a_{\pi_{R}(i)}^{(t+1)}$ is the same as minimization of the
average $V$-values for the corresponding row index $\pi_{R}(i)$ of the data
matrix $\boldsymbol{X}$. The other equations
$\lambda^{(t+1)}\tilde{b}_{j}^{(t+1)}=\lambda^{(t+1)}b_{\pi_{C}(j)}^{(t+1)}$
and $(\widetilde{\sigma}^{(t+1)})^{2}=(\sigma^{(t+1)})^{2}$ follows similarly
from (2.6) and (2.7) respectively. Again, restricting the parameter estimates
in $\boldsymbol{\Theta}=[0,\infty)\times S_{n}^{+}\times
S_{p}^{+}\times[0,\infty)$, the inductive step follows.

###  B.6 Proof of Theorem 3.7

First we show the existence of the minimizer of $\mathcal{H}_{n,p}$ for fixed
$n$ and $p$. Note that, as the component $\lambda$ of $\boldsymbol{\theta}$
increases to infinity, each term $V(x_{ij};\lambda a_{i},b_{j},\sigma^{2})$
decreases to $0$. On the other hand, for $\lambda=0$, each term
$V(x_{ij};0,b_{j})$ is finite and bounded below. Similarly, as
$\sigma\rightarrow\infty$, $\mathcal{H}_{n,p}\rightarrow 0$, and as
$\sigma\rightarrow 0$, each term $V(x_{ij};\lambda a_{i},b_{j},\sigma^{2})$
increases to infinity, hence is also bounded below. Therefore, there exists
compact set $\mathcal{C}=[0,M_{1}]\times S_{n}^{+}\times
S_{p}^{+}\times[\epsilon_{2},M_{2}]$ for some small $\epsilon_{2}>0$ such that
for any $\boldsymbol{\theta}\in\mathcal{C}$ the objective function
$\mathcal{H}_{n,p}(\boldsymbol{\theta})$ is bounded below. Since
$\mathcal{H}_{n,p}$ is continuous in its argument, it follows that it attains
a minimum at some $\boldsymbol{\theta}^{\ast}\in\mathcal{C}$.

Moving on, we shall employ Theorem 10.3.3 of Argyros (2008) to show that our
estimation process actually converges to the desired minimizer. Going by the
notation of Argyros (2008), let $F_{t}$ denotes one step of iteration. Since,
each of the iteration steps aims to identify a minima and there is a
possibility of existence of multiple minima, we have

$\left(\lambda^{(t+1)},\\{a_{i}^{(t+1)}\\}_{i=1}^{n},\\{b_{j}^{(t+1)}\\}_{j=1}^{p},(\sigma^{2})^{(t+1)}\right)\in
F_{t}\left(\lambda^{(t)},\\{a_{i}^{(t)}\\}_{i=1}^{n},\\{b_{j}^{(t)}\\}_{j=1}^{p},(\sigma^{2})^{(t)}\right).$

Before proceeding further, we identify the set of desired points $P$ as the
singleton set $\\{\theta^{\ast}\\}$. Note that this is essentially a one-step
process. The conditions (C1)-(C3) of Theorem 10.3.3 of Argyros (2008) follows

  1. 1.

Condition (C1) is easily satisfied by the definition of minima, as the
iteration process applied on $\boldsymbol{\theta}^{\ast}$ would retain
$\boldsymbol{\theta}^{\ast}$ as it is the minimum of $\mathcal{H}_{n,p}$
restricting our attention to $\boldsymbol{\Theta}$.

  2. 2.

Condition (C3) is implied by the existence of the compact set
$\boldsymbol{C}$.

  3. 3.

To show that condition (C2) follows, we choose the Lyapunov function as
$\mathcal{H}_{n,p}$ itself. Note that, $\mathcal{H}_{n,p}$ is necessarily
continuous in its arguments. Due to the specific iteration rules in
(2.5)-(2.7), all estimates are obtained in alternating fashion by minimizing
the parameters subject to the most recent value of the other parameters. Thus
it follows that

$\displaystyle\mathcal{H}_{n,p}(\lambda^{(t)},\\{a_{i}^{(t)}\\}_{i=1}^{n},\\{b_{j}^{(t)}\\}_{j=1}^{p},(\sigma^{2})^{(t)})$
$\displaystyle\geq\mathcal{H}_{n,p}(\lambda^{(t+1/2)},\\{a_{i}^{(t+1)}\\}_{i=1}^{n},\\{b_{j}^{(t)}\\}_{j=1}^{p},(\sigma^{2})^{(t)})$
$\displaystyle\geq\mathcal{H}_{n,p}(\lambda^{(t+1)},\\{a_{i}^{(t+1)}\\}_{i=1}^{n},\\{b_{j}^{(t+1)}\\}_{j=1}^{p},(\sigma^{2})^{(t)})$
$\displaystyle\geq\mathcal{H}_{n,p}(\lambda^{(t+1)},\\{a_{i}^{(t+1)}\\}_{i=1}^{n},\\{b_{j}^{(t+1)}\\}_{j=1}^{p},(\sigma^{2})^{(t+1)}).$

Therefore, a full step of iteration monotonically decreases (may not be
strictly) the function $\mathcal{H}_{n,p}$. This shows that the chosen
Lyapunov function is closed.

Now an application of Theorem 10.3.3 of Argyros (2008) completes the proof.

###  B.7 Proof of Theorem 3.8

First, we observe that, the stereographic transformation mentioned in the
discussion prior to Theorem 3.8 can be employed and would remain valid because
of Assumptions (A2) and (A3). Now, to prove the consistency, we shall take a
route similar to the one taken by Ghosh and Basu (2013) as in the case of
MDPDE for INH setup. Instead of showing that the rSVDdpd estimator i.e., the
minimizer $\boldsymbol{\theta}^{\ast}_{n,p}$ of $\mathcal{H}_{n,p}$ is
consistent for $\boldsymbol{\theta}^{g}_{n,p}$, we shall show instead that
$\boldsymbol{\eta}^{\ast}_{n,p}$ is consistent for
$\boldsymbol{\eta}^{g}_{n,p}$. By a slight abuse of notation throughout this
proof, we shall use $\mathcal{H}_{n,p}(\boldsymbol{\eta})$ to indicate
$\mathcal{H}_{n,p}(\mathcal{T}(\boldsymbol{\theta}))$. To prove that
$\boldsymbol{\eta}^{\ast}_{n,p}$ is consistent for
$\boldsymbol{\eta}^{g}_{n,p}$, we shall show that for any sufficiently small
$r>0$,
$\mathcal{H}_{n,p}(\boldsymbol{\eta}_{n,p})>\mathcal{H}_{n,p}(\boldsymbol{\eta}^{g}_{n,p})$
for sufficiently large $n$ and $p$, for any $\boldsymbol{\eta}_{n,p}$ with
$\|\boldsymbol{\eta}_{n,p}-\boldsymbol{\eta}^{g}_{n,p}\|_{2}=r$. This means
that the value of $\mathcal{H}_{n,p}$ at the surface of the ball of radius $r$
centered at $\boldsymbol{\eta}^{g}_{n,p}$ would be higher than its value at
$\boldsymbol{\eta}^{g}_{n,p}$, and hence by the smoothness of
$\mathcal{H}_{n,p}$, it is ensured that there will be a local minima strictly
inside that ball. Proceeding as in Ghosh and Basu (2013), we start with the
Taylor series expansion of $\mathcal{H}_{n,p}(\boldsymbol{\eta}_{n,p})$ about
$\boldsymbol{\eta}^{g}_{n,p}$, for any fixed $n$ and $p$. For notational
convenience, we suppress the subscripts from $\boldsymbol{\eta}$ and
$\boldsymbol{\eta}^{g}$ which should be obvious from the context.

$\displaystyle\mathcal{H}_{n,p}(\boldsymbol{\eta})-\mathcal{H}_{n,p}(\boldsymbol{\eta}^{g})$
$\displaystyle=\left.\frac{\partial\mathcal{H}_{n,p}}{\partial\lambda}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}(\lambda-\lambda^{g})+\sum_{i=1}^{(n-1)}\left.\frac{\partial\mathcal{H}_{n,p}}{\partial\alpha_{i}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}(\alpha_{i}-\alpha_{i}^{g})$
$\displaystyle\quad+\sum_{j=1}^{(p-1)}\left.\frac{\partial\mathcal{H}_{n,p}}{\partial\beta_{j}}\right|_{\boldsymbol{\eta}_{n,p}=\boldsymbol{\eta}^{g}_{n,p}}(\beta_{j}-\beta_{j}^{g})+\left.\frac{\partial\mathcal{H}_{n,p}}{\partial\sigma^{2}}\right|_{\boldsymbol{\eta}_{n,p}=\boldsymbol{\eta}^{g}_{n,p}}(\sigma^{2}-(\sigma^{g})^{2})$
$\displaystyle\quad+\frac{1}{2}\sum_{k_{1},k_{2}}\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\eta_{k_{1}}\partial\eta_{k_{2}}}\right|_{\boldsymbol{\eta}_{n,p}=\boldsymbol{\eta}^{g}_{n,p}}(\eta_{k_{1}}-\eta_{k_{1}}^{g})(\eta_{k_{2}}-\eta_{k_{2}}^{g})$
$\displaystyle\quad+\frac{1}{3!}\sum_{k_{1},k_{2},k_{3}}\left.\frac{\partial^{3}\mathcal{H}_{n,p}}{\partial\eta_{k_{1}}\partial\eta_{k_{2}}\partial\eta_{k_{3}}}\right|_{\boldsymbol{\eta}_{n,p}=\boldsymbol{\eta}^{g}_{n,p}}(\eta_{k_{1}}-\eta_{k_{1}}^{g})(\eta_{k_{2}}-\eta_{k_{2}}^{g})(\eta_{k_{3}}-\eta_{k_{3}}^{g})$
$\displaystyle=S_{1,1}+S_{1,2}+S_{1,3}+S_{1,4}+\frac{1}{2}S_{2}+\frac{1}{6}S_{3},$

where the quantities $S_{1,1},S_{1,2},S_{1,3},S_{1,4},S_{2}$ and $S_{3}$
respectively denote the summands they are replacing. Here, $\eta_{k}$ denotes
the $k$-th coordinate of the vector $\boldsymbol{\eta}_{n,p}$. Also
$\alpha_{i}$’s and $\beta_{j}$’s are the natural parametric representation of
the elements of left ($a_{i,n}$) and right singular vectors ($b_{j,p}$)
respectively, where the dimension subscripts ($n$ and $p$) have been
suppressed for notational convenience as indicated before.

Clearly, the smoothness of $\mathcal{H}_{n,p}$ along with Assumption (A1) on
the normality of the errors, indicates that $\int\mathcal{H}_{n,p}g_{ij}(x)dx$
can be differentiated thrice with respect to $\boldsymbol{\eta}_{n,p}$, and
the derivative can be taken under the integral sign. Thus, it follows that

$\mathbb{E}\left[\left.\frac{\partial\mathcal{H}_{n,p}}{\partial\eta_{k}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right]=\left.\frac{\partial\mathbb{E}\mathcal{H}_{n,p}}{\partial\eta_{k}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}=0,\
k=1,2,\dots(n+p),$ (B.5)

since the population version of the objective function
$\mathbb{E}\mathcal{H}_{n,p}$ is minimized at the true parameter
$\boldsymbol{\eta}^{g}$. Thus, by a generalized version of Khinchin’s Weak Law
of Large numbers, it follows that as $n$ and $p$ both increase to infinity,
each of the first order partial derivatives goes in probability to $0$.
However, the problem arises as there are potentially infinitely many terms (as
the parameter space increases in dimension). This jeopardizes any approach to
naturally extending the proof of Theorem 3.1 of Ghosh and Basu (2013).

Let us consider each of the sums $S_{1,1},S_{1,2},S_{1,3}$ and $S_{1,4}$
pertaining to the first order derivative separately. Since
$\left.\frac{\partial\mathcal{H}_{n,p}}{\partial\lambda}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}$
and
$\left.\frac{\partial\mathcal{H}_{n,p}}{\partial\sigma^{2}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}$
both converges in probability to $0$, hence for sufficiently large $n$ and
$p$, we have $|S_{1,1}|<r^{3}$ and $|S_{1,4}|<r^{3}$ with probability tending
to $1$. Now, to deal with an increasing number of summands in $S_{1,2}$ or
$S_{1,3}$ we apply Chebyshev’s inequality to bound these by bounding its
expectation and variance separately. By chain rule of differentiation,

$s_{n,p}=\sum_{i=1}^{(n-1)}\left.\frac{\partial\mathcal{H}_{n,p}}{\partial\alpha_{i}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}=\sum_{k=1}^{n}\left.\frac{\partial\mathcal{H}_{n,p}}{\partial
a_{k}}\right|_{\boldsymbol{\theta}=\boldsymbol{\theta}^{g}}\left(\sum_{i=1}^{(n-1)}\left.\frac{\partial
a_{k}}{\partial\alpha_{i}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right).$
(B.6)

As in the case of (B.5), one can verify that
$\mathbb{E}\left[\left.\frac{\partial\mathcal{H}_{n,p}}{\partial
a_{k}}\right|_{\boldsymbol{\theta}=\boldsymbol{\theta}^{g}}\right]=0$ for all
$k=1,2,\dots n$, and, therefore, (B.6) implies that $\mathbb{E}(s_{n,p})=0$.

Turning to its variance, since the part of $\mathcal{H}_{n,p}$ dependent on
$a_{k}$ would consist of only the $k$-th row of the data matrix
$\boldsymbol{X}$, which are assumed to be independently distributed in the
current setup, it follows that

$\mathrm{Cov}\left[\left.\frac{\partial\mathcal{H}_{n,p}}{\partial
a_{k}}\right|_{\boldsymbol{\theta}=\boldsymbol{\theta}^{g}},\left.\frac{\partial\mathcal{H}_{n,p}}{\partial
a_{l}}\right|_{\boldsymbol{\theta}=\boldsymbol{\theta}^{g}}\right]=0,\ 1\leq
k,l\leq n,\text{ with }k\neq l.$

Therefore, it follows that

$\displaystyle\mathrm{Var}(s_{n,p})$
$\displaystyle=\sum_{k=1}^{n}\left(\sum_{i=1}^{(n-1)}\left.\frac{\partial
a_{k}}{\partial\alpha_{i}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right)^{2}\mathrm{Var}\left(\left.\frac{\partial\mathcal{H}_{n,p}}{\partial
a_{k}}\right|_{\boldsymbol{\theta}=\boldsymbol{\theta}^{g}}\right)$
$\displaystyle=\sum_{k=1}^{n}\left(\sum_{i=1}^{(n-1)}\left.\frac{\partial
a_{k}}{\partial\alpha_{i}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right)^{2}\frac{1}{n^{2}p^{2}}(2\pi)^{-\alpha}\left(\sigma^{g}\right)^{-(2\alpha+2)}(1+2\alpha)^{-3/2}\sum_{j=1}^{p}(\lambda^{g})^{2}(b_{j}^{g})^{2}$
$\displaystyle=\sum_{k=1}^{n}\left(\sum_{i=1}^{(n-1)}\left.\frac{\partial
a_{k}}{\partial\alpha_{i}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right)^{2}\frac{1}{n^{2}p^{2}}(2\pi)^{-\alpha}\left(\sigma^{g}\right)^{-(2\alpha+2)}(1+2\alpha)^{-3/2}(\lambda^{g})^{2}.$

From Cauchy-Schwartz inequality, it follows that
$\sum_{i=1}^{(n-1)}a_{i}^{g}\leq\sqrt{n-1}$ since
$\sum_{i=1}^{n}(a_{i}^{g})^{2}=1$, and also $|a_{i}^{g}|\leq 1$. This yields

$\sum_{i=1}^{(n-1)}\left.\frac{\partial
a_{k}}{\partial\alpha_{i}}\right|_{\boldsymbol{\eta}}=\begin{cases}(1-a_{n})\sum_{i=1}^{(n-1)}a_{i}&\text{
if }k=n\\\ (1-a_{n})-a_{k}\sum_{i=1}^{(n-1)}a_{i}&\text{ if }k\neq
n\end{cases}$

both of which are less than of equal to $2\sqrt{n-1}$, for any $n\geq 5$.
Therefore, for sufficiently large choice of $n$ and $p$,

$\mathrm{Var}(s_{n,p})\leq
4(2\pi)^{-\alpha}(1+2\alpha)^{-3/2}(\lambda^{g})^{2}\frac{(\sigma^{g})^{-(2\alpha+2)}}{np^{2}}=\Theta\left((np)^{(2\alpha+2)/(\alpha+2)}/np^{2}\right).$

since, $(\sigma^{g})=\Theta\left((np)^{-1/(\alpha+2)}\right)$. Now as
$0\leq\alpha\leq 1<2$ implies $\frac{2(2\alpha+2)}{\alpha+2}<3$ and as
$\lim_{\begin{subarray}{c}n\rightarrow\infty\\\
p\rightarrow\infty\end{subarray}}\frac{n}{p}=c<\infty$, it follows that
$\mathrm{Var}(s_{n,p})\rightarrow 0$ as $n,p\rightarrow\infty$. Thus, the
quantity

$\left|\sum_{i=1}^{(n-1)}\left[\left.\frac{\partial\mathcal{H}_{n,p}}{\partial\alpha_{i}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right]\right|\rightarrow
0,$

with probability tending to one. Therefore, with
$|\alpha_{i}-\alpha_{i}^{g}|<r$, we have $|S_{1,2}|<r^{3}$, for sufficiently
large $n$ and $p$, with probability tending to $1$. Reversing the role of $n$
and $p$, and considering
$\sum_{j=1}^{(p-1)}\left[\left.\frac{\partial\mathcal{H}_{n,p}}{\partial\beta_{j}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right]$
instead, one can show that $|S_{1,3}|<r^{3}$ for sufficiently large $n$ and
$p$ with probability tending to $1$. Thus, combining everything we obtain
$|S_{1}|\leq|S_{1,1}|+|S_{1,2}|+|S_{1,3}|+|S_{1,4}|)<4r^{3}$ for sufficiently
large $n$ and $p$, with probability tending to $1$.

Now, turning our attention to the term $S_{2}$, we start by noting that

$\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\lambda^{2}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta^{g}}}\right]=(2\pi)^{-\alpha/2}(\sigma^{g})^{-(\alpha+2)}(1+\alpha)^{-1/2}\frac{\sum_{i}\sum_{j}(a_{i}^{g})^{2}(b_{j}^{g})^{2}}{np}=\Theta\left(\frac{(\sigma^{g})^{-(\alpha+2)}}{np}\right).$

Because $(\sigma^{g})=\Theta\left((np)^{-\frac{1}{(\alpha+2)}}\right)$ by
condition (A4), it then follows that asymptotically
$\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\lambda^{2}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta^{g}}}\right]$
is bounded away from zero, and, in particular,
$\inf_{n,p}\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\lambda^{2}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta^{g}}}\right]>0$.
The expectated value of the mixed order derivatives, on the other hand, is
equal to $0$ as follows

$\displaystyle\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\lambda\partial\alpha_{i}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta^{g}}}\right]$
$\displaystyle=\sum_{k=1}^{n}\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\lambda\partial
a_{k}}\right|_{\boldsymbol{\theta}=\boldsymbol{\theta^{g}}}\right]\left.\frac{\partial
a_{k}}{\partial\alpha_{i}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta^{g}}}$
$\displaystyle=\sum_{k=1}^{n}(2\pi)^{-\alpha/2}(\sigma^{g})^{-(\alpha+2)}(1+\alpha)^{-1/2}\frac{\sum_{j}\lambda^{g}a_{k}^{g}b_{j}^{g}}{np}\left.\frac{\partial
a_{k}}{\partial\alpha_{i}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta^{g}}}$
$\displaystyle=\frac{\lambda^{g}(2\pi)^{-\alpha/2}(\sigma^{g})^{-(\alpha+2)}(1+\alpha)^{-1/2}}{np}\sum_{k=1}^{n}a_{k}^{g}\left.\frac{\partial
a_{k}}{\partial\alpha_{i}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta^{g}}}=0,$

since $\sum_{i=1}^{n}a_{k}^{g}\left.\frac{\partial
a_{k}^{g}}{\partial\alpha_{i}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}=\frac{\partial}{\partial\alpha_{i}}\left(\sum_{k=1}^{n}a_{k}^{2}\right)=0$.
Similarly, we also have
$\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\lambda\partial\beta_{j}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta^{g}}}\right]=0$.
A chain rule of differentiation can be used to obtain the second order
derivatives of $\mathcal{H}_{n,p}$ with respect to $\alpha_{i}$’s as

$\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\alpha_{i}\partial\alpha_{j}}=\sum_{k=1}^{n}\frac{\partial\mathcal{H}_{n,p}}{\partial
a_{k}}\frac{\partial^{2}a_{k}}{\partial\alpha_{i}\partial\alpha_{j}}+\sum_{k=1}^{n}\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial
a_{k}\partial a_{l}}\frac{\partial a_{k}}{\partial\alpha_{i}}\frac{\partial
a_{k}}{\partial\alpha_{j}},$

as the first order derivative of $\mathcal{H}_{n,p}$ with respect to $a_{k}$
would only yield nonzero quantity corresponding to $k$-th row, and for $k\neq
l$, $\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial a_{k}\partial a_{l}}=0$. A
simple calculation by taking expectation on both sides reveals that

$\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\alpha_{i}\partial\alpha_{j}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right]=\begin{cases}\frac{(\lambda^{g})^{2}(2\pi)^{-\alpha/2}(\sigma^{g})^{-(\alpha+2)}(1+\alpha)^{-1/2}}{np}(1-a_{n}^{g})^{2}&\text{
if }i=j\\\ 0&\text{ if }i\neq j\end{cases}.$

Since $a_{n}^{g}\neq\pm 1$ by condition (B2), and
$\sigma^{g}=\Theta((np)^{-1/(\alpha+2)})$ by condition (A4), it again follows
that
$\inf_{n,p}\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\alpha_{i}^{2}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right]>0$.
A very similar conclusion can be said about the quantities
$\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\beta_{i}\partial\beta_{j}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right]$
as well. As shown in Ghosh and Basu (2013), for normal linear regression setup
the scale and the location estimator become asymptotically uncorrelated.
Similarly, here also we have

$\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\alpha_{i}\partial\sigma^{2}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right]=\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\beta_{j}\partial\sigma^{2}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right]=\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\lambda\partial\sigma^{2}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right]=0,$

and

$\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial(\sigma^{2})^{2}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right]=(2\pi)^{-\alpha/2}(\sigma^{g})^{-(\alpha+4)}\frac{(2+\alpha^{2})}{4(1+\alpha)^{3/2}}.$

Again Assumption (A4) implies that
$\inf_{n,p}\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial(\sigma^{2})^{2}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right]>0$.
Finally, we consider the $(n+p)\times(n+p)$ matrix $\boldsymbol{\Psi}_{n,p}$
whose $(k_{1},k_{2})$-th element is given by
$\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\eta_{k_{1}}\eta_{k_{2}}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right]$.
Then, the previous calculations imply that $\boldsymbol{\Psi}_{n,p}$ is of
full rank for any $n$ and $p$, and also
$\inf_{n,p}\text{tr}\left(\boldsymbol{\Psi}_{n,p}\right)>0$, where
$\text{tr}(\boldsymbol{A})$ denotes the trace of the matrix $\boldsymbol{A}$.
Also, due to Assumption (B1), the integrated likelihood can be differentiated
and the derivative can be obtained under the integral sign. Therefore,
$\mathbb{E}\left[\left.\frac{\partial^{2}\mathcal{H}_{n,p}}{\partial\eta_{k_{1}}\partial\eta_{k_{2}}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right]=\left.\frac{\partial\mathbb{E}(\mathcal{H}_{n,p})}{\partial\eta_{k_{1}}\partial\eta_{k_{2}}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}$,
and hence, $\boldsymbol{\Psi}_{n,p}$ is the Fisher’s information matrix at the
true parameter $\boldsymbol{\eta}^{g}$. Therefore, the matrix
$\boldsymbol{\Psi}_{n,p}$ is positive semi-definite for each $n$ and $p$.
Together with the result on full rank of $\boldsymbol{\Psi}_{n,p}$ and
$\inf_{n,p}\text{tr}(\boldsymbol{\Psi}_{n,p})>0$, it implies that
$\boldsymbol{\Psi}_{n,p}$ is positive definite for each $n$ and $p$, and
$\inf_{n,p}\text{minimum eigenvalue of
}\boldsymbol{\Psi}_{n,p}=\lambda_{0}>0$.

Now, we decompose $S_{2}$ by considering elements of $\boldsymbol{\Psi}_{n,p}$
as follows

$\sum_{k_{1},k_{2}}(\boldsymbol{\Psi}_{n,p})_{k_{1},k_{2}}(\eta_{k_{1}}-\eta_{k_{1}}^{g})(\eta_{k_{2}}-\eta_{k_{2}}^{g})=\sum_{k_{1},k_{2}}\left[(\boldsymbol{\Psi}_{n,p})_{k_{1},k_{2}}-E\left((\boldsymbol{\Psi}_{n,p})_{k_{1},k_{2}}\right)\right](\eta_{k_{1}}-\eta_{k_{1}}^{g})(\eta_{k_{2}}-\eta_{k_{2}}^{g})\\\
+\sum_{k_{1},k_{2}}\mathbb{E}\left((\boldsymbol{\Psi}_{n,p})_{k_{1},k_{2}}\right)(\eta_{k_{1}}-\eta_{k_{1}}^{g})(\eta_{k_{2}}-\eta_{k_{2}}^{g})$

Here, we can apply an orthogonal transformation on
$(\boldsymbol{\eta}-\boldsymbol{\eta}^{g})$ to express it as a linear
combination of the eigenvectors of $\Psi_{n,p}$, so that the second term can
be made greater than or equal to $\lambda_{0}r^{2}$. Also, it is evident that
the first summation has expected value equal to $0$. By a similar routine
calculation, one can show that the variance of the first term goes to $0$ as
well. Therefore, for sufficiently large $n$ and $p$, with probability tending
to $1$, $S_{2}>(-r^{3}+\lambda_{0}r^{2})$.

Finally, turning to $S_{3}$, we note that the expectated values of the third
order derivatives are bounded as shown below.

$\displaystyle\left|\mathbb{E}\left[\left.\frac{\partial^{3}\mathcal{H}_{n,p}}{\partial\lambda^{3}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}\right]\right|$
$\displaystyle=M_{1}\times\left|\frac{(\sigma^{g})^{-(\alpha+2)}}{np}\sum_{i,j}(a_{i}^{g}b_{j}^{g})^{3}\right|\leq
M_{1}\times\frac{(\sigma^{g})^{-(\alpha+2)}}{np},$ (B.7)
$\displaystyle\left|\mathbb{E}\left[\left.\frac{\partial^{3}\mathcal{H}_{n,p}}{\partial
a_{i}^{3}}\right|_{\boldsymbol{\theta}=\boldsymbol{\theta}^{g}}\right]\right|$
$\displaystyle=M_{2}\times\left|\frac{(\sigma^{g})^{-(\alpha+2)}}{np}\sum_{j}(\lambda^{g}b_{j}^{g})^{3}\right|\leq
M_{2}\times(\lambda^{g})^{3}\frac{(\sigma^{g})^{-(\alpha+2)}}{np},$ (B.8)
$\displaystyle\left|\mathbb{E}\left[\left.\frac{\partial^{3}\mathcal{H}_{n,p}}{\partial
b_{j}^{3}}\right|_{\boldsymbol{\theta}=\boldsymbol{\theta}^{g}}\right]\right|$
$\displaystyle=M_{3}\times\left|\frac{(\sigma^{g})^{-(\alpha+2)}}{np}\sum_{i}(\lambda^{g}a_{i}^{g})^{3}\right|\leq
M_{3}\times(\lambda^{g})^{3}\frac{(\sigma^{g})^{-(\alpha+2)}}{np},$ (B.9)

where $M_{1},M_{2},M_{3}$ are some finite constant. The inequalities follow
from Cauchy-Schwartz inequality and the observation that
$\sum_{i}(a_{i}^{g})^{2}=\sum_{j}(b_{j}^{g})^{2}=1$. Because of Assumption
(A4), the quantity $\frac{(\sigma^{g})^{-(\alpha+2)}}{np}=\Theta(1)$. Since
these bounds now do not depend on the particular choice of $i$ (or $j$), one
can show that these expectated values of third order derivatives are uniformly
bounded for any $n$ and $p$. This would imply that the third order derivatives
$\left.\frac{\partial\mathcal{H}_{n,p}}{\partial\eta_{k_{1}}\partial\eta_{k_{2}}\partial\eta_{k_{3}}}\right|_{\boldsymbol{\eta}=\boldsymbol{\eta}^{g}}$
as random variables remain uniformly bounded for
$\eta_{k_{1}},\eta_{k_{2}},\eta_{k_{3}}\in\\{\lambda,\alpha_{1},\dots\alpha_{(n-1)},\beta_{1},\dots\beta_{(p-1)}\\}$.
Let us denote this bound by $\widetilde{M}_{1}$. However, the expectation of
the third order derivative of $\mathcal{H}_{n,p}$ with respect to $\sigma^{2}$
is not finite in view of Assumption (A4) , since it is
$\Theta\left((\sigma^{g})^{-(\alpha+4)}/np\right)$. To circumvent this
problem, we use Assumption (A5) that the random entries $X_{ij}$ of the data
matrix are uniformly bounded for any sufficiently large $n$ and $p$ with
probability tending to $1$. Therefore,
$\frac{\partial^{3}\mathcal{H}_{n,p}}{\partial(\sigma^{2})^{3}}$, being a
continuous function of the entries $X_{ij}$, also remains uniformly bounded
for any sufficiently large $n$ and $p$ with probability tending to $1$. We
denote this bound by $\widetilde{M}_{2}$.

Combining these two bounds and taking $M=\max\\{M_{1},M_{2}\\}$, we obtain
that $S_{3}\leq
M\sum_{k_{1},k_{2},k_{3}}(\eta_{k_{1}}-\eta_{k_{1}}^{g})(\eta_{k_{2}}-\eta_{k_{2}}^{g})(\eta_{k_{3}}-\eta_{k_{3}}^{g})$
for sufficiently large $n$ and $p$, with probability tending to $1$. Since it
is well known that in finite-dimensional vector space, the usual $L_{p}$ norms
are equivalent, therefore,
$\sum_{k_{1},k_{2},k_{3}}(\eta_{k_{1}}-\eta_{k_{1}}^{g})(\eta_{k_{2}}-\eta_{k_{2}}^{g})(\eta_{k_{3}}-\eta_{k_{3}}^{g})\leq
Cr^{3}$ for some finite constant $C$. Therefore, for sufficiently large $n$
and $p$, $S_{3}<MCr^{3}$, with probability tending to $1$. Therefore, for
sufficiently large $n$ and $p$

$\mathcal{H}_{n,p}(\boldsymbol{\eta})-\mathcal{H}_{n,p}(\boldsymbol{\eta}^{g})>(-4r^{3}+\lambda_{0}r^{2}-r^{3}-MCr^{3})=\lambda_{0}r^{2}-(MC+5)r^{3}$

with probability tending to $1$. Now, a choice of $r<\lambda_{0}/(MC+5)$ would
ensure that
$\mathcal{H}_{n,p}(\boldsymbol{\eta})>\mathcal{H}_{n,p}(\boldsymbol{\eta}^{g})$
for any $\boldsymbol{\eta}$ satisfying
$\|\boldsymbol{\eta}-\boldsymbol{\eta}^{g}\|_{2}=r$ (where $\|\cdot\|_{2}$
denotes the Euclidean $L_{2}$ norm). This is exactly what we intended to show
at the beginning.

Finally, following the discussion at the beginning of the proof, we obtain
that $\boldsymbol{\eta}^{\ast}$ is consistent for $\boldsymbol{\eta}^{g}$, and
by continuous mapping theorem, we also have $\boldsymbol{\theta}^{\ast}$
consistent for $\boldsymbol{\theta}^{g}$. This completes the proof.

## References

  * Golub and Reinsch [1970] G. H. Golub and C. Reinsch.  Singular value decomposition and least squares solutions.  _Numerische Mathematik_ , 14(5):403–420, Apr 1970\.  ISSN 0945-3245.  doi: 10.1007/BF02163027.  URL https://doi.org/10.1007/BF02163027.
  * Greenacre [2017] Michael Greenacre.  _Correspondence analysis in practice_.  CRC press, 2017.
  * Hofmann [1999] Thomas Hofmann.  Probabilistic latent semantic indexing.  In _Proceedings of the 22nd Annual International ACM SIGIR Conference on Research and Development in Information Retrieval_ , SIGIR ’99, page 50–57, New York, NY, USA, 1999. Association for Computing Machinery.  ISBN 1581130961.  doi: 10.1145/312624.312649.  URL https://doi.org/10.1145/312624.312649.
  * Anandarajan et al. [2019] Murugan Anandarajan, Chelsey Hill, and Thomas Nolan.  _Semantic Space Representation and Latent Semantic Analysis_ , pages 77–91.  Springer International Publishing, Cham, 2019.  ISBN 978-3-319-95663-3.  doi: 10.1007/978-3-319-95663-3\\_6.  URL https://doi.org/10.1007/978-3-319-95663-3_6.
  * Drineas et al. [2004] Petros Drineas, Alan Frieze, Ravi Kannan, Santosh Vempala, and V Vinay.  Clustering large graphs via the singular value decomposition.  _Machine learning_ , 56(1-3):9–33, 2004.
  * Cheng et al. [2019] Hongchuan Cheng, Yimin Zhang, Wenjia Lu, and Zhou Yang.  A bearing fault diagnosis method based on vmd-svd and fuzzy clustering.  _International Journal of Pattern Recognition and Artificial Intelligence_ , 33(12):1950018, 2019.  doi: 10.1142/S0218001419500186.  URL https://doi.org/10.1142/S0218001419500186.
  * Dappuri et al. [2020] Bhasker Dappuri, M. Purnachandra Rao, and Madhu Babu Sikha.  Non-blind rgb watermarking approach using svd in translation invariant wavelet space with enhanced grey-wolf optimizer.  _Multimedia Tools Appl._ , 79(41–42):31103–31124, nov 2020.  ISSN 1380-7501.  doi: 10.1007/s11042-020-09433-0.  URL https://doi.org/10.1007/s11042-020-09433-0.
  * Zhao and Jia [2017] Ming Zhao and Xiaodong Jia.  A novel strategy for signal denoising using reweighted svd and its applications to weak fault feature enhancement of rotating machinery.  _Mechanical Systems and Signal Processing_ , 94:129–147, 2017.  ISSN 0888-3270.  doi: https://doi.org/10.1016/j.ymssp.2017.02.036.  URL https://www.sciencedirect.com/science/article/pii/S0888327017301061.
  * Rezaei and Khalili [2019] Azadeh Rezaei and Mehdi Khalili.  A robust blind audio watermarking scheme based on dct-dwt-svd.  In Shahram Montaser Kouhsari, editor, _Fundamental Research in Electrical Engineering_ , pages 101–113, Singapore, 2019. Springer Singapore.  ISBN 978-981-10-8672-4.
  * Grondin and Glass [2019] François Grondin and James Glass.  Svd-phat: A fast sound source localization method.  In _ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 4140–4144, 2019.  doi: 10.1109/ICASSP.2019.8683253.
  * Zhang et al. [2016] Dashan Zhang, Jie Guo, Xiujun Lei, and Chang’an Zhu.  Note: Sound recovery from video using svd-based information extraction.  _Review of Scientific Instruments_ , 87(8):086111, 2016.  doi: 10.1063/1.4961979.  URL https://doi.org/10.1063/1.4961979.
  * Franceschini et al. [2015] Andrea Franceschini, Jianyi Lin, Christian von Mering, and Lars Juhl Jensen.  SVD-phy: improved prediction of protein functional associations through singular value decomposition of phylogenetic profiles.  _Bioinformatics_ , 32(7):1085–1087, 11 2015.  ISSN 1367-4803.  doi: 10.1093/bioinformatics/btv696.  URL https://doi.org/10.1093/bioinformatics/btv696.
  * Bustamam et al. [2018] A. Bustamam, S. Formalidin, and T. Siswantining.  Clustering and analyzing microarray data of lymphoma using singular value decomposition (svd) and hybrid clustering.  _AIP Conference Proceedings_ , 2023(1):020220, 2018.  doi: 10.1063/1.5064217.  URL https://aip.scitation.org/doi/abs/10.1063/1.5064217.
  * Das et al. [2017] Lopamudra Das, J. K. Das, and Sarita Nanda.  Advanced protein coding region prediction applying robust svd algorithm.  In _2017 2nd International Conference on Man and Machine Interfacing (MAMI)_ , pages 1–6, 2017.  doi: 10.1109/MAMI.2017.8307887.
  * Kumar et al. [2011] Nishith Kumar, Mohammed Nasser, and Subaran Chandra Sarker.  A new singular value decomposition based robust graphical clustering technique and its application in climatic data.  _Journal of Geography and Geology_ , 3(1):227, 2011.
  * Hawkins et al. [2001] Douglas M. Hawkins, Li Liu, and Stanley Young.  Robust singular value decomposition.  Technical Report 122, National Institute of Statistical Sciences (NISS), 12 2001.
  * Liu et al. [2003] Li Liu, Douglas M. Hawkins, Sujoy Ghosh, and S. Stanley Young.  Robust singular value decomposition analysis of microarray data.  _Proceedings of the National Academy of Sciences_ , 100(23):13167–13172, 2003.  doi: 10.1073/pnas.1733249100.  URL https://www.pnas.org/doi/abs/10.1073/pnas.1733249100.
  * Ammann [1993] Larry P. Ammann.  Robust singular value decompositions: A new approach to projection pursuit.  _Journal of the American Statistical Association_ , 88(422):505–514, 1993.  doi: 10.1080/01621459.1993.10476301.  URL https://www.tandfonline.com/doi/abs/10.1080/01621459.1993.10476301.
  * Ke and Kanade [2005] Qifa Ke and Takeo Kanade.  Robust l1 norm factorization in the presence of outliers and missing data by alternative convex programming.  In _Proceedings of (CVPR) Computer Vision and Pattern Recognition_ , volume 1, pages 739 – 746, June 2005.  ISBN 0-7695-2372-2.  doi: 10.1109/CVPR.2005.309.
  * Jung [2010] Kang-Mo Jung.  Robust singular value decomposition balsed on weighted least absolute deviation regression.  _Communications for Statistical Applications and Methods_ , 17(6):803–810, 11 2010.
  * Rey [2007] William Rey.  Total singular value decomposition. robust svd, regression and location-scale.  2007\.
  * Markovsky and Van Huffel [2007] Ivan Markovsky and Sabine Van Huffel.  Overview of total least-squares methods.  _Signal Processing_ , 87(10):2283–2302, 2007\.  ISSN 0165-1684.  doi: https://doi.org/10.1016/j.sigpro.2007.04.004.  URL https://www.sciencedirect.com/science/article/pii/S0165168407001405.  Special Section: Total Least Squares and Errors-in-Variables Modeling.
  * Zhang et al. [2013] Lingsong Zhang, Haipeng Shen, and Jianhua Z. Huang.  Robust regularized singular value decomposition with application to mortality data.  _The Annals of Applied Statistics_ , 7(3):1540 – 1561, 2013.  doi: 10.1214/13-AOAS649.  URL https://doi.org/10.1214/13-AOAS649.
  * Wang [2017] Deshen Wang.  Adjustable robust singular value decomposition: Design, analysis and application to finance.  _Data_ , 2(3), 2017.  ISSN 2306-5729.  doi: 10.3390/data2030029.  URL https://www.mdpi.com/2306-5729/2/3/29.
  * Gabriel and Zamir [1979] K. Ruben Gabriel and S. Zamir.  Lower rank approximation of matrices by least squares with any choice of weights.  _Technometrics_ , 21(4):489–498, 1979.  ISSN 00401706.  URL http://www.jstor.org/stable/1268288.
  * Roy et al. [2021] Subhrajyoty Roy, Ayanendranath Basu, and Abhik Ghosh.  rsvddpd: A robust scalable video surveillance background modelling algorithm.  _arXiv preprint arXiv:2109.10680_ , 2021.
  * Basu et al. [1998] Ayanendranath Basu, Ian R. Harris, Nils L. Hjort, and M. C. Jones.  Robust and efficient estimation by minimising a density power divergence.  _Biometrika_ , 85(3):549–559, 1998.  ISSN 00063444.  URL http://www.jstor.org/stable/2337385.
  * Ghosh and Basu [2013] Abhik Ghosh and Ayanendranath Basu.  Robust estimation for independent non-homogeneous observations using density power divergence with applications to linear regression.  _Electronic Journal of Statistics_ , 7(none):2420 – 2456, 2013.  doi: 10.1214/13-EJS847.  URL https://doi.org/10.1214/13-EJS847.
  * Cichocki et al. [2011] Andrzej Cichocki, Sergio Cruces, and Shun-ichi Amari.  Generalized alpha-beta divergences and their application to robust nonnegative matrix factorization.  _Entropy_ , 13(1):134–170, 2011.  ISSN 1099-4300.  doi: 10.3390/e13010134.  URL https://www.mdpi.com/1099-4300/13/1/134.
  * Giraud et al. [2005] L. Giraud, J. Langou, and M. Rozloznik.  The loss of orthogonality in the gram-schmidt orthogonalization process.  _Computers & Mathematics with Applications_, 50(7):1069–1075, 2005.  ISSN 0898-1221.  doi: https://doi.org/10.1016/j.camwa.2005.08.009.  URL https://www.sciencedirect.com/science/article/pii/S0898122105003366.  Numerical Methods and Computational Mechanics.
  * Mehta [2004] Madan Lal Mehta.  _Random matrices_.  Elsevier, 2004.
  * Tracy and Widom [1993] Craig A. Tracy and Harold Widom.  Introduction to random matrices.  In G. F. Helminck, editor, _Geometric and Quantum Aspects of Integrable Systems_ , pages 103–130, Berlin, Heidelberg, 1993. Springer Berlin Heidelberg.  ISBN 978-3-540-48090-7.
  * Huber [1973] Peter J. Huber.  Robust Regression: Asymptotics, Conjectures and Monte Carlo.  _The Annals of Statistics_ , 1(5):799 – 821, 1973\.  doi: 10.1214/aos/1176342503.  URL https://doi.org/10.1214/aos/1176342503.
  * Portnoy [1984] Stephen Portnoy.  Asymptotic Behavior of $M$-Estimators of $p$ Regression Parameters when $p^{2}/n$ is Large. I. Consistency.  _The Annals of Statistics_ , 12(4):1298 – 1309, 1984.  doi: 10.1214/aos/1176346793.  URL https://doi.org/10.1214/aos/1176346793.
  * He and Shao [2000] Xuming He and Qi-Man Shao.  On parameters of increasing dimensions.  _Journal of Multivariate Analysis_ , 73(1):120–135, 2000.  ISSN 0047-259X.  doi: https://doi.org/10.1006/jmva.1999.1873.  URL https://www.sciencedirect.com/science/article/pii/S0047259X99918730.
  * Stacklies et al. [2007] Wolfram Stacklies, Henning Redestig, Matthias Scholz, Dirk Walther, and Joachim Selbig.  pcamethods – a bioconductor package providing pca methods for incomplete data.  _Bioinformatics_ , 23:1164–1167, 2007.
  * Zhang and Pan [2013] Lingsong Zhang and Chao Pan.  _RobRSVD: Robust Regularized Singular Value Decomposition_ , 2013\.  URL https://CRAN.R-project.org/package=RobRSVD.  R package version 1.0.
  * Rudelson and Vershynin [2010] Mark Rudelson and Roman Vershynin.  Non-asymptotic theory of random matrices: extreme singular values.  In _Proceedings of the International Congress of Mathematicians 2010 (ICM 2010) (In 4 Volumes) Vol. I: Plenary Lectures and Ceremonies Vols. II–IV: Invited Lectures_ , pages 1576–1602. World Scientific, 2010.
  * Alman and Williams [2021] Josh Alman and Virginia Vassilevska Williams.  _A Refined Laser Method and Faster Matrix Multiplication_ , pages 522–539.  2021\.  doi: 10.1137/1.9781611976465.32.  URL https://epubs.siam.org/doi/abs/10.1137/1.9781611976465.32.
  * Argyros [2008] I.K. Argyros.  _Convergence and Applications of Newton-type Iterations_.  Springer New York, 2008.  ISBN 9780387727431.  URL https://books.google.co.in/books?id=j5_vsZGAWLcC.