text:
# Deciding the Nature of the “Coarse Equation” through Microscopic Simulations: the Baby-Bathwater Scheme

## I Introduction

It is often the case that a microscopic or fine description of a physical system is available, while we are interested in its macroscopic or coarse behavior. Consider, as an example, a biased random walk model for which the particle density asymptotically evolves according to a macroscopic law such as the Burgers equation. Typically, the study of macroscopic behavior starts with obtaining a closed PDE-level description (a “coarse equation”) for the time evolution of the expected or ensemble-averaged fields of a few, low order moments of the micro-state phase space distribution. For our example, this would be the zero-th moment, the density field. Then, an array of mathematical and computational tools (numerical integration, fixed-point algorithms etc.) can be brought to bear on the coarse equation.

Over the last few years, we have been developing a class of numerical algorithms which attempt to analyze the coarse behavior without ever obtaining the coarse equation in closed form . The common character of these schemes is to use short, appropriately initialized bursts of microscopic simulations to estimate the quantities which, if the coarse equation was available, we would simply evaluate using the equation itself. Such quantities, estimated on-demand, include the time derivative of the evolving coarse fields, to be used in coarse projective integration , or the effect of the time-evolution operator for the implicit coarse Jacobian, to be used in Newton-Krylov type contraction mappings like the Recursive Projection Method (RPM)  or in eigenvalue/vector computations. These methods are based on matrix-free large scale scientific computing, and we sometimes collectively refer to them as “equation-free methods”. What makes these computations possible is the assumption of a separation of time scales in the dynamics of the evolving microscopic distribution. Typically, one finds that the hierarchy of coupled equations involving higher cumulants of the microscopic distribution constitutes a singularly perturbed problem: higher-order cumulants become, in the course of simulation, quickly slaved to (become deterministic functionals of) the lower-order cumulants. The consequences of slaving, realized in the computer as a black-box, embody the closures that allow us to solve for the coarse behavior. Fundamentally it is no different than if the closures are expressed in closed form first and then evaluated later. This way of thinking and newly developed computing technology can in practice exceed the traditional approach in both accuracy and total cost, especially if the constitutive relation is nonlinear and multi-dimensional. An example of this type is given in . An additional advantage of such methods is their ability to detect parametric regimes where the (number of moments used in the) present closure model is inadequate and hence appropriate refinements (including higher order moments) are necessary.

In the projective integration method  one takes advantage of the slow dynamics of the coarse variables to carry out only bursts of microscopic simulations connected via projections (in effect, extrapolations and/or interpolations) over gaps of time. In the same spirit of exploiting regularity, but now in space, we have developed the so-called gaptooth scheme  by evolving the full process only in an array of small spatial boxes (the teeth) separated by empty regions (the gaps). Clearly, the two methods are closely related by the physics of the problem. Indeed we can have a combined gaptooth-projective integration scheme; this is the focus of another paper . Here, we simply want to point out the fact that in the gaptooth method, the teeth communicate with each other via appropriate boundary conditions for the microscopic simulations performed inside them. And here lies the raison d’etre of this paper.

It is a well-known fact that certain features of a given equation affect the nature of the appropriate numerical solver. A Hamiltonian dynamics problem, for example, is best integrated by a symplectic integrator; often, finite difference solvers of partial differential equations (PDE) are built to respect certain properties of the PDE, such as conservation laws. Most importantly, the highest spatial order of an evolution equation critically affects the types of boundary conditions leading to a well-posed problem.

In a completely analogous manner, the way in which the microscopic model is solved seperately in each tooth in the gaptooth scheme, and the boundary conditions applied to the edges of each tooth, must respect the nature of the unavailable equations and their order. Furthermore, gaptooth algorithms compatible with conservation laws (e.g. using fluxes to estimate temporal derivatives, see for example ) are predicated upon knowing that the unavailable equation possesses certain conservation laws.

When the closed-form equation is available, some of these questions (e.g. the order of the highest spatial derivative in an evolution equation) can be answered by direct inspection. Other issues (such as the existence of conservation laws, or integrability) may, in the case of closed form equations, be relatively obvious, or may require a lot of work.

What we explore in this paper is the development of computer-assisted methodologies to answer the above questions when closed form equations are not available. The idea is that we can probe the consequences of these answers on the dynamics of the unavailable coarse equations using microscopic, particle- or agent-based, simulators by trying out large classes of appropriately chosen initial conditions. We will illustrate what we call the Baby-Bathwater algorithm on examples of particle systems realizing the Burgers and Korteweg-de Vries (KdV) equations, for the task of inferring the highest order of spatial derivative on the right-hand-side, and for answering questions concerning coarse conservation laws.

The paper is organized as follows: In Section II we briefly present our illustrative particle-based example. In Section III we discuss the determination of the highest order spatial derivative active in the “unavailable equation”. In Section IV we explore the possible existence of conservation laws. In the concluding Section we discuss the scope and limitations of the procedure, as well as additional questions that may be addressed through this approach. An interesting “twist” about reverse coarse integration arises in discussing the exploration of possible “coarse Hamiltonianity” of the unavailable coarse equation.

## II Numerical Experiment Setup

Our illustrative examples will be based on simple numerical experiments. We will first, as a sanity check, demonstrate the approach using a traditional numerical simulator of a known evolution equation as a “black box”. We will then substitute the simulator of the known continuum equation with a particle-based simulator, and repeat the procedure. Our first illustration will be the Burgers equation,

$$u_t+uu_x=\nu u_{xx},$$
(1)
as well as a particle based simulator constructed so that the evolution of its density resembles (at the appropriate limit, reproduces) the Burgers evolution. Since one of the issues to be explored is the number and type of boundary conditions in evolving the equation, our simulations must be possible without this a priori knowledge. We therefore use periodic boundary conditions (PBC) enforced on $`x[0,2\pi )`$ in all our exploratory simulations. One of the attractive features of the Burgers is that, for any initial profile $`u(x,t=0)`$, and even with PBC, the Cole-Hopf transformation  provides an analytical solution. The accuracy of the numerics can thus be checked directly.

A biased random walker-based particle simulator mimicking the Burgers dynamics was also constructed to demonstrate the direct application of our procedure on microscopic, particle based solvers. A detailed study of the features of this particle model is reported elsewhere . As a reference, the diffusion equation,

$$u_t=\nu u_{xx},$$
(2)
has the well-known microscopic realizations of Langevin dynamics or unbiased random walkers. It is not too difficult to conjure up a similar realization motivated by the Burgers equation using random walkers: a unit mass $`u(x,t)𝑑x=1`$ in the coarse description corresponds to $`Z`$ walkers, where $`Z`$ is a large integer constant. In the simulation, $`N`$ random walkers move on $`[0,2\pi )`$ at discrete timesteps $`t_n=nh`$. At each step, (a) the walkers’ positions $`\{x_i\}`$ are sorted, (b) each walker $`i`$ checks out the position of the walker $`m`$-places ahead, $`x_i^{m+}`$, and $`m`$-places behind, $`x_i^m`$ (properly accounting for PBC, of course). The difference $`x_i^{m+}x_i^m`$ is inversely proportional to the local density of walkers, therefore (c) every walker moves by $`\mathrm{\Delta }x_i`$ sampled from $`N(mh/Z(x_i^{m+}x_i^m),2\nu h)`$, a biased Gaussian distribution. The $`x_i`$’s are then wrapped around to $`[0,2\pi )`$, and the process repeats. This achieves a coarse-grained flux of $`ju^2/2\nu u_x`$ as motivated by the (1) by assigning each walker a drift speed of $`u/2`$. Quantifying the approximation of the Burgers evolution is an interesting subject that we take up separately in ; this is not, however, an important issue for this paper. It is only for benchmarking purposes that the relation to a known macroscopic equation is brought up. One can start by presenting a microscopic evolution law, without knowing anything about its corresponding coarse equation, and apply our algorithms on it directly.

Relating the fine with the coarse description requires the use of lifting and restriction operators. Lifting constructs particle distributions conditioned on some of their lower moment fields (here the zero-th moment, or coarse density field); it clearly is a one-many operator, and several microscopic “copies” of a given macroscopic initial condition are often required as discussed in detail in . The restriction (here computing moment fields of a given particle distribution) is a form of projection. Clearly the restriction (back to coarse variables) of a lifting of these coarse variables should be the identity (or close to it, due to noise effects). The lifting and restriction operators we constructed for this work, with $`u`$ interpreted as the coarse density on $`[0,2\pi )`$ PBC, are given in Appendix A. Fig. 1 shows a result of the “reversibility test”: we randomly generate a coarse density $`u(x)`$, lift it to a random walker distribution, then restrict back to $`\stackrel{~}{u}(x)`$, and observe the very good agreement between $`\stackrel{~}{u}(x)`$ with $`u(x)`$. Notice that the only point where this agreement may be less satisfactory is close to local maxima or minima, where the derivative changes sign.

While the proximity of our particle scheme to the Burgers evolution is not the issue in this paper, we briefly illustrate the correspondence of the evolution of an initial profile through the two approaches. Fig. 2 shows the analytical solution obtained through the Cole-Hopf transformation; it also contains the result of a $`21991`$-particle simulation, after the configurations have been processed by the $`\widehat{}`$ operator of Appendix A to extract the coarse density field estimate. A small value of viscosity $`\nu =0.1`$ is picked to accentuate the behavior of steepening wave-front with time. The microscopic simulation clearly captures the important features of the coarse behavior. Ensemble averaging with the same initial condition in coarse field $`u(x,0)`$ would reduce the error, but as we can see, even a single microscopic simulation using a reasonable number of particles may still perform quite efficiently.

Lastly, we mention the existence of particle methods to solve partial differential equations (PDEs) such as the Korteweg - de Vries (KdV) equation:

$$u_t=6uu_xu_{xxx}$$
(3)
which can be formulated as conservation laws . In this paper we use our construction above for the Burgers example. We should highlight once more that our ultimate purpose is not to construct particle solvers of given equations, but rather to decide on features of the unavailable coarse equations for given microscopic schemes.

## III Identifying The Highest Spatial Order of Coarse Variables

As we mentioned in Section I, system identification lies at the heart of the equation-free approach. Here, we suppose the coarse dynamics follows a certain time-evolution equation of the form:

$$u_t=f(u,u_x,u_{xx},..,u_x^{(N)}),$$
(4)
that is unavailable to us. We have already identified $`u`$ (the “coarse variable” for which we believe that a coarse deterministic equation exists in closed form), and constructed the lifting and restriction operators that connect macro/micro descriptions. We seek a general approach to decide qualitative questions, such as (a) what is $`N`$, and (b) whether $`f`$ can be written as $`𝐣`$, without having $`f`$ in closed form. One important motivation for this lies in that “production run” simulations of the problem via equation-free computation (for example through the gaptooth scheme) do not require knowledge of $`f`$, but are affected by the knowledge of $`N`$ (through “teeth” boundary conditions). What we do have is a microscopic simulator embodied in a computer code that can be initialized at will; the physical details of the microscopic code are both extremely important (that is where the “underlying physics” lies) and - for our purposes - irrelevant: we will use the microscopic simulation code as an “input-output” (I/O) black box. By probing the coarse I/O response of the black box, the question that we would like to address is whether we can decide on (a) and/or (b).

It may appear initially that we are trying to answer a circular question: in order to probe the coarse input-output response of a microscopic simulator we need to run it, and in order to run it we need well-posed boundary conditions, which - among other factors - depend on (a) and (b). To cut the knot, we use (for the decision stage exploratory runs) the Born-von Karman periodic boundary conditions. We are going to assume that the microscopic simulations can be carried out in PBC, which is an option prevalent among microscopic simulators. This enables us to probe the system’s response to only the initial $`u(x,0)`$ profile input.

The so-called baby-bathwater identification scheme works as follows:

1. Take an integer $`n`$, starting from $`1`$.
2. Pick a random point $`x_0`$ in the spatial periodic box.
3. Generate $`n`$ random numbers, designated as $`u(x_0,0)`$, $`u_x(x_0,0)`$, $`u_{xx}(x_0,0)`$, .., $`u_x^{(n1)}(x_0,0)`$ of $`u(x,0)`$.
4. Generate a conditionally random profile $`u(x,0)`$ compatible with the PBC and consistent with the above $`u(x_0,0)`$, $`u_x(x_0,0)`$, $`u_{xx}(x_0,0)`$, .., $`u_x^{(n1)}(x_0,0)`$ requirements. This can almost always be accomplished, for example, by summing $`2L`$ sine and cosine harmonics of the PBC:

