text:
# Conditional Expectations and Renormalization

## 1 Introduction

In the optimal prediction (OP) methods presented in earlier work by the author and others ,,,, an estimate of the future solution of an underresolved problem, or of a problem where the initial data are only partially known, was obtained by solving a reduced system of equations for the conditional expectation of the solution given the partial data. This system, closely related to a generalized Langevin equation of the Mori-Zwanzig type ,, is derived in detail in . Short-time estimates can be obtained by keeping only the first term on the right hand side of this system, obtaining a relation between the rate of change of a reduced set of variables and conditional expectations of the full rate of change; a simplified derivation of this relation is given below. Hald’s theorem  asserts that if one starts with a Hamiltonian system, then the reduced system obtained in this way is also Hamiltonian, with a Hamiltonian equal to a conditional free energy of the original Hamiltonian system.

Renormalization group (RNG) transformations ,, reduce the dimensionality of a system of equations by integrating out unwanted scales. That there is a qualitative resemblance between OP and RNG methods is quite clear, and has been pointed out in particular in the related work of Goldenfeld et al. ,. In the present paper we focus on the special case of Hamiltonian systems in thermal equilibrium, and show that in this case the RNG transformations of the Hamiltonian can be obtained by integrating conditional expectations of the derivatives of the Hamiltonian; loosely speaking, RNG transformations are integrals of OP reductions. This remark, based on Hald’s theorem, makes possible the efficient evaluation of the coefficients of the new Hamiltonians in RNG transformations (the “RNG parameter flow”) by simple Monte-Carlo methods, for example by Swendsen’s small-cell Monte-Carlo RNG ,. The coefficients in the new Hamiltonian define the reduced system of equations used to estimate the future in OP. To illustrate the construction, we apply it to spin systems and obtain explicitly the parameter flows in addition to critical points, critical exponents, and order parameters. We exhibit in detail a particular implementation that is a little awkward if viewed as an instance of a RNG but is particularly convenient for OP.

A little thought shows that what is offered in the present paper is a numerical short-cut. Suppose $`x=(x_1,x_2,\mathrm{})`$ is a set of $`n`$ random variables ($`n`$ may be infinite), and let $`m<n`$; partition $`x`$ so that $`x=(\widehat{x},\stackrel{~}{x})`$, $`\widehat{x}=(x_1,x_2,\mathrm{},x_m),\stackrel{~}{x}=(x_{m+1},x_{m+2},\mathrm{})`$. Let $`p=p(x)`$ be the joint probability density of all the variables, and consider the problem of finding a function $`\widehat{H}=\widehat{H}(\widehat{x})`$ such that

$$\mathrm{exp}(\widehat{H}(\widehat{x}))=p(x)𝑑\stackrel{~}{x},$$
$`(1)`$
where $`d\stackrel{~}{x}=dx_{m+1}dx_{m+2}\mathrm{}`$. There is no question that $`\widehat{H}`$ is well defined but the obvious ways of finding it can be costly. We are offering effective ways to do so. There are other situations where one wants to integrate out unwanted variables inside nonlinear functions, and our short-cut may serve there as well; in subsequent papers we shall apply it to problems in irreversible statistical mechanics and, equivalently, to problems involving the full long-time OP equations.

## 2 Conditional expectations and optimal prediction

Consider a set $`x`$ of random variables $`(x_1,x_2,\mathrm{},x_n)`$ with a joint probability density of the form $`Z^1e^{H(x)}`$, $`Z=e^{H(x)}𝑑x`$, $`dx=dx_1dx_2\mathrm{}dx_n`$. Consider the space $`L_2`$ of function $`u(x),v(x),\mathrm{}`$, with the inner product $`u,v=E[uv]=u(x)v(x)Z^1\mathrm{exp}(H)𝑑x`$, where $`E[]`$ denotes an expected value.

Partition the variables into two groups as above, $`x=(\widehat{x},\stackrel{~}{x})`$, $`\widehat{x}=(x_1,\mathrm{},x_m)`$, $`m<n`$. Given a function $`f(x)`$, its conditional expectation given $`\widehat{x}`$ is

$$E[f(x)|\widehat{x}]=\frac{{\displaystyle f(x)e^{H(x)}𝑑\stackrel{~}{x}}}{{\displaystyle e^{H(x)}𝑑\stackrel{~}{x}}};$$
$`(2)`$
it is the average of $`f`$ keeping $`\widehat{x}`$ fixed. The conditional expectation is a function of $`\widehat{x}`$ only, and it is the best approximation of $`f`$ in the mean square sense by a function of $`\widehat{x}`$:

