Appendix A Alias method#

The methods to sample from a discrete probability for n events, presented in Algorithms 1 and 2, become inefficient when n is large. To devise an efficient algorithm for large n, we must distinguish the case where the list of probabilities is constantly changing from the case where it remains unchanged. In the latter case, there are several ways to boost efficiency by performing a modestly expensive operation once and then using more efficient operations in subsequent samplings. One such approach is to sort the list of probabilities, an operation generally requiring (n log n) operations, and use a bisection method, requiring (log n) operations, to select the events from the list. What we really want, however, is a method requiring only (1) operations. Walker’s alias algorithm (Walker, 1977) has this property.

Suppose we need to repeatedly choose one of five colors randomly – red, blue, yellow, pink, and green – with weights 6, 1, 3, 2, and 8 (Fig. A.1). We first generate five sticks with lengths 6, 1, 3, 2, and 8 and paint each stick with the color that it represents. Next, we define a linked list of the sticks that are longer than the average (which is 4), and another for the sticks that are shorter than the average. In the present case, the long-stick list is $( {}^{} \mathrm{red}^{\ast \prime} {}^{} \mathrm{green}^{\ast \prime} )$ , and the short-stick list is $( ^{\ast} \mathrm{b l u e^{; \ast}}^{\ast \ast} \mathrm{y e l l o w^{; \ast}}^{\ast} \mathrm{p i n k^{; \ast}} )$ . We pick the first item from each list and cut the longer (red) stick in two pieces in such a way that if we join one of the two to the shorter (blue) stick, we obtain an average-length stick. As a result, we are left with a red stick of length 3 and a joint stick of length 4. Since the original red stick was made shorter than the average length, we remove it from the long-stick list and append it to the short-stick list. On the other hand, the original blue stick has become an average-length stick, so we remove it from the short-stick list. Then, we pick the first item from each list and repeat the same operations again and again. When finished, we have five sticks of average length, some with a single color and others with two colors.

At each step at least one of the lists is shortened by one. Because the length of the list is finite, the whole procedure must come to an end within a finite number of operations. While constructing the joint stick-lists takes ${\mathcal{O}} ( n )$ operations, this scaling does not matter in most cases, because once we have constructed the lists we can use them many times with little additional cost.

Figure A.1 Constructing aliases in Walker’s method.

Now, to use these lists efficiently in the present case, we generate a random number $r \in [ 0 , 20 ]$ (20 is the sum of the weights). Suppose it is 9.6. Instead of accumulating and comparing as in Algorithm 1, we simply divide r by the average stick length, that is, 4. This number means that our tentative choice is the third stick. But this time “third” does not necessarily mean “yellow.” It can be an alias of “green.” In the present case, the third stick, which was originally yellow, is now a composite of a yellow piece of length 3 and a green piece of length 1 so we have to choose among these two colors with probabilities 3/4 and 1/4. What is important is that this operation is (1).

Appendix B Rejection method#

Rejection methods sample a variety of probability densities in a quite different way from the methods discussed in Chapter 2. The methods in Section 2.3 always returned a sample for each random number drawn. In a rejection method some draws do not return a sample: They are rejected. There are a number of rejection methods, some of which target specific densities. Of the general methods, von Neumann’s is the most famous. It is quite general purposed, works in more than one dimension, and does not require the distribution of interest $f ( x )$ to be normalized.

Suppose $f ( x )$ equals $c p ( x )$ where $p ( x )$ is a probability density and c is the unknown normalization of $f ( x )$ , and $g ( x )$ is another probability density. The procedure described in Algorithm 48 enables the sampling of x from $p ( x )$ . The loop generates an x until it is accepted.The other samples are rejected.

It follows directly from the until statement that because $\zeta$ is drawn from a uniform distribution over [0, 1], the probability of accepting an x is

where C is a constant such that $C g ( x ) \geq f ( x )$

Algorithm 48 Rejection method.
Input: Function $f ( x )$ to be sampled, some density $g ( x )$ , and a constant C such that
$C g ( x ) \geq f ( x )$
repeat
Sample x from $g ( x )$
Compute the ratio $\begin{array} {r} {R \left( x \right) = \frac{f \left( x \right)} {C g \left( x \right)} \left( \leq 1 \right)} \end{array}$
Draw a uniform random number $\zeta \in [ 0 , 1 ]$
until $\zeta < R ( x )$
return x.

Next, we have

We have P(accept|x) but want P(x|accept). To find it we use Bayes’s theorem:

Because x is sampled independently of $\zeta$ , the two events in the joint probability are statistically independent, and we can write

which proves the validity of the algorithm.

The efficiency of the algorithm depends on the ratio $c / C$ and the choice of $g ( x )$ . The first condition follows from (B.1): A small value of $c / C$ implies low acceptance. The second is illustrated in Fig. B.1. Here, we assumed $f ( x ) = p ( x )$ and $g ( x )$ is uniform over the sampling interval $( a , b )$ . The algorithm is a variant of blindly throwing darts at a target and taking the ratio of the number of target hits relative to the number of times the bounding box is hit. For efficiency, $g ( x )$ should “cover” $f ( x )$ tightly and cost little to evaluate.

Figure B.1 The rejection algorithm. The function to be sampled is $f ( x ) = x \sin ( x )$ The covering function $g ( x )$ is uniform.

Appendix C Extended-ensemble methods#

While extending the state space with auxiliary graph variables, as done by the loop/cluster and worm algorithms, can be a powerful way of accelerating the simulation, there are also other strategies. The approach discussed here is a form of importance sampling (Section 10.4.3). In some cases, we can accelerate a slow relaxation by using a surrogate target weight $\tilde{W} ( C )$ instead of using the given target weight $W ( C )$ . At finite temperatures, for example, we simply carry out Markov-chain Monte Carlo simulations with the transition probabilities satisfying the detailed balance condition for $\tilde{W} ( C )$ instead of $W ( C )$ . To obtain the proper thermodynamic average with respect to the original weight, we calculate

where $\langle \cdots \rangle_{\tilde{W}}$ is the Monte Carlo average with respect to the surrogate $\tilde{W} . ~ \mathrm{At}$ finite temperatures, simulation techniques based on this idea are called extendedensemble Monte Carlo methods. Into this category fall various algorithms such as the reweighting method (Ferrenberg and Swendsen, 1988), multicanonical Monte Carlo (Berg and Neuhaus, 1991), the Wang-Landau method (Wang and Landau, 2001), simulated tempering (Marinari and Parisi, 1992), replica Monte Carlo (Geyer, 1991; Hukushima and Nemoto, 1996), and the like. They differ in how the importance function $\tilde{W} ( C )$ is constructed.

For an efficient sampling, the weight $\tilde{W} ( C )$ must have a sufficiently large overlap with the target weight $W ( C )$ . If not, states that dominate the ensemble of $W ( C )$ may seldomly appear in the simulation. For example, when applied to some classical systems, the reweighting method of Ferrenberg and Swendsen (1988) uses for $\tilde{W} ( C )$ the Boltzmann weight at a temperature different from the one at which the thermodynamic average is needed, that is, $\tilde{W} ( C ) \propto e^{- E ( C ) / T^{\prime}}$ and $W ( C ) \propto e^{- E ( C ) / T}$ In this case, $\vert T^{\prime} - T \vert$ must be smaller than some energy scale proportional to $1 / \sqrt{V}$ with V being the volume of the system. The proportionality has this form because the width of the energy distribution, or the density of states profile, is proportional to ${\sqrt{V}} ,$ whereas the distance between the mean values of the two distributions is proportional to $c V | T^{\prime} - T |$ when $T^{\prime}$ and $T$ are close to each other. Here, $c$ is the specific heat, which is of the order of unity unless the system is not close to a phase transition.

In most extended-ensemble methods, we do not fix $\tilde{W} ( C )$ but instead use it throughout the simulation as an adjustable function and gradually optimize it. One of the earliest successful examples of this adaptive strategy is the multicanonical method (Berg and Neuhaus, 1991, 1992). In the multicanonical method, the weight depends on the configuration C only through the energy $E ( C )$ . Therefore, we have

and may adjust $\tilde{w} ( E )$ to make the frequency of the event $E ( C ) = E$ independent of $E .$ . This adjustment ensures that the resulting distribution function has a significant overlap with the Boltzmann distribution at any temperature. In this way, we are able to obtain information for a wide range of temperatures with only a single Monte Carlo simulation. For a number of problems, the energy takes only a discrete set of values, while for others, it takes a continuum of values, at least bounded from below. In the continuum case, we bin the energy; that is, we effectively discretize it. The Monte Carlo sampling we now describe is a random walk on the discrete values (or bins) of the energy E.

In an extended ensemble method, we perform several sets of samplings. The weight $\tilde{w}$ is fixed during each set, and it is updated at the end of each set for use in the next set. For the first set, we make an initial guess of the appropriate weight, most often taking $\tilde{w} ( E ) = \mathrm{const} .$ . In every set, the histogram $h ( E )$ is recorded; that is, every time a new state C is generated, we update the histogram (which is initially set to $h ( E ) = 0$ for all $E )$ via

We update $C$ by Metropolis sampling with the weight $\tilde{W} ( C ) \equiv \tilde{w} ( E ( C ) )$ . At the end of each set, $\tilde{w} ( E )$ is updated using

and then the histogram is reset $( h ( E ) 0 )$

The next set with this new weight starts with the last configuration of the current set.

When properly normalized, the histogram must converge to $h^{*} ( E ) \propto g ( E ) \tilde{w} ( E )$ in the limit of infinite simulation length. Here

is the density of states and $S ( E ) = k_{B}$ ln $g ( E )$ is the microcanonical entropy. Therefore, if the simulation is sufficiently long, $h ( E ) / \tilde{w} ( E )$ is a good approximation to $g ( E )$ apart from a normalization constant. With $g ( E )$ thus obtained, we compute the canonical average of an arbitrary quantity $Q$ at an arbitrary inverse temperature $\beta$ as

In this expression, $Q ( E )$ is the microcanonical average of $Q$ at discretized values of the energy E. We compute $Q ( E )$ during the last set of samplings by classifying the sampled configurations according to the energy, and taking the average of each “microcanonical” set of configurations. In other words,

