Finding Relations in Polynomial Time

1 Introduction

Given a set of observations of several variables, to find efficiently parsimonious relations that they satisfy is a basic problem in the study of scientific discovery [Lan­ gley et al., 1987], [Simon et al., 1997] and, possibly, in its enhancement. More generally, one may seek relations between a subset of the given set of variables. Enumer­ ating all subsets and considering them in turn reduces this problem to the former one, but entails a combinato­ rial explosion. The purpose of this paper is to show that, for large classes of relations, the latter problem can be solved in polynomial time. More precisely, a polynomial algorithm is proposed for finding a basis of all affine rela­ tions between the variables (or powers, ratios, products, or logarithms of them). This algorithm uses basic results of linear algebra, also exploited in data analysis, and in particular in principal component analysis [Pontier et al.,

Pierre Hansen GERAD and

Ecole des Hautes Etudes Commerciales 3000, Chemin de la cote-Ste-Catherine Montreal H3T 2A7, Quebec, Canada e-mail: piP.rreh @crt. umontreal. ca

1990]. However, the purpose there, i.e., to explain differ­ ences between observations, is the reverse of ours, i.e., to find what properties the observations have in common.

The method is described in the next section. In Sec­ tion 3 it is applied to several classical discoveries al­ ready studied previously by the Bacon programs [Lang­ ley et al., 1987]. In Section 4, results of its embedding in the AutoGraphiX (AGX) system for computer-aided graph theory are presented. The resulting entirely au­ tomated system discovers five novel relations (or conjec­ tures) in graph theory, two of which pertain to math­ ematical chemistry. Three of these conjectures involve five invariants, which is more than in most propositions of graph theory. Conclusions are drawn in Section 5. Proofs of two among the five conjectures are given in the Appendix.

2 Basic principles

2.1 Description of the algorithm

Consider a set of m observations of n quantitative vari­ ablesx1,x2, ... ,xn. LetXmxn = (Xij) denotethecorre­ sponding data set. Our aim is to find a basis of all affine relations (i.e., relations of the form a1x1 + a2x2 + · · · + anXn = b, where the ai and b are constants) satisfied by the observations. Let p. = ;k t X 1m denote the vector of mean values for the columns of X. Then any affine relation between columns of X becomes a linear relation (i.e., a relation of the form a~ x~ + a2x2 + · · · + a~x~ = 0) between columns of X'= X- (1mtP.), the matrix of cen­ tered values. Furthermore, all coefficients (except b) are

and thus

holds, then

which means that if the linear relation (1) holds for the columns of X', it also holds for those of V. As Vis symmetriC, it still holds for its lines. The second step of our algorithm is thus to compute V.

The third and last step is to diagonalize V (with, how­ ever, some empty lines if there are relations) . This can be done, by Gaussian elimination. In the resulting ma­ trix V' , Dim(Im(V)) lines contain non-zero terms and correspond to independent variables. The remaining n- Dim(Im(V)) lines contain only zeroes and corre­ spond to dependent variables which may be expressed as linear combinations of the independent ones. These relations form a basis of the null-space of V. Using the initial data one can then compute the right-hand sides of the corresponding affine relations.

The case were a dependent variable is expressed as a monomial of several independent ones, instead of an affine expression is easily reduced to the previous one. Indeed, assume

where the et and f are constants. Taking logarithms, this relation becomes

So the method just described can be used to find ex­ ponents et (equal to 0 if x~ does not appear in the mono­ mial) and the initial data to find f.

The proposed algorithm is illustrated on several ex­ amples in Section 3.

Centering variables, the first step of our algorithm, thus transforms the problem of finding affine relations into the problem of finding linear ones.

Consider then the variance-covariance matrix V de­ fined by Vn xn = tX'X'. If the relation

2.2 Extension

The two classes of relations described in the previous subsection can be easily and substantially enlarged by computing additional terms from the initial data. The easiest case is to take, in addition to the given variables, all squares and products of pairs of them. Then one uses n + n 2 ±n = n 2 ± 3 n variables instead of n. The formulas have tfte form 2

Ratios, higher powers than two or products of three vari­ ables may also be considered. If Xj is expressed as a monomial, sums or differences of variables may be taken as additional terms.

2.3 Complexity

The first step of the method, computing centered values involves computing J.L, in O(mn) time and subtracting J.Li from each X i j, again in O(mn). The second step, comput­ ing the variance-covariance matrix V, takes 0( mn 2 ) time. The third step, diagonalizing V by Gaussian elimination requires O(n 3 ) time. In order to avoid spurious relations, one must have m > n. Thus the overall complexity is in O(mn 2 ) (or in O(n3) if one makes the reasonable assump­ tion that m is O(n)). If squares and products of pairs of variables are added the complexity rises to O(n 6 ). Such problems, with n moderate, are still solvable in a rea­ sonable time.