$$E\left[\left(f(x)E[f(x)|\widehat{x}]\right)^2\right]E\left[\left(f(x)h(\widehat{x})\right)^2\right]$$
$`(3)`$
for any function $`h=h(\widehat{x})`$. $`E[f|\widehat{x}]`$ is the orthogonal projection of $`f`$ onto the subspace $`\widehat{L}_2`$ of $`L`$ that contains functions of $`\widehat{x}`$ only. $`E[f(x)|\widehat{x}]`$ can be approximated by expansion in a basis of $`\widehat{L}_2`$; keeping only a suitable finite number $`\mathrm{}`$ of basis functions $`\phi _1(\widehat{x}),\phi _2(\widehat{x}),\mathrm{},\phi _{\mathrm{}}(\widehat{x})`$, and minimizing the distance between $`f`$ and the span of the $`\phi _i(\widehat{x})`$, one finds

$$E[f|\widehat{x}]=\underset{i=1}{\overset{\mathrm{}}{}}c_i\phi _i(\widehat{x}),$$

where the $`c_i`$ satisfy the equation

$$\mathrm{\Phi }c=r,$$
$`(4)`$
where $`\mathrm{\Phi }`$ is the matrix with elements $`\mathrm{\Phi }_{ij}=\phi _i,\phi _j`$, $`c=(c_1,\mathrm{},c_{\mathrm{}})`$, and $`r=(f,\phi _1,f,\phi _2\mathrm{},f,\phi _{\mathrm{}})`$. Usually the inner products can be calculated by Metropolis sampling.

Suppose you want to find a function $`\widehat{H}=\widehat{H}(\widehat{x})`$ such that

$$e^{\widehat{H}(\widehat{x})}=e^{H(\widehat{x},\stackrel{~}{x})}𝑑\stackrel{~}{x},$$
$`(5)`$
i.e., write the marginal probability density of the variables $`\widehat{x}`$ in exponential form. Suppose one can write

$$H(x)=\underset{i=1}{\overset{\mathrm{}}{}}\alpha _i\phi _i(x),$$

and and let $`im`$, where $`m`$ is the number of components of the vector $`\widehat{x}=`$ Then

$$\begin{array}{ccc}E\left[\frac{}{x_i}H(x)|\widehat{x}\right]& =& \frac{{\displaystyle \frac{}{x_i}H\left(x\right)e^{H(x)}𝑑\stackrel{~}{x}}}{{\displaystyle e^{H(x)}𝑑\stackrel{~}{x}}}\\ & & \\ & =& \frac{}{x_i}\left(\mathrm{log}e^{H(x)}𝑑\stackrel{~}{x}\right).\end{array}$$
$`(6)`$
An analogous relation between the derivative of a logarithm of a partially integrated density and a conditional expectation arises also in the context of expectation-maximization is statistics .

If one can find a basis for $`\widehat{L}_2`$ consisting of functions of the form $`\frac{}{x_1}\phi _j(\widehat{x})`$, $`j=1,\mathrm{},`$ and provided the set of variables is homogenous so that for all $`im`$ the coefficients $`c_j`$ in the expansions $`\left[\frac{}{x_i}H|\widehat{x}\right]=_{j=1}^{\mathrm{}}c_j\frac{}{x_i}\phi _j(\widehat{x})`$ are independent of $`i`$, then the expansion

$$\widehat{H}(\widehat{x})=c_j\phi _j(\widehat{x}).$$
$`(7)`$
follows immediately. This is our key observation.

This construction is just Hald’s theorem for OP : Suppose one has a system of differential equations (written as ordinary differential equations for simplicity) of the form

$$\frac{d}{dt}\phi (t)=R\left(\phi (t)\right),\phi (0)=x$$
$`(8)`$
where $`\phi ,R`$ and $`x`$ are $`n`$-vectors with components $`\phi _i,R_i,x_i,`$ $`i=1,\mathrm{},n`$ and $`t`$ is the time. Suppose we partition as above $`\phi =(\widehat{\phi },\stackrel{~}{\phi })`$, $`R=(\widehat{R},\stackrel{~}{R})`$, where $`\widehat{\phi }`$ contains the first $`m`$ components of $`\phi `$, etc. Suppose the system (8) is Hamiltonian, i.e., $`m,n`$ are even, $`R_i=\frac{}{x_{i1}}H`$ for $`i`$ even, $`R_i=\frac{}{x_{i+1}}H`$ for $`i`$ odd; $`H=H(x)`$ is the Hamiltonian and $`Z^1e^H`$ is then an invariant probability density for the system.

Suppose we can afford to solve only $`m<n`$ of the equations in (8) or have only $`m`$ data components $`\widehat{x}`$. We want to solve equations for $`\widehat{\phi }`$:

$$\frac{d\widehat{\phi }}{dt}=\widehat{R}(\phi ),\widehat{\phi }(0)=\widehat{x},$$

where $`im`$, but the argument of $`\widehat{R}`$ is the whole vector $`\phi `$. It is natural to approximate $`\widehat{R}_i(\phi )`$ by the closest function of $`\widehat{\phi }`$ for each $`im`$, i.e., solve

