Simulated annealing (sa) for baroque numbers--a warm up

The idea of applying the simulated annealing method to baroque numbers was with me for a long time. Finally I had invited others on pl.sci.matematyka (on 2007-05-19) to implement it in parallel with me (independently). I wanted to have some company, it'd be interesting. Somehow, after some positive (and negative) feedback, I ended up alone with a first version 2007-July. It was able to repeat the results obtained in 1990-ies. It couldn't do more because I used a 32-bit C++, and I have represented prime powers directly as actual integers, thus limiting them to the range p^e <>32 (I even stayed lazily withing p^e ≤ 2³¹). Then I barely started to work on a more advanced version 2007-Aug before I already had to stop (for reasons not realated directly to the project). Now, (starting near the end of 2008-August but really this September, I hope) I am trying to psyche myself up for another round. This time I am going to remove the said limitation, and I'd like to actually discover new baroque numbers, not known before.

First let me present here a very naive approach, just for the illustration (I would never actually code it; you may program this way only if you enjoy programming for its own sake). Select a natural number topNum, say topNum := 10000 and an array of primes, say:

P := (2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79}

Your space of "vertices" V consists now of natural numbers n, 1 < 28 =" 2²⋅7 is one of them. You want to find baroque numbers, which belong to V. For instance, 28 is one of them; indeed, the sum of divisors of 28 is:

sd(28) = 1+2+4+7+14+28 = 2⋅28

hence brq(28) = 2, which means that 28 is a perfect number.

Now define the penalty function pen(n) as the number of primes p ε P such that there exists natural exponent e for which p^e is a divisor of n but not of sd(n). Observe thatg n is baroque if and only if pen(n) = 0. Thus now we may say that our goal is to find many (preferably all) n ε V for which pen(n) = 0.

The algorithm may start with vertex v(0) := 43 (or any other prime from P). Let's assume that the algorithm has already reached vertex v(n). Selecting the next vertex, v(n+1), involves two stages: we select a candidate, then the candidate is accepted with probability prob(n), or rejected. If it is rejected, then we select another candidate, which is going to be accepted as v(n+1), with the same probability prob(n), or rejected, etc., until certain v(n+1) gets finally accepted.

In order to select v(n+1), the algorithm selects randomly a prime p ε P. The product v(n)⋅p is the first candidate for v(n+1). Then it checks condition v(n)⋅p ≤ topNum. If it is satisfied, then it checks penalty: if pen(v(n)⋅p) < pen(n); then v(n+1) := v(n)⋅p is accepted as the next vertex. If the penalty didn't decrease then we still allow a prob(n) := 1/(1 + 10⋅log(1+n)) chance that v(n)⋅p is accepted. If it is rejected, for one reason or another, but p is a divisor of n then v(n+1) := v(n)/n is tried, i.e. the algorithm checks pen(v(n+1)). If v(n+1) := v(n)/n t is accepted then the task of finding the next vertex is accomplished (and algorithm will look for v(n+2)). Otherwise a new prime p ε P is randomly selected and the described candidate process is repeated, etc, until finally a new v(n+1) gets accepted.

The describes algorithm can be discussed, modified, refined but my goal was just to provide an idea how a simulated annealing might work for baroque numbers. Among the shortcomings of the described version is the necessity of representing v(n) actually and directly (naively) as an integer. This either limits the program to a small range of numbers or forces us to use a multi-precision software library--indeed, baroque numbers ten to be huge. In the next posts I'll describe a bit less naive approach.

Remark Above, I have proposed

prob(n) := 1/(1 + 10⋅log(1+n))
One may also try one of the many other possibilities, e.g.

prob(n) := 10/(10 + n^1/4)

Simulated annealing algorithm

The goal is to obtain optimal or nearly optimal solutions for a given problem. Thus first of all we need to define a finite set V of vertices, which correspond to the solutions of the problem, and a non-negative real function, called penalty,

pen : V → R,
for which we would like to find a vertex v ε V such that the value pen(v) is minimal or nearly minimal, so that such optimal v will serve as a solution. This abstract statement sounds simple but in applications the vertices may represent complex configurations or situations.

Some pairs of the vertices are connected by directed edges, so that together, the vertices and the edges, they form a directed graf G := (V E), where E is the set of all edges. Edges are just one-way roads which allow us to travel directly from one vertex to another. An edge leading from vertex v to vertex w is also called (in simulated annealing) a move from v to w. Simulated annealing is applied only to connected graphs, meaning that one can reach any vertex from any other via a finite sequence of consecutive edges (moves).

It is up to the designer of the algorithm to define the graph G = (V E) and the penalty function pen : V → R. One may do it well or one may do it poorly. It's an art. The same goes for the selection of other elements of the algorithm.

In addition, we consider also a real parameter T ≥ 0, called temperature, or rather a decreasing sequence of temperatures T(0) ≥ T(1) ≥ ... ≥ 0, which approach (or even attains) the limit value 0.

The algorithm, i.e. the search for the optimal vertex, works as follows: we start a travel over graph G at a high initial temperature T(0) and in an initial vertex v(0). Let's assume that we are already in vertex v(n) at temperature T(n). Now algorithm selects randomly a candidate for a move from vertex v(n). Let v(n+1) be the destination of the candidate move. If the penalty has decreased, pen(v(n+1)) then v(n+1) is accepted. If not then it is accepted with the probability prob(n) which approaches zero together with the temperature. If the candidate is rejected then another candidate v(n+1) is selected, and the same acceptance procedure applies to it; and so on, until a candidate is accepted.

Remark 1 In addition to simulated annealing algorithm there are also other optimizing algorithms which serve similar purpose. One of them is the so-called genetic algorithm, which is like simulated annealing but richer: besides moves it has also procreation, where one gets a new vertex out of a pair of two previous ones. This provides more possibilities (while it is harder now to provide a scientific foundation as solid as in the case of simulated annealing).

Remark 2 The algorithm may look just for one solution. Then after finding it, it may stop. Or it may look for many solutions. Then after finding a consecutive solution it restarts itself.