$$u(x,0)=b_0+\underset{i=1}{\overset{L}{}}a_l\mathrm{sin}(lx)+b_l\mathrm{cos}(lx),x[0,2\pi ),$$
(5)
with $`L>n/2`$. Because we have $`2L+1`$ coefficients, even though there are $`n`$ constraints to satisfy, we still have some random degrees of freedom left in (5). In practice, this initialization can be accomplished by applying conjugate gradient minimization of the $`n`$-dimensional residual norm starting from a random $`\{a_l,b_l\}`$ vector.
5. Lift $`u(x,0)`$ of (5), run it in the microscopic simulator for time $`\mathrm{\Delta }`$, restrict it back to $`\stackrel{~}{u}(x,\mathrm{\Delta })`$, and estimate:

$$\stackrel{~}{u}_t(x_0,0)\frac{\stackrel{~}{u}(x_0,\mathrm{\Delta })\stackrel{~}{u}(x_0,0)}{\mathrm{\Delta }}.$$
(6)
Note that here $`\stackrel{~}{u}(x_0,0)`$ instead of $`u(x_0,0)`$ is used in the finite difference. This will cancel some internal noise from the lifting and restriction operations.
6. Repeat step (v) $`I`$ times to obtain an ensemble averaged $`\stackrel{~}{u}_t(x_0,0)`$ to reduce the microscopic noise.
7. Repeat step (iv) $`J`$ times, collect the $`\stackrel{~}{u}_t(x_0,0)`$ estimates:

$$(\stackrel{~}{u}_t^1(x_0,0),\stackrel{~}{u}_t^2(x_0,0),\mathrm{},\stackrel{~}{u}_t^J(x_0,0)),$$
(7)
compute the sample variance $`\sigma ^2(\stackrel{~}{u}_t(x_0,0))`$.
8. Repeat step (ii) $`K`$ times, compute the averaged sample variance $`\sigma ^2(\stackrel{~}{u}_t)_n`$.
9. Go back to step (i), $`nn+1`$. $`N`$ is identified when from $`n=N`$ to $`n=N+1`$, the averaged sample variance $`\sigma ^2(\stackrel{~}{u}_t)_{N+1}`$ decreases drastically to practically $`0`$.

Fig. 3 shows such families of constructed initial profiles with progressively more controlled initial derivatives. The basic idea is very simple: even though $`f`$ could have complicated functional dependencies on $`u,u_x,u_{xx},..,u_x^{(N)}`$, if they are all fixed, $`u_t`$ should have no dispersion even as $`u_x^{(N+1)},u_x^{(N+2)},\mathrm{}`$ are varied randomly. The “critical integer order” $`N`$ is identified when the variance at $`N`$ controlled derivatives $`u(x_0,0)`$, $`u_x(x_0,0)`$, $`u_{xx}(x_0,0)`$, .., $`u_x^{(N1)}(x_0,0)`$ jumps to a finite value; we then have already thrown out the “baby” (the highest relevant spatial derivative $`u_x^{(N)}`$) with the “bathwater” (the higher, non-relevant ones).

It is important to recognize that the time derivative estimation (6) does not occur instantaneously. A short “healing” period should elapse, during which the higher cumulants of the lifted phase space (micro-state) distribution become functionals of the lower order, slow governing cumulants. This separation of time scales, which fundamentally underlies the existence of a deterministic coarse equation closing with the lower cumulants, is discussed in more detail in .

As a sanity check, this algorithm is also applied to a traditional continuum PDE time-stepper “black box” first. Fig. 4(a) and 4(b) show the results of applying our decision scheme to forward Euler finite-difference PDE solvers of the Burgers and KdV equations, respectively. A spatial mesh of $`\mathrm{\Delta }x=2\pi /100`$ is adopted, and we define

$`u_x^{mk}`$ $`{\displaystyle \frac{u_{(m+1)\mathrm{\Delta }x,kh}u_{(m1)\mathrm{\Delta }x,kh}}{2\mathrm{\Delta }x}},`$ (8)
$`u_{xx}^{mk}`$ $`{\displaystyle \frac{u_{(m+1)\mathrm{\Delta }x,kh}+u_{(m1)\mathrm{\Delta }x,kh}2u_{m\mathrm{\Delta }x,kh}}{\mathrm{\Delta }x^2}},`$ (9)
$`u_{xxx}^{mk}`$ $`{\displaystyle \frac{u_x^{m+1,k}+u_x^{m1,k}2u_x^{mk}}{\mathrm{\Delta }x^2}},`$ (10)
$`u_{\mathrm{avg}}^{mk}`$ $`{\displaystyle \frac{u_{(m+1)\mathrm{\Delta }x,kh}+u_{m\mathrm{\Delta }x,kh}+u_{(m1)\mathrm{\Delta }x,kh}}{3}}.`$ (11)
We use,

$$\frac{u_{m\mathrm{\Delta }x,(k+1)h}u_{m\mathrm{\Delta }x,kh}}{h}=\nu u_{xx}^{mk}u_{m\mathrm{\Delta }x,kh}u_x^{mk},$$
(12)
to integrate the Burgers equation forward, and

$$\frac{u_{m\mathrm{\Delta }x,(k+1)h}u_{m\mathrm{\Delta }x,kh}}{h}=6u_{\mathrm{avg}}^{mk}u_x^{mk}u_{xxx}^{mk},$$
(13)
to integrate the KdV equation forward. $`u(x_0,kh)`$ is obtained by cubic spline over $`\{u_{m\mathrm{\Delta }x,kh}\}`$, and $`u_t(x_0,0)`$ is evaluated by finite differences, same as in (6). As can be seen in Fig. 4(a) and 4(b), $`N`$ is identified to be $`2`$ using the Burgers PDE time-stepper and $`3`$ using the KdV PDE time-stepper: the variances drop by more than four decades in both cases when going from $`N`$ to $`N+1`$ controlled derivatives. To see where the remaining “noise” comes from, note that,

$$u(x_0,\mathrm{\Delta })u(x_0,0)=u_t(x_0,0)\mathrm{\Delta }+u_{tt}(x_0,0)\frac{\mathrm{\Delta }^2}{2}+\mathrm{},$$
(14)
and clearly $`u_{tt}(x,0)`$ has higher-than-$`u_x^{(N)}`$ spatial derivative dependencies, which are, however, scaled by $`\mathrm{\Delta }`$ compared to the leading term. Thus the sample variances should drop by $`\mathrm{\Delta }^2`$ for $`n>N`$, which explains the observed magnitude of the four-decade decrease.

We then apply the identification scheme to the microscopic simulator of Section II, with the lifting and restriction operators constructed in Appendix A. The results are shown in Fig. 5. Under favorable conditions such as $`\nu =1`$ and $`m=100`$, it takes about $`10`$ minutes of computer time on a single 1GHz-CPU personal computer to obtain a reasonably good microscopic noise reduction so the variance drops by about $`2`$ decades going from $`n=2`$ to $`n=3`$. Under unfavorable conditions such as $`\nu =0.1`$ and $`m=1`$, it can take up to $`1,000`$ minutes of computer time to obtain the same $`2`$ decades drop. Compared with the deterministic finite-difference PDE time-steppers, identification of a microscopic simulator is undoubtedly much more computationally intensive, even though fundamentally there is no difference between the two “black boxes”. The problem of microscopic noise reduction is a persistent issue among all “equation-free” methods including bifurcation , projective / gaptooth integration , and identification, and calls for a unified treatment. This “one time” decision, performed at the beginning of studying a problem, will critically affect subsequent production runs of the microscopic simulator.

Here, one must pay special attention to the rank ($`M`$) of the restriction operator (see Appendix A). As can be seen in Figs. 1 and 2, our proposed restriction operator satisfies the constraint of reversibility and also accurately represents the profile’s long-time evolution. However, these merits do not guarantee automatically good short-time $`\stackrel{~}{u}_t`$ estimates by finite-difference. Special attention must be paid to the restriction operator $`\widehat{}`$: for example, if the highest harmonic in (5) for $`u(x,0)`$ is $`L`$, then with $`M=L`$, we can get good reversibility test of $`u(x,0)`$. Unless we use $`M=2L`$ for restricting the Burgers microscopic dynamics, however, we would not get a good estimate of $`u_t`$, because the nonlinear interaction $`uu_x`$ in (1) creates higher harmonics in $`u_t`$ up to $`2L`$. If $`M=L`$ is still used, it is equivalent to forcing a least-square projection of a $`4L+1`$ vector to a $`2L+1`$ subspace, which may work well enough in the long term, but is too inaccurate for short-term finite difference estimates. Unless this is taken care of, the $`u_t`$ estimate using our $`\widehat{}`$ is found to not even be superior to a crude bin-count density estimator with bin-width $`(2\pi /n)/8`$ about $`x_0`$, as $`2\pi /n`$ is the shortest wavelength in $`u(x,0)`$.

Lastly, we note that (4) represents a wide category of coarse dynamics; those with higher time-derivatives and mixed derivatives can be converted to a multi-variate version of (4) and the baby-bathwater identification scheme will still, in principle, work. A notable exception is the incompressible fluid dynamics case, where the sound-speed is infinite and the pressure plays the role of a global Lagrange multiplier. The incompressible fluid model is but a mathematical idealization of a certain physical limit. It is nonetheless useful and important enough, that the fact that it is not directly amenable to the baby-bathwater identification is worth mentioning. In general, the baby-bathwater identification presented here will not work for dynamics with instantaneous remote-action over macroscopic lengthscales, such as,

$$u_t(x,t)=𝑑\xi u(\xi ,t)K(x\xi ),$$
(15)
for which it is easy to show that $`u_t(x,t)`$ correlates with infinite number of local spatial derivatives $`\{u_x^{(n)}(x,t)\}`$.

## IV Identifying Conservation Laws

In section III above we address the concern of how to identify the highest spatial derivative of an unavailable coarse equation of the type (4). It is natural to try to decide other qualitative questions: for example, whether the coarse dynamics conserve a specific quantity,

$$Gg(u,u_x,u_{xx},..,u_x^{(N^{\prime \prime })})dx,$$
(16)
or not. In the simplest case, we ask whether $`gu`$ is conserved. We note that is equivalent to asking whether the RHS of (4) can be written as,

$$f(u,u_x,u_{xx},..,u_x^{(N)})=_xj(u,u_x,u_{xx},..,u_x^{(N^{})}),$$
(17)
or not. Alternatively, we ask whether there exists $`j(u,u_x,u_{xx},..,u_x^{(N^{})})`$ such that,

$$\frac{d}{dt}_{x_0}^{x_1}u(x,t)𝑑x=j(x_0,t)j(x_1,t),$$
(18)
for arbitrary $`x_0`$,$`x_1`$. Whereas in section III we try to identify features of $`f(u,u_x,u_{xx},..,u_x^{(N)})`$ through (4), here we can try to identify consequences of $`j(u,u_x,u_{xx},..,u_x^{(N^{})})`$ and its features through (18). The process of the baby-bathwater identification can be carried over; the only difference is that it is going to be a boundary scheme. In one dimension, the boundary sheme reduces to a two-point scheme as follows:

