Investigation of Finite-Size 2D Ising Model with a Noisy Matrix of Spin-Spin Interactions

We analyze changes in the thermodynamic properties of a spin system when it passes from the classical two-dimensional Ising model to the spin glass model, where spin-spin interactions are random in their values and signs. Formally, the transition reduces to a gradual change in the amplitude of the multiplicative noise (distributed uniformly with a mean equal to one) superimposed over the initial Ising matrix of interacting spins. Considering the noise, we obtain analytical expressions that are valid for lattices of finite sizes. We compare our results with the results of computer simulations performed for square N = L × L lattices with linear dimensions L = 50 ÷ 1000. We find experimentally the dependencies of the critical values (the critical temperature, the internal energy, entropy and the specific heat) as well as the dependencies of the energy of the ground state and its magnetization on the amplitude of the noise. We show that when the variance of the noise reaches one, there is a jump of the ground state from the fully correlated state to an uncorrelated state and its magnetization jumps from 1 to 0. In the same time, a phase transition that is present at a lower level of the noise disappears.


Introduction
Calculation of the partition function is an essential of statistical physics and informatics. A few conceptual models allow exact solutions [1][2][3][4][5][6]. Among these a 2D Ising model [7], though simple, deserves special attention because of its importance for investigating critical effects. Having contributed a lot to the development of the spin glass theory, the Edwards-Anderson model [8] and Sherrington-Kirkpatrick model [9] are also worth mentioning. However, there are not many models that permit exact solutions. This is the reason why numerical methods are mostly used for tackling complex systems. Of them, two methods are most suitable for our purpose. The first is the Monte-Carlo method [10,11]. It enables us to analyze a system and determine its critical parameters quite accurately [12][13][14][15][16]. The thorough consideration of the method can be found in papers [17,18]. Unfortunately, the method needs a great deal of computations and does not allow direct calculation of the free energy. The second method uses the approach [19,20], which has recently given rise to the fast algorithm [21,22] that finds the free energy by computing the determinant of a matrix. The algorithm is popular because it allows the user to compute the free energy quite accurately and at the same time determine the energy and configuration of the ground state of a system.
The methods of statistical physics help researchers to understand the behavior of complex neural nets and evaluate the capacity of neural-net storage systems [23][24][25][26][27][28]. The machine learning and computer-aided image processing need fast calculations of the partition function of specific interconnect matrices [29,30]. The realization of Hinton's ideas [31,32] gave rise to algorithms of deep learning and image processing [33][34][35][36]. Based on the optimization of the free energy of a spin (neuron) system, the algorithms, from the formal viewpoint, comes down to the optimization of the spin correlation in neighboring layers or within a single layer of a neural net. It should be understood that the system has a phase transition because the spin correlation grows abruptly at the critical point (the correlation length becomes nearly as great as the size of the whole system). In this case the optimization of the neural network becomes temperature dependent, which makes the learning algorithm almost impracticable.
The aim of the paper is to study the properties of a finite spin system whose Hamiltonian is defined as a quadratic Functional (1). The functional is often used in machine learning and image processing. Quantities s i = ±1 may stand for either pixel class (object/background) in an image [35], or the neuron activity indication in a Bayes neural network [36]. We will use the physical notation calling quantities s i = ±1 spins. The model under consideration has two limiting cases. The conventional 2D Ising model with regular interconnections presents the first case; the Edwards-Anderson model is the second case. The properties of our model lie somewhere in between. We introduce adjusting parameters in Functional (1), which allows us to go from the 2D Ising model to Edwards-Anderson model in a smooth manner and investigate the thermodynamic characteristics of the system in the transient state.
To avoid misunderstanding, let us point out two things. First, our interest is finite systems. For this reason, there is an expected discrepancy with Onsager results obtained at N → ∞ . Second, we cannot use the results of the spin glass theory to the full because the finite system under consideration is ergodic: it does not have multiple phase transitions caused by frustrations and provide self-averaging [37,38].

