text:
# Searching for Plannable Domains can Speed up Reinforcement Learning

## 1. Introduction

### 1.1. Reinforcement Learning

Reinforcement learning involves sequential decision making in uncertain environments. The sequential aspect of the decision problem reflects the fact that the immediate cost or benefit of any state of the environment may play only a small part in determining the true value of any state. The aim of the decision-making agent is to develop a *policy*, which maximizes the benefit of acting in its environment over an extended period of time.

An often used and efficient framework, which describes stochastic, sequential decision problems, is the *Markov decision process* (MDP) (for a review, see \[Puterman, 1994\]). When a problem description satisfies the requirements of the MDP framework, well-known algorithms can be used to determine an optimal policy, such as various forms of the dynamic programming \[Bellman, 1957, Bertsekas, 1987, Sutton, 1991\], Q-learning \[Watkins, 1989\], or SARSA \[Singh and Sutton, 1996\] methods.

These algorithms proceed by maintaining a *value function*, updating it according to the experiences and propagating the values of states. Under appropriate conditions, RL algorithms are shown to develop an optimal policy \[Bertsekas, 1987, Singh et al., 2000\]. However, the basic forms of these algorithms propagate the values step-by-step, therefore convergence may be very slow. Various methods have been developed to overcome this difficulty, such as prioritized sweeping \[Moore and Atkeson, 1993\], eligibility traces \[Sutton, 1988, Singh and Sutton, 1996\] and also planning methods to be described below.

### 1.2. Planning

Besides RL, another successful and well-studied approach to solving decision problems is *planning*. As a central problem of classical AI research, effective algorithms have been developed to solve planning problems. However, classical decision-theoretic methods usually assume that the world is deterministic. This is an appropriate approximation in some cases, but *not always*. There may be certain parts of the world, that are highly stochastic, making classical planning unreliable, or even useless. It is a plausible idea to integrate planning and reinforcement learning to gain more general and robust methods.

Sutton and colleagues as well as others \[Sutton, 1991, Sutton, 1990, Peng and Williams, 1993, Forbes and Andre, 2000\] have integrated planning and learning in the *Dyna* architecture: Planning is treated as being virtually identical to reinforcement learning. Learning updates the appropriate value function estimates according to experience as it actually occurs, whereas planning updates the same value function estimates for simulated transitions chosen from the world model.

Other researchers \[Sacerdoti, 1977, Korf, 1985\] used planning to develop *macro actions* (i.e., fixed sequences of actions) that could speed-up value propagation and learning. A macro action could be a complete sub-policy (such as ‘go to the door’, ‘search a wall’, etc.) \[Hauskrecht et al., 1998, Precup et al., 1998, McGovern and Sutton, 1998\]. The main difficulty with macro actions is how to construct them: they must be either handcrafted (see, e.g., \[Kaelbling, 1993, Kalmár et al., 1998\]) or one must try to generate them automatically (see, e.g., \[Dietterich, 2000\]). In this latter case a great number of useless macros might be generated, which might even deteriorate learning \[McGovern and Sutton, 1998, Kalmár and Szepesvári, 1999\].

### 1.3. Thesis of the Paper

We propose a method called plannable RL (pRL), that has similarities with both planning methods mentioned: we make focused value updates based on hypothetical experiences gained from a model. This is similar to the Dyna architecture. On the other hand, we maintain two separate (but interacting) value functions, one for learning and one for planning. We use a simple model to evaluate the value function for planning: values are updated when transitions are considered plannable (i.e., when those are close to deterministic). This planning-value function is then used to compute macro actions. We show that these macro actions are *near-optimal*. This means that we do not need to generate large numbers of (possibly bad) macro actions, and fast learning is still possible.

### 1.4. Structure of the Paper

