# 1 Introduction

## 1 Introduction

The set covering problem (SCP) is a well known NP-hard combinatorial optimization problem, which represents many real-world resource allocation problems. Exact solutions can be obtained by e.g. a branch-and-bound approach for modestly sized problems. For larger problems various approximative schemes have been suggested (see e.g. ). In this paper we develop a novel approach based on feedback Artificial Neural Networks (ANN), derived from the mean field approximation to the thermodynamics of spin systems.

ANN is a computer paradigm that has gained a lot of attention during the last 5-10 years. Most of the activities have been directed towards feed-forward architectures for pattern recognition or function approximation. ANN, in particular feedback networks, can also be used for difficult combinatorial optimization problems (e.g. -). Here ANN introduces a new method that, in contrast to most existing search and heuristics techniques, is not based on exploratory search to find the optimal configuration. Rather, the neural units find their way in a fuzzy manner through an interpolating, continuous space towards good solutions. There is a close connection between feedback ANN and spin systems in statistical physics. Consequently, many mathematical tools used for dealing with spin systems can be applied to feedback ANN. Two steps are involved when using ANN for combinatorial optimization:

1. Map the problem onto an energy function, e.g.

$$E\left(𝒮\right)=\frac{1}{2}\underset{ij}{}w_{ij}s_is_j$$
(1)   
where $`𝒮=\left\{s_i;i=1\mathrm{}N\right\}`$ is a set of binary spin variables $`s_i\{0,1\}`$, representing the elementary choices involved in minimizing $`E`$, while the weights $`w_{ij}`$ encode the costs and constraints.
2. To find configurations with low $`E`$, iterate the mean field (MF) equations

$$v_i=1/\left(1+\mathrm{exp}\left(\frac{1}{T}\frac{E\left(𝒱\right)}{v_i}\right)\right)$$
(2)   
where $`T`$ is a fictitious temperature while $`𝒱=\left\{v_i\right\}`$, where $`v_i[0,1]`$ represents the thermal average $`s_i_T`$, and allows for a probabilistic interpretation.

Eqs. (1,2) only represent one example. More elaborate encodings have been considered, e.g. based on Potts spins allowing for more general basic decisions elements than simple binary ones . A propagator formalism based on Potts neurons has been developed for handling topological complications in e.g. routing problems .

The ANN approach for SCP that we develop here differs from the one that was successfully used for the somewhat related knapsack problem in , in particular with respect to encoding the constraints. Whereas a non-linear step-function was used in , we will here use a multilinear penalty, which in addition to being theoretically more appealing, also appears to be very efficient. Furthermore, an automatic procedure for setting the relevant $`T`$-scale is devised.

The algorithm is extensively tested against a set of publicly available benchmark problems  with sizes (rows$`\times `$columns) ranging from 200$`\times `$1000 to 5000$`\times `$10<sup>6</sup>. The approach yields results, typically within a few percent from the exact optimal solutions, for sizes where these are available. Comparisons with other approximate methods also come out well. The algorithm is extremely rapid – the typical CPU demand is only a few seconds (on a 400MHz Pentium II).

A public domain WWW server has been set up, where arbitrary problems can be solved interactively.

This paper is organized as follows: In Sect. 2 we define the set covering problem, and in Sect. 3 we describe its encoding in terms of a neural network energy function and discuss the mean field treatment. Sect. 4 contains numerical explorations and comparisons. A brief summary is given in Sect. 5. Appendices A and B contain a derivation of the mean field equations, and some algorithmic implementation issues, respectively. Tables from the numerical explorations are found in Appendix C, while Appendix D contains pointers and instructions for the WWW server. It should be stressed that this paper is self-contained – no prior knowledge of feedback neural networks or the mean field approximation is necessary.

## 2 The Set Covering Problem

The set covering problem (SCP) is the problem of finding a subset of the columns of an $`M`$x$`N`$ zero-one matrix $`A=\{a_{ki}\{0,1\};k=1,\mathrm{},M;i=1,\mathrm{},N\}`$ that covers all rows at a minimum cost, based on a set of column costs $`\{c_i;i=1,\mathrm{},N\}`$. SCP is conveniently described using a set of binary variables, $`𝒮=\{s_i\{0,1\};i=1,\mathrm{},N\}`$. More precisely, SCP is defined as follows:

