Fractal Structure and Generalization Properties of Stochastic Optimization Algorithms
Year: 2022

Abstract:
Understanding generalization in deep learning has been one of the major challenges in statistical learning theory over the last decade. While recent work has illustrated that the dataset and the training algorithm must be taken into account in order to obtain meaningful generalization bounds, it is still theoretically not clear which properties of the data and the algorithm determine the generalization performance. In this study, we approach this problem from a dynamical systems theory perspective and represent stochastic optimization algorithms as random iterated function systems (IFS). Well studied in the dynamical systems literature, under mild assumptions, such IFSs can be shown to be ergodic with an invariant measure that is often supported on sets with a fractal structure. As our main contribution, we prove that the generalization error of a stochastic optimization algorithm can be bounded based on the ‘complexity’ of the fractal structure that underlies its invariant measure. Then, by leveraging results from dynamical systems theory, we show that the generalization error can be explicitly linked to the choice of the algorithm (e.g., stochastic gradient descent – SGD), algorithm hyperparameters (e.g., step-size, batch-size), and the geometry of the problem (e.g., Hessian of the loss). We further specialize our results to specific problems (e.g., linear/logistic regression, one hidden-layered neural networks) and algorithms (e.g., SGD and preconditioned variants), and obtain analytical estimates for our bound. For modern neural networks, we develop an efficient algorithm to compute the developed bound and support our theory with various experiments on neural networks.

