A methodology for comparing algorithms and what else it can come in handy for

Having read recently the article “Introduction to Optimization. Annealing simulation ” wanted to participate in a comparison of optimization algorithms. But they really would be nice to compare. And in the materials of the original article no quantitative data are given. So, I thought, we must first formulate the criteria for comparison. What I propose to do in this article.



So, by what criteria can you compare algorithms similar to "simulated annealing"? Most often they compare:
• efficiency
• speed of operation
• memory usage
And I would like to add
• repeatability (That is, check how stable the algorithm is in the results).

Criteria for evaluation

As the first criterion, we will use the length of the resulting route. I have no idea what the best route is on a randomly generated set of points. What to do? Let’s try, for simplicity, to run everything on one set. More pedantic ones, for greater reliability of the results, can then be expanded to several sets and average the results. In the meantime, you need to add the reading of the test set from the file to be able to repeat the experiment.
To determine the speed of operation, we calculate the time spent by the processor on the algorithm, since there is such a function in Octave. Of course, on different computers the result will also be different. Therefore, for comparison, all results can be scaled by any one algorithm that is always run. For example, according to the algorithm from the original article, since he got the first.
On the outlay of memory immediately came up with nothing. Therefore, we will try to figure out on the fingers. As I understand it, the main expenses of the author are two arrays of indices. One for the current route, the other for the candidate route. And, of course, an array of route points along two coordinates per point. We will not take into account the remaining costs, since they will be different in different runtime environments and do not depend on the size of the test data. (Perhaps here I am somewhat harsh, but quite open to suggestions). Total main expenses are N * 4, where N is the number of points to be circumvented.

