# Nonlinear Inductive Matrix Completion based on One-layer Neural Networks

Kai Zhong
<EMAIL_ADDRESS>
UT-Austin Supported in part by NSF grants CCF-1320746, IIS-1546452 and
CCF-1564000. Zhao Song
<EMAIL_ADDRESS>
Harvard University & UT-Austin Work done while hosted by Jelani Nelson.
Prateek Jain
<EMAIL_ADDRESS>
Microsoft Research, India  Inderjit S. Dhillon
<EMAIL_ADDRESS>
UT-Austin & Amazon/A9 Supported in part by NSF grants CCF-1320746, IIS-1546452
and CCF-1564000.

The goal of a recommendation system is to predict the interest of a user in a
given item by exploiting the existing set of ratings as well as certain
user/item features. A standard approach to modeling this problem is Inductive
Matrix Completion where the predicted rating is modeled as an inner product of
the user and the item features projected onto a latent space. In order to
learn the parameters effectively from a small number of observed ratings, the
latent space is constrained to be low-dimensional which implies that the
parameter matrix is constrained to be low-rank. However, such bilinear
modeling of the ratings can be limiting in practice and non-linear prediction
functions can lead to significant improvements. A natural approach to
introducing non-linearity in the prediction function is to apply a non-linear
activation function on top of the projected user/item features. Imposition of
non-linearities further complicates an already challenging problem that has
two sources of non-convexity: a) low-rank structure of the parameter matrix,
and b) non-linear activation function. We show that one can still solve the
non-linear Inductive Matrix Completion problem using gradient descent type
methods as long as the solution is initialized well. That is, close to the
optima, the optimization function is strongly convex and hence admits standard
optimization techniques, at least for certain activation functions, such as
Sigmoid and tanh. We also highlight the importance of the activation function
and show how ReLU can behave significantly differently than say a sigmoid
function. Finally, we apply our proposed technique to recommendation systems
and semi-supervised clustering, and show that our method can lead to much
better performance than standard linear Inductive Matrix Completion methods.

##  1 Introduction

Matrix Completion (MC) or Collaborative filtering [CR09, GUH16] is by now a
standard technique to model recommendation systems problems where a few user-
item ratings are available and the goal is to predict ratings for any user-
item pair. However, standard collaborative filtering suffers from two
drawbacks: 1) Cold-start problem: MC can’t give prediction for new users or
items, 2) Missing side-information: MC cannot leverage side-information that
is typically present in recommendation systems such as features for
users/items. Consequently, several methods [ABEV06, Ren10, XJZ13, JD13] have
been proposed to leverage the side information together with the ratings.
Inductive matrix completion (IMC) [ABEV06, JD13] is one of the most popular
methods in this class.

IMC models the ratings as the inner product between certain linear mapping of
the user/items’ features, i.e., $A(x,y)=\langle U^{\top}x,V^{\top}y\rangle$,
where $A(x,y)$ is the predicted rating of user $x$ for item $y$,
$x\in\mathbb{R}^{d_{1}},y\in\mathbb{R}^{d_{2}}$ are the feature vectors.
Parameters $U\in\mathbb{R}^{d_{1}\times k},V\in\mathbb{R}^{d_{2}\times k}$
($k\leq d_{1},k\leq d_{2}$) can typically be learned using a small number of
observed ratings [JD13].

However, the bilinear structure of IMC is fairly simplistic and limiting in
practice and might lead to fairly poor accuracy on real-world recommendation
problems. For example, consider the Youtube recommendation system [CAS16] that
requires predictions over videos. Naturally, a linear function over the pixels
of videos will lead to fairly inaccurate predictions and hence one needs to
model the videos using non-linear networks. The survey paper by [ZYS17]
presents many more such examples, where we need to design a non-linear ratings
prediction function for the input features, including [LLL+16] for image
recommendation, [WW14] for music recommendation and [ZYL+16] for
recommendation systems with multiple types of inputs.

We can introduce non-linearity in the prediction function using several
standard techniques, however, if our parameterization admits too many free
parameters then learning them might be challenging as the number of available
user-item ratings tend to be fairly small. Instead, we use a simple non-linear
extension of IMC that can control the number of parameters to be estimated.
Note that IMC based prediction function can be viewed as an inner product
between certain latent user-item features where the latent features are a
linear map of the raw user-item features. To introduce non-linearity, we can
use a non-linear mapping of the raw user-item features rather than the linear
mapping used by IMC. This leads to the following general framework that we
call non-linear inductive matrix completion (NIMC),

$A(x,y)=\langle\mathcal{U}(x),\mathcal{V}(y)\rangle,$ (1)

where $x\in\mathcal{X},y\in\mathcal{Y}$ are the feature vectors, $A(x,y)$ is
their rating and
$\mathcal{U}:\mathcal{X}\rightarrow\mathcal{S},\mathcal{V}:\mathcal{Y}\rightarrow\mathcal{S}$
are non-linear mappings from the raw feature space to the latent space.

The above general framework reduces to standard inductive matrix completion
when $\mathcal{U},\mathcal{V}$ are linear mappings and further reduces to
matrix completion when $x_{i},y_{j}$ are unit vectors $e_{i},e_{j}$ for $i$-th
item and $j$-th user respectively. When $[x_{i},e_{i}]$ is used as the feature
vector and $\mathcal{U}$ is restricted to be a two-block (one for $x_{i}$ and
the other for $e_{i}$) diagonal matrix, then the above framework reduces to
the dirtyIMC model [CHD15]. Similarly, $\mathcal{U}/\mathcal{V}$ can also be
neural networks (NNs), such as feedforward NNs [SCH+16, CAS16], convolutional
NNs for images and recurrent NNs for speech/text.

In this paper, we focus on a simple nonlinear activation based mapping for the
user-item features. That is, we set $\mathcal{U}(x)=\phi(U^{*\top}x)$ and
$\mathcal{V}(x)=\phi(V^{*\top}x)$ where $\phi$ is a nonlinear activation
function $\phi$. Note that if $\phi$ is ReLU then the latent space is
guaranteed to be in non-negative orthant which in itself can be a desirable
property for certain recommendation problems.

Note that parameter estimation in both IMC and NIMC models is hard due to non-
convexity of the corresponding optimization problem. However, for "nice" data,
several strong results are known for the linear models, such as [CR09, JNS13,
GJZ17] for MC and [JD13, XJZ13, CHD15] for IMC. However, non-linearity in NIMC
models adds to the complexity of an already challenging problem and has not
been studied extensively, despite its popularity in practice.

In this paper, we study a simple one-layer neural network style NIMC model
mentioned above. In particular, we formulate a squared-loss based optimization
problem for estimating parameters $U^{*}$ and $V^{*}$. We show that under a
realizable model and Gaussian input assumption, the objective function is
locally strongly convex within a "reasonably large" neighborhood of the ground
truth. Moreover, we show that the above strong convexity claim holds even if
the number of observed ratings is nearly-linear in dimension and polynomial in
the conditioning of the weight matrices. In particular, for well-conditioned
matrices, we can recover the underlying parameters using only
$\operatorname{poly}\log(d_{1}+d_{2})$ user-item ratings, which is critical
for practical recommendation systems as they tend to have very few ratings
available per user. Our analysis covers popular activation functions, e.g.,
sigmoid and ReLU, and discuss various subtleties that arise due to the
activation function. Finally we discuss how we can leverage standard tensor
decomposition techniques to initialize our parameters well. We would like to
stress that practitioners typically use random initialization itself, and
hence results studying random initialization for NIMC model would be of
significant interest.

As mentioned above, due to non-linearity of activation function along with
non-convexity of the parameter space, the existing proof techniques do not
apply directly to the problem. Moreover, we have to carefully argue about both
the optimization landscape as well as the sample complexity of the algorithm
which is not carefully studied for neural networks. Our proof establishes some
new techniques that might be of independent interest, e.g., how to handle the
redundancy in the parameters for ReLU activation. To the best of our
knowledge, this is one of the first theoretically rigorous study of neural-
network based recommendation systems and will hopefully be a stepping stone
for similar analysis for "deeper" neural networks based recommendation
systems. We would also like to highlight that our model can be viewed as a
strict generalization of a one-hidden layer neural network, hence our result
represents one of the few rigorous guarantees for models that are more
powerful than one-hidden layer neural networks [LY17, BGMSS18, ZSJ+17].

Finally, we apply our model on synthetic datasets and verify our theoretical
analysis. Further, we compare our NIMC model with standard linear IMC on
several real-world recommendation-type problems, including user-movie rating
prediction, gene-disease association prediction and semi-supervised
clustering. NIMC demonstrates significantly superior performance over IMC.

###  1.1 Related work

Collaborative filtering: Our model is a non-linear version of the standard
inductive matrix completion model [JD13]. Practically, IMC has been applied to
gene-disease prediction [ND14], matrix sensing [ZJD15], multi-label
classification[YJKD14], blog recommender system [SCLD15], link prediction
[CHD15] and semi-supervised clustering [CHD15, SCH+16]. However, IMC restricts
the latent space of users/items to be a linear transformation of the
user/item’s feature space. [SCH+16] extended the model to a three-layer neural
network and showed significantly better empirical performance for multi-
label/multi-class classification problem and semi-supervised problems.

Although standard IMC has linear mappings, it is still a non-convex problem
due to the bilinearity $UV^{\top}$. To deal with this non-convex problem,
[JD13, Har14] provided recovery guarantees using alternating minimization with
sample complexity linear in dimension. [XJZ13] relaxed this problem to a
nuclear-norm problem and also provided recovery guarantees. More general norms
have been studied [RSW16, SWZ17a, SWZ17b, SWZ18], e.g. weighted Frobenius
norm, entry-wise $\ell_{1}$ norm. More recently, [ZDG18] uses gradient-based
non-convex optimization and proves a better sample complexity. [CHD15] studied
dirtyIMC models and showed that the sample complexity can be improved if the
features are informative when compared to matrix completion. Several low-rank
matrix sensing problems [ZJD15, GJZ17] are also closely related to IMC models
where the observations are sampled only from the diagonal elements of the
rating matrix. [Ren10, LY16] introduced and studied an alternate framework for
ratings prediction with side-information but the prediction function is linear
in their case as well.

Neural networks: Nonlinear activation functions play an important role in
neural networks. Recently, several powerful results have been discovered for
learning one-hidden-layer feedforward neural networks [Tia17, ZSJ+17, JSA15,
LY17, BGMSS18, GKKT17], convolutional neural networks [BG17, ZSD17, DLT18a,
DLT+18b, GKM18]. However, our problem is a strict generalization of the one-
hidden layer neural network and is not covered by the above mentioned results.

Notations. For any function $f$, we define $\widetilde{O}(f)$ to be
$f\cdot\log^{O(1)}(f)$. For two functions $f,g$, we use the shorthand
$f\lesssim g$ (resp. $\gtrsim$) to indicate that $f\leq Cg$ (resp. $\geq$) for
an absolute constant $C$. We use $f\eqsim g$ to mean $cf\leq g\leq Cf$ for
constants $c,C$. We use $\operatorname{poly}(f)$ to denote $f^{O(1)}$.

Roadmap. We first present the formal model and the corresponding optimization
problem in Section 2. We then present the local strong convexity and local
linear convergence results in Section 3. Finally, we demonstrate the empirical
superiority of NIMC when compared to linear IMC (Section 4).

##  2 Problem Formulation

Consider a user-item recommender system, where we have $n_{1}$ users with
feature vectors $X:=\\{x_{i}\\}_{i\in[n_{1}]}\subseteq\mathbb{R}^{d_{1}}$,
$n_{2}$ items with feature vectors
$Y:=\\{y_{j}\\}_{j\in[n_{2}]}\subseteq\mathbb{R}^{d_{2}}$ and a collection of
partially-observed user-item ratings, ${\cal
A}_{\text{obs}}=\\{A(x,y)|(x,y)\in\Omega\subseteq X\times Y\\}$. That is
$A(x_{i},y_{j})$ is the rating that user $x_{i}$ gave for item $y_{j}$. For
simplicity, we assume $x_{i}$’s and $y_{j}$’s are sampled i.i.d. from
distribution $\mathcal{X}$ and $\mathcal{Y}$, respectively. Each element of
the index set $\Omega$ is also sampled independently and uniformly with
replacement from $S:=X\times Y$.

In this paper, our goal is to predict the rating for any user-item pair with
feature vectors $x$ and $y$, respectively. We model the user-item ratings as:

$A(x,y)=\phi(U^{*\top}x)^{\top}\phi(V^{*\top}y),$ (2)

where $U^{*}\in\mathbb{R}^{d_{1}\times k}$, $V^{*}\in\mathbb{R}^{d_{2}\times
k}$ and $\phi$ is a non-linear activation function. Under this realizable
model, our goal is to recover $U^{*},V^{*}$ from a collection of observed
entries, $\\{A(x,y)|(x,y)\in\Omega\\}$. Without loss of generality, we set
$d_{1}=d_{2}$. Also we treat $k$ as a constant throughout the paper. Our
analysis requires $U^{*},V^{*}$ to be full column rank, so we require $k\leq
d$. And w.l.o.g., we assume $\sigma_{k}(U^{*})=\sigma_{k}(V^{*})=1$, i.e., the
smallest singular value of both $U^{*}$ and $V^{*}$ is $1$.

Note that this model is similar to one-hidden layer feed-forward network
popular in standard classification/regression tasks. However, as there is an
inner product between the output of two non-linear layers, $\phi(U^{*}x)$ and
$\phi(V^{*}y)$, it cannot be modeled by a single hidden layer neural network
(with same number of nodes). Also, for linear activation function, the problem
reduces to inductive matrix completion [ABEV06, JD13].

Now, to solve for $U^{*}$, $V^{*}$, we optimize a simple squared-loss based
optimization problem, i.e.,

$\min_{U\in\operatorname{{\mathbb{R}}}^{d_{1}\times
k},V\in\operatorname{{\mathbb{R}}}^{d_{2}\times
k}}f_{\Omega}(U,V)=\sum_{(x,y)\in\Omega}(\phi(U^{\top}x)^{\top}\phi(V^{\top}y)-A(x,y))^{2}.$
(3)

Naturally, the above problem is a challenging non-convex optimization problem
that is strictly harder than two non-convex optimization problems which are
challenging in their own right: a) the linear inductive matrix completion
where non-convexity arises due to bilinearity of $U^{\top}V$, and b) the
standard one-hidden layer neural network (NN). In fact, recently a lot of
research has focused on understanding various properties of both the linear
inductive matrix completion problem [GJZ17, JD13] as well as one-hidden layer
NN [GLM18, ZSJ+17].

In this paper, we show that despite the non-convexity of Problem (3), it
behaves as a convex optimization problem close to the optima if the data is
sampled stochastically from a Gaussian distribution. This result combined with
standard tensor decomposition based initialization [ZSJ+17, KCL15, JSA15]
leads to a polynomial time algorithm for solving (3) optimally if the data
satisfies certain sampling assumptions in Theorem 2.1. Moreover, we also
discuss the effect of various activation functions, especially the difference
between a sigmoid activation function vs RELU activation (see Theorem 3.2 and
Theorem 3.4).

Informally, our recovery guarantee can be stated as follows,

######  Theorem 2.1 (Informal Recovery Guarantee).

Consider a recommender system with a realizable model Eq. (2) with sigmoid
activation, Assume the features $\\{x_{i}\\}_{i\in[n_{1}]}$ and
$\\{y_{j}\\}_{j\in[n_{2}]}$ are sampled i.i.d. from the normal distribution
and the observed pairs $\Omega$ are i.i.d. sampled from
$\\{x_{i}\\}_{i\in[n_{1}]}\times\\{y_{j}\\}_{j\in[n_{2}]}$ uniformly at
random. Then there exists an algorithm such that $U^{*},V^{*}$ can be
recovered to any precision $\epsilon$ with time complexity and sample
complexity (refers to $n_{1},n_{2},|\Omega|$) polynomial in the dimension and
the condition number of $U^{*},V^{*}$, and logarithmic in $1/\epsilon$.

##  3 Main Results

Our main result shows that when initialized properly, gradient-based
algorithms will be guaranteed to converge to the ground truth. We first study
the Hessian of empirical risk for different activation functions, then based
on the positive-definiteness of the Hessian for smooth activations, we show
local linear convergence of gradient descent. The proof sketch is provided in
Appendix C.

The positive definiteness of the Hessian does not hold for several activation
functions. Here we provide some examples. Counter Example 1) The Hessian at
the ground truth for linear activation is not positive definite because for
any full-rank matrix $R\in\mathbb{R}^{k\times k}$, $(U^{*}R,V^{*}R^{-1})$ is
also a global optimal. Counter Example 2) The Hessian at the ground truth for
ReLU activation is not positive definite because for any diagonal matrix
$D\in\mathbb{R}^{k\times k}$ with positive diagonal elements,
$U^{*}D,V^{*}D^{-1}$ is also a global optimal. These counter examples have a
common property: there is redundancy in the parameters. Surprisingly, for
sigmoid and tanh, the Hessian around the ground truth is positive definite.
More surprisingly, we will later show that for ReLU, if the parameter space is
constrained properly, its Hessian at a given point near the ground truth can
also be proved to be positive definite with high probability.

###  3.1 Local Geometry and Local Linear Convergence for Sigmoid and Tanh

We define two natural condition numbers for the problem that captures the
"hardness" of the problem:

######  Definition 3.1.

Define $\lambda:=\max\\{\lambda(U^{*}),\lambda(V^{*})\\}$ and
$\kappa:=\max\\{\kappa(U^{*}),\kappa(V^{*})\\}$, where
$\lambda(U)=\sigma_{1}^{k}(U)/(\Pi_{i=1}^{k}\sigma_{i}(U))$,
$\kappa(U)=\sigma_{1}(U)/\sigma_{k}(U)$, and $\sigma_{i}(U)$ denotes the
$i$-th singular value of $U$ with the ordering $\sigma_{i}\geq\sigma_{i+1}$.

First we show the result for sigmoid and tanh activations.

######  Theorem 3.2 (Positive Definiteness of Hessian for Sigmoid and Tanh).

Let the activation function $\phi$ in the NIMC model (2) be sigmoid or tanh
and let $\kappa,\lambda$ be as defined in Definition 3.1. Then for any $t>1$
and any given $U,V$, if

$\displaystyle\quad n_{1}\gtrsim t\lambda^{4}\kappa^{2}d\log^{2}d,\quad
n_{2}\gtrsim t\lambda^{4}\kappa^{2}d\log^{2}d,\quad|\Omega|\gtrsim
t\lambda^{4}\kappa^{2}d\log^{2}d,$ $\displaystyle\quad\text{and
}\quad\|U-U^{*}\|+\|V-V^{*}\|\lesssim 1/(\lambda^{2}\kappa),$

then with probability at least $1-d^{-t}$, the smallest eigenvalue of the
Hessian of Eq. (3) is lower bounded by:

$\displaystyle\lambda_{\min}(\nabla^{2}f_{\Omega}(U,V))\gtrsim
1/(\lambda^{2}\kappa).$

Remark. Theorem 3.2 shows that, given sufficiently large number of user-items
ratings and a sufficiently large number of users/items themselves, the Hessian
at a point close enough to the true parameters $U^{*}$, $V^{*}$, is positive
definite with high probability. The sample complexity, including $n_{1},n_{2}$
and $|\Omega|$, have a near-linear dependency on the dimension, which matches
the linear IMC analysis [JD13]. Strong convexity parameter as well as the
sample complexity depend on the condition number of $U^{*},V^{*}$ as defined
in Definition 3.1. Although we don’t explicitly show the dependence on $k$,
both sample complexity and the minimal eigenvalue scale as a polynomial of
$k$. The proofs can be found in Appendix C.

As the above theorem shows the Hessian is positive definite w.h.p. for a given
$U,V$ that is close to the optima. This result along with smoothness of the
activation function implies linear convergence of gradient descent that
samples a fresh batch of samples in each iteration as shown in the following,
whose proof is postponed to Appendix E.1.

######  Theorem 3.3.

Let $[U^{c},V^{c}]$ be the parameters in the $c$-th iteration. Assuming
$\|U^{c}-U^{*}\|+\|V^{c}-V^{*}\|\lesssim 1/(\lambda^{2}\kappa)$, then given a
fresh sample set, $\Omega$, that is independent of $[U^{c},V^{c}]$ and
satisfies the conditions in Theorem 3.2, the next iterate using one step of
gradient descent, i.e., $[U^{c+1},V^{c+1}]=[U^{c},V^{c}]-\eta\nabla
f_{\Omega}(U^{c},V^{c}),$ satisfies

$\displaystyle\|U^{c+1}-U^{*}\|_{F}^{2}+\|V^{c+1}-V^{*}\|_{F}^{2}\leq(1-M_{l}/M_{u})(\|U^{c}-U^{*}\|_{F}^{2}+\|V^{c}-V^{*}\|_{F}^{2})$

with probability $1-d^{-t}$, where $\eta=\Theta(1/M_{u})$ is the step size and
$M_{l}\gtrsim 1/(\lambda^{2}\kappa)$ is the lower bound on the eigenvalues of
the Hessian and $M_{u}\lesssim 1$ is the upper bound on the eigenvalues of the
Hessian.