$$\frac{d\widehat{\phi }}{dt}=E[\widehat{R}(\phi )\widehat{\phi }].$$
$`(9)`$
The approximation (9) is valid only for a short time, as one can see from the full equation for the evolution of $`\frac{d\widehat{\phi }}{dt}`$ in ,, see also below. Hald’s theorem asserts that the system (9) is Hamiltonian, with Hamiltonian $`\widehat{H}=\widehat{H}(\widehat{x})=\mathrm{log}e^{H(x)}𝑑\stackrel{~}{x}`$, a relation identical to equation (6). The existence of $`\widehat{H}`$ shows that the approximation (9) cannot be valid for long times: the predictive power of partial initial data decays at $`t\mathrm{}`$ for a nonlinear system, and the best estimate of $`\widehat{\phi }(t)`$ should decay to unconditional mean of $`\phi `$ (which is usually zero). The existence of a reduced Hamiltonian shows that this decay can happen only to a limited extent and thus the approximation can in general be valid only for short times. Equations (9) constitute the short time, or “first-order”, OP approximation.

Suppose however that instead of picking specific values for the initial data $`\widehat{x}`$ one samples them from the invariant density $`Z^1e^{\widehat{H}(\widehat{x})}`$. The distribution of the $`\widehat{x}`$’s is then invariant, and equal to their marginal distribution in the full system (8) when the data are sampled from the invariant distribution $`Z^1e^{H(x)}`$, as one can also see from the identities $`\mathrm{exp}(\widehat{H})=\mathrm{exp}\left(\mathrm{log}e^H𝑑\stackrel{~}{x}\right)=e^H𝑑\stackrel{~}{x}`$. The system (9) then generates the marginal probability density of part of the variables of a system at equilibrium. Thus OP at equilibrium is a way of reducing the number of variables without affecting the statistics of the variables that remain. One can make short-time predictions about the future from the reduced system with coefficients computed at equilibrium because it is self-consistent to assume for short times that unresolved degrees of freedom are in thermal equilibrium, as is explained in the OP papers cited above.

## 3 Renormalization

For simplicity, we work here with real-space renormalization applied to variables associated with specific sites in a plane, $`x^{(1)}=(x_{I_1},x_{I_2},\mathrm{})`$, where $`I_k=(i_k,j_k)`$, $`i_k,j_k`$ are integers, all the $`I_k`$ are inside a square $`D`$ of side $`N`$ with $`N`$ large, and the variables are listed in some convenient order. The Hamiltonian $`H=H^{(1)}`$ is a function of $`x^{(1)}`$, $`H^{(1)}=H^{(1)}(x^{(1)})`$. The need for the superscript (1) will appear shortly. We assume that the partition function $`Z=e^{H^{(1)}(x^{(1)})}𝑑x^{(1)}`$ is well defined, where $`dx^{(1)}=dx_{I_1}^{(1)}dx_{I_2}^{(1)}\mathrm{}`$.

Suppose we group the variables $`x_{I_1},x_{I_2},\mathrm{}`$ into groups of $`\mathrm{}`$ variables (for example, we could divide $`D`$ into squares each containing 4 variables). The variables can be referred to as “spins”in conformity with common usage in physics. Associate with each group a new variable $`x_{J_1}^{(2)},x_{J_2}^{(2)},\mathrm{},`$ where $`J_1,J_2,\mathrm{}`$ is some ordering of the new variables and $`x_{J_k}^{(2)}`$ is a function (not necessarily invertible) of the $`x_I^{(2)}`$ in the group labeled by $`J_k`$, for example $`x_{J_k}^{(2)}=g(x_{I_{m+1}}^{(1)},x_{I_{m+2}}^{(1)},\mathrm{},x_{I_{m+\mathrm{}}}^{(1)})`$ for the appropriate $`m`$. The vector $`x^{(2)}`$ is $`x^{(2)}=(x_{J_1}^{(2)},x_{J_2}^{(2)}\mathrm{})`$. We can write

$$\begin{array}{ccc}Z& =& e^{H^{(1)}(x^{(1)})}𝑑x^{(1)}\\ & & \\ & =& 𝑑x^{(2)}\delta \left(x^{(2)}g(x^{(1)})\right)e^{H^{(1)}(x^{(1)})}𝑑x^{(1)}\end{array}.$$

where $`dx^{(2)}=dx_{J_1}^{(2)}dx_{J_2}^{(2)}\mathrm{},`$ and the $`\delta `$ function is a product of delta functions, one per group. If one defines $`H^{(2)}(x^{(2)})`$ by the equation

