Beyond Scalar Probes: Exploiting Vector-Valued Outputs in ReLU Networks For Signature Extraction
arXiv:2609.35487v1 [cs.CR] 28 Sep 2026
Gorka Abad1 , Claude Carlet1,4 , Ermes Franch1 , Stjepan Picek1,2,3 , and Vincent Rijmen1,5 1
University of Bergen, Bergen, Norway {ermes.franch,gorka.abad}@uib.no 2 University of Zagreb, Zagreb, Croatia 3 Radboud University, Nijmegen, The Netherlands [email protected] 4 Université Paris 8, Saint-Denis, France [email protected] 5 KU Leuven, Leuven, Belgium [email protected]
Abstract. We revisit cryptanalytic extraction of ReLU networks from a geometric and algebraic perspective. Rather than restricting attention to a single output component, we study the full vector-valued behavior across adjacent linear regions. This leads to a rank-one characterization of Jacobian differences that recovers the usual row-signature information while also revealing complementary column-side information. Our experiments show how this additional structure can be used in extraction and improves numerical estimation under different numericalprecision regimes (float64, float32 and float16). We extend the analysis beyond the high-precision and output-rounding settings commonly considered in the literature towards the low-precision settings encountered in many practical settings. Code is available online6 .
Keywords: cryptanalytic model extraction; ReLU networks; Jacobian discontinuities; parameter identifiability; row signatures; column signatures
1
Introduction
Machine learning (ML) and, in particular, deep learning (DL) can represent valuable intellectual property. Their development may require considerable computational resources, large and carefully collected datasets, and substantial engineering effort. This makes the confidentiality of their trained parameters an important aspect of machine-learning security. Model extraction attacks [14] are a direct threat to this confidentiality. These attacks aim to imitate as faithfully as possible a target model, exploiting black-box access to the same. In this 6
https://anonymous.4open.science/r/relu-signature-leakage-BE51/README.md
framework, the attacker can query the target model on chosen inputs and observe the outputs without knowing any other information about internal states or parameters. Initial approaches aimed at task-accuracy extraction using queries to the original model to train a surrogate that will imitate the original model on a similar task or on a similar data distribution [14]. A stronger objective is fidelity extraction, where the goal is to recover a model implementing the same function as the target (i.e., for any input, the output of the stolen model will be exactly the same as the original). Cryptanalytic attacks [6] may go further and aim at parameter extraction, recovering the underlying parameters up to transformations that leave the realized function unchanged. These approaches leverage the piecewise linear structure of some neural networks [12, 9, 6]. Following previous work [6], in this work we consider the setting where the attacker can choose queries and have access to the raw, vector-valued output of the network. Previous works analyze the transition through the boundaries of linear cells to recover the weights and the biases of the model. We instead consider the complete affine maps associated with two adjacent linear cells. If two cells differ only in the activation of one neuron, their Jacobians satisfy: ∆ = J1 − J2 = cr ⊺ , hence ∆ has rank 1. The rows of this matrix encode the same neuron information used in [6] and follow-up works, while the columns add information lost when observing only a single output component. Since all nonzero rows of ∆ are scaled copies of the same direction, they provide redundant observations of the neuron signature. Rather than treating these observations independently, we exploit the known rank-one structure of ∆. Given the noisy observation S = ∆+E, where E represents the error introduced by finite precision arithmetic, its dominant singular component is the best rankone approximation of S in Frobenius norm. We use the dominant right singular vector of S as an estimate of the common signature direction, combining the information carried by all output components. We also investigate two refinements intended to improve numerical precision, at the cost of additional queries. The first maximizes the interval used to estimate the Jacobian in each direction, similarly to [13], while the second adjusts an initial estimate of the neuron weights by probing along approximate tangent directions in adjacent linear cells. These refinements can improve estimation when the cell-membership test is reliable, but may degrade accuracy at lower precision when the test accepts intervals crossing cell boundaries. A side advantage of using the full raw output vector instead of a single component is that, when neuron transitions occur in the last hidden layer, the matrix ∆ exposes column signatures of the final weight matrix, giving a leakage channel complementary to the usual forward extraction of row signatures. These leaks do not directly strengthen end-to-end extraction attacks. Its value is mainly structural, clarifying what information the full vector response contains that scalar-output methods discard. Our contributions can be summarized as: 2
1. We show that a simple transition, in which exactly one ReLU neuron changes activation, induces a rank-one Jacobian difference ∆ = cr ⊺ (Eq. (11)). Its right factor recovers the row-side information exploited by previous extraction attacks, while its left factor exposes complementary column-side information. We characterize when these column signatures correspond to W (ℓ) and discuss the limitations of recovering biases and signs from the column side alone. 2. We use the rank-one redundancy of ∆ to construct an SVD-based rowsignature estimator that combines all output components. Since each oracle query already returns the complete raw-output vector, this additional information can be exploited without increasing the number of oracle queries. 3. We evaluate the resulting estimator against a Carlini-style single-output baseline under float64, float32, and float16 arithmetic. Our experiments show improved parameter recovery under finite precision, particularly in the float32 regime, and identify the cell-membership test as the main limitation of the adaptive refinements.
2
Preliminaries
2.1
Notation
We denote by R the set of real numbers. The vector space of n-tuples of real numbers is denoted by Rn . Vectors are written in bold lowercase letters, such as v ∈ Rn , and are interpreted as column vectors unless otherwise stated. Similarly, the transpose v ⊺ will be interpreted as a row vector. For a vector v, its i-th component is denoted by vi where i goes from 1 to n. Matrices with real entries, m rows and n columns, are elements of Rm×n , and are denoted by bold uppercase letters, such as M . The vector formed by the entries of the i-th row of M is denoted by mi,∗ , and the j-th column of M is denoted by m∗,j . Both are column vectors; hence, the i-th row of M is written as m⊺i,∗ . The transpose of a matrix M ∈ Rm×n will be indicated by M ⊺ ∈ Rn×m . Given two functions f (x), g(x), we will indicate their composition as (f ◦ g)(x) = f (g(x)). 2.2
Deep Neural Networks
A deep feedforward neural network consists of several layers. Algebraically, each hidden layer is the composition of an affine map and a component-wise activation function. Let n0 , . . . , nℓ be positive integers denoting the layer widths. For k = 1, . . . , ℓ, let A(k) : Rnk−1 → Rnk , A(k) (h) = W (k) h + b(k) , where W (k) ∈ Rnk ×nk−1 is the weight matrix and b(k) ∈ Rnk is the bias (k) (k) vector. We denote the entries of W (k) by wi,j and those of b(k) by bi . 3
Let h(0) = x be the network input and let σ denote the component-wise application of an activation function σ : R → R. For k = 1, . . . , ℓ − 1, define the preactivation vector z (k) ∈ Rnk and hidden activation vector h(k) ∈ Rnk by (1) z (k) = W (k) h(k−1) + b(k) , h(k) = σ z (k) . When soft labels are exposed, a softmax is typically applied to the final affine output, while hard-label models return its argmax. In our setting, we assume access to the raw output before either operation. Consequently, the last layer is affine, and its raw network output is F (x) = A(ℓ) (h(ℓ−1) ) = W (ℓ) h(ℓ−1) + b(ℓ) .
(2)
Definition 1 (ℓ-layer deep neural network). An ℓ-layer deep neural network is a function F : Rn0 → Rnℓ of the form F (x) = A(ℓ) ◦ σ ◦ A(ℓ−1) ◦ · · · ◦ σ ◦ A(1) (x).
(3)
The input layer is not counted among the ℓ layers. From now on, we restrict to the ReLU activation σ(x) = ReLU(x) = max{0, x}. (k)⊺
(k)
For k < ℓ, let wi,∗ denote the i-th row of W (k) . We denote by ηi R the i-th neuron of layer k, defined by (k) (k)⊺ (k) ηi (h) = ReLU wi,∗ h + bi .
: Rnk−1 →
For a network input x, the input to this neuron is h(k−1) (x). (k)
Definition 2 (Critical point [6]). The neuron ηi is said to be in an active state, inactive state, or critical state at h ∈ Rnk−1 if (k)⊺
(k)
wi,∗ h + bi is positive, negative, or zero, respectively. A point h ∈ Rnk−1 satisfying (k)⊺
(k)
wi,∗ h + bi
=0
(4)
is called a critical point. The set of all critical points is the critical hyperplane n o (k) (k)⊺ (k) . Hi := h ∈ Rnk−1 | wi,∗ h = −bi (k)⊥
The critical hyperplane is parallel to the (nk−1 −1)-dimensional space wi,∗
:=
(k)⊺ (k) {h ∈ R | wi,∗ h = 0}. In particular, the vector wi,∗ is orthogonal to the (k) underlying space of the hyperplane Hi . This observation was exploited in [5] to (k)⊺ recover the row vectors wi,∗ up to some scalar. An equivalent way of extracting nk−1
these vectors was previously developed in [6] through differential techniques. 4
2.3
Invariant Transformation of Neural Network
In Definition 1, we defined F as a composition of affine functions A(k) (x) = W (k) x+b(k) interleaved by activation functions σ. Since σ acts component-wise, for any permutation matrix P ∈ Rnk ×nk we have σ(P h) = P σ(h). Therefore, for every hidden layer k < ℓ, the function F is preserved under the transformation F (x) = A(ℓ) ◦ · · · ◦ (A(k+1) ◦ P −1 ) ◦ σ ◦ (P ◦ A(k) ) ◦ · · · ◦ A(1) (x). When the activation function is ReLU, there is another transformation of the affine functions that preserves F . Indeed, ReLU(ax) = a ReLU(x) for every a > 0. More generally, let k < ℓ and let D ∈ Rnk ×nk be a diagonal matrix with strictly positive diagonal entries. Since ReLU acts component-wise, σ(Dh) = Dσ(h). Therefore, the function represented by the network is unchanged when the parameters of two consecutive layers are transformed as f (k) = DW (k) , W
e b(k) = Db(k) ,
f (k+1) = W (k+1) D −1 . W
Equivalently, the rows of W (k) and the corresponding entries of b(k) may be scaled by arbitrary positive factors, provided that the corresponding columns of W (k+1) are scaled by their reciprocals. This transformation does not apply to the output layer, because there is no subsequent layer in which to compensate for the rescaling. 2.4
Signatures
Row signatures are defined in [3] as follows: (k)
(k)⊺
Definition 3 (Row Signature [3, Definition 4]). Let ηi (h) = σ(wi,∗ h + (k)
(k)
bi ) be a neuron on the k-th layer. The row signature of ηi (k)
(k)
(k) (k)⊺ (wi,1 )−1 wi,∗ =
1,
wi,2
(k)
wi,1
,...,
wi,nk−1
is the row vector
!
(k)
.
wi,1
That is, the i-th row of W (k) is normalized on its first component. (k)
Notice this normalization requires the first component wi,1 to be nonzero. Since the set of such vectors has measure zero (or, rather, a small measure, because of finite precision) and it is extremely unlikely to happen, there is no need to address this particular issue. In the same way, one can define a column signature for the columns of W (k) normalized on their first components. In this paper, we will use a slightly different definition of row signature: 5
(k)
(k)⊺
(k)
Definition 4 (Row Signature (this work)). Let ηi (h) = σ(wi,∗ h + bi ) (k)
be a neuron on the k-th layer, and assume that wi,∗ ̸= 0. The row signature (k)
of ηi
is the row vector (k)
sign(wi,t ) (k)
∥wi,∗ ∥2
(k)⊺
wi,∗ ,
(k)
where t = min{r : wi,r ̸= 0} and sign is the function outputting +1 if the input is positive and −1 if it is negative. That is, the i-th row of W (k) is scaled to have norm 1 whose first nonzero coordinate is positive. (k)
In a similar way, for a nonzero column w∗,j , its column signature can be defined as the vector (k)
sign(wt,j ) (k) ∥w∗,j ∥2
(k)
w∗,j ,
(k)
t = min{r : wr,j ̸= 0}.
One advantage of the modified definitions is that they also work for rows or columns starting with a zero element. Additionally, their computation does not incur numerical issues when the rows or columns start with elements close to zero. Using the invariants described in the previous section, it is possible to construct different algebraic descriptions of the same neural network. Consequently, we may use positive rescaling to normalize either the columns of W (2) , . . . , W (ℓ) or the rows of W (1) , . . . , W (ℓ−1) . Positive rescaling does not reverse orientation. Therefore, after normalization, each row or column is equal to either its canonical signature or the negative of that signature. The signature of a row determines a unit normal direction for the associated critical hyperplane, but it does not determine the offset of that hyperplane; the offset also depends on the corresponding normalized bias. Moreover, even when the hyperplane is known, the canonical signature does not determine which side is active, because reversing both the weight vector and the bias preserves the hyperplane while exchanging its two sides. Recovering this orientation sign is therefore necessary for functional reconstruction. For a column signature, the missing sign analogously determines the direction of the neuron’s contribution to the next layer.
3
Related Work
Query-based model reverse engineering predates modern neural-network extraction. Lowd and Meek [11] studied adversarial query strategies for reverse engineering linear and Boolean classifiers. In 2016, Tramèr et al. pioneered the first attack methods based on synthetic data generation (i.e., interrogating the model on chosen inputs to generate pairs 6
(x, f (x)) that can be used to train a surrogate model) and simple query techniques. More precisely, early model extraction attacks typically relied on synthetic data generation or simple interrogation methods to mimic the model’s behavior [14]. However, these methods have proved to be insufficient to accurately replicate the actual model parameters. Jagielski et al. [9] followed a different approach that made significant progress towards the extraction goal with high accuracy and high fidelity, although their methods still lacked the query efficiency offered by subsequent cryptanalytic techniques. On the theoretical side, Fefferman [8] showed that, under suitable assumptions and with exact knowledge of the realized function, neural-network parameters can be reconstructed. In an orthogonal line of work, Batina et al. [1] used electromagnetic side-channel leakage to recover neural-network architecture and parameters. The attacks closest to ours exploit the piecewise linear structure of ReLU networks. Milli et al. [12] used gradient information exposed by model explanations to reconstruct network parameters, while Jagielski et al. [9] developed high-fidelity extraction from standard model outputs. A significant milestone in this field came with the 2020 paper by Carlini et al. [6], where the problem was framed as a cryptanalytic problem and showed how critical points and derivative discontinuities can reveal neuron signatures with high precision. Compared to [9], the attack in [6] achieved 220 times more accurate results with 100 times fewer queries. However, the difficulty of identifying the signs (positive or negative) of neurons in deep networks and the need for an excessive number of input-output pairs to obtain a highly accurate model limited the practical applicability of the attack. Canales-Martı́nez et al. developed techniques, among which one they named neuron wiggling [2] to solve the problem of identifying the signs of hidden neurons in deep networks and thus reduce the time complexity of the attack to a polynomial level. It is important to note that both attacks mentioned above assume the attacker has access to raw output. Hard label extraction was later studied by Chen et al., who proposed an extraction method that requires a polynomial number of queries but suffers from exponential runtime [7]. This theoretical challenge was subsequently overcome by Carlini et al. [4]. The authors demonstrated that parameters could be extracted by analyzing the bending of decision boundaries (i.e., points at which the model output switches from one class to another) in proximity to a neuron switch, thereby achieving polynomial-time complexity and making hard-label attacks practically feasible. More recently, Liu et al. focused on extracting signatures from networks [10] that are deeper than the ones achieved by [6] and [3]. Our work takes a complementary view of this cryptanalytic extraction line. Rather than restricting to a single output component around a critical point, we consider the full vector-valued affine behavior of adjacent linear cells. This reveals a rank-one Jacobian difference whose rows recover the usual neuron signature while its columns expose additional information about later layers. A second motivation is numerical: existing extraction techniques rely on highly accurate signature recovery, and finite-precision evaluation can become a signif7
icant source of error, particularly in reduced-precision settings such as float32 and float16. The full Jacobian difference provides several redundant observations of the same row direction, suggesting that rank-one approximation or averaging may improve the stability of signature estimation. We investigate both this potential numerical gain and the additional column-side information, while also identifying the limitations of this approach.
4
Two-Side Signature Recovery
An important property of DNNs with ReLU, or another piecewise-linear activation function, is that the function they realize is piecewise affine (often called piecewise linear in this context). We will see how this property allows us to extract the neuron signatures. Definition 5 (Piecewise linear function). A function F : Rn → Rm is said to be piecewise linear if there exists a finite collection C = {C1 , . . . , Ct } of St polyhedra,n each described by a finite system of affine inequalities, such that i=1 Ci = R and F|Ci (x) = J|Ci x + b|Ci , where F|Ci denotes the restriction of F to the set Ci and J|Ci ∈ Rm×n , b|Ci ∈ Rm . In other words, the restriction of F to the sets Ci is an affine function. The sets Ci are also known as linear cells. We define the rank of F|Ci as Rank(J|Ci ). If F is continuous and two full-dimensional cells Ci , Cj share a codimension-one boundary, their affine restrictions agree on that boundary. Consequently, every point x on the shared boundary satisfies (J|Ci − J|Cj )x + b|Ci − b|Cj = 0, For a ReLU network, such boundary pieces arise when a neuron preactivation is zero. In the network input space, they are preimages of the critical hyperplanes introduced in Definition 2 under the subnetwork preceding the corresponding neuron. We refer to these preimages as the cell boundary. The composition of affine functions and piecewise linear functions results in piecewise linear functions whose linear cells are delimited by portions of hyperplanes. Recall from Definition 1 that the function implemented by a neural network can be described as F (x) = A(ℓ) ◦ σ ◦ · · · ◦ A(2) ◦ σ ◦ A(1) (x),
(3)
where A(k) (x) = W (k) x + b(k) for some weight matrices W (k) ∈ Rnk ×nk−1 and bias vectors b(k) ∈ Rnk , hence they are affine functions. Let σ be the component-wise application of ReLU, which is a piecewise linear function; then the function F is also a piecewise linear function. The borders of the linear cells of F are determined by the activation function. 8
F|C1 F|C2 F|C3
F|C4
F|C5
F|C7
F|C6
Fig. 1. A subdivision of R2 into linear cells. On each cell Ci , the restriction F|Ci of F is affine.
Once we fix a particular input, this will determine a particular configuration of active-inactive neurons. Each of these patterns can be represented by substituting the function σ at layer k with a diagonal matrix whose diagonal (k)⊺ (k) entries are 1 if wi,∗ h(k−1) (x) + bi > 0 (i.e., the corresponding neuron is active) and 0 otherwise. For a linear cell C, we denote this diagonal matrix by (k) DC ∈ Rnk ×nk . Each point in the input space x ∈ Rn0 lies within a linear cell or on the border between two or more cells; inside the same linear cell, all the (k) diagonal matrices DC are constant. As the borders are constituted by portions of hyperplanes having measure zero, a random point x will be inside a linear cell with probability 1.7 We denote the linear cell containing the vector x0 by C0 , and the restriction of the function F to the cell C0 by F|C0 . Locally (if we restrict the domain to C0 ⊆ Rn0 ), the function F is a composition of affine functions, and it is therefore affine. Using this notation, we define: (ℓ−1)
J|C0 = W (ℓ) DC0 Li,C0 = W b|C0 =
(ℓ)
ℓ−1 X
(1)
· · · DC0 W (1) ∈ Rnℓ ×n0
(ℓ−1) DC0 · · · W (i+1) , i = 1, . . . , ℓ − 1 (i)
Li,C0 DC0 b(i) + b(ℓ)
(5) (6) (7)
i=1
the restricted function F|C0 can be expressed as F|C0 (x) = J|C0 x + b|C0 .
(8)
Informally, we will refer to the rank of the matrix J|C0 as the rank of F|C0 . The affine function F|C0 : Rn0 → Rnℓ , can be reconstructed choosing n0 + 1 affine independent points in C0 where n0 is the dimension of the input space. Provided we have access to the evaluation of the function at these points, and that all these points belong to the same linear cell C0 , it is always possible to 7
In practice, this probability will be only close to 1 as we cannot have infinite precision.
9
determine the function F|C0 and the matrix J|C0 . Another interpretation of the matrix J|C0 is that it is the Jacobian of the function F|C0 that is ∂F1|C , 0
∂e ∂F2|C1 0 , ∂e1
J|C0 =
∂F1|C0 ∂e2 ∂F2|C0 ∂e2
.. .
.. .
∂Fnℓ |C0 , ∂Fnℓ |C0 ∂e1 ∂e2
∂F1|C0 ∂en0 ∂F2|C0 ∂en0
··· ··· .. . ···
.. .
∂Fnℓ |C0 ∂en0
,
(9)
∂F
i|C0 , where ∂e is the partial derivative in direction ej of the i-th component of j the function F|C0 and we denote by ej the unit vector of the standard basis having j-th component 1 and being zero elsewhere. This simple fact becomes crucial when we compare the function in two adjacent linear cells, that is, two linear cells separated by a portion of a cell boundary. Let C0 and C1 be two adjacent linear cells that differ only in the activation of the (i) (i) j-th neuron of the i-th layer. The diagonal matrices satisfy DC0 − DC1 = ±Ej,j ni ×ni where Ej,j ∈ R is the matrix that is 0 in each entry except in the position (j, j) where it is 1. It is useful to introduce some notation for the matrix products before and after the activation at layer i. For a linear cell C0 , define
(i−1)
Ri,C0 = W (i) DC0
(ℓ−1)
Li,C0 = W (ℓ) DC0
(1)
W (i−1) · · · DC0 W (1) ∈ Rni ×n0 , W (ℓ−1) · · · W (i+1) ∈ Rnℓ ×ni .
We can use this to factor the matrix J|C0 as (i)
J|C0 = Li,C0 DC0 Ri,C0 .
(10)
The only difference between J|C0 and J|C1 is in a single activation of i-th layer. Without loss of generality, we label the cells so that this neuron is active in (i) (i) C0 and inactive in C1 . Therefore, DC0 − DC1 = Ej,j , and the difference between the two Jacobians is given by ⊺ ∆(C0 , C1 ) = J|C0 − J|C1 = Li,C0 Ej,j Ri,C0 = c∗,j rj,∗ ∈ Rnℓ ×n0 ,
(11)
⊺ where c∗,j is the j-th column of Li,C0 and rj,∗ is the j-th row of Ri,C0 . Each row of the matrix ∆(C0 , C1 ) is a scalar multiple of the row vector ⊺ rj,∗ , while each column is a scalar multiple of the column vector c∗,j . In terms of Definition 3, every row of ∆(C0 , C1 ) has the same row-signature and every column has the same column-signature. There are two special cases in which we can exploit this decomposition to recover the row-signatures or the columnsignatures of the matrices W (i) .
1. In the first case, the activation patterns for C0 and C1 differ for the jth neuron in the first layer. Since R1,C0 = W (1) , each row of the matrix 10
∆(C0 , C1 ) is a scalar multiple of the j-th row of W (1) . This is exactly the same case exploited in [6] and its follow-up works. The connection is even more evident if we recall that J|C0 is the Jacobian of the local affine function F|C0 (x). 2. In the second case (new), the activation patterns differ in the j-th neuron of the last hidden layer. Since Lℓ−1,C0 = W (ℓ) , the columns of ∆(C0 , C1 ) are scalar multiples of the j-th column of W (ℓ) . As already stated, the first case was exploited to collect all the signatures of the first layer and from there attack the second layer, and so on. Then, a natural question arises: Can we follow a similar approach, collecting the column signature of the last layer, and proceed backward towards the first? If successful, then the two approaches could be easily combined, reconstructing the function from both ends. 4.1
Finding adjacent linear cells
To determine signatures, it is essential to find adjacent linear cells. Consider an input x0 and a unit vector v, we can choose a step length δ > 0 and consider the sequence of points S(x0 , v) = {xt = x0 + tδv}t∈Z . We denote the finite differences of consecutive elements of this sequence as Dv,δ F (xt ) =
F (x0 + (t + 1)δv) − F (x0 + tδv) F (xt+1 ) − F (xt ) = . δ δ
(12)
If the interval [xt , xt+1 ] is contained in a single linear cell C0 of the piecewise linear function F , then Dv,δ F (xt ) coincides with the directional derivative of F along v at xt , which we denote by Dv F (xt ). In particular, Dv,δ F (xt ) = Dv F (xt ) = J|C0 v. If the interval crosses a cell boundary, Dv,δ F (xt ) remains a finite difference and does not need to coincide with the directional derivative at xt . Within the linear cell C0 the function F can be expressed as F|C0 (x) = J|C0 x+b|C0 , its finite difference along v is the constant vector Dv,δ F (x) = J|C0 v. Suppose x0 ∈ C0 , and let xt1 −1 ∈ C0 be the last element in the sequence S(x0 , v) that belongs to the same cell, we will have Dv,δ F (x0 ) = Dv,δ F (x1 ) = · · · = Dv,δ F (xt1 −2 ) = J|C0 v. As xt1 ∈ C1 while xt1 −1 ∈ C0 there will be some value of 0 ≤ λ < 1 for which xt1 −1 + λδv lies exactly on the cell boundary dividing C0 from C1 , then Dv,δ F (xt1 −1 ) = (λJ|C0 + (1 − λ)J|C1 )v ̸= J|C0 v,
(13)
except in the unlikely case where J|C0 v = J|C1 v. From xt1 onward, the finite difference Dv,δ F (xt1 ) will be again a constant vector given by J|C1 v as long as we stay within C1 . If the sequence leaves C1 at a 11
certain point xt2 , we will detect this transition by observing another discontinuity in the value of Dv,δ F . Denote by d the difference between the directional derivatives of the affine restrictions of F to the two cells: d = J|C0 v − J|C1 v = ∆(C0 , C1 )v. With a dichotomy, it is possible to identify the location of the critical point where this discontinuity is located [6]. Alternatively [3], this point can be located as xt1 −1 +λδv, solving for lambda in equation (13) noticing that the values of J|C0 v and J|C1 v are known as they can be calculated from the evaluations of F . Observe that [6] and subsequent works restrict the function F to a single component. In that case, the finite difference is a scalar, corresponding to the same component of the vector Dv,δ F (xt ). Transitions from one linear cell to another are detected when this scalar changes. Clearly, it is possible for two points xt0 and xt1 to be in different linear cells, where the vectors Dv,δ F (xt0 ) ̸= Dv,δ F (xt1 ), while the same vectors are equal in the component to which F has been restricted. In this case, the method in [6] fails to distinguish C0 , C1 as two distinct linear cells. This event has measure zero, but it might happen that the two scalars are so close that they cannot be reliably detected. 4.2
Extracting a row signature
In [6] and following works, signatures are extracted by comparing the derivative on both sides of a cell boundary on a single component q of the output vector Fq (x) = [F (x)]q ∈ R. Before describing our technique, we will briefly describe how signatures are recovered in [3]. Let ei be the standard basis unit vector that has 1 in its i-th entry and zero elsewhere, and let x∗ be a point on the cell boundary. They define αi,+ = Dei ,δ Fq (x∗ + δei ) and αi,− = Dei ,δ Fq (x∗ − δei ). For the first layer, they consider the signature αn0 ,+ − αn0 ,− α1,+ − α1,− α2,+ − α2,− , ,..., , αr,+ − αr,− αr,+ − αr,− αr,+ − αr,− and by choosing r = 1, they will obtain a signature whose first component is 1. It can be verified that, if no other neuron flips on the points used to estimate the values of αi,+ , αi,− , the signature obtained is a scalar multiple of the weight of the neuron. (k) Once the weights of a neuron ηi are recovered, its bias can be determined from Equation (4) in Definition 2. For convenience, we recall that (k)⊺
(k)
−wi,∗ h(k−1) (x) = bi ,
(14)
where h(k−1) (x) is a point that belongs to the critical hyperplane associated (k) with the neuron ηi and x is a point that belongs to a cell boundary of two (k) adjacent cells differing by the activation of ηi . Further details on the recovery of such a point, and hence of the bias, are given in Section 4.4. 12
Our method: extraction from the Jacobian In the discussion following (11), we already observed how the Jacobians of two adjacent linear cells C0 , C1 differ by a rank 1 matrix. If the neuron that flipped between C0 and C1 is in the first layer, the rows of ∆(C0 , C1 ) = J|C0 − J|C1 will all be scalar multiples of the row in the matrix W (1) corresponding to the weights of the neuron that flipped. To compute the Jacobians J|Ci we need to find a point x∗ ∈ Ci that is at a distance of at least δ from any cell boundaries. In this way, we can query the model on x∗ and on neighbor points x∗ + δvj , where v1 , . . . vn0 are linearly independent unit vectors, granted that all these points are still in the same linear cell. Assume that both x∗ , x∗ +δvj ∈ Ci and compute yj = F (x∗ +δvj )−F (x∗ ) = δJ|Ci vj , we call Y the matrix whose columns are given by the vectors yj and V the matrix whose columns are the vectors vj . Then J|Ci =
1 Y V −1 , δ
a natural choice is to choose vj = ej the standard basis, in this way the matrix V is the identity. It is crucial that all the points x∗ + δvj belong to the same cell Ci . In Section 4.1 we introduced the sequence S(x0 , u). Let x0 ∈ C0 , we name xt1 the first element in C1 , xt2 the first element in C2 and so on. A first good candidate is x +x given by x∗ = t1 2 t2 −1 , more generally, a good point to compute J|Ci is given by the midpoint between xti and xti+1 −1 . This point will be almost as far as possible from the boundaries of cell C1 along the line defined by x0 and the vector u. Unfortunately, this will not prevent the same point from being close to any other cell boundary. A practical consistency test for x∗ and x∗ + δvj is to compare finite differences along one or more probe directions us : F (x∗ + µus ) − F (x∗ ) F (x∗ + δvj + µus ) − F (x∗ + δvj ) − ≤ τ. µ µ In exact arithmetic τ = 0; using several probe directions reduces the chance of accepting points from different cells. Sign Similar to other techniques, our signatures will miss a global sign, in [6], [3], and [10], several techniques are proposed to recover this last missing information about each neuron. These techniques can be applied to our case in the same way as reconstructing the first layer. 4.3
Attacking deeper layers
Once the first layer is reconstructed, we can proceed to attack the second. The process will be the same for each subsequent layer, so we will assume that we have successfully reconstructed the first k − 1 layers and that we want to attack the k-th layer. 13
In our model, a signature is recovered as normalized rows and columns of ∆(C0 , C1 ) = J|C0 −J|C1 , where C0 , C1 are adjacent linear cells (i.e., cells differing only by one activation). Suppose the neuron whose activation determines the separation between C0 , C1 is the j-th neuron of the k-th layer and that we have already successfully reconstructed the first k − 1 layers. From (11) we have (k−1)
∆(C0 , C1 ) = J|C0 − J|C1 = Lk,C0 Ej,j Rk,C0 = Lk,C0 Ej,j W (k) DC0
Rk−1,C0 ,
Notice that, the matrix Lk,C0 Ej,j W (k) will be a rank 1 matrix whose rows are (k)⊺ all scalar multiples of wj,∗ , the j-th row of W (k) and whose columns are all scalar multiples of the j-th column of Lk,C0 . The matrix Rk−1,C0 is known, as we have already reconstructed the first k−1 layers. Then, if this matrix has rank nk−1 , it is possible to find an invertible (k−1) sub-matrix nk−1 × nk−1 and we can recover the matrix Lk,C0 Ej,j W (k) DC0 . Notice that this last operation hides two problems. One is that the matrix Rk−1,C0 might not have the desired rank, and the other is that the columns (k−1) of Lk,C0 Ej,j W (k) DC0 will be zero in correspondence with the inactive neurons of the (k − 1)-th layer. The first problem is known as insufficient rank signatures and the second as partial signatures. We will treat them separately in the following two paragraphs. Partial signatures Partial signatures emerged in a slightly different way in [6] and subsequent works. The solution proposed in [6] can be used for partial signatures extracted from the Jacobian. For the sake of clarity, we will briefly describe the method proposed in [6] to merge partial signatures in our framework. Under the assumption that Rk−1,C0 ∈ Rnk−1 ×n0 has rank nk−1 , the signa(k)⊺ ture we are recovering will be the j-th row of W (k) , denoted by wj,∗ , multiplied (k−1)
by the matrix DC0
. Unless all the neurons in the (k − 1)-th layer are active,
(k)⊺ (k−1) the row vector wj,∗ DC0 will have some zeros corresponding to the zeros in (k−1) DC0 . Notice that, if we observe a partial signature relative to the same j-th
neuron of the k-th layer from a different pair of adjacent cells Ci , Ci+1 , the matrix (k−1) DCi will likely reveal some coordinates that were previously zero and mask others with zeros. In [6], a technique was already proposed to identify and merge partial signatures coming from the same neuron. This method can be illustrated with a simple example. Consider two partial signatures v = (v1 , v2 , 0, v4 , 0)⊺ and u = (u1 , u2 , u3 , 0, u5 )⊺ that have two non-zero entries in the first two components. If the two vectors are a scaled and masked version of the same row (k)
vector wj,∗ , then uu12 =
(k)
wj,1
(k)
wj,2
= vv12 . This simple test gives a criterion to establish
whether two signatures are relative to the same neuron and can be merged. The implicit assumption is that different neurons will have different ratios between the overlapping components. To merge the two signatures, each vector is normalized by one of their common nonzero components, and the missing entries of one signature are filled using the corresponding nonzero entries of the other. In 14
the example above, normalizing with respect to the first component yields v1 v2 u3 v4 u5 , , , , . v1 v1 u1 v1 u1 Insufficient rank signatures The problem of insufficient rank signatures arises when the matrix Rk−1,C0 has rank lower than nk−1 . This problem is more severe than partial signatures as it can prevent the recovery of the signature of a neuron in the k-th layer. A solution to this problem was proposed in [10] combining linear equations coming from distinct portions (distinct pairs of adjacent cells) of cell boundaries relative to pairs of adjacent cells differing in the activation of the same neuron. Instead of explaining the method in [10], we will describe a similar method adapted to our framework. Let us denote by ∆1 , ∆2 , . . . , ∆m the matrices ∆i = ∆(Ci,0 , Ci,1 ) = J|Ci,0 − J|Ci,1 relative to m distinct pairs of adjacent cells Ci,0 , Ci,1 that differ only in the activation of the same neuron in the k-th layer. All these matrices will have rank 1, and each will satisfy an equation of the form (k)⊺
∆i = v|Ci,0 wj,∗ M|Ci,0 , (k−1)
where v|Ci,0 is the j-th column of Lk,Ci,0 and M|Ci,0 = DCi,0 Rk−1,Ci,0 ∈ (k)⊺
Rnk−1 ×n0 . Noticing that both ∆i and v|Ci,0 wj,∗
are rank 1 matrices, we can (k)⊺
write the equation above as a system of linear equations in the unknown wj,∗ , ignoring scalar multipliers as (k)
⊺ M|C wj,∗ = λdi i,0 (k)
where di is a row of the matrix ∆i up to scalar multiplication. Treating wj,∗ as the unknown, the solution to this system will be given by the vector space (k) ⊺ Vi = ⟨wi ⟩+ker(M|C ), where wi is a particular solution to the system. As wj,∗ i,0 is a solution to all the systems, it will belong to the intersection of all the vector ⊺ spaces Vi . Notice that the matrix M|C has zero columns corresponding to the i,0 ⊺ inactive neurons of the (k − 1)-th layer, so the kernel of M|C contains at least i,0 the space spanned by ej for each inactive neuron j of the (k − 1)-th layer. We ⊺ can recover a partial signature if the kernel of the matrix M|C from which we i,0 remove the zero columns is trivial. If this kernel is not trivial, we can intersect with a different solution, keeping in mind to remove the zero columns where the two matrices share zero columns and keep the zero columns where only one of the two matrices has a zero column. A priori, we do not know if ∆i and ∆j are relative to the same neuron, but we can test them in the same way proposed in [10]. In particular, the intersection of the two vector spaces Vi ∩Vj is expected to be trivial if the sum of their dimension is smaller than the number of variables (i.e., is smaller than the number of surviving columns in the two matrices). If the intersection has dimension 1, most likely we have recovered the signature even if the matrices M|Ci,0 and M|Cj,0 were rank deficient. 15
4.4
Extracting bias
Similar to what is done in [6], after we successfully recovered a complete signature (k) for a neuron ηi , it is possible to recover its bias from the following equation (deduced from Relation (14)): (k)⊺
(k)
−wi,∗ h(k−1) (x∗ ) = bi
(15)
where h(k−1) (x∗ ) ∈ Rnk−1 is a point that belongs to the critical hyperplane of (k) ηi . Since a signature determines the weight row only up to scale and sign, we normalize the recovered row as (k)
wi,∗
(k)
bi,∗ = w
(k)
∥wi,∗ ∥2
.
Consequently, Equation (15) recovers the normalized bias (k)
bb(k) = i
bi
(k)
∥wi,∗ ∥2
, (k)
up to the selected sign, rather than the physical parameter bi . Suppose we have already reconstructed the neural network up to layer k and let x1 , x2 be two points belonging to two adjacent cells differing in the activation (k) of neuron ηi . That is, the two cells are divided by the portion cell boundary (k) (k) (h(k−1) )−1 (Hi ), where Hi ⊆ Rnk−1 is the critical hyperplane relative to the (k) neuron ηi (see Definition 2). We assume that the segment [x1 , x2 ] crosses (k) exactly one cell boundary, at a point x∗ ∈ (h(k−1) )−1 (Hi ). Since x1 and x2 lie on opposite sides of the same critical boundary, one can search along the segment joining them for the point at which the local affine behavior of F changes. The continuity of F guarantees that the two affine pieces agree at the critical point. Let x∗ = x1 + t∗ (x2 − x1 ),
t∗ ∈ [0, 1],
(16)
and denote by J1 and J2 the Jacobians of F in the cells containing x1 and x2 , respectively. The continuity of F ensures that F (x1 ) + t∗ J1 (x2 − x1 ) = F (x2 ) − (1 − t∗ )J2 (x2 − x1 ). Rearranging, we get t∗ (J1 − J2 )(x2 − x1 ) = F (x2 ) − F (x1 ) − J2 (x2 − x1 ). Denoting by y the vector on the right-hand side of the equation and by a the vector (J1 − J2 )(x2 − x1 ), we get the simpler equation t∗ a = y. In [6], a single output coordinate (for example the first), is used to determine t∗ as t∗ = ay11 . This formula is correct if all computations are exact. Since we 16
typically work with finite numerical precision, the equation t∗ = ay11 might have no exact solution. In this case, it is usually better to determine the least-squares solution b t = arg min ∥ta − y∥22 , t
which is given by ⊤
a y b t= ⊤ . a a Hence, coordinates with larger |ar | contribute more strongly while coordinates for which ar is close to zero are naturally downweighted, since small numerical perturbations in yr are strongly amplified in the ratio yr /ar . 4.5
Numerical precision
Enhance precision Due to numerical precision, the recovery of the signatures will always have some error. These errors will accumulate when we explore deeper layers; this affects both the fidelity of the reconstructed model and the capacity to reliably extract deep layers. In previous works, signatures were collected as single row-signatures and grouped by similarity, fixing a threshold distance. In our case, the matrix ∆(Ci , Cj ) = J|Ci − J|Cj will contain rows (or columns) all referring to the same signature. Rather than extracting these row-signatures independently, we apply Singular Value Decomposition (SVD) to ∆(Ci , Cj ) and estimate their common direction from its dominant right singular vector. This allows the information contained in all output components to be combined into a single signature estimate, improving the precision of our extraction. If the entries correspond to non-adjacent cells, the matrix is no longer approximately rank-1; this is reflected in the presence of multiple significant singular values. We propose two further techniques to improve the numerical precision. Adaptive steps Inside a linear cell C, the Jacobian J|C is constant, so a directional derivative can be computed using any nonzero interval contained entirely in that cell. In particular, for a unit vector v and h+ , h− ≥ 0 with h+ + h− > 0, if [x − h− v, x + h+ v] ⊆ C, then in exact arithmetic F (x + h+ v) − F (x − h− v) = J|C v. h+ + h−
(17)
If each endpoint evaluation has an additive error of norm at most η, its contribution to the derivative error is bounded by 2η/(h+ + h− ), before accounting for rounding in the subtraction and division. This motivates using larger steps. This principle was also used by Noman et al. [13, Sec. 4.1] in a study of deterministic output rounding as a defense against cryptanalytic extraction. Their 17
step-spacing attack adaptively enlarges symmetric finite-difference intervals, subject to a local linearity test, to recover directional derivatives despite output rounding. The defense rounds only the returned outputs, while retaining fullprecision weights and forward computation. Our construction applies the same large-step principle to vector-valued Jacobian estimation, allowing asymmetric intervals and using the baseline Jacobian for consistency checks. We first obtain a baseline estimate Jb|C using the small-step procedure described above. Once this estimate is available, we can test whether a candidate point z is locally consistent with the same Jacobian. For a unit probe direction v and a small step δ, we check whether F (z + δv) − F (z) − δ Jb|C v ≤ τ,
(18)
where τ accounts for both evaluation errors and uncertainty in the baseline estimate. A single probe requires at most two additional oracle queries; several probe directions (and correspondingly more queries) can be used to reduce the risk of accepting a point with a different local Jacobian. This is a heuristic consistency test, not a certificate of cell membership: distinct cells can agree along the tested directions, and endpoint tests alone cannot exclude intervening boundaries. Its reliability also depends on the quality of the baseline estimate. Starting from a point x in C, we extend a segment along each coordinate direction ei . On the positive side, we test x + δei , x + 2δei , doubling the step at each iteration, stopping at the first failed consistency test or at a prescribed + search limit. We retain the last accepted displacement x+2ki δei in one direction and repeat the same procedure independently in the opposite direction, retaining − x − 2ki δei . The two points need not be equally distant from x. For convenience + − we rename hi,+ = 2ki δ and hi,− = 2ki δ, we then estimate the i-th column of the Jacobian by F (x + hi,+ ei ) − F (x − hi,− ei ) Jb|C ei = , hi,+ + hi,−
(19)
provided the total length hi,+ + hi,− is sufficiently large. A finite search limit allows the procedure to terminate when the cell is unbounded in a tested direction. We require each coordinate interval to have length at least a prescribed threshold hi,+ + hi,− > Lmin . We retain the estimates for coordinates whose intervals satisfy this criterion and retry only the remaining coordinates from other points consistent with the same cell, using at most T candidate starting points in total. If, after these attempts, any coordinate still lacks a sufficiently long interval, we discard this cell for the current signature-recovery attempt. This does not rule out recovering the same neuron from another pair of adjacent cells. The length threshold is a numerical safeguard rather than a guarantee of a particular estimation error. Rotating adjustment As a further refinement, we use directional probes to correct an initial estimate of the row signature. For two adjacent cells C0 , C1 with a 18
nonzero rank-one Jacobian difference, the matrix that would be obtained in the absence of numerical errors has the exact factorization ∆(C0 , C1 ) = J|C0 − J|C1 = λan⊺ ,
∥a∥ = ∥n∥ = 1,
λ > 0,
(20)
where n is the true row-signature direction and a is the direction of the corresponding change in the output derivative. The preceding estimation procedure, b a b , and n b of these quantities. For followed by SVD, provides initial estimates λ, a unit direction v, the difference between the directional derivatives in the two b can reveal a cells is ∆(C0 , C1 )v = λa(n⊺ v). Thus, a direction orthogonal to n residual component of the true signature that the initial estimate misses. To measure this difference, we choose a point xj in each cell Cj , j ∈ {0, 1}, and nonnegative displacements hj,+ (v), hj,− (v) such that [xj − hj,− (v)v, xj + hj,+ (v)v] ⊆ Cj ,
Lj (v) = hj,+ (v) + hj,− (v) > 0.
The centers and lengths may be chosen separately for each direction, using the maximal-interval search described in the previous paragraph. Define Dj (v) =
F (xj + hj,+ (v)v) − F (xj − hj,− (v)v) , Lj (v)
d(v) = D0 (v) − D1 (v) = λa(n⊺ v).
(21) (22)
The last equality holds in exact arithmetic when both intervals remain in their respective cells. Dividing each cell’s output increment by its own interval length allows unequal and asymmetric intervals to be used. If each endpoint evaluation in Cj has error of norm at most ηj , the evaluation-error contribution to the directional jump is bounded by 2η0 2η1 + , L0 (v) L1 (v) before accounting for arithmetic rounding. Thus, for fixed endpoint-error bounds, longer intervals improve this bound, although a short or noisy interval on one side can dominate it. b to an orthonormal basis {b We complete n n, t1 , . . . , tn0 −1 } of the input space. b , while an additional Probing along the ti measures the components missed by n b b supplies a reference component. Let d(v) probe along n denote the numerical estimate of (22). Consider the ratios rbi =
b i) b ⊺ d(t a , b n) b ⊺ d(b a
i = 1, . . . , n0 − 1.
(23)
If we had access to exact directional measurements d(ti ) = λa(n⊺ v), the comb ), provided both the mon factor λb a⊺ a would cancel, giving ri = (n⊺ ti )/(n⊺ n common factor and the reference component are nonzero. This would yield the corrected direction P 0 −1 b + ni=1 n rbi ti b corr = q (24) n Pn0 −1 2 . 1 + i=1 rbi 19
In exact arithmetic, this recovers n up to the usual global sign; in finite precision, it can still give a candidate correction to the initial estimate. Long intervals along the approximate tangent directions may improve these measurements, but the cell-consistency checks remain heuristic. When the initial b and a b are well aligned with their true directions, the exact reference estimates n signal b ⊺ d(b b) a n) = λ(b a⊺ a)(n⊺ n has magnitude close to λ. Under this assumption, angular misalignment is not expected to make the denominator in (23) small. The reference signal may nevertheless be weak relative to measurement uncertainty if the Jacobian jump itself is small. Even when the reference signal is reliable, the small tangential components needed for the correction may be obscured by numerical noise. 4.6
Layer determination for two-side signatures
When there is just one hidden layer, we are always in the case where two adjacent cells differ by a neuron in the last (and only) layer. This means that, from different pairs of adjacent cells, we can collect several signatures referring to the same (ℓ−1) (ℓ−1) layer. A naive heuristic is to assume that each neuron η1 , . . . , ηnℓ−1 , has a non-null probability p1 , . . . , pnℓ−1 > 0 to switch and that p = min{p1 , . . . , pnℓ−1 }. The hardest signature to get will be the one corresponding to the neuron that is the least likely to switch. The expectation for the number of attempts before observing the least likely signature is 1/p attempts. In the presence of several hidden layers, we can collect column signatures that do not correspond to the column of the matrix W (ℓ) . In general, there is no way to check if two different cells differ in the last layer or in an intermediate one. With a similar argument to [6] and [3], we expect that signatures corresponding to the columns of the last matrix will have a higher frequency than signatures observed as a consequence of a switch in some intermediate state. To justify this assumption, recall Eq. (11). For two linear cells C0 , C1 that differ only in the (i) activation of the neuron ηj we have ⊺ ∆(C0 , C1 ) = Li,C0 Ej,j Ri,C0 = c∗,j rj,∗ ∈ Rnℓ ×n0 ,
Suppose there exist two other cells C2 , C3 that differ in the activation of the (i) same neuron ηj . If i = 1, then we have R1,C0 = W (1) = R1,C2 , which means (1)
that the matrices ∆(C0 , C1 ) and ∆(C2 , C3 ) will all share the same row wj,∗ up to some multiplicative scalar. That is, we would observe the same signature twice. An identical argument can be made for the columns if the neuron that switches on both pairs C0 , C1 and C2 , C3 lies in the last hidden layer ℓ − 1. In (ℓ−1) particular, let ηj be that neuron, then we have Lℓ−1,C0 = W (ℓ) = Lℓ−1,C2 , which means that the matrices ∆(C0 , C1 ) and ∆(C2 , C3 ) will all share the same (ℓ) column w∗,j up to some scalar. In all the other cases, as there is no reason why the activation patterns relative to the layer i + 1 to ℓ are the same for two different cells, we expect that 20
in general Ri,C0 ̸= Ri,C2 and Li,C0 ̸= Li,C2 . Thus, while we find two independent transitions of the same neuron, we will observe two different column signatures. A challenge for the whole process is that there is no guarantee that the frequency argument will let us attribute the signature to the correct layer. In [6], it was observed that the preimages of critical hyperplanes relative to the nodes in the second and successive layers are hyperplanes bent by different activation patterns in the first layer. This means that the signatures relative to the neurons in the first layer are always the same throughout the measurements, whereas the signatures relative to deeper nodes are “distorted” in different ways. Another interesting observation was made in [10], noticing that, treating the weights of the row signatures as random variables, the variance inside row signatures coming from deeper layers is larger compared to the variance of the attacked layer. 4.7
Limits of column signatures
Working with Jacobians allows us to recover both row and column signatures. A tempting idea is to use column signatures to reconstruct the last layer and go backward. This would allow us to peel the neural network on both sides instead of being limited to working only from the first towards the last layer. Unfortunately, this approach presents two important limitations on the recovery of the biases and the signs. Bias recovery Starting their attack from the first layer, Carlini et al. [6] could determine a critical point x̂ for which wi⊺ x̂+bi = 0. From this equation, it is easy to recover the value of the i-th component of the bias vector b as bi = −wi⊺ x̂. In our case, once we collect all the column-signatures from W (ℓ) , and consider different inputs xi indexed by i ∈ I, we will have the equation W (ℓ) h(ℓ−1) (xi ) + b(ℓ) = F (xi ).
(25)
where h(ℓ−1) (xi ) is the output of the (ℓ−1)-th layer at input xi . We have no information about the vector h(ℓ−1) (xi ) besides that it is the output of a componentwise ReLU, meaning its coordinates are all positive or null. We also do not know exactly W (ℓ) but we know a matrix Ŵ (ℓ) whose columns ŵ1 , . . . , ŵm are the column signatures of the columns of W (ℓ) . The relation between the original weights W (ℓ) and the signature matrix is given by Ŵ (ℓ) = W (ℓ) M where M = P D is the product between a permutation P and a diagonal matrix D. If there is a solution yi = h(ℓ−1) (xi ) for Eq. (25), then ŷi = M −1 yi will be a solution of Ŵ (ℓ) ŷi + b(ℓ) = F (xi ).
(26)
For each of these equations, the weight matrix Ŵ (ℓ) and the outputs F (xi ) are known, but both b(ℓ) and ŷi are unknown. From Eq. (26), it follows that F (xi ) ∈ ColSpan(ŵ1 , . . . , ŵm , b(ℓ) ) for all i; then b(ℓ) ∈ ∩i∈I ColSpan(ŵ1 , . . . , ŵm , F (xi )). When Rank(Ŵ (ℓ) ) < nℓ − 1 the spaces ColSpan(ŵ1 , . . . , ŵm , F (xi )) ⊊ Rnℓ are 21
strictly contained in Rnℓ but, at best, we can determine b(ℓ) only up to an element in ColSpan(ŵ1 , . . . , ŵm ) that is contained in each subspace of the form ColSpan(ŵ1 , . . . , ŵm , F (xi )). Sign recovery Determining the correct sign for the column-signatures is also challenging. A simple sufficient case in which it would be easier is if there is no bias (or if we could assume to know the bias from another source) and that the the matrix W (ℓ) has full column rank Rank(W (ℓ) ) = nℓ . That is, if the equation W (ℓ) y + b(ℓ) = F (x) admits a unique solution. Suppose that we manage to recover and distinguish all the column-signatures of the matrix W (ℓ) and denote by Ŵ (ℓ) the matrix whose columns are those signatures. To determine the correct sign of these columns, we can start by setting them positive (i.e., having the first non-zero component positive). Observe that each component of h(ℓ−1) (x) is the output of a ReLU, then the solution y must have all its components ≥ 0. If some components of the solution y we found are negative, it means we must multiply the corresponding column in Ŵ by −1; in this way, we will force each component of the solution to be positive, while determining the correct sign of each signature. If a coordinate of y is zero because the corresponding neuron is inactive, that observation provides no information about the sign of the associated column, and several inputs may therefore be required. Moreover, if W (ℓ) does not have full column rank, the solution y is not unique and non-negativity alone may not determine the column signs. The scope of our contribution is limited to the signature-estimation primitive, which we evaluate experimentally in the following section, rather than a new endto-end extraction attack. As in prior cryptanalytic extraction work, the vector estimator assumes access to raw multi-dimensional outputs. Although column signatures reveal additional structural information, they do not by themselves enable backward extraction, since ambiguities in bias, orientation, and layer attribution remain. The adaptive refinements additionally rely on heuristic cellconsistency tests whose reliability decreases under low numerical precision.
5
Experiments
We call the network under attack the oracle, a black box we query with an input and read a raw output vector from, with no access to weights. We evaluate the two-sided estimator of Section 4 against Carlini et al.’s row-signature estimator [6], isolating the estimation step from critical-point search. Every method receives the same estimated critical point x∗ and estimates its own output row and coordinate signs from queries alone, never from the oracle’s true weights. 5.1
Setup
We attack single-hidden-layer ReLU networks F (x) = W2 ReLU(W1 x + b1 ) + b2 with 784-dimensional input and hidden widths h ∈ {8, 32}, under three simu22
lated oracle precisions (float64, float32, float16) applied to both the forward pass and every finite-difference step. We run every method against two weight sources: random Gaussian weights and weights after training on MNIST (reaching 92% test accuracy). We report angular error dangle = 1 − |û · u|/(∥û∥∥u∥) ∈ [0, 1] between the recovered and true row direction, the maximum coordinate-wise deviation d∞ between the normalized and sign-aligned parameter vectors, including the bias, and the number of oracle queries that each estimator spends past x∗ . Each query counts as a single oracle input returning the full output vector F (x) ∈ Rnℓ , not one per scalar component. We compare against two variants of Carlini et al.’s attack. Their original construction and an adaptive version of it, and evaluate three variants of our own attack. C-band Carlini et al.’s two-point safety-band construction: for each direction, a fixed offset is added before probing. C-band-adaptive the probing interval is maximized, constrained to remain inside the same cell, in each direction to improve measurement accuracy. J-first-row an ablation of J-vector-fixed that reads w directly off row 0 of the same fixed-step ∆ (Eq. (11)) instead of taking its leading right singular vector, isolating whether SVD’s pooling across rows is actually earning its keep. J-vector-fixed our method reads directly from ∆ (Eq. (11)) at a fixed probing interval δ (Eq. (12)), with no refinement. J-vector-adaptive the “Adaptive steps” refinement from Section 4: the probing interval δ is grown independently per side to the largest interval passing a local consistency check, before reading ∆. J-vector-rotate the “rotating adjustment” refinement from Section 4: a further correction pass that probes along directions perpendicular to the current signature estimate, to recover the residual component of the true neuron signature missed by the initial estimate. 5.2
Results
Table 1 reports the median angular error dangle at each precision, averaged over both widths. We report results separately for random and MNIST-trained weights. We also report the joint direction-and-bias error d∞ , computed similarly to [6]. We first rescale (w, b) and (ŵ, b̂) to the same unit-normal hyperplane form. The bias is divided by the same norm as its corresponding weight vector. We then sign-align the estimates using s = sign(w · ŵ) and compute the largest absolute coordinate-wise difference between the two augmented vectors. Thus, d∞ is an L∞ error, not an average over coordinates. For each precision, the table reports the median d∞ over trials and the median query cost. Values below 10−15 indicate recovery that is exact up to the numerical precision. 23
Finally, we include J-first-row. Instead of extracting w as the leading right singular vector, J-first-row reads it directly from row 0 of the same fixed-step ∆. This isolates the contribution of the SVD step. Table 1. Median angular error dangle , median joint direction-and-bias error d∞ , and median query cost, random vs. trained (MNIST) weights, by simulated oracle precision, with trials from widths {8, 32} pooled before taking a single median (n = 90 random / 90 trained at float64/float32, 44 random / 20 trained at float16). The valid-transition rate, the fraction of attempted clean single-neuron crossings that were usable, was 100%/99% (random/trained) at float64, 98%/97% at float32, and only 2.4%/1.1% at float16. Both weight sources estimate their own row, sign, and starting point from queries alone. Method C-band (Carlini) C-band-adaptive (Carlini) J-first-row (ours) J-vector-fixed (ours) J-vector-adaptive (ours) J-vector-rotate (ours) Method
dangle f64 (rand/trained)
dangle f32 (rand/trained)
dangle f16 (rand/trained)
< 10−15 / < 10−15 < 10−15 / < 10−15 < 10−15 / < 10−15 < 10−15 / < 10−15 < 10−15 / < 10−15 < 10−15 / < 10−15
1.5 × 10−4 / 7.5 × 10−5 5.3 × 10−3 / 7.0 × 10−3 7.1 × 10−4 / 4.0 × 10−5 5.4 × 10−5 / 2.2 × 10−6 5.5 × 10−3 / 1.1 × 10−2 4.3 × 10−3 / 9.0 × 10−3
2.6 × 10−1 / 3.3 × 10−1 8.6 × 10−1 / 7.2 × 10−1 4.1 × 10−1 / 3.4 × 10−1 4.9 × 10−2 / 3.8 × 10−1 2.4 × 10−2 / 7.9 × 10−1 8.2 × 10−3 / 7.6 × 10−1
d∞ f64 (rand/trained)
d∞ f32 (rand/trained)
d∞ f16 (rand/trained)
C-band (Carlini) 1.5 × 10−10 / 1.1 × 10−10 C-band-adaptive (Carlini) 1.6 × 10−10 / 1.0 × 10−10 J-first-row (ours) 1.9 × 10−9 / 5.0 × 10−10 J-vector-fixed (ours) 4.0 × 10−10 / 1.0 × 10−10 J-vector-adaptive (ours) 4.6 × 10−15 / 3.1 × 10−15 J-vector-rotate (ours) 3.7 × 10−15 / 1.9 × 10−15
7.5 × 10−3 / 4.1 × 10−3 4.5 × 10−2 / 5.6 × 10−2 1.5 × 10−2 / 3.9 × 10−3 2.8 × 10−3 / 8.2 × 10−4 4.0 × 10−2 / 4.9 × 10−2 3.8 × 10−2 / 3.9 × 10−2
2.1 × 10−1 / 2.2 × 10−1 2.0 × 10−1 / 2.6 × 10−1 2.0 × 10−1 / 2.6 × 10−1 1.3 × 10−1 / 2.2 × 10−1 6.4 × 10−2 / 2.5 × 10−1 3.5 × 10−2 / 3.0 × 10−1
Method C-band (Carlini) C-band-adaptive (Carlini) J-first-row (ours) J-vector-fixed (ours) J-vector-adaptive (ours) J-vector-rotate (ours)
Queries f64 (rand/trained) Queries f32 (rand/trained) Queries f16 (rand/trained) 6 268 / 6 268 18 804 / 18 804 1 570 / 1 570 1 570 / 1 570 55 348 / 55 445 201 738 / 201 510
6 208 / 6 226 18 772 / 20 428 1 570 / 1 570 1 570 / 1 570 41 918 / 42 462 146 110 / 145 050
5 758 / 5 692 27 290 / 26 651 1 570 / 1 570 1 570 / 1 570 14 588 / 17 412 52 172 / 58 888
Increasing the probing interval benefits J-vector-adaptive and J-vector-rotate only when the local consistency check can reliably distinguish true cell boundaries from numerical fluctuations. In float64, numerical fluctuations are much smaller than neuron activation jumps, allowing the interval to grow to the true boundary without significant error. At lower precisions, however, numerical fluctuations can mask boundary crossings, causing interval growth to overshoot into a neighboring cell and contaminate ∆. J-vector-fixed avoids this issue by using a short, fixed interval. We verified this against oracle ground truth. For every step accepted by (18), we compared the endpoint’s true activation pattern, computed from the target network’s privileged weights, with that of the base point. We sampled 20 coordinate directions, both signs, and both sides (J− , J+ ) per trial, using the population from Table 1. At float64, none of the 14, 400 accepted intervals crossed a cell boundary. At float32, 55% and 44% of accepted intervals crossed 24
a boundary for random and trained weights, respectively. At float16, these rates were 15% and 3.5%. This lower rate at float16 than at float32, is explained by the consistency threshold (below), which is more conservative at float16 and so causes the adaptive search to reject more intervals and terminate earlier, before they can overshoot. Thus, overshooting directly explains the degradation of interval growth at reduced precision. The consistency threshold is adapted to the numerical precision of the oracle. Reducing this threshold makes the test more conservative and can limit the acceptance of intervals that cross a cell boundary, but it can also cause clean intervals to be discarded. Its selection therefore involves a trade-off between accepting contaminated intervals and rejecting valid ones. At float64, all methods achieve essentially exact recovery, with dangle < 10−15 for all six estimators. At float32, J-vector-fixed gives the strongest accuracy on random weights, outperforming C-band-adaptive by nearly two orders of magnitude and C-band by approximately a factor of 3. On trained weights, its advantage over C-band-adaptive exceeds three orders of magnitude, while requiring approximately 4–12× fewer queries. J-first-row further shows the benefit of pooling information across rows. Replacing the single-row estimate with SVD pooling across all rows of ∆ reduces dangle by approximately 13× and 18× for random and trained weights at float32, respectively, at identical query cost. At float16, the improvement is approximately 8× for random weights. The trained-weight case is the only exception, where J-first-row slightly outperforms J-vector-fixed (3.4 × 10−1 versus 3.8 × 10−1 ). At float16, refined estimators are highly sensitive to the weight source. Jvector-rotate outperforms both Carlini variants on random weights, whereas J-vector-adaptive and J-vector-rotate degrade to approximately 0.7–0.8 error on trained weights. J-vector-fixed performs similarly to C-band and C-bandadaptive in both cases. The trained-weight sample is also smaller (n = 20 versus 44), reflecting the lower frequency of clean single-neuron crossings. In query cost, J-vector-fixed consistently requires fewer queries than both Carlini variants. J-vector-adaptive and J-vector-rotate are generally more expensive than C-band and C-band-adaptive, with the exception of float16, where J-vector-adaptive uses fewer queries than C-band-adaptive. Bias recovery. Recovering a neuron means recovering its bias along with its row direction, since together they place the critical hyperplane in input space. Table 2 reports the same comparison for the normalized bias error dbias . Throughout, we recover the normalized hyperplane offset b/∥w∥, not the physical parameter b, consistent with the d∞ metric above. Both Carlini variants recover the bias from a single output coordinate, b̂ = −ŵ⊤ x∗ . Our three methods instead solve the full-vector least-squares problem of Section 4, using every coordinate’s crossing to jointly estimate it. Query cost matches Table 1, since bias is computed directly from the neuron weights and the queries already used to detect the boundary between two cells. 25
Table 2. Median normalized bias error dbias , random vs. trained (MNIST) weights. Method
dbias f64 (rand/trained) dbias f32 (rand/trained) dbias f16 (rand/trained)
C-band (Carlini) 1.5 × 10−10 / 1.1 × 10−10 6.8 × 10−3 / 3.8 × 10−3 2.0 × 10−1 / 2.0 × 10−1 C-band-adaptive (Carlini) 1.6 × 10−10 / 1.0 × 10−10 2.5 × 10−2 / 3.2 × 10−2 1.8 × 10−1 / 2.7 × 10−1 J-vector-fixed (ours) 3.6 × 10−10 / 1.0 × 10−10 2.6 × 10−3 / 8.2 × 10−4 1.3 × 10−1 / 2.2 × 10−1 J-vector-adaptive (ours) 4.2 × 10−15 / 2.6 × 10−15 3.2 × 10−2 / 4.1 × 10−2 5.1 × 10−2 / 2.1 × 10−1 J-vector-rotate (ours) 3.6 × 10−15 / 1.6 × 10−15 3.0 × 10−2 / 3.4 × 10−2 2.6 × 10−2 / 2.2 × 10−1
J-vector-adaptive and J-vector-rotate recover the bias much more accurately than the Carlini methods at float64. Their error is around 10−15 , compared with about 10−10 for the Carlini methods. This improvement comes from estimating the boundary offset using all available output coordinates rather than relying on a single coordinate. At float32 and float16, however, this advantage shrinks, since the directions (the neuron weights) computed by the two methods also degrade in precision.
6
Conclusion
When a single neuron changes activation across two adjacent linear cells, the resulting Jacobian difference has rank one. Its right factor contains the rowside information exploited by previous extraction attacks, while its left factor provides complementary column-side information. By applying SVD to the full Jacobian difference, the information from all output components can be combined to obtain more accurate signature estimates. Our experiments show that this improves signature recovery under finite precision, particularly in the float32 regime, without increasing the number of oracle queries. We also considered adaptive refinements based on enlarging the intervals used to estimate the Jacobians. Their main limitation in float32 and float16 was the reliability of the cell-membership test: accepting points outside the intended cell contaminates the estimated Jacobian difference. A natural direction for future work is therefore to develop more reliable membership tests under reduced numerical precision. The column-side information also exposes a new leakage channel; however, using it to extract layers beyond the last and combining it with forward extraction into a complete bidirectional procedure remain open problems.
References 1. Batina, L., Bhasin, S., Jap, D., Picek, S.: CSI NN: Reverse engineering of neural network architectures through electromagnetic side channel. In: 28th USENIX Security Symposium (USENIX Security 19). pp. 515–532. USENIX Association, Santa Clara, CA (Aug 2019), https://www.usenix.org/conference/ usenixsecurity19/presentation/batina
26
2. Canales-Martı́nez, I.A., Chávez-Saab, J., Hambitzer, A., Rodrı́guez-Henrı́quez, F., Satpute, N., Shamir, A.: Polynomial time cryptanalytic extraction of neural network models. In: Joye, M., Leander, G. (eds.) Advances in Cryptology – EUROCRYPT 2024. pp. 3–33. Springer Nature Switzerland, Cham (2024) 3. Canales-Martı́nez, I.A., Chavez-Saab, J., Hambitzer, A., Rodrı́guez-Henrı́quez, F., Satpute, N., Shamir, A.: Polynomial time cryptanalytic extraction of neural network models. Cryptology ePrint Archive, Paper 2023/1526 (2023), https: //eprint.iacr.org/2023/1526 4. Carlini, N., Chávez-Saab, J., Hambitzer, A., Rodrı́guez-Henrı́quez, F., Shamir, A.: Polynomial time cryptanalytic extraction of deep neural networks in the hard-label setting. In: Fehr, S., Fouque, P.A. (eds.) Advances in Cryptology – EUROCRYPT 2025. pp. 364–396. Springer Nature Switzerland, Cham (2025) 5. Carlini, N., Chávez-Saab, J., Hambitzer, A., Rodrı́guez-Henrı́quez, F., Shamir, A.: Polynomial time cryptanalytic extraction of deep neural networks in the hard-label setting. Cryptology ePrint Archive, Paper 2024/1580 (2024), https://eprint. iacr.org/2024/1580 6. Carlini, N., Jagielski, M., Mironov, I.: Cryptanalytic extraction of neural network models (2020), https://arxiv.org/abs/2003.04884 7. Chen, Y., Dong, X., Guo, J., Shen, Y., Wang, A., Wang, X.: Hard-label cryptanalytic extraction of neural network models. In: International Conference on the Theory and Application of Cryptology and Information Security. pp. 207–236. Springer (2024) 8. Fefferman, C., et al.: Reconstructing a neural net from its output. Revista Matemática Iberoamericana 10(3), 507–556 (1994) 9. Jagielski, M., Carlini, N., Berthelot, D., Kurakin, A., Papernot, N.: High accuracy and high fidelity extraction of neural networks. In: 29th USENIX security symposium (USENIX Security 20). pp. 1345–1362 (2020) 10. Liu, H., Siproudhis, A., Experton, S., Lorenz, P., Boura, C., Peyrin, T.: Navigating the deep: End-to-end extraction on deep neural networks. In: Daemen, J., Thomé, E. (eds.) Advances in Cryptology – EUROCRYPT 2026. pp. 482–512. Springer Nature Switzerland, Cham (2026) 11. Lowd, D., Meek, C.: Adversarial learning. In: Proceedings of the eleventh ACM SIGKDD international conference on Knowledge discovery in data mining. pp. 641–647 (2005) 12. Milli, S., Schmidt, L., Dragan, A.D., Hardt, M.: Model reconstruction from model explanations. In: Proceedings of the Conference on Fairness, Accountability, and Transparency. p. 1–9. FAT* ’19, Association for Computing Machinery, New York, NY, USA (2019). https://doi.org/10.1145/3287560.3287562, https://doi.org/10.1145/3287560.3287562 13. Noman, I., Vukelic, M., Jasuja, V.: Output rounding is not a free defense against cryptanalytic neural network extraction. 6.5610 final project report, Massachusetts Institute of Technology (2026), https://65610.csail.mit.edu/2026/reports/ cryptanalytic_nn_extract.pdf, spring 2026 14. Tramer, F., Zhang, F., Juels, A., Reiter, M.K., Ristenpart, T.: Stealing machine learning models via prediction apis. In: 25th USENIX security symposium (USENIX Security 16). pp. 601–618 (2016)
27