- Email: [email protected]
A genetic algorithm for circulant Euclidean distance matrices
APPelLeED MAT~EMAT~C5 AND
CON NOTAT~ON
EI_SEVIER
Applied Mathematics and Computation 97 (1998) 1 9 7 2 0 8
A genetic algorithm for circulant Euclidean distance matrices Wilson Rivera-Gallego 1 NSF Engineering Research Center for Computational Field Simulation, Mississippi State University, Mississippi State, MS 39762, USA
Abstract
This paper presents a fast genetic algorithm to determine three-dimensional configurations of points that generate circulant Euclidean Distance Matrices (EDMs). A parallel implementation is possible by using the message passing interface (MPI) standard. In addition, theoretical results about the polyhedral structure of both the cone of circulant symmetric positive semidefinite matrices and the cone of circulant EDMs are introduced. © 1998 Elsevier Science Inc. All rights reserved.
1. Introduction
The first known references of Euclidean Distance Matrices (EDMs) appeared in multi-dimensional scaling [1]. Members of a multi-dimensional sample are represented by points in some geometrical space, usually Euclidean. It is possible to build up arrays of numbers, called ordering diagrams, where pairs of close points usually represent similar samples and faraway points represent dissimilar samples. The goal is to use this dissimilarity matrix in order to describe some characteristics of the sample. For example, in psychometric each element in the dissimilarity matrix may represent the dissimilarity of a pair of stimuli. It is required to represent the set of stimuli by a configuration of points in a low-dimensional Euclidean space [2]. However, it is in molecular conformation [3] where stronger applications of EDMs are found. For example, interatomic distances are measured by using
1 E-mail: [email protected]. 0096-3003/98/$ - see front matter © 1998 Elsevier Science Inc, All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 9 7 ) 1 0 1 43-6
198
w. Rivera-GallegoI Appl. Math. Comput. 97 (1998) 197-208
nuclear magnetic resonance (NMR) and they are located in a rectangular matrix. It is necessary to determine the possible three-dimensional configurations of points which would generate the matrix [4]. Also, it is interesting to observe the structure of the configurations from the properties of the matrix [5]. Schoenberg [6] found a natural parametrization of the EDMs cone from the symmetric positive semidefinite matrices. Gower [7,8] described some properties of the EDMs and pointed out some connections between these properties and the configurations of points that generate the matrices. Critchley [9] worked on the spectral decomposition of the EDMs. Hayden et al. [5] obtained Gower's results, but in a simple way offering us a more geometrical view of the EDMs cone. Rivera [10] studied some special structures of EDMs. Genetic Algorithms (GAs) are search algorithms based on the genetic processes observed in natural evolution. Theoretical concepts were first introduced by Holland [11] and posteriorly well described by Golberg [12]. Basically, GAs work with an initial population of individuals. Each individual (chromosome) consists of a string of characters (genes) and represents a possible solution in a search space. A fitness value is associated with each individual. By using genetic operators such as selection, crossover (mating between individuals), mutation (random modification on genes), offsprings with better fitness values are generated. This process keeps on in order to find an optimal solution to a particular problem. Section 2 introduces the definitions. In Section 3, theoretical results are presented. In Section 4, an outline of genetic algorithms is given and some details of its implementation for circulant EDMs are discussed. Section 5 deals with discussion of practical results. Finally, conclusions are given in Section 6. 2. Euclidean distance matrices
If there exist n points, pl,P2,... ,p,, in ~k such that [[Pi __pj[[2 = d2 then D = [d2.] E ~,×n is an EDM. The minimum value of k is called embedding dimension. A matrix symmetric A E R n×" is positive semidefinite if x'rAx >10 for all x E R n, x ¢ 0. In particular, if it is possible to write a matrix A c R "×" as a product A -----CC T then A is symmetric positive semidefinite. A matrix C E Rk×k such that the ith row corresponds to the coordinates of the point p~ is called Coordinates Matrix. Note that given a configuration of n points Pl,P2,... ,p, in Rk and the corresponding coordinates matrix C then the matrix B = CC T is symmetric positive semidefinite. Also
bij =pXpj and in particular
for i,j = 1 , 2 , . . . , n
W. Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197-208
b,i
= Ilp~ll 2
199
for i = 1 , 2 , . . . , n .
2.1. P a r a m e t r i z a t i o n
of EDM
Let us define: f2n = {B E Rn×n:B is a Symmetric positive semidefinite matrix}, An = {D E Rn×n:D is an Euclidean Distance Matrix}.
(1)
It is clear that d~ = IlPillz + [IPAI2 - 2p•pj - - bii -~- bjj - 2bij.
(2)
Because B = [bij] is a positive semidefinite matrix, the above expression is a natural parametrization of the EDMs ~c : (2 n --+ A n ,
[l£(n)]ij -~- bii -~- bjj - 2bij.
(3)
The map x can be expressed in matrix form by x ( B ) = b e ~ + e b T - 2B,
(4)
where b = diag(B) and e is the vector whose components are all equal to 1. Each E D M corresponds to different positive semidefinite matrices depending on the coordinates system. Thus, the function x is not injective. In order to get an inverse transformation of x, it is necessary to define the map Zs : An ---, On, z s ( O ) = - ½ ( I - e s T ) O ( I -- seT),
SE~"
(5)
such t h a t s T e = l ,
where s is a vector in R n which fixes the coordinates system. In case that s = e / n , the origin of the coordinates system corresponds to the centroid of the configuration of points [5,10].
2.2. C i r c u l a n t m a t r i c e s
A matrix A = [aij] E C n×n is called C i r c u l a n t M a t r i x if aij = ak
j - i + 1 - k (mod n).
We denote a circulant matrix by A = circ(al, a 2 , . . . , an).
W. Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197-208
200
Let F be the Fourier M a t r i x given by 1
f ~
1
1
1
...
1
1
w
w2
...
w(n-1)
1.
W2 . .
W4
"'"
w2(.-1)
1
w (~-1)
w 2(n-l/
...
w(n-l)(,-l)
.
(6)
where w = cos (2n/n) + i sin (2~/n). As pointed out in [13], a circulant matrix A E C "×" can be expressed as C = FAF*, where A = d i a g o n a l ( 2 1 , 2 2 , . . . , 2,), 2j = Pp(w/-'),
(7)
Pp(z) = cl + CzZ ~- " " " q- CnZ n - l .
This result and the s y m m e t r y of circulant E D M s can be used to generate a nice representation of this set o f matrices. We point out that the circulant distance matrices set is a polyhedral cone whose extreme points can be identified easily.
3. T h e o r e t i c a l results T h e o r e m 3.1. The cone o f circulant symmetric positive semidefinite matrices o f
order n is a polyhedral cone with extreme points given by
Aj=~T
j=I,...,L,
(8)
where
= ~fi--[cos 2 ( i - 1 ) ( j - 1 ) n n Vnl
sin 2 ( i - 1 ) ( j - 1 ) r c ] n
E ~nx2,
(9)
/
L = [n/2J i f n i s o d d a n d L = [n/2J + 1 i f n iseven. • = 1 i f j = [n/2J + 1, and = 2 otherwise.
Proof. Given a matrix circulant and symmetric C = circ(cl, c2, • • •, cn), we have that Cj = Cn_j+ 2
for j = 2 , 3 , . . . , n .
Hence, 2j=2,_j+2
for j = 2 , 3 , . . . , n .
(10)
W. Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197-208
201
F r o m (7) n
C----
Zj=l •'JY; ,
where fj = ~ 1 [1, w (j-x), w2~-1),... , w(,_~)0. 1)]T '
,
=
1 [ 1 , w - " - ' ) , w -2"-1> , .
. . , w_(n_X)(j_l)]
,
Thus, f f 7 = nL1[w(j_li(;_~) ]j (;,k) for i,k = 1,2, ... ,n.
(11)
Also f~_j+2 = f .
for j = 2 , 3 , . . . , n ,
f~-j+zf,*_j+2 = J~J~*
(12)
_ - -1 [w_(j_,lti_k) ] n
for i , k = 1 2, (i,k)
'
. n.
(13)
""
Because the vectors f,-j+2 and fj c o r r e s p o n d to the same eigenvalue, the matrix C can be d e c o m p o s e d as Ln/2J + 1
C = 21f, f ; + Z
2J(~f7 +f,-j+zf,*-j+2)
(14)
j=2
if n is odd, or tn/2J
C = Zlflf~ + )~kfkf; + Z 2 j ( f j f j * +fn-j+zf;-j+2)
(15)
j=2
if n is even, with k = [n/2J + 1. We define Aj = f j f j + f,-j+zf,*_j+2.
(16)
It is very easy to show that Aj = ~ / 0 f
j = 1,...,L.
(17)
On the other hand, it can be verified that if ss x is a circulant matrix o f rank 1, then there exists some a E ~ such that s = ~fl if n is odd, and s = ef~ or s = c~f, if n is even. In order to show this statement, it is sufficient to see that if s [S1,$2,... ,S,] x then for some ~ 6 R =
[si] = c~
i = 1 , 2 , . . . ,n.
202
w. Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197-208
Because ss x must have constant diagonal. Also, because C is circulant Cin:Ci+l, 1
for i = 1 , 2 , . . . , n -
1,
or equivalently SiS n =
Si+lS 1
for i = 1 , 2 , . . . , n -- 1.
By expanding the above expression, the statement can be verified.
[]
It is easy to show that B is a circulant positive semidefinite matrix if and only if x(B) is a circulant EDM. In addition, K is a linear transformation. As a consequence, we have the following corollary.
Corollary 3.2. The circulant E D M s set is a polyhedral cone whose extreme points are images f r o m the extreme points o f the circulant symmetric positive semidefinite matrices under the application ~. 3.1. Comments
In practice, it may be useful to determine substructures with special properties within molecular configurations. In particular, a subset of points within a configuration may generate a nearly circulant EDM. Thus, the characterization of the extreme points of the circulant symmetric positive semidefinite matrices becomes useful in order to determine the three-dimensional coordinates of the points. Our G A implementation for circulant EDMs makes good use of this characterization. Notice that, because 21flf; is a constant matrix and x(21flf~) = O, f l f ; does not have any effect on the corresponding EDM. Thus, for n odd the rank of the extreme points is 2 and for n even the rank is also 2 except for f~.fk* which has unitary rank. Hence, three-dimensional configurations that generate circulant E D M only there exist for an even number of points. Moreover, if we are looking for this sort of three-dimensional configurations, we only need to find a combination of one of the eigenvalues 2j for j = 2 , . . . , Ln/2J and 2t~/2j+l.
4. Genetic algorithms and circulant euclidean distance matrices A general G A is as follows: Generate r a n d o m l y an initial population. determine fitness of individuals. do
( S e l e c t i o n operator. C r o s s o v e r operator.
203
Rivera-GallegolAppL Math. Comput. 97(19~)197-2~
Mutation operator. Determine fitness of individuals. while (not convergence)
Individuals are selected from the population to mating processes by a selection operator. The probability of an individual being selected is proportional to its fitness. Offsprings are created by combining genes from the parents. This is made by crossover operators. For example, in one point crossover the parents' chromosomes are cut at a random crossover point and the sections are exchanged. In a uniform crossover each gene in the offspring is generated by copying randomly the gene from one or the other parent. In addition, the mutation operator changes genes randomly selected. Mutation is used to recover information and keep diversity in the population. Each individual is represented by the structure typedef struct
( int chrom[ARRAY]; double Fitness, values[2]; int parentl, parent2; ) individual;
Genetic information is stored into the binary array chrom. Recall that, in order to find a three-dimensional configuration, it is only necessary to know two eigenvalues: one of the eigenvalues 2j for j = 2 , . . . , [n/2J and the eigenvalue J,Ln/eJ+l. The decode function translates the information in chrom by dividing it into two parts, each one corresponding to one the eigenvalues, and stores the results into the array values. The length of the array chrom depends on the parameter called precision. The genetic information at each part, corresponding to an eigenvalue, is as follows: 2Precision
.
.
.
20
.
.
.
2Precision-1
Thus, the length of chrom is 4 × precision. Notice that the eigenvalues must be positive values. Hence, it is not necessary to worry about the sign. Each individual has a fitness value corresponding to the Frobenius norm of the difference between the circulant EDM generated by the individual and the input matrix. A tolerance parameter regulates the generation of individuals at each step. Usually, the initial value of tolerance is 2 Maxgen where Maxgen is the maximum number of allowed generations. At each generation, the tolerance is reduced to
204
w. Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197-208
the half guaranteeing those better individuals are selected. Thus, the G A stops when the Maxgen value is reached or after 10000 intents to select an individual that satisfies the tolerance value. Because there are Ln/2J feasible three-dimensional configurations, each one corresponding to a combination of one 2; for j = 2 , . . . , [n/2J and /~[n/2J+l the G A has a natural parallelization. The parallelized version is build up by using the Message Passing Interface (MPI) [14]. M P I has the particularity to be portable to different platforms.
5. Practical results
The code has been run on a SGI Onys Reality Engine Server with 12 processors. Different problems have been tested and the G A has shown to be a very robust algorithm. For example, by using crossover probability = 0.6, mutation
Table 1 Points
Time (min)
10
1.24
20 50 100
3.12 5.75 7.23
0.5-
-0.5.
i
-1
-1
Fig, 1.22 = 2.0, 25 = 4.125, + (real points), * (approximations).
205
w. Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197 208
probability = 0.05, Maxgen = 20, population size = 20 a n d t a k i n g different configurations with different n u m b e r o f points, the results are given in T a b l e 1. T h e time is the average a m o n g the times f r o m all o f the cases tested for a specific n u m b e r o f points. T h e time spent b y the a l g o r i t h m d e p e n d s on the c o m p l e x i t y o f the p r o b l e m . In p a r t i c u l a r , h o w m u c h closely circulant is the inp u t matrix.
E x a m p l e 5.1. F o r n = 8 the eigenvalues 2.2 = 2.0 a n d 25 = 4.0 generate the following c i r c u l a n t E D M : 0
2.2929
1.0000
3.7071
2.0000
3.7071
1.0000
2.2929'
2.2929
0
2.2929
1.0000
3.7071
2.0000
3.7071
1.0000
1.0000
2.2929
0
2.2929
1.0000
3.7071
2.0000
3.7071
3.7071
1.0000
2.2929
0
2.2929
1.0000
3.7071
2.0000
2.0000
3.7071
1.0000
2.2929
0
2.2929
1.0000
3.7071
3.7071
2.0000
3.7071
1.0000
2.2929
0
2.2929
1.0000
1.0000
3.7071
2.0000
3.7071
1.0000
2.2929
0
2.2929
2.2929
1.0000
3.7071
2.0000
3.7071
1.0000
2.2929
0
+
+
0.5. +
4-
-0.5,
+ -1:
-0.5 ~ - 0 . 5 -1 -1
Fig. 2. 23 = 1.75, 25 = 6.625, + (real points), * (approximations).
(18)
W. Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197-208
206
O.
~
+
1 "10 . 5 ~
+
~ ~
-0.5 " ~ . ~ ' ~ - 0 . 5 -1 -1
Fig. 3.24 = 1.875, 25 = 6.625, + (real points), * (approximations).
We introduce errors in the measure o f the distances. Thus, our input matrix is 0 2.2729 0.9600 3.6671 1.9400 3.6771 1.0300 2.2229
2.2729 0 2.2629 0.9700 3.7371 1.9700 3.657l 0.9200
0.9600 2.2629 0 2.3129 1.0500 3.7571 1.9800 3.6371
3.6671 0.9700 2.3129 0 2.2429 0.9200 3.4071 2.0800
1.9400 3.7371 1.0500 2.2429 0 2.1929 0.9400 3.6371
3.6771 1.9700 3.7571 0.9200 2.1929 0 2.3729 0.9000
1.0300 2.2229 3.6571 0.9200 1.9800 3.6371 3.4071 2.0800 0.9400 3.6371 2.3729 0.9000 0 2.2229 2.2229 0
(19)
R u n n i n g our algorithm with a population o f 20 individuals and Maxgen = 20, we obtain three configurations. The first one (Fig. 1) corresponds to 22 = 2.0, 25 = 4 . 1 2 5 and fitness=0.81. The second one (Fig. 2) corresponds to 23 = 1.75, 25 = 6.625 and fitness = 10.29. The last one (Fig. 3) is that corresponding to 24 = 1.875, 25 = 6.625 and fitness = 10.62. It is interesting to notice the structure o f three-dimensional configurations that generate circulant E D M s . Fig. 4 shows an example for a configuration with 20 points.
6. Conclusions A fast genetic algorithm to determine three-dimensional configurations that generate circulant E D M s has been developed. A parallel version o f the algo-
W. Rivera-Gallego I Appl. Math. Comput. 97 (1998) 197-208 "5
207
"8
0.5, "~9 ~7
~t3
4t5
a4 ~B
~o
-0.5> 0.5 .5
-0,5 .0.5
Fig. 4. n = 20, 22 = 2.0, and 211 = 4.0.
rithm was implemented by using MPI standard. Interesting theoretical results regarding the polyhedral structure of the cone of circulant symmetric positive semidefinite matrices and the cone of circulant EDMs were introduced. Recognition of substructures within molecular configurations is very important in order to determine more complex structures. The results obtained in this work may be useful to explore other configuration such as Toeplitz structures. Moreover, our genetic algorithm can be easily modified in order to solve other problems such as matrix completion. It would be interesting to use the knowledge from the structure of the set of EDMs to generate domain knowledge genetic algorithms. References [1] J. De Leuuw, W. Heiser, Theory of Multidimensional Scaling, in: P.R. Krishnaiah (Ed.), Handbook of Statistics, North-Holland, Amsterdam, 1982. [2] G. Young, A.S. Householder, Discussion of a set of points in terms of their mutual distances, Psychometrika 3 (1938) 19-22. [3] G.M. Crippen, T.F. Havel, Distance Geometry and Molecular Conformation, Wiley, New York, 1988. [4] W. Glunt, T.L. Hayden, M. Raydan, Molecular conformations from distance matrices, J. Comput. Chem. 1 (1993) 114-120. [5] T.L. Hayden, J. Wells, W. Liu, P. Tarazaga, The cone of distance matrices, Linear Algebra Appl. 144 (1991) 153-169. [6] I. Schoenberg, Remarks to Maurice Frechet's Article, Sur la definition axiomatique d' une classe d'espaces vectoriels distanci~s applicables vectoriellement sur L'espace de Hilbert, Ann. Math. 36 (1935) 724-732.
208
I,K Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197-208
[7] J.C. Gower, Euclidean distance geometry, Math. Scientist 7 (1982) 1 14. [8] J.C. Gower, Properties of euclidean and non-euclidean distance matrices, Linear Algebra Appl. 67 (1985) 81-97. [9] F. Critchley, On certain linear mappings between inner product and squared-distance matrices, Linear Algebra Appl. 105 (1988) 91 107. [10] W. Rivera, Matrices de distancia con estructuras especiales. M.Sc. Thesis, University of Puerto Rico, Mayaguez, 1994. [11] J.H. Holland, Adaptation of Natural and Artificial Systems. University of Michigan Press, Ann Arbor, MI, 1975. [12] D.E. Golberg, Genetic Algorithms in Search, Optimization and Machine Learning. AddisonWesley, Reading, MA, 1989. [13] P.J. Davis, Circulant Matrices, Wiley, New York, 1979. [14] Message Passing Interface Forum, MPI: A Message-Passing Interface Standard, March, 1994.
CON NOTAT~ON
EI_SEVIER
Applied Mathematics and Computation 97 (1998) 1 9 7 2 0 8
A genetic algorithm for circulant Euclidean distance matrices Wilson Rivera-Gallego 1 NSF Engineering Research Center for Computational Field Simulation, Mississippi State University, Mississippi State, MS 39762, USA
Abstract
This paper presents a fast genetic algorithm to determine three-dimensional configurations of points that generate circulant Euclidean Distance Matrices (EDMs). A parallel implementation is possible by using the message passing interface (MPI) standard. In addition, theoretical results about the polyhedral structure of both the cone of circulant symmetric positive semidefinite matrices and the cone of circulant EDMs are introduced. © 1998 Elsevier Science Inc. All rights reserved.
1. Introduction
The first known references of Euclidean Distance Matrices (EDMs) appeared in multi-dimensional scaling [1]. Members of a multi-dimensional sample are represented by points in some geometrical space, usually Euclidean. It is possible to build up arrays of numbers, called ordering diagrams, where pairs of close points usually represent similar samples and faraway points represent dissimilar samples. The goal is to use this dissimilarity matrix in order to describe some characteristics of the sample. For example, in psychometric each element in the dissimilarity matrix may represent the dissimilarity of a pair of stimuli. It is required to represent the set of stimuli by a configuration of points in a low-dimensional Euclidean space [2]. However, it is in molecular conformation [3] where stronger applications of EDMs are found. For example, interatomic distances are measured by using
1 E-mail: [email protected]. 0096-3003/98/$ - see front matter © 1998 Elsevier Science Inc, All rights reserved. PII: S 0 0 9 6 - 3 0 0 3 ( 9 7 ) 1 0 1 43-6
198
w. Rivera-GallegoI Appl. Math. Comput. 97 (1998) 197-208
nuclear magnetic resonance (NMR) and they are located in a rectangular matrix. It is necessary to determine the possible three-dimensional configurations of points which would generate the matrix [4]. Also, it is interesting to observe the structure of the configurations from the properties of the matrix [5]. Schoenberg [6] found a natural parametrization of the EDMs cone from the symmetric positive semidefinite matrices. Gower [7,8] described some properties of the EDMs and pointed out some connections between these properties and the configurations of points that generate the matrices. Critchley [9] worked on the spectral decomposition of the EDMs. Hayden et al. [5] obtained Gower's results, but in a simple way offering us a more geometrical view of the EDMs cone. Rivera [10] studied some special structures of EDMs. Genetic Algorithms (GAs) are search algorithms based on the genetic processes observed in natural evolution. Theoretical concepts were first introduced by Holland [11] and posteriorly well described by Golberg [12]. Basically, GAs work with an initial population of individuals. Each individual (chromosome) consists of a string of characters (genes) and represents a possible solution in a search space. A fitness value is associated with each individual. By using genetic operators such as selection, crossover (mating between individuals), mutation (random modification on genes), offsprings with better fitness values are generated. This process keeps on in order to find an optimal solution to a particular problem. Section 2 introduces the definitions. In Section 3, theoretical results are presented. In Section 4, an outline of genetic algorithms is given and some details of its implementation for circulant EDMs are discussed. Section 5 deals with discussion of practical results. Finally, conclusions are given in Section 6. 2. Euclidean distance matrices
If there exist n points, pl,P2,... ,p,, in ~k such that [[Pi __pj[[2 = d2 then D = [d2.] E ~,×n is an EDM. The minimum value of k is called embedding dimension. A matrix symmetric A E R n×" is positive semidefinite if x'rAx >10 for all x E R n, x ¢ 0. In particular, if it is possible to write a matrix A c R "×" as a product A -----CC T then A is symmetric positive semidefinite. A matrix C E Rk×k such that the ith row corresponds to the coordinates of the point p~ is called Coordinates Matrix. Note that given a configuration of n points Pl,P2,... ,p, in Rk and the corresponding coordinates matrix C then the matrix B = CC T is symmetric positive semidefinite. Also
bij =pXpj and in particular
for i,j = 1 , 2 , . . . , n
W. Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197-208
b,i
= Ilp~ll 2
199
for i = 1 , 2 , . . . , n .
2.1. P a r a m e t r i z a t i o n
of EDM
Let us define: f2n = {B E Rn×n:B is a Symmetric positive semidefinite matrix}, An = {D E Rn×n:D is an Euclidean Distance Matrix}.
(1)
It is clear that d~ = IlPillz + [IPAI2 - 2p•pj - - bii -~- bjj - 2bij.
(2)
Because B = [bij] is a positive semidefinite matrix, the above expression is a natural parametrization of the EDMs ~c : (2 n --+ A n ,
[l£(n)]ij -~- bii -~- bjj - 2bij.
(3)
The map x can be expressed in matrix form by x ( B ) = b e ~ + e b T - 2B,
(4)
where b = diag(B) and e is the vector whose components are all equal to 1. Each E D M corresponds to different positive semidefinite matrices depending on the coordinates system. Thus, the function x is not injective. In order to get an inverse transformation of x, it is necessary to define the map Zs : An ---, On, z s ( O ) = - ½ ( I - e s T ) O ( I -- seT),
SE~"
(5)
such t h a t s T e = l ,
where s is a vector in R n which fixes the coordinates system. In case that s = e / n , the origin of the coordinates system corresponds to the centroid of the configuration of points [5,10].
2.2. C i r c u l a n t m a t r i c e s
A matrix A = [aij] E C n×n is called C i r c u l a n t M a t r i x if aij = ak
j - i + 1 - k (mod n).
We denote a circulant matrix by A = circ(al, a 2 , . . . , an).
W. Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197-208
200
Let F be the Fourier M a t r i x given by 1
f ~
1
1
1
...
1
1
w
w2
...
w(n-1)
1.
W2 . .
W4
"'"
w2(.-1)
1
w (~-1)
w 2(n-l/
...
w(n-l)(,-l)
.
(6)
where w = cos (2n/n) + i sin (2~/n). As pointed out in [13], a circulant matrix A E C "×" can be expressed as C = FAF*, where A = d i a g o n a l ( 2 1 , 2 2 , . . . , 2,), 2j = Pp(w/-'),
(7)
Pp(z) = cl + CzZ ~- " " " q- CnZ n - l .
This result and the s y m m e t r y of circulant E D M s can be used to generate a nice representation of this set o f matrices. We point out that the circulant distance matrices set is a polyhedral cone whose extreme points can be identified easily.
3. T h e o r e t i c a l results T h e o r e m 3.1. The cone o f circulant symmetric positive semidefinite matrices o f
order n is a polyhedral cone with extreme points given by
Aj=~T
j=I,...,L,
(8)
where
= ~fi--[cos 2 ( i - 1 ) ( j - 1 ) n n Vnl
sin 2 ( i - 1 ) ( j - 1 ) r c ] n
E ~nx2,
(9)
/
L = [n/2J i f n i s o d d a n d L = [n/2J + 1 i f n iseven. • = 1 i f j = [n/2J + 1, and = 2 otherwise.
Proof. Given a matrix circulant and symmetric C = circ(cl, c2, • • •, cn), we have that Cj = Cn_j+ 2
for j = 2 , 3 , . . . , n .
Hence, 2j=2,_j+2
for j = 2 , 3 , . . . , n .
(10)
W. Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197-208
201
F r o m (7) n
C----
Zj=l •'JY; ,
where fj = ~ 1 [1, w (j-x), w2~-1),... , w(,_~)0. 1)]T '
,
=
1 [ 1 , w - " - ' ) , w -2"-1> , .
. . , w_(n_X)(j_l)]
,
Thus, f f 7 = nL1[w(j_li(;_~) ]j (;,k) for i,k = 1,2, ... ,n.
(11)
Also f~_j+2 = f .
for j = 2 , 3 , . . . , n ,
f~-j+zf,*_j+2 = J~J~*
(12)
_ - -1 [w_(j_,lti_k) ] n
for i , k = 1 2, (i,k)
'
. n.
(13)
""
Because the vectors f,-j+2 and fj c o r r e s p o n d to the same eigenvalue, the matrix C can be d e c o m p o s e d as Ln/2J + 1
C = 21f, f ; + Z
2J(~f7 +f,-j+zf,*-j+2)
(14)
j=2
if n is odd, or tn/2J
C = Zlflf~ + )~kfkf; + Z 2 j ( f j f j * +fn-j+zf;-j+2)
(15)
j=2
if n is even, with k = [n/2J + 1. We define Aj = f j f j + f,-j+zf,*_j+2.
(16)
It is very easy to show that Aj = ~ / 0 f
j = 1,...,L.
(17)
On the other hand, it can be verified that if ss x is a circulant matrix o f rank 1, then there exists some a E ~ such that s = ~fl if n is odd, and s = ef~ or s = c~f, if n is even. In order to show this statement, it is sufficient to see that if s [S1,$2,... ,S,] x then for some ~ 6 R =
[si] = c~
i = 1 , 2 , . . . ,n.
202
w. Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197-208
Because ss x must have constant diagonal. Also, because C is circulant Cin:Ci+l, 1
for i = 1 , 2 , . . . , n -
1,
or equivalently SiS n =
Si+lS 1
for i = 1 , 2 , . . . , n -- 1.
By expanding the above expression, the statement can be verified.
[]
It is easy to show that B is a circulant positive semidefinite matrix if and only if x(B) is a circulant EDM. In addition, K is a linear transformation. As a consequence, we have the following corollary.
Corollary 3.2. The circulant E D M s set is a polyhedral cone whose extreme points are images f r o m the extreme points o f the circulant symmetric positive semidefinite matrices under the application ~. 3.1. Comments
In practice, it may be useful to determine substructures with special properties within molecular configurations. In particular, a subset of points within a configuration may generate a nearly circulant EDM. Thus, the characterization of the extreme points of the circulant symmetric positive semidefinite matrices becomes useful in order to determine the three-dimensional coordinates of the points. Our G A implementation for circulant EDMs makes good use of this characterization. Notice that, because 21flf; is a constant matrix and x(21flf~) = O, f l f ; does not have any effect on the corresponding EDM. Thus, for n odd the rank of the extreme points is 2 and for n even the rank is also 2 except for f~.fk* which has unitary rank. Hence, three-dimensional configurations that generate circulant E D M only there exist for an even number of points. Moreover, if we are looking for this sort of three-dimensional configurations, we only need to find a combination of one of the eigenvalues 2j for j = 2 , . . . , Ln/2J and 2t~/2j+l.
4. Genetic algorithms and circulant euclidean distance matrices A general G A is as follows: Generate r a n d o m l y an initial population. determine fitness of individuals. do
( S e l e c t i o n operator. C r o s s o v e r operator.
203
Rivera-GallegolAppL Math. Comput. 97(19~)197-2~
Mutation operator. Determine fitness of individuals. while (not convergence)
Individuals are selected from the population to mating processes by a selection operator. The probability of an individual being selected is proportional to its fitness. Offsprings are created by combining genes from the parents. This is made by crossover operators. For example, in one point crossover the parents' chromosomes are cut at a random crossover point and the sections are exchanged. In a uniform crossover each gene in the offspring is generated by copying randomly the gene from one or the other parent. In addition, the mutation operator changes genes randomly selected. Mutation is used to recover information and keep diversity in the population. Each individual is represented by the structure typedef struct
( int chrom[ARRAY]; double Fitness, values[2]; int parentl, parent2; ) individual;
Genetic information is stored into the binary array chrom. Recall that, in order to find a three-dimensional configuration, it is only necessary to know two eigenvalues: one of the eigenvalues 2j for j = 2 , . . . , [n/2J and the eigenvalue J,Ln/eJ+l. The decode function translates the information in chrom by dividing it into two parts, each one corresponding to one the eigenvalues, and stores the results into the array values. The length of the array chrom depends on the parameter called precision. The genetic information at each part, corresponding to an eigenvalue, is as follows: 2Precision
.
.
.
20
.
.
.
2Precision-1
Thus, the length of chrom is 4 × precision. Notice that the eigenvalues must be positive values. Hence, it is not necessary to worry about the sign. Each individual has a fitness value corresponding to the Frobenius norm of the difference between the circulant EDM generated by the individual and the input matrix. A tolerance parameter regulates the generation of individuals at each step. Usually, the initial value of tolerance is 2 Maxgen where Maxgen is the maximum number of allowed generations. At each generation, the tolerance is reduced to
204
w. Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197-208
the half guaranteeing those better individuals are selected. Thus, the G A stops when the Maxgen value is reached or after 10000 intents to select an individual that satisfies the tolerance value. Because there are Ln/2J feasible three-dimensional configurations, each one corresponding to a combination of one 2; for j = 2 , . . . , [n/2J and /~[n/2J+l the G A has a natural parallelization. The parallelized version is build up by using the Message Passing Interface (MPI) [14]. M P I has the particularity to be portable to different platforms.
5. Practical results
The code has been run on a SGI Onys Reality Engine Server with 12 processors. Different problems have been tested and the G A has shown to be a very robust algorithm. For example, by using crossover probability = 0.6, mutation
Table 1 Points
Time (min)
10
1.24
20 50 100
3.12 5.75 7.23
0.5-
-0.5.
i
-1
-1
Fig, 1.22 = 2.0, 25 = 4.125, + (real points), * (approximations).
205
w. Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197 208
probability = 0.05, Maxgen = 20, population size = 20 a n d t a k i n g different configurations with different n u m b e r o f points, the results are given in T a b l e 1. T h e time is the average a m o n g the times f r o m all o f the cases tested for a specific n u m b e r o f points. T h e time spent b y the a l g o r i t h m d e p e n d s on the c o m p l e x i t y o f the p r o b l e m . In p a r t i c u l a r , h o w m u c h closely circulant is the inp u t matrix.
E x a m p l e 5.1. F o r n = 8 the eigenvalues 2.2 = 2.0 a n d 25 = 4.0 generate the following c i r c u l a n t E D M : 0
2.2929
1.0000
3.7071
2.0000
3.7071
1.0000
2.2929'
2.2929
0
2.2929
1.0000
3.7071
2.0000
3.7071
1.0000
1.0000
2.2929
0
2.2929
1.0000
3.7071
2.0000
3.7071
3.7071
1.0000
2.2929
0
2.2929
1.0000
3.7071
2.0000
2.0000
3.7071
1.0000
2.2929
0
2.2929
1.0000
3.7071
3.7071
2.0000
3.7071
1.0000
2.2929
0
2.2929
1.0000
1.0000
3.7071
2.0000
3.7071
1.0000
2.2929
0
2.2929
2.2929
1.0000
3.7071
2.0000
3.7071
1.0000
2.2929
0
+
+
0.5. +
4-
-0.5,
+ -1:
-0.5 ~ - 0 . 5 -1 -1
Fig. 2. 23 = 1.75, 25 = 6.625, + (real points), * (approximations).
(18)
W. Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197-208
206
O.
~
+
1 "10 . 5 ~
+
~ ~
-0.5 " ~ . ~ ' ~ - 0 . 5 -1 -1
Fig. 3.24 = 1.875, 25 = 6.625, + (real points), * (approximations).
We introduce errors in the measure o f the distances. Thus, our input matrix is 0 2.2729 0.9600 3.6671 1.9400 3.6771 1.0300 2.2229
2.2729 0 2.2629 0.9700 3.7371 1.9700 3.657l 0.9200
0.9600 2.2629 0 2.3129 1.0500 3.7571 1.9800 3.6371
3.6671 0.9700 2.3129 0 2.2429 0.9200 3.4071 2.0800
1.9400 3.7371 1.0500 2.2429 0 2.1929 0.9400 3.6371
3.6771 1.9700 3.7571 0.9200 2.1929 0 2.3729 0.9000
1.0300 2.2229 3.6571 0.9200 1.9800 3.6371 3.4071 2.0800 0.9400 3.6371 2.3729 0.9000 0 2.2229 2.2229 0
(19)
R u n n i n g our algorithm with a population o f 20 individuals and Maxgen = 20, we obtain three configurations. The first one (Fig. 1) corresponds to 22 = 2.0, 25 = 4 . 1 2 5 and fitness=0.81. The second one (Fig. 2) corresponds to 23 = 1.75, 25 = 6.625 and fitness = 10.29. The last one (Fig. 3) is that corresponding to 24 = 1.875, 25 = 6.625 and fitness = 10.62. It is interesting to notice the structure o f three-dimensional configurations that generate circulant E D M s . Fig. 4 shows an example for a configuration with 20 points.
6. Conclusions A fast genetic algorithm to determine three-dimensional configurations that generate circulant E D M s has been developed. A parallel version o f the algo-
W. Rivera-Gallego I Appl. Math. Comput. 97 (1998) 197-208 "5
207
"8
0.5, "~9 ~7
~t3
4t5
a4 ~B
~o
-0.5> 0.5 .5
-0,5 .0.5
Fig. 4. n = 20, 22 = 2.0, and 211 = 4.0.
rithm was implemented by using MPI standard. Interesting theoretical results regarding the polyhedral structure of the cone of circulant symmetric positive semidefinite matrices and the cone of circulant EDMs were introduced. Recognition of substructures within molecular configurations is very important in order to determine more complex structures. The results obtained in this work may be useful to explore other configuration such as Toeplitz structures. Moreover, our genetic algorithm can be easily modified in order to solve other problems such as matrix completion. It would be interesting to use the knowledge from the structure of the set of EDMs to generate domain knowledge genetic algorithms. References [1] J. De Leuuw, W. Heiser, Theory of Multidimensional Scaling, in: P.R. Krishnaiah (Ed.), Handbook of Statistics, North-Holland, Amsterdam, 1982. [2] G. Young, A.S. Householder, Discussion of a set of points in terms of their mutual distances, Psychometrika 3 (1938) 19-22. [3] G.M. Crippen, T.F. Havel, Distance Geometry and Molecular Conformation, Wiley, New York, 1988. [4] W. Glunt, T.L. Hayden, M. Raydan, Molecular conformations from distance matrices, J. Comput. Chem. 1 (1993) 114-120. [5] T.L. Hayden, J. Wells, W. Liu, P. Tarazaga, The cone of distance matrices, Linear Algebra Appl. 144 (1991) 153-169. [6] I. Schoenberg, Remarks to Maurice Frechet's Article, Sur la definition axiomatique d' une classe d'espaces vectoriels distanci~s applicables vectoriellement sur L'espace de Hilbert, Ann. Math. 36 (1935) 724-732.
208
I,K Rivera-Gallego / Appl. Math. Comput. 97 (1998) 197-208
[7] J.C. Gower, Euclidean distance geometry, Math. Scientist 7 (1982) 1 14. [8] J.C. Gower, Properties of euclidean and non-euclidean distance matrices, Linear Algebra Appl. 67 (1985) 81-97. [9] F. Critchley, On certain linear mappings between inner product and squared-distance matrices, Linear Algebra Appl. 105 (1988) 91 107. [10] W. Rivera, Matrices de distancia con estructuras especiales. M.Sc. Thesis, University of Puerto Rico, Mayaguez, 1994. [11] J.H. Holland, Adaptation of Natural and Artificial Systems. University of Michigan Press, Ann Arbor, MI, 1975. [12] D.E. Golberg, Genetic Algorithms in Search, Optimization and Machine Learning. AddisonWesley, Reading, MA, 1989. [13] P.J. Davis, Circulant Matrices, Wiley, New York, 1979. [14] Message Passing Interface Forum, MPI: A Message-Passing Interface Standard, March, 1994.