A Spin Glass Model for the Loss Surfaces of Generative Adversarial Networks.

We present a novel mathematical model that seeks to capture the key design feature of generative adversarial networks (GANs). Our model consists of two interacting spin glasses, and we conduct an extensive theoretical analysis of the complexity of the model's critical points using techniques from Random Matrix Theory. The result is insights into the loss surfaces of large GANs that build upon prior insights for simpler networks, but also reveal new structure unique to this setting which explains the greater difficulty of training GANs. © Crown 2022.

Introduction

By making various modeling assumptions about standard multi-layer perceptron neural networks, [1] argued heuristically that the training loss surfaces of large networks could be modelled by a spherical multi-spin glass. Using theoretical results of [2], they were able to arrive at quantitative asymptotic characterisations, in particular the existence of a favourable ‘banded structure’ of local-optima of the loss. There are clear and acknowledged deficiencies with their assumptions [3] and recent observations have shown that the Hessians of real-world deep neural networks do not behave like random matrices from the Gaussian Orthogonal Ensemble (GOE) of Random Matrix Theory at the macroscopic scale [4-6], despite this being implied by the spin-glass model of [1]. Moreover, there have been questions raised about whether the mean asymptotic properties of loss surfaces for deep neural networks (or energy surfaces of glassy objects) are even relevant practically for gradient-based optimisation in sub-exponential time [7-9], though interpretation of experiments with deep neural networks remains difficult and the discussion about the true shape of their loss surfaces and the implications thereof is far from settled. Nevertheless, spin-glass models present a tractable example of high-dimensional complex random functions that may well provide insights into aspects of deep learning. Rather than trying to improve or reduce the assumptions of [1], various authors have recently opted to skip the direct derivation from a neural network to a statistical physics model, instead proposing simple models designed to capture aspects of training dynamics and studying those directly. Examples include: the modified spin glass model of [10] with some explicitly added ‘signal’; the simple explicitly non-linear model of [11]; the spiked tensor ‘signal-in-noise’ model of [12]. In a slightly different direction, [13] removed one of the main assumptions from the [1] derivation, and in so doing arrived at a deformed spin-glass model. All of this recent activity sits in the context of earlier work connecting spin-glass objects with simple neural networks [14-16] and, more generally, with image reconstruction and other signal processing problems [17].

One area that has not been much explored in the line of the above-mentioned literature is the study of architectural variants. Modern deep learning contains a very large variety of different design choices in network architecture, such as convolutional networks for image and text data (among others) [18, 19], recurrent networks for sequence data [20] and self-attention transformer networks for natural language [21, 22]. Given the ubiquity of convolutional networks, one might seek to study those, presumably requiring consideration of local correlations in data. One could imagine some study of architectural quirks such as residual connections [23], and batch-norm has been considered to some extent by [24]. In this work, we propose a novel model for generative adversarial networks (GANs) [25] as two interacting spherical spin glasses. GANs have been the focus of intense research and development in recent years, with a large number of variants being proposed [26-32] and rapid progress particularly in the field of image generation. From the perspective of optimisation, GANs have much in common with other deep neural networks, being complicated high-dimensional functions optimised using local gradient-based methods such as stochastic gradient descent and variants. On the other hand, the adversarial training objective of GANs, with two deep networks competing, is clearly an important distinguishing feature, and GANs are known to be more challenging to train than single deep networks. Our objective is to capture the essential adversarial aspect of GANs in a tractable model of high-dimensional random complexity which, though being a significant simplification, has established connections to neural networks and high dimensional statistics.

Our model is inspired by [1, 12, 33, 34] with spherical multi-spin glasses being used in place of deep neural networks. We thus provide a complicated, random, high-dimensional model with the essential feature of GANs clearly reflected in its construction. By employing standard Kac-Rice complexity calculations [2, 35, 36] we are able to reduce the loss landscape complexity calculation to a random matrix theoretic calculation. We then employ various Random Matrix Theory techniques as in [13] to obtain rigorous, explicit leading order asymptotic results. Our calculations rely on the supersymmetric method in Random Matrix Theory, in particular the approach to calculating limiting spectral densities follows [37] and the calculation also follows [38, 39] in important ways. The greater complexity of the random matrix spectra encountered present some challenges over previous such calculations, which we overcome with a combination of analytical and numerical approaches. Using our complexity results, we are able to draw qualitative implications about GAN loss surfaces analogous to those of [1] and also investigate the effect of a few key design parameters included in the GAN. We compare the effect of these parameters on our spin glass model and also on the results of experiments training real GANs. Our calculations include some novel details, in particular, we use precise sub-leading terms for a limiting spectral density obtained from supersymmetric methods to prove a required concentration result to justify the use of the Coulomb gas approximation. We note that our complexity results could be also be obtained in principle using the methods developed in [40], however our work was completed several months before this pre-print appeared. Our approach for computing the limiting spectral density may nevertheless be the simplest and would be used as input to the results of [40].

