Product Distribution Field Theory
David H. Wolpert
This paper presents a novel way to approximate a distribution governing a system of coupled particles with a product of independent distributions. The approach is an extension of mean field theory that allows the independent distributions to live in a different space from the system, and thereby capture statistical dependencies in that system. It also allows different Hamiltonians for each independent distribution, to facilitate Monte Carlo estimation of those distributions. The approach leads to a novel energy-minimization algorithm in which each coordinate Monte Carlo estimates an associated spectrum, and then independently sets its state by sampling a Boltzmann distribution across that spectrum. It can also be used for high-dimensional numerical integration, (constrained) combinatorial optimization, and adaptive distributed control. This approach also provides a simple, physics-based derivation of the powerful approximate energy-minimization algorithms semi-formally derived in~\cite{wowh00, wotu02c, wolp03a}. In addition it suggests many improvements to those algorithms, and motivates a new (bounded rationality) game theory equilibrium concept.
cond-mat0307630.tex
% ****** Start of file apssamp.tex ******
%
%   This file is part of the APS files in the REVTeX 4 distribution.
%   Version 4.0 of REVTeX, August 2001
%
%   Copyright (c) 2001 The American Physical Society.
%
%   See the REVTeX 4 README file for restrictions and more information.
%
% TeX'ing this file requires that you have AMS-LaTeX 2.0 installed
% as well as the rest of the prerequisites for REVTeX 4.0
%
% See the REVTeX 4 README file
% It also requires running BibTeX. The commands are as follows:
%
%  1)  latex apssamp.tex
%  2)  bibtex apssamp
%  3)  latex apssamp.tex
%  4)  latex apssamp.tex
%
\documentclass[twocolumn,showpacs,preprintnumbers,amsmath,amssymb]{revtex4}
%\documentclass[preprint,showpacs,preprintnumbers,amsmath,amssymb]{revtex4}


\usepackage{psfig} 
\usepackage{stmaryrd} % for the Cartesian product
%\usepackage{mathabx} % for the Cartesian product
\usepackage{graphicx} % standard LaTeX graphics tool

\begin{document}

\preprint{APS/123-QED}

\title{Product Distribution Field Theory}

\author{David H. Wolpert}
\affiliation{NASA Ames Research Center, Moffett Field, CA, 94035, USA\\
{\tt dhw@email.arc.nasa.gov}}

\begin{abstract} This paper presents a novel way to approximate a
distribution governing a system of coupled particles with a product of
independent distributions.  The approach is an extension of mean field
theory that allows the independent distributions to live in a
different space from the system, and thereby capture statistical
dependencies in that system. It also allows different Hamiltonians for
each independent distribution, to facilitate Monte Carlo estimation of
those distributions. The approach leads to a novel energy-minimization
algorithm in which each coordinate Monte Carlo estimates an associated
spectrum, and then independently sets its state by sampling a
Boltzmann distribution across that spectrum.  It can also be used for
high-dimensional numerical integration, (constrained) combinatorial
optimization, and adaptive distributed control.  This approach also
provides a simple, physics-based derivation of the powerful
approximate energy-minimization algorithms semi-formally derived
in~\cite{wowh00, wotu02c, wolp03a}. In addition it suggests many
improvements to those algorithms, and motivates a new (bounded
rationality) game theory equilibrium concept.  \end{abstract}


\pacs{89.20.Ff,89.75.-k , 89.75.Fb}% PACS 

%\keywords{Suggested keywords}%Use showkeys class option if keyword
                              %display desired



\maketitle              % typesets the title of the contribution



\nocite{chjo02, kono03_short}