Essential Expressions, the Equation of State
Let us consider a system described by the Hamiltonian: This system consists of N Ising spins s i = ±1(i = 1, 2, . . . , N), positioned at the nodes of a planar grid, the nodes being numbered by index i. Only interactions with four nearest neighbors are considered. Spin-to-spin interactions J ij are random and defined as where ε ij is a random zero-mean variable uniformly distributed over the interval ε ij ∈ [−η, η]. We have chosen the uniform distribution of ε ij to be able to control J ij : when η ≤ 1, all interactions are positive J ij ≥ 0 . For the sake of simplicity, we assume that J = 1. Our interest is the free energy of the system: where the partition function Z = ∑ S e −NβE(S) is defined as a sum over all possible configurations S and β = 1/kT is the reverse temperature. The knowledge of the free energy makes it possible to compute the basic measurable parameters of the system: where free energy U = E is the ensemble average at given β, σ 2 = E 2 − E 2 is the variance of energy and C = β 2 σ 2 is the specific heat.
Along with that, we are interested in the configuration S 0 of the ground state, its energy E 0 = E(S 0 ) and the magnetization M 0 = 1 The properties of the system depend on the dimension of the system N and adjusting parameter η. Unfortunately, we cannot allow for the effect of the both parameters simultaneously, so we consider the contribution of each separately.

The Effect of the Finite Grid Dimension
Let us consider how the fact of the grid having a finite dimension affects its properties. Let us take η = 0 as the starting point. In this case the behavior of the system can be described by the expression (see reference [39]) which is true for finite systems with free boundary conditions: where Here K 1 = K 1 (κ) and K 2 = K 2 (κ) are full elliptical integrals of the first and second type correspondingly: Expressions (5)-(7) are the well-known Onsager solution [7], which is true for N → ∞ , modified for the case of finite N. Though true for N 1, the expressions agree well with the experimental data even at relatively small grid dimensions (L ∼ 25). As could be expected, when N → ∞ , Formula (6) give p → 1 , a 1,2 → 1 , ∆ → 0 , δ → 0 and Expression (5) turn into well-known ones [7]. Expression (5) agree excellently with experimental data: the relative error is less than 0.2% at L = 50. With the growing L, the error decreases rapidly and is within the limits of experimental error at L = 1000 10 −5 for σ 2 . By way of comparison Figures 3, 6 and 7 gives the plots of Function (5) for L = 400.
Expression (5) allow the N-dependences of the critical values of the reverse temperature, internal energy and energy variance of the system: where β ∞ = 1 2 ln

The Effect of Noise
Let us consider the random character of quantities J ij (η = 0). Let D(E) be the number of states of energy E. Then the sum of states can be presented as Z = ∑ E D(E)e −NβE . Passing from summation to integration, we get (to within an insignificant constant): where Ψ(E) = ln D(E)/N. Applying the saddle-point method to integral (9), we get The first expression in (10) defines the free energy, the second determines E at the saddle point where the derivative of function Ψ(E) − βE turns to zero.
The form of spectral function Ψ(E) is known only for the one-dimensional Ising model. That is why we turn to the so-called n-vicinity method [28] to calculate the spectral function. The idea of the method is to divide the whole space of 2 N states into N classes (n vicinities) and approximate the energy distribution in each class by a corresponding Gaussian. In brief, the approach is as follows: Let us denote the ground-state configuration as S 0 . Let class Ω n be a set of configurations S n that differs from S 0 in that they have n spins directed oppositely to the spins in S 0 . The number of configurations in the class is equal to the number of compositions of N in n, all configurations having the same (relative) magnetization m = N −1 · S m S T 0 = 1 − 2n/N. The distribution of state energies within the n-vicinity was shown [28] to follow the normal distribution D n (E): where Here E 0 is the ground state energy, σ 2 h0 is the variance of ground-state local fields. In this case we have σ 2 The sought-for distribution D(E) is found by summing D n (E) over all n. Using the Stirling formula and passing from summation to integration with respect to variable m = 1 − 2n/N, we get for D(E): where If we evaluate integral (13) by the saddle-point method, for the spectral function we get Ψ(E) = −F(m, E), where m is the solution of equation ∂F(m, E)/∂m = 0. Let us combine (13)- (14) and (9)- (10). Then the free energy can be written as where variables m = m(β) and E = E(β) are derived from the equations: It is easy to notice that set of Equation (16) is solvable when m = 0. Correspondingly, when the values of β are less than certain critical value β c , (16) and (12) gives us E m = 0, σ 2 m = 2 and E = −2β, the free energy taking the form f (β) = − ln N − β 2 . The phase transition occurs when β allows yet another solution to (16) at m = 0. Note that substituting the second equation from (16) into the first one allows us to eliminate variable E. Doing things this way and performing several transformations, we obtain the equation of state that holds only one variable m: where β = β/r. Here we introduced adjusting coefficient r to allow for the finite grid dimension: r = 1 when L → ∞ , r = 1.11 giving the excellent agreement with experiments at L ∼ 400. The critical temperature is defined as value β = β c at which there is a nontrivial solution of (17). This solution has to be found by a numerical method: when β > β c , we find m = 0 that satisfies (17) and compute the corresponding value of energy E = E m − βσ 2 m . Substitution of these values in (15) yields the corresponding value f (β).
Unfortunately, the n-vicinity method has an essential fault: it is applicable only when the condition ∑ J ij 2 / N∑ J 2 ij ≥ 4 ln 2 holds. In our case this condition works when (1 + σ 2 η ) · ln 2 ≤ 1, that is when η < 1.2. For such relatively small values of η Formulae (15)-(17) gives β c and f (β) that predict the experimental results well (see Figures 1 and 2).

