On the complexity of nonsmooth automatic differentiation
Authors: Jerome Bolte
Ryan Boustany
Edouard Pauwels
Béatrice Pesquet-Popescu
Venue: ICLR 2023 notable top 25% (2023)
Keywords: Automatic differentiation
nonsmooth derivatives
computational complexity
cheap derivatives
conservative gradients

Abstract: Using the notion of conservative gradient, we provide a simple model to estimate the computational costs of the backward and forward modes of algorithmic differentiation for a wide class of nonsmooth programs. The complexity overhead of the backward mode turns out to be independent of the dimension when using programs with locally Lipschitz semi-algebraic or definable elementary functions. This extends considerably the Baur-Strassen's smooth cheap gradient principle. We illustrate our results by establishing fast backpropagation results of conservative gradients through feedforward neural networks with standard activation and  loss functions. Nonsmooth backpropagation's cheapness contrasts with concurrent forward approaches, which have, to this day, dimensional-dependent worst case  overhead estimates. We provide further results suggesting the superiority of backward propagation of conservative gradients. Indeed, we relate the complexity of computing a large number of directional derivatives to that of matrix multiplication, and we show that finding two subgradients in the Clarke subdifferential of a function is a NP-hard problem.

Published as a conference paper at ICLR 2023
ON THE COMPLEXITY OF NONSMOOTH AUTOMATIC DIF- FERENTIATION
Jérôme Bolte!? Ryan Boustany Edouard Pauwels 2, & Béatrice Pesquet-Popeseu
Toulouse School of Economics
2 Université de Toulouse
IRIT, CNRS. Institut Universitaire de France (IUF).
Thales LAS France
{je rome bol te, ryan. "yan.bouste
fr, peoourdd.paueer sereelssrrii siirit. fr,
beatrice, pesuuettthaaesggrpp.co
ABSTRACT
Using the notion of conservative gradient, we provide a simple model to estimate the computational costs of the backward and forward modes of algorithmic dif- ferentiation for wide class of nonsmooth programs. The overhead complexity the backward mode turns out to be independent of the dimension when using programs with locally Lipschitz semi-algebraic or definable elementary functions. This considerably extends Stur-Strassen's 's smooth cheap gradient principle. We illustrate results by establishing fast backpropagation results of conservative gradients through feedforward neural networks with standard activation and loss functions. Nonsmooth backpropagation's cheapness contrasts with concurrent forward approaches, which have, to this day, dimensional -eepeddenn worst-case overhead estimates. We provide further results suggesting the superiority of back- ward propagation of conservative gradients. Indeed, we relate the complexity computing large number of directional derivatives to that of matrix multiplica- tion, and we show that finding two subgradients in the Clarke subdifferential of function an NP-hard problem.
INTRODUCTION
Automatic evaluation of derivatives: Algorithmic differentiation (AD) appeared around 60 years ago (Beda al. (1959); Wengert (1964)), and has been since then constantly developed and used in many contexts, see Griewank et al. (1989); Griewank and Walther (2008) for thorough discussion. Today, the core of modern learning architectures (Rumelhart et al., 1986; LeCun al., 2015; Baydin et al., 2018), to the point that training neural network (NN) is ultimately way combine the outputs of AD. There are many practical and thcoretical developments available nowadays: flexible and efficient numerical libraries (Abadi et al., 2016; Paszke al., 2019; Bradbury al., 2018), an implicit differentiation theory (Griewank and Faure, 2003; Griewank and Walther, 2008) and its extensions (Agrawal et al., 2019; Bai et al., 2019; Bolte al., 2021: Blondel al., 2021), the adjoint method (Farrell et al., 2013; Pearlmutter, 1995; Plessix, 2006) with application to neural style differentiation of optimization algorithms (Griewank and Faure, 2003; Mehmood and Ochs, 2020; Bertrand et al., 2020; Lorraine et al., 2020), or differentiation of conjugate gradient algorithms (Gratton et al., 2014).
ODEs (Chen al., 2018), "piggyback
Backward algorithmic differentiation, or backpropagation, plays particular role when smooth optimization tasks are at stake, as it evaluates the gradient of function with cost proportional to that of function evaluations, independently of dimension. This property, called the cheap gradient principle (Wolfe, 1982; Griewank and Walther, 2008), at the root the machine learning libraries revolution. According the key complexity theory version of this result due to Baur and Strassen (1983), arithmetic complexity of the evaluation the derivative times the complexity of function evaluation. Extensions exist for smooth differentiable functions Baur and Strassen (1983); Griewank and Walther (2008) but standard computational practice of AD consists of little known about the nonsmooth case.
rational function at most

---

Published as a conference paper at ICLR 2023
The objective of this paper is precisely to present simple, general, nonsmooth cheap conservative principle and to explore other complexity results for evaluating nonsmooth derivatives. This extends the cheap gradient principle of smooth AD to the path differentiablee world Bolte and Pauwels (2020b) which includes semi-algebraic and more generally definable functions Coste (2000a;b), a class that contains the vast majority of machine learning programs used in practice, see for example Bolte and Pauwels (2020b).
Nonsmooth AD & computational complexity: Sorting values, pooling data, thresholding func- tions, or determining closest points are some of the most essential numerical decision operations. They are ubiquitous in machine learning and modern optimization. All of them are nonsmooth, and most of them have very desirable feature: they cheap compute, much cheaper than smoothed surrogates. For instance, the famous ReLU activation in deep learning, whose role is to threshold to zero negative values allow for the inactivity of neurons, requires only one bit of encoding in theory. On the other hand, other nonlinear activations potentially require auxiliary algorithms for their evaluation, incurring a higher computational cost. This simplicity of use also comes with the issue of finding an adequate way of training models and, thus differentiating objects.
The standard computational practice of AD consists in applying differential calculus rules directly to nonsmooth objects, replacing gradients by surrogates, typically Clarke subgradients. This is how AD is performed within TensorFlow, PyTorch or Jax. This approach has shown tremendous success (LeCun et al., 2015) and has been massively used for the last 10 years. Yet, despite this empirical success, Barton al. claimed in Barton et al. (2018) that "there does not seem to exist /at this day/ true analogous reverse AD mode to compute generulized derivatives for nonsmooth functions", illustrating the difficulty of nonsmooth AD. Conservative gradients were introduced as a faithful mathematical model capturing the formal application of calculus rules to subdifferentials by Bolte and Pauwels (2020a;b); Bolte et al. (2021). The author unfamiliar with this notion may reduce, in a MI. context, conservative gradients and Jacobians. Our goal is to provide an adequate computational complexity theory for conservative calculus, theory that will therefore match standard practical approaches.
outputs of calculus rules formally applied Clarke subgradients
Among other possible first-order options offered by nonsmooth calculus, we also investigate the properties of directional derivatives and those of the Clarke subdifferential. For directional derivatives, our motivation comes from the fact that this nonsmooth operation has general calculus rules, while the Clarke subdifferential is central in terms of variational interpretation.
The main thesis of this work is that conservative gradients have computational properties similar to smooth derivatives, which are much more favorable than those of alternative nonsmooth oracles such subgradients or directional derivatives.
Contributions:
We provide a simple computational model for addressing the question of complexity theory of nonsmooth numerical programs.
For the backward mode, we prove a cheap conservative gradient principle à la Baur-Strassed generalizing state of the art nonsmooth programs modeling most NNs. We establish that, regardless of dimension, the computational cost of conservative gradient is of the order of that of function evaluation. Our results provide theoretical validation of the fact that the cost of backpropagation does not depend on the programs' smoothness.
For the forward mode, we relate the computational cost of p directional derivatives to that of p x p matrix multiplication. We provide lower complexity bounds that illustrate the limits to which this deficiency may be improved. This applies to existing nonsmooth AD frameworks (Khan and Barton, 2012: 2013).
We establish that computing two distinct elements in the Clarke subdifferential of given point is NP-hard for simple ReLU programs. This result also applies to the lexicographic subdifferential. In contrast, we show that the problem can be solved in polynomial time for conservative gradients. This reflects the computational difficulty of dealing with the Clarke subdifferential.
A result of independent interest: deciding differentiability of a ReLU program at a point is Na-hard.
Relation with existing work: Conservative gradients were introduced in Bolte and Pauwels (2020a;b) to model "formal subdifferentiation" used by practitioners and nonsmooth "backpropaga- tion". They were further studied in Lewis and Tian (2021); Davis and Drusvyatskiy (2021); Bolte al. (2021) and empirically investigated in Bertoin al. (2021). Computational complexity was
2

---

Published as a conference paper at ICLR 2023
only qualitatively considered. We provide rigorous description of this aspect based arithmetic computational cost framework capturing programming with nondifferentiable components. The quest for computationally cheap nonsmooth derivative has long history in AD literature. Existing works of Griewank (Griewank and Walther, 2008; Griewank, 2013; Griewank and Rojas, 2019; Griewank and Walther, 2020) essentially based piecewise smoothness structures (Scholtes, 2012). A cheap subgradient principle was also given in Kakade and Lee (2018), but it requires a very strong qualification condition. As illustrated Griewank and Rojas (2019), such qualification conditions can be computationally hard to check in practice.
In another research line, based on chain rules for directional derivatives, Khan-Barre (Khan and Barton, 2012; 2013; 2015; Barton et al., 2018) studied the vector forward mode AD. In particular, they investigated the forward AD framework to evaluate elements the lexicographic subdifferential (see Nesterov (2005)), which is contained in the Clarke subdifferential. In the worst case, the computational overhead ratio they obtain is proportional to the ambient dimension. This contrasts with our cheap gradient principle, whose constant is dimension-less. While these contributions are most relevant to nonsmooth AD, their applicability to large-scale learning models is limited, due to the central role of forward AD.
Organization of the paper: We introduce elements of nonsmooth analysis and, in particular, the notion of conservative gradient used throughout this work in Section 2. Section 3 describes general model of computation that allows one express the computational cost and complexity of programs, functions and their conservative gradients. This section also presents an abstract program algorithmic differentiation framework. These elements are gathered in Section which presents our extension of the Baur-Strassen result with the cheap conservative gradient principle and its illustrations. To conclude, Section 5, we describe computational lower bounds for evaluating directional derivatives and distinct subgradients for simple programs.
NONSMOOTH GENERALIZED GRADIENTS
They are fundamental to expressing variations of nonsmooth losses in Machine Learning. Given locally Lipschitz continuous function F TRP TR. the Clarke subdifferential of F is
a F=F(z) = conv
lim F((zz) E diff F, 1-++0000 k
(I)
where d the full measure where F is differentiable and V the standard gradient (Clarke, 1983). The subdifferential set-valued, which we write 70 F: R2 = RR. For each x E RP, elements of 2= F(x) are called Clarke subgradients of F. selection d in F, is function d: I2 RR such that for all I IRP, d(x)) Fz). If is C² then F [VF) everywhere so the only possible F.. We will manipulate derived dictionaries, which typically provide selection in either the Clarke subdifferential, or more general sel-velued maps.
selection
Example 1 For ReLU: t + max(0, t), we have (eReLU(t) is (0) ift < 0, 0, 0, ift 0 and [o, if O. We may define the function ReL.U' selection ReeLU
ReLU'(() = 1, if 0, ReLU'(() 0, otherwise.
The chain-rule, essential to AD, generally fails for Clarke subgradients. This is why we now consider the more flexible notion conservative gradients.
Definition (Conservative gradient) Let F: R R be a locally Lipschitz continuous function and Dr: R Rp locally bounded, nonempty and graph closed set-valued map. Then Dr is conservative gradient for F, if for any absolutely continuous curve y: [o, 1] IR,,
dt F(((t)) = (w, (())
Yoe Dr(r(t)),
for almost all [0.
In this case, F is called path differentiable. Conservative Jacobians are defined similarly. As in Section 2, d: R RP is selection of Dr if d(x)) Dr(x) for all TRP

---

Published as a conference paper at ICLR 2023
A rich class of path differentiable functions is given by locally Lipschitz continuous semi semi-algebraic functions with the Clarke subdifferential as conservative gradient. Actually, virtually all functions used in machine learning are path differentiable (Bolte and Pauwels, 2020a:b). The most salient facts about path differentiable functions and their conservative gradients are:
(Clarke subgradient), for all x: R2, F()) conv conv(Dr(zz)).
GGradient almost everywhere) Conservative gradients are gradients a.e (Bolte and Pauwels, 2020a).
(First-order oracle) Selection in conservative gradients can be used as surrogate gradients, while
preserving convergence guaranties (Bolte and Pauwels, 2020a;b; Bolte et al., 2021).
Conservative Jacobians can be composed while preserving conservativity (Bolte and Pauwels, 2020a), a feature which do not enjoy Clarke Jacobians: let F: E" - R" G: R - R' be locally Lipschitz continuous mappings, dF: RP I xx and de: RT R be selections in conservative Jacobians for F and G respectively. Then the product mapping x: - da (xx)) x dr((x) is a selection in conservative Jacobian for G F. The use of conservative Jacobians provides very convenient framework to model AD in the nonsmooth case, see Bolte and Pauwels (2020a:b).
fundamental theorem is the following:
Theorem (Path differentiable functions ubiquitous) (Bolte and Pauwels (2020a)) Locally Lipchitz semialgebraic (or definable) functions are path differentiable.
3 PROGRAMS, COMPLEXITY AND AUTOMATIC DIFFERENTIATION
3.1 CALCULUS MODEL, PROGRAMS, COMPUTATIONAL COST AND COMPLEXITY
dictionary is a finite set of real functions (e.g- (t,,,x,))) it is paired with (D), a set of elementary programs implementing them in real arithmetic. Starting from (D) we aim at capturing the notion of "program of programs" at any depth. As this is an inductive process, we call k E N program "level", which is simply an induction counter needed for consistency. Recursively, programs of level k + 1, in (D)+ consist of combinations of outputs of programs of level k, in (D)).. For example if P1 and Pg are elementary programs in (D))), then the program which sums the outputs of P1 and Po is of level 1. More precisely:
Algorithm 1:
Let p, q be input and output sizes respectively and m p + q a memory size. A predecessor relation a set valued map m) mm) such that for i = I...., for pr)()) < i.
Program data:
pr: {1, (I..... m] =
(p. q. m. pr, (g:)
(r( is empty if & p and nonempty otherwise.
Input: x (r), Tp)
An adapted program sequence (g:) in (D), is a set of programs such that Si has [pr(1)] input arguments and single output, for all i Given (p, q, m, pr, (+))=p++)) 1)). the program given in Algo- rithm 1 is level k program on D.
I: for i p 1 p + 2, do
2: I2 = g (r)()) where
3: Tpr(1) = (z) jepr(i)
..... m.
4: end for
Return: y := (-
The set of programs with dictionary D is P(D) = Ukao P*(D). We shall see however that Pk (D) pl (D) for all k, using modification of the computational graph.
A cost on dictionary D is nonnegative function on D, it extends additively by induction on programs on D through the rule cost(P) = >>= cost(g;) where P is a program on D as described in Algorithm 1. A direct example is the dictionary of arithmetic functions (+, x,/). J. together with addition or multiplication by fixed constants, denoted by + and x respectively' see also Section A.1. Throughout the paper, we assume that dictionaries contain at least operations + and x.
Each program on D may be represented by program in Pl (D) with the same cost, by expanding all subprograms until they reduce to an elementary program. Cost evaluation is thus well defined on such programs. As detailed in Appendix 1, this model of computation is equivalently expressed using directed acyclic graphs
Constants need to he distinguished from variables (for instance define polynomial)
4

---

Published as a conference paper at ICLR 2023
To sum up, have defined the set of programs P(D) on D, which includes programs of programs. The programs 9i in Algorithm may be taken in P(D). The cost of a program evaluated through the calls it makes to elementary programs in the dictionary.
Programs functions: A program P defines unique input-output function f: we say that P "computes" f, or "implements" J. and with a slight abuse of notation, we will identify and when there is no ambiguity (e.g. derivative of P). We use the equivalence relation to relate programs computing the same function. The equivalence classes correspond to functions expressible by programs with a given dictionary D. Given function f: IR. - R2 and program P on dictionary D, with p inputs and outputs, we write f = [P] to denote the fact that P is in the equivalence class of programs computing f. that is, implements f.
Complexity of function: The complexity of function f over dictionary D is the quan- tity comp(f, D) = inf (P)st((), s.t PE (PD), f = [[P]. the infimum being over all programs implementing f dictionary D. It could be infinite, if it is finite then it is attained.
3.2 AUTOMATIC DIFFERENTIATION
We pertain programs implementing functions, that is Algorithm 1 with single outputs q 1.
Given a dictionary D of locally Lipschitz path differentiable functions, derived dictionary is a set of functions D' = 'D which extends D and contains operations required to express at least an element in a conservative gradient for each of the functions in DD for example, an element in the Clarke subdifferential. We also consider a cost function on D', which we denote by cost and which extends to programs over D'. Given programs 9: on D. i = p 1,..., m. we define d; a derived program on D', with [pr(i)] inputs and outputs. which returns an element of conservative gradient for Ji (as for instance Clarke subgradient, or simply a gradient in the C' case). By gd,, we denote program on D' evaluating (z), d, (x)) jointly for a given æ. We denote by Algorithm I', an extension of Algorithm 1 which additionally returns w; dz((zz(1) for p + 1,... ...... by replacing line 2 in Algorithm 1 with call to gd; instead of 9:- The backward (resp. forward) AD program backprop( P) (resp. forprop P)) is defined as follows:
Algorithm 2: Algorithmic differentiation of as in Section 3.1
Input: variables (x:) Forward evaluation with derivatives: evaluate w; = de(Ipz(()) = p + m, with Algorithm I': Algorithm with gd, instead of 9; line 2.
1: Forward mode:
1: Backward mode:
2: Initialize: er e, I,...,P. P.
2: Initialize: v 6m
from canonical basis in RP.
3: for t m....p do
for = p i,...m .....m. do
4: for E pr(t) do
4:
5: Update coordinate j of
das Σ dx
dx
+ :- m[>] + w w.[i
(f)
6: end for
where x ".l,,..,,pp
5: end for
7: end for
Return: cx
drm. and eventually Tmm
Return: (x(s));1 and eventually Tm.
Note that Algorithm 2 starts with Algorithm i... Algorithm with gd, instead of gi on line Its computational cost, denoted cost(gd:), should be thought of as an exogenous parameter: it may model, for instance, the use of underlying software libraries or the hardware properties.
4 COMPUTATIONAL COMPLEXITY OF NONSMOOTH AD
We now evaluate the complexity of the forprop and backprop operations for conservative gradients in the path-differentiable case which encompasses, mentioned earlier, all semi-algebraic and
5

---

Published as a conference paper at ICLR 2023
definable locally Lipschitz functions. We show, in particular, that backpropagation with conservative gradients has a computational overhead ratio that independent of the dimension. This is in contrast with the best known algorithmic oracles for the Clarke subdifferential (see Khan and Barton (2012; 2013; 2015); Barton et al. (2018) and Appendix A.2), whose computational overhead ratio scales linearly with the dimension.
Theorem (Complexity of nonsmooth AD) Let be a program over dictionary 'D of path- differentiable functions with p inputs as in Algorithm 1 & 2. Then, the corresponding function [P] is path differentiable, there is a conservative gradient Dr for the function [P] such that:
(i) (Cost of backward mode) At each input point E IR,, the output of program backprop( P) is in (2(z) and we have cost cost(backprop(P P)) & wh cost(P); where max (cd))( + 2 max (+)s cost( (ost(x))lpr(i)) / cost (9:))
Wb
(3)
1-p+ 1
(ii) (Cost of forward mode) Ar each imput point x RP, the output of program forprop(P) is in D Dr(z) and we have cost (os((ooprro(()) & wy cost(P) where
W
max {(cost (gd,) + pipr(1) cost( x p((pr()) ()cost(+))) cost (g:)).
i-ptl,m
There a dissymmetry between the two modes since the constant wp is independent of the dimension p. This is why property (i) is sometimes called the "cheap conservative gradient principle" extending the classical smooth one which was derived by Baur and Strassen (1983) for real rational functions. Theorem 2 describes worst case upper bounds (maximum over i), which are tight, for example if pr()), costs of gi and gd, are independent of i.
We will consider several examples now.
The class of ReLU programs Let DRELU be the dictionary composed elementary arithmetic operations, logarithm, exponential and the ReLU function:
(4) ReLU program iss a program with dictionary DRLLU; can be expressed a compositional form (Section 3.1) with program sequences in DRELU- Note that this yields path differentiable functions.
DReLU := {+, x, tc, xc, inv, exp, log, ReLU).
Assumption (Computational Cost) In Algorithms (2), define the dictionary D' D'ELLU := DReLU' [ReLU) as in Example 1; then, all operations from D' D'RLLU have unit cost (see Remark 1).
Corollary (Backprop complexity of ReLU programs) Let P be ReLU program, under As- sumption 1, we have: cost (ootbbcckrop(p))) & 5cost(P), ()). This extends to more complex cost weighting schemes (Remark 1) and to selection functions which virtually capture all losses in ML (Remark 2).
Table 1:
Complexity constant of wh Theorem 2 for elementary g in DRALU and derived program with dictionary PRELU- This proves Corollary (more details in appendix B.
(+, x)
(+c, xc)
log
exp
inv
ReLU 3
(cost(ggd= 2cost x prl) cost (g)
3
4
5
Remark 1 (On refined cost systems) Unit cost in Assumption 1 gives simple interpretation to Corollary the cost program is the total number of numerical operations. This rough estimate of computational complexity, could be refined with different weighting schemes. However, the obtained constant 5 is robust to many different weighting choices, far beyond Assumption We detail an example in the Appendix B.2 for which the cost of all smooth nonlinear operations different from + or x Cuoulin > and we model the cost of sign branching in computation of ReLU and ReLU' with constant CRALUU 0. This yields the same constant as Corollary 1.
Remark (Beyond ReLU programs) Many other dictionaries could be considered. ReLU is an example chosen for its simplicity, but Corollary would hold similarly (with the same constant 5) for many different nonsmooth activations or components such as absolute value, max-pooling ELU function, & and lz norms. Similar results could be developed for the class of selection functions, which encompasses the vast majority of ML building blocks (see Bolte and Pauwels (2020b)). This is sketched in Appendix B.3.

---

Published as a conference paper at ICLR 2023
Chaining backpropagation derived programs Our approach is flexible enough to describe "programs of programs" and backpropagation chaining. Let P be a program as in Algorithm 1, with adapted ReLU program sequence {(g, ((x))=+++>> If cost(g:) >> |pr(i)], 9: is "long program", with many operations per input.We We may set gd, backprop(g.) using Algorithm 2, i p 1, ...... m. From Corollary I, we have cost(gd,) /cottt & 5, and for long programs wp = in Theorem 2. This illustrates the versatility of our approach it captures the complexity of chaining backprop operations, the resulting estimate being quite sharp in the regime of long programs.
Beyond backpropagation Programs may be differentiated by other means than backpropagation. Examples include, forward propagation, with applications in optimization and algorithmic unrolling (Mehmood and Ochs, 2020; Lorraine al., 2020; Maclaurin et al., 2015), implicit differentiation Agrawal et al. (2018); Winston and Kolter (2020); Bai et al. (2019); Bolte et al. (2021) with application in optimization and hyperparameter optimization (Bertrand et al., 2020), adjoint differentiation (Plessix. 2006) in programs with components involving ordinary differential equations (Courtier and Rabier, 1997: Chen et al., 2018), differentiation of conjugate gradient (Gratton et al., 2014), Cholesky algorithm (Smith, 1995). approximation of Jacobian matrices involving non-uniform FFT (Wang and Fessler, 2021).
Let P be program in Algorithm 1. Theorem 2 relates the complexity of combining derived programs in Algorithm 2 to the following quantities, for i = + 1...., m:
cost(gdd)/coss((g) ) the "computational overhead ratio". (xr(i))css(()))oool(
the ratio between multiplication cost and average cost per input argument.
The first quantity depends the technique used obtain gd,- The second quantity is typically less than 2 (at least one arithmetic operation per input) and becomes negligible for long programs (many operations per input).
For example in Mehmood and Ochs (2020); Lorraine et al. (2020); Maclaurin et al. (2015), for one i, the program gd, is optimization algorithm in RP, a long program differentiated using forward propagation. The corresponding overhead ratio is in this case 3p + 5 (Theorem 2). If combined with an outer backward pass, we obtain a cennensionddependeee backward differentiation. Our model provides computational cost estimates for mixed techniques, here combination of inner forward and outer backward propagation.
overhead ratio, in contrast with full
ON THE COMPUTATIONAL HARDNESS OF GENERALIZED GRADIENTS
Let P and DP be two programs such that DP evaluates jointly P and derivative of P. In the sequel, we use the term (computational) overhead ratio of DP to denote the quantity costt(D))) and cust(P)
costtPP) As established in Theorem 2, this ratio is dimensionless in the case of backpropagation with conservative gradients. Are there other ways to compute cheap nonsmooth gradients? Toward an answer to this question, we discuss this ratio for other nonsmooth differentiation oracles: directional derivatives (for which relate worst-case complexity to that of matrix multiplication), lexicographic derivatives with forward AD (with overhead ratio of order p Barton et al. (2018)). As for the Clarke subdifferential, prove the hardness of subgradients enumeration. Our motivation to estimate the complexity of these particular types of derivatives (directional, lexicographic and Clarke) is that they serve as a basis to alternative impiementable AD approaches (see Barton et al. (2018) and references therein), and are thus concurrent strategies of conservative gradient backpropagation. The results presented below do not provide a definitive answer, but they strongly suggest that backpropagation of conservative gradients has a much more favorable complexity.
computational overhead ratio of derivatives of P to denote the quantity comp(DP)
5.1 THE OVERHEAD RATIO FOR EVALUATING P DIRECTIONAL DERIVATIVES
Given G: RP - IR locally Lipschitz and x, d RP, the directional derivative of G at in direction d is given by jim, + td)- - G(z))/t when the limit exists. This section considers a family of functions with p inputs and q real parameters, represented by a locally Lipschitz function F: R x R3 + R. for which we investigate hardness of evaluation of p directional derivatives. The function F may describe, for instance, ReLU feedforward neural network empirical loss, parameterized by
7

---

Published as a conference paper at ICLR 2023
q real weights, with p inputs. For functions represented by ReLU programs, we prove overhead ratio of order 20-22-(11) where w is the matrix multiplication exponent (see definition below). In all rigor, it is not known whether w > or w = 2, so the derived ratio could be essentially dimensionless (if w = 2), though all practical evidences are against this far. The best known lower bound is w < 2.37, and in practice, the matrix multiplication exponent is closer to 2.7, both corresponding to dimension-depent dependent overhead, in contrast with the smooth case with essentially dimensionless overhead ratio to evaluate p directional derivatives (essentially a gradient).
Complexity of matrix multiplication: Throughout this section, we set D = (+, x, +c, xc), with unit costs (corresponding to polynomial functions). Denote by c(p) complexity of p x p matrix multiplication. More precisely, if f: Rxx x Px - "RE is such that f(A, B) AB for all, square matrices A, B E R we have c(p) comp(f, D), which we may write c(p) = pl where w is called the matrix multiplication exponent. Note that c(p) > p², as one needs at least one operation for each of the 2p² entries.
Directional derivatives: Given a function F: IP x R R. we denote by Ff: X RPxP I' the function which associates to ø E RR, y E R' and a matrix A E RPxp the directional derivatives with respect to variable, for fixed y, in directions given by the columns of A. The proof of the following theorem is given in Section C.
Theorem (Computational ratio for directional derívatives) There exists a function F: R. x RS - R and a program Pr implementing F on dictionary [+, x, ReLU,+c, xc} (all operations have unit cost), such that for any program P' implementing (y, A) + F1 y, A) on derived dictionary {+, x, ReLU, ReLU', +c, xc).
cost( (Pr) > (dp) - 5p)/(44p))))) -p 2
Theorem 3 has & parameters, parametric dependency is required to express hardness. Indeed, for some parameter values, computation may be trivial (e.g. null values). Alternatively, states that for some values of the q parameters, computing directional derivatives has cost as in (5).
The bound in (5) is sharp up to multiplicative constants for linear ReLU networks, Remark 5 in Appendix A.2.
Our overhead estimate is roughly p--2, it constitutes a bottleneck: "cheap nonsmooth p directional derivatives principle", would imply easy matrix multiplication, to the point that 2. Since the seminal work of Strassen et al. (1969), it is known that w & log2(7) = 2.81. Determining the precise exponent w has been an object of intense research Robinson (2005). Asymptotically, one has & w < 2.373, see Williams (2012); Le Gall (2014), the best known bound being given in Alman and Williams (2021). In this case, the estimate in (5) is roughly p.373
Consequences:
These estimates may involve non-constructive existence proofs, or suffer from the curse of recursion: meaningful efficiency occurs only for inaccessible sizes. According to Dumas and Pan (2016), for values p & 10th the most efficient practical algorithms have a complexity of the order 2 resulting in an overhead of order 7.77 in contrast with the constant overhead incurred by nonsmooth backpropagation. More discussion is given in Appendix A.2.
Comparison with the smooth case: If F C evaluating p directional derivatives is comparatively casier because F"(x, d) = (VF(x), d) for all x, d e IR?. Hence, one may first evaluate V (once), cost similar that of F (cheap gradient principle), and then evaluate p scalar products, at cost p². If the cost of F is of order p² at least (for example F is feedforward neural network with inputs and layer of p hidden neurons), then this is overall proportional to the cost of computing F.
5.2 COMPUTING CLARKE SUBGRADIENTS USING FORWARD AUTOMATIC DIFFERENTIATION
In Khan and Barton (2012; 2013; 2015), several automatic differentiation strategies proposed to evaluate elements of the Clarke subdifferential. These approaches are based on directional (Shapiro, 1990) and lexicographic derivatives (Nesterov, 2005) which satisfy chain rule under structural assumptions. The chain rule may be implemented using the vector forward mode of automatic differentiation (Barton al., 2018), which suffers from computational overhead scaling linearly
8