heading: 1 Introduction
text: In statistical learning, many problems can be naturally formulated as a risk minimization problem
minw2Rd n R(w) := Ez⇠⇡[`(w, z)] o , (1)
where z 2 Z denotes a data sample coming from an unknown distribution ⇡, and ` : Rd ⇥ Z ! R+ is the composition of a loss and a function from the hypothesis class parameterized by w 2 Rd. Since the distribution ⇡ is unknown, one needs to rely on empirical risk minimization as a surrogate to (1),
minw2Rd n R̂(w,Sn) := (1/n) Xn i=1 `(w, zi) o , (2)
where Sn := {z1, . . . , zn} denotes a training set of n points that are independently and identically distributed (i.i.d.) and sampled from ⇡, and model training often amounts to using an optimization algorithm to solve the above problem.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Statistical learning theory is mainly interested in understanding the behavior of the generalization error, i.e., R̂(w,Sn) R(w). While classical results suggest that models with large number of parameters should suffer from poor generalization [SSBD14, AB09], modern neural networks challenge this classical wisdom: they can fit the training data perfectly, yet manage to generalize well [ZBH+17, NBMS17]. Considering that the generalization error is influenced by many factors involved in the training process, the conventional algorithm- and data-agnostic uniform bounds are typically overly pessimistic in a deep learning setting. In order to obtain meaningful, non-vacuous bounds, the underlying data distribution and the choice of the optimization algorithm need to be incorporated in the generalization bounds [ZBH+17, DR17].
Our goal in this study is to develop novel generalization bounds that explicitly incorporate the data and the optimization dynamics, through the lens of dynamical systems theory. To motivate our approach, let us consider stochastic gradient descent (SGD), which has been one of the most popular optimization algorithms for training neural networks. It is defined by the following recursion:
wk = wk 1 ⌘rR̃k(wk 1), where rR̃k(w) := rR̃⌦k(w) := (1/b) X
i2⌦k r`(w, zi).
Here, k represents the iteration counter, ⌘ > 0 is the step-size (also called the learning-rate), rR̃k is the stochastic gradient, b is the batch-size, and ⌦k ⇢ {1, 2, . . . , n} is a random subset drawn with or without replacement with cardinality |⌦k| = b for all k.
Constant step-size SGD forms a Markov chain with a stationary distribution w1 ⇠ µ, which exists and is unique under mild conditions [DDB20, YBVE20], and intuitively we can expect that the generalization performance of the trained model to be intimately related to the behavior of the risk R(w) over this limit distribution µ. In particular, the Markov chain defined by the SGD recursion can be written by using random functions h⌦k at each SGD iteration k, i.e.,
wk = h⌦k(wk 1), with h⌦k(w) = w ⌘rR̃k(w). (3) Here, the randomness in h⌦k is due to the selection of the subset ⌦k. In fact, such a formulation is not specific to SGD; it can cover many other stochastic optimization algorithms if the random function h⌦k is selected accordingly, including second-order algorithms such as preconditioned SGD [Li17]. Such random recursions (3) and characteristics of their stationary distribution have been studied extensively under the names of iterated random functions [DF99] and iterated function system (IFS) [Fal04]. In this paper, from a high level, we relate the ‘complexity of the stationary distribution’ of a particular IFS to the generalization performance of the trained model.
We illustrate our context in two toy examples. In the first one, we consider a 1-dimensional quadratic problem with n = 2 and `(w, z1) = w2/2 and `(w, z2) = w2/2 w. We run SGD with constant step-size ⌘ to minimize the resulting empirical risk. We simply choose ⌦k ⇢ {1, 2} uniformly random with batch-size b = 1, and we plot the histograms of stationary distributions for different step-size choices ⌘ 2 {1/100, 1/3, 2/3} in Figure 1. We observe that the support of the stationary distribution of SGD depends on the step-size: As the step-size increases the support becomes less dense and a fractal structure in the stationary distribution can be clearly observed.
This behavior is not surprising, at least for this toy example. It is well-known that the set of points that is invariant under the resulting IFS (termed as the attractor of the IFS) for
the specific choice of ⌘ = 2/3 is the famous ‘middle-third Cantor set’ [FW09], which coincides with the support of the stationary distribution of the SGD.
As another example, we run SGD with constant step-size in order to train an ordinary linear regression model for a dataset of n = 5 samples and d = 2 dimensions, i.e., a>i w ⇡ yi, where for i = 1, . . . , 5, yi and each coordinate of ai are drawn uniformly at random from the interval [ 1, 1]. Figure 2 shows the heatmap of the resulting stationary distributions for different step-size choices ⌘ ranging from 0.1 to 0.9 (bright colors represent higher density). We observe that for small step-size choices, the stationary distribution is dense, whereas a fractal structure can be clearly observed as the step-size gets larger.
Fractals are complex patterns and the level of this complexity is typically measured by the Hausdorff dimension of the fractal, which is a notion of dimension that can take fractional values1, and can be much smaller than the ambient dimension d. Recently, assuming that SGD trajectories can be well-approximated by a certain type of stochastic differential equations (SDE) [ŞGN+19, ŞSG19, NcGR19, ŞZTG20], it is shown that the generalization error can be controlled by the Hausdorff dimension of the trajectories of the SDE, instead of their ambient dimension d [ŞSDE20]. That is, the ambient dimension that appears in classical learning theory bounds is replaced with the Hausdorff dimension.
The fractal geometric approach presented in [ŞSDE20] can capture the ‘low dimensional structure’ of fractal sets and provides an alternative perspective to the compression-based approaches that aim to understand why overparametrized networks do not overfit [AGNZ18, SAN20, SAM+20, HJTW21]. However, SDE approximations for SGD often serve as mere heuristics, and guaranteeing a good approximation typically requires unrealistically small step-sizes [LTE19]. For more realistic step-
sizes, theoretical concerns have been raised about the validity of conventional SDE approximations for SGD [LMA21, GŞZ21, Yai19]. Another drawback of the SDE approximation is that the bounds in [ŞSDE20] are implicit, in the sense that they cannot be related to algorithm hyperparameters, problem geometry, or data.
We address these issues and present a direct, discrete-time analysis by exploiting the connections between IFSs and stochastic optimization algorithms. Our contributions are summarized as follows:
• We extend [ŞSDE20] and show that the generalization error can be linked to the Hausdorff dimension of invariant measures (rather than the Hausdorff dimension of sets as in [ŞSDE20]). More precisely, under appropriate conditions, we establish a generalization bound for the stationary distribution of IFS w1 ⇠ µ. That is, with probability at least 1 2⇣,
|R̂(w1,Sn) R(w1)| .
s dimHµ log 2(n)
n +
log(1/⇣)
n , (4)
for n large enough, where dimHµ is the (upper) Hausdorff dimension of the measure µ. • By leveraging results from IFS theory, we further link dimHµ to (i) the form of the recursion (e.g.,
h⌦k in (3)), (ii) algorithm hyperparameters (e.g., ⌘, b), and (iii) problem geometry (e.g., Hessian of R̃k), through a single term, which encapsulates all these components and their interaction.
• We establish bounds on dimHµ for SGD and preconditioned SGD algorithms, when used to minimize various empirical risk minimization problems such as least squares, logistic regression, support vector machines. In all cases, we explicitly link the generalization performance of the model to the hyperparameters of the underlying training algorithm.
• Finally, we numerically compute key quantities that appear in our generalization bounds, and show empirically that they have a statistically significant correlation with the generalization error.
Notation and preliminaries. Bd(x, r) ⇢ Rd denotes the closed ball centered around x 2 Rd with radius r. A function f : Rd1 ! Rd2 is said to be (Fréchet) differentiable at x 2 Rd1 if there exists a d1⇥d2 matrix Jf (x) : Rd1 ! Rd2 such that limkhk!0 kf(x+h) f(x) Jf (x)hk/khk = 0. The matrix Jf (x) is called the differential of f , also known as the Jacobian matrix at x, and determinant of Jf (x) is called the Jacobian determinant [HS74]. For real-valued functions f, g, we define f(n) = !(g(n)) if limn!1 |f(n)|/g(n) = 1. For a set A, |A| denotes its cardinality. For a scalar-valued function f̃ : R ! R, we define kf̃k1 = maxr2R |f̃(r)|.
1The Hausdorff dimension of the middle-third Cantor set in Figure 1 is log3(2) ⇡ 0.63 whereas the ambient dimension is 1 [Fal04, Example 2.3].

heading: 2 Technical Background on Fractal Geometry
text: Fractal sets emerge virtually in all branches of science, and fractal-based techniques have been used in machine learning [SHTY13, MSS19, DSD20, ŞSDE20, AGZ21]. The inherent ‘complexity’ of a fractal set often plays an important role and it is typically measured by its fractal dimension, where several notions of dimension have been proposed [Fal04]. In this section, we briefly mention two important notions of fractal dimension, which will be used in our theoretical development.
Minkowski dimension of a set. The Minkowski dimension (also known as the box-counting dimension [Fal04]) is defined as follows. Let F ⇢ Rd be a set and for > 0, let N (F ) denote a collection of sets that contains the smallest number of closed balls of diameter at most which cover F . Then the upper-Minkowski dimension of F is defined as follows:
dimMF := lim sup !0
h log |N (F )| / log(1/ ) i . (5)
To visualize the upper-Minkowski dimension of a set F , consider the set F lying on an evenly spaced grid and count how many boxes are required to cover the set. The upper-Minkowski dimension measures how this number changes as the grid is made finer using a box-counting algorithm.
Hausdorff dimension of a set. An alternative to the purely geometric Minkowski dimension, the Hausdorff dimension [Hau18] is a measure theoretical notion of fractal dimension. It is based on the Hausdorff measure, which generalizes the traditional notions of area and volume to noninteger dimensions [Rog98]. More precisely, for s 0, let F ⇢ Rd and > 0, and denote H s (F ) := inf P1 i=1 diam(Ai)
s, where the infimum is taken over all the -coverings {Ai}i of F , that is, F ⇢ [iAi with diam(Ai) < for every i. The s-dimensional Hausdorff measure of F is defined as the monotone limit Hs(F ) := lim !0 Hs (F ). When s 2 N, Hs is the s-dimensional Lebesgue measure up to a constant factor; hence the generalization of ‘volume’ to fractional orders.
Based on the Hausdorff measure, the Hausdorff dimension of a set F ⇢ Rd is then defined as follows: dimH F := sup {s > 0 : H s(F ) > 0} = inf {s > 0 : Hs(F ) < 1} .
In other words, the Hausdorff dimension of F is the ‘moment’ s when Hs(F ) drops from 1 to 0, that is, Hr(F ) = 0 for all r > dimH F and Hr(F ) = 1 for all r < dimH F .
We always have 0  dimH F  d, and when F is bounded, we always have 0  dimH F  dimMF  d [Fal04]. Furthermore, the Hausdorff dimension of Rd equals d, and the Hausdorff dimension of smooth Riemannian manifolds correspond to their intrinsic dimension, e.g. dimH Sd 1 = d 1, where Sd 1 is the unit sphere in Rd. Hausdorff dimension of a probability measure. IFSs generate invariant measures as the number of iterates goes to infinity, and random fractals arise from such invariant measures. There has been a growing literature that studies the structure of such random fractals [Saz00, NSB02, MS02, Ram06, FST06, JR08], where the notion of fractal dimension can be extended to measures, and our theory will rely on the Hausdorff dimension of invariant measures associated with stochastic optimization algorithms. In particular, we will mainly use the upper Hausdorff dimension dimHµ of a Borel probability measure µ on Rd, which is defined as follows: dimHµ := inf {dimH A : µ(A) = 1}. In other words, dimHµ is the smallest Hausdorff dimension of all measurable sets with full measure.

heading: 3 Generalization Bounds for Stochastic Optimization Algorithms as IFSs
text: In this section, we will present our main theoretical results which relate the generalization error to the upper-Hausdorff dimension of the invariant measure associated with a stochastic optimization algorithm. We consider a standard supervised learning setting, where Z = X ⇥ Y , where X is the space of features and Y is the space of labels, and ⇡ is the unknown data distribution on Z .
For mathematical convenience, in order to construct the training set with n elements, we first consider an infinite sequence of i.i.d. data samples from the data distribution ⇡, then we take the first n elements from this infinite sequence. More precisely, we consider the (countable) product measure ⇡ 1 = (⇡ ⌦ ⇡ ⌦ . . . ) defined on the cylindrical sigma-algebra. Then, we consider the infinite i.i.d. data sequence as S ⇠ ⇡1, i.e., S = (zj)j 1 with zj i.i.d. ⇠ ⇡ for all j = 1, 2, . . . . Finally, we define the training set as Sn := (z1, . . . , zn), i.e., we take the first n elements of S. To avoid technical
complications, throughout the paper we will assume that all the encountered functions and sets are measurable. All the proofs are given in the supplement.
Given a dataset Sn, we represent the training algorithm as an IFS, which is based on the following recursion: wk = h⌦k(wk 1;Sn), where the mini-batch ⌦k is i.i.d. sampled according to some distribution (e.g., sampling without-replacement uniformly among all possible mini-batches). This compact representation enables us to cover a broad range of optimization algorithms with a unified notation, including SGD (see (3)), as well as preconditioned SGD, and stochastic Newton methods. For example, if we take h⌦k(w;Sn) = w ⌘Hk(w) 1rR̃k(w) where Hk(w) is an estimate of the Hessian of R̃k, we cover stochastic Newton methods [EM15]. Similar constructions can be made for other popular algorithms, such as SGD-momentum [Qia99], RMSProp [HSS12], or Adam [KB15].
Notice that there are only finitely many values that ⌦k can take. For example, in the case of withoutreplacement mini-batch sampling with batch-size b, there are in total mb = m wo-replacement b := n b many subsets of {1, 2, . . . , n} with cardinality b. Alternatively, another setup would be to divide the dataset into mb = mbatchb := n/b batches with each batch having b elements, and at each iteration k, we can randomly choose one of the batches. In both examples we can enumerate as S1, S2, . . . , Smb . If the probability of sampling the mini-batch ⌦k = Si is pi for every i, then, with a slight abuse of notation, we can rewrite the IFS recursion as:
wk = hUk(wk 1;Sn), (6) where Uk is a random variable taking values in {1, 2, . . . , mb} and pi := P(Uk = i). If the minibatch sampling is uniform (i.e., the default option in practice), we have pi = 1/mb; however, we are not restricted to this option, the sampling scheme is allowed to be more general. We finally call the triple (Rd, {hi(·;Sn)}mbi=1, {pi} mb i=1) an iterated function system (IFS).
Given a dataset Sn, we are interested in the limiting behavior of the training algorithm (6). We characterize this behavior by considering the invariant measure µW |Sn of the IFS (also called stationary distribution), that is a Borel probability measure on Rd, such that w1 ⇠ µW |Sn . To be able to work in this context, we first need to ensure that the recursion (6) admits an invariant measure, i.e., µW |Sn exists. Accordingly, we require the following mild conditions on the IFS (6). Let U be a random variable with the same distribution as Uk. If the recursion (6) is Lipschitz on average, i.e.,
E[LU | Sn] < 1, with LU := sup x,y2Rd
khU (x;Sn) hU (y;Sn)k
kx yk , (7)
and is contractive on average, i.e., if
E [log(LU ) | Sn] = Xmb
i=1 pi log(Li) < 0, with pi > 0, for any i = 1, . . . , mb, (8)
then it can be shown that the process is ergodic and admits a unique invariant measure where the limit ⇢ := limk!1(1/k) log hUkhUk 1 · · ·hU1 exists almost surely and is a constant [Elt90], where ⇢ is called the Lyapunov exponent. Furthermore, under further technical assumptions, it can be shown that (6) is geometrically ergodic [DF99]. We note that this condition for the existence of the invariant measure is only applicable to stochastic optimization algorithms with a constant stepsize in which case the random map hU is not time-varying. If decaying stepsize is used instead, then the limit may degenerate to be a singleton. For example, in the toy example illustrated in Figure 1 with quadratics in dimension one, if we use SGD with decaying stepsize ⌘k = c/k where the positive constant c is small enough, then the limit of the iterates will be a singleton as the iterates will converge to the global minimum of the optimization objective (see e.g. [GOP21, GOP19]).
Our goal will be to relate the generalization error to dimHµW |Sn . To achieve this goal, at first sight, it might seem tempting to extract a full-measure set A by using the definition of µW |Sn , such that µW |Sn(A) = 1 and dimH A ⇡ dimHµW |Sn , and then directly invoke the results from [ŞSDE20], which would link the generalization error to dimH A, hence, also to dimHµW |Sn . However, since [ŞSDE20] does not consider an IFS framework, the conditions they require (e.g., boundedness of A, dimM A = dimH A) are not suited to IFSs, and hence prevent us from directly using their results.
As a remedy, we make a detour and show that we can find almost full-measure sets A, such that µW |Sn(A) ⇡ 1 and dimMA ⇡ dimHµW |Sn (notice that in this case we directly use the Minkowski dimension of A, as opposed to its Hausdorff dimension). To achieve this goal, we require the following geometric regularity condition on the invariant measure.
H 1. For ⇡1-almost all S and all n 2 N+, the recursion (6) satisfies (7) and (8) and the limit limr!0 h log µW |Sn(Bd(w, r)) / log r i exists for µW |Sn -almost every w.
This is a common condition [Pes08] and is satisfied for a large class of measures. For instance, ‘sufficiently regular’ measures with the property that C1rs  µ(B(x, r))  C2rs for some constant s > 0 and positive constants C1, C2 will satisfy this assumption. Such measures are called Ahlforsregular (cf. [ŞSDE20, Assumption H4] for a related condition), and it is known that IFSs that satisfy certain ‘open set conditions’ lead to Ahlfors regular invariant measures (see [MT10, Section 8.3]). Yet, our assumption is more general and does not immediately require Ahlfors-regularity.
Under H1, we now formalize our key observation, which serves as the basis for our bounds. Proposition 1. Assume that H1 holds. Then for every " > 0, n 2 N+, and ⇡1-almost every S, there exists := (",Sn) 2 (0, 1] and a bounded measurable set ASn, ⇢ Rd, such that
µW |Sn(ASn, ) 1 , and dimMASn,  dimHµW |Sn + ", (9)
and (",Sn) ! 0 as " ! 0.
Thanks to this result, we can now leverage the proof technique presented in [ŞSDE20, Theorem 2], and link the generalization error to dimHµW |Sn through dimMASn, . We shall emphasize that, mainly due to the sets ASn, not being of full-measure, our framework introduces additional non-trivial technical difficulties that we need to tackle in our proof.
We now introduce our second assumption, which roughly corresponds to a ‘topological stability’ condition, and is adapted from [ŞSDE20, Assumption H5]. Formally, consider the (countably infinite) collection of closed balls of radius , whose centers are on the fixed grid N := n⇣ (2j1+1)
2 p d
, . . . , (2jd+1)
2 p d
⌘ : ji 2 Z, i = 1, . . . , d o , and for a set A ⇢ Rd, define N (A) :=
{x 2 N : Bd(x, ) \ A 6= ;}, which is the collection of the centers of the balls that intersect A. H 2. Let Z1 := (Z ⇥ Z ⇥ · · · ) denote the countable product endowed with the product topology and let B be the Borel -algebra generated by Z1. For a Borel set A ⇢ Rd, let F,G be the sub- - algebras of B generated by the collections of random variables given by {R̂(w,Sn) : w 2 Rd, n 1} and n 1 {w 2 N (ASn, )} , µW |Sn(ASn, ), dimHµW |Sn : , 2 Q>0, w 2 N , n 1 o , where
ASn, is given in Proposition 1. There exists a constant M 1 such that for any F 2 F, G 2 G, we have P (F \ G)  MP (F )P(G).
H2 simply ensures that the dependence between the training error and the topological properties of the support of µW |Sn can be controlled via M . Hence, it can be seen as a form of algorithmic stability [BE02], where M measures the level of stability of the topology of µW |Sn : a small M indicates that the geometrical structure of µW |Sn does not heavily depend on the particular value of Sn. The constant M is also related to the mutual information [XR17, AAV18, HŞKM21], but may be better behaved than the mutual information as it relies on very specific functions of the random variables. Similar to mutual information, a-priori there is no reason to expect M to be finite for general algorithms; intuitively, however, the more stochasticity an algorithm incorporates the more we expect the set ASn, and the loss landscape to decouple. For example, for a purely random algorithm (which ignores Sn) the two objects will be independent. In the other extreme, where the algorithm is deterministic given the sample may fail to be finite. Since we are controlling the generalization error on the support, which is itself a random set depending on the sample, we require H2 to be able to make progress.
We require one final assumption, which states that the loss ` is sub-exponential. H 3. ` is L-Lipschitz continuous in w, and when z ⇠ ⇡, for all w, `(w, z) is (⌫, )-sub-exponential, that is, for all | | < 1/, we have logEz⇠⇡ [exp ( (`(w, z) R(w)))]  ⌫2 2/.
Armed with these assumptions, we can now present our main result. Theorem 1. Assume that H1 to 3 hold and dimHµW |Sn = !(log log(n)/ log(n)), ⇡1-almost-surely. Then, the following bound holds for sufficiently large n:
|R̂(W,Sn) R(W )|  8⌫
s dimHµW |Sn log 2(nL2)
n +
log(13M/⇣)
n , (10)
with probability at least 1 2⇣ over the joint distribution of S ⇠ ⇡1, W ⇠ µW |Sn .
This theorem shows that the Hausdorff dimension of the invariant measure acts as a ‘capacity metric’ and the generalization error is therefore directly linked to this metric, i.e., the complexity of the underlying fractal structure has close links to the generalization performance. On the other hand, the condition dimHµW |Sn = !(log log(n)/ log(n)) is very mild and makes sure that the dimension of the IFS does not decrease very rapidly with increasing number of data points n. We shall mention that Theorem 1 has an asymptotic nature as we do not have an explicit control on how large n should be. This is due to the fact that the notions of Minkowski and Hausdorff dimensions are essentially asymptotic, which unfortunately prevents us from obtaining any truly non-asymptotic result. However, obtaining nonasymptotic results might be possible with further assumptions on the fractal dimension of the invariant measures and their supports.
Theorem 1 enables us to access the rich theory of IFSs, where bounds on the Hausdorff dimension are readily available, and connect them to statistical learning theory. The following result is a direct corollary to Theorem 1 and [Ram06, Theorem 2.1] (see Theorem S2 in the supplementary document). Corollary 1. Assume that the conditions of Theorem 1 hold. Furthermore, consider the recursion (6) and assume that hi are continuously differentiable with derivatives Jhi that are ↵-Hölder continuous for some ↵ > 0. Then, there exists a constant M > 1 such that for sufficiently large n:
|R̂(W,Sn) R(W )|  8⌫ s E log2(nL2)⇥Pmb
i=1 pi R Rd log(kJhi(w)k)dµW |Sn(w) ⇤ n
+ log(13M/⇣)
n , (11)
with probability 1 2⇣ over S ⇠ ⇡1, W ⇠ µW |Sn , where E := Pmb
i=1 pi log(pi) denotes the negative entropy of the mini-batch sampling scheme, k · k denotes the operator norm (with the `2-norm being the underlying norm), and Jhi is the Jacobian of hi defined in the notation section.
Note that the Hölder condition is mainly used to ensure that the invariant measure exists and the constant ↵ does not directly interact with the bound. However, it might affect the rate of convergence to the invariant measure. By this result, we discover an interesting quantity, Pmb
i=1 pi R Rd log(kJhi(w)k)dµW |Sn(w), which
simultaneously captures the effects of the data and the algorithm2. To see it more clearly, let us consider the SGD recursion (3), where kJhi(w)k = kI ⌘r2R̃Si(w)k and {Si} mb i=1 denotes the enumeration of the mini-batches. Then, the overall quantity becomes
EU,W ⇥ log kI ⌘r2R̃SU (W )k ⇤ , (12)
where the expectation is taken over the mini-batch index U 2 {1, . . . , mb} with P(U = i) = pi, and W ⇠ µW |Sn . We clearly observe that this term depends on (i) the algorithm choice through the form of hi, (ii) step-size ⌘, (iii) batch-size through mb, (iv) problem geometry through r2R̃, and (v) data distribution through µW |Sn . We believe that such a compact representation of all these constituents and their interaction is novel and will open up interesting future directions.

heading: 4 Analytical Estimates for the Hausdorff Dimension
text: The generalization bound presented in Theorem 1 applies to a number of stochastic optimization algorithms that can be represented with an IFS and to a large class of losses that can be non-convex or convex. It is controlled by the Hausdorff dimension of the invariant measure µW |Sn which needs to be estimated. In the numerical experiments section, we will discuss how this quantity can be estimated from the dataset Sn and the iterates of the underlying algorithm.
Corollary 1 shows that for smooth losses, the Hausdorff dimension can be controlled with the expectation of the norm of the logarithm of the Jacobian log(kJhi(w)k) with respect to the invariant measure µW |Sn . In general, an explicit characterization of the invariant measure is not known. Nevertheless, under additional appropriate assumptions that can hold in practice, such as boundedness of the data of the loss, we next discuss that it is possible to get uniform lower and upper bounds on the quantity kJhi(w)k which leads to analytical upper bounds on dimHµW |Sn .
2Note that thanks to [Ram06], we can allow state-dependent pi = pi(w); yet, we do not consider this option as its application is not immediately clear.
As illustrative examples; in the following, we will consider the setting where we divide Sn into mb = mbatchb = n/b batches with each batch having b elements, and then we discuss how analytical estimates on the (upper) Hausdorff dimension dimHµW |Sn can be obtained for some particular problems including least squares, regularized logistic regression, and one hidden-layer networks. In the supplementary document, we also discuss how similar bounds can be obtained for support vector machines and other algorithms such as preconditioned SGD and stochastic Newton methods.
Least squares. We consider the least squares problem, with data points zi = (ai, yi) and loss
`(w, zi) := a T i w yi 2 /2 + kwk2/2, (13)
where > 0 is a regularization parameter. Proposition 2 (Least squares). Consider the least squares problem (13). Assume the step-size ⌘ 2 (0, 1R2+ ), where R := maxi kaik is finite. Then, we have the following upper bound:
dimHµW |Sn  log (n/b)
log(1/(1 ⌘ )) . (14)
Note that here ` is only pseudo-Lipschitz |`(w, zi) `(w0, zi)|  Li(1 + kwk+ kw0k)kw w0k for some Li > 0, rather than globally Lipschitz. However; the conditions in Proposition 2 ensure that w will stay in a bounded region, in which case ` becomes Lipschitz. Also note that only the logarithm of the Lipschitz constant directly enters the bound.
We observe that, for fixed n, the upper bound for dimHµW |Sn is decreasing both in ⌘ and b. This behavior is not surprising: large ⌘ results in chaotic behaviors (cf. Figures 1,2), and in the extreme case where b = n, the algorithm becomes deterministic and hence converges to a single point, in which case the Hausdorff dimension becomes 0. However, the decrement due to b does not automatically grant good generalization performance: since the algorithm becomes deterministic, the stability constant M in H2 can get arbitrarily large, hence the bound in Theorem 1 could become vacuous. This outcome reveals an interesting tradeoff between the Hausdorff dimension and the constant M , through the batch-size b, and investigating this tradeoff is an interesting future direction.
We further notice that the numerator in (14) is logarithmically increasing with n, which is compensated by the factor 1/n in Theorem 1. Nevertheless, we can take the batch-size in proportion with n (i.e., setting mb to a constant value), in order avoid this logarithmic growth. We also note that the input dimension d potentially affects the term R, which forms the bound for the input data, and hence the input dimension will indirectly affect the generalization bound. Finally, regarding the remaining bounds in this section, even though their forms might differ from (14), similar remarks also apply. Hence, we will omit the discussion.
Regularized logistic regression. Given the data points zi = (ai, yi), consider regularized logistic regression with the loss:
`(w, zi) := log 1 + exp yia T i w + kwk2/2, (15)
where > 0 is a regularization parameter. We have the following result. Proposition 3 (Regularized logistic regression). Consider the regularized logistic regression (15). Assume the step-size ⌘ < 1/ and the input data is bounded, i.e. R := maxi kaik < 2 p . We have:
dimHµW |Sn  log (n/b)
log(1/(1 ⌘ + 1 4 ⌘R2))
. (16)
Next, we consider a non-convex formulation for robust regression (see e.g. [MBM18]), with zi = (ai, yi) and the loss
`(w, zi) := ⇢ (yi hw, aii) + kwk 2 /2, (17)
where > 0 is a regularization parameter and ⇢ is a non-convex function, assumed to be twice continously differentiable, where a standard choice is Tukey’s bisquare loss defined as ⇢Tukey(t) = 1 (1 (t/t0)2)3 for |t|  t0, and ⇢Tukey(t) = 1 for |t| t0 (see e.g. [MBM18]), and exponential squared loss: ⇢exp(t) = 1 e |t| 2/t0 , where t0 > 0 is a tuning parameter (see e.g. [WJHZ13]).
Proposition 4 (Non-convex formulation for robust regression). Consider the non-convex formulation for robust regression (17). Assume the step-size ⌘ < 1 +R2(2/t0) , where R = maxi kaik < p t0/2. Then, we have the following upper bound for the Hausdorff dimension:
dimHµW |Sn  log (n/b)
log(1/(1 ⌘ + ⌘R2 2t0 )) . (18)
One hidden-layer neural network. Given the data points zi = (ai, yi). Let ai 2 Rd be the input and yi 2 Rm be the corresponding output, and let wr 2 Rd be the weights of the r-th hidden neuron of a one hidden-layer network, and br 2 R is the output weight of hidden unit r. For simplicity of the presentation, following [ZMG19, DZPS19], we only optimize the weights of the hidden layer, i.e. w = ⇥ w
T 1 w T 2 . . . w T m
⇤ is the decision variable with the regularized squared loss:
`(w, zi) := kyi ŷik 2 + kwk2/2, ŷi := Xm r=1 br w T r ai , (19)
where the non-linearity : R ! R is smooth and > 0 is a regularization parameter. Proposition 5 (One hidden-layer network). Consider the one hidden-layer network (19). Assume the step-size ⌘ < 1
2 . Then, we have the following upper bound for the Hausdorff dimension:
dimHµW |Sn  log (n/b)
log(1/(1 ⌘( C))) , (20)
where C := Mykbk1k 00k1R2+(maxj kvjk1) 2
< , where R := maxi kaik, My := maxi kyi ŷik, and vi := ⇥ b1 0(wT 1 ai)ai b2 0(wT2 ai)ai · · · bm 0(wTmai)ai ⇤T .
In Proposition 4, we assumed that C is uniformly bounded and > C for mathematical convenience. However, in general, C by its definition might increase in m and n so that we do not expect that we can choose small for large data sets and wide networks.

heading: 5 Experiments
text: Our aim now is to empirically investigate whether the bound in Corollary 1 is informative, in that it is predictive of a neural network’s generalization error. As the second term of this bound cannot be evaluated, we will assume that it is negligible and focus our efforts on the first term. Further, because the denominator of the first term is the only term that depends on the invariant measure, we want to establish that the inverse of this denominator is predictive of a neural network’s generalization error. Note however that for complex models, such as modern neural networks, analytically bounding kJhi(x)k becomes highly non-trivial.
In our experiments, we fix the algorithm to SGD and we develop a numerical method for computing the term (12). Noting that Jhi(w) = I ⌘r2R̃i(w), for simplicity we denote the inverse of (12) as the ‘complexity’: R = 1/ ⇥Pmb i=1 pi R Rd log (kJhi(w)k) dµW |Sn(w) ⇤ 3.
To approximate the expectations, we propose the following simple Monte Carlo strategy:
R 1 ⇡ h 1/(NWNU ) iXNW i=1 XNU j=1 log ⇣ kJhUj (Wi)k ⌘ , (21)
where Uj denotes i.i.d. random mini-batch indices that are drawn without-replacement from {1, . . . , n} and Wi i.i.d. ⇠ µW |S . Assuming (8) is ergodic [DF99], we treat the iterates wk as i.i.d. samples from µW |Sn for large k, hence, log(kJhUj (Wi)k) can be computed on these iterates, and (21) can be computed accordingly. Our implementation for computing kJhUj (Wi)k for neural networks with millions of parameters is detailed in the supplementary document. Though the size of Jhi is very large in our experiments (⇡ 20M⇥20M on average), our algorithm can efficiently compute the norms without constructing JhUj (Wi), by extending the approach presented in [YGKM20]. Our implementation is available at https://github.com/umutsimsekli/fractal_generalization.
3Note that the first term of the bound in Corollary 1 suggests computing p R rather than R; however, both
choices yield very similar results with high statistical significance.
In Figure 3 we plot the estimates of R for a variety of convolutional (CONV) and fully connected network (FCN) architectures trained on CIFAR10, SVHN and Boston House Prices (BHP). For the full details of the models, the hardware used, their run-time, and the convergence criterion used, see Section S7 in the supplement. The plot demonstrates that R and generalization error are positively correlated and that this correlation is significant (p-value ⌧ 0.001) for all model architectures. This provides evidence that the bound on the generalization error in Corollary 1 is informative.
To support our findings in Section 4 that the bound for the Hausdorff dimension dimHµW |Sn is monotonically decreasing in the step-size ⌘, we plot R against ⌘ in Figure 4 for the networks trained on BHP in Figure 3. R decreases with increasing ⌘, clearly backing our findings.
We note that these results were inconclusive for classification models trained with a cross-entropy loss, in that we could not clearly observe a negative or positive correlation. Future work will further study this lack of correlation, particular to classification models.

heading: 6 Conclusion
text: In this work, we investigated stochastic optimization algorithms through the lens of IFSs and studied their generalization properties. Under some assumptions, we showed that the generalization error can be controlled based on the Hausdorff dimension of the invariant measure determined by the iterations, which can lead to tighter bounds than the standard bounds based on the ambient dimension. We proposed an efficient methodology to estimate the Hausdorff dimension in deep learning settings and supported our theory with several experiments on neural networks. We also illustrated our bounds on specific problems and algorithms such as SGD and its preconditioned variants, which unveil new links between generalization, algorithm parameters and the Hessian of the loss.
Our study does not have a direct societal impact as it is largely theoretical. The limitation of our study is its asymptotic nature due to operating on invariant measures. Future work will address (i) obtaining nonasymptotic bounds in terms of the number of iterations k, (ii) including the term of (12) as a regularizer to the optimization problem, which would be an alternative to the methods that aim to “decrease the intrinsic dimension” [ZQH+18, BLGŞ21].
Acknowledgements. MAE is partially funded by CIFAR AI Chairs program, and CIFAR AI Catalyst grant, NSERC Grant [2019-06167]. MG’s research is supported in part by the grants Office of Naval Research Award Number N00014-21-1-2244, National Science Foundation (NSF) CCF-1814888, NSF DMS-2053485, NSF DMS-1723085. UŞ’s research is supported by the French government under management of Agence Nationale de la Recherche as part of the “Investissements d’avenir” program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute). LZ is grateful to the partial support from NSF DMS-2053454 and a Simons Foundation Collaboration Grant.

Review summaries:
The paper gives novel generalization bounds for models trained with SGD (or variants of SGD) based on the Hausdorff dimension. The main contribution is Theorem 1 which shows, under certain assumptions, that for a general class of models, their generalization bound after training with SGD can be bounded by a term that depends on the Hausdorff dimension of the stationary distribution which the stochastic process converges to. Next, several interpretations of this bound are given for simple models like least squares, and 2-layer NN. Finally, an experiment shows on several datasets that there is a correlation between the generalization capacity of NN and a certain term in the given bound.

The authors give a general generalization bound for stochastic optimization algorithms, depending on the Hausdorff dimension of the invariant measure. This bound takes into account the dataset and training algorithm, and draws on theory on fractal dimensions and random iterated function systems in dynamical systems. The authors link this bound to quantities that can be estimated from running the algorithm in question - namely an expected logarithm of norm of Jacobian. They specialize the theory for least squares (SGD & stochastic Newton), logistic regression, SVMs, and 1-hidden-layer neural networks. They present experiments showing that this quantity does in fact correlate with with generalization performance for fully-connected and convolutional neural networks trained on CIFAR10 and SVHN.

The authors study the generalization error of iterative stochastic optimization algorithms whose iterates form an ergodic Markov chain. Particular examples include SGD and stochastic Newton methods for a large class of loss functions and models. The generalization error of these algorithms is upper-bounded by the Hausdorff of the invariant measure of the associated Markov chain. Bounds on this dimension are then given for some common example models and loss functions. Statistical correlation between the generalization error of deep neural networks and a quantity in the derived bounds is established numerically, but is only proven for a single layer model.