Year: 2022

Abstract:
Self-training, a semi-supervised learning algorithm, leverages a large amount of unlabeled data to improve learning when the labeled data are limited. Despite empirical successes, its theoretical characterization remains elusive. To the best of our knowledge, this work establishes the first theoretical analysis for the known iterative self-training paradigm and proves the benefits of unlabeled data in both training convergence and generalization ability. To make our theoretical analysis feasible, we focus on the case of one-hidden-layer neural networks. However, theoretical understanding of iterative self-training is non-trivial even for a shallow neural network. One of the key challenges is that existing neural network landscape analysis built upon supervised learning no longer holds in the (semisupervised) self-training paradigm. We address this challenge and prove that iterative self-training converges linearly with both convergence rate and generalization accuracy improved in the order of 1/ √ M , where M is the number of unlabeled samples. Experiments from shallow neural networks to deep neural networks are also provided to justify the correctness of our established theoretical insights on self-training.

text: √ M , where M is the number of unlabeled
samples. Experiments from shallow neural networks to deep neural networks are also provided to justify the correctness of our established theoretical insights on self-training.

heading: 1 INTRODUCTION
text: Self-training (Scudder, 1965; Yarowsky, 1995; Lee et al., 2013; Han et al., 2019), one of the most powerful semi-supervised learning (SemiSL) algorithms, augments a limited number of labeled data with unlabeled data so as to achieve improved generalization performance on test data, compared with the model trained by supervised learning using the labeled data only. Self-training has shown empirical success in diversified applications such as few-shot image classification (Su et al., 2020; Xie et al., 2020; Chen et al., 2020a; Yalniz et al., 2019; Zoph et al., 2020), objective detection (Rosenberg et al., 2005), robustness-aware model training against adversarial attacks (Carmon et al., 2019), continual lifelong learning (Lee et al., 2019), and natural language processing (He et al., 2019; Kahn et al., 2020). The terminology “self-training” has been used to describe various SemiSL
algorithms in the literature, while this paper is centered on the commonly-used iterative self-training method in particular. In this setup, an initial teacher model (learned from the labeled data) is applied to the unlabeled data to generate pseudo labels. One then trains a student model by minimizing the weighted empirical risk of both the labeled and unlabeled data. The student model is then used as the new teacher to update the pseudo labels of the unlabeled data. This process is repeated multiple times to improve the eventual student model. We refer readers to Section 2 for algorithmic details.
Despite the empirical achievement of self-training methods with neural networks, the theoretical justification of such success is very limited, even in the field of SemiSL. The majority of the theoretical results on general SemiSL are limited to linear networks (Chen et al., 2020b; Raghunathan et al., 2020; Oymak & Gulcu, 2020; Oneto et al., 2011). The authors in (Balcan & Blum, 2010) show that unlabeled data can improve the generalization bound if the unlabeled data distribution and target model are compatible. For instance, the unlabeled data need to be well-chosen such that the target function for labeled data can separate the unlabeled data clusters, which, however, may not be able to be verified ahead. Moreover, (Rigollet, 2007; Singh et al., 2008) proves that unlabeled data can improve the convergence rate and generalization error under a similar clustering assumption, where the data contains clusters that have homogeneous labels. A recent work by Wei et al. (2020) analyzes SemiSL on nonlinear neural networks and proves that an infinite number of unlabeled data can improve the generalization compared with training with labeled data only. However, Wei et al. (2020) considers single shot rather than iterative SemiSL, and the training problem aims to minimize the consistency regularization rather than the risk function in the conventional self-training method (Lee et al., 2013). Moreover, Wei et al. (2020) directly analyzes the global optimum of the nonconvex training problem without any discussion about how to achieve the global optimum. To the best of our knowledge, there exists no analytical characterization of how the unlabeled data affect the generalization of the learned model by iterative self-training on nonlinear neural networks.
Contributions. This paper provides the first theoretical study of iterative self-training on nonlinear neural networks. Focusing on one-hidden-layer neural networks, this paper provides a quantitative analysis of the generalization performance of iterative self-training as a function of the number of labeled and unlabeled samples. Specifically, our contributions include
1. Quantitative justification of generalization improvement by unlabeled data. Assuming the existence of a ground-truth model with weights W ∗ that maps the features to the corresponding labels, we prove that the learned model via iterative self-training moves closer to W ∗ as the number M of unlabeled data increases, indicating a better testing performance. Specifically, we prove that the Frobenius distance to W ∗, which is approximately linear in the generalization error, decreases in the order of 1/ √ M . As an example, Figure 1 shows that the proposed theoretical bound
matches the empirical self-training performance versus the number of unlabeled data for image classification; see details in Section 4.2.
2. Analytical justification of iterative self-training over single shot alternative. We prove that the student models returned by the iterative self-training method converges linearly to a model close toW ∗, with the rate improvement in the order of 1/ √ M .

heading: 3. Sample complexity analysis of labeled and unlabeled data for learning a proper model.
text: We quantify the impact of labeled and unlabeled data on the generalization of the learned model. In particular, we prove that the sample complexity of labeled data can be reduced compared with supervised learning.

heading: 1.1 RELATED WORKS
text: Semi-supervised learning. Besides self-training, many recent SemiSL algorithms exploit either consistency regularization or entropy minimization. Consistency regularization is based on the as-
sumption that the learned model will return same or similar output when the input is perturbed (Laine & Aila, 2016; Bachman et al., 2014; Sajjadi et al., 2016; Tarvainen & Valpola, 2017; Reed et al., 2015). (Grandvalet & Bengio, 2005) claims that the unlabeled data are more informative if the pseudo labels of the unlabeled data have lower entropy. Therefore, a line of works (Grandvalet & Bengio, 2005; Miyato et al., 2018) adds a regularization term that minimizes the entropy of the outputs of the unlabeled data. In addition, hybrid algorithms that unify both the above regularizations have been developed like (Berthelot et al., 2019a;b; Sohn et al., 2020).
Domain adaptation. Domain adaptation exploits abundant data in the source domain to learn a model for the target domain, where only limited training data are available (Liebelt & Schmid, 2010; Vazquez et al., 2013; Zhang et al., 2013; Long et al., 2015; Tzeng et al., 2014). Source and target domain are related but different. Unsupervised domain adaptation (Ganin & Lempitsky, 2015; Ganin et al., 2016; Gong et al., 2013; Bousmalis et al., 2016), where training data in target domain are unlabeled, is similar to SemiSL, and self-training methods have been used for analysis (Zou et al., 2018; Tang et al., 2012; French et al., 2018). However, self-training and unsupervised domain adaptation are fundamentally different. The former learns a model for the domain where there is limited labeled data, with the help of a large number of unlabeled data from a different domain. The latter learns a model for the domain where the training data are unlabeled, with the help of sufficient labeled data from a different domain.
Generalization analysis of supervised learning. In theory, the testing error is upper bounded by the training error plus the generalization gap between training and testing. These two quantities are often analyzed separately and cannot be proved to be small simultaneously for deep neural networks. For example, neural tangent kernel (NTK) method (Jacot et al., 2018; Du et al., 2018; Lee et al., 2018) shows the training error can be zero, and the Rademacher complexity in (Bartlett & Mendelson, 2002) bounds the generalization gap (Arora et al., 2019a). For one-hidden-layer neural networks (Safran & Shamir, 2018), the testing error can be proved to be zero under mild conditions. One common assumption is that the input data belongs to the Gaussian distribution (Zhong et al., 2017; Ge et al., 2018; Kalai et al., 2008; Bakshi et al., 2019; Zhang et al., 2016; Brutzkus & Globerson, 2017; Li & Yuan, 2017; Soltanolkotabi et al., 2018). Another line of approaches (Brutzkus et al., 2018; Li & Liang, 2018; Wang et al., 2019) consider linearly separable data.
The rest of this paper is organized as follows. Section 2 introduces the problem formulation and self-training algorithm. Major results are summarized in Section 3, and empirical evaluations are presented in Section 4. Section 5 concludes the whole paper. All the proofs are in the Appendix.