While this procedure is quite useful in studying various systems with slow dynamics, there is a drawback: The range of $E$ visited by the random walk in energy space widens very slowly, partially due to the poor initial guess of the weight $\tilde{w} ( E )$ and partially due to the diffusive nature of the random walk.

The Wang-Landau variant of the method (Wang and Landau, 2001) improves the situation. This method updates $\tilde{w}$ via

every time a new state is generated. However, the whole simulation may still be organized as multiple sets of samplings, and the reduction factor $0 ~ < ~ r ~ < ~ 1$ is fixed for each set. At the beginning of each sampling, the histogram is reset, that is, $h ( E ) = 0$ for all $E ;$ the reduction factor is updated as $r \gets r^{\alpha}$ where $0 < \alpha < 1$ ; while the weight $\tilde{w} ( E )$ is kept unchanged. Each set of samplings does not have a preset duration but is terminated when the histogram becomes approximately flat. Since the dynamic update (C.2) strongly penalizes random walkers staying at the same value of the energy, the visited region widens much faster than in the ordinary multicanonical Monte Carlo method. The efficiency of the algorithm depends on the choice of the initial value of $r ,$ the value of $\alpha ,$ and the condition for terminating each set of samplings. The original work (Wang and Landau, 2001) suggested $r = 1 / e$ $\alpha = 1 / 2$ , and each set was terminated when the maximum deviation from the average in the histogram became less than 20% of the average. Obviously, there is room for further optimization of the updating scheme for these parameters and the termination condition.

It is natural to consider the combination of the extended-ensemble method with the loop/cluster algorithms. However, in a typical quantum Monte Carlo method, the weight $\tilde{W} ( C )$ is not a function of the energy, and there is no reason to assign a special role to the energy. In addition, the value of the energy is not discretized or binned as conveniently. Janke and Kappler (1995), Yamaguchi and Kawashima (2002), Yamaguchi et al. (2002), and Troyer et al. (2003) used the extendedensemble method in the loop algorithm and in the directed-loop algorithm by histogramming the number of vertices (or graph elements). In other words, they replaced the factor $\beta^{n}$ in the high-temperature series expansion by a function f (n),

This replacement changes the vertex assignment in the directed-loop algorithm and the graph assignment in the loop/cluster algorithm. For the directed-loop algorithm, the starting expression (6.37) becomes

Here, $\begin{array} {r} {n ~ = ~ \sum_{k l} \gamma_{k l}} \end{array}$ is the number of vertices in $G \equiv \{\gamma_{k l} \}$ and $1 / M^{n}$ replaces $( M - n ) ! / M !$ . The replacement of $\beta^{n}$ by $f ( n )$ affects the weight only as a function of the graph G. It does not change the weight as a function of $C \equiv \{\psi ( k ) \}$ . Correspondingly, in the resulting extended-weight directed-loop algorithm, the worm updating cycle is exactly the same as before with only the graph assignment process modified.

With respect to the graph assignment, we need to modify the first two for-loops in Algorithm 20 to update the vertices for each uniform time-interval terminated by kinks at both ends. Let I be one such interval. From (C.3) it follows that the probability of having ν vertices in this interval is proportional to

where $n_{0}$ is the total number of vertices in the rest of the system. The binomial coefficient is the number of ways to assign ν vertices to the $I / ( \Delta \tau )$ cells. Taking the limit of $\Delta \tau = \beta / M \to 0$ , we obtain the unnormalized probability of having ν vertices as

Therefore, to assign vertices in I, we could in principle generate a random integer ν from this distribution and assign ν vertices to positions sampled from a uniform distribution over I. However, doing this is often unpractical, as computing the normalization constant is generally expensive. An alternative is to first choose with probability $\frac{1} {2}$ whether we propose to increase or decrease the number of vertices by one, and then update the number of vertices by the Metropolis algorithm; namely, we accept the proposal with the probability min(1, $, w_{I} ( \nu \pm 1 ) / w_{I} ( \nu ) )$ , with ν being the current number of the vertices in I.

Extended-ensemble methods are especially useful in the computation of the quantities directly related to the density of states, such as the free energy, entropy, and specific heat, which normally would require a series of simulations at different temperatures or the numerical integration over temperature if other methods are adopted. An example of the computation of these quantities can be found in Troyer et al. (2003). Yamaguchi et al. (2002) proposed further improvements by employing the broad-histogram method (de Oliveira et al., 1998a,b; Wang et al., 1999). This method exploits the exact relation between the expectation value of the transition matrix and the density of states,

where $N ( S , E \ \to \ E^{\prime} )$ is the number of possible paths to reach any of the states with energy $E^{\prime}$ when we are currently in a state with energy E. Generally, N is a macroscopic quantity proportional to the system’s volume. This means that by using the above equation we can estimate $g ( E )$ only from the expectation values of macroscopic quantities, which we can estimate with higher precision than each entry of the histogram itself.

Appendix D Loop/cluster algorithms: SU(N) model#

As mentioned at the beginning of Section 6.2.6, in some cases the loop algorithm (Algorithm 17) for $S ~ = ~ {\frac{1} {2}}$ is applicable to spin systems with $S \ > \ \frac 12$ without splitting spins. This is always possible if the Hamiltonian decomposes graphically. For example, let us consider the bilinear-biquadratic interaction model with $S = 1$ (Harada and Kawashima, 2001, 2002),

The model (D.1) obviously has SU(2) symmetry. At $\theta = \pm \pi / 2$ and $\theta = \pm \pi / 4$ however, it possesses a higher symmetry. Consider the case of $\theta = - \pi / 2$ for which

where the constant −1 has been added for convenience. The Hamiltonian (D.2), as well as the Hamiltonian at other special values of $\theta_{i}$ , has SU(3) symmetry. In the ordinary $S^{z .}$ basis, $S_{i}^{z} | {s_{i}} \rangle = s_{i} | s_{i} \rangle$ , the matrix elements of the pair Hamiltonian are

The function $\sigma ( S_{p} ) = \pm 1$ carries a sign: It is −1 if and only if either the initial state $( s_{i} , s_{j} )$ is (0, 0) and the final state $( s_{i}^{\prime} , s_{j}^{\prime} )$ is $( \pm 1 , \mp 1 )$ or vice versa. It is easy to see that this sign is irrelevant since the negative signs always occur an even number of times, leaving the sign of the whole system positive. Therefore, in (D.3) we neglect $\sigma$ . The function $\Delta ( S_{p} , g_{\mathrm{H}} )$ is

This quantity is the horizontal graph operator in Table 6.2, generalized to the present three-state model.

With the pair Hamiltonian expressed in terms of the graph operators, it is straightforward to apply the general prescription discussed in Section 6.2.2. The only difference between the resulting S = 1 algorithm and the $\begin{array} {r} {S = {\frac{1} {2}}} \end{array}$ case is that the density of the graph elements is J, not $J / 2$ , because of the larger diagonal matrix element and that the variables have three values, $s_{i} = \{- 1 , 0 , 1 \}$ , instead of two. Accordingly, loop-flipping is not a coin toss, but proceeds by choosing one of the three values randomly with equal probability.

We can construct a similar algorithm for cases with lower symmetry, for example, SU(2) symmetry, for the parameter region $- 3 \pi / 4 \le \theta \le - \pi / 2$ . The graphical decomposition of the pair Hamiltonian in this case is (Harada and Kawashima, 2001)

$$

• H_{i j} = J \left[ ( - \sin \theta + \cos \theta ) D_{i j} ( g_{\mathrm{h}} ) - ( \cos \theta ) D_{i j} ( g_{\mathrm{d}} ) \right] ,\tag{D.4} $$

where an irrelevant sign and an additive constant have been omitted. The symbol $D_{i j} ( g_{\mathrm{d}} )$ corresponds to the diagonal graph in Fig. 6.2. Its matrix elements are

The loop construction and the loop flipping are done as before.

The special case of the SU(N) bilinear-biquadratic model that possesses SU(3) symmetry generalizes to a SU(N) symmetric model. The local spin variable of the SU(N) Heisenberg model takes N possible values, that is, $\sigma_{i} = \{1 , 2 , \ldots , N \}$ Expressed in terms of the generators of the SU(N) algebra,

Here, $S_{i}^{\mu \nu} , \bar{S}_{j}^{\nu \mu}$ satisfy

The operators $\bar{S}_{j}^{\mu \nu}$ satisfy analogous conditions. We consider the model defined on a bipartite lattice. On one sublattice (say, the sublattice A) we adopt the fundamental representation that corresponds to the Young tableau with a single box. On the other sublattice (B), we adopt the conjugate representation. The bar on the operator ${\bar{S}}_{i}^{\mu \nu}$ is a reminder that we are using different representations for the spins on sublattices A and B.

More explicitly, the matrix element of $S_{i}^{\mu \nu}$ on the sublattice A is $\left. \alpha \left| S_{i}^{\mu \nu} \right| \beta \right. =$ $\delta_{\alpha \mu} \delta_{\beta \nu} - {\textstyle \frac{1} {N}} \delta_{\alpha \beta}$ , while on B $\begin{array} {r} {\langle \alpha | \bar{S}_{j}^{\mu \nu} | \beta \rangle = - \bar{\delta}_{\alpha \nu} \delta_{\beta \mu} + \frac{1} {N} \delta_{\alpha \beta}} \end{array}$ . Therefore, we are able to rewrite (D.5) as

If we define the operator $P_{i j}$ by its matrix element as

it is the projection operator to the color singlet-state on the bond $\langle i , j \rangle$

To convince ourselves that this is indeed a singlet state, we have to recall that we are using the conjugate representation to the fundamental one for one of the two sublattices. In other words, one of the two αs in $| \alpha , \alpha \rangle$ does not represent the same state as the other $\alpha .$ but it is the product of all the states but $\alpha .$ (It is easy to confirm that this state is invariant under $P_{i j} . )$ Using the SU(N) generators, we express the projection operator as

Therefore, apart from an additive constant,

Note also that the matrix elements of the projection operator are the product of Kronecker deltas, so they can be represented by a graph element, $g_{\mathrm{h}}$ , that is,