Remark. The linear convergence requires each iteration has a set of fresh
samples. However, since it converges linearly to the ground-truth, we only
need $\log(1/\epsilon)$ iterations, therefore the sample complexity is only
logarithmic in $1/\epsilon$. This dependency is better than directly using
Tensor decomposition method [JSA15], which requires $O(1/\epsilon^{2})$
samples. Note that we only use Tensor decomposition to initialize the
parameters. Therefore the sample complexity required in our tensor
initialization does not depend on $\epsilon$.

###  3.2 Empirical Hessian around the Ground Truth for ReLU

We now present our result for ReLU activation. As we see in Counter Example 2,
without any further modification, the Hessian for ReLU is not locally strongly
convex due to the redundancy in parameters. Therefore, we reduce the parameter
space by fixing one parameter for each $(u_{i},v_{i})$ pair, $i\in[k]$. In
particular, we fix $u_{1,i}=u_{1,i}^{*},\forall i\in[k]$ when minimizing the
objective function, Eq. (3), where $u_{1,i}$ is $i$-th element in the first
row of $U$. Note that as long as $u_{1,i}^{*}\neq 0$, $u_{1,i}$ can be fixed
to any other non-zero values. We set $u_{1,i}=u_{1,i}^{*}$ just for simplicity
of the proof. The new objective function can be represented as

$\displaystyle
f_{\Omega}^{\mathrm{ReLU}}(W,V)=\frac{1}{2|\Omega|}\sum_{(x,y)\in\Omega}(\phi(W^{\top}x_{2:d}+x_{1}(u^{*(1)})^{\top})^{\top}\phi(V^{\top}y)-A(x,y))^{2}.$
(4)

where $u^{*(1)}$ is the first row of $U^{*}$ and $W\in\mathbb{R}^{(d-1)\times
k}$.

Surprisingly, after fixing one parameter for each $(u_{i},v_{i})$ pair, the
Hessian using ReLU is also positive definite w.h.p. for a given $(U,V)$ around
the ground truth.

######  Theorem 3.4 (Positive Definiteness of Hessian for ReLU).

Define $u_{0}:=\min_{i\in[k]}\\{|u_{1,i}^{*}|\\}$. For any $t>1$ and any given
$U,V$, if

$\displaystyle\quad n_{1}\gtrsim
u_{0}^{-4}t\lambda^{4}\kappa^{12}d\log^{2}d,\quad n_{2}\gtrsim
u_{0}^{-4}t\lambda^{4}\kappa^{12}d\log^{2}d,\quad|\Omega|\gtrsim
u_{0}^{-4}t\lambda^{4}\kappa^{12}d\log^{2}d,$
$\displaystyle\quad\|W-W^{*}\|+\|V-V^{*}\|\lesssim
u_{0}^{4}/\lambda^{4}\kappa^{12},$

then with probability $1-d^{-t}$, the minimal eigenvalue of the objective for
ReLU activation function, Eq. (4), is lower bounded,

$\displaystyle\lambda_{\min}(\nabla^{2}f_{\Omega}^{\mathrm{ReLU}}(W,V))\gtrsim
u_{0}^{2}/\lambda^{2}\kappa^{4}.$

Remark. Similar to the sigmoid/tanh case, the sample complexity for ReLU case
also has a linear dependency on the dimension. However, here we have a worse
dependency on the condition number of the weight matrices. The scale of
$u_{0}$ can also be important and in practice one needs to set it carefully.
Note that although the activation function is not smooth, the Hessian at a
given point can still exist with probability $1$, since ReLU is smooth almost
everywhere and there are only a finite number of samples. However, owing to
the non-smoothness, a proof of convergence of gradient descent method for ReLU
is still an open problem.

###  3.3 Initialization

To achieve the ground truth, our algorithm needs a good initialization method
that can initialize the parameters to fall into the neighborhood of the ground
truth. Here we show that this is possible by using tensor method under the
Gaussian assumption.

In the following, we consider estimating $U^{*}$. Estimating $V^{*}$ is
similar.

Define $M_{3}:=\leavevmode\nobreak\ \mathbb{E}[A(x,y)\cdot(x^{\otimes
3}-x\widetilde{\otimes}I)]$, where
$x\widetilde{\otimes}I:=\sum_{j=1}^{d}[x\otimes e_{j}\otimes
e_{j}+e_{j}\otimes x\otimes e_{j}+e_{j}\otimes e_{j}\otimes x]$. Define
$\gamma_{j}(\sigma):=\mathbb{E}[\phi(\sigma\cdot z)z^{j}],\;\forall
j=0,1,2,3$. Then $M_{3}=\sum_{i=1}^{k}\alpha_{i}\overline{u}_{i}^{*\otimes
3},$ where $\overline{u}_{i}^{*}=u_{i}^{*}/\|u_{i}^{*}\|$ and
$\alpha_{i}=\gamma_{0}(\|v_{i}^{*}\|)\left(\gamma_{3}(\|u_{i}^{*}\|)-3\gamma_{1}(\|u_{i}^{*}\|)\right)$.
When $\alpha_{i}\neq 0$, we can approximately recover $\alpha_{i}$ and
$\overline{u}_{i}^{*}$ from the empirical version of $M_{3}$ using non-
orthogonal tensor decomposition [KCL15]. When $\phi$ is sigmoid,
$\gamma_{0}(\|v_{i}^{*}\|)=0.5$. Given $\alpha_{i}$, we can estimate
$\|u_{i}^{*}\|$, since $\alpha_{i}$ is a monotonic function w.r.t.
$\|u_{i}^{*}\|$. Applying Lemma B.7 in [ZSJ+17], we can bound the
approximation error of empirical $M_{3}$ and population $M_{3}$ using
polynomial number of samples. By [KCL15], we can bound the estimation error of
$\|u_{i}^{*}\|$ and $\overline{u}_{i}^{*}$. Finally combining Theorem 3.2, we
are able to show the recovery guarantee for sigmoid activation, i.e., Theorem
2.1.

Although tensor initialization has nice theoretical guarantees and sample
complexity, it heavily depends on Gaussian assumption and realizable model
assumption. In contrast, practitioners typically use random initialization.

##  4 Experiments

---

Dataset | $n$ | $d$ | $k$ | $|\Omega|$ | NIMC | IMC |  | NIMC-
---
RFF
| IMC-
---
RFF
mushroom | 8124 | 112 | 2 | 5$n$ | 0 | 0.0049 | 0 | 0
20$n$ | 0 | 0.0010 | 0 | 0
segment | 2310 | 19 | 7 | 5$n$ | 0.0543 | 0.0694 | 0.0197 | 0.0257
20$n$ | 0.0655 | 0.0768 | 0.0092 | 0.0183
covtype | 1708 | 54 | 7 | 5$n$ | 0.1671 | 0.1733 | 0.1548 | 0.1529
20$n$ | 0.1555 | 0.1600 | 0.1200 | 0.1307
letter | 15000 | 16 | 26 | 5$n$ | 0.0590 | 0.0704 | 0.0422 | 0.0430
20$n$ | 0.0664 | 0.0760 | 0.0321 | 0.0356
yalefaces | 2452 | 100 | 38 | 5$n$ | 0.0315 | 0.0329 | 0.0266 | 0.0273
20$n$ | 0.0212 | 0.0277 | 0.0064 | 0.0142
usps | 7291 | 256 | 10 | 5$n$ | 0.0211 | 0.0361 | 0.0301 | 0.0185
20$n$ | 0.0184 | 0.0320 | 0.0199 | 0.0152

Figure 1: The left two figures (top: sigmoid, bottom: ReLU) plot rate of
success of GD over synthetic data. White blocks denote $100\%$ success rate.
The right table presents error in semi-supervised clustering using NIMC and
IMC.

In this section, we show experimental results on both synthetic data and real-
world data. Our experiments on synthetic data are intended to verify our
theoretical analysis, while the real-world data shows the superior performance
of NIMC over IMC. We apply gradient descent with random initialization to both
NIMC and IMC.

###  4.1 Synthetic Data

We first generate some synthetic datasets to verify the sample complexity and
the convergence of gradient descent using random initialization. We fix
$k=5,d=10$. For sigmoid, set the number of samples $n_{1}=n_{2}=n=\\{10\cdot
i\\}_{i=1,2\cdots,10}$ and the number of observations $|\Omega|=m=\\{2kd\cdot
i\\}_{i=1,2,\cdots,10}$. For ReLU, set $n=\\{20\cdot i\\}_{i=1,2\cdots,10}$
and $m=\\{4kd\cdot i\\}_{i=1,2,\cdots,10}$. The sampling rule follows our
previous assumptions. For each $n,m$ pair, we make 5 trials and take the
average of the successful recovery times. We say a solution ($U,V$)
successfully recovers the ground truth parameters when the solution achieves
0.001 relative testing error, i.e.,
$\|\phi(X_{t}U)\phi(X_{t}U)^{\top}-\phi(X_{t}U^{*})\phi(X_{t}U^{*})^{\top}\|_{F}\leq
0.001\cdot\|\phi(X_{t}U^{*})\phi(X_{t}U^{*})^{\top}\|_{F},$ where
$X_{t}\in\mathbb{R}^{n\times d}$ is a newly sampled testing dataset. For both
ReLU and sigmoid, we minimize the original objective function (3). We
illustrate the recovery rate in left figures in Figure 1. As we can see, ReLU
requires more samples/observations than that for sigmoid for exact recovery
(note the scales of $n$ and $m/2kd$ are different in the two figures). This is
consistent with our theoretical results. Comparing Theorem 3.2 and Theorem
3.4, we can see the sample complexity for ReLU has a worse dependency on the
conditioning of $U^{*},V^{*}$ than sigmoid. We can also see that when $n$ is
sufficiently large, the number of observed ratings required remains the same
for both methods. This is also consistent with the theorems, where $|\Omega|$
is near-linear in $d$ and is independent of $n$.

###  4.2 Semi-supervised Clustering

We apply NIMC to semi-supervised clustering and follow the experimental
setting in GIMC [SCH+16]. In this problem, we are given a set of items with
their features, $X\in\mathbb{R}^{n\times d}$, where $n$ is the number of items
and $d$ is the feature dimension, and an incomplete similarity matrix $A$,
where $A_{i,j}=1$ if $i$-th item and $j$-th item are similar and $A_{i,j}=0$
if $i$-th item and $j$-th item are dissimilar. The goal is to do clustering
using both existing features and the partially observed similarity matrix. We
build the dataset from a classification dataset where the label of each item
is known and will be used as the ground truth cluster. We first compute the
similarity matrix from the labels and sample $|\Omega|$ entries uniformly as
the observed entries. Since there is only one features we set $y_{j}=x_{j}$ in
the objective function Eq. (3).

We initialize $U$ and $V$ to be the same Gaussian random matrix, then apply
gradient descent. This guarantees $U$ and $V$ to keep identical during the
optimization process. Once $U$ converges, we take the top $k$ left singular
vectors of $\phi(XU)$ to do k-means clustering. The clustering error is
defined as in [SCH+16]. Like [SCH+16], we define the clustering error as
follows,

$\text{error}=\frac{2}{n(n-1)}\left(\sum_{(i,j):\pi^{*}_{i}=\pi^{*}_{j}}1_{\pi_{i}\neq\pi_{j}}+\sum_{(i,j):\pi^{*}_{i}\neq\pi^{*}_{j}}1_{\pi_{i}=\pi_{j}}\right),$

where $\pi^{*}$ is the ground-truth clustering and $\pi$ is the predicted
clustering. We compare NIMC of a ReLU activation function with IMC on six
datasets using raw features and random Fourier features (RFF). The random
Fourier feature is
$r(x)=\frac{1}{\sqrt{q}}\cdot\begin{bmatrix}\sin(Qx)^{\top}&\cos(Qx)^{\top}\end{bmatrix}^{\top}\in\operatorname{{\mathbb{R}}}^{2q}$
and each entry of $Q\in\mathbb{R}^{q\times d}$ is i.i.d. sampled from
$\mathcal{N}(0,\sigma)$. We use Random Fourier features in order to see how
increasing the depth of the neural network changes the performance. However,
our analysis only works for one-layer neural networks, therefore, we use
Random Fourier features, which can be viewed as using two-layer neural
networks but with the first-layer parameters fixed.

$\sigma$ is chosen such that a linear classifier using these random features
achieves the best classification accuracy. $q$ is set as $100$ for all
datasets. Datasets mushroom, segment, letter,usps,covtype are downloaded from
libsvm website. We subsample covtype dataset to balance the samples from
different classes. We preprocess yalefaces dataset as described in [KTWA14].
As shown in the right table in Figure 1, when using raw features, NIMC
achieves better clustering results than IMC for all the cases. This is also
true for most cases when using Random Fourier features.

###  4.3 Recommendation Systems

Recommender systems are used in many real situations. Here we consider two
tasks.