The role that statistical physics models such as spherical multi-spin glasses are to ultimately play in the theory of deep learning is not yet clear, with arguments both for and against their usefulness and applicability. We provide a first attempt to model an important architectural feature of modern deep neural networks within the framework of spin glass models and provide a detailed analysis of properties of the resulting loss (energy) surface. Our analysis reveals potential explanations for observed properties of GANs and demonstrates that it may be possible to inform practical hyperparameter choices using models such as ours. Much of the advancement in practical deep learning has come from innovation in network architecture, so if deep learning theory based on simplified physics models like spin-glasses is to keep pace with practical advances in the field, then it will be necessary to account for architectural details within such models. Our work is a first step in that direction and the mathematical techniques used may prove more widely valuable.

The paper is structured as follows: in Sect. 2 we introduce the interacting spin glass model; in Sect. 3 we use a Kac-Rice formula to derive random matrix expressions for the asymptotic complexity of our model; in Sect. 4 we derive the limiting spectral density of the relevant random matrix ensemble; in Sect. 5 we use the Coulomb gas approximation to compute the asymptotic complexity, and legitimise its use by proving a concentration result; in Sect. 6 we derive some implications of our model for GAN training and compare to experimental results from real GANs; in Sect. 7 we conclude. All code used for numerical calculations of our model, training real GANs, analysing the results and generating plots is made available1.

An Interacting Spin Glass Model

We use multi-spin glasses in high dimensions as a toy model for neural network loss surfaces without any further justification, beyond that found in [1, 13]. GANs are composed of two networks: generator (G) and discriminator (D). G is a map and D is a map . G’s purpose is to generate synthetic data samples by transforming random input noise, while D’s is to distinguish between real data samples and those generated by G. Given some probability distribution on some , GANs have the following minimax training objectivewhere are the parameters of the discriminator and generator respectively. With , has some probability distribution . When successfully trained, the initially unstructured examples are easily distinguished by D, this in turn drives improvements in G, bring closer to . Ultimately, the process successfully terminates when is very close to and D performs little better than random at the distinguishing task. To construct our model, we introduce two spin glasses:where , all the are i.i.d. and are similarly i.i.d. . We then define the models for the discriminator and generator losses: plays the role of the loss of the discriminator network when trying to classify genuine examples as such. plays the role of loss of the discriminator when applied to samples produced by the generator, hence the sign difference between and . are the weights of the discriminator, and the weights of the generator. The are surrogates for the training data (i.e. samples from ) and the are surrogates for the noise distribution of the generator. For convenience, we have chosen to pull the scale outside of the and include it as a constant multiplier in (4)–(5). In reality, we should like to keep as i.i.d. but take to have some other more interesting distribution, e.g. normally or uniformly distributed on some manifold. Using [x] to denote the integer part of x, we take for fixed , , and study the regime . Note that there is no need to distinguish between and in the limit.

Remark 1

Our model is not supposed to have any direct relationship to GANs. Rather, we have used two spin glasses as models for high-dimensional random surfaces. The spin glasses are related by sharing some of their variables, namely the , just as the two training objectives in GANs share the discriminator weights. In prior work modeling neural network loss surfaces as spin glasses, the number of spins corresponds to the number of layers in the network, therefore we have chosen p spins for and for , corresponding to p layers in the discriminator and q layers in the generator, but the generator is only ever seen in the losses composed with the discriminator. One could make other choices of and to couple the two glasses and we consider one such example in the appendix Sect. 1.

Kac-Rice Formulae for Complexity

Training GANs involves jointly minimising the losses of the discriminator and the generator. Therefore, rather than being interested simply in upper-bounding a single spin-glass and counting its stationary points, the complexity of interest comes from jointly upper bounding both and and counting points where both are stationary. Using to denote the M-sphere2, we define the complexityfor some Borel sets and where denote the Riemannian covariate derivatives on the hyperspheres with respect to the discriminator and generator weights respectively. Note: So that the finer structure of local minima and saddle points can be probed, we also define the corresponding complexity with Hessian index prescriptionwhere i(M) is the index of M (i.e. the number of negative eigenvalues of M). We have chosen to consider the indices of the Hessians and separately, just as we chose to consider separately vanishing derivatives and . We believe this choice best reflects the standard training loop of GANs, where each iteration updates the discriminator and generator parameters in separate steps.

