Christoforos Charalambous
email: c.charalambous@frederick.ac.cy

Krzysztof Fleszar

Khalil S Hindi
email: khalil.hindi@aub.edu.lb

A Hybrid Searching Method for the Unrelated Parallel Machine Scheduling Problem

Keywords: parallel machine scheduling, hybrid optimization, variable neighborhood search, mixed-integer programming

The work addresses the NP-hard problem of scheduling a set of jobs to unrelated parallel machines with the overall objective of minimizing makespan. The solution presented proposes a greedy constructive algorithm followed by an application of a Variable Neighborhood Decent strategy that continually improves the incumbent solution until a local optimum is reached. The strength of the approach lies in the adoption of different objectives at various stages of the search to avoid early local optimum entrapment and, mainly, in the hybridization of heuristic methods and mathematical programming for the definition and exploration of neighborhood structures. Experimental results on a large set of benchmark problems attest to the efficacy of the proposed approach.

Introduction

The unrelated parallel machines scheduling problem is that of scheduling without preemption n jobs available at time zero on m unrelated machines to minimize the makespan (maximal machine completion time), C max . Henceforth, we will also use the term 'span' to signify the completion time of a machine.

The machines are unrelated in that the processing time of a job depends on the machine to which it is assigned. The practical importance of the problem stems from the fact that it is common in many industrial applications to have parallel resources with different capabilities, perhaps procured at different times. These resources would be capable of carrying out the same tasks, but the time taken to perform a task would depend on the resource on which it is performed. Such applications can be found in industries like painting, plastic, textile, glass, semiconductor, chemical, and paper manufacturing, as well as in some service industries [START_REF] Jp Arnaout | A two-stage ant colony optimisation algorithm to monomize the makespan of unrelated parallel machines with sequencedependent setup times[END_REF].

Literature on parallel machine scheduling problems is extensive; for general references and surveys, see [START_REF] Blazewicz | Scheduling computer and manufacturing processes[END_REF][START_REF] Brucker | Scheduling algorithms[END_REF][START_REF] Lee | Optimization and heuristics of scheduling[END_REF][START_REF] Mokotoff | Parallel machine scheduling problems: A survey[END_REF][START_REF] Pinedo | Scheduling: Theory, Algorithms, and Systems[END_REF]. However, most of the literature addresses identical machines, where the processing time of a job is the same regardless of the machine to which it is assigned, or uniform machines, where the processing time of a job is proportional to the speed of the machine. Of the relatively small number of works that address unrelated machines, most deal with the case without setup times. For minimization of the makespan, Martello et al offered exact and approximation algorithms [START_REF] Martello | Exact and approximation algorithms for makespan minimization on unrelated parallel machines[END_REF]; Mokotoff and Chretienne a cutting plane algorithm [START_REF] Mokotoff | A cutting plane algorithm for the unrelated parallel machine scheduling problem[END_REF]; Mokotoff and Jimeno [START_REF] Mokotoff | Heuristics based on partial enumeration for the unrelated parallel processor scheduling problem[END_REF] heuristics based on partial enumeration; Ghirardi and Potts [START_REF] Ghirardi | Makespan minimization for scheduling unrelated parallel machines: A recovering beam search approach[END_REF] a recovering beam search method; and Shchepin and Vakhania [START_REF] Shchepin | An absolute approximation algorithm for scheduling unrelated machines[END_REF][START_REF] Shchepin | An optimal rounding gives a better approximation for scheduling unrelated machines[END_REF] approximation algorithms.

Throughout, the following notation is used: N = {1, . . . , n} set of jobs, n being the total number of jobs. M = {1, . . . , m} set of machines, m being the total number of machines. p jk processing time of job j on machine k M (j) the machine to which job j is assigned to. J k set of jobs assigned to machine k. C k span (completion time) of machine k. C max = max k∈M C k makespan (maximum machine completion time).

Definition of objective criteria

The problem's main objective is to minimize the solution makespan (C max ). However, focusing solely on the main objective may inhibit the development of the search process as other desirable characteristics are likely to be ignored. In this work, in addition to the main objective, auxiliary criteria are employed to assist the development of high quality solutions. The criteria used are:

1. Sum of spans: To be able to minimize the makespan it is helpful to associate jobs to machines that process them efficiently. The lower the sum of spans (an average measurement on how efficiently the jobs have been assigned), the more slack can be expected on non-makespan inducing machines, leading to an increased capacity for solution improvement through job reassignments. 2. Number of makespan-inducing machines Neighborhood search structures operate by altering incumbent solutions through a small number of job reassignments aiming at an overall improvement. The fewer the number of machines that have a span equal to the makespan, the more probable the existence of a job reallocation pattern that would reduce the makespan.

The mode in which objective criteria are used depends on the neighborhood structure, as described in Section 4.

Initial solution

Numerous constructive algorithms have been proposed in the literature including the efficient APPROX method of Martello [START_REF] Martello | Exact and approximation algorithms for makespan minimization on unrelated parallel machines[END_REF]. These algorithms focus on finding techniques for developing solutions that minimize the objective function. In this work the primal goal of the solution construction is not to minimize makespan but rather to seed the searching phase of the algorithm with a solution that resides in a promising region. The rational against opting for makespan minimization is that fast constructive algorithms invariably adopt a greedy approach in the solution construction that often opts for short-sighted moves which worsen the objective function least. Such an approach leads to relatively good solutions but, due to the naive nature of choices, there is little room left for further improvements.

The initial solution is calculated based on the linear relaxation of the UPMSP min C max [START_REF] Jp Arnaout | A two-stage ant colony optimisation algorithm to monomize the makespan of unrelated parallel machines with sequencedependent setup times[END_REF] subject to

k∈M x jk = 1 ∀j ∈ N (2)
C max ≥ j∈N p jk x jk ∀k ∈ M (3)
x jk ∈ {0, 1} ∀j ∈ N, k ∈ M (4)
where x jk is a binary variable signifying whether job j is assigned to machine k.

When the integrality constraint in 4 is relaxed, the solution generated by the solver will consist of a large set of x jk variables set to zero, and the rest being either one or set to a fractional value. The initialization heuristic sets the upper bound of all x variables that have value zero to zero (i.e. forbidding the corresponding assignment). Subsequently, all fractional variables are sorted based on regret. The regret of a variable is defined as the added processing time that would incur in the solution if the corresponding variable was set to one and is calculated by p jk (1-x jk ). The variable with the maximum regret has its upper bound set to zero and the model is resolved until no variable has a fractional value. Note that if the model has no fractional variables, the corresponding solution is feasible.

Neighborhood Moves and Variable Neighborhood Decent

Improving an incumbent solution can often be achieved through the reassignment of jobs to machines. The means of these reassignments, henceforth called moves, is central to the definition of the neighborhood structures used in the searching approach. All moves are based on a series of additions (allocating a job to a machine) and removals (deallocating a job from a machine) and can be categorized as follows:

Transfer removes job j from its machine, M (j), and adds it to machine k = M (j) (one removal and one addition). Swap removes jobs i and j (i = j) from their respective machines (M (i) = M (j)), and adds i to M (j) and j to M (i) (two removals and two additions).

Closed ejection chain involves a closed chain of h ≥ 3 machines. A job is removed from each machine of the chain and added to the next machine in the chain (h removals and h additions). Open ejection chain involves an open chain of h ≥ 3 machines. A job is removed from each but the last machine of the chain and added to the next machine in the chain (h -1 removals and h -1 additions). Compound move is any collection of transfers, swaps, closed, and open ejection chains such that on each machine at most one addition (restriction 1) and at most one removal (restriction 2) is performed (h removals and h additions, 1 ≤ h ≤ m). Furthermore, an expanded compound move can be effected if restriction 1 or restriction 2 is removed, but not both.

The valid moves of a neighborhood define its structure and size. The simpler the moves that can be employed the faster the exploration of the neighborhood. On the other hand, simpler structures are less likely to be capable of identifying improving moves than more involved ones.

The proposed algorithm defines three neighborhood structures of increasing complexity named small (SNS), medium (MNS) and large (LNS). SNS is applied first and upon local optimum entrapment the search progresses to the MNS, altering also the objective criterion. When MNS also reaches a local optimum, LNS is employed. If one iteration of LNS leads to an improvement, the cycle is repeated otherwise the incumbent solution is the algorithm's proposed solution.

Small Neighbourhood Search

The initial solution is often quite poor but with high potential, given the high slacks. Improvements are obtained through a small neighborhood structure that is searched exhaustively, adopting the best improving move each time. In SNS, the neighborhood structure is defined as any transfer or swap that involves the machine inducing the makespan (in case more that one machines have spans equal to the makespan, one is randomly designated as the makespan machine).

A move is considered improving if, for all the machines affected by it, the new span is smaller than the incumbent makespan. Given the set of improving moves of a small neighborhood, the best is considered the one that minimizes the sum of spans (see Section 2). To further assist the guidance of the search to promising regions, in cases where no improving move can be found through transfers or swaps, or where the move increases the total span to more than 10% of the lower bound, the neighborhood structure is augmented to include ejection chains (open or close) of size 3. SNS terminates when no improving move can be found.

The adoption of a dynamic neighborhood structure proved necessary to balance the effectiveness and efficiency of the search. Exhaustively examining transfer and swap moves can be done very fast, given the simplicity of the moves. If a satisfactory move is identified, something likely in the early stages, it is adopted, thus allowing a speedy convergence without sacrificing the solution characteristics. Only upon non-satisfactory moves is the, more involved, ejection chain move examined.

Medium Neighbourhood Search

Focusing only on moves that involve the makespan machine inhibits the search from identifying more promising solutions. Clearly, moves in which the makespan machine is not involved cannot improve the overall objective of makespan minimization. However, moves involving the remaining machines can lead to improvements by increasing slack space (reducing the sum of spans) or reducing the number of makespan machines, both of which are desirable as explained in Section 2. In MNS, valid moves include any transfer or swap between any pair of machines. A move is considered undesirable if it worsens any of the two auxiliary objectives. MNS adopts a first-improving approach, committing the first move identified that reduces the sum of spans without increasing the number of makespan machines. Using a first-improving instead of a best-improving approach and limiting the moves to transfers and swaps was preferred to economize on computational effort, given the increase in the set of machines that could be involved in MNS moves.