The two spins bound by this graph element have the same value, whereas they were antiparallel previously. But this difference is removable by a gauge transformation and is not essential. Thus, we have the graphical decomposition of the SU(N) Heisenberg model Hamiltonian, from which a loop algorithm follows. The algorithm is almost identical to the one discussed in Section 6.2.3 for the $\begin{array} {r} {S = {\frac{1} {2}}} \end{array}$ antiferromagnet. The only difference is that each local variable takes N values instead of 2.

Appendix E Long-range interactions#

Most examples discussed in this book are models with on-site or nearest-neighbor interactions. Physical systems, however, often have interactions reaching beyond one lattice constant, and there are cases where such long-range interactions alter the essential properties of the system. A classic example is the dipolar interaction in magnets. Because of the long-range dipolar interaction, magnetic domains form, reducing the total magnetization. Another more recent example is the effect of dipolar interactions among ultra-cold atoms in an optical lattice. Because of the much larger lattice constant in such systems, the magnitude of the dipolar interaction relative to the nearest-neighbor interaction is enhanced compared with a conventional solid. Frustrated systems are other good examples, since for these the short-range interactions sometimes cancel, making the middle/long-range interactions more relevant.

From the computational point of view, long-range interactions generally increase the computational complexity of Monte Carlo simulations. Consider, for example, the Ising model

where $J ( R_{i j} )$ is some decreasing function of the distance between two sites $R_{i j}$ . If we adopt the single-spin update, every time we attempt flipping a spin we must compute the molecular field produced by all the other spins. This requires $\mathcal{O} ( N )$ operations for a lattice of N sites. Therefore, every Monte Carlo sweep takes a computational time proportional to $N^{2}$ instead of N, as in the case of the shortrange interaction models. A neat trick exists to reduce this $N^{2}$ dependence to NM, where M is the effective number of interacting spins defined by

If M is bounded from above by some constant that is independent of the system size, this trick reduces the computational complexity back to $\mathcal{O} ( N )$ . Even if M diverges in the thermodynamic limit, reducing the complexity from $N^{2}$ to NM can still be a major gain. For example, if the coupling constant $J ( R )$ is inversely proportional to R, M is proportional to $L^{d - 1}$ if L is the linear size and the dimension is $d \ge 3$ . In this case, $N M \propto L^{2 d - 1}$ , which is less than $N^{2}$ by a factor L.

As a first method and specific example, let us discuss the approach originally proposed by Luijten and Blote (1995), which is based on a look-up table. We¨ consider the Swendsen-Wang cluster algorithm for an Ising model with long-range interactions. The system has N sites. In the following, we first compute the probability for placing a bond assuming that all the spins are parallel, and then, before we actually place a bond, we examine whether the spins are parallel. (If not, we do not place a bond, of course.) The probability of having a bond between the p-th $( p = 1 , 2 , \ldots , N ( N - 1 ) / 2 )$ pair of spins can be written as

where $K_{p}$ describes the interaction strength for the pair $p .$ For example, if $p$ is the pair $( i , j ) , K_{p} \equiv \beta J_{i j}$ . Now consider a sequence of $N_{\mathrm{pair}} \equiv N ( N - 1 ) / 2$ boxes, each representing a pair of sites, and suppose that a ball is placed in each box with the corresponding probability (E.1). The probability that no ball is put in any box between the m-th and the n-th becomes

If we define an array (“look-up table”) A of length $N_{\mathrm{pair}}$ with elements

then (E.2) can be written as

The pairs that are to be connected can then be calculated as follows: A random number $\zeta_{1} \in [ 0 , 1 )$ is chosen and the array A is searched for the first index, say, $p_{1}$ such that $A_{p_{1}} \leq \zeta_{1}$ . For the calculation of the next bond, the values of A must be divided by $A_{p_{1}}$ , or equivalently, the random number $\zeta_{2} \in [ 0 , 1 )$ multiplied by $A_{p_{1}}$ The site number $p_{2}$ of the second connected spin is the first index, such that $A_{p_{2}} \leq$ $A_{p_{1}} \zeta_{2}$ . This process is continued until $A_{p_{k - 1}} \zeta_{k} < A_{N_{\mathrm{pair}}}$ . As was already mentioned, the outlined procedure assumes that all spins are aligned. Before actually inserting a proposed bond, we have to check whether the two spins are parallel. Because $A_{0} = 1$ , the above algorithm corresponds to the pseudocode in Algorithm 49.

430 Appendix E Algorithm 49 Cluster algorithm with long-range interactions: look-up table approach. Input: $N_{\mathrm{pair}} .$ , array A defined in (E.3). l , an initially empty list of connected pairs, $p = 0$ loop Draw random number $\zeta \in [ 0 , 1 )$ $r = A_{p} \zeta$ if $r < A_{N_{\mathrm{pair}}}$ then break ; end if Search the lowest index p such that $A_{p} \leq \zeta$ If the spins connected by pair p are parallel, add p to the list l ; end loop return the list l of connected pairs.

If the bisection method is used to search for the lowest index such that $A_{p} \leq \zeta$ the calculation time is of order (log $N_{\mathrm{pair}} ) = \mathcal{O} ( \log N )$ for each proposed bond. Since the number of bonds connecting to a given site is of order M, and we have N sites, the total calculation time for the bond placing is O(NM log N).

We can even get rid of the log N factor, although it may not matter in practical applications. To do this, let us regard a missing bond on a pair of interacting sites as a “pair survival.” If the survival probability of a pair p is $e^{- \Gamma_{p}}$ , we can interpret it as surviving the decay rate of unity for the duration $\Gamma_{p} . ~ ( \Gamma_{p} ~ = ~ 2 \beta J_{p}$ in the Swendsen-Wang algorithm.) In the case of many pairs, we have many intervals with various durations $\Gamma_{p}$ . Let us now consider a line segment of length $\begin{array} {r} {\dot{\Gamma} \equiv \sum_{p = 1}^{N_{\mathrm{pairs}}} \Gamma_{p} .} \end{array}$ Instead of judging for each individual interval whether pairs survive or not, we distribute uniformly and randomly decay events over this line segment with unit density. Specifically, we generate the total number of decay events by drawing an integer n from the Poisson distribution of mean and distribute n events uniformly and randomly on the line. This procedure is like throwing n darts onto the line. As a result, some darts may or may not fall into the interval $\Gamma_{p} . \mathrm{H}$ we have no dart in $\Gamma_{p} ,$ the pair p survives (no bond assigned); otherwise, it is killed (a bond assigned). Given that the expectation value of n is proportional to , the computational cost for generating n bonds is of order . In the case of the Swendsen-Wang algorithm this cost is proportional to $\sum_{i j} J_{i j}$ and hence O(NM).

Finding the intervals hit by any dart is less trivial. We carry out the task n times (for each of the randomly distributed darts). Generating a uniform random number is easy and costs $\mathcal{O} ( 1 )$ for each dart thrown. The difficult part is finding the index p of the interval hit. Using the technique discussed in Appendix A we can do this in a time (1) (Knuth, 1997; Walker, 1977; Fukui and Todo, 2009). As a result, throwing darts and judging survival takes a computational time of $\mathcal{O} ( \Gamma ) \sim \mathcal{O} ( N M )$ in total.

The application of these ideas to quantum simulations is rather straightforward. For example, in the directed-loop algorithm for a long-range interaction model, we generate the total number of vertices n by the Poisson distribution with mean $\Gamma .$ defined in this case as $\Gamma = \sum_{l} \Gamma_{l}$ with

In the placement of the vertices, the algorithm, however, differs from the general procedure described above. There, the number of darts $n_{l}$ on the interval $\Gamma_{l}$ did not matter as long as $n_{l} \neq 0$ . Here, it matters since $n_{l}$ is the number of the vertices we assign to the interacting spins involved in $H_{l}$ . For each of the $n_{l}$ vertices, we generate an imaginary time in the interval $\tau \in [ 0 , \beta ]$ and with probability

decide whether we place the vertex there. In this way, we compensate for the overestimate of the graph assignment density made in (E.4) by using the maximum value.

Appendix F Thouless’s theorem#

In this appendix, we prove Thouless’s theorem (Thouless, 1961; Blaizot and Ripka, 1986). It says that if the state $| \phi \rangle$ is the Slater determinant $| \phi \rangle = a_{N_{\sigma}}^{\dagger} \cdot \cdot \cdot a_{2}^{\dagger} a_{1}^{\dagger} | 0 \rangle$ where $\begin{array} {r} {a_{i}^{\dagger} = \sum_{j} c_{j}^{\dagger} \Phi_{j i}} \end{array}$ is defined by the $N \times N^{\sigma}$ matrix $\Phi$ , then the propagation of this state $e^{- \Delta \tau c^{\dagger} M c} \left| \phi \right.$ via a noninteracting Fermion Hamiltonian is to a new Slater determinant $\left| \phi^{\prime} \right.$ , where $\left| \phi^{\prime} \right. = b_{N_{\sigma}}^{\dagger} \cdot \cdot \cdot b_{2}^{\dagger} b_{1}^{\dagger} \left| 0 \right.$ , and $\begin{array} {r} {b_{i}^{\dagger} = \sum_{j} c_{j}^{\dagger} \Phi_{j i}^{\prime}} \end{array}$ is defined by the $N \times N^{\sigma}$ matrix $\Phi^{\prime} = \exp ( - \Delta \tau M ) \Phi$ . In other words, it says that if the Hamiltonian is noninteracting, a state of noninteracting Fermions evolves into another state of noninteracting Fermions and tells us how to compute the new state via the multiplication of two matrices.

This theorem underlies implicitly or explicitly many Fermion algorithms. Our use of the Hubbard-Stratonovich transformation replaces the exponential of the Hamiltonian of a system of interacting Fermions by the exponential of a Hamiltonian of Fermions that interact only with imaginary-time-dependent auxiliary fields. With this transformation, we were able to reduce the manybody formalism to numerical methods whose core computations are matrix-matrix multiplications. The role of the theorem is most apparent in the constrained path and phase algorithms (Chapter 11) where the projection to the ground state is cast as an open-ended imaginary-time propagation where each discrete imaginarytime step evolves one Slater determinant into another. It is implicit in the finitetemperature determinant method (Chapter 7). There, building the algorithm upon a generalization of the partition function of noninteracting Fermions was most natural. Our language was that of imaginary-time propagators, and our proof of the determinant form of basic single-particle propagators evoked the theorem. While we developed the zero-temperature determinant method by analogy to the finite-temperature one, we could have developed it from the point of view of the constrained path and phase ones by simply capping the length of the projection.