We have chosen to treat the parameters of each network as somewhat separate by placing them on their own hyper-spheres. This reflects the minimax nature of GAN training, where there really are 2 networks being optimised in an adversarial manner rather than one network with some peculiar structure.

We could have taken and required but, as in the previous comment, our choice is more in keeping with the adversarial set-up, with each network seeking to optimize separately its own parameters in spite of the other.

We will only be interested in the case and , for .

To calculate the complexities, we follow the well-trodden route of Kac-Rice formulae as pioneered by [35, 36]. For a fully rigorous treatment, we proceed as in [2, 13].

Lemma 1

and thereforewhere is the joint density of , the density of , and the density of , all implicitly evaluated at .

Proof

Routine application of a theorem of [41]. See appendix Sect. 1.

With Lemma 1 in place, we can now establish the following Kac-Rice expression specialised to our model:

Lemma 2

For GOE matrix M and independent GOE matrix , defineFor , defineDefine the constantwhere the variances areand is the surface area of the N sphere. The expected complexity is then

Define the matrixappearing in the expression for in Lemma 1. Note that takes the place of a Hessian (though it is not symmetric). We begin with the distribution ofNote that the integrand in (14) is jointly spherically symmetric in both and . It is therefore sufficient to consider in the region of a single point on each sphere. We choose the north poles and coordinate bases on both spheres in the region of their north poles. The remaining calculations are routine Gaussian manipulations which appear in the appendix Sect. 1. One findswhere are independent matrices, is an independent matrix and G is an independent Ginibre matrix. Note that the dimensions are and rather and . This is simply because the hypersphere is an dimensional manifold, and similarly .

We can simplify by summing independent Gaussians to obtainwhereand is a GOE matrix independent of and G.

There is an alternative reformulation of that will also be useful. Indeed, because , let us write aswhere is a GOE matrix of size , is an independent GOE matrix of size andIf follows thatNow define the constantsand then we arrive atThe variances of and derive from those of computed in appendix Sect. 1 (see (139), (143)):Similarly the density is found in (155): We have now collected all the inputs required for Lemma 1. The domain of integration B arises from the constraints and and the re-scaled variables (23). This completes the proof.

We will need the asymptotic behaviour of the constant , which we now record in a small lemma.

Lemma 3

As ,

By Stirling’s formulawhere we have used

Limiting Spectral Density of the Hessian

Our intention now is to compute the the expected complexity via the Coulomb gas method. The first step in this calculation is to obtain the limiting spectral density of the random matrixwhere, note, is just a shifted version of H as defined in Lemma 2. Here the upper-left block is of dimension , and the overall dimension is N. Let be the limiting spectral measure of and its density. The supersymmetric method provides a way of calculating the expected Stieltjes transforms of [37]:Recall that a density and its Stieltjes transform are related by the Stieltjes inversion formulaThe function Z(J) can be computed using a supersymmetric representation of the ratio of determinants. Firstly, we recall an elementary result from multivariate calculus, where M is a real matrix:By introducing the notion of Grassmann variables and Berezin integration, we obtain a complimentary expression:Here the are purely algebraic objects defined by the anti-commutation ruleand are separate objects, with the complex conjugation unary operator defined so that and Hermitian conjugation is then defined as usual by The set of variables generate a graded algebra over . Mixed vectors of commuting and anti-commuting variables are called supervectors, and they belong to a vector space called superspace. The integration symbol is defined as a formal algebraic linear operator by the propertiesFunctions of the the Grassmann variables are defined by their formal power series, e.g.where the termination of the series follows from , which is an immediate consequence of (33). From this it is apparent that (34), along with (33), is sufficient to define Berezin integration over arbitrary functions of arbitrary combinations of Grassmann variables. Finally we establish our notation for supersymmetric (or graded) traces of supermatrices. We will encounter supermatrices of the formwhere A, D are square block matrices of commuting variables and B, C are rectangular block matrices of Grassmann variables. In this case, the graded trace is given by . We refer the reader to [42] for a full introduction to supersymmetric methods.

