Okay, once I looked into simulated
annealing and here just read all the
comments on this thread. So, apparently
what I have to contribute is not here yet.

Here I outline non-linear duality theory
and Lagrangian relaxation and report a
war story where this approach badly beat
some efforts with simulated annealing.

Basically non-linear duality theory has
some surprising properties commonly just
too powerful to ignore.

Basically we're considering a case of
optimization or mathematical
programming.

So, let R denote the set of real numbers,
and suppose positive integers m and n are
given. Let R^n be the usual n-dimensional
real vector space of n-tuples. To use
matrix notation, we regard each n-tuple x
in R^n as a matrix with n rows and one
column, that is, as n x 1.

Similarly for R^m.

We are given set C a subset of R^n and
functions

f: R^n --> R

g: R^n --> R^m

We seek x in R^n to solve non-linear
program (NLP)

min z = f(x)

subject to

g(x) <= 0

x in C

where the 0 here is the m x 1 matrix of
zeros.

We let the feasible region R, a subset of
R^n, be

F = { x | x in C and g(x) <= 0 }

We regard g as m x 1 where for j = 1, 2,
..., m component j of g is

g_j: R^n --> R

where here g_j, borrowing TeX notation, is
g with subscript j.

Okay, given 1 x m l, we let the
Lagrangian

L( x, l ) = f(x) + l g(x)

Then for l >= 0 and x in F,

f(x) >= L( x, l )

and

z >= L( x, l )

Then the dual is to find l to maximize

M( l ) = min_{x in C} L( x, l )

Then M( l ) <= z to that M( l ) is a lower
bound on z.

In particular (H. Everett), if g(x) = 0, l
>= 0, and

z >= M( l ) = min_{x in C} L( x, l ) =

min_{x in C} f(x) + l g(x) = f(x)

and x in F so that x solves NLP.

A practical special case is essentially if
spend all the productive resources and do
so optimally, then are optimal.

Of course, the usual notation used for l
is lambda, and for a while around DC
Everett ran Lambda Corporation which did
resource allocation problems for the US
DoD.

Theorem: M is concave.

Proof: For t in the interval [0,1] and m
x 1 l_1 and l_2, we have

M( tl_1 + (1 - t)l_2) =

min_{x in C} L( x, t l_1 + (1 - t) l_2) =

min_{x in C} f(x) + ( t l_1 + (1 - t) l_2) g(x) =

min_{x in C} ( t + (1 - t)) f(x) + ( t l_1 + (1

- t) l_2) g(x) =

min_{x in C} t L( x, l_1 ) + (1 - t) L( x, l_2 )

>= t min_{x in C} L( x, l_1 ) + (1 - t) min_{x in C} L(
x, l_2 )

= t M( l_1 ) + (1 - t) M( l_2 )

so that M is concave. Done.

Immediately we have supporting hyperplanes
for our concave M:

Suppose y and l are given and

M( l ) = L( y, l )

Then for 1 x m u,

M( l + u ) = min_{x in C} L( x, l + u )

= min_{x in C} ( L( x, l ) + u g(x) )

<= L( y, l ) + u g(x)

= M( l ) + u g(x)

so that for 1 x m u

s( u ) = M( l ) + u g(x)

is a supporting hyperplane for concave
M( l ); that is,

M( l + u ) <= s( u ) = M( l ) + u g(x)

and s( 0 ) = M( l + 0 ) = M( l ).

So, here we have the basic theory of
non-linear duality as needed by the
iterative algorithmic approach of
primal-dual Lagrangian relaxation.

War Story:

Once some guys had a big resource
application problem having to do with a
suddenly legal case of cross selling for
banks. They had formulated their problem
as a 0-1 integer linear program with
~40,000 constraints and 600,000 variables.

They had tried simulated annealing, ran
for days, and quit with no idea how close
they were to optimality.

I looked at their formulation and found
that, except for 16 of the constraints,
the problem was comparatively easy.

So, I wrote out

L( x, l ) = f(x) + l g(x)

where l was 1 x 16 and g represented the
16 troublesome constraints. The set C
represented the 0-1 constraints and the
rest of the 40,000 constraints.

Then for each l, finding x to solve

M( l ) = min_{x in C} L( x, l )

was easy. Then I just had to find l to
solve

max M( l )

subject to

l >= 0

and that was just to maximize a concave
function.

So, I did some primal-dual iterations:

Set

k = 1

Pick l^k >= 0

Step 1 (Primal):

Find x = x^k to solve

min_{x in C} L( x, l^k )

Then

M( l^k ) = L( x^k, l^k )

Step 2 (Termination)

If x^k in F then

M( l^k ) <= z <= f( x^k )

so that if

M( l^k ) - f( x^k )

is sufficiently small, then
terminate. Else if k i
to large, terminate.

Step 3 (Dual):

Find 1 x m u and w in R to solve

max w

subject to

w <= M( l^p ) + u g(x^p)

p = 1, 2, ..., k

Of course, this is just a linear
program. Set

l^(k + 1) = l^K + u

k = k + 1

Go to Step 1.

There are refinements: E.g., for the
linear program to find l^{k + 1}, can hope
to get better numerical performance with
by using central cutting planes.

For the real problem, in 905 seconds on a
90 MHz processor, I ran 500 primal-dual
iterations and, then, had a feasible
solution, from the bound inequalities,
guaranteed to be within 0.025% of
optimality.

Worked fine! So, simulated annealing is
not the only tool in the tool box!
Sometimes Lagrangian relaxation can work
well, too!