Sampling the space of solutions of an artificial neural network

First, we trained the model with GD for approximately $\sim 10^{7}$ epochs and with a fixed learning rate $\eta=1.0$. The solutions $\boldsymbol{s}_{i}$ found with GD, from which we started the sampling, have mean energy $\langle U\rangle_{\boldsymbol{s}}=(8.5\pm 0.04)\cdot 10^{-2}$ and mean training error $\langle\epsilon\rangle_{\boldsymbol{s}}=0$ (Fig. 2a). The similarity between them is peaked at $q\approx 0.6$ (see upper panel in Fig. 2c), which means that GD finds solutions with an almost fixed mutual distance.

Starting from GD solutions, we explored the space of parameters of the network with the HMC algorithm at different temperatures (Fig. 2b). For $T<1.8\cdot 10^{-1}$ the average energy $\langle U\rangle_{T}$ starts to flatten, suggesting that here the system freezes in the lowest–energy available states. The average energy is comparable with that of the solutions $\boldsymbol{s}_{i}$ and also the mean training error remains $\langle\epsilon\rangle=0$. At these temperatures the sampling of the space of weights is difficult and the system gets trapped in local minima of the loss. The similarity distribution $\rho(q)$ (see Eq. (19)) calculated on weights sampled from HMC trajectories starting from the same initial condition (intra-state overlap distribution) is strongly peaked around 1, see Fig. 2c, middle panel. This is markedly different from the inter-state overlap distribution, which is obtained by measuring the overlap between weights along HMC trajectories starting from strictly different initial conditions and is peaked at $q\approx 0.63$. Thus, the system is not able to equilibrate at these temperatures but stays close to the initial condition.

In the range of temperature $1.8\cdot 10^{-1}<T<1.8\cdot 10^{2}$ the average energy increases almost linearly $\langle U\rangle_{T}\sim T^{\delta}$, with a scaling exponent $\delta\approx 0.989\pm 0.019$. This is typical of a glassy thermodynamic phase with a sub-exponential number of accessible states. In this intermediate temperature range, there is no marked difference between the inter and intra state overlap distribution (Fig. 2c, lower panel). This suggests that the system is not trapped in local minima as in the low–temperature phase. The average similarity between states displays a single peak centered at $q\approx 0.6$, similar to that observed at low temperature, but with a larger variance.

In the high temperature phase, the average energy deviates from the linear behavior and the training error further increases. The distribution of similarities $q$ is very broaden. However, here the HMC algorithm becomes inefficient because the high velocities extracted before each HMC move require a very small step $\delta t$ for a correct integration of the equations of motion and thus for an acceptable acceptance rate.

We summarize these results in Fig. 2d, where we show a sketch of the space of parameters at fixed norm. On this hypersphere, the solutions found by the GD define energy basins which are explored by the system at low and intermediate temperatures. These basins are disconnected at low temperature, since HMC remains confined near the initialization and cannot find pathways connecting the different basins while remaining at low energy and at zero training error.

To explore the connectivity of energy basins at low temperatures, we performed double-ratchet simulations, hoping to find paths which are not found by the HMC algorithm. This method is designed to identify low-barrier pathways between two weight configurations by following the minimal gradient in the direction that brings them closer together (Fig. 3a).

The double ratchets are initialized on different pairs of solutions $(\boldsymbol{s}_{i},\boldsymbol{s}_{j})$ found after a short thermalization of the system of approximately $10^{4}$ HMC steps. Along each simulation, the weights do not cross significant energy barriers, as compared to the average energy at this temperature, keeping the training error equal to zero. Conversely, the linear interpolation between pairs of points along the double–ratchet trajectory results in barriers that are substantially higher than the average energy (Fig. 3b). These results suggest that, although the solutions found by gradient descent are linearly mode disconnected, low-energy tortuous paths joining them still exist.

The anchor point $\boldsymbol{s}^{\star}$ obtained at the end of the double–ratchet trajectory has an energy that is slightly lower than the one of the initial points $(\boldsymbol{s}_{i},\boldsymbol{s}_{j})$ and $\langle U\rangle_{\boldsymbol{s}^{\star}}=(6.80\pm 0.09)\cdot 10^{-2}$ is compatible with the average $\langle U\rangle_{T_{L}}$, that is the constrained dynamic arising from the double–ratchet simulation helps the equilibration of the system. The anchor weights $\boldsymbol{s}^{\star}$ obtained from several double ratchet simulations initialized from different weight pairs $(\boldsymbol{s}_{i},\boldsymbol{s}_{j})$ exhibit larger similarity than the respective initial vectors and are separated by energy barriers along their linear interpolation (Fig. 3c). This suggests the presence of a compact and narrow region of the solution space that must be traversed to move from one gradient descent solution to another. The similarity between the anchor weights $\boldsymbol{s}^{\star}$ and the associated initial weights is also centered around $0.6$, as that between the starting solutions.