Large Neighbourhood Search

As explained, neighborhood structures that employ limited scope moves have limited capabilities of yielding improvements. This is particularly true when the incumbent solution is relatively tight, as solutions obtained after SNS and MNS are likely to be. To escape from such solutions, a large neighborhood structure is defined that explores all feasible compound moves.

Such structures can identify complicated moves that may yield improvements but the difficulty lies in the development of efficient mechanisms for facilitating the search. To achieve this a MILP model is developed as described below.

min C max + k∈M C k 1000|M | (5) subject to j∈N y jk ≤ 1 ∀k ∈ M (6) j∈J k l∈M y jl ≤ 1 ∀k ∈ M (7)
C k = C k + j∈J p jk y jk - j∈J k l∈M p jk y jl ∀k ∈ M (8) C max ≥ C k ∀k ∈ M (9) C max ≤ C max -1 ( 10
)
where

y jk -boolean variable set to 1 if job j is moved to machine k C k -incumbent span of machine k C max -incumbent makespan
The objective function (equation 5) minimizes the makespan using also the auxiliary objective of sum of spans minimization as a secondary objective. Constraint 6 ensures that at most one job can be added to any machine whereas constraint 7 limits the number of jobs that can be removed from a machine to one. Constraint 8 calculates the new span of each machine by adding to the incumbent span the processing times of all jobs moved to that machine and subtracting the processing times of all jobs that were allocated to the machine and have been moved to another machine. Constraint 9 links the makespan to the new spans and constraint 10 ensures that the makespan is improved.

LNS commences by restricting the compound moves to those that allow at most one insertion and one removal on each machine. If no such move is found, the constraint 6 is relaxed, allowing for more than one job to be added to a machine. If again no improving compound move is identified, constraint 6 is reinstated and constraint 7 is relaxed, allowing for more than one jobs to be removed from a machine but limiting the jobs to be added to a machine to one.

Given the large number of binary variables defined by the model, it is crucial to reduce as many as possible so that the model can be solved by a MP solver in realistic execution times. The following reductions are defined:

-y jk = 0 ∀j ∈ J k , k ∈ M -Assuming p max k represents the largest processing time on machine k of the jobs in J k , then y jk = 0 ∀j ∈ N if C k + p jk -p max k > C max -1 -Given two jobs i, j ∈ J k , y j,l = 0 ∀l ∈ M \{k} if p ik ≥ p jk ∧ p il ≤ p jl
The first reduction ensures that a job cannot be moved to its current allocation. The second reduction reduces a move if even assuming the removal of the job with the largest processing time from the destination machine, the move would still be bigger than or equal to the incumbent makespan. Finally, the third reduction states that if two jobs are currently assigned to the same machine, then if the first has bigger or equal processing time on the incumbent machine and for a given destination machine it has less or equal processing time, such a move would always be preferable therefore the variable corresponding to moving the second is reduced. Note that if constraint 7 is relaxed, the second and third reductions do not apply.

Since the incumbent solution is likely to be tight given the prior application of SNS and MNS, the amount of variables reduced is high and experimentation has shown that relevant models can be addressed within acceptable execution times. As a safety net, the time assigned to the MP solver for a model instance has been capped to 1 second.

Computational results

The proposed algorithm was implemented in C#.NET. The MIP problems were solved using ILOG CPLEX version 12.1, with the number of threads limited to 1 and the relative and absolute gaps set to zero. All tests were executed on a x5570 2.93GHz Intel Xeon processor using a single thread on a single core.

Three classes of random problem instances were generated, adopting the same pattern as that used in the literature. Specifically, for the first two classes processing times were taken randomly from uniform distributions U (10, 100) and U (10, 1000), respectively. For the third, speeds are assigned to machines and lengths to jobs and a weighted randomized processing time is calculated for each job-machine pair (see. [START_REF] Ghirardi | Makespan minimization for scheduling unrelated parallel machines: A recovering beam search approach[END_REF] for more details).

For each problem class, 20 different instances were generated with the machine number obtained from {10, 20, 30, 40, 50} and the job number from {100, 200, 500, 1000}. Computational results aggregated over classes, machines, and jobs are provided in Table 1. All presented results refer to the average over the 20 instances. The proposed algorithm was compared against the APPROX algorithm of Martello et al [START_REF] Martello | Exact and approximation algorithms for makespan minimization on unrelated parallel machines[END_REF] and the Recovering Beam Search method with beam width 3 of Ghirardi and Potts [START_REF] Ghirardi | Makespan minimization for scheduling unrelated parallel machines: A recovering beam search approach[END_REF] (denoted by RBS). For each algorithm the deviation from the lower bound (obtained through applying Langrangean relaxation on the UPMSP model) and the average and maximum execution time are provided. Finally, the last table column (denoted as 'Imp') shows the percentage of improvement (or deterioration) of the proposed solution against the best known deviation. For APPROX and RBS the deviation and times presented in [START_REF] Ghirardi | Makespan minimization for scheduling unrelated parallel machines: A recovering beam search approach[END_REF] were used (implemented in C++ and executed on a PC-Pentium III/866). As shown, the proposed VND greatly outperforms the other algorithms when compared by either class, machines or jobs. VND performs best in the third class which also has the highest time demands. Comparing by number of machines, it seems that VND consistently outperforms the other algorithms independent of machine size. On the contrary, when comparing by job number, it is evident that the performance improvement is linked with the problem size. Also, taking into account processor performance benchmarks, VND achieves comparable execution times with RBS.

Conclusion

A variable-neighborhood descent, hybridized with mathematical programming elements, that addresses the NP-hard problem of scheduling unrelated parallel machines with the objective of minimizing the makespan, has been presented. Neighborhood search structures, of increasing complexity, are investigated to guide the search to promising regions. Local optimum entrapment is avoided through the use of large neighborhoods which are explored with the use of mathematical programming. The performance of the proposed scheme has been evaluated by comparing its solutions on a set of benchmark problem instances to those generated by the best available algorithms. The results show that the proposed scheme provides considerably better solutions.

Table 1 .

 1 Comparative results on set of problem instances

	Aggr. by	APPROX			RBS			VND
			Dev(%) avT(s) maxT(s)	Dev(%) avT(s) maxT(s)	Dev(%) avT(s) maxT(s)	Imp(%)
	cl.	1	6.1	17.7	83.0	5.1	79.6	392.2	4.0	1.3	3.8	31.7
		2	12.0	13.5	61.9	10.7	74.0	351.5	10.1	1.7	5.4	25.2
		3	3.1	54.6	212.4	2.1	92.7	427.8	1.4	2.5	7.8	51.2
	m	10	1.9	22.8	174.4	1.2	63.8	262.1	0.9	0.8	2.1	41.5
		20	4.0	24.4	176.1	3.0	65.3	260.4	2.3	1.4	3.2	40.5
		30	6.5	28.0	181.8	5.2	75.1	291.2	4.7	2.0	5.2	33.6
		40	9.7	33.0	201.7	8.3	98.9	391.7	7.0	2.3	6.6	35.5
		50	13.2	34.8	212.4	12.1 107.4	427.8	11.0	2.6	7.8	29.0
	n 100	17.2	< 1	< 1	15.7	< 1	< 1	15.3	0.4	1.1	7.8
		200	6.6	1.4	3.1	5.2	1.8	3.7	4.6	1.0	2.4	22.2
		500	2.7	15.5	31.3	1.9	29.3	41.9	1.0	2.1	4.1	52.8
		1000	1.7	97.0	212.4	1.1 285.7	392.2	0.4	3.8	7.8	61.2
	overall	7.1 28.6	212.4	6.0 82.1	427.8	5.2	1.8	7.8	36.0


Jacek Kabziński
email: jacek.kabzinski@plodz.pl

One-Dimensional Linear Local Prototypes for Effective Selection of Neuro-Fuzzy Sugeno Model Initial Structure

Keywords:

We consider a Takagi-Sugeno-Kang (TSK) fuzzy rule based system used to model a memory-less nonlinearity from numerical data. We develop a simple and effective technique allowing to remove irrelevant inputs, choose a number of membership functions for each input, propose well estimated starting values of membership functions and consequent parameters. All this will make the fuzzy model more concise and transparent. The final training procedure will be shorter and more effective.

Introduction

We consider a Takagi-Sugeno-Kang (TSK) fuzzy rule based system used to model a memory-less nonlinearity from numerical data. It is typical for such application to perform a training of the system to improve modeling accuracy. One of the first (and very successful) approaches to the problem was Adaptive Neuro-Fuzzy Inference System (ANFIS) proposed by Jang [START_REF] Jang | ANFIS: Adaptive-network-based fuzzy inference system[END_REF]. The idea was to realize a fuzzy inference system as a neural network and to apply neural network learning techniques to train the fuzzy model. The ANFIS technique (where least squares optimization of consequent parameters is used and gradient decent error back propagation is applied to tune membership functions parameters) is offered as a standard tool of Matlab and is still very popular and employed to solve various modeling problems -for example [START_REF] Buyukbingola | Adaptive neuro-fuzzy inference system (ANFIS): A new approach to predictive modeling in QSAR applications: A study of neuro-fuzzy modeling of PCP-based NMDA receptor antagonists[END_REF][START_REF] Kovac | The ANFIS-based route preference estimation in sea navigation[END_REF]. Several other learning algorithms were developed to train fuzzy models including genetic, bacterial and PSO algorithms [START_REF] Alcalá | Local Identification of Prototypes for Genetic Learning of Accurate TSK Fuzzy Rule-Based Systems[END_REF][START_REF] Lin | The design of neuro-fuzzy networks using particle swarm optimization and recursive singular value decomposition[END_REF][START_REF] Gal | Modified bacterial memetic algorithm used for fuzzy rule base extraction[END_REF]. Before the training of parameters is started we have to propose initial structure of the fuzzy model. For a real-world problems we may consider a great number of potential inputs, but it is essential to select only the really important ones for system description. Also the number of membership functions for each input should be carefully chosen to compromise between model accuracy and complexity. The necessity of ANFIS input selection was noticed already by Jang [START_REF] Jang | Input selection for ANFIS learning[END_REF] and this aspect of fuzzy model complexity reduction was investigated by many researchers [START_REF] Hu | Input selection in learning systems: a brief review of some important issues and recent developments[END_REF], but the problem remains open. The most popular approaches to initial system architecture selection are based on several clustering techniques [START_REF] Hu | Input selection in learning systems: a brief review of some important issues and recent developments[END_REF].