$$e^{H^{(2)}(x^{(2)})}=\delta \left(x^{(2)}g(x^{(1)})\right)e^{H(x^{(1)})}𝑑x^{(1)},$$
$`(10)`$
then $`Z=e^{H^{(2)}(x^{(2)})}𝑑x^{(2)}`$.

The mapping $`x^{(1)}x^{(2)}`$, followed by a change of numbering of the remaining variables so that $`J_1,J_2\mathrm{}`$ (the indices of the new variables $`x^{(2)}`$) enumerate the new variables by going through all integer pairs in a reduced domain of side $`N/\sqrt{\mathrm{}}`$, is a real-space renormalization group transformation; it produces a new set of variables which has less spatial detail than the previous set and such that distances between the remaining spins have been scaled down by $`\sqrt{\mathrm{}}`$. If the calculation is set up so that the mapping $`x^{(1)}x^{(2)}`$, $`H^{(1)}H^{(2)}`$ can be repeated, for example, if the range of the variables $`x^{(1)}`$ is invariant and the Hamiltonians $`H^{(1)},H^{(2)}`$ can be represented in the same finite-dimensional basis, then one can produce in this way a sequence of Hamiltonians $`H^{(1)},H^{(2)},H^{(3)},\mathrm{}`$; the fixed points of the transformation $`H^{(n)}H^{(n+1)}`$ for a spin system of infinite spatial extent include the critical points of the system, see any discussion of the RNG, for example ,.

Consider the special case where $`x_J^{(2)}`$ is one of the $`x_I^{(1)}`$ in its group–i.e., replace a block of spins by one of the spins in the block. More general and widely used assignments of block variables will not be needed in the present paper and will be discussed elsewhere. We can identify the spins that remain with $`\widehat{x}`$ of the preceding section and the spins that disappear with $`\stackrel{~}{x}`$. Equation (10) becomes a special case of equation (5), and can be solved for $`H^{(2)}`$ by taking conditional expectations of the derivatives of $`H^{(1)}`$.

Note that the usual RNG representation of a renormalized Hamiltonian by means of additional couplings  is interpreted here as an expansion of a conditional expectation in a convergent series. The new interpretation may be useful both in understanding what is happening on the computer and in deriving error estimates. The relation between the RNG and conditional expectations shows that the latter can be calculated recursively, as we show in the example below.

We have written the RNG transformation above in notation suitable for spins with a continuous range. The case of discrete (e.g., Ising) spins is automatically included, even though it may seem odd to differentiate functions with a discrete domain and range. Indeed, add to the Hamiltonian $`H`$ a term of the form

$$\frac{1}{\epsilon }\underset{i}{}\left(\underset{j}{}\psi (x_ix_{0j})\right)$$

where $`\epsilon `$ is small, the sum is over all spins, the product is over a finite number of values $`x_{0j}`$, and $`\psi 0`$ has a minimum at 0 and is positive elsewhere. For small $`\epsilon `$ such a term will constrain the $`x_i`$ to take on the values $`x_{0j}`$, but since at the origin the derivative of $`\psi `$ is zero the calculation of the conditional expectations is unaffected by this term and the limit $`\epsilon 0`$ can be taken without actually doing anything on the computer. All one has to do is make sure that in the Monte-Carlo sampling only the values $`x_{0j}`$ are sampled. Indeed, results below will be given for Ising spins which take on the values $`+1`$ and $`1`$, with a “bare” (unrenormalized) Hamiltonian $`H^{(1)}=\beta x_Ix_J`$, with summation over locations $`I,J`$ that are neighbors on the lattice; $`\beta =1/T`$, where $`T`$ is the temperature.

## 4 A decimation RNG/OP scheme for a spin system

We consider in detail a RNG/recursive OP scheme where the number of variables is halved at each step. The spins are located on a square lattice with nodes $`I_k=(i_k,j_k)`$, $`i_k,j_k`$ integers, and at each step of the recursion those for which $`i_k+j_k`$ is odd are eliminated while those for which $`i_k+j_k`$ is even are kept. The spins with $`i_k+j_k`$ even constitute $`\widehat{x}`$ and the others $`\stackrel{~}{x}`$; the choice of which are even and which are odd is a matter of convention only (see Figure 1). The variables are labeled by $`I_k:x_{I_1},x_{I_2},\mathrm{}`$.

Given a location $`I=(i,j)`$, we group the other variables according to their distance from $`I`$: group 1 contains only $`x_I`$, the variable at $`I`$. Group 2 (relative to $`I`$) contains those variables whose distance from $`I`$ is 1, group 3 contains those variables whose distance to $`I`$ is $`\sqrt{2}`$, etc. We form the “collective” variables

$$X_{k,I}=\frac{1}{n_k}\underset{\mathrm{group}k}{}x_J$$