Minimize $`{\displaystyle \underset{i=1}{\overset{N}{}}}c_is_i`$ (3)   
subject to $`{\displaystyle \underset{i=1}{\overset{N}{}}}a_{ki}s_i1,k=1,\mathrm{},M`$ (4)   
with $`s_i\{0,1\}i=1,\mathrm{},N`$ (5)   
Eq. (4) states that at least one column must cover a particular row. The special case where all costs $`c_i`$ are equal is called the unicost SCP. There is subclass of SCPs that has a nice graphical interpretation: If the zero-one matrix $`A`$ has the property that each row contains exactly two 1’s then we can interpret $`A`$ as the vertex-edge matrix of a graph with $`N`$ vertices and $`M`$ edges. The unicost SCP is then a vertex covering problem, where the task is to find the minimal number of vertices that covers all edges of the graph. SCPs (including weighted vertex covering) are NP-hard combinatorial optimization problems. If the inequalities of Eq. (4) are replaced by equalities, one has the set partition problem (SPP). Both SCP and SPP have numerous resource allocation applications.

The following quantities will be used later:

Density: $`\rho ={\displaystyle \frac{1}{MN}}{\displaystyle \underset{ki}{}}a_{ki}`$ (6)   
Column Sums: $`\widehat{N}_i={\displaystyle \underset{k}{\overset{M}{}}}a_{ki}`$ (7)   
Row Sums: $`\widehat{M}_k={\displaystyle \underset{i}{\overset{N}{}}}a_{ki}`$ (8)   
## 3 The Mean Field Approach

### 3.1 The Energy Function for SCP

We start by mapping the SC problem of Eqs. (3, 4, 5) onto a spin energy function $`E(𝒮)`$ (step 1 in the introduction),

$$E(𝒮)=\underset{i=1}{\overset{N}{}}c_is_i+\alpha \underset{k=1}{\overset{M}{}}\underset{i=1}{\overset{N}{}}\left(1a_{ki}s_i\right)$$
(9)   
The first term yields the total cost and the second one represents the covering constraint of Eq. (4) by imposing a penalty if a row is not covered by any column.

The constraint term is a multilinear polynomial in the spin variables $`s_i`$, i.e. it is a linear combination of products $`s_1s_2\mathrm{}s_K`$ of distinct spins. This is attractive from a theoretical point of view, and the differentiability of $`E`$ enables a more quantitative analysis of the dynamics of the mean field algorithm.

An alternative would be to implement the inequality constraints using a piecewise linear function ,

$$\alpha \underset{k=1}{\overset{M}{}}\varphi \left(1\underset{i=1}{\overset{N}{}}a_{ki}s_i\right)$$
(10)   
with $`\varphi (x)=x\mathrm{\Theta }(x)=x`$ if $`x>0`$ and 0 otherwise. This yields a non-differentiable energy function and generally an inferior performance as compared to the polynomial representation (9).

### 3.2 The Statistical Mechanics Framework

The next step is to minimize $`E(𝒮)`$. Using some local updating rule will most often yield a local minimum close to the starting point, with poor solutions as a result. Simulated annealing (SA)  is one way of escaping from local minima since it allows for uphill moves in $`E`$. In SA a sequence of configurations $`𝒮`$ is generated according to a stochastic algorithm, such as to emulate the probability distribution

$$P(𝒮)=\frac{e^{E\left(𝒮\right)/T}}{{\displaystyle \underset{𝒮^{}}{}}e^{E\left(𝒮^{}\right)/T}}$$
(11)   
where the sum runs over all possible configurations $`𝒮^{}`$. The parameter $`T`$ (temperature) acts as a noise parameter. For large $`T`$ the system will fluctuate heavily since $`P(𝒮)`$ is very flat. For the SCP this implies that the sequence contains mostly poor and infeasible solutions. On the other hand, for a small $`T`$, $`P(𝒮)`$ will be narrow, and the sequence will be strongly dependent upon the initial configuration and contain configurations only from a small neighborhood around the initial point. In SA one generates configurations while lowering $`T`$ (annealing), thereby diminishing the risk of ending up in a suboptimal local minimum. This is quite CPU-consuming, since one has to generate many configurations for each temperature following a careful annealing schedule (typically $`T_k=T_0/\mathrm{log}(1+k)`$ for some $`T_0`$) in order to be certain to find the global minimum.