Using the integral results (31), (32) we can then writewhere the measure is is a vector of N complex commuting variables, and are vectors of N Grassmann variables, and we use the [t] notation to denote the splitting of each of the vectors into the first and last components, as seen in [38]:We then split the quadratic form expressions in (36)Taking the GOE averages is now simple [37, 43]:where the supersymmetric matrices are given byIntroducing the tensor notationandwe can compactly writeWe now perform two Hubbard-Stratonovich transformations [37]where and inherit their form from Q, Q[1]with real commuting variables, and Grassmanns; the factor i is introduced to ensure convergence. Integrating out over is now a straightforward Gaussian integral in superspace, givingRecalling the definition of , we haveand so one immediately obtainsTo obtain the limiting spectral density (LSD), or rather its Stieltjes transform, one must find the leading order term in the expansion for (50). This can be done by using the saddle point method on the manifolds. We know that the contents of the exponential must vanish at the saddle point, since the LSD is , so we in fact need only compute at the saddle point. We can diagonalise within the integrand of (50) and absorb the diagonalising graded U(1/1) matrix into . The resulting saddle point equations for the off-diagonal entries of the new (rotated) dummy variable are trivial and immediately give that is also diagonal at the saddle point. The saddle point equations are then(53) and (54) combine to give an explicit expression for :With a view to simplifying the numerical solution of the coming quartic, we define and then a line of manipulation with (54) and (55) givesBy solving (56) numerically for fixed values of , we can obtain the four solutions . These four solution functions arise from choices of branch for and determining the correct branch directly is highly non-trivial. However, for any , at most one of the will lead to a positive LSD, which gives a simple way to compute numerically using (30) and (50):Plots generated using (57) and eigendecompositions of matrices sampled from the distribution of are given in Fig. 1 and show good agreement between the two. Note the three different forms: single component support, two component support and the transition point between the two, according to the various parameters. In these plots, the larger lobes on the left correspond to the upper left block, which is much larger than the lower-right block (since here). One can see this by considering large , for which there must be a body of eigenvalues in the region of owing to the upper left block. Since only features in the upper-left block, not all of the eigenvalues can be located around , and the remainder are found in the other lobe of the density which is around 0 in Fig.  1.

Fig. 1

Example spectra of showing empirical spectra from 100 matrices and the corresponding LSDs computed from (56). Here , , =1 and is varied to give the three different behaviours

The Asymptotic Complexity

In the previous section, we have found the equilibrium measure, , of the ensemble of random matricesThe Coulomb gas approximation gives us a method of computing :We have access to the density of pointwise (in x and ) numerically, and so (59) is a matter of one-dimensional quadrature. Recalling (14), we then havewhereDue to Lemma 3, the constant term has asymptotic formWe then define the desired asand we haveUsing these numerical methods, we obtain the plot of in B and a plot of for some example values, shown in Figs. 2, 3. Numerically obtaining the maximum of on B is not as onerous as it may appear, since grows quadratically in at moderate distances from the origin.

Fig. 2

for . Red lines show the boundary of the integration region B

Fig. 3

and its cross-sections, fixing separately and . Here

We numerically verify the legitimacy of this Coulomb point approximation with Monte Carlo integrationwhere is the j-th eigenvalues of the i-th i.i.d. sample from the distribution of . The results, comparing at for a variety of are show in Fig. 4. Note the strong agreement even at such modest N, however to rigorously substantiate the Coulomb gas approximation in (59), we must prove a concentration result.

Fig. 4

Comparison of (59) and (65), verifying the Coulomb gas approximation numerically. Here . Sampled matrices for MC approximation are dimension , and MC samples have been used

Lemma 4

Let be a sequence of random matrices, where for each Nand , . Let be the empirical spectral measure of and say weakly almost surely. Then for any as .

We begin by establishing an upper bound. Take any , then Take also any , then trivially Overall we have, for any , Thence an application of Hölder’s inequality givesConsidering , we haveand soThe entries of are Gaussians with variance or and all the diagonal and upper diagonal entries are independent. All of these variances are , sowhere the are i.i.d. standard Gaussians for . It follows thatElementary calculations giveand sothus when we take , we have .

Considering , it is sufficient now to showwhere , a continuous and bounded function. For any , we haveThe entries of are Gaussian with variance and so obey a log-Sobolev inequality as required by Theorem 1.5 from [44]. The constant, c, in the inequality is independent of , so we need not compute it exactly. The theorem from [44] then givesWe have shownWe now need to establish a complimentary lower bound to complete the proof. By Jensen’s inequalityfor any . Convergence in law of to and the dominated convergence theorem givefor large enough , because has compact support. It remains to show that the expectation inside the exponent in the second term of (82) converges to zero uniformly in N in the limit .

