Year: 2020

Abstract:
The tremendous success of deep neural networks has motivated the need to better understand the fundamental properties of these networks, but many of the theoretical results proposed have only been for shallow networks. In this paper, we study an important primitive for understanding the meaningful input space of a deep network: span recovery. For k < n, let A ∈ Rk×n be the innermost weight matrix of an arbitrary feed forward neural network M : R → R, so M(x) can be written as M(x) = σ(Ax), for some network σ : R → R. The goal is then to recover the row span of A given only oracle access to the value of M(x). We show that if M is a multi-layered network with ReLU activation functions, then partial recovery is possible: namely, we can provably recover k/2 linearly independent vectors in the row span of A using poly(n) non-adaptive queries to M(x). Furthermore, if M has differentiable activation functions, we demonstrate that full span recovery is possible even when the output is first passed through a sign or 0/1 thresholding function; in this case our algorithm is adaptive. Empirically, we confirm that full span recovery is not always possible, but only for unrealistically thin layers. For reasonably wide networks, we obtain full span recovery on both random networks and networks trained on MNIST data. Furthermore, we demonstrate the utility of span recovery as an attack by inducing neural networks to misclassify data obfuscated by controlled random noise as sensical inputs.

heading: 1 INTRODUCTION
text: Consider the general framework in which we are given an unknown function f : Rn → R, and we want to learn properties about this function given only access to the value f(x) for different inputs x. There are many contexts where this framework is applicable, such as blackbox optimization in which we are learning to optimize f(x) (Djolonga et al., 2013), PAC learning in which we are learning to approximate f(x) (Denis, 1998), adversarial attacks in which we are trying to find adversarial inputs to f(x) (Szegedy et al., 2013), or structure recovery in which we are learning the structure of f(x). For example in the case when f(x) is a neural network, one might want to recover the underlying weights or architecture (Arora et al., 2014; Zhang et al., 2017). In this work, we consider the setting when f(x) = M(x) is a neural network that admits a latent low-dimensional structure, namely M(x) = σ(Ax) where A ∈ Rk×n is a rank k matrix for some k < n, and σ : Rk → R is some neural network. In this setting, we focus primarily on the goal of recovering the row-span of the weight matrix A. We remark that all our results generalize in a straightforward manner to the case when A is rank r < k.
Span recovery of general functions f(x) = g(Ax), where g is arbitrary, has been studied in some contexts, and is used to gain important information about the underlying function f . By learning Span(A), we in essence are capturing the relevant subspace of the input to f ; namely, f behaves identically on x as it does on the projection of x onto the row-span of A. In statistics, this is known as effective dimension reduction or the multi-index model Li (1991); Xia et al. (2002). Another important motivation for span recovery is for designing adversarial attacks. Given the span of A, we compute the kernel of A, which can be used to fool the function into behaving incorrectly on inputs which are perturbed by vectors in the kernel. Specifically, if x is a legitimate input correctly classified by f and y is a large random vector in the kernel of A, then x+ y will be indistinguishable from noise but we will have f(x) = f(x+ y).
Several works have considered the problem from an approximation-theoretic standpoint, where the goal is to output a hypothesis function f̃ which approximates f well on a bounded domain. For instance, in the case that A ∈ Rn is a rank 1 matrix and g(Ax) is a smooth function with bounded derivatives, Cohen et al. (2012) gives an adaptive algorithm to approximate f . Their results also give an approximation Ã to A, under the assumption that A is a stochastic vector (Ai ≥ 0 for each i and ∑ iAi = 1). Extending this result to more general rank k matrices A ∈ Rk×n, Tyagi & Cevher (2014) and Fornasier et al. (2012) give algorithms with polynomial sample complexity to find approximations f̃ to twice differentiable functions f . However, their results do not provide any guarantee that the original matrix A itself or a good approximation to its span will be recovered. Specifically, the matrix Ã used in the hypothesis function f̃(x) = g̃(Ãx) of Tyagi & Cevher (2014) only has moderate correlation with the true row span of A, and always admits some constant factor error (which can translate into very large error in any subspace approximation).
Furthermore, all aforementioned works require the strong assumption that the matrix of gradients is well-conditioned (and full rank) in order to obtain good approximations f̃ . In contrast, when f(x) is a non-differentiable ReLU deep network with only mild assumptions on the weight matrices, we prove that the gradient matrix has rank at least k/2, which significantly strengthens span recovery guarantees since we do not make any assumptions on the gradient matrix. Finally, Hardt & Woodruff (2013) gives an adaptive approximate span recovery algorithm with poly(n) samples under the assumption that the function g satisfies a norm-preserving condition, which is restrictive and need not (and does not) hold for the deep neural networks we consider here.
On the empirical side, the experimental results of Tyagi & Cevher (2014) for function approximation were only carried out for simple one-layer functions, such as the logistic function in one dimension (where A has rank k = 1), and on linear functions gT (Ax+ b), where g ∈ Rk has i.i.d. Gaussian entries. Moreover, their experiments only attempted to recover approximations Ã to A when A was orthonormal. In addition, Fornasier et al. (2012) experimentally considers the approximation problem for f when f(Ax) is a third degree polynomial of the input. This leaves an experimental gap in understanding the performance of span recovery algorithms on non-smooth, multi-layer deep neural networks.
When f(x) is a neural network, there have been many results that allow for weight or architecture recovery under additional assumptions; however nearly all such results are for shallow networks. Arora et al. (2014) shows that layerwise learning can recover the architecture of random sparse neural networks. Janzamin et al. (2015) applies tensor methods to recover the weights of a two-layer neural network with certain types of smooth activations and vectorvalued output, whereas Ge et al. (2019); Bakshi et al. (2019) obtain weight recovery for ReLU activations. Zhang et al. (2017) shows that SGD can learn the weights of two-layer neural networks with some specific activations. There is also a line of work for improperly two-layer networks, where the algorithm outputs an arbitrary hypothesis function which behaves similarly to the network under a fixed distribution (Goel et al., 2017; Goel & Klivans, 2019).
Learning properties of the network can also lead to so-called model extraction attacks or enhance classical adversarial attacks on neural networks (Jagielski et al., 2019). Adversarial attacks are often differentiated into two settings, the whitebox setting where the trained network weights are known or the blackbox setting where the network weights are unknown but attacks can still be achieved on external information, such as knowledge of the dataset, training algorithm, network architecture, or network predictions. Whitebox attacks are well-studied and usually use explicit gradients or optimization procedures to compute adversarial inputs for various tasks such as classification and reinforcement learning (Szegedy et al., 2013; Huang et al., 2017; Goodfellow et al., 2014). However, blackbox attacks are more realistic and it is clear that model recovery can enhance these attacks. The work of Papernot et al. (2017) attacks practical neural networks upon observations of predictions of adaptively chosen inputs, trains a substitute neural network on the observed data, and applies a whitebox attack on the substitute. This setting, nicknamed the practical blackbox setting (Chen et al., 2017), is what we work in, as we only observe adaptively chosen predictions without knowledge of the network architecture, dataset, or algorithms. We note that perhaps surprisingly, some of our algorithms are in fact entirely non-adaptive.