The purpose of this work is to develop a simple and effective technique allowing to remove irrelevant inputs, choose a number of membership functions for each input, propose well estimated starting values of membership functions and consequent parameters. All this will make the fuzzy model more concise and transparent. The resulting final training procedure will be shorter and more effective. The proposed approach is based on a concept of fuzzy projection of the data on each input space [START_REF] Lin | Using fuzzy partitions to create fuzzy systems from input-output data and set the initial weights in a fuzzy neural network[END_REF][START_REF] Kabziński | Effective selection of neuro-fuzzy Sugeno model architecture[END_REF], combined with piece-wise linear approximation allowing to choose automatically a number of membership functions according to the required accuracy and to estimate starting values of all fuzzy system parameters. We demonstrate the effectiveness of the proposed technique investigating several examples and comparing results with popular clustering methods.

Neuro-fuzzy inference system

We consider a neural network realisation of TSK fuzzy model [START_REF] Jang | ANFIS: Adaptive-network-based fuzzy inference system[END_REF]. The network possesses N inputs x 1 ,...,x N , one output y, m i membership functions

  i i j i m j i ,...., 1 , ,   associated with the i-th input   N i ,...., 1
 , N m R i i N   
       , , , 1 1 1 1 1 1
where

  i i m j ,...., 1

and r is the rule number. We can choose any smooth membership function, for example generalized bell-shape function:

 i j b x x c a , ( )          1 1 2 (2)
The rule firing strength is calculated as:

w x r i j i N i i     , ( ) 1 (3)
for the r-th rule and the normalised firing strength is:

w w w r r i i N R    1 . ( 4
)
The final output is calculated as:

y w f x x r r N r N R    1 1 ( ,..., ) . (5)
It is well known that this Sugeno fuzzy model may be realised as a five-layer neural network: the first layer calculates the membership functions values, the second layerfiring strengths, the third -normalised firing strengths, the fourth and the fifth calculates the final output [START_REF] Jang | ANFIS: Adaptive-network-based fuzzy inference system[END_REF]. The standard training procedure proposed for ANFIS combines gradient descent back-propagation method to tune membership functions parameters and least squares method to find optimal consequent parameters.

The system is supposed to model a memory-less nonlinearity given by input-output data with possible extraneous inputs. We assume that that the input-output data are given by

    m k N N y x x k k N k ,..., 1 , , , , , 1     . ( 6
)

Initial structure selection procedure

To test the significance of the i-th input,

  N i ,...., 1 
we consider a single-input fuzzy model (SIFM) described below:

 input -x i , output -c i ,  input linguistic categories:   m k x IS x k i i ,..., 1 ,  ,  membership functions:  i k i k b x x x a , , ( )          1 1 2 , ( 7
)
 rules:  . ,..., 1 0 , 0 0 0 , , , , , , , , , , , m k
x if y q p x if q x y p q x p c THEN x IS x IF k i k k i k i k i k i k i k k i k i i k i i k i i          . ( 8
)
The action curve given by the output c i of this system for the input data -x i,k generalises information coded by x y i k k ,  . The degree of this generalisation depends on parameter a . As we want to cover the whole range of the i-th input by a few (say 3 or 5) membership functions, it is reasonable to take such a that the set in which a membership function is 'active', say { :

( ) . } , x x i k   0 3 covers 10-30% of the interval min( ), max( ) , , k i k k i k x x      
. It may be proved that the shape of the action curve of SIFM is robust to single outliers in the measured data and to the measurement noise.

If the i-th input is inessential the curve generated by corresponding SIFM will be flat, if it is meaningful the curve will cover significant part of min( ), max( )

k k k k y y       .
Building SIFM for every input, plotting and testing it's action curve

    m k x c x A k i i k i i , , 1 )) ( , ( , ,
  
we are able to classify the importance of the input and to select N significant inputs with corresponding curves. We must remember that sometimes data symmetry may lead to false diagnosis that the input is irrelevant one.

In this case other techniques for input selection must be used [START_REF] Hu | Input selection in learning systems: a brief review of some important issues and recent developments[END_REF], or the modelling domain should be narrowed, or nonlinear data transformation may be applied.

Selection of membership functions for each input is based on piece-wise linear approximation of action curves derived above. Uniform or mean-square approach are both applicable. The input data have to be sorted in ascending order and repeated values must be removed -it is important especially for grid-type data. It is not necessary to obtain continuous piece-wise linear approximation. Mean square approximation means that starting from the left we are to approximate (by a linear function) maximal number of subsequent points

)) ( , ( , , k i i k i x c x
, under the condition that the mean-square error is limited by a given parameter  . Uniform approximation consists in placing the longest possible line segments ('sticks') inside a  -wide 'tube' surrounding the curve

    m k x c x A k i i k i i , , 1
)) ( , ( , ,    on the plane. Both approaches are easy to implement numerically and produces equivalent results. Of course design parameter  influences the number of intervals and so the number of membership functions. As the result of piece-wise linear approximation for the i-th significant input we obtain m i intervals

    i j i j i j i m j x x I , , 2 , 1 , , max , min ,
   . ( 9
)
and a linear polynomial

    i j i j i i j i i j i i j i m j x x x p x p x P , , 2 , 1 , , ) ( , max , min , 0 , 1 ,
     . ( 10
)
for each interval. The next step is to build a neuro-fuzzy (ANFIS) model to imitate curve

    m k x c x A k i i k i i , , 1 )) ( , ( , ,   
with small number of membership functions. The proposed membership functions for the i-th significant input will be m i generalized bell-shape functions

) ( , x j i  with parameters       i j i j i j i j i j i j i m j x x c x x a , , 2 , 1 , 2 1 , 2 1 , min , max , , min , max ,       . ( 11
) So the function ) ( , x j i  is spanned over j i I
, and centred at the middle point of for all membership functions. The rules for the proposed SISO neuro-fuzzy model will be:

 , ,..., 2 , 1 , , 0 , 1 , i j i i j i i j i i m j p x p c THEN IS x IF     (12)
where starting values of parameters are taken from piecewise linear approximation results and are given by ( 10) and (11). Each rule describes a linear, local fuzzy prototype of fuzzy-filtered one-dimensional data. Single input single output neuro-fuzzy models proposed for each significant input are trained using the data

  m k x c x k i i k i , , 1
j i i j i i j i i p x p c THEN IS x IF , 0 , 1 ,    with optimised parameters   j i j i p p , 0 , 1 , ˆ.
The multi-input model, which will be a starting point for the final training is equipped

with N m R i i N    1 rules:     i i j N j r r N j N j j N N j m j p p N q q x p x p y THEN IS x AND AND IS x IF N N N , , 2 , 1 , ... 1 ... ) ( ... ) ( , 0 , 1 0 , 1 1 , 1 1 , , 1 1 1 1 1
           (13)
and membership functions parameters are 𝑎 𝑖,𝑗 , 𝑏 𝑖,𝑗 , 𝑐 𝑖,𝑗 . The above neuro-fuzzy inference system is trained by a selected procedure to obtain the final model.

Examples

The results of the proposed approach were compared with three standard methods of generation initial structure of fuzzy inference system from numerical data applied in Matlab [11]:  genfis1 -generates a FIS structure from a training data set, using a grid partition on the data (no clustering), number and type of membership functions is given by the user;  genfis2 -generates Fuzzy Inference System structure from data using subtractive clustering, cluster radii is given by the user;  genfis3 -generates Fuzzy Inference System structure from data using FCM clustering, number of clusters is given by the user. User-defined parameters were chosen to obtain the same number of membership functions as generated automatically by the proposed method. For both functions the starting error as well as the whole learning performance was superior in case of the proposed method. The design procedure is illustrated in fig. 12345. The next task was modeling a section of a multidimensional ball with one 'dumb' input. 1000 points were generated on random in the unit multi-dimensional cube. The results are presented in fig. 6. Again the proposed method is better than the standard clustering procedures.

We propose simple and effective procedure for creating neuro-fuzzy TSK models of complex systems from the numerical input-output data. Significant inputs, number of membership functions for each input and initial values of all system parameters are set automatically based on accuracy requirements. Because the initial structure and parameters are set properly, we need a few training iterations for the neural network representation of our model to converge. The generated FIS may be used as a starting point to any training algorithms -also those incorporating reduction of fuzzy rules what is highly recommended, although was not discussed here. The proposed technique was compared with clustering methods implemented in Matlab. Results obtained from the proposed method provided several times smaller starting training error, leading to simpler models, smaller errors and reduction of necessary training time. Several important problems were not discussed here, for example elimination of relevant but dependent inputs.

  choice of the third parameter b i,j is arbitrary -it is reasonable to start with .

  We are to model two 2-dimmensional functions: a) 𝑦 = 𝑥 1 2 + 𝑥 2 2 -cos 18𝑥 1 -cos 18𝑥 2 (Riastragin function), 1][0,1]. We generate 400 points of grid data in [0,1][0,1] and calculate corresponding values of y. We make a false assumption that the third input influences the function -we consider the triples (x 1 , x 2 , x 3 ) where value of x 3 is generated on random from [0,1]. As we see in both cases the curve A 3 . is flat -we can eliminate input x 3 as an inessential one. Piece-wise linear mean square approximation was applied and for 001 . 0   numbers of membership functions was: for problem a) 𝑚 1 , 𝑚 2 = 6,6 , for problem b) 𝑚 1 , 𝑚 2 = 2,3 . Because of parameters generated by the linear, local prototype approach, the initial position of one-dimensional FIS is close to the desired one and the training is fast and effective. Finally the initial FIS was generated according to (13) and compared with those generated by genfis1-3.