---

Published as a conference paper at ICLR 2023
in p, contrary to the reverse mode in Theorem 2. Reducing this factor is open question. for compositional functions involving only univariate nonsmoothness such as absolute value (Khan, 2018). More details given in Appendix A.2.1.
5.3 COMPUTATIONAL HARDNESS SUBGRADIENT ENUMERATION
We investigate in this section the hardness finding subgradients for programs defined on the elementary dictionary Do- {+, ReLUU with unit costs. Let denote by P(Do) the set of such programs. We will, with slight abuse of notation. identify program P E Do = [+, ReLU] with the function it computes to state our complexity result (proof in Section D).
Theorem 4 (Clarke subgradients and NP- NP-Hardness)
(i) The problem of finding two distinct subgradients in the Clarke subdifferential of P E P(Do) at given imput (or one single subgradient it is reduced to singleton) is NP- (ii) Deciding if P P(Do) is not differentiable at some given input is NP-hard.
Remark In Theorem 4, numerical parameters and inputs are constrained to be {-1, O, 1), so that the hardness result does not depend on numerical representation and only involves program size (strong NP-hardness). See Appendix for more details.
The above problems (i)-(ii) enter the field of computational complexity as we consider programs P E P(Do) with natural notion of size, given by their cost, cost( P), the number of operations (recall that we assumed unit costs). Since the considered programs implement piecewise linear functions, it follows from (Barton et al., 2018, Proposition 2.7) that, our hardness result also holds for the lexicographic subdifferential Nesterov (2005), which reduces in this case to the set of neighboring gradients (see Section D).
The counterpart of the above problem for AD conservative gradients as in Definition 2 is tractable, illustrating major computational difference between Clarke subdifferential and AD conservative gradient. The proof is in Section D.4, by reduction to graph shortest path problem.
Proposition 1 (Finding two elements in autodiff conservative gradients is tractable) Given P PD)),, with conservative gradient Dp given by Theorem 2. finding two elements Dp(z) at a given input (or single element if Dp(x) is a singleton) is solvable in polynomial time.
CONCLUSION
We extended the "cheap gradient" principle to nonsmooth automatic differentiation with a flexible version of Haur-Strassen result: the overhead ratio of conservative gradients independent of the dimension. On the other hand, we showed that the potential gain in efficiency of forward AD for multiple directional derivatives limited due to an intrinsic connection to matrix multiplication. Finally, we have shown that for simple ReLU networks, the enumeration of Clarke subgradients is computationally hard, in contrast to the enumeration of conservative gradients.
The global picture is significantly different from the smooth case, with well understood "cheap gradient" principle that yields "cheap p directional derivatives", illustrating the specificities of nonsmoothness. Our results confirm the centrality of conservative gradients in nonsmooth AD and machine learning: they generalize gradients with a clear "cheap principle", contrary to concurrent notions. An important open question in this context is the complexity of subgradients, or, in other words, the existence of a "cheap subgradient principle". We conjecture negative answer in general.
ACKNOWLEDGMENTT AND DISCLOSURE OF FUNDING
The authors acknowledge the support of the AI Interdisciplinary Institute ANITI funding under the grant agreement AN P13A-10110004, The authors acknowledge the support of the Association nationale de la recherche et de la technologie (ANRT) and Thales LAS France, which contributed to Ryan B's 's grant. Jérome B. and Edouard P. acknowledge the financial support of Air Force Office of Scientific Research, Air Force Material Command, USAF, under grant numbers FA9550-19-1 7026 PA8655-22-1-700 and ANR MaSDOL 19-CE23-00119 7-01. Jérôme B. also acknowledges the