Movie recommendation for users. We use Movielens[Res97] dataset, which has not
only the ratings users give movies but also the users’ demographic information
and movies’ genre information. Our goal is to predict ratings that new users
will give the existing movies. We randomly split the users into existing users
(training data) and new users (testing data) with ratio 4:1. The user features
include 21 types of occupations, 7 different age ranges and one gender
information; the movie features include 18-19 (18 for ml-1m and 19 for
ml-100k) genre features and 20 features from the top 20 right singular values
of the training rating matrix (which has size #training users -by- #movies).
In our experiments, we set $k$ to be 50. Here are our results on datasets
ml-1m and ml-100k. For NIMC, we use ReLU activation. As shown in Table 1, NIMC
achieves much smaller RMSE than IMC for both ml-100k and ml-1m datasets.

Dataset | #movies | #users | # ratings | # movie feat. | # user feat. |  | RMSE
---
NIMC
| RMSE
---
IMC
ml-100k | 1682 | 943 | 100,000 | 39 | 29 | 1.034 | 1.321
ml-1m | 3883 | 6040 | 1,000,000 | 38 | 29 | 1.021 | 1.320
Table 1: Test RMSE for recommending new users with movies on Movielens dataset.  |  |  |
---|---|---|---
(a) | (b) | (c) | (d)
Figure 2: NIMC v.s. IMC on gene-disease association prediction task.

Gene-Disease association prediction.  We use the dataset collected by [ND14],
which has 300 gene features and 200 disease features. Our goal is to predict
associated genes for a new disease given its features. Since the dataset only
contains positive labels, this is a problem called positive-unlabeled learning
[HND15] or one-class matrix factorization [YHDL17]. We adapt our objective
function to the following objective,

$f(U,V)=\frac{1}{2}\left(\sum_{(i,j)\in\Omega}(\phi(U^{\top}x_{i})^{\top}\phi(V^{\top}y_{j})-A_{ij})^{2}+\beta\sum_{(i,j)\in\Omega^{c}}(\phi(U^{\top}x_{i})^{\top}\phi(V^{\top}y_{j}))^{2}\right),$
(5)

where $A$ is the association matrix, $\Omega$ is the set of indices for
observed associations, $\Omega^{c}$ is the complementary set of $\Omega$ and
$\beta$ is the penalty weight for unobserved associations. There are totally
12331 genes and 3209 diseases in the dataset. We randomly split the diseases
into training diseases and testing diseases with ratio 4:1. The results are
presented in Fig 2. We follow [ND14] and use the cumulative distribution of
the ranks as a measure for comparing the performances of different methods,
i.e., the probability that any ground-truth associated gene of a disease
appears in the retrieved top-$r$ genes for this disease.

In Fig 2(a), we show how $k$ changes the performance of NIMC and IMC. In
general, the higher $k$, the better the performance. The performance of IMC
becomes stable when $k$ is larger than 100, while the performance of NIMC is
still increasing. Although IMC performs better than NIMC when $k$ is small,
the performance of NIMC increases much faster than IMC when $k$ increases.
$\beta$ is fixed as $0.01$ and $r=100$ in the experiment for Fig 2(a). In Fig.
2(b), we present how $\beta$ in Eq. (5) affects the performance. We tried over
$\beta=[10^{-4},10^{-3},10^{-2},10^{-1},1]$ to check how the value of $\beta$
changes the performance. As we can see, $\beta=10^{-3}$ and $10^{-2}$ give the
best results. Fig. 2(c) shows the probability that any ground-truth associated
gene of a disease appears in the retrieved top-$r$ genes for this disease
w.r.t. different $r$’s. Here we fix $k=200$, and $\beta=0.01$. Fig. 2(d) shows
the precision-recall curves for different methods when $k=200$, and
$\beta=0.01$.

##  5 Conclusion

In this paper, we studied a nonlinear IMC model that represents one of the
simplest inductive model for neural-network-based recommender systems. We
study local geometry of the empirical risk function and show that, close to
the optima, the function is strongly convex for both ReLU and sigmoid
activations. Therefore, using a smooth activation function like sigmoid
activation along with standard tensor initialization, gradient descent
recovers the underlying model with polynomial sample complexity and time
complexity. Thus we provide the first theoretically rigorous result for the
non-linear recommendation system problem, which we hope will spur further
progress in the area of deep-learning based recommendation systems. Our
experimental results on synthetic data matches our analysis and the results on
real-world benchmarks for semi-supervised clustering and recommendation
systems show a superior performance over linear IMC.

## References

  * [ABEV06] Jacob Abernethy, Francis Bach, Theodoros Evgeniou, and Jean-Philippe Vert.  Low-rank matrix factorization with attributes.  arXiv preprint cs/0611124, 2006.
  * [BG17] Alon Brutzkus and Amir Globerson.  Globally optimal gradient descent for a convnet with Gaussian inputs.  In ICML. https://arxiv.org.pdf/1702.07966, 2017.
  * [BGMSS18] Alon Brutzkus, Amir Globerson, Eran Malach, and Shai Shalev-Shwartz.  SGD learns over-parameterized networks that provably generalize on linearly separable data.  In ICLR. https://arxiv.org/pdf/1710.10174, 2018.
  * [CAS16] Paul Covington, Jay Adams, and Emre Sargin.  Deep neural networks for YouTube recommendations.  In Proceedings of the 10th ACM Conference on Recommender Systems, pages 191–198. ACM, 2016.
  * [CHD15] Kai-Yang Chiang, Cho-Jui Hsieh, and Inderjit S Dhillon.  Matrix completion with noisy side information.  In Advances in Neural Information Processing Systems, pages 3447–3455, 2015.
  * [CR09] Emmanuel J. Candès and Benjamin Recht.  Exact matrix completion via convex optimization.  Foundations of Computational Mathematics, 9(6):717–772, December 2009.
  * [DLT18a] Simon S Du, Jason D Lee, and Yuandong Tian.  When is a convolutional filter easy to learn?  In ICLR. https://arxiv.org/pdf/1709.06129, 2018.
  * [DLT+18b] Simon S Du, Jason D Lee, Yuandong Tian, Barnabas Poczos, and Aarti Singh.  Gradient descent learns one-hidden-layer CNN: Don’t be afraid of spurious local minima.  In ICML. https://arxiv.org/pdf/1712.00779, 2018.
  * [GJZ17] Rong Ge, Chi Jin, and Yi Zheng.  No spurious local minima in nonconvex low rank problems: A unified geometric analysis.  In International Conference on Machine Learning, pages 1233–1242. https://arxiv.org/pdf/1704.00708, 2017.
  * [GKKT17] Surbhi Goel, Varun Kanade, Adam Klivans, and Justin Thaler.  Reliably learning the ReLU in polynomial time.  In 30th Annual Conference on Learning Theory (COLT). https://arxiv.org/pdf/1611.10258, 2017.
  * [GKM18] Surbhi Goel, Adam Klivans, and Reghu Meka.  Learning one convolutional layer with overlapping patches.  In ICML. https://arxiv.org/pdf/1802.02547, 2018.
  * [GLM18] Rong Ge, Jason D Lee, and Tengyu Ma.  Learning one-hidden-layer neural networks with landscape design.  In ICLR. https://arxiv.org/pdf/1711.00501, 2018.
  * [GUH16] Carlos A Gomez-Uribe and Neil Hunt.  The Netflix recommender system: Algorithms, business value, and innovation.  ACM Transactions on Management Information Systems (TMIS), 6(4):13, 2016.
  * [Har14] Moritz Hardt.  Understanding alternating minimization for matrix completion.  In Foundations of Computer Science (FOCS), 2014 IEEE 55th Annual Symposium on, pages 651–660. IEEE, 2014.
  * [HKZ12] Daniel Hsu, Sham M Kakade, and Tong Zhang.  A tail inequality for quadratic forms of subGaussian random vectors.  Electronic Communications in Probability, 17(52):1–6, 2012.
  * [HND15] Cho-Jui Hsieh, Nagarajan Natarajan, and Inderjit Dhillon.  Pu learning for matrix completion.  In International Conference on Machine Learning, pages 2445–2453, 2015.
  * [JD13] Prateek Jain and Inderjit S Dhillon.  Provable inductive matrix completion.  arXiv preprint arXiv:1306.0626, 2013.
  * [JNS13] Prateek Jain, Praneeth Netrapalli, and Sujay Sanghavi.  Low-rank matrix completion using alternating minimization.  In STOC, 2013.
  * [JSA15] Majid Janzamin, Hanie Sedghi, and Anima Anandkumar.  Beating the perils of non-convexity: Guaranteed training of neural networks using tensor methods.  arXiv preprint 1506.08473, 2015.
  * [KCL15] Volodymyr Kuleshov, Arun Chaganty, and Percy Liang.  Tensor factorization via matrix factorization.  In AISTATS, pages 507–516, 2015.
  * [KTWA14] Matt Kusner, Stephen Tyree, Kilian Weinberger, and Kunal Agrawal.  Stochastic neighbor compression.  In International Conference on Machine Learning, pages 622–630, 2014.
  * [LLL+16] Chenyi Lei, Dong Liu, Weiping Li, Zheng-Jun Zha, and Houqiang Li.  Comparative deep learning of hybrid representations for image recommendations.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 2545–2553, 2016.
  * [LY16] Ming Lin and Jieping Ye.  A non-convex one-pass framework for generalized factorization machine and rank-one matrix sensing.  In Advances in Neural Information Processing Systems, pages 1633–1641, 2016.
  * [LY17] Yuanzhi Li and Yang Yuan.  Convergence analysis of two-layer neural networks with ReLU activation.  In Advances in Neural Information Processing Systems, pages 597–607. https://arxiv.org/pdf/1705.09886, 2017.
  * [ND14] Nagarajan Natarajan and Inderjit S Dhillon.  Inductive matrix completion for predicting gene–disease associations.  Bioinformatics, 30(12):i60–i68, 2014.
  * [Ren10] Steffen Rendle.  Factorization machines.  In Data Mining (ICDM), 2010 IEEE 10th International Conference on, pages 995–1000. IEEE, 2010.
  * [Res97] GroupLens Research.  Movie lens dataset.  In University of Minnesota. http://www.grouplens.org/taxonomy/term/14, 1997.
  * [RSW16] Ilya Razenshteyn, Zhao Song, and David P Woodruff.  Weighted low rank approximations with provable guarantees.  In Proceedings of the forty-eighth annual ACM symposium on Theory of Computing (STOC), pages 250–263. ACM, 2016.
  * [SCH+16] Si Si, Kai-Yang Chiang, Cho-Jui Hsieh, Nikhil Rao, and Inderjit S Dhillon.  Goal-directed inductive matrix completion.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 1165–1174. ACM, 2016.
  * [SCLD15] Donghyuk Shin, Suleyman Cetintas, Kuang-Chih Lee, and Inderjit S Dhillon.  Tumblr blog recommendation with boosted inductive matrix completion.  In Proceedings of the 24th ACM International on Conference on Information and Knowledge Management, pages 203–212. ACM, 2015.
  * [SWZ17a] Zhao Song, David P Woodruff, and Peilin Zhong.  Low rank approximation with entrywise $\ell_{1}$-norm error.  In Proceedings of the 49th Annual Symposium on the Theory of Computing (STOC). ACM, https://arxiv.org/pdf/1611.00898, 2017.
  * [SWZ17b] Zhao Song, David P Woodruff, and Peilin Zhong.  Relative error tensor low rank approximation.  arXiv preprint arXiv:1704.08246, 2017.
  * [SWZ18] Zhao Song, David P Woodruff, and Peilin Zhong.  Towards a zero-one law for entrywise low rank approximation.  2018\.
  * [Tia17] Yuandong Tian.  An analytical formula of population gradient for two-layered ReLU network and its applications in convergence and critical point analysis.  In ICML. https://arxiv.org/pdf/1703.00560, 2017.
  * [Tro12] Joel A. Tropp.  User-friendly tail bounds for sums of random matrices.  Foundations of Computational Mathematics, 12(4):389–434, 2012.
  * [WW14] Xinxi Wang and Ye Wang.  Improving content-based and hybrid music recommendation using deep learning.  In Proceedings of the 22nd ACM international conference on Multimedia, pages 627–636. ACM, 2014.
  * [XJZ13] Miao Xu, Rong Jin, and Zhi-Hua Zhou.  Speedup matrix completion with side information: Application to multi-label learning.  In NIPS, pages 2301–2309, 2013.
  * [YHDL17] Hsiang-Fu Yu, Hsin-Yuan Huang, Inderjit S Dihillon, and Chih-Jen Lin.  A unified algorithm for one-class structured matrix factorization with side information.  In AAAI, 2017.
  * [YJKD14] Hsiang-Fu Yu, Prateek Jain, Purushottam Kar, and Inderjit Dhillon.  Large-scale multi-label learning with missing labels.  In ICML, pages 593–601, 2014.
  * [ZDG18] Xiao Zhang, Simon S Du, and Quanquan Gu.  Fast and sample efficient inductive matrix completion via multi-phase procrustes flow.  In ICML. https://arxiv.org/pdf/1803.01233, 2018.
  * [ZJD15] Kai Zhong, Prateek Jain, and Inderjit S. Dhillon.  Efficient matrix sensing using rank-1 Gaussian measurements.  In International Conference on Algorithmic Learning Theory, pages 3–18. Springer, 2015.
  * [ZSD17] Kai Zhong, Zhao Song, and Inderjit S Dhillon.  Learning non-overlapping convolutional neural networks with multiple kernels.  arXiv preprint arXiv:1711.03440, 2017.
  * [ZSJ+17] Kai Zhong, Zhao Song, Prateek Jain, Peter L Bartlett, and Inderjit S Dhillon.  Recovery guarantees for one-hidden-layer neural networks.  In ICML. https://arxiv.org/pdf/1706.03175, 2017.
  * [ZYL+16] Fuzheng Zhang, Nicholas Jing Yuan, Defu Lian, Xing Xie, and Wei-Ying Ma.  Collaborative knowledge base embedding for recommender systems.  In Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, pages 353–362. ACM, 2016.
  * [ZYS17] Shuai Zhang, Lina Yao, and Aixin Sun.  Deep learning based recommender system: A survey and new perspectives.  arXiv preprint arXiv:1707.07435, 2017.

## Appendix

##  Appendix A Notation

For any positive integer $n$, we use $[n]$ to denote the set
$\\{1,2,\cdots,n\\}$. For random variable $X$, let $\mathbb{E}[X]$ denote the
expectation of $X$ (if this quantity exists).

For any vector $x\in\mathbb{R}^{n}$, we use $\|x\|$ to denote its $\ell_{2}$
norm.

We provide several definitions related to matrix $A$. Let $\det(A)$ denote the
determinant of a square matrix $A$. Let $A^{\top}$ denote the transpose of
$A$. Let $A^{\dagger}$ denote the Moore-Penrose pseudoinverse of $A$. Let
$A^{-1}$ denote the inverse of a full rank square matrix. Let $\|A\|_{F}$
denote the Frobenius norm of matrix $A$. Let $\|A\|$ denote the spectral norm
of matrix $A$. Let $\sigma_{i}(A)$ to denote the $i$-th largest singular value
of $A$.

We use ${\mathbf{1}}_{f}$ to denote the indicator function, which is $1$ if
$f$ holds and $0$ otherwise. Let $I_{d}\in\mathbb{R}^{d\times d}$ denote the
identity matrix. We use $\phi(z)$ to denote an activation function. We use
${\cal D}$ to denote a Gaussian distribution ${\cal N}(0,I_{d})$. For integer
$k$, we use $\operatorname*{{\mathcal{D}}}_{k}$ to denote
$\operatorname*{{\mathcal{N}}}(0,I_{k})$.

For any function $f$, we define $\widetilde{O}(f)$ to be
$f\cdot\log^{O(1)}(f)$. In addition to $O(\cdot)$ notation, for two functions
$f,g$, we use the shorthand $f\lesssim g$ (resp. $\gtrsim$) to indicate that
$f\leq Cg$ (resp. $\geq$) for an absolute constant $C$. We use $f\eqsim g$ to
mean $cf\leq g\leq Cf$ for constants $c,C$.

##  Appendix B Preliminaries

We state some useful facts in this section.

######  Fact B.1.

Let $A=\begin{bmatrix}a_{1}&a_{2}&\cdots&a_{k}\end{bmatrix}$. Let
$\operatorname{diag}(A)\in\operatorname{{\mathbb{R}}}^{k}$ denote the vector
where the $i$-th entry is $A_{i,i}$, $\forall i\in[k]$. Let ${\bf
1}\in\operatorname{{\mathbb{R}}}^{k}$ denote the vector that the $i$-th entry
is $1$, $\forall i\in[k]$. We have the following properties,

$\displaystyle\mathrm{(\@slowromancap i@)}$ $\displaystyle\leavevmode\nobreak\
\sum_{i=1}^{k}(a_{i}^{\top}e_{i})^{2}=\|\operatorname{diag}(A)\|_{2}^{2},$
$\displaystyle\mathrm{(\@slowromancap ii@)}$
$\displaystyle\leavevmode\nobreak\
\sum_{i=1}^{k}(a_{i}^{\top}a_{i})^{2}=\|A\|_{F}^{2},$
$\displaystyle\mathrm{(\@slowromancap iii@)}$
$\displaystyle\leavevmode\nobreak\
\sum_{i=1}^{k}\sum_{j=1}^{k}(a_{i}^{\top}a_{j})=\|A\cdot{\bf 1}\|_{2}^{2},$
$\displaystyle\mathrm{(\@slowromancap iv@)}$
$\displaystyle\leavevmode\nobreak\ \sum_{i\neq
j}a_{i}^{\top}a_{j}=\|A\cdot{\bf 1}\|_{2}^{2}-\|A\|_{F}^{2}.$

###### Proof.

Using the definition, it is easy to see that (i@), (ii@) and (iii@) are
holding.

Proof of (iv@), we have

$\displaystyle\sum_{i\neq
j}a_{i}^{\top}a_{j}=\sum_{i,j}a_{i}^{\top}a_{j}-\sum_{i=1}^{k}a_{i}^{\top}a_{i}=\|A\cdot{\bf
1}\|_{2}^{2}-\|A\|_{F}^{2}.$

where the last step follows by (ii@) and (iii@). ∎

######  Fact B.2.

Let $A=\begin{bmatrix}a_{1}&a_{2}&\cdots&a_{k}\end{bmatrix}$. Let
$\operatorname{diag}(A)\in\operatorname{{\mathbb{R}}}^{k}$ denote the vector
where the $i$-th entry is $A_{i,i}$, $\forall i\in[k]$. Let ${\bf
1}\in\operatorname{{\mathbb{R}}}^{k}$ denote the vector that the $i$-th entry
is $1$, $\forall i\in[k]$. We have the following properties,

$\displaystyle\mathrm{(\@slowromancap i@)}$ $\displaystyle\leavevmode\nobreak\
\sum_{i\neq
j}a_{i}^{\top}e_{i}e_{i}^{\top}a_{j}=(\operatorname{diag}(A)^{\top}\cdot(A\cdot{\bf
1}))-\|\operatorname{diag}(A)\|_{2}^{2},$
$\displaystyle\mathrm{(\@slowromancap ii@)}$
$\displaystyle\leavevmode\nobreak\ \sum_{i\neq
j}a_{i}^{\top}e_{j}e_{j}^{\top}a_{j}=(\operatorname{diag}(A)^{\top}\cdot(A\cdot{\bf
1}))-\|\operatorname{diag}(A)\|_{2}^{2},$
$\displaystyle\mathrm{(\@slowromancap iii@)}$
$\displaystyle\leavevmode\nobreak\ \sum_{i\neq
j}a_{i}^{\top}e_{i}a_{j}^{\top}e_{j}=(\operatorname{diag}(A)^{\top}\cdot{\bf
1})^{2}-\|\operatorname{diag}(A)\|_{2}^{2},$
$\displaystyle\mathrm{(\@slowromancap iv@)}$
$\displaystyle\leavevmode\nobreak\ \sum_{i\neq
j}a_{i}^{\top}e_{j}a_{j}^{\top}e_{i}=\langle
A^{\top},A\rangle-\|\operatorname{diag}(A)\|_{2}^{2}.$

###### Proof.

Proof of (i@). We have

$\displaystyle\sum_{i\neq j}a_{i}^{\top}e_{i}e_{i}^{\top}a_{j}=$
$\displaystyle\leavevmode\nobreak\
\sum_{i,j}a_{i}^{\top}e_{i}e_{i}^{\top}a_{j}-\sum_{i=1}^{k}a_{i}^{\top}e_{i}e_{i}^{\top}a_{i}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\sum_{i,j}a_{i,i}e_{i}^{\top}a_{j}-\|\operatorname{diag}(A)\|_{2}^{2}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\sum_{i=1}^{k}a_{i,i}e_{i}^{\top}\sum_{j=1}^{k}a_{j}-\|\operatorname{diag}(A)\|_{2}^{2}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
(\operatorname{diag}(A)^{\top}\cdot(A\cdot{\bf
1}))-\|\operatorname{diag}(A)\|_{2}^{2}$

Proof of (ii@). It is similar to (i@).

Proof of (iii@). We have

$\displaystyle\sum_{i\neq j}a_{i}^{\top}e_{i}a_{j}^{\top}e_{j}=$
$\displaystyle\leavevmode\nobreak\
\sum_{i,j}a_{i}^{\top}e_{i}a_{j}^{\top}e_{j}-\sum_{i=1}a_{i}^{\top}e_{i}a_{i}^{\top}e_{i}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\sum_{i=1}^{k}a_{i}^{\top}e_{i}\cdot\sum_{j=1}^{k}a_{j}^{\top}e_{j}-\sum_{i=1}^{k}a_{i}^{\top}e_{i}a_{i}^{\top}e_{i}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\sum_{i=1}^{k}a_{i,i}\cdot\sum_{j=1}^{k}a_{j,j}-\sum_{i=1}^{k}a_{i,i}a_{i,i}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
(\operatorname{diag}(A)^{\top}\cdot{\bf
1})^{2}-\|\operatorname{diag}(A)\|_{2}^{2}$

Proof of (iv@). We have

$\displaystyle\sum_{i\neq j}a_{i}^{\top}e_{j}a_{j}^{\top}e_{i}=$
$\displaystyle\leavevmode\nobreak\ \sum_{i\neq
j}\operatorname{tr}[a_{i}^{\top}e_{j}a_{j}^{\top}e_{i}]$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\ \sum_{i\neq
j}\operatorname{tr}[e_{j}a_{j}^{\top}e_{i}a_{i}^{\top}]$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\ \sum_{i\neq j}\langle
e_{j}a_{j}^{\top},a_{i}e_{i}^{\top}\rangle$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\ \sum_{i,j}\langle
e_{j}a_{j}^{\top},a_{i}e_{i}^{\top}\rangle-\sum_{i=1}^{k}\langle
e_{i}a_{i}^{\top},a_{i}e_{i}^{\top}\rangle$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\ \langle
A^{\top},A\rangle-\|\operatorname{diag}(A)\|_{2}^{2}.$

where the second step follows by
$\operatorname{tr}[ABCD]=\operatorname{tr}[BCDA]$, the third step follows by
$\operatorname{tr}[AB]=\langle A,B^{\top}\rangle$. ∎

##  Appendix C Proof Sketch

At high level the proofs for Theorem 3.2 and Theorem 3.4 include the following
steps. 1) Show that the population Hessian at the ground truth is positive
definite. 2) Show that population Hessians near the ground truth are also
positive definite. 3) Employ matrix Bernstein inequality to bound the
population Hessian and the empirical Hessian.

We now formulate the Hessian. The Hessian of Eq. (3),
$\nabla^{2}f_{\Omega}(U,V)\in\mathbb{R}^{(2kd)\times(2kd)}$, can be decomposed
into two types of blocks, ($i\in[k],j\in[k]$),

$\frac{\partial^{2}f_{\Omega}(U,V)}{\partial u_{i}\partial
v_{j}},\frac{\partial^{2}f_{\Omega}(U,V)}{\partial u_{i}\partial u_{j}},$

where $u_{i}$($v_{j}$, resp.) is the $i$-th column of $U$ ($j$-th column of
$V$, resp.). Note that each of the above second-order derivatives is a
$d\times d$ matrix.

The first type of blocks are given by:

$\displaystyle\frac{\partial^{2}f_{\Omega}(U,V)}{\partial u_{i}\partial
v_{j}}=\underset{\Omega}{\widehat{\operatorname*{{\mathbb{E}}}}}\left[\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(v_{j}^{\top}y)xy^{\top}\phi(v_{i}^{\top}y)\phi(u_{j}^{\top}x)\right]+\delta_{ij}\underset{\Omega}{\widehat{\operatorname*{{\mathbb{E}}}}}\left[h_{x,y}(U,V)\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(v_{i}^{\top}y)xy^{\top}\right],$

where
$\widehat{\operatorname*{{\mathbb{E}}}}_{\Omega}[\cdot]=\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}[\cdot]$,
$\delta_{ij}=1_{i=j}$, and

$h_{x,y}(U,V)=\phi(U^{\top}x)^{\top}\phi(V^{\top}y)-\phi(U^{*\top}x)^{\top}\phi(V^{*\top}y).$

For sigmoid/tanh activation function, the second type of blocks are given by:

$\displaystyle\frac{\partial^{2}f_{\Omega}(U,V)}{\partial u_{i}\partial
u_{j}}=\underset{\Omega}{\widehat{\operatorname*{{\mathbb{E}}}}}\left[\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(u_{j}^{\top}x)xx^{\top}\phi(v_{i}^{\top}y)\phi(v_{j}^{\top}y)\right]+\delta_{ij}\underset{\Omega}{\widehat{\operatorname*{{\mathbb{E}}}}}\left[h_{x,y}(U,V)\phi^{\prime\prime}(u_{i}^{\top}x)\phi(v_{i}^{\top}y)xx^{\top}\right].$
(6)

For ReLU/leaky ReLU activation function, the second type of blocks are given
by:

$\displaystyle\frac{\partial^{2}f_{\Omega}(U,V)}{\partial u_{i}\partial
u_{j}}=\underset{\Omega}{\widehat{\operatorname*{{\mathbb{E}}}}}\left[\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(u_{j}^{\top}x)xx^{\top}\phi(v_{i}^{\top}y)\phi(v_{j}^{\top}y)\right].$

Note that the second term of Eq. (6) is missing here as $(U,V)$ are fixed, the
number of samples is finite and $\phi^{\prime\prime}(z)=0$ almost everywhere.

In this section, we will discuss important lemmas/theorems for Step 1 in
Appendix C.1 and Step 2,3 in Appendix C.3.

###  C.1 Positive definiteness of the population hessian

The corresponding population risk for Eq. (3) is given by:

$f_{\mathcal{D}}(U,V)=\frac{1}{2}\operatorname*{{\mathbb{E}}}_{(x,y)\sim\mathcal{D}}[(\phi(U^{\top}x)^{\top}\phi(V^{\top}y)-A(x,y))^{2}],$
(7)

where $\mathcal{D}:=\mathcal{X}\times\mathcal{Y}$. For simplicity, we also
assume $\mathcal{X}$ and $\mathcal{Y}$ are normal distributions.

Now we study the Hessian of the population risk at the ground truth. Let the
Hessian of $f_{\mathcal{D}}(U,V)$ at the ground-truth $(U,V)=(U^{*},V^{*})$ be
$H^{*}\in\mathbb{R}^{(2dk)\times(2dk)}$, which can be decomposed into the
following two types of blocks ($i\in[k],j\in[k]$),

$\displaystyle\frac{\partial^{2}f_{\mathcal{D}}(U^{*},V^{*})}{\partial
u_{i}\partial u_{j}}=$
$\displaystyle\operatorname*{{\mathbb{E}}}_{x,y}\left[\phi^{\prime}(u_{i}^{*\top}x)\phi^{\prime}(u_{j}^{*\top}x)xx^{\top}\phi(v_{i}^{*\top}y)\phi(v_{j}^{{}^{*}\top}y)\right],$
$\displaystyle\frac{\partial^{2}f_{\mathcal{D}}(U^{*},V^{*})}{\partial
u_{i}\partial v_{j}}=$
$\displaystyle\operatorname*{{\mathbb{E}}}_{x,y}\left[\phi^{\prime}(u_{i}^{*\top}x)\phi^{\prime}(v_{j}^{*\top}y)xy^{\top}\phi(v_{i}^{*\top}y)\phi(u_{j}^{{}^{*}\top}x)\right].$

To study the positive definiteness of $H^{*}$, we characterize the minimal
eigenvalue of $H^{*}$ by a constrained optimization problem,

$\displaystyle\lambda_{\min}(H^{*})=\min_{(a,b)\in\mathbb{B}}\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(u_{i}^{*\top}x)\phi(v_{i}^{*\top}y)x^{\top}a_{i}+\phi^{\prime}(v_{i}^{*\top}y)\phi(u_{i}^{*\top}x)y^{\top}b_{i}\right)^{2}\right],$
(8)

where $(a,b)\in\mathbb{B}$ denotes that
$\sum_{i=1}^{k}\|a_{i}\|^{2}+\|b_{i}\|^{2}=1$. Obviously,
$\lambda_{\min}(H^{*})\geq 0$ due to the squared loss and the realizable
assumption. However, this is not sufficient for the local convexity around the
ground truth, which requires the positive (semi-)definiteness for the
neighborhood around the ground truth. In other words, we need to show that
$\lambda_{\min}(H^{*})$ is strictly greater than $0$, so that we can
characterize an area in which the Hessian still preserves positive
definiteness (PD) despite the deviation from the ground truth.

Challenges. As we mentioned previously there are activation functions that
lead to redundancy in parameters. Hence one challenge is to distill properties
of the activation functions that preserve the PD. Another challenge is the
correlation introduced by $U^{*}$ when it is non-orthogonal. So we first study
the minimal eigenvalue for orthogonal $U^{*}$ and orthogonal $V^{*}$ and then
link the non-orthogonal case to the orthogonal case.

###  C.2 Warm up: orthogonal case