In fact, the code for any one of these four methods can be easily converted to one of the other three.

The proof of the theorem is simple and mainly uses the anticommutation relation $\{c_{i} , c_{j}^{\dagger} \} = \delta_{i j}$ . To evaluate $e^{- \Delta \tau c^{\dagger} M c} \left| \phi \right.$ , we commute $e^{- \Delta \tau c^{\dagger} M c}$ and $a_{N_{\sigma}}^{\dagger}$

so

Repeating this process until the operator exponential acts on the vacuum state yields

Next, we note that $e^{- \Delta \tau c^{\dagger} M c} \left| 0 \right. = \left[ I - \Delta \tau c^{\dagger} M c + {\textstyle \frac{1} {2}} \big ( \Delta \tau c^{\dagger} M c \big )^{2} + \cdot \cdot \cdot \right] \left| 0 \right. = \left| 0 \right.$ and hence

We could stop here, but we will add multiple propagator steps because they are required for the Trotterized path integrals in several algorithms. Each propagator step inherits the Slater determinant from the previous step. Multiple steps simply use multiple applications of the theorem.

Repeating this multiplication for each exponential produces

After defining $b_{i}^{\dagger} = \sum_{j k} c_{j}^{\dagger} \left[ e^{- \Delta \tau M_{N_{\tau}}} \cdot \cdot \cdot e^{- \Delta \tau M_{2}} e^{- \Delta \tau M_{1}} \right]_{j k} \Phi_{k i} .$ , we get

Appendix G Hubbard-Stratonovich transformations#

The original Hubbard-Stratonovich transformation converts the exponential of a two-body operator into a sum of exponentials of one-body operators. Each term in the sum depends on an auxiliary field. The models typically simulated have a small number of orbitals per lattice site and density-density Coulomb interactions acting only on site. This simplicity leads to the number of auxiliary fields being at most a few times the number of lattice sites. In these cases, Hirsch’s discrete transformation (Hirsch, 1983) is more efficient than the standard continuous transformation. Here, we consider simulations of more complex models and discuss Hubbard-Stratonovich transformations for models with extended-range interactions and more general Coulomb interactions. More specifically, we discuss the use of the continuous Hubbard-Stratonovich transformation when the interactions are long ranged and present a novel discrete Hubbard-Stratonovich transformation for interactions that are not simply products of number operators. In both cases, the auxiliary fields are real. We also mention in passing other Hubbard-Stratonovich transformations used in special applications. These transformations include ones that preserve the SU(2) symmetry (Assaad, 1998) and couple the auxiliary field to fluctuations in the pairing field (Batrouni and Scalettar, 1990).

G.1 Extended interactions#

We start by considering the exponentiation of the interaction

where the subscript i represents a combination of electron spin, orbital quantum numbers, and lattice position identifiers. The interactions, described by the elements $V_{i j}$ of a symmetric N ×N matrix, are not restricted to being on site or being between nearest neighbors. We have $n_{i}^{2} = n_{i}$

With $n = ( n_{1} , n_{2} , . . . , n_{N} )^{T}$ being a N-vector of occupation numbers, the operative Hubbard-Stratonovich transformation is

where $y = ( y_{1} , y_{2} , \ldots , y_{N} )^{T}$ is a vector of auxiliary fields. The convergence of this multi-variable Gaussian integral requires V to be positive definite. If it is not, we write

and choose $V_{0}$ such that the matrix $V^{\prime}$ obtained by shifting the diagonal elements of $V \ \mathrm{by} -2 V_{0}$ is positive definite. Since $\begin{array} {r} {n^{T} n = \sum_{i} n_{i} = \sum_{i} c_{i}^{\dagger} c_{i}^{}} \end{array}$ the additional exponential we introduce is that of a one-body interaction. We accommodate it by absorption into another one-body part of the Hamiltonian.

The generalized Hubbard-Stratonovich transformation (G.2) has had the same effect as the simpler transformations in Section 7.1.1, namely, converting the interacting problem into a noninteracting one. We can thus use it in a Monte Carlo simulation, instead of the others. What are the auxiliary fields and how do we sample them?

With V positive definite, a Cholesky factorization (Meyer, 2000) exists such that $V = R^{T} R .$ , where R is an upper triangular matrix with positive diagonal elements. Because of the factorization we can write $n^{T} V n = y^{T} y$ where $y ~ = ~ R n$ . After generating a vector $\boldsymbol{\zeta}^{T} = ( \zeta_{1} , \zeta_{2} , \ldots , \zeta_{N} )$ of Gaussian random variables of unit variance and zero mean, we then obtain a sample of the vector y by using the backsubstitution method to solve the linear system of equations $\zeta = R y$ . Note that we need only N auxiliary fields.

A general two-body interaction can be put into the form

With the definition of the operator $\rho_{i j} = c_{i}^{\dagger} c_{j}$ , and the mapping of the subscript pairs ij and kl to the integers m and n, (G.3) takes the quadratic form

After ensuring the $V_{m n}$ are elements of a positive-definite matrix, we sample the auxiliary field y in the same manner as we did for (G.1).

The Gaussian Hubbard-Stratonovich transformation is useful for a variety of applications other than those with long-range interactions. An important way to improve the efficiency in many of these simulations is to subtract background terms

prior to making the transformation, leading to sampling from a shifted Gaussian (Purwanto and Zhang, 2004).

G.2 General discrete transformations#

We ease into the general discrete Hubbard-Stratonovich transformation for interactions that are not simply products of number operators by first revisiting Hirsch’s transformation from a perspective different from those that are. The new perspective is an expression of the exponential of the operators as a quadratic function of operators.

With $V_{i j} > 0$ and the relation $n_{i}^{2} = n_{i}$ , Hirsch’s transformation for the exponential of each term in the potential energy (G.1) is

with cosh $\begin{array} {r} {( \Delta \tau J_{i j} ) = \exp ( \frac{1} {2} \Delta \tau V_{i j} )} \end{array}$ . If expanded as a power series, the left-hand side of this equation becomes

where we defined

From these results follows the trivial identity

if the product of $\alpha_{i j}$ and $\beta_{i j}$ equals $\lambda_{i j}$ . Noting that

unless $\alpha_{i j} = - \beta_{i j}$ , we choose¹

Next, we comment that because the action of exp[ $\begin{array} {r} {- \sigma \Delta \tau J_{i j} \left( n_{i} - n_{j} \right) + \frac{1} {2} \Delta \tau V_{i j}} \end{array}$ $\left( n_{i} + n_{j} \right) ]$ on a Slater determinant produces another one (Appendix F), then so must the action of the operator $1 + \sigma \alpha_{i j} n_{i} + \sigma \beta_{i j} n_{j} + \alpha_{i j} \beta_{i j} n_{i} n_{j}$ . If

where

the consequence of the exponential operators acting on $| \psi_{A} \rangle$ is $| \psi_{A^{\prime}} \rangle$ , where $A^{\prime} =$ $\exp ( \chi^{+} ) \exp ( \chi^{-} ) A , \chi^{+}$ is an $N \times N$ matrix that is null except for its i-th diagonal element equaling $- \Delta \tau ( {\scriptstyle{\frac{1} {2}}} V_{i j} - \sigma J_{i j} )$ , and $\chi^{-}$ is an $N {\times} N$ matrix that is null except for its j-th diagonal element equaling $- \Delta \tau ( {\textstyle{\frac{1} {2}}} V_{i j} + \sigma J_{i j} )$ . The effect of $\exp ( \chi^{-} )$ on A is to multiply the j-th row of A by exp $- \Delta \tau ( \frac{1} {2} V_{i j} + \sigma J_{i j} ) ]$ , and the effect of $\exp ( \chi^{+} )$ on this matrix product is to multiply its i-th row by exp[ $- \Delta \tau ( \textstyle{\frac{1} {2}} V_{i j} - \sigma J_{i j} ) ]$ . It is straightforward to show that these matrix-matrix multiplications are equivalent to

where $\alpha_{i j}$ and $\beta_{i j}$ are as defined by (G.4) and (G.6). For an arbitrary matrix $X , X_{\bullet i}$ is an N-vector formed from the i-th column of X, and $X_{i \bullet}$ is an N-vector formed from its i-th row. In particular, the outer product $I_{\bullet i} I_{i \bullet}^{T}$ is an $N \times N$ null matrix except for its i-th diagonal element equaling unity. Hence, the two different-looking Hubbard-Stratonovich transformations, one using an exponential of a product of operators and the other using a quadratic polynomial of operators, are equivalent.

We now write the two-body potential energy in still another form,²

with

We can show that

where $\lambda_{k l} = V_{k k l l} - V_{k l k l}$ . Hence

and

and choose

and then ask, Do the $\sigma$ -dependent operators in (G.7) yield a single Slater determinant $| \psi \rangle_{A^{\prime}}$ if they act on $| \psi \rangle_{A} 2$

Rombouts et al. (1999a) showed that the answer to this question is yes. More specifically, if $| \psi \rangle_{A}$ is a Slater determinant defined by the $N \times M$ matrix A, then the operator $I + \alpha c_{i}^{\dagger} c_{l} + \beta c_{j}^{\dagger} c_{k} + \alpha \beta c_{i}^{\dagger} c_{j}^{\dagger} c_{l} c_{k}$ acting on $| \psi_{A} \rangle$ produces a Slater determinant $\left| {\psi_{A^{\prime}}} \right.$ defined by the $N \times M$ matrix

Equations (G.7) and (G.8) define a discrete Hubbard-Stratonovich transformation for the general two-body interaction (G.2). After the Monte Carlo procedure selects an auxiliary field σ , (G.9) defines a simple procedure for updating the Slater determinant.

We now express this updating in a more familiar form. We start by letting $\Delta =$ $\alpha I_{\bullet i} I_{k \bullet}^{T} + \beta I_{\bullet j} I_{l \bullet}^{T}$ . If we are interested in the ratio of the overlap integrals of two Slater determinants after and before a proposed change, we have

where $G = R ( L^{T} R )^{- 1} L^{T}$ is the single-particle Green’s function matrix (Appendix I). Because of the sparsity of $\Delta$ , (G.10) reduces to a determinant of a $2 \times 2$ matrix