By (30), it is sufficient to consider , which is computed via (50). Let us define the function so thatHenceforth, are the solution to the saddle point equations (51–54) and are integration variables. Around the saddle pointfor some . We use the notation for and similarly . A superscript asterisk on or any of its derivatives is short hand for evaluation at the saddle point. While the Hessian of may not in general vanish at the saddle point,and so we must go to at least the cubic term in the expansion of around the saddle point, i.e.The bosonic (BB) and fermionic (FF) coordinates do not interact, so we can consider derivatives of as block tensors. Simple differentiation giveswhere follows similarly withBy the saddle point equations (51)–(54) we haveLet . Thenfor any C. Let and then choose and to givewith , , . The expressions for follow identically. We thus haveor perhaps with the the integration ranges reversed depending on the signs of . We havewhere we have definedThis last bound follows from a standard Cauchy rotation of integration contour if any of has vanishing real part. (97) is valid for , but if and , then the preceding calculations are simplified and we still obtain an upper bound but proportional to . Similarly with and and similarly for . The only remaining cases are or . But recall (93) and (53)–(54). We immediately see that if and only if , which occurs for no finite . Therefore, for fixed , and any where is positive and is decreasing in . Since is bounded, it follows that is bounded, and thereforeas uniformly in N, and so the lower bound is completed.

Equipped with this result, we can now prove the legitimacy of the Coulomb gas approximation in our complexity calculation. The proof will require an elementary intermediate result which has undoubtedly appeared in various places before, but we prove it here anyway for the avoidance of doubt.

Lemma 5

Let be a random symmetric real matrix with independent centred Gaussian upper-diagonal and diagonal entries. Suppose that the variances of the entries are bounded above by for some constant . Then there exists some constant such that

Let denote the variance of . ThenSimple integration with a change of variables givesand then, for large enough N, Stirling’s formula givesSo finallyso defining gives the result.

Theorem 1

For any , let be a random matrix distributed as in the statement of Lemma 4. Then as

Let be some constant, independent of N. Introduce the notation , and thenWe have the upper bound (81) of Lemma 4 but this cannot be directly applied to (107) since the bound relies on uniformity in which can only be established for bounded . We use a much cruder bound instead. First, letand then has centred Gaussian entries with variance , so Lemma 5 applies, and we findfor some constant which is independent of and N, but we need not compute it.

Now we haveBut, since is bounded and has compact support, we can choose R large enough (independent of N) so thatfor all with and for some fixed L independent of N. Whenceas . Finally, for in , the result of the Lemma 4 holds uniformly in , soThe result follows from (113), (114) and the triangle inequality.

Asymptotic Complexity with Prescribed Hessian Index

Recall the complexity defined in (7):The extra Hessian signature conditions in (7) enforce that both generator and discriminator are at low-index saddle points. Our method for computing the complexity in the previous subsection relies on the Coulomb gas approximation applied to the spectrum of . However, the Hessian index constraints are formulated in the natural Hessian matrix (16), but our spectral calculations proceed from the rewritten form (21). We find however that we can indeed proceed much as in [13]. Recall the key Hessian matrix given in (16) bywhere , , G is Ginibre, and all are independent. Note that we have used (23) to slightly rewrite (16). We must address the problem of computingIndeed, we introduce integration variables , being -vectors of commuting and anti-commuting variables respectively. Use [t] notation to split all vectors into the first and last components. LetWith these definitions, we have [13]where is the normalised measure of the and the ellipsis represents terms with no dependence on or , which we need not write down. The crux of the matter is that we must computebut [13] has performed exactly these calculations (see around (5.146) therein) and so there exist constants such thatandwhereHere is the rate function of the largest eigenvalue of the GOE as obtained in [45] and used in [2, 13]:Note that for and for we simply have . Note also that

We have successfully dealt with the Hessian index indicators inside the expectation, however we need some way of returning to the form of in (21) so the complexity calculations using the Coulomb gas approach can proceed as before. We can achieve this with inverse Fourier transforms:from which we obtainIt follows thatSo we see that the relevant exponent in this case is the same as for but with additional GOE eigenvalue large deviation terms, giving the complexity limitPlots of for a few values of are shown in Fig. 5.

Fig. 5

Contour plots of for a few values of . Here

Remark 2

