Year: 2021

Abstract:
We describe an algorithm that learns two-layer residual units using rectified linear unit (ReLU) activation: suppose the input x is from a distribution with support space R and the ground-truth generative model is a residual unit of this type, given by y = B∗ [ (A∗x) + + x ] , where ground-truth network parametersA∗ ∈ Rd×d represent a nonnegative full-rank matrix andB∗ ∈ Rm×d is full-rank with m ≥ d and for c ∈ R, [c]i = max{0, ci}. We design layer-wise objectives as functionals whose analytic minimizers express the exact ground-truth network in terms of its parameters and nonlinearities. Following this objective landscape, learning residual units from finite samples can be formulated using convex optimization of a nonparametric function: for each layer, we first formulate the corresponding empirical risk minimization (ERM) as a positive semi-definite quadratic program (QP), then we show the solution space of the QP can be equivalently determined by a set of linear inequalities, which can then be efficiently solved by linear programming (LP). We further prove the statistical strong consistency of our algorithm, and demonstrate its robustness and sample efficiency through experimental results.

text: [ (A∗x) + + x ] , where ground-truth network parametersA∗ ∈ Rd×d
represent a nonnegative full-rank matrix andB∗ ∈ Rm×d is full-rank with m ≥ d and for c ∈ Rd, [c+]i = max{0, ci}. We design layer-wise objectives as functionals whose analytic minimizers express the exact ground-truth network in terms of its parameters and nonlinearities. Following this objective landscape, learning residual units from finite samples can be formulated using convex optimization of a nonparametric function: for each layer, we first formulate the corresponding empirical risk minimization (ERM) as a positive semi-definite quadratic program (QP), then we show the solution space of the QP can be equivalently determined by a set of linear inequalities, which can then be efficiently solved by linear programming (LP). We further prove the statistical strong consistency of our algorithm, and demonstrate its robustness and sample efficiency through experimental results.

heading: 1 INTRODUCTION
text: Neural networks have achieved remarkable success in various fields such as computer vision (LeCun et al., 1998; Krizhevsky et al., 2012; He et al., 2016a) and natural language processing (Kim, 2014; Sutskever et al., 2014). This success is largely due to the strong expressive power of neural networks (Bengio & Delalleau, 2011), where nonlinear activation units, such as rectified linear units (ReLU) (Nair & Hinton, 2010) and hyperbolic tangents (tanh) play a vital role to ensure the large learning capacity of the networks (Maas et al., 2013). However, the nonlinearity of neural networks makes them significantly more difficult to train than linear models (Livni et al., 2014). Therefore, with the development of neural network applications, finding efficient algorithms with provable properties to train such nontrivial neural networks has become an important and a relatively new goal.
A
ReLU
B Residual networks, or ResNets (He et al., 2016a), are a class of deep neural networks that adopt skip connections to feed values between nonadjacent layers, where skipped layers may contain nonlinearities in between. Without loss of expressivity, ResNets avoid the vanishing gradient problem by directly passing gradient information from previous layers to current layers where otherwise gradients might vanish without skipping. In practice, ResNets have shown strong learning efficiency in several tasks, e.g. achieving at least 93% test accuracy on CIFAR-10 classification, lowering single-crop error to 20.1% on the 1000-class ImageNet dataset (Russakovsky et al., 2015; He et al., 2016b).
Common ResNets are often aggregated by many repeated shallow networks, where each network acts as a minimal unit with this kind of skip propagation, named a residual unit (He et al., 2016b). Given the flexibility and simplicity of residual units, much theoretical work has been devoted to study them and develop training algorithms for them in a way that sidesteps from the standard backpropagation regime and provides guarantees on the quality of estimation (see subsection 1.1). In this paper, we propose algorithms that can learn a general class of single-skip twolayer residual units with ReLU activation as shown on the right by the equation: y = B [ (Ax) + + x ] ,
where for a scalar c, c+ = max{0, c} (for a vector, this maximization is applied coordinate-wise), x
is a random vector as the network input with support space Rd, andA ∈ Rd×d andB ∈ Rm×d are weight matrices of layer 1 and layer 2, respectively.
Compared to previous work (Ge et al., 2018; 2017; Zhang et al., 2018; Wu et al., 2019; Tian, 2017; Du et al., 2017; Brutzkus & Globerson, 2017; Soltanolkotabi, 2017; Li & Yuan, 2017; Zhong et al., 2017), the introduction of residual connections simplifies the recovery of the network parameters by removing the permutation and scaling invariance. However, the naive mean square error minimization used in estimating the parameters remains nonconvex. Unlike most previous work, we do not assume a specific input distribution nor that the distribution is symmetric (Ge et al., 2018; Du & Goel, 2018).
We show that under regularity conditions on the weights of the residual unit, the problem of learning the unit can be formulated through quadratic programming (QP). We use nonparametric estimation (Guntuboyina et al., 2018) to estimate the ReLU function values in the networks. We further rewrite our constructed quadratic programs to linear programs (LPs). The LP formulation is simpler to optimize and has the same solution space as the QP for the network parameters.

heading: 1.1 RELATED WORK
text: Provable learning of neural networks has become an important topic of research. Arora et al. (2014) recover a multi-layer generative network with sparse connections and Livni et al. (2014) study the learning of multi-layer neural networks with polynomial activation. Goel et al. (2018) learns a one-layer convolution network with a perceptron-like rule. They prove the correctness of an iterative algorithm for exact recovery of the target network. Rabusseau et al. (2019) describe a spectral algorithm for two-layer linear networks. Recent work has connected optimization and two-layer network learning (Ergen & Pilanci, 2021; Sahiner et al., 2020; Ergen & Pilanci, 2020; Pilanci & Ergen, 2020) and has showed how to optimize networks layer by layer (Belilovsky et al., 2019).
For learning a one-layer ReLU network, Wu et al. (2019) optimize the norm and direction of neural network weight vectors separately, and Zhang et al. (2018) use gradient descent with a specific initialization. For learning a two-layer ReLU network, Ge et al. (2017) redesign the optimization landscape such that it is more amenable to theoretical analysis and Ge et al. (2018) use a momentbased method for estimating a neural network Janzamin et al. (2015). Many others have studied ReLU networks in various settings (Tian, 2017; Du et al., 2017; Brutzkus & Globerson, 2017; Soltanolkotabi, 2017; Li & Yuan, 2017; Zhong et al., 2017; Goel & Klivans, 2017).
The study of ReLU networks with two hidden layers has also been gaining attention (Goel & Klivans, 2017; Allen-Zhu et al., 2019), with a focus on PAC learning (Valiant, 1984). In relation to our work, Allen-Zhu & Li (2019) examined the PAC-learnable function of a specific three-layer neural network with residual connections. Their work differs from ours in two aspects. First, their learnable functions include a smaller (in comparison to the student network) three-layer residual network. Second, the assumptions they make on their three-layer model are rather different than ours.
In relation to nonparametric estimation, Guntuboyina et al. (2018) treat the final output of a shaperestricted regressor as a parameter, placing some restrictions on the type of function that can be estimated (such as convexity). They provide solutions for estimation with isotonic regression (Brunk, 1955; Ayer et al., 1955; van Eeden, 1956), convex regression (Seijo et al., 2011), shape-restricted additive models (Meyer, 2013; Chen & Samworth, 2016) and shape-restricted single index models (Kakade et al., 2011; Kuchibhotla et al., 2017).

heading: 1.2 MAIN RESULTS
text: We design quadratic objective functionals with a linear bounded feasible domain that take network parameters and functions as variables to estimate the ground-truth network parameters and nonlinearities. The values of the objectives are moments over the input distribution. Thm. 1.1 summarizes the landscapes of the objectives.
Theorem 1.1 (objective landscape, informal). Suppose a ground-truth residual unit has nondegenerate weights in both layers and nonnegative weights in layer 1. Then there exist quadratic functionals defined in linear-constrained domains whose minimizers a) are unique, and are the exact ground-truth, or b) are not unique, but can be adjusted to the exact ground-truth.
In practice, the exact moments over unobserved variables are not available. We can only construct the empirical risk minimization (ERM) of the moment-valued objectives by generated samples. With functions as variables in moment-valued objectives being optimized nonparametrically, the empirical objectives become quadratic functions with linear constraints (QP). We further show the convexity of our QP which guarantees its solution in polynomial time w.r.t. sample size and dimension. With the solution to the QP available, the strong consistency of our network learner is guaranteed: Theorem 1.2 (strong consistency, informal). Suppose samples are generated by a ground-truth residual unit that has nondegenerate weights in both layers and nonnegative weights in layer 1. Then there exists an algorithm that learns a network a.s.−−→ the exact ground-truth as sample size grows.
Roadmap: AssumeA andB are student network weights of layer 1 and 2. We first give a warm-up vanilla linear regression approach only knowing A is in Rd×d and B is in Rm×d (section 3). We then move to the details of our nonparametric learning for layer 2 (section 4), and similarly, layer 1 (section 5). We formalize Thm. 1.1, showing how nonparametric learning allows us to select the values of A and B from reduced spaces that are seeded by A∗ and B∗, under which LR is faster than vanilla LR. In section 6, we describe the strong consistency of our methods for respective layers, formalizing Thm. 1.2 using the continuous mapping theorem (Mann & Wald, 1943).

heading: 2 PRELIMINARIES
text: We describe the notations used in this paper, introduce the model and its underlying assumptions, and state conditions which simplify the problem but which can be removed without loss of learnability.
Notation The ReLU residual units we use take a vector x ∈ Rd as input and return a vector y ∈ Rm. We use A∗ ∈ Rd×d and B∗ ∈ Rm×d to denote the ground-truth network parameters for layer 1 and 2, respectively. We use circumflex to denote predicted terms (e.g. an empirical objective function f̂ , estimated layer 1 weights Â). We use n to denote the number of i.i.d. samples available for the training algorithm, {x(i),y(i)} to denote the i-th sample drawn. For an integer k, We define [k] to be {1, 2, . . . , k}, and e(j) as the standard basis vector with a 1 at position j. All scalar-based operators are element wise in the case of vectors or matrices unless specified otherwise. We use x and y to refer to the input and output vectors, respectively, as random variables.
Linear Regression: Linear regression (LR) in this paper refers to the problem of estimating L for unbiased model y = Lx. In this case, estimation is done by minimizing the empirical risk R̂ (L) = 12n ∑ i∈[n] ∥∥Lx(i) − y(i)∥∥2 – using linear least squares (LLS). Its solution has a closed form which we use as given. We refer the reader to Hamilton (1994) for more information.
Models, Assumptions and Conditions Following previous work (Livni et al., 2014), we assume a given neural network structure specifies a hypothesis class that contains all the networks conforming to this structure. Learning such a class means using training samples to find a set of weights such that the neural networks predictions generalize well to unseen samples, where both the training and the unseen samples are drawn from an unobserved ground-truth distribution. In this paper, the hypothesis class is given by ReLU residual units and we assume it has sufficient expressive power to fit the ground-truth model. More specifically, we discuss the realizable case of learning, in which the ground-truth model is set to be a residual unit taken from the hypothesis class with the form:
y = B∗ [ (A∗x) + + x ] , (1)
and is used to draw samples for learning. Unlike other multi-layer ReLU-affine models which do not apply skip connections, we cannot permute the weight matrices of the residual unit and retain the same function because of skip-adding x (it breaks symmetry). This helps us circumvent issues of identifiability1, and allows us to precisely estimate the ground-truth weight matricesA∗ andB∗.
Our general approach for residual unit layer 2 learns a scaled ground-truth weight matrix that also minimizes the layer 2 objective. The existence of such scaled equivalence of our layer 2 approach comes from what is defined below.
1Here, we are referring to the ability to identify the true model parameters using infinite samples.
Definition 2.1 (component-wise scale transformation). A matrix A ∈ Rd×d is said to be a scale transformation w.r.t. the j-th component if (Aj,:) > = Aj,j · e(j).
Additionally, estimation is more complex when the layer 2 weightsB∗ is nonsquare. For simplicity of our algorithm presentation for layer 2, we stick to Cond. 2.1 in the following sections2.
Condition 2.1 (layer 2 objective minimizer unique). A∗ is not a scale transformation w.r.t. any components andB∗ is a square matrix, i.e. m = d.