In the mean field (MF) approach the costly stochastic SA is approximated by a deterministic process. MF also contains an annealing procedure. The original binary variables $`s_i`$ are replaced by continuous mean field variables $`v_i[0,1]`$, with a dynamics given by iteratively solving of the MF equations for each $`T`$.

An additional advantage of the MF approach is that the continuous mean field variables can evolve in a space not accessible to the original variables. The intermediate configurations at non-zero $`T`$ have a natural probabilistic interpretation.

### 3.3 Mean Field Theory Equations

Our objective is now to minimize $`E`$ using the MF method. The binary spin variables $`s_i`$ are replaced by mean field variables $`v_i`$, representing solutions to the MF equations (a derivation of these is given in Appendix A).

$$v_i=1/\left(1+\mathrm{exp}\left(\frac{1}{T}\mathrm{\Delta }E_i\left(𝒱^{(i)}\right)\right)\right),i=1,\mathrm{},N$$
(12)   
where

$$\mathrm{\Delta }E_i\left(𝒱^{(i)}\right)=E\left(𝒱^{(i)},v_i=1\right)E\left(𝒱^{(i)},v_i=0\right)=E/v_i$$
(13)   
The set $`𝒱^{(i)}`$ denotes the complementary set $`\{v_j,ji\}`$. From Eq. (9) we get

$$\mathrm{\Delta }E_i\left(𝒱^{(i)}\right)=c_i\alpha \underset{k=1}{\overset{M}{}}a_{ki}\underset{ji}{\overset{N}{}}\left(1a_{kj}v_j\right)$$
(14)   
The MF equations are solved iteratively while annealing in $`T`$. It should be noted that the equation for $`v_i`$ does not contain any feedback, i.e no explicit dependence on $`v_i`$ itself. This yields a smooth convergence, and it is usually only necessary with a few iterations at each value of $`T`$. What remains to be specified are the parameters $`\alpha `$ and $`T`$. The latter will be discussed next and we return to the choice of $`\alpha `$ in connection with the numerical explorations in Sect. 4.

### 3.4 Critical Temperatures

In the limit of high temperatures, $`T\mathrm{}`$, the MF variables $`\{v_i\}`$ will, under the dynamics defined by iteration of Eq. (12), converge to a trivial symmetric fixed point with $`v_i=1/(1+\mathrm{exp}(0))=\frac{1}{2}`$, corresponding to no decision taken. At a finite but high $`T`$, the corresponding fixed point will typically deviate slightly from the symmetric point.

For many problems, a bifurcation occurs (indicative of a transition from a disordered phase to an ordered one) at a critical temperature $`T_c`$, where the trivial fixed point loses stability and other fixed points emerge, which as $`T0`$ converge towards definitive candidate solutions to the problem, in terms of $`v_i\{0,1\}`$. For some problems, a cascade of bifurcations occur, each at a distinct critical temperature, but in the typical case there is a single bifurcation.

It is then of interest to estimate of the position of $`T_c`$, which defines a suitable starting point for the MF algorithm. Such an estimate can be obtained by means of a linear stability analysis for the dynamics close to the fixed point. For the special case of a symmetric unicost SCP with constant row and column sums for the matrix $`A`$, there will be a sharp transition around

$$T_c\alpha \rho ^2M2^{\rho N}$$
(15)   
where $`\rho `$ is defined in (6).

For non-unicost problems, $`T_c`$ is harder to estimate, and there might even be no bifurcation at all. For such problems, a suitable initial $`T`$ is instead determined by means of a fast preliminary run of the algorithm (see below).

## 4 Numerical Explorations

### 4.1 Implementation Details

The annealing schedule for $`T`$ and the value of the constraint parameter $`\alpha `$ have to be determined before we can run the algorithm. The former is accomplished by a geometric decrease of $`T`$,