---

Published as a conference paper at ICLR 2023
support of ANR Chess, grant collaborators in the Thales LAS France, especially Andrei Purica, for helpful comments. We are grateful to Serge Gratton, Pierre Weiss and Pierre Boudier for useful reference suggestions.
TSE-P and the Centre Lagrange. We thank our
REFERENCES
Martin Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, Manjunath Kudlur, Josh Levenberg, Rajat Monga, Sherry Moore, Derek G. Murray, Benoit Steiner, Paul Tucker, Vijay Vasudevan, Pete Warden, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng. Tensorflow: A system for large-scale machine learning. In 12th USENIX Symposium on Operating Systems Design and Implementation (OSDI 16), pages 265-283, 2016. URL https: //ww. www.ix.eix..ommmmmmmmmm org/ system/ files/ conferenceeosdill osdi 1 abadi. pdf..
Akshay Agrawal, Robin Verschueren, Steven Diamond, and Stephen Boyd. A rewriting system for convex optimization problems. Journal of Control and Decision, 5(1):42-60, 2018.
Akshay Agrawal, Shane Barratt, Stephen Boyd, Enzo Busseti, and Walaa M Moursi. Differentiating through a cone program. Appl. Numer. Optim, 122):107-115, 2019.
Josh Alman and Virginia Vassilevska Williams. refined laser method and faster matrix multiplica- tion. In Proceedings of the 2021 ACM-SIAM Symposium Discrete Algorithms (SODA), pages 522-539. SIAM, 2021.
Raman Arora, Amitabh Basu, Poorya Mianjy, and Anirbit Mukherjee. Understanding deep neural networks with rectified linear units. International Conference on Learning Representations, Conference Track Proceedings, 2018.
Shaojie Bai, J Zico Kolter, and Vladlen Koltun. Deep equilibrium models. Advances in Neural Information Processing Systems, 32, 2019.
Paul I. Barton, Kamil A. Khan, Peter Stechlinski, and Harry A.J. Watson. Computationally relevant generalized derivatives: theory, evaluation and applications. Optimization Methods and Software, doi: 1111////156788.00171111 10. 10.1080/ 010.//105555878 2017 13774385.
33(4-6):1030-100 1072, 2018.
URL https: /ooi....... org/
Walter Baur and Volker Strassen. The complexity of partial derivatives. Theoretical Computer Science, 22:317-330, 1983.
Atilim Gunes Baydin, Barak A Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind. Automatic differentiation in machine learning: a survey. Journal of Marchine Learning Research, 18:1 13:1-43,22 2018.
L. M. Beda, L. N. Korolev, N. V. Sukkikh, and T. S. Frolova. Programs for automatic differentiation for the machine BESM. Technical report, Institute for Precise Mechanics and Computation Techniques, Academy of Science, Moscow, USSR, 1959.
David Bertoin, Jérôme Bolte, Sébastien Gerchinovitz, and Edouard Pauwels. Numerical influence of relu't (0) on backpropagation. Advances in Neural Information Processing Systems, 34, 2021.
Quentin Bertrand, Quentin Klopfenstein, Mathieu Blondel, Samuel Vaiter, Alexandre Gramfort, and Joseph Salmon. Implicit differentiation of lasso-type models for hyperparameter optimization. In International Conference on Machine Learning, pages 810-82 821. PMLR, 2020.
Mathieu Blondel, Quentin Berthet, Marco Cuturi, Roy Frostig, Stephan Hoyer, Felipe Llinares-López, Fabian Pedregosa, and Jean-Philippe Vert. Efficient and modular implicit differentiation. arXiv preprint arXiv:2105.15183 2021.
Jacek Bochnak. Michel Coste, and Marie-Françoise Roy. Real algebraic geometry, volume 36. Springer Science Business Media, 2013.
Jérôme Bolte and Edouard Pauwels. Conservative set valued fields, automatic differentiation. stochastic gradient methods and deep learning. Mathematical Programming, pages 1-33, 2020a.
10