heading: 3 WARM-UP: VANILLA LINEAR REGRESSION
text: Consider a ground-truth two-layer residual unit. If we assume that the inputs only contain vectors with negative entries, i.e. x < 0, the effect of the ReLU function in the residual unit then disappears because of the nonnegativity ofA∗. The residual unit turns into linear model y = B∗x. Thus, direct LR on samples with negative inputs can learn the exact ground-truth layer 2 parameterB∗ with at least d samples, when formulating the LR as a solvable full-rank linear equation system.
On the contrary, if the inputs only contain vectors with positive entries, i.e. x > 0, all the neurons in the hidden layer are then activated by the ReLU and the nonlinearity is eliminated. The residual unit in this case turns into y = B∗ (A∗ + Id) x. Taking the value that left-multiplies x as a single weight matrix D∗, this is also a linear model. Direct LR on at least d samples with positive inputs by the residual unit can learn the exactD∗. Since we have the access to the exactB∗, solving for the exact A∗ can be accomplished through solving a full-rank linear equation systemB∗ · Ã = D∗, where the unique solution Ã = A∗ + Id.
While simple, this vanilla LR approach requires a large number of redundant samples, since sampled inputs usually have a small proportion of fully negative/positive vectors. Taking random input vectors i.i.d. with respect to each entry as an example, the probability of sampling a vector with all negative entries is pd−, where p− is the probability of sampling a negative vector entry, then the expected number of samples to get one negative vector is 1/pd−. Denoting p+ similarly, 1/p d + samples are expected for a positive vector. Besides, each LR step in this approach requires d such samples respectively to make the linear equation system full-rank, which implies the expected sample size to be d-exponential d · ( 1/pd− + 1/p d + ) . For other common random vectors like Gaussian samples, the proportions of fully negative/positive vectors in sampled inputs are also expected to decrease exponentially as d grows, as high-dimensional random vectors like Gaussian are essentially concentrated uniformly in a sphere (Johnstone, 2006). Technical and experimental details about sample size expectations and the vanilla LR algorithm are further discussed in Appendix F.

heading: 4 NONPARAMETRIC LEARNING: LAYER 2
text: We present how we learn a residual unit layer 2 under Cond. 2.1 (estimatingB∗): We first design an objective functional with the arguments being a matrix and a function. The objective uses expectation of a loss over the true distribution generating the data and is uniquely minimized by [B∗]−1 and a rectifier function (ReLU). We then formulate its ERM using nonparametric estimation as a standard convex QP, further simplified as an LP that has the same capability as the QP to learn layer 2.

heading: 4.1 OBJECTIVE DESIGN AND LANDSCAPE
text: Consider the formulation of a residual unit as in Eq. 1. It is possible to rewrite the model as equation: C∗y = (A∗x) + + x, where the output of the hidden neuron with skip addition is on both sides of the equation, and C∗B∗ = Id. We aim to estimate the inverse ofB∗ by matrix variable C and the nonlinearity x 7→ (A∗x)+ by a function variable h. The objective is formulated as risk functional by the L2 error between values respectively computed by C and h
G2 (C, h) = 1 2 Ex [ ‖h (x) + x−Cy‖2 ] , (2)
2In Appendix G, we show that estimation of B∗ remains solvable without satisfying Cond. 2.1.
where the estimator C ∈ Rd×m, the domain of h is the nonnegative3 continuous4 Rd → Rd function space, written as C0≥0 in shorthand. This objective is quadratic because the forward mapping x 7→ h(x) + x and the backward mapping y 7→ Cy are both linear w.r.t. C and h, and the two are linearly combined in a L2 norm. The objective in Eq. 2 is minimized by the ground-truth, i.e. C∗ and x 7→ (A∗x)+, is one of its minimizers. However, it is not simple to describe other variable values that minimize the objective if any exist. We give that detail under Cond. 2.1, the minimizer of G2 is unique, as the exact ground-truth in the given domain. Theorem 4.1 (objective minimizer, layer 2). Define G2(C, h) as Eq. 2, whereC ∈ Rd×m, h ∈ C0≥0. Then under Cond. 2.1, G2(C, h) reaches its zero minimum iff C = [B∗] −1 and h : x 7→ (A∗x)+.
Technical details are in Appendix H, where we use Lem. 4.2 (in subsection 4.3) to prove a more general theorem which does not require Cond. 2.1 and is sufficient for Thm. 4.1. In the next subsection, we construct the ERM of G2 and present our convex QP formulation.

heading: 4.2 ERM WITH NONPARAMETRIC ESTIMATION IS CONVEX QP
text: Consider the second layer objective (Eq. 2) with nonnegative continuous function space as the domain of h. We follow Vapnik (1992) and define its standard empirical risk functional:
Ĝ2(C, h) = 1
2n ∑ i∈[n] ∥∥∥h(x(i)) + x(i) −Cy(i)∥∥∥2 . (3) The function variable h ∈ C0≥0 can be optimized either parametrically or nonparametrically. If we parameterize h, and show that the nonlinearity w.r.t. its parameters would make Ĝ2 lose its quadratic form. Instead, we estimate h nonparametrically: for each sample input x(i), we introduce variables ξ(i) that estimates mapped values by h. This avoids introducing nonlinearity to the objective and keeps Ĝ2 quadratic. On the other hand, the domain of h, i.e. nonnegative continuous function space, turns into a set of linear inequalities as constraints when optimizing nonparametrically. In this sense, learning the second layer of the residual unit can be formulated as the following QP:
min C,Ξ
ĜNPE2 (C,Ξ) := 1
2n ∑ i∈[n] ∥∥∥ξ(i) + x(i) −Cy(i)∥∥∥2 , s.t. ξ(i) ≥ 0, ∀ i ∈ [n]. (4) where Ξ = {ξ(i)}ni=1 is the nonparametric estimator of x 7→ (A∗x) +.
Nonparametric Estimation Validation: A solution to the ERM with nonparametric estimation, i.e. the QP, is guaranteed to be sufficient for minimizing the standard empirical risk functional (Eq. 3). More specifically, assuming C and Ξ = {ξ(i)}ni=1 are a solution to layer 2 QP (Eq. 4), it is clear that C and h ∈ C0≥0 such that h(x(i)) = ξ(i) minimize the empirical risk functional Eq. 3. Conversely, a minimizer of the standard empirical risk Eq. 3, C and h, corresponds to a solution to layer 2 QP as we set ξ(i) = h(x(i)). Therefore, minimizing Ĝ2 and solving layer 2 QP are empirically equivalent.
Convexity: The convexity of the QP: Eq. 4 is also guaranteed. First, constraints are linear. Second, for each sample with index i ∈ [n], the L2 norm wraps linearity w.r.t. C and ξ(i). Such formulation ensures the quadratic coefficient matrix is positive semidefinite. Thus, with the sum of convex functions still being convex, the QP objective (Eq. 4) becomes convex. Even without the knowledge of how samples are generated, this QP would be a convex program. Strict proofs are in Appendix I.

heading: 4.3 LP SIMPLIFICATION
text: Consider single-sample error written as g2 (C, h;x,y) = 12 ‖h (x) + x−Cy‖ 2. If there is a feasible C such that Cy − x ≥ 0 holds for all x ∈ Rd, then C and h : x 7→ Cy − x always minimize g2 as zero, and thereby minimize the layer 2 objective (Eq. 2). Thus, we obtain a condition that is equivalent to G2 reaching the minimum in Thm. 4.1 and avoids randomness (see Lem. 4.2).
3Setting h as nonnegative ensures that a) only ReLU nonlinearity minimizes G2 (see Thm. 4.1), and b) h’s nonparametric estimator is linearly constrained (explained in subsection 4.2).
4If h is not a continuous function, only a null set of discontinuities is possible to make G2 reach zero as its minimum. Setting h as continuous simplifies our theoretical results that still strictly support empirical discussion.
Algorithm 1 Learn a ReLU residual unit, layer 2. 1: Input: {(x(i),y(i))}ni=1, samples drawn by
Eq. 1. 2: Output: B̂, Ξ̂, a layer 2 and x 7→ (A∗x)+
estimate. 3: Go to line 4 if QP, line 5 if LP. 4: Solve QP: Eq. 4 and obtain a ĜNPE2 minimizer,
denoted by Ĉ, Ξ̂. Go to line 6. 5: Solve LP: Eq. 5 and obtain a minimizer Ĉ,
then assign ξ̂(i) ← Ĉy(i) − x(i). 6: return Ĉ−1, Ξ̂.
Algorithm 2 Learn a ReLU residual unit, layer 1. 1: Input: {(x(i),h(i))}ni=1, layer 1 samples. 2: Output: Â, a layer 1 estimate. 3: Solve QP: Eq. 7 or LP: Eq. 8 and obtain a ĜNPE1 minimizer, denoted by Â, Φ̂. {Φ̂ is no longer needed.}
4: for all j ∈ [d] do 5: k̂j ← LR{h(i)j , Âj,:x(i)}h(i)j >0. 6: Âj,: ← Âj,:/k̂j . {Rescale Â.} 7: end for 8: return Â.
Lemma 4.2. G2(C, h) reaches its zero minimum iff Cy − x ≥ 0 holds for any x ∈ Rd and its corresponding residual unit output y, and h : x 7→ Cy − x.
The pointwise satisfaction of the inequality in Lem. 4.2 describes the solution space for the minimization of G2. The sufficiency of satisfying this inequality to minimize G2 is directly obtained by assigning h : x 7→ Cy − x where C complies with Cy − x ≥ 0 for all x in the support space Rd. The necessity of satisfying this inequality comes from its contraposition: If a C violates the inequality, there must be a non-null set of x that yield the violation due to the continuity of Cy − x. In this sense, the resulting G2 value becomes nonzero.
Empirically speaking, we cannot solve an inequality that holds w.r.t. to every point in the support space if we only observe finite samples. We can only estimate C by solving the inequality that holds w.r.t. each sample. Following this, we formulate such estimation as to find a feasible point in the space defined by a set of linear inequalities, each of which corresponds to a sample: Cy(i)−x(i) ≥ 0. Each point in the feasibility defined by the inequalities has a one-to-one correspondence in the solution space to layer 2 QP (Eq. 4): C ↔ (C, {Cy(i) −x(i)}i∈[n]). The set of inequalities can be solved by a standard LP with a constant objective and with constraints that are inequalities, i.e.
min C
const, s.t. Cy(i) − x(i) ≥ 0, ∀ i ∈ [n]. (5)
With the one-to-one correspondences, LP: Eq. 5 and QP: Eq. 4 have equivalent solution spaces. Moreover, LP: Eq. 5 is also a convex program. Alg. 1 summarizes the layer 2 estimator: Simply solve the QP/LP5 and return the inverse of Ĉ as the estimate of layer 2 weights and Ξ̂ as x 7→ (A∗x)+ estimate. Regardless of time complexity, QP and LP in Alg. 1 work equivalently since their solution spaces are equivalent to each other.
Our nonparametric learning directly finds a unique layer 2 estimate under Cond. 2.1. In the general case without Cond. 2.1 (discussed in Appendix G), nonparametric learning essentially reduces the possible values ofB from Rm×d to aB∗ scale-equivalent matrix space, where LR uses sampled data much more efficiently than vanilla LR on layer 2 (section 3).

