{'fraction_non_alphanumeric': 0.030689329556185082, 'fraction_numerical': 0.011769863752866585, 'mean_word_length': 4.353493410362882, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Monte Carlo Determination of Multiple Extremal Eigenpairs

We present a Monte Carlo algorithm that allows the simultaneous determination of a few extremal eigenpairs of a very large matrix without the need to compute the inner product of two vectors or store all the components of any one vector. The new algorithm, a Monte Carlo implementation of a deterministic one we recently benchmarked, is an extension of the power method. In the implementation presented, we used a basic Monte Carlo splitting and termination method called the comb, incorporated the weight cancellation method of Arnow {\it et al.}, and exploited a new sampling method, the sewing method, that does a large state space sampling as a succession of small state space samplings. We illustrate the effectiveness of the algorithm by its determination of the two largest eigenvalues of the transfer matrices for variously-sized two-dimensional, zero field Ising models. While very likely useful for other transfer matrix problems, the algorithm is however quite general and should find application to a larger variety of problems requiring a few dominant eigenvalues of a matrix.


I. INTRODUCTION
A common problem in computational physics is computing the eigenpairs of large matrices. We will present a new Monte Carlo algorithm that allows the simultaneous determination of a few extremal eigenpairs of a very large matrix without the need to orthogonalize pairs of vectors to each other or store all the components of any vector. This algorithm is a Monte Carlo implementation of deterministic one we recently benchmarked [1]. The newly benchmarked algorithm is based on a refinement of the power method recently developed by Booth [2,3] and does not require, as does the standard Ritz estimator [4,5], the ex- Carlo splitting and termination technique called the comb [6,7], incorporates the weight cancellation method of Arnow et al. [8], and exploits a new sampling method, the sewing method [9], that does a large state space sampling as a succession of small state space samplings.
In refining the power method [1], we were targeting its use on matrices so large that they are unassailable deterministically because no single vector can be stored in memory. As the system size increases, finding a few extremal eigenpairs of the transfer matrix of the twodimensional Ising model becomes such a problem. We will illustrate the effectiveness of the algorithm by determining the two largest eigenvalues of the transfer matrices for variouslysized two-dimensional zero field Ising models, exploiting Onsager's exact results [10,11] for their values as benchmarks. We comment that two extremal eigenvalues of this matrix are of significant physical interest: the logarithm of λ 1 is proportional to the free energy, and the logarithm of the ratio λ 2 /λ 1 is proportional to the reciprocal of the correlation length that controls long range spin correlations near the critical point. Although our algorithm is extendable to finding more than two extremal eigenpairs, we will focus on finding just λ 1 and λ 2 [10,11].
In the next section, Section II, we summarize the basic features of the transfer matrix of the two-dimensional Ising model in a zero magnetic field. In subsequent sections we will reference these features to make our algorithm presentation more concrete. In Section III, we summarize our extension of the power method for the determination of the two eigenpairs corresponding to the two eigenvalues of largest absolute value. Then, in Section IV, we discuss the basics of our Monte Carlo implementation of this algorithm. We first use the algorithm for medium-sized matrices for which computer memory is adequate to store the eigenvectors and then use it for much larger-sized matrices for which it is not. For the Monte Carlo sampling of states in latter case, we use the sewing algorithm [9] that facilitates sampling of states from a large space from a smaller space. Results presented are for the determination of the two largest eigenvalues of the Ising transfer matrix for various lattice sizes. In the last section, we summarize our work and comment on its likely application to other systems.
We note that the intent of the present application is presenting and benchmarking a new and relatively general numerical method and not numerically studying the finite-size scaling of the eigenvalues of the matrix used for the benchmarking. Because we reproduce the exact eigenvalues to satisfactory accuracy, our estimates will enjoy the same scaling as the well known exact results for this problem. When the magnetic field is not zero and the eigenvalues are not exactly known, various researchers have calculated up to the four largest eigenvalues by a deterministic version of the power method and have made extensive studies of the scaling of these eigenvalues [12]. With the algorithm to be presented we have reproduced the first two of these eigenvalues and hence would reproduce the basic scaling. A more extensive study of the eigenvalues of the field dependent Ising model will be presented elsewhere [13]. The Monte Carlo approach increases the lattice sizes accessible. We note that these deterministic calculations required multiple processors. All our calculations were done on a single processor.

II. ISING MODEL
In his much celebrated work Onsager calculated many of the properties of the twodimensional Ising model exactly [10]. His calculations started with the expression of the partition function in terms of its transfer matrix [10,14,15]. He then found all the eigenvalues of this matrix analytically and showed that the scaling the dominant eigenvalue with the area in the thermodynamic limit (the area of the model approaching infinity) implied the onset of long-range ordering among the spin variables. The ratio of the second largest eigenvalue to the first is associated with the spatial behavior of this ordering. At the critical temperature, this ratio approaches unity as the area approaches infinity.
We will consider an m × n Ising model defined with periodic boundary conditions in one direction and open boundary conditions in the other. The two-dimensional Ising model's energy is Here, (i, j) are the coordinates of a lattice site. The Ising spin variable µ i,j on each site has the value of ±1, the exchange constant J > 0, and µ i,m+1 = µ i,1 . The symbol σ j = (µ 1,j , µ 2,j , . . . , µ m,j ) denotes a column configuration of Ising spins and there are 2 m possible configurations for each column.
The transfer matrix A(σ, σ ′ ) follows from a re-expression of the partition function [11] Z (m, where ν = J/k B T , k B is Boltzmann's constant, and T is the temperature, and A(σ, σ ′ ) is a 2 m × 2 m matrix whose elements are As customary for Ising model simulations, we represent a configuration σ by the first m bits of an integer between 0 and 2 m − 1, with a set bit being a +1 Ising spin and an unset bit being a −1 spin. This convention maps the matrix element A(σ, σ ′ ) between configurations to an element A ij between integers.
We note that A(σ, σ ′ ) is asymmetric and its elements are greater than zero so the matrix is maximally dense and hence irreducible. Because of the positivity and irreducibility the Perron-Frobenius Theorem [16] says that dominant eigenvalue is real, positive, and nondegenerate and all components of the corresponding eigenstate are real and have the same sign. The two largest eigenvalues of A, for finite m, are [17] λ 1 = (2 sinh 2ν) m/2 exp 1 2 (η 1 + η 3 + · · · + η 2m−1 ) where cosh η k = cosh 2ν coth 2ν − cos πk m The transfer matrix of the Ising model can be symmetrized [11] but we saw no computational advantage for using this form.

III. POWER METHOD
For some real-valued M × M matrix A, not necessarily symmetric, we will be concerned with the M eigenpairs (λ α , ψ α ) satisfying In the simplest application of the power method [4], an iteration is started with some ψ, normalized in a convenient, but otherwise relatively arbitrary, manner and consists of iterating the two steps φ = Aψ until some convergence criterion is met. If we write and if |λ 1 | > |λ 2 | ≥ |λ 3 | ≥ · · · ≥ |λ N |, then after n iterations Accordingly, as n → ∞, Thus, the dominant eigenpair is simultaneously determined. For the norm of the vector φ whose components are φ i , a frequent choice is For deterministic calculations of a few dominant eigenpairs, say N, one of two approaches are typically tried. One approach is to use the power method to determine the dominant eigenpair, use deflation to project out this state out of the matrix, and then reuse the power method on the deflated matrix. To determine several eigenpairs simultaneously, the power method can be generalized to where Φ and Ψ are M × N matrices whose columns are orthogonalized to each other.
This orthogonality needs maintenance throughout the computation or else all N vectors, represented by the columns of the initial Ψ, will converge to the one associated with the dominant eigenvalue [5]. This approach is called orthogonal [5] or simultaneous iteration [18].
For Monte Carlo calculations of a few dominant eigenpairs, we are proposing a quite different approach based on Booth's proposed refinement of the power method [2,3]. This refinement uses the observation that for any eigenpair (λ, ψ) and for each non-zero component of the eigenvector, the eigenvalue equation Aψ = λψ can be rewritten as and that similar equations can also be written for any number of groupings of components, where the R i are rules for different groupings. The groupings are quite flexible: they can overlap and their union need not cover the entire space. In addition, any two groupings, say From L groupings of the components, Booth constructed L estimators for the L largest eigenvalues and forced them to become equal by adjusting certain parameters at each iteration step. For the dominant two eigenvalues, we will do something similar.
First we note that for almost any starting point ψ = α ω α ψ α , the power method will converge to (λ 1 , ψ 1 ). To find the two extremal eigenvalues, we need two normalized, but not necessarily orthogonal, starting points . At each step, we will apply A to them individually. Without any intervention both will project the same dominant eigenfunction so at each step we adjust the relationship between their sum to direct one to the dominant state and the other to the next dominant one.
Formally, we start the iteration with ψ = ψ ′ + ηψ ′′ . Suppose at the n th step, ψ ′ and ψ ′′ have iterated toψ ′ andψ ′′ , then at the (n + 1) th step we invoke (15) and and (16) and require that which leads to The algorithm is to apply A repeatedly until two real solutions η 1 and η 2 for (18) exist. One solution will be then used to guide further iterations to (λ 1 , ψ 1 ); the other, to (λ 2 , ψ 2 ).
The power method becomes [1]: choose two starting points ψ ′ and ψ ′′ , which need not be orthogonal, then for each iteration step, compute and if the roots of (18) are real [19], update using Eigenvalues can be estimated from where η 1 and η 2 generate the largest and next largest eigenvalue estimates.

IV. MONTE CARLO IMPLEMENTATION
The basic operation of a power method is a matrix-vector multiplication. Here, we now describe how we used the Monte Carlo method to estimate the repetition of such multiplications.
In the basis defining the matrix elements of A, we write We will call the amplitudes ω ′ i and ω ′′ i weights even though they are not necessarily all positive nor are the sums of their absolute values unity. We will also assume that the elements A ij of the M × M matrix A are easily generated, as we are ultimately interested in cases where M is so large that this matrix must be generated on-the-fly as opposed to being stored. Next, we imagine we have N particles distributed over the M basis states defining A. Generally, N ≪ M. Then, at each iteration step, we interpret A ij as the weight of a particle arriving in state |i per unit weight of a particle in state |j . The action of A ij on a ψ thus causes all particles currently in state |j to jump to |i , carrying to |i their current weight ω j , modified by A ij .
The jumps will be executed probabilistically. To do this we let the total weight leaving state |j be and define the transition probability from |j to |i be The number W j is called the state weight multiplier. How we use these densities will depend on the size and types of the matrices under consideration.

A. Medium Matrices
If M is sufficiently small so we can store all components of our vectors, then a Monte Carlo procedure for jumping is easily constructed. Instead of always (i.e., with probability 1) moving weight A ij from state |j to state |i , we will instead sample a |i from T ij and multiply the transferred weight by the ratio of the true probability (1.0) to the sampled probability (T ij ); that is, if state |i is sampled, the weight arriving in state |i from |j is multiplied by As for many Monte Carlo simulations, as is the case for the transfer matrix of the Ising model, the particle weights defining the eigenvector associated with the largest eigenvalue can be made all positive. The second eigenfunction however must be represented by some particles of negative weight and some particles of positive weight. These negative and positive weights must for some jumps at least partially cancel to maintain a correct estimation of the second eigenfunction. When N ≪ M, this cancellation does not occur often enough in a Monte Carlo simulation without proper design: because the number of states vastly out numbers the number of particles, the probability that a negatively and a positively weighted particle randomly arrive in the same state becomes trivially small.
There are several ways to arrange the cancellation [2]. We found the Arnow et al. [8] algorithm effective and convenient. First, we consider two particles of weights w 1 and w 2 in states |j 1 and |j 2 , and then let T ij 1 and T ij 2 be the probabilities of their reaching state |i .
Next we use the weight multiplier W j for jumping from state |j to state |i and suppose that states |i 1 and |i 2 are sampled from the density T ij 1 +T ij 2 . This density can be sampled by sampling a |i 1 from T ij 1 and a |i 2 from T ij 2 . The true probability that particle 1 jumps to |i 1 is T i 1 j 1 and the true probability that particle 2 jumps to |i 1 is T i 1 j 2 . The ratio of the true density to the sampled density for particle 1 is so that the weight arriving at state |i 1 from particle 1 is The true probability that particle 2 arrives at |i 1 is T i 1 j 2 . The ratio of the true density to the sampled density for particle 2 is so that the weight arriving at |i 1 from particle 2 is Thus the total weight arriving at |i 1 is By similar arguments, the total weight arriving at |i 2 is We note that to get meaningful cancellation, say between particle 1 with weight w 1 < 0 and particle 2 with weight w 2 > 0, the transition probabilities must overlap somewhat. For example, if in (33) T i 2 j 1 = ǫ ≪ 1, then the total weight arriving at i 2 is essentially just the weight arriving of particle 2 alone.
How one arranges better overlap is problem dependent. For our Ising simulations, we sorted the particles into state order (a state is represented by the bits of an integer). Particles 1 and 2 were then sampled together according to the Arnow et al. scheme, then particles 3 and 4 are sampled together, and so forth. The fact that the list is ordered means that there are (typically) many nearby states |i that are accessible from both particles ℓ and ℓ + 1 with nontrivial transition probabilities.
As the iteration progresses, the absolute value of the weights of some particles becomes very large, and those of some others, very small. As standard for Monte Carlo methods with weighted particles, particles with weights of small magnitude are stochastically eliminated and those with large magnitudes are stochastically split. To do this we used a procedure called the comb [6,7]. It is described in Appendix A.
The steps of the algorithm are: First, we initialize the states and weights of two vectors.
For the Ising simulation, each vector had the same states but different weights. We selected the states uniformly and randomly over the interval (0, 2 m −1) and selected the ω ′ i uniformly and randomly over the interval (0,1) and the ω ′′ i uniformly and randomly over the interval (-0.5,0.5). Then, for a fixed number of times we iterate. At each iteration we execute the jump procedure for each particle, place the particle list in state order, effect cancellations, estimate the eigenvalues from (23), update ψ ′ and ψ ′′ , and then comb. R 1 consisted of the states for which more than half of its m bits were 0's and R 2 consisted of the states for which more than half of its m bits were 1's. Additional algorithmic details are given in Appendix B. Table I shows the computed λ 1 and λ 2 for m = 12. A computer program was written to make 20 independent (different random number seeds) calculations using N particles per iteration, for various values of N. Note that even when N = 100 << M = 2 12 = 4096 the method still separates the eigenfunctions. The simulations were run at the critical temperature of the infinite lattice; that is, ν = 0.4406867935097715 [11]. The last line in the table gives the eigenvalues deterministically obtained by using (5). From the last few lines, we also see that when the number of particles floods the number of states, that is, when all the basis states are being used multiple times, exceptional accuracy is obtained.
Also in the table is the timing for each run. The runs were done on a single 1.5 GHz Mac PowerPC processor. We note that the runs appear to scale sub-linearly with the number of particles. This scaling is deceiving. For a small number of particles, the time required to set up the transition matrix A ij is significant compared to the rest of the calculation. The Monte Carlo part is dominated by the state order sort which scales as N log N, so as N becomes very large, the run time scaling should eventually be slightly super-linear.

B. Large Matrices
For M ≤ 2 12 , sampling from the cumulative probability C i = i k=0 T kj works well. If the number of states gets too large (M = 2 12 was our limit), then C i cannot be sampled directly because it cannot fit in the computer's memory. In this case, we could just randomly pick from any state |j any state |i with probability 1/M instead of always picking a state |i (i.e., picking it with probability 1). Then, if state |i is sampled, the weight of particles arriving in state |i is The problem with this approach is that the A ij can have immense variation so that this simple sampling scheme is unlikely to work well as a Monte Carlo method. This situation is especially true for the Ising problem. A large part of such variations however can be removed by sampling the new state in stages and then sewing the stages together.
To explain the sewing procedure [9], we will first assume that we can write any state |i in our basis as a direct product of the states in a smaller basis, for example,

Instead of transferring weight
from state |j to state |i with probability we will use the a ij that would apply to the smaller set of states and then make an appropriate weight correction.
For each smaller set of states, we rewrite the analogous transition probability from state |j to state |i as and the analogous weight multiplier as We thus will sample |i 1 and |i 2 from the probability function Now, we define C ij to be the weight correction necessary to preserve the expected weight transfer from state |j = |j 2 |j 1 to state |i . It satisfies Using (38) in (41), we find Thus This sewing method generalizes easily. For k sets of states, (42) and (43) become and For the Ising problem, the weight correction is where νD i is the energy difference per k B T between calculating with the bits together and the bits separately [9]. Using this sewing algorithm for the sampling of states, we computed the first and second eigenvalues and their standard deviations for variously sized two-dimensional Ising models by sewing sets of 8 bits. The results are shown in Table II processor. ν = 0.4406867935097715, the bulk critical value. R 1 consisted of the states for which more than half of its m bits were 0's and R 2 consisted of the states for which more than half of its m bits were 1's. In Table III, we present the run averages for the M = 2 48 case. We note that 2 48 ≈ 2.8 × 10 14 .

V. CONCLUDING REMARKS
We have presented a Monte Carlo algorithm that enables the determination of two extremal eigenvalues of a very large matrix. The explicit demonstration of the power of our algorithm was the determination of the two largest eigenvalues of the transfer matrix of the two-dimensional Ising model in a zero magnetic field. The convenience of this matrix was the existence of exact expressions for its eigenvalues for any finite-sized system. We were able to reproduce the exact values for the two dominant eigenvalues to within the statistical error of our simulations.
Physics and chemistry provide numerous problems where obtaining several extremal eigenvalues is important. One problem class is finding the ground state and a few excited states of atoms, molecules, solids, and nuclei [20]. For these quantum problems, quantum Monte Carlo Monte Carlo (QMC) projector methods [21,22] are commonly used. In fact, it was for these methods that the Arnow et al . [8] method used here was first proposed as a method for taming the fermion sign problem [21]. We comment that our new algorithm requires sampling from the second eigenfuction which must have positive and negative components just as a fermion state must also have. Our successful cancellation of signed particles representing this eigenfunction and the slightly superlinear scaling of our computation time with system complexity, as opposed to an exponential scaling, demostrates that we have "solved" the sign problem [23] for our application.
We also note that several investigators [24] have adapted the technology of the diffusion QMC method [21]  a priori and was used to deflate the matrix so the projection was to the second largest eigenvalue. We comment that Booth [2] has presented an algorithm for determining a single eigenvalue without the need to know or concurrently determine any other one.
Our transfer matrix was positive, asymmetric, and dense. How is our algorithm changed if the matrix lacks one or more of these properties? For simple test cases, we have successfully constructed deterministic procedures for matrices whose elements are complex valued. Also for simple test cases, we have had success for real asymmetric matrices whose eigenvalues are complex valued. Devising Monte Carlo algorithms for real, symmetric, sparse matrices has however received more of our attention [28].
For an indefinite but symmetric matrix, the symmetry and reality assure the reality of the eigenpairs, but do not fix the sign of the dominant eigenvalue or the components of its eigenvector. Such matrices can also have the dominant eigenvalue being degenerate. In our algorithm a signature of the dominant eigenvalue being degenerate is the quantity q 1 (19) in our quadratic equation (18) becoming extremely small. In cases where we knew the degeneracy a priori , we saw our estimate of λ 1 being very close to λ 2 , that is, we saw our estimates of the dominant degenerate eigenvalue being slightly but artificially separated.
This type of difficulty is inherent to a power method. If the matrix elements are of mixed sign, the modifications to the algorithm are quite simple: we take the absolute value of the A ij and multiply the weight multiplier by the sign of A ij . We also have to account for the sign of the weights ω ′ i and ω ′′ i appropriately. Weight cancellation must be done with care. If the matrix is sparse, other options for selecting the state to jump to become available that avoid storing the cumulative distribution, improve efficiency, and reduce variance [28].
Here, we have focused principally on the determination of eigenvalues, noting the rapid convergence to their values, but what about the determination of the eigenvectors? If we can store the vectors in computer memory we can accurately determine them. For a problem in a continuum, for example, finding the dominant eigenvalues of an integral equation such as those occurring for the transport of neutrons, we have accurately estimated the two dominant eigenfunctions with an efficiency enhanced over the basic power method for determining just the dominant eigenfunction [29].
The dynamics of determining eigenvectors is as follows: Because the standard power method for the dominant eigenpair converges to ψ 1 as λ 2 /λ 1 (10), the dominant eigenvector becomes increasingly difficult to determine as λ 1 and λ 2 become very close in magnitude.
Because we are projecting to the first two dominant eigenpairs, our method converges to ψ 1 as λ 3 /λ 1 . This is usually a gain relative to the standard power method. Convergence to ψ 2 goes as λ 3 /λ 2 . This ratio illustrates the fact that in principle we can accelerate convergence to ψ i by seeking the L > i highest eigenfunctions and thus converging to ψ i as λ L+1 /λ i .
To do this, we need to use L groupings and L starting states and adjust the sum of the iterated states to guide them to convergence to L different eigenvalues. We have observed accelerated convergence for the L = 3 case on simple test cases. We comment that finding the eigenvectors of the transfer matrix of the Ising and other models is a problem for which finding the eigenvectors becomes increasingly difficult as the system size increases, if we could in fact store them. For this type of matrix, the eigenvalue spectra is expected to become quasi-continuous starting with the second eigenvalue and the first and second eigenvalues become asymptotically degenerate at a critical point [30].
In closing, we believe that the algorithm presented here is accurate, easy to implement, and applicable to many other problems. Different problems afford the opportunity to improve its efficiency by modifying some of its details. Wider use of the algorithm will define more crisply its strengths and limitations than is possible by just the present application.

APPENDIX A
The comb is a stochastic procedure for selecting M particles, not necessarily all different, from a list of N unevenly weighted particles and preserving the total weight. If ω i > 0 is the weight of an individual particle in the original list, its weight in the new list becomes i ω i is the total weight. To effect this procedure, we construct the cumulative sums of the original weights where j = 0, 1, 2, . . . , N and C 0 = 0 and C N = W T , draw a random number ξ uniformly distributed over (0,1), and then for k = 1, 2, . . . , M, select particle j for the new list if Dependent of the size of the difference C j−1 − C j relative to W T /M, particle j is selected zero, one, or more times. The procedure calls the random number generator only once, and the new particle list is always of a predetermined size.
In the present simulations we have a list of particles and to each particle i is associated a state |i and two weights, ω ′ i and ω ′′ i . If we were to comb the two lists of weights separately, we would in general produce two lists of states. We wanted one such list. Accordingly, we combed in the following manner: we formed and then generated the cumulative distribution.
Now we combed the C j as before, noting that we have chosen a normalization such that W T = 1. Instead of giving the selected particle j a weight W T /M for both w ′ j and w ′′ j , we instead assign The following are some useful particulars of our Monte Carlo implementation of our modified power method. First, instead of updating via (21) we update via This form is more symmetric and makes more explicit that we are trying to remove ψ ′′ from ψ ′ and vice versa. To avoid overcorrecting we have sometimes found it useful to update via where η ′ 1 = sign(η 1 ) min(α, |η 1 |) (B3) η ′ 2 = sign(η 2 ) min(α, |η 2 |) (B4) and α is some small positive number.
Besides the eigenvalue estimators (23), there are others (15) equally valid that provide multiple cross-checks to (23). These include all of which involve sums needed for (23) so they cost nothing extra to compute. Another use of these estimators is monitoring the utility of the regions R 1 and R 2 , as an estimator will perform inconsistently if too few particles occupy a region or if the sum of the weights in a given region is small. The most useful regions are those that are a major subset of the positive and negative regions of the second eigenvector.
For the Ising model we obtained some additional useful weight cancellation by exploiting the fact that most of the particles are in fully magnetized states and before using the comb replacing all the particles in the "up" state by a single particle whose weight is the sum of the "up" weights of these particles and doing similarly for the particles in the "down" state.
Because the ω ′ i should be positive, after combing we explicitly set each ω ′ i equal to its absolute value. THe step is useful because even though we initialize ψ ′ to have all positive components, the weight cancelation procedure can make some of them negative. Once so, they might project out rapidly and may adversely affect the eigenpair estimation until they do.