1. Take an integer $`n`$, starting from $`1`$.
2. Pick two random points $`x_0`$ and $`x_1`$ in the PBC.
3. Generate $`2n`$ random numbers, which are to be designated $`u(x_0,0)`$, $`u_x(x_0,0)`$, $`u_{xx}(x_0,0)`$, .., $`u_x^{(n1)}(x_0,0)`$ and $`u(x_1,0)`$, $`u_x(x_1,0)`$, $`u_{xx}(x_1,0)`$, .., $`u_x^{(n1)}(x_1,0)`$, of $`u(x,0)`$.
4. Generate a conditionally random profile $`u(x,0)`$ compatible with the PBC that is consistent with the above $`u(x_0,0)`$, $`u_x(x_0,0)`$, $`u_{xx}(x_0,0)`$, .., $`u_x^{(n1)}(x_0,0)`$ and $`u(x_1,0)`$, $`u_x(x_1,0)`$, $`u_{xx}(x_1,0)`$, .., $`u_x^{(n1)}(x_1,0)`$ requirements. This can always be done by (5) with $`L>n`$. As we discussed above, we have $`2L+1`$ coefficients, even though there are $`2n`$ constraints to satisfy, we still have some random degrees of freedom left in $`u(x,0)`$.
5. Lift $`u(x,0)`$ of (5), run it in the microscopic simulator for time $`\mathrm{\Delta }`$, restrict it back to $`\stackrel{~}{u}(x,\mathrm{\Delta })`$, estimate:

$$\stackrel{~}{U}_t(0)\frac{_{x_0}^{x_1}\stackrel{~}{u}(x^{},\mathrm{\Delta })𝑑x^{}_{x_0}^{x_1}\stackrel{~}{u}(x^{},0)𝑑x^{}}{\mathrm{\Delta }}.$$
(19)
6. Repeat step (v) $`I`$ times to obtain an ensemble averaged $`\stackrel{~}{U}_t(0)`$ to reduce the microscopic noise.
7. Repeat step (iv) $`J`$ times, collect the $`\stackrel{~}{U}_t(0)`$ estimates:

$$(\stackrel{~}{U}_t^1(0),\stackrel{~}{U}_t^2(0),\mathrm{},\stackrel{~}{U}_t^J(0)),$$
(20)
compute the sample variance $`\sigma ^2(\stackrel{~}{U}_t(0))`$.
8. Repeat step (ii) $`K`$ times, compute the averaged sample variance $`\sigma ^2(\stackrel{~}{U}_t)_n`$.
9. Go back to step (i), $`nn+1`$. A conservation law is positively identified when going from $`n=N^{}`$ to $`n=N^{}+1`$, the averaged sample variance $`\sigma ^2(\stackrel{~}{U}_t)_{N^{}+1}`$ decreases drastically to practically $`0`$.

Fig. 6 plots families of initial profiles thus constructed with progressively more controlled initial derivatives. Fig. 7(a) and 7(b) show the results of applying the identification scheme to the Burgers finite-difference PDE time-stepper (12) and the KdV finite-difference PDE time-stepper (13), respectively. $`N^{}`$ is identified to be $`1`$ for (12) and $`2`$ for (13).

Two comments are in order: first, we probe the consequences of conservation (i.e. that boundary fluxes are the only cause of change for the conserved quantity in a domain); second, we obtain (as a side-product) the highest spatial derivative of the conserved quantity $`u`$ in the constitutive equation for the flux. It is important to note that if the procedure progressively returns negative answers (e.g., if the sample variance is non-zero for a given number $`n`$ of controlled derivatives), this does not imply that a conservation law does not exist. It only implies that a conservation law of the class encompassed by our equation (4) with spatial derivatives up to the tested order $`n`$ does not exist. In that sense, our procedure provides sufficient confirmation, but its success is not necessary for a conservation law in a different class to prevail. Nevertheless, the class we consider is wide enough to encompass many known examples and problems of interest to applications.

Note that $`N^{}`$ is one order less than $`N`$ identified in section III in both cases. This is true in 1-D because of Eq. (17). So in 1-D, the baby-bathwater schemes give a definite answer to whether $`u`$ is conserved or not in a finite $`N1`$ steps, as long as $`N`$ is identified first.

## V Discussion

In section IV we proposed methods to check whether a coarse quantity (such as the mass corresponding to the coarse density $`g(u,u_x,u_{xx},..,u_x^{(N^{\prime \prime })})`$) is conserved, without knowledge of the coarse evolution equation. The obvious question that arises is how do we know which $`g`$ to check? The path that we suggest here, in the equation-free setting, is to examine the consequences of conservation laws. For example, consider the conservation of (linear) momentum. An equivalent statement, through Noether’s theorem  is the existence of translational invariance. If we numerically establish the latter, then we can claim the former. Let us then consider initial conditions to the available integrator which are shifts of an original profile e.g., $`u(xx_0)`$, $`u(x(x_0+ϵ))`$, $`u(xx_0+2ϵ)`$, etc. Then if we time evolve the problem, using our microscopic time-stepper, and the equation is translationally invariant, upon reaching the integration reporting horizon, we can back-shift the profile (by the original shift amount). If all back-shifts provide an identical profile, we can conclude translational invariance and hence linear momentum conservation. An additional note of caution is that the examination of such consequences is relevant when Noether’s theorem applies, hence when there is an underlying Lagrangian/Hamiltonian structure in the problem (we discuss separately the issue of Hamiltonian nature below). Notice, however the modulo the proviso of “Hamiltonianity”, this methodology can be used to establish additional dynamical invariants, e.g. the invariance with respect to phase of the evolution of a field can be related to norm invariance etc.

Establishing an underlying Hamiltonian structure in a sense proceeds in a similar fashion through its correlation with invariance with respect to time reversal. The crudest way to examine this is by simply running the integrator with a negative time-step (if that option is available). A more refined way to check the same symmetry is by examining computations of the spectrum (e.g., eigenvalues) of linearization of the coarse PDE. In particular, a straightforward consequence of the Hamiltonian nature is that all linearization eigenvalues should come in quartets, namely if $`\lambda `$ is an eigenvalue, then so are $`\lambda ,\lambda ^{},\lambda ^{}`$, where  denotes complex conjugation. It is fortunate that time-stepper based numerical analysis techniques for the numerical approximation of the leading spectrum of such a linearization are well-developed for the case of large scale continuum simulations (see for example ). If an eigenvalue $`\lambda `$ of the linearization is identified, matrix-free eigen-computations with shift can be used to explore the existence of the $`\lambda `$ eigenvalue (in general, real eigenvalues will come in pairs and complex conjugate eigenvalues in quartets).

While coarse time-reversibility can be answered by exploring the spectrum of the linearization, it raises the interesting question of how to integrate backward in time with the microscopic code. Consider a molecular dynamics configuration with a certain set of velocities at time zero, and the same molecular configuration with “flipped” velocities. A well known (and testable) consequence of microscopic reversibility and the “molecular chaos” ansatz is that, whether we integrate the molecular dynamics equations forward or backward in time starting from a randomly picked phase point, we will get “the same” forward in time evolution of the coarse macrosocpic observables. It is interesting, however, that the coarse projective integration techniques in an equation-free context can be used to attempt integration of the coarse variables backward in time (under appropriate conditions about the spectrum of the unavailable equation) as follows: Consider the lifting of a particular coarse initial condition to consistent molecular realizations; flipping the molecular velocities for these realizations will not affect the coarse procedure. We then evolve microscopically the molecular configurations (whether with the original or with flipped velocities) forward in time long enough for the higher moments to heal, say for a time $`\tau >\tau _{mol}`$. We now estimate the time derivative $`du/dt`$ of the healed coarse variables from the restriction of the “tail end” of the molecular trajectories, and then take a large, macroscopic Euler step backward in time for the coarse variables. We lift again, run molecular dynamics forward or backward microscopically, estimate the coarse forward time derivative, and take another coarse backward step. This procedure can of course be done in a much more sophisticated way as far as the coarse backward time step is concerned - algorithms like Runge-Kutta or Adams method can be combined with the MD computations to integrate density expectations backward in time for the coarse equations on the “slow manifold”. This “see-saw” forward-backward coarse integration procedure can also be used on stiff systems of ODEs and even dissipative PDEs under the appropriate conditions to evolve trajectories backwards on a slow manifold. The numerical analysis of these algorithms in the continuum case is an interesting subject in itself, and we are currently pursuing it . It is interesting that the technique, in the molecular dynamics case, can be used to coarsely integrate backward in time on a free energy surface, and thus help molecular simulations escape from free energy minima; we have already confirmed this in the case of Alanine dipeptide folding in water at room temperature through molecular dynamics simulations .

Finally, a more complicated question than checking the existence of one (a specific, and hence related to a specific invariance, in accordance with the above discussion) integral of the motion is the one of integrability. The latter necessitates infinite integrals of the motion, normally established by means of identifying Lax pairs and using the inverse scattering transformation machinery . However, one can also use in this case consequences of integrability to establish it. For instance, in recent work  it was qualitatively argued (and verified through numerical experiments in different settings) that a feature particular to integrable Hamiltonian systems is the presence of double continuous spectrum eigenvalues, when linearizing around a (coarse PDE) solitary wave under periodic boundary conditions. These as well as other criteria (such as the existence of point spectrum eigenvalues in the spectral gap ) can also be (conversely) used to potentially rule out the existence of integrable structure. In short, the spectral properties of the coarse PDE linearization can be used to establish or disprove not only the Hamiltonian (see above), but also potentially the integrable nature of the flow. While these are just initial thoughts towards attempting to decide vital questions about the nature of the unavailable closed equation, it is important to note that, what is computationally involved is a time-stepper based identification of facts about the spectrum of the linearization of an operator. This “computational technology” is quite mainstream in the case of large scale continuum simulators, and can be straightfowardly adapted to the case of coarse timesteppers in conjunction with the lifting-restriction steps. Variance reduction will clearly be the most significant step in the wide applicability of these and similar-spirited approaches.

Acknowledgements The authors would like to thank Frank Alexander for stimulating discussions. This work was partially supported by the National Science Foundation (IGK,PGK) AFOSR (CWG,IGK) and the Clay Institute (PGK).

## A Coarse Density Lifting / Restriction Operators

For a coarse field $`𝐮`$, the lifting operator $`\widehat{\mu }`$ generates a microstate $`𝐔`$: $`𝐔=\widehat{\mu }𝐮`$. Similarly, for a microstate $`𝐔`$, the restriction operator $`\widehat{}`$ returns a coarse field estimate $`\stackrel{~}{𝐮}`$: $`\stackrel{~}{𝐮}=\widehat{}𝐔`$. Both $`\widehat{\mu }`$ and $`\widehat{}`$ can be one-to-one or one-to-many operators, but we demand that $`\widehat{}\widehat{\mu }\widehat{I}`$ asymptotically when the wavelength of $`𝐮`$ is large enough compared to the microscopic length scale . For 1-D coarse density field $`u(x)`$ under $`x[0,2\pi )`$ PBC, we use the following $`\widehat{\mu }`$, $`\widehat{}`$ operators for the sake of definiteness in numerical experiments, even though their construction is not unique.

The lifting operator $`\widehat{\mu }`$, $`u(x)\{x_i\}`$:

1. Estimate $`u_{\mathrm{max}}^{\mathrm{approx}}u_{\mathrm{max}}\mathrm{max}_{x[0,2\pi )}u(x)`$. Pick $`u_{\mathrm{safe}}`$ that is “safely” greater than $`u_{\mathrm{max}}`$, for example $`u_{\mathrm{safe}}=1.1u_{\mathrm{max}}^{\mathrm{approx}}`$.
2. Define $`N^{}2\pi u_{\mathrm{safe}}Z`$. Create $`N^{}`$ particles $`\{x_i\}`$ with each $`x_i`$ independently drawn from uniform distribution on $`[0,2\pi )`$.
3. Go to each particle $`i`$, randomly decimate it with probability $`1\frac{2u(x_i)}{u_{\mathrm{max}}^{\mathrm{approx}}+u_{\mathrm{safe}}}`$. Count the total number of surviving particles $`N^{\prime \prime }`$.
4. Compute quadrature,

$$QZ_0^{2\pi }u(x)𝑑x$$
(A1)
randomly round to $`N=Q`$ or $`N=Q+1`$ such that $`N=Q`$. Randomly pick $`N^{\prime \prime }N`$ particles out of the $`N^{\prime \prime }`$ survivors and decimate them. We now have a set of particles $`\{x_i\}`$, totally numbered either $`Q`$ or $`Q+1`$.

The restriction operator $`\widehat{}`$, $`\{x_i\}\stackrel{~}{u}(x)`$:

1. Define microscopic density function,

$$a(x)\frac{1}{Z}\underset{i=1}{\overset{N}{}}\delta (xx_i^{})$$
(A2)
and corresponding cumulant function,