In this section, we consider the case when $U^{*},V^{*}$ are unitary matrices,
i.e., $U^{*\top}U^{*}=U^{*}U^{*\top}=I_{d}$. ($d=k$). This case is easier to
analyze because the dependency between different elements of $x$ or $y$ can be
disentangled. And we are able to provide lower bound for the Hessian. Before
we introduce the lower bound, let’s first define the following quantities for
an activation function $\phi$.

$\displaystyle\alpha_{i,j}:=$
$\displaystyle\operatorname*{{\mathbb{E}}}_{z\sim\mathcal{N}(0,1)}[(\phi(z))^{i}z^{j}],$
(9) $\displaystyle\;\beta_{i,j}:=$
$\displaystyle\operatorname*{{\mathbb{E}}}_{z\sim\mathcal{N}(0,1)}[(\phi^{\prime}(z))^{i}z^{j}],$
$\displaystyle\gamma:=$
$\displaystyle\operatorname*{{\mathbb{E}}}_{z\sim\mathcal{N}(0,1)}[\phi(z)\phi^{\prime}(z)z],$
$\displaystyle\rho:=$
$\displaystyle\min\\{(\alpha_{2,0}\beta_{2,0}-\alpha_{1,0}^{2}\beta_{1,0}^{2}-\beta_{1,0}^{2}\alpha_{1,1}^{2}),\leavevmode\nobreak\
(\alpha_{2,0}\beta_{2,2}-\alpha_{1,0}^{2}\beta_{1,2}^{2}-\gamma^{2})\\}.$

We now present a lower bound for general activation functions including
sigmoid and tanh.

######  Lemma C.1.

Let $(a,b)\in\mathbb{B}$ denote that
$\sum_{i=1}^{k}\|a_{i}\|^{2}+\|b_{i}\|^{2}=1$. Assume $d=k$ and $U^{*},V^{*}$
are unitary matrices, i.e.,
$U^{*\top}U^{*}=U^{*}U^{*\top}=V^{*}V^{*\top}=V^{*\top}V^{*}=I_{d}$, then the
minimal eigenvalue of the population Hessian in Eq. (8) can be simplified as,

$\min_{(a,b)\in\mathbb{B}}\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(x_{i})\phi(y_{i})x^{\top}a_{i}+\phi^{\prime}(y_{i})\phi(x_{i})y^{\top}b_{i}\right)^{2}\right].$

Let $\beta,\rho$ be defined as in Eq. (9). If the activation function $\phi$
satisfies $\beta_{1,1}=0$, then $\lambda_{\min}(H^{*})\geq\rho.$

Since sigmoid and tanh have symmetric derivatives w.r.t. $0$, they satisfy
$\beta_{1,1}=0$. Specifically, we have $\rho\approx 0.000658$ for sigmoid and
$\rho\approx 0.0095$ for tanh. Also for ReLU, $\beta_{1,1}=1/2$, so ReLU does
not fit in this lemma. The full proof of Lemma C.1, the lower bound of the
population Hessian for ReLU and the extension to non-orthogonal cases can be
found in Appendix D.

###  C.3 Error bound for the empirical Hessian near the ground truth

In the previous section, we have shown PD for the population Hessian at the
ground truth for the orthogonal cases. Based on that, we can characterize the
landscape around the ground truth for the empirical risk. In particular, we
bound the difference between the empirical Hessian near the ground truth and
the population Hessian at the ground truth. The theorem below provides the
error bound w.r.t. the number of samples $(n1,n2)$ and the number of
observations $|\Omega|$ for both sigmoid and ReLU activation functions.

######  Theorem C.2.

For any $\epsilon>0$, if

$\displaystyle
n_{1}\gtrsim\epsilon^{-2}td\log^{2}d,n_{2}\gtrsim\epsilon^{-2}td\log^{2}d,|\Omega|\gtrsim\epsilon^{-2}td\log^{2}d,$

then with probability at least $1-d^{-t}$, for sigmoid/tanh,

$\displaystyle\|\nabla^{2}f_{\Omega}(U,V)-\nabla^{2}f_{\mathcal{D}}(U^{*},V^{*})\|\lesssim\epsilon+\|U-U^{*}\|+\|V-V^{*}\|;$

for ReLU,

$\displaystyle\|\nabla^{2}f_{\Omega}(U,V)-\nabla^{2}f_{\mathcal{D}}(U^{*},V^{*})\|\lesssim\left(\|V-V^{*}\|^{1/2}+\|U-U^{*}\|^{1/2}+\epsilon\right)(\|U^{*}\|+\|V^{*}\|)^{2}.$

The key idea to prove this theorem is to use the population Hessian at $(U,V)$
as a bridge.

On one side, we bound the population Hessian at the ground truth and the
population Hessian at $(U,V)$. This would be easy if the second derivative of
the activation function is Lipschitz, which is the case of sigmoid and tanh.
But ReLU doesn’t have this property. However, we can utilize the condition
that the parameters are close enough to the ground truth and the piece-wise
linearity of ReLU to bound this term.

On the other side, we bound the empirical Hessian and the population Hessian.
A natural idea is to apply matrix Bernstein inequality. However, there are two
obstacles. First the Gaussian variables are not uniformly bounded. Therefore,
we instead use Lemma B.7 in [ZSJ+17], which is a loosely-bounded version of
matrix Bernstein inequality. The second obstacle is that each individual
Hessian calculated from one observation $(x,y)$ is not independent from
another observation $(x^{\prime},y^{\prime})$, since they may share the same
feature $x$ or $y$. The analyses for vanilla IMC and MC assume all the
items(users) are given and the observed entries are independently sampled from
the whole matrix. However, our observations are sampled from the joint
distribution of $\mathcal{X}$ and $\mathcal{Y}$.

To handle the dependency, our model assumes the following two-stage sampling
rule. First, the items/users are sampled from their distributions
independently, then given the items and users, the observations $\Omega$ are
sampled uniformly with replacement. The key question here is how to combine
the error bounds from these two stages. Fortunately, we found special
structures in the blocks of Hessian which enables us to separate $x,y$ for
each block, and bound the errors in stage separately. See Appendix E for
details.

##  Appendix D Positive Definiteness of Population Hessian

###  D.1 Orthogonal case

We first study the orthogonal case, where $d=k$ and $U^{*},V^{*}$ are unitary
matrices, i.e.,
$U^{*\top}U^{*}=U^{*}U^{*\top}=V^{*}V^{*\top}=V^{*\top}V^{*}=I_{d}$.

####  D.1.1 Lower bound on minimum eigenvalue

######  Lemma D.1 (Restatement of Lemma C.1).

Let $(a,b)\in\mathbb{B}$ denote that
$\sum_{i=1}^{k}\|a_{i}\|^{2}+\|b_{i}\|^{2}=1$. Assume $d=k$ and $U^{*},V^{*}$
are unitary matrices, i.e.,
$U^{*\top}U^{*}=U^{*}U^{*\top}=V^{*}V^{*\top}=V^{*\top}V^{*}=I_{d}$, then the
minimal eigenvalue of the population Hessian in Eq. (8) can be simplified as,

$\lambda_{\min}(H^{*})=\min_{(a,b)\in\mathbb{B}}\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(x_{i})\phi(y_{i})x^{\top}a_{i}+\phi^{\prime}(y_{i})\phi(x_{i})y^{\top}b_{i}\right)^{2}\right].$
(10)

Let $\beta,\rho$ be defined as in Eq. (9). If the activation function $\phi$
satisfies $\beta_{1,1}=0$, then $\lambda_{\min}(H^{*})\geq\rho.$

###### Proof.

In the orthogonal case, we can easily transform Eq. (8) to Eq. (10) since
$x,y$ are normal distribution. Now we can decompose Eq. (10) into the
following three terms.

$\displaystyle\leavevmode\nobreak\
\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(x_{i})\phi(y_{i})x^{\top}a_{i}+\phi^{\prime}(y_{i})\phi(x_{i})y^{\top}b_{i}\right)^{2}\right]$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\underbrace{\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(x_{i})\phi(y_{i})x^{\top}a_{i}\right)^{2}\right]}_{C}+\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(y_{i})\phi(x_{i})y^{\top}b_{i}\right)^{2}\right]$
$\displaystyle\leavevmode\nobreak\
+\underbrace{2\operatorname*{{\mathbb{E}}}_{x,y}\left[\sum_{i,j}\phi^{\prime}(x_{i})\phi(y_{i})x^{\top}a_{i}\phi^{\prime}(y_{j})\phi(x_{j})y^{\top}b_{j}\right]}_{D}.$

Note that the first term is similar to the second term, so we just lower bound
the first term and the third term. Define
$A=[a_{1},a_{2},\cdots,a_{k}],B=[b_{1},b_{2},\cdots,b_{k}]$. Let $A_{o}$ be
the off-diagonal part of $A$ and $A_{d}$ be the diagonal part of $A$, i.e.,
$A_{o}+A_{d}=A$. And let $g_{A}=\operatorname{diag}(A)$ be the vector of the
diagonal elements of A. We will bound $C$ and $D$ in the following.

For $C$, we have

$\displaystyle\leavevmode\nobreak\
\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(x_{i})\phi(y_{i})x^{\top}a_{i}\right)^{2}\right]$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\sum_{i=1}^{k}\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\phi^{\prime}(x_{i})\phi(y_{i})x^{\top}a_{i}\right)^{2}\right]+\sum_{i\neq
j}\operatorname*{{\mathbb{E}}}_{x,y}\left[\phi^{\prime}(x_{i})\phi(y_{i})x^{\top}a_{i}\cdot\phi^{\prime}(x_{j})\phi(y_{j})x^{\top}a_{j}\right]$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\sum_{i=1}^{k}\alpha_{2,0}\left[(a_{i}^{\top}e_{i})^{2}(\beta_{2,2}-\beta_{2,0})+\beta_{2,0}\|a_{i}\|^{2}\right]$
$\displaystyle\leavevmode\nobreak\ +\sum_{i\neq
j}\alpha_{1,0}^{2}\left[\beta_{1,0}^{2}a_{i}^{\top}a_{j}+(\beta_{1,2}\beta_{1,0}-\beta_{1,0}^{2})(a_{i}^{\top}e_{i}e_{i}^{\top}a_{j}+a_{i}^{\top}e_{j}a_{j}^{\top}e_{j})+\beta_{1,1}^{2}(a_{i}^{\top}e_{i}a_{j}^{\top}e_{j}+a_{i}^{\top}e_{j}a_{j}^{\top}e_{i})\right]$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\ C_{1}+C_{2}.$

where the last step follows by

$\displaystyle C_{1}=$ $\displaystyle\leavevmode\nobreak\
\sum_{i=1}^{k}\alpha_{2,0}\left[(a_{i}^{\top}e_{i})^{2}(\beta_{2,2}-\beta_{2,0})+\beta_{2,0}\|a_{i}\|^{2}\right]$
$\displaystyle C_{2}=$ $\displaystyle\leavevmode\nobreak\ \sum_{i\neq
j}\alpha_{1,0}^{2}\left[\beta_{1,0}^{2}a_{i}^{\top}a_{j}+(\beta_{1,2}\beta_{1,0}-\beta_{1,0}^{2})(a_{i}^{\top}e_{i}e_{i}^{\top}a_{j}+a_{i}^{\top}e_{j}a_{j}^{\top}e_{j})+\beta_{1,1}^{2}(a_{i}^{\top}e_{i}a_{j}^{\top}e_{j}+a_{i}^{\top}e_{j}a_{j}^{\top}e_{i})\right]$

First we can simplify $C_{1}$ in the following sense,

$\displaystyle C_{1}=$ $\displaystyle\leavevmode\nobreak\
\alpha_{2,0}(\beta_{2,2}-\beta_{2,0})\sum_{i=1}^{k}(a_{i}^{\top}e_{i})^{2}+\alpha_{2,0}\beta_{2,0}\sum_{i=1}^{k}\|a_{i}\|_{2}^{2}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\alpha_{2,0}(\beta_{2,2}-\beta_{2,0})\|\operatorname{diag}(A)\|_{2}^{2}+\alpha_{2,0}\beta_{2,0}\|A\|_{F}^{2},$

where the last step follows by Fact B.1.

We can rewrite $C_{2}$ in the following sense

$\displaystyle
C_{2}=\alpha_{1,0}^{2}(\beta_{1,0}^{2}C_{2,1}+(\beta_{1,2}\beta_{1,0}-\beta_{1,0}^{2})\cdot(C_{2,2}+C_{2,3})+\beta_{1,1}^{2}(C_{2,4}+C_{2,5})).$

where

$\displaystyle C_{2,1}=$ $\displaystyle\leavevmode\nobreak\ \sum_{i\neq
j}a_{i}^{\top}a_{j}$ $\displaystyle C_{2,2}=$
$\displaystyle\leavevmode\nobreak\ \sum_{i\neq
j}a_{i}^{\top}e_{i}e_{i}^{\top}a_{j}$ $\displaystyle C_{2,3}=$
$\displaystyle\leavevmode\nobreak\ \sum_{i\neq
j}a_{i}^{\top}e_{j}e_{j}^{\top}a_{j}$ $\displaystyle C_{2,4}=$
$\displaystyle\leavevmode\nobreak\ \sum_{i\neq
j}a_{i}^{\top}e_{i}a_{j}^{\top}e_{j}$ $\displaystyle C_{2,5}=$
$\displaystyle\leavevmode\nobreak\ \sum_{i\neq
j}a_{i}^{\top}e_{j}a_{j}^{\top}e_{i}$

Using Fact B.1, we have

$\displaystyle C_{2,1}=\|A\cdot{\bf 1}\|_{2}^{2}-\|A\|_{F}^{2}.$

Using Fact B.2, we have

$\displaystyle C_{2,2}=$ $\displaystyle\leavevmode\nobreak\
(\operatorname{diag}(A)^{\top}\cdot(A\cdot{\bf
1}))-\|\operatorname{diag}(A)\|_{2}^{2},$ $\displaystyle C_{2,3}=$
$\displaystyle\leavevmode\nobreak\
(\operatorname{diag}(A)^{\top}\cdot(A\cdot{\bf
1}))-\|\operatorname{diag}(A)\|_{2}^{2},$ $\displaystyle C_{2,4}=$
$\displaystyle\leavevmode\nobreak\ (\operatorname{diag}(A)^{\top}\cdot{\bf
1})^{2}-\|\operatorname{diag}(A)\|_{2}^{2},$ $\displaystyle C_{2,5}=$
$\displaystyle\leavevmode\nobreak\ \langle
A^{\top},A\rangle-\|\operatorname{diag}(A)\|_{2}^{2}.$

Thus,

$\displaystyle C_{2}=$ $\displaystyle\leavevmode\nobreak\
\alpha_{1,0}^{2}(\beta_{1,0}^{2}(\|A\cdot{\bf 1}\|_{2}^{2}-\|A\|_{F}^{2})$
$\displaystyle\leavevmode\nobreak\
+(\beta_{1,2}\beta_{1,0}-\beta_{1,0}^{2})2\cdot(\operatorname{diag}(A)^{\top}\cdot(A\cdot{\bf
1})-\|\operatorname{diag}(A)\|_{2}^{2})$ $\displaystyle\leavevmode\nobreak\
+\beta_{1,1}^{2}((\operatorname{diag}(A)^{\top}\cdot{\bf 1})^{2}+\langle
A^{\top},A\rangle-2\|\operatorname{diag}(A)\|_{2}^{2})).$

We consider $C_{1}+C_{2}$ by focusing different terms, for the
$\|A\|_{F}^{2}$(from $C_{1}$ and $C_{2}$), we have

$\displaystyle(\alpha_{2,0}\beta_{2,0}-\alpha_{1,0}^{2}\beta_{1,0}^{2})\|A\|_{F}^{2}.$

For the term $\langle A,A^{\top}\rangle$ (from $C_{2,5}$), we have

$\displaystyle\alpha_{1,0}^{2}\beta_{1,1}^{2}\langle A,A^{\top}\rangle.$

For the term $\|\operatorname{diag}(A)\|_{2}^{2}$ (from $C_{1}$ and $C_{2}$),
we have

$\displaystyle(\alpha_{2,0}(\beta_{2,2}-\beta_{2,0})-2\alpha_{1,0}^{2}(\beta_{1,2}\beta_{1,0}-\beta_{1,0}^{2})-2\alpha_{1,0}\beta_{1,1}^{2})\|\operatorname{diag}(A)\|_{2}^{2}$

For the term $\|A\cdot{\bf 1}\|_{2}^{2}$ (from $C_{2,1}$), we have

$\displaystyle\alpha_{1,0}^{2}\beta_{1,0}^{2}\|A\cdot{\bf 1}\|_{2}^{2}.$

For the term $\operatorname{diag}(A)^{\top}\cdot A\cdot{\bf 1}$ (from
$C_{2,2}$ and $C_{2,3}$), we have

$\displaystyle
2\alpha_{1,0}^{2}(\beta_{1,2}\beta_{1,0}-\beta_{1,0}^{2})\operatorname{diag}(A)^{\top}\cdot
A\cdot{\bf 1}.$

For the term $(\operatorname{diag}(A)^{\top}\cdot{\bf 1})^{2}$ (from
$C_{2,4}$), we have

$\displaystyle\alpha_{1,0}^{2}\beta_{1,1}^{2}(\operatorname{diag}(A)^{\top}\cdot{\bf
1})^{2}.$

Putting it all together, we have

$\displaystyle C_{1}+C_{2}=$ $\displaystyle\leavevmode\nobreak\
(\alpha_{2,0}\beta_{2,0}-\alpha_{1,0}^{2}\beta_{1,0}^{2})\|A\|_{F}^{2}+\alpha_{1,0}^{2}\beta_{1,1}^{2}\langle
A,A^{\top}\rangle$ $\displaystyle\leavevmode\nobreak\
+(\alpha_{2,0}(\beta_{2,2}-\beta_{2,0})-2\alpha_{1,0}^{2}(\beta_{1,2}\beta_{1,0}-\beta_{1,0}^{2})-2\alpha_{1,0}^{2}\beta_{1,1}^{2})\cdot\|\operatorname{diag}(A)\|^{2}$
$\displaystyle\leavevmode\nobreak\
+\alpha_{1,0}^{2}\beta_{1,0}^{2}\|A\cdot{\bf
1}\|^{2}+2\alpha_{1,0}^{2}(\beta_{1,2}\beta_{1,0}-\beta_{1,0}^{2})(\operatorname{diag}(A)^{\top}\cdot
A\cdot{\bf
1})+\alpha_{1,0}^{2}\beta_{1,1}^{2}(\operatorname{diag}(A)^{\top}\cdot{\bf
1})^{2}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\
(\alpha_{2,0}\beta_{2,0}-\alpha_{1,0}^{2}\beta_{1,0}^{2})(\|A_{o}\|_{F}^{2}+\|g_{A}\|^{2})+\alpha_{1,0}^{2}\beta_{1,1}^{2}(\langle
A_{o},A_{o}^{\top}\rangle+\|g_{A}\|^{2})$ $\displaystyle\leavevmode\nobreak\
+(\alpha_{2,0}\beta_{2,2}-\alpha_{2,0}\beta_{2,0}-2\alpha_{1,0}^{2}\beta_{1,2}\beta_{1,0}+2\alpha_{1,0}^{2}\beta_{1,0}^{2}-2\alpha_{1,0}^{2}\beta_{1,1}^{2})\cdot\|g_{A}\|^{2}$
$\displaystyle\leavevmode\nobreak\
+\alpha_{1,0}^{2}\beta_{1,0}^{2}(\|g_{A}\|^{2}+\|A_{o}\cdot{\bf
1}\|^{2}+2g_{A}^{\top}\cdot A_{o}\cdot{\bf 1})$
$\displaystyle\leavevmode\nobreak\
+2\alpha_{1,0}^{2}(\beta_{1,2}\beta_{1,0}-\beta_{1,0}^{2})(g_{A}^{\top}\cdot
A_{o}\cdot{\bf
1}+\|g_{A}\|^{2})+\alpha_{1,0}^{2}\beta_{1,1}^{2}(g_{A}^{\top}\cdot{\bf
1})^{2}$ $\displaystyle=$ $\displaystyle\leavevmode\nobreak\
(\alpha_{2,0}\beta_{2,0}-\alpha_{1,0}^{2}\beta_{1,0}^{2})\|A_{o}\|_{F}^{2}+\alpha_{1,0}^{2}\beta_{1,1}^{2}\langle
A_{o},A_{o}^{\top}\rangle+(\alpha_{2,0}\beta_{2,2}-\alpha_{1,0}^{2}\beta_{1,1}^{2})\cdot\|g_{A}\|^{2}$
$\displaystyle\leavevmode\nobreak\
+\alpha_{1,0}^{2}\beta_{1,0}^{2}(\|A_{o}\cdot{\bf
1}\|^{2})+2\alpha_{1,0}^{2}\beta_{1,2}\beta_{1,0}(g_{A}^{\top}\cdot
A_{o}\cdot{\bf 1})+\alpha_{1,0}^{2}\beta_{1,1}^{2}(g_{A}^{\top}\cdot{\bf
1})^{2}.$

By doing a series of equivalent transformations, we have removed the
expectation and the formula $C$ becomes a form of $A$ and the moments of
$\phi$. These equivalent transforms are mainly based on the fact that
$x_{i},x_{j},y_{i},y_{j}$ for any $i\neq j$ are independent on each other.

Similarly we can reformulate $D$,