---

Published as a conference paper at ICLR 2023
Jéröme Bolte, Tam Le, Edouard Pauwels, and Tony Silveti-Falls, Nonsmooth implicit differentiation for machine- leaming and optimization. Advances in Neural Information Processing Systems, 34, 2021.
Jérôme Bolte and Edouard Pauwels. A mathematical model for automatic differentiation in machine learning. In Conference on Neural Information Processing Systems, 2020b.
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL
http:// /ttp://gttthbb/t com/ comm///oll/////oooo jax.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Frank H Clarke. Optimization and nonsmooth analysis. SIAM, 1983.
Michel Coste. An introduction to o-minimal geometry. Istituti editoriali e poligrafici internazionali Pisa, 2000a.
Michel Coste. introduction to semialgebraic geometry, 2000b.
Courtier and F. Rabier. The use of adjoint equations numerical weather prediction. Ammospere- Ocean, 35(sup1):303-322) 1997. doi: 10. 1000000055000000000 1997.96899999998755 URL https: /doi. org/10. 080 000050000050000 997 .6877354.
Damek Davis and Dmitriy Drusvyatskiy. Conservative and semismooth derivatives are equivalent for semialgebraic maps. arXiv preprini arXiv:2102.08488 2021.
Jean-Guillaume Dumas and Victor Pan. Fast matrix multiplication and symbolic computation. arXiv preprint arXiv: 2016.
Patrick E Farrell, David Ham, Simon W Funke, and Marie Rognes. Automated derivation the adjoint of high-level transient finite element programs. SIAM Journal Scientific Computing, 35 (4):C369-C39) C39C3 2013.
Serge Gratton, David Titley-Pelquuin Philippe Toint, and Jean Tshimanga Iungaa Differentiating the method of conjugate gradients. SIAM Journal Matrix Analysis and Applicationn, 35(1):110-126 2014. doi: 10. 13/1137/108888848. URL https: /ooi .rwgorrgg/oooo 1tt.//...1.11111 11371.1201811999
A. Griewank and A. Rojas. Treating artificial neural net training as a nonsmooth global optimization problem. International Conference on Machine Learning, Optimization, and Data Science (pp. 759-770). Springer, Cham., 2019.
A. Griewank and A. Walther. Beyond the oracle: Opportunities of piecewise differentiation. In Numerical Nonsmooth Optimization (pp. 331-361). Springer. Cham., 2020.
Andreas Griewank. On stable piecewise linearization and generalized algorithmic differentiation. Optimization Methods and Software, 28, 07 2013. doi: 10. 10.1080/15556788.2111..999
Andreas Griewank and Christèle Faure. Piggyback differentiation and optimization. In Large-scale PDE-comstrained optimization, pages 148-164. Springer, 2003.
Andreas Griewank and Andrea Walther. Evaluating derivatives: principles and techniques of
algorithmic differentiation. SIAM, 2008
Andreas Griewank al. On automatic differentiation. Mathematical Programming: recent develop- ments and applications, 666):83-107 1989.
Sham M Kakade and Jason Lee. Provably correct automatic sub-differentiation for qualified programs. In Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
Kamil Khan. Branch-locking ad techniques for nonsmooth composite functions and nonsmooth implicit functions. Optimization Methods and Software, 33(4-6) 2018.
11

---

Published as a conference paper at ICLR 2023
Kamil A Khan and Paul I Barton. Evaluating an element of the clarke generalized jacobian of a piecewise differentiable function. In Recent Advances in Algorithmic Differentiation, pages 115-125. Springer, 2012.
Kamil A Khan and Paul Barton. Evaluating an element of the clarke generalized jacobian of composite piecewise differentiable function. ACM Transactions on Mathematical Software (TOMS), 39(4):1-25, 2013.
Kamil A Khan and Paul I Barton. A vector forward mode of automatic differentiation for generalized derivative evaluation. Optimization Methods and Software, 30(6):1 185-1212, 2015.
François Le Gall. Powers of tensors and fast matrix multiplication. In Proceedings of the 39th
international symposium on symbolic and algebraic computation, pages 296-303, 2014.
Yann LeCun, Yoshua Bengio, and Geoffrey Hinton. Deep learning, nature, 421 2015.
Adrian Lewis and Tonghua Tian. The structure of conservative gradient fields. arXiv preprint arXiv:2101.00699, 2021.
Jonathan Lorraine, Paul Vicol, and David Duvenaud. Optimizing millions of hyperparameters by implicit differentiation. In International Conference on Artificial Intelligence and Statistics, pages 1550-1552. PMLR, 2020.
Dougal Maclaurin, David Duvenaud, and Ryan Adams. Gradient-barl -aasee hyperparameter optimization through reversible learning. In International conference on machine learning, pages 2122- PMLR, 2015.
Sheheryar Mehmood and Peter Ochs. Automatic differentiation of some first-order methods in parametric optimization. In International Conference on Artificial Intelligence and Statistics, pages 1584-1594, PMLR, 2020.
Yu Nesterov. Lexicographic differentiation of nonsmooth functions. Mathematical programming. 104(2):699-700, 2005.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary Vetito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-peeformance deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, F. D'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems 32, pages 0124-8035. Curran Associates, Inc., 2019. URL http: /ttppp:pp/pppps...pptesp nwwries cc/psssrc///p////pp/pp/ 9015-pytorch- htwp::/www....aaaaa..eeeeeeeeeer pdf.
Barak A Pearlmutter. Gradient caleulations for dynamic recurrent neural networks: A survey. IEEE Transactions on Neural networks, 6(5):1212-1228, 1995.
R-E Plessix. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications. Geophysical Journal International, 167(2):455-503, 2006.
Maithra Raghu, Ben Poole, Jon Kleinberg, Surya Ganguli, and Jascha Sohl-Dickstein. On the expressive power of deep neural networks. In international conference on machine learning, pages 2847-2854. PMLR, 2017.
Sara Robinson. Toward optimal algorithm for matrix multiplication. SIAM news, 38(9):1-3, 2005.
David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams. Learning Representations by Back-propagating Errors. Nature, 323(6088):533--36, 1986. doi: 10.1038/323533a0, URL http::/www. nature. com/ art wttttceeee/33333aae
Stefan Scholtes. Introduction to piecewise differentiable equations. Springer Science Business Media, 2012.
12

---

Published as a conference paper at ICLR 2023
Alexander Schrijver. Theory of linear and integer programiming. John Wiley & Sons, 1998.
Alexander Shapiro. On concepts of directional differentiabilityy. Journal of optimization theory and applications, 66(3):477-487, 1990.
Stephen P Smith. Differentiation of the cholesky algorithm. Journal of Computational and Graphical Statistics, 4(2):134-144) 1995.
Volker Strassen al. Gaussian elimination not optimal. Numerische mathematik, 13(4):354-356, 1969.
Guanhua Wang and Jeffrey A. Fessler. Efficient approximation of jacobian matrices involving non-uniform fast fourier transform (nufft), 2021. URL httt::://arsii..rr/abbb////// 02912.
Robert Edwin Wengert. A simple automatic derivative evaluation program. Communications of the ACM, 7(8):463 464, 1964.
Virginia Vassilevska Williams. Multiplying matrices faster than coppersmiih-winnoradd In Pro- ceedings of the forty-fourth annual ACM symposium on Theory of computing, pages 887-898, 2012.
Ezra Winston and J Zico Kolter. Monotone operator equilibrium networks. Advances in neural information processing systems, 33:10718-10722 2020.
Philip Wolfe. Checking the calculation of gradients. ACM Transactions on Mathematical Software (TOMS), 844:3337-34, 1982.
13

---

