text:
# COLOR GRADIENTS IN EARLY-TYPE GALAXIES IN CLUSTERS AT THE REDSHIFT FROM 0.37 TO 0.56

## 1 INTRODUCTION

It has been known that nearby elliptical galaxies have color gradients; colors in an elliptical galaxy gradually become bluer with increasing radius (e.g., Vader et al. 1988; Franx, Illingworth, & Heckman 1989; Peletier et al. 1990a; Peletier, Valentijn, & Jameson 1990b, Goudfrooij et al. 1994; Michard 1999). Since many of elliptical galaxies show radial gradients in metal absorption line strengths such as Mg2, Fe1(5270 Å) and Fe2(5335 Å) (e.g., Carollo, Danziger, & Buson 1993; Davies, Sadler, & Peletier 1993; Gonzalez 1993; Kobayashi & Arimoto 1999), the origin of the color gradients has been naively interpreted to be the metallicity gradients.

However, such an interpretation for the origin of the color gradient is premature, because both metallicity gradient and age gradient in stellar population can cause the same color gradient, and we cannot distinguish the cause for the gradient. This is called age-metallicity degeneracy, which was originally pointed out by Worthey, Trager, & Faber (1996) in terms of the origin of the color-magnitude relation of nearby elliptical galaxies (see also Arimoto 1996). In order to break this degeneracy and to know the primary origin of the color gradients in elliptical galaxies, comparing the observed color gradients in distant ellipticals with predicted model gradients caused by either the metallicity gradient or the age gradient is a very effective approach, as was successful for examining the origin of the color-magnitude (CM) relation (Kodama & Arimoto, 1997). Tamura et al. (2000; hereafter called Paper I) constructed the two models both of which reproduce a typical color gradient of elliptical galaxies at $`z=0`$ using a population synthesis model. In one model, the mean metallicity of the stellar population decreases with increasing radius at a fixed old mean age. While in the other one, the mean age decreases with a radius at a fixed mean metallicity. These models were then made evolve back in time. The evolution of color gradients thus predicted are confronted with the observed ones in distant ($`z=0.11.0`$) ellipticals sampled from the Hubble Deep Field-North (HDF-N; Williams et al. 1996). As a result, Paper I found that the metallicity gradient is the primary origin of color gradients and the age gradient model cannot reproduce the observed gradient at such redshift.

The elliptical galaxies in the HDF-N, however, are only those in field environment. It has never been obvious that ellipticals in clusters evolve similarly as those in field. In rich clusters, it has been found that the color-magnitude relation still holds even at around $`z1`$ (e.g., Stanford, Eisenhardt, & Dickinson 1998) and these observational results seem to favor the classical monolithic collapse scenario associated with the galactic wind and high-$`z`$ formation (e.g., $`z>3`$) of elliptical galaxies (e.g., Kodama et al. 1998). However, this kind of evolution has not been established for ellipticals in lower density environment (but see Kodama, Bower, & Bell 1998). Some predictions either theoretically or observationally show that field ellipticals formed by recent (at $`z1`$) merging processes (e.g., Baugh, Cole, & Frenk 1996; Barger et al. 1999). An internal structure of a galaxy such as a metallicity gradient and an age gradient must depend on its formation process. If cluster ellipticals pass different formation histories from those for field ellipticals, their internal structures, thus the origin of the color gradients, may not be the same. Or some environmental effects on color gradients may exist. Thus, the same approach is needed for cluster ellipticals to clarify the origin of their color gradients.

It is noted that dust extinction in elliptical galaxies may also have some effects on the color gradients (Goudfrooij & de Jong 1995; Wise & Silva 1996; Silva & Wise 1996). However, about half of the detection towards ellipticals in far infrared with IRAS are around $`3\sigma `$ threshold and confirmation is needed to be definitive (Bregman et al. 1998). In addition, spatial distribution of dust in an elliptical galaxy as well as dust mass which could affect a color gradient are not established yet. These are still open problems and will be examined in detail in our forthcoming papers. Therefore, in this paper, we have chosen to focus on age and metallicity effects only.

This paper is organized as follows. The sample selection and data analysis of elliptical galaxies are described in § 2. Histograms of color gradients are presented in § 3 together with the representative color profiles of the sample ellipticals. Discussion is given in § 4. The cosmological parameters adopted throughout this paper are the same as those in Paper I; $`H_0=50`$ km s-1 Mpc-1, $`\mathrm{\Omega }_0=0.2`$ and $`\mathrm{\Lambda }=0`$.

## 2 DATA AND SAMPLE SELECTION

In order to examine color gradients in elliptical galaxies in distant clusters, deep imaging data with a high angular resolution are necessary in more than two bands. Thus, we choose to use the archival data taken with the WFPC2 on-board the HST. Smail et al. (1997) obtained deep imaging data of 11 distant rich clusters at the redshift from 0.37 to 0.56, most of which have exposure times more than 10,000 sec to examine their morphology in detail, and presented catalogs of the detected objects. In this paper, their reduced imaging data which are available at their web site111http://star-www.dur.ac.uk/~irs/morphs.html are used. From these data, we select six clusters whose images were taken in two bands to obtain galaxy colors. Cluster IDs and their basic properties taken from Smail et al. (1997) are given in Table 1. Among these clusters, Cl 0024$`+`$16, Cl 0016$`+`$16, and Cl 0054$``$27 are classified as high concentration clusters and Cl 0939$`+`$47 and Cl 0412$``$65 are low concentration clusters (Dressler et al. 1997). For A370 and Cl 0939$`+`$47, the images of their outer fields were taken. Therefore, environmental effects on color gradients may be examined.

Early type galaxies (E, E/S0, S0/E) in these clusters are sampled based on the catalog by Smail et al. (1997). Since our main interest is to examine their color gradients, high signal-to-noise ratio (S/N) is required and thus galaxies brighter than the apparent magnitude of 21 mag in $`I_{814}`$ band are selected in all of the sample clusters. This apparent magnitude roughly corresponds to the absolute magnitude of $`20`$ mag in $`V`$ band. Our sample galaxies are listed in Table 2 with their basic parameters, most of which are taken from Smail et al. (1997). In each cluster, the sample galaxies form the tight CM relations as shown in Figure 1, though the relations in A370, Cl 0939$`+`$47, and Cl 0412$``$65 are somewhat loose. Since Cl 0939$`+`$47 and Cl 0412$``$65 are classified as low concentration clusters, and the outer regions of the clusters were imaged for A370 and Cl 0939$`+`$47, some environmental effects on the CM relations might be seen. However, a discussion for this subject is beyond our scope in this paper. In the figure, a solid square with an open circle shows the object removed from the following analysis, because their color gradients could not be derived accurately enough to use our discussion due to their low signal-to-noise ratios (S/Ns) or due to their close positions to an edge of the image.

## 3 PROFILES AND COLOR GRADIENTS

In deriving color profiles of the sample galaxies, we should take care of a positional coincidence between a blue image and a red image, and a difference of the sizes of point spread functions (PSFs) between the two images. (One pixel corresponds to $`0_{}^{\prime \prime }10`$ for the Wide Field Cameras and $`0_{}^{\prime \prime }046`$ for the Planetary Camera.) Positional displacement between a blue image and a red one causes spurious asymmetry of a color distribution in a galaxy and has a serious effect on the color gradient. We examined the centroids of the sample galaxies in both blue and red images and corrected a displacement, if a systematic shift between the two images in each cluster was found. The sizes of PSFs should be also estimated and be adjusted between blue and red images. We simulated the PSF at each position of a sample galaxy in each band by using TinyTim v4.4 (Krist 1995; Krist & Hook 1997), and corrected the difference of the sizes between the two images, which were estimated by fitting a Gaussian to the PSFs. Sky value was determined by “mode” in an annulus with an inner radius of $`9^{\prime \prime }`$ and a width of 3′′ by using the “phot” task in the IRAF apphot package and the obtained sky was subtracted.

We next made azimuthally averaged radial surface brightness profiles in both blue and red images with a radial sampling of 0.′′1. These profiles in each galaxy are made along the ellipses fitted with a fixed position angle to the isophotes of the galaxy images in the red band. (The position angles are taken from Table 2A by Smail et al. (1997).) Thus the surface brightness profiles in the two bands are derived based on the same fitted elliptical isophotes. In Figure 2, representative surface brightness profiles of our sample galaxies are shown. The profiles of the brightest, intermediate, and the faintest sample galaxies in each cluster are shown in top, middle, and bottom panels, respectively. The profiles of almost all the galaxies in our sample are well represented by the $`r^{1/4}`$ law. To derive color profiles, the surface brightness profile in the red band is subtracted from the profile in the blue band. Figure 3 shows the resulting color profiles with the same arrangement as in Figure 2. An error bar to each data point includes a photometric error, a local sky subtraction error (1$`\%`$ of a sky value around each object is adopted), and a dispersion of colors along each elliptical isophote. It is important to describe here the two dimensional color distribution in the sample galaxies. The color maps were constructed by dividing the blue image by the red image after adjusting the positional displacement and difference of the sizes of the PSFs. Almost all the color maps show smooth color distribution and do not show any asymmetric or peculiar features. Thus the color profiles well represent the color distribution in each galaxy.

Finally, slopes of the color profiles, i.e., color gradients, in the sample galaxies are derived by the least square fitting within effective radii, which are obtained from the surface brightness profiles in the red band assuming the $`r^{1/4}`$ law. The fitting for the color gradients are done after removing data points with errors larger than 0.3 mag. We do not derive the gradients of the profiles with accepted data points fewer than 3. The average number of the data points used for the fitting is about 8 and the rejected objects are minor. Resulting color gradients in the sample galaxies are listed in Table 2 with the 1$`\sigma `$ fitting errors. Figure 3 indicates resulting slopes of the color gradients as well as the slopes with $`\pm 1\sigma `$. In this figure, abscissa refers to logarithm of a radius normalized by the outermost radius ($`r_f`$) used for the fitting in each color profile. For most of the sample galaxies, $`r_f`$ is roughly equal to $`r_e`$.

Figure 4 shows a histogram for the color gradients of the sample galaxies in each cluster. Each bin of the histograms is set to be 0.2 mag/dex which is comparable to the average value of the fitting error. It is found that the distributions of the gradients are very narrow except for a few outliers, which are #535 and #738 in Cl 0024$`+`$16, #2005 in Cl 0939$`+`$47, and #2050 in Cl 0016$`+`$16. The former three are significantly out of the CM relation towards the blue side, but the last one is almost on the relation. Considering that the range of the distribution of the color gradients is comparable with or only slightly larger than the estimated error for the slopes, the intrinsic dispersions of the color gradients must be considerably small. (The dispersion of the color gradients in nearby elliptical galaxies is about 0.04 mag/dex (Peletier et al. 1990a; Paper I).) It is intriguing that the color gradients of elliptical galaxies are uniform even at intermediate redshift. Furthermore, this encourages the comparison between model gradients and observed gradients in distant clusters, despite rather large errors for the observed slopes.

## 4 ORIGIN OF COLOR GRADIENTS IN CLUSTER ELLIPTICALS

### 4.1 Models

In order to examine whether the origin of the color gradient is the stellar metallicity or the age, we adopt the same approach as that in Paper I and the reader should refer to it in detail. We briefly summarize it here.