where $`n_k`$ is the number of variables in the group (1 for group 1, 4 for group 2, etc.). From these variables one can form a variety of translation-invariant polynomials in $`x`$ of various degrees: $`_Ix_IX_{k,I}=_IX_{1,I}X_{k,I},`$ $`_I(X_{k,I})^2(X_{k+1,I})^2,`$ $`_I(X_{k,I})^4,\mathrm{}`$. In practice the domain over which one sums must be finite, and it is natural to impose periodic boundary conditions at its edges to preserve the translation invariance. We wrote out explicitly only polynomials of even degrees because the Hamiltonians we consider are invariant under the transformation $`xx`$. The translation-invariant polynomials built up from the $`X_{k,I}`$ can be labeled $`\phi _1(x),\phi _2(x),\mathrm{}`$ in some order.

Expand the $`n`$-th renormalized Hamiltonian in a series and keep the first $`\mathrm{}`$ terms:

$$H^{(n)}=\underset{k=1}{\overset{\mathrm{}}{}}\alpha _k^{(n)}\phi _k(x).$$
$`(11)`$
The derivative of this series at the spin $`x_I`$ is

$$\frac{}{x_I}H^{(n)}=\underset{1}{\overset{\mathrm{}}{}}\alpha _k^{(n)}\phi _k^{}(x),\phi _k^{}=\frac{}{x_I}\phi _k.$$
$`(12)`$
The functions $`\phi _k^{}`$ are easily evaluated, for example:

$$\left(\underset{J}{}x_JX_{k,J}\right)^{}=2X_{k,I}$$

$$\left(\underset{J}{}(X_{k,J})^4\right)^{}=4\left(\underset{\mathrm{group}k}{}x_J^3\right)/n_k^2,$$

etc., where “group $`k`$” refers to distances from $`I`$, the variable with respect to which we are differentiating (see Figure 2).

Pick a variable $`x_I`$ in $`\widehat{x}`$ (i.e., $`I=(i,j),`$ $`i+j`$ even in our conventions). Some of the functions $`\phi _k^{}`$ in (12) will be functions of $`\widehat{x}`$ only and some will be functions of both $`\widehat{x}`$ and $`\stackrel{~}{x}`$ or of $`\stackrel{~}{x}`$ only. The task at hand is to project the latter on the former and then rearrange the series so as to shrink the scale of the physical domain. To explain the construction we consider a very special case.

Suppose one can write

$$H^{(n)}(x)=\alpha _2^{(n)}\phi _2+\alpha _3^{(n)}\phi _3+\alpha _4^{(n)}\phi _4,$$
$`(13)`$
where $`\phi _k(x)=x_IX_{k,I}`$ for $`k=2,3,4,\mathrm{}`$. Note that $`\phi _k^{}=\frac{}{x_I}\phi _k`$ is a function only of $`\widehat{x}`$ when $`k=3,4`$ (and when $`k=6`$, as we shall need to know shortly) but not when $`k=2`$ or 5 (see Figures 1, 2). We now calculate the conditional expectations of the derivatives of $`H^{(n)}`$ given $`\widehat{x}`$ by projecting them on the space of functions of $`\widehat{x}`$. First we project $`\phi _2^{}`$ on the span of $`\phi _3^{},\phi _4^{},\phi _6^{}`$ (note that $`\phi _6^{}`$ is not in the original expansion (13)). Form the matrix $`\mathrm{\Phi }`$ with rows $`\phi _k^{},\phi _3^{},\phi _k^{},\phi _4^{},\phi _k^{},\phi _6^{}`$ for $`k=3,4,6`$, the primes once again denoting differentiation with respect to $`x_I`$. Form the vector $`r`$ with component $`(\phi _2^{},\phi _3^{},\phi _2^{},\phi _4^{},\phi _2^{},\phi _6^{})`$. Let $`c=(c_1,c_2,c_3)`$ be the solution of $`\mathrm{\Phi }c=r`$ (see equation (4)). The coefficients $`c`$ are the coefficients of the orthogonal projection of $`\phi _2^{}`$ onto the span of $`\phi _3^{},\phi _4^{},\phi _6^{}`$ which is contained in $`\widehat{L}_2`$. After projection, the coefficients of $`\phi _3,\phi _4`$ in (13) become

$$\alpha _3^{\mathrm{new}}=\alpha _3^{(n)}+\alpha _2^{(n)}c_1,$$

$$\alpha _4^{\mathrm{new}}=\alpha _4^{(n)}+\alpha _2^{(n)}c_2,$$

and $`\phi _6`$ acquires the coefficient $`\alpha _2^{(n)}c_3`$.

If one relabels the remaining spins so that they occupy the lattice previously occupied by all the spins, group 3 becomes group 2, group 4 becomes group 3, and group 6 becomes group 4 (see Figure 2). The new Hamiltonian $`H^{(n+1)}`$ now has the representation