$$T_{t+1}=kT_t$$
(16)   
where $`k`$ is set to 0.80 and $`T_0`$ is determined by a fast prerun of the algorithm (see below). The number of iterations of Eqs. (12) for each value of $`T`$ is not fixed; $`T`$ is lowered only when all $`v_i`$ have converged.

In order to ensure a valid solution at low $`T`$, the size of the constraint term in Eq. (14) must be larger than the largest cost $`c_{max}`$ that is part of the solution. Using a too large $`\alpha `$ will however reduce the solution quality since $`E`$ is then dominated by the covering constraint. Our choice of $`\alpha `$ therefore depends on $`c_{max}`$. Ideally, $`M/(\rho M)=1/\rho `$ columns would suffice to cover each row of $`A`$. If we further optimistically assume that the $`1/\rho `$ smallest costs can be chosen for the solution, $`c_{max}`$ can easily be found. However, for most problems we need more than $`1/\rho `$ columns, which makes it difficult to estimate $`c_{max}`$ except for unicost problems where all column costs are equal.

To determine (an approximate) $`c_{max}`$ for a non-unicost SCP, we perform a fast prerun with a smaller annealing factor $`k=0.65`$, and with $`\alpha =1.01`$. From this prerun one can also obtain an estimate the critical temperature $`T_c`$ as the $`T`$ where the saturation (defined below) deviates from 0. The second run is then initated at $`T_0=2T_c`$, thereby avoiding unnecessary updates of $`v_i`$ at high $`T`$.

This procedure for setting $`\alpha `$ also requires rescaling of the costs; for all problems we set $`c_ic_i/\mathrm{max}_j(c_j)`$. See Table 1 for a summary of the parameters used.

The evolution of the MF variables $`\{v_i\}`$ is conveniently monitored by the saturation $`\mathrm{\Sigma }`$,

$$\mathrm{\Sigma }=\frac{4}{N}\underset{i=1}{\overset{N}{}}\left(v_i1/2\right)^2$$
(17)   
A completely “undecided” configuration, $`\mathrm{\Sigma }=0`$, means that every $`v_i`$ has the value $`1/2`$. During the annealing process, as the $`v_i`$ approach either 1 or 0, $`\mathrm{\Sigma }`$ converges to 1. The transition between $`\mathrm{\Sigma }=0`$ and $`\mathrm{\Sigma }=1`$ is usually smooth for a generic SCP. However, when a bifurcation is encountered (see above), $`\mathrm{\Sigma }`$ can change abruptly; this occurs e.g. for unicost SCP where there is no natural ordering among the costs $`\{c_i\}`$.

Fig. 1 shows the evolution of the mean field variables $`\{v_i\}`$ and the saturation $`\mathrm{\Sigma }`$ for the problems 4.1 and cyc09; the latter is a unicost problem (see Appendix C), and as can be seen from Fig. 1b, it clearly exhibits a bifurcation.

A summary of the algorithm can be found in Appendix B, while Appendix D gives the address of and instructions for a WWW server, where the MF algorithm is applied to user-defined SCPs.

### 4.2 Numerical Results

The performance of the algorithm is evaluated using 16 problem sets found in the OR-Library benchmark database . These 16 sets consist of 91 problems, out of which 19 are unicost SCP. The algorithm is coded in C and the computations are done on a 400MHz Pentium II PC. The details of the OR-library problems are given in Table C1 and our results can be found in Tables C2, C3 and C4 in Appendix C. The optimal or currently best known values are taken from refs. .

For each of the test SCPs, 10 trials of our algorithm are performed. In Tables C2-C4, the best and the average costs are listed for each problem. For 10 of the problems our method found the optimal solution. The MF results typically are within a few percent of the optimal solutions, as can be seen in Table 2. Large relative deviations from optimum are seen in the unicost CLR problems, which appear to be difficult for our approach. However, the optimal (integer) costs for these problems are low (23 - 26); this gives a large effect on the relative deviations even for a small change in the found costs. Decreasing the obtained costs by unity will change the mean relative deviation from $`16\%`$ to $`10\%`$. It is also important to notice that we use a common set of algorithm parameters ($`\alpha ,k,T_0`$) for all unicost problems, without any parameter optimization for each problem. Another set of parameters might be more advantageous for the CLR problems.