Evaluating the Spectral Density
The algorithm we use allows us to compute function     f f and its derivatives. In turn, this allows us to investigate how energy distribution amplitude. Indeed, it is easy to derive from Formulae (10)

Evaluating the Spectral Density
The algorithm we use allows us to compute function     f f and its derivatives. In turn, this allows us to investigate how energy distribution varies with the noise amplitude. Indeed, it is easy to derive from Formulae (10) the equation for the spectral function:

Evaluating the Spectral Density
The algorithm we use allows us to compute function f = f (β) and its derivatives. In turn, this allows us to investigate how energy distribution D(E) = exp[NΨ(E)] varies with the noise amplitude. Indeed, it is easy to derive from Formulae (10) the equation for the spectral function: and its derivatives dΨ Note that Ψ(E) is entropy up to a constant and Equation (18) are well-known Legendre transformations, which are applicable for analyzing the spectral density of finite-dimension models [40,41]. It follows from these equations that when β varies from β = 0 to β = ∞, E changes from 0 to E 0 and for each value of β we have a pair of values of E and Ψ(E). In so doing we determine the form of function Ψ(E) and its derivatives. The plots of function Ψ(E) and its derivatives presenting experimental data are given in Section 4.
The minimum of function d 2 Ψ/dE 2 at point E = 0 changes into the maximum as the noise amplitude grows. Let us find η at which it occurs. It can be noticed that with E → 0 the entropy can be presented as the series: where µ 4 = E 4 /σ 4 J is the fourth cumulant, which in our case is described by the expression [28]: From (20)- (21) it follows that in the center point of the curve (E = 0) quantity d 2 Ψ/dE 2 is determined by expression: and the fourth derivative d 4 Ψ/dE 4 E=0 = µ 4 /σ 4 J changes its sign at η = η c , when µ 4 = 0:

The Experiment Description
We make an intensive use of the Kasteleyn-Fisher algorithm [19,20] here to compute the free energy of the 2D square spin system. The algorithm gives exact results because the finding of the partition function is reduced to computation of the determinant of a matrix generated in accordance with the model under consideration. The algorithm permits us to exactly calculate the free energy of a spin system for an arbitrary planar graph with arbitrary links in a polynomial time. More information about the algorithm can be found in [21]. In this paper, we use the realization [22] of the algorithm that can give the same results in a shorter time. Using this algorithm, we were able to examine the behavior of free energy f = f (β; η) and its derivatives for a few lattices of different dimensions N = L × L. Additionally, paper [22] offers the algorithm for searching the ground state. This algorithm helped us to investigate the energy and magnetization of the ground state as functions of noise amplitude. For each value, we generated a great number of matrices but the results were practically the same when we changed one matrix to another.
Let us point out that the both algorithms we use are only applicable to planar lattices. It means that we considered only lattices with free boundary conditions because lattices with periodic boundary conditions do not belong to a planar graph. The length of the lattice varied from L = 25 to L = 10 3 . Most of the plots present the results for L = 400. The results for other sizes did not differ qualitatively.
The free energy is computed to 15-digit accuracy after the decimal point. Because we use the finite-difference method to compute the derivatives, the number of significant digits after the decimal point is about 7 for U(β) and 4 for σ 2 (β). With large grid dimensions (L ∼ 1000) and with β > 1 the computation error becomes too big and the plots of second derivatives start exhibiting oscillations. It is interesting to notice that introduction of little noise into grid interconnections allows us to decrease these oscillations.