heading: 5 NONPARAMETRIC LEARNING: LAYER 1
text: With layer 2 learned, outputs by hidden neurons become observable. The two-layer problem is thereby reduced to single-layer. Consider a ground-truth single-layer model: h = (A∗x)+. To construct a learning objective for this model, we rewrite the model as a nonlinearity plus a linear mapping byA∗: h = (−A∗x)+ +A∗x, where on both sides of the equation is the output of layer 1, and the nonlinearity is x 7→ (−A∗x)+. The objective of layer 1 is formulated as:
G1 (A, r) = 1 2 Ex [ ‖r(x) +Ax− h‖2 ] , (6)
5We use CVX (Grant & Boyd, 2014; 2008) that calls SDPT3 (Toh et al., 1999) (a free solver under GPLv3 license) and solves our convex QP/LP in polynomial time. See Appendix C for technical details.
whereA ∈ Rd×d is of the layer 1 weights estimator, the domain of r is also C0≥0. The minimizerA of the risk G1 falls into a matrix space such that for any matrixA in the space, each row ofA is a scaled-down version of the same row ofA∗ without changing the direction (Thm. 5.1). Theorem 5.1 (objective minimizer space, layer 1). Define G1(A, r) as Eq. 6, where A ∈ Rd×d, r ∈ C0≥0. Then G1(A, r) reaches its zero minimum iff for each j ∈ [d], Aj,: = kjA∗j,: where 0 ≤ kj ≤ 1, and r : x 7→ (A∗x)+ − diag(k) ·A∗x.
The scale equivalence in the solution space is derived from a ReLU inequality: (x)+ ≥ kx where 0 ≤ k ≤ 1. Let the j-th row of A be a scaled-down version of A∗, i.e. Aj,: = kjA∗j,: where 0 ≤ kj ≤ 1. According to the inequality, we have ( A∗j,:x
)+ ≥ kjA∗j,:x, which indicates (A∗x)+ ≥ diag(k) ·A∗x. Thus, r minimizing G1 in Thm. 5.1 lies in feasibility C0≥0 whenA = diag(k) ·A∗. Due to the existence of scale equivalence, we must compute the scale factor to obtain the ground-truth weightsA∗. The scale factor kj is sufficiently obtainable with (Aj,:x) + and ( A∗j,:x )+ observable: Conditioned on nonnegative ReLU input, we have a linear model Aj,:x = kjA∗j,:x where Aj,:x andA∗j,:x are observed. Thm. 5.2 summarizes layer 1 scale factor property, allowing us to correct a minimizer of G1 to the ground-truth weightsA∗ by computing a scalar for each row. Theorem 5.2 (scale factor, layer 1). Assume A is a minimizer of G1. Then for any j ∈ [d], (Aj,:x) + / ( A∗j,:x )+ is always equal to the scale factor kj given that (Aj,:x) + > 0.
The elimination of randomness for G1 follows the same pattern as layer 2. If there is a feasible A such that (A∗x)+ −Ax ≥ 0 holds for all x ∈ Rd, such A is a solution to our objective by Eq. 6. We can also have the following proposition that avoids randomness. Lemma 5.3. G1(A, r) reaches its zero minimum iff h − Ax ≥ 0 holds for any x ∈ Rd and its corresponding hidden output h, and r : x 7→ h−Ax.
We now turn into empirical discussion. Similar to layer 2, we formulate layer 1 QP by its empirical objective with nonparametric estimation and linear constraints representing the nonnegativity of r:
min A,Φ
ĜNPE1 (A,Φ) := 1
2n ∑ i∈[n] ∥∥∥φ(i) +Ax(i) − h(i)∥∥∥2 , s.t. φ(i) ≥ 0, ∀ i ∈ [n]. (7) where Φ = {φ(i)}ni=1 is the function estimator of x 7→ (A∗x)
+ − diag(k) ·A∗x where k refers to the scaling equivalence. Similarly, the solution space of QP: Eq. 7 can be represented by a set of linear inequalities as constraints of the following efficiently solvable LP
min A
const, s.t. h(i) −Ax(i) ≥ 0, ∀ i ∈ [n]. (8)
Alg. 2 describes how a layer 1 is learned: First, a G1 minimizer estimate Â is obtained by solving the QP/LP. Then for the j-th row, the scale factor kj is estimated by running LR on h (i) j and Âj,:x (i) s.t. h(i)j > 0 to correct Â, as Aj,:x = kjhj is an unbiased linear model given that hj > 0. By nonparametric learning, the value space ofA is reduced from Rd×d to {diag(k) ·A∗ | 0 ≤ kj ≤ 1}, where LR uses sampled data much more efficiently than vanilla LR layer 1 (section 3).

heading: 6 FULL ALGORITHM AND ANALYSIS
text: The full algorithm concatenates Alg. 1 and 2 in a layerwise fashion, with observations of input/output by the ground-truth network: a) Estimates layer 2 and nonlinearity: x 7→ (A∗x)+ by Alg. 1. b) Estimates layer 1 by running Alg. 2 on input samples and the nonlinearity estimate. It has provable guarantees. For empirical analysis, we use Ĉn to denote the estimation ofC∗ from n random samples. Similar notations are applied to other estimations. First of all, our methods to solve respective layers, Alg. 1 and 2, are strongly consistent if any convex QPs/LPs involved can be solved exactly.
Lemma 6.1 (layer 2 strong consistency). Under Cond. 2.1, Ĉn a.s.−−→ C∗ and B̂n a.s.−−→ B∗, n→∞.
For Ĉn: In Appendix J, we prove its more general a.s. convergence without satisfying Cond. 2.1, where the solution space to layer 2 objective is a non-compact continuous set where all the elements
are scaling equivalences. We use Hausdorff distance (Rockafellar & Wets, 2009) as metric and prove that the empirical solution space a.s. converges to the theoretical solution space. Then the more general a.s. convergence holds with the strong consistency of layer 2 scale factor estimator where we use LR to estimate the scale factors.
For B̂n: Using the continuous mapping theorem (Mann & Wald, 1943), we directly propagate the strong consistency of C∗ estimator to its inverseB∗’s estimator.
According to full algorithm description, layer 1 estimation uses Ĉny(i) − x(i) as the outputs where i ∈ [n]. Thus, the strong consistency of the hidden neuron estimator is also guaranteed by the continuous mapping theorem. Following the proof sketch of the Ĉn a.s. convergence, we obtain the strong consistency of layer 1 estimator (See Lem. 6.2).
Lemma 6.2 (layer 1 strong consistency). Ân a.s.−−→ A∗, n→∞.
By the continuous mapping theorem, the strong consistency of the full algorithm (Thm. 6.3, formal version of Thm. 1.2), which is commonly defined by a loss function that is continuous on network weights, is implied by the strong consistency of network weights estimators (Lem. 6.1 and 6.2). See Appendix J for proofs of strong consistency discussions in this section.
Theorem 6.3 (strong consistency, formal). Define L as L2 output loss. L(Ân, B̂n) a.s.−−→ 0, n→∞.

heading: 7 EXPERIMENTS
text: We provide experimental analysis to demonstrate the effectiveness and robustness of our approach in comparison to stochastic gradient descent (SGD) on L2 output loss: L(A,B) = 12Ex ‖ŷ − y‖ 2 ,
where we parameterize the output prediction by ŷ = B [ (Ax) + + x ] . Our proposed methods
outperform SGD in terms of sample efficiency and robustness to different network weights and noise strengths, which indicates a poor optimization landscape of L2 output loss for ReLU residual units.
Setup: The ground-truth weights are generated through i.i.d. folded standard Gaussian6 and standard Gaussian for layer 1 and 2 respectively, i.e. A∗ i.i.d.∼ |N | (0, 1), B∗ i.i.d.∼ N (0, 1). The input distribution is set to be an i.i.d. zero mean Gaussian-uniform equal mixture N (−0.1, 1) – U (−0.9, 1.1). SGD is conducted on mini-batch empirical losses of L(A,B) with batch size 32 for 256 epochs in each learning trial. We apply time-based learning rate decay η = η0/ (1 + γ · T ) with initial rate η0 = 10
−3 and decay rate γ = 10−5, where T is the epoch number. The above hyperparameters are tuned to outperform other hyperparameters in learning ReLU residual units in terms of output errors.
Evaluation: We use relative errors to measure the accuracy of our vector/matrix estimates: For a network with weights A and B and its teacher network with weights A∗ and B∗, a) layer 1 error refers to ‖A−A∗‖ / ‖A∗‖, similar to layer 2. b) output error refers to Ê [‖ŷ − y‖ / ‖y‖] by test data. Due to the equivalence between the solution spaces of our QP and LP without label noise, we choose LP in noiseless experiments, referred to as “ours”. In addition, to reduce variance, the results of learning the same ground-truths are computed as means across 16 trials.
Sample Efficiency: Consider Fig. 1. It depicts the variation in the prediction errors (warmer color, larger error) as a function of d (input dimension) and the number of samples the learning algorithm is using. We compare SGD against our algorithm. We observe that our approach to the estimating the neural network is more sample efficient. For SGD, the estimation is relatively easy with only up to 10 dimensions. As expected, once the dimension grows, the sample size required for the same level of error as our method is larger. Still, overall, our method is capable of learning robustly with small sample sizes and more efficiently than SGD even for larger sample sizes.
Network Weight Robustness: This experiment aims to verify whether our method can learn a broader class of residual units. In Tab. 1, our method shows a light-tailed distribution with nearly zero means and standard deviations for layers 1 and 2, and output errors across various ground-truth networks, whereas SGD is less robust in the same context. Our method shows strong robustness to network weight changes, indicating its applicability across the whole hypothesis class.
6A folded Gaussian is the absolute value of a Gaussian, with p.d.f. p(|x|) where x ∼ N , denoted as |N |. We use folded Gaussian to ensure layer 1 weights are nonnegative.
SGD
160 192 224 256 288 320 352 384
8
10
12
14
16
18
20
22
Ours
160 192 224 256 288 320 352 384
8
10
12
14
16
18
20
22
0.1
0.2
0.3
0.4
0.5
0.6
Figure 1: Output errors by SGD and our method for different dimensions and sample sizes. For the same d, we fix the ground-truth network as sample size grows.
Table 1: Means / Standard deviations of network estimate errors for different network weights. Values are computed from the process of learning 128 different ground-truth networks with d = 16. 512 training samples are drawn for each learning trial.
Layer 1 Layer 2 Output Mean Std Mean Std Mean Std
SGD 0.715 0.090 1.203 0.134 0.431 0.038 Ours 0.039 0.008 ≈ 0 ≈ 0 0.055 0.008
0 0.2 0.4 0.6 0.8 1
0
0.1
0.2
0.3
0.4
0.5
SGD LP QP
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
SGD LP QP
0 0.2 0.4 0.6 0.8 1
0
0.2
0.4
0.6
0.8
1
1.2
1.4
SGD LP QP
Figure 2: Respective errors of layers 1 and 2 and outputs of different label noise strengths for SGD, LP and QP. We fix the ground-truth weights with d = 10 and only the noise strength varies. 512 training samples are drawn for each learning trial.
Noise Robustness: Fig. 2 confirms the robustness of our methods when output noise exists. Samples are generated by a ground-truth residual unit with output noise being i.i.d. zero-mean Gaussian in different strengths (i.e. standard deviations). We try both QP and LP because in noisy setting the two approaches are not equivalent w.r.t. the solution space7. First, SGD always gives larger errors than our methods, even though it is hardly affected by tuning the noise strength. For QP/LP, all the errors for layer 1, 2 and output grow almost linearly as noise strength increases, indicating that both QP/LP learn the optima robustly when output noise is present, where QP slightly outperforms LP.

heading: 8 CONCLUSION
text: In this paper, we address the problem of learning a general class of two-layer residual units and propose an algorithm based on landscape design and convex optimization: We demonstrate firstly that minimizers of our objective functionals can express the exact ground-truth network. Then, we show that the corresponding ERM with nonparametric function estimation can be solved using convex QP/LP, which indicates polynomial-time solvability w.r.t. sample size and dimension. Moreover, our algorithms that are used to estimate both layers as well as the whole networks are strongly consistent, with very weak conditions on input distributions.
Our work opens the door to a variety of open problems to explore. We provide a strong consistency result without assuming a particular input distribution. It would be interesting to explore the sample complexity of our algorithm when stronger assumptions are placed on input distribution. Another extension could be to dispose of the limitations on ground-truth weights, e.g. m ≥ d, non-degeneration, nonnegativity ofA∗. In addition, given that our algorithm solves ReLU by using nonparametric estimation through convex optimization, this might provide an inspiration to solving learning problems with other nonlinearities using nonparametric methods that improve or change the optimization landscape.
7See Appendix B for noisy model discussion.

heading: REPRODUCIBILITY STATEMENT
text: The full set of assumptions of all theoretical results in this paper are stated - section 2, Appendix F and G. The complete proofs of all theoretical results are all included - Appendix F, G, H, I and J. Code in MATLAB is zipped as supplementary file. The numerical simulation details, including data sizes and splits, hyperparameters and how they were chosen, clear definitions of the specific measures or statistics used to report results, are discussed in detail - the Setup and Evaluation paragraphs in section 7. Results with central tendency (e.g. mean) and variation (e.g. stddev) are reported - Tab. 1. The use of external computing infrastructure, CVX, is clarified, and the consent obtaining and licensing of CVX and SDPT3 solver used in this work exactly follows the original author instructions: http://cvxr.com/cvx/doc/citing.html - footnote 5.