Recall that the limiting spectral measure of the Hessian displays a transition as the support splits from one component to two, as shown in Fig. 1. Let us comment on the relevance of this feature to the complexity. The spectral measure appears in one place in the above complexity calculations: the Coulomb gas integral . The effect of integrating against the measure is to smooth out the transition point. In other words, if has two components or is at the transition point, one expects to be able to construct another measure supported on a single component such that . We interpret this to mean that the Coulomb gas integral term does not display any features that can be unambiguously attributed to the transition behaviour of the spectral measure.

Implications

Structure of Low-Index Critical Points

We examine the fine structure of the low-index critical points for both spin glasses. [1] used the ‘banded structure’ of low-index critical points to explain the effectiveness of gradient descent in large multi-layer perceptron neural networks. We undertake to uncover the analogous structure in our dual spin-glass model and thence offer explanations for GAN training dynamics with gradient descent. For a range of values, starting at (0, 0), we compute on an appropriate domain. In the plane, we then find the maximum , and separately , such that . In the large N limit, this procedure reveals the regions in the plane where critical points of each index of the two spin glasses are found. Figure 6 plots these maximum values as contours on a shared plane. The grey region in the plot clearly shows the ‘ground state’ boundary beyond which no critical points exist. We use some fixed values of the various parameters: .

Fig. 6

Contours in the plane of the maximum and such that . results shown with a red colour red scheme, and with blue/green. The grey region on the left lies outside the domain of definition of . Here . The arrow indicates the approximate location of the contour discussed in the main text

These plots reveal, unsurprisingly perhaps, that something resembling the banded structure of [1] is present, with the higher index critical points being limited to higher loss values for each network. The 2-dimensional analogues of the boundary of [1] are evident in the bunching of the contours at higher values. There is, however further structure not present in the single spin-glass multi-layer perceptron model. Consider the contour of at the bottom of the full contour plot in Fig. 6. Imagine traversing a path near this contour from right to left (decreasing values); an example path is approximately indicated by a black arrow on the figure. At all points along such a path, the only critical points present are exact local minima for both networks, however the losses range over These three states correspond qualitatively to known GAN phenomena:

low generator loss, high discriminator loss;

some balance between generator and discriminator loss;

high generator loss, low discriminator loss.

discriminator collapses to predicting ‘real’ for all items;

successfully trained model;

generator collapses to producing garbage samples which the discriminator trivially identifies.

Overall, the analysis of our model reveals a loss surface that favours convergence to states of low loss for at least one of the networks, but not necessarily both. Moreover, our plots of and in Figs. 3, 5 demonstrate clearly the competition between the two networks, with the minimum attainable discriminator loss increasing as the generator loss decreases and vice-versa. We thus have a qualitative similarity between the minimax dynamics of real GANs and our model, but also a new two-dimensional banded critical points structure. We can further illuminate the structure by plotting, for each , the approximate proportion of minima with both and out of all points where at at least one of those conditions holds. The expression iswhich gives the log of the ratio in units of N. We show the plot in Fig. 7. Note that, for large N, any region of the plot away from a value of zero contains exponentially more bad minima – where one of the networks has collapsed – than good minima, with equilibrium between the networks. The model therefore predicts the existence of good local minima (in the bottom left of Fig. 7) that are effectively inaccessible due to their being exponentially outnumbered by bad local minima.

Fig. 7

Contour plot of the log ratio quantity given in (132). This is the approximate proportion of minima with both and out of all points where at at least one of those conditions holds

The structure revealed by our analysis offers the following explanation of large GAN training dynamics with gradient descent:

As with single feed-forward networks, the loss surface geometry encourages convergence to globally low values of at least one of the network losses.

The same favourable geometry encourages convergence to successful states, where both networks achieve reasonably low loss, but also encourages convergence to failure states, where the generator’s samples are too easily distinguished by the discriminator, or the discriminator has entirely failed thus providing no useful training signal to the generator.

Remark 3

A natural question in the context of our analysis of low-index critical points is: do such points reflect the points typically reached by gradient descent algorithms used to train real GANs? There has been much discussion in the literature of the analogous question for single networks and spin glasses [1, 7, 9]. It is not clear how to settle this question in our case, but we believe our model and its low-index critical points give a description of the baseline properties to be expected of high-dimensional adversarial optimisation problems late in the optimisation procedure. In addition, the unstructured random noise present in spin glasses may be more appropriate in our model for GANs than it is for single spin-glass models of single networks, as GAN generators do genuinely contain unstructured latent noise, rather than just the highly-structured data distributions seen on real data.