In our implementation of the MF algorithm, the time $`\tau `$ for a complete update of all variables $`\{v_i\}`$ scales approximately linearly with the number of non-zero entries in $`A`$,

$$\tau NM\rho $$
(18)   
This appealing property is feasible due to efficient calculations of $`\mathrm{\Delta }E_i`$ in Eq. (14) that utilizes the sparse nature of $`A`$. For the total solution time, $`\tau `$ should be multiplied by the number of iterations needed – empirically around 100, independently of problem size. Fig. 2 shows the mean solution time versus $`NM\rho `$.

Compared to other heuristic approaches, ours does not find the optimal cost as often as e.g the genetic algorithm . It is however very competitive with respect to speed. In ref.  nine different approximation algorithms were tested on a large number of unicost problems (including the set considered here); our approach is comparable to the top ones in both performance and speed.

## 5 Summary

We have developed a mean field feedback neural network approach for solving set covering problems. The method is applied to a standard set of benchmark problems available in the OR-library database. The method is also implemented in a public WWW server.

The inequality constraints involved are conveniently handled by means of a multi-linear penalty function that fits nicely into the mean field framework. The bifurcation structure of the mean field dynamics involved in the neural network approach is analyzed by means of a linearized dynamics. A simple and self-contained derivation of the mean field equations is provided.

High quality solutions are consistently found throughout a range of problem sizes ranging up to $`5\times 10^3`$ rows and $`10^6`$ columns for the OR-library problems without having to fine-tune the parameters, with a time consumption scaling as the number of non-zero matrix elements. The approach is extremely efficient, typically requiring a few seconds on a Pentium II 400 MHz computer.

The mean field approach to SCP can easily be modified to apply to the related, and more constrained, set partitioning problem. One simply has to replace the inequality constraint term with one that handles the equality constraint present in the set partitioning problem.

### Acknowledgments:

This work was in part supported by the Swedish Natural Science Research Council, the Swedish Board for Industrial and Technical Development and the Swedish Foundation for Strategic Research.

## Appendix Appendix A. Mean Field Approximation

Here follows for completeness a derivation of the mean field equations (see e.g. ). Let $`E\left(𝒮\right)`$ be an energy function of a set of binary decision variables (spins) $`𝒮=\{s_i|s_i\{0,1\},i=1,\mathrm{},N\}`$. If we assume a Boltzmann probability distribution for the spins, the average $`s_i_T`$ will be given by

$$s_i_T=\frac{{\displaystyle \underset{𝒮}{}}s_i\mathrm{exp}\left(E(𝒮)/T\right)}{{\displaystyle \underset{𝒮}{}}\mathrm{exp}\left(E(𝒮)/T\right)}$$
(A1)   
where the sums run over all possible configurations $`𝒮`$. We can manipulate this expression to obtain,

$$s_i_T=\frac{{\displaystyle \underset{𝒮}{}}\mathrm{exp}\left(E(𝒮)/T\right){\displaystyle \frac{{\displaystyle \underset{s_i=0,1}{}}s_i\mathrm{exp}\left(E(𝒮^{(i)},s_i)/T\right)}{{\displaystyle \underset{s_i=0,1}{}}\mathrm{exp}\left(E(𝒮^{(i)},s_i)/T\right)}}}{{\displaystyle \underset{𝒮}{}}\mathrm{exp}\left(E(S)/T\right)}$$
(A2)   
where $`𝒮^{(i)}`$ denotes the set $`\{s_j,ji\}`$ of all spins but $`s_i`$. If we now perform the sums over $`s_i`$ in the numerator, we get