$$c(x)_0^xa(x^{})𝑑x^{}.$$
(A3)
Clearly, at the the first, second, third particle positions $`x_{n_1}`$, $`x_{n_2}`$, $`x_{n_3}`$, $`c(x_{n_1})=1/Z`$, $`c(x_{n_2})=2/Z`$, $`c(x_{n_3})=3/Z`$, etc. And we have $`c(0)=0`$, $`c(2\pi )=N/Z`$.
2. Define a residual function $`r(x)`$,

$$r(x)c(x)\frac{Nx}{2\pi Z}$$
(A4)
which is the difference between $`c(x)`$ and the cumulant of a homogenized particle gas background. The idea is that $`r(0)=r(2\pi )=0`$, so it is a periodic function and can be approximated by,

$$r(x)\stackrel{~}{r}(x)=\underset{i=1}{\overset{M}{}}a_n(\mathrm{cos}(nx)1)+b_n\mathrm{sin}(nx).$$
(A5)
In fact, a sound strategy is to least-square fit $`\stackrel{~}{r}(x)`$ (its $`\{a_n\}`$,$`\{b_n\}`$ coefficients) to $`r(x)`$ at $`x=x_{n_i}`$’s, where $`\{x_{n_i}\}`$ is the sorted list of $`\{x_i\}`$. $`r(x)`$ can be easily evaluated at $`x_{n_i}`$’s, noting the last sentence of step 1.
3. The coarse density estimate can be obtained by taking the derivative of $`\stackrel{~}{c}(x)Nx/2\pi Z+\stackrel{~}{r}(x)`$,

$$\stackrel{~}{u}(x)=\frac{N}{2\pi Z}+\underset{i=1}{\overset{M}{}}na_n\mathrm{sin}(nx)+nb_n\mathrm{cos}(nx).$$
(A6)

It is worth noting that although the constructed $`\widehat{}`$ depends on $`M`$, it satisfies the particle number conservation exactly because the finite harmonics all integrate to zero and only the background contribution remains. In fact, $`\widehat{}\widehat{\mu }`$ also satisfies exact particle number conservation to the original $`u(x)`$ under probabilistic average. Further, one can show $`\widehat{}\widehat{\mu }=\widehat{I}`$ exactly for $`u(x)`$ in the first $`M`$ harmonics subspace.


format: text

text:
# On the survey-propagation equations for the random K-satisfiability problem

## 1 Introduction

Recently many progresses  have been done on the analytic and numerical study of the random K-satisfiability problem , using the approach of survey-propagation that generalizes the more old approach based on the “Min-Sum” 111The “Min-Sum” is the the zero temperature limit of the “Sum-Product” algorithm and sometimes is also called belief propagation. In the statistical mechanics language  the belief propagation equations are the extension of the TAP equations for spin glasses  and the survey-propagation equations are the TAP equations generalized to the broken replica case. algorithm  .

The aim of this note is to sketch a possible path to a proof of the existence an uniqueness of the survey-propagation equations for a given instance of the problem. Before presenting the main arguments, for reader convenience I will present an heuristic derivation of the survey-propagation equations in section 2; the full derivation can be found in the original papers . In section 3, I will present the main sequence of conjectures that may lead to the proof of the existence and uniqueness of the survey-propagation equations in the appropriate range of parameters. Finally I will present some conclusions.

## 2 A fast heuristic derivation of the survey equations

### 2.1 The random K-sat problem

In the random K-sat problem there are $`N`$ variable $`\sigma (i)`$ that may be true of false (the index $`i`$ will sometime called a node). An instance of the problem is given by a set of $`M\alpha N`$ clauses. Each clause is characterized by set of three nodes ($`i_1`$,$`i_2`$, $`i_3`$), that belong to the interval $`1N`$ and by three Boolean variables ($`b_1`$,$`b_2`$, $`b_3`$). In the random case the $`i`$ and $`b`$ variables are random with flat probability distribution. Each clause $`c`$ is true if the expression

$$E_c(\sigma (i_1^c)XORb_1^c)OR(\sigma (i_2^c)XORb_2^c)OR(\sigma (i_3^c)XORb_3^c)$$
(1)
is true. The problem is satisfiable iff we can find a set of the variables $`\sigma `$ such that all the clauses are true. The entropy  of a satisfiable problem is the logarithm of the number of the different sets of the $`\sigma `$ variables that make all the clauses true.

To a given problem we can associate a graph (the factor graph ) where the nodes are connected to the clauses ($`3\alpha `$ in average) and each clause is connected to three nodes. The properties of this graph play a very important role.

The goal of the analytic approach consists in finding for given $`\alpha `$ and for large values of $`N`$ the probability that a random problem (i.e. a problem with random chosen clauses) is satisfiable. The $`01`$ law  is supposed to be valid: for $`\alpha <\alpha _c`$ all systems (with probability one when $`N`$ goes to infinity) are satisfiable and their entropy is proportional to $`N`$ with a constant of proportionality that does not depend on the problem. On the other hand, for $`\alpha >\alpha _c`$ no random system (with probability one) is satisfiable. An heuristic argument has been given  that suggest that $`\alpha _c=\alpha ^{}4.27`$ where $`\alpha ^{}`$ can be computed using the survey-propagation equations defined later. There is already an incomplete proof  that $`\alpha ^{}`$ is a rigorous upper bound to $`\alpha _c`$.

### 2.2 The belief propagation equations

In the following analysis it will be crucial that in limit $`N\mathrm{}`$ the problem becomes locally trivial: this makes possible the computations of $`\alpha _c`$ and of the other properties of the system. Let us be more precise. We can define a distance among the $`N`$ nodes in the following way:

* Two nodes are at distance 1 if there is a clause that contains both of them. At large $`N`$ the number of nodes at distance 1 from a given node has a Poisson distribution with average $`6\alpha `$.
* Two nodes are at distance $`k`$ if they are not at a distance $`k1`$ and there is a chain of $`k`$ overlapping clauses that touch both of them (or equivalently the second node is at distance 1 from a node at distance $`k1`$ from the first node).

In the limit $`N\mathrm{}`$ the set of nodes at distance $`k`$ from a given node form a tree, without closed loops: locally the system looks like a tree. The solution of the K-sat problem on a tree (with given boundary conditions) can be trivially done: the belief propagation algorithm (defined later) is exact. When $`N`$ goes to infinity the systems is not a tree (with probability 1) and this fact makes hard to find an actual solution of the problem for a given instance; it may even destroy the global existence of a solution to the belief propagation equations. However the local treeness of the problem is enough to allow an analytic treatment.

Let us take a large system for $`\alpha <\alpha _c`$ and let us consider the set $`C`$ of all configurations that satisfies all the clauses. Our first task is to compute the probabilities $`p_T(i)`$ and $`p_F(i)`$ that the variable $`\sigma (i)`$ is true or false (obviously $`p_T(i)+p_F(i)=1`$), if $`\sigma (i)`$ belong to a random configuration in $`C`$.

The presence of a simple local structure allows us to write down simple local equations . In the case of belief propagation equation  one proceed as follows. For each clause that contains the node $`i`$ (we will use the notation $`ci`$ although it may be not the most appropriate) $`p_T(i,c)`$ is the probability that the variable $`\sigma (i)`$ would be true in absence of the clause $`c`$. If the node $`i_1^c`$ were contained in only one clause, we would have that

$`p_T(i_1^c)=u_T(p_T(i_2^c,c),p_T(i_3^c,c),b_1^c,b_2^c,b_3^c)u_T(i,c),`$
$`p_F(i_1^c)=1u_T(i,c),`$ (2)
where $`u_T`$ is an appropriate function that is defined by the previous relation. An easy computation shows that when all the $`b`$ are false, the variable $`\sigma (i_1^c)`$ must be true if both variable $`\sigma (i_2^c)`$ and $`\sigma (i_3^c)`$ are false, otherwise it can have any value. Therefore we have in this case that

$$u_T(i,c)=p_F(i_2^c,c)p_F(i_3^c,c)+\frac{1p_F(i_2^c,c)p_F(i_3^c,c)}{2}$$
(3)
In a similar way, if some of the $`b`$ variable are true, we should exchange the $`p_F`$ with the $`p_T`$ for the corresponding variables. Finally we have that

$`p_T(i,c)={\displaystyle \frac{\underset{di,dc}{}u_T(i,d)}{Z(i,c)}}`$
$`p_F(i,c)={\displaystyle \frac{\underset{di,dc}{}u_F(i,d)}{Z(i,c)}}`$
$`Z(i,c)={\displaystyle \underset{di,dc}{}}u_T(i,d)+{\displaystyle \underset{di,dc}{}}u_F(i,d).`$ (4)
In total there are $`3M`$ variables $`p_T(i,c)`$ and $`3M`$ equations. These equations are called in the literature under different name (e.g. belief propagation, TAP equations ) and we naively expect that these equations (belief propagations) are satisfied (or quasi-satisfied, i. e. corrections of order $`1/N`$ can be present ). For the time being let us suppose that such a solution exists and it is unique. In such case we expect that

$`p_T(i)={\displaystyle \frac{\underset{di}{}u_T(i,d)}{Z(i)}}`$
$`p_F(i)={\displaystyle \frac{\underset{di,}{}u_F(i)}{Z(i)}}`$
$`Z(i)={\displaystyle \underset{di,}{}}u_T(i,d)+{\displaystyle \underset{di}{}}u_F(i,d).`$ (5)
We note the previous formulae can be written in a more compact way if we introduce a two dimensional vector $`\stackrel{}{p}`$, with components $`p_T`$ and $`p_F`$. We define the product of these vector

$$c_T=a_Tb_Tc_F=a_Fb_F,$$
(6)
if $`\stackrel{}{c}=\stackrel{}{a}\stackrel{}{b}`$.

If the norm of a vector is defined by

$$|\stackrel{}{a}|=a_T+a_F,$$
(7)
we finally find that

$$p(i,c)=\frac{\underset{di,dc}{}u_T(i,d)}{|_{di,dc}u_T(i,d)|}.$$
(8)
For each clause $`c`$ we can define the probability that the clause would be satisfied in a system where the clause is not present. We will denote this probability by $`Z(c)`$. One finds that in the case where all the $`b`$ variables are false

$$Z(c)=1p_F(i_1^c)p_F(i_2^c)p_F(i_2^c).$$
(9)
Finally the total entropy (apart correction that are subleading when $`N`$ goes to infinity) is given by

$$S=\underset{i}{}\mathrm{log}(Z(i))+2\underset{c}{}\mathrm{log}(Z(c)).$$
(10)
### 2.3 The survey propagation equations

One can argue that the situation is more complex . The belief propagation algorithm works at low value of $`\alpha `$ but it fails when $`\alpha `$ becomes too large. The previous equations may have multiple solutions or quasi-solutions that are very different one from the other.

This fact has been interpreted  in the following way: the set of all configurations $`C`$ that satisfy all the formulae can be divided into many sets that are well separated one by the others (this sets are sometimes called states in statistical mechanics  or lumps ). The previous belief equations correspond to the probability restricted to one given set.

The new picture is the following: we have an exponential large number of solutions (or quasi-solutions) of the belief equations (we call this set $`B`$) and we would like to know this number (i.e. the exponential of the complexity $`\mathrm{\Sigma }(\alpha )`$). The total number of configurations that satisfies all the clauses is given by

$$\mathrm{exp}(S(\alpha ))=\mathrm{exp}(\mathrm{\Sigma }(\alpha )+S_s(\alpha ))$$
(11)
where $`\mathrm{exp}(S_s(\alpha ))`$ is the total number of configurations that satisfy all the clauses in a generic state.

One expects that $`\mathrm{\Sigma }(\alpha )`$ vanishes at $`\alpha _c`$, so that it computation is extremely important. In order to reach this goal we can mimic the steps that we have done for going from the variables $`\sigma `$ to the probability $`p_T(i)`$; however in this case we are going to play the game at an higher level of abstractions, where states (or quasi-solutions of the belief equations) play the same role of a single configuration of the $`\sigma (i)`$ in the previous approach.

