Robust and Resource Efficient Identification of Two Hidden Layer Neural Networks

1 Introduction

Deep learning is perhaps one of the most sensational scientific and technological developments in the industry of the last years. Despite the spectacular success of deep neural networks (NN) outperforming other pattern recognition methods, achieving even superhuman skills in some domains [12, 36, 58], and confirmations of empirical successes in other areas such as speech recognition [25], optical charachter recognition [8], games solution [44, 56]

, the mathematical understanding of the technology of machine learning is in its infancy. This is not only unsatisfactory from a scientific, especially mathematical point of view, but it also means that deep learning currently has the character of a black-box method and its success can not be ensured yet by a full theoretical explanation. This leads to lack of acceptance in many areas, where interpretability is a crucial issue (like security, cf.

[10]) or for those applications where one wants to extract new insights from data [61].

Several general mathematical results on neural networks have been available since the 90’s [2, 17, 38, 39, 46, 47, 48], but deep neural networks have special features and in particular superior properties in applications that still can not be fully explained from the known results. In recent years, new interesting mathematical insights have been derived for undestanding approximation properties (expressivity) [27, 54] and stability properties [9, 68] of deep neural networks. Several other crucial and challenging questions remain open.

A fundamental one is about the number of required training data to obtain a good neural network, i.e.

, achieving small generalization errors for future data. Classical statistical learning theory splits this error into bias and variance and gives general estimations by means of the so-called VC-dimension or Rademacher complexity of the used class of neural networks

[55]. However, the currently available estimates of these parameters [26] provide very pessimistic barriers in comparison to empirical success. In fact, the tradeoff between bias and variance is function of the complexity of a network, which should be estimated by the number of sampling points to identify it uniquely. Thus, on the one hand, it is of interest to know which neural networks can be uniquely determined in a stable way by finitely many training points. On the other hand, the unique identifiability is clearly a form of interpretability.

The motivating problem of this paper is the robust and resource efficient identification of feed forward neural networks. Unfortunately, it is known that identifying a very simple (but general enough) neural network is indeed NP-hard

[7, 33]. Even without invoking fully connected neural networks, recent work [20, 41]

showed that even the training of one single neuron (ridge function or single index model) can show any possible degree of intractability, depending on the distribution of the input. Recent results

[3, 34, 42, 57, 52], on the other hand, are more encouraging, and show that minimizing a square loss of a (deep) neural network does not have in general or asymptotically (for large number of neurons) poor local minima, although it may retain the presence of critical saddle points.
In this paper we present conditions for a fully nonlinear two-layer neural network to be provably identifiable with a number of samples, which is polynomially depending on the dimension of the network. Moreover, we prove that our procedure is robust to perturbations. Our result is clearly of theoretical nature, but also fully constructive and easily implementable. To our knowledge, this work is the first, which allows provable de-parametrization of the problem of deep network identification, beyond the simpler case of shallow (one hidden) layer neural networks already considered in very recent literature [3, 32, 21, 34, 42, 43, 57, 52]

. For the implementation we do not require black-box high dimensional optimization methods and no concerns about complex energy loss landscapes need to be addressed, but only classical and relatively simple calculus and linear algebra tools are used (mostly function differentiation and singular value decompositions). The results of this paper build upon the work

[20, 21], where the approximation from a finite number of sampling points have been already derived for the single neuron and one-layer neural networks. The generalization of the approach of the present paper to networks with more than two hidden layers is suprisingly simpler than one may expect, and it is in the course of finalization [22], see Section 5 (v) below for some details.

1.1 Notation

Let us collect here some notation used in this paper. Given any integer $m \in N$ , we use the symbol $[m] := {1, 2, \dots, m}$ for indicating the index set of the first $m$ integers. We denote $B_{1}^{d}$ the Euclidean unit ball in $R^{d}$ , $S^{d - 1}$ the Euclidean sphere, and $μ_{S^{d - 1}}$

is its uniform probability measure. We denote

ℓ_{q}^{d}

the

d

-dimensional Euclidean space endowed with the norm

∥ x ∥_{ℓ_{q}^{d}} = {(\sum_{j = 1}^{d} | x_{j} |^{q})}_{j}^{1 / q}

. For

q = 2

we often write indifferently

∥ x ∥ = ∥ x ∥_{2} = ∥ x ∥_{ℓ_{2}^{d}}

. For a matrix

M

we denote

σ_{k} (M)

its

k^{t h}

singular value. We denote

S

the sphere of symmetric matrices of unit Frobenius norm

∥ \cdot ∥_{F}

. The spectral norm of a matrix is denoted

∥ \cdot ∥

. Given a closed convex set

C

we denote

P_{C}

the orthogonal projection operator onto

C

(sometimes we use such operators to project onto subspaces of

R^{d}

or subspaces of symmetric matrices or onto balls of such spaces). For vectors

x_{1}, \dots, x_{k} \in R^{d}

we denote the tensor product

x_{1} \otimes \dots \otimes x_{k}

as the tensor of entries . For the case of

k = 2

the tensor product

x \otimes y

of two vectors

x, y \in R^{d}

equals the matrix

x y^{T} = (x_{i} y_{j})_{i j}

. For any matrix

M \in R^{m \times n}

vec (M) := (m_{11}, m_{21}, \dots, m_{m 1}, m_{12}, m_{22}, \dots, m_{m n})^{T} \in R^{m n},

(1)

is its vectorization, which is the vector created by the stacked columns of $M$ .

1.2 From one artificial neuron to shallow, and deeper networks

1.2.1 Meet the neuron