heading: APPENDICES
text: We include an outline of the appendices:
• Appendix A: An empirical analysis of the running time of our algorithm and an analysis of the effect of condition number on the estimation accuracy.
• Appendix B: A brief discussion of our learning algorithms in a noisy context, referring to the LP slack variable technique that is used in our experiments with noise in the main paper. Also, this section provides an insight to potential future direction of this work.
• Appendix C: A detailed explanation of the computational complexity and the methods through which the QPs/LPs in this paper are solved.
• Appendix D: A preliminarily proposed general optimization formula for the extension of learning multi-layer ResNets, to show the scalability of our algorithm.
• Appendix E: A justification of how our QP/LP approaches get rid of the exponentiality that vanilla LR approach (section 3) has.
• Appendix F: A formal result and an empirical validation of the example given in the main paper in the context of the vanilla linear regression method (section 3).
• Appendix G: A generalization of layer 2 learning to the case where Cond. 2.1 is not satisfied.
• Appendix H: Proofs regarding the minimizers of the objective functions we use.
• Appendix I: Proofs that justify the convexity of our QPs.
• Appendix J: Proofs that show our estimation algorithm is strongly consistent.

heading: A RUNNING TIME AND SENSITIVITY TO CONDITION NUMBER
text: Following the setup in section 7, we report additional experiments testing the running time of our algorithm in comparison to SGD and the feasibility of estimating the inverse of B∗ when ill-conditioned.

heading: A.1 RUNNING TIME EFFICIENCY VERSUS SGD
text: In Fig. 3, we compare the running time of our algorithm against the running time of SGD for ground-truth networks sampled at random. The running time of our algorithm is significantly lower, with an error level which is also significantly lower. Furthermore, this is demonstrated across input dimensions, where the difference between the running time of SGD until convergence and our algorithm’s running time increases as d increases. In addition, the running time of our algorithm is fixed, and does not depend on the ground-truth network, while the running time of SGD varies based on the network sampled.

heading: A.2 LAYER 2 WEIGHTS CONDITION NUMBER ROBUSTNESS
text: Our algorithm depends on a matrix-inversion step to obtain the estimated parameters. This affects the estimation ofB∗. Due to this, we further investigate the effect of the condition number ofB∗ in the ground-truth parameters on the accuracy of estimation. For this experiment, to obtain ground-truth network parameters with different condition number, we follow the procedure of Ge et al. (2018), and multiply a diagonal matrix with exponentially-dropping values on the diagonal (λ−i for the i-th element) by two random orthonormal matrices, U and V>, on the left and on the right respectively: B = U diag ( λ−1, . . . , λ−d ) V>.
In Fig. 4, we compare the effect of the condition number of B∗ on SGD and on our algorithm. SGD turns out to be quite sensitive to the condition number, especially for lower dimensions. Our algorithm, on the other hand, is almost not affected by the condition number, with a final estimation error consistently close to 0, recovering the matrixB∗.

heading: B DISCUSSION: NOISY MODEL
text: Here, we discuss our learning methods in the noisy case. We introduce output noise to our model, namely
y = B∗ [ (A∗x) + + x ] + z, (9)
where the label noise z ∈ Rd is an i.i.d. random vector with respect to each component, satisfying E [z] = 0. In addition z and x are statistically independent.
Taking layer 2 as an example, our original objective functional does not reach zero by substituting the ground-truth: G2(C∗,x 7→ (A∗x)+) = Ez[‖C∗z‖2] = σ2 TrC∗C∗> where σ is the noise strength. However, if layer 2 is well conditioned, the ground-truth will still assign a value close to zero to the objective. In this sense, with G2’s continuity, the ground-truth can approximately minimize G2, which validates our QP approach in terms of learning noisy residual units.
Our original LP (Eq. 5) fails to give feasible solutions due to possible violation of the inequality C∗y−x ≥ 0, sinceC∗y−x = (A∗x)+ +C∗z is not necessarily an entrywise nonnegative vector because the term C∗z might have negative entries. So we introduce slack variables ζ(i) to soften the constraints:
min C,Z
1
n ∑ i∈[n] 1> · ζ(i), (10)
s.t. Cy(i) − x(i) ≥ −ζ(i), ζ(i) ≥ 0. (11)
For a sample (x(i),y(i)) with noise z(i), if ( A∗x(i) )+ +C∗z(i) < 0, then its L1 norm would be added to the objective. This method remedies violations to the inequality C∗y − x ≥ 0. With access to a sufficiently large sample and with the “stability” assumption, our solution Ĉ would be close to C∗ since large deviations seldom rise and their penalties are diluted in the objective.

heading: C DISCUSSION: SOLVING CONVEX PROGRAMS
text: The theoretical foundation of solving convex QP and LP has been driven to maturity in terms of computational complexity (Kozlov et al., 1980; Murty, 1980) and convergence analysis (Jarre, 1990; Zhang et al., 1992). The time complexity for a convex QP/LP is analyzed in terms of the number of scalar variables N and the number of bits L in the input (Ye & Tse, 1989). For example, under Cond. 2.1, N = d2 + nd for the QP because the matrix variables have d2 scalars and the function estimators have nd scalars. Similarly we have N = d2 for the noiseless LP and N = d2 + nd for the noisy LP. It is guaranteed that primal-dual interior point methods can solve the convex QP/LP in a polynomial number of iterations O( √ NL), where each iteration costs at worst O(N2.5) arithmetic operations from the Cholesky decomposition for the needed matrix inversion (Monteiro & Adler, 1989; Karmarkar, 1984). This indicates that our QPs/LPs are guaranteed to be solvable in O(N3L) arithmetic operations, which is generally an O(poly(n,m, d, L)) complexity. To solve convex QPs/LPs in this paper, we use CVX, a commonly used package for specifying and solving convex programs (Grant & Boyd, 2014; 2008). For both QPs and LPs, CVX calls a solver, SDPT3 (Toh et al., 1999), which is specified for semidefinite-quadratic-linear programming and applies interior-point methods with the computational complexity mentioned above. Experimentally, SDPT3 indeed specifies and solves our programs fast and makes our numerical results robust and stable.

heading: D DISCUSSION: EXTENSION TO MULTI-LAYER
text: In this appendix, we propose a generalized optimization formula of nonparametric learning of multilayer residual networks to show that our algorithm is scalable. Assume we have an L-layer network of the form
r0 = x rl = ( W ∗l−1rl−1 )+ + rl−1, l = 1, . . . , L− 1
y = W ∗L−1rL−1.
We can scale the objective to use nonparameter variables {ξ(i)l }l∈[L−1] to estimate every ReLU activation (W ∗l rl)
+ in this network. Generally, for the l-th ReLU we can have the outer layer objective
min ξ (i)
[l] ,Vl ∥∥∥ξ(i)l + · · ·+ ξ(i)1 + x(i) − Vlĥ(i)l ∥∥∥ , s.t. ξ(i)[l] ≥ 0 where hat denotes estimated values and ĥ(i)L−1 = y (i), and its corresponding inner layer objective
min φ
(i) l ,Wl−1,ξ (i) l ∥∥∥φ(i)l +Wl−1 (ξ̂(i)l−1 + · · ·+ ξ̂(i)1 + x(i))− ξ(i)l ∥∥∥ , s.t. φ(i)l ≥ 0, ξ(i)l ≥ 0. We can apply techniques, e.g., tune the combinations of those objectives with respect to layers or whether to have a hat or not, or even get them weighted, to keep its convex quadratic form and at the same time improve its optimization landscape. Specifically when L = 2, this becomes the two-layer ResNet learning discussed in the main paper.

heading: E SAMPLE EFFICIENCY: VANILLA LR VERSUS OURS
text: In this appendix, we describe an example to justify how our core QP/LP approaches eliminate the exponential sample efficiency compared to the vanilla LR approach. In layer 2 learning, C’s optimization is to find a feasible point in the space determined by n inequalities each of which corresponds to a sample. Taking the j-th row of the inequality, i.e. Cj,:y ≥ xj . Every time we get a new sample x(i), y(i) in,
• The inequality Cj,:y(i) ≥ x(i)j eliminates the solution space of Cj,: by one of the spaces divided by plane Cj,:y(i) = x (i) j . This property guarantees fast convergence speed at early
phase since Cj,:’s feasibility starts from Rd.
• In fact, when ReLU does not activateA∗x(i) at the j-th row, i.e.A∗x(i) ≤ 0, the theoretical solution of Cj,: that Cj,:B∗ = (0, . . . , 0, 1, 0, . . . , 0) (a 1 at index j) directly lies on the divisive plane Cj,:y(i) = x (i) j :
Cj,:y (i) = Cj,:B ∗ [( A∗x(i) )+ + x(i) ] = x (i) j .
This property remarkably speeds up further constraints on the solution space of Cj,: to the correct estimate and is with high probability, since it directly depends on the sign ofA∗x(i).
With the vanilla LR approach mentioned in section 3 all the dimensions need to have the correct sign at the same time. This is the reason for the exponential complexity of this approach. Our main algorithm only requires one dimension to get the correct sign for a single sample, which avoids the exponential sample size expectation. The convergence speed is determined by the actual probability of each dimension not getting activated by ReLU, and such probability is determined by specified input distribution.
Algorithm 3 Learn a ReLU residual unit by LR. 1: Input: {(x(i),y(i))}ni=1, n samples drawn by Eq. 1. 2: Output: Â, B̂, estimated weight matrices. 3: B̂ ← LR{(x,y) ∈ {(x(i),y(i))}bn/2ci=1 | x < 0}. 4: D̂ ← LR{(x,y) ∈ {(x(i),y(i))}ni=bn/2c+1 | x > 0}. {the estimation ofB ∗ (A∗ + Id)}
5: Solve full-rank linear equation system B̂ · Ã = D̂ 6: return Ã− Id, B̂.
F VANILLA LINEAR REGRESSION: MORE DETAILS
Alg. 3 gives the full list of steps of learning two-layer residual units by LR, where we first split the drawn samples into two halves for the respective two key steps, then for both halves we filter negative and positive vectors, respectively. By running LR on both filtered sets we obtain an estimation of the ground-truth network.
In the following sections, we formalize this intuition and describe the exponential sample complexity of this approach under entry-wise i.i.d. setting as an example of the inefficiency of this method. Theorem F.1 (exponential sample complexity, vanilla LR). Assume the input vectors are i.i.d. with respect to each component. Then Alg. 3 learns a neural network from a ground-truth residual unit with at least O ( d · 2d+1 ) expected number of samples.
Proof. For each j ∈ [d], let Pj be the marginal probability of xj being positive, i.e. Pj = P (xj > 0) > 0. Thus, the probability that a sample is positive P (x > 0) = ∏ j∈[d] Pj . Then the expected
number of sampling trials to obtain d positive samples is d/ ∏
j∈[d] Pj . Similarly, the expected number of sampling trials to obtain d negative samples is d/ ∏ j∈[d] [1− Pj − P (xj = 0)]. Let n be a random natural number s.t. with n samples the network has learned. The expected number of samples that guarantees successful learning is the sum of the two expectations
E[n] = d { 1∏
j∈[d] Pj + 1∏ j∈[d] [1− Pj − P (xj = 0)]
}
≥ d [ 1∏
j∈[d] Pj + 1∏ j∈[d] (1− Pj)
]
≥ d · 2 ∏ j∈[d] 1√ Pj (1− Pj)
≥ d · 2 ∏ j∈[d]
1
(Pj + 1− Pj) /2
= d · 2d+1.
Thus E[n] ≥ O(d · 2d+1).
With an exponential sample complexity, the time complexity of the LR approach is thereby also exponential because filtering through all the samples costs time with the same complexity as the number of samples, even though the LR itself only costs polynomial time in the number of samples. To summarize, this vanilla LR approach learns exact ground-truth residual units but with exponential complexity in terms of both computational cost and sample size. While this approach to learning a residual unit such as ours is simple and intuitive, we aim to improve on this approach, making full use of all samples available.