In summary (Fig. 3d), the basins identified by GD solutions are connected by non–linear paths pointing towards a denser, more compact region with slightly lower energy. The low–energy manifold of the space of weights has thus a spiky shape, with the GD solutions lying on its spikes.

We then studied the geometry of the center of the spiky low energy manifold in more detail using coupled replica simulations, starting with five GD solutions coupled to their barycenter and increasing the harmonic couplings until all the running points converge to a single vector $\boldsymbol{c}$ (Fig. 4a).

We compared the properties of the centers found in different coupled replica simulations with the solutions $\boldsymbol{s}$ found by GD and the anchor points $\boldsymbol{s}^{\star}$ found by double ratchets initialized on different GD solutions at $T=T_{L}$. The final configuration $\boldsymbol{c}$ always has an energy that is lower by several standard deviations than the mean energy at the simulation temperature (Fig. 4b,c). The centers obtained from the coupled replica simulations are confined to a very narrow region of the solution space (Fig. 4b), even narrower than the region spanned by the anchor points (middle panel of Fig. 3c).

In addition, the similarity distribution between centers and anchor points is larger than that between centers and GD solutions. This ordering suggests that the centers lie deep within the bulk of the solution manifold, with anchor points positioned more peripherally and GD solutions even farther away. This “nested overlap” structure has already been observed in the negative perceptron problem, see Annesi et al. (2024).

Despite the high similarity between the $\boldsymbol{c}$ vectors ($q\approx 0.90\pm 0.03$), the latter are not linearly connected at fixed norm (Fig. 4b, upper panel, green curves), even though the barrier height is significantly lower than the one between $\boldsymbol{c}$ and GD solutions $\boldsymbol{s}$, as well as $\boldsymbol{s}^{*}$ points (Fig. 4b, upper panel, blue and red curves, respectively).

Finally, we performed double–ratchet simulations between pairs of $\boldsymbol{c}$ vectors and between centers and GD solutions. In the latter case, at $T=T_{L}$ the weight initialized on the center $\boldsymbol{c}$ stays very close to it along the double ratchet dynamics (and its energy within two standard deviations from the average at the simulation temperature); the weight starting from the GD solution instead slowly lowers its energy until the coupled $\boldsymbol{c}$ vector is reached (Fig. 4c, blue curves).

On the contrary, when the same simulation is performed between pairs of center points at $T<T_{L}$, the latter are able to connect to each other through paths whose energy is lower than the average $\langle U\rangle_{T_{L}}$ (Fig. 4c, green curves).

From this last computational analysis, we conclude that the center of the spiky low–energy manifold is not barrier–free. Although there are low–energy paths connecting all low–energy regions of the structure, these are simple as linear paths. Moreover, there is a small slope of the energy towards the center.

By increasing the temperature to intermediate values (cf. Fig. 2), the system reaches a phase where the HMC algorithm can easily sample the space of weights. In this phase the training error is still small but non–negligible (up to $\epsilon\approx 0.3$) but the system can still learn some of the data we present.

The similarity distribution between states shows a single peak centered at $q\approx 0.6$, indicating that most of the sample states are equivalent to each other. This is different from what is found at low temperatures, where there is a difference between points at the periphery of the structure and weights in the center of the low–energy manifold. Such a single peak in $\rho(q)$ is compatible with what is known in the language of complex systems as replica-symmetric behavior of the system Mézard et al. (1987).

The states sampled during simulations in this intermediate regime ($T=1.8\cdot 10^{1}$) show an average similarity $\langle q\rangle=0.54\pm 0.03$ with respect to GD solutions $\boldsymbol{s}_{i}$ and $\langle q\rangle=0.58\pm 0.04$ with respect to both the $s^{*}$ and $\boldsymbol{c}$ vectors found in the low–temperature phase. Consequently, here the system is never close to any specific region of the subspace explored at low temperatures. In particular, the large mean distance between the center of the manifold and the typical states sampled at intermediate temperatures suggests that the former is entropically unfavorable. This fact, along with the unimodal shape of the cosine–similarity distribution, suggests that the visited manifold at intermediate temperatures still retains a (less rugged) spiky shape, but with a “hollow” center characterized by high free energy.

We then compared the properties of the space of weights that we found at the interpolation threshold with that of the overparametrized regime.

We performed various HMC simulations at $\alpha=0.2$ for a wide range of temperatures, all of which never indicate a frozen behavior as, similarly to the intermediate–temperature case at $\alpha=1.8$, the intra e inter overlap distribution coincide, see Fig. 5a. Moreover, the energy profile along the linear paths between the sampled states shows a central free energy barrier, both at $\alpha=1.8$ and at most temperatures at $\alpha=0.2$. In the last case, the central barrier disappears for vanishing temperatures and the energy profile assumes a convex shape (Fig. 5b).