Experimental Results
In the experiments, we calculate the free energy and its derivatives and find the ground-state configuration and energy. The accent is given to the finding of the critical point and corresponding quantities. The location of the maximum of curve σ 2 = σ 2 (β) is used to find the critical temperature. Most important experimental data are presented in Figures 1-7

Experimental Results
In the experiments, we calculate the free energy and its derivatives and find the ground-state configuration and energy. The accent is given to the finding of the critical point and corresponding quantities. The location of the maximum of curve     2 2 ( ) is used to find the critical temperature.
Most important experimental data are presented in Figures 1-7 and Table 1.

Experimental Results
In the experiments, we calculate the free energy and its derivatives and find the ground-state configuration and energy. The accent is given to the finding of the critical point and corresponding quantities. The location of the maximum of curve     2 2 ( ) is used to find the critical temperature.
Most important experimental data are presented in Figures 1-7 and Table 1.

The Free and Internal Energy
Experimental dependencies f = f (β) and U = U(β) are shown in Figures 1 and 2. It is seen from Figure 1 that the curves go down with η because the ground-state energy grows. When noise is small (η < 1.2), the curves of free energy f (β) and internal energy U(β) almost merge (Figures 1 and 2). When η < 1.7 the curves U(β) demonstrate a cusp ( Figure 2) which corresponds to the phase transition. When η ∼ 1.7, the cusp disappears and the further increase of noise changes only the asymptotic behavior of curves f (β) and U(β) according to (26).

The Energy Variance
Curves σ 2 = σ 2 (β) are shown in Figure 3. It should be noted that because the n-vicinity method gives a piecewise-linear approximation of the energy variance, the red marks in Figure 3 indicates values obtained by using the generalization of Onsager solution to a finite-dimension case according to Formula (5). The formula gives a perfect agreement with experimental data, yet it is applicable only in a zero-noise case.
The behavior of curves σ 2 = σ 2 (β) near point β = 0 is quite expected for any η: when β = 0, the energy variance is equal to σ 2 0 and, according to (20), grows gradually in proportion to noise variance σ 2 η = η 2 /3. With great β the behavior of curves σ = σ(β) is much dependent on η. It is seen in Figure 3 that the energy variance peaks corresponding to the phase transition are observed only at η < 1.7. The peaks become lower with growing η and move to the right at the same time. When η > 1.8, the peaks disappear at all, only the maximum at β = 0 remains.
It is interesting that all the curves in Figure 3a have the common intersection point near β ≈ 0.29. We could not find out why it is so. The intersection moves to the right slowly with the growing noise amplitude.

The Critical Temperature
The critical temperature is defined by the location of the maximum of curve σ = σ(β) or by the presence of a cusp on it. Figure 4 shows how the variance peak location and height vary with the growing noise. Holding true only for η < 1.2, the numerical solution of the equation of state (17) gives β s that agrees with the experimental data perfectly. For greater η it is possible to use the approximate expression resulting from the experiment: where β c0 is the zero-noise critical value resulted from (8). The peak height lowers linearly with the growing noise amplitude: where σ 2 s0 is the variance at η = 0 defined in (8). It follows that if η ≈ √ 3, σ 2 s falls to zero. It means that when η > √ 3, the variance peak disappears and we can say that the critical temperature is zero.

The Ground State
The results we obtained testify that when the noise amplitude η ≈ 1.7 (at σ η ≈ 1), the quality of the system changes. The ground state configuration experiences the most noticeable changes (see Figure 5). Clear that with zero noise the ground state is fully correlated, that is, all spins are the same s i = 1. The situation keeps as long as all matrix elements J ij > 0, that is, η < 1. However, (see Figure 5) the ground-state energy proved to remain almost the same for σ η as big as σ η ≈ 1. Then it starts decreasing gradually and comes to an asymptotic value [42]: corresponding to the energy of the ground state in the Edwards-Anderson model. The ground-state magnetization changes stepwise from 1 to 0 when the noise deviation comes close to unit σ η ≈ 1.