a result we could have obtained from (G.10) by using Wick’s theorem. Using the Sherman-Morrison formula (Meyer, 2000; Press et al., 2007)

twice yields an efficient procedure for updating the Green’s function: First, we execute

and then

Appendix H Multi-electron propagator#

In this appendix, we derive the expressions for the multi-electron propagator heavily used in Chapters 7 and 11 and in Appendix I. We choose to give the details for the derivations to introduce various techniques that are useful for other analyses and to engender a keener appreciation of the physical content of the multipropagator relation.

Our main objective is proving the identity

which expresses the multielectron propagator as a determinant of a matrix of singleelectron propagators

As a preliminary step, we prove a related result for electrons in nonoverlapping states, whose creation and destruction operators satisfy

where $\begin{array} {r} {B_{i j} \equiv \langle 0 | c_{i} c_{j}^{\dagger} | 0 \rangle} \end{array}$ is the overlap. The result we first prove is

To demystify this result let us remark that if the electrons were in orthonormal states, $B_{i j} = \delta_{i j}$ , the determinant would equal a product of δ-functions and would simply express the normalization of the many-electron basis.

We first note that (H.4) is trivially true for $N_{\sigma} = 1$ . We next assume that it is true for $N_{\sigma} - 1$ . Then, in the inner product on the left-hand side of (H.3), we commute electron operators until we move $c_{j_{N_{\sigma}}}$ all the way to the right, so it operates directly on the vacuum state |0 and produces 0. The anticommutation relation (H.3) says this sequence of commutations generates

By assumption, the inner products on the right-hand side of the above equation are determinants for $N_{\sigma} \mathrm{~ - ~} 1$ electrons, and accordingly the right side is the cofactor expansion of the determinant we would obtain by expanding down the far righthand column of the matrix for the $N_{\sigma}$ electron system. Thus, we have proven the identity (H.4).

To start the proof of (H.1), we consider the imaginary-time propagation of an $N_{\sigma}$ -state piecewise in $\Delta \tau$ steps:

Because of the Hubbard-Stratonovich transformation, the electron operator $H ( \ell \Delta \tau )$ is a quadratic form in the electron creation and destruction operators. Designating the $N \times N$ matrix that defines this form as $M_{\ell}$ and defining vectors of creation and destruction operators as

yields

Our task now is propagating the right-hand $N_{\sigma} {\mathrm{- b o d y}}$ state $| i_{N_{\sigma}} , \dots , i_{2} , i_{1} \rangle =$ $c_{i_{N_{\sigma}}}^{\dagger} \cdots c_{i_{2}}^{\dagger} c_{i_{1}}^{\dagger} | 0 \rangle$ in (H.1) by the product of these exponentials. We first prove that

where $\begin{array} {r} {c^{\dagger} M c = \sum_{i j} c_{i}^{\dagger} M_{i j} c_{j}} \end{array}$ . For this, we note that

or more compactly that $\left( c^{\dagger} M c \right) c^{\dagger} = c^{\dagger} \left[ M + I \left( c^{\dagger} M c \right) \right]$ . Using this commutation relation n times, we obtain

Then, using the power series representation of the exponential of an operator, $\begin{array} {r} {\exp ( A ) = I + A + \frac{1} {2 !} A^{2} + \cdot \cdot \cdot} \end{array}$ , we find $e^{- \Delta \tau c^{\dagger} M c} c^{\dagger} = c^{\dagger} e^{- \Delta \tau M - \Delta \tau \left( c^{\dagger} M c \right) I}$ and, noting that M and $c^{\dagger} M c$ commute, the relation (H.6).

To get (H.1), we make a proof by induction similar to the one used to go from (H.3) to (H.4). We first note that after a Trotter breakup (H.5) and an application of Thouless’s theorem (Appendix F)

equals (F.1) and

We now define

With this definition, we see that the inner product in (H.1) is trivially true for $N_{\sigma} = 1$ :

We next assume that it is true for $N_{\sigma} - 1$ . Then, in the inner product on the left-hand side, we commute electron operators until we move $c_{j_{N_{\sigma}}}$ all the way to the right so

it operates directly on the vacuum state |0 and produces 0. The relation (H.7) says that this sequence of commutations generates

where we left out the τ dependence of the B matrix elements for notational simplicity. By assumption, the inner products on the right-hand side are determinants for $N_{\sigma} - 1$ electrons, and accordingly, the right side is the cofactor expansion of the determinant we would obtain by expanding down the far right-hand column of the matrix for the $N_{\sigma} {\mathrm{- e l e c t r o n}}$ system. Thus, we have proven the identity (H.1).

Appendix I Zero-temperature determinant method#

While the zero-temperature determinant algorithm projects to the ground state from a trial state, it does so in a way different from that of the power methods discussed in Chapter 10. It focuses on (7.1) for a large but fixed value of $\beta$ and is a procedure to estimate expectation values

Usually, $\begin{array} {r} {\tau = \frac{1} {2} \beta} \end{array}$ . One may also calculate

that is, average the expectation values over a small interval around ${\scriptstyle{\frac{1} {2}}} \beta$ to reduce the variance. It is also possible to extend the integration so it ranges from 0 to $\beta$ . Our purposes are served if we take $\begin{array} {r} {\tau = \frac{1} {2} \beta} \end{array}$

The simulations are done for $N_{\sigma}$ fixed. Almost always, $| \psi_{L} \rangle = | \psi_{R} \rangle$ . These states are usually chosen to be a superposition of basis states for fixed $N^{\sigma}$ , for example, the noninteracting or the Hartree-Fock solution for the system being simulated.

Finite-sized systems of electrons display shell effects caused by degenerate states being separated by energy gaps. For closed shells, electrons occupy all states below some gap and hence are in a nondegenerate state. In this case, the $| \psi_{L} \rangle$ and $| \psi_{R} \rangle$ are often adequately chosen to be the noninteracting state. For open shells, not all states below the gap are occupied, and hence the state is degenerate. In this case it is better to choose $| \psi_{L} \rangle$ and $| \psi_{R} \rangle$ as a symmetry-adapted linear combination $\sum_{\phi} c_{\phi} | \phi \rangle$ of the noninteracting degenerate states. Then, even if $| \psi_{L} \rangle = | \psi_{R} \rangle$ , we need the expectation values of operators between a $| \phi^{\prime} \rangle$ and $| \phi \rangle$ that differ. In the following, we assume that the $| \psi_{L} \rangle = | \psi_{L}^{\uparrow} \rangle | \psi_{L}^{\downarrow} \rangle$ and $| \psi_{R} \rangle = | \psi_{R}^{\uparrow} \rangle | \psi_{R}^{\downarrow} \rangle$ are represented by single but different Slater determinants defined by $N_{\sigma} \times N$ and $N \times N_{\sigma}$ matrices $\Phi_{L}^{\sigma}$ and $\Phi_{R}^{\sigma}$ . For example,

where

As is standard, we express the denominator of (I.1) as a path-integral (7.1), make the Trotter approximation and the Hubbard-Stratonovich transformation (7.9), and then recycle the analysis of Section 7.1.2 to find

With the definitions

we can show that

which reduces to (7.43). In zero-temperature fixed $N^{\sigma}$ calculations, $L^{\sigma}$ and $R^{\sigma}$ are $N^{\sigma} \times N$ and $N \times N^{\sigma}$ rectangular matrices. We still update the Green’s function matrix by (7.45).

For the measurement of the equal-time Green’s function, an analysis similar to that in Section 7.2.2 yields the formula¹

Here, in contrast to the finite-temperature algorithm, τ should be interpreted not as an imaginary time, but rather as a Trotter index. For “equal time” measurements nothing changes, while unequal-time measurements are unphysical. If we move from one Trotter step to a later one, we can use relations such as

At any Trotter index, the Green’s function matrix is primarily composed of two strings of $B^{\sigma}$ matrices, one whose number of factors is decreasing from the left and another whose number of factors is increasing from the right as imaginary time advances. Creating and storing the sequences of partial matrix products (and their factorizations) accelerates the simulation. We can store the partial products for each spin $\sigma$ as indicated in the following equation:

The partial products generated in the computation of $L^{\sigma} ( \tau )$ and $R^{\sigma} ( \tau )$ are stored. To move from time τ to time $\tau + \Delta \tau$ we need the factors $L^{\sigma} ( \tau + \Delta \tau )$ and $R^{\sigma} ( \tau + \Delta \tau )$ . The left factor $L^{\sigma} ( \tau + \Delta \tau )$ is already available. The right factor $R^{\sigma} ( \tau + \Delta \tau ) = B^{\sigma} ( \tau + \Delta \tau , \tau ) R^{\sigma} ( \tau )$ is easily computed. $R ( \tau )$ may be stored in lieu of $L^{\sigma} ( \tau )$ , which is no longer needed. With this storage scheme, we can create the matrix needing inversion just with a few matrix multiplications. The cost of repeatedly stabilizing the matrix inversion becomes minor. Storing partial products in this manner has the added advantage of reducing the accumulation of round-off errors. Instead of storing the partial products for each Trotter index, doing the storage every so many steps is still advantageous. For example, in the above we replace $\Delta \tau \to m \Delta \tau$

The matrix factorization procedure for doing the matrix inversion is a bit different. Using the factored forms

in the expression for the Green’s function

we find that

Algorithm 50 Modified Gram-Schmidt method: matrix column orthonormaliza tion. Input: $\overline{{M \times N}}$ matrix U. for k = 1 to N do $d = 0$ ; for i = 1 to M do $d \gets d + U ( i , k ) U ( i , k )$ end for $d \sqrt{d}$ ; for i = 1 to M do $U ( i , k ) \gets U ( i , k ) / d$ end for for $j = k + 1$ to N do v = 0 ; for i = 1 to N do $\nu \nu + U ( i , k ) U ( i , j )$ end for for i = 1 to M do $U ( i , j ) \gets U ( i , j ) - \nu U ( i , k )$ end for end for end for return the column-orthonormalized matrix U.