In conclusion, the intermediate–temperature configurations visited by the system are similarly distributed both near the interpolation threshold and in the overparametrized regime, where in both scenarios the manifold exhibits a symmetric shape, in the sense that the system populates a single kind of state, belonging to the spiky periphery of the space. Only at low temperatures do the two cases differ, since in the former the explored subspace maintains a spiky shape, where part of the center of the manifold and its periphery are populated, whilst in the latter it becomes convex.

In the overparametrized regime, we can compute analytically both the energy (loss) and the training error along the geodesic path connecting two weights $\boldsymbol{w}^{1}$, $\boldsymbol{w}^{2}$ in the thermodynamic limit, extending the previous analysis performed in ref. Annesi et al. (2023) to the tree committee machine case. In full generality we will assume to sample the two weights from two different Boltzmann distributions $\boldsymbol{w}^{1}\sim p_{1}(\boldsymbol{w};\mathcal{D})$, $\boldsymbol{w}^{2}\sim p_{2}(\boldsymbol{w};\mathcal{D})$ whose individual form are the same as in equation (7) but each of them differing from the choice of the loss function and the prior, i.e.

Notice, moreover, that in equations (22) the training data $\mathcal{D}$ is the same for both Boltzmann distributions. As in Section II we choose a Gaussian distribution as a prior $p_{1}(\boldsymbol{w})$, $p_{2}(\boldsymbol{w})$, with a $L2$ regularization parameter that we denote respectively by $\lambda_{1}$ and $\lambda_{2}$. In the large $N$ limit the choice of the regularization parameter will induce a non-trivial value of the norm of the weights $\boldsymbol{w}^{1}$ and $\boldsymbol{w}^{2}$. In the following we will suppose that $\lambda_{1}$ and $\lambda_{2}$ are chosen such that the norm of $\boldsymbol{w}_{1}$ is the same as the one of $\boldsymbol{w}^{2}$, i.e. $|\boldsymbol{w}^{1}|^{2}=|\boldsymbol{w}^{2}|^{2}\equiv Q$.

We are interested in computing the average training error and training loss landscape on the geodesic path joining $\boldsymbol{w}^{1}$ and $\boldsymbol{w}^{2}$ at fixed squared norm $Q$. This can be obtained as follows. Given $\boldsymbol{w}^{1}$ and $\boldsymbol{w}^{2}$ we define the interpolating weight $\boldsymbol{w}(\gamma)$ with parameter $\gamma\in[0,1]$ as the vector whose components are the linear interpolation of the components of $\boldsymbol{w}^{1}$ and $\boldsymbol{w}^{2}$

In order to compute the geodesic path at the same squared norm $Q$, we finally need to project the whole straight path on the hypersphere of radius $\sqrt{Q}$. This defines the weight

where

In the previous expression we have introduced the quantity

which corresponds to the overlap between $\boldsymbol{w}^{1}$ and $\boldsymbol{w}^{2}$. In the previous equation we have denoted by $\langle\cdot\rangle$ the average over the Boltzmann distribution (7) and by $\mathbb{E}_{\mathcal{D}}$ the average over the dataset.

We are interested in finding what is the average value of a loss $\widetilde{\mathcal{L}}$ of the projected interpolated weight $\widetilde{\boldsymbol{w}}(\gamma)$ as a function of $\gamma$, and how this profile depends on the choice of the endpoints $\boldsymbol{w}^{1}$ and $\boldsymbol{w}^{2}$ via the probability density functions in (22). In formulas, what we compute is

Using a general loss $\widetilde{\mathcal{L}}$ allows us also to both have access to training loss and training error profiles. Equation (27) can be computed by the replica method in the large $N$ limit. In the following we will focus only on the infinite width limit $K\to\infty$ with the ratio $K/N\to 0$, as done in Baldassi et al. (2019). The full calculation is reported in Appendix B. Here we only mention that the result of the calculation, assuming no replica symmetry breaking of the solution space (i.e. in the so called Replica Symmetric ansatz), will depend on simple geometrical quantities. Those are the typical overlap between two weights $\boldsymbol{w}^{a}$, $\boldsymbol{w}^{b}$ extracted from $p_{1}$, $\boldsymbol{w}^{a},\boldsymbol{w}^{b}\sim p_{1}(\boldsymbol{w};\mathcal{D})$ (respectively from $p_{2}$, i.e. $\boldsymbol{w}^{a},\boldsymbol{w}^{b}\sim p_{2}(\boldsymbol{w};\mathcal{D})$), i.e.