The Entropy
The change of the ground-state configuration and energy results in a change of energy distribution density Ψ(E). The curves of Ψ(E) and its derivatives are shown in Figures 6 and 7.
The disappearance of the phase transition is easy to notice if we look at the curve of the second derivative d 2 Ψ/dE 2 . It is seen in Figure 7a that the sink in the middle of the curve (E = 0) rises with growing η and, according to (23) the minimum of d 2 Ψ/dE 2 at E = 0 turns into a maximum when η ≈ 1.5. The peaks at points E = ±U c separate with growing η ( U c → E 0 ) and become lower like d 2 Ψ/dE 2 = −σ −2 c until full disappearance at η ≈ 1.7. When η > 1.7, curve d 2 Ψ/dE 2 has a noticeably convex shape and the phase transition peaks disappear. Moreover, in this case function d 2 Ψ/dE 2 is well described by the expression: Formula (27) gives good approximation of experimental data (accurate to 0.5% over the energy interval 0 ≤ |E| ≤ 0.91|E 0 |).

Conclusions
In this paper, we have considered the Ising model on a two-dimensional grid with noise-polluted interconnections. In the limiting case N → ∞ such system demonstrates the following properties: with low noise the system have all characteristics of conventional Ising model, with high noise it turns into the Edwards-Anderson spin glass model. The goal of our experiments was to observe the transition between these two limiting cases in the finite-dimension system N ≤ 10 6 . It proved that when the noise is weak (σ η < 1), the behavior of the system is much like the behavior of the conventional Ising model. We expected that with heavy noise (σ η >> 1), the behavior of the system would be like that of the Edwards-Anderson model. However, the experimental results are significantly different from the expectation. It turned out that even when the noise is relatively weak ( σ η ∼ 1), the system undergoes considerable changes.
First, when σ η ∼ 1, the energy spectrum D(E) changes radically (it is clearly seen in Figure 7): the curves of d 2 Ψ/dE 2 has a two-humped form at σ η < 1 and with σ η > 1 become simply convex. Moreover, the ground-state magnetization changes to zero when σ η > 1. It means that when the threshold value η = √ 3 is surpassed, the ground-state configuration goes off the initial state by distance of 1 2 N in the Hamming's terms. In other words, the system undergoes a zero-temperature phase transition. The transition is followed by the change of the ground-state energy from E 0 = −2J to asymptotic value (26).
Second, the experimental relation between the critical temperature and noise divergence differs greatly from the well-known [8] expression kT c = 2 9 ∑ α J 2 iα 1/2 , which in our terms takes the form: We can see that the classical theory predicts that β c should fall with the growing deviation of noise. Moreover, Expression (28) predicts finite values of β c for any large η. The experiment yields the opposite result: in accordance with (24) β c grows in proportion with σ 2 η . The experiment also shows that β c grows with η and when η → √ 3 it reaches its maximum β c = 0.625, the phase transition disappears at η > √ 3 (σ η > 1). It can be said conceptually that when the threshold value η = √ 3 is surpassed, the jump T c → 0 occurs.
In our opinion, the difference between the experiment and theoretical predictions is due to the finite dimension of the system. First, the finite system is ergodic and even at low temperatures does not have spontaneous magnetization, which can be tested easily with the help of Monte-Carlo algorithm. Second, the self-averaging principle used for building the theory for N → ∞ is not realizable for finite N. Additionally, the use of terms "critical temperature" and "phase transition" is not quite correct in description of finite-dimension systems. For our estimates, we use approximate Expressions (5) and (6), which are valid for a special case of the free boundaries conditions and finite L. More general and more accurate estimates can be obtained using the results of papers [45,46], where the authors analyzed the Ising random-bond model with a tunable fraction of negative bonds and the paper [47], where the finite size of the lattice was taken into account accurately.
Finite-dimension grids are of interest in image processing and machine learning. In our paper, the grid dimensions were N = L × L with L = 25 ÷ 1000. If we consider a planar grid as a model of a flat pixel image, such dimensions are very popular. The main conclusion that can be drawn from our results is that the learning algorithms based on free energy optimization are temperature insensitive in the most popular condition of η >> 1 because there is no observable phase transition in this case.
Author Contributions: Authors contributed equally. All authors participated in the design of the survey, its realization, and in the writing of the manuscript.
Funding: This research received no external funding.