The $U_{L}^{\sigma}$ and $U_{R}^{\sigma}$ are $N^{\sigma} \times N$ and $N \times N^{\sigma}$ rectangular matrices, while the $D_{L}^{\sigma} , V_{L}^{\sigma}$ $D_{L}^{\sigma}$ , and $V_{L}^{\sigma}$ matrices are $N^{\sigma} \times N^{\sigma}$ . The Green’s function, however, depends only on the rectangular $U_{L}^{\sigma}$ and $U_{R}^{\sigma}$ matrices, so we can stabilize the calculation (Algorithm 50) just by orthonormalizing the rows of $L^{\sigma}$ and columns of $R^{\sigma}$ . Sugiyama and Koonin (1986) and Sorella et al. (1989) were among the first to note the importance of matrix column orthonormalization.

The zero-temperature algorithm is an elegant complement to the finitetemperature one. Unfortunately, it too has a sign problem. The more recent lattice fixed-node, constrained-path, and constrained-phase methods (Chapter 11) are useful in taming it.

Appendix J#

Anderson impurity model: chain representation#

This appendix explains the mapping of the Anderson impurity model (8.8)–(8.11) onto a semi-infinite chain. For simplicity, we may assume that the hybridization parameters $V_{p}$ are real. The goal is to transform the impurity and bath annihilation operators $\{d , c_{p_{1}} , c_{p_{2}} , \ldots \}$ to new operators $\{d , c_{1} , c_{2} , \ldots \}$ such that $H_{0} + H_{\mathrm{bath}} + H_{\mathrm{mix}}$ becomes tridiagonal:

The symmetric and orthogonal matrix U associated with this tranformation can be constructed as a product of Householder transformations

which describe mirror operations. Let us assume that the bath has N sites. In the first step, we choose v such that the vector $\vec{V} = ( V_{p_{1}} \ : V_{p_{2}} \ : V_{p_{3}} \ : \cdot \cdot \cdot V_{p_{N}} )^{T}$ is mapped onto the first element ${\vec{e}}_{1} \colon{\vec{\nu}}_{1} = {\vec{V}} \mp | V | {\vec{e}}_{1} \Rightarrow h_{\nu_{1}} {\vec{V}} = \pm | V | {\vec{e}}_{1}$ , with $\begin{array} {r} {| V | ^{2} = \sum_{p = 1}^{N} | V_{p} | ^{2}} \end{array}$ Defining the block matrix

with diagonal blocks of dimension 1 and N, respectively, we obtain the transformed Hamiltonian

We can now repeat the procedure, and define a transformation $h_{\nu_{2}}$ that maps the vector defined by the off-diagonal elements in the first column of the $N \times N$ matrix $A_{2}$ onto $\vec{e}_{2}$ . The Householder matrix

with diagonal blocks of dimension 2 and $N - 1$ may then be used to bring the Hamiltonian into tridiagonal form up to the second row and column. The transformation U for an impurity model with N bath sites eventually becomes

In the case of multi-orbital systems with diagonal baths, similar transformations map the system onto a star-shaped geometry, with the impurity site at the center.

Appendix K#

In this appendix, we derive the effective action for the Anderson impurity model (8.8)–(8.11). We denote the impurity annihilation operators by d, the bath annihilation operators by c, and an arbitrary type of operator by a. The partition function of the impurity model is

To derive the action, we introduce the coherent states

which are eigenstates of the destruction operators: $a_{\alpha} | \xi \rangle = \xi_{\alpha} | \xi \rangle$ . The $\xi_{\alpha}$ are elements of a Grassmann algebra, which allows us to account for the anticommutation relations of the Fermionic operators in a path-integral formulation. Specifically, the $\xi_{\alpha}$ satisfy

Three important relations of a Grassman algebra, which we will need in the following, are the scalar product, the closure relation, and the formula for Gaussian integrals:

Using the coherent states and the anticommutation of Grassmann numbers, we express the partition function as

where $\{| n \rangle \}$ is some basis of the impurity model Hilbert space. In the last step, we introduced vectors of Grassmann variables to simplify the notation.

By repeatedly inserting the closure relation (K.2) and introducing antiperiodic boundary conditions, we find

with $\Delta \tau = \beta / M$ . In the limit $\Delta \tau \to 0$ the exponent becomes

Hence, we obtain the expression for the partition function as a path integral over Grassmann variables:

To formally define the effective action for the impurity, we distinguish between the Grassmann variables $\xi_{d}$ associated with the impurity and $\xi_{c}$ associated with the bath and write

The Hamiltonian of our impurity model is of the form

so the action explicitly reads

where in the second line we have introduced the matrix $M_{i j} = \delta_{i j} ( - \partial_{\tau} - \varepsilon_{i} )$ and the vectors $\begin{array} {r} {\vec{J}_{\sigma}^{*} ( \tau ) = \xi_{d \sigma}^{*} ( \tau ) \vec{V}_{\sigma}} \end{array}$ and $\vec{\cal J}_{\sigma} ( \tau ) = \vec{V}_{\sigma}^{*} \xi_{d \sigma} ( \tau )$

We now use (K.3) to integrate out the bath degrees of freedom:

The coupling to the bath leads to a retarded term in the effective action, which is given by the hybridization function $\Delta_{\sigma} ( \tau )$ defined in (8.14). In fact, the Matsubara coefficients of this function are

Inserting the definition of $H_{\mathrm{loc}} = U n_{\uparrow} n_{\downarrow} - \mu ( n_{\uparrow} + n_{\downarrow} )$ , we finally obtain

$\mathcal{G}_{0 \sigma}$ is the “Weiss Green’s function” (8.15), which is related to the hybridization function by

The Grassmann path integral expression for the partition function is therefore

Alternatively, we can express the partition function in terms of the operators $d , d^{\dag}$ , and the hybridization function $\Delta$ :

Appendix L Continuous-time auxiliary-field algorithm#

This appendix presents a continuous-time impurity solver that combines the weakcoupling expansion with an auxiliary-field decomposition of the interaction term (Gull et al., 2008). This continuous-time auxiliary-field method (CTAUX) is an adaptation of an algorithm by Rombouts et al. (1999b) for Fermion lattice models. It is suitable for models with density-density interactions and in this case is equivalent to the weak-coupling continuous-time method of Section 8.4.

We present the algorithm for the Anderson impurity model (8.8)–(8.11). Following Rombouts and collaborators, we define $H_{2} = H_{U} - \textstyle{\frac{1} {7}} U ( n_{\uparrow} + n_{\downarrow} ) - K / \beta$ and $\begin{array} {r} {H_{1} = H - H_{2} = H_{\mu} + \frac{1} {2} U ( n_{\uparrow} + n_{\downarrow} ) + H_{\mathrm{bath}} + H_{\mathrm{mix}} + K / \beta} \end{array}$ , with K being some arbitrary nonzero constant. After the expansion of the partition function in $H_{2} , ( 8.30 )$ gives the statistical weight of a configuration of n “interaction vertices.” We now extend our configuration space by decoupling each interaction vertex using the discrete Hubbard-Stratonovich decoupling formula (Hirsch, 1983)

$$

• H_{2} = K / \beta - U ( n_{\uparrow} n_{\downarrow} - \frac{1} {2} ( n_{\uparrow} + n_{\downarrow} ) ) = \frac{K} {2 \beta} \sum_{s = - 1 , 1} e^{\gamma s ( n_{\uparrow} - n_{\downarrow} )} , $$

with

The configuration space becomes the collection of all possible auxiliary spin configurations $C = \{( \tau_{1} , s_{1} ) , \ldots , ( \tau_{n} , s_{n} ) \} , n = 0 , 1 , \ldots , \tau_{i} \in [ 0 , \beta ) , s_{i} = \pm 1$ on the imaginary-time interval $[ 0 , \beta )$ . These configurations have the weight

We first separate the electron spin components and then proceed to the analytical calculation of the trace. Introducing $\begin{array} {r} {H_{1}^{\sigma} = - ( \mu - \frac{1} {2} U ) n_{\sigma} + \sum_{p} \varepsilon_{p} c_{p \sigma}^{\dagger} c_{p \sigma} +} \end{array}$ $\begin{array} {r} {\sum_{p} ( V_{p \sigma} d_{\sigma}^{\dagger} c_{p \sigma} + V_{p \sigma}^{*} c_{p \sigma}^{\dagger} d_{\sigma} )} \end{array}$ , we factorize the trace in (L.2) as

Then, using the identity $e^{\gamma s \sigma n_{\sigma}} = e^{\gamma s \sigma} d_{\sigma}^{\dagger} d_{\sigma} + d_{\sigma} d_{\sigma}^{\dagger} = e^{\gamma s \sigma} - ( e^{\gamma s \sigma} - 1 ) d_{\sigma} d_{\sigma}^{\dagger}$ , we express the trace factors in terms of noninteracting impurity Green’s functions $\tilde{\mathcal{G}}_{0}$ (with chemical potential shifted as $\mu \to \tilde{\mu} = \mu - {\textstyle \frac{1} {2}} U )$

where $\tilde{Z}_{0}^{\sigma} = \mathrm{Tr} [ e^{- \beta H_{1}^{\sigma}} ]$ . Here, $N_{\sigma}$ represents an $n \times n$ matrix defined by the location of the decoupled interaction vertices, the spin orientations, and the noninteracting Green’s functions:

where we used the notations $e^{\Gamma_{\sigma}} \equiv \mathrm{diag} ( e^{\gamma \sigma s_{1}} , \ldots , e^{\gamma \sigma s_{n}} ) , [ \tilde{\mathcal{G}}_{0}^{\sigma} ]_{i j} = \tilde{\mathcal{G}}_{0}^{\sigma} ( \tau_{i} - \tau_{j} )$ for $i \neq j , [ \tilde{\mathcal{G}}_{0}^{\sigma} ]_{i i} = \tilde{\mathcal{G}}_{0}^{\sigma} ( 0^{+} ) . ^{1}$ Combining (L.2), (L.3), and (L.4) yields

for the weight of configuration $C = \{( \tau_{1} , s_{1} ) , . . . , ( \tau_{n} , s_{n} ) \}$

For ergodicity, it is sufficient to insert and remove spins with random orientations at random times. The sampling procedure with fast matrix updates is thus completely analogous to that of the weak-coupling algorithm discussed in Section 8.4.