As an example for comparison, using linear regression allows at most one relation to be found at a time while our method finds a basis of all the relations underlying the data. Furthermore, linear regression needs one vari­ able to be explained by the others, which means that a regression must be computed for each of the n - 1 first variables to be explained by the next ones. To ensure that all relations are found, one cannot escape the "best-subset" criterion for each regression, thus try­ ing all combinations of variables, which has exponential complexity. Treating a 20 variables and 44 observations example with the SAS statistical software needs 10 sec­ onds while our system needs 0.0013 sec (on a Sun Ultra I, 140 MHz). Using 25 variables instead of 20 and still 44 observations, the time required by SAS increases to 4 minutes and a half while our method only needs 0.002 seconds (130 000 times faster).

3 Rediscovering known laws

3.1 Kepler's third law

Kepler's third law is that p 2 / d 3 = ct where p denotes the period of a planet's orbit and d its average distance to the sun. As all those considered in this section, this law was rediscovered by one of the Bacon programs [Langley et al., 1987]. The computing time they use is much larger than with our method. Note, however, that the aim of these programs is different from that pursued here: i.e., to understand the reasoning which led to the discovery rather than to find an efficient but computer-dependent way to discover the corresponding relation.

A generator was used to produce the following four simulated observations of p and d (taking units such that ct = 1):

or, switching to logarithms,

The corresponding variance-covariance matrix is

Diagonalization leads to

which shows that a relation of the form log(p) = log( d)ct' holds. Then using Xtog, we find ct' = 1.5, and an equiv­ alent form of the original relation, p = dl. 5 . Computing time is too small to be measured.

3.2 Ideal-gas law

The ideal gas law is PV = nRT where P is the pres­ sure (in Pascals), V the volume (in cubic meters), n the number of moles and T the temperature (in Kelvins). R = 8.32 is the universal ideal gas constant. The genera­ tor gave simulated observations of P, V, n and T. These were then converted to the older units, to simulate the situation prevalent at the time of discovery, i.e., before absolute temperatures were used. Then the pressure p is measured in atmospheres, the volume v in liters and the t temperature in Celsius degrees. The conversion rules are P = 101325p, V = v/1000 and T = t+273.15. Using 25 observations and the extended set of terms including squares and products of pairs of variables led to the re­ lation nt = - 273.15n + 12.1785pv in less than 1 ~ 0 of a second. Converting back to the SI System units gives nT = 0.120192PV, which is an equivalent form of the law presented above.

Note that this relation is even easier to find, using the SI System units, if we switch to the logarithms of the values; then an extended data set is not needed.

Sixteen observations are enough and the variance­ covariance matrix for these data is

Diagonalizing this matrix becomes

which shows that a relation of the form log(T) = log( P) + log(V) -log( n) + ct holds. From the logarithms of the data, we find ct = 0.120192. So the equivalent form T = 0.120192PVn- 1 of the ideal gas law is found instantaneously.

3.3 Newton's gravity law

Newton's Law is that F = G"}}f where m and Mare the masses of two bodies, D the distance between them, F their force of attraction and G = 6.67e- 11 the universal gravitational constant. The problem is formally very similar to the previous one, as one seeks a monomial involving three variables.

The generator provided 25 simulated observations of F m M and D, then the system obtained from their , , 2 67 11 h. h . logarithms the relations F = mM D- 6. e- w 1c 1s equivalent to the law given above. Again the derivation was instantaneous.

In addition to this result, the system gave spurious laws, i.e., m + F = m, M + F = M and D + F =D. It is easy to recognize that they are incorrect, and due to numerical errors. A more careful scaling of units before computation avoids this error.

3.4 Ohm's law

Ohms' law is that I L = r I + v where I denotes cur­ rent intensity, L inductance, r resistance and v tension. The generator provided 25 simulated observations of I, L, r and v. Using an extended data set with squares and products of pairs of variables, i.e., 14 terms (or 10 when seeking for expression of a dependent variable as a monomial, because in that case sums of twice the same variable are redundant) the relation r I = -v + I L was found in less than 1 ~ 0 of a second.

Note that using only the original observations would not have allowed to find this relation.

4 Finding new relations in graph theory

4.1 The AutoGraphiX system