\section{Introduction} \label{sec:intro} Mean Field Approximations
(MFA's) have many applications outside of physics, combinatorial
optimization by mean-field annealing being perhaps the most prominent
of them~\cite{pean87,nish01}. One problem with MFA's is that because
they treat all particles as independent, they can be a poor
approximation of the desired distribution. Another is that especially
in large systems with complex Hamiltonians, they can be expensive to
evaluate.

We extend MFA's by introducing coordinate transformations that allow
us to convert a coupled distribution into a decoupled one. This allows
us to improve the accuracy of the MFA. It also allows us to employ
Monte Carlo mean-field annealing for optimal adaptive control. We also
introduce separate Hamiltonians for each particle. This allows us to
exploit recent results in collectives theory~\cite{wolp03a} to speed
up MFA evaluations. This extended MFA justifies heuristics recently
found to beat simulated annealing in adaptive control problems by
orders of magnitude~\cite{wowh00, wotu02c, wotu03a}. It also provides
novel ways to do (constrained) combinatorial optimization, a novel
(bounded rational) game theory equilibrium
concept~\cite{auha92,pisl02}, and a way to calculate such equilibria.





\section{Self-consistent Distributions}

Let $\zeta$ be a space with elements $z$ and cardinality $|\zeta|$.  A
{\it semi-coordinate system} of $\zeta$, $\Xi$, is a space $\xi$
having elements $x = (x_1 \in \xi_1, \ldots, x_{M} \in \xi_{M})$
together with a surjective map from $\xi$ to $\zeta$, written as
$\zeta(.)$.  $\zeta(.)$ need not be invertible; the (semi-)coordinates
of a point $x \in \xi$ may not be unique. No a priori restriction is
made on whether $\xi$ and $\zeta$ are countable, uncountable,
time-extended, etc.

The space of possible probability distributions (or density functions,
as the case may be) over $\xi$ is $\cal{P}^\dagger$. Any $p \in
{\cal{P}}^\dagger$ induces a distribution over $\zeta$, $p(z) = \int
dx \; p(x) \delta(z - \zeta(x))$. (We implicitly assume integral
measures and delta functions that match $\xi$.)  Expressions like
``$p(\xi) = p'(\xi)$'', ``$p_i = p'_i$" and ``$p_{(i)} = p'_{(i)}$'',
where $p, p' \in {\cal{P}}^\dagger$, mean $\forall x \in \xi$, $p(x) =
p'(x)$, $p(x_i) = p'(x_i)$, and $p(x_j) = p'(x_j) \; \forall j \ne i$,
respectively.  $\cal{P}$ is the set of all product distributions over
$\xi$, i.e., the submanifold of all $p \in {\cal{P}}^\dagger$ obeying
$p(\xi) = \prod_i p(\xi_i)$. We will use the obvious parameterization
of elements of $\cal{P}$ as the vectors of their marginal
distributions, written $q = (q_1(\xi_1), \ldots, q_M(\xi_M))$. Note
that changing $\zeta(.)$ with $\xi$ fixed will change the manifold
$\cal{P}$ in general but won't affect $\cal{Q}$, the space of all $q$.

We are interested in {\it self-consistent} $q$ that obey $q_i(\xi_i) =
A_i(\xi_i, q_{_{(i)}}) \; \forall i$ for some functions $\{A_i\}$. By
Brouwer's theorem, for any smooth \{$A_i$\}, the map $q \rightarrow
\{A_i(\xi_i, q_{_{(i)}})\}$ (which we call {\it parallel Brouwer
updating}) has at least one fixed point. Here we consider such $A_i$
that set $q_i$ by minimizing a functional of $(q_i, q_{_{(i)}})$. As
an example, that functional might be the same for all $i$, for each
one measuring the error of $(q_i, q_{_{(i)}})$ as an approximation of
some $p \in {\cal{P}}^\dagger$. In this case the fixed point(s) of the
\{$A_i$\} are locally optimal approximations in $\cal{P}$ of that $p$.

Write the cross entropy from $p$ to $p'$ as $S(p \; || \; p') \equiv
-\int dx \; p(x) {\mbox{ln}}[\frac{p'(x)}{\mu(x)}]$, where as usual
$\mu$ is a prior probability over $\xi$ that ensures the argument of
the logarithm is unitless. So the entropy of $p$ is just $S(p \; || \;
p) \equiv S(p)$, and the Kullback-Leibler distance from $p$ to $p'$ is
$KL(p \; || \; p') \equiv S(p \; || \; p') - S(p)$~\cite{coth91}.
Next define {\it private Hamiltonians} $h_{i,j}$ where $\zeta(x) =
\zeta(x') \Rightarrow h_{i,j}(x) = h_{i,j}(x') \; \forall i, j, x$ and
$x'$, so expected values of a Hamiltonian have the same value in
$\zeta$ and $\xi$.  We write $h_i(\xi) \equiv \sum_{j=1}^{m_i}
\beta_{i,j} h_{i,j}(\xi)$ where the temperatures $\{\beta_{i,j}\}$ are
usually non-negative and write the associated free energies as
${F}_{h_i}(p) \equiv \int dx \; p(x) h_{i}(x) \; - \; S(p)$.  So
argmin$_{q'_i \in {\cal{Q}}_i} {F}_{h_i}(q'_i, q_{_{(i)}})$ is the
$q'_i$ maximizing the entropy of $q'_i(\xi_i) \prod_{k \ne i}
q_{k}(\xi_{k})$, subject to expected values of $i$'s private
Hamiltonians implicit in $i$'s temperatures. The associated fixed
points obey

$ $

$
q_i^{g}(\xi_i) \; = \;
\frac{\mu_i(\xi_i)}{ N_i({[h_i]}_{i,q^g})} 
e^{-{[h_i]}_{i,q^g}(\xi_{i})}
$

$ $

\noindent where $N_i([U]_{i,p})$ is the partition function and the
bracket notation indicates expectation conditioned on the value of the
subscripted coordinate: for any $U$ and $p \in {\cal{P}}^\dagger$,
$[U]_{i,p}(\xi_{i}) \equiv \int d{x_{(i)}} U(\xi_i, x_{(i)}) p(x_{(i)}
\mid \xi_i)$.

We can choose \{$A_i$\} based on objective functions other than the
free energy; the important thing is to incorporate a concave function
like entropy, to avoid the expensive process of checking the borders
of $\cal{Q}$ for fixed points. For example, all the local minima of
${F}_{h_i}(q)$ over ${\cal{Q}}$ are interior to ${\cal{Q}}$.
Moreover, for any fixed $q_{_{(i)}}$, there is only one local minimum
over ${\cal{Q}}_i$.

The following six examples all have uniform $\mu$.  $\xi = \zeta$ (as
in conventional MFA's) in all but the fifth example:

$ $

{\it Example 1}: Let $m_i = 1 \; \forall i$. Then at $q^g$
each coordinate $i$ obeys its own canonical ensemble and is
independent of the values of the other coordinates. It is coupled to
its own heat bath, with its own effective Hamiltonian set by the
distributions over the other coordinates, ${[h_i]}_{i,q^g}(\xi_i)$.
In contrast, if $\forall i, m_i = 2$ and $h_{i,2} = (h_{i,1})^2$, each
$q^{g}_i(\xi_i)$ is a Gaussian over the values ${[h_i]}_{i,1}(\xi_i)$
rather than a Boltzmann distribution over them.  In this case
$\beta_{i,1}$ and $\beta_{i,2}$ specify the mean and variance of
$h_{i,1}$.

{\it Example 2}: We have a {\it team game} when $\forall i,
i'$, and $j$, $m_i = m_{i'}$, and $h_{i} = h_{i'}$.  In a team game we
have a single shared objective function trading off entropy and
expected values of the Hamiltonians.  For $m = 1$ we define the {\it
world Hamiltonian} $H$ as the single shared private Hamiltonian. For
this the minimizer of the free energy over all ${\cal{P}}^\dagger$ is
the joint canonical ensemble distribution, ${\bf{p}}^{\beta
H}$~\cite{jayn57}. In contrast, $q^{\beta H}$, the minimizer over
${\cal{P}}$, is an MFA for ${\bf{p}}^{\beta H}$.  Since $\zeta$ is
arbitrary, by gradually increasing $\beta$ one can converge on a $q$
that is a delta function about the $z$ that minimizes $H(z)$. Since
$\zeta$ is arbitrary, this optimization algorithm can be applied for
any underlying space one wishes to search. In particular, parallel
Brouwer updating with gradually increasing $\beta$ is equivalent to
mean-field annealing~\cite{pean87}.

\noindent In general $q^{\beta H}(\xi_i) \ne {\bf{p}}^{\beta
H}(\xi_i)$, the marginal of ${\bf{p}}^{\beta H}(\xi_i)$. However the
$q^{\beta{H}}$ with the lowest free energy is the $q$ that best
approximates ${\bf{p}}^{\beta{H}}$, as measured by $KL(q^{\beta H},
{\bf{p}}^{\beta H})$ To have the fixed points equal the marginals we
must instead have the $\{A_i\}$ minimize $KL({\bf{p}}^{\beta H},
q^{\beta H})$.)

{\it Example 3:} One way to enforce a constraint $f(z) = 0$ is via a
generalization of penalty functions.  Say we have an $m=1$ team game
with world Hamiltonian $H$.  Choose the objective function $\int dx
p(x) [\beta H(x) + \beta_f [f(\zeta(x))]^2] - S(p)$. Then for $\beta_f
\rightarrow \infty$ our constraint will be enforced, if possible.
Often that solution is on ${\cal{Q}}$'s border though, which slows the
search --- finite $\beta_f$ moves the solution(s) to the interior of
${\cal{Q}}$, by weakening the constraint. Alternatively, choose the
objective function $\int dx p(x) (\beta H(x) + \lambda
[f(\zeta(x))]^2) + \gamma \lambda^2 - S(p)$, where $\gamma \ge 0 $ is
a constant.  Minimizing over both $p \in {\cal{Q}}$ and $\lambda \in
\mathbb{R}$ forces $f(z) = 0$ if $\gamma = 0$, while nonzero
$\gamma_j$ weakens the constraint but forces the fixed point $q$(s)
inside ${\cal{Q}}$.

{\it Example 4}: Consider a team game with $m = m'+1$ and
$h_{j}(x) = h_{j}(x_j) \in {\mathbb{N}}^+ \; \forall j \le m'$. Here
$p^{g}$ is an MFA for the grand canonical ensemble of a system with
$m'$ particle types and Hamiltonian $h_{m'+1}$: $x_i$ encodes the
states of all type $i$ particles, and $h_{i}$ counts their
number~\cite{jayn57}.

{\it Example 5}: Say $\zeta(.)$ is bijective but $\xi \ne
\zeta$. Then $\zeta(.)$ is akin to a rotation, in that each $\xi_i$
couples multiple components of $\zeta$. More generally, $|\xi| \ne
|\zeta|$ means $\xi$ is an $M$-dimensional representation of an $M' <
M$-dimension system.  This allows arbitrary distributions over $\zeta$
to be expressed as a product distribution.

\noindent As an illustration, say $\zeta$ is two-dimensional with
elements $(z_1, z_2)$. Take $\xi_1 = \zeta_1$, and have an additional
$\xi$ coordinate for each separate value of $z_1$, which we write as
$\xi_{z_1}$. So there are $1+|\zeta_{1}|$ coordinates altogether. Have
each $\xi_{z_{1}}$ contain $|\zeta_2|$ elements, i.e., $|\zeta_2|$
separate subsets of $\zeta$, and label those subsets by the values
$z_2 \in \zeta_2$.  Have $\zeta(x_1, \ldots, x_{1 + |\zeta_1|}) =
(x_1, x_{x_1})$, i.e., $z_1 = x_1$ and $z_2 = x_{z_1}$. So we can
write $p(z_1) = \sum_{z_2} p(z_1, z_2) = \sum_{z_2} \sum_{x: x_1 =
z_1, x_{z_1} = z_2} q(x) = q_{x_1}(z_1)$. Therefore $p(z_2 \mid z_{1})
= q(\xi_{z_{1}} = z_2 \mid \xi_1 = z_1)$. So $p(z_2 \mid z_{1}) =
q_{z_1}(z_2)$.

\noindent This allows us to do optimal distributed adaptive control in
$\zeta$ with an MFA over $\xi$. Let $\zeta_2$ be the state of the
plant one wants to control, $\zeta_{1}$ the control variable you can
set, and ${\cal{H}}(\zeta_2)$ the objective function one wants to
minimize. The goal in optimal control is to find argmin$_{z_{1}}
E({\cal{H}} \mid z_{1})$. Have the $q_{i \ge 2}$ fixed to $P(\zeta_2
\mid \zeta_{1})$, the distribution relating the plant variable and the
control variables. Take $m_1 = 1$, $h_{1,1}(\xi) =
{\cal{H}}(\zeta(\xi))$, and redefine $F_{h_1}$ to involve the entropy
of just $\xi_1$ rather than all of $\xi$.  Anneal $\beta \rightarrow
\infty$, so that $q_1(\xi_1)$ becomes a single delta function. Write
the resultant $q^{\beta H}(\xi)$ as $ \delta_{\xi_1, x_1} \prod_{i \ge
2} q_i(\xi_i)$. This distribution minimizes $\sum_{x'} h_{1,1}(x')
\delta_{x'_1, x_1} \prod_{i \ge 2} q_i(x'_i) = \sum_{x'} h_{1,1}(x')
q_{x_1}(x'_{x_1}) = \sum_{z} {\cal{H}}(z) p(z_2 \mid z_1 = x_1)$. So
$z_1 = x_1$ solves our optimal control problem.

{\it Example 6}: Formally, when $m_i = 1 \; \forall i$ and $\xi =
\zeta$, each coordinate is a player in a non-cooperative game
~\cite{auha92,pisl02}, with $x$ the players' joint move and
\{$-h_{i,1}$\} their payoff functions. As the $\beta_i \rightarrow
\infty$, $q^{\{\beta_i h_{i,1}\}}$ becomes a (normal form) mixed
strategy Nash equilibria.

\noindent Now consider the general bounded rational (i.e., non-Nash equilibrium)
scenario. Say the players are told the entropy of the joint system,
and each finds its best possible mixed strategy subject to the other
players' strategies and entropy value. Then the system is at a $q^g$.

\noindent As an alternative interpretation of $q^g$, say we have a {\it
rationality} functional $R(U, q_i)$ that measures how peaked any $q_i$
is about argmax$_{x_i} U(x_i)$ for any $U : \xi_i \rightarrow
\mathbb{R}$. We require that $R(U, q_i) = \beta$ if $q_i(\xi_i)
\propto \mu(\xi_i) e^{-\beta U(\xi_i)}$, and that the $q_i$ satisfying
$R(U, q_i) = \beta$ that has maximal entropy is $\mu(\xi_i) e^{-\beta
U(\xi_i)} / N_i(U)$.  Now say we are told the \{$\beta^*_i$\}\}, the
rationalities of the $M$ players for their associated effective
Hamiltonians \{${[h_{i,1}]}_{i,q}$\}. Then the information-theoretic
optimal estimate for the associated $q$ is the minimizer over $q$ and
the \{$\lambda_i$\} of the free energy $\sum_i \lambda_i
[f_i(R({[h_{i,1}]}_{i,q}, q_i)) - f_i(\beta^*_i)] - S(q)$ for any
monotonically increasing functions $\{f_i\}$~\cite{jayn57}. At any
local minimum of this free energy $q = q^g$ with $\beta_i = \beta^*_i
\; \forall i$.

\noindent This use of rationality functionals expresses the
bounded-rational solution to any non-team game as the minimizer of a
single objective function whose local minima are all interior to
$\cal{Q}$.  As an illustration, choose $f_i(\beta)$ to be the ideal
expected Hamiltonian $\partial_{\beta} {\mbox{ln}}(N_i(\beta
{[h_{i,1}]}_{i,q}))$, and $R(U, q_i) \equiv {\mbox{argmin}}_\beta \;
KL(q_i \; || \; \frac{e^{-\beta U}}{N_i(\beta U)})$. For such a choice
$f_i(R({[h_{i,1}]}_{i,q}, q_i))$ is the actual expected Hamiltonian,
$\int dx_i \; q_i(x_i) {[h_{i,1}]}_{i,q}(x_i)$.









\section{Finding Fixed Points}

When we have an $m=1$ team game we can use variants of gradient
descent to find minima of the (single) free energy. Another approach
for such scenarios is parallel Brouwer updating. More generally,
define the {\it free energy gap} at $q$ for coordinate $i$ as
ln$[N_i({[h_i]}_{i,q^g})] + \int dx_i q_i(x_i) {[h_i]}_{i,q^g}(x_i) +
\int dx_i q_i(x_i) {\mbox{ln}}\frac{q_i(x_i)}{\mu_i(x_i)}$.  This is
how much ${F}_{h_i}$ is reduced if only $q_i$ undergoes the Brouwer
update.  Define {\it serial} Brouwer updating as only updating one
$q_i$ at a time. In an $m=1$ team game, any such update must reduce
${F}_{\beta{H}}(p)$, in contrast to the case with parallel Brouwer
updating.  In {\it greedy} serial Brouwer updating, instead of cycling
through all $i$, at each iteration we update only the coordinate with
the largest gap; this maximizes the free energy drop in that update.

A practical difficulty with these schemes for finding fixed points is
that evaluating ${[h_i]}_{i,q^g}$ can be very difficult in large
systems.  An alternative is to use Monte Carlo simple sampling to get
estimates of the effective Hamiltonians, and use those to update $p$.
In this scheme, given a $q$ at iteration $t$, each $x_i$ is separately
set by randomly sampling $q_i(t)$, thereby generating $x(t)$. Next the
pair $(x_i(t), h_i(x(t)))$ is combined with previous pairs and the
update rule to set $q_i(t+1)$, and then the process repeats.

To simplify the analysis, consider simple gradient descent of the free
energy for the $m=1$ team game. Have $p$ be constant through each
successive block of $\tau$ timesteps, updating only when we go from
some block $m$ to block $m + 1$, with the update based on observations
during block $m$.  Say we have a team game and are at a
block-transition, $t = \tau m + 1$, and let $n_i(t) \in \nu_i(t)$ be
all information the algorithm controlling $q_i$ has at that time,
including the associated Monte Carlo samples. So we have a posterior
conditional distribution of possible gradient descent directions
$P(\vec{F}^i_{{H},\beta}(q(t)) \mid n_i(t))$, where $\vec{F}^i_{H,
\beta}(q(t))$ is the components involving coordinate $i$ of the
projection onto ${\cal{P}}$ of the free energy gradient $\nabla
F_{{H}, \beta}(q(t))$.

In gradient descent updating, that distribution should set the vector
to add to $q_i(t)$ to get $q_i(t+1)$.  More precisely, since agent $i$
knows $q_i(t)$, presuming quadratic loss reflects quality of the
update, the Bayes-optimal estimate of the gradient is the posterior
expected gradient, $\int d[q_{_{(i)}}(t)] \; P(q_{_{(i)}}(t) \mid n_i)
\times \vec{F}^i_{H, \beta}(q(t))$. Expanding, the $x_i$ component of
$\vec{F}^i_{H, \beta}(q)$ is $u_i(x_i) - \sum_{x_i} u_i(x_i) /
|{\xi}_i|$, where $u_i(x_i) \equiv \beta {[H]}_{i,q^g}(x_i) \;+$
ln$[q(x_i)]$.  Rather than evaluate the integral of
${[H]}_{i,q^g}(\xi_i)$, we can use a maximum likelihood estimator,
i.e., replace that integral with $({\hat{H}})_{i,{n_i}}(\xi_i)$, the
average of the observed $H$ values over the $\tau_{x_i}$ instances (of
the just-completed block) when $\xi_i = x_i$.

Unfortunately, often in very large systems the convergence of
$(\hat{H})_{i,n_i}(\xi_i)$ is very slow, since the distribution sampled by
the Monte Carlo to produce $n_i$ is very broad.  To address this,
posit that the differences \{$(\hat{H})_{i,q^g}(x_i) -
(\hat{H})_{i,q^g}(x'_i), \; x_i, x'_i \in \xi_i$\} are unchanged when
one replaces $H$ with some $h_i$.  This means that $q_i$(t)
is unchanged by that replacement. The set of all $h_i$ guaranteed to
have this character, regardless of the form of $q(x_{i},
\xi_{_{(i)}})$, is the set of all {\it difference Hamiltonians},
$h_i({\xi}) = {H}({\xi}) - D_i({\xi}_{_{(i)}})$ for some function $D_i$.
Now across block $m$ we are sampling $P(h_i(\xi))$ to generate $n_i$,
and then evaluating ${(\hat{h_i})}_{i,n_i}(x_i)$.  The associated
variances of values (one for each of the $x_i$),
Var$((\hat{h_i})_{i,\nu_i}(x_i))$, govern the accuracy of the estimate of
the free energy gradient.  For well-chosen $D_i$ these variances may
be far smaller than when $h_i = H$. In particular, if the number of
coordinates coupled to $\xi_i$ through ${H}$ does not grow as the
system does, often such difference Hamiltonian variances will not grow
much with system size, whereas the variances
Var$(\hat{H}_{i,\nu_i}(x_i))$ will grow greatly. Furthermore, very
often such a difference Hamiltonian is far easier to evaluate than is
$H$, due to cancellation in subtracting $D_i$.

More precisely, for practical reasons we want $i$'s update algorithm
to be robust against misperception of $q_{_{(i)}}$.  So for quadratic
loss, assuming no $\tau_{x_i} = 0$, consider the $h_i$ minimizing
$\int d{q_{_{(i)}}}$ $P(q_{_{(i)}}) [\int dn^u_i P(n^u_i \mid
n^{\xi_i}_i, q_{_{(i)}}, h_i)$ $\{{\vec{F}}^i_{H}(q(\xi)) -
{\hat{F}}^i_{n_i}(q_i)\}^2$], where $P(q_{_{(i)}})$ reflects any prior
information we might have concerning $q_{_{(i)}}$ (e.g., that it is
likely that the associated ${\vec{F}}^i_{H}(q)$ is close to that
estimated for the previous block of $\tau$ steps). Here
${\hat{F}}^i_{n_i}(q_i)$ is our estimator for ${\vec{F}}^i_{\beta
H}(q)$, and $n^u_i$ is the Hamiltonian values contained in $n_i$, the
associated $\xi_i$ values, $n_i^{\xi_i}$, being independent of $h_i$
and $q_{_{(i)}}$ and therefore fixed.

The inner integral is a sum, of the (square of the) bias
${\hat{F}}^i(q_i) - {\vec{F}}^i_{H}(q)$ with the variance, $\int
dn^u_i P(n^u_i \mid n^{\xi_i}_i, q_{_{(i)}}, h_i) \;
\{{\hat{F}}^i_{n_i}(q_i) - {\hat{F}}^i(q_i)\}^2$, where
${\hat{F}}^i(q_i) \equiv \int dn^u_i P(n^u_i \mid n^{\xi_i}_i,
q_{_{(i)}}, h_i) {\hat{F}}^i_{n_i}(q_i)$.  By only considering
difference utilities we guarantee that the bias equals 0. Now expand
our variance of vectors as a sum of variances of scalars, one for each
value $x_i$.  Since $n_i$ is IID generated that sum is $\frac{|\xi_i|
- 1}{|\xi_i|} \sum_{x_i} {\mbox{Var}}((\hat{h_i})_{i,\nu_i}(x_i)) /
\tau_{x_i}$.  The difference Hamiltonian minimizing this is ${H}(s) -
\sum_{x'_i} \frac{\tau_{x'_i}^{-1}}{\sum_{x''_i} \tau_{x''_i}^{-1}}
{H}(x'_i, {x}_{_{(i)}})$~\cite{wolp03a}.  (If all the terms in this
sum cannot be stored because $|\xi_i|$ is too large, then
$n_i^{\xi_i}$ must be averaged over as well, which can be approximated
by replacing the $\tau_{x_i}$ with $q(x_i)\tau$.)  Being independent
of $q_{_{(i)}}$, this Hamiltonian minimizes our $q_{_{(i)}}$ integral,
regardless of $P(q_{_{(i)}})$. For the same reason it is optimal if if
the integral is replaced by a worst-case bound over $q_{_{(i)}}$.

An extensive series of experiments have been conducted with this
optimal Hamiltonian under the approximation of uniform $q(x_i)$, and
under the approximation that $\tau_{x_i} = 0$ for one and only one
$x_i$ value.  These experiments compared this Hamiltonian with the
team game, for many different $H$, under a variant of the parallel
Brouwer update rule~\cite{wowh00, wotu02c} that crudely corrected for
non-stationarity. In that work these algorithms were semi-formally
justified as ways to minimize $H$, with no appreciation for free
energy, self-consistency, or the like. Indeed, due to the nonlinearity
of the Brouwer update rule's dependence on ${[h_i]}_{i,q^g}$, bad
distortions may arise with the update used in that previous work,
$q(x_i) \rightarrow e^{-\beta E(h_i \mid n_i, \xi_i = x_i)} /
\sum_{x'_i} e^{-\beta E(h_i \mid n_i, \xi_i = x_i)}$.  Also, there are
no known assurances that small Var$(h_i \mid n_i, \xi_i = x_i)$ result
in more accurate updating with parallel Brouwer updating, even when
that updating is done correctly.  Despite these shortcomings, in the
experiments in~\cite{wowh00, wotu02c} the approximations to the
optimal Hamiltonian minimized $H$ up to orders of magnitude faster
than team games, with the improvement growing with problem size.


\section{AVOIDING LOCAL MINIMA}

Say we have an $m=1$ team game, are ultimately interested in annealing
$\beta$, and are currently at a local minimum $q^{H}
\in \cal{Q}$ of the shared free energy. Then to break out of that
minimum we can simply raise $\beta$ and restart the updating, since we
want to raise $\beta$ anyway, and in general doing so will change the
${F}_{H}$ so that the Lagrange gaps become nonzero.

A way to break free without changing the free energies is to switch to a
coordinate system $\Xi^2$ for $\xi$, and thereby change $\cal{P}$.  As
an example, $\Xi^2$ can join two components of $\xi$ into an aggregate
coordinate. Since we can now have statistical dependencies between
those two components, the $\xi^2$ space product distributions map to a
superset of ${\cal{P}}$. In general the local minima of that superset
do not coincide with local minima of ${\cal{P}}$.

Less trivially, say $\xi^2 = \xi$, and $\xi(.)$ is the identity map
for all but a few components, indicated as indices $1 \rightarrow
n$. Have $\xi(.)$ be a bijection, so that for any fixed $x^2_{n+1
\rightarrow M} = x_{n+1 \rightarrow M}$, the effect of the coordinate
transformation is merely to ``shuffle'' the associated mapping taking
coordinates $1\rightarrow n$ to $\zeta$. Say we have a $m=1$ team
game, and set $q^{\xi^2}_{n+1 \rightarrow M} = q^{\xi}_{n+1
\rightarrow M}$. This means we can estimate the expectations of $\beta
H$ conditioned on possible $x^2_{1\rightarrow n}$ from the Monte Carlo
samples conditioned on $\xi(x^2_{1\rightarrow n})$. So for any
$\xi(.)$ we can estimate $E(H)$ as $\int dx^2_{1\rightarrow n}
p^{\xi^2}(x^2_{1\rightarrow n}) E(H \mid \xi(x^2_{1, \ldots,
n}))$. Now entropy is the sum of the entropy of coordinates $n+1
\rightarrow M$ plus that of coordinates $1 \rightarrow
n$. Accordingly, for any choice of $\xi(.)$ and
$q^{\xi^2}_{1\rightarrow n}$, we can approximate ${L}^{\xi}_{H}$
as (our associated estimate of) $E(H)$ minus the entropy of
$p^{\xi^2}_{1\rightarrow n}$, minus a constant unaffected by choice of
$\xi(.)$.

So for finite and small enough $|\xi_{1\rightarrow n}|$, we can use
our estimates $E(H \mid \xi(x^2_{1 \rightarrow n}))$ to search
for the ``shuffling'' $\xi(.)$ and distribution
$q^{\xi^2}_{1\rightarrow n}$ that minimizes ${L}^{\xi}_{H}$
(penalizing by the bias$^2$ plus variance expression if we intend to
do more Monte Carlo). The search can involve a series of free energy
descents over $\xi^2_{1\rightarrow n}$ for each possible $\xi(.)$, or
use cruder heuristics, like having $q^{\xi^2}_{1\rightarrow n} =
q^{\xi}_{1\rightarrow n}$, and only varying $\xi(.)$. Not only should
this coordinate transformation lower the free energy, it should also
result in a new surface through ${\cal{P}}^\dagger$ that is no longer
at a local minimum. More generally, for arbitrary $\xi^2$ we can bound
$F_{\xi^2} \ge F_{\xi^1} \ge F_{\xi^2} - {\mbox{max}}_{x^1}$
[ln$(\sum_{x^2} \delta_{\xi^1(x^2), x^1})]$ (and similarly for
uncountable $\xi^2$). So if after switching to $\xi^2$ we can then
reduce $F_{\xi^2}$ to less than the value $F_{\xi^1}$ had when we made
the switch, then we know we have also reduced $F_{\xi^1}$ to below
that pre-switch value. An upper bound on the how large that drop in
$F_{\xi^2}$ needs to be is ${\mbox{max}}_{x^1}$ [ln$(\sum_{x^2}
\delta_{\xi^1(x^2), x^1})]$.








$ $

\noindent
I would like to thank Chiu Fan Lee, Bill Macready, and Stefan
Bieniawski for helpful discussion.



\begin{thebibliography}{10}

\bibitem{auha92}
R.J. Aumann and S.~Hart.
\newblock {\em Handbook of Game Theory with Economic Applications}.
\newblock North-Holland Press, 1992.

\bibitem{chjo02}
D.~Challet and N.~F. Johnson.
\newblock Optimal combinations of imperfect objects.
\newblock {\em Phys. Rev. Let.}, 89:028701, 2002.

\bibitem{kono03_short}
G.~Korniss et~al.
\newblock {\em Science}, 299:677, 2003.

\bibitem{jayn57}
E.~T. Jaynes.
\newblock Information theory and statistical mechanics.
\newblock {\em Physical Review}, 106:620, 1957.

\bibitem{nish01}
H.~Nishimori.
\newblock {\em Statistical physics of spin glasses and information processing}.
\newblock Oxford University Press, 2001.

\bibitem{pean87}
C.~Peterson and J.~R. Anderson.
\newblock A mean-field theory learning algorithm for neural networks.
\newblock {\em Complex Systems}, 1:995, 1987.

\bibitem{pisl02}
E.~W. Piotrowski and J.~Sladkowski.
\newblock An invitation to quantum game theory.
\newblock /quant-ph/0211191, 2002.

\bibitem{wotu02c}
D.~Wolpert and K.~Tumer.
\newblock Beyond mechanism design.
\newblock In H.~Gao et~al., editor, {\em International Congress of
  Mathematicians 2002 Proceedings}. Qingdao Publishing, 2002.

\bibitem{wolp03a}
D.~H. Wolpert.
\newblock Theory of collective intelligence.
\newblock In K.~Tumer and D.~H. Wolpert, editors, {\em Collectives and the
  Design of Complex Systems}, New York, 2003. Springer.

\bibitem{wotu03a}
D.~H. Wolpert, K.~Tumer, and E.~Bandari.
\newblock Improving search by using intelligent coordinates.
\newblock 2003.
\newblock submitted.

\bibitem{wowh00}
D.~H. Wolpert, K.~Wheeler, and K.~Tumer.
\newblock Collective intelligence for control of distributed dynamical systems.
\newblock {\em Europhysics Letters}, 49(6), March 2000.

\end{thebibliography}




\end{document}













\section{The Mathematics of Designing Collectives}
\label{sec:math}
Our uncertainty concerning the state of the system is reflected in a
probability distribution over $\zeta$. Our ability to control the
system consists of setting the value of some characteristic of the
collective, e.g., setting the functions $\u{g}_i$.  Indicating that
value of the {\bf global coordinate} by $s$, our analysis revolves
around the following {\bf central equation} for $P(G \mid s)$, which
follows from Bayes' theorem:

\begin{eqnarray}
& & P(G \mid s) = \\ \nonumber
& &  \int  d{N}_G P(G \mid {N}_G, s) \int
d{N}_g P({N}_G \mid {N}_g, s) P({N}_g \mid s) \; ,
\label{eq:central}
\end{eqnarray}

\noindent
where ${N}_g$ and ${N}_G$ are the {\bf intelligence} vectors
of the particles with respect to $h_\eta$ and $G$,
respectively. Intelligence, defined as the ``standardization'' of
utility functions so that the numeric value they assign to a $z$ only
reflects their ranking of $z$ relative to certain other elements of
$\zeta$, is given by:
\begin{eqnarray}
N_{\eta,U}(z ) \equiv \int d\mu_{z_{\;\hat{}\;\eta}} (z')
	\Theta[U(z) - U(z')] \; ,
\end{eqnarray}

\noindent
where $\Theta$ is the Heaviside function, and where the subscript on
the (normalized) measure $d\mu$ indicates it is restricted to $z'$
sharing the same non-$\eta$ components as $z$. 

Note that $N_{\eta,h_\eta}(z) = 1$ means that particle $\eta$ is fully
rational at $z$, in that its move maximizes the value of its utility,
given the moves of the particles. In other words, a point $z$ where
$N_{\eta,h_\eta}(z) = 1$ for all particles $\eta$ is one that meets the
definition of a game-theory Nash equilibrium.  On the other hand, a
$z$ at which all components of ${N}_G = 1$ is a maximum $G$ along
all coordinates of $z$ (which of course does not mean it is a maximum
along an off-axis direction). So if we can get these two points to be
identical, then if the particles do well enough at maximizing their
impulse functions we are assured we will be near an axis-maximizing
point for $G$.

To formalize this, consider our decomposition of $P(G \mid s)$. If we
can choose $s$ so that the third conditional probability in the
integrand is peaked around vectors ${N}_g$ all of whose components
are close to 1, then we have likely induced large (impulse
function) intelligences.  Intuitively, this ensures that the 
impulse functions have high ``signal-to-noise''.  If we can also have
the second term be peaked about ${N}_G$ equal to ${N}_g$, then
${N}_G$ will also be large.  It is in the second term that the
requirement that the impulse functions be ``aligned with $G$''
arises. Note that our desired form for the second term in
Equation~\ref{eq:central} is assured if we have chosen private
utilities such that ${N}_g$ equals ${N}_G$ exactly for all
$z$. Such a system is said to be {\bf factored}.  Finally, if the
first term in the integrand is peaked about high $G$ when ${N}_G$
is large, then our choice of $s$ will likely result in high $G$, as
desired.  In this letter we concentrate on the second and third terms,
and show how to simultaneously set them to have the desired form.

As an example, any ``team game'' in which all impulse
functions equal $G$ is factored~\cite{crba96}.  However team games
often have very poor forms for term 3 in Equation~\ref{eq:central},
forms which get progressively worse as the size of the collective
grows. This is because for such impulse functions each particle
$\eta$ will usually confront a very poor ``signal-to-noise'' ratio in
trying to discern how its state affects its utility $h_\eta = G$,
since so many other particle's states also affect $G$ and therefore
dilute $\eta$'s effect on its own impulse function.

We now focus on algorithms based on impulse functions
\{$h_\eta$\} that optimize the signal/noise ratio reflected in the
third term, subject to the requirement that the system be factored.
To understand how these algorithms work, say we are given an arbitrary
function $f(z_\eta)$ over particle $\eta$'s moves, two such moves
${z_\eta}^1$ and ${z_\eta}^2$, a utility $U$, a value $s$ of
the global coordinate, and a move by all particles other than $\eta$,
$z_{\;\hat{}\;\eta}$. Define the associated {\bf learnability} by
\begin{eqnarray}
& & \Lambda_f(U ; z_{\;\hat{}\;\eta}, s, {z_\eta}^1, {z_\eta}^2)
\equiv \\ \nonumber
& & \sqrt{ \frac{ [E(U ; z_{\;\hat{}\;\eta}, {z_\eta}^1) -
 E(U ; z_{\;\hat{}\;\eta}, {z_\eta}^2)]^2}
{\int dz_\eta [f(z_\eta)Var(U ; z_{\;\hat{}\;\eta}, z_\eta)] }} \;
 .
\label{eq:learnability}
\end{eqnarray}
\noindent
The expectation values in the numerator are formed by averaging over
the training set of the learning algorithm used by particle $\eta$,
$n_\eta$. Those two averages are evaluated according to the two
distributions $P(U | n_\eta) P(n_\eta | z_{\;\hat{}\;\eta},
{z_\eta}^1)$ and $P(U | n_\eta) P(n_\eta | z_{\;\hat{}\;\eta},
{z_\eta}^2)$, respectively. (That is the meaning of the semicolon
notation.)  Similarly the variance being averaged in the denominator is
over $n_\eta$ according to the distribution $P(U | n_\eta) P(n_\eta |
z_{\;\hat{}\;\eta}, {z_\eta})$.

The denominator in Equation~\ref{eq:learnability} reflects how
sensitive $U(z)$ is to changing $z_{\;\hat{}\;\eta}$. In
contrast, the numerator reflects how sensitive $U(z)$ is to
changing $z_{\eta}$. So the greater the learnability of a impulse function $h_\eta$, the more $h_\eta(z)$ depends only on
the move of particle $\eta$, i.e., the better the associated
signal-to-noise ratio for $\eta$. Intuitively then, so long as it does
not come at the expense of decreasing the signal, increasing the
signal-to-noise ratio specified in the learnability will make it
easier for $\eta$ to achieve a large value of its intelligence. This
can be established formally: if appropriately scaled, $g'_\eta$ will
result in better expected intelligence for particle $\eta$ than will
$h_\eta$ whenever $\Lambda_f(g'_\eta; z_{\;\hat{}\;\eta}, s,
{z_\eta}^1, {z_\eta}^2) \; > \; \Lambda_f(h_\eta;
z_{\;\hat{}\;\eta}, s, {z_\eta}^1, {z_\eta}^2)$ for all
pairs of moves ${z_\eta}^1,
{z_\eta}^2$ \cite{wolp03aa}. 

One can solve for the set of all impulse functions that are factored
with respect to a particular world utility. Unfortunately though, in
general a collective cannot both be factored and have infinite
learnability for all of its particles~\cite{wolp03a}.  However consider
{\bf difference} utilities, of the form
\begin{eqnarray}
U(z) = \beta[G(z) - D(z_{\;\hat{}\;\eta})]
\end{eqnarray}
\noindent
Any difference utility is factored~\cite{wowh00}. In addition, for
all pairs ${z_\eta}^1, {z_\eta}^2$, under benign approximations the
difference utility maximizing $\Lambda_f(U ; z_{\;\hat{}\;\eta}, s,
{z_\eta}^1, {z_\eta}^2)$ by
\begin{eqnarray}
D(z_{\;\hat{}\;\eta}) = 
  E_f(G(z) \mid z_{\;\hat{}\;\eta}, s) \;,
\end{eqnarray}
\noindent
up to an overall additive constant, where the expectation value is
over $z_\eta$. We call the resultant difference utility the {\bf
Aristocrat} utility ($AU$), loosely reflecting the fact that it measures
the difference between a particle's actual state and the average
state.  If each particle $\eta$ uses an appropriately rescaled version
of the associated $AU$ as its impulse function, then we have
ensured good form for both terms 2 and 3 in Equation~\ref{eq:central}.

Using $AU$ in practice is sometimes difficult, due to the need to
evaluate the expectation value. Fortunately there are other utility
functions that, while being easier to evaluate than $AU$, still are both
factored and possess superior learnability to the team game utility,
$h_\eta = G$. One such impulse function is the {\bf Wonderful
Life} Utility (WLU).  The $WLU$ for particle $\eta$ is parameterized by a
pre-fixed {\bf clamping parameter} $CL_\eta$ chosen from among
$\eta$'s possible moves:
\begin{eqnarray}
WLU_\eta \equiv G(z) - G(z_{\;\hat{}\;\eta}, CL_\eta) \; .
\end{eqnarray}
\noindent
WLU is factored no matter what the choice of clamping
parameter. Furthermore, while not matching the high learnability of
$AU$, $WLU$ usually has far better learnability than does a team game, and
therefore (when appropriately scaled) results in better expected
intelligence \cite{tuwo00,wotu99a,wotu01a}.


***


At this point we have transformed the problem to that of finding the
optimal \{$D_i$\}. To address it in this new form, first note that the
iterative algorithm governing particle $i$'s interactions with the
rest of the system is a standard algorithm from the field of
Reinforcement Learning (RL) ~\cite{suba98}. That field studies
algorithms that enable an ``agent'' (here a particle $i$) to use past
experience to choose its ``actions'' in such a way as to maximize the
the resultant value of its (private) ``utility'' function. In RL, the
algorithm outlined above is called a ``Boltzmann learning algorithm''.
The use of finite $\beta$ in such an algorithm is traditionally
motivated as a way to ``trade off exploration and exploitation'' in
the agent's search for how best to increase its rewards. (Accordingly,
in the absence of annealing, in general $\beta$ should not be
infinite.)

Next, consider a system containing many agents each of which can be
viewed as using an RL algorithm to maximize an associated private
utility function. Combined with a {\it world utility} function $G$
grading the behavior of such a system as a whole, that system is a
{\it collective} ~\cite{wowh00, wotu01a, wolp03a}. So our problem is
how to modify the private utilities of the agents in a collective to
speed up convergence of the world utility --- a problem that has been
extensively studied recently.

To illustrate this connection, identify the world utility with the
full Hamiltonian of the system.  So we are not concerned directly with
minimizing free energy but rather with minimizing
${G}$. Accordingly, identify the action of particle $i$ with its
state $s_i$ and its reward with $h_i({s})$. (Especially if we were
directly concerned with minimizing free energy, it might make sense to
consider an alternative identification where the action of each agent
is the Euclidean vector of its marginal distribution.) In such a
context, $G = -{G}$ and we call $G - D_i$ the {\it difference
utility} for agent $i$.  (The minus signs are due to the fact that
game theory, RL and the study of collectives traditionally casts
problems as maximization of utilities rather than minimizations of
Hamiltonians).




****

\section{Combination of Imperfect Objects}

We illustrate the foregoing with the problem of combining imperfect
objects ~\cite{chjo02}. This problem involves a set of noisy
observational devices each producing a single real number by sampling
a Gaussian of fixed width, where all the Gaussians are centered on the
same value.  The problem is to choose the subset of those devices
having smallest average (over the subset) error.

Formally, the problem is to minimize $\epsilon \equiv \frac{|
\sum_{j=1}^m x_j a_j |} {\sum_{k=1}^m x_k}$
where $x_j \in \{0, 1\}$ is whether device $j$ is selected, and there
are $m$ devices in the collection, having associated distortions
$\{a_j\}$.  We identify $\epsilon$ with the world utility (so that
here the goal is to minimize $G$, not maximize it) and assigned one
agent to set each $x_j$.

The agents' estimates of their private utilities were formed in an
initial 20 time-step period in which all agent's moves were chosen
uniformly randomly.  (All private utilities result in the same average
performance during that period, since none are actually used then.)
Following this, with each successive $t$ a fixed fraction of the
agents still moving randomly switch to using their Boltzmann learning
algorithms instead. This gradual introduction of learning ``softens
the discontinuity'' in each agent's environment when the other agents
start using their learning algorithms and therefore change their
moves.  For $m=100$, four agents turned on their learning algorithms
at each time step (i.e., by $t=26$, all agents were using their
learners). For simplicity, given how short our runs were, $\beta$ was
fixed throughout any run.

\begin{figure} [htb]
%\vspace*{-.2in}
\centerline{\psfig{figure={100.ps},height=2.2in,width=3.4in}}
\caption{Combination of Imperfect Objects Problem, $m=100$.}
\label{fig:n=100}
%\vspace*{-.2in}
\end{figure}

Figure~\ref{fig:n=100} shows the convergence properties of different
private utilities in a system with 100 agents.  The results reported
are based on 20 different $\{a_j\}$ configurations, each performed 50
times (i.e., each point on the figure is the average of $20 \times 50
= 1000$ runs).  $G$, $AU$, and $WLU$ show the performance of agents
using Boltzmann learners with private utilities provided by $G$, $AU$,
and $WLU$, respectively.  $S$ shows the performance of greedy search
where new $x_j$'s are generated at each step and selected if the
solution is better than the current best solution. Because the runs
are only 200 timesteps long, algorithms such as simulated annealing
(SA) do not outperform simple search: there is simply no time for
annealing.

Because of the environmental discontinuity experienced by each agent
as the other agents turn on their learning algorithms, $G$ during the
turn-on period may degrade temporarily, before the system settles
(e.g., this is the case with $WLU$).  After all settling, systems in
which agents use the $G$ utility perform poorly, due to large noise in
their $\hat{h}_i$, while systems in which agents use $AU$ perform
best, and the search algorithm's performance falls between the two.

\begin{figure} [htb]
\centerline{\psfig{figure={scaleW.ps},height=2.0in,width=3.4in}}
\caption{Scaling in the Combination of Imperfect Objects Problem.}
\label{fig:scaleW}
%\vspace*{-.2in}
\end{figure}

We also investigated the scaling properties with each private utility.
Figure~\ref{fig:scaleW} shows the results (the $t=200$ average
performance over 1000 runs).  As $m$ grows two competing factors come
into play.  On the one hand, there are more degrees of freedom to use
to minimize $G$. On the other hand, the problem becomes more
difficult: the search space gets larger for $S$, and there is more
noise in the system for the learning algorithms. To account for these
effects and calibrate the performance values as $m$ varies, in the
figure we also provide the performance of the algorithm that randomly
selects its state (``Ran'').  Note that the difference between the
performances of $S$ and $AU$ increases when the system size increases,
up to a factor of twenty for $m=1000$.

All algorithms but $AU$ have slopes similar to that of
``Ran'', demonstrating that they cannot use the additional
degrees of freedom provided by the larger $m$.  Only $AU$ effectively
uses the new degrees of freedom, providing gains that are
proportionally higher than the other algorithms (i.e., the rate at
which $AU$'s performance improves outpaces what is ``expected'' based
on the random algorithm's performance).


%\nocite{joja98,zhan98,huho88,wada92}





Any coordinate system \{$\xi_i$\} of the lattice ${\xi}^N$ is a set
of partitions of that lattice, one partition for each $\xi_i$, where
specifying the joint value ($x_i \in \xi_i \; \forall i$) uniquely
fixes $s \in \xi^N$~\cite{wolp03a}. One of the simplest
transformations of a coordinate system is to ``fracture a particle
into multiple particles'', i.e., replace one of the partitions $\xi_i$
with multiple partitions \{$\xi'_j$\} such that $x_i \in \xi_i
\rightarrow \{x'_j \in \xi'_j\}$ is invertible.  Another is to
aggregate multiple particles into a composite particle (i.e., replace
a set of partitions \{$\xi_i$\} with a new partition $\xi'_j$ where
each $x'_j \in \xi'_j$ is the intersection of (the $s$'s in) one
partition element from each of the $\xi_i$. More sophisticated
transformations allow each new partition not to be aligned with some
subset of the $N$ axes of ${\xi}^N$, as they are in
aggregating/splitting. The first-order self-consistent distribution is
then redefined to be with respect to the variables $x_i$ specifying
the elements of the new partitions.  Some future work investigates the
effects of changing the coordinate system on the free energy of the
associated first-order self-consistent distribution and on convergence
to that distribution.


\footnote{Formally, define the vector $q_w$ to 
be a probability distribution over possible values of $s_{(i)}$, 
evaluated for $s_{(i)} = w$. Note that we don't know $p$, only 
$n_i$. Accordingly write $h_i(s_i, p) = \sum_w {G}(s_i, w) 
q_w$. So $E(h_i \mid n_i) = \sum_{w} {G}(s_i, w) [\int dp P(p 
\mid n_i) q_{w}] = \sum_{w} {G}(s_i, w) P(s_{(i)} = w \mid 
n_i)$.}


To use non-difference utilities this way, we need a way to find a
guaranteed associated self-consistent distribution, just like free
energy descent does this for team games. To that end first define
particle $i$'s free energy by $F_{i,g_i,\beta}(P) \equiv E(g_i) -
\beta^{-1}S(q_i)$, and the associated {\it free energy floor} of
particle $i$ by $F^*_{i,g_i,\beta}(P) \equiv
-\beta^{-1}$ln$[\sum_{s_i} e^{-\beta E(g_i \mid s_i)}]$. These are the
actual free energy for particle $i$ for private Hamiltonian $E(g_i
\mid s_i)$ and the lowest free energy for that Hamiltonian that could
arise by varying $P(s_i)$ while keeping $P(s_{(i)})$ fixed.

$F_{i,g_i,\beta}(P) \ge F^*_{i,g_i,\beta}(P) \; \forall i$, with
equality $\forall i$ iff we are at a self-consistent Boltzmann
distribution.  Therefore we can always find such a distribution by
gradient descent in any function of $P$ that has a local minimum only
at $P$ where $F_{i,g_i,\beta}(P) = F^*_{i,g_i,\beta}(P) \; \forall
i$. One example of such a function is $\sum_i F_{i,g_i,\beta}(P) -
F^*_{i,g_i,\beta}(P)$.



%As another example, say that there are constraints coupling the
%components of $s$, where we can express those constraints by
%identifying each $s_i$ with a different set of components $z_{[j
%\rightarrow j+K]}$ of some vector $z$, the overlaps between those sets
%imposing the constraints. Then the $K$-th order self-consistent
%distribution with minimal free energy has $F_i(z_{j \rightarrow j+K})
%\propto e^{}.


This leads to a novel energy-minimization
algorithm where each particle Monte-Carlo estimates its own spectrum,
and then independently samples a Boltzmann distribution across that
spectrum to set its state.  This algorithm can be used even if the
functional form of $G$ is unknown (i.e., we only have Monte
Carlo samples of it).

Recent results in the field of collectives~\cite{wolp03a} can be used
to improve the convergence of this algorithm.  The resultant hybrid
framework, intimately related to bounded rationality non-cooperative
game theory, provides novel explanations for many previously observed
experimental results concerning collective systems~\cite{wowh00,
wotu01a, wotu02c}. It also makes predictions concerning such systems,
in particular concerning how to improve their performance, in
particular how to approximate optimal distributed control, even when
the the separate controllers cannot fully communicate with one
another.


Additional benefits of minimizing variances arise if the distributions
$P(g_i(s_i, \xi_{(i)}) \mid n_i))$ and $P({G}(s_i, \xi_{(i)}
\mid n_i))$ for particle $i$ are all related by affine
transformations, i.e., that $\forall \; s_i, P(g_i(s_i, \xi_{(i)}) =
y \mid n_i) = P(K{G}(s_i, \xi_{(i)}) + h_i(s_i) = y \mid
n_i))$ for some constant $K$ and function $h_i$
\footnote{This is a generalization of assuming that the distributions
$P({G}(s_i, \xi_{(i)} \mid n_i))$ are all Gaussians, and that
for the underlying distribution $P(s_{(i)})$ replacing ${G}
\rightarrow g_i$ rescales those Gaussians by the same factor $K$ and
then shifts them by potentially varying amounts, $h_i(s_i)$.}.
In~\cite{wolp03a} it is proven that such affineness means that
shrinking the variances for particle $i$ cannot change
$\frac{\exp{-\beta <g_{i,n_i}(s_i)>}} {\sum_{s'_i} \exp{-\beta
<g_{i,n_i}(s'_i)>}}$ to de-emphasize the preferable $s_i$ having small
${[g_i]}_{i,q^g}(s_i)$. 

The only assumption about $p(\xi_{(i)})$ used to derive this result
is the implicit one in the affine property; no more detailed
information concerning $p(\xi_{(i)})$ and its relation to $n_i$ is
needed. While this assumption often will not hold exactly, it often
holds to a good approximation in many large highly coupled systems due
to the central limit theorem. 




Now by the weak law of large numbers, for any fixed positive real
number $\epsilon$, $P(n_i : |{[g_i]}_{i,q^g}(s_i) -
\hat{g}_{i,\nu_i}(s_i)| > \epsilon) <
\frac{{\mbox{Var}}(\hat{g}_{i,n_i}(s_i))}{\epsilon^2 L_{s_i}}$.  This
suggests that if particle $i$ replaces its estimates of $G$
values with those of an appropriately chosen difference utility, and
uses those estimates in its free energy descent update rule, the
algorithm will converge to the same $p$ but far more quickly.


We are interested in vectors each of whose $N$ components is a
distribution over $\xi$, ${\bar{p}}(\xi) \equiv
({\bar{p}}_1(\xi), \ldots, {\bar{p}}_N(\xi))$.  Let
${\bar{{\cal{P}}}}$ be a space of such vectors. Define the associated
space $\cal{P} \subseteq {\cal{P}}^\dagger$ having elements $p$ as the
image of ${\bar{{\cal{P}}}}$ under the Hadamard product: $p(s) =
\prod_{i=1}^N {\bar{p}}_i(s) \; \forall s$. Say
that in addition to this map from ${\bar{ {\cal{P}} }}$ to $\cal{P}$ ,
we have one in the other direction.  Any $p$ at which these two maps
are inverses is {\it self-consistent}. When ${\bar{\cal{P}}}$ is the
set of vectors $\bar{p}$ such that $\forall i, {\bar{p}}_i(s)$ depends
only on the value $s_i$, we say a self-consistent $p$ is {\it first
order self-consistent}. Under any first-order self-consistent $p$ the
particle's states are independent.

The map from ${\bar{ {\cal{P}} }}$ to $\cal{P}$ is continuous, and
therefore for any smooth \{$A_i$\} so is the compound endomap
${\cal{P}} \rightarrow {\bar{ {\cal{P}} }} \rightarrow {\cal{P}}$. So
by Brouwer's fixed point theorem that endomap has at least one fixed
point, i.e., at least one self-consistent distribution must exist for
any such \{$A_i$\}.  We call this endomap the {\it parallel Brouwer
update rule}.

Furthermore, many search algorithms (e.g., simulated annealing
(SA)) intersperse ``exploitation'' and ``exploration'' steps.
Typically such algorithms make random moves in the exploration step.
A natural alternative is to sample (an estimate of) the
self-consistent distribution to explore. (From a game-theoretic
perspective, this is ``agentizing'' the variables of a search
problem.) In recent work this led to substantially faster minimization
of $G$ ~\cite{wotu03a} than SA, with the improvement growing
with problem size.


Let ${\bar{p}} \equiv ({\bar{p}}_1, \ldots, {\bar{p}}_M)$ be a vector
of distributions over $\xi$, and let ${\bar{\cal{P}}}$ be a space
of such vectors. Let $\cal{P}$ be the associated space of
distributions given by $p(s) = \prod_{j=1}^M {\bar{p}}_j(s)$.  Say
that in addition to this product map from ${\bar{\cal{P}}}$ to
$\cal{P}$ we have a smooth map from $\cal{P}$ into ${\bar{\cal{P}}},
which we label $B$.  Any $p$ at which these two maps are inverses is
{\it self-consistent}. Now every $B$ induces an endomap $p \rightarrow
{\bar{p}} \rightarrow p'$, which we call the {\it parallel Brouwer
update rule}. Applying Brouwer's fixed point theorem to this rule
shows that there is at least one self-consistent point for any $B$.

We say $\bar{p}$ and the $p$ formed from it are {\it first order} if
$\forall j, {\bar{p}}_j$ lives in $\cal{P}^\dagger_j$. By restricting
${\bar{\cal{P}}}$ to only contain first-order points, we see that for
any associated $B$ there must be first-order self-consistent $p$'s.
Such a $p$ obeys $p(s) = \prod_j {\bar{p}}_j(s_j)$, i.e., the
particles are independent under such a $p$.

Define the entropies $S(p) \equiv -\sum_{s} p(s) {\mbox{ln}}[(p(s)) /
\mu(s)]$ and $S({\bar{p}}_i(\xi)) \equiv -\sum_{s} p(s)
{\mbox{ln}}[({\bar{p}}_i(s)) / \mu(s)]$.  Let the \{$g_{i,j}$\} be
arbitrary real-valued functions of ${s}$, and define $E(g_{i,j}(s_i,
\xi_{(i)})) \equiv \sum_{(i)} \; g_{i,j}(s_i, s_{(i)}) \;
p(s_{(i)})$. Then $p(s_i) \propto \mu(s_i) \exp{-[ \sum_{j=1}^{m_i}
\beta_{i,j} E(g_{i,j}(s_i, \xi_{(i)}))]}$ is the $p \in {\cal{P}}$
that maximizes $S(p(\xi_i))$ subject to a fixed $p(\xi_{(i)}) \in
{\cal{P}}_{(i)}$ and provided values of the $E(g_{i,j}(\xi_i,
\xi_{(i)}))$.  and similarly for $p(s)$.  


Then we define $b_{\{a_i\},s^*}( \{
{\bar{p}}^i \}) \equiv \prod_i \frac{{\bar{p}}^i(s_i \mid
A(\{s^*_{j>i}\}, \{s_{j < i}\}))} {{\bar{p}}^i(s^*_i \mid A(\{s^*_{j >
i}\}, \{s_{j<i}\}))}$ where $s^* \in \xi$ and each ${\bar{p}}^i \in
{\cal{P}}^\dagger$. Next define ${\cal{P}}_{\{a_i\},s^*} \subset
{\cal{P}}^\dagger$ as the set of all $p$ that are proportional to
$b_{\{a_i\},s^*}(\{{\bar{p}}^i\})$ for some set \{${\bar{p}}^i$\}.  If
$a_i(s) = s_{(i)} \; \forall i$, then up to the same overall
proportionality constant $b_{\{a_i\},s^*}(p, p, \ldots, p) = p \;
\forall p \in {\cal{P}}^\dagger$, i.e., then ${\cal{P}}_{\{a_i\},s^*}
= {\cal{P}}^\dagger$.  We will use the shorthand $\cal{P} \equiv
{\cal{P}}_{\{a_i\},s^*}$, $b \equiv b_{\{a_i\},s^*}$, for any $p \in
{\cal{P}}^\dagger$, $q_{i}(s) \equiv p(s_i \mid a_i(s_{(i)}))$ and
$q_{(i)} \equiv (q_1, \ldots, p_{i-1}, p_{i+1}, \ldots, p_N)$.

By Brouwer's theorem, any smooth map $B$ from ${\cal{P}}^\dagger$ into
${\cal{P}}$ has a fixed point; we call application of $B$ the
associated {\it Brouwer update rule}. Here we consider the rule $p
\rightarrow b(\{ {\mbox{argmin}}_{ {\bar{p}}_i \in
{\cal{P}}^\dagger_i} {F}_{g_i}({\bar{p}}_i, p_{(i)}) \})$,
where the free energies ${F}_{g_i}(p') \equiv
\sum_{j=1}^{m_i}
\beta_{i,j} g_{i,j}(s) p'(s) - H(p')$ for arbitrary {\it private
utility functions} $g_{i,j}$, $H(p')$ is the entropy, $-\sum_s p'(s)
{\mbox{ln}}[\frac{p'(s)}{\mu(s)}]$, and  usually the
calibration distribution $\mu \in {\cal{P}}$ and the
\{$\beta_{i,j}$\} are all non-negative.

To solve the minimization, at each separate possible value of $a_i(s)$
form the gradient of ${F}_{g_i}({\bar{p}}_i, p_{(i)})$ over
the $|\xi_i|$ components of ${\bar{p}}_i$, introducing an extra
Lagrange parameter to force ${\bar{p}}_i \in
{\cal{P}}^\dagger_i$. This gives ${\bar{p}}_i(s_i \mid a_i(s)) \propto
\mu_i(s) e^{-\sum_j \beta_{i,j} {[g_i]}_{i,j}(s_i,
a_i(s))}$, where ${[g_i]}_{i,j}(s_i, a_i(s)) \equiv \sum_{s'}
g_{i,j}(s') P_{s_i, a_i(s)}(s')$ and $P_{s_i, a_i(s)}(s')$ is given by
normalizing $\delta_{s'_i, s_i} \delta_{a_i(s'), a_i(s)} 
\frac{\prod_{k \ne i} p_k(s'_{k} \mid A(\{s^*_{j>k}\}, \{s_{j \le k}\})}
{\prod_{k} p_k(s^*_{k} \mid A(\{s^*_{j> \ge k}\},
\{s_{j<k}\})}$ over all $s'$.
We indicate such a fixed point of $B$ by $p^*_{g,\beta}$.


***

As usual we can view a set of lists \{${\omega}_i$\} as a directed
graph connecting the particles, where there is an edge pointing from
$i$ to $j$ iff $j \in \omega_i$~\cite{graphical models for
physicists}.  Given a set of distributions \{$p^i \in
{\cal{P}}^\dagger$\}, $p(\xi) \propto \prod_i p^i(\xi_i \mid
\xi_{{\omega}_i})$ need not imply that $p(\xi_i \mid
\xi_{\omega_i}) \propto p^i(\xi_i \mid \xi_{\omega_i})$; this
$p$ may not even be a Markov Random Field with neighborhood structure
\{$\omega_i$\}.  However if no path leading from $k$ along the
directed edges of \{$\omega_i$\} returns to $k$, then $p(\xi_k \mid
\xi_{\omega_k})$ equals $p^k(\xi_k \mid \xi_{\omega_k})$
exactly. In particular, this is true for all particles if $\forall i
\; {\omega}_i$ points only to particles with index $j > i$, If in fact
$\omega_i$ is all $j > i$ for every particle $i$, then $\prod_i
p'(\xi_i \mid \xi_{\omega_i}) = p'(\xi) \; \forall p' \in
{\cal{P}}^\dagger$, i.e., this $\{\omega_i\}$ supports an exact
approximation of any $p \in {\cal{P}}^\dagger$.

Write the entropy of $p$ as $H(p) \equiv -\sum_s p(s)
{\mbox{ln}}[\frac{p(s)}{\mu(s)}]$. Next specify free energies
${F}_{g_i}(p) \equiv \sum_{j=1}^{m_i} \beta{\omega}_{i,j}
g_{i,j}(s) p(s) - H(p)$. Usually the Lagrange parameters
\{$\beta_{i,j}$\} are all non-negative and the calibration
distribution $\mu \in {\cal{P}}$.

In general we use the symbol $q$ to refer to vectors of $N$
conditional distributions $(p^1(\xi_1 \mid \xi_{\omega_1}),
p^2(\xi_2 \mid \xi_{\omega_2}), \ldots)$. We refer to the space
of all possible $q$ as $\cal{Q}$.  So for any $q \in {\cal{Q}}$, $q_i
\in {\cal{Q}}_i$ is the set of distributions over values $s_i$
conditioned on values $s_{\omega_i}$.  Also define $A(q)$ as $\prod_i
q_i(\xi_i \mid \xi_{{\omega}_i})$ multiplied by a normalizing
constant and define ${\cal{P}} \subset {\cal{P}}^\dagger$ as the range
of $A$, i.e., the set of all $p$ that equal $A(q)$ for some $q$.  In
addition, write ${F}_{g_i}(q) \equiv
{F}_{g_i}(A(q))$. For each $i$ we call the associated
$g_{i,j}$ the {\it private utilities} of $i$. So argmin$_{q'_i \in
{\cal{Q}}_i} {F}_{g_i}(q'_i, q_{(i)})$ is the $q'_i$ that
maximizes the entropy of $A(q'_i, q_{(i)})$, subject both to the
specified $q_{(i)}$, and to expected values of $i$'s private utilities
implicitly given by the values of $i$'s Lagrange parameters.

By Brouwer's theorem, any smooth map $B$ from ${\cal{Q}}$ into itself
has a fixed point.  Here we consider the map $q \rightarrow A[\{
\; {\mbox{argmin}}_{q'_i \in {\cal{Q}}_i}
{F}_{g_i}(q'_i, q_{(i)}) \; \}]$.  To solve the
minimization for one particular $i$, separately for each value of
$s_{{\omega}_i}$ form the gradient of ${F}_{g_i}(q'_i,
q_{(i)})$ over the $|\xi_i|$ components of $q'_i(\xi_i \mid
s_{\omega_i})$, introducing an extra Lagrange parameter to force
$q'_i$ to be normalized. The result is

$ $

$
q'_i(\xi_i \mid \xi_{{\omega}_i}) \propto 
\mu(\xi_i \mid \xi_{{\omega}_i}) e^{-\sum_j
\beta_{i,j} {[g_i]}_{i,j}(\xi_{i}, \xi_{\omega_i}) - 
H({\bar{q}}_{\xi_{i}, \xi_{\omega_i} } )},
$

$ $

\noindent where ${[g_i]}_{i,j}(\xi_{i}, \xi_{\omega_i}) \equiv
\sum_{s'_{(i,{\omega}_i)}} g_{i,j}(\xi') {\bar{q}}_{\xi_{i},
\xi_{\omega_i}} (s'_{(i,{\omega}_i)})$, and $\forall s$,
${\bar{q}}_{s_{i}, s_{\omega_i}} (s'_{(i,{\omega}_i)})$ is
given by normalizing $\delta_{s_i,s'_i} \delta_{s_{\omega_i},
s'_{\omega_i}} \prod_{k \ne i} q_k(s'_{k} \mid s'_{{\omega}_k}) $ over
all $s'$.  We indicate the associated fixed point of $B$ by
$q^{g,\beta}$, and $A[q^{g,\beta}]$ by $p^{g,\beta}$.

To illustrate this result, consider the case where (as in conventional
MFA's) $\mu$ is uniform and the particles are independent, i.e.,
${\omega}_i = \O \; \forall i$. In this scenario $\forall i \;
{\bar{p}}_{s_{i}, s_{\omega_i}} (s_{(i, {\omega}_i)}) =
{\bar{p}}_{s_i} (s_{(i)})$ and is independent of $s_i$. Furthermore
there are no closed paths in \{$\omega_i$\}, so the formula for
$q^i(s_i \mid s_{\omega_i})$ gives $p^{g,\beta}(s_i \mid
s_{\omega_i})\; \forall i$. We consider several special cases of this
situation:

\noindent
{\it Example 1}: Let $m_i = 1 \; \forall i$. Then at
$p^{g,\beta}$ each particle obeys its own single-particle canonical
ensemble, being coupled to its own heat bath, and having its own
Hamiltonian ${[g_i]}_{i,q^g}(s_i) = \sum_{s'_{(i)}} g_i(s_i, s'_{(i)})
\prod_{j \ne i} p^{g,\beta}(s_j)$. In contrast, if $\forall i, m_i =
2$, and $g_{i,2} = (g_{i,1})^2$, then $\beta_{i,1}$ and $\beta_{i,2}$
specify the mean and variance of $g_{i,1}$. In this case
$p^{g,\beta}(s_i)$ is given by a Gaussian over the values
${[g_i]}_{i,1}(s_i)$ rather than a Boltzmann distribution over them.

\noindent
{\it Example 2}: We have a {\it team game} when $\forall i, i', j, m_i
= m_{i'}$, $\beta_{i,j} = \beta_{i',j}$, and $g_{i,j} = g_{i',j}$. In
this situation the free energies of all the particles are the same.  For
$m = 1$ define the {\it world utility} $G$ as the single shared
private utility. For this case the global free energy is just the
system's free energy for Hamiltonian $G$, and its minimizer over all
${\cal{P}}^\dagger$ is the canonical ensemble distribution, $p_{G,
\beta}$~\cite{jayn57}. In contrast, the minimizer over ${\cal{P}}$ is
an MFA for $p_{G, \beta}$.

\noindent
{\it Example 3}: Say we have a team game with $m = N+1$, where $g_{j}(s) =
g_{j}(s_j) \in {\mathbb{N}}^+ \; \forall j$. Then $p^{g,\beta}$ is
an ``MFA'' for the grand canonical ensemble of a system with $N$ types
of particles and Hamiltonian $g_{N+1}$: $s_i$ encodes the states of the
particles of type $i$, and $g_{i}$ is the number of such particles.


An important modification to the foregoing is to add a particle $N+1$,
still require $p(s_{-[N+1]}) = p(s_{1}, \ldots, s_N) = \prod_{i=1}^N
A_i(p(\xi_{(i)}, s_{i}))$, and make no restriction concerning
$s_{N+1}$. For this modification redefine $A_i(p(s_i, \xi_{(i)}))
\propto e^{-\beta E(g_i \mid s_i)}$ (so we implicitly require $p(s_i)
\ne 0$, which may mean $\beta < \infty$). Under this modification all
of the preceding conclusions still hold, for the vectors
$p(\xi_{-[N+1]})$ (rather than the full $p(\xi)$), provided we
take the free energy to be $F_{-[N+1],{G},\beta}(p) \equiv E({G}) -
\beta^{-1}H(p(\xi_{-[N+1]}))$.

{\bf {dhw: Change the following --- it should just amount to having
one of the distributions in the argument of $A$ in the Brouwer update
rule be fixed, rather than set by optimizing an associated
free energy. The associated $s_i$ is the plant. More generally, when
several of the $p_i$ are fixed, with their associated $\omega_i$ they
express a priori constraints on the optimization.}}

As an example, the problem of optimal control for objective $G$,
control variables $s_{-[N+1]}$, and plant variable $s_{N+1}$
is to find lim$_{\beta \rightarrow
\infty}$argmin[$F_{{G}', \beta}(p(\xi_{-[N+1]}))$]
where ${G}'(s_{-[N+1]}) \equiv E({G} \mid
s_{-[N+1]})$. However $E({G}') = \sum_{s_{-[N+1]}} E({G}
\mid s_{-[N+1]}) p(s_{-[N+1]}) = E({G})$. So
$F_{{G}', \beta}(p(\xi_{-[N+1]})) =
F_{-[N+1],{G},\beta}(p)$. Accordingly, by finding the 1st-order
Boltzmann self-consistent distribution over the control variables
$s_{-[N+1]}$ for large $\beta$ we can approximate optimal control.

Each coordinate $\xi_i$ is an $x_i$-indexed set, its elements being
subsets of $\zeta$, $\cup_{x'_{(i)}} \zeta(x_i, x'_{(i)})$. (The
special case where each coordinate is a partition of $\zeta$ is
investigated in~\cite{wolp03a}.)

The set of all $q \omega$ is defined as ${\cal{P}}_\xi \subset
{\cal{P}}^\dagger_\xi$. When $\xi$ is implicit we suppress the
subscript on ${\cal{P}}_\xi$ and ${\cal{P}}^\dagger_\xi$.


In general the mapping of $p^{\xi^2}$ onto ${\cal{P}}^{\xi^2}$ will
change the free energy. So the starting value of the free energy after
the change of coordinates may be worse than that of the original
$p^\xi$. Furthermore, if $|\xi^2| \ne |\xi|$,
${L}^{\xi^2}_{\beta G} \ne f([q^{\xi^2}])$ in general; we must
include a distortion factor to relate $H^{\xi^2}$ to $H^\xi$, which
may complicate the analysis. As an example of how to get around these
problems, take $\xi^2 = \xi \vartimes \mathbb{B}$, $\zeta(\xi, 2) =
\zeta(\xi)$, and $\zeta(\xi, 2) = \zeta^2(\xi)$ for some desired
replacement of $\zeta(.)$, $\zeta^2(.)$. In performing the
transformation, impose $p(\xi^2_{M+1}) = a \delta_{\xi^2_{M+1},1} + (1
- a)\delta_{\xi^2_{M+1},1}$ with $a = 1$, and $p(\xi_{-(M+1)}) =
p(\xi)$. (Note this is a product distribution.)  At this point
${L}^{\xi^2}_{\beta G}(p^{\xi^2}) = {L}^{\xi}_{\beta
G}(p^{\xi})$. Now decrease ${L}^{\xi^2}_{\beta G}(q^{\xi^2})$,
in the usual Monte Carlo way. However in successive blocks ``anneal"
$a$, ending at $a = 0$. At the end of this process we again have
${L}^{\xi^2}_{\beta G}(p^{\xi^2}) = {L}^{\xi}_{\beta
G}(p^{\xi})$, but now we have replaced $\zeta(.)$ with $\zeta^2(.)$
as the mapping taking $(\xi^2_1, \xi^2_2) \rightarrow \zeta$.



The potential drop from that initial value of the free energy should be
roughly indicated by the sum of the free energy gaps in the new
coordinate system, or alternatively by the magnitude of the
free energy's gradient in the new system --- such indicators can be
used to choose the transformation.


and all we know about the players is their expected payoffs. Then the
information-theoretic optimal estimate for the joint strategy of the
players is $q^{g}$ where the $\beta_i$ are given by those expected
payoffs~\cite{jayn57}. Alternatively, say player $i$ is told that the
players are statistically independent (the most common scenario
studied in game theory) and is also told (or estimates) the expected
values of the \{$g_{j \ne i, 1}$\}. Then the optimal estimate for $i$
to make of the joint strategy of the other players is $\prod_{k \ne i}
q^{g}_k(x_k)$.  So the players are ``coordinated'' by all making the
most conservative assumption they can about each other, given what
they are told about expected payoffs. An alternative view is that

cond-mat0307630.extracted_bib
\begin{thebibliography}{10}

\bibitem{auha92}
R.J. Aumann and S.~Hart.
\newblock {\em Handbook of Game Theory with Economic Applications}.
\newblock North-Holland Press, 1992.

\bibitem{chjo02}
D.~Challet and N.~F. Johnson.
\newblock Optimal combinations of imperfect objects.
\newblock {\em Phys. Rev. Let.}, 89:028701, 2002.

\bibitem{kono03_short}
G.~Korniss et~al.
\newblock {\em Science}, 299:677, 2003.

\bibitem{jayn57}
E.~T. Jaynes.
\newblock Information theory and statistical mechanics.
\newblock {\em Physical Review}, 106:620, 1957.

\bibitem{nish01}
H.~Nishimori.
\newblock {\em Statistical physics of spin glasses and information processing}.
\newblock Oxford University Press, 2001.

\bibitem{pean87}
C.~Peterson and J.~R. Anderson.
\newblock A mean-field theory learning algorithm for neural networks.
\newblock {\em Complex Systems}, 1:995, 1987.

\bibitem{pisl02}
E.~W. Piotrowski and J.~Sladkowski.
\newblock An invitation to quantum game theory.
\newblock /quant-ph/0211191, 2002.

\bibitem{wotu02c}
D.~Wolpert and K.~Tumer.
\newblock Beyond mechanism design.
\newblock In H.~Gao et~al., editor, {\em International Congress of
  Mathematicians 2002 Proceedings}. Qingdao Publishing, 2002.

\bibitem{wolp03a}
D.~H. Wolpert.
\newblock Theory of collective intelligence.
\newblock In K.~Tumer and D.~H. Wolpert, editors, {\em Collectives and the
  Design of Complex Systems}, New York, 2003. Springer.

\bibitem{wotu03a}
D.~H. Wolpert, K.~Tumer, and E.~Bandari.
\newblock Improving search by using intelligent coordinates.
\newblock 2003.
\newblock submitted.

\bibitem{wowh00}
D.~H. Wolpert, K.~Wheeler, and K.~Tumer.
\newblock Collective intelligence for control of distributed dynamical systems.
\newblock {\em Europhysics Letters}, 49(6), March 2000.

\end{thebibliography}