To compute the contribution of a configuration C to the Green’s function, $G_{C}^{\sigma} ( \tau )$ we insert in (L.2) a creation operator $d_{\sigma}^{\dagger}$ at time 0 and an annihilation operator $d_{\sigma}$ at time τ and divide by $w_{C}$ . Wick’s theorem then leads to

As in the weak-coupling algorithm, it is advantageous to accumulate the quantity

by binning the time points $\tilde{\tau}$ on a fine grid. Then, after the simulation, the Green’s function is computed from the convolution integral

Let us briefly address the role of the parameter K and its effect on the expansion order. It follows from (8.69) that $\begin{array} {r} {\langle - H_{2} \rangle = \frac{1} {\beta} \int_{0}^{\beta} d \tau \langle - H ( \tau ) \rangle = \frac{1} {\beta} \langle n \rangle} \end{array}$ , and because $\begin{array} {r} {\langle - H_{2} \rangle = K / \beta - U \langle n_{\uparrow} n_{\downarrow} - \frac{1} {2} ( n_{\uparrow} + n_{\downarrow} ) \rangle} \end{array}$ , the average perturbation order $\langle n \rangle$ is related to the parameter K and the potential energy $E_{\mathrm{pot}} = U \langle n_{\uparrow} n_{\downarrow} \rangle$ by

Increasing K leads to a higher perturbation order (and thus slower matrix updates), but through (L.1) also to a smaller value of $\gamma$ and thus to less polarization of the auxiliary spins. The average perturbation order grows proportionally to $U$ and the inverse temperature.

For models with density-density interaction, the weak-coupling algorithm discussed in Section 8.4 is equivalent to the continuous-time auxiliary-field method (Mikelsons et al., 2009). The weak-coupling algorithm is based on the following representation of the interaction term (see (8.31)):

The second term on the right-hand side can be absorbed into a shift of the chemical potential $\mu \tilde{\mu} = \mu - {\textstyle \frac{1} {2}} U$ (with corresponding redefinition of the noninteracting Green’s function $\mathcal{G}_{0} \tilde{\mathcal{G}}_{0}$ and of the noninteracting partition function $Z_{0} \tilde{Z}_{0} )$ , and for future convenience we write the third term as $K^{*} / \beta$ , so that

With this definition, the weight of a weak-coupling configuration becomes

with $[ M_{\sigma}^{- 1} ]_{i j} = \tilde{\mathcal{G}}_{0}^{\sigma} ( \tau_{i} - \tau_{j} ) - \textstyle \frac 12 \delta_{i j} - ( \textstyle \frac 12 + \delta ) s_{i} \sigma \delta_{i j}$ . This expression should be compared with the weight of a corresponding configuration in the CTAUX formalism, which reads

with $[ N_{\sigma}^{- 1} ]_{i j} = [ e^{\Gamma_{\sigma}} - \tilde{\mathcal{G}}_{0}^{\sigma} ( e^{\Gamma_{\sigma}} - I ) ]_{i j} = e^{\gamma \sigma s_{i}} \delta_{i j} - \tilde{\mathcal{G}}_{0}^{\sigma} ( \tau_{i} - \tau_{j} ) ( e^{\gamma \sigma s_{j}} - 1 )$ . The parameters $\gamma$ and K are related by cosh $( \gamma ) = 1 + \beta U / 2 K$

To derive the relationship between the weak-coupling parameter δ and the CTAUX parameter $\gamma_{:}$ , which makes the two methods identical, we manipulate the matrix $N_{\sigma}^{- 1}$ with the goal of bringing it into a form analogous to that of $M_{\sigma}^{- 1}$ . The first step is to rewrite

The expression in the curly brackets is a diagonal matrix with elements $1 / ( 2$ tanh $( \gamma \sigma s / 2 ) ) = \sigma s / ( 2 \operatorname{tanh} ( \gamma / 2 ) )$ . The last identity is true because $\sigma s = \pm 1$ . We furthermore note that $\textstyle \prod_{\sigma} ( I - e^{\Gamma_{\sigma}} )$ is a diagonal matrix with elements $2 - 2 \cosh ( \gamma ) =$ $- \beta U / K$ . This means that the last factor in $\left( \mathrm{L} . 5 \right)$ contributes a factor $( - \frac{\beta U} {K} )^{n}$ to the product of determinants. Hence, we obtain

By comparing the right-hand side with the expression for $w_{C}^{\mathrm{w e a k - c o u p l i n g}} .$ we conclude that the two algorithms produce configurations $C = \{( \tau_{1} , s_{1} ) , \dots ( \tau_{n} , s_{n} ) \}$ with identical weights if the parameters $\delta$ and $\gamma$ satisfy the condition

or, equivalently, if ${\cal{K}} = \beta U [ ( \textstyle{\frac{1} {2}} + \delta )^{2} - \textstyle{\frac{1} {4}} ] = {\cal{K}}^{*}$ . Note that $\delta \ge 0$ implies $K \geq 0$

Appendix M Continuous-time determinant algorithm#

In this appendix, we describe the continuous-time version of the determinant method, which has been developed by Iazzi and Troyer (2014). This algorithm combines the efficient linear-in- $\boldsymbol{\cdot} \beta$ scaling of the determinant algorithm (Chapter 7) with a continuous-time sampling of interaction vertices analogous to the weakcoupling impurity solver described in Section 8.4 of Chapter 8.

First, let us recall that in the discrete-time determinant method we transformed the partition function via a Trotter approximation and Hubbard-Stratonovich transformation into

with

and

a string of multielectron propagators for some noninteracting Hamiltonian. The continuous-time formulation is based on a similar determinant, but avoids a fixed time-discretization and Hubbard-Stratonovich decoupling. Instead the times in the multi-electron propagators are generated by a weak-coupling expansion and are sampled stochastically. We now describe the method for finite-temperature simulations, following the presentation in the original paper (Iazzi and Troyer, 2014).

As in the case of the continuous-time impurity solvers, we switch to an interaction representation. The Hamiltonian of the lattice model is split into $H = H_{0} + H_{1}$ ， where $H_{0}$ is the noninteracting part of the Hamiltonian and $H_{1}$ is the interacting part. In this representation,

where $H_{1} = \textstyle \sum_{i} V_{i}$ and $V_{i}$ is some site-dependent two-body interaction in a lattice of N sites. To be specific, let us consider the Hubbard model, for which $V_{i} ~ =$ $\begin{array} {r} {U ( n_{i \uparrow} - \frac{1} {2} ) ( n_{i \downarrow} - \frac{1} {2} )} \end{array}$ . For future convenience, we also introduce the dimensionless vertex $\nu_{i} = 4 ( n_{i \uparrow} - \textstyle{\frac{1} {2}} ) ( n_{i \downarrow} - \textstyle{\frac{1} {2}} )$ . As in the case of the weak-coupling continuous-time solver, we now expand the partition function in powers of the interaction term

where

and $C = \{( i_{1} , \tau_{1} ) , . . . , ( i_{k} , \tau_{k} ) \}$ is a configuration of k interaction vertices at spacetime coordinates $( i_{\ell} , \tau_{\ell} )$

Since the propagation from one vertex to the next is given by the noninteracting part of the Hamiltonian, the weight can (in a manner analogous to the method described in Chapter 7) be expressed as a determinant,

where $\nu ( i_{\ell} , \tau ) = e^{H_{0} \tau} \nu_{i_{\ell}} e^{- H_{0} \tau}$ , and the $2 \times 2$ matrices $H_{0}$ and $\nu_{i_{\ell}} = e^{i \pi ( n_{i_{\ell} \uparrow} + n_{i_{\ell} \downarrow} )}$ have the elements

The fact that the weights are such determinants allows us to use measurement formulas for Green’s functions and other observables that are similar to those used in the discrete-time determinant methods. The Monte Carlo sampling and updating, however, are analogous to the procedures used in the weak-coupling impurity solver.

The factor $( - U / 4 )^{k}$ in (M.1) introduces a sign problem, which, as in the case of the weak-coupling impurity solver, we can avoid by introducing auxiliary fields that shift the chemical potential for spin up and down or by expanding in a shifted interaction term, as in the CTAUX method (Appendix L). Following Iazzi and Troyer (2014), we choose the latter approach and expand the partition function in powers of $\begin{array} {r} {V_{i} = - U ( 1 - n_{i \uparrow} n_{i \downarrow} ) \equiv - \frac{1} {2} U \nu_{i}} \end{array}$ . Introducing the auxiliary Ising variable α (Appendix G), we then write

The configuration $C = \{( i_{1} , \tau_{1} , \alpha_{1} ) , . . . , ( i_{k} , \tau_{k} , \alpha_{k} ) \}$ becomes a collection of k auxiliary spin variables, and the matrix for the auxiliary spin dependent vertex becomes

In complete analogy to the sampling discussed in Section 8.4, the Monte Carlo sampling of the configurations C proceeds by random insertion and removal of vertices. For example, if the initial configuration C has k vertices and the new configuration with $k + 1$ vertices is $C^{\prime}$ , then

The computational bottleneck is the updating of the matrix $B ( C , \beta )$ . To do this calculation in a time $\mathcal{O} ( \beta N^{3} )$ , Iazzi and Troyer (2014) propose to switch to the eigenbasis of $H_{0} \colon H_{0} = O E O^{\dagger}$ , with E a diagonal matrix and O being an orthonormal transformation. In this basis, the time-evolution operators become diagonal, while the vertex operators become an identity matrix plus an outer product of two vectors: $\delta_{k k^{\prime}} + \gamma O_{k i_{\ell}}^{\dagger} O_{i_{\ell} k^{\prime}}$ . As in the case of the finite-temperature determinant algorithm (Chapter 7), the calculation of the matrix products requires a stabilization procedure.

Equal-time observables are expressible as functions of the matrix B. For example, the equal-time Green’s function $\langle c_{i \sigma} c_{j \sigma^{\prime}}^{\dagger} \rangle$ is obtained by computing the matrix elements of

The time-dependent Green’s function at time $\tau$ is measured using the partial propagator $B ( C , \tau ) ~ = ~ e^{- \tau H_{0}} \prod_{\tau_{i} < \tau} \nu ( i_{\ell} , \tau_{\ell} )$ . In this case, we accumulate the matrix elements of

The continuous-time version of the zero-temperature determinant algorithm (Appendix I) has been discussed by Wang et al. (2015).

Appendix N Correlated sampling#