$\displaystyle\leavevmode\nobreak\
\operatorname*{{\mathbb{E}}}_{x,y}\left[\sum_{i,j}\phi^{\prime}(x_{i})\phi(y_{i})x^{\top}a_{i}\phi^{\prime}(y_{j})\phi(x_{j})y^{\top}b_{j}\right]$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\sum_{i}\operatorname*{{\mathbb{E}}}_{x,y}\left[\phi^{\prime}(x_{i})\phi(y_{i})x^{\top}a_{i}\phi^{\prime}(y_{j}i)\phi(x_{i})y^{\top}b_{i}\right]+\sum_{i\neq
j}\operatorname*{{\mathbb{E}}}_{x,y}\left[\phi^{\prime}(x_{i})\phi(y_{i})x^{\top}a_{i}\phi^{\prime}(y_{j})\phi(x_{j})y^{\top}b_{j}\right]$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\sum_{i}\gamma^{2}a_{i}^{\top}e_{i}b_{i}^{\top}e_{i}+\sum_{i\neq
j}\alpha_{1,1}^{2}a_{i}^{\top}e_{j}b_{j}^{\top}e_{i}+\alpha_{1,1}\beta_{1,1}(a_{i}^{\top}e_{j}b_{j}^{\top}e_{j}+a_{i}^{\top}e_{i}b_{j}^{\top}e_{i})+\beta_{1,1}^{2}a_{i}^{\top}e_{i}b_{j}^{\top}e_{j}$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
(\gamma^{2}-\beta_{1,0}^{2}\alpha_{1,1}^{2}-2\alpha_{1,0}\alpha_{1,1}\beta_{1,0}\beta_{1,1}-\alpha_{1,0}^{2}\beta_{1,1}^{2})g_{A}^{\top}g_{B}$
$\displaystyle\leavevmode\nobreak\ +\beta_{1,0}^{2}\alpha_{1,1}^{2}\langle
A,B^{\top}\rangle+\alpha_{1,0}^{2}\beta_{1,1}^{2}(g_{A}^{\top}1)(g_{B}^{\top}1)$
$\displaystyle\leavevmode\nobreak\
+\alpha_{1,0}\alpha_{1,1}\beta_{1,0}\beta_{1,1}[(A1)^{\top}g_{B}+(B1)^{\top}g_{A}]$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
(\gamma^{2}-\alpha_{1,0}^{2}\beta_{1,1}^{2})g_{A}^{\top}g_{B}+\beta_{1,0}^{2}\alpha_{1,1}^{2}\langle
A_{o},B_{o}^{\top}\rangle+\alpha_{1,0}^{2}\beta_{1,1}^{2}(g_{A}^{\top}1)(g_{B}^{\top}1)$
$\displaystyle\leavevmode\nobreak\
+\alpha_{1,0}\alpha_{1,1}\beta_{1,0}\beta_{1,1}[(A_{o}1)^{\top}g_{B}+(B_{o}1)^{\top}g_{A}].$

Combining the above results, we have

$\displaystyle\lambda_{\min}(H^{*})=$ $\displaystyle\leavevmode\nobreak\
\min_{\|A\|_{F}^{2}+\|B\|_{F}^{2}=1}\bigg{(}\beta_{1,0}^{2}\alpha_{1,1}^{2}\|A_{o}+B_{o}^{\top}\|_{F}^{2}$
(11) $\displaystyle\leavevmode\nobreak\
+\|\alpha_{1,0}\beta_{1,0}A_{o}1+\alpha_{1,0}\beta_{1,2}g_{A}+\alpha_{1,1}\beta_{1,1}g_{B}\|^{2}$
$\displaystyle\leavevmode\nobreak\
+\|\alpha_{1,0}\beta_{1,0}B_{o}1+\alpha_{1,0}\beta_{1,2}g_{B}+\alpha_{1,1}\beta_{1,1}g_{A}\|^{2}$
$\displaystyle\leavevmode\nobreak\
+(\alpha_{2,0}\beta_{2,0}-\alpha_{1,0}^{2}\beta_{1,0}^{2}-\beta_{1,0}^{2}\alpha_{1,1}^{2}-\alpha_{1,0}^{2}\beta_{1,1}^{2})(\|A_{o}\|_{F}^{2}+\|B_{o}\|_{F}^{2})$
$\displaystyle\leavevmode\nobreak\
+1/2\cdot\alpha_{1,0}^{2}\beta_{1,1}^{2}(\|A_{o}+A_{o}^{\top}\|_{F}^{2}+\|B_{o}+B_{o}^{\top}\|_{F}^{2})$
$\displaystyle\leavevmode\nobreak\
+[\alpha_{2,0}\beta_{2,2}-\alpha_{1,0}^{2}\beta_{1,1}^{2}-\alpha_{1,0}^{2}\beta_{1,2}^{2}-\alpha_{1,1}^{2}\beta_{1,1}^{2}]\cdot(\|g_{A}\|^{2}+\|g_{B}\|^{2})$
$\displaystyle\leavevmode\nobreak\
+2(\gamma^{2}-\alpha_{1,0}^{2}\beta_{1,1}^{2}-2\alpha_{1,0}\alpha_{1,1}\beta_{1,1}\beta_{1,2})g_{A}^{\top}g_{B}$
$\displaystyle\leavevmode\nobreak\
+\alpha_{1,0}^{2}\beta_{1,1}^{2}(g_{A}^{\top}1+g_{B}^{\top}1)^{2}\bigg{)}.$

The final output of the above formula has a clear form: most non-negative
terms are extracted. $A,B$ are separated into the off-diagonal elements and
off-diagonal elements and these two terms can be dealt with independently. Now
we consider the activation functions that satisfy $\beta_{1,1}=0$, which
further simplifies the equation. Note that Sigmoid and $\tanh$ satisfy this
condition.

Finally, for $\beta_{1,1}=0$, we obtain

$\displaystyle\lambda_{\min}(H^{*})=$ $\displaystyle\leavevmode\nobreak\
\min_{\sum_{i=1}^{k}\|a_{i}\|^{2}+\|b_{i}\|^{2}=1}\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(x_{i})\phi(y_{i})x^{\top}a_{i}+\phi^{\prime}(y_{i})\phi(x_{i})y^{\top}b_{i}\right)^{2}\right]$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\min_{\|A\|_{F}^{2}+\|B\|_{F}^{2}=1}(\alpha_{2,0}\beta_{2,0}-\alpha_{1,0}^{2}\beta_{1,0}^{2}-\beta_{1,0}^{2}\alpha_{1,1}^{2})(\|A_{o}\|_{F}^{2}+\|B_{o}\|_{F}^{2})$
$\displaystyle\leavevmode\nobreak\
+(\alpha_{2,0}\beta_{2,2}-\alpha_{1,0}^{2}\beta_{1,2}^{2}-\gamma^{2})(\|g_{A}\|^{2}+\|g_{B}\|^{2})$
$\displaystyle\leavevmode\nobreak\
+\beta_{1,0}^{2}\alpha_{1,1}^{2}\|A_{o}+B_{o}^{\top}\|_{F}^{2}+\gamma^{2}\|g_{A}+g_{B}\|^{2}$
$\displaystyle\leavevmode\nobreak\
+\alpha_{1,0}^{2}(\|\beta_{1,0}g_{A}+\beta_{1,2}A_{o}1\|^{2}+\alpha_{1,0}^{2}\|\beta_{1,0}g_{A}+\beta_{1,2}B_{o}1\|^{2})$
$\displaystyle\geq$ $\displaystyle\leavevmode\nobreak\
\underbrace{\min\\{(\alpha_{2,0}\beta_{2,0}-\alpha_{1,0}^{2}\beta_{1,0}^{2}-\beta_{1,0}^{2}\alpha_{1,1}^{2}),(\alpha_{2,0}\beta_{2,2}-\alpha_{1,0}^{2}\beta_{1,2}^{2}-\gamma^{2})\\}}_{:=\rho}.$

For sigmoid, we have $\rho=0.000658$; for tanh, we have $\rho=0.0095$.

∎

The following lemma will be used when transforming non-orthogonal cases to
orthogonal cases.

######  Lemma D.2.

For any $A=[a_{1},a_{2},\cdots,a_{k}]\in\operatorname{{\mathbb{R}}}^{d\times
k}$, we have,

$\displaystyle\operatorname*{{\mathbb{E}}}_{x,y\sim\mathcal{D}_{k}}\left[\left\|\sum_{i=1}^{k}\phi^{\prime}(x_{i})\phi(y_{i})a_{i}\right\|^{2}\right]\geq\leavevmode\nobreak\
(\alpha_{2,0}\beta_{2,0}-\alpha_{1,0}^{2}\beta_{1,0}^{2})\|A\|_{F}^{2}.$

###### Proof.

Recall ${\bf 1}\in\operatorname{{\mathbb{R}}}^{d}$ denote the all ones vector.

$\displaystyle\leavevmode\nobreak\
\operatorname*{{\mathbb{E}}}_{x,y\sim\mathcal{D}_{k}}\left[\left\|\sum_{i=1}^{k}\phi^{\prime}(x_{i})\phi(y_{i})a_{i}\right\|^{2}\right]$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\operatorname*{{\mathbb{E}}}_{x,y\sim\mathcal{D}_{k}}\left[\sum_{i=1}^{k}(\phi^{\prime}(x_{i})\phi(y_{i}))^{2}\|a_{i}\|^{2}\right]+\operatorname*{{\mathbb{E}}}_{x,y\sim\mathcal{D}_{k}}\left[\sum_{i\neq
j}\phi^{\prime}(x_{i})\phi(y_{i})\phi^{\prime}(x_{j})\phi(y_{j})a_{i}^{\top}a_{j}\right]$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
(\alpha_{2,0}\beta_{2,0}-\alpha_{1,0}^{2}\beta_{1,0}^{2})\|A\|_{F}^{2}+\alpha_{1,0}^{2}\beta_{1,0}^{2}\|A\cdot
1\|^{2}$ $\displaystyle\geq$ $\displaystyle\leavevmode\nobreak\
(\alpha_{2,0}\beta_{2,0}-\alpha_{1,0}^{2}\beta_{1,0}^{2})\|A\|_{F}^{2}.$

Thus, we complete the proof. ∎

Now let’s show the PD of the population Hessian of Eq. (4) for the ReLU case.
where $u^{*(1)}$ is the first row of $U^{*}$ and $W\in\mathbb{R}^{(d-1)\times
k}$.

######  Lemma D.3.

Consider the activation function to be ReLU. Assume $k=d$, $U^{*},V^{*}$ are
unitary matrices and $u^{*}_{1,i}\neq 0,\forall i\in[k]$. Then the minimal
eigenvalue of the corresponding population Hessian of Eq. (4) is lower
bounded,

$\lambda_{\min}(\nabla^{2}f_{\mathcal{D}}^{\mathrm{ReLU}}(W^{*},V^{*}))\gtrsim\min_{i\in[k]}\\{u_{1,i}^{*2}\\},$

where $W^{*}=U^{*}_{2:d,:}$ is the last $d-1$ rows of $U^{*}$ and

$\displaystyle
f_{\mathcal{D}}^{\mathrm{ReLU}}(W,V):=\operatorname*{{\mathbb{E}}}_{x,y}\left[(\phi(W^{\top}x_{2:d}+x_{1}(u^{*(1)})^{\top})^{\top}\phi(V^{\top}y)-A(x,y))^{2}\right],$
(12)

###### Proof.

By fixing $u_{i,1}=u^{*}_{i,1},\forall i\in[k]$, we can rewrite the minimal
eigenvalue of the Hessian as follows. For simplicity, we denote
$\lambda_{\min}(H):=\lambda_{\min}(\nabla^{2}f_{\mathcal{D}}^{\mathrm{ReLU}}(W^{*},V^{*}))$.
First we observe that

$\lambda_{\min}(H)=\min_{\begin{subarray}{c}\sum_{i=1}^{k}\|a_{i}\|^{2}+\|b_{i}\|^{2}=1\\\
a_{i,1}=0,\forall
i\in[k]\end{subarray}}\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(u_{i}^{*\top}x)\phi(v_{i}^{*\top}y)x^{\top}a_{i}+\phi^{\prime}(v_{i}^{*\top}y)\phi(u_{i}^{*\top}x)y^{\top}b_{i}\right)^{2}\right].$
(13)

Without loss of generality, we assume $V^{*}=I$. Set $x=U^{*}s$, then we have

$\displaystyle\lambda_{\min}(H)$
$\displaystyle=\min_{\begin{subarray}{c}\sum_{i=1}^{k}\|a_{i}\|^{2}+\|b_{i}\|^{2}=1\\\
a_{i,1}=0,\forall
i\in[k]\end{subarray}}\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(s_{i})\phi(y_{i})s^{\top}U^{*\top}a_{i}+\phi^{\prime}(y_{i})\phi(x_{i})y^{\top}b_{i}\right)^{2}\right]$
$\displaystyle=\min_{\begin{subarray}{c}\sum_{i=1}^{k}\|a_{i}\|^{2}+\|b_{i}\|^{2}=1\\\
u^{*(1)}a_{i}=0,\forall
i\in[k]\end{subarray}}\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(s_{i})\phi(y_{i})s^{\top}a_{i}+\phi^{\prime}(y_{i})\phi(x_{i})y^{\top}b_{i}\right)^{2}\right],$

where $u^{*(1)}$ is the first row of $U^{*}$ and the second equality is
because we replace $U^{*\top}a_{i}$ by $a_{i}$. In the ReLU case, we have

$\alpha_{1,0}=\alpha_{1,1}=\alpha_{2,0}=\beta_{1,0}=\beta_{1,1}=\beta_{1,1}=\beta_{2,0}=\beta_{2,2}=\gamma=1/2.$

According to Eq. (11), we have

$\displaystyle\lambda_{\min}(H)\geq\min_{\|A\|_{F}^{2}+\|B\|_{F}^{2}=1,u^{*(1)}A=0}C_{0}($
$\displaystyle\|A_{o}\|_{F}^{2}+\|B_{o}\|_{F}^{2}+\|A_{o}+A_{o}^{\top}\|_{F}^{2}/2+\|B_{o}+B_{o}^{\top}\|_{F}^{2}/2$
$\displaystyle+\|A_{o}+B_{o}^{\top}\|_{F}^{2}+\|g_{A}+g_{B}\|^{2}$
$\displaystyle+\|A_{o}1+g_{A}+g_{B}\|^{2}+\|B_{o}1+g_{A}+g_{B}\|^{2}+(g_{A}^{\top}1+g_{B}^{\top}1)^{2}),$

where $C_{0}$ is a universal constant. Now we show that there exists a
positive number $c_{0}$ such that $\lambda_{\min}(H)\geq c_{0}$. If there is
no such number, i.e., $\lambda_{\min}(H)=0$, then we have $A_{o}=B_{o}=0$,
$g_{A}=-g_{B}$. By the assumption that $u_{1,i}^{*}\neq 0$ and the condition
$u^{*(1)}A=0$, we have $g_{A}=g_{B}=0$, which violates
$\|A\|_{F}^{2}+\|B\|_{F}^{2}=1$. So $\lambda_{\min}(H)>0$. An exact value for
$c_{0}$ is postponed to Theorem D.6, which gives the lower bound for the non-
orthogonal case. ∎

###  D.2 Non-orthogonal Case

The restriction of orthogonality on $U,V$ is too strong. We need to consider
general non-orthogonal cases. With Gaussian assumption, the non-orthogonal
case can be transformed to the orthogonal case according to the following
relationship.

######  Lemma D.4.

Let $U\in\mathbb{R}^{d\times k}$ be a full-column rank matrix. Let
$g:\mathbb{R}^{k}\rightarrow[0,\infty)$. Define
$\lambda(U)=\sigma_{1}^{k}(U)/(\prod_{i=1}^{k}\sigma_{i}(U))$. Let ${\cal D}$
denote the normal distribution. Then

$\operatorname*{{\mathbb{E}}}_{x\sim\mathcal{D}_{d}}\left[g(U^{\top}x)\right]\geq\frac{1}{\lambda(U)}\operatorname*{{\mathbb{E}}}_{z\sim\mathcal{D}_{k}}\left[g(\sigma_{k}(U)z)\right].$
(14)

Remark This lemma transforms $U^{\top}x$, where the elements of $x$ are mixed,
to $\sigma_{k}(U)z$, where all the elements are independently fed into $g$
with the sacrifices of a condition number of $U$. Using Lemma D.4, we are able
to show the PD for non-orthogonal $U^{*},V^{*}$.

###### Proof.

Let $P\in\mathbb{R}^{d\times k}$ be the orthonormal basis of $U$, and let
$W=[w_{1},w_{2},\cdots,w_{k}]=P^{\top}U\in\mathbb{R}^{k\times k}$.

$\displaystyle\leavevmode\nobreak\
\operatorname*{{\mathbb{E}}}_{x\sim\mathcal{D}_{d}}[g(U^{\top}x)]$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\operatorname*{{\mathbb{E}}}_{z\sim\mathcal{D}_{k}}[g(W^{\top}z)]$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\int(2\pi)^{-k/2}g(W^{\top}z)e^{-\|z\|^{2}/2}\mathrm{d}z$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\
\int(2\pi)^{-k/2}g(s)e^{-\|W^{\dagger\top}s\|^{2}/2}|\det(W^{\dagger})|\mathrm{d}s$
$\displaystyle\geq$ $\displaystyle\leavevmode\nobreak\
\int(2\pi)^{-k/2}g(s)e^{-\sigma_{1}^{2}(W^{\dagger})\|s\|^{2}/2}|\det(W^{\dagger})|\mathrm{d}s$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\int(2\pi)^{-k/2}g\left(\frac{1}{\sigma_{1}(W^{\dagger})}t\right)e^{-\|t\|^{2}/2}|\det(W^{\dagger})|/\sigma_{1}^{k}(W^{\dagger})\mathrm{d}t$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\frac{1}{\lambda(W)}\int(2\pi)^{-k/2}g(\sigma_{k}(W)t)e^{-\|t\|^{2}/2}\mathrm{d}t$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\frac{1}{\lambda(U)}\operatorname*{{\mathbb{E}}}_{z\sim\mathcal{D}_{k}}[g(\sigma_{k}(U)z)],$

where the third step follows by replacing $z$ by $z=W^{\dagger\top}s$, the
fourth step follows by the fact that
$\|W^{\dagger\top}s\|\leq\sigma_{1}(W^{\dagger})\|s\|$, and the fifth step
follows replacing $s$ by $s=\frac{1}{\sigma_{1}(W^{\dagger})}t$. ∎

Using Lemma D.4, we are able to provide the lower bound for the minimal
eigenvalue for sigmoid and tanh.

######  Theorem D.5.

Assume $\sigma_{k}(U^{*})=\sigma_{k}(V^{*})=1$. Assume $\beta_{1,1}$ defined
in Eq. (9) is $0$. Then the minimal eigenvalue of Hessian defined in Eq. (8)
can be lower bounded by,

$\displaystyle\lambda_{\min}(H^{*})\geq\frac{\rho}{\lambda(U^{*})\lambda(V^{*})\max\\{\kappa(U^{*}),\kappa(V^{*})\\}}$

where

$\displaystyle\lambda(U)=\sigma_{1}^{k}(U)/(\Pi_{i=1}^{k}\sigma_{i}(U)),\kappa(U)=\sigma_{1}(U)/\sigma_{k}(U).$

###### Proof.

Let $P\in\mathbb{R}^{d\times k},Q\in\mathbb{R}^{d\times k}$ be the orthonormal
basis of $U^{*},V^{*}$ respectively. Let
$R\in\operatorname{{\mathbb{R}}}^{k\times
k},S\in\operatorname{{\mathbb{R}}}^{k\times k}$ satisfy that $U^{*}=P\cdot R$
and $V^{*}=Q\cdot S$. Let
$P_{\perp}\in\mathbb{R}^{d\times(d-k)},Q_{\perp}\in\mathbb{R}^{d\times(d-k)}$
be the orthogonal complement of $P,Q$ respectively. Set $a_{i}=P\cdot
s_{i}+P_{\perp}\cdot t_{i}$ and $b_{i}=Q\cdot p_{i}+Q_{\perp}\cdot q_{i}$.
Then we can decompose the minimal eigenvalue problem into three terms.

$\displaystyle\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(u_{i}^{*\top}x)\phi(v_{i}^{*\top}y)x^{\top}a_{i}+\phi^{\prime}(v_{i}^{*\top}y)\phi(u_{i}^{*\top}x)y^{\top}b_{i}\right)^{2}\right]$
$\displaystyle=$
$\displaystyle\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(u_{i}^{*\top}x)\phi(v_{i}^{*\top}y)x^{\top}(Ps_{i}+P_{\perp}t_{i})+\phi^{\prime}(v_{i}^{*\top}y)\phi(u_{i}^{*\top}x)y^{\top}(Qp_{i}+Q_{\perp}q_{i})\right)^{2}\right]$
$\displaystyle=$
$\displaystyle\underbrace{\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(u_{i}^{*\top}x)\phi(v_{i}^{*\top}y)x^{\top}Ps_{i}+\phi^{\prime}(v_{i}^{*\top}y)\phi(u_{i}^{*\top}x)y^{\top}Qp_{i}\right)^{2}\right]}_{C_{1}}$
$\displaystyle+\underbrace{\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(u_{i}^{*\top}x)\phi(v_{i}^{*\top}y)x^{\top}P_{\perp}t_{i}\right)^{2}\right]}_{C_{2}}+\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(v_{i}^{*\top}y)\phi(u_{i}^{*\top}x)y^{\top}Q_{\perp}q_{i}\right)^{2}\right],$