Remark 4

The issue of meta-stability is also worth mentioning. In single spin glasses, the boundary between fixed index and unbounded index critical points is meta-stable [46, 47]. From the random matrix theory perspective, the boundary corresponds to the left edge of the Wigner semi-circle [2]. There are O(N) eigenvalues in any finite interval at the left of the Wigner semi-circle, corresponding to O(N) Hessian eigenvalues in any neighbourhood around zero. The 2D analogue of the boundary in our double spin-glass model is expected to possess the same meta-stability: the Wigner semi-circle is replaced by the measure studied in Sect. 4, to which the preceding arguments apply. In the context of deep neural networks, there is a related discussion concerning “wide and flat local optima” of the loss surface, i.e. local optima for which many of the Hessian eigenvalues are close to zero. There are strong indications that deep neural networks converge under gradient-based optimisation to such optima [48-53] and that they are perhaps better for generalisation (i.e. test set loss) than other local optima, however some authors have challenged this view [54-58]. It is beyond the scope of the present work to analyse the role of meta-stability further, however we note that the indications from machine learning are that it is most significant when considering generalisation, however our work simplifies to the case of a single loss rather than separately considering training and test loss.

Hyperparameter Effects

Our proposed model for GANs includes a few fixed hyperparameters that we expect to control features of the model, namely and . Based on the results of [1, 2, 13], and the form of our analytical results above, we do not expect p and q (the number of layers in the discriminator and generator) to have any interesting effect beyond ; this is clearly a limitation of the model. We would expect there to exist an optimal value of that would result in minimum loss, in some sense. The effect of is less clear, though we guess that, in the studied limit, all are effectively equivalent. Intuitively, choosing corresponds to one network having a negligible number of parameters when compared with the other and we would expect the much larger network to prevail in the minimax game, however our theoretical results above are valid strictly for .

In the following two subsections we examine effect of and in our theoretical and in real experiments with a DCGAN [26]. Additional supporting plots are given in the appendix.

Effect of Variance Ratio

In the definition of complexity, and are upper bounds on the loss of the discriminator and generator, respectively. We are interested in the region of the plane such that , this being the region where gradient descent algorithms are expected to become trapped. We therefore investigate the minimum loss such that , this being, for a given , the theoretical minimum loss attainable by the GAN. We consider two natural notions of loss: We vary over a range of values in and compute .

;

.

To compare the theoretical predictions of the effect of to real GANs, we perform a simple set of experiments. We use a DCGAN architecture [26] with 5 layers in each network, using the reference PyTorch implementation from [59], however we introduce the generator noise scale . That is, the latent input noise vector for the generator is sampled from . For a given , we train the GANs for 10 epochs on CIFAR10 [60] and record the generator and discriminator losses. For each , we repeat the experiment 30 times and average the minimum attained generator and discriminator losses to account for random variations between runs with the same . We note that the sample variances of the loss were typically very high, despite the PyTorch random seed being fixed across all runs. We plot the sample means, smoothed with rolling averaging over a short window, in the interest of clearly visualising whatever trends are present. The results are shown in Fig. 8.

Fig. 8

The effect of . Comparison of theoretical predictions of minimum possible discriminator and generator losses to observed minimum losses when training DCGAN on CIFAR10. The blue cross-dashed lines show the experimental DCGAN results, and solid red lines show the theoretical results . and are used in the theoretical calculations, to best match the DCGAN architecture. is shown on a log-scale

There is a striking similarity between the generator plots, with a sharp decline between and around , after which the minimum loss is approximately constant. The picture for the discriminator is less clear. Focusing on the sections , both plots show a clear minimum, at around in experiments and in theory. Note that the scales on the y-axes of these plots should not be considered meaningful. Though there is not precise correspondence between the discriminator curves, we claim that both theory and experiment tell the same qualitative story: increasing to at least around gives the lowest theoretical generator loss, and then further increasing to, tentatively, some value in gives the lowest possible discriminator loss at no detriment to the generator.

We are not aware of tuning being widely used in practice for real GANs, rather it is typically taken to be unity. We have chosen this parameter, as it can be directly paralleled in our spin glass model, therefore allowing for the above experimental comparison. Naturally there are other parameters of real GANs that one might wish to study (such as learning rates and batch sizes) however these are much less readily mirrored in the spin glass model and complexity analysis, precluding comparisons between theory and experiment. Nevertheless, the experimental results in Fig. 8 do demonstrate that tuning in real GANs could be of benefit, as does not appear to be the optimal value.