The simplest artificial neural network $f : Ω \subset R^{d} \to R$ is a network consisting of exactly one artificial neuron, which is modeled by a ridge-function (or single-index model) $f$ as

f (x) = ϕ (a^{T} x + θ) = g (a^{T} x),

(2)

where $g : R \to R$ is the shifted activation function $ϕ (\cdot + θ)$ and the vector $a \in R^{d}$ expresses the weight of the neuron. Since the beginning of the 90’s [31, 30], there is a vast mathematical statistics literature about single-index models, which addresses the problem of approximating $a$ and possibly also $g$ from a finite number of samples of $f$ to yield an expected least-squares approximation of $f$ on a bounded domain $Ω \subset R^{d}$

. Now assume for the moment that we can evaluate the network

f

at any point in its domain; we refer to this setting as active sampling. As we aim at uniform approximations, we adhere here to the language of recent results about the sampling complexity of ridge functions from the approximation theory literature, e.g., [13, 20, 41]. In those papers, the identification of the neuron is performed by using approximate differentiation. Let us clarify how this method works as it will be of inspiration for the further developments below. For any

ϵ > 0

, points

x_{i}

i = 1, \dots m_{X}

, and differentiation directions

φ_{j}

j = 1, \dots m_{Φ}

we have

\frac{f (x_{i} + ϵ φ_{j}) - f (x_{i})}{ϵ} \approx \frac{\partial f (x_{i})}{\partial φ_{j}} = g^{'} (a^{T} x_{i}) a^{T} φ_{j} .

(3)

Hence, differentiation exposes the weight of a neuron and allows to test it against test vectors $φ_{j}$ . The approximate relationship (3) forms for every fixed index $i$ a linear system of dimensions $m_{Φ} \times d$ , whose unknown is $x_{i}^{*} = g^{'} (a^{T} x_{i}) a$ . Solving approximately and independently the systems for $i = 1, \dots m_{X}$ yields multiple approximations $^a = x_{i}^{*} / ∥ x_{i}^{*} ∥_{2} \approx a$ of the weight, the most stable of them with respect to the approximation error in (3) is the one for which $∥ x_{i}^{*} ∥_{2}$ is maximal. Once $^a \approx a$ is learned then one can easily construct a function $^f (x) =^g ({^a}^{T} x)$ by approximating $^g (t) \approx f (^a t)$ on further sampling points. Under assumptions of smoothness of the activation function $g \in C^{s} ([0, 1])$ , for $s > 1$ , $g^{'} (0) \neq 0$ , and compressibility of the weight, i.e., $∥ a ∥_{ℓ_{q}^{d}}$ is small for $0 < q \leq 1$ , then by using $L$ sampling points of the function $f$ and the approach sketched above, one can construct a function $^f (x) =^g ({^a}^{T} x)$ such that

∥ f -^f ∥_{C (Ω)} \leq C ∥ a ∥_{ℓ_{q}^{d}} ⎧ ⎨ ⎩ L^{- s} + ∥ g ∥_{C^{s} ([0, 1])} {(\frac{1 + log (d / L)}{L})}^{1 / q - 1 / 2} ⎫ ⎬ ⎭ .

(4)

In particular, the result constructs the approximation of the neuron with an error, which has polynomial rate with respect to the number of samples, depending on the smoothness of the activation function and the compressibility of the weight vector $a$ . The dependence on the input dimension is only logarithmical. To take advantage of the compressibility of the weight, compressive sensing [23] is a key tool to solve the linear systems (3). In [13] such an approximation result was obtained by active and deterministic choice of the input points $x_{i}$ . In order to relax a bit the usage of active sampling, in the paper [20] a random sampling of the points $x_{i}$ has been proposed and the resulting error estimate would hold with high probability. The assumption $g^{'} (0) \neq 0$ is somehow crucial, since it was pointed out in [20, 41] that any level of tractability (polynomial complexity) and intractability (super-polynomial complexity) of the problem may be exhibited otherwise.

1.2.2 Shallow networks: the one-layer case

Combining several neurons leads to richer function classes [38, 39, 46, 47, 48]. A neural network with one hidden layer and one output is simply a weighted sum of neurons whose activation function only differs by a shift, i.e.,

f (x) = m \sum i = 1 b_{i} ϕ (a_{i}^{T} x + θ_{i}) = m \sum i = 1 g_{i} (a_{i}^{T} x),

(5)

where $a_{i} \in R^{m}$ and $b_{i}, θ_{i} \in R$ for all $i = 1, \dots, m$ . Sometimes, it may be convenient below the more compact writing $f (x) = 1^{T} g (A^{T} x)$ where $g = (g_{1}, \dots, g_{m})$ and $A = [a_{1} | \dots | a_{m}] \in R^{d \times m}$ ¹¹1Below, with slight abuse of notation, we may use the symbol $A$ also for the span of the weights ${a_{i} : i = 1, \dots, m}$ .. Differently from the case of the single neuron, the use of first order differentiation

\nabla f (x) = m \sum i = 1 g_{i}^{'} (a_{i}^{T} x) a_{i} \in A = span {a_{1}, \dots, a_{m}},

(6)

may furnish information about $A = span {a_{1}, \dots, a_{m}}$ (active subspace identification [14, 15], see also [20, Lemma 2.1]), but it does not allow yet to extract information about the single weights $a_{i}$ . For that higher order information is needed. Recent work shows that the identification of a network (5) can be related to tensor decompositions [1, 32, 21, 43]. As pointed out in Section 1.2.1 differentiation exposes the weights. In fact, one way to relate the network to tensors and tensor decompositions is given by higher order differentiation. In this case the tensor takes the form

D^{k} f (x) = m \sum i = 1 g_{i}^{(k)} (x) a_{i} \otimes \dots \otimes a_{i}      k - times,

which requires that the $g_{i}$ ’s are sufficiently smooth. In a setting where the samples are actively chosen, it is generally possible to approximate these derivatives by finite differences. However, even for passive sampling there are ways to construct similar tensors [32, 21], which rely on Stein’s lemma [59] or differentiation by parts or weak differentiation. Let us explain how passive sampling in this setting may be used for obtaining tensor representations of the network. If the probability measure of the sampling points $x_{i}$ ’s is $μ_{X}$ with known (or approximately known [18]) density $p (x)$ with respect to the Lebesgue measure, i.e., $d μ_{X} (x) = p (x) d x$ , then we can approximate the expected value of higher order derivatives by using exclusively point evatuations of $f$ . This follows from

	$\frac{1}{N} N \sum i = 1 f (x_{i}) (- 1)^{k} \frac{\nabla^{k} p (x_{i})}{p (x_{i})}$	$\approx \int_{R^{d}} f (x) (- 1)^{k} \frac{\nabla^{k} p (x)}{p (x)} p (x) d x$
		$= \int_{R^{d}} \nabla^{k} f (x) d μ_{X} (x) = E_{x \sim μ_{X}} [\nabla^{k} f (x)]$
		$= m \sum i = 1 (\int_{R^{d}} g^{(k)} (a_{i}^{T} x) d μ_{X} (x)) a_{i} \otimes \dots \otimes a_{i}      k - times .$

In the work [32] decompositions of third order symmetric tensors ( $k = 3$ ) [1, 35, 51] have been used for the weights identification of one hidden layer neural networks. Instead, beyond the classical results about principal Hessian directions [37], in [21] it is shown that using second derivatives $(k = 2)$ actually suffices and the corresponding error estimates reflect positively the lower order and potential of improved stability, see e.g., [16, 28, 29]. The main part of the present work is an extension of the latter approach and therefore we will give a short summary of it with emphasis on active sampling, which will be assumed in this paper as the sampling method. The first step of the approach in [21] is taking advantage of (6) to reduce the dimensionality of the problem from $d$ to $m$ .

Reduction to the active subspace.

Before stating the core procedure, we want to introduce a simple and optional method, which can help to reduce the problem complexity in practice. Assume $f : Ω \subset R^{d} \to R$ takes the form (5), where $d \geq m$ and that $a_{1}, \dots, a_{m} \in R^{d}$ are linearly independent. From a numerical perspective the input dimension $d$ of the network plays a relevant role in terms of complexity of the procedure. For this reason in [21] the input dimension is effectively reduced to the number of neurons in the first hidden layer. With this reasoning, in the sections that follow we also consider networks where the input dimension matches the number of neurons of the first hidden layer.

Assume for the moment that the active subspace $A = span {a_{1}, \dots, a_{m}}$ is known. Let us choose any orthonormal basis of $A$ and arrange it as the columns of a matrix $^A \in R^{d \times m}$ . Then

f (x) = f (P_{A} x) = f (^A {^A}^{T} x),

which can be used to define a new network

^f (y) := f (^A y) : R^{m} \to R .

(7)

whose weights are $α_{1} = {^A}^{T} a_{1}, \dots, α_{m} = {^A}^{T} a_{m}$ , all the other parameters remain unchanged. Note that $^A α_{i} = P_{A} a_{i} = a_{i}$ , and therefore $a_{i}$ can be recovered from $α_{i}$ . In summary, if the active subspace of $f$ is approximately known, then we can construct $^f$ , such that the identification of $f$ and $^f$ are equivalent. This allows us to reduce the problem to the identification of $^f$ instead of $f$ , under the condition that we approximate $P_{A}$ well enough [21, Theorem 1.1]. As recalled in (6) we can produce easily approximations to vectors in $A$ by approximate first order differentiation of the original network $f$ and, in an ideal setting, generating $m$ linear independent gradients would suffices to approximate $A$ . However, in general, there is no way to ensure a priori such linear independence and we have to account for the error caused by approximating gradients by finite differences. By suitable assumptions on $f$ (see the full rank condition on the matrix $J [f]$ defined in (8) below) and using Algorithm 1 we obtain the following approximation result.

Input: Given a shallow neural network

f

as in (5), step-size of finite differences

ϵ > 0

, number of samples

m_{X}

1 begin

2 Draw

x_{1}, \dots, x_{m_{X}}

uniformly from the unit sphere

S^{d - 1}

3 Calculate the estimated gradients

Δ_{ϵ} f (x_{1}), \dots, Δ_{ϵ} f (x_{m_{X}})

by first order finite differences with stepsize

ϵ

4 Compute the singular value decomposition

(Δ_{ϵ} f (x_{1}) ∣ ∣ \dots ∣ ∣ Δ_{ϵ} f (x_{m_{X}})) = (\begin{matrix} {^U}_{1} & {^U}_{2} \end{matrix}) (\begin{matrix} {^Σ}_{1} & 0 0 & {^Σ}_{2} \end{matrix}) ⎛ ⎝ \begin{matrix} {^V}_{1}^{T} {^V}_{2}^{T} \end{matrix} ⎞ ⎠,

where

{^Σ}_{1}

contains the

m

largest singular values. Set

P_{^A} = {^U}_{1} {^U}_{1}^{T}

5 end

Output:

P_{^A}

Algorithm 1 Active subspace identification [21]

Theorem 1 ([21],Theorem 2.2).

Assume the vectors $(a_{i})_{i = 1}^{m}$ are linear independent and of unit norm. Additionally, assume that the $g_{i}$ ’s are smooth enough. Let $P_{^A}$ be constructed as described in Algorithm 1 by sampling $m_{X} (d + 1)$ values of $f$ . Let $0 < s < 1$ , and assume that the matrix

	$J [f]$	$:=$	$E_{X \sim μ_{S^{d - 1}}} \nabla f (X) \otimes \nabla f (X)$		(8)
		$=$			(8)

has full rank, i.e., its $m$ -th singular value fulfills $σ_{m} (J [f]) \geq α > 0$ . Then

∥ P_{A} - P_{^A} ∥_{F} \leq \frac{2 C_{1} ϵ m}{\sqrt{α (1 - s)} - C_{1} ϵ m},

with probability at least $1 - m exp (- \frac{m_{X} α s^{2}}{2 m^{2} C_{2}^{2}})$ , where $C_{1}, C_{2} > 0$ are absolute constants depending on the smoothness of $g_{i}$ ’s.

Identifying the weights.

As clarified in the previous section we can assume from now on that $d = m$ without loss of generality. Let $f$ be a network of the type (5), with three times differentiable activation functions $(g_{i})_{i = 1, \dots, m}$ , and independent weights $(a_{i})_{i = 1, \dots m} \in R^{m}$ of unit norm. Then $f$ has second derivative

\nabla^{2} f (x) = m \sum i = 1 g_{i}^{''} (a_{i}^{T} x) a_{i} \otimes a_{i} \in A = span {a_{1} \otimes a_{1}, \dots, a_{m} \otimes a_{m}},

(9)

whose expression represents a non-orthogonal rank-1 decomposition of the Hessian. The idea is, first of all, to modify the network by an ad hoclinear transformation (withening) of the input

f (W^{T} x) = m \sum i = 1 g_{i} (a_{i}^{T} W^{T} x)

(10)

in such a way that $(W a_{i} / ∥ W a_{i} ∥_{2})_{i = 1, \dots, m}$ forms an orthonormal system. The computation of $W$ can be performed by spectral decomposition of any positive definite matrix

G \in^A \approx A, γ I ≼ G .

In fact, from the spectral decomposition of $G = U D U^{T}$ , we define $W = D^{- \frac{1}{2}} U^{T}$ (see [21, Theorem 3.7]). This procedure is called whitening and allows to reduce the problem to networks with nearly-orthogonal weights, and presupposes to have obtained $^A≈A=span{a1⊗a1,…,am⊗am}$ . By using (9) and a similar approach as Algorithm 1 (one simply substitutes there the approximate gradients with vectorized approximate Hessians), one can compute $^A$ under the assumption that also the second order matrix

	$H [f]$	$:=$	$E_{X \sim μ_{S^{m - 1}}} vec (\nabla^{2} f (X)) \otimes vec (\nabla^{2} f (X))$
		$=$	$\int_{S^{m - 1}} vec (\nabla^{2} f (x)) \otimes vec (\nabla^{2} f (x)) d μ_{S^{m - 1}} (x)$

is of full rank, where $vec (\nabla^{2} f (x))$ is the vectorization of the Hessian $\nabla^{2} f (x)$ .

After whitening one could assume without loss of generality that the vectors $(a_{i})_{i = 1, \dots m} \in R^{m}$ are nearly orthonormal in the first place. Hence the representation (9) would be a near spectral decomposition of the Hessian and the components $a_{i} \otimes a_{i}$

would represent the approximate eigenvectors. However, the numerical stability of spectral decompositions is ensured only under spectral gaps

[50, 4]. In order to maximally stabilize the approximation of the

a_{i}

’s, one seeks for matrices

M \in^A

with the maximal spectral gap between the first and second largest eigenvalues. This is achieved by the maximizers of the following nonconvex program

M = arg max ∥ M ∥ s.t. M \in^A, {∥ M ∥}_{F} \leq 1,

(11)

where $∥ \cdot ∥$ and $∥ \cdot ∥_{F}$ are the spectral and Frobenius norms respectively. This program can be solved by a suitable projected gradient ascent, see for instance [21, Algorithm 3.4] and Algorithm 3 below, and any resulting maximizer has the eigenvector associated to the largest eigenvalue in absolute value close to one of the $a_{i}$ ’s. Once approximations ${^a}_{i}$ to all the $a_{i}$ ’s are retrieved, then it is not difficult to perform the identification of the activation functions $g_{i}$ , see [21, Algorithm 4.1, Theorem 4.1]. The recovery of the network resulting from this algorithmic pipeline is summarized by the following statement.

Theorem 2 ([21],Theorem 1.2).

Let $f$ be a real-valued function defined on the neighborhood of $Ω = B_{1}^{d}$ , which takes the form

f (x) = m \sum i = 1 g_{i} (a_{i} \cdot x),

for $m \leq d$ . Let $g_{i}$ be three times continuously differentiable on a neighborhood of $[- 1, 1]$ for all $i = 1, \dots, m$ , and let ${a_{1}, \dots, a_{m}}$ be linearly independent. We additionally assume both $J [f]$ and $H [f]$ of maximal rank $m$ . Then, for all $ϵ > 0$ (stepsize employed in the computation of finite differences), using at most $m_{X} [(d + 1) + (m + 1) (m + 2) / 2]$ random exact point evaluations of $f$ , the nonconvex program (11) constructs approximations ${^a1,…,^am}$ of the weights ${a_{1}, \dots, a_{m}}$ up to a sign change for which

(m \sum i = 1 ∥ {^a}_{i} - a_{i} ∥_{2}^{2})^{1 / 2} ≲ ε,

(12)

with probability at least $1 - m exp (- \frac{m_{X} c}{2 max {C_{1}, C_{2}}^{2} m^{2}})$ , for a suitable constant $c > 0$ intervening (together with some fixed power of $m$ ) in the asymptotical constant of the approximation (12). Moreover, once the weights are retrieved one constructs an approximating function $^f : B_{1}^{d} \to R$ of the form

^f (x) = m \sum i = 1 {^g}_{i} ({^a}_{i} \cdot x),

such that

∥ f -^f ∥_{C (Ω)} ≲ ϵ .

(13)

While this result have been generalized to the case of passive sampling in [21] and through whitening allows for the identification of non-orthogonal weights, it is restricted to the case of $m \leq d$ and linearly independent weights ${a_{i} : i = 1, \dots, m}$ .

The main goal of this paper is generalizing this approach to account for both the identification of two fully nonlinear hidden layer neural networks and the case where $m > d$ and the weights are not necessarily nearly orthogonal or even linearly independent (see Remark 2 below).

1.2.3 Deeper networks: the two layer case

What follows further extends the theory discussed in the previous sections to a wider class of functions, namely neural networks with two hidden layers. By doing so, we will also address a relevant open problem that was stated in [21], which deals with the identification of shallow neural networks where the number of neurons is larger than the input dimension. First, we need a precise definition of the architecture of the neural networks we intend to consider.

Definition 3.

Let $0 < m_{1} \leq m_{0} \leq d$ , and ${a_{1}, \dots, a_{m_{0}}} \subset R^{d}$ , ${b_{1}, \dots, b_{m_{1}}} \subset R^{m_{0}}$ be sets of unit vectors, and denote $A := [a_{1} | \dots | a_{m_{0}}] \in R^{d \times m_{0}}$ , $B := [b_{1} | \dots | b_{m_{1}}] \in R^{m_{0} \times m_{1}}$ . Let $g_{1}, \dots, g_{m_{0}}$ and $h_{1}, \dots, h_{m_{1}}$ be univariate functions, and denote $G_{0} = diag (g_{1}^{'} (0), \dots, g_{m_{0}}^{'} (0))$ . We define

F (d, m_{0}, m_{1}) := {f : R^{d} \to R : f (x) = m_{1} \sum ℓ = 1 h_{ℓ} (m_{0} \sum i = 1 b_{i ℓ} g_{i} (a_{i}^{T} x))},

(14)

with ${a_{1}, \dots, a_{m_{0}}} \subset R^{d}$ , ${b_{1}, \dots, b_{m_{1}}} \subset R^{m_{0}}$ , $g_{1}, \dots, g_{m_{0}}$ and $h_{1}, \dots, h_{m_{1}}$ satisfying

[label=(A0)]
$g_{i}^{'} (0) \neq 0 \forall i = 1, \dots, m_{0}$ ,
a frame condition for the system ${a_{1}, \dots, a_{m_{0}}, v_{1}, \dots, v_{m_{1}}}$ with $v_{ℓ} := \frac{A G_{0} b_{ℓ}}{∥ A G_{0} b_{ℓ} ∥}$ , i.e. there exist constants $c_{f}, C_{F} > 0$ such that

(15)

for all $x \in R^{d}$ ,
the derivatives of $g_{i}$ and $h_{ℓ}$ are uniformly bounded according to

$max i = 1, \dots, m_{0} sup t \in R ∣ ∣ g_{i}^{(k)} (t) ∣ ∣ \leq κ_{k}, and max i = 1, \dots, m_{1} sup t \in R ∣ ∣ h_{ℓ}^{(k)} (t) ∣ ∣ \leq η_{k}, k = 0, 1, 2, 3.$ (16)

Sometimes it may be convenient below the more compact writing $f (x) = 1^{T} h (B^{T} g (A^{T} x))$ where $g = (g_{1}, \dots, g_{m_{0}})$ , $h = (h_{1}, \dots, h_{m_{1}})$ . In the previous section we presented a dimension reduction that can be applied to one layer neural networks, and which can be useful to reduce the dimensionality from the input dimension to the number of neurons of the first layer. The same approach can be applied to networks defined by the class $F (d, m_{0}, m_{1})$ . For the approximation error of the active subspace, we end up with the following corollary of Theorem 1.

Corollary 1 (cf. Theorem 1).

Assume that $f \in F (d, m_{0}, m_{1})$ and let $P_{^A}$ be constructed as described in Algorithm 1 by sampling $m_{X} (d + 1)$ values of $f$ . Let $0 < s < 1$ , and assume that the $m_{0}$ -th singular value of $J [f]$ fulfills $σ_{m_{0}} (J [f]) \geq α > 0$ . Then we have

∥ P_{A} - P_{^A} ∥_{F} \leq \frac{2 C_{3} ϵ m_{0} m_{1}}{\sqrt{(1 - s) α} - C_{3} ϵ m_{0} m_{1}},

with probability at least $1 - m_{0} exp (- \frac{s^{2} m_{x} α}{2 C_{4} m_{1}})$ and constants $C_{3}, C_{4} > 0$ that depend only on $κ_{j}, η_{j}$ for $j = 0, \dots, 3$ .

From now on we assume $d = m_{0}$ .

2 Approximating the span of tensors of weights

In the one layer case, which was described earlier, the unique identification of the weights is made possible by constructing a matrix space whose rank-1 basis elements are outer products of the weight profiles of the network. This section illustrates the extension of this approach beyond shallow neural networks. Once again, we will make use of differentiation and overall there will be many parallels to the approach in [21]. However, the intuition behind the matrix space will be less straightforward, because we can not anymore directly express the second derivative of a two layer network as a linear combination of symmetric rank-1 matrices. This is due to the fact that the Hessian matrix of a network $f \in F (m_{0}, m_{0}, m_{1})$ has the form

	$\nabla^{2} f (x) =$	$m_{1} \sum ℓ = 1 h_{ℓ}^{'} (b_{ℓ}^{T} g (A^{T} x)) m_{0} \sum i = 1 b_{i ℓ} g_{i}^{''} (a_{i}^{T} x) a_{i} \otimes a_{i}$
	$+$	$m_{1} \sum ℓ = 1 m_{0} \sum i, j = 1 h_{ℓ}^{''} (b_{ℓ}^{T} g (A^{T} x)) b_{i ℓ} b_{j ℓ} g_{i}^{'} (a_{i}^{T} x) g_{j}^{'} (a_{j}^{T} x) (a_{i} \otimes a_{j} + a_{j} \otimes a_{i})$

Therefore, $\nabla^{2} f (x) \in span {a_{i} \otimes a_{j} + a_{j} \otimes a_{i} | i, j = 1, \dots, m_{0}}$ , which has dimension $\frac{m_{0} (m_{0} + 1)}{2}$

and is in general not spanned by symmetric rank-1 matrices. This expression is indeed quite complicated, due to the chain rule and the mixed tensor contributions, which are consequently appearing. At a first look, it would seem impossible to use a similar approach as the one for shallow neural networks recalled in the previous section. Nevertheless a relatively simple algebraic manipulation allows to recognize some useful structure: For a fixed

x \in R^{m_{0}}

we rearrange the expression as

	$\nabla^{2} f (x)$	$= m_{1} \sum ℓ = 1 h_{ℓ}^{'} (b_{ℓ}^{T} g (A^{T} x)) m_{0} \sum i = 1 b_{i ℓ} g_{i}^{''} (a_{i}^{T} x) a_{i} \otimes a_{i}$
		$+ m_{1} \sum ℓ = 1 h_{ℓ}^{''} (b_{ℓ}^{T} g (A^{T} x)) [m_{0} \sum i = 1 b_{i ℓ} g_{i}^{'} (a_{i}^{T} x) a_{i}] \otimes [m_{0} \sum j = 1 b_{j ℓ} g_{j}^{'} (a_{j}^{T} x) a_{j}],$

which is a combination of symmetric rank-1 matrices since $\sum_{j = 1}^{m_{0}} b_{j ℓ} g_{j}^{'} (a_{j}^{T} x) a_{j} \in R^{m_{0}}$ . We write the latter expression more compactly by introducing the notation

\nabla^{2} f (x) = m_{0} \sum i = 1 γ_{i} (x) a_{i} \otimes a_{i} + m_{1} \sum ℓ = 1 τ_{ℓ} (x) v_{ℓ} (x) \otimes v_{ℓ} (x),

(17)

where $G_{x} = diag (g_{1}^{'} (a_{1}^{T} x), \dots, g_{m_{0}}^{'} (a_{m_{0}}^{T} x)) \in R^{m_{0} \times m_{0}}$ and

$v_{ℓ} (x)$	$= A G_{x} b_{ℓ} \in R^{m_{0}}$	$for ℓ \in [m_{1}],$	(18)
$γ_{i} (x)$	$= g_{i}^{''} (a_{i}^{T} x) m_{1} \sum ℓ = 1 h_{ℓ}^{'} (b_{ℓ}^{T} g (A^{T} x)) b_{i ℓ} \in R$	$for i \in [m_{0}],$	(19)
$τ_{ℓ} (x)$	$= h_{ℓ}^{''} (b_{ℓ}^{T} g (A^{T} x)) \in R$	$for ℓ \in [m_{1}] .$	(20)

Figure 1: Illustration of the relationship between

W

(black line) and

span {\nabla^{2} f (x) ∣ ∣ x \in R^{m_{0}}}

(light blue region) given by two non-linear cones that fan out from

\nabla^{2} f (0)

. There is no reason to believe that the these cones are symmetric around

W

. The gray cones show the maximal deviation of

^W

from

W

Let us now introduce the fundamental matrix space

W = W (f) := span {a_{1} \otimes a_{1}, \dots, a_{m_{0}} \otimes a_{m_{0}}, v_{1} \otimes v_{1}, \dots, v_{m_{1}} \otimes v_{m_{1}}},

(21)

where $a_{1}, \dots, a_{m_{0}}$ are the weight profiles of the first layer and

v_{ℓ} := v_{ℓ} (0) / {∥ v_{ℓ} (0) ∥}_{2} = A G_{0} b_{ℓ} / ∥ A G_{0} b_{ℓ} ∥_{2},

for all $ℓ = 1, \dots, m_{1}$ encode entangled information about $b_{1}, \dots, b_{m_{1}}$ . For this reason, we call the $v_{ℓ}$ ’s entangled weights. Let us stress at this point that the definition and the constructive approximtion of the space $W$ is perhaps the most crucial and relevant contribution of this paper. In fact, by inspecting carefully the expression (17), we immediately notice that $\nabla^{2} f (0) \in W$ , and also that the first sum in (17), namely $\sum_{i = 1}^{m_{0}} β_{i} (x) a_{i} \otimes a_{i}$ , lies in $W$ for all $x \in R^{m_{0}}$ . Moreover, for arbitrary sampling points $x$ , deviations of $\nabla^{2} f (x)$ from $W$ are only due to the second term in (17). The intuition is that for suitable centered distributions of sampling points $x_{i}$ ’s so that $a_{j}^{T} x_{i} \approx 0$ and $G_{x_{i}} \approx G_{0}$ , the Hessians $(\nabla^{2} f (x_{i}))_{i = 1, \dots, m_{X}}$ will distribute themselves somehow around the space $W$ , see Figure 1 for a two dimensional sketch of the geometrical situation. Hence, we would attempt an approximation of $W$ by PCA of a collection of such approximate Hessians. Practically, by active sampling (targeted evaluations of the network $f$ ) we first construct estimates $(Δ_{ϵ}^{2} f (x_{i}))_{i = 1, \dots, m_{X}}$ by finite differences of the Hessian matrices $(\nabla^{2} f (x_{i}))_{i = 1, \dots, m_{X}}$ (see Section 2.1), at sampling points $x_{1}, \dots, x_{m_{X}} \in R^{m}$ drawn independently from a suitable distribution $μ_{X}$ . Next, we define the matrix

^W = (vec (Δ_{ϵ}^{2} f (x_{1})), \dots, vec (Δ_{ϵ}^{2} f (x_{m_{X}}))),

whose columns are the vectorization of the approximate Hessians. Finally, we produce the approximation $^W$ to $W$ as the span of the first $m_{0} + m_{1}$ left singular vectors of the matrix $^W$ . The whole procedure of calculating $^W$ is given in Algorithm 2. It should be clear that the choice of $μ_{X}$ plays a crucial role for the quality of this method. In the analysis that follows, we focus on distributions that are centered and concentrated. Figure 1 helps to form a better geometrical intuition of the result of the procedure. It shows the region covered by the Hessians, indicated by the light blue area, which envelopes the space $W$ in a sort of nonlinear/nonconvex cone originating from $\nabla^{2} f (0)$ . In general, the Hessians do not concentrate around $W$ in a symmetric way, which means that the “center of mass” of the Hessians can never be perfectly aligned with the space $W$ , regardless of the number of samples. In this analogy, the center of mass is equivalent to the space estimated by Algorithm 2

, which essentially is a non-centered principal component analysis of observed Hessian matrices. The primary result of this section is Theorem

4, which provides an estimate of the approximation error of Algorithm 2 depending on the subgaussian norm of the sample distribution

μ_{X}

and the number of neurons in the respective layers. More precisely, this result gives a precise worst case estimate of the error caused by the imbalance of mass. For reasons mentioned above, the error does not necessarily vanish with an increasing number of samples, but the probability under which the statement holds will tend to

1

. In Figure 1, the estimated region is illustrated by the gray cones that show the maximal, worst case deviation of

^W

. One crucial condition for Theorem 4 to hold is that there exists an

α > 0

such that

σ_{m_{0} + m_{1}} (E_{X \sim μ_{X}} v e c (\nabla^{2} f (X)) \otimes vec (\nabla^{2} f (X))) \geq α .

(22)

This assumption makes sure that the space spanned by the observed Hessians has, in expectation, at least dimension $m_{0} + m_{1}$ . Aside from this technical aspect this condition implicitly helps to avoid network configurations, which are reducible, for certain weights can not be recovered. For example, we can define a network in $F (2, 2, 1)$ with weights given by

a_{1} = ⎛ ⎜ ⎝ \begin{matrix} \frac{1}{\sqrt{2}} \frac{1}{\sqrt{2}} \end{matrix} ⎞ ⎟ ⎠, a_{2} = ⎛ ⎜ ⎝ \begin{matrix} \frac{1}{\sqrt{2}} - \frac{1}{\sqrt{2}} \end{matrix} ⎞ ⎟ ⎠, b_{1} = (\begin{matrix} 10 \end{matrix}) .

It is easy to see that $a_{2}$ will never be used during a forward pass through the network, which makes it impossible to recover $a_{2}$ from the output of the network.
In the theorem below and in the proofs that follow we will make use of the subgaussian norm ${∥ \cdot ∥}_{ψ_{2}}$

of a random variable. This quantity measures how fast the tails of a distribution decay and such a decay plays an important role in several concentration inequalities. More in general, for

p \geq 1

, the

ψ_{p}

-norm of a scalar random variable

Z

is defined as

{∥ Z ∥}_{ψ_{p}} = inf {t > 0 : E exp (| Z / t |^{p}) \leq 2} .

For a random vector $X$ on $R^{d}$ the $ψ_{p}$ -norm is given by

{∥ X ∥}_{ψ_{p}} = sup x \in S^{d - 1} {∥ | ⟨ X, x ⟩ | ∥}_{ψ_{p}} .

The random variables for which ${∥ X ∥}_{ψ_{1}} < \infty$ are called subexponential and those for which ${∥ X ∥}_{ψ_{2}} < \infty$ are called subgaussian. More in general, the Orlicz space $L_{ψ_{p}} = L_{ψ_{p}} (Ω, Σ, P)$ consists of all real random variables $X$ on the probabillity space $(Ω, Σ, P)$ with finite ${∥ X ∥}_{ψ_{p}}$ norm and its elements are called $p$ -subexponential random variagles. Below, we mainly focus on subgaussian random variables. In particular, every bounded random variable is subgaussian, which covers all the cases we discuss in this work. We refer to [66]

for more details. One example of a subgaussian distribution is the uniform distribution on the unit-sphere, which has subgaussian norm

{∥ X ∥}_{ψ_{2}} = \frac{1}{\sqrt{d}}, X \sim Unif (S^{d - 1})

Input: Neural network

f

, number of estimated Hessians

m_{X}

, step-size of the finite difference approximation

ϵ > 0

, probability distribution

μ_{X}

1 begin

2 Draw

x_{1}, \dots x_{m_{X}}

independently from

μ_{X}

3 Calculate the matrix

M = (vec (Δ_{ϵ}^{2} f (x_{1})) | \dots | v e c (Δ_{ϵ}^{2} f (x_{m_{X}})))

4 Set

U Σ V^{T} = SVD (M)

5 Denote by

U_{1}

the first

m_{0} + m_{1}

columns of

U

P_{^W} = (u_{1} | \dots | u_{m_{0} + m_{1}}) (u_{1} | \dots | u_{m_{0} + m_{1}})^{T}

7 end

Output:

P_{^W}

Algorithm 2 Approximating

W

Theorem 4.

Let $f \in F (m_{0}, m_{0}, m_{1})$ be a neural network within the class described in Definition 3 and consider the space $W$ as defined in (21). Assume that $μ_{X}$ is a probability measure with $s u p p (μ_{X}) \subset B_{1}^{m_{0}}$ , $E X = 0$ , and that there exists an $α > 0$ such that

σ_{m_{0} + m_{1}} (E_{X \sim μ_{X}} v e c (\nabla^{2} f (X)) \otimes vec (\nabla^{2} f (X))) \geq α .

(23)

Then, for any $ϵ > 0$ , Algorithm 2 returns a projection $P_{^W}$ that fulfills

(24)

for a suitable subspace $W^{*} \subset W$ (we can actually assume that $W^{*} = W$ according to Remark 1 below) with probability at least

1 - 2 e^{- c_{1} L_{1} m_{X} min {L_{1} m_{X}, 1}} - (m_{0} + m_{1}) e^{- L_{2} m_{X}}

where

L_{1}

:= {∥ X ∥}_{ψ_{2}}^{2} log (m_{0} + 1) m_{1}, L_{2} := α (8 C_{1} {∥ A ∥}^{4} {∥ B ∥}^{4} m_{1})^{- 1},

$c_{1}, c_{2}$ are absolute constants and $C, C_{1}, C_{Δ} > 0$ are constants depending on the constants $κ_{j}, η_{j}$ for $j = 0, \dots, 3$ .

Remark 1.

If $ϵ > 0$ is sufficiently small, due to (23) the space $^W$ returned by Algorithm 2 has dimension $m_{0} + m_{1}$ . If the error bound (24) in Theorem 4 is such that ${∥ ∥ P_{W^{*}} - P_{^W} ∥ ∥}_{F} < 1$ , then $^W$ and $W^{*}$ must have the same dimension. Moreover, $W^{*} \subset W$ and $dim (W) = m_{0} + m_{1}$ would necessarily imply that $W = W^{*}$ . Hence, for ${∥ ∥ P_{W^{*}} - P_{^W} ∥ ∥}_{F} < 1$ and $ϵ > 0$ sufficiently small, we have $W^{*} = W$ .

As already mentioned above, for $μ_{X} = Unif (S^{m_{0} - 1})$ we have ${∥ X ∥}_{ψ_{2}} = \frac{1}{\sqrt{m_{0}}}$ . In this case the error bound 24 behaves like

∥

Robust and Resource Efficient Identification of Two Hidden Layer Neural Networks

Authors

Stable Recovery of Entangled Weights: Towards Robust Identification of Deep Neural Networks from Minimal Samples

Neural networks with trainable matrix activation functions

Jacobi-type algorithm for low rank orthogonal approximation of symmetric tensors and its convergence analysis

Identification of Shallow Neural Networks by Fewest Samples

On the Complexity of Learning Neural Networks

A Unified and Constructive Framework for the Universality of Neural Networks

1 Introduction

1.1 Notation

1.2 From one artificial neuron to shallow, and deeper networks

1.2.1 Meet the neuron

1.2.2 Shallow networks: the one-layer case

Reduction to the active subspace.

Theorem 1 ([21],Theorem 2.2).

Identifying the weights.

Theorem 2 ([21],Theorem 1.2).

1.2.3 Deeper networks: the two layer case

Definition 3.

Corollary 1 (cf. Theorem 1).

2 Approximating the span of tensors of weights

Theorem 4.

Remark 1.

Robust and Resource Efficient Identification of Two Hidden Layer Neural Networks

Authors

Related Research

Stable Recovery of Entangled Weights: Towards Robust Identification of Deep Neural Networks from Minimal Samples

Neural networks with trainable matrix activation functions

Jacobi-type algorithm for low rank orthogonal approximation of symmetric tensors and its convergence analysis

Identification of Shallow Neural Networks by Fewest Samples

On the Complexity of Learning Neural Networks

A Unified and Constructive Framework for the Universality of Neural Networks

This week in AI

1 Introduction

1.1 Notation

1.2 From one artificial neuron to shallow, and deeper networks

1.2.1 Meet the neuron

1.2.2 Shallow networks: the one-layer case

Reduction to the active subspace.

Theorem 1 ([21],Theorem 2.2).

Identifying the weights.

Theorem 2 ([21],Theorem 1.2).

1.2.3 Deeper networks: the two layer case

Definition 3.

Corollary 1 (cf. Theorem 1).

2 Approximating the span of tensors of weights

Theorem 4.

Remark 1.