heading: F.1 EXPERIMENTS: SAMPLE EFFICIENCY
text: We present the results of the vanilla LR in learning the residual units. As discussed in section 3, LR requires full-rank linear systems parameterized by samples to learn the exact ground-truth parameters.
However, with degenerate linear systems, LR is completely incapable of learning the parameters. Therefore, there exists a hard threshold for the number of samples required by LR that makes the linear equation system full-rank. With such a dichotomous constraint on LR learning, we take the learning success rate among 1000 trials as the metric to evaluate the performance of LR with different sample sizes and number of dimensions.
Tab. 2 shows the learning success rates with zero-mean Gaussian inputs. For each fixed number of dimensions, it appears there is a hard threshold that switches the learnability of LR. The exponential sample complexity is also reflected there as the number of dimensions grows linearly. Tab. 3 shows the rates with input mean non-zero, where an overall decline in success rates happens. This observation is explainable because the bottleneck of LR learning is the lower value found between the probability of sampling a positive and a negative vector. A positive mean reduces the probability of the latter, and thereby increases the sample size required.

heading: G GENERALIZATION OF LAYER 2 LEARNING
text: In this appendix, we discuss how layer 2 is learned without satisfying Cond. 2.1. It is a fact that A∗ being a scale transformation w.r.t. some components causes scaling equivalence to our layer 2 objective functional minimizers just as in layer 1. To be more specific, we give a general version of Thm. 4.1 that handles the case without Cond. 2.1 and describes the scaling equivalence of G2 minimizers. See Thm. G.1 for a formal description.
Theorem G.1 (objective minimizer space, layer 2, general). Let G2 be 12Ex [ ‖h (x) + x−Cy‖2 ] as a functional, where C ∈ Rd×m, h ∈ C0≥0. Then G2 (C, h) reaches its zero minimum iff CB∗ = diag(k) where for each j ∈ [d],
a) 11+A∗j,j ≤ kj ≤ 1, ifA ∗ is a scale transformation w.r.t. the j-th component. b) kj = 1, ifA∗ is not a scale transformation w.r.t. the j-th component.
As in layer 1, the scaling equivalence in layer 2 can also be obtained by assigning A∗ as a scale transformation w.r.t. some component and using the properties of the ReLU nonlinearity. See subsection H.1 for a detailed explanation of Thm. G.1. To obtain the scale factors k and correct a scale-equivalent G2 minimizerC to the ground-truth, we observe linear models parameterized by the scale factors, where for each j ∈ [d], [Cy]j = kjxj given that xj < 0. See Thm. G.2 and its proof for justifications of the linear models that are used to compute k.
Theorem G.2 (scale factor, layer 2, general). Assume C is a minimizer of G1 in the context of Thm. G.1. Then for any j ∈ [d], the following three propositions are equivalent:
a) A∗ is a scale transformation w.r.t. the j-th component. b) [Cy]j /xj is a constant c n j given that xj < 0. c) [Cy]j /xj is a constant c p j given that xj > 0.
Additionally, if one of the above three propositions is true, then kj = cnj . Proof. For j ∈ [d], we first prove a) =⇒ b), c): From ( A∗j,: )> = A∗j,j · e(j) we have
[Cy]j xj =
[ kje (j) ]> [ (A∗x) + + x ]
xj
= kj
[( A∗j,:x )+ + xj ] xj = kj [( A∗j,jxj )+ + xj ] xj
= { kj , xj < 0 kj ( A∗j,j + 1 ) , xj > 0 .
Then we prove ¬ a) =⇒ ¬ b), ¬ c): ¬ a) =⇒ ∃ j′ ∈ [d] and j′ 6= j s.t. A∗j,j′ 6= 0. Recall
[Cy]j xj = kj
[( A∗j,:x )+ + xj ] xj .
The value of xj′ affects the value of [Cy]j /xj because A ∗ j,j′ 6= 0. In fact, from supp p(x) = Rd we have supp p(xj′ | xj) = R, indicating that the value of [Cy]j /xj can never be kept as a constant when given both xj < 0 and xj > 0 because xj′ can be any real number.
With the exact derivations above, we are able to obtain a left inverse ofB∗, namely C, that satisfies CB∗ = Id. Consider Eq. 1 left multiplied byB∗C
B∗Cy = y. (12)
Eq. 12 is also a noiseless unbiased linear model where Cy and y are observable, and therebyB∗ is computable due to the easy solvability of its linearity.
Alg. 5 describes how a residual unit second layer is learned without satisfying Cond. 2.1: We first solve the QP/LP and obtain a scaling equivalence to an estimated left inverse ofB∗ and the function estimate, namely Ĉ and Ξ̂. Then we compute the scale factor estimate k̂ by running Alg. 4, where for each component index j ∈ [d], we first use a tolerance parameter as a threshold to determine whether ground-truth layer 1 is a scale transformation, and if so, we run LR to estimate the model [Cy]j = kjxj , otherwise the scale factor is directly assigned by 1. Upon correcting Ĉ and Ξ̂ by k̂, we obtain a layer 2 estimate B̂ by running LR to estimate the linear model Eq. 12. The strong consistency of our results in this appendix is justified in Appendix J.

heading: H EXACT DERIVATION OF OBJECTIVE FUNCTIONAL MINIMIZERS

heading: H.1 LAYER 2
text: Lem. 4.2 is proved as follows.
Algorithm 4 Rescale a ĜNPE2 minimizer. 1: Parameters: εtol > 0, LR objective tolerance. 2: Input: {(x(i),y(i))}ni=1, samples drawn by Eq. 1; Ĉ, a ĜNPE2 minimizer. 3: Output: k̂, a layer 2 scale factor estimate. 4: for each j ∈ [d] do 5: k̂j ← LR { x (i) j , [ Ĉy(i) ] j } x (i) j <0 .
6: if the LR objective optimal R̂j(k̂j) > εtol then 7: k̂j ← 1. 8: end if 9: end for
10: return k̂.
Algorithm 5 Learn a ReLU residual unit layer 2. 1: Input: {(x(i),y(i))}ni=1, samples drawn by Eq. 1. 2: Output: B̂, Ξ̂, a layer 2 and x 7→ (A∗x)+ estimate. 3: Go to line 4 if QP, line 5 if LP. 4: Solve QP: Eq. 4 and obtain a ĜNPE2 minimizer, denoted by Ĉ, Ξ̂. Go to line 6. 5: Solve LP: Eq. 5 and obtain a minimizer Ĉ, then assign ξ̂(i) ← Ĉy(i) − x(i) for each i ∈ [n]. 6: Run Alg. 4 on {(x(i),y(i))}i∈[n], Ĉ and obtain k̂. 7: B̂ ← LR { diag−1(k̂) · Ĉy(i),y(i) } i∈[n] .
8: ξ̂(i) ← diag−1(k) [ ξ̂(i) + x(i) ] − x(i) for each i ∈ [n]. {Correct Ξ̂ to the function estimation
of (A∗x)+.} 9: return B̂, Ξ̂.
Proof. “⇐= ”: Since h(x) = Cy − x ≥ 0, random vector h(x) + x−Cy is always a zero vector, which implies G2(C, h) = 0. Hence “⇐= ” holds.
“ =⇒ ”: Since r.v. ‖h(x) + x−Cy‖2 ≥ 0 we have G2(C, h) = 0 =⇒ λ ({ x ∈ Rd ∣∣∣ ‖h(x) + x−Cy‖2 > 0}) = 0 (13) where λ is the Lebesgue measure on Rd.
Proof by contradiction: Assume that ∃x′ ∈ Rd, ∃ i ∈ [d], s.t. (Cy′ − x′)i < 0, where y′ is the corresponding network output. Let f(x) = (Cy − x)i which is continuous on Rd. Therefore for = −f(x′) > 0, ∃ δ > 0, ∀x ∈ B(x; δ) i.e. ‖x− x′‖ < δ,
|f(x)− f(x′)| < =⇒ 2f(x′) < f(x) < 0 =⇒ [x−Cy]i > 0 =⇒ [h(x) + x−Cy]i > 0 =⇒ ‖h(x) + x−Cy‖ 2 > 0.
Since λ (B(x; δ)) = πd/2
Γ (d/2 + 1) δd > 0 and B(x; δ) is a subset of the measured set in Eq. 13,
we have a contradiction. Thus Cy − x ≥ 0 holds in a pointwise manner in Rd, indicating{ x ∈ Rd | h(x) 6= Cy − x } must be a null set. With h’s continuity, h must be x 7→ Cy − x. Therefore “ =⇒ ” holds.
Lemma H.1. Cy − x ≥ 0 holds for any x ∈ Rd and its corresponding output y only if CB∗ is a diagonal matrix.
Proof. LetD = CB∗. We rewrite Cy − x ≥ 0 as D [ (A∗x) + + x ] − x ≥ 0. (14)
For further use, we substitute x by −x and the resulting inequality D [ (−A∗x)+ − x ] + x ≥ 0
still holds. Added by Eq. 14 we have D (∣∣A∗k,: · x∣∣)d×1 ≥ 0. (15)
Proof by contradiction: Assume D is not diagonal, then ∃ i 6= j, such that Di,j 6= 0. Consider the following two cases:
a) Di,j > 0. We take the i-th row of Eq. 14 as follows
d∑ k=1 Di,k [( A∗k,: · x )+ + xk ] ≥ xi. (16)
Let x−j = 0 and xj < 0. Then ∑d k=1Di,k [( A∗k,: · x )+ + xk ] = Di,j · xj < 0 =⇒ ⊥.
b) Di,j < 0. We take the i-th row of Eq. 15 as follows
d∑ k=1 Di,k ∣∣A∗k,: · x∣∣ ≥ 0. (17)
Let x = (A∗)−1 v, where v = e(j). Then ∑d
k=1Di,k ∣∣∣A∗k,: · x∣∣∣ = Di,j < 0 =⇒ ⊥. With Lem. 4.2 and Lem. H.1, we prove Thm. G.1 as follows.
Proof. With Lem. 4.2, we only need to prove ∀x ∈ Rd, Cy − x ≥ 0 ⇐⇒ CB∗ = diag(k). “⇐= ”: We have
Cy − x = CB∗ [ (A∗x) + + x ] − x = diag(k) [ (A∗x) + + x ] − x. (18)
For the i-th row of Eq. 18, consider the following two cases:
a) A∗i,: ·x = Ai,ixi, i.e.A∗ is a scale transformation w.r.t. the i-th row. With 11+A∗i,i ≤ ki ≤ 1 we have
[Cy − x]i = ki [( A∗i,ixi )+ + xi ] − xi = ki ( A∗i,ixi )+ + (ki − 1)xi.
For xi ≥ 0, [Cy − x]i = ( kiA ∗ i,i + ki − 1 ) xi ≥ 0.
For xi < 0, [Cy − x]i = (ki − 1)xi ≥ 0. b) A∗i,: · x 6= Ai,ixi. With ki = 1 we have
[Cy − x]i = ki [( A∗i,:x )+ + xi ] − xi = ( A∗i,:x )+ ≥ 0. Hence “⇐= ” holds. “ =⇒ ”: LetD = CB∗. With Lem. H.1,D is diagonal. Consider the following two cases:
a) A∗i,: · x = Ai,ixi. The i-th inequality can be written as
[Cy − x]i = Di,i [( A∗i,ixi )+ + xi ] − xi
= Di,i ( A∗i,ixi )+ + (Di,i − 1)xi ≥ 0. (19)
Proof by contradiction: we need to find x ∈ Rd which contradicts with Eq. 19 in the following three cases:
a) If Di,i ≤ 0, then, xi > 0 =⇒ Di,i ( A∗i,ixi )+ + (Di,i − 1)xi < 0 =⇒ ⊥. b) If Di,i > 1, then, xi < 0 =⇒ ⊥.
c) If 0 < Di,i < 11+A∗i,i , then ∃ a > 0, s.t. Di,i = 1 1+A∗i,i+a . Letting xi > 0, we have
Di,i ( A∗i,ixi )+ + (Di,i − 1)xi =
( A∗i,ixi )+ 1 +A∗i,i + a − ( A∗i,i + a ) xi 1 +A∗i,i + a < 0 =⇒ ⊥.
Hence 11+A∗i,i ≤ Di,i ≤ 1. b) A∗i,: · x 6= Ai,ixi, i.e. ∃ j 6= i, s.t.A∗i,j > 0. The i-th inequality can be written as
[Cy − x]i = Di,i [( A∗i,:x )+ + xi ] − xi
= Di,i ( A∗i,:x )+ + (Di,i − 1)xi ≥ 0. (20)
Proof by contradiction: we need to find x ∈ Rd which contradicts with Eq. 20 in the following three cases:
a) If Di,i ≤ 0, then, xi > 0 =⇒ Di,i ( A∗i,:x )+ + (Di,i − 1)xi < 0 =⇒ ⊥. b) If Di,i > 1, then, xi < 0 ∧ x−i ≤ 0 =⇒ ⊥. c) If 0 < Di,i < 1, then, xi > 0 ∧ xj ≤ −
A∗i,i A∗i,j
xi ∧ xk ≤ 0 =⇒ ⊥, where k 6= i, j. Hence Di,i = 1.
HenceD = CB∗ = diag(k), and thereby “ =⇒ ” holds.