where we omit the terms containing a single independent Gaussian variable,
whose expectation is zero. Using Lemma D.4, we can lower bound the term
$C_{1}$ as follows,

$\displaystyle C_{1}=$ $\displaystyle\leavevmode\nobreak\
\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(u_{i}^{*\top}x)\phi(v_{i}^{*\top}y)x^{\top}U^{*}R^{-1}s_{i}+\phi^{\prime}(v_{i}^{*\top}y)\phi(u_{i}^{*\top}x)y^{\top}V^{*}S^{-1}p_{i}\right)^{2}\right]$
$\displaystyle\geq$ $\displaystyle\leavevmode\nobreak\
\frac{1}{\lambda(U^{*})\lambda(V^{*})}\cdot\operatorname*{{\mathbb{E}}}_{x,y\sim\mathcal{D}_{k}}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(\sigma_{k}(U^{*})x_{i}))\phi(y_{i})x^{\top}R^{-1}s_{i}\sigma_{k}(U^{*})\right.\right.$
$\displaystyle\leavevmode\nobreak\
\left.\left.+\phi^{\prime}(\sigma_{k}(V^{*})y_{i})\phi(\sigma_{k}(U^{*})x_{i})y^{\top}S^{-1}p_{i}\sigma_{k}(V^{*})\right)^{2}\right].$

And

$\displaystyle C_{2}\geq$
$\displaystyle\operatorname*{{\mathbb{E}}}_{x,y}\left[\left\|\sum_{i=1}^{k}\phi^{\prime}(u_{i}^{*\top}x)\phi(v_{i}^{*\top}y)t_{i}\right\|^{2}\right]$
$\displaystyle\geq$
$\displaystyle\frac{1}{\lambda(U^{*})\lambda(V^{*})}\operatorname*{{\mathbb{E}}}_{x,y\sim\mathcal{D}_{k}}\left[\left\|\sum_{i=1}^{k}\phi^{\prime}(\sigma_{k}(U^{*})x_{i})\phi(\sigma_{k}(V^{*})y_{i})t_{i}\right\|^{2}\right].$

Without loss of generality, we assume $\sigma_{k}(U^{*})=\sigma_{k}(V^{*})=1$.
Then according to Lemma D.1 and Lemma D.2, we have

$\displaystyle\lambda_{\min}(H)\geq$ $\displaystyle\leavevmode\nobreak\
\frac{1}{\lambda(U^{*})\lambda(V^{*})\max\\{\kappa(U^{*}),\kappa(V^{*})\\}}$
$\displaystyle\leavevmode\nobreak\
\cdot\min\\{(\alpha_{2,0}\beta_{2,0}-\alpha_{1,0}^{2}\beta_{1,0}^{2}-\beta_{1,0}^{2}\alpha_{1,1}^{2}),(\alpha_{2,0}\beta_{2,2}-\alpha_{1,0}^{2}\beta_{1,2}^{2}-\gamma^{2})\\}.$

Considering the definition of $\rho$ in Eq. (9), we complete the proof. ∎

For the ReLU case, we lower bound the minimal eigenvalue of the Hessian for
non-orthogonal cases.

######  Theorem D.6.

Consider the activation to be ReLU. Assume $U^{*},V^{*}$ are full-column-rank
matrices and $u^{*}_{1,i}\neq 0,\forall i\in[k]$. Then the minimal eigenvalue
of the Hessian of Eq. (12) is lower bounded,

$\displaystyle\lambda_{\min}(\nabla^{2}f_{\mathcal{D}}^{\mathrm{ReLU}}(W^{*},V^{*}))\gtrsim\frac{1}{\lambda(U^{*})\lambda(V^{*})}\left(\frac{\min_{i\in[k]}\\{|u_{1,i}^{*}|\\}}{(1+\|u^{*(1)}\|)\max\\{\|U^{*}\|,\|V^{*}\|\\}}\right)^{2},$

where $u^{*(1)}$ is the first row of $U^{*}$.

###### Proof.

Let $P\in\mathbb{R}^{d\times k},Q\in\mathbb{R}^{d\times k}$ be the orthonormal
basis of $U^{*},V^{*}$ respectively. Let
$R\in\operatorname{{\mathbb{R}}}^{k\times
k},S\in\operatorname{{\mathbb{R}}}^{k\times k}$ satisfy that $U^{*}=P\cdot R$
and $V^{*}=Q\cdot S$. Let
$P_{\perp}\in\mathbb{R}^{d\times(d-k)},Q_{\perp}\in\mathbb{R}^{d\times(d-k)}$
be the orthogonal complement of $P,Q$ respectively. Set $a_{i}=P\cdot
s_{i}+P_{\perp}\cdot t_{i}$ and $b_{i}=Q\cdot p_{i}+Q_{\perp}\cdot q_{i}$.
Similar to the proof of Theorem D.5, Lemma D.2 and Lemma D.3, we have the
following.

$\displaystyle\leavevmode\nobreak\
\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(u_{i}^{*\top}x)\phi(v_{i}^{*\top}y)x^{\top}a_{i}+\phi^{\prime}(v_{i}^{*\top}y)\phi(u_{i}^{*\top}x)y^{\top}b_{i}\right)^{2}\right]$
$\displaystyle\geq$ $\displaystyle\leavevmode\nobreak\
\frac{1}{\lambda(U^{*})\lambda(V^{*})}\operatorname*{{\mathbb{E}}}_{x,y\sim\mathcal{D}_{k}}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(\sigma_{k}(U^{*})x_{i}))\phi(y_{i})x^{\top}R^{-1}s_{i}\sigma_{k}(U^{*})\right.\right.$
$\displaystyle\leavevmode\nobreak\
\left.\left.+\phi^{\prime}(\sigma_{k}(V^{*})y_{i})\phi(\sigma_{k}(U^{*})x_{i})y^{\top}S^{-1}p_{i}\sigma_{k}(V^{*})\right)^{2}\right]$
$\displaystyle\leavevmode\nobreak\
+\frac{1}{\lambda(U^{*})\lambda(V^{*})}\operatorname*{{\mathbb{E}}}_{x,y\sim\mathcal{D}_{k}}\left[\left\|\sum_{i=1}^{k}\phi^{\prime}(\sigma_{k}(U^{*})x_{i})\phi(\sigma_{k}(V^{*})y_{i})t_{i}\right\|^{2}\right]$
$\displaystyle\leavevmode\nobreak\
+\frac{1}{\lambda(U^{*})\lambda(V^{*})}\operatorname*{{\mathbb{E}}}_{x,y\sim\mathcal{D}_{k}}\left[\left\|\sum_{i=1}^{k}\phi^{\prime}(\sigma_{k}(U^{*})x_{i})\phi(\sigma_{k}(V^{*})y_{i})q_{i}\right\|^{2}\right]$
$\displaystyle\geq$ $\displaystyle\leavevmode\nobreak\
\frac{1}{16\lambda(U^{*})\lambda(V^{*})}(\|\widehat{A}_{o}\|_{F}^{2}+\|\widehat{B}_{o}\|_{F}^{2}+\|g_{\widehat{A}}+g_{\widehat{B}}\|^{2}+3(\|\widehat{T}\|_{F}^{2}+\|\widehat{Q}\|_{F}^{2})),$

where $\widehat{A}=[R^{-1}s_{1},R^{-1}s_{2},\cdots,R^{-1}s_{k}]$,
$\widehat{B}=[S^{-1}p_{1},S^{-1}p_{2},\cdots,S^{-1}p_{k}]$,
$\widehat{T}=[t_{1},t_{2},\cdots,t_{k}]$,
$\widehat{Q}=[q_{1},q_{2},\cdots,q_{k}]$.

Similar to Eq. (13), we can find the minimal eigenvalue of the Hessian by the
following constrained minimization problem.

$\displaystyle\lambda_{\min}(H)$
$\displaystyle=\min_{\begin{subarray}{c}\sum_{i=1}^{k}\|a_{i}\|^{2}+\|b_{i}\|^{2}=1\\\
a_{i,1}=0,\forall
i\in[k]\end{subarray}}\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\sum_{i=1}^{k}\phi^{\prime}(u_{i}^{*\top}x)\phi(v_{i}^{*\top}y)x^{\top}a_{i}+\phi^{\prime}(v_{i}^{*\top}y)\phi(u_{i}^{*\top}x)y^{\top}b_{i}\right)^{2}\right],$

which is lower bounded by the following formula.

$\displaystyle\underset{\widehat{A},\widehat{B},\widehat{T},\widehat{P}}{\min}$
$\displaystyle\leavevmode\nobreak\
\frac{1}{16\lambda(U^{*})\lambda(V^{*})}(\|\widehat{A}_{o}\|_{F}^{2}+\|\widehat{B}_{o}\|_{F}^{2}+\|g_{\widehat{A}}+g_{\widehat{B}}\|^{2}+3(\|\widehat{T}\|_{F}^{2}+\|\widehat{Q}\|_{F}^{2}))$
(15) s.t. $\displaystyle\leavevmode\nobreak\
\|R\widehat{A}\|_{F}^{2}+\|S\widehat{B}\|_{F}^{2}+\|\widehat{T}\|_{F}^{2}+\|\widehat{Q}\|_{F}^{2}=1$
$\displaystyle\leavevmode\nobreak\
e_{1}^{\top}PR\widehat{A}+e_{1}^{\top}P_{\perp}\widehat{T}=0$

If we assume the minimum of the above formula is $c_{1}$. We show that
$c_{1}>0$ by contradiction. If $c_{1}=0$, then $\widehat{T}=\widehat{Q}=0$,
$\widehat{A}_{o}=\widehat{B}_{o}=0$, $g_{\widehat{A}}=-g_{\widehat{B}}$. Since
$\widehat{T}=0$, we have
$e_{1}^{\top}PR\widehat{A}=e_{1}^{\top}U^{*}\widehat{A}=0$. Assuming
$(e_{1}^{\top}U^{*})_{i}\neq 0,\forall i$, we have
$g_{\widehat{A}}=g_{\widehat{B}}=0$. This violates the condition that
$\|R\widehat{A}\|_{F}^{2}+\|S\widehat{B}\|_{F}^{2}+\|\widehat{T}\|_{F}^{2}+\|\widehat{Q}\|_{F}^{2}=1$.

Now we give a lower bound for $c_{1}$. First we note,

$\displaystyle\|R\widehat{A}\|_{F}^{2}+\|S\widehat{B}\|_{F}^{2}+\|\widehat{T}\|_{F}^{2}+\|\widehat{Q}\|_{F}^{2}\leq\|R\|^{2}\|\widehat{A}\|^{2}_{F}+\|S\|^{2}\|\widehat{B}\|_{F}^{2}+\|\widehat{T}\|_{F}^{2}+\|\widehat{Q}\|_{F}^{2}.$

Therefore,

$\displaystyle\|\widehat{A}\|^{2}_{F}+\|\widehat{B}\|_{F}^{2}+\|\widehat{T}\|_{F}^{2}+\|\widehat{Q}\|_{F}^{2}\geq\frac{1}{\max\\{\|U^{*}\|^{2},\|V^{*}\|^{2}\\}}.$

Also, as $e_{1}^{\top}U^{*}\widehat{A}_{o}+(e_{1}^{\top}U^{*})\odot
g_{\widehat{A}}^{\top}+e_{1}^{\top}P_{\perp}\widehat{T}=0$, where $\odot$ is
the element-wise product, we have

$\displaystyle\|g_{\widehat{A}}\|^{2}$
$\displaystyle\leq(\frac{1}{\min\\{|u_{1,i}^{*}|\\}}(\|u^{*(1)}\|\|\widehat{A}_{o}\|+\|\widehat{T}\|)^{2}$
$\displaystyle\leq\left(\frac{1+\|u^{*(1)}\|}{\min\\{|u_{1,i}^{*}|\\}}\right)^{2}2(\|\widehat{A}_{o}\|_{F}^{2}+\|\widehat{T}\|_{F}^{2}).$

Note that
$\|g_{\widehat{A}}\|^{2}+\|g_{\widehat{A}}+g_{\widehat{B}}\|^{2}\geq\frac{1}{2}\|g_{\widehat{B}}\|^{2}$.
Now let’s return to the main part of objective function Eq. (15).

$\displaystyle\leavevmode\nobreak\
\|\widehat{A}_{o}\|_{F}^{2}+\|\widehat{B}_{o}\|_{F}^{2}+\|g_{\widehat{A}}+g_{\widehat{B}}\|^{2}+3(\|\widehat{T}\|_{F}^{2}+\|\widehat{Q}\|_{F}^{2})$
$\displaystyle\geq$ $\displaystyle\leavevmode\nobreak\
\frac{2}{3}(\|\widehat{A}_{o}\|_{F}^{2}+\|\widehat{T}\|_{F}^{2})+\frac{1}{3}\|\widehat{A}_{o}\|_{F}^{2}+\|\widehat{B}_{o}\|_{F}^{2}+\|g_{\widehat{A}}+g_{\widehat{B}}\|^{2}+\|\widehat{T}\|_{F}^{2}+\|\widehat{Q}\|_{F}^{2}$
$\displaystyle\geq$ $\displaystyle\leavevmode\nobreak\
\frac{1}{3}\left(\frac{\min\\{|u_{1,i}^{*}|\\}}{1+\|u^{*(1)}\|}\right)^{2}\|g_{\widehat{A}}\|^{2}+\frac{1}{3}\|\widehat{A}_{o}\|_{F}^{2}+\|\widehat{B}_{o}\|_{F}^{2}+\|g_{\widehat{A}}+g_{\widehat{B}}\|^{2}+\|\widehat{T}\|_{F}^{2}+\|\widehat{Q}\|_{F}^{2}$
$\displaystyle\geq$ $\displaystyle\leavevmode\nobreak\
\frac{1}{12}\left(\frac{\min\\{|u_{1,i}^{*}|\\}}{1+\|u^{*(1)}\|}\right)^{2}(\|g_{\widehat{A}}\|^{2}+\|g_{\widehat{B}}\|^{2})+\frac{1}{3}\|\widehat{A}_{o}\|_{F}^{2}+\|\widehat{B}_{o}\|_{F}^{2}+\|\widehat{T}\|_{F}^{2}+\|\widehat{Q}\|_{F}^{2}$
$\displaystyle\geq$ $\displaystyle\leavevmode\nobreak\
\frac{1}{12}\left(\frac{\min\\{|u_{1,i}^{*}|\\}}{1+\|u^{*(1)}\|}\right)^{2}\left(\|g_{\widehat{A}}\|^{2}+\|g_{\widehat{B}}\|^{2}+\|\widehat{A}_{o}\|_{F}^{2}+\|\widehat{B}_{o}\|_{F}^{2}+\|\widehat{T}\|_{F}^{2}+\|\widehat{Q}\|_{F}^{2}\right)$
$\displaystyle\geq$ $\displaystyle\leavevmode\nobreak\
\frac{1}{12}\left(\frac{\min\\{|u_{1,i}^{*}|\\}}{(1+\|u^{*(1)}\|)\max\\{\|U^{*}\|,\|V^{*}\|\\}}\right)^{2}.$

Therefore,

$\displaystyle c_{1}\geq$
$\displaystyle\frac{1}{200\lambda(U^{*})\lambda(V^{*})}\left(\frac{\min\\{|u_{1,i}^{*}|\\}}{(1+\|u^{*(1)}\|)\max\\{\|U^{*}\|,\|V^{*}\|\\}}\right)^{2}.$

∎

##  Appendix E Positive Definiteness of the Empirical Hessian

For any $(U,V)$, the population Hessian can be decomposed into the following
$2k\times 2k$ blocks ($i\in[k],j\in[k]$),

$\displaystyle\frac{\partial^{2}f_{\mathcal{D}}(U,V)}{\partial u_{i}\partial
u_{j}}=$ $\displaystyle\leavevmode\nobreak\
\operatorname*{{\mathbb{E}}}_{x,y}\left[\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(u_{j}^{\top}x)xx^{\top}\phi(v_{i}^{\top}y)\phi(v_{j}^{\top}y)\right]$
(16)
$\displaystyle+\delta_{ij}\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\phi(U^{\top}x)^{\top}\phi(V^{\top}y)-\phi(U^{*\top}x)^{\top}\phi(V^{*\top}y)\right)\phi^{\prime\prime}(u_{i}^{\top}x)\phi(v_{i}^{\top}y)xx^{\top}\right]$
$\displaystyle\frac{\partial^{2}f_{\mathcal{D}}(U,V)}{\partial u_{i}\partial
v_{j}}=$ $\displaystyle\leavevmode\nobreak\
\operatorname*{{\mathbb{E}}}_{x,y}\left[\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(v_{j}^{\top}y)xy^{\top}\phi(v_{i}^{\top}y)\phi(u_{j}^{\top}x)\right]$
$\displaystyle+\delta_{ij}\operatorname*{{\mathbb{E}}}_{x,y}\left[\left(\phi(U^{\top}x)^{\top}\phi(V^{\top}y)-\phi(U^{*\top}x)^{\top}\phi(V^{*\top}y)\right)\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(v_{i}^{\top}y)xy^{\top}\right],$

where $\delta_{ij}=1$ if $i=j$, otherwise $\delta_{ij}=0$. Similarly we can
write the formula for $\frac{\partial^{2}f_{\mathcal{D}}(U,V)}{\partial
v_{i}\partial v_{j}}$ and $\frac{\partial^{2}f_{\mathcal{D}}(U,V)}{\partial
v_{i}\partial u_{j}}$.

Replacing $\operatorname*{{\mathbb{E}}}_{x,y}$ by
$\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}$ in the above formula, we can obtain
the formula for the corresponding empirical Hessian,
$\nabla^{2}f_{\Omega}(U,V)$.

We now bound the difference between $\nabla^{2}f_{\Omega}(U,V)$ and
$\nabla^{2}f_{\mathcal{D}}(U^{*},V^{*})$.

######  Theorem E.1 (Restatement of Theorem C.2).

For any $\epsilon>0$, if

$\displaystyle n_{1}\gtrsim\epsilon^{-2}td\log^{2}d,\quad
n_{2}\gtrsim\epsilon^{-2}t\log
d,\quad|\Omega|\gtrsim\epsilon^{-2}td\log^{2}d,$

then with probability $1-d^{-t}$, for sigmoid/tanh,

$\|\nabla^{2}f_{\Omega}(U,V)-\nabla^{2}f_{\mathcal{D}}(U^{*},V^{*})\|\lesssim\epsilon+\|U-U^{*}\|+\|V-V^{*}\|,$

for ReLU,

$\|\nabla^{2}f_{\Omega}(U,V)-\nabla^{2}f_{\mathcal{D}}(U^{*},V^{*})\|\lesssim\left(\left(\frac{\|V-V^{*}\|}{\sigma_{k}(V^{*})}\right)^{1/2}+\left(\frac{\|U-U^{*}\|}{\sigma_{k}(U^{*})}\right)^{1/2}+\epsilon\right)(\|U^{*}\|+\|V^{*}\|)^{2}.$

###### Proof.

Define $H(U,V)\in\mathbb{R}^{(2kd)\times(2kd)}$ as a symmetric matrix, whose
blocks are represented as

$\displaystyle H_{u_{i},u_{j}}=$ $\displaystyle\leavevmode\nobreak\
\operatorname*{{\mathbb{E}}}_{x,y}\left[\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(u_{j}^{\top}x)xx^{\top}\phi(v_{i}^{\top}y)\phi(v_{j}^{\top}y)\right],$
(17) $\displaystyle H_{u_{i},v_{j}}=$ $\displaystyle\leavevmode\nobreak\
\operatorname*{{\mathbb{E}}}_{x,y}\left[\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(v_{j}^{\top}y)xy^{\top}\phi(v_{i}^{\top}y)\phi(u_{j}^{\top}x)\right].$

where $H_{u_{i},u_{j}}\in\mathbb{R}^{d\times
d},H_{u_{i},v_{j}}\in\mathbb{R}^{d\times d}$ correspond to
$\frac{\partial^{2}f_{\mathcal{D}}(U,V)}{\partial u_{i}\partial
u_{j}},\frac{\partial^{2}f_{\mathcal{D}}(U,V)}{\partial u_{i}\partial v_{j}}$
respectively.

We decompose the difference into

$\displaystyle\|\nabla^{2}f_{\Omega}(U,V)-\nabla^{2}f_{\mathcal{D}}(U^{*},V^{*})\|$
$\displaystyle\leq\|\nabla^{2}f_{\Omega}(U,V)-H(U,V)\|+\|H(U,V)-\nabla^{2}f_{\mathcal{D}}(U^{*},V^{*})\|.$

Combining Lemma E.2, E.14, we complete the proof. ∎

######  Lemma E.2.

For any $\epsilon>0$, if

$\displaystyle n_{1}\gtrsim\epsilon^{-2}td\log^{2}d,\quad
n_{2}\gtrsim\epsilon^{-2}t\log
d,\quad|\Omega|\gtrsim\epsilon^{-2}td\log^{2}d,$

then with probability $1-d^{-t}$, for sigmoid/tanh,

$\|\nabla^{2}f_{\Omega}(U,V)-H(U,V)\|\lesssim\epsilon+\|U-U^{*}\|+\|V-V^{*}\|,$

for ReLU,