Published as a conference paper at ICLR 2023
This is the appendix for "On the complexity of nonsmooth automatic differentiation".
CONTENTS
A Further comments, discussion and technical elements
14
B Proofs related to Section
15
C Proofs of Section 5.1
20
D Proofs of Section 5.3
21
A FURTHER COMMENTS, DISCUSSION AND TECHNICAL ELEMENTS
A.1 COMMENTS ON SECTION 3
A.1.1 COMPUTATIONAL MODEL IN SECTION 3.1
DAG representation and examples 3.1: We start with a remark regarding representations of
programs as directed acyelic graphs and use them illustrate the model of computation proposed in the main text. It reduces to that of arithmetic circuit complexity for dictionary composed of elementary arithmetic operations.
Remark 4 (Programs as directed graphs) A predecessor relation trivially describes directed acyclic graph (DAG). Therefore, a program is equivalently represented as a DAG, nodes corre- sponding either to input variables (empty predecessor) computation (nonempty predecessor). Directed edges connect predecessor nodes to their successors. Each computation node contains lower-level program (with single output), with the number of input edges being coherent with the number of arguments. The cost of a node is that of the underlying program and the cost of is the sum of the costs of its nodes. Nodes without outer edges output nodes. See examples in Appendix A.
We represent programs using the DAG representation as in Remark Let us define simple dictionary D = {+, x and introduce level O elementary program Po such that Pa(a, b) a + b meaning that Po computes the quantity a + b. Po identified with from the dictionary. We also introduce a level 1 program Pi such that (a b) a x (b + c). We can construct an equivalent level program, Q1 such that Qi b) ] = a + a x c, in this case, we have Pi QI, or [P]] [Qi] since they compute the same quantity. The level 2 program P2 is such that P2(a, b, c, d) (a + b) x (c+ d) Q1 c, d) + Pi c, d) and uses level programs Q1 and Pi in its computation nodes. The Directed Acyelic Graphs (DAGs) representing these programs are given in Figure 1. Assuming cost( +) = cost(x) = 1, we have (Po) 1, cost P1) = cost Q1) and cost(P2) (ost(Q1) cost(P) + cost(x) 6.
P
a
d
(a) Po
(b) P
(c) Q1
(d) P
Figure I: DAG illustrating different programs with dictionary x). (a) Po b) a
of level 0 which is identified with from the dictionary, (b) (a a(b c), of level 1.
(c) (a b) = ab + ac, level 1 and equivalent to Pt. (d) b2(a,b,c,d) = (a +b( d) = Qi c, Pa(b, c, d), of level
14

---

Published as a conference paper at ICLR 2023
A.2 COMMENTS ON SECTION 5
A.2.1 FORWARD AD AND CLARKE
Nesterov (2005) introduced the notion of lexicographic subdifferential, denoted here F for Lipschitz function F: R R. The construction of ê based on successive local approximations of F with directional derivatives, and one has ((r) 2 Fxx) for all such that the first term is well defined.
It is known that automatic differentiation can be used to compute directional derivatives, particu- larly the forward mode of automatic diffcrentiation (Griewank and Walther, 2008). Based on this observation, Khan and Barton developed several algorithms to evaluate elements of 6 F. based on directional derivatives (Khan and Barton, 2012: 2013; 2015). They concentrate on piecewise C functions, for example Scholtes (2012), and propose handle compositional structures with different restrictions on the function class considered, such as functions in abs- normal forms (Khan and Barton, 2012), broader classes (Khan and Barton, 2013; Barton et al., 2018).
All these procedures either require to evaluate directional derivatives (Khan and Barton, 2012; 2013), rely on forward chain rule propagation for lexicographic derivatives (Khan and Barton, 2015; Barton et al., 2018), which also require to maintain directional derivatives. For this reason, all these methods suffer from a multiplicative computational overhead ratio of the order of in the worst case, and it is not known if this could be improved (Barton et al., 2018), although efforts have been made in this direction (Khan, 2018).
A.2.2 MATRIX MULTIPLICATIONS
Remark 5 The lower bound described in Theorem 3 is sharp for linear ReLU network F as in (11) involving only square p x p matrices. Indeed, directional derivatives F in directions a1.... Ap- can be computed with roughly Lc(p) operations, using matrix multiplication algorithm realizing the c(p) bound, for example using the forward mode of AD Khan and Barton (2012; 2013). The naive Pr algorithm for forward evaluation performs roughly 2.p² operations which results in the bound (neglecting terms of order one in numerator and denominator).
comp(Fa, (ReLU, ReLU"))
c(p)
(FF)
2p
for this class of networks, to be compared with (5). Finally, remark that in the smooth case such complexity estimates reduce to gradient computation which can be done using backward algorithmic differentiation with constant multiplicative overhead ratio.
We denote by Fa the function Fa: (y, A) - (0, y, A) which computes p directional derivatives given point. Setting implementing Fa, we have shown that asymptotically, for any <>
lim supp- log(c(p)) log(p), since P' is an arbitrary program
comp(Fa, D [ReLU, ReLU')) 2--tte 400,
sup p. D.F=(1P)).reeP(PPttR
cost (Pr)
where the supremum is taken over all p and all functions F; IRPX4 - R implemented by a program Pp with dictionary D J [RcLU). is not known whether w > 2.
B PROOFS RELATED TO SECTION
Proof of Theorem 2:
Given program P in Section 3.1, the path differentiability of [P] is immediate by composition and the chain rule property. The associated conservative gradient Dy is constructed in Bolte and Pauwels (2020a).
We have the following cost estimates which can be deduced from the definition of the cost of a program in Section 3.1.
A Aggorihhm forward evaluation:
cost(P) = (elgoriith 1) = Σ
cost (g:)
(6)
i --++1
15

---

Published as a conference paper at ICLR 2023
Algorithm 1 forward evaluation with derivatives: Algorithm 1' with gd; instead of 9s line
cost(Algorithm 1')
cost (gd:)
i=p+1
Algorithm 2 backward AD cost:
cosk(backprop(P) P)) = cost Algoorithm 1') +
Σ
[p(())(coss(+t) + cost(x))
1=p+1
cost (gd:) + (co(t)(cost((t) + cost(x)) x)).
(8)
Algorithm 2 forward AD cost:
cost(forprop(P)) cost(Algorithm 1')
Σ
plpr(i)]cos((x) + p((or(t) - 1)cost(+)
p
cost (gd:) + plpr(i))cos((>> x) + p([pr(t) 1))oss(++)
(9)
Let us derive the complexity bound of Algorithm according to Algorithm 2.
Backward AD complexity result:
Using (8) and the fact that cost has value in RT, have
cost(backprop( P)) = Σ
cost (gd) + [pr(i) + cost(x))
m
cost. (gd:) [p](i)(cos((+) + coot(()))
cost(g:)) x
cost(g:)
i=p+1
cost (gd.) + [pr(t) + cost(x))
m
max imp- 1,m
Σ i=p+1
cost(gi),
cost(g))
where the inequality is due to factorization by the maximal value. Using (6), we obtain
cost cost(backprop(P2 P)) wa x cost(P)
where wa is given in (3). This proves point (i).
Forward AD complexity result:
Using (9) and the fact that cost has value in R$, have
cost(forprop(P)) cost (gd,) + plpx(i)) cost)( x) p((pr()) - 1)cost(+)
cg:) x cost (gd:) plp plp=( + p((pr(i) 1) 1)
Σ
cost(g:)
=
cost (gd,) + p|pz(i)) x) p(ppr(i) - ()coos(()) x
max 1mp+1m
cost(g;)
cost(g:).
<-pt1
where the inequality is due to factorization by the maximal value. Using (6), obtain
cost(forprop(P)) & wy x cost( P)
where wy is given in (3).
16

---

Published as a conference paper at ICLR 2023
B.1 JUSTIFICATION OF THE COMPLEXITY TABLE 1 OF THE DICUU-DICI
The proof of Corollary follows from Theorem by computing the relevant constants. They are shown in Table 1, let justify the proposed numbers.
Case (cost( x). cost(+)) Let us define g(a, b) = a b. To evaluate g. we need one operation from DRLLU. The derived program d related to 9. should satisfy d(a, (b, a) which does not require additional operation. Therefore, from Assumption we can deduce that cost(g) and cost(gd) L. We get the same result for cost(+) ) by applying identical reasoning.
Case 2 (cost(xe (o0,,cos((())) Let us define g(a) = C a. To evaluate we need one operation from DReLU. The derived program d related to g. should satisfy d(a) = which does not require additional operation from ReLUU- Therefore, from Assumption I we can deduce that cost = and cost(gd) 1. We get the same result for cost(+c) by applying identical reasoning.
Case 3 (log) Let us define g(a) = log(a). To evaluate 9, we need one operation from DRELU- The derived program d related 9, should satisfy d(a) 1/a, which requires the inverse operation from RELL. Therefore, from Assumption we can deduce that cost(g)
and cost(gd) 2.
Case 4 (cost(exp)) Let us define g(a) exp(a). To evaluate g. we need operation from DRELU- The derived program d related to 9, should satisfy d(a) g(a) which does not require operation from DRELU. Finally, from Assumption 1 can deduce that cost(g) = 1 and cost(gd) =
Case 5 (cost(inw)) Let us define g(a) = To evaluate g. we need one operation from DRELU- The derived program d related to g, should satisfy dl(a) = which requires one additional multiplication compute the square and one (-1) mudiplicatiion operation from D DRLLU- Finally, from Assumption I we can deduce that cost(g) = 1 and cost(gd) = 3.
Case 6 (cost (Re)U) Let us define g(x)) (eLU(w) max(x, 0). To evaluate g. need to evaluate the sign of x. The derived program ReLU can be computed also from the sign of x without further operation. We have cost(g) by hypothesis, but it is also reasonable to consider cost(gd) = 1 as both operations only require sign evaluation of the same object.
Remark 6 Since DRELU dictionary contains the ReLU function, we can build other non-smooth functions such as the maximum and the absolute value. For example, max(x, y} = ReLU y) + y [e y) ReLU(w) ReL.U-w).
B.2 AN EXTENSION OF TABLE 1
The justifications of the following are similar to Section B. simply taking into consideration different types of operations. Taking Cponlin CReLU
recover table We replace ReLU by ReLU which corresponds to its usage in practice and allows us to balance the cost of ReLU operations and that of multiplications.
The justification is the same as in Section B.1 taking into consideration different types of operations. For the ReL.U operation, the justification as follows.
Case 7 (x (xcost The operation has two argument and requires one sign evaluation and one multiplication in the worst case, so we assign it the cost + CReLU The differentiated program should compute the function (a, - (ReL.U()) a x ReLU'(b)). One can write a program to compute jointly g and d as follows: return (a x b,b, if > and (0, O, 0) if < O. This only requires bit sign check which cost is CReLU and multiplication. We therefore model this operation such that cos) cost(g)
+ CReLU-
Further refinements could be considered including various type of computational operations, such as memory moves, these are beyond the scope of the present paper.
B.3 ADDITIONAL ELEMENTARY NONSMOOTH PROGRAMS AND COST EXAMPLES
For simplicity, we do not discuss the dictionary and its related derived dictionary as there are many possibilities, one of them being DReLU and D' PROLU as all the considered operations be equivalently
17

---

Published as a conference paper at ICLR 2023
Table 2: Extension cost table, Chonlin 2 1 is the cost of nonlinear operations and CReLU 2 0 the cost of sign evaluation for ReLU or ReLU".
(+, x)
(c,, XC)
log
cxp
inv
ReLU
cost g
1
Cponlin
Cuoulin
Cnonlin
CReLU
p
2
1
1
1
costi
1
Cnuonlin
Cuonli
Cnoulin + 2
CReLU
cost
2
Suonlia
cost (g)
Cpantan
cost p
4
2
cost (g)
Cronlin
Chonltn
Caanlin
1+enala
cost(ggd) + 2cost x pr
5
3
$
4
$
3
5
cost. (g)
expressed with R.LU, We use the same framework as in B.2 and we identify the cost of comparing two real numbers with CReLU > For each program, g and associated derived program d, we let
cost(gd) + 2cost( x |prl
cost (g)
Table 3: Extension cost table. Cuoulin > 1 is the cost nonlinear operations and CReLU 2 0 the cost of sign evaluation for ReI.U or ReLU'. For simplicity CReLU is abbreviated on and Gnoulin is abbreviated Cal
g
(+,x) x
ELU
3-max-pool
cost(g g
1+ CR
2 CR + Cnl
153 + 8CR
+ 2nCR
n22 CR
|
2
1
9
n
cost (dt,dd. g
1
+ CR
2 + CR + Cnl
153 + SCR
n + 2ncR
n(2 + CR
cost((gd]
1
1
cost (g)
cost x ]|rr
4
cost (g)
ItCR
2+oR tua
ScR < 1.12
(2+cRR) )
w
5
<
$
$ 3
S 2
Case 8 (Absolute value and Leaky-ReLU) Recall that |x]] x O and -x otherwise. Simi- larly (ea = x if x > O and az otherwise, for some parameter a E (0, 1) so that both cases exactly the same. The reasoning and result exactly the same for both operations treat the absolute value. The construction is similar as what was proposed for xcost (ReLU) treated in the previous section.
Let be a program to evaluate multiplication that cost(g) + CReLU- Similarly it is possible built program which returns (x, 1) if x > 0 and (-x,- -1) otherwise, this computes (gd) and require the exact same operations so that cost(gd) (g)
| 1. in the worst case it requires one sign evaluation and one
CReLU.
Case 9 (ELU)
f(0)=
ale" ((e--)(
if x 2 with a > 0.
if
Let g be program to evaluate the ELU function, it requires a sign evaluation and in the worst case one nonlinear operation to evaluate e" one multiplication to evaluate ae", and one substraction to evaluate ac" Therefore, cost(g) CReLU Cnoulin + 2. The derived program d requires the same sign and returns 1 or ac" depending on the sign. This does not require additional operation and therefore the joint computation of and d satisfies cost(gd) cost(g).
18

---

Published as a conference paper at ICLR 2023
Case 10 (max-m-linear) Set number of inputs and m 2 2 number of linear functions which are parameters, represented by a matrix A and a fixed input vector of size n represented by x E IR". Setting maxm R to I. the function which evaluates the maximum of m numbers, we consider g program which evaluates the function A -- maxm (Az). Recall that x is fixed so that the number of imputs is m The multiplication requires m. (2n 1) multiplications and additions and the evaluation of maxm requires (m - 1)CROLLU as it requires m - pairwise comparisons. We therefore have gost(g) = m (2n 1) + (m 1)CReLU.
As for the derived program d, setting M, = 0 except for row number i which attains the maximum in g which is set to x, we have an element of a conservative gradient for g. It is possible to jointly compute g(A) and d(A) by invoking a program which returns (i) M.) where i is any index realizing the max and Mj is as discussed. This does not require more operations and we have therefore cost(gd) = cgst(gg) = m. x 2m + (m C)
Case 11 (Two dimensional max-pooling (3 x 3-max-pool) We consider a kernel of size 3 3 for simplicity. The goal is to differentiate with respect to the kernel weights for fixed input. Let denote a program implementing such a function, is of the same form as max-m-linear except that the matrix A is of size 9 x 25 (padding values at the boundary of the 3 x 3 patch, this gives 5 x 5 = 25 inputs and 9 outputs), but is sparse and can be parametrized by only 9 values, and the evaluation of the linear function for fixed 5 x 5 input only requires 9 x (9 8) 153 addition and multiplications. We then take the maximum of these 9 outputs so that and cost g) = 153 + SCaoulin- For the same reason max-m-linear, we have cost(gd) cost(g) 153 8Cnonlin-
Case 12 (1 (n-mm,,, ||- (|>)) Denote by g a program which evaluate the 1 norm on R". It has n inputs. addition, n. multiplication by -1 and
In the worst case, its evaluation can be done with pairwise comparisons. Therefore we have costíg) = 22n + nCRELU - 1. For the same reasons as all examples before, it is possible identify derived program d without requiring additional operation 1.
so that cost(gd) = cost(g) 2n + RCReLU
Case 13 (Median of numbers) Denote by a program that evaluates the median of n. numbers. This can be done by sorting the n numbers and outputting the value corresponding to ( which requires roughly log(n) operations, depending on the algorithm used. The sorting operation is permutation, one could apply the same permutation to the vector (1, (1.2..... n) without additional operation required. The number at position 12], call it is the index of the value associated with the median. Setting d to be the null vector in Rº with value 1 at position only, we have a selection in conservative gradient for the median with no additional operation required. Therefore in this case cost(g) = cost(gd).
Case 14 (Selection functions) This example encompasses virtually all examples used in machine learning and extends the median example above. Assume that . IRP - R is locally Lipschitz, given in the form
f(x) - fa(e)
where s: R (1....,mi) is an index selection function, and for each i 1,..., m. f:: TR. R is 2 function. Let g be program computing f. one possibility is to first evaluate s(x) at cost C2 and then evaluate [e(z)(2) at cost ct. As shown in Bolte and Pauwels (2020b), under very mild restrictions on and f (which should be expressed with logarithums, polynomials, exponentials etc. ...). the function
E. V (z)
is a conservative gradient for f. It can be seen that it is possible to evaluate jointly (g.gd) by first computing s, at a cost Ca, then evaluate fa and fs jointly at a cost ev.
cost(g) = C. + Cj cost(gd) C, cv
cost(
+ Cv
& C. + 5cs
cost(g) C. c
Co c
where we used exp & 5cj, the cheap gradient principle for smooth programs. This ratio is close to 5 if c, is negligible, we recover the usual ratio for smooth programs. It is close to 1 c. dominates, which is the case in the median example where f. just corresponds to coordinate number of the input and has constant derivative.
19

---

Published as a conference paper at ICLR 2023
C PROOFS OF SECTION 5.1
C.I PROOF OF THE MAIN RESULT
Proof of Theorem 3: Let U e RexP be an orthogonal matrix with entries in [-1,1] which columns denoted by u1,..., Up (with squared norm p). Assume that we have variables a matrix M E Rexp and two matrices A, B E RpxP with columns .1...., Up and b,,...,b respectively.
Consider the function
F: (2, 1 P Σ CI p
The pair (M, B) will be identified as y in the statement of the theorem. Considering the dictionary of elementary functions {+, x, RRRLUU, +c, xc). F has representation as a program Pp using the identity |d| = (t)U(t) + ReLU(-t) for all e R. We may construct Pp such that cost (Pp) 6p² + 2p 8y² where count 2p² p operation for each matrix vector multiplication to evaluate UB" Mæ (there are three of them), p multiplication by - to evaluate BU Mx. 2p application of ReLU (on U Mr and -UB" ( Mx). p additions of ReLU outputs evaluate p applications of the absolute value, p for the outer sum and for the division. Now consider the constraints
SUg Ma.) = u4, == P.
(10)
The set of matrices A, B, M satisfying this constraint is an open set, call it S. We now restrict our attention to this open set and argue that cost(P") does not change if the input variables are constrained to be in S.
We have for all i 1, ,,...,P and BA M) E S, the following directional derivatives with respect to variable
Fi, B, a, p
(o.) Ma, UUB Ma: by. Ma.
p
Setting the function G: (A, B)
- Ei,, M, = TM we have that is a polynomial and G(A, B, M) E b.y = BAT. Note that this does not depend on M.
Fix P' any program implementing the directional derivatives function (y, A) - (0, y, A) of F
described above, with dictionary {+, x, ReLU, ReLU', +c, x c), as the statement of the theorem.
Claim There program P2 on dictionary D {+, x, +C, xe) such that G [Pz] (on the whole space) and cost(P2)= cost((P) + p.
We use the DAG representation of programs as in Remark 4. Therefore P' is described by a DAG which node are either input nodes computation nodes implementing functions from PRELU- We will modify the program by simple modifications of the computation nodes. We may obtain program Po implementing on S with dictionary PRELU with cost & cost(P") + p by summing the outputs of P'. The ReLU' nodes P. represent semialgebraie function Coste (2000a;b) with values in finite set. Therefore, there is dense open semialgebraic set on which all ReLU nodes in Po are locally constant (Coste, 2000a, Theorem 6.7). Reducing S if necessary, we obtain program P1 on dictionary DRELU such that P1 Po on S by replacing each ReLU' node in Po by the corresponding constants. We have cost(Pi) P1) & cost(Po) (we replace computing nodes by constants). By Lemma cost(P1 & cost (Po) & cost(P") + p and G = [2]] (on the whole space). This proves the claim.
1, there is program P2 D such that cost(P2)
We may obtain program D2 implementing V G with dictionary D by backward algorithmic differentiation on Pz, that is D2 backprop( P2). we have therefore
comp(BA' BAT D) cost cost(PP, D2)
5cost( P2)
5p 5cost PP),
20

---

Published as a conference paper at ICLR 2023
where the first incquality is because D2 is program computing BA" for all A, B on dictionary D, the second is because adding computation increases the cost, the third is property of backward algorith- mic differentiation D and the last one is by construction of Pz. Note that Comp(BA', D) c(p) by definition, therefore we have the claimed lower bound
cost(P") cost cost(Pr)
dp) 5p c(p) 5p 5cost(Pp)
8p²
C.2 AN ADDITIONAL LEMMA
Lemma Let Q: IRP - be a polynomial and P1 be a program (without loss of generality of level
1) the dictionary D {+, x, ReI.U, U.UUU++ +c, xc), such that Q - [PP] for all inputs restricted to an open set S = RP. Then there is a level program P2 on the dictionary DD D2 (RReL[]
{+, +c, xc} such that Q2 [P&] (for all inputs in IP), Furthermore, if cost cost(xe) then, (2)) = (I))
Proof: : We use the DAG representation of programs as in Remark 4. Therefore P1 is described by DAG which node are either input nodes computation nodes implementing functions from D1. The function computed by Pi as well as each of its nodes are semi alii-alebbracc Bochnak et al. (2013); Coste (2000a;b). For each ReLU node in the graph representing P1 (assume that there N of them) we associate a number: the function ReLUU evaluated on its input (with the convention that Re)U'(0) finite set, by semialgebraicity. there is an open subset of S' c S such that G is constant on S (Coste, 2000a, Theorem 6.7). Consider P2 which computation graph the same as that of P1 except that each absolute value node is replaced by multiplication by the corresponding ReLU value (which is constant on S'). Then Q = [Pi] = [Pa] for all inputs in the open set S', All computation nodes of programs on D are multivariate polynomials and two polynomials which agree on open set are equal globally. This concludes the proof.
0). This defines semialgebraic function G: E2 {0, 1).N. As it has values in
D PROOFS OF SECTION 5.3
We investigate in this section the hardness of finding a Clarke subgradient for programs defined on the elementary dictionary Do (+:-, ReLU). We start with an equivalent representation of these programs as linear ReLU networks with skip connections and specific weight matrices. This equivalence preserve representation size up to polynomial factors. We will then prove hardness result on such ReLU networks. This will provide proof arguments for Theorem by the polynomial time equivalence of the two representation. We proceed similarly prove Proposition 1, using the equivalence with the two representations.
D.1 POLYNOMIAL TIME EQUIVALENCE WITH LINEAR RELU NETWORKS WITH SKIP CONNECTIONS
Given set of matrices Mj E (-1,0,1)Mxx, M2 e (-1,0,1)"2xm, ML-1
{-1, 0, 1ppe Ixpl 2, Mr, (-1, 0, 1,101.11222.... we consider the function F: IR R.
F M2
(M
(11)
where i;: IRP - R are given functions which apply to each coordinate, an activation function which either the identity the ReLU function. There is an obvious notion of size for this representation, corresponding to the number of free parameters (matrix entries and coordinates on which ReLU identity is applied), the size of the representation is PL-1 Σ/-1 P: x pi-1 + Pi-
function F given in (II) can be represented by a program on Do of equivalent size, this correspond to naive implementation. Similarly, any program P P(Dn) on p inputs and with a single output can be represented by a network as (11) which size is at most 18cost (P)3. Indeed, we may assume that cost(P) p/2 without loss of generality, otherwise, the program would not perform operations on some of the input variables and it could be simplified by removing variables which do not affect the output. Recall that m in Algorithm 1 is the memory footprint of P. in case,
21

---

Published as a conference paper at ICLR 2023
is m = = (o) the number of inputs plus the total number of operations. Note that have m & 3cost(P). Each operation +. or ReLU in the program can be represented by a m x m matrix composed with certain c: R R" which contribution to the Relu network size is at most (m3 + m) & 2m & 18cost( since m is integer and m & Bc)st(P).. There are cost P) such operations
that a program can be represented equivalently by linear Relu network, with L = cost layers which contribution to the network size is at most 18cost(PP) (P) so that the size of the resulting network at most 18cost(P) which is the desired bound since.
We have shown that working with functions represented as in equation (I1) is equivalent to work with programs in (Do) as possible to switch from one to the other at a cost of an increase of the representation size which is only cubic. Therefore we will from now on work with functions represented as linear relu networks with skip connections as in (11), and NP-hardness or polynomial time results on such function will be valid on PDD)) by the construction above.
D.2 FURTHER PROPERTIES OF LINEAR ReLU NETWORKS
Throughout this section F denotes a with representation in (11). This function is positively homogeneous, it satisfies F(0) 0 and it. By piecewise linearity, its Clarke subdifferential is polyhedron (see e.g., Arora al. (2018); Raghu et al. (2017)). The Clarke subdifferential is conservative gradient for this function, and we will associate to it a different conservative gradient, associated to Algorithm 2
Definition 2 (Autodiff conservative gradient) We consider a specific conservative gradient for F. it is given by (x) {M] DM D2... 12....7__. - diagonal matrix which entries respects the sign pattern of the corresponding activation function: the activation is identity, 0 if the activation is ReLU and the input is negative, if the input is positive and all elements in [o, 1J if the input is null. We have in particular
where for i = ,,..., L 1, D; is
(0( (Mf D M, Di
(12)
where in this case, diagonal entries of matrices D, corresponding to ReLU activations are arbitrary in [0, 1] and the remaining diagonal entries are (corresponding to identity activations).
The autodiff conservative gradient is associated with the algorithmic differentiation of natural numerical program implementing F as in Subsection 3.2. Furthermore, one can check that given a program P E P(Do), after the transformation outlined in Section D.I, have that Dj coincides with Dp in Theorem 2. In the following definition, Dp could fer example, the Clarke subdifferential of F or the algorithmic differentiation conservative gradient DT.
We consider the following problem.
Problem (Conservative gradient enumeration) Given matrices M1 e RRIXP, M2 E R22XP1 ... MM.-1 E PP-- 1 xp2-2, M, E RixpL-1, and functions L1...., $ f-11, consider F: E. R
the associated linear ReLU network with skip connections in (11), x E R² and Dp: IP = R
conservative gradient for F. Compute two distinct elements in Dp(a) or one element if it is singleton.
This problem enters the field of computational complexity as we have associated to it a representation size corresponding to the number of "free parameters" to be chosen: each matrix entry and the activation (ReLU or identity) corresponding to each coordinate, resulting in a number of parameters Pr.-1 + Σ P; X 1i + p:- In what follows, we will consider integral or rational entries for matrices and input x with the common notion of bit size. Schrijver (1998).
D.2.1 CLARKE ENUMERATION IS NP-IARD FOR RELU NETWORKS
The decision version of Problem 1, under the same assumptions, is to decide if there exists two distinct elements in Dy (z:), that is, decide if Dr(x)) not reduced singleton.
Theorem (Finding two Clarke subgradients is NP-Hard) Decision version of problem (1) with matrix and vector entries in (-,00,1] and Dr 2 F is NP-hard.
22

---

Published as a conference paper at ICLR 2023
Sketch of proof: We encode a boolean formula p boolean variable, in a linear ReLU network with p inputs, size proportional to that of и. We do so by replacing "or" operations by maxima, "and" operations by minima, negation by multiplication by -1 and adding ReLU operations to the result. Using Lemma 3 in appendix D.5, the resulting F is represented by a linear ReL.U network. By construction, is global minimum of F so 0 8 F(), and F takes positive values if and only if is satisfiable if and only if F(O) * {0}. We detail this proof in coming sections.
Theorem illustrates the hardness enumerating Clarke subgradients of linear ReLU networks. For F as in (11) and x E RP, 2 Fx)) is not a singleton if and only if F is not differentiable at x, therefore:
Corollary (Deciding non-differentiability of a NN is NP-Hard) Given a linear ReLU network as in (11) with matrices as in Theorem 5 and z E R,, deciding F is not differentiable at is NP-hard.
In the coming section, we will provide proof for Theorem 5 and Corollary 2. By the polynomial time equivalence of the representation of programs in P(Db) and functions as in (11) detailed in Section D.I, this proves Theorem
We add remark on lexicographic subdifferential. It follows from (Barton et al., 2018, Proposition 2.7) that, for linear ReLU network as in (11), the lexicographic subdifferential Nesterov (2005) is the set of neighboring gradients and is contained in Clarke subdifferential.
Corollary (Finding two lexicographic subgradients is NP-Hard) Theorem 5 remains true if DF is the lexicographic subdifferential.
D.3 PROOF OF THE MAIN HARDNESS RESULT
Preliminary on 3-SAT We will use reduction to 3-SAT problem which is among the most well known NP-complete problems. Recall that a boolean formula is built from boolean variables, and operators: AND (conjunction, denoted ^) OR (disjunction, v) and NOT (negation, -). A literal, is either a variable or the negation of a variable. clause a disjunction of literals (or a single literal). A formula is in conjunctive normal form (CNF), if it is conjunction of clauses clause. 3-SAT is the decidability problem associated to CNF formulas with clauses containing 3 literals, such formulas called 3-CNF formulas.
Example 2 The formula (b) v b2 v -b3) (bi by v -bs) (-b2 v -b v be) is 3-CNF with boolean variables b1,..., by and 3 clauses.
Problem 2 (3-SAT) Given N N and a boolean function with boolean arguments bi..... bp represented by a 3-CNF formula with n clauses, decide if there exists an assignment (b1...., bp)) {0, 1 such that π (b.,...,b) bp) 1.
Proof of Theorem 5:
The reduction is 3-SAT.
Consider a 3-CNF function in p variables bi....,b, bp with n clauses of size 3. We may assume without loss of generality that n is of the form 2* for & E N by adding clauses which are always true and increasing the number of clauses by a factor at most 2. We will consider p real variables .I..... Zp. Consider the first clause of æ, say for example (&z by -by). We associate to each literal the corresponding variable æ if the literal is equal to a variable, and - the negation of the corresponding variable, for example æ1, 22, T33 These are combined using ReLU max resulting in the expression ReLU [11 T2, 106.-230 x3))-
We proceed similarly for each clause, obtain the max is over three numbers. The max numbers the same as the max of 4 numbers (by copying one of the inputs) and, according to Lemma can be represented by ReLU network with 2 ReLU layers of size at most x 2 6 with weight matrices in (,1,0,1).
2 expressions involving ReLU max where
We may therefore represent the ReLU max expressions with network with p inputs and n outputs, with 3 ReLU layers (2 for each max and one for the outer ReLU) size at most 6n (6 nodes for each max) involving matrices with entries in (-,,0,1). These expressions
23

---

Published as a conference paper at ICLR 2023
combined using the operator min applied to the n 2 clause. Thanks to Lemma 3 again, using min(a, a. b) max a. - b]. the max over the 2 numbers can be expressed with layers of size at most x 2 n
We call the resulting network F. has representation as in (11), with matrices with entries in
Zg 1, 1} in Problem 1. It contains (n))(n) 3 ReLU layers of size at most 6n and it has therefore description which size polynomially bounded in n which is proportional to the bit size representation of the 3-CNF formula ...
Example 3 If the 3-CNF formula is given by (b) v b2 v -ba) a (br v by v -bs) (-b2 v by v b6) (b2 -62 v be) with p boolean variables and 4 clauses, we will obtain a network computing the following expression in 6 real variables ;
F(x1....,26))
mmin(ReLUmaxx()). (..88..))))) -2a)), Reeeeeemmmmm.eaaa -I5)) RcLU(max( -22, -3, m6)), ReL.U(max(ns))) - 16))))
We have the following rules for min and max over real numbers a, b, (we use the convention sign(0) 0).
max(a, b, c)
(a 0) (b> (c> 0).
max(a, 0
(iggg))
0.
min(a,b,c) > 0
+
(a> 0) ^ (b 0) (c 0).
min(a,b,c) >
min(sign( sign (igg 0.
a 0
(-a < 0)
sigu(a) > 0.
. ReL.U(maxxssgnn(0) sign(b), sign ())) 10.
Because of the min oReLU structure, we have F(x) > 0 for all x, furthermore, F(0) 0, so that O a global minimum of F and E 20 F(O). For any x, we have F(x) if and only the output of each max is positive, if and only if each max clause contains a positive argument. We therefore have that F(x) > 0 if and only if (sign(z)) > where sign is the coordinatewise application of the sign, taking value 0 at O.
We have the following chain of equivalence
F()) # [0}
3x e RP, F(x) #
B E RP, F(x) >
Dx e R2, i; 0 1, .p......) F(x) > F(sign(x)) > 0
3x e RP,
* 0(Vi ..)
Br E 13", F(x) >
1x (-1,1)", T(b) 1, b. b- = 11=1...p)
where outputs if the boolean argument is true, and otherwise. The first equivalence is by Lemma 2, the second is because F > 0, the third is because F is continuous, the fourth is by the discussion above and the fifth is obvious because all possible (-1, 1} patterns can be described coordinatewise sign applied vectors in R2' with nonzero entries. For the last equivalence, for ar; E (-1, we set x; 1. Each ReL.U max applied to the sign vector corresponds to clause and its output is in {0, 1). The output of each ReLU max clause is if and only if at least one of its argument is 1, if and only one of the litteral of the corresponding disjunction is and only if the disjunction applied to the corresponding boolean variables is true. Otherwise, it is O. Similarly, the min combination has positive output if and only if all max outputs are 1 if and only all the disjunctions applied to variables b, true.
b; 0 if x;
and by
This shows that Problem is NP-hard, because 0 E 2" F(F( and F*F(()) * {0) if and only if there exists two distinct elements in 8 FOF()).
24

---

Published as a conference paper at ICLR 2023
D.4 PROOF OF FEASIBILITY FOR AUTODIFF CONSERVATIVE GRADIENT
The counterpart of Problem for AD conservative gradient in Definition 2 is tractable, illustrating a major computational difference between Clarke subdifferential and AD conservative gradient. The proof is in Section D.4, by reduction to a graph shortest path problem. By the polynomial time equivalence between linear ReLU network and programs on { -, ReLU) proved in Section D.I. this proves Proposition I.
Proposition 2 Problem (1) with matrix entries in Q and Dr = Di is polynomial time solvable.
Proof of Proposition 2: Consider the following polynomial expression:
Mf (Q1 + (0-
(13)
where we decomposed D; Q: + Qi in Definition 2, such that Q₂ is constant, diagonal, with zero entries except for the 1 entries which are enforced by the network activation and sign pattern: strictly positive activation before application of ReLU when network is evaluated at æ, or identity
activations. Furthermore, Q: contains qi Pi diagonal variables to be chosen in [0, 1J corresponding to the zero activation pattern before application of ReLU, for = 1...., L - 1. The strictly negative values before application of ReLU do not play additional role, they correspond diagonal entries constrained to 0 in both Q and Q. i == I..... L - 1. Note that polynomial is constant on a box and only if it is constant so the polynomial expression in (13) is constant when diagonal entries constrained in [0. 1]. if and only if it is constant. So the problem reduces to decide if the polynomial expression in (13) is non constant, with respect variables Q1, QL- We show that this reduces to graph connectivity problem over 2 E q vertices and edge weight given by partial products in (13).
First, the problem can be reduced to finding a non-zero value in the expression in (13). Indeed, one can substract the value obtained choosing Qi 0, i I...., L - 1 and use the following block representation:
(MF -MP)
(
Q1
Q1 0
)))
MM 0
O MT_1) ) ( M ....MMM-..AAA
(Q (QL1+ + QL-1
)im1)
0
M MY
= (Q + Q1) ... M (Q, + QI-
(14)
Therefore, expression (13) is nonconstant if and only if expression in (14) takes nonzero value for some assignment of Q...... Qr.-1. The number of variables in (13) and (14) is the same and they have exactly the same form. Therefore we assume without loss of generality that the problem is to
decide if the polynomial expression in (13) is not equal the null polynomial.
Expression (13) is a vector function each of its coordinates being a polynomial function. It is not uniformly null if and only if and only if there exists a coordinate which not the null polynomial, so we may add diagonal matrix Qo with po p diagonal entries in [0, 1] (and Qo 0 for the sake of symmetry) and Mo E R I the vector of all ones and find a nonzero value for the product
Mf (Qu Q6) (Q) + ... M.....M/.. + QL-1)ME.
(15)
1. denote by dj E [0. 1]st. the vector containing the diagonal entries of matrix Q. this corresponds exactly to the variable diagonal elements of D; in Definition 2. Denote by P(do.... di) the obtained polynomial, P is multilinear in do.... dz-1. that is, it has an affine dependency for one block vector if the others are fixed. Therefore the hessian of P has zero diagonal blocks and the function is harmonic (hessian has zero trace), as a consequence, the maximum principle for harmonic functions entails that its maximum and minimum on any polytope are attained at vertices.
Expression (15) is now real valued and therefore defines polynomial. For each 0=
For 0,..., L. denote by A, = TR&:, the convex hull of the origin and the canonical basis vectors, this is a qi dimensional simplex with nonempty interior. The polynomial P in (15) is identically zero and only if it vanishes on the product of simplices Ao x x... x Å (which has non empty interior). if and only if it vanishes maximum principle. In other words, P is not identically zero, if and only if it contains a nonzero element when each d; is restricted to be an element of the canonical basis (zero vector with exactly one nonzero entry) or the null vector.
the product set of the edges of these simplices by the
Define graph with layer structure:
25

---

Published as a conference paper at ICLR 2023
The source layer V-1 contains single source node, V--1,1
The zero-th layer Vo contains qo nodes Vo,1 V0.00
Recursively, the i-th layer Vi contains y nodes Vi ... Vi.go for 1... L
The sink layer VE contains single node node VL
We connect nodes between consecutive layers, respecting the order induced by the layer structure. -1. L- and
For
o,..., T, with we connect layers V; and Vj as follows
Compute the quantity
R
MTQm x MT.
where if 4 1 the product reduces to the identity (R Mz))
For k = 1.... and = 1........ 9j. add an edge with between Ug and Ej if Rg,1 # O.
The resulting graph has number of nodes equal the number of ReLU functions in F plus p additional nodes and the source and sink nodes. Computation of edges can be done in polynomial time: it requires most 4(L 1)² matrix products involving most 2L 1 matrices. Indeed the product of m matrices has polynomial time complexity in the representation bit size of the m input matrices.
In this graph, directed path from the source to the sink visits each layer at most once, and in that case it visits single node. Each such path corresponds to monomial with nonzero coefficient appearing in the polynomial P in (15) by construction of the graph structure. Conversely each given monomial in (15) is uniquely associated to path the graph which corresponds to the nodes associated to variables in the monomial. Therefore, the source is connected to the sink if and only if there a nonzero monomial in (15). and only the corresponding polynomial is nonzero. Furthermore, each path corresponds to the evaluation of the program edge of the product Ag x ... x A Therefore finding a path connecting the source to the sink allows to compute nonzero element in the product and if no such path exists, the polynomial is identically zero
nonzero coefficient
So have shown that the truth value of problem with DF Dg. is the same the source being connected to the sink by directed path in the graph we defined, which has size polynomialy bounded compared network size. Connectivity can be solved, for example using Dijkstra's algorithm, in time ((V|3) where [V] is the number of nodes (or vertices). A path represents nonzero element of (0) and no such path exists, we conclude that (0( {0). This shows that the problem is solvable in polynomial time and concludes the proof.
D.5 ADDITIONAL LEMMAS
The following lemma provides a characterization of singleton subgradient for linear ReL.U networks.
Lemma Let F be a linear ReLU network, then F=F( {0} if and only F constant.
Proof: If F is constant, the result is immediate because F = 0. Now, suppose that 6 F(0) {0)- We know that F is piecewise linear and there exists finite set of polyhedron whose union is RP. where F is affine linear over each polyhedron. Furthermore, F is positively homogeneous, therefore for each æ I2, F=)(() d0 F(xx) with À O. Setting R IRP, the full measure where F is differentiable, one has that for all x R and
F()) = conv v IRP, 3yk k
0 with yk E R, U/
F -
(0)-
k-
Therefore, each affine part has zero derivative on each polyhedra and by continuity we conclude that F is constant.
The next lemma describes an explicit representation of maximum of finitely many numbers using a ReLU network with weights in (-,,0,
26

---

Published as a conference paper at ICLR 2023
Lemma 3 Civen k E N, k O, there exists F, ReLU network with k ReLU layers size most
3 x 2 and weight matrices with entries in (- 0, 1}. with = 2 inputs such that for any x RR,
F())
max Fz-
i=
Proof: : We proceed by recursion on k. Note that for any 21,42 R
11 x2) = ReLU((11 - z2) + x2 = ReLU(<1 -22) + ReLUU(2)) - ReLU(-=2).
matrices
-1 0
B (1 1 -1)
-1
The function Fi: IR² - - given by
(x)x) BReLU(Ar)
satisfies (i( = max(#1, x2)] This proves the result for it
Now assume that for & > 1, we have network with ReLU layers of size most 3 x 2 represented by matrices "1....,,MM. with entries in {-1, 0, 1}. such that the corresponding ReLU network, as in (11) Fx: R _ IR satisfies for all x E R
Fx(x)
max Tz.
i=l...,," 2
Set Fs the concatenation of two copies of Fk, that is Fx: R2 - RR? such that for all æ, y R22
F2 y) = - 2i
-I.,, 26
The matrices representing F. can be described in block form
[Me
M. a.)
M,
R (2pp))
2* and pı: 1. This network made of k layers of size at most 3 x 2 it has 2 inputs and two outputs and its weight matrices have elements in {-1, 0, 1). The block representation the last matrix of this network is the form
for 1,..., + where Po
(Me+1
O
E 2xl
MK+1
where is the size of the row vector Mk+1- We have
A x
Mc+1
00
1
x
MM411
Mc+1
M
0
Mc+1
-M&+1
3x1121
0
-1
-M++1/
Fi(Fa(x,, v) for all x, e 22 In matrix notation We set FL+1(2,)) v) Fx(l) have
FR+1(x,))) BReLU(AF,(Z,.))
The involved matrices are Mr-2 B, A x Ma41 and Mr M,. They all have entries in (-1, and the corresponding network has layers of size at most 3 2 The result then holds by recursion.
27


Reviews:
[{"review_id": "pvStxi3Z9js", "reviewer": "Reviewer_vqhQ", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper utilizes the concept of a conservative gradient to bound the complexity of non-smooth derivatives, and shows that the computational complexity for backprop with conservative gradients has a computational overhead ratio independent of dimension (Theorem 2), and hence generalizes Baur-Strassen to non-smooth programs. \n\nThe authors show that in general for other differentiation oracles (like directional derivatives, etc) for the forward mode the overhead depends upon the dimension p (for example is O(p^\\omega) for directional derivatives), and therefore can be limited.", "summary_of_the_review": "The authors extend the cheap gradient principle and show that it is independent of the dimension. They further show that bound the forward auto-diff complexity will have an overhead due to matrix multiplication. ", "strength_and_weaknesses": "The paper points out the different complexity for computing derivatives for non-smooth functions, and shows a variety of results regarding the complexity ratio of different derivative methods. \n\n- the paper provides and rigorously prove bounds for the computational hardness non-smooth AD as well as for different generalized gradient methods.\n- The authors show that calculating distinct elements in a clarke differential of a given point is NP-hard.\n- The proofs (as much as i could follow them) seem correct.\n\nHowever, i think the proof sketches could be made more clearer in the main text, since it was a bit hard to follow.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}, {"review_id": "PmsIcnR38o", "reviewer": "Reviewer_LBds", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper studies the computational cost of forward and backward automatic differentiation (AD° for *nonsmooth* programs, based on the notion of conservative gradient introduced by Bolte and Pauwels (2020a, b).\nThe relative cost of backward mode AD/backpropagation (compared to the original program/function) is proved to be independent of the dimension for path-differentiable functions (that encompass semialgebraic and definable locally Lipschitz functions). Theorem 2 is a nice result that extends results that were known in the smooth case under the name of \"cheap gradient principle\" (Baur and Strassen, 1983).\nLastly, the authors show that such results on conservative gradients do not hold for other nonsmooth alternatives: the cost ratio for computing $p$ directional derivatives depends on the dimension $p$ unless the matrix exponent is equal to 2 (Theorem 3), and finding two distinct Clarke subgradients in the Clarke subdifferential is NP hard (Theorem 4, by reduction to 3-SAT).", "summary_of_the_review": "Well-written and theoretical sound paper on important recent development in automatic differentiation for non smooth functions, in particular neural networks.", "strength_and_weaknesses": "Strengths:\n- Nonsmooth operations/layers are ubiquitous in model neural networks (ReLU activation, sorting operations, max pooling, etc), and backpropagation is the default algorithm to train these networks. Having a rigorous theory of what backpropagation actually computes, and at which cost, is of paramount importance.\n- The concept of conservative gradients fills this need, and thus studying it in more detail is important for our understanding of theory and practice. \n- The paper is well-written, with a wealth of examples, a motivated impact.\n\nWeaknesses:\n- no important ones (see minor clarifications below)\n\n### Minor remarks and question\n- In Theorem 2, ii, an upper bound on the ratio is derived where $\\omega_f$ depends on the dimension $p$. For completeness, one may add that this upper bound is reached, for example in the case where $|pr(i)|$, and the costs of $g_i$ and $gd_i$ are independent of $i$. Otherwise it could just be a bad upper bound.\n- In the first paragraph of Section 5.1: is it known than the matrix multiplication exponent $\\omega$ is $> 2$? If $\\omega$ were to be equal to $2$, then the overhead ratio of order $p^{\\omega - 2 + o(1)}$ could in fact be dimensionless. This is later discussed in the \"Consequences\" part (and at the end of Section A.2.2) but it could be mentioned earlier: in all rigor, there is currently no result that says $\\omega > 2$, and so the derived ratio could be dimensionless (though all practical evidences are against this so far), which contradicts \"in contrast with the smooth case with essentially dimensionless overhead.\".\n- The ReLU requires only 2 bits of encoding > is it not a single bit (a binary value encoding the sign of the input)?\n\n\n### Very minor cosmetic:\n- the authors may consider configuring hyperref through \\hypersetup in order to avoid the flashy green and red boxes (see 2nd answer here: https://tex.stackexchange.com/questions/823/remove-ugly-borders-around-clickable-cross-references-and-hyperlinks)\n- when the citation is part of the sentence, it should not use parentheses: \"see Griewank et al (1989)\". When it is not part of the citation, it should use parentheses: \"It is at the core of modern learning architectures (Rumelhart et al., 1986).\". Use \\citet and \\citep respectively.\n- , little is known about the nonsmooth case. > \"but\" little is known about the nonsmooth case?\n- Standard computational practice of AD consists of: \"consists of\" means \"is composed of\". You mean \"consists in\"\n- in an ML > in a ML\n- In particular,they > missing space after comma\n- To conculde\n- Let $F, G$ be locally Lipschitz continuous mapping > mappings\n- Lipchitz\n- Path differentiable functions are ubiquituous > ubiquitous\n- consist of combinations outputs of > consists of combinations *of* outputs of\n- is an optimization algorithms > algorithm\n- The function F may describe, for instance, a ReLU feedforward neural network empirical loss, parameterized by $p$ real weights, with $q$ inputs: isn't $p$ the number of inputs and $q$ the number of parameters in the previous sentence?\n- then evaluate $p$ scalar product > products\n- for example $F$ is feedforward neural network > $F$ is a feedforward\n- Case 13: This can be done by sorting the n numbers and output the value value > and outputting + \"value\" duplicated\n- Therefore, the each affine part has zero derivative > remove \"the\"\n- fintiely", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}, {"review_id": "C_Gjj87wpK", "reviewer": "Reviewer_9fsL", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "  In this paper, the authors consider the problem of computing subgradients of non-smooth non-convex functions that are given as programs that are based on the small set of elementary functions. Their main contributions are the following:\n  1. They analyze the computational complexity of the algorithm suggested by Bolte and Pauwels 2020a for forward and backward evaluations of a relaxation of subgradients called conservative gradients. An important corollary of their analysis is that when we apply backward evaluation, the cost of evaluating the gradient is only a dimension-independent constant time larger than the cost of evaluating the function value. The authors also show that this is not the case for the forward evaluation of the gradient.\n  2. They show that given a program P that computes a function $f$ using only the operations $\\{ +, -, \\mathrm{ReLu} \\}$, and a point $x$, deciding whether $f$ is differentiable at $x$ is NP-hard. In contrast deciding whether $f$ is path-differentiable at $x$, which means that there is a unique conservative gradient, is possible in polynomial time.", "summary_of_the_review": "  The contribution of this paper is interesting but I am not sure if it is enough to guarantee acceptance at ICLR. For this reason I am slightly leaning towards rejection.", "strength_and_weaknesses": "Strengths\n======================\n  - The problem of theoretically understanding the complexity of computing subgradients of non-smooth functions is very important and interesting.\n\nWeaknesses - Comments\n========================\n\n  1. The main problem with the paper is that most of the material is a restatement of results provided by Bolte and Pauwels 2020a. In particular, the algorithms for forward and backward computation of conservative gradients are provided in exactly the same form in Bolte and Pauwels 2020a, at least in their arxiv version. Given that, the novelty of this part of the paper is the proof of Theorem 2, which can be found in page 15 and mostly follows a very simple induction argument similar to the existing proofs for the smoothed case.\n\n  2. I am not sure I understand the motivation of Theorem 3. Why would we need to compute the directional derivatives since we can compute the whole subgradient?\n\n  3. Related to Theorem 4: The problem of given a function $f$, a vector $v$, and a point $x$ verify whether $v$ is a subgradient of $f$ at $x$ is solvable in polynomial time or not?\n\n  4. If the answer to question 3 is positive then I am puzzled with Theorem 4 (b). The authors state that deciding the differentiability of a non-smooth function at a point $x$ is NP-hard. But I believe that the correct statement is that deciding the non-differentiability of $f$ is NP-hard. In fact, if the answer to 3 is positive, deciding the non-differentiability should be NP-complete, since the existence of two subgradients is a witness of the non-differentiability of $f$ at $x$. This means that deciding the differentiability is coNP-complete and we do not expect it to either be in NP or being NP-hard.\n\n  5. There is a small issue with Theorem 4 that usually is swept under the rug, but when proving computational hardness results it needs attention. What is the representation of the numbers that you assume in Theorem 4?  \n  **(a)** If you use binary representation, then the hardness result of Theorem 4 implies that there is no algorithm with complexity polynomial in $p$ (the dimension) and $\\log(1/\\epsilon)$, where $\\epsilon$ is the accuracy that we choose to represent the numbers.  \n  **(b)** If you can show the hardness result even when you use unary representation, then the Theorem 4 implies that there is no algorithm with complexity polynomial in $p$ and $1/\\epsilon$.  \n  If the case is (a) then this would be a less interesting result because the ML community focuses on algorithms with running time $\\mathrm{poly}(1/\\epsilon)$. In any case, a discussion about this is missing and is important.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: notable-top-25%