$$H^{(n+1)}=\alpha _2^{(n+1)}\phi _2+\alpha _3^{(n+1)}\phi _3+\alpha _4^{(n+1)}\phi _4,$$

with

$$\begin{array}{ccc}\alpha _2^{(n+1)}& =& \alpha _3^{(n)}+\alpha _2^{(n)}c_1,\\ & & \\ \alpha _3^{(n+1)}& =& \alpha _4^{(n)}+\alpha _2^{(n)}c_2,\\ & & \\ \alpha _4^{(n+1)}& =& \alpha _2^{(n)}c_3.\end{array}$$

More generally, if $`H^{(n)}`$ is expressed as a truncated series, partition the terms in the series for $`\frac{}{x_I}H^{(n)}`$ into functions of $`\widehat{x}`$ and functions of both $`\widehat{x}`$ and $`\stackrel{~}{x}`$. Add to the terms which are functions of $`\widehat{x}`$ additional terms which are functions of $`\widehat{x}`$ and are chosen so that after relabelling they acquire the form of terms already in the series (just as above, terms that depend on $`X_{3,J}`$, for example, become terms that depend on $`X_{2,J}`$ after relabelling). Project the functions of $`x`$ on the span of the expanded set of functions of $`\widehat{x}`$, collect terms and relabel. This is a renormalization step, and it can be repeated.

Note that it if one wants to reduce the number of variables by a given factor, one can in principle use an analogous RNG/conditional expectation construction and get there in one iteration; the recursive construction is easier to do and the intermediate Hamiltonians, whose coefficients constitute the parameter flow in the renormalization, contain useful information.

The discussion so far may suggest that one sample the Hamiltonians recursively, i.e., start with $`H^{(1)}`$, find $`H^{(2)}`$, use Monte-Carlo to sample the density $`Z^1e^{H^{(2)}}`$ and find $`H^{(3)}`$ etc. The disadvantages of this approach are: (i) The sampling of the densities $`Z^1e^{H^{(n)}}`$ can be much more expensive for $`n>1`$ than for $`n=1`$ because each proposed Monte-Carlo move may require that the full series for $`\frac{}{x_I}H^{(n)}`$ be summed twice; and (ii) each evaluation of a new Hamiltonian is only approximate because the series are truncated, and, more important, the Monte-Carlo evaluation of the coefficients may have limited accuracy. These errors accumulate from step to step and may produce false fixed points and other artifacts.

The remedy lies in Swendsen’s observation , that the successive Hamiltonians can be sampled without being known explicitly. Sample the original Hamiltonian, remove the unwanted spins and relabel the remaining spins so as to cover the original lattice, as in the relabelling step in the renormalization; the probability density of the remaining spins is $`Z^1e^{H^{(2)}}`$; repeating $`n`$ times yields samples of $`Z^1e^{H^{(n+1)}}`$. The price one pays is that to get an $`m`$ by $`m`$ sample of $`Z^1e^{H^{(n)}}`$ one has to start by sampling a $`2^qm`$ by $`2^qm`$ array of non-renormalized spins, where $`q`$ is either $`(n+1)/2`$ or $`n/2`$ depending on the parity of $`n`$ and on programming choices; the trade-off is in general very worthwhile. What has been added to Swendsen’s calculation is an effective evaluation of the coefficients of the expansion of $`H^{(n)}`$ from the samples.

The programming here requires some care. With the decimation scheme as in Figure 1, after one removes the unwanted spins in $`x^{(n)}`$ the remaining spins, the variables $`x^{(n+1)}`$, live on a lattice with a mesh size $`\sqrt{2}`$ larger than before; after relabelling they find themselves on a lattice with the same mesh size as before but arranged at a $`\pi /4`$ angle with respect to the previous lattice. To extract a square array from at this set of spins one has to make the size of the box that includes all the spins half the size of the previous box. At the next renormalization one obtains $`x^{(n+2)}`$ which can be extracted from $`x^{(n)}`$ by taking one spin in four and the resulting box size is the same as the size of the box that contains $`x^{(n+1)}`$. One may worry a little about boundary conditions for $`x^{(n+1)}`$: the periodicity of $`x^{(n)}`$ is not the same as the periodicity one has to assume for $`x^{(n+1)}`$ because of the rotation; the resulting error is too small to be detected in our calculations.

## 5 Some numerical results

We now present some numerical results obtained with the RNG/conditional expectation scheme. The problem we apply the construction to is Ising spins; more interesting applications will be presented elsewhere. The point being made is that the construction can be effectively implemented. The results are presented for Ising spins.

In table I we list the coefficients $`\alpha _k^{(n)}`$ in the expansion of $`H^{(n)}`$ for $`n=1,\mathrm{}7`$ and $`T=2.27`$. The functions $`\phi _k`$ are as follows:

$$\phi _k=x_JX_{k,J}\mathrm{for}k=1,2,3,4,5,6$$