The new quantity is the full survey probability $`s(i|\stackrel{}{p})`$, i.e. a function of $`\stackrel{}{p}`$ that is defined at each node that is equal to the probability of finding a solution (or a quasi-solution) of the belief equations with $`\stackrel{}{p}(i)=\stackrel{}{p}`$. In other words $`s(i|\stackrel{}{p})`$ an indirect probability, i.e. a probability of a probability.

We introduce the quantity $`s(i,c|\stackrel{}{p})`$ that is the distribution probability for the probability (i.e. $`s(i|\stackrel{}{p})`$) when the clause $`c`$ is removed. The equations for the full survey probability $`s(i,c|\stackrel{}{p})`$ can be obtained using the techniques of , however we will not consider them here. Indeed we need them if our aim is to compute the total entropy (or $`S_s(\alpha )`$). If our aim is more modest and we want to compute only $`\mathrm{\Sigma }(\alpha )`$ and consequently $`\alpha _c`$, a simple approach is possible .

The crucial observation is that a given solution of the belief equations may have $`p_T(i,c)=1`$ or $`p_T(i,c)=0`$ or $`0<p_T(i,c)<1`$.

We can coarse grain the probability distribution of the beliefs by introducing the quantity $`s_T(i,c)`$ that is defined as the probability of finding $`p_T(i,c)=1`$, in the same way $`s_F(i,c)`$ as the probability of finding $`p_T(i,c)=0`$ and $`s_I(i,c)`$ is the probability of finding $`0<p_T(i,c)<1`$. As discussed in , by considering the equations for only these coarse grained surveys, it is possible to compute the complexity $`\mathrm{\Sigma }`$ and the value of $`\alpha _c`$.

We can use a more compact notation by introducing a three dimensional vector $`\stackrel{}{s}`$ given by

$$\stackrel{}{s}=\{s_T,s_I,s_F\}.$$
(12)
Everything works as before with the only difference that we have a three component vector instead of a two component vector. Generalizing the previous arguments one can introduce the quantities $`\stackrel{}{u}(i,c)`$ that is the value that the survey at $`i`$ would take if only the clause $`c`$ would be present in $`i`$. In the case where all the $`b`$ are false, a simple computation gives

$$\stackrel{}{u}(i,c)=\{s_F(i_2^c,c),s_F(i_3^c,c),1s_F(i_2^c,c),s_F(i_3^c,c),0\}.$$
(13)
The formula can generalized as before 222It always happens that the vector $`\stackrel{}{u}`$ has only one zero component ($`u_Tu_F=0`$). This fact may be used to further simplify the analysis. to the case of different values of $`b`$. One finally finds

$$\stackrel{}{p}(i,c)=\frac{\underset{di,dc}{}\stackrel{}{u}(i,d)}{|_{di,dc}\stackrel{}{u}(i,d)|},$$
(14)
where we have defined product in such a way that

$$\stackrel{}{a}\stackrel{}{b}=\{a_Tb_T+a_Ib_T+a_Tb_I,a_Ib_I,a_Fb_F+a_Ib_F+a_Fb_I\}.$$
(15)
The previous equations are the survey propagation equations (as defined in ) and the reader can find there the details of the derivation.

If one finds a solution to the survey equation one can compute the survey probability as

$$\stackrel{}{s}(i)=\frac{\underset{di,}{}\stackrel{}{u}(i,d)}{|_{di}\stackrel{}{u}(i,d)|}$$
(16)
and the total complexity is given by

$$\mathrm{\Sigma }=\underset{i}{}\mathrm{log}(Z(i)+2\underset{c}{}\mathrm{log}(Z(c))$$
(17)
where now

$$Z(i)=\mathrm{ln}(|\underset{di}{}\stackrel{}{u}(i,d)|),Z(c)=\mathrm{ln}(|\stackrel{}{s}(i,c)\stackrel{}{u}(i,c)|)$$
(18)
## 3 The infinite rooted tree and some conjectures

In the previous section we have sketched an heuristic derivation of the survey equations: it should be clear that these equations should be taken as they are and there is no warranty of any kind, either expressed or implied, that the belief equations (that have been used heuristically to construct the survey equations) do have a global solution.

The aim of this note is not to present a more precise derivation of the survey equation, but to discuss the problem of finding solutions to the survey propagation equations: this is a well posed mathematical problem independently from the origin of these equations. If the survey equations would have no solutions, the previous arguments would be empty or if survey equations had an exponential number of solutions, we should go up to an higher level of abstraction.

Numerical experiments on systems of size up to $`N=10^6`$ show that in the interval $`\alpha _L<\alpha <\alpha _U`$ ($`\alpha _L`$ and $`\alpha _U`$ asymptotically do not depend on the system, they are near to 3.9 and 4.36 for large $`N`$), the survey equations do have one solution that is obtained by iterations. For $`\alpha <\alpha _L`$ the survey equation converge to the trivial solution $`s_I(i,c)=1`$. On the other end for $`\alpha >\alpha _U`$ the iterative equations do not converge. Moreover the difference among a solution for the survey propagation equations and those for a perturbed survey propagation equation (e.g. by adding or removing a clause, or by fixing a survey to an arbitrary value) diverges when we approach $`\alpha _U`$ from below. The complexity change sign at a value $`\alpha ^{}`$ such that $`\alpha _U>\alpha ^{}>\alpha _L`$.

### 3.1 The infinite tree

These results call for an analytic proof. The aim of this note is to suggest a possible approach.

We propose to generalize the construction that Aldous has recently used in the study of the random matching problem . For any node $`i`$ of a given problem for finite $`N`$ we associate an infinite tree rooted in $`i`$ that is constructed in the following way 333Sometimes the tree is not infinite, e.g. if the site $`i`$ has no neighbour, however for not too small $`\alpha `$ the tree is infinite in most of the cases. If no loops were present (a rather unlikely possibility for large $`N`$) in the original graph the infinite rooted tree would be identical to the original graph.. Let us denote by $`\gamma `$ a node of the infinite tree. We assume that there is a function $`M(\gamma )`$ that maps the nodes of the infinite tree onto the nodes of the original problems in such a way that if $`\gamma _1`$, $`\gamma _2`$ and $`\gamma _3`$ belong to the same clause, also the three nodes $`i_k=M(\gamma _k)`$ belong to the clause, (the variables $`b_k`$ of the two clauses have the same values). We can further impose that the number of nodes at distance one from $`\gamma `$ is equal to the number of nodes at distance one from $`M(\gamma )`$. The construction is simpler that it may looks. In the cases of a problem with four variables and clauses that involve only two sites (i.e. a 2-sat, not 3-sat for graphical convenience) the original graph is shown in figure 1, while the center of the infinite graph, rooted in 0, is shown in fig. (2).

The construction of the infinite rooted tree is the simplest way to project a finite graph on a tree preserving as much as possible the structure of the original graph. The fact that for large $`N`$ the original problem does not contain small loops implies that the original problem is locally very similar to a subset of the infinite tree. This fact suggest that the solution on the K-sat problem on the infinite tree may give us information on the solution of the K-sat problem on the original problem.

An other object that we can construct is a random infinite tree. It is a random tree where the number of nearest nodes has a Poisson distribution. It is evident that for large $`N`$ the infinite tree associated to a given random problem becomes locally very near to a random infinite tree: the first contains correlations that vanish when $`N`$ goes to infinity. Moreover the properties the infinite random tree can be studied analytically. One would like to compare the properties of a given problem with the properties of the associated infinite tree; the hope is that the infinite tree associated to a random problem with $`N`$ clauses should become similar to the random infinite tree when $`N`$ goes to infinity.

We will say that a problem on the infinite tree has a unique solution iff, when we impose a generic boundary conditions one the $`k+1`$ shell, the behaviour at the center of the tree does not depend on the boundary conditions with probability one when $`k`$ goes to infinity.

The problem of computing a solution that satisfies all the clauses cannot have a unique solution in the previous sense: the entropy has a term proportional to $`N`$ and there is an exponentially large number of different solutions.

The belief propagation equations have an unique solution for a given belief on the boundary (quasi solutions fade away if closed loops are absent). An explicit computation show that when $`\alpha `$ is greater than a critical value (around 3.9)  the solution of the belief does depend on the boundary. One could also argue that if the solution of the belief equations would be independent from the boundary conditions when $`k`$ goes to infinity, there should be an essential unique solution of the belief equation for a given problem and for sufficiently large $`\alpha `$ this does not happen  (i.e. Aldous essential uniqueness property fails for the belief equations).

### 3.2 Three conjectures

We have now to exclude the possibility that the survey propagation equations do have a stable solution on the rooted tree near $`\alpha _c`$. This is an highly non-trivial requirement that fails in other models or for other forms of survey equations in the K-sat problem.

If the approach of  is correct the following conjectures should be true:

1. For $`\alpha <\alpha _U`$ and large $`N`$ the infinite tree associated to random problem has one stable solution with probability one. If this happens, the corresponding survey probabilities will be denoted by $`\stackrel{}{s}_{\mathrm{}}(i,c)`$.
2. For large $`N`$ the survey probabilities $`\stackrel{}{s}_{\mathrm{}}(i,c)`$, should be an approximate solution of the survey probability equations of the original problem. In other words for a given sample and $`\alpha <\alpha _U`$ there is solution of the survey probability equations near to $`\stackrel{}{s}_{\mathrm{}}(i,c)`$.
3. The value of $`\alpha _U`$ is determined by the following condition: for $`\alpha <\alpha _U`$ the infinite random tree has only one stable solution, while this does not happens $`\alpha >\alpha _U`$.

A direct numerical test of these conjectures is not easy especially in the interesting region where $`\alpha `$ is near to 4. In principle it is possible to construct in an explicit way the first $`k`$ shells of the infinite tree and to study what happens for large $`k`$. Unfortunately the number of first neighbour nodes is $`6\alpha `$ so that for $`\alpha =4`$ the $`k^{th}`$ shell contains of the order $`24^k`$ nodes, that is a very large number for numerical analysis already for $`k=5`$. If $`N`$ is not much greater than $`24^5`$, there will be many repetitions of the same nodes in the first five shells of the tree and on such scale the original problem does not looks very tree like. I have done studied numerically problems up to $`N=10^6`$ and $`k=4`$ and the data (e.g. at $`\alpha =4`$) are consistent with the first two conjectures although it is difficult to arrive to convincing evidence.

If we accept also the last conjecture, we can compute the value of $`\alpha _U`$ on the infinite random lattice and we can compare it with the numerical estimates for a given problem, i.e. $`\alpha _U4.36`$. At this end we must study the survey equations on the infinite random tree. In this case it is natural to suppose a translational invariance property .

Let us call $`𝒫_k(\stackrel{}{s})`$ the probability distribution of the survey in a generic node at distance $`k`$ form the origine 444For each node $`i,`$ we have to consider all the different surveys with one of the clauses removed, for simplicity we will not indicate the $`c`$ dependence and in the following $`\stackrel{}{s}(i)`$ will be nickname for all the $`\stackrel{}{s}(i,c`$).. Translational invariance implies that $`𝒫_k(\stackrel{}{s})`$ does not depend on $`k`$: it will be denoted by $`𝒫(\stackrel{}{s})`$. This quantity plays a crucial role in the approach 555It is convenient to recall the heuristic definition.of $`𝒫(\stackrel{}{s})`$. We decompose into states the set of the configuration that satisfies all the conditions. For each state we compute the belief probability that a given variable is true. The survey probability characterizes the distribution probability of the belief at a given site: the survey is a probability of a probability (an indirect probability). Finally $`𝒫(\stackrel{}{s})`$ is the probability of finding a site with that particular survey probability. In other words $`𝒫(\stackrel{}{s})`$ is a probability of a probability of a probability; in the simplifying case we are studying here is a function (we only care if a belief is equal to $`\pm 1`$ or not), while in the more general case it would be a functional..

### 3.3 A consistency check

In principle it possible that this construction fails: the equations for the survey may have a solution that depends on the boundary condition. In such a case a more complex construction should be done  and the aim of this note is to exclude that this happens.

The properties of $`𝒫(\stackrel{}{s})`$ can be well studied numerically and hopefully analytically. Indeed the surveys of the nodes on a shell at distance $`k`$ from the origine can be expressed in terms of the surveys at the nodes at distance $`k+1`$; using the supposed translational invariance of the probability distribution we get a consistency equation.