$`s_i_T`$ $`=`$ $`{\displaystyle \frac{{\displaystyle \underset{𝒮}{}}\mathrm{exp}\left(E(𝒮)/T\right){\displaystyle \frac{1}{1+\mathrm{exp}\left(\mathrm{\Delta }E_i/T\right)}}}{{\displaystyle \underset{𝒮}{}}\mathrm{exp}\left(E(𝒮)/T\right)}}{\displaystyle \frac{1}{1+\mathrm{exp}\left(\mathrm{\Delta }E_i/T\right)}}_T`$ (A3)   
$`\text{where}\mathrm{\Delta }E_i(𝒮^{(i)})`$ $`=`$ $`E(𝒮^{(i)},s_i=1)E(𝒮^{(i)},s_i=0)`$ (A4)   
So far there are no approximations, the expression for $`s_i_T`$ has just been rewritten. Eq. (A3) states that the expectation value of $`s_i`$ is equal to the expectation value of a nonlinear function $`f`$ of all the other spins. The mean field approximation consists of approximating the expectation value $`f\left(𝒮^{(i)}\right)`$ by $`f\left(𝒮^{(i)}\right)`$. With $`v_i`$ denoting $`s_i`$, and $`𝒱^{(i)}`$ the complementary set $`\{v_j,ji\}`$, this amounts to making the replacement

$$\frac{1}{1+\mathrm{exp}\left(\mathrm{\Delta }E_i\left(𝒮^{(i)}\right)/T\right)}_T\frac{1}{1+\mathrm{exp}\left(\mathrm{\Delta }E_i\left(𝒱^{(i)}\right)/T\right)}$$
(A5)   
in Eq. (A3). This results in a set of self-consistency equations for $`𝒱`$, the mean field equations

$$v_i=\frac{1}{1+\mathrm{exp}\left(\mathrm{\Delta }E_i\left(𝒱^{(i)}\right)/T\right)},i=1\mathrm{}N$$
(A6)   
which in general must be solved numerically, e.g. by iteration.

The mean field approximation often becomes exact in the limit of infinite range interactions where each spin variable interacts with all the others. This can be seen from $`\mathrm{\Delta }E_i(𝒮^{(i)})`$ which then becomes a sum of many (approximately) independent random numbers, and a central limit theory can be applied.

## Appendix Appendix B. Algorithm Details

Here follows a summary of the MF annealing algorithm for finding approximate solutions to set covering problems. The procedure presented below is used for all problems in this study except for the large rail problems where numerical problems caused by limited machine precision comes into play. The problem arises when calculating the product in Eq. (14), which for large problems can contain many factors. A work-around is implemented by a simple truncation,

$$v_i:=\{\begin{array}{cc}0\hfill & \text{if }v_i<0.05\hfill \\ v_i\hfill & \text{otherwise}\hfill \end{array}$$
(B1)   
This numerical fix is only used when calculating the $`\mathrm{\Delta }E_i\left(𝒱^{(i)}\right)`$ in Eq. (14).

Algorithmic outline of our approach:

## Appendix Appendix C. Benchmark Results

## Appendix Appendix D. WWW Server

A program executing the mean field algorithm for set covering problems as presented in this paper can be publicly used by means of a World Wide Web server. A user can interactively submit a file defining an instance of SCP, and obtain the found solution. The URL of the WWW server is:

http://www.thep.lu.se/complex/mf\_server.html

An instance of SCP is defined by specifying the costs $`c_i`$ and the matrix $`A`$. Two formats, row and column ordering, are supported for the file that lists $`c_i`$ and the none-zero entries of $`A`$; they are defined as follows:

The column sums $`\widehat{N}_i`$ and row sums $`\widehat{M}_k`$ are defined in Eqs. (7-8). As an example, consider the SCP instance defined by

$$\stackrel{}{c}=(1,2,3,4,5),A=\left(\begin{array}{ccccc}1& 0& 1& 0& 1\\ 0& 1& 0& 1& 0\\ 1& 1& 0& 0& 1\\ 0& 0& 1& 1& 1\end{array}\right)$$
(D1)   
for which the column and row ordering formats read:

Row Ordering

4 5

1 2 3 4 5

3 1 3 5

2 2 4

3 1 2 5

3 3 4 5

Column Ordering

4 5

1 2 1 3

2 2 2 3

3 2 1 4

4 2 2 4

5 3 1 3 4

The solver returns the found cost (energy) $`E`$ (Eq.(3)), together with a characterization of the problem. Upon request, a file that lists the columns used in the solution is also provided.

There is a limitation ($`B`$) on the size of the problems that can be submitted to the server. Instances of SCP with $`MN\rho >B`$ will not be considered. Presently $`B`$, which is limited by the available memory of the server, is given by $`3x10^6`$.