heading: H.2 LAYER 1
text: The proof of Lem. 5.3 is similar to that of Lem. 4.2 because the two lemmas follow the same idea, which is to link objective functional minimization with always-hold inequalities. With Lem. 5.3 , we prove Thm. 5.1 as follows.
Proof. With Lem. 5.3, we only need to prove ∀x ∈ Rd, (A∗x)+ − Ax ≥ 0 ⇐⇒ ∀ i ∈ [d],Ai,: = kiA ∗ i,:, where 0 ≤ ki ≤ 1. This is equivalent to what it is for a single row, i.e. ∀x ∈
Rd, ( a∗>x )+ − a>x ≥ 0 ⇐⇒ a = ka∗, where 0 ≤ k ≤ 1.
“⇐= ”: The case where a∗ = 0 is clear. If a∗ is not a zero vector and a = ka∗, we have a) If a∗>x ≥ 0, ( a∗>x )+ − a>x = a∗>x− ka∗>x = (1− k)a∗>x ≥ 0.
b) If a∗>x < 0, ( a∗>x )+ − a>x = −ka∗>x ≥ 0.
Hence “⇐= ” holds. “ =⇒ ”: If a∗ = 0, a must be a zero vector, otherwise let x = a, then −a>x < 0 =⇒ ⊥. If a∗ is not a zero vector, consider two cases below:
a) If a 6= ka∗ where k ≥ 0. Let x = ‖a ∗‖ ‖a‖ · a− a ∗, then x is not a zero vector, and
a∗>x = ‖a∗‖2 (cos θ − 1) < 0 , a>x = ‖a‖ ‖a∗‖ (1− cos θ) > 0
} =⇒ ( a∗>x )+ − a>x < 0 =⇒ ⊥
where θ denotes the angle between a and a∗. b) If a = ka∗ where k > 1, then let x = a∗ we have(
a∗>x )+ − a>x = (1− k) ‖a∗‖2 < 0 =⇒ ⊥.
Hence a = ka∗ where 0 ≤ k ≤ 1 if a∗ is not zero, and thereby “ =⇒ ” holds.

heading: I QP CONVEXITY
text: The LPs in the main paper are trivially convex. So in this appendix, we only justify the convexity of our QPs: We first prove the convexity of single-sample objectives, then the convexity of the empirical objectives with nonparametric estimation, i.e. ĜNPE1 and Ĝ NPE 2 is obtained by the convexity of convex function summations.
Lemma I.1. Suppose f(u) = 12 ‖Tu− b‖ 2 where u is a real matrix. Then f is convex w.r.t. u.
Lem. I.1 is easily obtained since the Hessian f ′′(T ) = T>T is positive semidefinite. In the following, we demonstrate and justify the convexity of the QPs of both layers by rewriting their single-sample objectives into the formulation of f and summing them without loss of convexity. Theorem I.2. QP: Eq. 4, and QP: Eq. 7 are convex optimization.
Proof. First of all, constraints of both QPs are trivially linear and convex. Thus, we only need to justify the convexity of the two empirical objectives, ĜNPE1 and Ĝ NPE 2 (see Eq. 7 and 4). Consider the single-sample version of ĜNPE1 , namely g NPE 1 (A,φ;x,h) = 1 2 ‖φ+Ax− h‖
2, which, in the formulation of f , can be rewritten with
T =  x> 1 x> 1 . . . . . .
x> 1
 , u =  A>1,: A>2,:
... A>d,: φ  , b = h (21) which guarantees the convexity of gNPE1 w.r.t. A and φ by Lem. I.1. For g NPE 2 (C, ξ;x,y) = 1 2 ‖ξ + x−Cy‖ 2, we have
T =  −y> 1 −y> 1 . . . . . .
−y> 1
 , u =  C>1,: C>2,:
... C>m,: ξ  , b = −x (22) which guarantees the convexity of gNPE2 w.r.t. C and ξ by Lem. I.1. Now we consider the summation. Taking layer 1 as an example, by definition we have
ĜNPE1 (A,Φ) = ∑ i∈[n] gNPE1 (A,φ (i);x(i),h(i)). (23)
For each i ∈ [n], equivalently, we take Φ as a variable instead of φ(i) in gNPE1 , but with only φ(i) ∈ Φ determining the value of gNPE1 . In this sense, g NPE 1 is convex w.r.t. A and Φ for each i ∈ [n]. Thus, the sum ĜNPE1 is convex w.r.t.A and Φ. Similarly, Ĝ NPE 2 is convex w.r.t. C and Ξ.

heading: J STRONG CONSISTENCY
text: In this appendix, we justify the strong consistency of our estimators for the residual unit layer 1/2 learning and the whole network.
According to Thm. G.1 and Thm. 5.1, the solutions to our objective functionals are continuous sets. In addition, there are theoretical intermediate results that are also represented as continuous sets, e.g. possible left-inverse matrices forB∗ in the results for layer 2. Thus, to analyze the consistency of our learning algorithm, we define distances between sets, so that the convergence of sets can be well defined. Further point convergence results, i.e. layer 1/2 estimator strong consistency, are based on the set convergence we define. Definition J.1 (deviation). Let A ⊆ M and B ⊆ M be two non-empty sets from a metric space (M, d). The deviation of set A from the set B, denoted by D(A,B), is
D(A,B) = sup a∈A d(a,B) = sup a∈A inf b∈B d(a, b), (24)
where sup and inf represent supremum and infimum, respectively. Definition J.2 (Hausdorff distance). Let A ⊆M and B ⊆M be two non-empty sets from a metric space (M, d). The Hausdorff distance between A and B, denoted by DH(A,B), is
DH(A,B) = max{D(A,B), D(B,A)}. (25) Remark. In the following, we use the Frobenius norm to define the distance between matrices, i.e. d(X,Y ) = ‖X − Y ‖F.

heading: J.1 LAYER 2
text: In this subsection, we prove the strong consistency of the layer 2 estimator.

heading: J.1.1 OBJECTIVE MINIMIZER SPACE ESTIMATOR
text: For simplicity of notation, we use Ŝn to denote our layer 2 QP/LP solution space by n random samples8 as a random set
Ŝn := {C ∈ Rd×m : Cy(i) − x(i) ≥ 0, ∀ i ∈ [n]} (26)
and S∗ to denote the value space of C in Lem. 4.2 which minimizes layer 2 objective functional
S∗ := {C ∈ Rd×m : Cy − x ≥ 0, ∀x ∈ Rd and its corresponding output y}. (27)
For further use, we name P̂n as the set of n sampled inputs which define Ŝn, i.e. P̂n := {x(i)}i∈[n] where each x(i) is the same random variable in Eq. 26. By the definitions above, we describe the strong consistency of our QP/LP as Lem. J.1.
Lemma J.1 (QP/LP strong consistency, layer 2). DH(Ŝn,S∗) a.s.−−→ 0 as n→∞.
Proof. First we prove that DH(Ŝn,S∗) p−→ 0 as n → ∞. Recall Thm. G.1, Cy − x ≥ 0 only if CB∗ = diag(k). We inherit the notation as definingD = CB∗. Theorem G.1 is based on Lem. H.1, and we prove them by raising points that show contradiction, i.e. violate the inequality that holds in a pointwise fashion:
a) In the proof of Lem. H.1, we use d points: −e(i), for i ∈ [d], to make∑d k=1Di,k [( A∗k,: · x )+ + xk ] < xi, so that Di,j (i 6= j) cannot be positive; and another
d points: (A∗)−1e(i), to show ∑d
k=1Di,k ∣∣∣A∗k,: · x∣∣∣ < 0, so that Di,j (i 6= j) cannot be negative. For each −e(i) we use here, since the violations follow strict inequalities, we know there exists a neighborhood of −e(i), Ni = N(−e(i)), such that ∀z ∈ Ni,∑d
k=1Di,k [( A∗k,: · z )+ + zk ] < zi. We can similarly find such neighborhood of each
(A∗)−1e(i) that the strict inequality holds within the neighborhood respectively. We index them as Nd+1 to N2d. b) In the proof of Thm. G.1, we further construct d points: for each i ∈ [d], we take a point x such that xi > 0 ∧ xj ≤ −
A∗i,i A∗i,j
xi ∧ xk ≤ 0, where k 6= i, j. This counterexample shows [Cy − x]i < 0, and eliminates the possibility of 0 < Di,i < 1 when A∗i,: is not a scale transformation. We can similarly find neighborhood of each point and index them as N2d+1 to N3d. Note that we omit some cases in the proof of Thm. G.1, because the first 2d points are sufficient to use in those cases to show contradiction.
In the sampling procedure, if we sample at least one point in each neighborhood Ni, Thm. G.1 assures the solution we get Ĉn would lie in the true optimal set S∗. The probability that the sampling procedure “omits” any of the neighborhoods is
P ( P̂n ⋂ N1 = ∅ or P̂n ⋂ N2 = ∅ or . . . or P̂n ⋂ N3d = ∅ ) ≤ 3d∑ i=1 P ( P̂n ⋂ Ni = ∅ )
≤ 3d[1− min i∈[3d]
P (Ni)]n (28)
Since the measure on each neighborhood P (Ni) = ∫ x∈Ni p(x) > 0,
P ( P̂n ⋂ N1 = ∅ or P̂n ⋂ N2 = ∅ or . . . or P̂n ⋂ N3d = ∅ ) → 0, as n→∞. (29)
8Here, we take samples as random variables for empirical analysis.
Here we obtain DH(Ŝn,S∗) p−→ 0 as n → ∞. Now we take the infinite sum over the both sides of Eq. 28 ∑ n∈[∞] P ( P̂n ⋂ N1 = ∅ or P̂n ⋂ N2 = ∅ or . . . or P̂n ⋂ N3d = ∅ )
≤ 3d ∑
n∈[∞]
[ 1− min
i∈[3d] P (Ni)
]n = 3d [ 1 mini∈[3d] P (Ni) − 1 ] < +∞.
By the Borel-Cantelli lemma (Borel, 1909), DH(Ŝn,S∗) a.s.−−→ 0 as n→∞.