heading: 1.1 OUR CONTRIBUTIONS
text: In this paper, we provably show that span recovery for deep neural networks with high precision can be efficiently accomplished with poly(n) function evaluations, even when the networks have poly(n) layers and the output of the network is a scalar in some finite set. Specifically, for deep networks M(x) : Rn → R with ReLU activation functions, we prove that we can recover a subspace V ⊂ Span(A) of dimension at least k/2 with polynomially many non-adaptive queries.1 First, we use a volume bounding technique to show that a ReLU network has sufficiently large piece-wise linear sections and that gradient information can be derived from function evaluations. Next, by using a
1We switch to the notation M(x) instead of f(x) to illustrate that M(x) is a neural network, whereas f(x) was used to represent a general function with low-rank structure.
novel combinatorial analysis of the sign patterns of the ReLU network along with facts in polynomial algebra, we show that the gradient matrix has sufficient rank to allow for partial span recovery.
Theorem 3.4 (informal) Suppose we have the network M(x) = wTφ(W1φ(W2φ(. . .Wdφ(Ax)) . . . ), where A ∈ Rk×n is rank k, φ is the ReLU and Wi ∈ Rki×ki+1 are weight matrices, with ki possibly much smaller than k. Then, under mild assumptions, there is a non-adaptive algorithm that makesO(kn log k) queries toM(x) and returns in poly(n, k)-time a subspace V ⊆ span(A) of dimension at least k2 with probability 1− δ. We remark that span recovery of the first weight layer is provably feasible even in the surprising case when the neural network has many “bottleneck" layers with small O(log(n)) width. Because this does not hold in the linear case, this implies that the non-linearities introduced by activations in deep learning allow for much more information to be captured by the model. Moreover, our algorithm is non-adaptive, which means that the points xi at whichM(xi) needs to be evaluated can be chosen in advance and span recovery will succeed with high probability. This has the benefit of being parallelizable, and possibly more difficult to detect when being used for an adversarial attack. In addition, we note that this result generalize to the case when A is rank r < k, in which setting our guarantee will instead be that we recover a subspace of dimension at least r2 contained within the span of A.
In contrast with previous papers, we do not assume that the gradient matrix has large rank; rather our main focus and novelty is to prove this statement under minimal assumptions. We require only two mild assumptions on the weight matrices. The first assumption is on the orthant probabilities of the matrix A, namely that the distribution of sign patterns of a vector Ag, where g ∼ N (0, In), is not too far from uniform. Two examples of matrices which satisfy this property are random matrices and matrices with nearly orthogonal rows. The second assumption is a non-degeneracy condition on the matrices Wi, which enforces that products of rows of the matrices Wi result in vectors with non-zero coordinates.
Our next result is to show that full span recovery is possible for thresholded networks M(x) with twice differentiable activation functions in the inner layers, when the network has a 0/1 threshold function in the last layer and becomes therefore non-differentiable, i.e., M(x) ∈ {0, 1}. Since the activation functions can be arbitrarily non-linear, our algorithm only provides an approximation of the true subspace Span(A), although the distance between the subspace we output and Span(A) can be made exponentially small. We need only assume bounds on the first and second derivatives of the activation functions, as well as the fact that we can find inputs x ∈ Rn such that M(x) 6= 0 with good probability, and that the gradients of the network near certain points where the threshold evaluates to one are not arbitrarily small. We refer the reader to Section 4 for further details on these assumptions. Under these assumptions, we can apply a novel gradient-estimation scheme to approximately recover the gradient of M(x) and the span of A.
Theorem 4.3 (informal) Suppose we have the network M(x) = τ(σ(Ax)), where τ : R → {0, 1} is a threshold function and σ : Rk → R is a neural network with twice differentiable activation functions, and such that M satisfies the conditions sketched above (formally defined in Section 4). Then there is an algorithm that runs in poly(N) time, making at most poly(N) queries to M(x), where N = poly(n, k, log( 1 ), log(
1 δ )), and returns with probability 1 − δ
a subspace V ⊂ Rn of dimension k such that for any x ∈ V , we have
‖PSpan(A)x‖2 ≥ (1− )‖x‖2
where PSpan(A) is the orthogonal projection onto the span of A.
Empirically, we verify our theoretical findings by running our span recovery algorithms on randomly generated networks and trained networks. First, we confirm that full recovery is not possible for all architectures when the network layer sizes are small. This implies that the standard assumption that the gradient matrix is full rank does not always hold. However, we see that realistic network architectures lend themselves easily to full span recovery on both random and trained instances. We emphasize that this holds even when the network has many small layers, for example a ReLU network that has 6 hidden layers with [784, 80, 40, 30, 20, 10] nodes, in that order, can still admit full span recovery of the rank 80 weight matrix.
Furthermore, we observe that we can effortlessly apply input obfuscation attacks after a successful span recovery and cause misclassifications by tricking the network into classifying noise as normal inputs with high confidence. Specifically, we can inject large amounts of noise in the null space of A to arbitrarily obfuscate the input without changing the output of the network. We demonstrate the utility of this attack on MNIST data, where we use span recovery to generate noisy images that are classified by the network as normal digits with high confidence. We note that this veers away from traditional adversarial attacks, which aim to drastically change the network output with humanly-undetectable changes in the input. In our case, we attempt the arguably more challenging problem of drastically changing the input without affecting the output of the network.

heading: 2 PRELIMINARIES
text: Notation For a vector x ∈ Rk, the sign pattern of x, denoted sign(x) ∈ {0, 1}k, is the indicator vector for the nonzero coordinates of x. Namely, sign(x)i = 1 if xi 6= 0 and sign(x)i = 0 otherwise. Given a matrix A ∈ Rn×m, we denote its singular values as σmin(A) = σmin{n,m}, . . . , σ1(A) = σmax(A). The condition number of A is denoted κ(A) = σmax(A)/σmin(A). We let In ∈ Rn×n denote the n× n identity matrix. For a subspace V ⊂ Rn, we write PV ∈ Rn×n to denote the orthogonal projection matrix onto V . If µ ∈ Rn and Σ ∈ Rn×n is a PSD matrix, we write N (µ,Σ) to denote the multi-variate Gaussian distribution with mean µ and covariance Σ.
Gradient Information For any function f(x) = g(Ax), note that ∇f(x) = A>g(Ax) must be a vector in the row span of A. Therefore, span recovery boils down to understanding the span of the gradient matrix as x varies. Specifically, note that if we can find points x1, .., xk such that {∇f(xi)} are linearly independent, then the full span of A can be recovered using the span of the gradients. To our knowledge, previous span recovery algorithms heavily rely on the assumption that the gradient matrix is full rank and in fact well-conditioned. Specifically, for some distribution D, it is assumed that Hf = ∫ x∼D∇f(x)∇f(x)
> dx is a rank k matrix with a minimum non-zero singular value bounded below by α and the number of gradient or function evaluations needed depends inverse polynomially in α. In contrast, in this paper, when f(x) is a neural network, we provably show that Hf is a matrix of sufficiently high rank or large minimum non-zero singular value under mild assumptions, using tools in polynomial algebra.

heading: 3 DEEP NETWORKS WITH RELU ACTIVATIONS
text: In this section, we demonstrate that partial span recovery is possible for deep ReLU networks. Specifically, we consider neural networks M(x) : Rn → R of the form
M(x) = wTφ(W1φ(W2φ(. . .Wdφ(Ax)) . . . )
where φ(x)i = max{xi, 0} is the RELU (applied coordinate-wise to each of its inputs), and Wi ∈ Rki×ki+1 , and w ∈ Rkd , and A has rank k. We note that ki can be much smaller than k. In order to obtain partial span recovery, we make the following assumptions parameterized by a value γ > 0 (our algorithms will by polynomial in 1/γ):
• Assumption 1: For every sign pattern S ∈ {0, 1}k, we have Prg∼N (0,In) [sign(φ(Ag)) = S] ≥ γ/2k. • Assumption 2: For any S1, . . . , Sd 6= ∅ where Si ⊆ [ki], we have wT (∏d i=1(Wi)Si ) ∈ Rk is entry-wise
non-zero. Here (Wi)Si is the matrix with the rows j /∈ Si set equal to 0. Moreover, we assume
Prg∼N (0,Ik) [M(g) = 0] ≤ γ
8 .
Our first assumption is an assumption on the orthant probabilities of the distribution Ag. Specifically, observe that Ag ∈ Rk follows a multi-variate Gaussian distribution with covariance matrix AAT . Assumption 1 then states that the probability that a random vector x ∼ N (0,AAT ) lies in a certain orthant of Rk is not too far from uniform. We remark that orthant probabilities of multivariate Gaussian distributions are well-studied (see e.g., Miwa et al. (2003); Bacon (1963); Abrahamson et al. (1964)), and thus may allow for the application of this assumption to a larger class of matrices. In particular, we show it is satisfied by both random matrices and orthogonal matrices. Our second assumption is a non-degeneracy condition on the weight matrices Wi – namely, that products of wT with non-empty sets of rows of the Wi result in entry-wise non-zero vectors. In addition, Assumption 2 requires that the network is non-zero with probability that is not arbitrarily small, otherwise we cannot hope to find even a single xwithM(x) 6= 0. In the following lemma, we demonstrate that these conditions are satisfied by randomly initialized networks, even when the entries of the Wi are not identically distributed. Lemma 3.1. If A ∈ Rk×n is an arbitrary matrix with orthogonal rows, or if n > Ω(k3) and A has entries that are drawn i.i.d. from some sub-Gaussian distribution D with expectation 0, unit variance, and constant sub-Gaussian norm ‖D‖ψ2 = supp≥1 p−1/2 (EX∼D|X|p)
1/p then with probability at least 1− e−k2 , A satisfies Assumption 1 with γ ≥ 1/2. Moreover, if the weight matrices w,W1,W2, ...,Wd with Wi ∈ Rki×ki+1 have entries that are drawn independently (and possibly non-identically) from continuous symmetric distributions, and if ki ≥ log( 16dδγ ) for each i ∈ [d], then Assumption 2 holds with probability 1− δ.

heading: 3.1 ALGORITHM FOR SPAN RECOVERY
text: The algorithm for recovery is given in Algorithm 1. Our algorithm computes the gradient ∇M(gi) for different Gaussian vectors gi ∼ N (0, Ik), and returns the subspace spanned by these gradients. To implement this procedure,
Algorithm 1: Span Recovery with Non-Adaptive Gradients Input: function M(x) : Rn → R, k: latent dimension , γ: probability parameter
1 Set r = O(k log(k)/γ) 2 for i = 0, . . . , r do 3 Generate random Gaussian vector gi ∼ N (0, In). 4 Compute Gradient: zi = ∇M(gi) . Lemma 3.2 5 end 6 return Span ( z1, z2, . . . , zr
) we must show that it is possible to compute gradients via the perturbational method (i.e. finite differences), given only oracle queries to the network M . Namely, we firstly must show that if g ∼ N (0, In) then ∇M(g) exists, and moreover, that ∇M(x) exists for all x ∈ B (g), where B (g) is a ball of radius centered at g, and is some value with polynomial bit complexity which we can bound. To demonstrate this, we show that for any fixing of the sign patterns of the network, we can write the region of Rn which satisfies this sign pattern and is -close to one of the O(dk) ReLU thresholds of the network as a linear program. We then show that the feasible polytope of this linear program is contained inside a Euclidean box in Rn, which has one side of length . Using this containment, we upper bound the volume of the polytope in Rn which is close to each ReLU, and union bound over all sign patterns and ReLUs to show that the probability that a Gaussian lands in one of these polytopes is exponentially small. Lemma 3.2. There is an algorithm which, given g ∼ N (0, Ik), with probability 1− exp(−nc) for any constant c > 1 (over the randomness in g), computes∇M(g) ∈ Rn with O(n) queries to the network, and in poly(n) runtime.
Now observe that the gradients of the network lie in the row-span of A. To see this, for a given input x ∈ Rn, let S0(x) ∈ Rk be the sign pattern of φ(Ax) ∈ Rk, and more generally define Si(x) ∈ Rki via
Si(x) = sign ( φ(Wiφ(Wi+1φ(. . .Wdφ(Ax)) . . . ) ) Then ∇M(x) = (wT · ( ∏d i=1(Wi)Si)))AS0 , which demonstrates the claim that the gradients lie in the row-span of A. Now define zi = ∇M(gi) where gi ∼ N (0, In), and let Z be the matrix where the i-th row is equal to zi. We will prove that Z has rank at least k/2. To see this, first note that we can write Z = VA, where V is some matrix such that the non-zero entries in the i-th row are precisely the coordinates in the set Si0, where S i j = Sj(g
i) for any j = 0, 1, 2, . . . , d and i = 1, 2, . . . , r. We first show that V has rank at least ck for a constant c > 0. To see this, suppose we have computed r gradients so far, and the rank of V is less than ck for some 0 < c < 1/2. Now V ∈ Rr×k is a fixed rank-ck matrix, so the span of the matrix can be expressed as a linear combination of some fixed subset of ck of its rows. We use this fact to show in the following lemma that the set of all possible sign patterns obtainable in the row span of V is much smaller than 2k. Thus a gradient zr+1 with a uniform (or nearly uniform) sign pattern will land outside this set with good probability, and thus will increase the rank of Z when appended. Lemma 3.3. Let V ∈ Rr×k be a fixed at most rank ck matrix for c ≤ 1/2. Then the number of sign patterns S ⊂ [k] with at most k/2 non-zeros spanned by the rows of V is at most 2
k √ k . In other words, the set S(V) = {sign(w) | w ∈ span(V), nnz(w) ≤ k2} has size at most 2k√ k .
Theorem 3.4. Suppose the network M(x) = wTφ(W1φ(W2φ(. . .Wdφ(Ax)) . . . ), where φ is the ReLU, satisfies Assumptions 1 and 2. Then the algorithm given in Figure 1 makes O(kn log(k/δ)/γ) queries to M(x) and returns in poly(n, k, 1/γ, log(1/δ)) time a subspace V ⊆ span(A) of dimension at least k2 with probability 1− δ.

heading: 4 NETWORKS WITH THRESHOLDING ON DIFFERENTIABLE ACTIVATIONS
text: In this section, we consider networks that have a threshold function at the output node, as is done often for classification. Specifically, let τ : R → {0, 1} be the threshold function: τ(x) = 1 if x ≥ 1, and τ(x) = 0 otherwise. Again, we let A ∈ Rk×n where k < n, be the innermost weight matrix. The networks M : Rn → R we consider are then of the form:
M(x) = τ(W1φ1(W2φ2(. . . φdAx)) . . . )
where Wi ∈ Rki×ki+1 and each φi is a continuous, differentiable activation function applied entrywise to its input. We will demonstrate that even for such functions with a binary threshold placed at the end, giving us minimal information about the network, we can still achieve full span recovery of the weight matrix A, albeit with the cost of an - approximation to the subspace. Note that the latter fact is inherent, since the gradient of any function that is not linear
in some ball around each point cannot be obtained exactly without infinitely small perturbations of the input, which we do not allow in our model.
We can simplify the above notation, and write σ(x) = W1φ1(W2φ2(. . . φdAx)) . . . ), and thus M(x) = τ(σ(x)). Our algorithm will involve building a subspace V ⊂ Rn which is a good approximation to the span of A. At each step, we attempt to recover a new vector which is very close to a vector in Span(A), but which is nearly orthogonal to the vectors in V . Specifically, after building V , on an input x ∈ Rn, we will query M for inputs M((In −PV )x). Recall that PV is the projection matrix onto V , and PV ⊥ is the projection matrix onto the subspace orthogonal to V . Thus, it will help here to think of the functions M,σ as being functions of x and not (In − PV )x, and so we define σV (x) = σ(A(In−PV )x), and similarly MV (x) = τ(σV (x)). For the results of this section, we make the following assumptions on the activation functions.
Assumptions:
1. The function φi : R→ R is continuous and twice differentiable, and φi(0) = 0. 2. φi and φ′i are Li-Lipschitz, meaning:
sup x∈R ∣∣∣∣ ddxφi(x) ∣∣∣∣ ≤ Li, sup x∈R ∣∣∣∣ d2d2xφi(x) ∣∣∣∣ ≤ Li
3. The network is non-zero with bounded probability: for every subspace V ⊂ Rn of dimension dim(V ) < k, we have that Prg∼N (0,In) [σV (g) ≥ 1] ≥ γ for some value γ > 0. 4. Gradients are not arbitrarily small near the boundary: for every subspace V ⊂ Rn of dimension dim(V ) < k
Prg∼N (0,In) [ |∇gσV (cg)| ≥ η,∀c > 0 such that σV (cg) = 1, and σV (g) ≥ 1] ≥ γ
for some values η, γ > 0, where ∇gσV (cg) is the directional derivative of σV in the direction g.
The first two conditions are standard and straightforward, namely φi is differentiable, and has bounded first and second derivatives (note that for our purposes, they need only be bounded in a ball of radius poly(n)). SinceM(x) is a threshold applied to σ(x), the third condition states that it is possible to find inputs x with non-zero network evaluation M(x). Our condition is slightly stronger, in that we would like this to be possible even when x is projected away from any k′ < k dimensional subspace (note that this ensures that Ax is non-zero, since A has rank k).
The last condition simply states that if we pick a random direction g where the network is non-zero, then the gradients of the network are not arbitrarily small along that direction at the threshold points where σ(c · g) = 1. Observe that if the gradients at such points are vanishingly small, then we cannot hope to recover them. Moreover, since M only changes value at these points, these points are the only points where information about σ can be learned. Thus, the gradients at these points are the only gradients which could possibly be learned. We note that the running time of our algorithms will be polynomial in log(1/η), and thus we can even allow the gradient size η to be exponentially small.

heading: 4.1 THE APPROXIMATE SPAN RECOVERY ALGORITHM
text: We now formally describe and analyze our span recovery algorithm for networks with differentiable activation functions and 0/1 thresholding. Let κi be the condition number of the i-th weight matrix Wi, and let δ > 0 be a failure probability, and let > 0 be a precision parameter which will affect the how well the subspace we output will approximate Span(A). Now fix N = poly(n, k, 1γ , ∑d i=1 log(Li), ∑d i=1 log(κi), log( 1 η ), log( 1 ), log( 1 δ )). The running time and query complexity of our algorithm will be polynomial in N . Our algorithm for approximate span recovery is given formally in Algorithm 2. Proposition 4.1. Let V ⊂ Rn be a subspace of dimension k′ < k, and fix any 0 > 0. Then we can find a vector x with 0 ≤ σV (x) − 1 ≤ 2 0 in expected O(1/γ + N log(1/ 0)) time. Moreover, with probability γ/2 we have that ∇xσV (x) > η/4 and the tighter bound of 0η2−N ≤ σV (x)− 1 ≤ 2 0.
We will apply the above proposition as input to the following Lemma 4.2, which is the main technical result of this section. Our approach involves first taking the point x from Proposition 4.1 such that σV (x) is close but bounded away from the boundary, and generating n perturbations at this point MV (x+ ui) for carefully chosen ui. While we do not know the value of σV (x+ui), we can tell for a given scaling c > 0 if σV (x+ cui) has crossed the boundary, since we will then have MV (x+ cui) = 0. Thus, we can estimate the directional derivative ∇uiσ(x) by finding a value ci via a binary search such that σV (x + ciui) is exponentially closer to the boundary than σV (x). In order for our estimate to be accurate, we must carefully upper and lower bound the gradients and Hessian of σv near x, and demonstrate that the linear approximation of σv at x is still accurate at the point x + ciui where the boundary is crossed. Since each
Algorithm 2: Span Recovery With Adaptive Gradients 1 V ← ∅, N = poly ( n, k, 1γ , ∑d i=1 log(Li), ∑d i=1 log(κi), log( 1 η ), log( 1 ), log( 1 δ ) ) , 0 ← 2−poly(N) 2 for i = 0, . . . , k do 3 Generate g ∼ N (0, In) until M((In −PV )g) = 1 4 Find a scaling α > 0 via binary search on values τ(σV (αg)) such that x = αg satisfies
0η2 −N ≤ σV (x)− 1 ≤ 2 0. . Proposition 4.1 5 Generate g1, . . . , gn ∼ N (0, In), and set ui = ( gi2 −N − x/‖x‖2 ) . 6 For each i ∈ [n], binary search over values c to find ci such that
1− β ≤ σV (x+ ciui) ≤ 1
where β = 2−N 2
20. 7 If any ci satisfies |ci| ≥ (10 · 2−N 0/η), restart from line 5 (regenerate the Gaussian g). 8 Otherwise, define B ∈ Rn×n via B∗,i = ui, where B∗,i is the i-th column of B. Define b ∈ Rn by bi = 1/ci.
9 Let y∗ be the solution to: min y∈Rn ‖yTB− bT ‖22
Set vi = y∗ and V ← Span(V, vi). 10 end 11 return V
value of 1/ci is precisely proportional to∇uiσ(x) = 〈∇σ(x), ui〉, we can then set up a linear system to approximately solve for the gradient∇σ(x) (lines 8 and 9 of Algorithm 2). Lemma 4.2. Fix any , δ > 0, and let N be defined as above. Then given any subspace V ⊂ Rn with dimension dim(V ) < k, and given x ∈ Rn, such that 0η2−N ≤ σV (x) − 1 ≤ 2 0 where 0 = Θ(2−N C
/ ) for a sufficiently large constant C = O(1), and such that ∇xσV (x) > η/2, then with probability 1 − 2−N/n 2
, we can find a vector v ∈ Rn in expected poly(N) time, such that ‖PSpan(A)v‖2 ≥ (1− )‖v‖2, and such that ‖PV v‖2 ≤ ‖v‖2.
Theorem 4.3. Suppose the network M(x) = τ(σ(Ax)) satisfies the conditions described at the beginning of this section. Then Algorithm 2 runs in poly(N) time, making at most poly(N) queries to M(x), where N = poly(n, k, 1γ , ∑d i=1 log(Li), ∑d i=1 log(κi), log( 1 η ), log( 1 ), log( 1 δ )), and returns with probability 1 − δ a subspace V ⊂ Rn of dimension k such that for any x ∈ V , we have ‖PSpan(A)x‖2 ≥ (1− )‖x‖2.

heading: 5 EXPERIMENTS
text: When applying span recovery for a given network, we first calculate the gradients analytically via auto-differentiation at a fixed number of sample points distributed according to a standard Gaussian. Our networks are feedforward, fully-connected with ReLU units; therefore, as mentioned above, using analytic gradients is as precise as using finite differences due to piecewise linearity. Then, we compute the rank of the resulting gradient matrix, where the rank is defined to be the number of singular values that are above 1e-5 of the maximum singular value. In our experiments, we attempt to recover the full span of a 784-by-80 matrix with decreasing layer sizes for varying sample complexity, as specified in the figures. For the MNIST dataset, we use a size 10 vector output and train according to the softmax cross entropy loss, but we only calculate the gradient with respect to the first output node.
Our recovery algorithms are GradientsRandom (Algorithm 1), GradientsRandomAda (Algorithm 2), and GradientsMNIST. GradientsRandom is a direct application of our first span recovery algorithm and calculates gradients via perturbations at random points for a random network. GradientsRandomAda uses our adaptive span recovery algorithm for a random network. Finally, GradientsMNIST is an application of GradientsRandom on a network with weights trained on MNIST data. In general, we note that the experimental outcomes are very similar among all three scenarios.
For networks with very small widths and multiple layers, we see that span recovery deteriorates as depth increases, supporting our theory (see Figure 1). This holds both in the case when the networks are randomly initialized with Gaussian weights or trained on a real dataset (MNIST) and whether we use adaptive or non-adaptive recovery algorithms. However, we note that these small networks have unrealistically small widths (less than 10) and when trained on MNIST, these networks fail to achieve high accuracy, all falling below 80 percent. The small width case is therefore only used to support, with empirical evidence, why our theory cannot possibly guarantee full span recovery under every network architecture.
For more realistic networks with moderate or high widths, however, full span recovery seems easy and implies a real possibility for attack (see Figure 2). Although we tried a variety of widths and depths, the results are robust to reasonable settings of layer sizes and depths. Therefore, we only present experimental results with sub-networks of a
network with layer sizes [784, 80, 40, 30, 20, 10]. Note that full span recovery of the first-layer weight matrix with rank 80 is achieved almost immediately in all cases, with less than 100 samples.
On the real dataset MNIST, we demonstrate the utility of span recovery algorithms as an attack to fool neural networks to misclassify noisy inputs (see Figure 3). We train a ReLU network (to around 95 percent accuracy) and recover its span by computing the span of the resulting gradient matrix. Then, we recover the null space of the matrix and generate random Gaussian noise projected onto the null space. We see that our attack successfully converts images into noisy versions without changing the output of the network, implying that allowing a full (or even partial) span recovery on a classification network can lead to various adversarial attacks despite not knowing the exact weights of the network.

heading: ACKNOWLEDGMENTS
text: The authors Rajesh Jayaram and David Woodruff would like to thank the partial support by the National Science Foundation under Grant No. CCF-1815840.

heading: A MISSING PROOFS FROM SECTION 3
text: We first restate the results which have had their proofs omitted, and include their proofs subsequently.
Lemma 3.1 If A ∈ Rk×n has orthogonal rows, or if n > Ω(k3) and A has entries that are drawn i.i.d. from some sub-Gaussian distribution D with expectation 0, unit variance, and constant sub-gaussian norm ‖D‖ψ2 = supp≥1 p −1/2 (EX∼D|X|p)1/p then with probability at least 1− e−k 2 , A satisfies Assumption 1 with γ ≥ 1/2.
Moreover, if the weight matrices w,W1,W2, ...,Wd with Wi ∈ Rki×ki+1 have entries that are drawn independently (and possibly non-identically) from continuous symmetric distributions, and if ki ≥ log( 16dδγ ) for each i ∈ [d], then Assumption 2 holds with probability 1− δ.
Proof. By Theorem 5.58 of Vershynin (2010), if the entries A are drawn i.i.d. from some sub-Gaussian isotropic distribution D over Rn such that ‖Aj‖2 = √ n almost surely, then √ n − C √ k − t ≤ σmin(A) ≤ σmax(A) ≤√
n + C √ k + t with probability at least 1 − 2e−ct2 , for some constants c, C > 0 depending only on ‖D‖ψ2 . Since the entries are i.i.d. with variance 1, it follows that the rows of A are isotropic. Moreover, we can always condition on the rows having norm exactly √ n, and pulling out a positive diagonal scaling through the first Relu of M(x), and absorbing this scaling into Wd. It follows that the conditions of the theorem hold, and we have √ n − C √ k ≤ σmin(A) ≤ σmax(A) ≤ √ n + C √ k with probability at least 1 − e−k2 for a suitably large re scaling of the constant C. Setting n > Ω(k3), it follows that κ(A) < (1 + 1/(100k)), which holds immediately if A has orthogonal rows.
Now observe that Ag is distributed as a multi-variate Gaussian with co-variance ATA, and is therefore given by the probability density function (pdf)
p′(x) = 1
(2π)k/2det(ATA)1/2 exp
( −1
2 xTATAx ) Let p(x) = 1
(2π)k/2 exp ( − 12x Tx )
be the pdf of an identity covariance GaussianN (0, Ik). We lower bound p′(x)/p(x) for x with ‖x‖22 ≤ 16k. In this case, we have
p′(x) p(x) =
1
det(ATA)1/2 exp
( −1
2
( xTATAx− ‖x‖22 )) ≥κ(A)−k exp ( −1
2
( ‖x‖22(1 + 1/(100k))− ‖x‖22 )) ≥κ(A)−k exp ( −1
2 |x‖22/(100k) ) ≥(1 + 1/(100k))−k exp ( −1
2 (16/100) ) ≥(1 + 1/(100k))−k exp ( −1
2 (16/100) ) ≥1/2
(1)
Thus for any sign pattern S, Pr[sign(Ag) = S : ‖Ag‖22 ≤ k] ≥ 12Pr[sign(g) = S : ‖g‖ 2 2 ≤ k]. Now Pr[sign(g) = S] = 2−k, and spherical symmetry of Gaussians, Pr[sign(g) = S : ‖g‖22 ≤ k] = 2−k, and thus Pr[sign(Ag) = S : ‖Ag‖22 ≤ k] ≥ 2−k−1. Now ‖Ag‖22 ≤ 2‖g‖22 which is distributed as a χ2 random variable with k-degrees of freedom. By standard concentration results for χ2 distributions (Lemma 1 Laurent & Massart (2000)), we have Pr[‖g‖22 ≥ 8k] ≤ e−2k, so Pr[‖Ag‖22 ≥ 16k] ≤ Pr[‖g‖22 ≥ 8k] ≤ e−2k. By a union bound, Pr[sign(Ag) = S : ‖Ag‖22 ≤ k] ≥ 2−k−1 − e−2k ≥ 142
−k, which completes the proof of the first claim with γ = 1/4.
For the second claim, by an inductive argument, the entries in the rows i ∈ Sj of the product WSj (∏d i=j+1(Wi)Si ) are drawn from a continuous distribution. Thus each column of WSj (∏d i=j+1(Wi)Si ) is non-zero with probability
1. It follows that 〈w, (∏d
i=1(Wi)Si ) ∗,j 〉 is the inner product of a non-zero vector with a vector w with continuous,
independent entries, and is thus non-zero with probability 1. By a union bound over all possible non-empty sets Sj , the desired result follows.
We now show that the second part of Assumption 2 holds. To do so, first let g ∼ N (0, In). We demonstrate that PrW1,W2,...,Wd,g [M(x) = 0] ≤ 1−γδ/100. Here the entries of the Wi’s are drawn independently but not necessarily identically from a continuous symmetric distribution. To see this, note that we can condition on the value of g, and condition at each step on the non-zero value of yi = φ(Wi+1φ(Wi+2φ(. . . φ(Ag) . . . ). Then, over the randomness of Wi, note that the inner product of a row of Wi and yi is strictly positive with probability at least 1/2, and so each coordinate of Wiyi is strictly positive independently with probability ≥ 1/2. It follows that φ(Wiyi) is non-zero with probability at least 1− 2−ki . Thus
PrW1,W2,...,Wd,g [M(g) 6= 0] ≥ d∏ i=1 (1− 2−ki)
≥ (
1− δγ 16d )d ≥ 1− δγ/8
(2)
where the second inequality is by assumption. It follows by our first part that
PrW1,W2,...,Wd,g [M(g) = 0] ≤ δγ/8
So by Markov’s inequality, PrW1,W2,...,Wd [Prg [M(g) = 0] ≥ γ/8] ≤ δ
Thus with probability 1− δ over the choice of W1, . . . ,Wd, we have Prg [M(g) = 0] ≤ γ/8 as desired.
Lemma 3.2 There is an algorithm which, given g ∼ N (0, Ik), with probability 1−exp(−nc) for any constant c > 1 (over the randomness in g), computes∇M(g) ∈ Rn with O(n) queries to the network, and in poly(n) running time.
Proof. Let Mi(x) = φ(Wiφ(Wi+1φ(. . . φ(Ax)) . . . ) where φ is the ReLU. If ∇M(g) exists, there is an > 0 such that M is differentiable on B (g). We show that with good probability, if g ∼ N (0, In) (or in fact, almost any continuous distribution), then M(g) is continuous in the ball B (g) = {x ∈ Rn | ‖x− g‖2 < } for some which we will now compute.
First, we can condition on the event that ‖g‖22 ≤ (nd)10c, which occurs with probability at least 1− exp(−(nd)5c) by concentration results for χ2 distributions Laurent & Massart (2000). Now, fix any sign pattern Si ⊆ [ki] for the i-th layer Mi(x) = φ(Wi(φ(. . . (φ(Ax) . . . ), and let S = (S1, S2, . . . , Sd+1). We note that we can enforce the constraint that for an input x ∈ Rn, the sign pattern of Mi(x) is precisely Si. To see this, note that after conditioning on a sign pattern for each layer, the entire network becomes linear. Thus each constraint that 〈(Wi)j,∗ ,Mi+1(x)〉 ≥ 0 or 〈(Wi)j,∗ ,Mi+1(x)〉 ≤ 0 can be enforced as a linear combination of the coordinates of x. Now fix any layer i ∈ [d + 1] and neuron j ∈ [ki], WLOG j ∈ Si. We now add the additional constraint that 〈(Wi)j,∗ ,Mi+1(x)〉 ≤ η, where η = exp(−poly(nd)) is a value we will later choose. Thus, we obtain a linear program with k + ∑d i=1 ki constraints and n variables. The feasible polytope P represents the set of input points which satisfy the activation patterns S and are η-close to the discontinuity given by the j-th neuron in the i-th layer.
We can now introduce the following non-linear constraint on the input that ‖x‖2 ≤ (nd)10c. Let B = B(nd)10c(~0) be the feasible region of this last constraint, and let P∗ = P ∩B. We now bound the Lesbegue measure (volume) V (P∗) of the region P . First note that V (P∗) ≤ V (P ′), where V (P ′) is the region defined by the set of points which satisfy:
〈y, x〉 ≥ 0 〈y, x〉 ≤ η
‖x‖22 ≤ (nd)10c (3)
where each coordinate of the vector y ∈ Rn is a linear combination of products of the weight matrices W`, ` ≥ i. One can see that the first two constraints for P ′ are also constraints for P∗, and the last constraint is precisely B, thus P∗ ⊂ P ′ which completes the claim of the measure of the latter being larger. Now we can rotate P ′ by the rotation
which sends y → ‖y‖2 ·e1 ∈ Rn without changing the volume of the feasible region. The resulting region is contained in the region P ′′ given by
0 ≤ x1 ≤ η‖x‖2∞ ≤ (nd)10c (4)
Finally, note that P ′′ ⊂ Rn is a Eucledian box with n− 1 side lengths equal to (nd)10c and one side length of ‖y‖2η, and thus V (P ′′) ≤ ‖y‖2η(nd)10nc. Now note we can assume that the entries of the weight matrices A,W1, . . . ,Wd are specified in polynomially many (in n) bits, as if this were not the case the output M(x) of the network would not have polynomial bit complexity, and could not even be read in poly(n) time. Equivalently, we can assume that our running time is allowed to be polynomial in the number of bits in any value of M(x), since this is the size of the input to the problem. Given this, since the coordinates of y were linear combinations of products of the coordinates of the weight matrices, and note that each of which is at most 2n C
for some constant C (since the matrices have polynomial bit-complexity), we have that P∗ ≤ η2nC (nd)10nc as needed.
Now the pdf of a multi-variate Gaussian is upper bounded by 1, so Prg∼N (0,In) [g ∈ P∗] ≤ V (P∗) ≤ η2n C (nd)10nc. It follows that the probability that a multivariate Gaussian g ∼ N (0, In) satisfies the sign pattern S and is η close to the boundary for the j-th neuron in the i-th layer. Now since there are at most 2k · ∏d i=1 2
ki ≤ 2nd possible combinations of sign patterns S, it follows that the the probability that a multivariate Gaussian g ∈ N (0, In) is η close to the boundary for the j-th neuron in the i-th layer is at most η2n C
(nd)10nc2nd. Union bounding over each of the ki neurons in layer i, and then each of the d layers, it follows that g ∈ N (0, In) is η close to the boundary for any discontinuity in M(x) is at most η2n C
(nd)10nc+12nd. Setting η ≤ 2(nd)20c+1)2−nC , it follows that with probability at least 1− exp(−(nd)c), the network evaluated at g ∈ N (0, In) is at least η close from all boundaries (note that C is known to us by assumption).
Now we must show that perturbing the point g by any vector with norm at most results in a new point g′ which still has not hit one of the boundaries. Note that M(g) is linear in an open ball around g, so the change that can occur in any intermediate neuron after perturbing g by some v ∈ Rn is at most ‖A‖2 ∏d i=1 ‖Wi‖2, where ‖ · ‖2 is the spectral norm. Now since each entry in the weight matrix can be specified in polynomially many bits, the Frobenius norm of each matrix (and therefore the spectral norm), is bounded by n22n C for some constant C. Thus
‖A‖2 d∏ i=1 ‖Wi‖2 ≤ (n22n C )d+1 = β
and setting = η/β, it follows that M(x) is differentiable in the ball B (x) as needed. We now generate u1, u2, . . . , un ∼ N (0, In), which are linearly independent almost surely. We set vi = ui2‖ui‖2 . SinceM(g) is a ReLU network which is differentiable on B (g), it follows thatM(g) is a linear function on B (g), and moreover vi ∈ B (g) for each i ∈ [n]. Thus for any c < 1 we have M(g)−M(g+cvi)c = ∇viM(x), thus we can compute ∇viM(x) for each i ∈ [n]. Finally, since the directional derivative is given by ∇viM(x) = 〈∇M(x), vi/‖vi‖2〉, and since v1, . . . , vn are linearly independent, we can set up a linear system to solve for ∇M(x) exactly in polynomial time, which completes the proof.
Lemma 3.3 Let V ∈ Rr×k be a fixed at most rank ck matrix for c ≤ 1/2. Then the number of sign patterns S ⊂ [k] with at most k/2 non-zeros spanned by the rows of V is at most 2
k √ k . In other words, the set S(V) = {sign(w) | w ∈ span(V), nnz(w) ≤ k2} has size at most 2k√ k .
Proof. Any vector w in the span of the rows of V can be expressed as a linear combination of at most ck rows of V. So create a variable xi for each coefficient i ∈ [ck] in this linear combination, and let fj(x) be the linear function of the x′is which gives the value of the j-th coordinate of w. Then f(x) = (f1(x), . . . , fk(x)) is a k-tuple of polynomials, each in ck-variables, where each polynomial has degree 1. By Theorem 4.1 of Hall et al. (2010), it follows that the number of sign patterns which contain at most k/2 non-zero entries is at most ( ck+k/2 ck ) . Setting c ≤ 1/2, this is at
most ( k k/2 ) ≤ 2 k √ k .
Theorem 3.4 Suppose the network M(x) = wTφ(W1φ(W2φ(. . .Wdφ(Ax)) . . . ), where φ is the ReLU, satisfies Assumption 1 and 2. Then the algorithm given in Figure 1 makes O(kn log(k/δ)/γ) queries to M(x) and returns in poly(n, k, 1/γ, log(1/δ))-time a subspace V ⊆ span(A) of dimension at least k2 with probability 1− δ.
Proof. First note that by Lemma 3.2, we can efficiently compute each gradient ∇M(gi) using O(n) queries to the network. After querying for the gradient zi = ∇M(gi) for r′ ≤ r independent Gaussian vectors gi ∈ Rn, we obtain the vector of gradients Z = V · A ∈ Rr′×k. Now suppose that V had rank ck for some c ≤ 1/2. Now consider the gradient ∇(gr′+1), which can be written as zr′+1 = wT (∏d i=1(Wi)Sr′+1i ) Diag(sign(Agr ′+1))A. Thus we can write zr ′+1 = Vr′+1A, where Vr′+1 ∈ Rk is a row vector which will be appended to the matrix V to form a new Z = VA ∈ Rr′+1×k after the (r′ + 1)-th gradient is computed. Specifically, for any j ∈ [r], we have: Vj = w T (∏d
i=1(Wi)Sji
) Diag(sign(Agj)).
Let Er′+1 be the event that M(gr ′+1) 6= 0. It follows necessarily that, conditioned on Er′+1, we have that sign(Vr′+1) = Sr ′+1 0 . The reason is as follows: if M(g r′+1) 6= 0, then we could not have Sr ′+1 j = ∅ for any j ∈ {0, 1, 2, . . . , d}, since this would result in M(gr′+1) = 0. It follows that the conditions for Assumption 2 to apply hold, and we have that wT (∏d i=1(Wi)Sr′+1i ) ∈ Rk is entry-wise non-zero. Given this, it follows that the sign pattern of
Vr′+1 = w T ( d∏ i=1 (Wi)Sr′+1i ) Diag(sign(Agj))
will be precisely Sr ′+1
0 , as needed.
Now by Lemma 3.3, the number of sign patterns of k-dimensional vectors with at most k/2 non-zero entries which are contained in the span of V is at most 2
k √ k . Let S be the set of sign patterns with at most k/2 non-zeros and such that for every S ∈ S, S is not a sign pattern which can be realized in the row span of V. It follows that |S| ≥ 122 k − 2 k √ k ≥ 142 k, so by Assumption 1, we have that Pr [ sign(Agr ′+1) ∈ S ] ≥ γ/4. By a union bound, we
have Pr [ sign(Agr ′+1) ∈ S, Er′+1 ] ≥ γ/8
Conditioned on sign(Agr ′+1) ∈ S) and Er′+1 simultaneously, it follows that adding the row vector Vr′+1 to the matrix V will increase its rank by 1. Thus after O(log(k/δ)/γ) repetitions, the rank of V will be increased by at least 1 with probability 1 − δ/k. By a union bound, after after r = O(k log(k/δ)/γ), V ∈ Rr×k will have rank at least k/2 with probability at least 1 − δ, which implies that the same will hold for Z since rank(Z) ≥ rank(V), which is the desired result.

heading: B MISSING PROOFS FROM SECTION 4
text: Proposition 4.1 Let V ⊂ Rn be a subspace of dimension k′ < k, and fix any 0 > 0. Then we can find a vector x with
0 ≤ σV (x)− 1 ≤ 2 0 in expectedO(1/γ+N log(1/ 0)) time. Moreover, with probability γ/2 we have that∇xσV (x) > η/4 and the tighter bound of
0η2 −N ≤ σV (x)− 1 ≤ 2 0
Proof. We begin by generating Gaussians g1, . . . and computing MV (gi). By Property 3 of the network assumptions, we need only repeat the process 1/γ times until we obtain an input y = (In − PV )gi with M(y) = MV (gi) = 1. Since all activation functions satisfy φi(0) = 0, we know that σV (0 · gi) = 0, and σV (gi) = 1. Since σ is continuous, it follows that ψgi(c) = σV (c · gi) is a continuous function ψ : R 7→ R. By the intermediate value theorem, there exists a c∗ ≤ 1 such that σV (c∗gi) = 1. We argue we can find a c with |c− c∗| ≤ 02−N in time O(N log(1/ 0)). To find c, we can perform a binary search. We first try c0 = 1/2, and if ψgi(c0) = 0, we recurse into [1/2, 1], otherwise if ψgi(c0) = 1 we recurse into [0, 1/2]. Thus, we always recurse into an interval where ψgi switches values. It follows that we can find a c with |c − c∗| ≤ 0‖gi‖22−N in time O(N log(‖gi‖2/ 0)) for some c∗ with σV (c
∗gi) = 1. Now observe that it suffices to binary search a total of O(N log(‖gi‖2/ 0)) times, since 2N is an upper bound on the Lipschitz constant of ψx, which gives 0 ≤ σV (cgi)− 1 ≤ 0. Now the expected running time to do this is O(1/γ+N log(‖gi‖2/ 0)), but since ‖gi‖2 has Gaussian tails, the expectation of the maximum value of ‖gi‖2 over
1/γ repetitions is O(log(1/γ) √ n), and thus the expected running time reduces to the stated bound, which completes the first claim of the Proposition.
For the second claim, note that ∇giσV (c∗gi) > 0 by construction of the binary search, and since σV (c∗gi) > 0 = 1, by Property 4 with probability γ we have that ∇giσV (gi) > η. Now with probability 1 − γ/2, we have that ‖gi‖22 ≤ O(n log(1/γ)) (see Lemma 1 Laurent & Massart (2000)), so by a union bound both of these occur with probability γ/2. Now since ‖(c∗ − c)x‖2 ≤ 02−N (after rescaling N by a factor of log(‖gi‖2) = O(log(n))), and since 2N is also an upper bound on the spectral norm of the Hessian of σ by construction, it follows that∇giσV (cgi) > η/2.
Now we set x← cgi+c 02−Ngi/(‖cgi‖2). First note that this increases σV (cx)−1 by at most 0, so σ(cx)−1 ≤ 2 0, so this does not affect the first claim of the Proposition. But in addition, note that conditioned on the event in the prior paragraph, we now have that σV (x) > 1 + η 02−N . The above facts can be seen by the fact that 2N is polynomially larger than the spectral norm of the Hessian of σ, thus perturbing x by 02−N additive in the direction of x will result in a positive change of at least 12 (η/4)( 02
−N ) in σ. Moreover, by applying a similar argument as in the last paragraph, we will have ∇xσV (cx) > η/4 still after this update to x. Lemma 4.2 Fix any , δ > 0, and letN = poly(n, k, 1γ , ∑d i=1 log(Li), ∑d i=1 log(κi), log( 1 η ), log( 1 ), log( 1 δ )). Then given any subspace V ⊂ Rn with dimension dim(V ) < k, and given x ∈ Rn, such that 0η2−N ≤ σV (x)− 1 ≤ 2 0 where 0 = Θ(2−N C
/ ) for a sufficiently large constant C = O(1), and ∇xσV (x) > η/2, then with probability 1− 2−N/n2 , we can find a vector v ∈ Rn in expected poly(N) time, such that
‖PSpan(A)v‖2 ≥ (1− )‖v‖2
and such that ‖PV v‖2 ≤ ‖v‖2.
Proof. We generate g1, g2, . . . , gn ∼ N (0, In), and set ui = gi2−N − x/‖x‖2. We first condition on the event that ‖gi‖2 ≤ N for all i, which occurs with probability 1 − 2−N/n 2
. Note ∇ui [σV (x)] = wTA(In − PV )ui, where wT = ∇ [σ(A(In −PV )x)]T ∈ Rk, which does not depend on ui. Thus wTA(In −PV ) is a vector in the row span of A(In −PV ). We can write
MV (x+ cui) = τ ( (σV (x) + cw TA(In −PV )ui + Ξ(cui) )
where Ξ(cui) = O(‖c(In − PV )ui‖222N ) = O(‖cui‖222N ) is the error term for the linear approximation. Note that the factor of N comes from the fact that the spectral norm of the Hessian of σ : Rn → R can be bounded by∏ i κiLi ≤ 2N . Fix some β > 0. We can now binary search again, as in Proposition 4.1, with O(log(N/β)) iterations over c, querying values MV (x+ cui) to find a value c = ci > 0 such that
1− β ≤MV (x+ cui) ≤ 1
so 1− β ≤ ( σV (x) + ciw TA(In −PV )ui + Ξ(ciui) ) ≤ 1
We first claim that the ci which achieves this value satisfies ‖ciui‖2 ≤ (10 · 2−N 0/η). To see this, first note that by Proposition 4.1, we have ∇xσV (x) > η/4 with probability γ. We will condition on this occurring, and if it fails to occur we argue that we can detect this and regenerate x. Now conditioned on the above, we first claim that ∇uiσV (x) ≥ η/8, which follows from the fact that we can bound the angle between the unit vectors in the directions of ui and x by
cos (angle(ui, x)) = 〈 ui ‖u‖2 , x ‖x‖2 〉 ≥ (1− n/2−N ) > (1− η/2−N/2)
along with the fact that we have∇xσV (x) > η/4. Since |σV (x)− 1| < 2 0 < 2−N C , and since 2N is an upper bound on the spectral norm of the Hessian of σ, we have that ∇uiσV (x + cui) > η/8 + 2−N > η/10 for all c < 2−2N . In other words, if H is the hessian of σ, then perturbing x by a point with norm O(c) ≤ 2−2N can change the value of the gradient by a vector of norm at most 22N‖H‖2 ≤ 2−N , where ‖H‖2 is the spectral norm of the Hessian. It follows that setting c = (10 · 2−N 0/η) is sufficient for σV (x+ cui) < 1, which completes the above claim. Now observe that if after binary searching, the property that c ≤ (10 · 2−N 0/η) does not hold, then this implies that we did not have ∇σ(x) > η/4 to begin with, so we can throw away this x and repeat until this condition does hold. By Assumption 4, we must only repeat O(1/γ) times in expectation in order for the assumption to hold.
Next, also note that we can bound ci ≥ 0η2−N/N , since 2N again is an upper bound on the norm of the gradient of σ and we know that σV (x) − 1 > 0η2−N . Altogether, we now have that |Ξ(ciui)| ≤ c2i 2N ≤ (10 · 2−N 0/η)22N .
We can repeat this binary search to find ci for n different perturbations u1, . . . , un, and obtain the resulting c1, . . . , cn, such that for each i ∈ [n] we have
1− σV (x)− Ξ(ciui)− βi = ciwTA(In −PV )ui
where βi is the error obtained from the binary seach on ci, and therefore satisfies |βi| ≤ 2−N 2 20 taking poly(N) iterations in the search. Now we know ci and ui, so we can set up a linear system
min y ‖yTB− bT ‖22
for an unknown y ∈ Rk where the i-th column of B is given by Bi = ui ∈ Rn, and bi = 1/ci for each i ∈ [n]. First note that the set {u1, u2, . . . , un} is linearly independent with probability 1, since {u1 − x, u2 − x, . . . , un − x} is just a set of Gaussian vectors. Thus B has rank n, and the above system has a unique solution y∗.
Next, observe that setting ŷ = w TA(In−PV ) (1−σV (x)) , we obtain that for each i ∈ [n]:
(ŷB)i = 1 ci − 1 ci Ξ(ciui)− βi 1− σV (x)
Thus
‖ŷTB − bT ‖2 ≤ ( n∑ i=1 ( 1 ci · Ξ(ciui)− βi 1− σV (x) )2)1/2
≤ n( 1 ci · (10 · 2 −N 0/η) 22N − 2−N2 20 1− σV (x) )21/2
≤ (2−N )O(1) 2 0
ci(1− σ(x))
(5)
Now setting y∗ such that (y∗)TB = bT , we have that cost of the optimal solution is 0, and ‖ŷB − bT ‖2 ≤ (2−N )O(1)
20 ci(1−σ(x)) , so ‖(y ∗ − ŷ)B‖2 ≤ (2−N )O(1) 2 0 ci(1−σ(x)) . By definition of the minimum singular value of
B, it follows that ‖y∗ − ŷ‖2 ≤ 1σmin(B) (2 −N )O(1) 20 ci(1−σ(x)) . Now using the fact that B = 2 −NG + X where G is a Gaussian matrix and X is the matrix with each row equal to x, we can apply Theorem 1.6 of Cook et al. (2018), we have σmin(B) ≥ 1/(2−N )O(1), So ‖y∗−ŷ‖2 ≤ (2−N )O(1) 2 0
ci(1−σ(x)) . Thus we have ‖y ∗‖2 ≤ ‖ŷ‖2+(2−N )O(1)
2 0
ci(1−σ(x)) , and moreover note that ‖ŷ‖2 ≥ ‖∇σ(x)‖2 11−σ(x) > η 8(1−σ(x)) . So we have that
‖y∗ − ŷ‖2 ‖ŷ‖2 ≤ (2−N )O(1) 20 ci(1−σ(x)) ‖y∗‖2
≤ (2−N )O(1) 2 0
ci
≤ (2−N )O(1) 0 ≤ (2−N )O(1) 0 ≤ 2−N/2
(6)
Where the last inequality holds by taking C larger than some constant in the definition of 0. Thus ‖y∗ − ŷ‖2 ≤ 2−N/2‖ŷ‖2, thus ‖y∗ − ŷ‖2 ≤ 2 · 2−N/2‖y∗‖2 ≤ ‖y∗‖2 after scaling N up by a factor of log2(1/ ). Thus by setting v = y∗, and observing that ŷ is in the span of A, we ensure ‖PSpan(A)v‖2 ≥ (1 − )‖v‖2 as desired. For the final claim, note that if we had y∗ = ŷ exactly, we would have ‖PV y∗‖2 = 0, since ŷ is orthogonal to the subspace V . It follows that since ‖y∗ − ŷ‖22 ≤ 2‖y∗‖22, we have ‖(In −PV )y∗‖22 ≥ (1− 2)‖y∗‖22, so by the Pythagorean theorem, we have ‖PV y∗‖22 = ‖y∗‖22 − ‖(In −PV )y∗‖22 ≤ 2‖y∗‖22 as desired.
Theorem 4.3 Suppose the network M(x) = τ(σ(Ax)) satisfies the conditions described at the beginning of this section. Then Algorithm 2 runs in poly(N) time, making at most poly(N) queries to M(x), where N =
poly(n, k, 1γ , ∑d i=1 log(Li), ∑d i=1 log(κi), log( 1 η ), log( 1 ), log( 1 δ )), and returns with probability 1 − δ a subspace V ⊂ Rn of dimension k such that for any x ∈ V , we have ‖PSpan(A)x‖2 ≥ (1− )‖x‖2
Proof. We iteratively apply Lemma 4.2, each time appending the output v ∈ Rn of the proposition to the subspace V ⊂ Rn constructed so far. WLOG we can assume v is a unit vector by scaling it. Note that we have the property at any given point in time k′ < k that V = Span(v1, . . . , vk′) where each vi satisfies that ‖PSpan{v1,...,vi−1}vi‖2 ≤ . Note that the latter fact implies that v1, . . . vk′ are linearly independent. Thus at the end, we recover a rank k subspace V = Span(v1, . . . , vk), with the property that ‖PSpan(A)vi‖22 ≥ (1− )‖vi‖22 for each i ∈ [k].
Now let V ∈ Rn×n be the matrix with i-th column equal to vi. Fix any unit vector x = Va ∈ V , where a ∈ Rn is uniquely determined by x. Let V = V+ +V− where V+ = PSpan(A)V and V− = V−V+ Then x = V+a+V−a, and
‖(In −PSpan(A))x‖2 ≤ ‖(In −PSpan(A))V+a‖2 + ‖(In −PSpan(A))V−a‖2 ≤ ‖(In −PSpan(A))V−a‖2 ≤ ‖V−‖2‖a‖2
(7)
First note that by the construction of the vi’s, each column of V− has norm O( ), thus ‖V−‖2 ≤ O( √ n ). Moreover, since ‖x‖2 = 1, it follows that ‖a‖2 ≤ 1σmin(V) , which we now bound. Since ‖PSpan{v1,...,vi−1}vi‖2 ≤ for each i ∈ [k], we have
‖Va‖2 ≥ ‖ n∑ i=1 viai‖2
≥ ‖ n∑ i=1 (In −PSpan{v1,...,vi−1})viai + n∑ i=1 PSpan{v1,...,vi−1}viai‖2
≥ ‖ n∑ i=1 (In −PSpan{v1,...,vi−1})viai‖2 − ‖ n∑ i=1 PSpan{v1,...,vi−1}viai‖2
≥ ‖ n∑ i=1 (In −PSpan{v1,...,vi−1})viai‖2 −O( )‖a‖2
= ( ‖
n∑ i=1 (In −PSpan{v1,...,vi−1})viai‖ 2 2 −O( )‖a‖2‖ n∑ i=1 (In −PSpan{v1,...,vi−1})viai‖2 +O( 2)‖a‖22
)1/2
= ( ‖
n∑ i=1 (In −PSpan{v1,...,vi−1})viai‖ 2 2 −O( )‖a‖22
)1/2
= ( n∑ i=1 ‖(In −PSpan{v1,...,vi−1})viai‖ 2 2 −O( )‖a‖22 )1/2
≥ ( n∑ i=1 (1−O( 2))|ai|2 −O( )‖a‖22 )1/2 ≥ ( ‖a‖22 −O( )‖a‖22
)1/2 ≥ (1−O( ))‖a‖2
(8)
Thus σmin(V) ≥ (1 − O( )), so we have ‖(In − PSpan(A))x‖2 ≤ ‖V−‖2 1σmin(V) ≤ 2 √ n . By the Pythagorean
theorem: ‖PSpan(A))x‖22 = 1− ‖(In −PSpan(A))x‖22 ≥ 1−O(n 2). Thus we can scale by a factor of Θ(1/ √ n) in the call to Lemma 4.2, which gives the desired result of ‖PSpan(A)x‖2 ≥ 1− .

Review summaries:
The paper studies how to recover the span of a NN from a limited number of queries. The problem belongs to the general question of how to reconstruct functions from black-box interaction and it may find application in obfuscation attacks where very large perturbations of the input do not affect the output. The main contribution of the paper is on the theoretical analysis of a simple non-adaptive and a more sophisticated adaptive algorithm. The main finding is that under mild conditions on the structure of the NN partial recovery is possible. The empirical validation show that in practice, it is often the case the full span recover is actually possible, as the structure and weights of common NN are "friendly" enough.

The paper considers the problem of recovering the span of the latent variables of a neural network with various activation functions. More precisely, if we write a neural network as M(x) = f(Ax), for some neural network f, and a matrix A: R^{k x d} -> R^{k}, k << d, we wish to recover the row span of A. The authors consider ReLU activations -- in which case they can recover at least "half" of the row span with ~kd queries to M(x), and smooth activations -- in which case they can approximately recover the row span in poly(k,d) queries. The authors also consider (empirically) applications of these algorithms to "input obfuscation": namely generating samples which are effectively noise, but the network classifies them as "structure" (e.g. digits on MNIST).