$\|\nabla^{2}f_{\Omega}(U,V)-H(U,V)\|\lesssim\epsilon\|U^{*}\|\|V^{*}\|.$

###### Proof.

We can bound $\|\nabla^{2}f_{\Omega}(U,V)-H(U,V)\|$ if we bound each block.

We can show that if

$\displaystyle n_{1}\gtrsim\epsilon^{-2}td\log^{2}d,\quad
n_{2}\gtrsim\epsilon^{-2}t\log
d,\quad|\Omega|\gtrsim\epsilon^{-2}td\log^{2}d,$

then with probability $1-d^{-t}$,

$\displaystyle\leavevmode\nobreak\
\left\|\left(\operatorname*{{\mathbb{E}}}_{x,y}-\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}\right)\left[\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(u_{j}^{\top}x)xx^{\top}\phi(v_{i}^{\top}y)\phi(v_{j}^{\top}y)\right]\right\|$
$\displaystyle\lesssim$ $\displaystyle\leavevmode\nobreak\
\epsilon\|U^{*}\|^{p}\|V^{*}\|^{p}$ Lemma E.3
$\displaystyle\leavevmode\nobreak\
\left\|\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}\left[\left(\phi(U^{\top}x)^{\top}\phi(V^{\top}y)-\phi(U^{*\top}x)^{\top}\phi(V^{*\top}y)\right)\phi^{\prime\prime}(u_{i}^{\top}x)\phi(v_{i}^{\top}y)xx^{\top}\right]\right\|$
$\displaystyle\lesssim$ $\displaystyle\leavevmode\nobreak\
\|U-U^{*}\|+\|V-V^{*}\|$ Lemma E.6 $\displaystyle\leavevmode\nobreak\
\left\|\left(\operatorname*{{\mathbb{E}}}_{x,y}-\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}\right)\left[\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(v_{j}^{\top}y)xy^{\top}\phi(v_{i}^{\top}y)\phi(u_{j}^{\top}x)\right]\right\|$
$\displaystyle\lesssim$ $\displaystyle\leavevmode\nobreak\
\epsilon\|U^{*}\|^{p}\|V^{*}\|^{p}$ Lemma E.7
$\displaystyle\left\|\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}\left[\left(\phi(U^{\top}x)^{\top}\phi(V^{\top}y)-\phi(U^{*\top}x)^{\top}\phi(V^{*\top}y)\right)\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(v_{i}^{\top}y)xy^{\top}\right]\right\|$
$\displaystyle\lesssim$ $\displaystyle\leavevmode\nobreak\
\|U-U^{*}\|+\|V-V^{*}\|,$ Lemma E.9

where $p=1$ if $\phi$ is ReLU, $p=0$ if $\phi$ is sigmoid/tanh.

Note that for ReLU activation, for any given $U,V$, the second term is $0$
because $\phi^{\prime\prime}(z)=0$ almost everywhere. ∎

######  Lemma E.3.

If

$\displaystyle n_{1}\gtrsim\epsilon^{-2}td\log^{2}d,\quad
n_{2}\gtrsim\epsilon^{-2}t\log
d,\quad|\Omega|\gtrsim\epsilon^{-2}td\log^{2}d,$

then with probability at least $1-d^{-t}$,

$\displaystyle\left\|\left(\operatorname*{{\mathbb{E}}}_{x,y}-\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}\right)\left[\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(u_{j}^{\top}x)xx^{\top}\phi(v_{i}^{\top}y)\phi(v_{j}^{\top}y)\right]\right\|\leq\epsilon\|v_{i}\|^{p}\|v_{j}\|^{p}$

where $p=1$ if $\phi$ is $\mathrm{ReLU}$, $p=0$ if $\phi$ is sigmoid/tanh.

###### Proof.

Let
$B(x,y)=\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(u_{j}^{\top}x)xx^{\top}\phi(v_{i}^{\top}y)\phi(v_{j}^{\top}y)$.
By applying Lemma E.11 and Property $\mathrm{(\@slowromancap
i@)}-\mathrm{(\@slowromancap iii@)},\mathrm{(\@slowromancap vi@)}$ in Lemma
E.4 and Lemma E.5, we have for any $\epsilon>0$ if

$\displaystyle n_{1}\gtrsim\epsilon^{-2}td\log^{2}d,\quad
n_{2}\gtrsim\epsilon^{-2}t\log d,$

then with probability at least $1-d^{-2t}$,

$\displaystyle\left\|\operatorname*{{\mathbb{E}}}_{x,y}[B(x,y)]-\frac{1}{|S|}\sum_{(x,y)\in
S}B(x,y)\right\|\leq\epsilon\|v_{i}\|^{p}\|v_{j}\|^{p}.$ (18)

By applying Lemma E.12 and Property $\mathrm{(\@slowromancap
i@)},\mathrm{(\@slowromancap iii@)}-\mathrm{(\@slowromancap v@)}$ in Lemma E.4
and Lemma E.5, we have for any $\epsilon>0$ if

$\displaystyle n_{1}\gtrsim\epsilon^{-1}td\log^{2}d,\quad
n_{2}\gtrsim\epsilon^{-2}t\log d,$

then

$\displaystyle\left\|\frac{1}{n_{1}}\sum_{l\in[n_{1}]}(\phi^{\prime}(u_{i}^{\top}x_{l})\phi^{\prime}(u_{j}^{\top}x_{l}))^{2}\|x_{l}\|^{2}x_{l}x_{l}^{\top}\right\|\lesssim
d,$

and

$\displaystyle\left\|\frac{1}{n_{2}}\sum_{l\in[n_{2}]}(\phi(v_{i}^{\top}y_{l})\phi(v_{j}^{\top}y_{l}))^{2}\right\|\lesssim\|v_{i}\|^{2p}\|v_{j}\|^{2p}.$

Therefore,

$\max\left(\left\|\frac{1}{|S|}\sum_{(x,y)\in
S}B(x,y)B(x,y)^{\top}\right\|,\left\|\frac{1}{|S|}\sum_{(x,y)\in
S}B(x,y)^{\top}B(x,y)\right\|\right)\lesssim\epsilon
d\|v_{i}\|^{2p}\|v_{j}\|^{2p}.$ (19)

We can apply Lemma E.13 and use Eq. (19) and Property $\mathrm{(\@slowromancap
i@)}$ in Lemma E.4 and Lemma E.5 to obtain the following result. If

$\displaystyle|\Omega|\gtrsim\epsilon^{-2}td\log^{2}d,$

then with probability at least $1-d^{-2t}$,

$\left\|\frac{1}{|S|}\sum_{(x,y)\in
S}B(x,y)-\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}B(x,y)\right\|\lesssim\epsilon\|v_{i}\|^{p}\|v_{j}\|^{p}.$
(20)

Combining Eq. (18) and (20), we finish the proof. ∎

######  Lemma E.4.

Define
$T(z)=\phi^{\prime}(u_{i}^{\top}z)\phi^{\prime}(u_{j}^{\top}z)zz^{\top}$. If
$z\sim\mathcal{Z}$, $\mathcal{Z}=\mathcal{N}(0,I_{d})$ and $\phi$ is ReLU or
sigmoid/tanh, the following holds for $T(z)$ and any $t>1$,

$\displaystyle\mathrm{(\@slowromancap i@)}\quad$
$\displaystyle\quad\underset{z\sim{\cal Z}}{\Pr}\left[\left\|T(z)\right\|\leq
5td\log n\right]\geq 1-n^{-1}d^{-t};$ $\displaystyle\mathrm{(\@slowromancap
ii@)}\quad$ $\displaystyle\quad\max_{\|a\|=\|b\|=1}\left(\underset{z\sim{\cal
Z}}{\mathbb{E}}\left[\left(a^{\top}T(z)b\right)^{2}\right]\right)^{1/2}\lesssim
1;$ $\displaystyle\mathrm{(\@slowromancap iii@)}\quad$
$\displaystyle\quad\max\left(\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)T(z)^{\top}]\right\|,\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)^{\top}T(z)]\right\|\right)\lesssim d;$
$\displaystyle\mathrm{(\@slowromancap iv@)}\quad$
$\displaystyle\quad\max_{\|a\|=1}\left(\underset{z\sim{\cal
Z}}{\mathbb{E}}\left[\left(a^{\top}T(z)T(z)^{\top}a\right)^{2}\right]\right)^{1/2}\lesssim
d;$ $\displaystyle\mathrm{(\@slowromancap v@)}\quad$
$\displaystyle\quad\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)T(z)^{\top}T(z)T(z)^{\top}]\right\|\lesssim d^{3};$
$\displaystyle\mathrm{(\@slowromancap vi@)}\quad$
$\displaystyle\quad\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)]\right\|\lesssim 1.$

###### Proof.

Note that $0\leq\phi^{\prime}(z)\leq 1$, therefore $\mathrm{(\@slowromancap
i@)}$ can be proved by Proposition 1 of [HKZ12]. $\mathrm{(\@slowromancap
ii@)}-\mathrm{(\@slowromancap vi@)}$ can be proved by Hölder’s inequality. ∎

######  Lemma E.5.

Define $T(z)=\phi(v_{i}^{\top}z)\phi(v_{j}^{\top}z)$. If $z\sim\mathcal{Z}$,
$\mathcal{Z}=\mathcal{N}(0,I_{d})$ and $\phi$ is ReLU or sigmoid/tanh, the
following holds for $T(z)$ and any $t>1$,

$\displaystyle\mathrm{(\@slowromancap i@)}\quad$
$\displaystyle\quad\underset{z\sim{\cal Z}}{\Pr}\left[\left\|T(z)\right\|\leq
5t\|v_{i}\|^{p}\|v_{j}\|^{p}\log n\right]\geq 1-n^{-1}d^{-t};$
$\displaystyle\mathrm{(\@slowromancap ii@)}\quad$
$\displaystyle\quad\max_{\|a\|=\|b\|=1}\left(\underset{z\sim{\cal
Z}}{\mathbb{E}}\left[\left(a^{\top}T(z)b\right)^{2}\right]\right)^{1/2}\lesssim\|v_{i}\|^{p}\|v_{j}\|^{p};$
$\displaystyle\mathrm{(\@slowromancap iii@)}\quad$
$\displaystyle\quad\max\left(\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)T(z)^{\top}]\right\|,\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)^{\top}T(z)]\right\|\right)\lesssim\|v_{i}\|^{2p}\|v_{j}\|^{2p};$
$\displaystyle\mathrm{(\@slowromancap iv@)}\quad$
$\displaystyle\quad\max_{\|a\|=1}\left(\underset{z\sim{\cal
Z}}{\mathbb{E}}\left[\left(a^{\top}T(z)T(z)^{\top}a\right)^{2}\right]\right)^{1/2}\lesssim\|v_{i}\|^{2p}\|v_{j}\|^{2p};$
$\displaystyle\mathrm{(\@slowromancap v@)}\quad$
$\displaystyle\quad\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)T(z)^{\top}T(z)T(z)^{\top}]\right\|\lesssim\|v_{i}\|^{4p}\|v_{j}\|^{4p};$
$\displaystyle\mathrm{(\@slowromancap vi@)}\quad$
$\displaystyle\quad\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)]\right\|\lesssim\|v_{j}\|^{p}\|v_{i}\|^{p}.$

where $p=1$ if $\phi$ is ReLU, $p=0$ if $\phi$ is sigmoid/tanh.

###### Proof.

Note that $|\phi(z)|\leq|z|^{p}$, therefore $\mathrm{(\@slowromancap i@)}$ can
be proved by Proposition 1 of [HKZ12]. $\mathrm{(\@slowromancap
ii@)}-\mathrm{(\@slowromancap vi@)}$ can be proved by Hölder’s inequality ∎

######  Lemma E.6.

If

$\displaystyle n_{1}\gtrsim\epsilon^{-2}td\log^{2}d,\quad
n_{2}\gtrsim\epsilon^{-2}t\log
d,\quad|\Omega|\gtrsim\epsilon^{-2}td\log^{2}d,$

then with probability at least $1-d^{-t}$,

$\displaystyle\leavevmode\nobreak\
\left\|\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}\left[\left(\phi(U^{\top}x)^{\top}\phi(V^{\top}y)-\phi(U^{*\top}x)^{\top}\phi(V^{*\top}y)\right)\phi^{\prime\prime}(u_{i}^{\top}x)\phi(v_{i}^{\top}y)xx^{\top}\right]\right\|$
$\displaystyle\lesssim$ $\displaystyle\leavevmode\nobreak\
(\|U-U^{*}\|+\|V-V^{*}\|).$

###### Proof.

We consider the following formula first,

$\displaystyle\left\|\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}\left[\left((\phi(u_{j}^{\top}x)-\phi(u_{j}^{*\top}x))\phi(v_{j}^{*\top}y)\right)\phi^{\prime\prime}(u_{i}^{\top}x)\phi(v_{i}^{\top}y)xx^{\top}\right]\right\|$
$\displaystyle\leq$
$\displaystyle\left\|\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}\left[\left|(u_{j}-u_{j}^{*})^{\top}x\right|xx^{\top}\phi(v_{j}^{*\top}y)\phi(v_{i}^{\top}y)\right]\right\|.$

Similar to Lemma E.3, we are able to show

$\displaystyle\leavevmode\nobreak\
\left\|\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}\left[\left|(u_{j}-u_{j}^{*})^{\top}x\right|xx^{\top}\phi(v_{j}^{*\top}y)\phi(v_{i}^{\top}y)\right]-\operatorname*{{\mathbb{E}}}_{(x,y)}\left[\left|(u_{j}-u_{j}^{*})^{\top}x\right|xx^{\top}\phi(v_{j}^{*\top}y)\phi(v_{i}^{\top}y)\right]\right\|$
$\displaystyle\lesssim$ $\displaystyle\leavevmode\nobreak\ \|U-U^{*}\|.$

Note that by Hölder’s inequality, we have,

$\displaystyle\left\|\operatorname*{{\mathbb{E}}}_{(x,y)}\left[\left|(u_{j}-u_{j}^{*})^{\top}x\right|xx^{\top}\phi(v_{j}^{*\top}y)\phi(v_{i}^{\top}y)\right]\right\|\lesssim\|U-U^{*}\|.$

So we complete the proof. ∎

######  Lemma E.7.

If

$\displaystyle n_{1}\gtrsim\epsilon^{-2}td\log^{2}d,\quad
n_{2}\gtrsim\epsilon^{-2}t\log
d,\quad|\Omega|\gtrsim\epsilon^{-2}td\log^{2}d,$

then with probability at least $1-d^{-t}$,

$\displaystyle\left\|\left(\operatorname*{{\mathbb{E}}}_{x,y}-\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}\right)\left[\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(v_{j}^{\top}y)xy^{\top}\phi(v_{i}^{\top}y)\phi(u_{j}^{\top}x)\right]\right\|\lesssim\epsilon\|v_{i}\|^{p}\|u_{j}\|^{p}.$

###### Proof.

Let $B(x,y)=M(x)N(y)$, where
$M(x)=\phi^{\prime}(u_{i}^{\top}x)\phi(u_{j}^{\top}x)x$ and
$N(y)=\phi^{\prime}(v_{j}^{\top}y)\phi(v_{i}^{\top}y)y^{\top}$. By applying
Lemma E.11 and Property $\mathrm{(\@slowromancap i@)}-\mathrm{(\@slowromancap
iii@)},\mathrm{(\@slowromancap vi@)}$ in Lemma E.8 , we have for any
$\epsilon>0$ if

$\displaystyle n_{1}\gtrsim\epsilon^{-2}td\log^{2}d,\quad
n_{2}\gtrsim\epsilon^{-2}td\log^{2}d,$

then with probability at least $1-d^{-2t}$,

$\left\|\operatorname*{{\mathbb{E}}}_{x,y}B(x,y)-\frac{1}{|S|}\sum_{(x,y)\in
S}B(x,y)\right\|\lesssim\epsilon\|u_{j}\|^{p}\|v_{i}\|^{p}.$ (21)

By applying Lemma E.12 and Property $\mathrm{(\@slowromancap
i@)},\mathrm{(\@slowromancap iv@)}-\mathrm{(\@slowromancap vi@)}$ in Lemma
E.8, we have for any $\epsilon>0$ if

$\displaystyle
n_{1}\gtrsim\epsilon^{-2}td\log^{2}d,n_{2}\gtrsim\epsilon^{-2}td\log^{2}d,$

then

$\displaystyle\left\|\frac{1}{n_{1}}\sum_{l\in[n_{1}]}M(x_{l})M(x_{l})^{\top}\right\|\lesssim\|u_{j}\|^{2p},\quad\left\|\frac{1}{n_{2}}\sum_{l\in[n_{2}]}N(y_{l})^{\top}N(y_{l})\right\|\lesssim\|v_{i}\|^{2p}.$

By applying Lemma E.12 and Property $\mathrm{(\@slowromancap
i@)},\mathrm{(\@slowromancap iv@)},\mathrm{(\@slowromancap
vii@)},\mathrm{(\@slowromancap viii@)}$ in Lemma E.8, we have for any
$\epsilon>0$ if

$\displaystyle
n_{1}\gtrsim\epsilon^{-2}td\log^{2}d,n_{2}\gtrsim\epsilon^{-2}td\log^{2}d,$

then

$\displaystyle\left\|\frac{1}{n_{1}}\sum_{l\in[n_{1}]}M(x_{l})^{\top}M(x_{l})\right\|\lesssim
d\|u_{j}\|^{2p},\quad\left\|\frac{1}{n_{2}}\sum_{l\in[n_{2}]}N(y_{l})N(y_{l})^{\top}\right\|\lesssim
d\|v_{i}\|^{2p}.$

Therefore,

$\max\left(\left\|\frac{1}{|S|}\sum_{(x,y)\in
S}B(x,y)B(x,y)^{\top}\right\|,\left\|\frac{1}{|S|}\sum_{(x,y)\in
S}B(x,y)^{\top}B(x,y)\right\|\right)\lesssim\epsilon
d\|v_{i}\|^{2p}\|u_{j}\|^{2p}$ (22)

We can apply Lemma E.13 and Eq. (22) and Property $\mathrm{(\@slowromancap
i@)}$ in Lemma E.8 to obtain the following result. If

$\displaystyle|\Omega|\gtrsim\epsilon^{-2}td\log^{2}d,$

then with probability at least $1-d^{-2t}$,

$\left\|\frac{1}{|S|}\sum_{(x,y)\in
S}B(x,y)-\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}B(x,y)\right\|\leq\epsilon\|v_{i}\|^{p}\|u_{j}\|^{p}.$
(23)

Combining Eq. (21) and (23), we finish the proof. ∎

######  Lemma E.8.

Define $T(z)=\phi^{\prime}(u_{i}^{\top}z)\phi(u_{j}^{\top}z)z$. If
$z\sim\mathcal{Z}$, $\mathcal{Z}=\mathcal{N}(0,I_{d})$ and $\phi$ is ReLU or
sigmoid/tanh, the following holds for $T(z)$ and any $t>1$,

$\displaystyle\mathrm{(\@slowromancap i@)}\quad$
$\displaystyle\quad\underset{z\sim{\cal Z}}{\Pr}\left[\left\|T(z)\right\|\leq
5td^{1/2}\|u_{j}\|^{p}\log n\right]\geq 1-n^{-1}d^{-t};$
$\displaystyle\mathrm{(\@slowromancap ii@)}\quad$
$\displaystyle\quad\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)]\right\|\lesssim\|u_{j}\|^{p};$
$\displaystyle\mathrm{(\@slowromancap iii@)}\quad$
$\displaystyle\quad\max_{\|a\|=\|b\|=1}\left(\underset{z\sim{\cal
Z}}{\mathbb{E}}\left[\left(a^{\top}T(z)b\right)^{2}\right]\right)^{1/2}\lesssim\|u_{j}\|^{p};$
$\displaystyle\mathrm{(\@slowromancap iv@)}\quad$
$\displaystyle\quad\max\left\\{\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)T(z)^{\top}]\right\|,\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)^{\top}T(z)]\right\|\right\\}\lesssim d\|u_{j}\|^{2p};$
$\displaystyle\mathrm{(\@slowromancap v@)}\quad$
$\displaystyle\quad\max_{\|a\|=1}\left(\underset{z\sim{\cal
Z}}{\mathbb{E}}\left[\left(a^{\top}T(z)T(z)^{\top}a\right)^{2}\right]\right)^{1/2}\lesssim\|u_{j}\|^{2p};$
$\displaystyle\mathrm{(\@slowromancap vi@)}\quad$
$\displaystyle\quad\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)T(z)^{\top}T(z)T(z)^{\top}]\right\|\lesssim
d\|u_{j}\|^{4p};$ $\displaystyle\mathrm{(\@slowromancap vii@)}\quad$
$\displaystyle\quad\max_{\|a\|=1}\left(\underset{z\sim{\cal
Z}}{\mathbb{E}}\left[\left(a^{\top}T(z)^{\top}T(z)a\right)^{2}\right]\right)^{1/2}\lesssim
d\|u_{j}\|^{2p};$ $\displaystyle\mathrm{(\@slowromancap viii@)}\quad$
$\displaystyle\quad\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)^{\top}T(z)T(z)^{\top}T(z)]\right\|\lesssim
d^{2}\|u_{j}\|^{4p}.$

###### Proof.

Note that $0\leq\phi^{\prime}(z)\leq 1$,$|\phi(z)|\leq|z|^{p}$, therefore
$\mathrm{(\@slowromancap i@)}$ can be proved by Proposition 1 of [HKZ12].
$\mathrm{(\@slowromancap ii@)}-\mathrm{(\@slowromancap viii@)}$ can be proved
by Hölder’s inequality. ∎