The procedure is the following. We consider a node with $`z`$ clauses, where $`z`$ has a Poisson distribution with average $`3\alpha `$ and we assume that the nearest nodes have the survey probability distributed according to $`𝒫(\stackrel{}{s})`$ and we compute the survey probabilities for the new node. If we average on $`z`$ and on the random clauses we get a new survey probability $`\stackrel{~}{𝒫}(\stackrel{}{s})`$, that obviously depends on $`𝒫(\stackrel{}{s})`$. The equation for $`𝒫(\stackrel{}{s})`$ are simply

$$\stackrel{~}{𝒫}(\stackrel{}{s})=𝒫(\stackrel{}{s}).$$
(19)
Using the techniques of  one can also construct a functional $`\mathrm{\Phi }(𝒫)`$ such that the actual solution of the equation (19) can be found by minimizing this functional (however we will not discuss this point).

It is a standard conjecture  that equations of the form (19) can be solved by iteration (and this likely follows from the convexity properties of $`\mathrm{\Phi }(𝒫)`$). If this is the case, we can use the method of population dynamics to find the solution of equation (19). The method is very simple and can be trivially implemented on a computer. It consists in describing a probability $`P`$ by an ensample of $`L`$ elements distributed according to this probability; the method becomes exact when $`L`$ goes to infinity.

We consider a set of $`L`$ surveys. Starting from this set we generate a new ensemble of $`L`$ surveys by using the standard procedure described in . The construction of an element of the new ensample is done as follow. We extract a Poisson distributed integer $`z`$ with average $`3\alpha `$. We take $`6z`$ surveys extracted randomly from the $`L`$ surveys and we extract random the $`b`$ of the corresponding $`k`$ clauses. Using equations (13,15) we compute one of the surveys that belong the new ensemble. Finally by repeating this operation $`L`$ times we obtain the new ensample.

By iterating this procedure $`I`$ time we find for large $`I`$ a probability distribution on the surveys, that is $`I`$ dependent. This procedure can be done also for large values of $`L`$ (e.g. $`L=10^6`$) and the convergence is rather fast (the corrections seem to be proportional to $`1/L`$ for generic $`\alpha `$). If the limits $`L`$ to $`\mathrm{}`$ and $`I`$ are smooth (the first should be done firstly) the resulting probability distribution is a solution of the equation (19).

In the same approach we can ask what happens in the population dynamics if we start from two different sets of surveys at the initial step. Let us indicate the $`i^{th}`$ survey at the iteration $`t`$ with $`\stackrel{}{s}(i,t)`$. Let assume to run twice the population algorithm (with the same random number generator) but taking two different sets as starting points: $`\stackrel{}{s}_1(i,0)`$ and $`\stackrel{}{s}_2(i,0)`$. We expect that for large $`i`$ the probability distribution of the survey should be the same; however it is a not evident if for a given $`i`$

$$\stackrel{}{s}_1(i,t)\stackrel{}{s}_2(i,t)_t\mathrm{}0.$$
(20)
It is natural to conjecture that if this happens, the survey equations have an unique solution on a random infinite tree. Indeed the computation of the surveys on the shell $`MI`$ as function of the survey on the shell $`MI+1`$, can be done exactly using the same algorithm we use in the population dynamics, with the only difference that the total number of surveys is constant in the population dynamics (i.e. it is equal to $`L`$) and it decreases with $`I`$ on a three (the number is proportional to $`(6\alpha )^{MI})`$). In the limits $`L`$ and $`M`$ going to infinity this difference should not be relevant.

In order verify if equation (20) is true is convenient to define a distance $`D(t)`$ as

$$D(t)=\frac{\underset{i=1,L}{}|\stackrel{}{s}_1(i,t)\stackrel{}{s}_2(i,t)|^2}{L}.$$
(21)
I have numerical studied the properties of $`D(t)`$ for large $`L`$ (up to $`L=10^6`$) finding little $`L`$ dependence (as expected). I find that for $`\alpha <\alpha _U4.36`$ for large $`t`$:

$$D(t)\mathrm{exp}(\lambda (\alpha )t),$$
(22)
where $`\lambda (\alpha )`$ is positive for $`\alpha <\alpha _U`$ and vanishes at $`\alpha =\alpha _U`$. For $`\alpha >\alpha _U`$, but not too large both $`vecs_1`$ \- and $`\stackrel{}{s}_2`$ go to the same probability distribution, but $`D(t)`$ does not go to zero. In other words $`\lambda (\alpha )`$ is maximum Liapunov exponent: when it is negative the iteration converges to a fixed point, while when it is positive chaos is present.

The estimated value of $`\alpha _U`$, i.e. 4.36, is larger that the value where the complexity vanishes (i.e. 4.27), so that the previous conjectures should be correct in the interesting region of positive complexity, where there should be solutions of the satisfiability conditions that correspond to the solutions of the survey propagation equations.

The reader should notice that the survey propagation equations do have a solution also for $`\alpha >\alpha _c`$ and the fact that the complexity turns out be negative is a warning that original satisfiability problem does not have any solution. In this case the heuristic construction of the surveys is empty.

## 4 Conclusions.

The main propose of this note is the construction (following Aldous ) of the infinite rooted tree associated to given satisfiability problem. This infinite rooted tree plays the role of a bridge among a finite instance of the problems and the infinite random tree where analytic computations  are done. It is argued that existence of an unique stable solution on the infinite tree (that apparently holds for $`\alpha <\alpha _U4.36`$) implies the existence of an unique stable solution of the survey equations on a large system in the same range of $`\alpha `$.

This result implies that the survey equations that have been used  in an algorithm to find an actual solution of an instance of the K-sat problem do have a solution. However independently from the interest of this application of the survey equations, I believe that the conjectures that have been put forward have a mathematical interest in their own because they clarify the fundamental hypothesis behind the approach of ; eventually they can be used to prove similar results in other problems (like the $`p`$-spin model).

## Acknowledgements

I thank Marc Mézard and Riccardo Zecchina for useful discussions and encouragement.


text:
# Multivector Functionalspublished: Advances in Applied Clifford Algebras 11(S3), 93-103 (2001).

## 1 Introduction

This is the last paper (VII) of series of papers dealing with the theory of multivector and extensor functions and multivector functionals. It is dedicated to the introduction of a key concept, that of *multivector functionals* and the study of their properties. Particularly important is the concept of *induced multivector functionals.* Several kinds of derivatives of multivector functionals, such as $`A`$-directional derivative and generalized concepts of curl, divergence and gradient are defined. Since the subject of the present paper has not been explored in the literature111For the best of our knowldge the only place where the concept has been rudimentary used was in . The concept has been used also in (,)., we present in section 3 several examples worked in detail of calculations of different types of derivatives for multivector functionals. Multivector functionals are fundamental for the formulation of the Lagrangian field theory of multivector and extensor fields on an arbitrary manifold, a subject that will be studied in a new series of papers.

## 2 Multivector Functionals

Any mapping which sends general extensors over $`V`$ into multivectors over $`V`$ will be called a *general multivector functional over* $`V`$*.*

In particular, the general functionals with image-values belonging to $`^rV`$ are said to be $`r`$*-vector functionals of general extensor.* For the cases $`r=0,`$ $`r=1,`$ $`r=2,\mathrm{}`$ and $`r=n`$ we speak about *scalar, vector, bivector,*$`\mathrm{}`$and *pseudoscalar functionals*, respectively.

For the applications we have in mind we shall need only some particular cases of these general functionals for which we will give special names.

Any mapping $`:ext_p^q(V)^rV`$ will be called a $`r`$*-vector functional of a* $`(p,q)`$*-extensor.* In accordance to what was said above, the cases for which $`[t]`$ belongs to $`,`$ $`V,`$ $`^2V,\mathrm{}`$ and $`^nV`$ will be named respectively as *scalar, vector, bivector,*$`\mathrm{}`$ and *pseudoscalar functionals of a* $`(p,q)`$*-extensor.*

### 2.1 Induced Multivector Functionals

Let $`F:\underset{k\text{ factors}}{\underset{}{^qV\times \mathrm{}\times ^qV}}^rV`$ be any $`r`$-vector function of $`k`$ $`q`$-vector variables. Take some $`k`$-uple of $`p`$-vectors $`(A^1,\mathrm{},A^k).`$

Associated to $`F`$ and with respect to $`(A^1,\mathrm{},A^k)`$ it is possible to construct a $`r`$-vector functional of a $`(p,q)`$-extensor, say $`_{(A^1,\mathrm{},A^k)},`$ given by

$$ext_p^q(V)t_{(A^1,\mathrm{},A^k)}[t]^rV\text{ such that}$$

$$_{(A^1,\mathrm{},A^k)}[t]=F[t(A^1),\mathrm{},t(A^k)].$$
(1)
It will be called the $`r`$*-vector functional of a* $`(p,q)`$*-extensor induced by* $`F,`$ *relative to* $`(A^1,\mathrm{},A^k)`$*.*

If $`F`$ is differentiable on $`\underset{k\text{ factors}}{\underset{}{^qV\times \mathrm{}\times ^qV}},`$ then $`_{(A^1,\mathrm{},A^k)}`$ is said to be *differentially-induced by* $`F`$ *with respect to* $`(A^1,\mathrm{},A^k).`$

In this way if $`_{(A^1,\mathrm{},A^k)}`$ is differentially-induced, then there must exist the standard derivatives of $`F`$ with respect to each $`p`$-vector variable $`X^1,\mathrm{}`$ and $`X^k`$ (the so-called *partial derivatives of* $`F`$), i.e., $`_{X^1}F,\mathrm{}`$ and $`_{X^k}F,`$ (see ).

Associated to $`_{X^1}F,\mathrm{}`$ and $`_{X^k}F`$ with respect to $`(A^1,\mathrm{},A^k)`$ we can define the following multivector functionals of a $`(p,q)`$-extensor:

$`ext_p^q(V)t`$ $`_{X^1}F[t(A^1),\mathrm{},t(A^k)]{\displaystyle V},`$
$`\mathrm{}`$
and

$$ext_p^q(V)t_{X^k}F[t(A^1),\mathrm{},t(A^k)]V.$$
(2)
We see that they are induced by the partial derivatives of $`F`$ with respect to $`(A^1,\mathrm{},A^k).`$

#### 2.1.1 Directional Derivative

Take an arbitrary $`p`$-vector $`A.`$ We introduce the$`A`$*-directional derivative* of the differentially-induced $`r`$-vector functional $`_{(A^1,\mathrm{},A^k)}`$ as being the multivector functional $`_{(A^1,\mathrm{},A^k)A}^{}`$ given by

$$ext_p^q(V)t_{(A^1,\mathrm{},A^k)A}^{}[t]\mathrm{\Lambda }V\text{ such that}$$

$$_{(A^1,\mathrm{},A^k)A}^{}[t]=\underset{i=1}{\overset{𝑘}{}}AA^i_{X^i}F[t(A^1),\mathrm{},t(A^k)].$$
(3)
Note that the algebraic object just defined associated to $`_{(A^1,\mathrm{},A^k)}`$ has the property of linearity with respect to the direction, i.e., for any $`\alpha ,\beta `$ and $`A,B^pV`$

$$_{(A^1,\mathrm{},A^k)\alpha A+\beta B}^{}[t]=\alpha _{(A^1,\mathrm{},A^k)A}^{}[t]+\beta _{(A^1,\mathrm{},A^k)B}^{}[t],$$
(4)
as expected to hold for a well-defined $`A`$-directional derivative of $`_{(A^1,\mathrm{},A^k)}.`$

#### 2.1.2 Derivatives

Let $`(\{e_k\},\{e^k\})`$ be a pair of arbitrary reciprocal bases of $`V.`$ It is also possible to introduce *four* derivatives-like operators for the differentially-induced $`r`$-vector functional $`_{(A^1,\mathrm{},A^k)}`$ as the following multivector functionals $`_{(A^1,\mathrm{},A^k)}^{}`$ defined by