$$\phi _{6+k}=(X_{k+1,J})^4,\mathrm{for}k=1,2,3$$

$$\phi _{10}=X_{2,J}^2X_{3,J}^2.$$

Note that as a result of the numbering of the $`\phi `$’s the last coefficient is not necessarily the smallest coefficient. This table represents the parameter flow and if the functions $`\phi _k`$ are written in terms of the variables $`x_J`$ the table defines the new system of equations for the reduced set of variables. Remember that in the projection on $`\widehat{L}_2`$ additional functions are used so that after relabelling the series has the same terms , but maybe with different coefficients, as before the renormalization. In $`H^{(1)}`$, $`\alpha _2`$ is the sole non-zero coefficient, and its value is determined by $`T`$ and the definition of $`X_{2,J}`$, in particular the presence of the coefficient $`n_2`$ (see above).

It is instructive to use the parameter flow to identify the critical temperature $`T_c`$. For $`T<T_c`$ the renormalization couples ever more distant spins while for $`T>T_c`$ the spins become increasingly decoupled. One can measure the increasing or decreasing coupling by considering the quadratic terms in the Hamiltonian (the terms of the form $`x_JX_{k,J}`$) and calculating the “second moments” $`M_2`$ of their coefficients $`\alpha _k^{(n)}`$:

$$M_2^{(n)}=\underset{k=2}{\overset{\mathrm{}}{}}d_k^2\alpha _k^{(n)}$$

where $`d_k`$ is the distance from $`J`$ of the spins in the group $`k`$ (see the definition of $`X_{k,J}`$), $`\alpha _k^{(n)}`$ is the coefficient of $`x_JX_{k,J}`$ in the expansion of $`H^{(n)}`$, and $`\mathrm{}`$ is the number of quadratic terms in this expansion. In Figure 3 we show the evolution of $`M_2^{(n)}`$ with $`n`$ for various values of $`T`$ (with $`\mathrm{}=5`$ and 7 functions over-all in the expansion, including non-quadratic functions).

In Figure 4 we show the evolution of $`M_2^{(n)}`$ near $`T_c=2.269\mathrm{}`$ with $`\mathrm{}=6`$ and 10 terms in the expansion. The non-uniform behavior of $`M_2`$ is not a surprise (it is related to the non-uniform convergence of critical exponents already observed by Swendsen). Each step in the renormalization used $`10^5`$ Monte-Carlo steps per spin. From these graphs one would conclude that $`T_c2.26`$, an error of .5%. The accuracy depends on the number of terms in the expansion and on the choice of terms; with only 6 terms (4 quadratic and 2 quartic), the error in the location of $`T_c`$ increases to about 3%. The point is not that this is a good way to find $`T_c`$ but that it is a check on the accuracy of the parameter flow. From the Table one can see that the system first approaches the neighborhood of a fixed point and then diverges from it, as one should expect in a discrete sequence of transformations.

| Table 1 | | | | | | | |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Parameter flow for the Ising model $`T=2.26`$, 10 basis functions | | | | | | | |
| iteration | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
| $`\alpha _1`$ | 0 | .26 | .35 | .44 | .48 | .52 | .54 |
| $`\alpha _2`$ | .893 | .47 | .47 | .35 | .30 | .25 | .21 |
| $`\alpha _3`$ | 0 | .32 | .20 | .23 | .21 | .20 | .18 |
| $`\alpha _4`$ | 0 | .04 | .08 | .11 | .12 | .13 | .13 |
| $`\alpha _5`$ | 0 | .07 | .11 | .13 | .13 | .12 | .12 |
| $`\alpha _6`$ | 0 | $``$.01 | .01 | .01 | .02 | .03 | .02 |
| $`\alpha _7`$ | 0 | $``$.08 | $``$.07 | $``$.10 | $``$.09 | $``$.09 | $``$.08 |
| $`\alpha _8`$ | 0 | .04 | .02 | .02 | .01 | .00 | $``$.10 |
| $`\alpha _9`$ | 0 | $``$.00 | $``$.01 | $``$.00 | $``$.00 | .00 | .00 |
| $`\alpha _{10}`$ | 0 | $``$.12 | $``$.17 | $``$.18 | $``$.18 | $``$.17 | $``$.16 |