Several systems have been proposed for computer-aided or automated discovery of relations in graph theory, prominent among which are Graph [Cvetkovic et al., 1981]; [Cvetkovic and Simic, 1994], and Graffiti [Fajt­ lowicz, 1988]; (Fajtlowicz, 1998]. They work with graph invariants, i.e., quantities such as the number of ver­ tices, number of edges, stability number, chromatic num­ ber, radius, diameter, etc., (see [Berge, 1973] for unde­ fined terms) which do not depend on labeling of ver­ tices or edges. Recently, a new system, called Au­ toGraphiX (AGX) has been developed [Caporossi and Hansen, 1997]. Its basic idea is to consider the problem of determining extremal graphs for some graph invari­ ant, possibly subject to constraints, and related problem such as refuting conjectures, as problems of combinato­ rial optimization defined on an infinite family of graphs. Then, powerful heuristics may be brought to bear to find presumably extremal graphs, of moderate size, from which much information may be gathered. More pre­ cisely, AGX addresses the following problems:

• (a) Find a graph satisfying given constraints. Let i1 (G), i2(G), ... , it(G) denote invariants of a graph G and P1, P2, ... , Pt the values they should take. Consider the problem

where 9n denotes the set of graphs with n vertices, and n may be varied. Any graph for which /(G) = 0 satisfies the given constraints.

• (b) Find optimal or near optimal values for an invariant io(G) subject to constraints. Consider the problem

where M is a sufficiently large constant.

• (c) Refute a conjecture, expressed as h(G) ~ g(G) where h(G) and g(G) are formulas depending on one or more invariants of G. Consider the problem

IT a graph G for which /(G) < 0 is found, the con­ jecture is refuted.

• (d) Suggest a conjecture. This can be done in various ways, which use parameterization on the number n of vertices and m of edges of G or other invariants. Surfaces of values for the invariant or formula un­ der study can be displayed graphically, which often suffices to obtain conjectures, and presumably ex­ tremal graphs can be studied, which gives further information.
• (e) Suggest a proof. Checking if presumably extremal graphs found are also obtained with simplified heuristics, such as e.g. local search with a single type of move may suggest ways to prove conjectures for all graphs, or for some classes of them, e.g., by showing that such a move can always be performed on any non-extremal graph and is an improving one.

The heuristics used in AGX fit into the Variable Neigh­ borhood Search metaheuristic framework [Hansen and Mladenovic, 1998]. A local search routine is applied to G, which employs several types of moves: add or re­ move an edge, displace one or two edges, add or remove a pending vertex, insert a vertex in an edge, etc. Once a locally optimal graph G1 has been found, one considers a set of neighborhoods Nl(Gl), N2(G1), . .. , Nkmax(Gl)

increasingly far from G1 and defined by addition or re­ moval of 1, 2, . . . , kmax edges. A graph G2 in the neigh­ borhood N1(G1), N2(G1), etc., is drawn at random and a local search performed leading to a graph G3 . IT G3 has a better value than G1 the search is recentered there. AGX has been applied to a series of problems of graph theory, some of which pertain to mathematical chem­ istry [Caporossi and Hansen, 1997]; [Caporossi et al., 1998a]; [Cvetkovic et al., 1998]; [Caporossi et al., 1998b]; [Caporossi et al., 1999]. 24 conjectures have been ob­ tained up to now, 11 of which have been proved and none of which has been refuted. In addition, AGX alone refuted 3 conjectures of Graffiti [ Caporossi and Hansen, 1997] and, together with some enumeration routines 6 more [Pujol, 1998]. As described, AGX is a tool for computer-assisted discovery. However, addition of a rou­ tine implementing the algorithm for finding relations de­ scribed in Section 2 transforms it into an entirely auto­ mated system. Results obtained when doing so are next presented.

4.2 A property of color-constrained trees with minimum index

A tree is a bipartite graph, whose vertices may be col­ ored in black and white such that no pair of adjacent vertices have the same color. IT the numbers of black and of white vertices are fixed the tree is color-constrained. In [Cvetkovic et al., 1998] AGX is used to study color­ constrained trees extremal with respect to their index, or largest eigenvalue of their adjacency matrix. Six con­ jectures are obtained, and four of them proved.

To explore further the extremal trees found, values of the following fifteen invariants were recorded: number n of vertices, number n1 of pending vertices, number m of edges, diameter D, radius r, stability number a (i.e., maximum number of pairwise non-adjacent vertices), av­ erage degree 8, average distance and sum of distances between pairs of vertices, energy (sum of the absolute values of the eigenvalues of the adjacency matrix), max­ imum degree, largest eigenvalue of the adjacency matrix, Hyperwiener index (sum of the squares of distances be­ tween pairs of vertices), Ran die index x (see below), and Chromatic number 'Y.