An observed color gradient can be reproduced by either a metallicity gradient or an age gradient of stellar populations in an elliptical galaxy at $`z=0`$. However, since the color gradient caused by a metallicity gradient is expected to follow a different evolution from that by an age gradient, the origin of the color gradients can be found by comparing the observed ones at high redshift with those predicted by model. For this purpose, using the population synthesis model (Kodama & Arimoto 1997; Kobayashi, Tsujimoto, & Nomoto 2000), we construct the two model galaxies; one model galaxy has the color gradient made by the pure metallicity gradient (thereafter called metallicity gradient model) without age gradient, and the other made by the pure age gradient (age gradient model) without metallicity gradient. In the metallicity gradient model, the metallicity gradient is produced by assuming that a galactic wind blowed later in the inner region in an elliptical galaxy; star formation continued longer and thus the mean stellar metallicity became higher at the inner region. For the age gradient model, star formation started earlier in the inner region and thus the mean age of stellar populations is older than that in the outer region. The stellar population in each region in an elliptical galaxy is assumed to be made by a single burst and to evolve independently of other regions. Model parameters used here are set to the same as those in Paper I, which are chosen so as to reproduce the typical color gradient at $`z=0`$. The mean value of $`\mathrm{\Delta }(BR)/\mathrm{\Delta }\mathrm{log}r`$ of $`0.09`$ mag$`/`$dex obtained by Peletier et al. (1990a) is adopted as the typical color gradient at $`z=0`$. Note that these model galaxies must be old (8 $``$ 15 Gyr) to reproduce colors in nearby elliptical galaxies. Then, we calculate the spectral evolution in each region of the model galaxies and their color gradients at any redshifts using the response functions including each filter on the HST. It should be emphasized that we do not intend to study physical formation process of elliptical galaxies in this paper, but aim at depicting the evolution of the color gradient caused by either metallicity gradient or age gradient to be compared with the observed ones. Actual physical process that made the metallicity/age gradient may be different from our brief ideas presented in the model description. However it is not a problem here, because once such gradient formed, subsequent evolution of the stellar population is unique and does not depend on the formation process.

The two lines in each panel of Figure 5 show the evolutionary track of the model color gradients; the solid curve indicates the evolution for the metallicity gradient model and the dotted curve for the age gradient model. The model color gradient by the metallicity gradient is almost constant with a redshift within $`z1`$, while that by the age gradient changes abruptly and shows a quite steep gradient even at $`z=0.3`$. We will compare the model gradients with the observed ones in the next subsection.

### 4.2 Model vs observation

The mean values of the color gradients in each cluster sample are plotted at their redshifts in Figure 5. An error bar attached to each point indicates a mean error of the gradients in each cluster. As clearly shown, the metallicity gradient is much more favorable as the origin of the color gradients. This result does not depend on cosmological parameters or parameters for an evolutionary model of galaxy within a reasonable range, and does not change even if we consider the dispersion of the color gradients in the sample galaxies (see Figure 4) and that in nearby ellipticals ($`0.04`$ mag/dex). Although the sample galaxies of which memberships in the clusters are spectroscopically confirmed are minor (Dressler et al. 1999), background or foreground contaminations are not expected to affect the result for the origin of the color gradients, because the result does not change even if we remove the galaxies which significantly deviate from the CM relation in each cluster. The color gradients in several sample galaxies may be affected by other galaxies close to them, and the color profile of a galaxy which locates close to an edge of the chip or on a joint between the cameras may be somewhat spurious. However, our result still holds even after removing the galaxies which may suffer from these effects.

Considering the result in Paper I, in both cluster and field, the primary origin of the color gradients in elliptical galaxies is considered to be the stellar metallicity. However, it is interesting to point out that the mean values of the color gradients seem to deviate upwards from the line for the metallicity gradient model. Our models are calibrated by the color gradients seen in nearby ellipticals by Peletier et al. (1990a), in which most of the sample ellipticals reside in field or group environment. Therefore, the upward deviation might indicate an environmental effect on the color gradients of elliptical galaxies between in rich clusters and in field. However, the correlation between the mean value and the degree of the concentration in each cluster is not seen. In addition, the mean gradients of the clusters of which outer field images were taken do not show larger values than others. Further detailed study on the color gradients in cluster ellipticals and field ones at $`z=0`$ as well as at high redshift should be done in the future.

We would like to thank C. Kobayashi, N. Arimoto, and T. Kodama for fruitful collaboration in Paper I. This work was financially supported in part by Grant-in-Aid for the Scientific Research (No. 11740123) by the Ministry of Education, Science, Sports and Culture of Japan.

Figure Caption

Figure. 1 — Color-magnitude (CM) diagrams for the sample galaxies in the clusters. A solid square with an open circle indicates the object whose color gradient cannot be obtained due to low S/N or the close position to an edge of the image. “#2” in the cluster ID refers to the outer field of the cluster.

Figure. 2 — Azimuthally averaged radial surface brightness profiles of representative galaxies in the sample are presented. The profiles of the brightest, intermediate, and the faintest sample galaxies in each cluster are shown in top, middle, and bottom panels, respectively. Solid squares are the profile in a blue band and open circles in a red band. An object ID is shown at the upper left in each panel.

Figure. 3 — Color profiles and fitted slopes of the representative galaxies in the sample are indicated. The galaxies presented and the arrangement of the panels are the same as those in Figure 2. Among the three solid lines in each panel, the middle one shows the best fit slope. Other two lines show the slopes with $`\pm 1\sigma `$ of the best fit. Abscissa refers to logarithm of a radius normalized by the outermost radius ($`r_f`$) used for the fitting in each color profile.

Figure. 4 — Histograms of the color gradients in the sample galaxies in each sample cluster. A total number of the objects in each histogram is shown at the upper left in each panel and a number in parentheses indicates a total number of the sample galaxies in each cluster. Cluster ID and its redshift are shown on the top of each panel.

Figure. 5 — A mean color gradient in each cluster versus redshift is shown. A solid curve in each panel represents the evolutionary track of the color gradient caused by the metallicity gradient and a dotted curve shows the track by the age gradient (see text in detail).


format: text

text:
# Continuous time dynamics of the Thermal Minority Game

## Abstract

We study the continuous time dynamics of the Thermal Minority Game. We find that the dynamical equations of the model reduce to a set of stochastic differential equations for an interacting disordered system with non-trivial random diffusion. This is the simplest microscopic description which accounts for all the features of the system. Within this framework, we study the phase structure of the model and find that its macroscopic properties strongly depend on the initial conditions.

Many of the current challenges for statistical physics have their origins in problems in biology  and economics . In particular, the application of ideas and techniques of the statistical mechanics of disordered systems can prove useful in the study of systems of adaptive and competitive agents, which are relevant, for example, to the microscopic modeling of financial markets; and, conversely, such problems can raise new issues for statistical physics. One of these systems is the minority game (MG) , a simple model based on Arthur’s “El Farol” bar problem  which describes the behaviour of a group of competing heterogeneous agents subject to the economic law of supply and demand. Despite its simplicity, the MG is very non-trivial, and although much progress has been made in the qualitative  and quantitative  understanding of its features, a full analytic solution of the MG is still missing.

The main hurdles in the way of an analytical study of the MG in its original formulation were its non-locality in time due to the dependence on the game history, its discrete kinematics and dynamics, and the “best-strategy” rule (see, however, ). The first of these obstacles was overcome in , where it was shown numerically that the macroscopic behaviour of the MG was unchanged if the real history was replaced by a random one. This allowed the study of a simpler stochastic Markovian problem instead of the original deterministic non-Markovian one.

In  a natural continuous generalization of the MG was presented. The “information” to which the agents reacted was taken as an external input to the system and it was shown that all the macroscopic features of the MG were preserved, as long as the external information was ergodic in time, the simplest choice being just noise. To handle the problem of the ‘best-strategy’ rule, the Thermal Minority Game (TMG) was introduced, in which a certain degree of stochasticity in the choice of the strategies by the agents was allowed, controlled by a parameter $`T`$, the “temperature”, the limits $`T=0`$ and $`T=\mathrm{}`$ corresponding to the original deterministic MG and the case of completely random strategy choices, respectively. The TMG displayed extra non-trivial structure as a function of $`T`$, notably that in the region where the MG performs worse than random, the system can be made to perform better than random by allowing a certain degree of individual stochastic error.

In the present paper we carefully study the continuous time limit of the TMG, in order to obtain the simplest microscopic description which accounts for all the macroscopic features of the system, and as a further step towards an analytical solution of the model. We confirm that the external information ‘integrates-out’ providing simply an effective coupling between agents. We also show the crucial dependence of the macroscopic properties of the model on the initial conditions. We find that the microscopic equations of the TMG reduce to a set of disordered stochastic differential equations with a non-trivial random diffusion matrix, and study the phase structure of the model in the temperature/dimension plane.

The setup of the TMG is as follows . The system consists of $`N`$ agents playing the game. At each time step $`t`$, each agent reacts to a common piece of “information” $`\stackrel{}{I}(t)`$, by making a “bid” $`b_i(t)`$ ($`i=1,\mathrm{},N`$). The information, defined as a unit-length vector in $`^D`$, is taken to be a random noise, $`\delta `$-correlated in time and uniformly distributed on the unit sphere . The bid $`b(t)`$ is defined to be a real number, which can be interpreted as placing an order in a market, of size $`|b(t)|`$ and positive/negative meaning buy/sell. Bids are prescribed by “strategies”: maps from information to bid, $`^D`$. For simplicity the strategy space $`\mathrm{\Gamma }`$ of the model is restricted to the subspace of homogeneous linear mappings. A strategy $`\stackrel{}{R}`$ is defined as a vector in $`^D`$, subject to the constraint $`\stackrel{}{R}=\sqrt{D}`$, and the prescribed bid is just the scalar product $`\stackrel{}{R}\stackrel{}{I}(t)`$. Each agent has $`S`$ strategies, drawn randomly and independently from $`\mathrm{\Gamma }`$ (with uniform distribution) remaining fixed throughout the game. In what follows we will restrict for simplicity to the case of two strategies per agent $`S=2`$, the generalization to $`S>2`$ being straightforward. At time step $`t`$ each agent $`i`$ chooses one of his/her strategies $`\stackrel{}{R}_i^{}(t)`$ to play with. The “total bid” (or “excess demand”) is then $`A(t)_jb_j(t)=_j\stackrel{}{R}_j^{}(t)\stackrel{}{I}(t)`$. The agents keep track of the potential success of the strategies by assigning points to them, which are updated according to $`P(\stackrel{}{R},t+1)=P(\stackrel{}{R},t)A(t)b(\stackrel{}{R})/N`$, where $`P(\stackrel{}{R},t)`$ represents the points of strategy $`\stackrel{}{R}`$ at time $`t`$.

In the original formulation of the MG the agents played in a deterministic fashion using their ‘best’ strategies, the ones with the highest number of points. In the TMG the natural generalization to non-deterministic behaviour is allowed. At time step $`t`$, each agent $`i`$ chooses $`\stackrel{}{R}_i^{}(t)`$ randomly from his/her $`\{\stackrel{}{R}_i^1,\stackrel{}{R}_i^2\}`$ with probabilities $`\{\pi _i^1(t),\pi _i^2(t)\}`$. The probabilities $`\pi _i^a(t)`$ are functions of the points parameterized by a temperature $`T`$, defined so as to interpolate between the MG case at $`T=0`$, all the way up to the totally random case $`\pi _i^1=\pi _i^2=1/2`$ at $`T=\mathrm{}`$. The qualitative behaviour of the system does not depend on the specific functional form of the $`\pi _i^a(t)`$. In  the probabilities were defined as $`\pi _i^{1,2}(t)\mathrm{exp}[\beta P(\stackrel{}{R}_i^{1,2},t)]`$ (with $`\pi _i^1(t)+\pi _i^2(t)=1`$ and $`\beta =1/T`$), while an alternative convenient form is given by $`\pi _i^{1,2}(t)\mathrm{exp}[\pm \beta z_i(t)]`$, where $`z_i(t)\mathrm{s}gn(p_i(t))`$, with $`p_i(t)[P(\stackrel{}{R}_i^1,t)P(\stackrel{}{R}_i^2,t)]/2`$. The consequential difference between these two definitions will be discussed below.