Effect of Size Ratio

Similarly to the previous section, we can investigate the effect of using while varying over (0, 1). To achieve this variation in the DCGAN, we vary the number of convolutional filters in each network. The generator and discriminator are essentially mirror images of each other and the number of filters in each intermediate layer are defined as increasing functions3 of some positive integers . We fix and vary to obtain a range of values, with . The results are shown in Fig. 9.

Fig. 9

The effect of . Comparison of theoretical predictions of minimum possible discriminator and generator losses to observed minimum losses when training DCGAN on CIFAR10. The blue cross-dashed lines show the experimental DCGAN results, and the solid red show the theoretical results . and are used in the theoretical calculations, to best match the DCGAN architecture

The theoretical model predicts a a broad range of equivalently optimal values centred on from the perspective of the discriminator loss, and no effect of on the generator loss. The experimental results similarly show a broad range of equivalently optimal centred around , however there appear to be deficiencies in our model, particularly for higher values. The results of the experiments are intuitively sensible: the generator loss deteriorates for closer to 1, i.e. when the discriminator has very many more parameters than the generator, and vice-versa for small .

Conclusions and Outlook

We have contributed a novel model for the study of large neural network gradient descent dynamics with statistical physics techniques, namely an interacting spin-glass model for generative adversarial neural networks. We believe this is the first attempt in the literature to incorporate advanced architectural features of modern neural networks, beyond basic single network multi-layer perceptrons, into such statistical physics style models. We have conducted an asymptotic complexity analysis via Kac-Rice formulae and Random Matrix Theory calculations of the energy surface of this model, acting as a proxy for GAN training loss surfaces of large networks. Our analysis has revealed a banded critical point structure as seen previously for simpler models, explaining the surprising success of gradient descent in such complicated loss surfaces, but with added structural features that offer explanations for the greater difficulty of training GANs compared to single networks. We have used our model to study the effect of some elementary GAN hyper-parameters and compared with experiments training real GANs on a standard computer vision dataset. We believe that the interesting features of our model, and their correspondence with real GANs, are yet further compelling evidence for the role of statistical physics effects in deep learning and the value of studying such models as proxies for real deep learning models, and in particular the value of concocting more sophisticated models that reflect aspects of modern neural network design and practice.

Our analysis has focused on the annealed complexity of our spin glass model (i.e. taking the logarithm after the expectation) rather than the quenched complexity (i.e. taking the expectation after the logarithm). Ideally one would compute both, as the quenched complexity is often considered to reflect the typical number of stationary points and is bounded above by the annealed complexity. Computing the quenched complexity is typically more challenging than the annealed and such a calculation for our model could be the subject of a further work requiring considerable technical innovations. Even the elegant and very general methods presented recently in [40] are restricted only to the annealed case. Agreement between annealed and quenched is known only in a few special cases closely related to spherical spin glasses [61-63] and is not expected in general [33]. It is conceivable that quenched and annealed complexity agree in the case of our model, as it closely related to spin glasses and possesses no distinguished directions (i.e. spikes) such as are present in [33]. Establishing agreement by existing methods requires analysis of pairs of correlated GOE-like matrices. Such an approach for our model may well require analysis of at least 4 correlated matrices (2 per diagonal block), and quite possibly more, including correlations between blocks. We leave this considerable challenge for future work.

From a mathematical perspective, we have extensively studied the limiting spectral density of a novel random matrix ensemble using supersymmetric methods. In the preparation of this paper, we made considerable efforts to complete the average absolute value determinant calculations directly using a supersymmetric representation, as seen in [13], however this was found to be analytically intractable (as expected), but also extremely troublesome numerically (essentially due to analytically intractable and highly complicated Riemann sheet structure in ). We were able to sidestep these issues by instead using a Coulomb gas approximation, whose validity we have rigorously proved using a novel combination of concentration arguments and supersymmetric asymptotic expansions. We have verified with numerical simulations our derived mean spectral density for the relevant Random Matrix Theory ensemble and also the accuracy of the Coulomb gas approximation.

We hope that future work will be inspired to further study models of neural networks such as we have considered here. Practically, it would be exciting to explore the possibility of using our insights into GAN loss surfaces to devise algorithmic methods of avoiding training failure. Mathematically, the local spectral statistics of our random matrix ensemble may be interesting to study, particularly around the cusp where the two disjoint components of the limiting spectral density merge.