heading: J.1.2 SCALE FACTOR ESTIMATOR
text: To avoid ambiguity, we use nsf to denote the number of samples used in Alg. 4. The samples pairs are {(x(i),y(i))}nsfi=1. Without loss of generality, the following discussion focuses on some fixed index j ∈ [d]. In Alg. 4, we plug in our estimator Ĉn and use LR to estimate kj given that x(i)j < 0
knsf(Ĉn) = arg min k
1
2nsf ∑ i∈[nsf] ∥∥∥[Ĉny(i)]j − kx(i)j ∥∥∥2 = ∑ i∈[nsf] x (i) j [Ĉny (i)]j∑ i∈[nsf] ( x(i)j )2 . (30)
We first give the strong consistency of layer 2 scale factor estimator for as nsf →∞, as described in Lem. J.2. Lemma J.2 (scale factor estimator strong consistency, layer 2). SupposeA∗ is a scale transformation w.r.t. the j-th component, and knsf(Ĉn) is the nsf-sample estimator of kj via LR: Eq. 30 given that x (i) j < 0. Define sets
Unsf,n := {knsf(C) : C ∈ Ŝn}, and U∗ :=
[ 1
1 +A∗j,j , 1
] . (31)
Then lim nsf→∞ lim n→∞
DH(Unsf,n,U∗) a.s. == 0.
Proof. Following our notation, Thm. G.1 and Thm. G.2 ensure that ifA∗ is a scale transformation w.r.t. the j-th component, for any C belonging to the true optimal set S∗, knsf(C) ∈
[ 1 1+A∗j,j , 1 ] .
And the “iff” statement strengthens that U∗ = {knsf(C) : C ∈ S∗} for any nsf ∈ Z+. Note that since S∗ ⊂ Ŝn, we have U∗ ⊂ Unsf,n. We only need to prove D(Unsf,n,U∗) a.s.−−→ 0 as n→∞.
∀ Ĉn ∈ Ŝn and ∀C ∈ S∗,
∣∣∣knsf(Ĉn)− knsf(C)∣∣∣ =
1∑ i∈[nsf] ( x(i)j )2 ∣∣∣∣∣∣ ∑ i∈[nsf] x(i)j [ (e(j))> ( Ĉn −C ) B∗ [( A∗x(i) )+ + x(i) ]]∣∣∣∣∣∣ ≤ 1∑
i∈[nsf]
( x(i)j )2  ∑ i∈[nsf] ∣∣∣x(i)j ∣∣∣ ∥∥∥(e(j))> (Ĉn −C)B∗∥∥∥ 2 (∥∥∥A∗x(i)∥∥∥ 2 + ∥∥∥x(i)∥∥∥ 2 ) ≤ 1∑
i∈[nsf]
( x(i)j )2  ∑ i∈[nsf] ∣∣∣x(i)j ∣∣∣ ∥∥∥Ĉn −C∥∥∥ F ‖B∗‖F (‖A ∗‖F + 1) ∥∥∥x(i)∥∥∥ 2 
=
∑ i∈[nsf] [∣∣∣x(i)j ∣∣∣ ‖B∗‖F (‖A∗‖F + 1) ∥∥x(i)∥∥2]∑ i∈[nsf] ( x(i)j )2 ∥∥∥Ĉn −C∥∥∥ F .
Then
sup Ĉn∈Ŝn inf C∈S∗ ∣∣∣knsf(Ĉn)− knsf(C)∣∣∣ (32) ≤ ∑ i∈[nsf]
[∣∣∣x(i)j ∣∣∣ ‖B∗‖F (‖A∗‖F + 1) ∥∥x(i)∥∥2]∑ i∈[nsf] ( x(i)j )2 sup Ĉn∈Ŝn inf C∈S∗ ∥∥∥Ĉn −C∥∥∥ F
(33)
which implies
D(Unsf,n,U) ≤
∑ i∈[nsf] [∣∣∣x(i)j ∣∣∣ ‖B∗‖F (‖A∗‖F + 1) ∥∥x(i)∥∥2]∑ i∈[nsf] ( x(i)j )2 D(Ŝn,S∗). (34)
Take the nsf →∞ limit over both sides of Eq. 34. With the strong law of large numbers9 we have
lim nsf→∞
D(Unsf,n,U) ≤ E [∣∣∣x(i)j ∣∣∣ ‖B∗‖F (‖A∗‖F + 1) ∥∥x(i)∥∥2]
E [ x(i)j ]2 D(Ŝn,S∗), w.p. 1. (35)

heading: Since D(Ŝn,S∗)
text: a.s.−−→ 0 as n→∞, we have D(Unsf,n,U∗) a.s.−−→ 0 as n→∞ then nsf →∞.

heading: J.1.3 LAYER 2 WEIGHTS ESTIMATOR
text: In Alg. 5, we solve B via LR. Let ẑ(i) = diag−1(k̂) · Ĉny(i) ∈ Rd, where k̂ is obtained through Alg. 4 with input Ĉn. Assume we are using sample size of nw to do the LR. The optimization problem is
min B ∑ i∈[nw] ∥∥∥y(i) −Bẑ(i)∥∥∥2 (36) Now we present the strong consistency of layer 2 estimator, as described in Thm. J.3.
Theorem J.3 (strong consistency, layer 2). Suppose B̂nsf is the solution to Eq. 36. Then B̂nw a.s.−−→ B∗ as n, nsf, nw →∞.
Proof. Let β denote vecB (flattening B into a vector), then Bẑ(i) = ( ẑ(i) )> ⊗ Im. Here the operation ⊗ denotes the Kronecker product. Then we can define an equivalent optimization problem
min β
1
2nw ∑ i∈[nw] ∥∥∥∥y(i) − [(ẑ(i))> ⊗ Im]β∥∥∥∥2 . (37) Take the derivatives of β, we obtain
− 2 ∑
i∈[nw]
[[( ẑ(i) )> ⊗ Im ]>( y(i) − [( ẑ(i) )> ⊗ Im ] β )] = 0 (38)
9Here we assume that the Kolmogorov’s strong law assumption on moments (Sen & Singer, 1994) is met as is commonly done in empirical analysis.
Then the optimal solution β̂nw of this optimization can be written in closed form:
β̂nw =  ∑ i∈[nw] [( ẑ(i) )> ⊗ Im ]> [( ẑ(i) )> ⊗ Im ]−1 ∑ i∈[nw] [( ẑ(i) )> ⊗ Im ]> y(i)  =
 ∑ i∈[nw] [ ẑ(i) ⊗ Im ] [( ẑ(i) )> ⊗ Im ]−1 ∑ i∈[nw] [ ẑ(i) ⊗ Im ] y(i)  =
 ∑ i∈[nw] [ ẑ(i) ( ẑ(i) )>] ⊗ Im −1 ∑ i∈[nw] ẑ(i) ⊗ Im y(i) 
=  ∑ i∈[nw] ẑ(i) ( ẑ(i) )>⊗ Im −1 ∑ i∈[nw] ẑ(i) ⊗ Im y(i) 
We inherit the notation from the last two subsections. By Lem. J.1, d(Ĉn,S∗) a.s.−−→ 0 as n → ∞. Thus ∀ ε > 0, ∃N such that ∀n ≥ N , d(Ĉn,S∗) ≤ ε w.p. 1, i.e. ∃Cn ∈ S∗ s.t. d(Ĉn,Cn) ≤ ε w.p. 1. Then by Lem. J.2, ∃K > 0, for ε that is small enough, ∃Nsf such that ∀nsf ≥ Nsf, we have
∣∣∣knsf(Ĉn)− knsf(Cn)∣∣∣ ≤ Kε w.p.1. For simplicity, we omit the under-script nsf of knsf in the following discussion. In fact,
∣∣∣ẑ(i)j − z(i)j ∣∣∣ = ∣∣∣∣∣ 1k(Ĉn) ( e(i) )> Ĉny (i) − 1 k(Cn) ( e(i) )> Cny (i) ∣∣∣∣∣ =
∣∣∣∣∣ 1k(Ĉn)k(Cn) [ k(Cn) ( e(i) )> Ĉn − k(Ĉn) ( e(i) )> Cn ] y(i) ∣∣∣∣∣ ≤
∣∣∣∣∣ 1k(Ĉn)k(Cn) ∣∣∣∣∣ ∣∣∣∣k(Cn)(e(i))> Ĉny(i) − k(Cn)(e(i))>Cny(i)∣∣∣∣
+ ∣∣∣∣∣ 1k(Ĉn)k(Cn) ∣∣∣∣∣ ∣∣∣∣k(Cn)(e(i))>Cny(i) − k(Ĉn)(e(i))>Cny(i)∣∣∣∣
≤ ( 1 +A∗j,j )2 1−Kε ( 1 +A∗j,j ) (∥∥∥y(i)∥∥∥∥∥∥Ĉn −Cn∥∥∥ F + ∣∣∣k(Ĉn)− k(Cn)∣∣∣ ‖CnB∗‖F (‖A∗‖F + 1) ∥∥∥x(i)∥∥∥)
≤ ( 1 +A∗j,j )2 1−Kε ( 1 +A∗j,j ) (∥∥∥y(i)∥∥∥+ dK (‖A∗‖F + 1) ∥∥∥x(i)∥∥∥) ε (39) Then ∥∥∥∥ẑ(i) (ẑ(i))> − z(i) (z(i))>∥∥∥∥
F
= ∥∥∥∥ẑ(i) (ẑ(i))> − z(i) (ẑ(i))> + z(i) (ẑ(i))> − z(i) (z(i))>∥∥∥∥ F
≤ ∥∥∥∥[ẑ(i) − z(i)] (ẑ(i))>∥∥∥∥
F
+ ∥∥∥∥z(i) [(ẑ(i))> − (z(i))>]∥∥∥∥ F
≤ ∥∥∥∥[ẑ(i) − z(i)] [ẑ(i) − z(i)]>∥∥∥∥
F
+ 2 ∥∥∥∥z(i) [(ẑ(i))> − (z(i))>]∥∥∥∥ F
≤ ∑ j∈[d] ( ẑ (i) j − z (i) j )2 + 2 ∥∥∥z(i)∥∥∥∥∥∥ẑ(i) − z(i)∥∥∥
It follows that ∥∥∥ẑ(i) (ẑ(i))> − z(i) (z(i))>∥∥∥
F is also bounded by O(ε). With similar techniques we
can prove ∥∥∥∥∥∥ ∑
i∈[nw]
ẑ(i) ( ẑ(i) )> − ∑
i∈[nw]
z(i) ( z(i) )>∥∥∥∥∥∥
F
≤ O(ε) (40)
and ∥∥∥ẑ(i) ⊗ Im − z(i) ⊗ Im∥∥∥ F ≤ O(ε) (41)
Denote [∑ i∈[nw] z (i) ( z(i) )>] as P and ∑ i∈[nw] z (i) as Q. Substitute z(i) with ẑ(i) in the above
expression we have P̂ and Q̂. Hence,∥∥∥β̂nw − β∗∥∥∥ F = ∥∥∥∥[P̂ ⊗ Im]−1 ([Q̂⊗ Im]y(i))− [P ⊗ Im]−1 ([Q⊗ Im] y(i))∥∥∥∥ F
≤ ∥∥∥∥[P̂ ⊗ Im]−1 ([Q̂⊗ Im]y(i))− [P̂ ⊗ Im]−1 ([Q⊗ Im] y(i))∥∥∥∥
F
+ ∥∥∥∥[P̂ ⊗ Im]−1 ([Q⊗ Im] y(i))− [P ⊗ Im]−1 ([Q⊗ Im] y(i))∥∥∥∥ F
= ∥∥∥∥[P̂ ⊗ Im]−1∥∥∥∥ F ∥∥∥[(Q̂−Q)⊗ Im]y(i)∥∥∥ F + ∥∥∥∥[P̂ ⊗ Im]−1 − [P ⊗ Im]−1∥∥∥∥ F∥∥∥[Q⊗ Im] y(i)∥∥∥
F
In the first part, by triangle inequality,∥∥∥∥[P̂ ⊗ Im]−1∥∥∥∥ F ≤ ∥∥∥∥[P̂ ⊗ Im]−1 − [P ⊗ Im]−1∥∥∥∥ F + ∥∥∥[P ⊗ Im]−1∥∥∥ F (42)
So we only need to prove ∥∥∥∥[P̂ ⊗ Im]−1 − [P ⊗ Im]−1∥∥∥∥
F
≤ O(ε) to claim ∥∥∥β̂nw − β∗∥∥∥
F ≤ O(ε).
Denote P̂ − P = ∆P . From Eq. 39, we know every entry of ∆P ⊗ Im can be bounded by O(ε). By simple calculation we have
(P + ∆P )−1 = P−1 − P−1∆PP−1 +O(ε2). (43) Then we have
P−1 − P̂−1 = P−1∆PP−1 +O(ε2) = O(ε). (44)