The set of unconstrained degrees of freedom of the TMG is given by the difference $`p_i(t)`$ of the points of the two strategies of each agent. The choice of strategies can then be defined by $`\stackrel{}{R}_i^{}(t)=\stackrel{}{h}_i+\stackrel{}{\xi }_i\mathrm{sgn}\left[s_i(t)+\mu (t)\right]`$, where $`s_i(t)\pi _i^1(t)\pi _i^2(t)`$, $`\stackrel{}{\omega }_i\left(R_i^1+R_i^2\right)/2`$, $`\stackrel{}{\xi }_i\left(R_i^1R_i^2\right)/2`$, and $`\mu (t)`$ is a stochastic random variable uniformly distributed between $`1`$ and $`1`$ and independently distributed in time. The equations for the point differences then read,

$$p_i(t+1)=p_i(t)\stackrel{}{a}(t)\stackrel{}{I}(t)\stackrel{}{\xi }_i\stackrel{}{I}(t),$$
(1)
where $`\stackrel{}{a}(t)_i\stackrel{}{R}_i^{}(t)/N`$. Eqs. (1), together with the random processes for $`\stackrel{}{I}(t)`$ and $`R_i^{}(t)`$, define the dynamics of the TMG.

The average of the total bid $`A(t)`$ over time and quenched disorder is zero, so the first relevant macroscopic observable of the TMG is its normalized standard deviation $`\sigma `$ (or “volatility”) $`\sigma ^2N^1\overline{A^2(t)}`$, where the overline means disorder average, and $`lim_t\mathrm{}\frac{1}{t}_{t_0}^{t_0+t}()𝑑t^{}`$. In  it was found that $`\sigma `$ had a non-trivial structure both as a function of the reduced dimension of the strategy space $`d=D/N`$ and of the temperature. The second important observable is the fraction $`\varphi `$ of “frozen” agents, defined as those for which one of the strategies always outperforms the other, $`\varphi N^1_i\delta \left(\left|z_i(t)\right|1\right)`$, with the normalization $`\delta (0)=1`$. It was introduced in  as an order parameter for the MG, where it was found that $`\varphi (d)`$ changed from zero to a finite value at $`d=d_c`$.

We now consider the continuous time limit of Eqs. (1) in such a way as to preserve all the macroscopic features of the TMG. To this end we introduce an arbitrary time step $`\mathrm{\Delta }t`$. We deal first with the information $`\stackrel{}{I}(t)`$. Let us assume that $`\stackrel{}{I}(t)`$ is a differential random motion in the space of strategies, i.e., $`\stackrel{}{I}(t)=\mathrm{\Delta }\stackrel{}{W}(t)`$, with zero mean and variance $`\mathrm{\Delta }t`$. Replacing in Eqs. (1) we obtain $`p_i(t+\mathrm{\Delta }t)=p_i(t)\stackrel{}{a}(t)\mathrm{\Delta }\stackrel{}{W}(t)\stackrel{}{\xi }_i\mathrm{\Delta }\stackrel{}{W}(t)`$. In the limit $`\mathrm{\Delta }t0`$, and using the Kramers-Moyal expansion , we get

$$dp_i(t)=\stackrel{}{a}(t)\stackrel{}{\xi }_idt/D+𝒪(dt^2).$$
(2)
Note that to $`𝒪(dt)`$ the noise has been eliminated in favour of an effective interaction among the agents, and the $`\sigma `$ becomes $`\sigma ^2=(ND)^1\overline{\stackrel{}{a}(t)\stackrel{}{a}(t)}`$.

At $`T=0`$, corresponding to the MG, Eqs. (2) are completely deterministic. To first order in $`dt`$ we have

$$dp_i(t)=[h_i+\underset{j}{}J_{ij}z_j(t)]dt,$$
(3)
where $`h_i_j\stackrel{}{\omega }_j\stackrel{}{\xi }_i/ND`$ and $`J_{ij}\stackrel{}{\xi }_j\stackrel{}{\xi }_i/ND`$, while the volatility reads

$$\sigma ^2=\overline{\mathrm{\Omega }}+2\underset{i}{}\overline{h_iz_i(t)}+\underset{ij}{}\overline{J_{ij}z_i(t)z_j(t)},$$
(4)
and $`\mathrm{\Omega }_{ij}\stackrel{}{\omega }_i\stackrel{}{\omega }_j/ND`$. In order to check the continuous time approximation at $`T=0`$ we have simulated Eqs. (3). Results are presented in Fig. 1. We can see that this approximation reproduces all the features of the original MG. Note that in Eqs. (3) all stochasticity coming from the information has dropped out, the only effect being a small quantitative deviation in the low $`d`$ region.

Eqs. (3) can be rewritten as $`d𝐩/dt=_𝐬`$, where $`𝐩(p_1,\mathrm{},p_N)`$, similarly for $`𝐬`$, and

$$=\frac{1}{2}\mathrm{\Omega }+\underset{i}{}h_is_i+\frac{1}{2}\underset{ij}{}J_{ij}s_is_j.$$
(5)
This is similar to what was done in  for the time and information averages of $`z_i`$. There the value of $`\sigma `$ was related to the average extrema of $``$ by assuming that the system equilibrated. A good agreement with the numerics was found in the phase $`d>d_c`$, but this method failed to reproduce the behaviour in the $`d<d_c`$ phase (see Fig. 1). This disagreement was speculated as due to the need for terms with higher order time derivatives in the continuous time equations. This is clearly wrong, since, as we have just shown, Eqs. (3) describe correctly the dynamics of the model for all values of $`d`$ (see Fig. 1).

The phase $`d<d_c`$ of the MG is very sensitive to the initial conditions. In Fig. 2 we show the results of simulating both the original dynamics Eqs. (1) and the continuous time approximation Eqs. (3) starting from random initial conditions $`|p_i(0)|=𝒪(1)`$. From Fig. 2 see that the behaviour of both $`\sigma `$ and $`\varphi `$ is different from that of Fig. 1 in the region $`d<d_c`$: the system stays in the better-than-random phase for all values of $`d`$. Again the continuous time dynamics is very close to the original discrete one. This sensitivity of the results to the initial conditions is a clear indication that the system does not equilibrate for $`d<d_c`$, and raises a question on the existence and character of the “phase transition” in the MG .

When the temperature is different from zero the TMG Eqs. (2) still depend on the stochastic choice of strategies $`R_i^{}(t)`$, even at leading order. At each time step, $`R_i^{}`$ takes one of the two possible values $`R_i^{1,2}`$, defining a stochastic jump process. In order to write the corresponding Master Equation we need to know the transition probabilities. The r.h.s. of Eq. (2), which we denote $`\mathrm{\Delta }_i`$, is a normalized sum of $`N`$ random numbers $`\stackrel{}{\xi }_iR_j^{}(t)`$, each with mean $`\stackrel{}{\xi }_i\stackrel{}{\omega }_j+\stackrel{}{\xi }_i\stackrel{}{\xi }_js_j(t)`$, and variance $`(\stackrel{}{\xi }_i\stackrel{}{\xi }_j)^2[1s_j^2(t)]`$. By the central limit theorem, we know that for $`N`$ large $`\mathrm{\Delta }_i`$ will tend to be normally distributed with mean $`\mathrm{\Delta }_i=/s_i`$, and variance $`\mathrm{\Delta }_i^2=_jJ_{ij}^2[1s_j^2(t)]`$, where $``$ stands for average over realizations of the random process $`\mu (t)`$. Moreover, $`\mathrm{\Delta }_i`$ and $`\mathrm{\Delta }_{ji}`$ are correlated, the covariance matrix given by

$`M_{ij}[𝐩(t)]`$ $``$ $`\mathrm{\Delta }_i\mathrm{\Delta }_j\mathrm{\Delta }_i\mathrm{\Delta }_j`$ (6)
$`=`$ $`{\displaystyle \underset{k}{}}J_{ik}J_{jk}\left[1s_k^2(t)\right].`$ (7)
Collecting these results, we obtain the transition probabilities in the large $`N`$ limit, $`W(𝐩^{}|𝐩)=\mathrm{\Phi }(_𝐬;𝕄)`$, where $`\mathrm{\Phi }`$ corresponds to the normal distribution with mean $`_𝐬`$ and covariance matrix $`𝕄\{M_{ij}\}`$. Note that $`/s_i𝒪(1)`$, and $`M_{ij}𝒪(1/N)`$, so that fluctuations are also of $`𝒪(1)`$ and thus are not suppressed when $`N\mathrm{}`$.

The $`R_i^{}(t)`$ are chosen independently at each time. If we make the natural assumption that in the limit $`dt0`$ their correlation at different times is a $`\delta `$-function, the Master Equation becomes a Fokker-Planck equation by means of Kramers-Moyal expansion

$$\frac{𝒫}{t}=\underset{i}{}\frac{}{p_i}\left(\frac{}{s_i}𝒫\right)+\frac{1}{2}\underset{ij}{}\frac{^2}{p_ip_j}\left(M_{ij}𝒫\right).$$
(8)
We therefore conclude that the dynamics of the TMG is effectively described by a set of stochastic differential equations for the point differences

$$d𝐩=_𝐬dt+𝕄d𝐖,$$
(9)
where $`𝐖(t)`$ is an $`N`$-dimensional Wiener process, and the volatility is given by $`\sigma ^2=2\overline{}+_i\overline{J_{ii}}_i\overline{J_{ii}s_i^2}`$.

We have checked by means of extensive numerical simulations that Eqs. (9) give the same results as Eqs. (2), up to statistical errors. Figs. 3 and 4 present the results from the continuous time dynamics Eqs. (9). For these simulations we have chosen for the strategy-use probabilities the form $`\pi _i^{1,2}(t)\mathrm{exp}\pm \beta z_i(t)`$ which makes the numerics simpler. Similar results can be obtained with the form $`\pi _i^{1,2}(t)\mathrm{exp}\beta P(\stackrel{}{R}_i^{1,2},t)`$, but a small-$`p`$ cutoff of $`𝒪(dt)`$ is required to avoid the system getting trapped in the $`p=0`$ region. In Fig. 3 we plot the volatility as a function of the temperature, showing that the behaviour is the same as the one found in : for $`d<d_c`$, as $`T`$ is increased $`\sigma `$ first drops to a minimum, and then increases towards the random case $`\sigma =1`$; for $`d>d_c`$, the optimum value is the MG one, and $`\sigma `$ simply grows monotonically to $`1`$. In the inset we give $`\sigma `$ as a function of $`d`$ for different temperatures.

Fig. 4 shows how the fraction of frozen agents varies as the temperature is increased. For all values of $`d`$ there is a clear jump at $`T=𝒪(1)`$ from the MG value to $`\varphi =1`$ . Figs. 3 and 4 determine the phase diagram of the TMG in the $`(d,T)`$ plane. It is schematically depicted in the inset of Fig. 4. For low $`d`$ and $`T`$ the system performs worse than random, while for large enough values of $`T`$ the system becomes random, independently of $`d`$. Better than random performance is achieved between these two regions. It is important to note that, as in the case of the MG, the phase structure of the TMG depends strongly on the initial conditions. In particular, the low $`(d,T)`$ part of the phase diagram shrinks to zero for finite random initial conditions.

In the case where the probabilities are defined as $`\pi _i^{1,2}(t)\mathrm{exp}\beta P(\stackrel{}{R}_i^{1,2},t)`$, the monotonic increase of $`\sigma `$ to the random value at large $`T`$ reported in  is due to finite waiting times, as pointed out in : for $`tNT`$ the volatility stays at the minimum value for any finite $`T1`$. This phenomenon is easily understood from Eqs. (9). For large values of $`T`$, there is first a transient in which $`s_i=\mathrm{tanh}\beta p_i`$ are close to zero, and Eqs. (9) reduce to $`dp_ih_idt+_{jk}J_{ik}J_{jk}dW_j`$, i.e., the point differences of all the agents do a randomly biased Brownian motion, and the system performs as in the random case. Eventually, however, $`p_i`$ become of $`𝒪(T)`$ and $`s_i`$ finite, and the system effectively behaves as for $`T𝒪(1)`$. Note that this cannot happen when $`s_i=\mathrm{tanh}\beta z_i`$.