Correlated sampling is a technique for reducing variance. It is sometimes called differential Monte Carlo because its principal use is for reducing the variance of the average difference between two quantities. A related use is estimating derivatives. The method exploits the fact that by promoting a statistical dependence between variables, often simply by having estimates share all or part of a sequence of random numbers, the total variance is less than the sum of the variances.

To explain the basic idea, suppose we seek a Monte Carlo estimate of the integral

where $p_{X} ( x )$ is the probability distribution of a random variable X. In other words, $F = \langle f ( x ) \rangle$ . We now consider another random variable Y, which is not necessarily independent of X, and construct a third random variable

where for the time being $g ( Y )$ is arbitrary. If $G = \langle g ( y ) \rangle$ , then we can express the variance of $\xi$ as

where $\rho$ is the correlation coefficient

As noted in Section $3.6 , - 1 \leq \rho \leq 1$ , and if X and Y are statistically independent, that is, uncorrelated, $\rho = 0$ and the variances add.

If G is known, then since

we could estimate F by estimating ξ. From (N.1), we see that the variance of $\xi$ is small if $\rho$ is close to unity. In fact, if $X = Y$ and $g = f$ , that is, the expectation values are perfectly correlated, then $\rho ~ = ~ 1$ and the variance of $\zeta$ equals zero. Consequently, when $X = Y$ and $g \approx f$ , it is more advantageous to sample $\xi$ than $X .$

We comment that correlated sampling is often used in the following way: Suppose the random variable Y is correlated with the random variable X and we can compute cov(X, Y) (3.23). If we want to estimate X, know Y, and compute $\operatorname{var} ( Y )$ , then we construct the random variable $X ( a ) = X + a ( Y - \langle Y \rangle )$ . Noting that X(a) and X have the same mean, we write

If we take $a = \operatorname{cov} ( X , Y ) / \operatorname{var} ( Y )$ , then

Finally, we note that if $\left. X \right. = \left. Y \right.$ , then we can always choose an a such that the variance of $X ( a ) = a X + ( 1 - a ) Y$ is smaller that the variance of X.

In both the energy and energy-variance minimization variational Monte Carlo methods (Section 9.3), the basic step is generating a sequence of trial states $| \psi_{T}^{1} \rangle , \psi_{T}^{2} \rangle , \dots , | \psi_{T}^{i} \rangle , . . .$ . from a base state $| \psi_{T}^{0} \rangle$ , where the parameters of each $| \psi_{T}^{i} \rangle$ differ only slightly from those of $| \psi_{T}^{0} \rangle$ . If we were to compute the energy of each state by drawing independent configurations, the statistical errors of the differences between any pair of estimates would add and their confidence intervals could overlap. By drawing configurations from the base trial state and computing all energies with this common set of configurations, the statistical error of the difference between any pair is less than the sum because the pair’s estimate is correlated.

Appendix O The Bryan algorithm#

There are multiple ways to maximize a nonlinear function of N variables. Here we describe the Levenberg-Marquardt method (Golub and Loan, 1989; Press et al., 2007), which is a variation of Newton’s method. After its description, we discuss how to adjust it for the specific problem of maximizing Q as stated in (12.39).

Our problem here is to maximize some general function

with respect to A. Starting from some initial guess for A, we iteratively improve it until the procedure converges. Suppose our current guess is A and we want to find a better solution A-. We start by considering the Taylor series expansion of Q about some point A,

where the vector δA is the deviation $\delta A = A^{\prime} - A$

are the components of the gradient of Q evaluated at A, and

are the elements of the Hessian matrix of $Q$ evaluated at A. Near the maximum, the Hessian is positive definite.

We can regard (O.1) as defining a local approximation for the gradient

Because we want to make the gradient at $A^{\prime}$ vanish, within the second-order approximation (O.1), we should make δA satisfy

so the relation between $A^{\prime}$ and A is

If the function $Q$ were in fact quadratic with a positive-definite Hessian, then the above analysis shows that we could move from $A$ to the maximum $\hat{\boldsymbol A}$ in a single step, that is, $\hat{A} = A + H^{- 1} g$ . When (O.1) is not exact, we can execute this stepping procedure iteratively,

until convergence (as measured, for example, by $0 < \sqrt{g^{T} \cdot g} \ll 1 )$ . The method will converge provided we can keep stepping uphill. Adding the matrix $\mu I$ to the Hessian shifts its eigenvalues, and adjusting $\mu$ at each step to maintain the local positive definiteness of $\mu I + H$ keeps us stepping uphill. This is the Levenberg-Marquardt strategy, summarized in Algorithm 51.

In Algorithm 51, what do we mean by solving for δA? Clearly, $\left[ \mu I + H \right] \cdot \delta A = g$ is a simultaneous system of linear equations for δA. For a given value of $\mu$ , we use the standard Gaussian elimination method with partial pivoting to solve them for δA, and check whether the new position $A + \delta A$ moves the new value of $Q$ above the old. If not, we could then repeat the process with new values of $\mu$ until it does. When we finally step uphill, we accept the new position and continue iterating.

Algorithm 51 Levenberg-Marquardt algorithm.
Input: Starting point A and a convergence parameter .
Set ν = 10, μ- = 0.01, A- = A, and Q- = Q(A) ;
repeat
$\mu \mu^{\prime} / \nu$ ;
repeat
Solve [μI + H(A-)] · δA = g(A-) for δA ;
if Q(A + δA) < Q- then
μ ← μν ;
end if
until Q(A + δA) > Q-.
μ- ← μ ;
A- ← A + δA ;
Q- ← Q(A + δA) ;
until $\left\| g ( A^{\prime} ) \right\| _{2} < \epsilon .$
return $A^{\prime} .$

Because H is symmetric, we alternatively could just as easily diagonalize it, writing the diagonalization as $H = Z \cdot \Xi \cdot Z^{T}$ , where $Z$ is an orthonormal matrix and $\Xi$ is a diagonal matrix whose diagonal elements are the eigenvalues $\xi_{i}$ of H. We now transform the system of linear equations

and thereby transform the coupled system of linear equations into a set of independent linear equations. Ensuring we step uphill becomes equivalent to adjusting μ so that $\mu + \xi_{i} > 0$ for all i. Because an orthonormal transformation preserves norms, our stopping criterion

remains unaltered by the coordinate transformation.

For our specific problem discussed in Section 12.3, the gradient is

with

where $\delta G$ is a vector whose components are $G_{i} - \bar{G}_{i}$ . The gradient of S follows from (12.43).

When $\nabla Q = 0$

The expression $C^{- 1}$ · δG is just a vector. What the above says is that the solution is a linear combination of the vectors representing the columns of $K^{T}$ . With the singular value decomposition $K^{T} = U \cdot D \cdot V^{T}$ , these vectors are actually the orthogonal columns of $U .$ Based on this observation, Bryan (1990) suggested parameterizing the solution as

The condition for the maximum then becomes

Hence,

is the gradient in the space defined by the parameters u. The Hessian in this space is

We now proceed to the maximization procedure, invoking $H \delta u = g$ for the Newton step using the just-defined expressions for g and H.

At each step of the iteration, we must compute g and H for the current estimate of u. We now recall, however, that before the start of the maximization procedure, we transformed the Green’s function and the kernel into a space where the covariance matrix C is diagonal, so $C^{- 1}$ is some diagonal matrix  (12.37). We also have $K = V \cdot D \cdot U^{T}$ . Accordingly,

We note that

where is a diagonal matrix whose diagonal elements are $A_{i}$ . Finally, we obtain

When we restrict the search space to the s-dimensional sub-dominant space of the kernel K, each of these matrices and their product becomes a symmetric $s \times s$ matrix. The solution to this reduced system of linear equations becomes a numerical task that we could handle by Gaussian elimination with partial pivoting. We could also proceed as before by diagonalizing the symmetric matrix $X \cdot Y$ . However, another, more preferable strategy is available.

The preferred strategy transforms H into a form similar to one we would get with the diagonalization. First, we note that X depends on the matrices we obtain from singular-value decomposition of the kernel K and the diagonalization of the covariance C. Because these matrix factorizations are done prior to starting the maximization procedure, they are done only once, and accordingly X is invariant during the maximization procedure. On the other hand, Y changes at each step as the values of $u_{i}$ and hence $A_{i}$ change. We also note that the matrix U, completing the definition of Y, does not change.

To transform H, we first diagonalize Y and write $\boldsymbol{Y} = \boldsymbol{W} \cdot \boldsymbol{\Omega} \cdot \boldsymbol{W}^{T}$ , where W is an orthonormal matrix and  is a diagonal matrix whose diagonal elements are the eigenvalues of Y. Next, we define the symmetric matrix $P = \sqrt{\Omega} \cdot W \cdot X \cdot W^{T} \cdot \sqrt{\Omega}$ and then diagonalize it, writing ${\cal P} = {\cal Z} \cdot \Xi \cdot{\cal Z}^{T}$ . If we define $F = W \cdot \sqrt{\Omega^{- 1}} \cdot Z ,$ then direct substitution shows that $F^{- T} \cdot F^{- 1} = Y$ . Now, we perform the series of transformations

and again reduce the simultaneous system of linear equations to a set of independent linear equations but in this case for the step $\delta u^{\prime} .$ . The Levenberg-Marquardt phase of the procedure chooses $\mu$ such that $\alpha + \xi_{i} + \mu > 0$ for all i.

The transformation just made is not norm preserving. The stopping condition $\sqrt{[ \delta u^{\prime} ]^{T} \cdot \delta u^{\prime}} < \epsilon$ is, however, equivalent to $\sqrt{\delta u^{T} \cdot Y \cdot \delta u} < \epsilon$ . Before parametrizing A by u and being mindful of the metric for flat curvature of the entropy, we would have used the stopping criterion $\sqrt{\sum_{i} ( \delta A_{i} )^{2} / A_{i}} ~ < ~ \epsilon$ . Remarkably, $\begin{array} {r} {\sum_{i} ( \delta A_{i} )^{2} / A_{i} = \sum_{i j} \delta u_{i} Y_{i j} \delta u_{j}} \end{array}$ ; that is, Y is the metric in u-space.

The readily available LAPACK and BLAS software facilitates the numerical implementation of this algorithm, which is due to Bryan (1990).