Fig. 1 .Fig. 2 .

 12 Fig.1. Single-input fuzzy model curves left: function a), right: function b)

Fig. 3 .

 3 Fig.3. Training of a single-input fuzzy model to the action curve, function b).

Fig. 4 .

 4 Fig.4. Action surface of an initial FIS (before the final training) and numerical data +, function b).

Fig. 5 .

 5 Fig.5. Final model training: left function a), right function b).

Fig. 6 .

 6 Fig.6. Modeling a multidimensional ball: left -single-input fuzzy curves, right -final model training

  with optimised parameters 𝑎 𝑖,𝑗 , 𝑏 𝑖,𝑗 , 𝑐 𝑖,𝑗 ,

	)), parameters (10, 11) are carefully chosen and as we train one-dimensional model, the ( , ( , ,   . As the starting values of training is fast and the results are accurate. Finally, after the training we get N one-dimensional neuro-fuzzy models, the i-th one is equipped with  m i membership functions ) ( ,  j i 
		m i rules


Gerrit K Janssens
email: gerrit.janssens@uhasselt.be

José Maria Pangilinan
email: joey.pangilinan@slu.edu.ph

Multiple criteria performance analysis of nondominated sets obtained by multi-objective evolutionary algorithms for optimisation

Keywords: evolutionary algorithms, multi-objective optimization, multicriteria analysis

The paper shows the importance of a multi-criteria performance analysis in evaluating the quality of non-dominated sets. The sets are generated by the use of evolutionary algorithms, more specifically through SPEA2 or NSGA-II. Problem examples from different problem domains are analyzed on four criteria of quality. These four criteria namely cardinality of the nondominated set, spread of the solutions, hyper-volume, and set coverage do not favour any algorithm along the problem examples. In the Multiple Shortest Path Problem (MSPP) examples, the spread of solutions is the decisive factor for the 2S|1M configuration, and the cardinality and set coverage for the 3S configuration. The differences in set coverage values between SPEA2 and NSGA-II in the MSPP are small since both algorithms have almost identical non-dominated solutions. In the Decision Tree examples, the decisive factors are set coverage and hyper-volume. The computations show that the decisive criterion or criteria vary in all examples except for the set coverage criterion. This shows the importance of a binary measure in evaluating the quality of nondominated sets, as the measure itself tests for dominance. The various criteria are confronted by means of a multi-criteria decision tool.

Introduction

Many real-world optimization problems are multi-objective by nature and with objectives that are in conflict. Mathematical techniques are available to find bestcompromise solutions by aggregating multiple objectives into a single function [START_REF] Ehrgott | Multicriteria Optimization[END_REF]. They have their drawbacks as they have difficulty dealing with concave and discontinuous Pareto fronts. Stochastic local or global search algorithms or population-based algorithms are often used when exact methods are infeasible to be applied and are used to solve difficult optimization problems. Among the many types of these algorithms, Evolutionary Algorithms (EA) seem particularly suitable to solve multi-objective optimization problems. An increasing number of research papers report comparative findings of several evolutionary algorithms in terms of computing speed and Pareto optimality as tested on various multi-objective problem instances or applications with known Pareto sets. In most practical or experimental cases however the Pareto sets are unknown.

The objective of the study is to describe the performance of evolutionary algorithms in terms of stability, computational complexity, diversity and optimality of solutions in different multi-objective optimization problems (MOOP). Multi-objective evolutionary algorithm (MOEA) experiments generate a variety of non-dominated sets in each problem domain. Comparisons of the quality of the non-dominated sets have not yet been presented in each of the cases. The stability of an algorithm is concerned with the sensitivity of the results to changes in the MOEA parameters settings. Computational complexity refers to the solution run-time complexity in terms of the size of the problem. Diversity measures the spread of solutions in the non-dominated set in order to provide the decision maker a true picture of trade-off solutions. Optimality measures the proximity of the best non-dominated set to the Pareto-optimal set.

Bosman and Thierens [START_REF] Bosman | A balance between proximity and diversity in multiobjective evolutionary algorithms[END_REF] argue that the quest for finding the components that result in the best EAs for multi-objective optimization is not likely to converge to a single, specific MOEA. They stated that the trade-off between the goals of proximity and diversity preservation plays an important role in the exploitation and exploration phases of any MOEA.

A comprehensive discussion of multi-objective evolutionary algorithms (MOEA) can be found in [START_REF] Deb | MultiObjective Optimization Using Evolutionary Algorithms[END_REF]. In addition, Coello [START_REF] Coello | An updated survey of GA-based multiobjective optimization techniques[END_REF] gives a summary of current approaches in MOEA and emphasizes the importance of new approaches in exploiting the capabilities of evolutionary algorithms in multi-objective optimization. Zitzler and Thiele [START_REF] Zitzler | Multiobjective evolutionary algorithms: a comparative study and the strength Pareto approach[END_REF] performed a comparative analysis of existing evolutionary algorithms in multi-objective optimization by means of well-defined quantitative performance measures. In this research a comparison is made using two well-known techniques: SPEA2 and NSGA-II.

Zitzler et al. [START_REF] Zitzler | SPEA2: Improving the strength Pareto evolutionary algorithm for multiobjective optimization[END_REF] introduced the Strength Pareto Evolutionary Algorithm 2 (SPEA2), which is an extension and improvement of the original work by Zitzler and Thiele [START_REF] Zitzler | Multiobjective evolutionary algorithms: a comparative study and the strength Pareto approach[END_REF]. SPEA2 integrates a fitness assignment strategy, which considers the number of individuals that an individual dominates and the number of its dominators. It uses a nearest-neighbor density estimation technique that guides the search more efficiently and avoids the formation of new solutions in only a few clusters in the search space. SPEA2 has a truncation procedure that preserves the best solutions when the number of non-dominated individuals exceeds the external population size.

Deb et al. [START_REF] Deb | A fast elitist multiobjective genetic algorithm: NSGA-II[END_REF] introduced an elitist non-dominated genetic algorithm (NSGA-II) that uses not only an elite-preserving strategy but also an explicit-diversity preserving mechanism. Initially, NSGA-II creates a random parent population, sorts the population based on non-domination and assigns each solution a fitness value equal to its non-domination level.

Consider first a case where an optimal algorithm exists to obtain the Pareto set. For example, in the Competitive Facility Location Problem (CFLP) [START_REF] Eiselt | Competitive location models: a framework and bibliography[END_REF], the quality of the non-dominated sets generated by the MOEA can be calculated and compared to the Pareto-optimal set generated by an algorithm by Carrizosa and Plastria [START_REF] Carrizosa | On minquantile and maxcovering optimization[END_REF]. The error ratio metric, which measures the closeness of the non-dominated set to the Pareto front in terms of set membership, can be used to measure the quality of solutions in the CFLP. However, if such an algorithm does not exist or is not available, like for two other problems under study here, the Multi-Objective Shortest Path Problem (MSPP) and the Decision Tree (DT) experiments, only approximations to their Pareto-optimal sets are available for performance analysis.

The next sections present some performance metrics that are useful in measuring the quality of non-dominated sets when the Pareto-optimal set is unknown, and utilizes a multi-criteria tool to determine the best non-dominated set for MSPP, and DT problems based on the performance metrics. Several performance metrics exist in literature, and several comparative studies have been conducted that evaluate them.

The studies presented above show a variety of results and no single MOEA performs better in the different performance metrics but most studies compare their algorithms with either NSGA-II or SPEA2 or both. The studies above mostly evaluate the performance of the selection operators of each MOEA without investigating the effect of the parameter settings on its performance. The current study investigates the performance of NSGA-II and SPEA2 in selected multi-objective optimization problems by means of a multi-criteria method in which various performance metrics of the Pareto front are presented to the decision-maker.

Performance metrics

Deb [START_REF] Deb | MultiObjective Optimization Using Evolutionary Algorithms[END_REF] states that there are two orthogonal goals for any multi-objective optimisation algorithm: (1) to identify solutions as close as possible to the true Pareto-optimal set and (2) to identify a diverse sets of solutions distributed evenly across the entire Pareto-optimal surface. This has led to several metrics that characterise either closeness, or diversity, or both. Examples of metrics which measure the closeness to the Pareto surface are the Error ratio [START_REF] Van Veldhuizen | Multiobjective Evolutionary Algorithms: Classifications, Analyses and New Innovations[END_REF] and the Set coverage [START_REF] Zitzler | Evolutionary Algorithms for Multiobjective Optimization: Methods and Applications[END_REF]. Examples of metrics which measure the diversity across the Pareto surface are the Spacing [START_REF] Schott | Fault tolerant design using single and multicriteria genetic algorithm optimization[END_REF] and the Spread [START_REF] Deb | A fast elitist multiobjective genetic algorithm: NSGA-II[END_REF]. A measure like the Hypervolume measures both closeness and diversity [START_REF] Zitzler | Multiobjective evolutionary algorithms: a comparative study and the strength Pareto approach[END_REF]. Most measures are unary quality measures, i.e. the measure assigns to each Pareto set approximation a number that reflects a certain quality aspect.

Hypervolume

The hypervolume metric [START_REF] Zitzler | Multiobjective evolutionary algorithms: a comparative study and the strength Pareto approach[END_REF] calculates the volume in the objective space covered by the members of the non-dominated set Q. For each solution i ∈ Q, a hypercube v i is computed from a reference point and the solution i as the diagonal corners of the hypercube, The reference point can be found by constructing a vector of worst objective function values. The hypervolume (HV) is calculated as:

( ) i | Q | i v volume HV 1 = = U . ( 1
)
The hypervolume metric is interesting because it is a single metric which is sensitive both to the overall advancement of the non-dominated set and to the distribution of individual points across the set. The placement of the reference point is critical and determines the sense and the magnitude of the hypervolume. Problems may appear if objectives have dissimilar scales or if some objectives are unbounded.

Spacing

Schott [START_REF] Schott | Fault tolerant design using single and multicriteria genetic algorithm optimization[END_REF] introduced a metric, which is a measure of the relative distances between consecutive solutions in the non-dominated set Q is calculated as:

( ) Q ∑ = - = Q i i d d S 1 (2) ∑ = ≠ ∧ ∈ - = M m k m i m i k Q k i f f min d 1 (3) Q ∑ - = Q i i d d 1 (4)
Schott's metric measures the diversity of a non-dominated set.

Set coverage metric

This metric is based on Zitzler [START_REF] Zitzler | Evolutionary Algorithms for Multiobjective Optimization: Methods and Applications[END_REF]. The metric computes the relative spread of solutions between two non-dominated sets A and B. The set coverage metric C(A, B) calculates the proportion of solutions in B that are weakly dominated by solutions of A:

( ) { } B b a A a B b B A p : ∈ ∃ ∈ = , C (5)
The metric value C(A, B) = 1 means all members of B are weakly dominated by A. On the other hand, C(A, B) = 0 means that no member of B is weakly dominated by A. This operator is not symmetric, thus it is necessary to calculate C(B, A). The set coverage metric measures convergence based on the concept of dominance relations.

Cardinality

This metric counts the number of solutions in the non-dominated set. It measures neither diversity nor convergence.

In order not to limit the description of the quality of a non-dominated set by using only a single metric, a multi-criteria evaluation seems appropriate. Hence, the computations in this paper evaluate the quality of non-dominated sets according to four criteria, which are mentioned above. The multi-criteria tool employed in the computations is Decision Lab [START_REF]Visual Decision: Decision Lab and Decision Lab[END_REF]. The Decision Lab software is a multi-criteria decision making software, which is based on the Preference Ranking Organization Method for Enrichment Evaluations (PROMETHEE) and the Graphical Analysis for Interactive Assistance (GAIA). The details of the PROMETHEE method are found in Brans et al. [START_REF] Brans | How to select and how to rank projects: the PROMETHEE method[END_REF]. GAIA, which makes use of principal component analysis, is a descriptive complement to the PROMETHEE methods [START_REF] Brans | The PROMCALC and GAIA Decision Support System for Multicriteria Decision Aid[END_REF].

Computations and results

Table 1 shows the computed values for each performance criterion in the different multi-objective optimization problems. There are four criteria. A smaller set cardinality is preferred, especially in the case of a continuous decision space. A spread that has smaller value means that the solutions on the non-dominated front are uniformly spaced therefore this criterion is minimized. Hypervolume and set coverage are maximized. Two non-dominated sets are compared in each MOOP. One set is generated by NSGA-II and the other by SPEA2. Decision Lab can rank more nondominated sets but since the set coverage is a binary quality measure, only two nondominated sets can be evaluated each time.

The Multi-objective Shortest Path Problem (MSPP) is an extension of the traditional shortest path problem and is concerned with finding a set of efficient paths with respect to two or more objectives that are usually in conflict. A variety of algorithms and methods such as dynamic programming, label selecting, label correcting, interactive methods, and approximation algorithms have been implemented and investigated with respect to the MSPP [START_REF] Ehrgott | A survey and annotated bibliography of multiobjective combinatorial optimization[END_REF]. The problem is known to be NP-complete. For fitness and selection, two objective configurations are considered for finding efficient paths: (3-S) and (2-S|1-M). S-type objectives are sum problems that are to be minimised and M-type objectives are max-min problems that are to be minimised. A 50-node of 10% density is the basis for the computed values in both MSPP (Multi-objective Shortest Path Problem) configurations. The hypervolume values for the MSPP are blank since they cannot be computed. This reduces the number of MSPP criteria to three.

An evolutionary algorithm is a promising technique to build oblique decision trees (see [START_REF] Cantu-Paz | Inducint oblique decision trees with evolutionary algorithms[END_REF] for a list of advantages). Benchmarking in this field is done by finding oblique partitions on a variety of datasets from the UCI (University of California at Irvine) machine learning repository. The Housing and the Optical Digits datasets are used as the non-dominated sets for the Decision Tree problem.

Decision Lab has a visualization tool that shows the relation between the criteria and the non-dominated set, and shows a preferred solution if it exists. Figure 1 shows a GAIA diagram for the MSPP that shows how each criterion relates to each action. The GAIA plane corresponds to the first principal components of the data, which ensures that a maximum quantity of information is available on the plane. An action in this case refers to a non-dominated set, generated by a specific algorithm. The orientation of the criteria axes indicates which criteria are in agreement with each other. The orientation of the position of an action indicates its strong features. The length of the axis correspond to a criterion's observed deviations between actions, the longer the axis the higher the deviation. The orientation of the pi-vector on the figure, which is the decision axis, points to the preferred action or solution considering all the criteria. In this example, pi does not point towards any action, which means that there is no compromise solution. This due to the condition that NSGA-II has a strong feature on cardinality, SPEA2 has a strong feature on spread, and there is no difference regarding set coverage. All the criteria have the same weight. Adding different weights to each criterion obviously changes the orientation of the pi decision axis.

The criterion weight is independent from the scale of the criterion which means the larger the value the more important the criterion. In order to compare the different criteria independently from their measurement units, the PROMETHEE method provides six preference functions. Five out of six preference functions need some parameters specified by the decision-maker (only the 'usual criterion' does not have parameters). The parameters are used to computer the 'level of preference' of one action over another.

The preference function translates the deviation d between the values of two actions on a single criterion in terms of a preference degree. The preference degree is an increasing function of the deviation, defined on an interval [0,1]. A value 0 has to be interpreted an indifference, while a value of 1 has to be interpreted as strict preference. This research makes use of the 'linear preference function', which is defined by Brans et al. [START_REF] Brans | How to select and how to rank projects: the PROMETHEE method[END_REF] as follows:

( )     > - < ≤ ≤ - = p d p d , p d p p / d d H or if 1 if (6)
As long as d is less than a threshold value p, the preference of the decision maker increases linearly with d. If d becomes larger than p, a strict preference exists. The linear preference function is used in the following computations and is associated with all criteria. The linear preference function is chosen as the function takes into account even the smallest difference in the scores between two alternatives. Table 2 contains some parameter values to be set by the decision-maker. The threshold values are expressed in a relative way. It should be understood that the choice of these values might influence the final decision. That means the values are to be set either on an objective consensus, or on a subjective but also consensus basis, or -in case these options are not available -are values which should be subject to sensitivity analysis.

Figures 1 and2 are used, for the MSPP, to show which non-dominated set has better quality under two preference functions: the 'usual criterion' and the 'linear preference function'. They show the plane for the 2S|1M configuration. The number of criteria has reduced to three as mentioned previously. The values for set coverage in both sets are zero, which means that there are no weakly dominated solutions from each set or that their solutions are similar. Figure 2 shows that SPEA2 is the preferred solution after incorporating the thresholds from Table 2, and the decisive criterion should be interpreted as the 'spread' of solutions. The set coverage criterion is not a factor since both sets do not cover any weakly dominated solutions between them.

Fig. 2. GAIA diagram for a 2S|1M configuration with preferences

With respect to the application, decision trees or classifiers, the fronts shown in Figure 3 are the non-dominated sets from the ensembles of trees generated for one of the UCI datasets, the Housing Set. The Housing set deals with data regarding housing values in the suburbs of Boston (12 attributes, 506 observations). The OC1 (Oblique Classifier 1 [START_REF] Murty | A system for induction of oblique decision trees[END_REF]) solutions are dominated by the solutions of either SPEA2, or NSGA-II, or AP (axis-parallel tests). The non-dominated solutions of the AP classifier are dominated by either SPEA2 or NSGA-II. Most of the non-dominated solutions in NSGA are dominated by solutions in SPEA2. It seems that SPEA2 produces the better non-dominated set through the projection of their non-dominated sets but needs to be validated using the PROMETHEE method. The non-dominated sets of AP and OC1 need not be tested for performance quality as their solutions are dominated by both MOEAs solutions. Figure 4 shows the GAIA plane of the Housing dataset options, and validates that SPEA2 is the preferred solution. The factors that favor SPEA2 are the hypervolume, the spread, and the set coverage. The result does not change when preference thresholds are added. In fact, the pi decision axis leans more to the direction of SPEA2 when preferences are added than it does when without any preferences.

Conclusions

Evolutionary algorithms have a great power of generating non-dominated sets in multi-objective optimization problems. Various strategies deliver various nondominated sets. In order to evaluate the power of such a strategy, several aspects are to be taken into consideration like diversity and convergence to optimality. In literature various metrics have been proposed to measure both performance characteristics. Some metrics measure only one of both, some measure a mix of both. To decide which strategy is better reduces to a multi-criteria problem. In this paper it is illustrated how a plain and clear method, the Promethee method, can be used to decide on the overall quality of an EA strategy. The illustration is made for a combinatorial problem, called the Multi-Objective Shortest Path Problem, and a data mining problem, called the Decision Tree problem.

Fig. 1 .

 1 Fig. 1. GAIA diagram for a 2S|1M configuration using the 'usual criterion'

Fig. 3 .

 3 Fig. 3. Non-dominated fronts for the Housing dataset

Fig. 4 .

 4 Fig.4. GAIA plane for the Housing dataset using the usual criterion

Table 1 .

 1 Computed criteria values of non-dominated sets

	Criterion

Table 2 .

 2 The linear preference function threshold values

	Criterion	Threshold
	Cardinality	75%
	Spread	50%
	Hypervolume	75%
	Set Coverage	50%


Maria Papaioannou

Costas Neocleous
email: costas.neocleous@cut.ac.cy

Anastasis Sofokleous

Nicos Mateou
email: nmateou@mod.gov.cy

Andreas Andreou
email: aandreou@ucy.ac.cy

Christos N Schizas
email: schizas@ucy.ac.cy

Kallipoleos Str

Α Generic Tool for Building Fuzzy Cognitive Map Systems

Keywords: Fuzzy Cognitive Maps, Soft Computing, Intelligent Systems

Α generic system for simulating complex dynamical systems along the paradigm of fuzzy cognitive maps (FCM) has been created and tested. The proposed system enables a user to design appropriate FCM structures, by specifying the desired concepts and the various parameters such as sensitivities, as well as a variety of shaping functions. The user is able to see the results, change the parameters, modify the functions, and rerun the system using an alteration of the final results and make new conclusions. The system is introduced and demonstrated using a simple real case. The results of a usability test of the system suggest that the system is capable of simulating complicated FCM structures in an effective manner, helping the user to reduce the degree of risks during decision making.

Introduction and Background Theory

Fuzzy Cognitive Maps (FCMs) constitute a powerful soft computing modeling method that emerged from the combination of Fuzzy Logic and Neural Networks as an extension of the Cognitive Maps [START_REF] Kosko | A dynamic systems approach to Machine Intelligence[END_REF]. FCMs were introduced first by Kosko [START_REF] Kosko | Fuzzy cognitive maps[END_REF] and since then a wide range of applications in modeling complex dynamic systems have been reported such as medical [START_REF] Papageorgiou | Fuzzy Cognitive Map based decision support system for thyroid diagnosis management[END_REF], environmental [START_REF] Kok | The potential of Fuzzy Cognitive Maps for semi-quantitative scenario development, with an example from Brazil[END_REF], supervisory [START_REF] Stylios | Fuzzy Cognitive Maps in modeling supervisory control systems[END_REF] and political [START_REF] Andreou | The Cyprus Puzzle and the Greek-Turkish Arms Race: Forecasting Developments Using Genetically Evolved Fuzzy Cognitive Maps[END_REF][START_REF] Neocleous | Application of Fuzzy Cognitive Maps to the Political-Economic Problem of Cyprus[END_REF][START_REF] Neocleous | Fuzzy Cognitive Models in Studying Political Dynamics -The case of the Cyprus problem[END_REF] systems. Essentially, a Fuzzy Cognitive Map is developed by integrating the existing experience and knowledge regarding a system. This can be achieved by using a group of human experts to describe the system's structure and behavior in different conditions.

In a graphical form, the FCMs are typically signed fuzzy weighted graphs, usually involving feedbacks, consisting of nodes and directed links connecting them. The nodes represent descriptive behavioral concepts of the system and the links represent cause-effect relations between the concepts [START_REF] Stylios | Modeling Complex Systems Using Fuzzy Cognitive Maps[END_REF]. In the context of FCM theory, the fuzzy value of a concept denotes the degree in which the specific concept is active in the general system, usually bounded in a normalized range of [0, 1]. Furthermore, the weights (or sensitivity) of the system's interrelations reflect the degree of influence between two concepts and they are usually made to vary in the interval [-1,+1] [START_REF] Stylios | Modeling Complex Systems Using Fuzzy Cognitive Maps[END_REF].

Considering the mathematical formulation of the FCM, an FCM system can be represented by a concept state vector A and an adjacency matrix W. The dimension of the concept vector is equal to the total number of distinct concepts n that compose the modeled system. Moreover, the values of its elements correspond to the states of the concepts. The nxn adjacency matrix represents the cause -effect interrelations among the concepts. The value of each cell [i,j] corresponds to the weight of the causeeffect relation between the i th concept and the j th concept.

The basic FCM inference is an iterative process of estimating the new values of the concept vector until they converge. The change of the value of a concept depends on the values of the connected concepts along with the corresponding weights of their interrelations. In practice, at each iteration, the produced concept state vector A t+1 emanates from a vector -matrix multiplication, using the previous concept vector A t and the adjacency matrix W.

After the system interacts for a number of iterations there are 3 possible system behaviors to appear: 1. Fixed point attractor, where the activation levels of all concepts become steady and don't alternate through iterations. 2. Limit Cycle, where an identical sequence of activation levels appears every fixed number of iterations. 3. Chaotic attractor where the activation levels keep changing for ever through iterations.

Additionally, the variation of a concept's impact on another concept may be defined through a manner expressed by a variety of functions (e.g. linear, stepwise, sigmoid, Gaussian and exponential) [START_REF] Koulouriotis | Efficiently modeling and controlling complex dynamic systems using evolutionary fuzzy cognitive maps (invited paper)[END_REF].

A wide range of real life dynamical systems are characterized by causality through their interrelations. Most of the times, experienced and knowledgeable users can identify the parameters and the conditions of their system. However, the complexity of these systems acts as a prohibitive factor on prognosis of their future states. Nevertheless, it is crucial for the users to make a prognosis on the cost of changing a state of a concept, in taking a decision whether they will take such a risk of change or not. FCM modeling constitutes an alternative way of building intelligent systems providing exactly this opportunity to the users; predicting the final states of the system caused by a change on initial concept states. The proposed tool implements the FCM inference mechanism in a simple manner by combining in a simple, easy to use and friendly environment the mechanics of FCM and a user's knowledge. This leads to the modeling of the system in an intelligent framework. An attempt has also been reported in [14].

The main purpose of designing this tool was to create a user-friendly software for building and simulating complex dynamic systems via FCM, with little computer literacy or experience from the part of a user. The only requirement from the user's part is to go through the fuzzification process for the concept values and the weights of their interrelations. The main goal is to help users to design and simulate FCM systems even when they have only the basic knowledge of FCM modeling, provided that they have a good and clear understanding of the parameters and conditions of their system, as well as the corresponding fuzzy values. The long-term objective is to develop a multi-function tool, incorporating certain modeling procedures and algorithms, such as genetic algorithms and neural networks, as well as fuzzification and defuzzification functions. The first edition of the tool implementation, which is limited to simple simulation and testing, is presented in this paper.

A Generic Tool for Running FCM systems

In FCM model, the causality propagation is accomplished through the configuration of the new concept values at each iteration. The proposed system allows the user to choose out of three functions for computing the new concept values.

A ୧ ୲ ൌ ݂ ൭ ൫ A ୨ ୲ିଵ W ୨୧ ൯ ୀଵ ஷ ൱ . ( 1
)[4] A ୧ ୲ ൌ ݂ ൭ ൫ A ୨ ୲ିଵ W ୨୧ A ୧ ୲ିଵ ൯ ୀଵ ஷ ൱ . ( 2
)[9] A ୧ ୲ ൌ ݂ ൭ ൫݇ ଵ A ୨ ୲ିଵ W ୨୧ ݇ ଶ A ୧ ୲ିଵ ൯ ୀଵ ஷ ൱ . ( 3
)[9]
In all the above equations ܣ ௧ is the concept value of concept C i at iteration t, n is the total number of the concepts affecting concept C i , ܹ is the value of the sesnsitivity connecting concept C j with (to) concept C i , k 1 is the coefficient which controls the influence of the interconnected concepts in the calculation of the new concept value ܣ ௧ , k 2 is the coefficient that expresses the contribution of the previous concept value in the calculation of the new one, f(.) is a transformation function which maps each concept value into the space of either [-1,1] or [0,1]. The user is able to choose between sigmoid, hyperbolic, bivalent or trivalent function f(.).

Interaction with the system

The graphical environment of the tool is created using the available interface features that MATLAB 7.9 (R2009 b) provides. The software tool incorporates Microsoft Excel for facilitating the storing and retrieving of data from files. By invoking the tool, a window with a tabbed menu navigation system appears as shown in Figure 1.

Fig. 1. The initial screen of the main window

The tabs divide the window into two screens. By selecting the firs tab "Set FCM Model Parameters" the user is lead to a screen where one can set or modify the FCM model's concepts, their initial values, and the weights of their interrelations. The second tab named Results brings out a screen with information about the final values of the concepts. In the case where no results are available for any simulation of any model, the tab Results is disabled. The user has the flexibility to switch between FCM model's parameters settings and the Results screens whenever one desires. A button named Run is provided for executing the FCM model

Model Building and editing

The user can either build an FCM model from scratch or modify an existing one through the "Set FCM mode's Parameters" screen. The names of the concepts are listed in a list-box, while the initial values of the concepts and their weights are exposed in tables. The first step of building a model is the definition of the names of the concepts. The user can insert, delete, edit the name of a concept or re-order the list of concepts. For each new concept added to the model, the system adds a new row to the "Initial Values of Concepts" and the "Weights between Concepts" table and initializes their numeric values to 0. Each numbered row of the table corresponds to the same numbered concept of the "Concepts" list. The user can edit the numeric values of the tables. There is also a pop-up menu next to each initial concept value that allows the user to lock the corresponding value, by choosing the option "constant" (i.e. execute this concept out of any calculations that alter its initial value).

Running the system and Results

The user can run the FCM model by using the button "Run" located at the right top of the main window. When the algorithm stops executing, the "Results" screen appears (Figure 2).

Fig. 2. The results screen

The final activation levels of the concepts appear on the "Final Concepts Values" list. Each concept in that list is colored and so do the values of the corresponding plot on the graph right next to the list. In the "Results" screen, the main graph presents the variation of the states of the model concepts. It is also possible to select the desired concepts to plot, from the "Final Concept Values" list and then click on the "Plot selected concepts" pushbutton right next to the graph. There is also an edit-box where the user can write comments and conclusions about a simulation. There is also a table displaying information about a simulation (functions used and the number of required iterations). There is also another list on this screen, named "Simulations History". Each time the user runs a model, the system saves all the information about the simulation and the "Simulations History" list adds one more "Simulation" at the bottom of the list. When the user selects a saved Simulation, all the corresponding information is retrieved and displayed. One can delete one or more simulations from the list. The user has also the ability to save a group of simulations in excel file format. The user is thus able to visit the results of each simulation of an FCM model easily and quickly, while one can transfer this simulations file to any computer running WINDOWS, without the need of the FCM tool.

Evaluation of the Tool

The credibility of the tool was tested by modeling and simulating a real process control problem taken from [START_REF] Stylios | Fuzzy Cognitive Maps in modeling supervisory control systems[END_REF]. The results of the FCM model built with the tool agree with the results reported in [START_REF] Stylios | Fuzzy Cognitive Maps in modeling supervisory control systems[END_REF].

For evaluating the usability of the tool, an evaluation questionnaire was developed. The evaluators had to carry out 5 tasks. The tasks cover all the functions a user is able to execute with the tool. After finishing the tasks the user had to answer 10 usability questions suggested in [START_REF] John | SUS -A quick and dirty usability scale[END_REF] and one extra question. Twenty senior computer science, mathematics and computer engineering students were invited to evaluate the tool. Seventeen of the selected students had never used tools for simulating graph models. All of the students had no knowledge about FCM modeling so an essential description of the FCM theory took place just before the evaluation.

Analytically, the 10 questions taken from [START_REF] John | SUS -A quick and dirty usability scale[END_REF] were: 1. I think that I would like to use this system frequently 2. I found the system unnecessarily complex 3. I though the system was easy to use 4. I think I would need the support of a technical person to be able to use this system 5. I found the various functions in this system were well integrated 6. I thought there was too much inconsistency in this system 7. I would imagine that most people would learn to use this system very quickly 8. I needed to learn a lot of things before I could get going with this system 9. I would easily learn how to user this tool 10. I found the system flexible to use And the extra question included in the questionnaires was: I think I would save time by using the tool instead of developing the FCM model without the tool.

The user had to express the degree of agreement or disagreement with each statement on a seven point Likert scale. Each questionnaire was first analyzed by calculating the score from each item (statement). For items 1, 3, 5, 7, 9, 10, 11, as listed above, the score was the scale position minus one. On the other hand, for items 2, 4, 6, 8 the score was seven minus the scale position [START_REF] John | SUS -A quick and dirty usability scale[END_REF]. As a result, there is a mapping between the Likert style user's choice, which ranges in [START_REF] Andreou | The Cyprus Puzzle and the Greek-Turkish Arms Race: Forecasting Developments Using Genetically Evolved Fuzzy Cognitive Maps[END_REF][START_REF] Papageorgiou | Fuzzy Cognitive Map based decision support system for thyroid diagnosis management[END_REF], to a score ranging in [0,6] where, {0,1,2,3,4,5,6} stand for {"very bad", "bad", "tend to be bad", "moderate", "tend to be good", "good", "very good"} system in terms of usability, respectively. As shown in Figure 3, none of the questions was scored as "very bad". Only one question had 5% scores in "bad" category. Only the 5% of four questions were scored below "moderate". Considering the worst scores, Question 6 Fig. 3. Histogram presenting the percentage of the answers distribution per each question. The grey dashed lines divide the histogram into the 7 answer categories.had 10% scores is the only one that had 10% scores falling in both "bad" and "tend to be bad" categories. After the first analysis of the questionnaires the students were asked why they thought there wasn't satisfactory consistency in the system. The users stated that they couldn't understand the word "consistency" and they got confused. Question 2, had also answers in the group of "tend to be bad","moderate" and "tend to be good" categories. A possible reason for that could be that the users felt they had to memorize many ways for utilizing data, since the first tab screen is about editing tables and a listbox, while the second one is more about pushing buttons and utilizing plots. 25% of the answers of Question 4 were scored as "tend to be bad", "moderate" and "tend to be good". A possible reason why users felt that way is that they didn't have any experience with graph models simulation tools or FCM. The same reason applies for the "moderate" and "tend to be good" scores of Question 3. Hence, it is absolutely reasonable that some of the users would face some difficulties using for the first time the tool and they would ask for help from a technical person. Besides, 95% of the answers of Question 9 fall in "good" and "very good" categories, stating that 95% of the users believe they would easily learn how to use the tool. This implies that the users understand that the difficulties they had to resolve during the evaluation came from their inexperience and lack of knowledge. Furthermore, it is noteworthy that all eleven questions had 70% -100% of their answers falling into categories "good" and "very good" indicating that most of the users thought the system was indeed successful in terms of usability. This fact is also reflected in Figure 4, whereas all question medians range between five and six (good and very good) and approximately 50% of them equal to very good which is the optimum score.

Conclusions and Future Work

Dynamical complex systems are met in most human decision making activities. Decision making problems concerning such systems can be solved by utilizing FCM modeling methodology. The proposed tool enables a user, who is not familiar with computer literacy, to design and simulate an FCM modeled system and benefit from the advantageous features of FCM modeling.

The tool and its capabilities were presented and described in detail in the presented paper. The friendly graphical interface helps the user to easily insert, save, retrieve and display data. The main advantages of the tool is the preservation of the information of multiple simulations of different FCM models at the same time, the display of the plot of each model final states in different windows and the feature of saving information of multiple simulations in an spreadsheet file so that the user is able to go back to the results of any simulation of any model without the need to rebuild and re-simulate the system on the tool.

The usability of the system was evaluated by developing an evaluation questionnaire. The conclusions drawn from the results of the evaluation are (1) the tool meets the usability standards as defined in [START_REF] John | SUS -A quick and dirty usability scale[END_REF] since the majority of the students gave very good scores for the usability of the tool through their answers (2) users having the basic knowledge of FCM theory are capable of using the tool successfully in designing and simulating FCM systems even if they had never interacted with FCM modeling systems in the past (3) the tool is easy to learn and easy to use.

The tool is expected to introduce common people to the challenging world of FCMs, by offering them a user-friendly functional environment which allows them to create and simulate FCM models.

In the future, the tool will be enhanced with pre-processing functions, helping and enabling the user to extract the fuzzy rules which characterize a system and embed this knowledge to the tool. The user will be also able to "run" a demonstration of built examples of the tool. In addition, the notion of time, in terms of time dependencies, as well as an algorithm inspired by the theory of evaluation of species will be introduced in the model.

Fig. 4 .

 4 Fig. 4. Histogram presenting the median of the answers of each question


Charalampos Doukas
email: doukas@aegean.gr

Ilias Maglogiannis

A Fast Mobile Face Recognition System for Android OS based on Eigenfaces Decomposition

Keywords: Mobile face recognition, Android OS, Face Detection, Eigenfaces decomposition

This paper presents a speed-optimized face recognition system designed for mobile devices. Such applications may be used in the context of pervasive and assistive computing for the support of elderly suffering from dementia in recognizing persons or for the development of cognitive memory games Eigenfaces decomposition and Mahalanobis distance calculation have been utilized whereas the recognition application has been developed for Android OS. The initial implementation and the corresponding results have proven the feasibility and value of the proposed system.

Introduction

Face recognition refers to automatically identifying or verifying a person from a digital image or a video frame from a video source. The process involves the face detection in images or image sequences, proper feature extraction, the creation of a reference set based on known faces and the detection of new faces based on the latter. Mobile face recognition can be utilized either for user authentication or for cognitive assistance in identifying individuals. The latter can be particularly useful for the elderly whose face recognition ability declines while aging. A plethora of face recognition algorithms have been proposed in literature ( [START_REF] Chellappa | Humain and machine recognition of faces: A survey[END_REF] - [START_REF] Turk | Eigenfaces for recognition[END_REF]). However, the application of the well known and established techniques to mobile devices, such as mobile phones, raises some challenges related to two major issues: the inefficiency of the algorithms and the limited processing capabilities of the mobile devices.

This paper presents a speed-optimized method for performing face recognition on Android mobile phones. Both training and recognition phases are performed on the device. The Android mobile operating system provides built-in functionality for detecting faces on still images. Eigenface decomposition is used in order to extract important features from the images and a training vector is created based on an initial image training set. The Mahalanobis cosine distance is utilized for measuring the similarity of the given train set images and the examined image in the Eigenface projection space. The latter methods have been selected for implementation on a mobile device due to lower complexity and memory requirements compared to other methods available in literature. In addition to the face recognition process performed thoroughly on the device, a web service has been developed that allows the feature extraction and classification using more advanced dace recognition techniques. The latter method requires network connectivity from the mobile device but can provide more accurate results.

Initial evaluation of the system has been performed on Google's G1 mobile phone and results are quite promising regarding the system's performance and accuracy. The paper presents the aforementioned, discusses related work and includes technical information and implementation details.

Related Work and Background Information

In general, the whole process of face recognition consists of detecting and selecting faces from the current image or video frame utilizing techniques like skin detection, process the selected image region and finally applying a recognition algorithm.

Several classifiers have been proposed in the literature e.g. minimum distance classification in the eigenspace ( [START_REF] Sirovitch | Low-dimensional procedure for the characterization of human faces[END_REF], [START_REF] Turk | Eigenfaces for recognition[END_REF]), Fisher's discriminant analysis [START_REF] Belhumeur | Eigenfaces vs fisherfaces: Recognition using class specific linear projection[END_REF], and neural networks [START_REF] Fleming | Categorization of faces using unsupervised feature extraction[END_REF]. Global techniques work well for classifying frontal views of faces, however they are not robust against pose changes since global features are highly sensitive to translation and rotation of the face. To avoid this problem an alignment stage can be added before classifying the face. Aligning an input face image with a reference face image requires computing correspondence between the two face images. The correspondence is usually determined for a small number of prominent points in the face like the center of the eye, the nostrils, or the corners of the mouth. The image region is also converted to grey levels and equalized, then cropped with an elliptic mask and normalized. Finally, the recognition algorithm compares the processed image to a set of prerecorded faces, all belonging to the same subject, to find the best match. This algorithm usually returns also a measure of similarity between the current face and the reference template.

Face Detection and Preprocessing

One of the most crucial parts in the presented face recognition system is the preprocessing of the considered image region. To align a face horizontally, a common technique consists in locating the eyes and then rotating the selected image so that their inclination is null. Also, the distance between the eyes is used to resize the same image to a given value of width and height.

Once the positions (uR, vR) and (uL, vL), of the right and left eye respectively, have been determined, we simply calculate the angle and the distance between them as follows:

  1 tan vL vR RL uL uR a     (1) 2 2 (
) ( )

RL L R L R d u u v v     (2)
The face is then rotated of an angle -a RL in order to align the eyes, and then resized to 24x24 pixels, i.e. the size of the templates used for comparison. Rotation and scaling, centered on the right eye (u R , v R ), are performed with a simple affine transformation as follows:

(1 ) ' '

(1 )

1 R R R R u a b a u bv u v v b a bu a v                             cos( ) RL a s a   sin( ) RL b s a   (3)
where (u, v) are the pixel coordinates of the source image and (u′, v′) those of the destination. In order to avoid possible outliers in the final image after rotation, the affine transformation is actually applied to a sub-region slightly bigger than the original face bounding box. The final step of the image processing consists in cropping the face's area with an elliptical mask. This reduces the influence of hair and background pixels at the four corners of the rectangular region.

Eigenfaces Decomposition and Similarity Detection

After the image processing described above, the most popular techniques for face recognition, called Eigenfaces [START_REF] Turk | Eigenfaces for recognition[END_REF] is applied. Eigenfaces decomposition and similarity detection relies on measuring the similarity between a new face image and a reference one, projecting both the images into an eigenspace, previously created by training, and calculating the distance between the projections. The idea behind the eigenface technique is to extract the relevant information contained in a facial image and represent it as efficiently as possible. Rather than manipulating and comparing faces directly, one manipulates and compares their representations.

To create a set of eigenfaces, one must prepare a training set of face images. The pictures constituting the training set should have been taken under the same lighting conditions, and must be normalized to have the eyes and mouths aligned across all images. They must also be all resampled to the same pixel resolution. Each image is treated as one vector, simply by concatenating the rows of pixels in the original image, resulting in a single row with r × c elements. For this implementation, it is assumed that all images of the training set are stored in a single matrix T, where each row of the matrix is an image. The mean is then subtracted. The average image a has to be calculated and then subtracted from each original image in T. The eigenvectors and eigenvalues of the covariance matrix S are thereafter calculated. Each eigenvector has the same dimensionality (number of components) as the original images, and thus can itself be seen as an image. The eigenvectors of this covariance matrix are therefore called eigenfaces (see Fig. 1). They are the directions in which the images differ from the mean image. The D x D covariance matrix will result in D eigenvectors, each representing a direction in the r × c dimensional image space. The eigenvectors (eigenfaces) with largest associated eigenvalue are kept. The latter can now be used to represent both existing and new faces: we can project a new (meansubtracted) image on the eigenfaces and thereby record how that new face differs from the mean face. The eigenvalues associated with each eigenface represent how much the images in the training set vary from the mean image in that direction.

The computation of eigenvectors is performed as follows:

Let T be the matrix of preprocessed training examples, where each row contains one mean-subtracted image. The covariance matrix can then be computed as S = T T T and the eigenvector decomposition of S is given by:

T i i ii Sv T Tv v    (4)
However TTT is a large matrix, and if instead we take the eigenvalue decomposition of:

T i i i TT u u   (5)
then we notice that by pre-multiplying both sides of the equation with TT, we obtain:

T T T i i T TT ui T u   (6)
Meaning that, if u i is an eigenvector of TT T , then v i =T T u i is an eigenvector of S.

Regarding measuring the similarity between different eigenvectors and performing image recognition, several methods have been proposed. The most common one, which is considered to give the best performance, is the Mahalanobis Cosine distance ( [START_REF] Saleh Aly | Face recognition under varying illumination using Mahalanobis self-organizing map[END_REF], [START_REF] Beveridge | The CSU Face Identification Evaluation System User's Guide: Version 5.0[END_REF]). As the name suggests, this is given by the cosine of the angle θij between the two face projections fi and fj in the Mahalanobis space:

( , ) cos The mobile face recognition system has been developed utilizing Google's Android mobile Operating System (OS) [START_REF]The Android mobile OS by GoogleTM[END_REF] and the appropriate software development kit (sdk). Android is a mobile operating system running on the Linux kernel. Several mobile device vendors already support it. The platform is adaptable to larger and traditional smart phone layouts and supports a variety of connectivity technologies (CDMA, EV-DO, UMTS, Bluetooth, and Wi-Fi). It supports a great variety of audio, video and still image format, making it suitable for capturing and displaying images. Finally, it supports native multi-touch technology, which increases the application's usability. In addition to the face recognition process implemented locally on the mobile device, a web service has bee also developed that enables a more advanced feature extraction and classification of the face images. An appropriate module enables the communication between the Android application and the Web Service utilizing SOAP messages. The Web Service has been developed using Java and Apache Axis [START_REF]The Apache Axis Framework[END_REF]. The classification module is based on the same java code developed for the Android for pre-processing the images and extracting the eigenvectors. Classification of the latter is performed using the Weka tool [START_REF] Hall | The WEKA Data Mining Software: An Update[END_REF] that contains a great variety of tools for pattern recognition. The data exchange between the mobile device and the web service involves the exchange of images and the classification results.

fi fj fi fj fi fj ij        (7) (a) (b) (c) (d)
Fig. 3 illustrates the main steps for performing face recognition on an Androidenabled mobile device. When faces are detected in an image, a validation of the faces found and the corresponding eye coordinates and eye distances is performed in order to minimize erroneous cases.

Initial Evaluation Results

In order to evaluate the proposed system a HTC G1 [START_REF]HTC G1 mobile phone[END_REF] mobile phone running Android OS version 1.6 has been used in order to host the developed face recognition application. The initial training set consisted of 126 images containing 7 different faces (4 males and 3 males, 18 sample images per face). All faces have been cropped from the background as described in Section 2 and resized to a 100x100 pixels resolution. All images reside into the mobile phone's memory card and the developed face recognition application has been initially used in order to generate the corresponding eigenfaces. The first 50 eigenvalues have been for each image and the corresponding vectors have been saved into the mobile's phone memory card. For the evaluation, 5 different images (i.e. not included in the training set) from each individual have been used. The creation of the training set has taken an average of 65 seconds for all 126 images on the device, whereas the processing, eigenfaces decomposition and mahalanobis distance calculation and comparison for each evaluation image has taken approximately only 1.2 seconds. The overall system accuracy in detecting faces in images is about 96%, whereas the total accuracy in recognizing faces that do not belong to the training set is 80%. The latter are quite promising considering the size of the evaluation set and the fact that both the training and recognition phases are performed on the mobile device at very reasonable time delays, allowing the usage of the proposed face recognition system in real applications. Regarding the ability of the system to utilize more advanced classifiers through the Weka tool for more accurate face recognition, a second experiment has been conducted. The same training set has been used, whereas for the evaluation process, the latter set along with new image samples have been used through 10-fold cross validation. Several classifiers like Bayes Networks, Support Vector Machines (SVM), and Neural Networks (NN) have been validated. The best performance has been achieved by SVM approaching a total accuracy of 98%, followed by 96.2% achieved by NN and 95% by the Bayes classifier respectively. The total time to process and classify an image is less than 1 sec at this case, but the overall network transmission of data introduces further delay until the result is displayed at the mobile device.

Conclusions

The advances of state of the art smart phones and PDAs have enabled the creation of fast and reasonably robust face detection/recognition applications like the one presented in this work. Such applications may be used in the context of pervasive and assistive computing for the support of elderly suffering from dementia in recognizing persons or for the development of cognitive memory games.

An important issue to be handled in face recognition is the "Face Variation" problem. This is closely related to changes of illumination, pose, and expression. As non uniform illumination affects the acquired face image, it is important to know how the recognition algorithm works under varying illumination. Pose variations include approaches on how to model faces seen from different points of view. Expression variations refer to how we can design algorithms for the recognition of emotions and other facial expressions [START_REF] Maglogiannis | Face Detection and Recognition of Human Emotion Using Markov Random Fields[END_REF]. Future research will focus on addressing such issues, on improvement of detection precision, the addition of online training module and other novel applications.

Fig. 1 .

 1 Fig. 1. Example illustrating eigenfaces decomposition: a) original face image, b) first eigenface, c) tenth eigenface, d) hundredth eigenface.