Eqs. (9) are much simpler than the original ones for the microscopic description of the TMG. The external information has been replaced by an interaction among the agents and the random strategy choice has given rise to the diffusive term. They describe a dynamics which is different from the relaxation of disordered systems usually found in physical problems: the random force $`_𝐬`$ can be written as the gradient of a potential function only up to a factor, which amounts to a metric in the space of $`𝐩`$, and the non-trivial diffusion matrix $`𝕄`$ depends both on the variables $`𝐩`$ and on the quenched disorder of the problem. Finding adequate analytic asymptotic solutions to this dynamics is the next challenging task.

We thank J.-P. Bouchaud, A. Cavagna, I. Giardina, G. Lythe and D. Williams for useful discussions. This work was supported by EC Grant No. ARG/B7-3011/94/27 and EPSRC Grant No. GR/M04426.


text:
# A Combinatorial, Strongly Polynomial-Time Algorithm for Minimizing Submodular Functions

## 1 Introduction

Grötschel, Lovász, and Schrijver  revealed the polynomial-time equivalence between the optimization and separation problems in combinatorial optimization via the ellipsoid method. Since then, many combinatorial problems have been shown to be polynomial-time solvable by means of their framework. The problem of minimizing submodular (set) functions is among these problems. Since the ellipsoid method is far from being efficient in practice and is not combinatorial, efficient combinatorial algorithms for submodular function minimization have been desired for a long time.

A function $`f`$ on all the subsets of a finite set $`V`$ is called submodular if it satisfies

$$f(X)+f(Y)f(XY)+f(XY)X,YV.$$

We suppose that $`f(\mathrm{})=0`$ without loss of generality throughout this paper.

Submodular functions arise in various branches of mathematical engineering such as combinatorial optimization and information theory. There are also close connections between submodularity and convexity . Examples include the matroid rank function, the cut capacity function, and the entropy function. In each of these and other applications, the base polyhedron defined by

$$\mathrm{B}(f)=\{xx𝐑^V,x(V)=f(V),XV:x(X)f(X)\}$$
(1.1)
often plays an important role, where $`x(X)=_{vX}x(v)`$ for any $`XV`$.

Linear optimization problems over base polyhedra are efficiently solvable by the greedy algorithm of Edmonds . Thus Grötschel, Lovász, and Schrijver  assert that the submodular function minimization, which is equivalent to the separation problem, is solvable in polynomial time by the ellipsoid method. Later, they also devise a strongly polynomial-time algorithm within their framework using the ellipsoid method .

A first step towards a combinatorial strongly polynomial-time algorithm was taken by Cunningham , who devised a strongly polynomial-time algorithm for testing membership in matroid polyhedra as well as a pseudopolynomial-time algorithm for minimizing submodular functions. Recently, Narayanan  improved the running time bounds of these combinatorial algorithms by a rounding technique. Based on the minimum-norm base characterization of minimizers due to Fujishige , Sohoni  gave another combinatorial pseudopolynomial-time algorithm for submodular function minimization.

For the problem of minimizing a symmetric submodular function over proper nonempty subsets, Queyranne  presented a combinatorial strongly polynomial-time algorithm, extending the undirected minimum cut algorithm of Nagamochi and Ibaraki .

In this paper, we present a combinatorial polynomial-time algorithm for submodular function minimization. Our algorithm uses an augmenting path approach with reference to a convex combination of extreme points of the base polyhedron. Such an approach was first introduced by Cunningham for minimizing submodular functions that arise from the separation problem for matroid polyhedra . This was adapted for general submodular function minimization by Bixby, Cunningham, and Topkis  and improved by Cunningham  to obtain a pseudopolynomial-time algorithm.

A fundamental tool in these algorithms is to move from one extreme point of the base polyhedron to an adjacent extreme point via an exchange operation that increases one coordinate and decreases another coordinate by the same quantity. This quantity is called the exchange capacity. These previous methods maintain a directed graph on the underlying set that represents the possible exchange operations. They are inefficient since the lower bound on the size of each augmentation is too small. In traditional network flow problems, it is possible to surmount this difficulty by augmenting only on paths of sufficiently large capacity . However, it has proved difficult to adapt this scaling approach to work in the setting of submodular function minimization, mainly because the amount of augmentation is determined by exchange capacities multiplied by the convex combination coefficients. These coefficients can be as small as the reciprocal of the maximum absolute value of the submodular function.

To overcome this difficulty, we augment the directed graph corresponding to allowable exchanges with the complete directed graph on the underlying set, letting the capacity of this additional arc set depend directly on our scaling parameter. This technique was first introduced by Iwata , who used it to develop the first polynomial-time capacity-scaling algorithm for the submodular flow problem of Edmonds and Giles . This algorithm was later refined by Fleischer, Iwata, and McCormick  into one of the fastest algorithms for submodular flow. Our work in this paper builds on ideas in this latter paper to develop a capacity-scaling, augmenting-path algorithm for submodular function minimization. The running time of the resulting algorithm is weakly polynomial, i.e., bounded by a polynomial in the size of the underlying set and the largest length of the function value. Even under the similarity assumption that the largest function value is bounded by a polynomial in the size of the underlying set, our algorithm is faster than the best previous combinatorial, pseudopolynomial-time algorithm .

We then modify our scaling algorithm to run in strongly polynomial time, i.e., in time bounded by a polynomial in the size of the underlying set, independently of the largest length of the function value. To make a weakly polynomial-time algorithm run in strongly polynomial time, Frank and Tardos  developed a generic preprocessing technique that is applicable to a fairly wide class of combinatorial optimization problems including the submodular flow problem and testing membership in matroid polyhedra. However, this framework does not apply to submodular function minimization. Instead, we devise a combinatorial algorithm that repeatedly detects an element that belongs to every minimizer or an ordered pair of elements with the property that if the first belongs to a minimizer then the second does.

There are some practical problems, in dynamic flows , facility location , and multi-terminal source coding , where the polynomial-time solvability relies on a submodular function minimization routine. Goemans and Ramakrishnan  discussed a class of submodular function minimization problems over restricted families of subsets. Their solution is combinatorial modulo an oracle for submodular function minimization on distributive lattices. Our algorithm can be used to provide combinatorial, strongly polynomial-time algorithms for these problems.

This paper is organized as follows. Section 2 provides preliminaries on submodular functions. Section 3 presents a scaling algorithm for submodular function minimization, which runs in weakly polynomial time. Section 4 is devoted to the strongly polynomial-time algorithm. Finally, we discuss extensions in Section 5.

## 2 Preliminaries

We denote by $`𝐙`$ and $`𝐑`$ the set of integers and the set of reals, respectively. Let $`V`$ be a finite nonempty set of cardinality $`|V|=n`$. For a vector $`x𝐑^V`$ we define a modular function $`x:2^V𝐑`$ by $`x(X)=_{vV}x(v)`$. For each $`uV`$, we denote by $`\chi _u`$ the unit vector in $`𝐑^V`$ such that $`\chi _u(v)=1`$ if $`v=u`$ and $`=0`$ otherwise.

Given a submodular function $`f`$ with $`f(\mathrm{})=0`$ and its associated base polyhedron $`\mathrm{B}(f)`$ as defined in (1.1), we call a vector $`x\mathrm{B}(f)`$ a base. An extreme point of $`\mathrm{B}(f)`$ is called an extreme base. A fundamental step in submodular function minimization algorithms is to move from one base $`x`$ to another base $`x^{}`$ via an exchange operation that increases one coordinate while decreasing another coordinate by the same amount. The maximum amount of increase that ensures $`x^{}\mathrm{B}(f)`$ is called the exchange capacity. More precisely, for any base $`x\mathrm{B}(f)`$ and any distinct $`u,vV`$ the exchange capacity is

$$\stackrel{~}{\mathrm{c}}(x,u,v)=\mathrm{max}\{\alpha \alpha 𝐑,x+\alpha (\chi _u\chi _v)\mathrm{B}(f)\}.$$
(2.1)
The exchange capacity can also be expressed as

$$\stackrel{~}{\mathrm{c}}(x,u,v)=\mathrm{min}\{f(X)x(X)uXV\backslash \{v\}\}.$$
(2.2)
In general, computing $`\stackrel{~}{\mathrm{c}}(x,u,v)`$ is as hard as submodular function minimization, even when $`x`$ is an extreme base. However, if $`x`$ is an extreme base, then for special pairs of vertices $`u`$ and $`v`$, the exchange capacity $`\stackrel{~}{\mathrm{c}}(x,u,v)`$ can be computed with one function evaluation as follows.

Let $`L=(v_1,v_2,\mathrm{},v_n)`$ be a linear ordering of $`V`$. For any $`k\{1,2,\mathrm{},n\}`$, we define $`L(v_k)=\{v_1,v_2,\mathrm{},v_k\}`$. Given such a linear ordering, the greedy algorithm of Edmonds  computes

$$y(v_i)=f(L(v_i))f(L(v_{i1}))(i=1,2,\mathrm{},n),$$
(2.3)
where $`L(v_0)=\mathrm{}`$. The resulting vector $`y𝐑^V`$ is an extreme base $`y\mathrm{B}(f)`$. Conversely, any extreme base can be generated by applying the greedy algorithm to an appropriate linear ordering. Note that a linear ordering $`L=(v_1,v_2,\mathrm{},v_n)`$ generates an extreme base $`y`$ if and only if $`y(L(v_i))=f(L(v_i))`$ for $`i=1,2,\mathrm{},n`$. For any base $`y\mathrm{B}(f)`$, a set $`XV`$ is called $`y`$-tight if $`y(X)=f(X)`$. A pair $`(u,v)`$ is called eligible for $`y`$ if $`u`$ immediately succeeds $`v`$ in some linear ordering that generates $`y`$. The following lemma enables us to compute an exchange capacity $`\stackrel{~}{\mathrm{c}}(y,u,v)`$ if $`(u,v)`$ is eligible for $`y`$.

###### Lemma 2.1

Let $`L`$ be a linear ordering of $`V`$ that generates an extreme base $`y\mathrm{B}(f)`$. Let $`L^{}`$ be the linear ordering obtained by interchanging $`u`$ and $`v`$ that are consecutive in $`L`$. Then the extreme base $`y^{}`$ generated by $`L^{}`$ satisfies

$$y^{}=y+\beta (\chi _u\chi _v)$$
(2.4)
with

$$\beta =f(L(u)\backslash \{v\})f(L(u))+y(v).$$
(2.5)
Moreover, we have $`\stackrel{~}{\mathrm{c}}(y,u,v)=\beta `$.

Proof. Equations (2.4) and (2.5) follow from the greedy algorithm (see (2.3)). By the definition (2.1) of the exchange capacity, we have $`\beta \stackrel{~}{\mathrm{c}}(y,u,v)`$. Since $`y(L(u))=f(L(u))`$, it follows from (2.2) and (2.5) that $`\beta \stackrel{~}{\mathrm{c}}(y,u,v)`$. Thus we obtain $`\beta =\stackrel{~}{\mathrm{c}}(y,u,v)`$.

We will use Lemma 2.1 to transform one extreme base into another and to update the corresponding linear ordering.

For any vector $`x𝐑^V`$, we denote by $`x^{}`$ the vector in $`𝐑^V`$ defined by $`x^{}(v)=\mathrm{min}\{0,x(v)\}`$ for $`vV`$. The following fundamental lemma easily follows from a theorem of Edmonds  on the vector reduction of polymatroids (see \[12, Corollaries 3.4 and 3.5\]).

###### Lemma 2.2

For a submodular function $`f:2^V𝐑`$ we have

$$\mathrm{max}\{x^{}(V)x\mathrm{B}(f)\}=\mathrm{min}\{f(X)XV\}.$$