We now use the renormalized system to calculate the magnetization $`m=E[x_I/n^2]`$. To get the correct non-zero $`m`$ for $`T<T_c`$ on a small lattice the symmetry must be broken, and we do this by imposing on all the arrays the boundary condition $`x_{boundary}=1`$ rather than the periodic boundary conditions used elsewhere in this paper. In Figure 5 we display $`m`$ computer with the bare (unrenormalized) Hamiltonian $`H^{(1)}`$ on 3 lattices: $`20`$ by $`20`$, $`40`$ by $`40`$, $`60`$ by $`60`$, as well as the results obtained on a $`20`$ by $`20`$ lattice by sampling the density defined by the renormalized Hamiltonian $`H^{(5)}`$ which corresponds in principle to an $`80`$ by $`80`$ bare calculation. We also display the exact Onsager values of $`m`$. The calculations focus on values of $`T`$ in the neighborhood of $`T_c`$ where the size of the lattice matters; one cannot expect the results to agree perfectly with the Onsager results on a finite lattice with periodic boundary conditions for any $`n`$; all one can expect is to have the values of the small renormalized calculation be consistent with results of a larger bare calculation. We observe that they do, up to the shift in $`T_c`$ already pointed out and due to the choice of basis functions.

The determination of the critical exponents for a spin model is independent of the determination of the coefficients in the expansion of $`H^{(n)}`$, and is mentioned here only because it does provide a sanity check on the constructions, in particular on the adequacy of the basis functions. For comparable earlier calculations, see in particular Swendsen’s chapter in . As is well known, if $`A`$ is the matrix of derivatives $`\alpha _i^{(n+1)}/\alpha _j^{(n)}`$ at $`T=T_c`$, those of its eigenvalues that are larger than 1 are the critical exponents of the spin system . The matrix $`A`$ can be found from the chain rule ,

$$\frac{}{\alpha _j^{(n)}}E\left[\phi _k\left(x^{(n+1)}\right)\right]=\underset{i}{}\frac{\alpha _i^{(n+1)}}{\alpha _j^{(n)}}\frac{E[\phi _k(x^{(n+1)})]}{\alpha _i^{(n+1)}}$$

and the sum is over all the coefficients that enter the expansion. The derivatives of the expectations are given by correlations as follows:

$$\frac{E[\phi _k(x^{(n+1)})]}{\alpha _j^{(n)}}=E\left[\phi _k(x^{(n+1)})\phi _j(x^{(n)})\right]E[\phi _k(x^{(n+1)})]E[\phi _j(x^{(n)})],$$

$$\frac{E[\phi _k(x^{(n+1)})]}{\alpha _i^{(n+1)}}=E\left[\phi _k(x^{(n+1)})\phi _i(x^{(n+1)})\right]E[\phi _k(x^{(n+1)})]E[\phi _i(x^{(n+1)})],$$

see . In most of the literature on real-space renormalization for Ising spins the variables $`x^{(n+1)}`$ are obtained from $`x^{(n)}`$ by “majority rule”, i.e., by assigning to the group that defines $`x^{(n+1)}`$ the value $`+1`$ if most of the members of the group are $`+1`$, the value $`1`$ if most of the members of the group are $`1`$, with ties resolved at random. For the decimation scheme described above our “pick one” rule ($`x^{(n+1)}`$ is one of the members of the group) is identical to the majority rule. There is an apparent difficulty in the decimation because at each recursion the number of terms in the summation that defines the basis functions is reduced by half while the square root of an integer is not in general an integer, so that one has to perform Swendsen sampling on rectangles so designed that the ratio of the areas of two successive rectangles is 1/2. This has not turned out to be harmful, and the value of $`\nu `$, the correlation exponent, was found to be $`1`$ (the exact value) $`\pm .01`$ with $`10^6`$ Monte-Carlo moves per spin, the error depending mainly on the number of Monte-Carlo moves which has to be very large, in line with previous experience . We also checked that in a renormalization scheme where a $`2\times 2`$ block of spins is replaced at each iteration by a single spin, the “majority rule” and our “pick one” rule for $`x^{(n+1)}`$ yield similar results. One needs fewer terms in the expansion of the Hamiltonian to get accurate values of the exponents than to get an accurate parameter flow, but a larger number of Monte-Carlo moves.

## 6 Conclusions

We have presented a simple relation between conditional expectations for systems at equilibrium on one hand and the RNG on the other, which makes it possible to find efficiently the coefficients in a reduced systems of equations for a subset of variables whose distribution as given by reduced system equals their marginal distribution in the original system. The numerical results above emphasized the neighborhood of the critical point in the simple example because this is where the variables are strongly coupled without separation of scales and a reduction in system size requires non-trivial tools. The next steps will be the application of these ideas to time-dependent problems and to finite-difference approximations of underresolved partial differential equations, along the lines suggested in ; this work will be presented elsewhere.

Acknowledgments. I would like to thank Prof. G.I. Barenblatt, Prof. N. Goldenfeld, Prof. O. Hald, Prof. R. Kupferman, Mr. K. Lin, and Mr. P. Stinis for very helpful discussions and comments. This work was supported in part by the Office of Science, Office of Advanced Scientific Computing Research, Mathematical, Information, and Computational Sciences Division, Applied Mathematical Sciences Subprogram, of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098 and in part by the National Science Foundation under grant number DMS89-19074.