Fig. 2 .

 2 Fig. 2. Architecture of the proposed mobile face recognition system.

Fig. 2

 2 Fig.2presents an illustration of the proposed system architecture. The mobile application utilizes native Android resources for capturing images using the mobile devices' image recording capabilities. The images can also be allocated on a remote or local repository. Android also provides native face detection functionality. Appropriate java methods can be used directly in any custom application and allow the detection of faces within images. The methods provide the coordinates of the eyes. Appropriate java classes (referred as modules) have been developed for performing the initial image processing (i.e., face cropping, rotation and scaling), the decomposition of the initial image into eigenfaces and the projection of the eigenvactors into the Mahalanobis domain and distance calculation between a given new image instance and the ones used as a training set. The images and the computed eigenvectors of the training set reside within the mobile device's storage repository (e.g., a memory card).In addition to the face recognition process implemented locally on the mobile device, a web service has bee also developed that enables a more advanced feature extraction and classification of the face images. An appropriate module enables the communication between the Android application and the Web Service utilizing SOAP messages. The Web Service has been developed using Java and Apache Axis[START_REF]The Apache Axis Framework[END_REF]. The classification module is based on the same java code developed for the Android

Fig. 3 .

 3 Fig. 3. Data flow diagram illustrating the basic steps for face detection and recognition using eigenface decomposition

Table 1 .

 1 Initial evaluation results for the proposed mobile face recognition system. System accuracy in face detection within images and overall face recognition accuracy are presented. Both training and recognition phases have been performed on the mobile device.

		Accuracy	in	Accuracy in Face
		detecting faces in	Recognition (%)
		images in training set
		(%)
	Male 1	94.4		80
	Male 2	100.0		80
	Male 3	100.0		100
	Male 4	88.8		60
	Female 1	88.8		80
	Female 2	100.0		80
	Female 3	100.0		80