If $`f`$ is integer-valued, then the maximizer $`x`$ can be chosen from among integral bases.

We will not use the integrality property indicated in the latter half of this lemma. Lemma 2.2 shows a min-max relation of strong duality. A weak duality is described as follows: For any base $`x\mathrm{B}(f)`$ and any $`XV`$ we have $`x^{}(V)f(X)`$. We call the difference $`f(X)x^{}(V)`$ a duality gap. Note that, if $`f`$ is integer-valued and the duality gap $`f(X)x^{}(V)`$ is less than one for some $`x\mathrm{B}(f)`$ and $`XV`$, then $`X`$ minimizes $`f`$.

## 3 A Scaling Algorithm

In this section, we describe a combinatorial algorithm for minimizing an integer-valued submodular function $`f:2^V𝐙`$ with $`f(\mathrm{})=0`$. We assume an evaluation oracle for the function value of $`f`$. Let $`M`$ denote an upper bound on $`|f(X)|`$ among $`XV`$. Note that we can easily compute $`M`$ by $`\mathrm{O}(n)`$ calls for the evaluation oracle as follows. Let $`y`$ be an extreme base generated by a linear ordering $`L`$. For any $`XV`$, we have $`y^{}(V)y(X)f(X){\displaystyle \underset{vV}{}}\mathrm{max}\{0,f(\{v\})\}`$. Thus we obtain $`M=\mathrm{max}\{|y^{}(V)|,{\displaystyle \underset{vV}{}}\mathrm{max}\{0,f(\{v\})\}\}`$.

### 3.1 Algorithm Outline

As indicated earlier, our algorithm uses an augmenting path approach to submodular function minimization . As with previous algorithms, we maintain a base $`x\mathrm{B}(f)`$ as a convex combination of extreme bases $`y_i\mathrm{B}(f)`$ indexed by $`iI`$, so that $`x=_{iI}\lambda _iy_i`$. Roughly speaking, these previous algorithms use a directed graph with the arc set defined by the pairs of vertices that are eligible for some $`y_i`$, $`iI`$. They seek to increase $`x^{}(V)`$ by performing exchange operations along a path of arcs from vertices $`s`$ with $`x(s)<0`$ to vertices $`t`$ with $`x(t)>0`$. The algorithms stop with an optimal $`x`$ when there are no more augmenting paths. The corresponding minimizer $`X`$ is determined by the set of vertices reachable from vertices $`s`$ with $`x(s)<0`$.

To adapt this procedure to a scaling framework, we use a complete directed graph on $`V`$ with arc capacities that depend directly on our scaling parameter $`\delta `$, an idea first introduced for submodular flows in . Let $`\phi :V\times V𝐑`$ be skew-symmetric, i.e., $`\phi (u,v)+\phi (v,u)=0`$ for $`u,vV`$, and $`\delta `$-feasible in that it satisfies capacity constraints $`\delta \phi (u,v)\delta `$ for every $`u,vV`$. The function $`\phi `$ can be regarded as a flow in the complete directed graph $`G=(V,E)`$ with the vertex set $`V`$ and the arc set $`E=V\times V`$. The boundary $`\phi :V𝐑`$ of $`\phi `$ is defined by

$$\phi (v)=\underset{uV}{}\phi (u,v)(vV).$$
(3.1)
Instead of trying to maximize $`x^{}(V)`$ directly, we define $`z=x\phi `$. Our algorithm seeks to maximize $`z^{}(V)`$ and thereby increases $`x^{}(V)`$.

We also maintain linear orderings $`L_i`$ for $`iI`$ and extreme bases $`y_i`$ generated by them. We start with an arbitrary linear ordering $`L`$ on $`V`$ and the extreme base $`x\mathrm{B}(f)`$ generated by $`L`$. In addition, we start with the zero flow $`\phi =\mathrm{𝟎}`$. Thus, initially $`z^{}(V)=x^{}(V)nM`$. We seek to increase $`z^{}(V)`$, and in doing so, obtain improvements in $`x^{}(V)`$, via the $`\delta `$-feasibilty of $`\phi `$.

The algorithm consists of scaling phases with a positive parameter $`\delta `$. It starts with $`\delta =M`$, cuts $`\delta `$ in half at the beginning of each scaling phase, and ends with $`\delta <1/n^2`$. Each $`\delta `$-scaling phase maintains a $`\delta `$-feasible flow $`\phi `$, and uses the residual graph $`G(\phi )=(V,E(\phi ))`$ with the arc set

$$E(\phi )=\{(u,v)|u,vV,uv,\phi (u,v)0\}.$$
(3.2)
Intuitively, $`E(\phi )`$ consists of the arcs through which we can augment the flow $`\phi `$ by $`\delta `$ without violating the capacity constraints.

A $`\delta `$-scaling phase starts by preprocessing $`\phi `$ to make it $`\delta `$-feasible, and then repeatedly searches to send flow along augmenting paths in $`G(\phi )`$ from $`S:=\{vvV,z(v)\delta \}`$ to $`T:=\{vvV,z(v)\delta \}`$. Such a directed path is called a $`\delta `$-augmenting path.

If there are no $`\delta `$-augmenting paths, then the algorithm checks whether there is a pair $`(u,v)`$ of vertices such that $`u`$ is reachable from $`S`$ by a path with residual capacity $`\delta `$, $`v`$ is not, and $`u`$ immediately succeeds $`v`$ in a linear ordering that generates $`y_i`$ for some $`iI`$. We perform the appropriate exchange operation, and modify $`\phi `$ by creating residual capacity on $`(u,v)`$ so that $`z=x\phi `$ is invariant. This operation may increase the set of vertices reachable from $`S`$ on paths of residual capacity $`\delta `$. Once a $`\delta `$-augmenting path is found, the algorithm augments the flow $`\phi `$ by $`\delta `$ through the path without changing $`x`$. As a consequence, $`z^{}(V)`$ increases by $`\delta `$ in one iteration. This is an extension of a technique for handling exchange capacity arcs in submodular flows first developed in .

### 3.2 Algorithm Details

We now describe the scaling algorithm more precisely. Figure 1 provides a formal description.

At the beginning of the $`\delta `$-scaling phase, after $`\delta `$ is cut in half, the current flow $`\phi `$ is $`2\delta `$-feasible. Then the algorithm modifies each $`\phi (u,v)`$ to the nearest value within the interval $`[\delta ,\delta ]`$ to make it $`\delta `$-feasible. This may decrease $`z^{}(V)`$ for $`z=x\phi `$ by at most $`\left(\genfrac{}{}{0pt}{}{n}{2}\right)\delta `$. The rest of the $`\delta `$-scaling phase aims at increasing $`z^{}(V)`$ by augmenting flow along $`\delta `$-augmenting paths.

Let $`W`$ denote the set of vertices reachable by directed paths from $`S`$ in $`G(\phi )`$. For each $`iI`$ we keep a linear ordering $`L_i`$ that generates $`y_i`$. We call a vertex $`vV\backslash W`$ active if $`v`$ is the last vertex in $`L_i`$ among vertices in $`V\backslash W`$ that satisfies $`W\backslash L_i(v)\mathrm{}`$. If $`v`$ is active in $`L_i`$, we call $`(i,v)`$ an active pair. We denote by $`Z`$ the set of the current active pairs.

If $`WT=\mathrm{}`$, there is no $`\delta `$-augmenting path in $`G(\phi )`$. Then, as long as there is an active pair $`(i,v)`$, i.e., $`Z\mathrm{}`$, the algorithm repeatedly picks an active pair $`(i,v)Z`$ and applies $`\mathrm{𝖯𝗎𝗌𝗁}(i,u,v)`$ to $`u`$ that succeeds $`v`$ in $`L_i`$. Note that $`v`$ active implies that $`uW`$.

The operation $`\mathrm{𝖯𝗎𝗌𝗁}(i,u,v)`$ starts with reducing the flow through $`(u,v)`$ by $`\alpha =\mathrm{min}\{\delta ,\lambda _i\stackrel{~}{\mathrm{c}}(y_i,u,v)\}`$. The boundary $`\phi `$ moves to $`\phi +\alpha (\chi _u\chi _v)`$. The operation $`\mathrm{𝖯𝗎𝗌𝗁}(i,u,v)`$ is called saturating if $`\alpha =\lambda _i\stackrel{~}{\mathrm{c}}(y_i,u,v)`$. Otherwise, it is called nonsaturating. A nonsaturating $`\mathrm{𝖯𝗎𝗌𝗁}(i,u,v)`$ adds to $`I`$ a new index $`k`$ with $`y_k:=y_i`$, $`\lambda _k:=\lambda _i\alpha /\stackrel{~}{\mathrm{c}}(y_i,u,v)`$, and $`L_k:=L_i`$. Whether the $`\mathrm{𝖯𝗎𝗌𝗁}(i,u,v)`$ is saturating or not, it updates $`y_i`$ as $`y_i:=y_i+\stackrel{~}{\mathrm{c}}(y_i,u,v)(\chi _u\chi _v)`$, $`\lambda _i:=\alpha /\stackrel{~}{\mathrm{c}}(y_i,u,v)`$ if $`\stackrel{~}{\mathrm{c}}(y_i,u,v)>0`$, and $`L_i`$ by interchanging $`u`$ and $`v`$. Then the current base $`x`$ moves to $`x+\alpha (\chi _u\chi _v)`$. Thus $`z=x\phi `$ is invariant.

Each time the algorithm applies the push operation, it updates the set $`W`$ of vertices reachable from $`S`$ in $`G(\phi )`$ and the set $`Z`$ of active pairs. If $`\mathrm{𝖯𝗎𝗌𝗁}(i,u,v)`$ is nonsaturating, it makes $`v`$ reachable from $`S`$ in $`G(\phi )`$, and hence $`W`$ is enlarged. Once $`v`$ becomes reachable from $`S`$ in $`G(\phi )`$, it will never become active again for any $`iI`$ until the algorithm finds a $`\delta `$-augmenting path or all the active pairs disappear. Note that we encounter at most $`n`$ nonsaturating pushes before we find a $`\delta `$-augmenting path or all the active pairs disappear. Each time the algorithm picks an active pair $`(i,v)`$ and applies $`\mathrm{𝖯𝗎𝗌𝗁}(i,u,v)`$, the vertex $`v`$ shifts towards the end of $`L_i`$. Hence the algorithm picks an active pair $`(i,v)`$ at most $`n`$ times before $`v`$ enters $`W`$ or $`W\backslash L_i(v)=\mathrm{}`$. At this point $`v`$ becomes inactive, and remains inactive until the next augmentation. Hence, for each $`iI`$ the total time required for processing active vertices in $`L_i`$ is $`\mathrm{O}(n^2)`$.

We note that we could relax the definition of an active vertex to include any vertex $`vV\backslash W`$ whose immediate successor in $`L_i`$ belongs to $`W`$. The correctness argument would apply without modifications. However, care is needed to obtain an efficient implementation.

If we find a $`\delta `$-augmenting path, the algorithm augments $`\delta `$ units of flow along the path, which effectively increases $`z^{}(V)`$ by $`\delta `$. We also compute an expression for $`x`$ as a convex combination of at most $`n`$ affinely independent extreme bases $`y_i`$, chosen from the current $`y_i`$’s. This computation is a standard linear programming technique of transforming feasible solutions into basic feasible solutions. If the set of extreme points are not affinely independent, there is a set of coefficients $`\mu _i`$ for $`iI`$ that is not identically zero and satisfies $`\mu _iy_i=0`$ and $`\mu _i=0`$. Using Gaussian elimination, we can start computing such $`\mu _i`$ until a dependency is detected. At this point, we eliminate the dependency by computing $`\theta :=\mathrm{min}\{\lambda _i/\mu _i|\mu _i>0\}`$ and update $`\lambda _i:=\lambda _i\theta \mu _i`$ for $`iI`$. At least one $`iI`$ satisfies $`\lambda _i=0`$. Delete such $`i`$ from $`I`$. We continue this procedure until we eventually obtain affine independence. Since a new index $`k`$ is added to $`I`$ only as a result of a nonsaturating push, $`|I|2n`$ after finding an augmenting path. The bottleneck in this procedure is the time spent computing the coefficients $`\mu _i`$, which is $`O(n^3)`$ overall.

