Non-Hermitian Localization in Biological Networks
NNon-Hermitian Localization in Biological Networks
Ariel Amir, Naomichi Hatano, and David R. Nelson
1, 3 School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA Institute of Industrial Science, University of Tokyo, Komaba, Meguro, Tokyo 153-8505, Japan Department of Physics, Harvard University, Cambridge, MA 02138, USA
We explore the spectra and localization properties of the N -site banded one-dimensional non-Hermitian random matrices that arise naturally in sparse neural networks. Approximately equalnumbers of random excitatory and inhibitory connections lead to spatially localized eigenfunctions,and an intricate eigenvalue spectrum in the complex plane that controls the spontaneous activity andinduced response. A finite fraction of the eigenvalues condense onto the real or imaginary axes. Forlarge N , the spectrum has remarkable symmetries not only with respect to reflections across the realand imaginary axes, but also with respect to 90 ◦ rotations, with an unusual anisotropic divergencein the localization length near the origin. When chains with periodic boundary conditions becomedirected, with a systematic directional bias superimposed on the randomness, a hole centered onthe origin opens up in the density-of-states in the complex plane. All states are extended on therim of this hole, while the localized eigenvalues outside the hole are unchanged. The bias dependentshape of this hole tracks the bias independent contours of constant localization length. We treat thelarge- N limit by a combination of direct numerical diagonalization and using transfer matrices, anapproach that allows us to exploit an electrostatic analogy connecting the “charges” embodied inthe eigenvalue distribution with the contours of constant localization length. We show that similarresults are obtained for more realistic neural networks that obey “Dale’s Law” (each site is purelyexcitatory or inhibitory), and conclude with perturbation theory results that describe the limit oflarge bias g, when all states are extended. Related problems arise in random ecological networksand in chains of artificial cells with randomly coupled gene expression patterns. PACS numbers: 87.18.Sn; 02.10.Yn; 63.20.Pw
I. INTRODUCTION
The simplest models of neural networks assumelong-range connectivity between individual neurons inthe brain, leading to synaptic matrices M ( i, j ) with con-nection strengths approximately independent of the sep-aration r ij = | (cid:126)r i − (cid:126)r j | in three dimensions. The eigen-value spectrum of M ( i, j ) controls the spontaneous ac-tivity and induced response of the network, and muchis known when its elements are chosen from simple ran-dom matrix ensembles in the limit of large matrix rank N. For example, classic treatments of the spectra of realsymmetric random matrices leading to the Wigner semi-circular density-of-states describing the distribution ofreal eigenvalues have been generalized by Sommers etal . to allow for a tunable asymmetry in Gaussian prob-ability distributions for the matrix elements M ( i, j ) and M ( j, i ). These authors introduce a parameter that inter-polates between the Hermitian limit ( M ( i, j ) = M ( j, i ))studied by Wigner , and the case of fully asymmetric ma-trices where M ( i, j ) and M ( j, i ) are independent ran-dom variables. In the latter, non-Hermitian limit, thesemi-circular eigenvalue distribution on the real axis isreplaced by the “Circular Law” , where the eigenval-ues are now uniformly distributed inside a circle in thecomplex plane, with a vanishing fraction lying outsidethe circle in the limit N → ∞ . For the general case,Sommers et al . found that the eigenvalue distribution isuniform inside an ellipse , whose aspect ratio along thereal and imaginary axes varies with the amount of non- Hermiticity .As pointed out by Rajan and Abbott , typical ap-plications to neuroscience require that each node in asynaptic conductivity network be either purely excitatoryor inhibitory (Dale’s Law), which leads to constraintson the signs of the matrix elements M ( i, j ): all entriesin a row describing an excitatory neuron must be posi-tive or zero, and all entries in an inhibitory row must benegative or zero. These authors then studied eigenvaluespectra of random matrices with long range connectiv-ity, with excitatory and inhibitory networks drawn fromdistributions with different means and with equal or dif-ferent standard deviations. When the strengths of theexcitatory and inhibitory connections are appropriatelybalanced, with equal standard deviations, the eigenvaluedistributions can be made to obey the Circular Law byimposing a mild constraint. However, when the standarddeviations are different, the eigenvalue density becomesnon-uniform within a circle in the complex plane.Less is known for N -site banded random matriceswith signed matrix elements, which might be an approx-imate model for neural networks such that M ( i, j ) ∼ exp[ −| (cid:126)r i − (cid:126)r j | /l ], where l d (in d-dimensions) is often alarge volume containing as many as 10 neurons. On spa-tial scales larger than l , the synaptic connectivity matrixbecomes sparse, with the largest elements concentratedalong the diagonal. Banded Hermitian random matricesin d dimensions, frequently studied in the context of solid-state physics, have long been known to have eigenvaluespectra characterized by a large number of spatially lo-calized eigenfunctions , and it is this phenomenon that a r X i v : . [ c ond - m a t . d i s - nn ] J a n we wish to study here. To focus on an extreme exam-ple of bandedness, consider a matrix describing a one-dimensional chain of N sites, where only the elements M ( j, j ), M ( j, j + 1) and M ( j + 1 , j ) describing on-siteand nearest-neighbor couplings can be nonzero. If wewish to impose periodic boundary conditions, we will set M ( i, j ) = M ( i, j + N ) = M ( i + N, j ) , ∀ i, j . If the latticespacing a ≈ l , this model is a rough approximation to thedenser neural networks discussed above, coarse-grainedout to a scale of order l , with each site representing thespatially averaged firing rates of many actual neurons.We concentrate here on off-diagonal randomness, and as-sume that all M ( j, j ) are identical, and describe, say,a site-independent damping to a background firing rate.For the Hermitian case, with M ( j, j + 1) = M ( j + 1 , j )chosen from some probability distribution, nearly allstates are localized in the limit of large N , with thelongest localization lengths occurring near the band cen-ter and the shortest localization lengths near the bandedges . See Appendix A for a brief review and numeri-cal illustration of this solid-state physics example, whichprovides a useful benchmark for the more intricate prob-lem with complex eigenvalues we study here. Chaudhari et al . have studied a related problem, with Hermitiancoupling strengths falling off exponentially in space andrandom self-couplings (diagonal randomness), in the con-text of one-dimensional neural networks, as well as lo-calization of the eigenmodes in a non-Hermitian matrixarising not from disorder but from a slow gradient in thediagonal elements. Here, we study sparse non-Hermitian matrices and the localization properties of their eigen-modes. An important feature of our model is the under-lying spatial structure (the connections are between near-est neighbors in real space), which distinguishes our workfrom recent, interesting studies of sparse non-Hermitianmatrices without such structure . A. From neural networks to random matrices
As stated above, we focus here on off-diagonal ran-domness in the neural connections, which is both non-Hermitian ( M ( j, j + 1) (cid:54) = M ( j + 1 , j )) and, importantly,also allows for M ( j, j + 1) and M ( j + 1 , j ) to be of op-posite sign roughly 50% of the time. We thus model aset of approximately balanced excitatory and inhibitorynearest-neighbor neural connections in one dimension,and study the localization properties of the intricate com-plex eigenvalue spectrum that results. To put our inves-tigations in context, consider first (using a convenientDirac notation | j (cid:105) to describe a neuron at site j ) thespectrum of a simple Hermitian
1d tridiagonal matrix S e n s o r y I npu t s P r o ce ss i n g FIG. 1. Schematic of a 1d neural network problem with pe-riodic boundary conditions. Sensory inputs, possibly after aprocessing step, are sent via feed-forward couplings into a cir-cular ring of N neurons | j (cid:105) , j = 1 , ...N , with nearest-neighborexcitatory and inhibitory connections. M (1 ,
2) and M (2 , s . with random connections, namely M = − N (cid:88) j =1 (cid:2) s + j | j (cid:105)(cid:104) j + 1 | + s − j | j + 1 (cid:105)(cid:104) j | (cid:3) ,s + j = s − j = s + s j > ,s j ∈ [ − ∆ , ∆] , ∆ = 0 . s (1)Here, all eigenvalues are real, and the symmetrical con-nections s + j = s − j between neighboring sites are guar-anteed to be positive by our choice of a relatively nar-row (∆ < s ) box distribution for the bond-to-bondfluctuations in the connection strengths relative to thebackground level s , and we have subtracted off a diag-onal contribution, assumed to be site-independent. Asshown in Appendix A, the localization length ξ ( ε ) of theeigenfunctions diverges near the band center at energy ε = 0. The quantity ξ ( ε ) describes the spatial scale overwhich an eigenfunction with energy ε is nonzero. If theeigenfunction φ ε ( j ) is large near a “center of localiza-tion” j ∗ , then roughly speaking its envelope decays likeexp[ −| j − j ∗ | /ξ ( ε )]. The localization length ξ ( ε ) is knownto diverge logarithmically , ξ ( ε ) ∼ log(1 /ε ), as ε → ρ ( ε ) to the localizationlength ξ ( ε ), known as the Thouless relation . In thiscase, the Thouless relation implies a strongly divergingdensity-of-states, ρ ( ε ) ∼ / [ | ε | log (1 /ε )], near the ori-gin. We shall see echoes of these results later in thispaper.We study here a generalization of Eq. (1) thatarises in one-dimensional neural networks with randomexcitatory and inhibitory nearest-neighbor connections.Following Chapter 7 of Ref. [16], consider the sparse“recurrent neural network” shown in Fig. 1, a chainof nodes with asymmetric connections between nearestneighbors and with periodic boundary conditions. Sen-sory inputs, possibly after a processing step, are sent viafeed-forward couplings into a circular ring of N neurons | j (cid:105) , j = 1 , ..., N . The nearest-neighbor excitatory and in-hibitory couplings M ( j, j + 1) and M ( j + 1 , j ) can notonly be unequal, but also of opposite sign, if one directionis excitatory and the other inhibitory. Consider a modelwhere the average firing rates v i and v j in neurons i and j (a coarse-grained description of the temporal densityof discrete spikes in these neurons) are coupled together,and obey τ dv i dt = − v i + F ( M ij v j + h i ) . (2)Here M ij ≡ M ( i, j ), τ is a characteristic neuron timeconstant (assumed for simplicity to be the same for allneurons) and we use the summation convention. Inputsto an animal brain from the outside, due to whiskers,retinal cells, olfaction, etc . (after a possible processingstep), are represented by h i = W ij u j , where the con-nection matrix W ij ≡ W ( i, j ) and the input firing rates { u i } represent the feed-forward part of this network. Theactivation function F ( w ) (often taken to have a nonlin-ear sigmoidal shape ) insures that the firing rates arebounded above (when inhibitory connections are present,additional constraints can be imposed to insure that thefiring rates can never be negative ). Here we assumefor simplicity that the activation function is the same forboth excitatory and inhibitory connections. The eigen-values and eigenvectors of the matrix M ij are clearlyimportant for understanding the behavior of a linearizedversion of Eq. (2), τ dv i dt = − v i + M ij v j + h i , (3)where we assume without loss of generality that F (0) = 0and F (cid:48) (0) = 1. This linear recurrent network is capableof both selective amplification and input integration .More generally, knowledge of the eigenvalues and eigen-functions of M ij is useful for studying spontaneous activ-ity and evoked responses . Spontaneous activity de-pends on whether the real parts of any of the eigenvaluesare large enough to destabilize the silent state in a lin-ear analysis, and the spectrum of eigenvalues with largereal parts provides valuable information about the spon-taneous activity in the full, nonlinear models, and aboutthe localization volume determining the size of the ac-tive clusters carrying out computations. Moreover, simi-lar matrices arise when nonlinear problems are linearizedabout a steady state.To see why random neural connections might berelevant, note that these can be formed during develop- ment, with many random attachments of axons and den-drites to other neurons. Then, over time, pruning (lossof connections) and adaptation (strengthening and weak-ening of various excitatory and inhibitory connections)occur as neural circuits “learn” various functions. Thelikely result is a mixture of structured and random com-ponents. The spectra and eigenfunctions of completelyrandom sparse neural network chains, with a mixture ofinhibitory and excitatory connections, could provide adescription of neural activity during the early stages ofdevelopment, and is also a useful reference model. Sim-ilar justifications have been advanced for studying thedense neural networks that obey Dale’s law treated inRef. [7]. B. Model and density-of-states
With this motivation, we now discuss the spectra ofnon-Hermitian matrices that generalize Eq. (1), namely M = N (cid:88) j =1 (cid:2) s + j e g | j + 1 (cid:105)(cid:104) j | + s − j e − g | j (cid:105)(cid:104) j + 1 | (cid:3) , (4)where for most of this paper we impose periodic bound-ary connections, | j + N (cid:105) = | j (cid:105) . The constant diagonalcontribution associated with Eq. (3) has again been sub-tracted off. The connection strengths s + j and s − j are in-dependent and identically distributed random variableschosen from a probability distribution P ( s ± j ), given by(see inset to Fig. 1), P ( s )= f − u for u < s < , − f − u for − < s < − u, , otherwise . (5)The parameter u , 0 < u ≤
1, controls the width of thepositive and negative parts of the distribution, while f ,0 < f <
1, determines the ratio of inhibitory to exci-tatory connections. This functional form excludes con-nections that are very close to zero, which would biasthe 1d network towards falling apart into disjoint pieces.The coupling g in Eq. (4) controls the strength of a sys-tematic clockwise ( g >
0) or counterclockwise ( g < g canhave a remarkable effect on the spectrum and localizationproperties. In this paper, we concentrate on the spectraand localization properties of eigenfunctions in the ap-proximately balanced case, f ≈ /
2, which representsthe greatest departure from conventional Hermitian lo-calization problems in one dimension . For now,we suppress the neuroscience constraints associated withDale’s Law, as might be appropriate if each node in thechain describes a large number of strongly coupled neu-rons randomly chosen to be excitatory and inhibitory.However, we shall later argue (Sec. II C) that a straight-forward modification of Eq. (4) that respects Dale’s Lawproduces negligible changes in the spectra and localiza-tion properties in the limit of large N .The case of f = 1 with 0 < u < and in the population dynamics of heteroge-neous 1d environments . When g >
0, this problem issometimes referred to as “directed localization” , ter-minology we adopt here as well. The spectrum of mod-els with f = 1 / u = 1 ( i.e ., M ( j, j + 1) = ± M ( j + 1 , j ) = ±
1, excitatory/inhibitory connec-tions chosen at random with equal probabilities) has alsobeen studied before , and has been shown to havean extremely rich structure. Here, we explore the local-ization properties of the eigenfunctions associated withthese spectra for a range of u values in some detail.Figures 2(a) and (b) exhibit the remarkable spectraassociated with Eqs. (4) and (5) when f = 1 / u = 1and for g = 0 and g = 0 . , whomentioned that this model might have interesting local-ization properties. Although the eigenvalues are in gen-eral complex, when g = 0 a significant fraction of them(about 20%) have condensed onto the real axis, see Fig. 3;similarly, about 20% have condensed onto the imaginaryaxis. In Appendix B this numerical analysis is extendedto the case of u (cid:54) = 1, with similar results. For large N ,the density-of-states is symmetric under reflections acrossthe real and imaginary axes, as well as across ± ◦ linesin the complex plane, as we shall show in Sec. II. Theremaining eigenvalues (approximately 60%) form an in-tricate, diamond-shaped structure. When u is near 1, thedensity-of-states appears to acquire a fractal-like bound-ary. See Appendix B for a summary of the density-of-states for the more general probability distribution ofEq. (5) for arbitrary u and f = 1 / C. Main Results
We are now in a position to summarize our mainresults. In Sec. II we discuss various symmetries asso-ciated with the density-of-states of the models studiedhere. Sec. III we show that almost all eigenfunctions arelocalized (similar to the 1d Hermitian case of AppendixA), with the smallest localization lengths near the bound-ary of the spectrum in the complex plane, and a diverginglocalization length near the origin. Our analysis of local-ization in this model has been guided by work of Derrida et al . on a related problem (with unimodular complex random couplings between sites), who derive an elegant − − − − (a) − − − − (b) FIG. 2. (a) Density-of-states (DOS) of the complex eigenvaluespectrum for the sparse random matrix defined by f = and u = 1 with g = 0 obtained through exact diagonalizationof 10,000 matrices of dimension 5000 × M ( j, j + 1) = ± M ( j + 1 , j ) = ±
1. (b)Density-of-states for the case g = 0 .
1, all other parametersare identical to (a). generalization of the Thouless formula for eigenvalues inthe complex plane: The inverse localization length is thetwo-dimensional electrostatic potential associated with acollection of charges at the eigenvalue locations in thecomplex plane. Our numerical analysis strongly suggeststhat the localization length diverges as the modulus ofthe eigenvalues tends to zero. Indeed, if the eigenvaluesare written λ = λ + iλ , where λ and λ are real, anumerical study of the inverse localization length definedvia the product of N random 2 × − − − − FIG. 3. DOS on the real axis for the parameters of Fig. 2(a),comprising close to 20% of all eigenvalues. for u near 1 leads to the following ansatz: ξ ( λ , λ ) ∝ | λ | + | λ | ) (cid:112) λ + λ , (6)which should be compared to the much weaker logarith-mic divergence discussed in Appendix A for 1d Hermitianhopping randomness. Although the divergence shownin Eq. (6) only holds for u near 1, the infinite localiza-tion length at the origin is more general, as discussed inSec. III E.As shown in Fig. 2(b), a hole surrounding the originwith angular corners opens up in the complex plane whenthese calculations are repeated for a clockwise bias pa-rameter g = 0 .
1. A similar hole opens up in the Feinberg-Zee model of Ref. [26], with complex unimodular hoppingmatrix elements . As we demonstrate in Sec. III, a largenumber of extended states now occupies the rim of hole.For large N , the eigenvalues of the localized states out-side the hole are unchanged, a spectral rigidity propertythat can be derived from a simple exponential “gaugetransformation” acting on the corresponding eigenfunc-tions for g = 0, see Ref. [20]. A corollary is that the g -dependent shape of holes like that in Fig. 2(b) tracksthe g = 0 contours of constant localization length, witha diverging localization length as the rim is approachedfrom the outside.What happens to directed localization for neuralnetworks that obey Dale’s Law, as studied for spatiallyextended neural connections in Ref. [7]? In Sec. II C,we argue that the above results should be unchanged forlarge N . We will do this by replacing the matrix M witha modified connectivity matrix, G = N (cid:88) k =1 σ k (cid:2) e g | k + 1 (cid:105)(cid:104) k | + e − g | k − (cid:105)(cid:104) k | (cid:3) , | k + N (cid:105) = | k (cid:105) (7) where the N real random numbers { σ k } are chosenfrom the probability distribution P ( σ k ), again given byEq. (5). Figure 7 illustrates a particular example of theDale’s Law connectivity matrix for N = 5. Note thatthe nonzero connections in the same row have the samesign. Equation (7) has site randomness, as opposed tothe bond randomness displayed in Eq. (4). Despite thefact that 2 N random numbers are necessary to specify M and only N random numbers specify G , we show viasimilarity transformations in Sec. II C that the spectraand localization properties of M and G are essentiallyidentical, a result which we have also confirmed numeri-cally. The underlying reason is that the spectral proper-ties are determined in both cases by above/below diag-onal products such as M ( j, j + 1) · M ( j + 1 , j ) = s + j s − j and G ( j, j + 1) · G ( j + 1 , j ) = σ j σ j +1 . These quantitieshave identical statistical properties.In Sec. IV we discuss large- g perturbation theory,that focuses on the changes in the eigenvalues and eigen-functions which are all extended in this limit. This anal-ysis can be carried out for arbitrary f and u , although itcannot capture the localization that results as the eigen-values move toward the origin with decreasing g .Section V contains a summary and outlook, includ-ing a brief discussion of the effect of diagonal randomness.Appendices A-C describe, respectively, a Hermitian ran-dom hopping model, the density-of-states for arbitrary u , and second-order perturbation theory for large g . II. SYMMETRIES ASSOCIATED WITH THEDENSITY-OF-STATES
In order to discuss spectral symmetries, we first in-troduce a similarity transformation which is applicable tothe present model with open boundary conditions. Con-sider an N × N tridiagonal matrix of the form A = d b a d b a d b . . . . . . . . . a N − d N − b N − a N − d N , (8)where { a x } , { b x } and { d x } are arbitrary real numbers,and all other entries vanish. We can symmetrize thismatrix by a diagonal similarity transformation, whose j th matrix element reads S jj = j − (cid:89) k =1 (cid:114) a k b k , (9)which we may call a generalized gauge transformation .The result of this symmetrization reads A (cid:48) = S − AS = d c c d c c d c . . . . . . . . . c N − d N − c N − c N − d N , (10)where c j = (cid:112) a j b j . (11)The matrices (8) and (10) are isospectral, due to theproperties of similarity transformations.Note that the spectrum of the matrix A (cid:48) dependsonly on the product of the opposing off-diagonal elements a j and b j , not independently on each of them (as alsofollows from calculating the characteristic polynomial ofthe matrix A ). Another important observation is thatthe matrix (10) is real and symmetric if a j and b j havethe same sign, but is non-Hermitian otherwise. The non-Hermitian Anderson chains proposed in Refs. [20 and38], in which all a j and b j are negative, therefore wouldhave only real eigenvalues unless we introduced periodicboundary conditions. If the matrix (8) has non-zero cor-ner elements a N for A N and b N for A N (thus couplingthe chain into a ring), the resulting matrix A (cid:48) = S − AS has non-zero corner matrix elements A (cid:48) N = (cid:0) S − AS (cid:1) N = (cid:112) a N b N N (cid:89) j =1 (cid:114) a j b j , (12) A (cid:48) N = (cid:0) S − AS (cid:1) N = (cid:112) a N b N N (cid:89) j =1 (cid:115) b j a j , (13)which make the matrix non-Hermitian and allows thepossibility of complex eigenvalues unless N (cid:89) j =1 a j = N (cid:89) j =1 b j . (14)Although this similarity transformation leaves thediagonal randomness intact, it packs all effects of therandom, non-Hermitian hopping terms into a single pairof corner matrix elements. This perspective is useful al-ready for simple cases where the elements { a j } and { b j } in Eq. (8) can be different, but are both constrained tobe of the same sign. Let us take a j = e − g s − j > b j = e g s + j >
0, consistent with Eq. (4), so that the corner matrix element A (cid:48) N takes the form A (cid:48) N = (cid:113) s − N s + N e Ng N (cid:89) j =1 (cid:118)(cid:117)(cid:117)(cid:116) s + j s − j . (15)If we now choose the elements (cid:8) s ± j (cid:9) according to theprobability distribution of Eq. (5) with f = 1 and0 < u <
1, all s ± j will be positive, and A (cid:48) N is real anddescribed by a log-normal distribution. It is simpler tostudy log A (cid:48) N , which behaves like a random walk. Asdiscussed in Sec. III, a closely related quantity determinesthe localization properties of eigenfunctions as functionof g . Focusing for simplicity on the case u = 0 + , wereadily find (cid:104) log A (cid:48) N (cid:105) = N g − , (16)and (cid:68) (cid:0) log A (cid:48) N − (cid:10) log A (cid:48) N (cid:11)(cid:1) (cid:69) = N O (1 /N ) , (17)where (cid:104)•(cid:105) represents an average over the disorder andsimilar results obtain for 0 < u <
1. Upon defining aneffective directional bias parameter g eff ≡ (cid:10) log A (cid:48) N (cid:11) /N ,we see that if microscopic bias is g = 0, then the hoppingrandomness represented by the elements (cid:8) s ± j (cid:9) leads to a g eff = O (1 /N / ), which vanishes in the limit large N .Thus, the hopping disorder is effectively undirected as N → ∞ in this case. When diagonal randomness is alsopresent, we expect that the localized states will remainlocalized with real eigenvalues, unless g exceeds a criti-cal value g c given by the minimum inverse localizationlength when g = 0.If a j and b j can have different signs, the matrix (8)is inherently non-Hermitian and can have complex eigen-values with or without periodic boundary terms. It isthis interesting case we focus on in the present paper.As discussed below, coupling the chain into a ring is cru-cial when g > A. Spectrum of sign-random model
Let us apply the above considerations to the sign-random non-Hermitian tight-binding chain given by the N × N matrix corresponding to Eq. (4) with g = 0: M = s +1 s − s +2 s − s +3 s − . . .. . . . . .. . . s + N − s − N − , (18) − − − − (a) − − − − (b) − − − − (c) FIG. 4. The spectra on complex planes of the matrix (18) for(a) N = 100, (b) N = 1000 and (c) N = 10 , where all remaining elements, including diagonal ones,vanish, and each of s ± j are randomly set to be ± /
2. Spectra found by numerical diagonaliza-tion of a random sample are shown in Fig. 4 for N = 10,1000 and 10 , g = 0 we can neglect the corner ma-trix elements.We can describe the symmetries of the spectrum inthe following way: First, since the matrix (18) is real,namely, M = M ∗ , if there is an eigenvalue λ n , theremust be another eigenvalue λ ∗ n . In other words, the spec-trum is symmetric with respect to reflections across thereal axis. Second, since the spectrum depends only onthe product of s + x and s − x , the matrices H and − H areisospectral, and hence if there is an eigenvalue λ n , theremust be another eigenvalue − λ n . In other words, thespectrum is symmetric under inversion in the complexplane λ → − λ . Upon combining the two symmetries,we see that the matrices M and − M ∗ are isospectral, − − − − FIG. 5. The spectrum on a complex plane of the matrix (18),this time with additional boundary elements; compare it withFig. 4(a). The system size is N = 100. and hence the spectrum is symmetric with respect to re-flections around the imaginary axis, too. This argumentholds in a statistical sense even if we add zero-mean di-agonal randomness into Eq. (18).Finally, we argue that the spectrum has statisticalsymmetry with respect to the reflections across the 45 ◦ lines Re E = ± Im E . According to the argument in thebeginning of this section, the spectrum depends only onwhether the product s + x s − x is +1 or −
1. In other words,the randomness of the matrix (18) is caused by indepen-dent probability distributions of N − { s + x s − x } , instead of 2 N −
1. Let us thenconsider the spectrum of the matrix i M . By multiplyingevery matrix element by i = √−
1, we flip the sign ofthe product of the opposing off-diagonal elements which,however, does not change the binomial distribution of the N − f = 1 /
2. There-fore, the matrices M and i M are statistically isospectral.Since the spectrum of i M is given by the 90 ◦ rotationof that of M , the spectrum is statistically symmetricwith respect to this operation. Combining this symmetrywith the other symmetries, we conclude that it is statis-tically symmetric with respect to reflections around the45 ◦ lines. The fact that the symmetry becomes betteras we increase the system size underlines the observationthat the symmetry is indeed statistical.Adding the boundary elements M N = s + N and M N = s − N does not change the spectrum in an es-sential way when g = 0; their first-order perturbationto an eigenvalue λ n with its normalized left- and right-eigenvectors (cid:104) ˜ ψ n | and | ψ n (cid:105) (we use the tilde symbol toemphasize that they are not Hermitian conjugate to eachother) is of order 1 /N at most (and is exponentially smallif the eigenfunctions are localized). Indeed, comparisonof numerical results of Figs. 4(a) and 5 with and with-out the boundary elements suggests that they are notonly statistically the same but also almost identical withoccasional differences, even for N = 100. B. Asymmetric amplitudes
Let us next introduce asymmetric amplitudes to thesign-random tight-binding model. Following Refs. [20,22, and 38], we express the asymmetry in the form (equiv-alent to Eq. (4)) M ( g ) = e g s +1 e − g s − N e − g s − e g s +2 e − g s − . . . . . .. . . . . .. . . e g s + N − e g s + N e − g s − N − , (19)where we assume g > e g s + N and e − g s − N ; if not, the spectrum would be g -independent be-cause it would depend only on the product of the oppos-ing off-diagonal elements. The diagonal similarity trans-formation T ( g ) xx = e − g ( x − (20)changes the matrix M ( g ) into M (cid:48) ( g ) = T ( g ) − H ( g ) T ( g ) = s +1 e − Ng s − N s − s +2 s − . . . . . . s + N − e Ng s + N s − N − , (21)which shows that the boundary elements are essential inhaving a strong dependence on g .As was discussed in Refs. [20 and 38], the spectrumof M ( g ) can in fact be an indicator of the localization ofthe eigenfunctions. Suppose that the eigenfunction ψ n ofan eigenvalue λ n of the original Hamiltonian M = M (0)is localized around a site x and behaves approximatelyas ψ n ( x ) ∼ e − κ n | x − x | (22)except for a phase factor. This quantity is also an ap-proximate eigenfunction of T ( g ) − M ( g ) T ( g ), becausethe first-order perturbative corrections due to the bound-ary elements are exponentially small, of order e − N ( κ n − g ) ,when g < κ n . Thus, the corresponding eigenfunction of M ( g ) is given by ψ n ( x ; g ) ∼ e − gx − κ n | x − x | (23) − − − − (a) − − − − (b) − − − − (c) FIG. 6. The spectra on complex planes of (a) the directedmatrix (19) with random signs for g = 0 . N = 1000;(b) the matrix (24) obeying Dale’s law for N = 1000; Notethe similarity to the spectrum in Fig. 4(b), that does not havethis restriction. (c) the matrix (25) for g = 0 . N = 1000.Note the close similarity with Fig. 6(a). except for a phase factor. Indeed, the periodic boundaryconditions are almost precisely satisfied for large N if g < κ n ; the discrepancy at the boundary is exponentiallysmall, of order e − N ( κ n − g ) . Therefore, the eigenvalue λ n of M (0) remains to be an eigenvalue of M ( g ) when g <κ n . This argument breaks down when g > κ n , for whichthe eigenvalue now moves as a function of g , with motionstarting when g = κ n . The numerical diagonalization ofa random sample with g = 0 . κ = 0 . g = 0, and vanishing κ for g = 0 . FIG. 7. Ring with N = 5 coupled neurons obeying Dale’sLaw: each neuron couples in a purely excitatory or purelyinhibitory manner to its two neighbors; solid arrows representpositive, excitatory connections, and dashed lines negative,inhibitory ones. The 5 × G corresponding to thisparticular choices of signs is indicated on the right, where onlythe sign on the nonzero matrix elements is indicated. Notethat non-zero connections in the same row have the same sign. C. Spectrum of models obeying Dale’s law
Figure 7 shows a network with N = 5 that respectsDale’s law, and the signs of the non-zero elements of thecorresponding matrix. We now argue that the resultspresented in this paper are readily extended also to thisscenario, which is more realistic for neural networks.To take this situation into account, we consider(taking g = 0 for now) G = σ σ σ σ σ σ . . .. . . . . .. . . σ N − σ N (24)instead of M in Eq. (18), where each of σ j randomlytakes ± /
2, although similar consid-erations apply to the more general probability distribu-tion of Eq. (5). The value of σ j indicates whether thetwo connections out of the j th neuron are excitatory orinhibitory.According to the previous argument, the spectrumdepends only on the product of opposing off-diagonal el-ements. In the case of the matrix (24), we can regard the N − { σ j σ j +1 = ± (cid:107) j = 1 , , . . . , N − } as independent random variables, just as for the ma-trix of Eq. (18) we can regard the N − { s + x s − x = ± (cid:107) x = 1 , , . . . , N − } as independent. There-fore, the matrices (18) and (24) are statistically isospec-tral; see Fig. 6(b) for the spectrum for one random sam-ple, obeying Dale’s law, to be compared with Fig. 4(b).The statistical isospectrality does not change muchwhen we introduce the boundary terms G N = σ N and G N = σ , because the perturbation of these terms tothe spectrum is of order 1 /N at most (and exponentially small if the states are localized). The only difference inthe statistics is the fact that the product of all super- andsub-diagonal elements of the matrix (18), including theboundary terms M N = s + N and H N = s − N , is randomand can take ±
1, but that of the matrix (24) is always+1. Finally, introduction of the asymmetry parameter g to G as in G ( g ) = e g σ e − g σ e − g σ e g σ e − g σ e g σ e − g σ . . .. . . . . .. . . e g σ N − e g σ N e − g σ N , (25)has the same effect as we showed in Sec. II B for the ma-trix (19); see Fig. 6(c) for a numerical illustration for onerandom sample. Note the close similarity with Fig. 6(a),for the matrix M ( g ), which is unconstrained by Dale’slaw. III. LOCALIZATION PROPERTIES
We now investigate the localization properties ofthe model via three different and complimentary routes:(i) By calculating the participation ratio of eigen-modes obtained via exact diagonalization;(ii) By using the transfer-matrix approach, and theequivalence between the Lyapunov exponents andthe inverse localization length;(iii) By numerically calculating the density-of-states(DOS) via exact diagonalization, and inferring thelocalization length via the Thouless relation , asgeneralized to localization with complex eigenval-ues in Ref. [33].We find analytically that for g = 0 and for any u that the localization length is infinite at λ = 0 ( i.e ., theinverse localization length κ vanishes at the origin), sug-gesting a diverging localization length as | λ | →
0. Sucha divergence is strongly supported by our numerical re-sults. Interestingly, we find that in contrast to the resultsof Ref. [33], the dependence κ on λ = x + iy in the vicin-ity of the origin is not isotropic. Through the Thoulessrelation, which we elaborate on below, we will show thatthis property is connected to the vanishing DOS at theorigin. In the following, we elaborate on the differentmethods and compare the results.0 A. Localization properties from numericaldiagonalization
A useful measure of the localization of an eigenvec-tor is its participation ratio, defined as P ≡ (cid:88) j | ψ j | / (cid:88) j | ψ j | . (26)Indeed, a perfectly localized eigenvector with supportonly at a single site would have P = 1, while a per-fectly delocalized one (with ψ j = 1 / √ N for every j ) has P = N . By averaging the participation ratio, or its in-verse, we may gain insights into the localization prop-erties of the system. Figure 8 shows the results of nu-merical diagonalization of 10,000 matrices of dimension5000 × f = 1 / u = 1.Figure 8(a) shows the tendency of states to be de-localized near the origin (vanishing inverse localizationlength (IPR)), becoming more localized away from theorigin. However, while this is a direct and straightfor-ward method, in the next subsections III B and III C, wewill find the localization length more accurately; we willshow that while it diverges near the origin, it does nothave a radial symmetry, and only achieves radial symme-try away from the origin.Upon repeating the analysis for g = 0 . g .We now comment briefly on the effect of periodicboundary vs . open conditions. The arguments given inSec. II suggest (and numerical diagonalizations confirm)that the g = 0 spectrum is nearly identical when peri-odic instead of open boundary conditions are employedin Fig. 8(a). In contrast, the hole and extended states inthe g = 0 . s + N → , s − N →
0, which breaksthe chain. Nevertheless the hole has a physical inter-pretation even for open chains: Although all eigenvaluesretain their g = 0 values, eigenfunctions inside the holebecome edge states, piled up on one side of the brokenchain. B. Transfer matrix approach
A well-established method for finding the localiza-tion length of a one-dimensional system calculates theLyapunov exponent via the transfer matrix technique .If ψ n is the eigenfunction amplitude on the n th site, the − − − − − − − − − − − − (a) − − − − − − − − − − − − (b) FIG. 8. (a) Inverse participation ratio (IPR) as a function ofeigenvalue, obtained via numerical diagonalization of 10,000matrices of dimension 5000 × g = 0 and periodicboundary conditions. With g = 0, the results with openboundary conditions would be nearly identical. The colorbars indicate the inverse localization length on a logarithmicscale. Note that the background (where there are no states) iswhite; the fractal nature of the spectrum implies that the IPRis not evaluated everywhere, but only on the support of theDOS. (b) IPR for the same parameters, but with g = 0 . × (cid:18) ψ n ψ n +1 (cid:19) tothe vector (cid:18) ψ n − ψ n (cid:19) with eigenvalue λ is given by T n = (cid:18) − e − g s − n − /s + n λe − g /s + n (cid:19) , (27)where we do not include diagonal disorder and s − n − and s + n are independent random variables representing the1off-diagonal randomness.The Lyapunov exponent can be extracted by takingthe limit κ ≡ lim N →∞ (cid:104) log( || T N · T N − ... T || ) (cid:105) /N, (28)where || .. || denotes the norm of the matrix, and (cid:104) .. (cid:105) en-semble averaging over the quenched disorder. It can beproven that under quite general conditions the limit ex-ists, and κ equals the inverse of the localization length ,which we identify (up to constants of order unity) withthe inverse participation ratio of Sec. IIIA.This procedure provides a numerically attractiveroute to finding the localization length, without havingto diagonalize large matrices. However, in practice N has to be large in order for the method to be accurate,which implies that the product will result in a matrixwith a large norm, imposing computational difficulties.We resolved this problem by working with the recursiverelation for the quantity r n ≡ ψ n +1 /ψ n (note that unlikeRef. [33], in our definition r is a complex number). FromEq. (27) we immediately find that r n +1 = − ( s − n − /s + n ) e − g /r n + λe − g /s + n . (29)In this case, the values of r n are well-behaved alsofor large n , leading to robust numerics. Upon evaluationof r ...r N = ψ n +1 /ψ , the Lyapunov exponent can befound in a similar fashion as κ ( λ ) = 1 N N (cid:88) j =1 log | r j | . (30)It is beneficial to omit the values of r j at the beginningof the sequence, to reduce the effects of the initial condi-tions, though in the limit of large N this is not strictlynecessary.Using this method, we obtained Fig. 9, which cor-roborates and complements the results of the exact nu-merical diagonalization. While Fig. 8 is computationallyexpensive, generating Fig. 9 takes several minutes on aPC, a testimony to the power of this technique. Note,however, that Eqs. (29) and (30) always deliver a valuefor κ ( λ ), regardless of whether there is actually a nor-malizable eigenfunction at that particular value of λ . C. Connection to the density-of-states via theThouless relation
A classic result in the theory of Anderson localiza-tion in one dimension is an elegant relation connectingthe density-of-states to the localization length, due toThouless . This relation can readily be generalized tothe non-Hermitian case , where it states that ∇ κ ( x, y ) = ρ ( x, y ) . (31) − − − − (a) − − − − (b) FIG. 9. (a) Inverse localization length κ ( λ ) as a functionof eigenvalue, obtained via Eq. (30) and the recursion equa-tion (29). Each point is obtained via 10,000 iterations. Inthis case κ does not take negative values anywhere. (b) Thecolored map shows κ ( λ ) for the same parameters, but with g = 0 .
1. There is a finite region in the vicinity of the originwith negative values of κ , which was given a white color (notincluded in the color map). This region corresponds to thegap of Fig. 8(b); the black dots superimposed on the plot arethe result of the diagonalization of a single 1000 × g = 0 .
1. Note that the boundary of the white holecorresponds almost exactly to the rim of the extended statesfound via the exact numerical diagonalization.
Here, the complex eigenvalue is λ = x + iy , and ρ ( x, y )is the density-of-states. This equation can be inverted,using the well-known analogy with 2d electrostatics,whereby ρ ( x, y ) represents a collection of infinite, chargedwires, perpendicular to the complex plane, each associ-ated with a logarithmic potential. Therefore we have κ ( x, y ) = (cid:90) ρ ( x (cid:48) , y (cid:48) ) log( | r − r (cid:48) | ) dxdy + C. (32)2 − − − − FIG. 10. The Lyapunov exponent κ ( λ ) was extracted fromthe diagonalization of a single 1000 × g = 0,using the electrostatic relation (32). The result is very similarto Fig. 9(a), which was obtained via the Ricatti recursionrelations. In the case g = 0, the constant is given by C = (cid:104) log( | s + j | ) (cid:105) = (cid:104) log( | s − j | ) (cid:105) , (33) i.e ., the average of the logarithm of the random matrixelements. Hence, in the case we are focusing on where | s + j | = 1, we find that C = 0. In the next section weshall show how the results for κ ( x, y ) for finite g can bemapped to the g = 0 behavior, which will show that inthe more general case we have C = 12 (cid:104) log( | s − j /s + j | ) (cid:105) = − g, (34)which follows from Eq. (37).This remarkable relation allows us to go back-and-forth between the two very different numerical proce-dures: obtaining κ via the recursion relation and ob-taining ρ via exact numerical diagonalization. Indeed,using a single realization of 1000 × g = 0. D. Hole in spectrum corresponds to contours ofLyapunov exponent
Consider the recursion relation of Eq. (29). It iseasy to “gauge away” the effect of g , by making the trans-formation y n ≡ r n e g , (35) − − − − FIG. 11. Constant κ contours of the heatmap of Fig. 9(a). upon which the equation takes the form y n +1 = − (cid:32) s − n − s + n (cid:33) /y n + λ/s + n . (36)This representation implies that for any complex eigen-value λ = x + iy , the effect of g is to decrease the Lya-punov exponent by an amount g : κ ( x, y ; g ) = κ ( x, y ; 0) − g (37)Hence, consistent with the gauge transformation resultof Eq. (23), for any g > κ < g will acquire a negative Lyapunov exponent. Sinceall states must be normalizable, the region with negative κ will not support any eigenfunctions, and corresponds tothe “hole” or gap seen in Figs. 2 and 6. This argumentimplies that the hole boundary corresponds to contourwhere κ = g , consistent with Fig. 9(b), where the resultsof exact diagonalization are superimposed on top of theLyapunov exponent heatmap.In Fig. 11, the contours of constant κ are shown. Asexpected from the electrostatic Thouless relation, awayfrom the effective support of the DOS the contours be-come circular, κ ( x, y ) ∝ log( x + y ), since all “charges”associated with the complex eigenvalues act as if theywere concentrated at the origin. Close to the origin thecontours obtain a roughly diamond or octagonal shape.This behavior is consistent with the vanishing DOS sug-gested by Fig. 2(a); via the Thouless relation, Eq. (31),if κ had a perturbative expansion such as κ ∼ x + y (see Ref. [33] for such result in a related model), thenthe DOS at the origin would have a finite, non-vanishingvalue.Furthermore, Fig. 12 shows the results for the y component of the gradient of the Lyapunov exponent,suggesting a δ -function contribution to the DOS alongthe x axis. Similar results can be obtained for the y axis.3 − − − − − − − − − FIG. 12. The y component of the gradient of the Lyapunovexponent. The abrupt change along the x is consistent witha condensation of eigenvalues along the real axis. For u = 1, the strength of the singular DOS along the x and y axis decays linearly close to the origin, as shownin Fig. 3. These results lead us to the following ansatzfor the behavior of κ near the origin for u = 1: κ ( x, y ) ∼ ( | x | + | y | ) (cid:112) x + y , (38) i.e ., it is a product of L and L norms. This ansatzis consistent with the eigenvalue condensations onto the x and y axis, and their linear density shown in Fig. 3.When appropriate higher-order cubic terms are addedto the ansatz of Eq. (38), the function becomes har-monic away from the x and y axes ( e.g .: by replacing d = (cid:112) x + y with [1 − e − d ]), consistent with the van-ishing DOS at the origin. The good agreement of theequipotential contours of the ansatz of Eq. (38) and ofthe numerically evaluated Lyapunov exponent is shownin Fig. 13. E. Vanishing of the Lyapunov exponent at theorigin
It is easy to see that for any distribution of thehopping matrix elements, the Lyapunov exponent mustvanish at the origin (in contrast to the behavior of mod-els with additional diagonal disorder ). To see this, con-sider the transfer matrix Eq. (27) for λ = 0. The productof two adjacent transfer matrices is in this case diagonal: S n = T n T n − = e − g − s − n − s + n − s − n − s + n − (39)with the elements on the diagonal being the ratio of tworandom variables. Therefore the Lyapunov exponent is NumericalAnalytic Ansatz − − − − FIG. 13. Numerically extracted contours of constant Lya-punov exponent near the origin, compared with those ofEq. (38). given by κ = 12 (cid:104) log( | s − j /s + j | ) (cid:105) = 0 , (40)a result that holds provided that s − j and s + j +1 are cho-sen from identical, independent probability distributionfunctions, as in Eq. (5). The vanishing value of κ atthe origin is numerically corroborated in Fig. 9. In fact,for λ = 0 (and for any probability distribution function P ( s )), there is an extended eigenfunction of Eq. (4) thatreads | λ = 0 , g = 0 (cid:105) = | (cid:105) + ( N − / (cid:88) m =1 ( − m s − s − ...s − m − s +2 s +4 ...s +2 m | m + 1 (cid:105)≡ | (cid:105) + ( N − / (cid:88) m =1 ψ m | m + 1 (cid:105) ,ψ m = ( − m s − s − ...s − m − s +2 s +4 ...s +2 m (41)where we have assumed N is odd and the amplitudes onall even sites vanish. A similar state can be constructedfor an even number of sites, with a mild restriction on s + N and s − N in both cases when periodic boundary condi-tions are imposed. That this zero energy state is indeedextended follows from the definition of the inverse lo-calization length within the transfer matrix method,[26-29] κ = lim N →∞ N (cid:104) log | ψ N |(cid:105) , where the average is overthe probability distributions of the matrix elements inEq. (41).4 − − − − FIG. 14. The eigenvalues of an exact numerical diagonaliza-tion of a matrix with N = 1000, g = 0 .
01 (blue dots). Fromeach point a line emanates, in the direction of the eigenvaluevelocity in the complex plane ( i.e . dλdg ), and with length pro-portional to the velocity. For eigenvalues away from the rim ofthe hole, no line is visible; spectral rigidity implies vanishingvelocities in this regime. F. Spectral rigidity outside the gap
Consider the model for g = 0, and “ramp up” g . Aswe argued in Sec. III D and as was discussed in Ref. [37]in the context of a related model, this results in a holethat tracks the contours of constant Lyapunov exponent.Thus, as g increases, the hole widens and “sweeps away”the eigenvalues in its vicinity. The hole hence acquires afinite fraction of the spectrum, concentrated on its one-dimensional rim. Since the rim of the hole correspondsto diverging Lyapunov exponent, these states have allbecome delocalized by the finite value of g , while thestates outside the hole are still localized, as explainedin Sec. II B. These states are insensitive to the boundaryconditions, and their eigenvalues will not be modified by g . This spectral rigidity is illustrated in Fig. 14. Tocalculate the eigenvalue velocity dλ/dg , we used first-order perturbation theory, which states that this deriva-tive is given by dλdg = dxdg + i dydg = ( (cid:126)v Lλ , B (cid:126)v Rλ ) , (42)where (cid:126)v Rλ > and (cid:126)v Lλ are the right and left eigenvectorsrespectively of the non-Hermitian matrix M ( g ) , ( .., ... )is the scalar product, and the matrix B is the matrixderivative of the matrix M ( g ) with respect to g . IV. PERTURBATION THEORY FOR LARGE g Our problem simplifies for large g . In this limitwe first neglect all terms of order e − g , and the remainingmatrix, with periodic boundary conditions, is of the form(illustrated for N = 4) M = s +1 e g s +2 e g
00 0 0 s +3 e g s +4 e g (43)with s + j = ± f = 1 / u = 1. We can attemptto “gauge out” the signs of { s + j } by applying a similaritytransformation H = Q − M Q with Q = c c c
00 0 0 c . (44)Choosing c c = s +1 , c c = s +3 ... c n c = s + n (with each c i = ±
1) results in a matrix of the form H = e g . (45)Note that H is proportional to the translational operatorfor a clockwise rotation of one lattice constant around thering. This procedure can only be applied when the prod-uct of the odd elements s +1 s +3 ... equals that of the evenelements s +2 s +4 ... (which occurs with probability 1 / { c i } in thiscase) leading to similar results.This matrix is readily diagonalized by plane waves, i.e ., right eigenvectors v Rj = e ikj , where the periodicboundary conditions imply that the allowed values of k must be k = 2 πn/N , n = 0 , ... ( N − v Lj = e − ikj . The resultingeigenvalues are then λ k = e g + ik , (46) i.e ., except for their magnitude e g they are the N rootsof unity. The eigenvectors of the original matrix M areplane waves modulated by random sign changes deter-mined by the elements s + j .So far we concluded that to zeroth order the eigen-values will sit at regular intervals on a circle . We maynow introduce the terms with the factor e − g as a pertur-bation, and calculate the shift of the eigenvalues to thefirst order in perturbation theory. The perturbation ma-trix is of the form (both before and after the similarity5 − − − − − − Numerical diagonalizationFirst-order perturbation theory FIG. 15. Comparison of first-order perturbation theory andexact numerical diagonalization, for g = 1, N = 10. The redcircle has radius e g . transformation) B = e − g s − s − s − s − . (47)Within first-order perturbation theory the shift in the k th eigenvalue is δλ k = ( (cid:126)v Lk , B (cid:126)v Rk ) , (48)and upon inserting the plane-wave eigenfunctions wehave δλ k = e − g [ e − ik ( s − + s − ... + s − n − + s − N ) /N ] , (49)Upon invoking the central-limit theorem, for large N we can replace the sum by a Gaussian variable withvariance N . Hence the eigenvalue will be shifted in the di-rection ± /λ , with a magnitude of order e − g / √ N . Notethat the magnitude of the shift is identical for all eigen-values. These results are illustrated in Fig. 15.It is straightforward to repeat these calculationsfor hopping elements s + j and s − j governed for the moregeneral probability distribution of Eq. (5). After a sim-ilarity transformation, M → M (cid:48) = P − M P , with P = diag { , /s +1 , / ( s +1 s +2 ) , ... } and up to corner ma-trix elements that do not affect our results as N → ∞ ,we have M (cid:48) = − N (cid:88) j =1 (cid:2) e g | j + 1 (cid:105)(cid:104) j | + s + j s − j e − g | j (cid:105)(cid:104) j + 1 | (cid:3) ≡ H + B , (50)where the periodic boundary conditions imply | j + N (cid:105) = | j (cid:105) . We recover the plane-wave eigenvectors discussedabove for H , and find from first-order perturbation the-ory (cid:104) λ k (cid:105) = e g + ik + e − g − ik N N (cid:88) j =1 s + j s − j . (51)With the help of the probability distribution of Eq. (5),we can now carry out a disorder average, with the result λ k = e g + ik + e − g − ik (1 + u ) ( f − / ≡ e g + ik + αe − g − ik , α = (1 + u ) ( f − / , (52)so that the eigenvalues will lie on an ellipse with majoraxis e g + αe − g and minor axis e g − αe − g . It is straight-forward to show for this generalized problem that thefluctuation of the k th eigenvalue about its mean valuesis O ( e − g / √ N ) h ( u, f ), where h ( u, f ) is a dimensionlessfunction of order unity.In Appendix C, we go to second -order in perturba-tion theory, and show that it leads to a similar picturequalitatively. V. SUMMARY AND OUTLOOK
In this paper we studied a coarse-grained, simplifiedmodel for the dynamics of neural networks, which uponlinearization close to a steady state leads to the studyof the eigenvector spectrum of an ensemble of sparse,non-Hermitian matrices. In contrast to most previousstudies in this context, here the connections were onlybetween neighboring neurons, i.e ., the model includeda spatial structure. For concreteness and simplicity, wefocused on a ring topology, which is realized in several in-stances in neuroscience . An additional parameter inour model, g , controlled the directional bias in the neu-ral network, i.e ., favoring clockwise over counterclockwiseconnections.Despite the deceptive simplicity of the model, itexhibits surprisingly rich behavior both in terms of theeigenvalue spectrum and in terms of the localizationproperties of the eigenvectors. Figure 16 shows the tra-jectories of eigenvalues for a particular instance N = 100,and for a value of g decreasing from one down to zero.The eigenvalues “flow” in the complex plane, until theirmotion ultimately ceases once the corresponding eigen-vectors become localized. For large values of g , we usedperturbation theory to show that the eigenvectors areapproximately plane waves (up to a similarity transfor-mation) and that the eigenvalues form a circle (or anellipse, more generally) in the complex plane. As g de-creases, eigenvalues move in the complex plane until theylocalize, after which “spectral rigidity” will take over andthe motion of the localized eigenvalue stops. The finalpositions of the eigenvalues for g = 0, when this game of“musical chairs” has ended, showed a remarkably intri-6 − − − − − − − − FIG. 16. Plot of the eigenvalues of a particular N = 100 ma-trix, where g is varied continuously from 1 (corresponding tothe outer circle) down to 0. These eigenvalues stop changingwith g when their eigenfunctions localize. (b) is an enlargedview of a part of (a). cate, fractal-like pattern . For any intermediate valueof g , the spectrum will show a pronounced “hole” or gapsurrounding the origin, with the eigenvalues which willultimately end inside the hole lying on its boundary, andwith localized states outside it.The spectra of conventional, highly-connected ran-dom matrices for large N can be grouped into univer-sality classes, such as those of the Gaussian orthogonalensemble and the Gaussian unitary ensemble, and thoseobeying the Circular Law . It is natural to ask aboutthe universality of the spectra and eigenfunctions of theone-dimensional sparse random matrices studied here.Because of its beautiful fractal-like spectrum, we havefocused here on directed localization in the bimodal non-Hermitian random hopping model of Feinberg and Zee .However, many of our conclusions also apply to the moregeneral model defined by Eqs. (4) and (5). For exam-ple, the symmetries under reflections across the real andimaginary axes and under 90 ◦ rotations in the complex plane discussed in Sec. II are preserved for arbitrary u when f = 1 /
2. As discussed in Sec. III E, there is alwaysa divergent localization length at the origin in this model.As summarized in Appendix B, when f = 1 /
2, approxi-mately equal numbers of the eigenvalues ( ∼ g = 0,as u varies from a bimodal distribution ( u = 1) to asymmetric double box distribution (0 < u <
1) to a sym-metric single box distribution ( u = 0). As f moves awayfrom 1 /
2, we expect that the spectrum becomes more el-liptical, consistent with the eigenvalue spectrum derivedin the large g limit in Eq. (52). Another aspect of uni-versality, that connected with Dale’s law in neuroscience,was addressed in Sec. II C.The gauge transformation argument leading toEq. (37) is quite general. Provided the localization lengthincreases monotonically at the origin, it predicts for 1drings a gap or hole bounded by a rim of extended states inthe spectrum for g > g c . Because the localization lengthdiverges at the origin for the model defined by Eqs. (4)and (5), we expect that g c = 0 in this case. When diag-onal disorder is present, the localization length for g = 0remains finite even at the origin and now g c > . Toillustrate the universal nature of the gap, Fig. 17 shows asingle box ( u = 0) spectrum with N = 1000 and g = 0 . d j chosen from a uni-form distribution with support [ − , N = 1000and g = 0 .
1. Although the single box spectrum inFig. 17(a) no longer has the fractal-like eigenvalue spec-trum shown in Fig. 2, a diamond-shape gap centered onthe origin with an enhanced density of states is clearlypresent. In Fig. 17(b) , we see that diagonal random-ness added to the bimodal model destroys the symmetryunder 90 ◦ rotations by removing the eigenvalue conden-sation onto the imaginary axis. Nevertheless, a hole inthe spectrum with an enhanced density of states on itsrim survives the imposition of diagonal randomness forthis value of g = 0 . > g c >
0. The large g perturbationtheory of Sec. IV can be used to show that all states aredelocalized (being plane-wave like) as g → ∞ for a wideclass of models, including those with diagonal random-ness. Hence, we expect that there exists another criticalvalue g c , such that for g > g c all states are delocalized.Localized eigenfunctions in neuroscience could be helpfulfor avoiding crosstalk between different neural computa-tion centers, and the extended states on the rim of thehole when g > areadapted to allow for spatial structure, with predator andprey species are localized to an array of lattice sites, butallowed to interact with their neighbors. For example,a site dominated by foxes would have an inhibitory ef-fect on neighboring sites occupied by rabbits, whereasrabbits would have an excitatory effect on nearby foxes.7t! − − − − − − − − − FIG. 17. (a) Single box ( u = 0) spectrum with N = 1000and g = 0 .
5. A macroscopic fraction of eigenvalues still con-densed onto the real and imaginary axes (compare Fig. 2), butthe remaining eigenvalues now have a diamond-like bound-ary. A diamond-shaped gap now replaces the octagonal holein Fig. 2(b). (b) Spectrum for the bimodal model with sym-metrical diagonal randomness with elements d j chosen froma uniformly from the interval [ − ,
1] , with N = 1000 and g = 0 .
1. Although only the real axis now has a macroscopicfraction of the eigenvalues, a lens-shaped gap replaces the oc-tagonal hole in Fig. 2(b).
Random excitatory and inhibitory connections in one di-mension could also be studied in chains of artificial cellswith spatially coupled gene expression patterns . VI. ACKNOWLEDGMENTS
We would like to thank F. Dyson, J. Hertz, Y. Lue,V. Murthy, T. Rogers and H. Sompolinsky for usefuldiscussions. Work by NH was supported by the JSPSKAKENHI Grants 15K05200, 15K05207, and 26400409.Work by DRN was supported by the NSF, through grantDMR13063667 and through the Harvard Materials Re-search Science and Engineering via grant DMR1420570.
Appendix A: Spectrum of the Hermitianrandom-hopping model
It is interesting to contrast our model of non-Hermitian localization with its Hermitian analogue,which also has a diverging localization length at the ori-gin and a connection between the density-of-states andthe inverse localization length. The Hermitian randomhopping model we consider is a reformulation of Eq. (1) H Herm = − N (cid:88) j =1 t j ( | j + 1 (cid:105)(cid:104) j | + | j (cid:105)(cid:104) j + 1 | ) , (A1)where { t j } is a set of mutually independent random vari-ables taking the values in the range [1 − ∆ , ≤ ∆ <
1. Although this is a standard one-dimensionalversion of the Anderson model , dominated by local-ized eigenstates, it is well established that the stateat λ = 0 is delocalized with both the localization lengthand the density-of-states diverging as | λ | → ρ ( λ ) andthe inverse localization length κ ( λ ) for ∆ = 0 .
85. We (a) − − (b) − − FIG. 18. (a) The density-of-states ρ ( λ ) and (b) the inverselocalization length κ ( λ ) for a Hermitian ring of length 1000with hopping randomness. We computed the former fromthe histogram of the eigenvalues that we obtained from ex-act diagonalization of 5000 random samples, while the latterfrom a new algorithm exploiting the Chebyshev-polynomialexpansion applied to 100 random samples. TABLE I. The fraction of the states on the real and imaginaryaxes as well as of the states with the zero eigenvalue. The datain the row “binomial” are for the binomial distribution ± − ,
1] of 500 samples of length 1000. distribution on on zeroreal axis imaginary axis eigenvaluesbinomial ( u = 1) 19.9% 19.9% 0 u = 0 .
95 19.8% 20.0% 0 u = 0 .
75 20.4% 20.6% 0 u = 0 . u = 0 .
25 24.7% 24.8% 66 (0.013%)one box (u=0) 33.7% 33.8% 792 (0.16%)numerically confirmed that the Hermitian version of theThouless formula connects these quantities : κ ( λ ) = P (cid:90) ∞−∞ dλ (cid:48) ρ ( λ (cid:48) ) log | λ − λ (cid:48) | , (A2)where P denotes the principal part. We can indeed seeevidence for singularities around λ = 0 in both figures.These singularities are expected to take the forms ρ ( λ ) ∼ | λ | log (1 / | λ | ) ; (A3) κ ( λ ) ∼ / log(1 / | λ | ) . (A4) Appendix B: Density-of-states on the real andimaginary axes for f = 1 / and arbitrary u In this Appendix we study the density of eigenstatesthat have condensed on the real and imaginary axes forthe model defined by Eqs. (4) and (5) with f = 1 / g = 0 as a function of u . For numerical purposes, wehere define the states to lie on the real and imaginaryaxes provided | Im λ n | < − , (B1) | Re λ n | < − . (B2)We list the fraction of the eigenvalues that satisfy theseconditions in Table I, where the zero eigenvalues are thestates that satisfy both criteria. As discussed in Sec. III,these states play an important role in determining how κ ( x, y ) vanishes near the origin.In all cases, the density-of-states (see Fig. 19) is sta-tistically the same on the real and imaginary axes. Thezero eigenvalues are absent until the two-box distribution(Eq. (5), 0 < u <
1) becomes close to the one-box dis-tribution ( u (cid:46) u ≤ . λ = 0. Appendix C: Perturbation theory of thesign-random tight-binding chain
We summarize here second-order perturbation the-ory applied to our model with f = 1 / u = 1, for large g .Upon adopting the similarity transformation of Eq. (9), T jj = j − (cid:89) k =1 b k (C1)we can bring the tridiagonal hopping matrix M ( g ) = s +1 e g s − N e − g s − e − g s +2 e g s − e − g . . . . . . s + N − e g s + N e g s − N − e − g , (C2)into the form M (cid:48) ( g ) = T − M ( g ) T = e g r N e − g r e − g e g r e − g . . . . . . e g e g r N − e − g , (C3)where all remaining matrix elements in Eqs. (C2)and (C3) are zero, r j = s + j s − j , j = 1 , .., N, (C4)and we have assumed (cid:81) j =1 ..N s + j = 1 in order to getsimplified corner matrix elements. The elements r j arepositive or negative random numbers if s + j and s − j arerandom and both positive and negative; In particular,when s + j and s − j are ±
1, we have r j = ± M ( g ) and M (cid:48) ( g ) are then isospectral.We now split the matrix M (cid:48) ( g ) into the matrix withelements proportional to e g and the matrix with the ele-ments proportional to e − g , and formulate the perturba-tion of the spectrum of the former with respect to the9 − − − − (a) − − − − (b) − − − − (d) − − − − (c) − − − − (f) − − − − (e) FIG. 19. The density-of-states for f = 1 / g = 0 on the real axis (filled light blue bars) and on the imaginary axis (open redbars) plotted as histograms N ( λ ): (a) binomial distribution ( u = 1); (b) u = 0 .
95; (c) u = 0 .
75; (d) u = 0 .
5; (e) u = 0 .
25; (f)one-box distribution ( u = 0). The system size is 1000 for all data and the number of samples is 500 for all data except for (a),where it is 1000. latter. We thus set M (cid:48) ( g ) = M + M , where M = e g , (C5) M = e − g r N r r
0. . . . . . r N − r N − . (C6) The zeroth-order eigenvalues and eigenvectors of M are given by λ (0) k n = e g + ik n , (C7) (cid:104) x | ψ (0) k n (cid:105) = 1 √ N e ik n x , (C8)where k n := 2 πnN , n = 0 , , , · · · , N − . (C9)0By setting ξ = Re λ (0) k n = e g cos k, (C10) η = Im λ (0) k n = e g sin k, (C11)we see that the eigenvalues are equidistantly aligned ona circle of radius e g in the complex λ plane: ξ + η = e g . (C12)Similarly to Sec. IV, we find the first-order eigenvalueperturbative shift in eigenvalues, λ (1) k n = (cid:104) ψ (0) k n | A | ψ (0) k n (cid:105) = e − g − ik n ¯ r (0) , (C13)where ¯ r (0) is the component at k = 0 of the Fouriertransform of the random variable r x :¯ r ( k ) = 1 N N − (cid:88) x =0 r x +1 e ikx . (C14)We can cast it in the following way: ξ := Re λ (1) k n = e − g ¯ r (0) cos k n , (C15) η := Im λ (1) k n = − e − g ¯ r (0) sin k n , (C16)Note that the first-order perturbation does not depend on the details of the random numbers but only on theaverage. Since we use the random numbers with a sym-metric probability distribution, ¯ r (0) vanishes in the limit N → ∞ . For a finite value of N ( N = 16), we find themovement of the eigenvalues as illustrated in Fig. 20(a).The second -order eigenvalue perturbation, obtainedby similar techniques, is given by λ (2) k n := − N − (cid:88) m =0 m (cid:54) = n (cid:104) ψ (0) k n | A | ψ (0) k m (cid:105)(cid:104) ψ (0) k m | A | ψ (0) k n (cid:105) λ (0) k m − λ (0) k n (C17)= − e − g N − (cid:88) m =0 m (cid:54) = n e − i ( k m + k n ) e ik m − e ik n | ¯ r ( k n − k m ) | , (C18)which we illustrate in Fig. 20(b). The third-order eigen-value perturbation is illustrated in Fig. 20(c) too. Moreanalyses reveal that the k th-order corrections generallybehave as λ (2) k ∝ e − g − ik , λ (3) k ∝ e − g − ik . (C19)Although large g perturbation theory is useful for captur-ing eigenvalue trends, it does not seem capable of deter-mining when eigenvalues stop moving with increasing g ;the corresponding eigenvalues remain delocalized withinthis approach. M. L. Mehta,
Random matrices, 3rd edition , Vol. 142 (Aca-demic press, New York, 2004). T. A. Brody, J. Flores, J. B. French, P. Mello, A. Pandey,and S. S. Wong, Reviews of Modern Physics , 385 (1981). H. Sommers, A. Crisanti, H. Sompolinsky, and Y. Stein,Physical Review Letters , 1895 (1988). E. P. Wigner, Ann. of Math., 62, 548 (1955); Ann. of Math.67, 324 (1958). J. Ginibre, Journal of Mathematical Physics , 440 (1965). V. L. Girko, Teoriya Veroyatnostei I Ee Primeneniya ,669 (1984). K. Rajan and L. Abbott, Physical Review Letters ,188104 (2006). P. W. Anderson, Physical Review , 1492 (1958). B. Shklovskii and A. Efros,
Electronic properties of dopedsemiconductors (Springer-Verlag, Berlin, 1984). R. Chaudhuri, A. Bernacchia, and X.-J. Wang, Elife ,e01239 (2014). F. L. Metz, I. Neri, and D. Boll´e, Physical Review E ,055101 (2011). I. Neri and F. Metz, Physical Review Letters , 030602(2012). H. Rouault, and S. Druckmann, Spectrum density of largesparse random matrices associated to neural networks,arXiv:1509.01893 (2015). See, e.g., T. A. L. Ziman, Physical Review Letters 49, 337(1982), and references therein. D. Thouless, Journal of Physics C: Solid State Physics , 77 (1972). P. Dayan and L. F. Abbott,
Theoretical neuroscience (Cambridge, MA: MIT Press, 2001). J. Hertz, A. Krogh, and R. G. Palmer,
Introduction to thetheory of neural computation (Basic Books, 1991). O. Shriki, D. Hansel, and H. Sompolinsky, Neural compu-tation , 1809 (2003). T. P. Vogels, K. Rajan, and L. Abbott, Annu. Rev. Neu-rosci. , 357 (2005). N. Hatano and D. R. Nelson, Physical Review B , 8651(1997). I. Y. Goldsheid and B. A. Khoruzhenko, Physical ReviewLetters , 2897 (1998). D. R. Nelson and N. M. Shnerb, Physical Review E ,1383 (1998). N. M. Shnerb and D. R. Nelson, Physical Review Letters80, 5172 (1998); See also K. A. Dahmen, D. R. Nelson, andN. M. Shnerb. ”Population dynamics and non-hermitianlocalization”, in Statistical Mechanics of Biocomplexity(Springer, Berlin,1999) pp. 124-151. K. B. Efetov, Physical Review Letters , 491 (1997). P. Brouwer, P. Silvestrov, and C. Beenakker, PhysicalReview B , R4333 (1997). J. Feinberg and A. Zee, Physical Review E , 6433 (1999). D. E. Holz, H. Orland, and A. Zee, Journal of Physics A:Mathematical and General , 3385 (2003). S. Chandler-Wilde, R. Chonchaiya, and M. Lindner, Op-erators and Matrices , 739 (2013). (a) − − − − − − (b) − − − − − − (c) − − − − − − FIG. 20. The directions of eigenvalue perturbations with N = 16 and g = 1 for a random sample out of a binomialdistribution. The blue dots indicate the positions of zeroth-order eigenvalues (C8) in all panels. The arrows in each panelindicate: (a) the first-order perturbations (C13) magnified 50times; (b) the second-order perturbations (C18) magnified 20times; (c) The arrows indicate the third-order perturbationsmagnified 100 times. S. Chandler-Wilde and E. Davies, Journal of Spectral The-ory , 147 (2012). S. Chandler-Wilde, R. Chonchaiya, and M. Lindner, Op-erators and Matrices , 633 (2011). R. Hagger, On the Spectrum and Numerical Range ofTridiagonal Random Operators, to appear in Journal ofSpectral Theory, arXiv:1407.5486 (2014). R. Hagger, Symmetries of the Feinberg-Zee Random Hop-ping Matrix, to appear in Random Matrices: Theory andApplications, arXiv:1412.1937 (2014). B. Derrida, J. L. Jacobsen, and R. Zeitak, Journal of Sta-tistical Physics , 31 (2000). H. Furstenberg, Transactions of the American Mathemat-ical Society , 377 (1963). K. Ishii, Progress of Theoretical Physics Supplement ,77 (1973). A. Crisanti, G. Paladin, and A. Vulpiani,
Products of ran-dom matrices (Springer, 1993). L. Molinari and G. Lacagnina, Journal of Physics A: Math-ematical and Theoretical , 395204 (2009). N. Hatano and D. R. Nelson, Physical Review Letters ,570 (1996). B. Doiron and A. Litwin-Kumar, Frontiers in computa-tional neuroscience (2014). J. D. Seelig and V. Jayaraman, Nature , 186 (2015). R. M. May, Mathematical Biosciences , 59 (1971). R. M. May, Nature , 413 (1972). M. W. McCoy, B. M. Bolker, K. M. Warkentin, and J. R.Vonesh, The American Naturalist , 752 (2011). E. Karzbrun, A. M. Tayar, V. Noireaux, and R. H. Bar-Ziv, Science , 829 (2014). G. Theodorou and M. Cohen, Physical Review B , 4597(1976). T. Eggarter and R. Riedinger, Physical Review B , 569(1978). N. Hatano and J. Feinberg, in preparation. R. Silver and H. R¨oder, Int. J. Mod. Phys. C , 735 (1994). A. V. R.N. Silver, H. Roeder and J. Kress, J. Comp. Phys. , 115 (1996). R. Silver and H. R¨oder, Physical Review E56