Then the algorithm of Section 2 was used to find a basis of affine relations on those invariants. In addition to well-known relations, such as m = n - 1 and 'Y = 2 which are valid for all trees, the unexpected following one was obtained

and the remaining 8 invariants are proved to be linearly independant from any of the considered ones.

This result can also be expressed as follows.

Conjecture 1. For all color-constrained trees with minimal index

Furthermore, AGX was used to see whether this con­ jecture could be extended to all trees. Minimizing and maximizing the left-hand side of (4) then led to the fol­ lowing results.

Conjecture 2. For all trees

Let La J denote the largest integer smaller or equal to a.

Conjecture 3. For all trees

Conjecture 2 is proved in the Appendix.

4.3 Two bounds for the connectivity index of chemical graphs

Chemical graphs are extensively used to represent hydro­ carbons. Vertices correspond to carbon atoms and edges to bonds between them (hydrogen atoms are usually not represented as their number and bonds are uniquely de­ fined). Consequently the maximum vertex degree is 4. Moreover, chemical graphs are connected. The connec­ tivity index (or Randic index) [Randic, 1975] has been introduced as a measure of molecular branching. It is defined as follows. Let (vi, v;) denote an edge of a graph G, d; and d; the degrees of its end-vertices. Then the weight of that edge is (di d;)- 1 1 2 and the connectivity index x( G) is the sum of the weights of all edges of G. This index is one of the most frequently used molecu­ lar structure-descriptors, and the only one to which two books [Kier and Hall, 1976]; [Kier and Hall, 1986] and numerous articles have been devoted.

In [Caporossi et al., 1998b] chemical trees with ex­ tremal connectivity index are studied, based on results of AGX, and fully characterized. Applying AGX, with the relation-finding routine, to chemical graphs with pre­ sumably minimal connectivity index, did at first give no relations. However, considering the set of local minima among these graphs gave two new relations:

Conjecture 4. For all chemical graphs

Conjecture 5. For all chemical graphs

Conjecture 4 is proved in the Appendix.

5 Concluding remarks

An algorithm for finding automatically from observa­ tions of a set of variables a basis of affine relations be­ tween these variables, or products, ratios or powers of them, or between logarithms of such terms has been presented. As it is polynomial in the number of vari­ ables considered, the combinatorial explosion encoun­ tered when each potential relation is considered in turn is avoided.

Moreover, adding a routine implementing this algo­ rithm to the system AGX led to an entirely automated system for finding relations in graph theory. Its use is il­ lustrated by the derivation of five novel conjectures, two of which pertain to mathematical chemistry. Many more applications are expected in the near future.

A Proofs of two conjectures

Theorem 1 {Conjecture 2). For all trees with stabil­ ity number a , m edges, n1 pending vertices, diameter D and radius r

Proof. Recall that in a tree D = 2r if D is even and D = 2r -1 if Dis odd. Hence, D- 2r = -(D mod (2)). Assume first that D is even. Define the star of a vertex Vi to be the set of edges incident with Vi· Clearly, stars of vertices of a stable set are disjoint. Consider a maxi­ mum stable set S with a 1 pending vertices and a2 inner vertices. Then, m ~ a1 + 2a2 = 2a - a1 ~ 2a - n1 and a~ Hm +nt).

Assume then that D is odd. Consider a path P of length D. Either only one end vertex of P belongs to S or there is a pair of consecutive vertices Vi and v; of P which do not belong to S, and the edge (vi, v;) does not belong to the star of any vertex of S. In both cases m ~ 2a- n1 + 1 and a ~ Hm + n1 - 1). D

Theorem 2 (Conjecture 4). For all chemical graphs G with connectivity index x( G), n1 pending vertices and m edges

Proof. Let ni, i = 1, ... , 4 denote the numbers of ver­ tices of degree i of G and Yii, i, j = 1, . . . , 4; i ~ j denote the numbers of edges of G with end vertices of degree i and j. Then consider the 6 equations:

obtained by counting vertices and degrees. These equations are clearly independent. Solving for n11 n2, n3, n4, Y14 and Y44 and substituting Y14 and Y44 in the objective function

leads to the relation

As Y12, Yl3 , Y22, Y23 , Y24, Y33 and Y34 are non-negative and have positive coefficients in {11) the result follows. 0

Conjecture 5 can be proved in a similar way.

Aknowledgments

The authors thank Franc;ois Bellavance and Robert Latour from Ecole des Hautes Etudes Commerciales de Montreal for their advices on statistical meth­ ods, and Ivan Gutman from university of Kragujevac (Yougoslavia) for discussions about AGX and its re­ sults. This work was supported by NSERC grant OG­ POO 39682, FCAR grant 32EQ 1048 and a grant from CETAI-HEC.