A $`\delta `$-scaling phase ends when either $`S=\mathrm{}`$, $`T=\mathrm{}`$, or $`Z=\mathrm{}`$. In the last case, we have a set of vertices $`WV`$ that are reachable from $`S`$ in $`G(\phi )`$ such that $`WT=\mathrm{}`$.

###### Lemma 3.1

If $`Z=\mathrm{}`$, then $`W`$ is tight for $`x`$.

Proof. If $`Z=\mathrm{}`$, for each $`iI`$ the first $`|W|`$ vertices in $`L_i`$ must belong to $`W`$. Then it follows from (2.3) that $`y_i(W)=f(W)`$. Since $`x=_{iI}\lambda _iy_i`$ and $`_{iI}\lambda _i=1`$, this implies $`x(W)=_{iI}\lambda _iy_i(W)=f(W)`$.

### 3.3 Correctness and Complexity

We now investigate the number of iterations in each $`\delta `$-scaling phase. To do this, we prove relaxed weak and strong dualities. The next lemma shows a relaxed weak duality.

###### Lemma 3.2

For any base $`x\mathrm{B}(f)`$ and any $`\delta `$-feasible flow $`\phi `$, the vector $`z=x\phi `$ satisfies $`z^{}(V)f(X)+\left(\genfrac{}{}{0pt}{}{n}{2}\right)\delta `$ for any $`XV`$.

Proof. For any $`XV`$ we have $`x(X)f(X)`$ and $`\phi (X)\left(\genfrac{}{}{0pt}{}{n}{2}\right)\delta `$, and hence $`z^{}(V)z(X)f(X)+\left(\genfrac{}{}{0pt}{}{n}{2}\right)\delta `$.

A relaxed strong duality is given as follows.

###### Lemma 3.3

At the end of each $`\delta `$-scaling phase, the following (i)(iii) hold for $`x`$ and $`z=x\phi `$.

If $`S=\mathrm{}`$, then $`x^{}(V)f(\mathrm{})n^2\delta `$ and $`z^{}(V)f(\mathrm{})n\delta `$.

If $`T=\mathrm{}`$, then $`x^{}(V)f(V)n^2\delta `$ and $`z^{}(V)f(V)n\delta `$.

If $`W`$ is tight for $`x`$, then $`x^{}(V)f(W)n^2\delta `$ and $`z^{}(V)f(W)n\delta `$.

Proof. When the $`\delta `$-scaling phase finishes with $`S=\mathrm{}`$, we have $`x(v)>\phi (v)\delta n\delta `$ for every $`vV`$, which implies $`x^{}(V)f(\mathrm{})n^2\delta `$ as well as $`z^{}(V)f(\mathrm{})n\delta `$. Similarly, when the $`\delta `$-scaling phase finishes with $`T=\mathrm{}`$, we have $`x(v)<\phi (v)+\delta n\delta `$ for every $`vV`$, which implies $`x^{}(V)x(V)n^2\delta =f(V)n^2\delta `$ as well as $`z^{}(V)x(V)n\delta `$.

When the $`\delta `$-scaling phase ends with $`x(W)=f(W)`$ due to Lemma 3.1, then $`SWV\backslash T`$ and $`\phi (W)<0`$. By the definitions of $`S`$ and $`T`$, we also have $`x(v)>\phi (v)\delta n\delta `$ for every $`vV\backslash W`$ and $`x(v)<\phi (v)+\delta n\delta `$ for every $`vW`$. Therefore we have $`x^{}(V)=x^{}(W)+x^{}(V\backslash W)x(W)n\delta |W|n\delta |V\backslash W|=f(W)n^2\delta `$ as well as $`z^{}(V)=z^{}(W)+z^{}(V\backslash W)x(W)\phi (W)|W|\delta \delta |V\backslash W|f(W)n\delta `$.

Lemma 3.3 implies that at the beginning of the $`\delta `$-scaling phase, after $`\delta `$ is cut in half, $`z^{}(V)`$ is at least $`f(X)2n\delta `$ for some $`XV`$. Making the current flow $`\delta `$-feasible decreases $`z^{}(V)`$ by at most $`\left(\genfrac{}{}{0pt}{}{n}{2}\right)\delta `$. Each $`\delta `$-augmentation increases $`z^{}(V)`$ by $`\delta `$. Since $`z^{}(V)`$ is at most $`f(X)+\left(\genfrac{}{}{0pt}{}{n}{2}\right)\delta `$ at the end of a $`\delta `$-phase by Lemma 3.2 the number of $`\delta `$-augmentations per phase is at most $`n^2+n`$ for all phases after the first. Since $`z^{}(V)=x^{}(V)nM`$ at the start of the algorithm, setting the initial $`\delta =M`$ is more than sufficient obtain a similar bound on the number of augmentations in the first phase.

As an immediate consequence of Lemmas 2.2 and 3.3, we also obtain the following.

###### Theorem 3.4

The algorithm obtains a minimizer of $`f`$ at the end of the $`\delta `$-scaling phase with $`\delta <1/n^2`$.

Proof. By Lemma 3.3, the output $`X`$ of the algorithm satisfies $`x^{}(V)f(X)n^2\delta >f(X)1`$. For any $`YV`$, the weak duality in Lemma 2.2 asserts $`x^{}(V)f(Y)`$. Thus we have $`f(X)1<f(Y)`$, which implies by the integrality of $`f`$ that $`X`$ minimizes $`f`$.

###### Theorem 3.5

Algorithm $`\mathrm{𝖲𝖥𝖬}`$ runs in $`\mathrm{O}(n^5\mathrm{log}(nM))`$ time.

Proof. The algorithm starts with $`\delta =M`$ and ends with $`\delta <1/n^2`$, so the algorithm consists of $`\mathrm{O}(\mathrm{log}(nM))`$ scaling phases. Each scaling phase finds $`\mathrm{O}(n^2)`$ $`\delta `$-augmenting paths. To find an augmenting path, we perform at most $`O(n^2)`$ pushes per extreme base. A saturating push requires $`\mathrm{O}(1)`$ time while a nonsaturating one $`\mathrm{O}(n)`$ time. Here, note that there are less than $`n`$ nonsaturating pushes per augmenting path. Hence, the time spent in pushes per augmenting path is $`\mathrm{O}(n^3)`$. After each augmentation, we also update the expression $`x=_{iI}\lambda _iy_i`$, which also takes $`\mathrm{O}(n^3)`$ time per augmentation. Thus the overall complexity of $`\mathrm{𝖲𝖥𝖬}`$ is $`\mathrm{O}(n^5\mathrm{log}(nM))`$.

In this section, we have shown a weakly polynomial-time algorithm for minimizing integer-valued submodular functions. The integrality of a submodular function $`f`$ guarantees that if we have a base $`x\mathrm{B}(f)`$ and a subset $`X`$ of $`V`$ such that the duality gap $`f(X)x^{}(V)`$ is less than one, $`X`$ is a minimizer of $`f`$. Except for this we have not used the integrality of $`f`$. It follows that for any real-valued submodular function $`f:2^V𝐑`$, if we are given a positive lower bound $`ϵ`$ for the difference between the second minimum and the minimum value of $`f`$, the present algorithm works for the submodular function $`(1/ϵ)f`$ and runs in $`\mathrm{O}(n^5\mathrm{log}(nM/ϵ))`$ time, where $`M`$ is an upper bound on $`|f(X)|`$ among $`XV`$.

## 4 A Strongly Polynomial-Time Algorithm

This section presents a strongly polynomial-time algorithm for minimizing submodular functions using the scaling algorithm in Section 3. The new algorithm exploits the following proximity lemma.

###### Lemma 4.1

At the end of the $`\delta `$-scaling phase, if $`x(w)<n^2\delta `$, then $`w`$ belongs to every minimizer of $`f`$.

Proof. Let $`X`$ be any minimizer of $`f`$. There exists a vector $`y\mathrm{B}(f)`$ with $`x^{}y^{}`$ such that $`y^{}(V)=f(X)`$. Note that $`y(v)0`$ for each $`vV\backslash X`$. By Lemma 3.3, there exists a subset $`YV`$ such that $`x^{}(V)f(Y)n^2\delta `$. Then we have $`y^{}(w)x^{}(w)y^{}(V)x^{}(V)f(X)f(Y)+n^2\delta n^2\delta `$. This implies $`y(w)<0`$ due to the assumption, and hence $`wX`$.

Let $`f:2^V𝐑`$ be a submodular function and $`x\mathrm{B}(f)`$ an extreme base whose components are bounded from above by $`\eta >0`$. Assume that there exists a subset $`YV`$ such that $`f(Y)\kappa `$ for some positive parameter $`\kappa `$, which will be specified later as $`\eta /2`$. We then apply the scaling algorithm starting with $`\delta =\eta `$ and the extreme base $`x\mathrm{B}(f)`$. After $`\mathrm{log}_2(n^3\eta /\kappa )`$ scaling phases, $`\delta `$ becomes less than $`\kappa /n^3`$. Since $`x(Y)f(Y)\kappa `$, at least one element $`wY`$ satisfies $`x(w)<n^2\delta `$. By Lemma 4.1, such an element $`w`$ belongs to every minimizer of $`f`$. We denote this procedure by $`\mathrm{𝖥𝗂𝗑}(f,x,\eta )`$.

We now discuss how to apply this procedure to design a strongly polynomial-time algorithm for minimizing a submodular function $`f`$. If $`f(V)>0`$, we replace the value $`f(V)`$ by zero. The set of minimizers remains the same unless the minimum value is zero, in which case we may assert that $`\mathrm{}`$ minimizes $`f`$.

An ordered pair $`(u,v)`$ of distinct vertices $`u,vV`$ is said to be compatible with $`f`$ if $`uX`$ implies $`vX`$ for every minimizer $`X`$ of $`f`$. Our algorithm keeps a directed acyclic graph $`D=(V,F)`$ whose arcs are compatible with $`f`$. Initially, the arc set $`F`$ is empty. Each time the algorithm finds a compatible pair $`(u,v)`$ with $`f`$, it adds $`(u,v)`$ to $`F`$. When this gives rise to a cycle in $`D`$, the algorithm contracts the strongly connected component $`UV`$ to a single vertex and modifies the submodular function $`f`$ by regarding $`U`$ as a singleton.

For each $`vV`$, let $`R(v)`$ denote the set of vertices reachable from $`v`$ in $`D`$ and $`f_v`$ the submodular function on the subsets of $`V\backslash R(v)`$ defined by

$$f_v(X)=f(XR(v))f(R(v))(XV\backslash R(v)).$$

A linear ordering $`(v_1,\mathrm{},v_n)`$ of $`V`$ is called consistent with $`D`$ if $`i<j`$ implies $`(v_i,v_j)F`$. Consider an extreme base $`x\mathrm{B}(f)`$ generated by a linear ordering $`(v_1,v_2,\mathrm{},v_n)`$ consistent with $`D`$. The extreme base generated by a consistent linear ordering is also called consistent. For any consistent extreme base $`x\mathrm{B}(f)`$, the greedy algorithm defines $`x(v)`$ by $`x(v)=f(U)f(U\backslash \{v\})`$ for some $`UV`$ with $`R(v)U`$. It then follows from the submodularity of $`f`$ that $`x`$ satisfies $`x(v)f(R(v))f(R(v)\backslash \{v\})`$ for each $`vV`$.

In each iteration, the algorithm computes