heading: 2 FORMALIZING SELF-TRAINING: NOTATION, FORMULATION, AND ALGORITHM
text: Problem formulation. Given N labeled data sampled from distribution Pl, denoted by D = {xn, yn}Nn=1, and M unlabeled data drawn from distribution Pu, denoted by D̃ = {x̃m}Mm=1. The aim is to find a neural network model g(W ), whereW denotes the trainable weights, that minimizes the testing error on data sampled from Pl.
Iterative self-training. In each iteration, given the current teacher predictor g(W (`)), the pseudolabels for the unlabeled data in D̃ are computed as ỹm = g(W (`); x̃m). The method then minimizes the weighted empirical risk f̂D,D̃(W ) of both labeled and unlabeled data through stochastic gradient
descent, where
f̂D,D̃(W ) = λ
2N N∑ n=1 ( yn − g(W ;xn) )2 + λ̃ 2M M∑ m=1 ( ỹm − g(W ; x̃m) )2 , (1)
and λ+ λ̃ = 1. The learned student model g(W (`+1)) is used as the teacher model in the next iteration. The initial model g(W (0)) is learned from labeled data. The formal algorithm is summarized as in Table 1.
Model and assumptions. This paper considers regression1, where g is a one-hidden-layer fully connected neural network equipped with K neurons. Namely, given the input x ∈ Rd and weights W = [w1,w2, · · · ,wK ] ∈ Rd×K , we have
g(W ;x) := 1
K K∑ j=1 φ(wTj x), (2)
where φ is the ReLU activation function2, and φ(z) = max{z, 0} for any input z ∈ R. Here, we fix the top layer weights as 1 for simplicity, and the equivalence of such a simplification is discussed in Appendix K.
Moreover, we assume an unknown ground-truth model with weightsW ∗ that maps all the features to the corresponding labels drawn from Pl, i.e., y = g(W ∗;x), where (x, y) ∼ Pl. The generalization function (GF) with respect to g(W ) is defined as
I ( g(W ) ) = E(x,y)∼Pl ( y − g(W ;x) )2 = E(x,y)∼Pl ( g(W ∗;x)− g(W ;x) )2 . (3)
By definition I ( g(W ∗) ) is zero. Clearly, W ∗ is not unique because any column permutation of W ∗, which corresponds to permuting neurons, represents the same function as W ∗ and minimizes GF in (3) too. To simplify the representation, we follow the convention and abuse the notation that the distance from W to W ∗, denoted by ‖W −W ∗‖F , means the smallest distance from W to any permutation ofW ∗. Additionally, some important notations are summarized in Table 2.
We assume the inputs of both the labeled and unlabeled data belong to the zero mean Gaussian distribution, i.e., x ∼ N (0, δ2Id), and x̃ ∼ N (0, δ̃2Id). The Gaussian assumption is motivated by the data whitening (LeCun et al., 2012) and batch normalization techniques (Ioffe & Szegedy, 2015) that are commonly used in practice to improve learning performance. Moreover, training onehidden-layer neural network with multiple neurons is NP-Complete (Blum & Rivest, 1992) without any assumption.
The focus of this paper. This paper will analyze three aspects about self-training: (1) the generalization performance of W (L), the returned model by self-training after L iterations, measured by ‖W (L) −W ∗‖F 3; (2) the influence of parameter λ in (1) on the training performance; and (3) the impact of unlabeled data on the training and generalization performance.

heading: 3 THEORETICAL RESULTS
text: Beyond supervised learning: Challenge of self-training. The existing theoretical works such as (Zhong et al., 2017; Zhang et al., 2020a;b;c) verify that for one-hidden-layer neural networks, if only labeled data are available, and x are drawn from the standard Gaussian distribution, then supervised learning by minimizing (1) with λ = 1 can return a model with ground-truth weights W ∗ (up to column permutation), as long as the number of labeled data N is at least N∗, which depends on κ,K and d. In contrast, this paper focuses on the low labeled-data regime when N is less than N∗. Specifically,
N∗/4 < N ≤ N∗. (4)
Intuitively, ifN < N∗, the landscape of the empirical risk of the labeled data becomes highly nonconvex, even in a neighborhood of W ∗, thus, the existing analyses for supervised learning do not hold in this region. With additional unlabeled data, the landscape of the weighted empirical risk becomes smoother near W ∗. Moreover, as M increases, and starting from a nearby initialization, the returned modelW (L) by iterative self-training can converge to a local minimum that is closer to W ∗ (see illustration in Figure 2).
Compared with supervised learning, the formal analyses of self-training need to handle new technical challenges from two aspects. First, the existing analyses of supervised learning exploit the fact that the GF and the empirical risk have the same minimizer, i.e., W ∗. This property does not hold for self-training asW ∗ no longer minimizes the weighted empirical risk in (1). Second, the iterative manner of self-training complicates the analyses. Specifically, the empirical risk in each iteration is different and depends on the model trained in the previous iteration through the pseudo labels.
In what follows, we provide theoretical insights and the formal theorems. Some important quantities λ̂ and µ are defined below
λ̂ := λδ2
λδ2 + λ̃δ̃2 , and µ = µ(δ, δ̃) :=
√ λδ2 + λ̃δ̃2
λρ(δ) + λ̃ρ(δ̃) , (5)
where ρ is a positive function defined in (73). λ̂ is an increasing function of λ. Also, from Lemma 11 (in Appendix), ρ(δ) is in the order of δ2 when δ ≤ 1 for ReLU activation functions. Thus, µ is a fixed constant, denoted by µ∗, for all δ, δ̃ ≤ 1. When δ and δ̃ are large, µ increases as they increase. The formal definition of N∗ in (4) is c(κ)µ∗2K3d log q, where c(κ) is some polynomial function of κ and can be viewed as constant.
3.1 INFORMAL KEY THEORETICAL FINDINGS
To the best of our knowledge, Theorems 1 and 2 provide the first theoretical characterization of iterative self-training on nonlinear neural networks. Before formally presenting them, we summarize the highlights as follows.

heading: 1. Linear convergence of the learned models.
text: The learned models converge linearly to a model close to W ∗. Thus, the iterative approach returns a model with better generalization than that by the single-shot method. Moreover, the convergence rate is a constant term plus a term in the order of 1/ √ M (see ∆1 in Figure 3), indicating a faster convergence with more unlabeled data.
2. Returning a model with guaranteed generalization in the low labeled-data regime. Even when the number of labeled data is much less than the required sample complexity to obtain W ∗ in supervised learning, we prove that with the help of unlabeled data, the iterative self-training can return a model in the neighborhood of W [λ̂], where W [λ̂] is in the line segment of W (0) (λ̂ = 0) and ground truth W ∗ (λ̂ = 1). Moreover, λ̂ is upper bounded by √ N/N∗. Thus W (L) moves closer to W ∗ as N increases (E0 in Figure 3), indicating a better generalization performance with more labeled data.
3. Guaranteed generalization improvement by unlabeled data. The distance between W (L) and W [λ̂] (E1 in Figure 3) scales in the order of 1/ √ M . With a larger number of unlabeled data M , W (L) moves closer to W [λ̂] and thus W ∗, indicating an improved generalization performance (Theorem 1). When N is close to N∗ but still smaller as defined in (12), both W (L) and W [λ̂] converge toW ∗, and thus the learned model achieves zero generalization error (Theorem 2).

heading: 3.2 FORMAL THEORY IN LOW LABELED-DATA REGIME
text: Takeaways of Theorem 1: Theorem 1 characterizes the convergence rate of the proposed algorithm and the accuracy of the learned model W (L) in a low labeled-data regime. Specifically, the iterates converge linearly, and the learned model is close to W [λ̂] and guaranteed to outperform the initial modelW (0). Theorem 1. Suppose the initializationW (0) and the number of labeled data satisfy
‖W (0) −W ∗‖F ≤ p−1 · ‖W ∗‖F
c(κ)µ2K3/2 with p ∈ (1 2 , 1 ] , (6)
and max { 1 K , p− 2p− 1 µ √ K }2 ·N∗ ≤ N ≤ N∗. (7)
If the value of λ̂ in (5) and unlabeled data amount M satisfy max { 1 K , p− 2p− 1 µ √ K } ≤ λ̂ ≤ min {√ N N∗ , p+ 2p− 1 µ √ K } , (8)
and M ≥ (2p− 1)−2c(κ)µ2 ( 1− λ̂ )2 K3d log q. (9)
Then, when the number T of SGD iterations is large enough in each loop `, with probability at least 1− q−d, the iterates {W (`)}L`=0 converge toW [λ̂] as
‖W (L) −W [λ̂]‖F ≤ (( 1 + Θ (µ(1−λ̂)√
M
)) · λ̂ )L · ‖W (0) −W [λ̂]‖2 + ( 1 + Θ (µ(1−λ̂)√ M )) · ‖W ∗ −W [λ̂]‖F ,
(10) whereW [λ̂] = λ̂W ∗ + (1− λ̂)W (0). Typically, when the iteration number L is sufficient large, we have
‖W (L) −W ∗‖F ≤ ( 1 + Θ (µ(1− λ̂)√
M
)) · 2(1− λ̂) · ‖W ∗ −W (0)‖F . (11)
The accuracy of the learned model W (L) with respect to W ∗ is characterized as (10), and the learning model is better than initial model as in (11) if the following conditions hold. First, the weights λ in (1) are properly chosen as in (8). Second, the number of unlabeled data is sufficiently large as in (9).
Selection of λ in self-training algorithms. When λ̂ increases, the required number of unlabeled data is reduced from (9), and the convergence point W (L) becomes closer toW ∗ from (11), which indicates a smaller generalization error. Thus, a large λ̂ within its feasible range (8) is desirable. When the initial model W (0) is closer to W ∗ (corresponding to a larger p), and the number of labeled data N increases, the upper bound in (8) increases, and thus, one can select a larger λ̂.
The initial model W (0). The tensor initialization from (Zhong et al., 2017) can return a W (0) that satisfies (6) when the number of labeled data isN = p2N∗ (see Lemma 3 in Appendix). Combining with the requirement in (7), Theorem 1 applies to the case that N is at least N∗/4.

heading: 3.3 FORMAL THEORY OF ACHIEVING ZERO GENERALIZATION ERROR
text: Takeaways of Theorem 2: Theorem 2 indicates the model returned by the proposed algorithm converges linearly to the ground truthW ∗. Thus the distance between the learned model and the ground truth can be arbitrarily small with the ability to achieve zero generalization error. The required sample complexity is reduced by a constant factor compared with supervised learning. Theorem 2. Consider the number of unlabeled data satisfies(
1− 1/(µ √ K) )2 ·N∗ ≤ N ≤ N∗, (12)
we choose λ̂ such that 1− 1/(µ √ K) ≤ λ̂ ≤ √ N/N∗. (13)
Suppose the initial modelW (0) and the number of unlabeled data M satisfy
‖W (0) −W ∗‖F ≤ ‖W ∗‖F
c(κ)µ2K3/2 and M ≥ c(κ)µ2(1− λ̂)2K3d log q, (14)
the iterates {W (`)}L`=0 converge to the ground truthW ∗ as follows,
‖W (L) −W ∗‖F ≤ [(
1 + c(κ)λ̂√ N + c(κ)(1− λ̂)√ M
) · µ √ K(1− λ̂) ]L · ‖W (0) −W ∗‖F . (15)
The models W (`)’s converge linearly to the ground truth W ∗ as (15) when the number of labeled data satisfies (12). In contrast, supervised learning requires at least N∗ labeled samples to estimate W ∗ accurately without unlabeled data, which suggests self-training at least saves a constant fraction of labeled data.

heading: 3.4 THE MAIN PROOF IDEA
text: Our proof builds upon and extends one recent line of works on supervised learning such as (Zhong et al., 2017; Zhang et al., 2020b;c; 2021). The standard framework of these works is first to show that the generalization function I(g(W )) in (3) is locally convex near W ∗, which is its global minimizer. Then, when M = 0 and N is sufficiently large, the empirical risk function using labeled data only can approximate I(g(W )) well in the neighborhood of W ∗. Thus, if initialized in this local convex region, the iterations, returned by applying gradient descent approach on the empirical risk function, converge toW ∗ linearly.
The technical challenge here is that in self-training, when unlabeled data are paired with pseudo labels, W ∗ is no longer a global minimizer of the empirical risk f̂D,D̃ in (1), and f̂D,D̃ does not approach I(g(W )) even when M and N increase to infinity. Our new idea is to design a population risk function f(W ; λ̂) in (17) (see appendix), which is a lower bound of f̂D,D̃ when M and N are infinity. f(W ; λ̂) is locally convex around its minimizer W [λ̂], and W [λ̂] approaches W ∗ as λ̂ increases. Then we show the iterates generated by f̂D,D̃ stay close to f(W ; λ̂), and the returned
modelW (L) is close toW [λ̂]. New technical tools are developed to bound the distance between the functions f̂D,D̃ and f(W ; λ̂).

heading: 4 EMPIRICAL RESULTS

heading: 4.1 SYNTHETIC DATA EXPERIMENTS
text: We generate a ground-truth neural network with the width K = 10. Each entry ofW ∗ is uniformly selected from [−2.5, 2.5]. The input of labeled data xn are generated from Gaussian distribution N (0, Id) independently, and the corresponding label yn is generated through (2) using W ∗. The unlabeled data x̃m are generated from N (0, δ̃2Id) independently with δ̃ = 1 except in Figure 7. d is set as 50 except in Figure 9. The value of λ is selected as √ N/(2Kd) except in Figure 8.
We consider one-hidden-layer except in Figure 4. The initial teacher model W (0) in self-training is randomly selected from {W |‖W −W ∗‖F /‖W ∗‖F ≤ 0.5} to reduce the computation. In
each iteration, the maximum number of SGD steps T is 10. Self-training terminates if ‖W (`+1) − W (`)‖F /‖W (`)‖F ≤ 10−4 or reaching 1000 iterations. In Figures 5 to 8, all the points on the curves are averaged over 1000 independent trials, and the regions in lower transparency indicate the corresponding one-standard-deviation error bars. Our empirical observations are summarized below.
(a) GF (testing performance) proportional to ‖W −W ∗‖F . Figure 4 illustrates the GF in (3) against the distance to the ground truth W ∗. To visualize results for different networks together, GF is normalized in [0, 1], divided by its largest value for each network architecture. All the results are averaged over 100 independent choice of W . One can see that for one-hidden-layer neural networks, in a large region near W ∗, GF is almost linear in ‖W −W ∗‖F . When the number of hidden layers increases, this region decreases, but the linear dependence still holds locally. This is an empirical justification of using ‖W −W ∗‖F to evaluate the GF and, thus, the testing error in Theorems 1 and 2. (b) ‖W (L) −W ∗‖F as a linear function of 1/ √ M . Figure 5 shows the relative error ‖W (L) − W ∗‖F /‖W ∗‖F when the number of unlabeled data and labeled data changes. One can see that the relative error decreases when either M or N increases. Additionally, the dash-dotted lines represent the best fitting of the linear functions of 1/ √ M using the least square method. Therefore, the relative error is indeed a linear function of 1/ √ M , as predicted by our results in (11) and (15).
(c) Convergence rate as a linear function of 1/ √ M . Figure 6 illustrates the convergence rate when M and N change. We can see that the convergence rate is a linear function of 1/ √ M , as predicted by our results (11) and (15). When M increases, the convergence rate is improved, and the method converges faster.
(d) Increase of δ̃ slows down convergence. Figure 7 shows that the convergence rate becomes worse when the variance of the unlabeled data δ̃ increases from 1. When δ̃ is less than 1, the convergence rate almost remains the same, which is consistent with our characterization in (10) that the convergence rate is linear in µ. From the discussion after (5), µ increases as δ̃ increases from 1 and stays constant when δ̃ is less than 1.
(e) ‖W (L)−W ∗‖F /‖W ∗‖F is improved as a linear function of λ̂. Figure 8 shows that the relative errors of W (L) with respect to W ∗ decrease almost linearly when λ̂ increases, which is consistent with the theoretical result in (11). Moreover, when λ exceeds a certain threshold positively correlated with N , the relative error increases rather than decreases. That is consistent with the analysis in (8) that λ̂ has an upper limit, and such a limit increases as N increases.
(f) Unlabeled data reduce the sample complexity to learn W ∗. Figure 9 depicts the phase transition of returning W (L). For every pair of d and N , we construct 100 independent trials, and each trial is said to be successful if ‖W (L) −W ∗‖F /‖W ∗‖F ≤ 10−2. The white blocks correspond to the successful trials, while the block in black indicates all failures. When d increases, the required number of labeled data to learn W ∗ is linear in d. Thus, the sample complexity bound in (12) is order-wise optimal for d. Moreover, the phase transition line when M = 1000 is below the one when M = 0. Therefore, with unlabeled data, the required sample complexity of N is reduced.
1 2 3 4 5
0.99
0.995
1
C on
ve rg
en ce
r at
e
Figure 7: Convergence rate with different δ̂.
when λ̂ and N change.

heading: 4.2 IMAGE CLASSIFICATION ON AUGMENTED CIFAR-10 DATASET
text: We evaluate self-training on the augmented CIFAR-10 dataset, which has 50K labeled data. The unlabeled data are mined from 80 Million Tiny Images following the setup in (Carmon et al., 2019)4, and additional 50K images are selected for each class, which is a total of 500K images, to form the unlabeled data. The self-training method is the same implementation as that in (Carmon et al., 2019). λ and λ̃ is selected as N/(M +N) and M/(N +M), respectively, and the algorithm stops after 200 epochs. In Figure 10, the dash lines stand for the best fitting of the linear functions of 1/ √ M via the least square method. One can see that the test accuracy is improved by up to 7% using unlabeled data, and the empirical evaluations match the theoretical predictions. Figure 11 shows the convergence rate calculated based on the first 50 epochs, and the convergence rate is almost a linear function of 1/ √ M , as predicted by (10).

heading: 5 CONCLUSION
text: This paper provides new theoretical insights into understanding the influence of unlabeled data in the iterative self-training algorithm. We show that the improved generalization error and convergence rate is a linear function of 1/ √ M , where M is the number of unlabeled data. Moreover, compared with supervised learning, using unlabeled data reduces the required sample complexity of labeled data for achieving zero generalization error. Future directions include generalizing the analysis to multi-layer neural networks and other semi-supervised learning problems such as domain adaptation.

heading: ACKNOWLEDGEMENT
text: This work was supported by AFOSR FA9550-20-1-0122, ARO W911NF-21-1-0255, NSF 1932196 and the Rensselaer-IBM AI Research Collaboration (http://airc.rpi.edu), part of the IBM AI Horizons Network (http://ibm.biz/AIHorizons).
4The codes are downloaded from https://github.com/yaircarmon/semisup-adv

heading: A OVERVIEW OF THE PROOF TECHNIQUES
text: We first provide an overview of the techniques used in proving Theorems 1 and 2.
1. Characterization of a proper population risk function. To characterize the performance of the iterative self-training algorithm via the stochastic gradient descent method, we need first to define a population risk function such that the following two properties hold. First, the landscape of the population risk function should be analyzable near {W (`)}L`=0. Second, the distance between the empirical risk function in (1) and the population risk function should be bounded near {W (`)}L`=0. The generalization function defined in (3), which is widely used in the supervised learning problem with a sufficient number of samples, failed the second requirement. To this end, we turn to find a new population risk function defined in (17), and the illustrations of the population risk function and objection function are included in Figure 12.
2. Local convex region of the population risk function. The purpose is to characterize the iterations via the stochastic gradient descent method in the population risk function. To obtain the local convex region of the population risk function, we first bound the Hessian of the population risk function at its global optimal. Then, we utilize Lemma 12 in Appendix H.1 to obtain the Hessian of the population risk function near the global optimal. The local convex region of the population risk function is summarized in Lemma 1, and the proof of Lemma 1 is included in Appendix H.1.
3. Bound between the population risk and empirical risk functions. After the characterization of the iterations via the stochastic gradient descent method in the population risk function, we need to bound the distance between the population risk function and empirical risk function. Therefore, the behaviors of the iterations via the stochastic gradient descent method in the empirical risk function can be described by the ones in the population risk function and the distance between these two. The key lemma is summarized in Lemma 2 (see Appendix H.2), and the proof is included in Appendix H.2.
In the following contexts, the details of the iterative self-training algorithm are included in Appendix B. We then first provide the proof of Theorem 2 in Appendix E, which can be viewed as a special case of Theorem 1. Then, with the preliminary knowledge from proving Theorem 2, we turn to present the full proof of a more general statement summarized in Theorem 3 (see Appendix F), which is related to Theorem 1. The definition and relative proofs of µ and ρ are all included in Appendix G. The proofs of preliminary lemmas are included in Appendix H.
B ITERATIVE SELF-TRAINING ALGORITHM
In this section, we implement the details of the mini-batch stochastic gradient descent used in each stage of the iterative self-training algorithm. After t number of iterations via mini-batch stochastic
gradient descent at `-th stage of self-training algorithm, the learned model is denoted asW (`,t). One can easily check thatW (`) in the main context is denoted asW (`,0) in this section and the following proofs. Last, the pseudo-code of the iterative self-training algorithm is summarized in Algorithm 1.
Algorithm 1 Iterative Self-Training Algorithm
Input: labeled D = {(xn, yn)}Nn=1, unlabeled data D̃ = {x̃m}Mm=1, and gradient step size η;
Initialization: preliminary teacher model with weightsW (0,0);
Partition: randomly and independently pick data from D and D̃ to form T subsets {Dt}T−1t=0 and {D̃t}T−1t=0 , respectively;
for ` = 0, 1, · · · , L− 1 do
ym = g(W (`,0); x̃m) for m = 1, 2, · · · ,M
for t = 0, 1, · · · , T − 1 do
W (`,t+1) = W (`,t) − η · ∇f̂Dt,D̃t(W (`,t)) + β ·
( W (`,t) −W (`,t−1) ) end for
W (`+1,0) = W (`,T )
end for

heading: C NOTATIONS
text: In this section, we first introduce some important notations that will be used in the following proofs, and the notations are summarized in Table 1.
As shown in Algorithm 1, W (`,t) denotes the learned model after t number of iterations via minibatch stochastic gradient descent at `-th stage of the iterative self-training algorithm. Given a student model W̃ , the pseudo label for x̃ ∈ D̃ is generated as
ỹ = g(W̃ ; x̃). (16)
Further, letW [p] = pW ∗ + (1− p)W (0,0), we then define the population risk function as
f(W ; p) = λ 2 Ex ( y∗(p)− g(W ;x) )2 + λ̃ 2 Ex̃ ( ỹ∗(p)− g(W ; x̃) )2 , (17)
where y∗(p) = g(W [p];x) with x ∼ N (0, δ2I) and ỹ∗(p) = g(W [p]; x̃) with x̃ ∼ N (0, δ̃2I). When p = 1, we haveW [p] = W ∗ and y∗(p) = y for data in D. Moreover, we use σi to denote the i-th largest singular value of W ∗. Then, κ is defined as σ1/σK , and γ = ∏K i=1 σi/σK . Additionally, to avoid high dimensional tensors, the first order derivative of the empirical risk function is defined in the form of vectorizedW as
∇f̂(W ) = [ ∂f ∂w1 T , ∂f ∂w2 T , · · · , ∂f ∂wK T ]T ∈ RdK (18)
with W = [w1,w2, · · · ,wK ] ∈ Rd×K . Therefore, the second order derivative of the empirical risk function is in Rdk×dk. Similar to (18), the high order derivatives of the population risk functions are defined based on vectorized W as well. In addition, without special descriptions, α = [αT1 ,α T 2 , · · · ,αTK ]T stands for any unit vector that in RdK with αj ∈ Rd. Therefore, we have
‖∇2f̂‖2 = max α ‖αT∇2f̂α‖2 = max α ( K∑ j=1 αTj ∂f̂ ∂wj )2 . (19)
Finally, since we focus on order-wise analysis, some constant numbers will be ignored in the majority of the steps. In particular, we use h1(z) & (or .,h)h2(z) to denote there exists some positive constant C such that h1(z) ≥ (or ≤,=)C · h2(z) when z ∈ R is sufficiently large.

heading: D PRELIMINARY LEMMAS
text: We will first start with some preliminary lemmas. As outlined at the beginning of the supplementary material, Lemma 1 illustrates the local convex region of the population risk function, and Lemma 2 explains the error bound between the population risk and empirical risk functions. Then, Lemma 3 describes the returned initial model W (0,0) via tensor initialization method (Zhong et al., 2017) purely using labeled data. Next, Lemma 4 is the well known Weyl’s inequality in the matrix setting. Moreover, Lemma 5 is the concentration theorem for independent random matrices. The definitions of the sub-Gaussian and sub-exponential variables are summarized in Definitions 1 and 2. Lemmas 6 and 7 serve as the technical tools in bounding matrix norms under the framework of the confidence interval.
Lemma 1. Given anyW ∈ Rd×K , let p satisfy
p . σK
µ2K · ‖W −W ∗‖F . (20)
Then, we have
λρ(δ) + λ̃ρ(δ̃) 12κ2γK2 ∇2f(W ; p) 7(λδ
2 + λ̃δ̃2)
K . (21)
Lemma 2. Let f and f̂ be the functions defined in (17) and (1), respectively. Suppose the pseudo label is generated through (16) with weights W̃ . Then, we have
‖∇f(W )−∇f̂(W )‖2 . λδ2
K
√ d log q
N · ‖W −W ∗‖+ λ̃δ̃
2
K
√ d log q
M · ‖W − W̃ ‖2
+
∥∥λδ2 · (W̃ −W [p])+ λ̃δ̃2 · (W ∗ −W [p])∥∥ 2
2K
(22)
with probability at least 1− q−d. Lemma 3 (Initialization, (Zhong et al., 2017)). Assuming the number of labeled data satisfies
N ≥ p2N∗ (23)
for some large constant q and p ∈ [ 1K , 1], the tensor initialization method, which is summarized in Appendix I, outputsW (0,0) such that
‖W (0,0) −W ∗‖F ≤ σK
p · c(κ)µ2K (24)
with probability at least 1− q−d. Lemma 4 (Weyl’s inequality, (Bhatia, 2013)). LetB = A+E be a matrix with dimension m×m. Let λi(B) and λi(A) be the i-th largest eigenvalues ofB andA, respectively. Then, we have
|λi(B)− λi(A)| ≤ ‖E‖2, ∀ i ∈ [m]. (25)
Lemma 5 ((Tropp, 2012), Theorem 1.6). Consider a finite sequence {Zk} of independent, random matrices with dimensions d1 × d2. Assume that such random matrix satisfies
E(Zk) = 0 and ‖Zk‖ ≤ R almost surely. Define
δ2 := max {∥∥∥∑
k
E(ZkZ∗k) ∥∥∥,∥∥∥∑
k
E(Z∗kZk) ∥∥∥}.
Then for all t ≥ 0, we have
Prob {∥∥∥∥∥∑ k Zk ∥∥∥∥∥ ≥ t } ≤ (d1 + d2) exp ( −t2/2 δ2 +Rt/3 ) .
Definition 1 (Definition 5.7, (Vershynin, 2010)). A random variable X is called a sub-Gaussian random variable if it satisfies
(E|X|p)1/p ≤ c1 √ p (26)
for all p ≥ 1 and some constant c1 > 0. In addition, we have
Ees(X−EX) ≤ ec2‖X‖ 2 ψ2 s2 (27)
for all s ∈ R and some constant c2 > 0, where ‖X‖φ2 is the sub-Gaussian norm of X defined as ‖X‖ψ2 = supp≥1 p−1/2(E|X|p)1/p.
Moreover, a random vector X ∈ Rd belongs to the sub-Gaussian distribution if one-dimensional marginal αTX is sub-Gaussian for any α ∈ Rd, and the sub-Gaussian norm of X is defined as ‖X‖ψ2 = sup‖α‖2=1 ‖α
TX‖ψ2 . Definition 2 (Definition 5.13, (Vershynin, 2010)). A random variableX is called a sub-exponential random variable if it satisfies (E|X|p)1/p ≤ c3p (28) for all p ≥ 1 and some constant c3 > 0. In addition, we have
Ees(X−EX) ≤ ec4‖X‖ 2 ψ1 s2 (29)
for s ≤ 1/‖X‖ψ1 and some constant c4 > 0, where ‖X‖ψ1 is the sub-exponential norm of X defined as ‖X‖ψ1 = supp≥1 p−1(E|X|p)1/p.
Lemma 6 (Lemma 5.2, (Vershynin, 2010)). Let B(0, 1) ∈ {α ∣∣‖α‖2 = 1,α ∈ Rd} denote a unit ball in Rd. Then, a subset Sξ is called a ξ-net of B(0, 1) if every point z ∈ B(0, 1) can be approximated to within ξ by some point α ∈ B(0, 1), i.e., ‖z − α‖2 ≤ ξ. Then the minimal cardinality of a ξ-net Sξ satisfies |Sξ| ≤ (1 + 2/ξ)d. (30) Lemma 7 (Lemma 5.3, (Vershynin, 2010)). Let A be an d1 × d2 matrix, and let Sξ(d) be a ξ-net of B(0, 1) in Rd for some ξ ∈ (0, 1). Then
‖A‖2 ≤ (1− ξ)−1 max α1∈Sξ(d1),α2∈Sξ(d2) |αT1Aα2|. (31)
Lemma 8 (Mean Value Theorem). Let U ⊂ Rn1 be open and f : U −→ Rn2 be continuously differentiable, and x ∈ U , h ∈ Rn1 vectors such that the line segment x + th, 0 ≤ t ≤ 1 remains in U . Then we have:
f(x+ h)− f(x) = (∫ 1
0
∇f(x+ th)dt ) · h,
where∇f denotes the Jacobian matrix of f .

heading: E PROOF OF THEOREM 2
text: With p = 1 in (17), the population risk function is reduced as
f(W ) = λ
2 Ex(y − g(W ;x)) +
λ̃ 2 Ex̃(ỹ∗ − g(W ; x̃)), (32)
where y = g(W ∗;x) with x ∼ N (0, δ2I) and ỹ∗ = g(W ∗; x̃) with x̃ ∼ N (0, δ̃2I). In fact, (32) can be viewed as the expectation of the empirical risk function in (1) given ỹm = g(W ∗; x̃m). Moreover, the ground-truth model W ∗ is the global optimal to (32) as well. Lemmas 9 and 10 are the special case of Lemmas 1 and 2 with p = 1. The proof of Theorem 2 is followed by the presentation of the two lemmas.
The main idea in proving Theorem 2 is to characterize the gradient descent term by the MVT in Lemma 8 as shown in (36) and (37). The IVT is not directly applied in the empirical risk function because of its non-smoothness. However, the population risk functions defined in (17) and (32), which are the expectations over the Gaussian variables, are smooth. Then, as the distance ‖∇f(W )−∇f(W ∗)‖F is upper bounded by a linear function of ‖W −W ∗‖F as shown in (47), we can establish the connection between ‖W (`,t+1) −W ∗‖F and ‖W (`,t) −W ∗‖F as shown in (50). Finally, by mathematical induction over ` and t, one can characterize ‖W (L,0) −W ∗‖F by ‖W (0,0) −W ∗‖F as shown in (52), which completes the whole proof. Lemma 9 (Lemma 1 with p = 1). Let f and f̂ are the functions defined in (32) and (1), respectively. Then, for anyW that satisfies,
‖W −W ∗‖F ≤ σK µ2K , (33)
we have λρ(δ) + λ̃ρ(δ̃)
12κ2γK2 ∇2f(W ) 7(λδ
2 + λ̃δ̃2)
K . (34)
Lemma 10 (Lemma 2 with p = 1). Let f and f̂ be the functions defined in (32) and (1), respectively. Suppose the pseudo label is generated through (16) with weights W̃ . Then, we have
‖∇f(W )−∇f̂(W )‖2 . (λδ2 K √ d log q N + (1− λ)δ̃2 K √ d log q M ) · ‖W −W ∗‖2
+ (1− λ)δ̃2
K
(√d log q M + 1 2 ) · ‖W̃ −W ∗‖2
(35)
with probability at least 1− q−d.
Proof of Theorem 2. From Algorithm 1, in the `-th outer loop, we have
W (`,t+1) =W (`,t) − η∇f̂Dt,D̃t(W (`,t)) + β(W (`,t) −W (`,t−1))
=W (`,t) − η∇f(W (`,t)) + β(W (`,t) −W (`,t−1)) + η · ( ∇f(W (`,t))−∇f̂Dt,D̃t(W (`,t)) ) .
(36)
Since∇f is a smooth function andW ∗ is a local (global) optimal to f , then we have
∇f(W (`,t)) =∇f(W (`,t))−∇f(W ∗)
= ∫ 1 0 ∇2f ( W (`,t) + u · (W (`,t) −W ∗) ) du · (W (`,t) −W ∗),
(37)
where the last equality comes from MVT in Lemma 8. For notational convenience, we use H(`,t) to denote the integration as
H(`,t) := ∫ 1 0 ∇2f ( W (`,t) + u · (W (`,t) −W ∗) ) du. (38)
Then, we have [ W (`,t+1) −W ∗
W (`,t) −W ∗
] = [ I − ηH(`,t) βI
I 0
][ W (`,t) −W ∗
W (`,t−1) −W ∗
]
+ η
[ ∇f(W (`,t))−∇f̂Dt,D̃t(W (`,t))
0
] .
(39)
LetH(`,t) = SΛST be the eigen-decomposition ofH(`,t). Then, we define
A(β) :=
[ ST 0
0 ST
] A(β) [ S 0
0 S
] = [ I − ηΛ + βI βI
I 0
] . (40)
Since
[ S 0
0 S
][ ST 0
0 ST
] = [ I 0
0 I
] , we know A(β) and [ I − ηΛ + βI βI
I 0
] share the same
eigenvalues. Let γ(Λ)i be the i-th eigenvalue of ∇2f(ŵ(t)), then the corresponding i-th eigenvalue of (40), denoted by γ(A)i , satisfies
(γ (A) i (β)) 2 − (1− ηγ(Λ)i + β)γ (A) i (β) + β = 0. (41)
By simple calculation, we have
|γ(A)i (β)| =  √ β, if β ≥ ( 1− √ ηγ (Λ) i )2 ,
1 2 ∣∣∣∣(1− ηγ(Λ)i + β) +√(1− ηγ(Λ)i + β)2 − 4β∣∣∣∣ , otherwise. (42) Specifically, we have
γ (A) i (0) > γ (A) i (β), for ∀β ∈ ( 0, (1− ηγ(Λ)i ) 2 ) , (43)
and γ(A)i achieves the minimum γ (A)∗ i = ∣∣∣1−√ηγ(Λ)i ∣∣∣ when β = (1−√ηγ(Λ)i )2. From Lemma 9, for any a ∈ Rd with ‖a‖2 = 1, we have
aT∇f(W (`,t))a = ∫ 1 0 aT∇2f ( W (`,t) + u · (W (`,t) −W ∗) ) adu ≤ ∫ 1 0 γmax‖a‖22du = γmax,
aT∇f(W (`,t))a = ∫ 1 0 aT∇2f ( W (`,t) + u · (W (`,t) −W ∗) ) adu ≥ ∫ 1 0
γmin‖a‖22du = γmin, (44)
where γmax = 7(λδ2+λ̃δ̃2) K , and γmin = λρ(δ)+λ̃ρ(δ̃) 12κ2γK2 . Therefore, we have
γ (Λ) min =
λρ(δ) + λ̃ρ(δ̃)
12κ2γK2 , and γ(Λ)max =
7(λδ2 + λ̃δ̃2)
K . (45)
Thus, we can select η = (
1√ γ (Λ) max+ √ γ (Λ) min
)2 , and ‖A(β)‖2 can be bounded by
min β ‖A(β)‖2 ≤1− √(λρ(δ) + λ̃ρ(δ̃) 12κ2γK2 ) / ( 2 · 7(λδ 2 + λ̃δ̃2) K
) =1− µ(δ, δ̃)√
168κ2γK ,
(46)
where µ(δ, δ̃) = ( λρ(δ)+λ̃ρ(δ̃)
λδ2+λ̃δ̃2
)1/2 .
From Lemma 10, we have
‖∇f(W (`,t))−∇f̂(W (`,t))‖2 = (λδ2 K √ d log q Nt + λ̃δ̃2 K √ d log q Mt ) · ‖W (`,t) −W ∗‖2
+ λ̃δ̃2
K (√d log q Mt + 1 2 ) · ‖W (`,0) −W ∗‖2.
(47)
Given ε > 0 and ε̃ > 0 with ε+ ε̃ < 1, let
η · λδ 2
K
√ d log q
Nt ≤ εµ(δ, δ̃)√ 168κ2γK ,
and η · λ̃δ̃ 2
K
√ d log q
Mt ≤ ε̃µ(δ, δ̃)√ 168κ2γK ,
(48)
where we need
Nt ≥ ε−2µ−2 ( λδ2 λδ2 + λ̃δ̃2 )2 κ2γK3d log q,
and Mt ≥ ε̃−2µ−2 ( λ̃δ̃2 λδ2 + λ̃δ̃2 )2 κ2γK3d log q.
(49)
Therefore, from (46), (47) and (48), we have
‖W (`,t+1) −W ∗‖2
≤ (
1− (1− ε− ε̃)µ(δ, δ̃)√ 168κ2γK
) ‖W (`,t) −W ∗‖2 + η · λ̃δ̃2
K (√d log q Mt + 1 2 ) · ‖W (`,0) −W ∗‖2
≤ (
1− (1− ε− ε̃)µ(δ, δ̃)√ 168κ2γK
) ‖W (`,t) −W ∗‖2 + η · λ̃δ̃2
K ‖W (`,0) −W ∗‖2
(50)
when M ≥ 4d log q. By mathematical induction on (50) over t, we have
‖W (`,t) −W ∗‖2 ≤ (
1− (1− ε− ε̃)µ√ 168κ2γK
)t · ‖W (`,0) −W ∗‖2
+
√ 168κ2γK (1− ε− ε̃)µ ·
√ K
14(λδ2 + λ̃δ̃2) · λ̃δ̃
2
K ‖W (`,0) −W ∗‖2
≤ [(
1− (1− ε− ε̃)µ√ 168κ2γK
)t +
√ κ2γλ̃δ̃2
(1− ε− ε̃)µ(λδ2 + λ̃δ̃2)
] · ‖W (`,0) −W ∗‖2
(51)
By mathematical induction on (51) over `, we have
‖W (`,T ) −W ∗‖2 ≤ [(
1− (1− ε− ε̃)µ√ 168κ2γK
)T +
√ κ2γλ̃δ̃2
(1− ε− ε̃)µ(λδ2 + λ̃δ̃2)
]` · ‖W (0,0) −W ∗‖2
(52)

heading: F PROOF OF THEOREM 1
text: Instead of proving Theorem 1, we turn to prove a stronger version, as shown in Theorem 3. One can verify that Theorem 1 is a special case of Theorem 3 by selecting λ̂ in the order of p and ε̃ is in the order of (2p− 1). The major idea in proving Theorem 3 is similar to that of Theorem 2. The first step is to characterize the gradient descent term on the population risk function by the MVT in Lemma 8 as shown in (58) and (59). Then, the connection between ‖W (`+1,0) −W [p]‖F and ‖W (`,0) −W [p]‖F are characterized in (64). Compared with proving Theorem 2, where the induction over ` holds naturally with large size of labeled data, the induction over ` requires a proper value of p as shown in (69). By induction over ` on (64), the relative error ‖W (L,0) −W [p]‖F can be characterized by ‖W (0,0) − W [p]‖F as shown in (71). Theorem 3. Suppose the initializationW (0,0) satisfies with
|p− λ̂| ≤ 2(1− ε̃)p− 1 µ √ K
(53)
for some constant ε̃ ∈ (0, 1/2), where
λ̂ := λδ2 λδ2 + λ̃δ̃2 = ( N κ2γK3µ2d log q ) 1 2 (54)
and
µ = µ(δ, δ̃) = λδ2 + λ̃δ̃2
λρ(δ) + λ̃ρ(δ̃) . (55)
Then, if the number of samples in D̃ further satisfies M & ε̃−2κ2γµ2 ( 1− λ̂ )2 K3d log q, (56)
the iterates {W (`,t)}L,T`,t=0 converge toW [p] with p satisfies (53) as
lim T→∞
‖W (`,T ) −W [p]‖2
≤ 1 1− ε̃
· ( 1− p∗ + µ √ K ∣∣(λ̂− p∗)∣∣) · ‖W (0,0) −W ∗‖2 + ε̃
(1− ε̃) · ‖W (`,0) −W [p]‖2,
(57)
with probability at least 1− q−d.
Proof of Theorem 3. From Algorithm 1, in the `-th outer loop, we have
W (`,t+1) =W (`,t) − η∇f̂Dt,D̃t(W (`,t)) + β(W (`,t) −W (`,t−1))
=W (`,t) − η∇f(W (`,t)) + β(W (`,t) −W (`,t−1)) + η · ( ∇f(W (`,t))−∇f̂Dt,D̃t(W (`,t)) ) (58)
Since∇f is a smooth function andW [p] is a local (global) optimal to f , then we have
∇f(W (`,t)) =∇f(W (`,t))−∇f(W [p])
= ∫ 1 0 ∇2f ( W (`,t) + u · (W (`,t) −W [p]) ) du · (W (`,t) −W [p]),
(59)
where the last equality comes from Lemma 8.
Similar to the proof of Theorem 2, we have
‖W (`,t+1)−W [p]‖2 ≤ ‖A(β)‖2·‖W (`,t)−W [p]‖2+η·‖∇f(W (`,t))−∇f̂Dt,D̃t(W (`,t))‖2. (60)
From Lemma 2, we have
‖∇f(W (`,t))−∇f̂(W (`,t))‖2
. λδ2
K
√ d log q
Nt · ‖W (`,t) −W ∗‖+ λ̃δ̃
2
K
√ d log q
Mt · ‖W (`,t) −W (`,0)‖2
+ ∣∣λδ2 · (W (0,0) −W [p])− λ̃δ̃2 · (W ∗ −W [p])∣∣ K
(61)
When ` = 0, following the similar steps from (41) to (46), we have
‖∇f(W (`,t))−∇f̂(W (`,t))‖2
. λδ2
K
√ d log q
Nt · ‖W (`,t) −W [p]‖+ λ̃δ̃
2
K
√ d log q
Mt · ‖W (`,t) −W [p]‖2
+ λδ2
K
√ d log q
Nt · ‖W ∗ −W [p]‖+ λ̃δ̃
2
K
√ d log q
Mt · ‖W (0,0) −W [p]‖2
+ ∣∣λδ2 · (1− p)− λ̃δ̃2 · p∣∣ K · ‖W (0,0) −W ∗‖2
(62)
and
‖W (`,t+1) −W [p]‖2 ≤ (
1− 1− ε̃ µ(δ, δ̃) √ 154κ2γK
) · ‖W (`,t) −W [p]‖2
+ η · (λδ2(1− p)
K
√ d log q
Nt + ∣∣λδ2 · (1− p)− λ̃δ̃2 · p∣∣ K ) · ‖W (0,0) −W ∗‖2
+ η · ε̃λ̃δ̃ 2 · p K
· √ d log q
Mt ‖W (0,0) −W ∗‖2.
(63)
Therefore, we have
lim T→∞
‖W (`,T ) −W [p]‖2
≤µ √ 154κ2γK
1− ε̃ · η · [(λδ2(1− p) K √ d log q Nt + ∣∣λδ2 · (1− p)− λ̃δ̃2 · p∣∣ K ) · ‖W (0,0) −W ∗‖2
+ ε̃λ̃δ̃2 · p K
· √ d log q
Mt · ‖W (0,0) −W ∗‖2 ] ≤µ √ 154κ2γK
1− ε̃ · K 14(λδ2 + λ̃δ̃2) · [(λδ2(1− p) K √ d log q Nt + ∣∣λδ2 · (1− p)− λ̃δ̃2 · p∣∣ K ) · ‖W (0,0) −W ∗‖2 +
ε̃λ̃δ̃2 · p K
· √ d log q
Mt · ‖W (0,0) −W ∗‖2 ] ' 1 1− ε̃ · ( 1− p+ √ K ·
∣∣(1− p)µλ̂− pµ(1− λ̂)∣∣) · ‖W (0,0) −W ∗‖2 + ε̃p
(1− ε̃) · ‖W (0,0) −W ∗‖2
= 1 1− ε̃ · ( 1− p+ µ √ K ∣∣λ̂− p∣∣) · ‖W (0,0) −W ∗‖2 + ε̃p (1− ε̃) · ‖W (0,0) −W ∗‖2,
(64)
where λ̂ = λδ 2
λδ2+λ̃δ̃2 .
To guarantee the convergence in the outer loop, we require
lim T→∞
‖W (`,T ) −W [p]‖2 ≤ ‖W (0,0) −W [p]‖2 = p‖W (0,0) −W ∗‖2,
and lim T→∞
‖W (`,T ) −W ∗‖2 ≤ ‖W (0,0) −W ∗‖2. (65)
Since we have
‖W (`,T ) −W [p]‖2 ≤‖W (`,T ) −W ∗‖2 + ‖W ∗ −W [p]‖2 =‖W (`,T ) −W ∗‖2 + (1− p) · ‖W ∗ −W (0,0)‖2,
(66)
it is clear that (65) holds if and only if
1 1− ε̃ · ( 1− p+ ε̃p+ µ √ K ∣∣λ̂− p∣∣)+ 1− p ≤ 1. (67)
To guarantee the iterates strictly converges to the desired point, we let
1 1− ε̃ · ( 1− p+ ε̃p+ µ √ K ∣∣λ̂− p∣∣)+ 1− p ≤ 1− 1 C (68)
for some larger constant C, which is equivalent to
|p− λ̂| ≤ 2(1− ε̃)p− 1 µ √ K . (69)
To make the bound in (69) meaningful, we need
p ≥ 1 2(1− ε̃) . (70)
When ` > 1, following similar steps in (64), we have
lim T→∞
‖W (`,T ) −W [p]‖2
≤ 1 1− ε̃
· ( 1− p+ µ √ K ∣∣(λ̂− p)∣∣) · ‖W (0,0) −W ∗‖2 + ε̃p
1− ε̃ · ‖W (`,0) −W [p]‖2,
(71)
Given (69) holds, from (71), we have
lim L→∞,T→∞
‖W (L,T ) −W [p]‖2
≤ 1 1− ε̃
· ( 1− p+ µ √ K ∣∣λ̂− p∣∣) · ‖W (0,0) −W ∗‖2
≤ 1 1− ε̃
· ( 1− p+ µ √ K ∣∣λ̂− p∣∣) · ‖W (0,0) −W ∗‖2.
(72)
G DEFINITION AND RELATIVE PROOFS OF ρ
In this section, the formal definition of ρ is included in Definition 3, and a corresponding claim about ρ is summarized in Lemma 11. One can quickly check that the ReLU activation function satisfies the conditions in Lemma 11.
The major idea in proving Lemma 11 is to show Hr(δ) and Jr(δ) in Definition 3 are in the order of δr when δ is small. Definition 3. Let Hr(δ) = Ez∼N (0,δ2) ( φ′(σKz)z r ) and Jr(δ) = Ez∼N (0,δ2) ( φ′2(σKz)z r ) . Then, ρ = ρ(δ) is defined as
ρ(δ) = min { J0(δ)−H20 (δ)−H21 (δ), J2(δ)−H21 (δ)−H22 (δ), H0(δ) ·H2(δ)−H21 (δ) } , (73)
where σK is the minimal singular value ofW ∗.
Lemma 11 (Order analysis of ρ). If ρ(δ) > 0 for δ ∈ (0, ξ) for some positive constant ξ and the sub-gradient of ρ(δ) at 0 can be non-zero, then ρ(δ) = Θ(δ2) when δ → 0+. Typically, for ReLU activation function, µ in (5) is a fixed constant for all δ, δ̃ ≤ 1.
Proof of Lemma 11 . From Definition 3, we know that Hr(δ) = Ez∼N (0,δ2)φ′(σKz)zr. Suppose we have Hr(δ) = Θ(δr) and Jr(δ) = Θ(δr), then from (73) we have
J0(δ)−H20 (δ)−H21 (δ) ∈ Θ(1)−Θ(δ2), J2(δ)−H21 (δ)−H22 (δ) ∈ Θ(δ2), H0(δ) ·H2(δ)−H21 (δ) ∈ Θ(δ2)−Θ(δ4).
(74)
Because ρ is a continuous function with ρ(z) > 0 for some z > 0. Therefore, ρ 6= J0(δ)−H20 (δ)− H21 (δ) when δ → 0+, otherwise ρ(z) < 0 for any z > 0. When δ → 0+, both J2(δ) − H21 (δ) − H22 (δ) and H0(δ) ·H2(δ)−H21 (δ) are in the order of δ2, which indicates that µ is a fixed constant when both δ and δ̃ are close to 0. In addition, J2(δ) − H21 (δ) − H22 (δ) goes to +∞ while both J0(δ)−H20 (δ)−H21 (δ) and H0(δ) ·H2(δ)−H21 (δ) go to −∞ when δ → +∞. Therefore, with a large enough δ, we have
ρ(δ) ∈ Θ(δ2)−Θ(δ4) or Θ(1)−Θ(δ2), (75)
which indicates that µ is a strictly decreasing function when δ and δ̃ are large enough.
Next, we provide the conditions that guaranteeHr(δ) = Θ(δr) hold, and the relative proof for Jr(δ) can be derived accordingly following the similar steps as well. From Definition 3, we have
lim δ→0+
Hr(δ)
δr = lim δ→0+ ∫ +∞ −∞ φ′(σKz) (z δ )r 1√ 2πδ e− z2 δ2 dz
(a) = lim
δ→0+ ∫ +∞ −∞ φ′(σKδt) tr√ 2π e−t 2 dt
= lim δ→0+ ∫ 0− −∞ φ′(σKδt) tr√ 2π e−t 2 dt+ lim δ→0+ ∫ +∞ 0+ φ′(σKδt) tr√ 2π e−t 2 dt
=φ′(0−) ∫ 0− −∞ tr√ 2π e−t 2 dt+ φ′(0+) ∫ +∞ 0+ tr√ 2π e−t 2 dt,
(76)
where equality (a) holds by letting t = zδ . It is easy to verify that∫ +∞ 0+ tr√ 2π e−t 2 dt = (−1)r ∫ 0− −∞ tr√ 2π e−t 2 dt,
and both are bounded for a fixed r. Thus, as long as either φ′(0−) or φ′(0+) is non-zero, we have Hr(δ) = Θ(δ
r) when δ → 0+. If φ has bounded gradient as |φ′| ≤ Cφ for some positive constant Cφ. Then, we have∣∣Hr(δ)
δr ∣∣ =∣∣∣ ∫ +∞ −∞ φ′(σKz) (z δ )r 1√ 2πδ e− z2 δ2 dz ∣∣∣
= ∣∣∣ ∫ +∞ −∞ φ′(σKδt) tr√ 2π e−t 2 dt ∣∣∣
≤Cφ · ∣∣∣ ∫ +∞ −∞ tr√ 2π e−t 2 dt ∣∣∣
(77)
Therefore, we have Hr(δ) = O(δr) for all δ > 0 when φ has bounded gradient. Typiclly, for ReLU function, one can directly calculate that Hr(δ) = δr for δ ∈ R, and ρ(δ) = Cδ2 when δ ≤ 1 for some constant C = 0.091. Then, it is easy to check that µ is a constant when δ, δ̃ ≤ 1.

heading: H PROOF OF PRELIMINARY LEMMAS
text: H.1 PROOF OF LEMMA 1
The eigenvalues of ∇2f(·; p) at any fixed point W can be bounded in the form of (80) by Weyl’s inequality (Lemma 4). Therefore, the primary technical challenge lies in bounding ‖∇2f(W ; p)− ∇2f(W [p]; p)‖2, which is summarized in Lemma 12. Lemma 13 provides the exact calulation of the lower bound of Ex (∑K j=1α T j xφ ′(w [p]T j x) )2 when x belongs to Gaussian distribution with zero mean, which is used in proving the lower bound of the Hessian matrix in (81).
Lemma 12. Let f(W ; p) be the population risk function defined in (17) with p and W satisfying (20). Then, we have
‖∇2f(W [p]; p)−∇2f(W ; p)‖2 . λδ2 + (1− λ)δ̃2 K · ‖W [p] −W ‖2 σK . (78)
Lemma 13 (Lemma D.6, (Zhong et al., 2017)). For any {wj}Kj=1 ∈ Rd, let α ∈ RdK be the unit vector defined in (19). When the φ is ReLU function, we have
min ‖α‖2=1 Ex∼N (0,σ2) ( K∑ j=1 αTj xφ ′(wTj x) )2 & ρ(σ), (79)
where ρ(σ) is defined in Definition 3.
Proof of Lemma 1. Let λmax(W ) and λmin(W ) denote the largest and smallest eigenvalues of ∇2f(W ; p) at pointW , respectively. Then, from Lemma 4, we have
λmax(W ) ≤ λmax(W [p]) + ‖∇2f(W ; p)−∇2f(W [p]; p)‖2, λmin(W ) ≥ λmin(W [p])−‖∇2f(W ; p)−∇2f(W [p]; p)‖2.
(80)
Then, we provide the lower bound of the Hessian matrix of the population function at W [p]. For any α ∈ RdK defined in (19) with ‖α‖2 = 1, we have
min ‖α‖2=1
αT∇2f(W [p]; p)α
= 1
K2 min ‖α‖2=1
[ λEx ( K∑ j=1 αTj xφ ′(w [p]T j x) )2 + λ̃Ex̃ ( K∑ j=1 αTj x̃φ ′(w [p]T j x̃) )2]
≥ 1 K2 min ‖α‖2=1 λEx ( K∑ j=1 αTj xφ ′(w [p]T j x) )2 + min ‖α‖2=1 λ̃Ex̃ ( K∑ j=1 αTj x̃φ ′(w [p]T j x̃) )2 ≥λρ(δ) + λ̃ρ(δ̃)
11κ2γK2 ,
(81)
where the last inequality comes from Lemma 13.
Next, the upper bound can be bounded as
max ‖α‖2=1
αT∇2f(W [p]; p)α
= 1
K2 max ‖α‖2=1
[ λEx ( K∑ j=1 αTj xφ ′(w [p]T j x) )2 + λ̃Ex̃ ( K∑ j=1 αTj x̃φ ′(w [p]T j x̃) )2]
≤ 1 K2 max ‖α‖2=1 λEx ( K∑ j=1 αTj xφ ′(w [p]T j x) )2 + max ‖α‖2=1 λ̃Ex̃ ( K∑ j=1 αTj x̃φ ′(w [p]T j x̃) )2 .
(82)
For Ex (∑K j=1α T j xφ ′(w [p]T j x) )2 , we have
Ex ( K∑ j=1 αTj xφ ′(w [p]T j x) )2 =Ex
K∑ j1=1 K∑ j2=1 αTj1xφ ′(w [p]T j1 x)αTj2xφ ′(w [p]T j2 x)
= K∑ j1=1 K∑ j2=1 ExαTj1xφ ′(w [p]T j1 x)αTj2xφ ′(w [p]T j2 x)
≤ K∑ j1=1 K∑ j2=1 [ Ex(αTj1x) 4Ex(φ′(w[p]Tj1 x)) 4Ex(αTj2x) 4Ex(φ′(w[p]Tj2 x)) 4 ]1/4
≤ K∑ j1=1 K∑ j2=1 3δ2‖αj1‖2‖αj2‖2
≤6δ2 K∑ j1=1 K∑ j2=1 1 2 (‖αj1‖22 + ‖αj2‖22)
=6Kδ2
(83)
Therefore, we have
max ‖α‖2=1
αT∇2f(W [p]; p)α
≤ 1 K2 max ‖α‖2=1 λEx ( K∑ j=1 αTj xφ ′(w [p]T j x) )2 + max ‖α‖2=1 λ̃Ex̃ ( K∑ j=1 αTj x̃φ ′(w [p]T j x̃) )2 ≤6(λδ 2 + λ̃δ̃2)
K .
(84)
Then, given (20), we have
‖W (0,0) −W [p]‖F = p‖W (0,0) −W ∗‖F . σK µ2K . (85)
Combining (85) and Lemma 12, we have
‖∇2f(W ; p)−∇2f(W [p]; p)‖2 . λρ(δ) + λ̃ρ(δ̃)
132κ2γK2 . (86)
Therefore, (86) and (80) completes the whole proof.
H.2 PROOF OF LEMMA 2
The task of bounding of the quantity between ‖∇f̂ −∇f‖2 is dividing into bounding I1, I2, I3 and I4 as shown in (89). I1 and I3 represent the deviation of the mean of several random variables to their expectation, which can be bounded through concentration inequality, i.e, Chernoff bound. I2 and I4 come from the inconsistency of the output label y and pseudo label ỹ in the empirical risk function in (1) and population risk function in (17). The major challenge lies in characterizing the upper bound of I2 and I4 as the linear function of W̃ −W [p] andW [p]−W ∗, which is summarized in (96).
Proof of Lemma 2. From (1), we know that
∂f̂
∂wk (W ) =
λ
N N∑ n=1 ( 1 K K∑ j=1 φ(wTj xn)− yn ) xn + 1− λ M M∑ m=1 ( 1 K K∑ j=1 φ(wTj x̃m)− ỹn ) x̃m
= λ
K2N N∑ n=1 K∑ j=1 ( φ(wTj xn)− φ(w∗Tj xn) ) xn
+ 1− λ K2M M∑ m=1 K∑ j=1 ( φ(wTj x̃m)− φ(w̃Tj x̃m) ) x̃m.
(87) From (32), we know that
∂f̂
∂wk (W ) =
λ
K2 Ex K∑ j=1 ( φ(wTj x)− φ(w∗Tj x) ) x+ 1− λ K2 Ex̃ K∑ j=1 ( φ(wTj x̃)− φ(w̃Tj x̃) ) x̃.
(88) Then, from (17), we have
∂f̂ ∂wk (W )− ∂f ∂wk (W ; p)
= λ
K2N K∑ j=1 [ N∑ n=1 ( φ(wTj xn)− φ(w∗Tj xn) ) xn − Ex ( φ(wTj x)− φ(w [p]T j x) ) x ]
+ 1− λ K2M K∑ j=1 [ M∑ m=1 ( φ(wTj x̃m)− φ(w̃Tj x̃m) ) x̃m − Ex̃ ( φ(wTj x̃)− φ(w [p]T j x̃) ) x̃ ]
= λ
K2 K∑ j=1 [ 1 N N∑ n=1 ( φ(wTj xn)− φ(w∗Tj xn) ) xn − Ex ( φ(wTj x)− φ(w∗Tj x) ) x ]
+ λ
K2 K∑ j=1 Ex [( φ(w∗Tj x)− φ(w [p]T j x) ) x ]
+ 1− λ K2 K∑ j=1 [ 1 M M∑ m=1 ( φ(wTj x̃m)− φ(w̃Tj x̃m) ) x̃m − Ex̃ ( φ(wTj x̃)− φ(w̃Tj x̃) ) x̃ ]
+ 1− λ K2 K∑ j=1 Ex̃ [( φ(w̃Tj x̃)− φ(w [p]T j (p)x̃) ) x̃ ]
:=I1 + I2 + I3 + I4.
(89)
For anyαj ∈ Rd with ‖αj‖2 ≤ 1, we define a random variableZ(j) = ( φ(wTj x)−φ(w∗Tj x) ) αTj x
and Zn(j) = ( φ(wTj xn)−φ(w∗Tj xn) ) αTj xn as the realization of Z(j) for n = 1, 2 · · · , N . Then, for any p ∈ N+, we have( E|Z|p )1/p = ( E|φ(wTj x)− φ(w∗Tj x)|p · |αTj x|p
)1/p ≤ ( E|(wj −w∗j )Tx|p · |αTj x|p
)1/p ≤C · δ2‖wj −w∗j ‖2 · p,
(90)
where C is a positive constant and the last inequality holds since x ∼ N (0, δ2). From Definition 2, we know that Z belongs to sub-exponential distribution with ‖Z‖ψ1 . δ2‖wj −w∗j ‖2. Therefore, by Chernoff inequality, we have
P {∣∣∣ 1 N N∑ n=1 Zn(j)− EZ(j) ∣∣∣ < t} ≤ 1− e−C(δ2‖wj−w∗j ‖2)2·Ns2 eNst (91)
for some positive constant C and any s ∈ R. Let t = δ2‖wj −w∗j ‖2 √ d log q N and s = 2 Cδ2‖wj−w∗j ‖2 · t for some large constant q > 0, we have
∣∣∣ 1 N N∑ n=1 Zn(j)− EZ(j) ∣∣∣ . δ2‖wj −w∗j ‖2 ·√d log qN (92)
with probability at least 1− q−d. From Lemma 7, we have∥∥∥ 1 N N∑ n=1 ( φ(wTj xn)− φ(w∗Tj xn) ) xn − Ex ( φ(wTj x)− φ(w∗Tj x) ) x ∥∥∥ 2
≤2δ2‖wj −w∗j ‖2 · √ d log q
N
(93)
with probability at least 1 − (q/5)−d. Since q is a large constant, we release the probability as 1− q−d for simplification. Similar to Z, we have∥∥∥ 1
M M∑ m=1 ( φ(wTj x̃m)− φ(w̃Tj x̃m) ) x̃m − Ex̃ ( φ(wTj x̃)− φ(w̃Tj x̃) ) x̃ ∥∥∥ 2
.δ̃2‖wj − w̃j‖2 · √ d log q
M
(94)
with probability at least 1− q−d.
For term Ex [(
φ(w∗Tj x)−φ(w [p]T j x)
)
x
]
, let us define the angle betweenw∗j andw [p] j as θ1. Figure 13 shows the subspace spanned by the vectorw∗j and w̃j . We divide the subspace by 4 pieces, where the gray region denotes area I, and the blue area denotes area II. Areas III and IV are the symmetries of II and I from the origin, respectively. Hence, we have
Ex [( φ(w∗Tj x)− φ(w [p]T j x) ) x ]
=Ex∈ area I [( φ(w∗Tj x)− φ(w [p]T j x) ) x ] + Ex∈ area II [( φ(w∗Tj x)− φ(w [p]T j x) ) x ]
+ Ex∈ area III [( φ(w∗Tj x)− φ(w [p]T j x) ) x ] + Ex∈ area IV [( φ(w∗Tj x)− φ(w [p]T j x) ) x ]
=Ex∈ area I [( φ(w∗Tj x)− φ(w [p]T j x) ) x ] + Ex∈ area II [ w∗Tj x̃x̃ ] − Ex∈ area III [ w [p]T j xx ] =Ex∈ area I [ (w∗j −w [p] j ) Txx ] + Ex∈ area II [ (w∗j −w [p] j ) Txx ]
= 1 2 Ex [ (w∗j −w [p] j ) Txx ]
(95)
Therefore, we have∥∥∥∥ λK2Ex[(φ(w∗Tj x)− φ(w[p]Tj x))x]+ 1− λK2 Ex̃[(φ(w̃Tj x̃)− φ(w[p]Tj x̃))x̃] ∥∥∥∥ 2
= ∥∥∥∥ λ2K2Ex[(w∗j −w[p]j )Txx]+ 1− λ2K2 Ex̃[(w̃j −w[p]j )Txx] ∥∥∥∥ 2
= ∥∥λδ2 · (w̃j −w[p]j )+ (1− λ)δ̃2 · (w∗j −w[p]j )∥∥2 2K2 .
(96)
From (93), (94) and (96), we have∥∥∥ ∂f̂ ∂wk (W ; p)− ∂f ∂wk (W ) ∥∥∥ 2
≤ λ K2 K∑ j=1 ∥∥∥ 1 N N∑ n=1 ( φ(wTj xn)− φ(w∗Tj xn) ) xn − Ex ( φ(wTj x)− φ(w∗Tj x) ) x ∥∥∥ 2
+ 1− λ K2 K∑ j=1 ∥∥∥ 1 M M∑ m=1 ( φ(wTj x̃m)− φ(w̃Tj x̃m) ) x̃m − Ex̃ ( φ(wTj x̃)− φ(w̃Tj x̃) ) x̃ ∥∥∥ 2
+ K∑ j=1 ∥∥∥∥ λK2Ex[(φ(w∗Tj x)− φ(w[p]Tj x))x]+ 1− λK2 Ex̃[(φ(w̃Tj x̃)− φ(w[p]Tj x̃))x̃] ∥∥∥∥ 2
≤ λ K2
δ2 √ d log q N · K∑ j=1 ‖wj −w∗j ‖2 + 1− λ K2 · δ̃2 √ d log q M · K∑ j=1 ‖wj − w̃j‖2
+ 1 2K2 · K∑ j=1 ∥∥λδ2 · (w̃j −w[p]j )+ λ̃δ̃2 · (w∗j −w[p]j )∥∥2 ≤ λ K3/2 δ2 √ d log q N · ‖W −W ∗‖2 + 1− λ K3/2 · δ̃2 √ d log q M · ‖W − W̃ ‖2
+ 1
2K3/2
∥∥λδ2 · (W̃ −W [p])+ (1− λ)δ̃2 · (W ∗ −W [p])∥∥ 2
(97)
with probability at least 1− q−d. In conclusion, let α ∈ RKd and αj ∈ Rd with α = [αT1 ,αT2 , · · · ,αTK ]T , we have
‖∇f(W )−∇f̂(W )‖2 = ∣∣∣αT (∇f(W )−∇f̂(W ))∣∣∣ ≤
K∑ k=1 ∣∣∣αTk ( ∂f̂∂wk (W )− ∂f∂wk (W )) ∣∣∣
. K∑ k=1 ∥∥∥ ∂f̂ ∂wk (W )− ∂f ∂wk (W ) ∥∥∥ 2 · ‖αk‖2
. λ K δ2 √ d log q N · ‖W −W ∗‖2 + 1− λ K · δ̃2 √ d log q M · ‖W − W̃ ‖2
+ 1
2K
∥∥λδ2 · (W̃ −W [p])+ (1− λ)δ̃2 · (W ∗ −W [p])∥∥ 2
(98)
with probability at least 1− q−d.
H.3 PROOF OF LEMMA 12
The distance of the second order derivatives of the population risk function f(·; p) at point W and W [p] can be converted into bounding P1, P2, P3 and P4, which are defined in (101). The major
idea in proving P1 is to connect the error bound to the angle between W and W [p]. Similar ideas apply in bounding the other three items as well.
Proof of Lemma 12. From (17), we have
∂2f
∂wj1∂wj2 (W [p]; p) =
λ
K2 Exφ′(w[p]Tj1 x)φ ′(w [p]T j2 x)xxT + 1− λ K2 Ex̃φ′(w [p]T j1 x̃)φ′(w [p]T j2 x̃)x̃x̃T ,
(99)
and ∂2f
∂wj1∂wj2 (W ; p) =
λ
K2 Exφ′(wTj1x)φ ′(wTj2x)xx T + 1− λ K2 Ex̃φ′(wTj1 x̃)φ ′(wTj2 x̃)x̃x̃ T ,
(100)
where w[p]j is the j-th column ofW [p]. Then, we have
∂2f
∂wj1∂wj2 (W ∗)− ∂
2f
∂wj1∂wj2 (W )
= λ K2 Ex [ φ′(w [p]T j1 x)φ′(w [p]T j2 x)− φ′(wTj1x)φ ′(wTj2x) ] xxT
+ 1− λ K2
Ex̃ [ φ′(w [p]T j1 x̃)φ′(w [p]T j2 x̃)− φ′(wTj1 x̃)φ ′(wTj2 x̃) ] x̃x̃T
= λ K2 Ex [ φ′(w [p]T j1 x) ( φ′(w [p]T j2 x)− φ′(wTj2x) ) + φ′(wTj2x) ( φ′(w [p]T j1 x)− φ′(wTj1x) )] xxT
+ 1− λ K2
Ex̃ [ φ′(w [p]T j1 x̃) ( φ′(w [p]T j2 x̃)− φ′(wTj2 x̃) ) + φ′(wTj2 x̃) ( φ′(w [p]T j1 x̃)− φ′(wTj1 x̃) )] x̃x̃T
= λ
K2
[ Exφ′(w[p]Tj1 x) ( φ′(w [p]T j2 x)− φ′(wTj2x) ) xxT + Exφ′(wTj2x) ( φ′(w [p]T j1 x)− φ′(wTj1x) ) xxT ] +
1− λ K2
[ Ex̃φ′(w [p]T j1 x̃) ( φ′(w [p]T j2 x̃)− φ′(wTj2 x̃) ) x̃x̃T + Ex̃φ′(wTj2 x̃) ( φ′(w [p]T j1 x̃)− φ′(wTj1 x̃) ) x̃x̃T ] := λ
K2 (P1 + P2) + 1− λ K2 (P3 + P4).
(101)
For any a ∈ Rd with ‖a‖2 = 1, we have aTP1a =Exφ′(w[p]Tj1 x) ( φ′(w [p]T j2 x)− φ′(wTj2x) ) (aTx)2 (102)
where a ∈ Rd. Let I = φ′(w[p]Tj1 x) ( φ′(w [p]T j2 x)−φ′(wTj2x) ) ·(aTx)2. It is easy to verify there exists a group of orthonormal vectors such that B = {a, b, c,a⊥4 , · · · ,a⊥d } with {a, b, c} spans a sub-
space that contains a,wj2 andw ∗ j2 . Then, for any x, we have a unique z = [ z1, z2, · · · , zd ]T such that x = z1a+ z2b+ z3c+ · · ·+ zda⊥d . Also, since x ∼ N (0, δ2Id), we have z ∼ N (0, δ2Id). Then, we have
I =Ez1,z2,z3 |φ′ ( wTj2x ) − φ′ ( w [p]T j2 x ) | · |aTx|2
= ∫ |φ′ ( wTj2x ) − φ′ ( w [p]T j2 x ) | · |aTx|2 · fZ(z1, z2, z3)dz1dz2dz3,
where x = z1a+ z2b+ z3c and fZ(z1, z2, z3) is probability density function of (z1, z2, z3). Next, we consider spherical coordinates with z1 = Rcosφ1, z2 = Rsinφ1sinφ2, z3 = Rsinφ1cosφ2. Hence,
I = ∫ |φ′ ( wTj2x ) − φ′ ( w [p]T j2 x ) | · |R cosφ1|2 · fZ(R,φ1, φ2)R2 sinφ1dRdφ1dφ2. (103)
It is easy to verify that φ′ ( wTj2x ) only depends on the direction of x and
fZ(R,φ1, φ2) = 1
(2πδ2) 3 2
e− z21+z 2 2+z 2 3 2δ2 = 1
(2πδ2) 3 2
e− R2 2δ2
only depends on R. Then, we have
I(i2, j2)
= ∫ |φ′ ( wTj2(x/R) ) − φ′ ( w [p]T j2 (x/R) ) | · |R cosφ1|2 · fZ(R)R2 sinφ1dRdφ1dφ2
= ∫ ∞ 0 R4fz(R)dR ∫ π 0 ∫ 2π 0 | cosφ1|2 · sinφ1 · |φ′ ( wTj2(x/R) ) − φ′ ( w [p]T j2 (x/R) ) |dφ1dφ2
(a) ≤3δ2 · ∫ ∞ 0 R2fz(R)dR ∫ π 0 ∫ 2π 0 sinφ1 · |φ′ ( wTj2(x/R) ) − φ′ ( w [p]T j2 (x/R) ) |dφ1dφ2
=3δ2 · Ez1,z2,z3 ∣∣φ′(wTj2x)− φ′(w[p]Tj2 x)|
≤3δ2 · Ex ∣∣φ′(wTj2x)− φ′(w[p]Tj2 x)|,
(104)
where the inequality (a) is derived from the fact that |cosφ1| ≤ 1 and∫ ∞ 0 R4 1 (2πδ2) 3 2 e− R2 2δ2 dR = ∫ ∞ 0 − R 3δ2 (2πδ2) 3 2 d(e− R2 2δ2 )
= ∫ ∞ 0 e− R2 2δ2 d R3δ2 (2πδ2) 3 2
=3δ2 ∫ ∞ 0 R2 1 (2πδ2) 3 2 e− R2 2δ2 dR.
(105)
Define a set A1 = {x|(w[p]Tj2 x)(w T j2 x) < 0}. If x ∈ A1, then w[p]Tj2 x and w T j2 x have different signs, which means the value of φ′(wTj2x) and φ ′(w [p]T j2
x) are different. This is equivalent to say that
|φ′(wTj2x)− φ ′(w [p]T j2
x)| = {
1, if x ∈ A1 0, if x ∈ Ac1 . (106)
Moreover, if x ∈ A1, then we have
|w[p]Tj2 x| ≤|w [p]T j2 x−wTj2x| ≤ ‖w [p] j2 −wj2‖2 · ‖x‖2. (107)
Let us define a set A2 such that
A2 = { x ∣∣∣ |w[p]Tj2 x|‖w∗j2‖2‖x‖2 ≤ ‖w ∗ j2 −wj2‖2 ‖w∗j2‖2 } = { θx,w∗j2 ∣∣∣| cos θx,w[p]j2 | ≤ ‖w [p] j2 −wj2‖2 ‖w[p]j2 ‖2 } . (108)
Hence, we have that
Ex|φ′(wTj2x)− φ ′(w [p]T j2 x)|2 =Ex|φ′(wTj2x)− φ ′(w [p]T j2 x)| =Prob(x ∈ A1) ≤Prob(x ∈ A2).
(109)
Since x ∼ N (0, δ2‖a‖22I), θx,w[p]j2 belongs to the uniform distribution on [−π, π], we have
Prob(x ∈ A2) = π − arccos
‖w[p]j2 −wj2‖2 ‖w[p]j2 ‖2
π ≤ 1 π tan(π − arccos ‖w[p]j2 −wj2‖2 ‖w[p]j2 ‖2 )
= 1
π cot(arccos ‖w[p]j2 −wj2‖2 ‖w[p]j2 ‖2 )
≤ 2 π ‖w[p]j2 −wj2‖2 ‖w[p]j2 ‖2 .
(110)
Hence, (104) and (110) suggest that
I ≤ 6δ 2
π
‖wj2 −w [p] j2 ‖2
σK · ‖a‖22. (111)
The same bound that shown in (111) holds for P2 as well.
P3 and P4 satisfy (111) except for changing δ2 to δ̃2.
Therefore, we have
‖∇2f(W [p]; p)−∇2f(W ; p)‖2
= max ‖α‖2≤1 ∣∣∣αT (∇2f(W [p]; p)−∇2f(W ; p))α∣∣∣ ≤
K∑ j1=1 K∑ j2=1 ∣∣∣∣∣αTj1 ( ∂2f ∂wj1∂wj2 (W [p]; p)− ∂ 2f ∂wj1∂wj2 (W ; p) ) αj2 ∣∣∣∣∣ ≤ 1 K2 K∑ j1=1 K∑ j2=1 ( λ‖P1 + P2‖2 + (1− λ)‖P3 + P4‖2 ) ‖αj1‖2‖αj2‖2
≤ 1 K2 K∑ j1=1 K∑ j2=1 4(λδ2 + (1− λ)δ̃2) ‖w[p]j2 −wj2‖2 σK ‖αj1‖2‖αj2‖2
≤ 4 K
( λδ2 + (1− λ)δ̃2 ) · ‖W
[p] −W ‖2 σK ,
(112)
where α ∈ RKd and αj ∈ Rd with α = [αT1 ,αT2 , · · · ,αTK ]T .
I INITIALIZATION VIA TENSOR METHOD
In this section, we briefly summarize the tensor initialization in (Zhong et al., 2017) by studying the target function class as
y = 1
K K∑ j=1 v∗jφ(w ∗T j x), (113)
where v∗j ∈ R. Note that for ReLU function, we have v∗jφ(w∗Tj x) = sign(v∗j )φ(|v∗j |w∗Tj x). Without loss of generalization, we can assume v∗j ∈ {+1,−1}. Additionally, it is clear that the function studied in (2) is the special case of (113) when v∗j = 1 for all j. In addition, Theorem 5.6 in (Zhong et al., 2017) show that the sign of v∗j can be directly recovered using tensor initialization, which indicates the the equivalence of (2) and (113) when using tensor initialization.
We first define some high order momenta in the following way:
M1 = Ex{yx} ∈ Rd, (114)
M2 = Ex [ y ( x⊗ x− δ2I )] ∈ Rd×d, (115)
M3 = Ex [ y ( x⊗3 − x⊗̃δ2I )] ∈ Rd×d×d, (116)
where Ex is the expectation over x and z⊗3 := z ⊗ z ⊗ z. The operator ⊗̃ is defined as
v⊗̃Z = d2∑ i=1 (v ⊗ zi ⊗ zi + zi ⊗ v ⊗ zi + zi ⊗ zi ⊗ v), (117)
for any vector v ∈ Rd1 and Z ∈ Rd1×d2 .
Following the same calculation formulas in the Claim 5.2 (Zhong et al., 2017), there exist some known constants ψi, i = 1, 2, 3, such that
M1 = K∑ j=1 ψ1 · ‖w∗j ‖2 ·w∗j , (118)
M2 = K∑ j=1 ψ2 · ‖w∗j ‖2 ·w∗jw∗Tj , (119)
M3 = K∑ j=1 ψ3 · ‖w∗j ‖2 ·w∗⊗3j , (120)
where w∗j = w ∗ j /‖w∗j ‖2 in (114)-(116) is the normalization of w∗j . Therefore, we can see that the information of {w∗j }Kj=1 are separated as the direction of wj and the magnitude of wj in M1, M2 and M3.
M1, M2 and M3 can be estimated through the samples { (xn, yn) }N n=1
, and let M̂1, M̂2, M̂3 denote the corresponding estimates. First, we will decompose the rank-K tensorM3 and obtain the {w∗j}Kj=1. By applying the tensor decomposition method (Kuleshov et al., 2015) to M̂3, the outputs, denoted by ŵ ∗ j , are the estimations of {sjw∗j}Kj=1, where sj is an unknown sign. Second, we will estimate sj , v∗j and ‖w∗j ‖2 throughM1 andM2. Note thatM2 does not contain the information of sj because s2j is always 1. Then, through solving the following two optimization problem:
α̂1 = arg min α1∈RK
: ∣∣∣M̂1 − K∑
j=1
ψ1α1,jŵ ∗ j ∣∣∣, α̂2 = arg min α2∈RK : ∣∣∣M̂2 − K∑
j=1
ψ2α2,jŵ ∗ j ŵ ∗T j ∣∣∣, (121) The estimation of sj can be given as
ŝj = sign(α̂1,j/α̂2,j).
Also, we know that |α̂1,j | is the estimation of ‖w∗j ‖ and v̂j = sign(α̂1,j/sj) = sign(α̂2,j). Thus,W (0) is given as [ sign(α̂2,1)α̂1,1ŵ ∗ 1, · · · , sign(α̂2,K)α̂1,Kŵ ∗ K ] .
Subroutine 1 Tensor Initialization Method 1: Input: labeled data D = {(xn, yn)}Nn=1; 2: Partition D into three disjoint subsets D1, D2, D3; 3: Calculate M̂1, M̂2 following (114), (115) using D1, D2, respectively; 4: Obtain the estimate subspace V̂ of M̂2; 5: Calculate M̂3(V̂ , V̂ , V̂ ) through D3; 6: Obtain {ŝj}Kj=1 via tensor decomposition method (Kuleshov et al., 2015) on M̂3(V̂ , V̂ , V̂ ); 7: Obtain α̂1, α̂2 by solving optimization problem (121); 8: Return: w(0)j = sign(α̂2,j)α̂1,jV̂ ûj and v (0) j = sign(α̂2,j), j = 1, ...,K.
To reduce the computational complexity of tensor decomposition, one can project M̂3 to a lowerdimensional tensor (Zhong et al., 2017). The idea is to first estimate the subspace spanned by {w∗j }Kj=1, and let V̂ denote the estimated subspace. Moreover, we have
M3(V̂ , V̂ , V̂ ) = Ex [ y ( (V̂ Tx)⊗3 − (V̂ Tx)⊗̃Ex(V̂ Tx)(V̂ Tx)T )] ∈ RK×K×K , (122)
Then, one can decompose the estimate M̂3(V̂ , V̂ , V̂ ) to obtain unit vectors {ŝj}Kj=1 ∈ RK . Since w∗ lies in the subspace V , we have V V Tw∗j = w ∗ j . Then, V̂ ŝj is an estimate of w ∗ j . The initialization process is summarized in Subroutine 1.

heading: J CLASSIFICATION PROBLEMS
text: The framework in this paper is extendable to binary classification problem. For binary classification problem, the output y given input x is defined as
Prob{y = 1} = g(W ∗;x) (123)
with some ground truth parameter W ∗. To guarantee the output is within [0, 1], the activation function is often used as sigmoid. For classification, the loss function is cross-entropy, and the objective function over labeled data D is defined as
fD(W ) = 1
N ∑ (xn,yn)∈D −yn log g(W ;xn)− (1− yn) log(1− g(W ;xn)). (124)
The expectation of objective function can be written as
EDfD(W ) =E(x,y) − y log(g(W ;xn))− (1− y) log(1− g(W ;x)) =ExE(y|x) − y log(g(W ;xn))− (1− y) log(1− g(W ;x))
=Ex [ − g(W ∗;x) log(g(W ;xn))− (1− g(W ∗;x)) log(1− g(W ;xn)) ] (125) Please note that (125) is exactly the same as (32) with λ = 1 when the loss function is squared loss.
For cross entropy loss function, the second order derivative of (125) is calculated as
∂fD(W ) ∂wj∂wk = 1 N [ yn g2(W ;x) + 1− yn (1− g(W ;x))2 ] · φ′(wTj x)φ′(wTk x)xxT . (126)
when j 6= k. Refer to (88) in (Fu et al., 2020) or (132) in (Zhang et al., 2020b), we have∥∥∥yn(φ′(wTj x)φ′(wTk x)) g2(W ;x) ∥∥∥ 2 ≤ ∥∥∥φ′(wTj x)φ′(wTk x) g2(W ;x) ∥∥∥ 2 ≤ K2. (127)
Following similar steps in (90), from Defintion 2, we know that αTj ∂fD(W ) ∂wj∂wk αk belongs to the subexponential distribution. Therefore, similiar results for objective function with cross-entropy loss can be established as well. One can check (Fu et al., 2020) or (Zhang et al., 2020b) for details.

heading: K ONE-HIDDEN LAYER NEURAL NETWORK WITH TOP LAYER WEIGHTS
text: For a general one-hidden layer neural network, the output of the neural network is defined as
g(W ,v;x) = 1
K K∑ j=1 vjφ(w T j x), (128)
where v = [v1, v2, · · · , vK ] ∈ RK . Then, the target function can be defined as
y = g(W ∗,v∗;x) = 1
K K∑ j=1 v∗jφ(w ∗T j x) (129)
for some unknown weightsW ∗ and v∗.
In the following paragraphs, we will provide a short description for the equivalence of (129) and (2) in theoretical analysis. Note that for ReLU functions, we have vjφ(wTj x) = sign(vj)φ(|vj |wTj x).
Without loss of generalization, we can assume vj , v∗j ∈ {+1,−1} for all j ∈ [K]5. From Appendix I, we know that the sign of v∗j can exactly estimated through tensor initialization. There, we can focus on analysis the neural network in the form as
g(W ;x) = 1
K K∑ j=1 v∗jφ(w T j x). (130)
Considering the objective function in (1) and population risk function in (17), we have∥∥∥ ∂f̂ ∂wk (W )− ∂f ∂wk (W ; p) ∥∥∥ 2
= ∥∥∥ λ K2N K∑ j=1 v∗j [ N∑ n=1 ( φ(wTj xn)− φ(w∗Tj xn) ) xn − Ex ( φ(wTj x)− φ(w [p]T j x) ) x ]
+ 1− λ K2M v∗j K∑ j=1 [ M∑ m=1 ( φ(wTj x̃m)− φ(w̃Tj x̃m) ) x̃m − Ex̃ ( φ(wTj x̃)− φ(w [p]T j x̃) ) x̃ ]∥∥∥ 2
≤ K∑ j=1 ·|v∗j | · ∥∥∥ λ K2N [ N∑ n=1 ( φ(wTj xn)− φ(w∗Tj xn) ) xn − Ex ( φ(wTj x)− φ(w [p]T j x) ) x ]
+ 1− λ K2M [ M∑ m=1 ( φ(wTj x̃m)− φ(w̃Tj x̃m) ) x̃m − Ex̃ ( φ(wTj x̃)− φ(w [p]T j x̃) ) x̃ ]∥∥∥ 2
= K∑ j=1 ∥∥∥ λ K2N [ N∑ n=1 ( φ(wTj xn)− φ(w∗Tj xn) ) xn − Ex ( φ(wTj x)− φ(w [p]T j x) ) x ]
+ 1− λ K2M [ M∑ m=1 ( φ(wTj x̃m)− φ(w̃Tj x̃m) ) x̃m − Ex̃ ( φ(wTj x̃)− φ(w [p]T j x̃) ) x̃ ]∥∥∥ 2 ,
(131)
which is exact the same as (89). Similar results can be derived for Lemma 12. Therefore, the conclusions and proofs of Lemma 1 and Lemma 2 does not change at all.
Additionally, fixing the second-layer weights and only training the hidden layer is the state-of-theart practice in analyzing two-layer neural networks (Arora et al., 2019b;a; Allen-Zhu et al., 2019; Safran & Shamir, 2018; Li & Liang, 2018; Brutzkus & Globerson, 2017; Oymak & Soltanolkotabi, 2018; Zhang et al., 2019). Additionally, as indicated in (Safran & Shamir, 2018), training a onehidden-layer neural network with all vj fixed as 1 has intractable many spurious local minima, which indicates that training problem is not trivial.
5To see this, one can view |v∗j |w∗j as the new ground truth weights, and the goal for this paper is to recover the new ground truth weights.

Review summaries:
This paper presents the first theoretical analysis of self-training on one hidden-layer neural network with gaussian input. The paper shows that under certain conditions (e.g., the initialization is neither too far or too close to the optimal), iterative self-training on the aformentioned setting can provably recover the ground truth labeling matrix and offers improvement in sample complexity over supervised learning with only labeled data. Experiments validate the theoretical results.

In this paper, the authors make first forrays into the study of the theoretical properties of the following iterative self-training algorithm in a semi-supervised regression setting where we have N labeled samples and M unlabeled samples (and the ground truth is realisable):  (1)  train the model on the labeled data, obtain pseudo labels for the unlabeled data by feeding them through the trained model  (2)  train an auxiliary loss (see equation (1)) on the augmented dataset the includes the unlabeled data with the corresponding pseudolabels (the loss is a convex combination of the empirical loss on the labelled and unlabeled datasets), and  (3) iteratively apply step (2) until convergence.   It is assumed that the ground truth is realisable (though without overparametrization assumptions), the labels are observed exactly (without noise), the model is a simple two-layer neural network with the second layer weights all fixed to one, and the labeled and unlabeled datapoints come from two isotropic Gaussian distributions of variances $\delta$ and $\tilde{\delta}$.   In theorem 2, the following surprising result is shown: under some reasonable assumptions on the initial point and the number of unlabeled datapoints, for a suitable choice of the convex combination parameter $\lambda$, the model converges towards the ground truth with an assumption on the number of samples which is weaker (by a constant factor at best) than what is required in the case where no unlabeled data is provided (which is a particular case of the case studied in [1])  In the more general Theorem 1, it is shown that under much weaker condition (without a requirement on the number of labelled samples), the model still converges towards a given convex combination of the ground truth and the initial point, and is guaranteed to outperform the initial model.    The main idea of the proofs is to define an auxiliary loss referred to as the population risk function (which estimates the risk taking as "labels" the images of the points by a convex combination of the ground truth and the initial point, cf. equation (17)) and to show (1) that the optimization lanscape of this functions is mild  (mainly proved in Lemmas 1 and 11), and (2) that this auxiliary function is close to the empirical risk which is minimized by the model (mainly proved in Lemma 2). The proofs of the techincal content of lemmas 1 and 2 rely both on computational geometry arguments and on existing results from [1].    Experiments are provided which show:  (1) an excellent match between the dependence of the bounds and that of the observed generalisation gaps on various parameters and (2) that the unlabeled data indeed improves performance in data-sparse regimes on some real data.     ======================= Reference: =======================   [1] Kai Zhong, Zhao Song, Prateek Jain, Peter L. Bartlett, Inderjit S. Dhillon. "Recovery Guarantees for One-hidden-layer Neural Networks", JMLR 2017 

The paper theoretically analyzes self-training in neural networks with one hidden layer. Self-training is a popular semisupervised learning algorithm where a model uses large unlabeled data by training on pseudolabels generated by a teacher model trained on some small set of labels. The paper studies the generalization performance of self-training on a two-layer network under a Gaussian assumption on the inputs, and has some synthetic and CIFAR-10 experiments to show that the predicted theoretical rates (1/\sqrt{ size unlabeled data}) roughly match empirical performance. 

This paper theoretically analyzes the iterative self-training algorithm with 1-hidden layer neural network. When the labeled and unlabeled data is generated from zero-mean Gaussian distributions with variances of different scales, it is shown that both the convergence rate and the generalization error decreases at a rate of $1/\sqrt{M}$ where $M$ is the number of unlabeled data. Experiments on synthetic data and augmented CIFAR-10 corroborate theoretical findings.