######  Lemma E.9.

If

$\displaystyle n_{1}\gtrsim td\log^{2}d,\quad n_{2}\gtrsim t\log
d,\quad|\Omega|\gtrsim td\log^{2}d,$

then with probability at least $1-d^{-t}$,

$\displaystyle\leavevmode\nobreak\
\left\|\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}\left[\left(\phi(U^{\top}x)^{\top}\phi(V^{\top}y)-\phi(U^{*\top}x)^{\top}\phi(V^{*\top}y)\right)\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(v_{i}^{\top}y)xy^{\top}\right]\right\|$
$\displaystyle\lesssim$ $\displaystyle\leavevmode\nobreak\
\|U-U^{*}\|+\|V-V^{*}\|.$

###### Proof.

We consider the following formula first,

$\displaystyle\leavevmode\nobreak\
\left\|\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}\left[\left((\phi(u_{j}^{\top}x)-\phi(u_{j}^{*\top}x))\phi(v_{j}^{*\top}y)\right)\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(v_{i}^{\top}y)xy^{\top}\right]\right\|$

Set
$M(x)=(\phi(u_{j}^{\top}x)-\phi(u_{j}^{*\top}x))\phi^{\prime}(u_{i}^{\top}x)x$
and $N(y)=\phi(v_{j}^{*\top}y)\phi^{\prime}(v_{i}^{\top}y)y^{\top}$ and follow
the proof for Lemma E.7. Also note that $\phi$ is Lipschitz, i.e.,
$|\phi(u_{j}^{\top}x)-\phi(u_{j}^{*\top}x)|\leq|u_{j}^{\top}x-u_{j}^{*\top}x|$.
We can show the following. If

$\displaystyle n_{1}\gtrsim td\log^{2}d,\quad n_{2}\gtrsim t\log
d,\quad|\Omega|\gtrsim td\log^{2}d,$

then with probability at least $1-d^{-t}$,

$\displaystyle\left\|\left(\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}-\operatorname*{{\mathbb{E}}}_{x,y}\right)\left[M(x)N(y)\right]\right\|\lesssim\|u_{j}-u_{j}^{*}\|.$

Note that by Hölder’s inequality, we have,

$\displaystyle\|\operatorname*{{\mathbb{E}}}_{x,y}\left[M(x)N(y)\right]\|\lesssim\|u_{j}-u_{j}^{*}\|.$

So we complete the proof.

∎

We provide a variation of Lemma B.7 in [ZSJ+17]. Note that the Lemma B.7
[ZSJ+17] requires four properties, we simplify it into three properties.

######  Lemma E.10 (Matrix Bernstein for unbounded case (A modified version of
bounded case, Theorem 6.1 in [Tro12], A variation of Lemma B.7 in [ZSJ+17])).

Let ${\cal B}$ denote a distribution over $\mathbb{R}^{d_{1}\times d_{2}}$.
Let $d=d_{1}+d_{2}$. Let $B_{1},B_{2},\cdots B_{n}$ be i.i.d. random matrices
sampled from ${\cal B}$. Let $\overline{B}=\mathbb{E}_{B\sim{\cal B}}[B]$ and
$\widehat{B}=\frac{1}{n}\sum_{i=1}^{n}B_{i}$. For parameters $m\geq
0,\gamma\in(0,1),\nu>0,L>0$, if the distribution ${\cal B}$ satisfies the
following four properties,

$\displaystyle\mathrm{(\@slowromancap i@)}\quad$
$\displaystyle\quad\underset{B\sim{\cal B}}{\Pr}\left[\left\|B\right\|\leq
m\right]\geq 1-\gamma;$ $\displaystyle\mathrm{(\@slowromancap ii@)}\quad$
$\displaystyle\quad\max\left(\left\|\underset{B\sim{\cal
B}}{\mathbb{E}}[BB^{\top}]\right\|,\left\|\underset{B\sim{\cal
B}}{\mathbb{E}}[B^{\top}B]\right\|\right)\leq\nu;$
$\displaystyle\mathrm{(\@slowromancap iii@)}\quad$
$\displaystyle\quad\max_{\|a\|=\|b\|=1}\left(\underset{B\sim{\cal
B}}{\mathbb{E}}\left[\left(a^{\top}Bb\right)^{2}\right]\right)^{1/2}\leq L.$

Then we have for any $\epsilon>0$ and $t\geq 1$, if

$n\geq(18t\log
d)\cdot((\epsilon+\|\overline{B}\|)^{2}+m\epsilon+\nu)/\epsilon^{2}\quad\text{\leavevmode\nobreak\
and\leavevmode\nobreak\ }\quad\gamma\leq(\epsilon/(2L))^{2}$

with probability at least $1-d^{-2t}-n\gamma$,

$\left\|\frac{1}{n}\sum_{i=1}^{n}B_{i}-\underset{B\sim{\cal
B}}{\mathbb{E}}[B]\right\|\leq\epsilon.$

###### Proof.

Define the event

$\xi_{i}=\\{\|B_{i}\|\leq m\\},\forall i\in[n].$

Define $M_{i}={\mathbf{1}}_{\|B_{i}\|\leq m}B_{i}$. Let
$\overline{M}=\mathbb{E}_{B\sim{\cal B}}[{\mathbf{1}}_{\|B\|\leq m}B]$ and
$\widehat{M}=\frac{1}{n}\sum_{i=1}^{n}M_{i}$. By triangle inequality, we have

$\displaystyle\|\widehat{B}-\overline{B}\|\leq\|\widehat{B}-\widehat{M}\|+\|\widehat{M}-\overline{M}\|+\|\overline{M}-\overline{B}\|.$
(24)

In the next a few paragraphs, we will upper bound the above three terms.

The first term in Eq. (24). For each $i$, let $\overline{\xi}_{i}$ denote the
complementary set of $\xi_{i}$, i.e.
$\overline{\xi}_{i}=[n]\backslash\xi_{i}$. Thus
$\Pr[\overline{\xi}_{i}]\leq\gamma$. By a union bound over $i\in[n]$, with
probability $1-n\gamma$, $\|B_{i}\|\leq m$ for all $i\in[n]$. Thus
$\widehat{M}=\widehat{B}$.

The second term in Eq. (24). For a matrix $B$ sampled from ${\cal B}$, we use
$\xi$ to denote the event that $\xi=\\{\|B\|\leq m\\}$. Then, we can upper
bound $\|\overline{M}-\overline{B}\|$ in the following way,

$\displaystyle\leavevmode\nobreak\ \|\overline{M}-\overline{B}\|$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\left\|\underset{B\sim{\cal B}}{\mathbb{E}}[{\mathbf{1}}_{\|B\|\leq m}\cdot
B]-\underset{B\sim{\cal B}}{\mathbb{E}}[B]\right\|$ $\displaystyle=$
$\displaystyle\leavevmode\nobreak\ \left\|\underset{B\sim{\cal
B}}{\mathbb{E}}\left[B\cdot{\mathbf{1}}_{\overline{\xi}}\right]\right\|$
$\displaystyle=$ $\displaystyle\leavevmode\nobreak\
\max_{\|a\|=\|b\|=1}\underset{B\sim{\cal
B}}{\mathbb{E}}\left[a^{\top}Bb{\mathbf{1}}_{\overline{\xi}}\right]$
$\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
\max_{\|a\|=\|b\|=1}\underset{B\sim{\cal
B}}{\mathbb{E}}[(a^{\top}Bb)^{2}]^{1/2}\cdot\underset{B\sim{\cal
B}}{\mathbb{E}}\left[{\mathbf{1}}_{\overline{\xi}}\right]^{1/2}$ by Hölder’s
inequality $\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\
L\underset{B\sim{\cal
B}}{\mathbb{E}}\left[{\mathbf{1}}_{\overline{\xi}}\right]^{1/2}$ by Property
(iv@) $\displaystyle\leq$ $\displaystyle\leavevmode\nobreak\ L\gamma^{1/2},$
$\displaystyle\text{\leavevmode\nobreak\ by\leavevmode\nobreak\
}\Pr[\overline{\xi}]\leq\gamma$ $\displaystyle\leq$
$\displaystyle\leavevmode\nobreak\ \frac{1}{2}\epsilon,$
$\displaystyle\text{\leavevmode\nobreak\ by\leavevmode\nobreak\
}\gamma\leq(\epsilon/(2L))^{2},$

which is

$\displaystyle\|\overline{M}-\overline{B}\|\leq\frac{\epsilon}{2}.$

Therefore, $\|\overline{M}\|\leq\|\overline{B}\|+\frac{\epsilon}{2}$.

The third term in Eq. (24). We can bound $\|\widehat{M}-\overline{M}\|$ by
Matrix Bernstein’s inequality [Tro12].

We define $Z_{i}=M_{i}-\overline{M}$. Thus we have $\underset{B_{i}\sim{\cal
B}}{\mathbb{E}}[Z_{i}]=0$, $\|Z_{i}\|\leq 2m$, and

$\displaystyle\left\|\underset{B_{i}\sim{\cal
B}}{\mathbb{E}}[Z_{i}Z_{i}^{\top}]\right\|=\left\|\underset{B_{i}\sim{\cal
B}}{\mathbb{E}}[M_{i}M_{i}^{\top}]-\overline{M}\cdot\overline{M}^{\top}\right\|\leq\nu+\|\overline{M}\|^{2}\leq\nu+\|\overline{B}\|^{2}+\epsilon^{2}+\epsilon\|\overline{B}\|.$

Similarly, we have $\left\|\underset{B_{i}\sim{\cal
B}}{\mathbb{E}}[Z_{i}^{\top}Z_{i}]\right\|\leq\nu+\|\overline{B}\|^{2}+\epsilon^{2}+\epsilon\|\overline{B}\|$.
Using matrix Bernstein’s inequality, for any $\epsilon>0$,

$\displaystyle\underset{B_{1},\cdots,B_{n}\sim{\cal
B}}{\Pr}\left[\frac{1}{n}\left\|\sum_{i=1}^{n}Z_{i}\right\|\geq\epsilon\right]\leq
d\exp\left(-\frac{\epsilon^{2}n/2}{\nu+\|\overline{B}\|^{2}+\epsilon^{2}+\epsilon\|\overline{B}\|+2m\epsilon/3}\right).$

By choosing

$\displaystyle n\geq(3t\log
d)\cdot\frac{\nu+\|\overline{B}\|^{2}+\epsilon^{2}+\epsilon\|\overline{B}\|+2m\epsilon/3}{\epsilon^{2}/2},$

for $t\geq 1$, we have with probability at least $1-d^{-2t}$,

$\displaystyle\left\|\frac{1}{n}\sum_{i=1}^{n}M_{i}-\overline{M}\right\|\leq\frac{\epsilon}{2}.$

Putting it all together, we have for $\epsilon>0$, if

$n\geq(18t\log
d)\cdot((\epsilon+\|\overline{B}\|)^{2}+m\epsilon+\nu)/(\epsilon^{2})\quad\text{\leavevmode\nobreak\
and\leavevmode\nobreak\ }\quad\gamma\leq(\epsilon/(2L))^{2}$

with probability at least $1-d^{-2t}-n\gamma$,

$\left\|\frac{1}{n}\sum_{i=1}^{n}B_{i}-\underset{B\sim{\cal
B}}{\mathbb{E}}[B]\right\|\leq\epsilon.$

∎

######  Lemma E.11 (Tail Bound for fully-observed rating matrix).

Let $\\{x_{i}\\}_{i\in[n_{1}]}$ be independent samples from distribution
$\mathcal{X}$ and $\\{y_{j}\\}_{j\in[n_{2}]}$ be independent samples from
distribution $\mathcal{Y}$. Denote
$S:=\\{(x_{i},y_{j})\\}_{i\in[n_{1}],j\in[n_{2}]}$ as the collection of all
the $(x_{i},y_{j})$ pairs. Let $B(x,y)$ be a random matrix of $x,y$, which can
be represented as the product of two matrices $M(x),N(y)$, i.e.,
$B(x,y)=M(x)N(y)$. Let $\overline{M}=\operatorname*{{\mathbb{E}}}_{x}M(x)$ and
$\overline{N}=\operatorname*{{\mathbb{E}}}_{y}N(y)$. Let $d_{x}$ be the sum of
the two dimensions of $M(x)$ and $d_{y}$ be the sum of the two dimensions of
$N(y)$. Suppose both $M(x)$ and $N(y)$ satisfy the following properties ($z$
is a representative for $x,y$, and $T(z)$ is a representative for
$M(x),N(y)$),

$\displaystyle\mathrm{(\@slowromancap i@)}\quad$
$\displaystyle\quad\underset{z\sim{\cal Z}}{\Pr}\left[\left\|T(z)\right\|\leq
m_{z}\right]\geq 1-\gamma_{z};$ $\displaystyle\mathrm{(\@slowromancap
ii@)}\quad$ $\displaystyle\quad\max_{\|a\|=\|b\|=1}\left(\underset{z\sim{\cal
Z}}{\mathbb{E}}\left[\left(a^{\top}T(z)b\right)^{2}\right]\right)^{1/2}\leq
L_{z};$ $\displaystyle\mathrm{(\@slowromancap iii@)}\quad$
$\displaystyle\quad\max\left(\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)T(z)^{\top}]\right\|,\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)^{\top}T(z)]\right\|\right)\leq\nu_{z}.$

Then for any $\epsilon_{1}>0,\epsilon_{2}>0$ if

$\displaystyle n_{1}\geq(18t\log
d_{x})\cdot(\nu_{x}+(\|\overline{M}\|+\epsilon_{1})^{2}+m_{x}\epsilon_{1})/\epsilon_{1}^{2}\quad\text{\leavevmode\nobreak\
and\leavevmode\nobreak\ }\quad\gamma_{x}\leq(\epsilon_{1}/(2L_{x}))^{2}$
$\displaystyle n_{2}\geq(18t\log
d_{y})\cdot(\nu_{y}+(\epsilon_{2}+\|\overline{N}\|)^{2}+m_{y}\epsilon_{2})/\epsilon_{2}^{2}\quad\text{\leavevmode\nobreak\
and\leavevmode\nobreak\ }\quad\gamma_{y}\leq(\epsilon_{2}/(2L_{y}))^{2}$

with probability at least
$1-d_{x}^{-2t}-d_{y}^{-2t}-n_{1}\gamma_{x}-n_{2}\gamma_{y}$,

$\left\|\operatorname*{{\mathbb{E}}}_{x,y}B(x,y)-\frac{1}{|S|}\sum_{(x,y)\in
S}B(x,y)\right\|\leq\epsilon_{2}\|\overline{M}\|+\epsilon_{1}\|\overline{N}\|+\epsilon_{1}\epsilon_{2}.$
(25)

###### Proof.

First we note that,

$\frac{1}{|S|}\sum_{(x,y)\in
S}B(x,y)=\frac{1}{n_{1}n_{2}}\left(\sum_{i\in[n_{1}]}M(x_{i})\right)\cdot\left(\sum_{j\in[n_{2}]}N(y_{j})\right),$

and

$\operatorname*{{\mathbb{E}}}_{x,y}[B(x,y)]=\left(\operatorname*{{\mathbb{E}}}_{x}[M(x)]\right)\left(\operatorname*{{\mathbb{E}}}_{y}[N(y)]\right).$

Therefore, if we can bound
$\|\operatorname*{{\mathbb{E}}}_{x}[M(x)]-\frac{1}{n_{1}}\sum_{i\in[n_{1}]}M(x_{i})\|$
and the corresponding term for $y$, we are able to prove this lemma.

By the conditions of $M(x)$, the three conditions in Lemma E.10 are satisfied,
which completes the proof.

∎

######  Lemma E.12 (Upper bound for the second-order moment).

Let $\\{z_{i}\\}_{i\in[n]}$ be independent samples from distribution
$\mathcal{Z}$. Let $T(z)$ be a matrix of $z$. Let $d$ be the sum of the two
dimensions of $T(z)$ and $\overline{T}:=\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)T(z)^{\top}]$. Suppose $T(z)$ satisfies the following
properties.

$\displaystyle\mathrm{(\@slowromancap i@)}\quad$
$\displaystyle\quad\underset{z\sim{\cal Z}}{\Pr}\left[\left\|T(z)\right\|\leq
m_{z}\right]\geq 1-\gamma_{z};$ $\displaystyle\mathrm{(\@slowromancap
ii@)}\quad$ $\displaystyle\quad\max_{\|a\|=1}\left(\underset{z\sim{\cal
Z}}{\mathbb{E}}\left[\left(a^{\top}T(z)T(z)^{\top}a\right)^{2}\right]\right)^{1/2}\leq
L_{z};$ $\displaystyle\mathrm{(\@slowromancap iii@)}\quad$
$\displaystyle\quad\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)T(z)^{\top}T(z)T(z)^{\top}]\right\|\leq\nu_{z},$

Then for any $t>1$, if

$\displaystyle n\geq(18t\log
d)\cdot(\nu_{z}+(\|\overline{T}\|+\epsilon)^{2}+m^{2}_{z})/\epsilon^{2}\quad\text{\leavevmode\nobreak\
and\leavevmode\nobreak\ }\quad\gamma_{z}\leq(\epsilon/(2L_{z}))^{2},$

we have with probability at least $1-d^{-2t}-n\gamma_{z}$,

$\displaystyle\left\|\frac{1}{n}\sum_{i\in[n]}T(z_{i})T(z_{i})^{\top}\right\|\leq\left\|\underset{z\sim{\cal
Z}}{\mathbb{E}}[T(z)T(z)^{\top}]\right\|+\epsilon.$

###### Proof.

The proof directly follows by applying Lemma E.10. ∎

######  Lemma E.13 (Tail Bound for partially-observed rating matrix).

Given $\\{x_{i}\\}_{i\in[n_{1}]}$ and $\\{y_{j}\\}_{j\in[n_{2}]}$, let’s
denote $S:=\\{(x_{i},y_{j})\\}_{i\in[n_{1}],j\in[n_{2}]}$ as the collection of
all the $(x_{i},y_{j})$ pairs. Let $\Omega$ also be a collection of
$(x_{i},y_{j})$ pairs, where each pair is sampled from $S$ independently and
uniformly. Let $B(x,y)$ be a matrix of $x,y$. Let $d_{B}$ be the sum of the
two dimensions of $B(x,y)$. Define
$\overline{B}_{S}=\frac{1}{|S|}\sum_{(x,y)\in S}B(x,y)$. Assume the following,

$\displaystyle\mathrm{(\@slowromancap i@)}\quad$
$\displaystyle\quad\|B(x,y)\|\leq m_{B},\forall(x,y)\in S,$
$\displaystyle\mathrm{(\@slowromancap ii@)}\quad$
$\displaystyle\quad\max\left(\left\|\frac{1}{|S|}\sum_{(x,y)\in
S}B(x,y)B(x,y)^{\top}\right\|,\;\left\|\frac{1}{|S|}\sum_{(x,y)\in
S}B(x,y)^{\top}B(x,y)\right\|\right)\leq\nu_{B}.$

Then we have for any $\epsilon>0$ and $t\geq 1$, if

$\displaystyle|\Omega|\geq(18t\log
d_{B})\cdot(\nu_{B}+\|\overline{B}_{S}\|^{2}+m_{B}\epsilon)/\epsilon^{2},$

with probability at least $1-d_{B}^{-2t}$,

$\displaystyle\left\|\overline{B}_{S}-\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}B(x,y)\right\|\leq\epsilon.$

###### Proof.

Since each entry in $\Omega$ is sampled from $S$ uniformly and independently,
we have

$\displaystyle\operatorname*{{\mathbb{E}}}_{\Omega}\left[\frac{1}{|\Omega|}\sum_{(x,y)\in\Omega}B(x,y)\right]=\frac{1}{|S|}\sum_{(x,y)\in
S}B(x,y).$

Applying the matrix Bernstein inequality Theorem 6.1 in [Tro12], we prove this
lemma. ∎

######  Lemma E.14.

For sigmoid/tanh activation function,

$\displaystyle\|H(U,V)-\nabla^{2}f_{\mathcal{D}}(U^{*},V^{*})\|\lesssim(\|V-V^{*}\|+\|U-U^{*}\|),$

where $H(U,V)$ is defined as in Eq. (17).

For ReLU activation function,

$\displaystyle\|H(U,V)-\nabla^{2}f_{\mathcal{D}}(U^{*},V^{*})\|\lesssim\left(\left(\frac{\|V-V^{*}\|}{\sigma_{k}(V^{*})}\right)^{1/2}\|U^{*}\|+\left(\frac{\|U-U^{*}\|}{\sigma_{k}(U^{*})}\right)^{1/2}\|V^{*}\|\right)(\|U^{*}\|+\|V^{*}\|).$

###### Proof.

We can bound each block, i.e.,

$\displaystyle\operatorname*{{\mathbb{E}}}_{x,y}\left[\phi^{\prime}(u_{i}^{\top}x)\phi^{\prime}(u_{j}^{\top}x)xx^{\top}\phi(v_{i}^{\top}y)\phi(v_{j}^{\top}y)-\phi^{\prime}(u_{i}^{*\top}x)\phi^{\prime}(u_{j}^{*\top}x)xx^{\top}\phi(v_{i}^{*\top}y)\phi(v_{j}^{*\top}y)\right].$
(26)