$$\eta =\mathrm{max}\{f(R(v))f(R(v)\backslash \{v\})vV\}.$$
(4.1)
If $`\eta 0`$, then an extreme base $`x\mathrm{B}(f)`$ consistent with $`D`$ satisfies $`x(v)0`$ for each $`vV`$. In this case $`x^{}(V)=x(V)=f(V)`$, which implies that $`V`$ minimizes $`f`$ by Lemma 2.2. If in addition $`f(V)=0`$, then the original function may have had a positive value of $`f(V)`$. Therefore, the algorithm returns $`\mathrm{}`$ or $`V`$ as a minimizer, according to whether $`f(V)=0`$ or $`f(V)<0`$.

If $`\eta >0`$, let $`u`$ be an element that attains the maximum in the right-hand side of (4.1). Then we have $`f(R(u))=f(R(u)\backslash \{u\})+\eta `$, which implies either $`f(R(u))\eta /2>0`$ or $`f(R(u)\backslash \{u\})<\eta /2<0`$ holds.

In the former case ($`f(R(u))\eta /2`$), we have $`f_u(V\backslash R(u))=f(V)f(R(u))\eta /2`$. The algorithm finds a consistent extreme base $`x\mathrm{B}(f_u)`$ generated by a linear ordering $`(v_1,\mathrm{},v_k)`$ consistent with $`D`$, where $`k=|V\backslash R(u)|`$. That is, let $`x(v_1)=f_u(\{v_1\})`$ and $`x(v_j)=f_u(\{v_1,v_2,\mathrm{},v_j\})f_u(\{v_1,v_2,\mathrm{},v_{j1}\})`$ for $`j=2,\mathrm{},k`$. Then the extreme base $`x\mathrm{B}(f_u)`$ satisfies $`x(v)f(R(v))f(R(v)\backslash \{v\})\eta `$. Thus we may apply the procedure $`\mathrm{𝖥𝗂𝗑}(f_u,x,\eta )`$ to find an element $`wV\backslash R(u)`$ that belongs to every minimizer of $`f_u`$. Since $`\kappa =\eta /2`$, the procedure terminates within $`\mathrm{O}(\mathrm{log}n)`$ scaling phases. Consequently, we obtain a new pair $`(u,w)`$ that is compatible with $`f`$. Hence the algorithm adds the arc $`(u,w)`$ to $`F`$.

In the latter case ($`f(R(u)\backslash \{u\})<\eta /2`$), we compute an extreme base $`x\mathrm{B}(f)`$ consistent with $`D`$ by the greedy algorithm, and then apply the procedure $`\mathrm{𝖥𝗂𝗑}(f,x,\eta )`$ to find an element $`wR(u)`$ that belongs to every minimizer of $`f`$. Since $`x(v)\eta `$ for every $`vV`$ and $`\kappa =\eta /2`$, the procedure terminates within $`\mathrm{O}(\mathrm{log}n)`$ scaling phases. Note that every minimizer of $`f`$ includes $`R(w)`$. Thus it suffices to minimize the submodular function $`f_w`$, which is now defined on a smaller underlying set. Figure 3 provides a formal description of the strongly polynomial-time algorithm.

###### Theorem 4.2

The algorithm in Figure 3 computes the minimizer of a submodular function in $`\mathrm{O}(n^7\mathrm{log}n)`$ time, which is strongly polynomial.

Proof. Each time we call the procedure $`\mathrm{𝖥𝗂𝗑}`$, the algorithm adds a new arc to $`D`$ or deletes a set of vertices. This can happen at most $`n^2`$ times. Thus the overall running time of the algorithm is $`\mathrm{O}(n^7\mathrm{log}n)`$, which is strongly polynomial.

## 5 Concluding Remarks

This paper presents a strongly polynomial-time algorithm for minimizing submodular functions defined on Boolean lattices. We now briefly discuss minimizing submodular functions defined on more general lattices.

Consider a submodular function $`f:𝒟𝐑`$ defined on a distributive lattice $`𝒟`$ represented by a poset $`𝒫`$ on $`V`$. Then the associated base polyhedron is unbounded in general (see ).

An easy way to minimize such a function $`f`$ is to consider the reduction of $`f`$ by a sufficiently large vector. As described in \[12, p. 56\], we can compute an upper bound $`\widehat{M}`$ on $`|f(X)|`$ $`(X𝒟)`$. Let $`f^{}`$ be the rank function of the reduction by a vector with each component being equal to $`\widehat{M}`$. The submodular function $`f^{}`$ is defined on $`2^V`$ and the set of minimizers of $`f^{}`$ coincides with that of $`f`$. Thus, we may apply our algorithms. However, each evaluation of the function value of $`f^{}`$ requires $`\mathrm{O}(n^2)`$ elementary operations in addition to a single call for the evaluation of $`f`$. Consequently, this approach takes $`\mathrm{O}(n^7\mathrm{min}\{\mathrm{log}(n\widehat{M}),n^2\mathrm{log}n\})`$ time.

Alternatively, we can slightly extend the algorithms in Sections 3 and 4 by keeping the base $`x\mathrm{B}(f)`$ as a convex combination of extreme bases $`y_i`$’s plus a vector in the characteristic cone of $`\mathrm{B}(f)`$. The latter can be represented as a boundary of a nonnegative flow in the Hasse diagram of $`𝒫`$. This extension enables us to minimize $`f`$ in $`\mathrm{O}(n^5\mathrm{min}\{\mathrm{log}(n\widehat{M}),n^2\mathrm{log}n\})`$ time.

Submodular functions defined on modular lattices naturally arise in linear algebra. Minimization of such functions has a significant application to computing canonical forms of partitioned matrices . It remains an interesting open problem to develop an efficient algorithm for minimizing submodular functions on modular lattices, even for those specific functions that arise from partitioned matrices.

Independently of this work, and almost simultaneously, Schrijver has also developed a combinatorial, strongly polynomial-time algorithm for submodular function minimization . His algorithm also extends Cunningham’s approach. However, the resulting algorithm is quite different from ours.

## Acknowledgments

We are grateful to Bill Cunningham, Michel Goemans, and Maiko Shigeno for their useful comments.


text:
# Untitled Document

1 Introduction
The interacting dynamical nature of any economic sector composed of many cooperatively interacting agents, has many features in common with the statistical physics of interacting systems. In fact, economists like Pareto investigated the power law properties of the wealth distributions more than a century ago and Stigler studied the market fluctuations by employing Monte Carlo methods more than thirty years back . Motivated by these investigations, many efforts are being made recently to apply the statistical physics methods to various economic problems. A major part of the recent efforts has gone in investigating the nature of fluctuations and their distributions in the stock markets . We believe, however, a fundamental dynamics occurs in the money market, which affects strongly the dynamics of other sectors in the economy. An understanding of the statistical mechanics of the money market is therefore essential, and some of its features are very intriguing . Dragulescu and Yakovenko  have shown very recently that for any arbitrary and random sharing but locally conserving money transaction between any two agents in a market, the money distribution goes to the equilibrium Gibb’s distribution of statistical mechanics: $`P(m)=(1/T)\mathrm{exp}(m/T)`$ where $`T=M/N`$, the average money per agent in the market ($`M`$ is the total money of $`N`$ agents in the market). This equilibrium distribution is extremely robust and various kinds of inter-agent monetary transactions, which locally conserve the total money, lead to the Gibb’s probability distribution satisfying $`P(m_1)P(m_2)=P(m_1+m_2)`$. A major achievement of this study  has been the precise identification of the temperature $`T`$ as the average money per agent in the market. This is due to the fact that the probability distribution is normalized and the total money is conserved. This precise identification had been missing in many of the previous attempts  though the identification from fiscal policy considerations was indeed very close . They had considered a market consisting of $`N`$ agents, where initially each one gets any arbitrary share $`m_i`$ of the total money $`M`$ in the market ($`_{i=1}^Nm_i=M`$; $`M`$ and $`N`$ fixed). The “trade” dynamics goes as follows. Select any two arbitrary agents $`i`$ and $`j`$ with money $`m_i`$ and $`m_j`$, respectively. These two agents then exchange their money through some trade, keeping their total amount of money $`m_i+m_j`$ conserved and no debt is allowed ($`m_i0`$ at any stage of the trade). There is no other restriction in the trade. Extensive numerical simulations show that this and various modifications of trade, like multi-agent transactions, etc., all lead to the robust Gibb’s distribution, independent of the initial distribution the market starts with. So, most of the agents end-up in this market with very little money! Supply of more money in the market (increasing $`T`$) can increase the width of the distribution, but the most probable money for any agent in the market remains zero. It may be mentioned that Ispolatov et al. had studied earlier the (non-equilibrium) wealth distributions in various asset exchange models  where the trade dynamics do not have time-reversal symmetry. Of course, the study of Levy and Solomon  indicates a power law distribution for the wealth in realistic (possibly non-equilibrium) markets or societies.

Although the model of Dragulescu and Yakovenko  is a simple and interesting one to start with, a very natural ingredient for any realistic economic agent is missing in the model: no economic agent trades with the entire money he or she possesses without saving a part of it; saving propensity is too natural a tendency for any economic agent . This saving propensity of course varies from agent to agent and even with the accumulation of wealth of a single agent. There are also country-wise variations of this saving propensity. We note that the fraction of savings $`\lambda `$, called the “marginal propensity to save” by economists , remains fairly constant, independent of the agents. We have taken it to be a constant in the model considered here. We show that in presence of this saving factor $`\lambda (0)`$ the money market becomes truly cooperative in nature and “critical behaviour”  sets in. The multiplicative property of the Gibb’s distribution $`P(m)`$, discussed above, was responsible for the absence of any cooperative feature of the statistics for $`\lambda =0`$. Once the local or individual measure of saving propensity is introduced ($`\lambda 0`$), a global order emerges in the entire money market, giving a non-vanishing most probable money for each agent in the market and nontrivial critical behaviour of the resulting statistics in the equilibrium.
2 Model and simulation results for money distribution
We again consider a simple model of the closed economic system where the total amount of money $`M`$ is conserved and the number of economic agents $`N`$ is fixed. Each economic agent $`i`$, which may be an individual or a corporate, possesses money $`m_i`$. An economic agent can exchange money with any other agent through some trade, keeping the total amount of money of both the agents conserved. We assume that each economic agent saves a fraction $`\lambda `$ of its money $`m_i`$ before trading. We again assume that an agent’s money must always be non-negative and therefore no debt is permitted. Let us now consider that an arbitrary pair of agents $`i`$ and $`j`$ get engaged in a trade so that their money $`m_i`$ and $`m_j`$ change by amounts $`\mathrm{\Delta }m_i`$ and $`\mathrm{\Delta }m_j`$ to become $`m_i^{}`$ and $`m_j^{}`$, where $`\mathrm{\Delta }m_i`$ is a random fraction of $`(1\lambda )(m_i+m_j)`$ and $`\mathrm{\Delta }m_j`$ is the rest of it. Conservation of the total money in each trade is ensured, as earlier.