heading: J.2 LAYER 1
text: In this subsection, we justify the strong consistency of layer 1 objective functional minimizer estimator in detail, i.e. the layer 1 QP/LP solution space. We will omit the detailed proof of Alg. 2 line 4 to 7 strong consistency since it is similar to the proof of Lem. J.2. Additionally, we also omit the strong consistency of the x 7→ (A∗x)+ function estimator because it can be directly obtained by Lem. J.1 and J.2 and the continuous mapping theorem.
We use a new optimization problem equivalent to the optimization of G1. Before that, we first define the equivalence between two optimization problems as follows. Definition J.3. Let opt1 and opt2 be two optimization problems, with f1, f2 as the respective objective functions. Then opt1 and opt2 are said to be equivalent if given a feasible solution to opt1, namely x1, a feasible solution to opt2 is uniquely corresponded, namely x2, such that f1(x1) = f2(x2), and vice versa.
The new optimization problem and its equivalence to the optimization of G1 is described in Lem. J.4. Lemma J.4. The optimization of G1 (Eq. 6) is equivalent to
min A
f(A) = 1 2 Ex [∥∥∥(Ax− h)+∥∥∥2] . (45)
Proof. To see this, suppose A1 is one optimal solution to Eq. 45, then we can construct r1(x) = (A1x− h)+ so that G1 (A1, r1) = f(A1) and the optimality implies
min A, r G1 (A, r) ≤ min A f(A).
On the other hand, suppose (A2, r2) is an optimum of G1. Let r3(x) = (h−A2x)+, then ∀x ∈ Rd and h be the corresponding hidden output, if [h − A2x]j ≥ 0, then [r3(x) + A2x − h]j = 0, otherwise [r3(x) +A2x− h]2j = [A2x− h]2 ≤ [h2 (x) +A2x− h]2j since h2 is nonnegative. So that we know that
min A, r G1(A, r) = G1(A2, r2) = G1(A2, r3) = f(A2)
From the optimality, we further have
min A, r G1(A, r) ≥ min A f(A)
From the simple calculation above, we can see that one optimal solution to Eq. 45 has a one-to-one correspondence to an optimal solution to G1.
Similarly, the empirical version of the two problems are equivalent, which indicates their consistency in the empirical estimation being equivalent. In the following, we justify the strong consistency of empirical Eq. 45 instead of G1
min A
f̂n (A) = 1
2n ∑ i∈[n] ∥∥∥∥(Ax(i) − h(i))+∥∥∥∥2 . (46) Denote T∗ := {diag(k) · A∗ | 0 ≤ kj ≤ 1, j ∈ [d]} as the true optimal solution set, and T̂n as the optimal solution set corresponding to the n-sample problem. In the following, we justify four conditions in a row that hold for f to derive the strong consistency of its optimal solution estimator.
Lemma J.5. Let T̂n′ be the layer 1 QP/LP solution space by n′ samples. Then there exists a compact set C determined byA∗, namely C(A∗), s.t. T̂n′ ⊂ C(A∗) w.p. 1 as n→∞.
Proof. ∀l ∈ [d], let a∗ be the l-th row ofA∗, and âl be the l-th row of Ân′ . We’d like first to prove that the set T̂ln′ = {al : al is the l-th row ofA, whereA ∈ T̂n′} (47) is compact w.p.1.
Suppose n′ > d, and among the n′ samples, we classify them into two folds. To avoid ambiguity, let u(i) be the points such that (a∗)>u(i) > 0, i ∈ [q]; and v(j) be the points such that (a∗)>v(j) < 0, j ∈ [n− q]. From the analysis of Thm. 5.1, we have (a∗)>u(i) ≥ â>u(i), ∀i ∈ [q] and â>v(j) ≤ 0, ∀j ∈ [n− q]. It follows that we can rewrite T̂ln′ as a polyhedron
T̂ln′ = {a ∈ Rd : a>u(i) ≤ (a∗)>u(i),a>v(j) ≤ 0} (48)
We are going to show the polyhedron T̂ln′ is bounded by contradiction. If it is not bounded, then ∃d ∈ Rd, d 6= 0 and ã ∈ T̂ln′ , such that ∀λ > 0, ã+ λd ∈ T̂ln′ . Then
(ã+ λd) > u(i) = ã>u(i) +λd>u(i) ≤ (a∗)> u(i) ⇐⇒ λd>u(i) ≤ (a∗)> u(i)− ã>u(i) (49)
similarly, (ã+ λd) > v(i) = ã>v(j) + λd>v(j) ≤ 0 (50)
From the definition, λ can be arbitrarily big, then d>u(i) ≤ 0, ∀i ∈ [q], and d>v(j) ≤ 0, ∀j ∈ [n−q].
Since we know span{u(i)} = Rd w.p. 1, then there ∃ some i∗ such that d>u(i∗) < 0, w.p. 1. (Otherwise if d>u(i) = 0 for all i ∈ [q], then either span{u(i)} 6= Rd or d = 0.) Under our assumption that, T̂ln′ is not bounded, we know the following system (w.r.t x) has a feasible solution w.p. 1 [
−U −V
] x ≥ 0, x>u(i ∗) < 0 (I)
where every row of U and V is ( u(i) )> and ( v(j) )> respectively. By Farkas’ lemma (Farkas, 1902), the system [−U>,−V >] · x = u(i
∗), x ≥ 0 (II) is not feasible (w.p. 1). We claim that u(i
∗) lies in the conic hull of −v(j)’s w.p. 1. So that the second system actually has a feasible solution and thus it raises the contradiction.
Denote the conic hull as
H = t ∈ Rd : t = ∑ j∈[n−q] λj ( −v(j) ) , λj ≥ 0 for ∀j ∈ [n− q]  Now suppose u(i
∗) /∈ H, by the supporting hyperplane theorem (Luenberger, 1997), ∃b ∈ Rd, b 6= 0, such that b>u(i
∗) ≤ b>t for ∀t ∈ H. Then by definition, −v(j) ∈ {t : b>t ≤ b>u(i
∗)}, for ∀j ∈ [n− q] Denote the hyperplane J = {t : b>t ≤ b>u(i∗)}, then
P ( u(i ∗) /∈ H ) ≤ P ( −v(j) ∈ J, for ∀j ∈ [n− q] ) = P ( −v(j) ∈ J )n−q Since we have in this case a geometric sequence, we know its infinite sum is bounded. By BorelCantelli lemma (Borel, 1909), we conclude that u(i
∗) ∈ H w.p. 1. Then system II is feasible w.p. 1. So that T̂ln′ is compact w.p. 1.
Now we prove that there exists a compact set C(A∗), s.t. T̂n′ ⊂ C(A∗) w.p. 1 as n→∞. Similarly, we focus on the analysis of one row. As discussed above, T̂ln′ is compact w.p. 1. Let Ŵ1n′ be the set of all u(i) sampled in estimating T̂n′ , and Ŵ2n′ be the set of all v(j) sampled. Now for another set T̂n′′ , similarly define sample point sets Ŵ1n′′ and Ŵ2n′′ . We claim that ∀u(i) ∈ Ŵ1n′ , u(i) lies in the conic hull of Ŵ1n′′ . Actually, this part of the proof is very similar to the way we prove u(i
∗) ∈ H w.p. 1, so we will omit it here.
∀a ∈ T̂n′′ , let x(1) and x(2) be two different points in Ŵ1n′′ . Then a>(λ1x(1) + λ2x(2)) ≤ (a∗)>(λ1x (1) + λ2x (2)) for ∀λ1 ≥ 0 and λ2 ≥ 0. This simple calculation reveals a>u(i) ≤ (a∗)>u(i) for ∀u(i) ∈ Ŵ1n′ (from the claim we made). This implies that T̂n′′ ⊂ C(A∗) w.p. 1 as n→∞, too.
Lemma J.6. The minimizer space of f(A), i.e. T∗ = {diag(k) · A∗ | 0 ≤ kj ≤ 1, j ∈ [d]}, is contained in C(A∗). Lemma J.7. f(A) is finite valued and continuous on C(A∗).
Lem. J.6 and J.7 are easily obtained by the formulation of f (see Eq. 45) and Lem. J.5.
Lemma J.8 (uniform a.s. convergence). f̂n(A) a.s.−−→ f(A) as n→∞, uniformly inA ∈ C(A∗). Proof. Name single-sample objective g(x,A) = 12 ∥∥∥(Ax− h)+∥∥∥2. The uniform a.s. convergence is guaranteed by the uniform law of large numbers (Jennrich, 1969):
a) By Lem. J.5, C(A∗) is a compact set. b) g is continuous w.r.t.A by its formulation and measurable over x at eachA ∈ C(A∗). c) In fact,
g(x,A) = 1
2 ∥∥∥(Ax− h)+∥∥∥2 ≤ ‖Ax‖2 + ‖A∗x‖2
≤ (‖A‖F + ‖A ∗‖F) ‖x‖ 2 .
SinceA ∈ C(A∗) is in a compact set,
g(x,A) ≤ [ sup
A∈C(A∗) ‖A‖F + ‖A ∗‖F
] ‖x‖2 . (51)
Thus the dominating function exists10.
By Lem. J.5, J.6, J.7 and J.8, all of the conditions are satisfied in (Shapiro et al., 2014, Thm. 5.3). Thus, we have the strong consistency of layer 1 objective optima estimator as described in Lem. J.9.
Lemma J.9 (QP/LP strong consistency, layer 1). DH(T̂n′ ,T∗) a.s.−−→ 0 as n→∞.
Similar to Lem. J.2, we have the strong consistency of the layer 1 scale factor estimator as described in Lem. J.10. Lemma J.10 (scale factor estimator strong consistency, layer 1). Let kn′sf(Ĉn) be the n ′ sf-sample estimator of kj via LR: Alg. 2 line 5. given that h (i) j > 0. Define sets
Vn′sf,n := {kn′sf(A) : A ∈ T̂n′}, and V ∗ := [0, 1] . (52)
Then lim n′sf→∞ lim n′→∞ DH(Vn′sf,n′ ,V ∗) a.s. == 0.
Remark. In case kj = 0, suppose the algorithm finds a solution over a continuous distribution with [0, 1] as support and the probability that it finds a solution with scale factor 0 is 0.
With Thm. 5.2 and the continuous mapping theorem, the strong consistency of layer 1 estimation is guaranteed.
Theorem J.11 (strong consistency, layer 1). Suppose Ân′ is scaled by k̂n′sf . Then Ân′ a.s.−−→ A∗ as n′, n′sf →∞.
By Thm. J.3 and J.11, Thm. 6.3 is guaranteed by the continuous mapping theorem.
10Here, we assume E[‖x‖2] < +∞ as is commonly done in empirical analysis.

Review summaries:
This paper gives an algorithm for learning two-layer neural networks with ReLU activation (realizable case). The network also consists of skip connections to avoid the problems of identifiability. As a warm up, the authors first give a simple algorithms based on linear regression, whenever enough samples could be observed with all-positive and all-negative entries.   They start with an algorithm for recovering the layer two parameters. For this step, they use the $\ell_2$ loss and show that its minimizer is unique and matches with the ground truth parameters under certain condition. Subsequently, the authors replace the expectation in the objective with an average over the samples (based on the empirical risk minimization principle), which they could write as a convex quadratic program without making any parametric assumption. Finally, they show that the quadratic program can be equivalently written as a linear program, which lets them recover the parameters of the second layer (in fact the inverse matrix of this weight matrix).  Once layer 2 is learnt, the original problem reduces to that of learning the single layer non-linearity. They write the objective as a sum of weight estimator and a function estimator. They show that the minimizer of this loss is achieved when the estimated weight matrix is a diagonal matrix times the ground truth and the function measures the difference between the true nonlinear function and the scaled wright matrix. Again, they learn the function non-parametrically using a quadratic program, which they equivalently write as a linear program to get the scaled weight matrix. The scaling factors in the diagonal matrix can then be learnt using a linear regression.  Finally, the authors conduct extensive experimentation to validate that their algorithms performs superior to the vanilla SGD-based algorithm. They also show that their algorithm, especially the QP algorithm is robust to noise.

This paper considers the problem of learning the parameters of a two-layer ReLU network with a residual unit: given samples of the form $(x_i, y_i)$, with $y_i = B^*[(A^* x_i)^+ + x_i]$, the authors provide an algorithm that provably learns $A^*$ and $B^*$, as the sample size grows. It is crucially assumed that the entries of $A^*$ are non-negative and, by exploiting this fact, the learning problem is formulated via quadratic programming (QP): the authors show that there exist quadratic functionals whose minimisers can be used to deduce the ground-truth parameters. Furthermore, the QPs are re-written as linear programs (LPs), which are simpler to optimise and have the same solution space. Numerical results demonstrate the superiority of the proposed approach with respect to gradient descent methods.