First, the basis of RL methods are reviewed (Section 2. The pRL algorithm is described and a pseudo-code of the algorithm is provided in Section 3. The near optimality of the algorithm is proven in Section 4. Computational demonstration are provided in the next section, i.e., in Section 5. The paper is finished by a discussion (Section 6).

## 2. Preliminaries

### 2.1. Reinforcement Learning (RL) and the MDP framework

Let us recall the definition of a Markov decision process (MDP) \[Puterman, 1994\]. A (finite) MDP is defined by the tuple $`X,A,R,P`$, where $`X`$ and $`A`$ denotes the finite set of states and actions, respectively. $`P:X\times A\times X[0,1]`$ is called the transition function; $`P(x,a,y)`$ gives the probability of arriving at state $`y`$ after executing action $`a`$ in state $`x`$. Finally, $`R:X\times A`$ is the reward function, $`R(x,a)`$ gives the immediate reward for choosing action $`a`$ in state $`x`$.

At each sequence of discrete time steps, $`t=0,1,\mathrm{}`$, the problem-solving agent observes its world to be in state $`s_tS`$ and executes an action $`a_tA`$. After this action, the agent receives reward $`r_{t+1}`$ from the world and observes the next state, $`s_{t+1}`$ according to the functions $`R`$ and $`P`$.

The agent’s objective is to choose each action $`a_t`$ so as to maximize the expected discounted return, $`_t`$, i.e., $`E\left(_t=r_{t+1}+\gamma r_{t+2}+\gamma ^2r_{t+3}+\mathrm{}\right)`$, where $`0\gamma <1`$ is the discount factor and $`E(.)`$ denotes the expected value of the argument.

A Markovian policy $`\pi `$ is a $`S\times A[0,1]`$ mapping, where $`\pi (x,a)`$ is the probability that in state $`x`$ the agent selects action $`a`$. The value of a state $`x`$ under policy $`\pi `$ is the expected value of the total discounted reward of starting from $`x`$ and then following policy $`\pi `$. Formally, $`V^\pi (x)=E_{P,\pi }(_{t=0}^{\mathrm{}}\gamma ^tr_t|s_0=x)`$ (i.e., the expected value depends on the transition probabilities and the policy applied). These values satisfy the recursive equations

$$V^\pi (x)=\underset{aA}{}\pi (x,a)\left(R(x,a)+\gamma \underset{yS}{}P(x,a,y)V^\pi (y)\right)\text{ for all }x.$$

A policy $`\pi ^{}`$ is optimal, if for all policies $`\pi `$, $`V^\pi ^{}(x)V^\pi (x)`$ for all $`xS`$ (it is easy to show that such a policy exists). The corresponding value function satisfies

(2.1)
$$V^{}(x)=\underset{aA}{\mathrm{max}}\left(R(x,a)+\gamma \underset{yS}{}P(x,a,y)V^{}(y)\right).$$

A standard way to find an optimal policy is to compute the optimal value function $`V^{}:X`$, which gives the value (i.e., the expected accumulated discounted reward) for each starting state. From $`V^{}`$, the optimal policy can be derived: the ‘greedy’ policy with respect to the optimal value function is an optimal policy (see, e.g., \[Sutton and Barto, 1998\] and references therein). This is the well-known *value iteration approach* \[Bellman, 1957\]. Functions $`V^\pi `$, $`V^{}`$ are called value functions, which are associated with a fixed policy $`\pi `$ and with the optimal policy, respectively.

We can also define the action value function $`Q^\pi (x,a)`$ as the expected value of the total discounted reward of starting from state $`x`$, choosing action $`a`$ and then following policy $`\pi `$. This value function is often more useful, because it provides the values of individual actions. It is easy to see that

$$Q^\pi (x,a)=R(x,a)+\gamma \underset{yS}{}P(x,a,y)V^\pi (y),$$

$$V^\pi (x)=\underset{aA}{}\pi (x,a)Q^\pi (x,a)$$

and thus

$$Q^\pi (x,a)=R(x,a)+\gamma \underset{yS}{}P(x,a,y)\underset{a^{}A}{}\pi (y,a^{})Q^\pi (y,a^{}).$$

The optimal action-value function $`Q^{}`$ satisfies an equation analogous to $`V^{}`$:

$$Q^{}(x,a)=R(x,a)+\gamma \underset{yS}{}P(x,a,y)\underset{a^{}A}{\mathrm{max}}Q^{}(y,a^{}).$$

### 2.2. Dynamic Programming and Two Basic RL Algorithms

Equation 2.1 can be used as an iteration:

(2.2)
$$V_{t+1}(x)=\underset{aA}{\mathrm{max}}\left(R(x,a)+\gamma \underset{yS}{}P(x,a,y)V_t(y)\right)$$

for every state $`x`$, and for an arbitrary $`V_0`$. The iteration is known to converge to $`V^{}`$. The method is called *value iteration*.

If the model of the environment (i.e. $`R`$ and $`P`$) is not known, or the state space is too large for solving Eq. (2.2), then sampling methods can be used. One such algorithm is called *Q-learning* \[Watkins, 1989\]. Q-learning uses the following iteration:

(2.3)
$$Q_{t+1}(s_t,a_t)=(1\alpha _t)Q_t(s_t,a_t)+\alpha _t\left(r_t+\gamma \underset{aA}{\mathrm{max}}Q_t(s_{t+1},a)\right),$$

where $`s_t`$, $`a_t`$ and $`r_t`$ are the state of the system, the selected action, and the immediate reward at time step $`t`$, respectively, $`Q_0`$ is arbitrary, and $`0\alpha _t1`$ is the learning rate.

It has been shown that if $`\alpha _t`$ converges to 0 properly ($`_t\alpha _t`$ is divergent, but $`_t\alpha _t^2`$ is convergent), and if every $`(x,a)`$ pair is updated infinitely often, then $`Q_t`$ converges to $`Q^{}`$ with probability 1. The proof can be found, e.g., in \[Singh et al., 2000\].

Another RL method is SARSA, which takes sampling to the extreme. It has an update rule similar to (2.3):

(2.4)
$$Q_{t+1}(s_t,a_t)=(1\alpha _t)Q_t(s_t,a_t)+\alpha _t\left(r_t+\gamma Q_t(s_{t+1},a_{t+1})\right),$$

where $`Q`$-values of the iteration are action value estimates of two state-action pairs, the one which just occurred and its predecessor. SARSA is convergent under the same assumptions as Q-learning \[Singh et al., 2000\]. A comprehensive overview of various reinforcement learning methods can be found in \[Sutton and Barto, 1998\] and in the references therein.

### 2.3. Planning with Dyna

Reinforcement learning often requires a large number of experiences (i.e. $`(s_t,a_t,r_t,s_{t+1})`$ tuples) to develop an appropriate policy. This is efficient when experience can be collected quickly, provided that one can afford the cost of explorations. *Dyna* offers a solution to improve the exploitation of previous experiences by integrating planning into the learning process.

Informally, planning is a process of computing a (near-)optimal policy for the existing (possibly inaccurate) model of the environment. Planning is an off-line method, it improves the policy without invoking additional interactions with the environment. DP, for example, executed on the available model, is a planning algorithm. Limitation arises if the model is not given: experience is needed to build a model. Another problem appears when DP is computationally intensive: performing a single and *complete* DP iteration requires $`O(|S|^2)`$ computation steps. Efforts have been made to overcome this drawback, yielding solutions like *prioritized sweeping* \[Moore and Atkeson, 1993\], or *real-time dynamic programming* \[Barto et al., 1995\].

Another successful approach is to combine reinforcement learning and dynamic programming in a single algorithm. This approach was first suggested by Sutton, who called the algorithm *Dyna*. The basic idea is that one DP update on one single state can be interpreted as a reinforcement learning step. In the Dyna architecture the agent repeats the following steps: (1) obtain experience from the environment, (2) use this to update the value function, (3) use the experience to improve the model of the environment (e.g. to approximate the transition probabilities), (4) obtain *hypothetical experience* from the updated model, and (5) use the hypothetical experience to update the value function. Note that Dyna focuses DP iterations on previously visited (and thus presumably important) regions of the state space, and therefore it might reduce the computational demand of DP significantly.

Dyna allows the adjustment of planning relative to collecting experiences: if there is no time for planning, steps (4) and (5) can be omitted. Conversely, if real experience is slow or costly, then multiple planning steps may be accomplished.

### 2.4. Planning with Macro Actions

Another possibility for utilizing additional computational capacity is to compile compound actions (macro actions) from basic ones. In its simplest form, macros are fixed sequences of actions (e.g. ‘go forward, turn right, go forward’). Many works have dealt with different versions of the macro concept. In particular, one might consider policies of sub-problems (e.g., policy for ‘finding a wall’ or policy for ‘going to the door’) as macros. In this case separate value functions on separate parts of the state-space will arise.

In order to integrate macros into reinforcement learning, two additional steps must be implemented: (6) the generation of macro actions and (7) the evaluation of macro actions. The latter one can be accomplished, e.g., by computing $`Q(x,a_{\text{macro}})`$ in all necessary states. Generating appropriate macro actions is a more difficult problem. One commonly used approach is to pre-wire them by hand \[Kaelbling, 1993\]. This is a straightforward way to encode prior knowledge into the learning problem if such prior knowledge is available and if the environment is steady. Efforts have been made also to generate macro actions in an automated fashion. The general approach is to use some heuristics to compile a number of macros from basic actions, evaluate the macros, and then keep the ones which are the most useful.

A useful set of macro actions can substantially speed up learning, because the agent can make larger steps toward its goal. On the other hand, evaluating bad macros can deteriorate performance (because trying a bad macro takes a large step in a wrong direction).

## 3. The pRL algorithm

In this section a novel algorithm is proposed, which combines DP and macros. First, a model is built for model-based value-function updates, just like in Dyna. These updates are used directly to find useful macro actions and to calculate their utility. A parallel approximate value function is maintained, which encodes the macros and their values according to the model. Greedy policy with respect to this approximate value function is applied. The resulting policy is then used to choose a macro action.

An important question is the choice of the model. One might attempt to learn an approximation of the $`P(x,a,y)`$ transition probabilities and the corresponding rewards. However, estimating and maintaining a table of size $`|S|^2|A|`$ may not be feasible. Finally, we think that planning is more effective in near-deterministic domains of the state space, and is less advantageous when little is known about the outcome of an action. Such problems frequently arise in real life. For example,

* when one has to pass a city, one might choose going through the city directly or taking the ring around the city. The first option could be less expensive in many cases, but large durations on travel time may arise if a traffic jam occurs in the inner areas. The second option may allow for more accurate planning and may have a competitive cost on the average.
* Requests can be distributed for mobile agents in many different ways. Planning becomes crucial when part of the system may be(come) unreliable.

Our planning algorithm (the planner) stores only those transitions that are almost deterministic, which we shall call *plannable transitions*.111Note that we allow for almost deterministic transitions but we do not assume an almost deterministic environment here. The algorithm will be called plannable-RL, or *pRL*, for short.

### 3.1. Basic Definitions, the Value Functions and their Update Rules

In the following we specify the pRL algorithm. First of all, we define plannable transitions formally:

###### Definition 3.1.

A transition $`(x,y)`$ is called *plannable with accuracy $`\kappa `$*, if it can be realized with probability $`\kappa `$, i.e. there exists an action $`a(x,y)`$ such that $`P(x,a(x,y),y)\kappa `$.

Let us denote the set of states that are plannable from state $`x`$ by $`T(x)`$, i.e.

$$T(x):=\{y:(x,y)\text{is }\kappa \text{-plannable}\}=\{ya(x,y):P(x,a(x,y),y)\kappa \}$$

We assume that we are given an inverse dynamics $`\varphi :X\times XA`$ such that $`\varphi (x,y)=a(x,y)`$, if $`(x,y)`$ is plannable, and it is arbitrary otherwise. This is a reasonable assumption; either the inverse dynamics is known in advance, or, the ‘learning by doing’ method can be sufficient to approximate the inverse dynamics for the plannable domains.

To compute $`T(.)`$, we approximate $`P(x,\varphi (x,y),y)`$ with $`\widehat{P}_t(x,y)`$ using the following approximation scheme:

(3.1)
$$\widehat{P}_{t+1}(x,y)=\{\begin{array}{cc}(1\alpha _t)\widehat{P}_t(x,y)+\alpha _t1\hfill & \text{if }s_t=x\text{}a_t=\varphi (x,y)\text{ and }s_{t+1}=y,\hfill \\ (1\alpha _t)\widehat{P}_t(x,y)+\alpha _t0\hfill & \text{if }s_t=x\text{}a_t=\varphi (x,y)\text{ and }s_{t+1}y,\hfill \\ \widehat{P}_t(x,y)\hfill & \text{otherwise},\hfill \end{array}$$

$$\widehat{P}_01.$$

Note that this iteration approximates the exact transition probabilities. A transition $`(x,y)`$ is considered plannable, iff $`\widehat{P}_t(x,y)\kappa `$. These almost-sure transitions will be considered as *sure* by our planning algorithm. The approximation simplifies the required computations, and – as we shall show later – provides near-optimal solutions.

###### Definition 3.2.

Connected components of plannable transitions are called *plannable domains*.

Note that plannable domains are disjoint. In turn, the number of plannable domains is smaller than the number of states.

The immediate rewards of plannable transitions are similarly approximated by $`\widehat{R}(x,y)`$.

In our algorithm, two value functions are maintained: a standard (basic) action-value function $`Q(x,a)`$, $`aA`$, which is updated by real experiences, and another value function called the planning value function. The planning value function will make use of a Dyna-like algorithm above the simplified model described in Eq. (3.1). Both functions suggest – possibly different – policies. In every time step, pRL can switch between these policies by examining which one seems better.

### 3.2. Updating the basic value function

For updating the basic value function, any traditional RL update rule can be applied. For example, one may use DP, Q-learning, SARSA, etc. We use SARSA because of its simplicity:

$$Q_{t+1}(s_t,a_t)=(1\alpha _t)Q_t(s_t,a_t)+\alpha _t\left(r_t+\gamma Q_t(s_{t+1},a_{t+1})\right),$$

where $`s_t`$, $`a_t`$ and $`r_t`$ are the state of the system, the selected action and the reward at time step $`t`$, respectively. $`Q_0`$ is arbitrary, and $`0\alpha _t1`$ denotes the learning rate. Note that the basic value of state $`x`$ at time $`t`$ is given by $`V_t(x)=\mathrm{max}_aQ_t(x,a)`$.

### 3.3. Updating the planning value function

Several update rules could be used for the planning value function. However, there is an important difference: an (approximate) model and an inverse dynamics is known for this case. Therefore We do not need to maintain an action-value table, instead a simple state value function can be used. This function will be denoted by $`\widehat{V}(x)`$. The corresponding policy can be determined, e.g., by the inverse dynamics. Note that this policy may differ from the policy suggested by the basic value function.

We chose the value iteration update rule in the following form:

$$\widehat{V}_{t+1}(x)=\underset{a}{\mathrm{max}}\underset{y}{}(P(x,a,y)(R(x,a,y)+\gamma ^{}\widehat{V}_t(y))).$$

The number of non-vanishing terms of this update may be considerably smaller than the number of all possible terms, because in our approximate model all transition probabilities are either 0 or 1. Transitions with low probability are omitted, whereas transitions with high probability (determined by $`\kappa `$) are considered as sure transitions. Note that this simplification could mislead action selection. The error depends on the degree of the simplification, which is determined by $`\kappa `$: lower $`\kappa `$ results in a coarser model. Assuming that immediate reward $`R`$ is approximated by $`\widehat{R}_t`$, the update rule can be rewritten as

(3.2)
$$\widehat{V}_{t+1}(x)=\mathrm{max}\{\underset{yT(x)}{\mathrm{max}}(\widehat{R}_t(x,y)+\gamma ^{}\widehat{V}_t(y)),V_t(x)\}.$$

Here we took into consideration that sometimes it is better to select actions according to the basic value function $`V_t(x)`$ (i.e. continue sampling or quit planning and start sampling) than to select action according to the planning value function $`\widehat{V}_t(x)`$) (i.e., continue planning, or quit sampling and start planning). Note that when no plannable states are available, one has to return to the basic value function. One may think of this technique that the space of plannable actions is extended by a new pseudo-action, which could be called ‘*stop planning*’: The choice of this action retreats to sampling and to action selection by means of the basic value function.

The update rule is applied for the planning value action in the plannable area around the current state $`x`$ in every time step. This area can be determined by a limited breadth-first search on the graph of plannable transitions, starting from state $`x`$. The search phase and the update phase require no interactions with the real system, these are ‘off-line’ evaluations. Due the to limited search and the limited number of DP updates, they take only $`𝒪(C)`$ steps.

### 3.4. Action selection

In classical RL methods (Q-learning, SARSA, etc.), $`\epsilon `$-soft or $`\epsilon `$-greedy policies with respect to the actual value function can serve to support exploration and exploitation (see, e.g., \[Sutton and Barto, 1998\]). Now we have two different value functions which may suggest different policies. We shall make use of the $`\epsilon `$-greedy policy when $`Q_t`$ is used, but if the value of planning from the actual state is higher, i.e. $`\widehat{V}_t(x)>V_t(x)`$, then action will be generated by means of the planning value function $`\widehat{V}_t`$. This decision is made by pRL in every state.

### 3.5. Plannable transitions as macros

Macros are not represented explicitly in pRL, only through the planning value function. This is advantageous because we do not need much space to store the macros - and we can still store the most useful macro with its value for each state.

The macro encoded by $`\widehat{V}`$ is *”choose the greedy action according to $`\widehat{V}_t`$, while $`\widehat{V}_t>V_t`$. Stop macro, if an action leads to a state other than the planned one.”*

More formally, let

$$S(x):=\mathrm{arg}\underset{y}{\mathrm{max}}(\widehat{R}_t(x,y)+\gamma ^{}\widehat{V}_t(y))$$

be the greedy selection according to $`\widehat{V}`$, and let

$$a(x):=\varphi (x,S(x)).$$

Then the macro of state $`x`$ is the action sequence $`a(x),a(S(x)),a(S(S(x))),\mathrm{},a(S^n(x)),\mathrm{}`$, with the stopping condition ”stop at the $`n`$th step, if either $`\widehat{V}_t(S^{n+1}(x))<V_t(S^{n+1}(x))`$, or taking action $`a(S^n(x))`$ in state $`S^n(x)`$ leads to a state other than $`S^{n+1}(x)`$.

Naturally, such listing of the macros can be included into the pRL algorithm, but is not necessary.

The pseudo-code of the algorithm is summarized in Fig. 1.

## 4. Near-optimality of pRL Macros

It is easy to see that the basic value function of pRL converges to the optimal value function. This follows from the fact that the convergence theorem for SARSA and Q-learning \[Singh et al., 2000\] requires (i) sufficient exploration (every $`(x,a)`$ pair should be visited infinitely often) and (ii) setting the learning rate properly ($`_t\alpha _t=\mathrm{}`$, but $`_t\alpha _t^2<\mathrm{}`$). Clearly, if these criteria are satisfied with standard RL algorithms, then they are also satisfied in pRL. However, convergence to the optimal policy is guaranteed only if macros are valued correctly, i.e., if $`\kappa =1`$. The rate of convergence can also be seriously affected by $`\kappa `$: $`\kappa <1`$ allows for a larger set of macros and may converge much faster to (possibly suboptimal) solutions. Therefore, it is an important issue how the learning and the utilization of macros will influence performance. For the analysis of this problem, an extension of the classical MDP model is necessary. For this reason, we briefly introduce $`\epsilon `$-MDPs \[Szita et al., 2002b\], and review related theorems. The $`\epsilon `$-MDP-theory will be applied to pRL to show that macros are near-optimal.

### 4.1. $`\epsilon `$-MDPs

The RL problem can be extended so that the environment is no longer required to be an MDP, it is only required to remain ‘near’ to an MDP, i.e. the environment is allowed to change over time, even in a non-Markovian manner.

The closeness of two environments (which have the same state- and action-sets) is measured by the distance of their transition functions. We say that the distance of two transition functions $`P`$ and $`P^{}`$ is $`\epsilon `$-small ($`\epsilon >0`$), if $`P(x,a,.)P^{}(x,a,.)_{L_1}\epsilon `$ for all $`(x,a)`$, i.e. $`_y|P(x,a,y)P^{}(x,a,y)|\epsilon `$ for all $`(x,a)`$. (Note that for a given state $`x`$ and action $`a`$, $`P(x,a,y)`$ is a probability distribution over $`yX`$.)

A tuple $`X,A,\{P_t\},R`$ is an $`\epsilon `$-MDP with $`\epsilon >0`$, if there exists an MDP $`X,A,P,R`$ (called the base MDP) such that the difference of the transition functions $`P`$ and $`P_t`$ is $`\epsilon `$-small for all $`t=1,2,3,\mathrm{}`$.

As a simple example for an $`\epsilon `$-MDP, consider an ordinary MDP perturbed by a small noise in each time step.

### 4.2. Convergence of Learning Algorithms in $`\epsilon `$-MDPs

One expects that such small perturbations in the environment may not disturb the performance of the learning algorithms very much. Nevertheless, one cannot expect that any algorithm finds an optimal value function in an $`\epsilon `$-MDP. (Such a solution may not even exist because of the perturbations of the environment). However, we shall guarantee that these algorithms find near-optimal value functions. Formally, the task is to show that $`limsup_t\mathrm{}V_tV^{}\text{const}\epsilon `$ (or $`limsup_t\mathrm{}Q_tQ^{}\text{const}\epsilon `$),222Unless otherwise noted, $`.`$ denotes the max-norm. where $`V^{}`$ and $`Q^{}`$ are the optimal value functions of the base MDP.

First, consider Q-learning. Recall the update:

(4.1)
$$Q_{t+1}(x_t,a_t)=(1\alpha _t(x_t,a_t))Q_t(x_t,a_t)+\alpha _t(x_t,a_t)(r_t+\gamma \underset{a}{\mathrm{max}}Q_t(\stackrel{~}{y}_t,a)),$$

where $`\stackrel{~}{y}_t`$ is selected by sampling, i.e., according to the probability distribution $`P_t(x_t,a_t,.)`$.

###### Theorem 4.1.

Let $`Q^{}`$ be the optimal value function of the base MDP of the $`\epsilon `$-MDP, and let $`M=\mathrm{max}_{x,a}Q^{}(x,a)\mathrm{min}_{x,a}Q^{}(x,a)`$. If

1. every state-action pair $`(x,a)`$ is updated infinitely often,
2. the learning rates satisfy $`_{t=0}^{\mathrm{}}\chi (x_t=x,a_t=a)\alpha _t(x,a)=\mathrm{}`$ and $`_{t=0}^{\mathrm{}}\chi (x_t=x,a_t=a)\alpha _t(x,a)^2<\mathrm{}`$ uniformly w.p.1,333$`\chi `$ denotes the indicator function, $`\chi (\text{condition})`$ is 1 if the condition is true, and 0 otherwise.

then $`lim\; sup_t\mathrm{}Q_tQ^{}\frac{2}{1\gamma }\gamma M\epsilon `$ w.p.1.

The proof can be found in \[Szita et al., 2002a\].

The case of *SARSA.* The proof of the near-optimality of the SARSA update is similar to the proof of Q-learning. Here it suffices to refer to \[Singh et al., 2000\]: It has been proven for MDPs in \[Singh et al., 2000\] that if Q-learning converges to the optimal value function, then – under the same conditions – SARSA converges, too. The proof in \[Singh et al., 2000\] carries over directly to $`\epsilon `$-MDPs \[Szita et al., 2002a\].

### 4.3. pRL in the $`\epsilon `$-MDP Framework

Recall that for finding macro actions, we have applied learning methods that find optimal policies and value functions *in the modified environment*, where almost sure transitions (with probabilities greater than $`\kappa `$) are treated as sure (probability 1). It may be worth noting that according to Eq. 3.2 the value function of the model (i.e., $`\widehat{V}`$) is inherited from the original value function. $`\widehat{V}(x)`$ can be different from $`V(x)`$ iff $`P(x,y)`$ is modified by the model. Note that modification of $`P(x,y)`$ is at most $`\epsilon `$. Therefore the modified environment is an $`\epsilon `$-MDP of the original one, with $`\epsilon =1\kappa `$, so one can apply the $`\epsilon `$-MDP-theory to prove the following:

###### Corollary 4.2.

The pRL algorithm that uses either DP, Q-learning or SARSA for the learning of the macro-value function, produces approximations $`\widehat{V}_t`$ such that $`limsup_t\mathrm{}\widehat{V}_tV^{}\text{const}(1\kappa )`$.

Consequently, the macro actions found by pRL are asymptotically near-optimal.

## 5. Experiments

To test pRL, a $`40\times 40`$ maze was generated. The agent starts in the upper left corner and its goal is in the lower right corner. In every state, the agent observes its position and can take four different actions (N(=north),S(=south),E(=east), and W(=west)). In every state $`s`$, an action is successful with probability $`P^{succ}(s)`$ if the attempted direction is successfully executed, and it fails with probability $`1P^{succ}(s)`$ when the agent steps to a random, wrong direction. A lower bound for $`P^{succ}(s)`$ was set to 0.7. Regions with higher $`P^{succ}(s)`$ were generated. The probability of $`P^{succ}(s)`$ was gradually increased to 1 within these regions. The agent received a small $`(0.1)`$ negative reward in each state, except for the goal state, where it received the reward of $`+200`$. In addition, some pitfall domains were generated randomly. Every domain contributed to the local reward by $`1`$. In case of overlaps amongst the domains rewards were cumulated. After reaching the goal, a new episode was started. An example problem is shown in Fig. 2.

In the experiments, SARSA was utilized with $`\epsilon `$-greedy action selection, with eligibility traces (\[Sutton, 1988, Singh and Sutton, 1996\]). The following parameters were used: learning rate was held constant at $`\alpha =0.001`$. The eligibility decay $`\lambda `$ was set to $`0.95`$, discount factor $`\gamma `$ was equal to $`0.98`$, the probability of random action selection was $`0.1`$. For pRL, the same setting was used and several $`\kappa `$ values were tried. The number of updates was set to 10.

### 5.1. Convergence speed

In theory, by setting $`\kappa =1`$, pRL becomes identical to SARSA, and by setting $`\kappa =0`$, pRL becomes equivalent to the simplest Dyna algorithm (see 3.1). By choosing an appropriate $`\kappa `$, pRL learns at least as quickly as the better candidate, thus it approximates the convergence rate of Dyna – which usually converges in a few trials in this case – even at low $`\kappa `$ values (Fig. 3). Problems may arise if it is not known whether the best solution corresponds to total planning (Dyna) or if no planning is possible at all (SARSA). By choosing an appropriate $`\kappa `$, (or a set of $`\kappa `$’s) the convergence rate of pRL can always achieve these boundaries. The method is most useful if computations for different $`\kappa `$’s can be afforded.

### 5.2. Optimality

The second experiment demonstrates the near-optimality of pRL. The performance of the resulting policy for different $`\kappa `$ values were computed (Fig. 4). The question is the value of the constant multiplier of Corollary 4.2; whether it is too high or not. Here, the assumption about deterministic transitions does not significantly influence the performance for high $`\kappa `$ values, as expected. It was also found that the resulting policies perform quite well for some fairly low $`\kappa `$ values; for example, $`\kappa 0.5`$ performed still reasonably well. This suggests that better (stronger) estimations might exist for this problem. However, there is a $`\kappa `$ region, which exhibits poor performance. This region is around $`\kappa =0.7`$. The poor performance is the consequence of the special properties of our toy problem. In this $`\kappa `$ region, learning is compromised by the fact that most transition probabilities had a value of 0.7 (except in the plannable domains), and thus the algorithm had large uncertainties whether a particular transition is plannable or not.

## 6. Discussion

In this work, we introduced a new algorithm called pRL, which integrates planning into reinforcement learning in a novel way. An attractive property of our algorithm is that it can find near-optimal value function and near-optimal policy under certain conditions. Furthermore, the macros created by pRL are near-optimal. Near-optimality is controlled by the certainty of planning: our algorithm deals with $`kappa`$-plannable transitions, which have transition probabilities $`0<\kappa 1`$.

Our algorithm was illustrated by computer simulations on a toy problem. It was found that SARSA is slow, but it leads to optimal solutions. On the other hand, Dyna may be fast but the resulting policies could be poor. Our algorithm incorporates both SARSA and Dyna. These are the extremes, SARSA corresponds to $`\kappa =1`$, whereas the original form of Dyna \[Sutton, 1991\] is recovered when $`\kappa =0`$. By tuning $`\kappa `$ in pRL and by using the arising computationally inexpensive model of pRL, one can quickly find almost optimal solutions. This feature is most advantageous when off-line computations are inexpensive. Comparisons amongst the different methods are provided in Table 1.

pRL fits well into existing planning paradigms: it can be seen as a prioritized sweeping method where only the plannable states (i.e. states which can be reached almost deterministically from the current state) are updated. In turn, our method is a particular (extended) form of prioritized sweeping (see, e.g. \[Moore and Atkeson, 1993\]). One can decide to select updates according to state values \[Moore and Atkeson, 1993\], to select updates ordered by increasing distance from the current state or by prediction difference \[Peng and Williams, 1993\]. Prioritized sweeping can dramatically improve the performance of planning. Using our method, one cane also have performance bounds.

pRL can support approaches aiming to solve partially observable Markov decision processes (POMDPs). Most practical environments are inevitably POMDPs, which appear to the Markovian agent in the form of highly non-deterministic transitions. POMDPs are generally solved by extending the state representation of the agent (e.g. by estimating the hidden state variables). pRL can be used as a ‘pre-filtering’ technique. One can decide to use pRL first and to decouple the almost deterministic domains of the state space, where – as a first approximation – no further extension of the representation is necessary.

Another promising area for pRL is seen in the novel *E-learning* method \[Lőrincz et al., 2002, Szita et al., 2002b, Szita et al., 2002a\]. E-learning is a natural candidate for pRL when trajectory planning is desirable. E-learning is a modification of traditional reinforcement schemes. In E-learning, the selection of an ‘action’ is equivalent to selecting a desired next state \[Lőrincz et al., 2002\]. After this selection is made, the desired state is passed to an underlying controller. E-learning views the controller (the policy, or an inverse dynamics, or both) as part of the environment. pRL using E-learning develops estimations of state-state transition probabilities. In turn, the set of plannable states is directly available in E-learning. Thus pRL is greatly simplified within the framework of E-learning.

In our simulations, the discount rate for the planning ($`\gamma ^{}`$) was equal to the discount rate used in SARSA ($`\gamma `$). A simple but useful extension arises if one allows for different discount rates. This could be used for a number of purposes. If $`\gamma ^{}\gamma `$ is our choice, then planning may cover larger domains. On the other hand, $`\gamma ^{}\gamma `$ means that planned actions will achieve larger short-term rewards. Thus one can tune the interaction of planning and reinforcement learning by setting the discount rates, and determine whether to use planning or direct reinforcement learning to achieve short-term and long-term goals.

## 7. Acknowledgments.

This work was supported by the Hungarian National Science Foundation (Grant OTKA 32487) and by EOARD (Grant F61775-00-WE065). Any opinions, findings and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the European Office of Aerospace Research and Development, Air Force Office of Scientific Research, Air Force Research Laboratory.


format: text

text:
# Optimization and Physics: On the satisfiability of random Boolean formulae

(March 18, 2024)

## Abstract

LECTURE GIVEN AT TH2002: Given a set of Boolean variables, and some constraints between them, is it possible to find a configuration of the variables which satisfies all constraints? This problem, which is at the heart of combinatorial optimization and computational complexity theory, is used as a guide to show the convergence between these fields and the statistical physics of disordered systems. New results on satisfiability, both on the theoretical and practical side, can be obtained thanks to the use of physics concepts and methods.

Combinatorial optimization aims at finding, in a finite set of possible configurations, the one which minimizes a certain cost function. The famous example of the traveling salesman problem (TSP) can serve as an illustration: A salesman must make a tour through $`N`$ cities, visiting each city only once, and coming back at its starting point. The cost function is a symmetric matrix $`c_{ij}`$, where $`c_{ij}`$ is the cost for the travel between cities $`i`$ and $`j`$. A permutation $`\pi `$ of the $`N`$ cities gives a tour $`\pi (1)\pi (2)\pi (3)\mathrm{}`$. Taking into account the equivalence between various starting points and the direction of the tour, one sees that the number of distinct tour is $`(N1)!/2`$. For each tour $`\pi `$, the total cost is $`C=c_{\pi (N)\pi (1)}+_{r=1}^{N1}c_{\pi (r)\pi (r+1)}`$, which can be computed in $`N`$ operations. The problem is to find the tour $`\pi `$ with lowest cost.

As can be seen on this example, the basic ingredients of the optimization problems which will interest us are the following:

* An integer $`N`$ giving the size of the problem (in the TSP, it is the number of cities).
* A set of configurations, which typically scales like $`\mathrm{exp}(N^\alpha )`$.
* A cost function, which one can compute in polynomial time $`O(N^\gamma )`$.

Let me mention a few examples, beside the TSP papad .

In the assignment problem, one is given a set of $`N`$ persons $`i=1,\mathrm{},N`$, a set of $`N`$ tasks $`a=1,\mathrm{},N`$, a $`N\times N`$ cost matrix $`c`$ where $`c_{ia}`$ is the cost for having task $`a`$ performed by person $`i`$. An assignment is a permutation $`\pi S_N`$ assigning task $`a=\pi (i)`$ to person $`i`$, and the problem is to find the lowest cost assignment, i.e. the permutation which minimizes $`C=_{i=1}^Nc_{i\pi (i)}`$.

In the spin glass problemMPV , one is given a set of $`N`$ spins, which could be for instance Ising variables $`\sigma _i\{\pm 1\}`$, the total energy of the configuration is a sum of exchange interaction energies between all pairs of spins, $`E(\{\sigma \})=_{1i<jN}J_{ij}\sigma _i\sigma _j`$, and one seeks the ground state of the problem, i.e. the configuration (among the $`2^N`$ possible ones) which minimizes the energy.

In physical terms, optimization problems consist in finding ground states. This task can be non trivial if a system is frustrated, which means that one cannot get the ground state simply by minimizing the energy locally. This is typically what happens in a spin glass. In some sense, statistical physics addresses a more general question. Assuming that the system is at thermal equilibrium at temperature $`T`$, every configuration $`𝒞`$ is assigned a Boltzmann probability, $`P(𝒞)=\frac{1}{Z}e^{\beta E(𝒞)}`$. Beside finding the ground state, one can ask also interesting questions about which are the typical configurations at a given temperature, like counting them (which leads to the introduction of entropy), or trying to know if they are located in one single region of phase space, or if they build up well separated clusters, as often happens in situation where ergodicity is broken. The generalization introduced by using a finite temperature, beside leading to interesting questions, can also be useful for optimization, both from the algorithmic point of view (for instance this is the essence of the simulated annealing algorithmKGV , which is a general purpose heuristic algorithm that can be used in many optimization problems), but also from an analytic point of view MPV . Conversely, smart optimization algorithms turn out to be very useful in the study of frustrated physical systems like spin glasses or random field models, and the cross-fertilization between these two fields (and also with the related domain of error correcting codes for information transmission codes ) has been very fruitful in the recent years TCS\_issue .

Before proceeding with one such example, let us briefly mention a few important results in optimization which will provide the necessary background and motivation. One of the great achievements of computer science is the theory of computational complexity. It is impossible to present it in any details here and I will just sketch a few main ideas, the interested reader can study it for instance in papad2 .

Within the general framework explained above, we can define three types of optimization problems: the optimization problem in which one wants to find the optimal configuration, the evaluation problem in which one wants to compute the optimal cost (i.e. the ground state energy), the decision problem in which one wants to know, given a threshold cost $`C_0`$, if there exists a configuration of cost less than $`C_0`$.

One classification of decision problems is based on the scaling with $`N`$ of the computer time required to solve them in the worst case. There are two main complexity classes:

* Class P = polynomial problems: they can be solved in a time $`t<N^\alpha `$. The assignment is an example of a polynomial problem, as is the spin glass problem in 2 dimensions.
* Class NP = non-deterministic polynomial: Given a ’yes’ solution to a NP problem, it can be checked in polynomial time. Roughly speaking this means that the energy is computable in polynomial time, so this class contains a wide variety of problems, including most of the ones of interest in physics. All problems mentioned above are in NP.

One nice aspect of focusing on polynomiality is that it allows to forget about the details of the definition of $`N`$, the implementation, language, machine, etc…: any ’reasonable’ such change (for instance one could have used the number of possible links in the assignment) will change the exponent of $`N`$ appearing in the computer time of a problem in P, but not transform it into an exponential behavior. A problem $`A`$ is said to be at least as hard as a problem $`B`$ if there exists a polynomial algorithm which transforms any instance of $`B`$ into an instance of $`A`$.

This allows to define the very important class:

* NP-complete A problem is NP complete if it is at least as hard as any problem in NP.

So the NP-complete are the hardest problems in NP. If one such problem can be solved in polynomial time, then all the problems in NP are solved in polynomial time. Clearly P is contained in NP, but it is not yet known whether P = NP , and this is considered as a major challenge.

A great result was obtained in this field by Cook in 1971 Cook : The satisfiability problem, which we shall describe below, is NP-complete. Since then hundreds of problems have been shown to belong to this class, among which the decision TSP or the spin glass in dimension larger or equal to 3.

The fact that 3-d spin glass is NP-complete while 2-d spin glass is P might induce one to infer that NP-completeness is equivalent with the existence of a glass transition. This reasoning is too naive and wrong; the reason is that the complexity classification relies on a worth-case analysis, while physicists study the behavior of a typical sample. The development of a typical case complexity theory has become a major goal average , also motivated by the experimental observations that many instances of NP-complete problems are easy to solveTCS\_issue .

One way of addressing this issue of a typical case complexity is to define a probability measure on the instances (= the ’samples’) of the optimization problem which one is considering. Typical examples are:

* TSP with independent random points uniformly distributed in the unit square
* assignment with independent affinities uniformly distributed on $`[0,1]`$
* CuMn spin glass at one percent Mn

Once this measure has been defined, one is interested in the properties of the generic sample in the $`N\mathrm{}`$ limit. In most cases, global properties (e.g. extensive thermodynamic quantities, among which the ground state energy), turn out to be self averaging. This means for instance that the distribution of the ground state energy density becomes more and more peaked around an asymptotic value in the large $`N`$ limit: almost all samples have the same ground state energy density. In this situation, a statistical physics analysis is appropriate. Early examples of the use of statistical physics in such a context are the derivation of bounds for the optimal length of a TSPMVan , the exact prediction of the ground state energy in the random assignment problem defined above, where the result $`E_0=\frac{\pi ^2}{6}`$, derived in 1985 through a replica analysisMP\_match , was recently confirmed rigorously by Aldous Aldous , or the link between spin glasses and graph partitioning FuAnd .

As statistical physics can be quite powerful at understanding the generic structure of an optimization problem, one may also hope that it can help finding better optimization algorithms. A successful example which was developed recently is the satisfiability problem, to which we now turn.

As we have seen, satisfiability is a core problem in optimization and complexity theory. It is defined as follows Hayes : A configuration is a set of $`N`$ Boolean variables $`x_i\{0,1\}i=1,\mathrm{},N`$. One is given $`M`$ constraints which are clauses, meaning that they are in the form of an OR function of some variables or their negations. For instance: $`x_1x_{27}\overline{x}_3`$, $`\overline{x}_{11}x_2`$, are clauses (notation: $`\overline{x}_3`$ is the negation of $`x_3`$). So the clause $`x_1x_{27}\overline{x}_3`$ is satisfied if either $`x_1=1`$, or $`x_{27}=1`$, or $`x_3=0`$ (these events do not exclude each other). The satisfiability problem is a decision problem. It asks whether there exists a choice of the Boolean variables such that all constraints are satisfied (we will call it a SAT configuration). This is a very generic problem, because one can see it as finding a configuration of the $`N`$ Boolean variables which satisfies the logical proposition built from the AND of all the clauses (in our example $`(x_1x_{27}\overline{x}_3)(\overline{x}_{11}x_2)\mathrm{}`$), and any logical proposition can be written in such a ’conjunctive normal form’ .

Satisfiability is known to be NP complete if it contains clauses of length $`3`$, but common sense and experience show that the problem can often be easy; for instance if the number of constraints per variable $`\alpha =\frac{M}{N}`$ is small, the problem is often SAT, if it is large, the problem is often UNSAT.

This behavior can be characterized quantitatively by looking at the typical complexity of the random 3-SAT problem, defined as follows. Each clause involves exactly three variables, chosen randomly in $`\{x_1,..,x_N\}`$; a variable appearing in the clause is negated randomly with probability $`1/2`$. This defines the probability measure on instances for the random 3-SAT problem. The control parameter is the ratio Constraints/Variables $`\alpha =\frac{M}{N}`$ .

The properties of random 3-SAT have been first investigated numerically KirkSel ; Crawford , and exhibit a very interesting threshold phenomenon at $`\alpha _c4.26`$: a randomly chosen sample is generically SAT for $`\alpha <\alpha _c`$ (meaning that it is SAT with probability $`1`$ when $`N\mathrm{}`$), generically UNSAT for $`\alpha >\alpha _c`$. The time used by the best available algorithms (which have an exponential complexity) to study the problem also displays an interesting behavior: For $`\alpha `$ well below $`\alpha _c`$, it is easy to find a SAT configuration; for $`\alpha `$ well above $`\alpha _c`$, it is relatively easy to find a contradiction in the constraints, showing that the problem is UNSAT. The really difficult region is the intermediate one where $`\alpha \alpha _c`$, where the computer time requested to solve the problem is much larger and increases very fast with system size. A lot of important work has been done on this problem, to establish the existence of a threshold phenomenon, give upper and lower bounds on $`\alpha _c`$, show the existence of finite size effects around $`\alpha \alpha _c`$ with scaling exponents. We refer the reader to the literature KirkSel ; Crawford ; MZKST ; bounds\_1 ; bounds\_2 ; bounds\_3 ; bounds\_4 , and just here extract a few crucial aspects for our discussion. The threshold phenomenon is a phase transition, and the neighborhood of the transition is the region where the algorithmic problem becomes really hard.

This relationship between phase transition and complexity makes a statistical physics analysis of this problem particularly interesting. Monasson and Zecchina have been the first ones to recognize this importance and to use statistical physics methods for an analytic study of the random 3-SAT problem MoZe\_prl ; MoZe\_pre . Basically this problem falls into the broad category of spin glass physics. This is immediately seen through the following formulation. To make things look more familiar, physicists like to introduce for each Boolean variable $`x_i\{0,1\}`$ an Ising spin $`\sigma _i\{1,1\}`$ through $`x_i=\frac{1+\sigma _i}{2}`$. A satisfiability problem like

$$(x_1x_{27}\overline{x}_3)(\overline{x}_{11}x_3x_2)\mathrm{}$$
(1)
can be written in terms of an energy function, where the energy is equal to one for each violated clause. Explicitly, in the above example, one would have

$$E=\frac{1+\sigma _1}{2}\frac{1+\sigma _{27}}{2}\frac{1\sigma _3}{2}+\frac{1\sigma _{11}}{2}\frac{1+\sigma _3}{2}\frac{1+\sigma _2}{2}+\mathrm{}$$
(2)
This is clearly a problem of interacting spins, with 1,2, and 3 spin interactions, disorder (in the choice of which variable interacts with which), and frustration (because some constraints are antagonist). More technically, the problem has a special type of three-spin interactions on a random hyper-graph.

Using the replica method, Monasson and Zecchina first showed the existence of a phase transition within the replica symmetric approximation, at $`\alpha _c=5.18`$, then showed that replica symmetry must be broken in this problem. Some variational approximation to describe the replica symmetry breaking effects were developed in particular in BMW ; FLRZ .

Recently, in a collaboration with G.Parisi and R. Zecchina MEPAZE ; MZ\_pre , we have developed a new approach to the statistical physics of the random 3-SAT problem using the cavity method. While the cavity method had been originally invented to deal with the SK model where the interactions are of infinite range MPV , it was later adapted to problems with ’finite connectivity’, in which a given variable interacts with a finite set of other variables. While this is easily done for systems which are replica symmetric (like the assignment, or the random TSP with independent links), it turned out to be considerably more difficult to develop the corresponding formalism and turn it into a practical method, in the case where replica symmetry is broken. This has been done in the last two years in joint works with G.Parisi MP\_Bethe , and has opened the road to the study of finite connectivity optimization problems with replica symmetry breaking like random K-sat. Curiously, while the cavity method is in principle equivalent to the replica method (although it proceeds through direct probabilistic analysis instead of using an analytic continuation in the number of replicas), it turns out that it is easier to solve this problem with the cavity method.

The analytic study of MEPAZE ; MZ\_pre  for the random 3-SAT problem at zero temperature shows the following phase diagram (see fig.1):

* For $`\alpha <\alpha _d=3.921`$, the problem is generically SAT; the solution can be found relatively easily, because the space of SAT configurations builds up a single big connected cluster. A $`T=0`$ Metropolis algorithm, in which one proposes to flip a randomly chosen variable, and accepts the change iff the number of violated constraints in the new configuration is less or equal to the old one, is able to find a SAT configuration. We call this the EASY-SAT phase
* For $`\alpha _d<\alpha <\alpha _c=4.267`$, the problem is still generically SAT, but now it becomes very difficult to find a solution (we call this the HARD-SAT phase). The reason is that the configurations of variables which satisfy all constraints build up some clusters. Each such cluster, which we call a ’state’, is well separated from the other clusters (passing from one to the other requires flipping an extensive number of variables). But there also exist many “metastable states”: starting from a random configuration, a local descent algorithm will get trapped in some cluster of configurations with a given finite energy (they all have the same number of violated clauses, and it is impossible to get out of this cluster towards lower energy configurations through any descending sequence of one spin flip moves). The number of SAT clusters $`𝒩`$ is exponentially large in $`N`$, it behaves as $`𝒩\mathrm{exp}(N\mathrm{\Sigma })`$; but the number of metastable clusters is also exponentially large with a larger growth rate, behaving like $`\mathrm{exp}(N\mathrm{\Sigma }_{ms})`$ with $`\mathrm{\Sigma }_{ms}>\mathrm{\Sigma }`$. The most numerous metastable clusters, which have an energy density $`e_{th}`$, will trap all local descent algorithms (zero temperature Metropolis of course, but also simulated annealing, unless it is run for an exponential time).
* For $`\alpha >\alpha _c`$, the problem is typically UNSAT. The ground state energy density $`e_0`$ is positive. Finding a configuration with lowest energy is also very difficult because of the proliferation of metastable states.

A more quantitative description of the thermodynamic quantities in the various phases is shown in fig.2. The most striking result is the existence of an intermediate SAT phase where clustering occurs. A hint of such a behavior had been first found in a sophisticated variational one step replica symmetry breaking approximation developed in BMW ; however this approximation predicted a second order phase transition (clusters separating continuously), while we now think that the transition is discontinuous: an exponentially large number of macroscopically separated clusters appears discontinuously at $`\alpha =\alpha _d`$. Another point which should be noticed is the fact that the complexity, and the energy $`e_{th}`$ in the HARD-SAT phase are rather small: around $`e_{th}\mathrm{3\; 10}^3`$ violated clause per variable for $`\alpha =4.2`$. This shows that until one reaches problems with at least a few thousands variables, one cannot feel the true complexity of the problem. This can explain why the existence of the intermediate phase went unnoticed in previous simulations.

The second type of results found in MEPAZE ; MZ\_pre ; BMZ  is a new class of algorithms dealing with the many states structure. Basically this algorithm amounts to using the cavity equations on one given sample. Originally the cavity method was developed to handle a statistical distribution of samples. For instance in the random 3-SAT case its basic strategy is to add one extra variable and connect it randomly to a number of new clauses which has the correct statistical distribution. In the large $`N`$ limit, the statistics of the local field on the new variable should be identical to the statistics of the local fields on any other variable in the absence of the new one. It turns out that this strategy can be adapted to study a single sample: one considers a given clause $`a`$, which involves three variables $`\sigma _1,\sigma _2,\sigma _3`$. The cavity field on $`\sigma _1`$ is the field felt by $`\sigma _1`$ in the absence of $`a`$. In the case where there exist many states, there is one such field for each possible state of the system, and the order parameter is the survey of these fields, in all the states of fixed energy density $`e`$:

$$P_0^e(h)=C^t\underset{\alpha }{}\delta \left(h_0^\alpha h\right)\delta \left(\frac{E^\alpha }{N}e\right)$$
(3)
One can then write a recursion recursion between these surveys. Looking for instance at the structure of fig.3, one gets the following iteration equation:

$`P_0^e(h)`$ $`=`$ $`C^t{\displaystyle P_{\sigma _1}^e(g_1)𝑑g_1P_{\tau _1}^e(h_1)𝑑h_1P_{\sigma _2}^e(g_2)𝑑g_2P_{\tau _2}^e(h_2)𝑑h_2}`$ (4)
$`\delta (hf(g_1,h_1,g_2,h_2))\mathrm{exp}\left({\displaystyle \frac{d\mathrm{\Sigma }}{de}}w(g_1,h_1,g_2,h_2)\right)`$
The function $`f`$ just computes the value of the new cavity field on $`\sigma _0`$ in terms of the four cavity fields $`g_1,h_1,g_2,h_2`$. It is easily computed by considering the statistical mechanics problem of the five-spin system $`\{\sigma _0,\sigma _1,\tau _1,\sigma _2,\tau _2\}`$ and summing over the 16 possible states of the spins $`\sigma _1,\tau _1,\sigma _2,\tau _2`$. The function $`w`$ computes the free energy shift induced by the addition of $`\sigma _0`$ to the system with the four spins $`\sigma _1,\tau _1,\sigma _2,\tau _2`$. The exponential reweighting term in (4) is the crucial piece of survey propagation: it appears because one considers the survey at a fixed energy density $`e`$. As the number of states at energy $`E=Ne+\delta E`$ increases in $`\mathrm{exp}(y\delta E)`$, where $`y=\frac{d\mathrm{\Sigma }}{de}`$, this favors the states with a large negative energy shift.

The usual cavity method for the random 3-SAT problem determines the probability, when a new variable is added at random, that its survey $`P_0^e(h)`$ is a given function $`P(h)`$: the order parameter is thus a functional, the probability of a probability. Because the fields are distributed on integers, this object is not as terrible as it may look and it is possible to solve the equation and compute the ’complexity curve’ $`\mathrm{\Sigma }(E)`$, giving the phase diagram described above.

The algorithm for one given sample basically iterates the survey propagation equation on one given graph. It is a message passing algorithm which can be seen as a generalization of the belief propagation algorithm familiar to computer scientistspearl . The belief propagation is a propagation of local magnetic fields, which is equivalent to using a Bethe approximation yed . Unfortunately, it does not converge in the hard-SAT region because various subparts of the graph tend to settle in distinct states, and there is no way to globally choose a state. In this region, the survey propagation, which propagates the information on the whole set of states (in the form of a histogram), does converge.

Based on the surveys, one can detect some strongly biased spins, which are fixed to one in almost all SAT configurations. The “Survey Inspired Decimation” (SID) algorithm fixes the spin which is most biased, then it re-runs the survey propagation on the smaller sample so obtained, and then iterates… An example of the evolution of the complexity as a function of the decimation is shown in 4. This algorithm has been tested in the hard SAT phase. It easily solves the ’large’ benchmarks of random 3sat at $`\alpha =4.2`$ with $`N=1000,2000`$ available at SATLIB . It turns out to be able to solve typical random 3-SAT problems with up to $`N=10^7`$ at $`\alpha =4.2`$ in a few hours on a PC, which makes it much better than available algorithms. The main point is that the set of surveys contains a lot of detailed information on a given sample and can probably be used to find many new algorithms, of which SID is one example.

To summarize, the recent statistical physics approaches to the random 3-SAT problem give the following results:

* An analytic result for the phase diagram of the generic samples
* An explanation for the slowdown of algorithms near $`\alpha _c=4.267`$: this is due to the existence of a HARD-SAT phase at $`\alpha [3.921,4.267]`$, with exponentially many metastable states.
* An algorithm for single sample analysisweb : Survey propagation converges and yields very non trivial information on the sample. It can be used for instance to decimate the problem and get an efficient SAT-solver in the hard-SAT region.

This whole set of results calls for a lot of developments in many directions.

On the analytical side, the cavity method results quoted above are believed to be correct, but they are not proven rigorously. It would be very interesting to turn these computations into a rigorous proof. A very interesting step in this direction was taken recently by Franz and Leone who showed that the result for the critical threshold $`\alpha _c`$ obtained by the cavity method on random K-SAT with $`K`$ even actually give a rigorous upper bound to the critical $`\alpha _c`$ FraLeo . The whole construction of the cavity method with the clustering structure, the many states and the resulting reweighting, has actually been checked versus rigorous computations on a variant of the SAT problem, the random XORSAT problem, where rigorous computations xorsat  have confirmed the validity of the approach.

On the numerical side, one needs to develop convergence proofs for survey propagation, and to derive the generalization of the algorithm for the case in which there exists local structures in the interaction graphs. This will amount to generalizing from a Bethe like approximation (with many states) to a cluster variational method with larger clusters (and with many states).

The techniques can also be extended to other optimization problems like coloring colouring . Beside the problems mentioned here, there exist many other fascinating problems of joint interests to physicists and computer scientists, like e.g the dynamical behavior of algorithms in optimization or error correcting codes, which I cannot survey here. Let me just quote a few recent references to help the readers through the corresponding bibliography cocmon ; global\_algorithm ; FLMR ; MonZec .

Acknowledgements It is a pleasure to thank all my collaborators over the years on various topics at the interface between statistical physics and optimization. I am particularly thankful to G. Parisi and R. Zecchina for the wonderful collaborations over the last three years which lead to the works on the cavity method and the satisfiability problem.


text:
# QCD Saturation and Deuteron-Nucleus Collisions

## I Erratum added on March 8, 2004

After the publication of our paper, the experimental data on hadron multiplicities in $`dAu`$ collisions at RHIC have been presented phobos-e ; brahms-e ; star-e . While the predicted multiplicity distributions are consistent with the data within $`30\%`$, a detailed comparison exhibits discrepancies both in the shape of rapidity distributions and in centrality dependence.

Should we consider this discrepancy as an argument against the Color Glass Condensate approach, or does it stem from something that we overlooked in our calculations? Our analysis presented below shows that the disagreement with the data originates mainly in the centrality determination procedure. In our paper, as well as in the previous publications dedicated to $`AuAu`$ collisions, we used the optical Glauber model to evaluate the numbers of participants $`N_{part}`$ and collisions $`N_{coll}`$ corresponding to a particular centrality cut imposed on the multiplicity distribution. Meanwhile, the experiments use Monte Carlo realizations of the Glauber model, which also take into account the geometry and acceptance of the detector.

Already for the $`AuAu`$ system it has been noted that in peripheral collisions, when $`N_{part}`$ and $`N_{coll}`$ become small, the discrepancy between the optical and Monte Carlo realizations becomes significant. Nevertheless, for $`N_{part}70`$ the agreement between the two approaches was reasonable, on the order of $`10\%`$. However in $`dAu`$ collisions, where the numbers of participants and collisions are always relatively small, the differences in centrality determination strongly affect the comparison of our predictions to the data at all centralities: indeed, the discrepancies between the number of participants given in our Table 1 and presented by the experimental Collaborations phobos-e ; brahms-e ; star-e  sometimes are as big as $`50\%`$. Part of the problem stems from the fact that in the optical approach in peripheral collisions $`N_{part}`$ can be smaller than two – this is because the overlap integral in this case has a meaning of the probability to have $`N_{part}=2`$; in the Monte Carlo approach, one instead triggers on the inelastic interaction event, so the number of participants is always $`N_{part}2`$. Note that since the shape of our rapidity distributions (see Eq.(12)) depends on the number and density of participants in $`Au`$ and deuterium separately, it is also affected by the uncertainty in centrality determination.

To investigate the influence of these differences in centrality determination, we have repeated the calculations according to our Eq.(12), but using experimental numbers of $`N_{part}`$ and $`N_{coll}`$. This can be done by multiplying our $`\rho _{part}^D`$ and $`\rho _{part}^A`$ by the factor $`\frac{N_{part}^{exp}}{N_{part}^{exp,MB}}\frac{N_{part}^{our,MB}}{N_{part}^{our}}`$. We found that this alone removes almost all of the discrepancy with the data. Moreover, as we discussed in the paper, our master equation Eq.(12) is based on the assumption that the gluon distribution inside the proton can be effectively described as a saturated one. As we emphasized, the value of the effective proton saturation momentum is somewhat uncertain, and we used $`Q_s^2(p;W=200\mathrm{GeV};y)=\mathrm{\hspace{0.17em}0.25}\mathrm{GeV}^2`$. We find that a larger value $`Q_s^2(p;W=200\mathrm{GeV};y)=\mathrm{\hspace{0.17em}0.34}\mathrm{GeV}^2`$ improves the fit. The results are shown in Fig.1 by the dashed lines. One can see that the agreement is now quite good, with an exception of the gold fragmentation region. The Color Glass Condensate approach cannot be justified in this region where the nuclear distribution is probed at large Bjorken $`x`$; the fragmentation of the $`Au`$ participants is expected to dominate there. If we assume that $`dN/dy`$ is equal to $`N_{part}^{Au}\times dN(pp)/dy`$ in the region of $`5<y<3`$, we get the results shown by the solid curves in Fig.1.


text:
# Nonlinear ac response of colloidal suspension with an intrinsic dispersion

## I Introduction

When an intense electric field is applied to a suspension containing highly polarized particled embedded in a host fluid, the particles will form chains due to the induced dipole moments inside the particles, and hence complex anisotropic structures occur. Also, the apparent viscosity of the suspension will be enhanced greatly. These structures can have anisotropic physical properties, such as the effective conductivity and permittivity . In this case, optical nonlinearity enhancement was found as well .

A convenient method of probing the nonlinear characteristics of the field-induced structures is to measure the harmonics of the nonlinear polarization under the application of a sinusoidal (ac) electric field . Therefore it is possible to perform a real-time monitoring of the field-induced aggregation process in such a suspension by measuring the optical nonlinear ac responses both parallel and perpendicular to the uniaxial anisotropic axis. Recently, an anisotropic Maxwell-Garnett theory  was applied to study the nonlinear polarization of these anisotropic structures . Anisotropic optical nonlinear ac response was indeed found.

When a nonlinear composite with nonlinear dielectric particles embedded in a host medium, or with a nonlinear host medium is subjected to a sinusoidal field, the electrical response in the composite will in general be a superposition of many sinusoidal functions . It is natural to investigate the effects of a nonlinear characteristics on the optical nonlinear ac response in a colloidal suspension which can be regarded as a nonlinear composite medium . The strength of the nonlinear polarization should be reflected in the magnitude of the harmonics. For these field-induced anisotropic structures, we found that as more and more chains are formed the harmonics of local field and induced dipole moment increase (decrease) for longitudinal (transverse) field case .

On the other hand, the effect of nonlinear characteristics on the inter-particle force has been analyzed in a colloidal suspension of nonlinear particles  and further extended to a nonlinear host medium .

Actually, the intrinsic dielectric dispersion often occurs due to the surface conductivity or the inhomogeneous structure (e.g., coating effect ) which often exist in colloidal suspensions. Hence it is of particular importance to discuss the effect of an intrinsic dielectric dispersion on the nonlinear polarization. In this regard, we will use a perturbation expansion method to calculate the nonlinear ac response of a nonlinear composite, i.e., harmonics of local field as well as induced dipole moment.

## II Nonlinear polarization and its higher harmonics

Let us consider the system in which the suspended particles of radius $`a`$ with a nonlinear dielectric constant are suspended in a host medium of linear dielectric constant $`ϵ_2`$. The nonlinear characteristic gives rise to a field-dependent dielectric coefficient . In this case, the electric displacement-electric field relation inside the spheres is given by:

$`𝐃_1`$ $`=`$ $`ϵ_1𝐄_1+\chi _1|𝐄_1|^2𝐄_1ϵ_1𝐄_1+\chi _1|𝐄_1|^2𝐄_1`$ (1)
$`=`$ $`(ϵ_{1H}+{\displaystyle \frac{\mathrm{\Delta }ϵ_1}{1+i\omega /\omega _0}})𝐄_1+\chi _1|𝐄_1|^2𝐄_1\stackrel{~}{ϵ}_1𝐄_1,`$ (2)
with $`i=\sqrt{1}`$, where $`ϵ_1`$ and $`\chi _1`$ are the linear coefficient and the weak nonlinear coefficient of the suspended particles, respectively. Here $`ϵ_1`$ is expressed as a Debye expression, $`ϵ_{1H}`$ is the dielectric constant at high frequency, $`\mathrm{\Delta }ϵ_1`$ is the intrinsic dispersion of the particle, $`\omega _0=2\pi /\tau _0`$ is the intrinsic characteristic frequency, with $`\tau _0`$ being the intrinsic relaxation time, and $`\omega `$ is the angular frequency of the external field. In Eq.(2), we have adopted an approximation: the local field inside the particles is assumed to be uniform. This assumption is called the decoupling approximation . It has been shown that such an approximation yields a lower bound for the accurate result for the local field . We further assumed that $`\chi _1`$ is independent of frequency, which is a valid assumption for low-frequency processes in colloidal suspensions. As a result, the induced dipole moment under an applied field $`𝐄(t)=E(t)\widehat{𝐳}`$ along $`z`$axis, with $`t`$ being time, is given by :

$$\stackrel{~}{𝐩}=\stackrel{~}{ϵ}_ea^3\stackrel{~}{b}𝐄(t),$$
(3)
where $`\stackrel{~}{ϵ}_e`$ is the effective dielectric constant of the suspension, and $`\stackrel{~}{b}`$ is the field-dependent dipole factor:

$$\stackrel{~}{b}=\frac{\stackrel{~}{ϵ_1}ϵ_2}{\stackrel{~}{ϵ_1}+2ϵ_2}=\frac{ϵ_1+\chi _1|𝐄_1|^2ϵ_2}{ϵ_1+\chi _1|𝐄_1|^2+2ϵ_2}.$$
(4)
In order to obtain the effective dielectric constant $`\stackrel{~}{ϵ}_e`$, we invoke the Maxwell-Garnett approximation (MGA) for anisotropic composites . For the longitudinal field case when the ac field is applied along the uniaxial anisotropic axis, the MGA has the form

$$\frac{\stackrel{~}{ϵ}_eϵ_2}{\beta _L\stackrel{~}{ϵ}_e+(3\beta _L)ϵ_2}=f\frac{\stackrel{~}{ϵ}_1ϵ_2}{\stackrel{~}{ϵ}_1+2ϵ_2}.$$
(5)
with $`f`$ being the volume fraction of particles, whereas for an transverse field case when the ac field is applied perpendicular to the uniaxial anisotropic axis, the MGA expression is obtained by replacing $`\beta _L`$ with $`\beta _T=(3\beta _L)/2`$ . Here $`\beta _L`$ and $`\beta _T`$ denote the local field factors parallel and perpendicular to the uniaxial anisotropic axis. These local field factors are defined as the ratio of the local field in the particles to the Lorentz cavity field . For isotropic composites, $`\beta _L=\beta _T=1`$, while both $`\beta _L`$ and $`\beta _T`$ will deviate from unity for an anisotropic distribution of particles in composites. These $`\beta `$ factors have been evaluated in a tetragonal lattice of dipole moments  and in various field-structured composites .

We resort to the spectral representation approach  to represent the local field inside the particle :

$$|𝐄_1|^2=\frac{𝐄^2(t)}{f}_0^1\frac{|s|^2\mu (x)}{|sx|^2}dx,$$
(6)
with the material parameter $`s=(1\stackrel{~}{ϵ}_1/ϵ_2)^1`$ as well as the spectral density function $`\mu (x)=f\delta [x(1f\beta )/3]`$, which satisfies the sum rule $`_0^1\mu (x)dx=f`$. Hence Eq.(6) has the form

$$|𝐄_1|^2=\frac{9|ϵ_2|^2}{|\stackrel{~}{ϵ}_1(1f\beta )+ϵ_2(2+f\beta )|^2}𝐄^2(t).$$
(7)
Actually, this is a self-consistent equation for $`|𝐄_1|^2`$, which can be solved at least numerically . In what follows, we will use a perturbation method instead to obtain analytic expressions for harmonics of the local electric field and the induced dipole moment.

On the other hand, due to the uniform local field in each particle, Eq.(7) can also be directly derived from :

$$|𝐄_1|^2=\frac{1}{f}𝐄^2(t)\left|\frac{\stackrel{~}{ϵ}_e}{\stackrel{~}{ϵ}_1}\right|.$$
(8)
If we apply a sinusoidal electric field, i.e. $`𝐄(t)=𝐄\mathrm{sin}(\omega t)`$, the induced dipole moment $`\stackrel{~}{p}`$ will depend on time sinusoidally, too. By virtue of the inversion symmetry, $`\stackrel{~}{p}`$ is a superposition of odd-order harmonics such that

$$\stackrel{~}{p}=p_\omega \mathrm{sin}\omega t+p_{3\omega }\mathrm{sin}3\omega t+\mathrm{}.$$
(9)
Also, the local electric field contains similar harmonics

$$\sqrt{|𝐄_1|^2}=E_\omega \mathrm{sin}\omega t+E_{3\omega }\mathrm{sin}3\omega t+\mathrm{}.$$
(10)
These harmonic coefficients can be extracted from the time dependence of the solution of $`\stackrel{~}{p}`$ and $`E_1(t)`$.

## III Analytic Solutions

In this section, we will apply the perturbation expansion  to extract the harmonics of the local electric field and the induced dipole moment. It is known that the perturbation expansion method is applicable to weak nonlinearity only, limited by the convergence of the series expansion.

We expand $`\stackrel{~}{p}`$ and $`\sqrt{\chi _1|𝐄_1|^2}`$ into a Taylor expansion, taking $`\chi _1|𝐄_1|^2`$ as a perturbative quantity:

$$\stackrel{~}{p}=\underset{n=0}{\overset{\mathrm{}}{}}\frac{a^3E(t)}{n!}\frac{^n}{\stackrel{~}{ϵ}_1^n}\left[ϵ_2\stackrel{~}{b}\frac{\stackrel{~}{ϵ}_1(1+3ff\beta )+ϵ_2(23f+f\beta )}{\stackrel{~}{ϵ}_1(1f\beta )+ϵ_2(2+f\beta )}\right]_{\stackrel{~}{ϵ}_1=ϵ_1}(\chi _1|𝐄_1|^2)^n,$$
(11)
$$\sqrt{\chi _1|𝐄_1|^2}=\underset{m=0}{\overset{\mathrm{}}{}}\frac{1}{m!}\frac{^m}{\stackrel{~}{ϵ}_1^m}\left[\frac{3|ϵ_2|\sqrt{\chi _1𝐄^2(t)}}{|(1f\beta )\stackrel{~}{ϵ}_1+(2+f\beta )ϵ_2|}\right]_{\stackrel{~}{ϵ}_1=ϵ_1}(\chi _1|𝐄_1|^2)^m.$$
(12)
In view of weak nonlinearity, we can rewrite Eqs.(10) and (11), keeping the lowest orders of $`\chi _1|𝐄_1|^2`$ and $`\chi _1𝐄^2(t)`$:

$`\stackrel{~}{p}`$ $``$ $`h_1E(t)+h_3\chi _1E^3(t),`$ (13)
$`\sqrt{\chi _1|𝐄_1^2|}`$ $``$ $`j_1\sqrt{\chi _1𝐄^2(t)}+j_3\left(\sqrt{\chi _1𝐄^2(t)}\right)^3,`$ (14)
where

$`h_1`$ $`=`$ $`a^3ϵ_2{\displaystyle \frac{ϵ_1ϵ_2}{ϵ_1+2ϵ_2}}{\displaystyle \frac{\rho +3f(ϵ_1ϵ_2)}{\rho }},`$ (15)
$`h_3`$ $`=`$ $`{\displaystyle \frac{27a^3|ϵ_2|^2ϵ_2^2}{\rho |\rho |^2(ϵ_1+2ϵ_2)}}\left[{\displaystyle \frac{3f(ϵ_1ϵ_2)}{\rho }}+{\displaystyle \frac{\rho +3f(ϵ_1ϵ_2)}{ϵ_1+2ϵ_2}}\right],`$ (16)
$`j_1`$ $`=`$ $`{\displaystyle \frac{3|ϵ_2|}{\rho }},j_3={\displaystyle \frac{27}{2}}|ϵ_2|^3{\displaystyle \frac{f\beta 1}{\rho |\rho |}}({\displaystyle \frac{1}{|\rho |^2}}+{\displaystyle \frac{1}{(\rho ^2)^{}}}),`$ (17)
with $`\rho =(1f\beta )ϵ_1+(2+f\beta )ϵ_2`$.

In the case of a sinusoidal field, we can expand $`E^3(t)`$ in terms of the first and the third harmonics. The comparison with Eq.(8) yields the harmonics of the induced dipole moment:

$$p_\omega =h_1E+\frac{3}{4}h_3\chi _1E^3,p_{3\omega }=\frac{1}{4}h_3\chi _1E^3.$$
(18)
Similarly, we find the harmonics of the local electric field

$$\chi _1^{1/2}E_\omega =j_1\chi _1^{1/2}E+\frac{3}{4}j_3(\chi _1^{1/2}E)^3,\chi _1^{1/2}E_{3\omega }=\frac{1}{4}j_3(\chi _1^{1/2}E)^3.$$
(19)
In the above analysis, we have used the identity $`\mathrm{sin}^3\omega t=(3/4)\mathrm{sin}\omega t(1/4)\mathrm{sin}3\omega t`$ to obtain the first and the third harmonics. Similar analysis can be used to extract the higher-order harmonics, by retaining more terms in the series expansion. For instance, to obtain fifth harmonics, we may take into account the identity $`\mathrm{sin}^5\omega t=(5/8)\mathrm{sin}\omega t(5/16)\mathrm{sin}3\omega t+(1/16)\mathrm{sin}5\omega t`$, and then obtain the harmonics of induced dipole moment and local electric field, respectively:

$`p_\omega `$ $`=`$ $`h_1E+{\displaystyle \frac{3}{4}}h_3\chi _1E^3+{\displaystyle \frac{5}{8}}h_5\chi _1^2E^5,`$ (20)
$`p_{3\omega }`$ $`=`$ $`{\displaystyle \frac{1}{4}}h_3\chi _1E^3{\displaystyle \frac{5}{16}}h_5\chi _1^2E^5,`$ (21)
$`p_{5\omega }`$ $`=`$ $`{\displaystyle \frac{1}{16}}h_5\chi _1^2E^5,`$ (22)
where

$`h_5`$ $`=`$ $`{\displaystyle \frac{81}{2}}{\displaystyle \frac{a^3|ϵ_2|^4}{|\rho |^4}}\{{\displaystyle \frac{9fϵ_2^2}{\rho ^2(ϵ_1+2ϵ_2)}}[{\displaystyle \frac{2(1f\beta )(ϵ_1ϵ_2)}{\rho }}+{\displaystyle \frac{3ϵ_2}{ϵ_1+2ϵ_2}}]+{\displaystyle \frac{3ϵ_2^2}{\rho (ϵ_1+2ϵ_2)^2}}`$ (24)
$`[{\displaystyle \frac{(1+3ff\beta )\rho (1f\beta )(\rho +3fϵ_13fϵ_2)}{\rho }}{\displaystyle \frac{2\rho +6f(ϵ_1ϵ_2)}{ϵ_1+2ϵ_2}}]\},`$
and

$`\chi _1^{1/2}E_\omega `$ $`=`$ $`j_1\chi _1^{1/2}E+{\displaystyle \frac{3}{4}}j_3(\chi _1^{1/2}E)^3+{\displaystyle \frac{5}{8}}j_5(\chi _1^{1/2}E)^5,`$ (25)
$`\chi _1^{1/2}E_{3\omega }`$ $`=`$ $`{\displaystyle \frac{1}{4}}j_3(\chi _1^{1/2}E)^3{\displaystyle \frac{5}{16}}j_5(\chi _1^{1/2}E)^5,`$ (26)
$`\chi _1^{1/2}E_{5\omega }`$ $`=`$ $`{\displaystyle \frac{1}{16}}j_5(\chi _1^{1/2}E)^5,`$ (27)
where

$$j_5=\frac{243}{4}\frac{(1f\beta )^2|ϵ_2|^5}{|\rho |^4}(\frac{3}{2}\frac{1}{\rho ^2|\rho |}+\frac{1}{|\rho |^3}+\frac{3}{2}\frac{\rho ^2}{|\rho |^5}).$$

Because of the weak nonlinearity under consideration, the effect of high order harmonics is quite small. Thus, without loss of generality, we only consider the Taylor expansion up to second term. In this regard, we will use Eqs.(18) and (19) to investigate the nonlinear ac response in the numerical calculations in next section.

## IV Numerical Results

We are now in a position to perform numerical calculations to investigate the dispersion effect on the harmonics of local electric field as well as induced dipole moment. We let $`f=0.15`$, $`ϵ_{1H}=6ϵ_0`$, $`ϵ_2=75ϵ_0`$, where $`ϵ_0`$ is the dielectric constant of free space. Additionally, set $`𝐩_0=ϵ_ea^3(ϵ_1ϵ_2)/(ϵ_1+2ϵ_2)𝐄`$ to normalize the induced dipole moment, where $`ϵ_e`$ is given by the anisotropic MGA (Eq.(5)) without nonlinear characteristic. Obviously, the expressions for the normalized local electric field and the normalized induced dipole moment are functions of a single variable $`\chi _1^{1/2}E`$ . In view of the weak nonlinearity, we set $`\chi _1^{1/2}E=0.9\sqrt{ϵ_0}`$ for our calculation.

In Fig.1(a)$``$(c), we investigate the harmonics of local field and induced dipole moment versus frequency for different intrinsic dispersion strengths $`\mathrm{\Delta }ϵ_1`$. It is obvious that increasing the dispersion strength reduces the harmonics of both the local field and the induced dipole moment. We find that, for increasing frequency, the harmonics of local field always increases. However, the harmonics of induced dipole moment exhibits a more complicated behavior. For small $`\mathrm{\Delta }ϵ_1`$, the ratio of the third to first harmonics of induced dipole moment decreases as the frequency increases, which is qualitatively in good agreement with the experimental data . However, for large $`\mathrm{\Delta }ϵ_1`$, its harmonic ratio also increases as the frequency increases. We expect this behavior can be demonstrated by experiment in the future.

In Fig.2, we investigate the harmonics of local field and induced dipole moment versus frequency for different intrinsic relaxation times $`\tau _0`$. A longitudinal field case ($`\beta _L=2`$) is discussed in Fig.2(a), while a transverse field case ($`\beta _T=0.5`$) in Fig.2(b). It is shown that, in either case, increasing relaxation time $`\tau _0`$ may enhance the harmonics of both local field and induced dipole moment. In addition, the transverse field case may predict stronger harmonics than the longitudinal field case within high-frequency region, which is actually valid for different dispersion strengths as well (no figures shown herein).

On the other hand, Klingenberg claimed that increasing external electric field leads to increasing harmonics of electric current . In this sense, we have also calculated the effect of external electric field on the nonlinear ac response of the induced dipole moment (no figures shown herein). Obviously, our calculations qualitatively demonstrate the experimental results  as well.

## V Discussion and conclusion

Here a few comments on our results are in order. In the present study, we have examined the case of nonlinear particles suspending in a linear host. Our considerations may be extended to a nonlinear host medium .

We have studied the effect of an intrinsic dispersion on the nonlinear ac response of a colloidal suspension. Actually, our formalism can be readily generalized to discuss two or more intrinsic dielectric dispersions of particles (or host fluid). In view of the fact that the intrinsic dielectric dispersion often occurs due to the surface conductivity or the inhomogeneous structure (e.g., coating effect ), we believe such effect also plays an important role in the dielectric behavior of biological cells, such as the electrorotation, dielectrophoresis and electro-orientation.

The perturbation method is used in the present work, but we believe similar results may also be predicted by a self-consistent theory .

In summary, a perturbation method has been employed to compute the local electric field and the induced dipole moment for suspensions in which the suspended particles have a nonlinear characteristic and an intrinsic dispersion, in an attempt to investigate the frequency effect on the nonlinear ac response. The results showed, for weak intrinsic dispersion strength, the ratio of the third to first harmonics of the induced dipole moment decreases as the frequency increases, which is qualitatively in agreement with experimental result. However, for a strong dispersion strength, the harmonics ratio increases as the frequency increases. Moreover, an increase in the intrinsic relaxation time may increase the strength of harmonics.

## Acknowledgments

This work was supported by the Research Grants Council of the Hong Kong SAR Government under project numbers CUHK 4284/00P and CUHK 4245/01P. K.W.Y. thanks Prof. A. R. Day for his interest in our work and his suggestion to examine the effects of inhomogeneities on the dielectric response of biological cell suspensions.

Fig.1 (a)

Fig.1 (b)

Fig.1 (c)

Fig.2 (a)

Fig.2 (b)