We performed computer simulations with fixed number of agents $`N`$ and total money $`M=NT`$. Most of our simulation results are for $`N=500`$ and $`T=100`$. However, we checked results for different $`N`$ values ($`250N1000`$) to check finite size effects, etc., and with different $`T`$ values (upto $`10000`$). Initially we divided the total money $`M`$ amongst $`N`$ agents equally so that $`m_i=M/N=T`$ for all $`i`$. We chose a fixed value of $`\lambda `$ between zero and unity and investigated its effects on the equilibrium distribution of money $`P(m)`$ in the market, giving the (normalized) number of agents $`P`$ with money $`m`$. We choose randomly two agents $`i`$ and $`j`$ having money $`m_i`$ and $`m_j`$, respectively. Then $`\mathrm{\Delta }m_i=ϵ(1\lambda )(m_i+m_j)`$ and $`\mathrm{\Delta }m_j=(1ϵ)(1\lambda )(m_i+m_j)`$, where $`ϵ`$ is a random number between zero and unity. Then $`m_i^{}=\lambda m_i+\mathrm{\Delta }m_i`$ and $`m_j^{}=\lambda m_j+\mathrm{\Delta }m_j`$ after the trade. Alternatively, this trade can also be viewed as $`m_im_{i}^{}{}_{}{}^{}`$, $`m_jm_{j}^{}{}_{}{}^{}`$ where $`m_{i}^{}{}_{}{}^{}=m_i\mathrm{\Delta }m`$, $`m_{j}^{}{}_{}{}^{}=m_j+\mathrm{\Delta }m`$ with $`\mathrm{\Delta }m=(1\lambda )[m_iϵ(m_i+m_j)]`$. This can be checked by straightforward substitution. These trades were repeated for large number of arbitrary choices of pairs of agents; each of these random choices of pairs is considered as one trade and consequently one time unit. The typical time upto which we run this algorithm is above $`5000`$ (sometimes upto $`50000`$ for larger $`\lambda `$). The probability distribution was determined after every $`50`$ time steps till a stationary distribution was obtained. We then took an average over $`2000`$ such stationary distributions to obtain a smooth distribution. We checked that the equilibrium distribution obtained is again extremely robust and does not depend at all on the initial distribution of money in the market. We investigated the nature of the equilibrium distribution $`P(m)`$ for various values of $`\lambda `$ ($`0\lambda <1`$). Apart from the stationary distribution, we also investigated the “relaxation” behaviour  of the distributions and obtained the time variations of $`P_1P(m=T)`$ till a steady behaviour was found. We define the relaxation time $`\tau _R`$ as the earliest time where $`P_1(t)`$ becomes practically independent of time.

The results for the equilibrium distribution $`P(m)`$ are shown in Fig. 1, for some values of $`\lambda `$. The inset shows that the equilibrium distribution is independent of the market size $`N`$ and the average money in the market or temperature $`T`$. The real money exchanged randomly in any trade is less than the total money, because of the saving by each agent. This destroys the multiplicative property of the distribution $`P(m)`$ (seen earlier for $`\lambda =0`$) and $`P(m)`$ changes from the Gibb’s form to the asymmetric Gaussian-like form as soon as a finite $`\lambda `$ is introduced. The $`\lambda =0`$ case was practically a random noise dominating one and therefore effectively a non-interacting market. Introduction of a finite amount of saving ($`\lambda 0`$), dictated by individual self-interest, immediately makes the money dynamics cooperative and the global ordering (in the distribution) is achieved. This kind of self-organization in the money market, coming out of pure self-interest of each agent, is reminiscent of the “invisible hand” effect  in the “free market” suggested originally by Adam Smith in 1776.

The relaxation behaviour of the distribution is shown in Fig. 2, where the time variation of $`P_1(t)`$ for different values of $`\lambda `$ is shown. Since $`m_i/T=1`$ for all $`i`$ to start with, $`P_1(t)`$ starts falling from unity in all the cases (not shown in Fig. 2). After some initial rapid decay, we see an extremely slow spin-glass type  ($`\mathrm{ln}t`$) relaxation behaviour. The inset shows the typical variation of the relaxation time $`\tau _R`$ with $`\lambda `$. The dynamics obviously becomes slower with increasing $`\lambda `$ and $`\tau _R`$ seems to diverge as $`\lambda `$ approaches unity. Precisely at $`\lambda =1`$, the dynamics of course stops.

An important feature of this humped distribution $`P(m)`$ at any non-vanishing $`\lambda `$ is the variation of the most probable money $`m_p(\lambda )`$ ( where $`P(m)`$ becomes maximum) of the agents. As discussed before, $`m_p=0`$ for $`\lambda =0`$ (Gibb’s distribution) and most of the economic agents in the market end-up losing most of their money. However, even with pure self-interest of each agent for saving a factor $`\lambda `$ of its own money at any trade, a global feature emerges: the entire market ends-up with a most-probable money $`m_p(\lambda )`$. This $`m_p(\lambda )`$ shifts in an interesting manner from $`m_p=0`$ (for $`\lambda =0`$) to $`m_pT`$ (for $`\lambda 1`$). We find that initially, for small $`\lambda `$, $`m_p(\lambda )`$ varies very closely as $`\lambda ^{1/2}`$ and then it crosses over at $`\lambda =\lambda _c0.45`$ to $`\lambda ^{1/3}`$ variation. This is shown in Fig. 3, where the curves for $`\lambda ^{1/2}`$ and $`\lambda ^{1/3}`$ are also indicated. We do not have any idea why this crossover in the exponent for $`\lambda `$ occurs at any finite $`\lambda (0.45)`$. Also the exponents $`1/2`$ and $`1/3`$ are very intriguing. We checked that this crossover point is perhaps also the point where the probability of an agent having zero money $`P_0P(m=0)`$ just disappears, as shown in the inset of Fig. 3. The half-width $`\mathrm{\Delta }m_p`$ and the peak height $`P_{m_p}`$ of the equilibrium distribution, which scale practically as $`(1\lambda )^{1/2}`$ and $`(1\lambda )^{1/2}`$ respectively (see Fig. 4), do not have any irregularity there. In fact, no other property of the equilibrium distribution $`P(m)`$ has any irregularity at $`\lambda _c`$. It may be mentioned that Ispolatov et al. also observed a singularity in the power law of non-equilibrium growth distribution of wealth in their multiplicative asset exchange model having broken time-reversal symmetric trade dynamics.

We noted that although the total distribution assumes some global cooperative feature, each individual’s money $`m_i`$ fluctuates randomly. Fig. 5 shows the time variation of the money of an arbitrarily chosen individual in the market for two different values of $`\lambda `$. The inset shows the variation of the time-averaged money $`m_i`$ of the agent and its fluctuation $`\mathrm{\Delta }m_i\sqrt{(m_im_i)^2}`$ with $`\lambda `$ after relaxation ($`t>\tau _R`$). Since the total money is conserved, $`m_i`$ remains constant ($`=T`$) here, while $`\mathrm{\Delta }m_i`$ goes down with $`\lambda `$ as $`(1\lambda )`$. This is because at any time the agents keep a fixed fraction of their individual money and receive a random fraction of the money traded proportional to $`(1\lambda )`$.

3 Summary and discussion
We thus considered here a very simple model of money market where the total money $`M`$ and the number of agents $`N`$ is fixed. Each pair of arbitrarily chosen agents in the market exchange money with each other through a trade (each trade considered as one time unit). During the trade, each agent saves a fraction $`\lambda `$ of its own money at that time and exchanges randomly out of the remaining money, conserving the total amount of money and not allowing any debt. We find that for $`\lambda =0`$, the market effectively becomes non-interacting and, no matter what initial distribution the agents start with, the resulting money distribution becomes the equilibrium Gibb’s distribution , where most of the agents in the market end-up with very little money. For any non-vanishing $`\lambda `$, the equilibrium distribution becomes the asymmetric Gaussian-like with the most probable money $`m_p`$ (corresponding to the peak in $`P(m)`$) shifting away from $`m=0`$ with increasing $`\lambda `$. This global feature, coming out of the individual self-interest of saving a part of its own money, can be considered to be a demonstration of the self-organization in the market suggested long ago by Adam Smith (“invisible hand” ). Apart from this we find intriguing singularities appearing in the equilibrium distribution $`P(m)`$: $`m_p\lambda ^\nu `$ where $`\nu =1/2`$ for $`\lambda <\lambda _c(0.45)`$ and $`\nu =1/3`$ for $`\lambda _c<\lambda <1`$. Also, $`\mathrm{\Delta }m_p(1\lambda )^{1/2}`$ and $`P_{m_p}(1\lambda )^{1/2}`$.

These singularities in the equilibrium distribution come obviously from the cooperative nature of the market interactions induced by the saving propensity of the agents. It may be mentioned that while such singular behaviour in the equilibrium money distribution is very natural here, in the corresponding physical (gas) system of Newtonian particles one gets regular distributions, e.g., the Gibb’s distribution (or for that matter, Bose or Fermi distributions for quantum particles) and never any singularity, because of the absence of any physical equivalent of the “saving” factor there. We thus believe, while regular distributions are common for minimally interacting physical many-body systems in equilibrium, singular distributions are typical of any working model of the markets.
Acknowledgements:  We are grateful to A. Dutta, S. S. Manna, S. Pradhan and P. Sen for useful discussions, and S. Redner for some useful communications. We also thank the anonymous referee for suggesting important corrections and improvement of the manuscript.

References
e-mail addresses :

(1)anirban@cmp.saha.ernet.in

(2)bikas@cmp.saha.ernet.in

1. V. Pareto, Cours d’Economie Politique, Lausanne and Paris (1897); G. J. Stigler, J. Business, 37, 117 (1964).
2. See e.g., R. N. Mantegna and H. E. Stanley, An Introduction to Econophysics, Cambridge University Press, Cambridge (2000).
3. B. K. Chakrabarti and S. Marjit, Ind. J. Phys., 69 B, 681 (1995).
4. S. Moss de Oliveira, P. M. C. de Oliveira and D. Stauffer, Evolution, Money, War and Computers, B. G. Teubner, Stuttgart (1999); P. W. Anderson, K. J. Arrow and D. Pines (Eds), The Economy as an Evolving Complex system, Addison-Wesley, Redwood City (1988).
5. A. Dragulescu and V. M. Yakovenko, cond-mat/0001432.
6. S. Ispolatov, P. L. Krapivsky and S. Redner, Euro. Phys. J. B, 2, 267-276 (1998).
7. M. Levy and S. Solomon, Physica A, 242, 90 (1997).
8. P. A. Samuelson, Economics, $`11`$th Edition, McGraw-Hill Inc., Auckland, 195-208 (1980).
9. M. Plischke and B. Bergersen, Equilibrium Statistical Physics, 2nd Edition, World Scientific, Singapore (1994); S. K. Ma, Statistical Mechanics, World Scientific, Singapore (1985).
10. See e.g., J. M. Keynes, The General Theory of Employment, Interest and Money, The Royal Economic Society, Macmillan Press, London (1973).

Figure captions
Fig. 1 : The stationary money distribution $`P(m)`$ versus money $`m`$ for different saving propensity factor $`\lambda `$ ($`N=500,T=100`$). The inset shows the same for different $`N`$ and $`T`$ at two fixed values of $`\lambda `$. The superposition of data indicates the absence of finite size effects and the temperature independence.
Fig. 2 : Relaxation of the distribution: $`P_1(t)`$ versus the time $`t`$ (number of trades) for different $`\lambda `$. Because of extremely slow relaxation, we show in the $`\mathrm{ln}t`$ scale for the later stage only. The vertical arrows indicate the relaxation time $`\tau _R`$. The inset shows the typical variation of $`\tau _R`$ with $`\lambda `$; $`\tau _R`$ diverges as $`\lambda 1`$.
Fig. 3 : The variation of the most probable money $`m_p`$ of an agent (peak position of the equilibrium distribution $`P(m)`$) as a function of $`\lambda `$. The dotted and dashed curves correspond to $`\lambda ^{1/2}`$ and $`\lambda ^{1/3}`$ respectively. The crossover point ($`\lambda _c`$) is indicated by the vertical arrow. The inset shows that the probability $`P_0`$ of agents with zero money also first disappears at almost the same point ($`\lambda _c`$) indicated by the vertical arrow.
Fig. 4 : The variation of the distribution half-width $`\mathrm{\Delta }m_p`$ as function of $`\lambda `$. The dashed curve corresponds to $`\sqrt{2.0(1\lambda )}`$. The inset shows the peak height $`P_{m_p}`$ as function of $`\lambda `$. The dotted curve here corresponds to $`1/\sqrt{2.4(1\lambda )}`$.
Fig. 5 : Time variation of money of an individual agent $`m_i`$ for two different values of $`\lambda `$. The inset shows the variations of $`m_i`$ and its fluctuations $`\mathrm{\Delta }m_i`$ as functions of $`\lambda `$.