Particle Swarm Optimisation

Having decided that I'm going to use my mfe_mae indicator as a target for my neural net, over the last couple of months I've been doing some research on what would make good features for this target. In the course of this I've also decided that Particle swarm optimization would be a useful tool to use.

Thanks to the pseudo-code on this wiki and the code on in this stackoverflow thread I've been able to write some Octave code to perform pso over one dimension, which is shown in the code box below:

clear all

function [out] = evaluate_fitness( R , val )
    out = R .* R .- val .* R ;
end

val = input( 'Enter test val: ' ) ;

% Initialize the particle positions and their velocities
n_particles = 100 ;
upper_limit = val ;
n_iterations = 50 ;

w = 0.9 .* ones( n_particles , 1 ) ; % is an inertial constant. Good values are usually slightly less than 1. Or it could be randomly initialized for each particle.

% c1 and c2 are constants that say how much the particle is directed towards good positions. They represent a "cognitive" and a "social" component, respectively, 
% in that they affect how much the particle's personal best and the global best (respectively) influence its movement. Usually we take c_1, c_2  approx = 2. 
% Or they could be randomly initialized for each particle.
C1 = 2.0 ;
C2 = 2.0 ;

X = upper_limit * rand( n_particles , 1 ) ; % == particle position vector containing lambda values
V = zeros( size(X) ) ; % particle velocity vector
X_lbest = X ; % == particle position vector containing lambda values for local best
 
% Initialize the global and local fitness to the worst possible. Fitness == LOOCV "press" statistic
fitness_gbest = Inf ; % _gbest == global best
fitness_lbest = fitness_gbest * ones( n_particles , 1 ) ; % _lbest == local best
 
% Loop until convergence, in this example a finite number of iterations chosen
for ii = 1 : n_iterations

    % evaluate the fitness of each particle, i.e. do the linear regression and get
    % the LOOCV statistic
    fitness_X = evaluate_fitness( X , val ) ;
 
    % Update the local bests and their fitness 
    ix = find( fitness_X < fitness_lbest ) ; % if local LOOCV "press" statistic improves
    fitness_lbest( ix ) = fitness_X( ix ) ;  % record this better LOOCV statistic value 
    X_lbest( ix ) = X( ix ) ;                % and the lambda value that gave rise to it    
 
    % Update the global best and its fitness 
    [ min_fitness min_fitness_index ] = min( fitness_X ) ;
    
    if ( min_fitness < fitness_gbest )       % if global LOOCV "press" statistic improves
        fitness_gbest = min_fitness ;        % record this better LOOCV statistic value
        X_gbest = X( min_fitness_index ) ;   % and the lambda value that gave rise to it
    end % end if    
 
    % Update the particle velocity and position vectors
    R1 = rand( n_particles , 1 ) ;
    R2 = rand( n_particles , 1 ) ;
    V = w .* V + C1 .* R1 .* ( X_lbest .- X ) .+ C2 .* R2 .* ( X_gbest .* ones( n_particles , 1 ) .- X ) ;
    X = X .+ V ;
    
end % end main ii loop

and which is vectorised as much as I can make it at the moment. The evaluate_fitness function I've chosen to use is a Quadratic function of the form 

f(x) = x^2 - bx 

which, when a positive value for "val" is input as the test value, ensures the function curve goes through the origin and has a minimum y-axis value at a point on the x-axis that is half the input test value. This make it easy to quickly verify that the pso code is performing as expected, with the global minimum x_axis value found by the algorithm being given by the variable X_gbest. My reasons for choosing a test function of this form, and for looking at pso in general, will be given in my next post.