$`_{(A^1,\mathrm{},A^k)}^{}[t]`$ $`={\displaystyle \frac{1}{p!}}(e^{j_1}\mathrm{}e^{j_p})_{(A^1,\mathrm{},A^k)e_{j_1}\mathrm{}e_{j_p}}^{}[t]`$ (5)
$`={\displaystyle \frac{1}{p!}}(e_{j_1}\mathrm{}e_{j_p})_{(A^1,\mathrm{},A^k)e^{j_1}\mathrm{}e^{j_p}}^{}[t],`$ (6)
where $``$ means either $`(),`$ $`(),`$ $`(\mathrm{})`$ or $`(`$*Clifford product*$`).`$

It should be noted that $`_{(A^1,\mathrm{},A^k)}^{}`$ are well-defined multivector functionals of $`(p,q)`$-extensor *only* associated with $`_{(A^1,\mathrm{},A^k)}`$ since, by taking into account eq.(4), $`_{(A^1,\mathrm{},A^k)}^{}[t]`$ are multivectors which do not depend on the choice of $`(\{e_k\},\{e^k\})`$.

Recall also that a straightforward calculation gives with the use eq.(3) that

$`_{(A^1,\mathrm{},A^k)}^{}[t]`$ $`={\displaystyle \frac{1}{p!}}(e^{j_1}\mathrm{}e^{j_p})(\underset{i=1}{\overset{𝑘}{{\displaystyle }}}(e_{j_1}\mathrm{}e_{j_p})A^i_{X^i}F[\mathrm{}])`$
$`=(\underset{i=1}{\overset{𝑘}{{\displaystyle }}}{\displaystyle \frac{1}{p!}}(e_{j_1}\mathrm{}e_{j_p})A^ie^{j_1}\mathrm{}e^{j_p})_{X^i}F[\mathrm{}]`$
$`_{(A^1,\mathrm{},A^k)}^{}[t]`$ $`=\underset{i=1}{\overset{𝑘}{{\displaystyle }}}A^i_{X^i}F[t(A^1),\mathrm{},t(A^k)].`$ (7)
Eq.(7) shows explicitly that $`_{(A^1,\mathrm{},A^k)}^{}`$ can be intrinsically defined without using any pair of reciprocal bases of $`V`$.

The special cases: $`_{(A^1,\mathrm{},A^k)}^{},`$ $`_{(A^1,\mathrm{},A^k)}^{},`$ $`\mathrm{}_{(A^1,\mathrm{},A^k)}^{}`$ and $`_{(A^1,\mathrm{},A^k)}^{}`$ (i.e., $``$*Clifford product*) will be called respectively the *curl, scalar divergence, left contracted divergence* and *gradient* of $`_{(A^1,\mathrm{},A^k)}.`$ Sometimes, $`_{(A^1,\mathrm{},A^k)}^{}`$ will be called the *standard derivative* of $`_{(A^1,\mathrm{},A^k)}.`$

We introduce now on the *real vector space of differentially-induced* $`r`$*-vector functionals of* $`(p,q)`$*-extensor* the following *four* derivative-like operators $`_t`$ as follows

$$_t_{(A^1,\mathrm{},A^k)}[t]=_{(A^1,\mathrm{},A^k)}^{}[t],$$
(8)
i.e., by eq.(7)

$$_t_{(A^1,\mathrm{},A^k)}[t]=\underset{i=1}{\overset{𝑘}{}}A^i_{X^i}F[t(A^1),\mathrm{},t(A^k)].$$
(9)
The special cases: $`_t,`$ $`_t,`$ $`_t\mathrm{}`$ and $`_t`$ (i.e., $``$*Clifford product*) will be called respectively the (functional) *curl, scalar divergence, left contracted divergence* and *gradient operator.* Sometimes, we will say that $`_t`$ is the *standard derivative operator with respect to* $`t`$*.*

$`_t_{(A^1,\mathrm{},A^k)}[t],`$ $`_t_{(A^1,\mathrm{},A^k)}[t],`$ $`_t\mathrm{}_{(A^1,\mathrm{},A^k)}[t]`$ and $`_t_{(A^1,\mathrm{},A^k)}[t]`$ (i.e., $``$*Clifford product*) will be named respectively as the *curl, scalar divergence, left contracted divergence* and *gradient* of$`_{(A^1,\mathrm{},A^k)}.`$ The gradient of $`_{(A^1,\mathrm{},A^k)}`$ will be called the *standard derivative* of $`_{(A^1,\mathrm{},A^k)}`$ *with respect to* $`t`$

It is still possible to define the noticeable derivative-like operator $`A_t`$ as follows

$`A_t_{(A^1,\mathrm{},A^k)}[t]`$ $`=(A{\displaystyle \frac{1}{p!}}e^{j_1}\mathrm{}e^{j_p})_{(A^1,\mathrm{},A^k)e_{j_1}\mathrm{}e_{j_p}}^{}[t]`$ (10)
$`=(A{\displaystyle \frac{1}{p!}}e_{j_1}\mathrm{}e_{j_p})_{(A^1,\mathrm{},A^k)e^{j_1}\mathrm{}e^{j_p}}^{}[t],`$ (11)
i.e., by eq.(4)

$$A_t_{(A^1,\mathrm{},A^k)}[t]=_{(A^1,\mathrm{},A^k)A}^{}[t].$$
(12)
Eq.(12) means that $`A_t`$ is the $`A`$*-directional derivative operator* which maps $`_{(A^1,\mathrm{},A^k)}`$ $``$ $`_{(A^1,\mathrm{},A^k)A}^{}.`$

It is often convenient when doing calculations to employ some abuses of notation for simplifying the handle of the fundamental formulas. Thus, eqs.(3) and (9) will be usually written

$`A_tF[t(A^1),\mathrm{},t(A^k)]`$ $`=\underset{i=1}{\overset{𝑘}{{\displaystyle }}}AA^i_{t(A^i)}F[t(A^1),\mathrm{},t(A^k)],`$ (13)
$`_tF[t(A^1),\mathrm{},t(A^k)]`$ $`=\underset{i=1}{\overset{𝑘}{{\displaystyle }}}A^i_{t(A^i)}F[t(A^1),\mathrm{},t(A^k)].`$ (14)
No confusion arises since $`A_t`$ and $`_t`$ denote derivation of $`r`$-vector functional with respect to $`(p,q)`$-extensor $`t,`$ and $`_{t(A^i)}`$ holds for derivation of $`r`$-vector function with respect to $`q`$-vector $`t(A^i)`$.

It should be noted that by employing the abused notation we can re-write eqs.(5) and (6) as

$`_tF[\mathrm{}]`$ $`={\displaystyle \frac{1}{p!}}(e^{j_1}\mathrm{}e^{j_p})(e_{j_1}\mathrm{}e_{j_p})_tF[\mathrm{}]`$ (15)
$`={\displaystyle \frac{1}{p!}}(e_{j_1}\mathrm{}e_{j_p})(e^{j_1}\mathrm{}e^{j_p})_tF[\mathrm{}].`$ (16)
#### 2.1.3 $`A`$-Directional Derivation Rules

###### Proposition 1.

Take a real $`\lambda `$ and a multivector $`M.`$ If $`tF[t(A^1),\mathrm{},t(A^k)]`$ is any differentially-induced $`r`$-vector functional of a $`(p,q)`$-extensor, then

$`A_t(\lambda F[\mathrm{}])`$ $`=\lambda A_tF[\mathrm{}],`$ (17)
$`A_t(F[\mathrm{}]M)`$ $`=(A_tF[\mathrm{}])M.`$ (18)
###### Proof.

It follows directly from eq.(13) by using the derivation formulas: $`_{X^i}(\lambda F(\mathrm{}))=\lambda _{X^i}F(\mathrm{})`$ and $`_{X^i}(F(\mathrm{})M)=(_{X^i}F(\mathrm{}))M.`$

###### Theorem 1.

Let $`tF[t(A^1),\mathrm{},t(A^k)]`$ and $`tG[t(A^1),\mathrm{},t(A^k)]`$ be any two differentially-induced $`r`$-vector functionals of a $`(p,q)`$-extensor.

The addition $`t(F+G)[t(A^1),\mathrm{},t(A^k)]`$ is a differentially-induced $`r`$-vector functional of a $`(p,q)`$-extensor and the following rule holds

$$A_t(F+G)[\mathrm{}]=A_tF[\mathrm{}]+A_tG[\mathrm{}].$$
(19)
###### Proof.

As we can see, it is an immediate consequence of the derivation rule $`_{X^i}(F+G)(\mathrm{})=_{X^i}(F)(\mathrm{})+_{X^i}G(\mathrm{}).`$

###### Theorem 2.

Let $`t\mathrm{\Phi }[t(A^1),\mathrm{},t(A^k)]`$ and $`tG[t(A^1),\mathrm{},t(A^k)]`$ be any differentially-induced scalar and $`r`$-vector functional of a $`(p,q)`$-extensor, respectively.

The scalar multiplication $`t(\mathrm{\Phi }G)[t(A^1),\mathrm{},t(A^k)]`$ is also a differentially-induced $`r`$-vector functional of a $`(p,q)`$-extensor and we have

$$A_t(\mathrm{\Phi }G)[\mathrm{}]=(A_t\mathrm{\Phi }[\mathrm{}])G[\mathrm{}]+\mathrm{\Phi }[\mathrm{}]A_tG[\mathrm{}].$$
(20)
It is rightly a Leibnitz-like rule.

###### Proof.

As the reader can easily prove, eq.(20) is an immediate consequence of the derivation rule $`_{X^i}(\mathrm{\Phi }G)(\mathrm{})=(_{X^i}\mathrm{\Phi }(\mathrm{}))G(\mathrm{})+\mathrm{\Phi }(\mathrm{})_{X^i}G(\mathrm{}).`$

###### Theorem 3.

Let $`t\mathrm{\Psi }[t(A^1),\mathrm{},t(A^k)]`$ and $`\lambda \varphi (\lambda )`$ be any differentially-induced scalar functional and a derivable ordinary real function, respectively. Then, $`t\varphi (\mathrm{\Psi }[t(A^1),\mathrm{},t(A^k)])`$ is a differentially-induced scalar functional and the following rule holds

$$A_t\varphi (\mathrm{\Psi }[\mathrm{}])=\varphi ^{}(\mathrm{\Psi }[\mathrm{}])A_t\mathrm{\Psi }[\mathrm{}].$$
(21)
It is an interesting and useful chain-like rule for $`A`$-directional derivation of a special type of scalar functionals.

###### Proof.

Eq.(21) follows easily from eq.(13) by taking into account the derivation rule $`_{X^i}\varphi \mathrm{\Psi }(\mathrm{})=\varphi ^{}\mathrm{\Psi }(\mathrm{})_{X^i}\mathrm{\Psi }(\mathrm{}).`$

## 3 Examples

###### Example 4.

Let $`hext_1^1(V)`$ and take $`a,b,cV`$. Then,

$`a_h(h(b)h(c))`$ $`=ab_{h(b)}(h(b)h(c))+ac_{h(c)}(h(b)h(c))`$
$`=abh(c)+ach(b),`$
$`a_h(h(b)h(c))`$ $`=h(abc+acb).`$ (22)
Also,

$`a_h(h(b)h(c))`$ $`=ab_{h(b)}(h(b)h(c))+ac_{h(c)}(h(b)h(c))`$
$`=ab(n1)h(c)ac(n1)h(b)`$
$`=(n1)h(abcacb),`$
$`a_h(h(b)h(c))`$ $`=(n1)h(a\mathrm{}(bc)).`$ (23)
In eqs.(22) and (23) we have used the derivative formulas $`_x(xy)=y`$ and $`_x(xy)=(n1)y`$, where $`n`$ is the dimension of $`V.`$

The second formula developed in this example has an interesting and useful generalization, which is:

The $`a`$-derivative of the $`k`$-vector functional $`ext_1^1(V)h\underset{¯}{h}(a^1\mathrm{}a^k)^kV,`$ with $`a^1,\mathrm{},a^kV,`$ is given by

$$a_h\underset{¯}{h}(a^1\mathrm{}a^k)=(nk+1)\underset{¯}{h}(a\mathrm{}(a^1\mathrm{}a^k)).$$
(24)
###### Example 5.

Let $`hext_1^1(V)`$ and take $`bV.`$

We shall calculate $`a_hh(b)`$ and $`a_hh^{}(b).`$ And, also $`_hh(b)`$ and $`_hh^{}(b).`$

First, we have

$`a_hh(b)`$ $`=ab_{h(b)}h(b)=(ab)n,`$
$`a_hh(b)`$ $`=n(ab),`$ (25)
were we used the derivative formula $`_xx=n.`$ Thus,

$$_hh(b)=e^je_j_hh(b)=e^jn(e_jb)=bn,$$

i.e.,

$`_hh(b)`$ $`=_hh(b)=nb,`$
$`_hh(b)`$ $`=_h\mathrm{}h(b)=0.`$
Now, by employing a trick we have

$$a_hh^{}(b)=a_h(h^{}(b)e^je_j)=a_h(bh(e^j)e_j).$$

Thus, by using eq.(18)

$`a_hh^{}(b)`$ $`=(\underset{i=1}{\overset{𝑛}{{\displaystyle }}}ae^i_{h(e^i)}bh(e^j))e_j=\underset{i=1}{\overset{𝑛}{{\displaystyle }}}ae^ib\delta _i^je_j,`$
$`a_hh^{}(b)`$ $`=ba,`$ (26)
were we used the derivative formula $`_x(bx)=b.`$ Thus,

$$_hh^{}(b)=e^je_j_hh^{}(b)=e^j(be_j).$$

It follows that

$`_hh^{}(b)`$ $`=e^j(be_j)+e^j(be_j)=b.`$
$`_hh^{}(b)`$ $`=e^j(be_j)+e^j(be_j)=0.`$
$`_h\mathrm{}h^{}(b)`$ $`=e^j\mathrm{}(be_j)+e^j\mathrm{}(be_j)=(e^jb)e_j(e^je_j)b=(1n)b.`$
$`_hh^{}(b)`$ $`=e^j(2e_jbe_jb)=(2n)b.`$
###### Example 6.

Let $`text_1^1(V).`$ The trace of $`t,`$ i.e., $`ttr[t]=t(e^j)e_j,`$ is a scalar functional and the bivector of $`t,`$ i.e., $`tbiv[t]=t(e^j)e_j,`$ is a bivector functional, both of them associated to $`t`$. We shall calculate $`a_ttr[t]`$ and $`a_tbiv[t].`$ And, also $`_ttr[t]`$ and $`_tbiv[t].`$

First, we have

$`a_ttr[t]`$ $`=\underset{i=1}{\overset{𝑛}{{\displaystyle }}}ae^i_{t(e^i)}(t(e^j)e_j)=\underset{i=1}{\overset{𝑛}{{\displaystyle }}}ae^i\delta _i^je_j,`$
$`a_ttr[t]`$ $`=a.`$ (27)
We have used once again the derivative formula $`_x(xy)=y.`$ Hence,

$$_ttr[t]=e^je_j_ttr[t]=e^je_j,$$

i.e.,

$`_ttr[t]`$ $`=0,`$
$`_ttr[t]`$ $`=_t\mathrm{}tr[t]=_ttr[t]=n.`$
Now, we have also

$`a_tbiv[t]`$ $`=\underset{i=1}{\overset{𝑛}{{\displaystyle }}}ae^i_{t(e^i)}(t(e^j)e_j)=\underset{i=1}{\overset{𝑛}{{\displaystyle }}}ae^i(n1)\delta _i^je_j,`$
$`a_tbiv[t]`$ $`=(n1)a,`$ (28)
were we have used once again the derivative formula $`_x(xy)=(n1)y.`$ Hence,

$$_tbiv[t]=(n1)e^je_j,$$

i.e.,

$`_tbiv[t]`$ $`=0,`$
$`_tbiv[t]`$ $`=_t\mathrm{}biv[t]=_tbiv[t]=(n1)n.`$
###### Example 7.

Let $`hext_1^1(V)`$ and take a non-zero $`I^nV.`$ We shall calculate the $`a`$-directional derivative of the pseudoscalar functional $`h\underset{¯}{h}(I),`$ i.e., $`a_h\underset{¯}{h}(I).`$

By employing one of the expansion formulas for pseudoscalars (see ), eq.(17) and eq.(24) we have

$`a_h\underset{¯}{h}(I)`$ $`=a_hI(e_1\mathrm{}e_n)\underset{¯}{h}(e^1\mathrm{}e^n)`$
$`=I(e_1\mathrm{}e_n)a_h\underset{¯}{h}(e^1\mathrm{}e^n)`$
$`=I(e_1\mathrm{}e_n)\underset{¯}{h}(a\mathrm{}(e^1\mathrm{}e^n)),`$
$`a_h\underset{¯}{h}(I)`$ $`=\underset{¯}{h}(a\mathrm{}I)=h(aI).`$ (29)
###### Example 8.

Let $`hext_1^1(V)`$ and take a non-zero $`I\mathrm{\Lambda }^nV.`$ The determinant of $`h,`$ i.e., $`hdet[h]`$ such that $`\underset{¯}{h}(I)=det[h]I,`$ is a characteristic scalar functional of $`h.`$ We shall calculate $`a_hdet[h]`$ and $`_hdet[h].`$

By employing eq.(18) and eq.(29) we have

$$a_hdet[h]=(a_h\underset{¯}{h}(I))I^1=\underset{¯}{h}(aI)I^1.$$

But, by taking into account the extensor formula $`h^1(a)=det^1[h]\underset{¯}{h}^{}(aI)I^1`$ (see) and recalling that $`det[h^{}]=det[h]`$ and $`h^{}=(h^{})^1=(h^1)^{}`$ we get

$$a_hdet[h]=det[h]h^{}(a).$$
(30)
Hence, it follows that

$$_hdet[h]=e^je_j_hdet[h]=det[h]e^jh^{}(e_j),$$

i.e.,

$`_hdet[h]`$ $`=det[h]h^{}(e_j)e^j=det[h]biv[h^1].`$
$`_hdet[h]`$ $`=_h\mathrm{}det[h]=det[h]h^1(e^j)e_j=det[h]tr[h^1].`$
$`_hdet[h]`$ $`=det[h]e^jh^{}(e_j)=det[h](tr[h^1]+biv[h^1]).`$
### 3.1 An Enlightening Discussion

Let us consider for example a differentially-induced scalar functional of $`(1,1)`$-extensor $`t\mathrm{\Phi }[t(a^1)].`$ We have the possibility for constructing a differentiable scalar function of $`n\times n`$ real variables $`(t_{11},\mathrm{},t_{1n},\mathrm{},t_{n1},\mathrm{},t_{nn})\widehat{\mathrm{\Phi }}(t_{11},\mathrm{},t_{1n},\mathrm{},t_{n1},\mathrm{},t_{nn}),`$ defined by

$$\widehat{\mathrm{\Phi }}(t_{11},\mathrm{},t_{1n},\mathrm{},t_{n1},\mathrm{},t_{nn})=\mathrm{\Phi }[t_{ij}(a^1e^i)e^j]$$
(31)
where $`t_{ij}=t(e_i)e_j`$ are the $`n\times n`$ matrix elements of $`t`$ with respect to $`\{e_k\}.`$

Eq.(31) shows that all information just contained into the *classical real function* $`(t_{11},\mathrm{},t_{nn})\widehat{\mathrm{\Phi }}(t_{11},\mathrm{},t_{nn})`$ whose real variables are $`t_{pq},`$ is also codified into the scalar functional $`t\mathrm{\Phi }[t(a^1)].`$

We shall search for the relationship which exists between the *ordinary partial derivatives* of $`\widehat{\mathrm{\Phi }}(\mathrm{})`$ with respect to each *tensor covariant component*222They are the $`n\times n`$ covariant components of a $`2`$-tensor $`T`$ in biunivocal correspondence with the $`(1,1)`$-extensor $`t,`$ see , i.e., $`T_{pq}T(e_p,e_q)=t(e_p)e_qt_{pq}.`$ $`t_{pq}`$ and the $`a`$-directional derivative of $`\mathrm{\Phi }[\mathrm{}].`$

By using $`_{\lambda _i}\mathrm{\Phi }(x(\lambda _1,\mathrm{},\lambda _k))=_{\lambda _i}x(\lambda _1,\mathrm{},\lambda _k)_x\mathrm{\Phi }(x(\lambda _1,\mathrm{},\lambda _k)),`$ a chain-like derivation rule, we may write

$`{\displaystyle \frac{\widehat{\mathrm{\Phi }}}{t_{pq}}}(t_{11},\mathrm{},t_{nn})`$ $`=_{t_{pq}}(t_{ij}(a^1e^i)e^j)_x\mathrm{\Phi }[t_{ij}(a^1e^i)e^j]`$
$`=\delta _{ij}^{pq}(a^1e^i)e^j_{t(a^1)}\mathrm{\Phi }[t(a^1)],`$
$`{\displaystyle \frac{\widehat{\mathrm{\Phi }}}{t_{pq}}}(t_{11},\mathrm{},t_{nn})`$ $`=(a^1e^p)e^q_{t(a^1)}\mathrm{\Phi }[t(a^1)].`$ (32)
Now, Clifford multiplication by $`(ae_p)e_q`$ (and summing over $`p,q`$) on both sides of eq.(32) yields

$`(ae_p)e_q{\displaystyle \frac{\widehat{\mathrm{\Phi }}}{t_{pq}}}(t_{11},\mathrm{},t_{nn})`$ $`=(aa^1)e_qe^q_{t(a^1)}\mathrm{\Phi }[t(a^1)]`$
$`=aa^1_{t(a^1)}\mathrm{\Phi }[t(a^1)],`$
$`(ae_p)e_q{\displaystyle \frac{\widehat{\mathrm{\Phi }}}{t_{pq}}}(t_{11},\mathrm{},t_{nn})`$ $`=a_t\mathrm{\Phi }[t(a^1)].`$ (33)
That is the required result relating both $`{\displaystyle \frac{\widehat{\mathrm{\Phi }}}{t_{pq}}}(t_{11},\mathrm{},t_{nn})`$ and $`a_t\mathrm{\Phi }[t(a^1)].`$

It is still possible to find a relationship between $`{\displaystyle \frac{\widehat{\mathrm{\Phi }}}{t_{pq}}}(t_{11},\mathrm{},t_{nn})`$ and the $``$-derivatives of $`\mathrm{\Phi }[t(a^1)].`$ From eq.(32) we have

$`e_p(e_q{\displaystyle \frac{\widehat{\mathrm{\Phi }}}{t_{pq}}}(t_{11},\mathrm{},t_{nn}))`$ $`=a^1(e_qe^q_{t(a^1)}\mathrm{\Phi }[t(a^1)])`$
$`=a^1_{t(a^1)}\mathrm{\Phi }[t(a^1)],`$
$`e_p(e_q{\displaystyle \frac{\widehat{\mathrm{\Phi }}}{t_{pq}}}(t_{11},\mathrm{},t_{nn}))`$ $`=_t\mathrm{\Phi }[t(a^1)].`$ (34)
That is the expected identity which relates both $`{\displaystyle \frac{\widehat{\mathrm{\Phi }}}{t_{pq}}}(t_{11},\mathrm{},t_{nn})`$ and $`_t\mathrm{\Phi }[t(a^1)].`$

## 4 Conclusions

In this paper we introduced the key concepts of a theory of multivector functionals. We studied several aspects of the notion of derivative that can be applied to these objects, as e.g., the $`A`$-directional derivatives and the generalized concepts of curl, divergence and gradient. We worked in details several examples where we calculate different types of derivatives for multivector functionals. It is worth to said once again that these objects play a decisive role in the development of a Lagrangian formalism for extensor fields as it will be seen in two future series of papers: *geometric theories of gravitation and Lagrangian formulation of the multivector and extensor fields theory*.

Acknowledgement: V. V. Fernández is grateful to FAPESP for a posdoctoral fellowship. W. A. Rodrigues Jr. is grateful to CNPq for a senior research fellowship (contract 201560/82-8) and to the Department of Mathematics of the University of Liverpool for the hospitality. Authors are also grateful to Drs. P. Lounesto, I. Porteous, and J. Vaz, Jr. for their interest on our research and useful discussions.