A small digression: everything described here was tested not in MatLab, as in the original article, but in Octave [http://www.gnu.org/software/octave/]. There are several reasons for this - firstly, I have it, and secondly, its programming language is almost completely compatible with MatLab, and finally, everyone can have it because it is free and public.

Analyzing ...

To get started, make sure compatibility for the existing code (download the link at the end of the original article). Everything works, but it turns out that there is no code that generates beautiful graphics in it. It doesn’t matter - we write our function from which we call the original algorithm and draw the result. All the same, it will be needed to kit the whole collection of statistics.

So the Analyze function:
function [timeUsed, distTotal] = Analyze( )
	% считывание массива городов из файла для повторяемости результата
	load points_data.dat cities
	% засечка времени выполнения 
	timeStart = cputime();
	% собственно поиск маршрута
	result = SimulatedAnnealing(cities, 10, 0.00001);
	% определение затраченного процессором времени
	timeUsed = cputime - timeStart;
	% вычисление общей длины маршрута
	distTotal = PathLengthMat(result);
	% рисование итогового маршрута	
	plot(result(:,1), result(:,2), '-*');
	drawnow ();
end


I did not write a function for generating a file with an array of city coordinates. You can do this directly from the command line:
	cities = rand(100, 2)*10;
	save points_data.dat cities;


Yes, if someone noticed that the CalculateEnergy function from the original algorithm is not used for the assessment, then there are two reasons for this. Firstly, there is an idea to try work on other metrics, but it is better to evaluate the result nevertheless by Euclidean distance. And secondly, although the author of the original article described in detail the calculation of the Euclidean distance, the code in the sources computes Manhattan. Euclidean for two points correctly calculated as follows:
function [dist] = distEuclid(A, B)
    dist = (A - B).^2;
    dist = sum(dist);
    dist = sqrt(dist);
end


Well, to determine how stable the result is from time to time, we will run all this in a cycle many times, collect the results and evaluate their average and the degree of spread of the results:
function [] = Statistics
	testNum = 1000;
	% инициализация массивов для сбора статистики
	timeUsed = zeros(testNum,1);
	distTotal = zeros(testNum,1);
	for n=1:testNum
		[timeUsed(n), distTotal(n)] = Analyze;
	end
	% определение средних значений и разброса
	time = [mean(timeUsed), std(timeUsed)]
	dist = [mean(distTotal), std(distTotal)]
end


The first disappointment is how long it takes! On a computer core i5-3450 3.1GHz 8Gb RAM, one pass is an average of 578 seconds. So I wanted to start and average 1000 cycles for greater reliability, but this is ... 6.69 days. I had to run for 100 cycles and after about 16 hours 4 minutes (that is, in the morning) we have the result: the
average calculation time is 578.1926 (standard deviation in the sample is 3.056), the average length of the resulting path is 91.0844 (standard deviation 2.49). That is, everything is more than stable, and, therefore, suitable for use. But how long. I could not resist and was tempted to do optimization.

Optimization

The first candidate for acceleration is distance estimation. Octave, like MatLab, is optimized for matrix computing - we’ll rewrite this as well. The function will take two arrays: A - the beginning of the segments, B - their ends, respectively:
function [dist] = distEuclid(A, B)
    dist = (A - B).^2;
    dist = sum(dist, 2);
    dist = sqrt(dist);
    dist = sum(dist);
end


To call this function, we introduce another one that will prepare the data by cyclically shifting the initial array of points to obtain an array of segment ends:
function [length] = PathLengthMat(cities)
	citiesShifted = circshift(cities, 1);
	length = distEuclid(cities, citiesShifted);
end


Well, of course, you need to rewrite the function SimulatedAnnealing and GenerateStateCandidate to work directly with point arrays. I will not give the full code, at the end there will be a link to the archive with all the source codes for those who are interested.
We start, look ... the average calculation time is 70.554 (standard deviation in the sample is 0.096), the average length of the resulting path is 90.403 (standard deviation 2.204).
The increase in speed is more than eight times! It’s no longer interesting to try further - you won’t get such an increase. True, memory consumption (according to the proposed methodology) now has two arrays of points in two coordinates (current solution and new candidate), that is ... N * 4. Nothing changed.
But, since we are evaluating the algorithm, we will check the temperature data. The initial “energy” of the generated point arrays ranges from 500 to 700, and the final output routes fluctuate around 90 plus / minus 5. Therefore, the temperatures chosen by the author practically provide only a limit of 100,000 iterations (by the way, it is also hard-coded in the code. To prevent looping, I’ve raised it to ten million, why bother myself when checking the algorithm?). Indeed, having experimented a little, we get approximately the same results for Tmax = 300 and Tmin = 0.001. At the same time, the execution time was reduced to 21 seconds.

Exploring options

Now try to compare with something. Initially, there was an idea to compare with the Manhattan metric, but it turned out that it was implemented. Then there is this optimization. In short - add the code for the matrix calculation of the Manhattan distance and try for it. It is also interesting how another temperature change function will show itself. Those who are deceived by the linearity of the temperature change function promised in the original article are forced to disappoint - it is not. To understand this, just look at its graph:


Probably, the linearity of the denominator of the function was meant, because the denominators of all the other options used in science really look worse.
For comparison, we will also try the temperature function T = Tmax * (A ^ k) where A is chosen, usually within 0.7 - 0.999. (When implemented, it is not greedy at all - Ti = A * Ti-1). Options for "annealing" with its use are called "quenching". It drops much faster to the minimum temperature, but more gentle at the initial stage. That is, the algorithm experiments with it at the beginning more, but finishes work faster. So - a summary table of results:
AlgorithmLength (Spread)Time (scatter)Memory
From the article91.0844 (2.4900)578.1926 (3.0560)N * 4
Dot Matrix (10-0.00001)90.4037 (2.2038)70.554593 (0.096166)N * 4
Dot Matrix (300-0.001)90.7956 (2.1967)20.94024 (0.16235)N * 4
Matrix Manhattan (10-0.00001)90.0506 (3.2397)70.96909 (0.7807)N * 4
Matrix "quenching" (300-0.001)92.0963 (2.3129)22.59591 (0.39581)N * 4


Conclusions.

  1. The ability to quantify the results of the algorithm is useful even if there is no one to compare it with.
  2. If you take into account the features of the tool when writing a program - you can save a lot of time.
  3. If you select the parameters for a specific task - you can save a little more.
  4. It is pointless to compare algorithms by speed without optimization - the execution time will be almost a random variable.
  5. For modern processors, raising and extracting a root is almost as simple as taking a module - you can not distort and write math with almost no simplifications.
  6. Many variations of one algorithm give very close results, if not a graphomaniac - use what is more convenient (or more familiar).


The promised sources can be found on the link.
Well, if you liked it, you can compare it with something else after some time.

Also popular now: