Valeri Beck
email: v.beck@alfa-carbon.de

The origin of the Solar System in the field of a standing sound wave

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The classical theory of the origin of the Solar System is based on the Kant-Laplace nebular hypothesis [START_REF] Kant | Kant's Cosmology, Universal Natural History and Theories of the Heavens[END_REF][START_REF] De Laplace | Expositions du système du monde[END_REF] suggesting that the Sun and the planets condensed out of a spinning nebula of gas and dust. During the condensation process, the spin of the nebula accelerated (the 'pirouette' effect), and the resulting distribution of angular momentum caused it to form a disc. The centre of the nebula evolved into a hot, highly compressed gas region -the protosun. Concentration and coalescence of dust particles in the remainder of the spinning disc led to the formation of planets orbiting the Sun.

This theory accounts for the general nature of the origin of the Solar System, but it cannot explain many of the observed facts. One problem is the angular momentum distribution: the Sun has more than 99.8% of the entire system mass, but only about 0.5% of the total angular momentum, with the remaining 99.5% residing in the orbiting planets. The hypothesis also fails to explain the 7 • tilt of the Sun's equatorial plane relative to the average orbital plane of the planets. Another serious problem is the distinction between small solid-surface inner planets and outer gas giants. Moreover, the observed regularity in planetary distances from the Sun -the so-called Titius-Bode law -has no explanation.

It is difficult to explain the existence of the Oort cloud beyond the Solar System's planets. It consists of trillions of small objects composed of dust and water, ammonia and methane ice and it is believed that these objects were scattered outwards by the gas giants at the planetary formation stage and then acquired distant circular orbits (out to about one light year) as a result of gravitational forces due to nearby stars. Such an Oort cloud emergence scenario seems very unlikely for such a large number of bodies.

Neptune is the most distant gas planet. Based on the decreasing series of giant planet masses -Jupiter: 318 earth masses (M ⊕ ); Saturn: 95.3 M ⊕ ; Uranus: 14.5 M ⊕ -we could expect Neptune's mass to be several times smaller than that of Uranus. Such a mass distribution can be explained by decreasing density of the gas nebula from the centre to the periphery, so that each planet has less gaseous substance than the previous one. In reality, Neptune has a mass of 17.5 M ⊕ , which is greater than the mass of Uranus.

Authors such as Chamberlin (1901), Moulton (1905), Schmidt (1944), von Weizsaecker (1944), McCrea (1960), Woolfson (1964) and Safronov (1972) [START_REF] Chamberlin | On a Possible Function of Disruptive Approach in the Formation of Meteorites, Comets, and Nebulae[END_REF][START_REF] Moulton | On the Evolution of the Solar System[END_REF][START_REF] Schmidt | Meteoritic theory of the origin of the Earth and planets[END_REF][START_REF] Von Weizsäcker | Über die Entstehung des Planetensystems[END_REF][START_REF] Mccrea | The origin of the solar system[END_REF][START_REF] Woolfson | A capture theory of the origin of the solar system[END_REF][START_REF] Safronov | Evolution of the Protoplanetary Cloud and Formation of the Earth and Planets[END_REF] have offered a variety of scenarios for the Solar System's formation. However, none of the existing theories is able to give a comprehensive picture of the emergence and development of the planetary system. A large number of recently discovered exoplanets (1883 confirmed planets by August 2015, NASA Exoplanet Archive1 ) allows us to estimate the number of planets in our galaxy as many billions, so that a comprehensive understanding of the processes leading to the emergence of planetary systems becomes critically important. In this article we present a new conceptual model of the Solar System's formation from a pre-solar gas and dust nebula in the field of a standing sound wave.

THEORETICAL CONSTRUCTION 2.1 A spherical standing sound wave

It is generally accepted that during a sufficiently strong process of heating, gravitational contraction in the centre of the pre-solar nebula started thermonuclear fusion of hydrogen into helium. A large amount of energy was emitted and radiation pressure prevented the further contraction of the gas nebula. However thermonuclear fusion could not start quietly, the process resulted in an explosion that caused a spherical shock wave originating from the central region. The specific explosion mechanism accompanied by the thermonuclear fusion processes is beyond the scope of this article; we note here only that the power of the explosion would have been large enough for the shock wave to spread all over the pre-solar nebula. A long time observed active processes, such as jets and outflows associated with star formation (e.g. [START_REF] Bally | Twin Herbig-Haro jets and molecular outflows in L1228[END_REF][START_REF] Bally | Observations of Jets and Outflows from Young Stars[END_REF][START_REF] Eislöffel | Parsec-scale molecular H 2 outflows from young stars[END_REF][START_REF] Frank | Jets and outflows from star to cloud: Observations confront theory[END_REF][START_REF] Busquet | The CHESS survey of the L1157-B1 bow-shock: high and low excitation water vapor[END_REF][START_REF] Nisini | Mapping water in protostellar outflows with Herschel-PACS and HIFI observations of L1448-C[END_REF][START_REF] Santangelo | Herschel-PACS observations of shocked gas associated with the jets of L1448 and L1157[END_REF][START_REF] Safron | HOPS 383: An Outbursting Class 0 Protostar in Orion[END_REF][START_REF] Stanke | Observations of molecular jets in Orion A[END_REF] ), suggests that the explosions fundamentally inherent in the birth of stars and can start the emergence of planetary systems as described in our model.

As the wave propagated, the gas particles in the pre-solar nebula oscillated radially, for two reasons: first, the difference in gas pressure; and secondly, the gravitational pull toward the centre of the nebula. The second factor began to play a crucial role at large distances from the central attracting mass concentrations: at some point, the accelerated gas particles at the bow shock would no longer return under the influence of gravity. Gas density at this distance from the centre would have been so low that the pressure difference could no longer cause the return movement of gas particles. The peripheral part separated, and the pre-solar nebula nebula was thus divided into a central spherical region and an expanding spherical shell of gas (Fig. 1).

Gas particles at the boundary of the central spherical nebula fell under gravity towards the centre then stopped when the pressure of the lower gas layers exceeded the gravitational pull, and began to move in the opposite direction. As a result of the interaction between gravity and contracting gas backpressure, the boundary of the spherical nebula began to oscillate radially and a reverse spherical wave propagated from the periphery towards the centre. The sound wave would have appeared to reflect from the gravitational pause -a spherical shell where the gas particles' velocities were not high enough to overcome the gravitational pull from the centre of the nebula and the gas particles that had moved outwards returned due to oscillations; that is, their velocities were lower than the escape velocity: where G is the gravitational constant, M is the mass of the attraction centre and r is the distance of the gravitational pause from the centre. The centre of the pre-solar nebula became quiescent after the first massive explosion: the central part, having been heated due to the emission of large amounts of energy, expanded, since the gravitational field was not yet intense enough to resist the much increased gas pressure, and the thermonuclear fusion of hydrogen combustion diminished. A large mass of gas was ejected into the surrounding space while the compact region of compressed and heated gas -the protosun -remained at the nebula's core. Several hundred years after this, the wave returning from the boundary of the spherical nebula reached the protosun, concentrated in the centre, and then began to propagate towards the periphery again: a rapid pressure increase resulted in a dramatic rise in temperature at the centre of the protosun, generating a new hydrogen explosion, much weaker than the first one but still strong enough to give extra energy to the wave reflected from the central region. The reflected wave travelled all the way from the centre to the periphery of the spherical nebula, was reflected from the gravitational pause and then returned to the centre again, causing another explosion. This process, repeated several times, eventually established regular oscillations: the wave propagated from the centre to the edge, was reflected from the boundary of the spherical nebula, and returned, causing another explosion that compensated for wave energy loss. Thus, the wave caused explosions while acquiring the energy it needed, establishing a self-sustaining process whose period, defined by the free oscillations of the spherical nebula boundaries, equalled hundreds of years.

ν e = 2GM r , (1)
The acoustic radiation pressure prevented gravitational contraction of the pre-solar nebula and compensated for deviations from the spherical shape that might result, for example, from the gas turbulence. Waves travelling from the centre and from the periphery interfered with each other to form a giant spherical 3/13 acoustic cavity resonator the size of the modern Solar System, including the Kuiper belt and the scattereddisc (Fig. 2). This resonator contained a standing wave with nodes and antinodes, the number of which would have been at least 10 for the pre-solar system: 8 planets + the asteroid belt + the scattered-disc objects. This article does not include a calculation of exact distances between the antinodes of the standing wave; we may simply state that the wavelength increases with increasing distance from the centre of attraction. Section [3.4] outlines some regularities that allow us to make some observations about the physical reasons underlying planetary positions in the Solar System.

Dust concentration in the antinodes

In our model, the pre-solar system can keep 'sounding' for many millions of years, as the periodic central explosions significantly retard the gravitational contraction of the protosun and the acoustic radiation pressure stabilizes the gas nebula, preventing its collapse. The dust present in the pre-solar nebula gradually concentrates in the antinodes of the acoustic oscillating system. Apart from gas viscosity, the process of solid particle concentration in the standing wave also relies on attraction from the large gas masses that periodically emerge in the antinodes of the standing wave. This attractive force makes the dust particles move towards the gas clusters in the antinodes and collide with each other, causing redistribution of their velocity vectors in such a way that the dust particles come to rest in the centre of the antinode.

During the gas expansion phase, dust particles in each antinode are gravitationally attracted towards two neighbouring antinodes. The two gravitational sources in these antinodes largely cancel each other 4/13 out and, moreover, are remote, so that the dust particles in the gas expansion phase mostly remain in the same place. The next phase of compression attracts new dust particles to the antinode, which also come to rest due to collisions with each other and to the viscosity of the compressed gas.

Migration of solid matter and formation of rings

A certain period of time after the standing wave was established, most of the dust would have been concentrated in the regions of spherical antinodes. Due to increased dust concentration, the number of particle collisions significantly increased, causing their coalescence and increase in size and mass. The gravitational field within the antinodes and the centrifugal force together caused dust to concentrate and form equatorial rings, with radii corresponding to the antinodes of the standing sound wave (Fig. 3). The model predicts that the increase in mass of the solid particles clusters also makes it increasingly difficult for them to stay at the antinodes of the standing sound wave as at the same time that they are affected by attraction from the central gas masses, principally the protosun. The largest clusters of dust fall toward the centre; this results in migration of solid matter from the outer regions of the pre-solar nebula towards the internal zones, while the centrifugal force does not allow the clusters to fall directly onto the protosun. The gravitational attraction of the gas masses in the antinodes directs the matter falling towards the rings forming in the equatorial plane, where the solid matter comes to rest. A significant amount of 5/13 solid matter collected in the region of the orbits of Venus and Earth, as it was concentrated here from most of the pre-solar nebula volume.

Fc Fg

This dust, gathered in compact rings, is thousands of times more concentrated here than in the primary gas and dust nebula, greatly accelerating the coalescence of particles and leading to the emergence of increasingly massive dust clusters. Agglomerated clusters of dust in narrow spinning rings move in almost identical circular orbits around the protosun and have very small relative velocities. Since their collisions do not result in fragmentation, they form massive planetesimals relatively quickly, which in turn agglomerate to form planets and their satellites. In this model, are no catastrophic collisions between the nascent planets, their satellites or other massive objects; the planets keep the circular orbits and planes of the equatorial dust rings.

Birth of the Sun

Gravitational contraction of the protosun, which had been significantly slowed as a result of regular explosions in the centre, still continued, and millions of years after the establishment of the standing sound wave gas temperature at the centre of the protosun became sufficiently high for thermonuclear fusion of hydrogen into helium to continue between explosions, resulting in the birth of a new star -the Sun. The newborn Sun stopped augmenting the sound wave by periodic explosions, and the standing wave diminished. The gas shell previously supported by acoustic pressure began to shrink, forming gas-giant planets around already existing solid nuclei. Jupiter acquired the greatest share of mass, as it was located in the region of highest gas density. All the other gas planets were situated in a lower density environment and thus obtained smaller masses. Jupiter also received large amounts of gas from the inner-planet region, blown outwards by the strong solar wind in the first few millions of years after the birth of the Sun.

The Oort cloud

The spherical expanding shell that separated from the spherical part of the nebula during the passage of the shock wave from the first explosion (see Fig. 1) moved faster than the escape velocity ν e and could not return to the centre as the pull of gravity from the central masses was too weak at such distance. With its expansion, this spherical shell accumulated increasing quantities of highly rarefied gas from the primary nebula at its inner edge, while its expansion rate gradually slowed and eventually stopped at a distance of about 1 light year from the centre. This formed a giant spherical region containing dust, ice and frozen gases particles in addition to gaseous hydrogen and helium. Over time the gas component of the shell dissipated in space, while solid particles were concentrated in increasingly large chunks of ice and dust -the cometary bodies -to form the Oort cloud, [START_REF] Oort | The Structure of the Cloud of Comets Surrounding the Solar System and a Hypothesis Concerning its Origin[END_REF] which has a weak gravitational connection with the central part of the system.

Neptune, the Kuiper belt and scattered-disc objects

After expansion of the gas shell ceased in the Oort cloud, a weak reverse wave formed within the shell and moved towards the central spherical nebula of gas, where planets had already emerged. Hundreds of thousands or millions of years later the reverse wave from the Oort cloud collided with the outer boundary of the central spherical nebula, causing a redistribution of matter on the edge of the Solar System. Large gas and ice masses from the outer antinodes of the standing sound wave were shifted, and the planet Neptune formed a little closer to the Sun relative the original position of the 9th antinode, while its mass increased several times by capture of gas and ice from the reverse wave from the Oort cloud. Scattered-disc objects such as the minor planet Eris acquired highly elongated orbits, as there was a long period under the influence of gravitational pull from the gas masses transported by the reverse wave from the Oort cloud. Minor planets in the Kuiper belt also gained significant eccentricity. The orbits of the low-mass objects in 6/ 13 the Kuiper belt and in the scattered disc have been changed so dramatically that they had become to reach Neptune and Uranus, and even began to go into the central region of the solar system and possibly caused the late heavy bombardment of the inner planets, in which was also attended cometary bodies from the Oort cloud, that have had unstable orbits in the first time after the formation.

The periphery processes were very slow, developing over many millions of years, and were relatively weak in their effect on the central region of the pre-solar system. The reverse wave from the Oort cloud only caused the formation of the Kuiper belt, an offset of Neptune's formation and its mass increase. The rest of the Solar System was and is still affected by the region of the Oort cloud only through comets.

DISCUSSION

Terrestrial and giant planets

The difference in chemical composition of the inner and outer planets is explained in our model by the migration of solid matter from the outer spherical shells to the internal ones (see Fig. 3). Before concentrating in stable spinning dust rings, the heavy chemical elements, such as iron and silicate travelled farther towards the centre of the system and were incorporated by the inner planets, with Earth and Venus having more matter, and Mercury and Mars having less. The composition of the closest to the Sun Mercury at the same time has obtained a significant amount of iron. The region between Mars and Jupiter had insufficient solid matter left for a proper planet, so only the minor planets of the asteroid belt were formed. Lighter and more volatile chemical compounds, such as water, methane or ammonia remained in significant quantities in the colder regions beyond the asteroid belt, and became the compositional basis of the giant planets.

Angular momentum distribution and dust structure

If we assume that the first powerful explosion in the centre of the protosun was not absolutely symmetrical, slight asymmetry in the explosion led to a redistribution of the angular momentum in the pre-solar nebula, with the less massive peripheral zone of the pre-solar nebula beginning to rotate faster. In addition, an initial asymmetry could produce the observed 7 • equatorial plane tilt relative to the Sun's plane of rotation.

Currently the Sun has more than 99.8% of the entire Solar System mass but only about 0.5% of the total angular momentum. If the rotational kinetic energy were evenly distributed over the original pre-solar nebula, the resulting rotation (taking into account the pirouette effect) would be thousands or millions of times slower than the current rotation of the planets. The angular momentum was initially distributed evenly in the pre-solar nebula, and this means that such weak rotation cannot result in the process of dust disc formation as described by the nebular hypothesis.

However, these discs are visible in images of some young stars (Fig. 4), and in some cases their internal structure can be identified. If a spherical dust cluster such as described in this article were observed from one side, it would be impossible to observe any individual spheres as they would overlap in the line of sight (Fig. 5).

The stage at which a significant amount of dust has already accumulated in the equatorial plane would give such a system the appearance of a disc, despite the fact that there are significant gaps between individual dust clusters. Fig. 6 shows images from the NASA Hubble Space Telescope, showing the complex three-dimensional structure of the dust around young stars that are similar to the structure, which should be observed according to the proposed model; for more images see also: Schneider et al. [START_REF] Schneider | Probing for Exoplanets Hiding in Dusty Debris Disks: Disk Imaging, Characterization, and Exploration with HST/STIS Multi-Roll Coronagraphy[END_REF] On these images we see the light reflected from the inner part the near side of the cup-shaped dust structure. In contrast to the images in visible light the image of HL Tauri on the Fig. 4 received in submillimeter radiation can show the overlapped spherical dust structures around the less than 100000 years young star.

Exoplanetary systems

Generally speaking, the explosion in a protostellar nebula can either accelerate the gas shell rotation around the future star, slow it down or even give it a spin in the opposite direction. The latter case has been observed in some exoplanetary systems (e.g. [START_REF] Narita | First evidence of a retrograde orbit of a transiting exoplanet HAT-P-7b[END_REF][START_REF] Winn | HAT-P-7: a retrograde or polar orbit, and a third body[END_REF][START_REF] Winn | ORBITAL ORIENTATIONS OF EXOPLANETS: HAT-P-4b IS PROGRADE AND HAT-P-14b IS RETROGRADE[END_REF][START_REF] Triaud | Spin-orbit angle measurements for six southern transiting planets-New insights into the dynamical origins of hot Jupiters[END_REF][START_REF] Queloz | WASP-8b: a retrograde transiting planet in a multiple system[END_REF][START_REF] Bayliss | Confirmation of a retrograde orbit for exoplanet WASP-17b[END_REF][START_REF] Hébrard | The retrograde orbit of the HAT-P-6b exoplanet[END_REF] ) and the standard nebular hypothesis cannot explain it, since it implies that the central star and its planets should always rotate and revolve in the same direction, following the rotation of the protostellar nebula. An asymmetry in the initial explosion in the centre of the protostar can also lead to a very strong tilt of the equatorial plane of the planetary system, with tilts of 45 • or even 90 • not impossible (e.g. [START_REF] Pont | The spin-orbit angle of the transiting hot Jupiter CoRoT-1b[END_REF][START_REF] Simpson | The spin-orbit angles of the transiting exoplanets WASP-1b, WASP-24b, WASP-38b and HAT-P-8b from Rossiter-McLaughlin observations[END_REF][START_REF] Winn | The oblique orbit of the super-Neptune HAT-P-11b[END_REF][START_REF] Hirano | A Possible Tilted Orbit of the Super-Neptune HAT-P-11b[END_REF] ).

Background for a physical and mathematical model of the Solar System

Oscillations occur in an environment where gravitational field intensity and gas pressure change substantially within a wavelength. It can be assumed that the wavelength is inversely proportional to the product of average values of gravitational acceleration g, which is inversely proportional to the square of the distance to the centre of attraction, and gas pressure p, which in turn depends directly on g:

λ ∼ 1 g • p (2)
or, supposing that p ∼ g:

λ ∼ r 4 , (3)
where λ is the acoustic wavelength, g and p are the average values of gravitational acceleration and gas pressure as function of r, the distance from the centre. and the scattered disk are represented by Ceres and Eris. The planets originated in the antinodes, so the figure also shows the distribution of the standing wave in the pre-solar gas nebula. We can see that the points on the graph may be approximated by two straight lines, supporting our hypothesis that locations of the standing wave antinodes depend on r 4 . Closer to the Sun, the constant of proportionality determining the wavelength increases, which may be explained by the fact that the speed of the sound wave increases with rising gas temperature: c ∼ √ T ; consequently, with the frequency f = c λ of the sound wave constant,

λ (T ) ∼ √ T , (4)
where c is the speed of the sound wave, λ is its wavelength and T is the temperature. The dependence hypothesis λ ∼ r 4 provides a background for developing a more complex physical and mathematical model for the Solar System's origin, as well as for planetary systems around other stars.

CONCLUSIONS

According to the proposed model of planetary origin from the gas and dust nebula in the field of a standing sound wave, only single stars or multiple stars with large distances between individual components can have planetary systems. Systems with closely located multiple star components cannot sustain the reverse wave after the explosion in the centre of one of the stars causes further periodic explosions due to the motion of the stars, so that it will fade out, while a gas-dust nebula that is not stabilized by the field of a standing sound wave will relatively swiftly accrete to the central stars without any planet formation.

Jensen and Akeson 32 describe the case of a binary young stars with components separated by about 386 AU, each of which is surrounded by a cloud of dust, and the plane of rotation of these clouds are strongly tilted relative to each other. The classical theory can not explain this tilt, the proposed model explains this by saying that the explosions in the centre of each of protostars occur independently of each other, and the dust clouds get their individual tilts.

The classical theory suggests that rotation of the protostellar cloud is essential for the formation of planets around the central star. The proposed model does not require such an assumption. If a gas-dust nebula does not rotate, after the first asymmetric explosion the central protostar and the peripheral gas masses will gain equal angular momentum, but opposite in sign. The gas-dust nebula will start spinning and at the same time a standing sound wave will emerge, beginning the process of planet formation. The resulting planetary system will have the central star and planets rotating in opposite directions.

1 Figure 1 .

 11 Figure 1. Separation of the pre-solar nebula into (1) a central spherical gas nebula and (2) an expanding spherical shell.

Figure 2 .

 2 Figure 2.The standing wave in the pre-solar nebula: (1) nodes, (2) antinodes (gas compression and expansion), (3) protosun, (4) real-scale wave. Nodes and antinodes are evenly spaced for illustration purposes; in reality, separations would increase dramatically with increasing distance from the centre.

Figure 3 .

 3 Figure 3. Migration of solid matter in the spherical nebula. Fg represents the gravitational pull toward the centre and Fc represents the centrifugal force. (1) spherical dust clusters in the antinodes, (2) dust clusters in the equatorial rings, (3) protosun, (4) nebula rotational axis, (5) direction of dust migration in the antinodes, (6) direction of massive dust cluster migration across the antinodes under the influence of the resultant force F.

7 / 13 Figure 4 .

 7134 Figure 4. An image of a protoplanetary dust nebula around HL Tauri received by ALMA (The Atacama Large Millimeter/submillimeter Array). Credit: ALMA (ESO/NAOJ/NRAO).

Figure 5 . 13 Figure 6 .

 5136 Figure 5. The averaged in the line of sight distribution of dust in the spherical antinodes. Closer to centre the dust has already reached the equatorial plane, on the periphery it need to pass a much greater way, so formed a cup-shaped dust distribution structure around the protostar.

Figure 7 .

 7 Figure 7. Distribution of planetary distances from the Sun. The vertical axis shows the fourth root of distance in kilometres.

http://exoplanetarchive.ipac.caltech.edu

2/13

Acknowledgements

This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.


Katerina Gkirtzou

Matthew B Blaschko
email: matthew.blaschko@inria.fr

The Pyramid Quantized Weisfeiler-Lehman Graph Representation

Keywords: Weisfeiler-Lehman algorithm, graph kernels, regularization, vector labels

Introduction

Graphs are a powerful and natural way to represent complex data with integrated structure. Graphs have been used in numerous applications in a number of different fields, such as (i) computer vision [START_REF] Greig | Exact maximum a posteriori estimation for binary images[END_REF][START_REF] Torresani | Feature correspondence via graph matching: Models and global optimization[END_REF] and biomedical imaging analysis [START_REF] Ng | Connectivity-informed sparse classifiers for fMRI brain decoding[END_REF][START_REF] Rao | Adaptive regularization for image segmentation using local image curvature cues[END_REF][START_REF] Goldstein | The genetic graph: A representation for the evolution of procedural knowledge[END_REF], (ii) bioinformatics [START_REF] Canutescu | A graphtheory algorithm for rapid protein side-chain prediction[END_REF][START_REF] Wagner | The small world inside large metabolic networks[END_REF], (iii) social networks analysis [START_REF] Scott | Social Network Analysis: A Handbook[END_REF] and (iv) chemoinformatics [START_REF]Chemoinformatics: A Textbook, 1st Edition[END_REF]. In many applications, the exploration of the data requires the ability to efficiently compare graphs and to provide a similarity measurement, a problem known as graph comparison.

Graph kernels are the most recent approaches that tackle both the problem of graph representation and graph comparison and have been found to be useful across a wide range of applications. They exploit the graph topology by decomposing the graph into sub-structures and aggregating statistics over these substructures. This strategy considers a measure of similarity between the graphs as a form of inner product. Graph kernels are commonly derived as an instance of the family of the R-convolution kernels [START_REF] Haussler | Convolution kernels on discrete structures[END_REF], which are a generic way of constructing kernels of complex objects by decomposing them into discrete structures and comparing all pairs of decompositions. Every new decomposition would yield a new kernel. A first approach would be to decompose the graphs into all possible subgraphs. However, calculating all subgraphs is at least as hard as deciding whether two graphs are isomorphic [START_REF] Gärtner | On graph kernels: Hardness results and efficient alternatives[END_REF]. So it is necessary to limit the decomposition of the graphs only into specific types of subgraphs that are computable in polynomial time [START_REF] Vishwanathan | Graph kernels[END_REF][START_REF] Shervashidze | Weisfeiler-Lehman graph kernels[END_REF]. There are three main categories of graph kernels, graph kernels based on (i) walks and paths, (ii) subgraphs of limited size, and (iii) subtree patterns. Table 1 shows a summary of the state of the art of graph kernels grouped per category, while we review extensively each of these categories in the following Section 1.2.

Although efficient graph kernels have been developed that have good performance on discretely labeled graphs, the literature on continuously or vector labeled graphs is still relatively undeveloped, as is shown in Table 1. As a result, existing graph kernels are poorly suited to applications that are naturally represented by continuous node labeled graphs, such as in graph representations for fMRI analysis and shape classification based on 3D mesh representations. We help to address this shortcoming of existing methods by proposing a novel framework for kernel construction, the pyramid quantized Weisfeiler-Lehman graph representation. Our method makes use of the efficiency of the Weisfeiler-Lehman kernel [START_REF] Shervashidze | Weisfeiler-Lehman graph kernels[END_REF] for discrete labels by converting continuous vector labeled graphs into a sequence of discretely labeled graphs using a pyramid quantization strategy. This results in a kernel with a multiplicative error bound on an optimal matching measure.

Firstly, we exploit the key concepts of the Weisfeiler-Lehman test of isomorphism and then see how it is used to produce subtree patterns for comparing discrete labeled graphs in Section 2. In Section 3 we present our pyramid quantization scheme and we then learn different strategies for combing the pyramid levels in Section 4. Secondly, we perform an experimental comparison with a graph matching technique in Section 5 and we evaluate our proposed method in two general tasks, an fMRI analysis task in Section 6 and 3D mesh shape classification in Section 7. Finally, we conclude with a discussion in Section 8.

Terminology

Before reviewing the graph kernels literature, we clarify our terminology to provide a unified presentation of previous work. We define a graph as a triplet G = (V, E, L), where V is the vertex set, E is the edge set and L : V → Σ is a function assigning a label from an alphabet Σ to each vertex in the graph. The neighborhood N(v) = {v |(v, v ) ∈ E} of a vertex v is the set of all vertices adjacent to v, i.e. all vertices connected with a single edge. The degree d(v) of a vertex v is the number of edges incident with v. Every graph has at most v vertices, e edges and a maximum degree of d. A walk in a graph is a sequence of adjacent vertices. A path is a walk that contains only distinct vertices, while a cycle is a closed walk. A rooted tree is an acyclic graph with a specified root vertex. The height of a tree is the maximum distance between the root vertex and any other vertex in the tree. Subtree patterns are labeled trees extracted from a labeled graph G for a given depth h and a given vertex v. Repetition of the same vertex is allowed in subtree pattern, but it is treated as distinct vertices, allowing a cycle-free pattern.

Related Work

The first class of kernels are based on walks (paths) and compute the number of matching pairs of random walks (paths). The standard formulation of random walk kernel counts the number of nodes in the walk which have the same label and requires O(v 6 ) computation for a pair of graphs [START_REF] Gärtner | On graph kernels: Hardness results and efficient alternatives[END_REF]. Restating the same problem in terms of Kronecker products reduces the runtime complexity to O(v 3 ) for a pair of graphs [START_REF] Vishwanathan | Graph kernels[END_REF]. Another approach based on dynamic programming speeds up the computations of the random walk kernel [START_REF] Harchaoui | Image classification with segmentation graph kernels[END_REF] at the cost of considering only walks of fixed size. The marginalized graph kernel [START_REF] Mahé | Extensions of marginalized graph kernels[END_REF] modifies the label of each vertex with the use of the Morgan index [START_REF] Morgan | The Generation of a Unique Machine Description for Chemical Structures -A Technique Developed at Chemical Abstracts Service[END_REF], which increases the specificity of labels by adding information on the number of walks starting from that node, and therefore reduces the run-time as fewer vertices will match. A graph kernel that compares the length of shortest path between pairs of nodes with matching source and sink labels in two graphs requires O(v 4 ) [START_REF] Borgwardt | Shortest-path kernels on graphs[END_REF]. A specialized graph kernel for chemoinformatics [START_REF] Ralaivola | Graph kernels for chemical informatics[END_REF] uses molecular fingerprinting techniques and counts labeled paths of length p that can be retrieved by depthfirst search from each vertex, an efficient approach for graphs with an average node degree of 2 or 3.

The second class of graph kernels are based on limited-size subgraph structures called graphlets. A naive computation of all graphlets of a graph, without considering labels, requires O(v k ) computations, where k ∈ {3, 4, 5} is the size of subgraphs. Since enumerating all graphlets is prohibitively expensive, even for small values of k, [START_REF] Shervashidze | Efficient graphlet kernels for large graph comparison[END_REF] showed that sampling a fixed number of graphlets suffices to bound the deviation of the empirical estimates of the graphlet distribution from the true distribution and for graphs of degree bounded by d, the exact number of all graphlets of size k can be determined in time O(vd k-1 ). Another kernel is the cyclic pattern kernels [START_REF] Horváth | Cyclic pattern kernels for predictive graph mining[END_REF], which count pairs of matching cyclic and tree patterns in two graphs. In the general case, the cyclic pattern kernel is NP-hard, but in specific cases the kernel can be computed efficiently. Finally, the neighborhood subgraph pairwise distance kernel [START_REF] Costa | Fast neighborhood subgraph pairwise distance kernel[END_REF] decomposes a graph into all pairs of neighborhood subgraphs of small radius r at increasing distances.

The third class of graph kernels is based on subtree patterns. The Ramon-Gärtner subtree kernel [START_REF] Ramon | Expressivity versus efficiency of graph kernels[END_REF] compares all pairs of nodes from two labeled graphs by iteratively comparing their subtree patterns of height h. Although the subtree kernel is more expressive than ker-nels based on walks, unfortunately it is computationally expensive, since for a set of n graphs it requires O(n 2 v 2 h4 d ). Mahé and Vert [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF] and Bach [START_REF] Bach | Graph kernels between point clouds[END_REF] refine the Ramon-Gärtner subtree kernel for applications in chemoinformatics with a parameter to control the complexity of the subtree features and hand-written digit recognition considering α-ary subtrees with at most α children per node. Unfortunately, the complexity of both kernels is still exponential in the smoothing and α parameter, and both kernels are feasible on small sized graphs only. More recently an efficient kernel, the Weisfeiler-Lehman subtree kernel, was introduced in [START_REF] Shervashidze | Weisfeiler-Lehman graph kernels[END_REF]. It uses the Weisfeiler-Lehman test of isomorphism [START_REF] Weisfeiler | A reduction of a graph to a canonical form and an algebra arising during this reduction[END_REF] to efficiently compute subtree patterns up to height h for discretely labeled graphs. For n pairs of discretely labeled graphs, the Weisfeiler-Lehman subtree kernel requires O(nhe + n 2 hv). At its core, the Weisfeiler-Lehman kernel achieves linear complexity in the graph size due to an efficient hashing scheme for discrete labeled graphs, a remarkable improvement over previous approaches.

The majority of these methods focus on either unlabeled or discretely labeled graphs, while an efficient and expressive representation and comparison of graphs with continuous high-dimensional vectors remains an open research problem. Initial results in this direction have been developed by [START_REF] Neumann | Efficient graph kernels by randomization[END_REF], which address the special case of scalar valued node labels with missing data, and by [START_REF] Borgwardt | Scalable kernels for graphs with continuous attributes[END_REF] based on sub-path similarities. We are unaware of previous work that has achieved computational complexity linear in the size of the graph, as in our method.

The Weisfeiler-Lehman kernel

The Weisfeiler-Lehman test of isomorphism

Our proposed algorithm uses the discretely labeled subtree pattern features introduced by the Weisfeiler-Lehman kernel [START_REF] Shervashidze | Weisfeiler-Lehman graph kernels[END_REF], which exploits the key concepts from the Weisfeiler-Lehman test of isomorphism [START_REF] Weisfeiler | A reduction of a graph to a canonical form and an algebra arising during this reduction[END_REF]. The key idea, given two graphs, is the construction of augmented node labels from the labels of all the neighbor nodes and their compression into new short labels. This process is repeated until either the label sets of the two graphs differ or the maximum number of iterations has been reached. In the former case, the two graphs are not isomorphic, while in the latter the test was not able to determine whether they are not isomorphic. Algorithm 1 provides pseudocode for the Weisfeiler-Lehman algorithm, while Figure 1 shows its first iteration. The algorithm has fast run-time O(he), where h is the maximum number of iterations of the test and e the maximum number of edges [START_REF] Shervashidze | Weisfeiler-Lehman graph kernels[END_REF].

The linear Weisfeiler-Lehman subtree kernel

For each iteration i of the Weisfeiler-Lehman algorithm, we obtain new compressed labels l i (v) for all nodes v in all graphs (see Line 12 in Algorithm 1). Each compressed label l i (v) is a subtree pattern of height i rooted at node v and histograms of their occurrences have been recently employed in a kernel for comparing graphs with discrete labels, the linear Weisfeiler-Lehman subtree kernel [START_REF] Shervashidze | Weisfeiler-Lehman graph kernels[END_REF]Definition 4]. Table 2 shows the subtree patterns up to depth h = 1 produced by the Weisfeiler-Lehman algorithm and their respective occurrences between the two graphs G, G shown in Figure 1(a).

The linear Weisfeiler-Lehman subtree kernel is computed as an inner product of a vector of statistics φ (h) (G) explicitly computed for each graph (Table 2). A key advantage of the kernel is its complexity being linear in the graph size O(nhe + n 2 hv) for n graphs, where e is the maximal number of edges, v the maximal number of vertices and h the depth of subtree patterns used [START_REF] Shervashidze | Weisfeiler-Lehman graph kernels[END_REF]Theorem 4]. In addition to these computational benefits, linear Weisfeiler-Lehman graph kernels have been shown to perform comparably to or better than a number of more computationally complex kernels [START_REF] Shervashidze | Weisfeiler-Lehman graph kernels[END_REF]. Finally, we note that the linear Weisfeiler-Lehman algorithm at depth 0 computes exactly the bag of words representation commonly used in natural language processing [START_REF] Harris | Distributional structure[END_REF][START_REF] Ko | A study of term weighting schemes using class information for text classification[END_REF] and computer vision [START_REF] Qiu | Indexing chromatic and achromatic patterns for contentbased colour image retrieval[END_REF][START_REF] Li | A bayesian hierarchical model for learning natural scene categories[END_REF].

The pyramid quantization strategy for continuous vector labels

The Weisfeiler-Lehman algorithm is efficient precisely because it makes use of a discrete labeling over nodes, which enables an efficient hashing scheme. A problem occurs when extending this method to continuous labeled graphs: we no longer have a notion of an exact match of a discrete label, and a hash function that counts approximate matches would implicitly define a single quantization of the vector space to a discrete set of labels. A single quantization is inexact, and gives only a weak relationship to the potentially rich geometry of the original label space. It is also not clear what the resolution of the quantization should be to maximize performance. To overcome this, we propose a pyramid quantization strategy similar to the one used by [START_REF] Grauman | The pyramid match kernel: Efficient learning with sets of features[END_REF][START_REF] Grauman | Approximate correspondences in high dimensions[END_REF] to determine a logarithmic number of quantizations, ending with a sequence of graphs with discrete labels of increasing granularity, for which we run the Weisfeiler-Lehman algorithm.

Table 1: An overview of graph kernel methods. From left to right we show the type of subgraphs used, the algorithm, its complexity when that is known and whether the kernel works on unlabeled, discretely, continuous or vector labeled graphs. Note that n is the number of graphs under comparison, v is the maximal number of nodes, e is the maximal number of edges, h is the height of subtree patterns, d is the maximum degree, k is the size of graphlets, δ is the graph diameter and f is the dimension of the node labels. O(n 2 v 6 ) Mahe et al. [START_REF] Mahé | Extensions of marginalized graph kernels[END_REF] Vishwanathan et al [START_REF] Vishwanathan | Graph kernels[END_REF] O(n 2 v 3 )

Paths

Borgwardt et al. [START_REF] Borgwardt | Shortest-path kernels on graphs[END_REF] O(n 2 v 4 ) Ralaivola et al. [START_REF] Ralaivola | Graph kernels for chemical informatics[END_REF] Feragen et al. [START_REF] Borgwardt | Scalable kernels for graphs with continuous attributes[END_REF] O(n

2 v 2 (m + log v + δ 2 + f ))
Graphlets

Horvath et al. [20] Shervashidze et al. [19]

O(vd k-1 ) Costa et al. [START_REF] Costa | Fast neighborhood subgraph pairwise distance kernel[END_REF] Subtree Patterns Ramon et al. [START_REF] Ramon | Expressivity versus efficiency of graph kernels[END_REF] O(n 2 v 2 h4 d ) Bach et al. [START_REF] Bach | Graph kernels between point clouds[END_REF] Mahe et al [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF] Shervashidze et al. [START_REF] Shervashidze | Weisfeiler-Lehman graph kernels[END_REF] O(nhe + n 2 hv) Neumann et al. [START_REF] Neumann | Efficient graph kernels by randomization[END_REF] Table 2: An example of the subtree patterns up to depth h = 1 for the graphs from Figure 1(a). The first row are the labels Σ 0 ∪ Σ 1 encountered up to depth h = 1, while the second and the third row the occurrences of subtrees for graph G and G respectively

Labels Σ 0 Labels Σ 1 Labels Σ 0 ∪ Σ 1 = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11} φ (1) (G) = (3, 1, 3, 0, 3, 0, 1, 2, 1, 0, 0) φ (1) (G ) = (2, 2, 3, 1, 1, 1, 1, 1, 0, 1, 1)
3.1. The pyramid quantization strategy Given a vector labeled graph G = (V, E, L), where L : V → R d is the function assigning a d-dimensional vector label to each vertex, we want to derive a hierarchical decomposition of R d as multi-resolution quantizations that will then be used to determine the discrete labeling of increasing granularity. First we construct a set of quantization functions Q (l) : R d → Σ (l) 0 that will encode the continuous vector labels into a quantization of a given resolution |Σ (l) 0 | = 2 l . The quantization function Q (l) is generated for multiple resolutions l ∈ {0, . . . , L}, where L = log 2 D , D is the number of unique values in the image of the vertex set V under the label function L. Note that the single quantization bin for Q (0) is big enough so that all vertices receive the same discrete label, while the coarser quantization Q (L) contains quantization bins that are small enough so each unique continuous vector label falls into its own quantization bin. To achieve this hierarchical quantization in the experiments performed here, we have used an agglomerative hierarchical clustering with Ward's minimum variance method [START_REF] Ward | Hierarchical Grouping to Optimize an Objective Function[END_REF].

Secondly, we compose the quantization function Q (l) with the labeling function L, ∀l ∈ {0, . . . , L}, and we approximate our initial vector labeled graph G as a sequence of graphs with discrete labels of increasing granularity:

G = (V, E, L) Q (l) • L → G (0) , . . . , G (L) (1) = (V, E, L (0) ), . . . , (V, E, L (L) ) ,
where L (l) : V → Σ (l) 0 is defined to be Q (l) • L, and Σ (l) 0 is the discrete label alphabet for a given level l of quantization. Note that the topology of the graph does not change in the sequence of graphs, only the continuous vector labels are discretized. This type of quantization scheme, when paired with a histogram intersection kernel and an appropriately weighted linear combination of kernel values across Algorithm 1 The one dimensional Weisfeiler-Lehman test of graph isomorphism Require: Two graphs, G = (V, E, L), G = (V , E , L ), with discrete labelings L : V → Σ and L : V → Σ over vertices, where Σ is a vertex label set and the maximum number of iterations h.

1: i ← 0 2: repeat 3: if i = 0 then { Multiset-label initialization } 4: M i (v) := l 0 (v) = L(v).
5:

else if i ≥ 1 then { Multiset-label determination } 6: Assign a multiset-label M i (v) to each node v in G and G which consists of the multiset {l i-1 (u)|u ∈ N(v)}, where N(v) denotes the neighbor set of v. { Sorting each multiset } 7:
Sort the elements in M i (v) in ascending order.

8:

Concatenate the elements in M i (v) into a string s i (v).

9:

Add l i-1 (v) as a prefix to s i (v). { Sorting the set of multisets } 10:
Sort all of the strings s i (v) for all v from G and G in ascending order. { Label compression via hashing } 11:

Map each string s i (v) to a new compressed label using a function f :

Σ * → Σ such that f (s i (v)) = f (s i (w)) ⇐⇒ s i (v) = s i (w). { Relabeling } 12: Set l i (v) := f (s i (v)
) for all nodes in G and G .

13:

end if

14: i ← i + 1 15: until {l i (v)|v ∈ V} {l i (v )|v ∈ V } or i > h
quantization levels, results in a multiplicative error bound on the optimal graph matching [32, Proposition 3]. We may therefore interpret the pyramid quantized Weisfeiler-Lehman graph representation as a function space that enables tight approximations to cost of the optimal matching over vector representations of subtree patterns. 1

The intersection Weisfeiler-Lehman subtree kernel

Each graph with vector valued labels after the pyramid quantization step described in Section 3.1 is represented as a sequence of graphs with nested quantizations of increasing granularity as described in Equation 1. We construct a kernel, using the following definition:

Definition 3.1 (The intersection Weisfeiler-Lehman subtree kernel). Let G (l) = (V, E, L (l) ) and G (l) = (V , E , L (l)
) be two graphs of the same quantization level l, where L (l) : V → Σ (l) 0 and L (l) : V → Σ (l) 0 , of two vector labeled graphs G = (V, E, L) and G = (V , E , L ), where L : V → R d and L : V → R d . Define Σ (l) i ⊆ Σ (l) as the set of symbols that occur as node labels at least once in G (l) or G (l) at the end of the i-th iteration of the Weisfeiler-Lehman algorithm. Let Σ (l) 0 be the set of original node labels of G (l) and G (l) . Assume all Σ (l) i are pairwise disjoint. Without loss of generality, assume that every

Σ (l) i = {σ (l) i1 , . . . , σ (l) i|Σ (l) i | } is ordered. Define a map η i : {G (l) , G (l) } × Σ (l) i → N such that η i (G (l) , σ (l) i j
) is the number of occurrences of the letter σ (l) i j in the graph G (l) . The intersection Weisfeiler-Lehman subtree kernel on two graphs G and G with h iterations is defined as

K (h) i-WLsubtree (G (l) , G (l) ) := I φ (l) (h) (G (l) ), φ (l) (h) (G (l) ) (2)
where Figure 1: An illustration of the Weisfeiler-Lehman algorithm for the iteration i = 1. Note that there is a one-to-one correspondence between the compressed labels and subtree patterns. For example, the compressed label 9 corresponds to a subtree pattern of height 1 where the root node has label 3 and its children nodes have labels 1, 2 and 3 respectively. and

φ (l) (h) (G (l) ) = η 0 (G (l) , σ (l) 01 ), . . . , η h (G (l) , σ (l) h|Σ (l) h | ) (3)
φ (l) (h) (G (l) ) = η 0 (G (l) , σ (l) 01 ), . . . , η h (G (l) , σ (l) h|Σ (l) h | ) (4)
I φ (l) (h) (G (l) ), φ (l) (h) (G (l) ) = (5) h i=0 |Σ (l) i | j=1 min η i (G (l) , σ (l) i j ), η i (G (l) , σ (l) i j )
Note that the intersection Weisfeiler-Lehman subtree kernel counts the overlap of subtree features φ (l) (h) between two graphs at the same quantization level l and for a given binning resolution l the intersection kernel is a positive-definite similarity function [START_REF] Odone | Building kernels from binary strings for image matching[END_REF]. We define and propose the use of the intersection kernel instead of the dot-product kernel, as it can give us better control over feature matches per level due to the overlapping label representation. Moreover, the histogram intersection kernel has been shown to give good results in a number of applications [START_REF] Odone | Building kernels from binary strings for image matching[END_REF].

For the pyramid quantized Weisfeiler-Lehman kernel described above the following monotonicity property holds: Theorem 3.2 (Monotonicity property of the pyramid quantized Weisfeiler-Lehman kernel with the granularity of the node labeling). The Weisfeiler-Lehman algorithm for a given height h of subtree patterns produces histograms whose intersection are monotonically decreasing in the granularity of the graph node labeling:

I φ l (h) (G (l) ), φ l (h) (G (l) ) ≥I φ l+1 (h) (G (l+1) ), φ l+1 (h) (G (l+1) ) ∀l, G (l) , G (l) , (6)
where φ l (h) (G (l) ) is the histogram of subtree patterns of height h computed at pyramid level l, and level l + 1 is more granular.

Proof. We first note that the number of subtree patterns of a given depth, h, of the Weisfeiler-Lehman algorithm is dependent only on the topology of the graph, and not on the graph labeling:

φ l (h) (G (l) ) 1 = v • (h + 1) ∀l ( 7
)
where v is the number of vertices in the graphs. We next note that the number of vertex labels is strictly monotonic in the pyramid level, |Σ l | < |Σ l+1 |, and that for each label σ ∈ Σ l at level l, there exist a nonempty set of labels S ⊂ Σ l+1 at level l + 1 such that L l+1 (u) ∈ S ⇐⇒ L l (u) = σ. To complete the proof, we observe that

∀G (l) ∀l φ l (h) (G (l) ) 1 = φ l+1 (h) (G (l+1) ) 1 ∧ φ l (h) (G (l) ) 0 < φ l+1 (h) (G (l+1) ) 0 =⇒ ∀G (l) , G (l) ∀l I φ l (h) (G (l) ), φ l (h) (G (l) ) ≥ I φ l+1 (h) (G (l+1) ), φ l+1 (h) (G (l+1) ) , (8)
where • 0 is the 0 pseudo norm.

Proposition 3.3. The intersection Weisfeiler-Lehman subtree kernel (Definition 3.1) specializes to a bag-ofwords intersection kernel for h = 0.

Proof. For h = 0, η l maps to Σ l 0 , a single quantization. The outer summation in Equation ( 5) is over one term, and the inner summation computes a histogram intersection over a single fixed quantization.

As the Weisfeiler-Lehman algorithm with h = 0 specializes to the bag of words model as a result of Proposition 3.3, we have as a result that our graph kernel for continuous vector valued node labels strictly generalizes the pyramid match kernel [START_REF] Grauman | The pyramid match kernel: Efficient learning with sets of features[END_REF]. Theorem 3.2 generalizes the monotonicity exploited in the construction of the pyramid match kernel to guarantee that Mercer's condition holds [32, Equation ( 4) & Proposition 1].

The pyramid quantized Weisfeiler-Lehman kernel complexity

The main steps of the pyramid quantized Weisfeiler-Lehman kernel are the calculation of the multiresolution hierarchical quantization of the continuous vector labels of the vertices, the assignment of the initial continuous vector label of each vertex to the appropriate path of the hierarchical quantization and the calculations of the subtree pattern features. Calculating a multi-resolution hierarchical quantization of m unique values requires O(m 3 ) calculations. In the ideal case, we would like to use all unique values in the image of the vertex set V under the label function L : V → R d , but since in some applications this number can be prohibitively large, a representative subset of size m D -where D is the number of unique values in the image of the vertex set V under the label function L -can be sufficient. Independent of whether all unique label values or a subset is used to calculate the multi-resolution hierarchical tree, this process should be distinguished as a separate step that can be done once per application domain, and it is not required during test time.

During test time, the steps of the pyramid quantized Weisfeiler-Lehman graph representation are the assignment of the continuous vector label of each vertex to the appropriate path of the multi-resolution hierarchical quantization and the calculation of the subtree pattern features. Note that since we want to approximate our initial vector labeled graph G(V, E, L) as a sequence of L = log 2 D graphs with discrete labels of increasing granularity, that means we need to calculate L Weisfeiler-Lehman kernels. For the assignment of a given vector label to the hierarchical quantization, it should be emphasized that the multi-resolution hierarchical quantization can be seen as a complete binary tree of height L. So the label's placement in the hierarchical quantization tree is determined by comparing it to the appropriate two bin centers at each of the L levels. So we start at the second coarser quantization level Q (1) where all labels receive the only two discrete labels, (the coarser level Q (0) is a trivial one where all labels get the same discrete label) and we compare it with the two bin centers. The label gets the assigned discrete label of the bin center that is nearest to it in terms of the distance used to create the hierarchical quantization tree. Then, we push the label down to hierarchical quantization tree and continue to finer levels only along the branch that is chosen at each level. So a label is first assigned to one of the coarser levels and then it is assigned to one of its children, and so on recursively. This amounts to a total of 2L comparisons in order to find all multiresolution discrete labels for a given vector label and O(vL) for a given graph. The complexity of label assignment is dominated by the complexity of the kernel computation for L quantization levels, and consequently the overall complexity for n graphs is O(L(nhe + n 2 hv)).

Learning to combine pyramid quantized Weisfeiler-Lehman features

Given a set of vector labeled graphs G = {G i = (V i , E i , L i )} 1≤i≤n where L : V → R d and a classification label y i ∈ Y for each graph G i , we classify graphs using the representation created by the Weisfeiler-Lehman algorithm after their quantization as described in Section 3. In order to maximize their classification performance we examine two different approaches, one through multiple kernel learning (Section 4.1) and another through direct sparse regularization of individual subtree patterns (Section 4.2).

Multiple kernel learning

Applying the intersection Weisfeiler-Lehman subtree kernel for each pair of graphs G (l) , G (l) for all the pyramid levels from Equation 1, we end up with a sequence of intersection Weisfeiler-Lehman kernels for a given height h of subtree patterns : ) ) . Taking also into consideration the observation by [START_REF] Lanckriet | A statistical framework for genomic data fusion[END_REF] that using multiple kernels instead of a single one can improve performance, we would like to combine our sequence of kernels. A common approach considers that the kernel K(G, G ) is actually a convex combination of basis kernels:

K (0) (h) (G (0) , G (0) ), . . . , K (L) (h) (G (l) , G (L
K (h) (G, G ) = L l=0 d l K (l) (h) (G (l) , G (l) ), d l ≥ 0, L l=0 d l = 1.
(9) For determining the weights d l we consider two different approaches: an automatic way through the framework of multiple kernel learning and another through a fixed weight scheme.

The automatic determination of the weights d l of the linear combination of our multiple kernels K (l) (h) (G (l) , G (l) ) as well as the coefficients α i , β in a single optimization problem is known as the multiple kernel learning (MKL) problem [START_REF] Zien | Multiclass multiple kernel learning[END_REF][START_REF] Sonnenburg | Large scale multiple kernel learning[END_REF][START_REF] Lanckriet | A statistical framework for genomic data fusion[END_REF][START_REF] Rakotomamonjy | Sim-pleMKL[END_REF]. The multiple kernel learning approach addresses the problem through a weighted 2 norm regularization formula. In addition, a 1 norm is posed as a constraint on the kernel weights a i . This additional constraint encourages sparse set of basis kernels as an inherited property from the 1 norm. The primal MKL problem is defined as

min f,β,ξ,d 1 2 L l=0 1 d l || f l || 2 H l + C n i=1 ξ i s.t. y i L l=0 f l (G (l) i ) + y i β ≥ 1 -ξ i ∀i (10)
ξ i ≥ 0 ∀i, L l=0 d l = 1, d l ≥ 0 ∀l
where each function f l belongs to a different RKHS H l associated with a kernel K (l) (k) . The constraint of the 1 norm on the weighted leads to sparse solutions, which may be undesirable if all pyramid levels are informative. We therefore additionally consider taking their empirical mean:

K (h) (G, G ) = L l=0 d l K (l) (h) (G (l) , G (l) ), where d l = 1 L + 1 . (11)
Visualization. When we consider a linear combination of our kernels and since the intersection kernel can be considered as a "quasi-linear" kernel [START_REF] Vedaldi | Multiple kernels for object detection[END_REF], we can develop visualizations that approximate the learned discriminant functions. The discriminant function for class c has the form:

f c (G) =b + L l=0 d l j α l j I φ (l) (h) (G (l) j ), φ (l) (h) (G (l) ) (12) ≈b + L l=0 d l w l , φ (l) (h) (G (l) ) ( 13
)
where l indexes the levels of the pyramid and j indexes over the samples in the training set. At each pyramid level l there is exactly one subtree pattern of a given height h rooted at each vertex of the graph. We may generate for each vertex v a visualization by coloring each vertex by the weight in w i corresponding to the subgraph pattern rooted at that node.

Elastic Net on the pyramid quantized Weisfeiler-Lehman subtree features

We also consider a direct sparse regularization on the subtree patterns using the Elastic Net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]. The Elastic Net linearly combines 1 with 2 regularization in order to appropriately trade off sparsity with a low variance estimator in the case of correlated signals. Formally, if φ (l) (h) (G (l) i ) is a feature vector of subtree pattern up to height h for a given quantization level l for a graph G i , the Elastic Net computes β = arg min

β∈R d λ 2 β 2 2 +λ 1 β 1 + 1 n n i=1 β, φ (l) (h) (G (l) i ) -y i 2 , (14)
where λ 1 , λ 2 ≥ 0 are scalar regularization parameters.

Experimental Comparison with Graph Matching

Graph matching has been proposed to construct graph kernels. The matching energy is used as a kernel value [START_REF] Duchenne | A graph-matching kernel for object categorization[END_REF][START_REF] Liu | Graph matching by simplified convex-concave relaxation procedure[END_REF][START_REF] Yang | Outlier robust point correspondence based on GNCCP[END_REF]. Although this does not lead to a strictly positive definite function, it can nevertheless yield reasonable performance in some domains. In this section, we construct synthetic data sets that show the relative advantages of the pyramid quantized Weisfeiler-Lehman kernel vs. graph matching based approaches.

We perform an experimental comparison using the factorized graph matching technique [START_REF] Zhou | Factorized graph matching[END_REF]. We perform a classification experiment with synthetic node labeled graphs from two classes.

In order to create the graphs, we consider a 2D unit square plane, two labeling schemes and two Gaussian distribution for labeling the nodes. Each labeling scheme represents a different graph class. In the first labeling scheme, i.e. for the graphs from class 0, we split the 2D plane in half on the x-axis and for the nodes that are placed on the left side are labeled from Distribution A, otherwise they are labeled from Distribution B (see Figure 2(a)). In the second labeling scheme, i.e. for the graphs from class 1, the 2D plane is split by a square in the middle of the plane occupying the half area of the 2D plane and for the nodes that are placed within the square are labeled from Distribution A, otherwise they are labeled from Distribution B (Figure 2(b)).

We use two different edge construction methods. In the first variant, we ensure that nodes lie on a regular grid, and use a four-connected graph topology. In the second variant, we randomly sample the locations of the nodes, and construct the topology using a knearest neighbor graph construction. We performed 10fold cross validation experiments with a fixed c value (c = 0.01) and accuracies are shown in Table 3.

Experiments on an fMRI analysis problem

We evaluate the pyramid quantized Weisfeiler-Lehman graph representation using real data from two different domains. The first evaluation described comes from the fMRI analysis area, while the second described in Section 7 focus on a 3D shape mesh classification problem. We approach the fMRI analysis by representing fMRI recordings as graphs, and we use our proposed method to learn from the interconnections between voxels. Our approach has an enriched capacity to model such dependencies by considering interconnections between voxels which may be functionally important. 2

Cocaine addiction dataset

We evaluate the approach on a dataset [START_REF] Goldstein | Anterior cingulate cortex hypoactivations to an emotionally salient task in cocaine addiction[END_REF][START_REF] Honorio | Can a single brain region predict a disorder?[END_REF][START_REF] Gkirtzou | fMRI analysis of cocaine addiction using k-support sparsity[END_REF][START_REF] Gkirtzou | fmri analysis with sparse weisfeiler-lehman graph statistics[END_REF][START_REF] Belilovsky | Predictive sparse modeling of fMRI data for improved classification, regression, and visualization using the k-support norm[END_REF] that contains an approximately equal number of co-caine addicted individuals and control subjects performing a neuropsychological experiment of block design, called a drug Stroop experiment. The classification task is to discriminate cocaine from control subjects. The data were preprocessed using SPM2 [START_REF] Friston | Statistical parametric maps in functional imaging: A general linear approach[END_REF] and a contrast map for each subject was produced. Only the subjects that complied to motion < 2mm translation, < 2 • rotation and at least 50% performance of the subject in an unrelated task [START_REF] Goldstein | Anterior cingulate cortex hypoactivations to an emotionally salient task in cocaine addiction[END_REF] were kept.

Graph construction

We build a graph representation from the contrast maps by determining a subset of voxels with the use of the Elastic Net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] on the raw voxel values. This method is particularly appropriate in fMRI where nearby voxels are likely to be correlated, and regions responsible for a given function or behavior distributes across multiple voxels. We have made use of k-nearest neighbor graphs on the voxels that were selected by an initial training of the Elastic Net. We symmetrize the k-nn relationship by considering the edges to indicate an undirected graph structure. While other models of connectivity are of interest [START_REF] Sporns | Networks of the Brain[END_REF][START_REF] Wee | Enriched white matter connectivity networks for accurate identification of MCI patients[END_REF], we have found that the use of k-nearest neighbors to determine the graph topology yields good performance in general and illustrates the advantages of the pyramid Weisfeiler-Lehman approach. Furthermore, the subtree statistics considered here implicitly account for longer distance connections for sufficiently deep subtree patterns. We set the number of neighbors k = 5 in all experiments.

To enrich our graph representations of the fMRI contrast maps, we take advantage of the activation information. At each voxel selected by the Elastic Net for the construction of the graph, we label it with its activation. Since the activation has continuous values, our graph representation is transformed to a continuous labeled graph. We subsequently compute subgraph statistics over this graph to generate a feature vector, φ (l) (h) (G (l) ) for a given height h of subtree patterns and a given quantization level l for a graph G. Finally, we use the Elastic Net on these subgraph statistics over all quantization levels, in order to determine our final prediction function, with a model selection step to determine appropriate values for λ 1 and λ 2 .

Results

We use the same experimental setup, a random splitting scheme with 50 trials, to estimate the classification performance of pyramid quantized Weisfeiler-Lehman graph representation and the baseline methods. In each trial, a random selection of 80% of the data are used for training, while the remaining 20% are used to estimate the performance. In Table 4 we show the performance of the pyramid quantized Weisfeiler-Lehman graph representation for four different depths of subtree patterns.

Our approach achieves a mean accuracy of 64.28% for subtree patterns up to depth two, a significant improvement over the bag of words kernel (h = 0). We also compare our proposed technique with three other methods on the same dataset: (i) Gaussian kernel ridge re-gression, (ii) the Elastic Net with raw voxels as features, and (iii) the Elastic Net with raw voxels and pyramid quantized Weisfeiler-Lehman subtree features concatenated in a joint feature vector. In Figure 4 we show the mean accuracy of the final system and the standard error. The pyramid quantized Weisfeiler-Lehman graph representation outperforms the rest of the methods and with a Wilcoxon signed rank test we determine that our proposed method is statistically significantly better (p = 0.02). Additionally, a reduction of over 14% in classification error is recorded between the Elastic Net on the raw voxels and our method. Figure 3(a) shows the areas selected by the Elastic Net, while Figure 3(b) and Figure 3(c) show the visualizations of the learned functions for control and cocaine addicted subjects, respectively. Note that Elastic Net on the raw voxels selected the rostral anterior cingulate cortex, an important region for addictive behavior [START_REF] Gkirtzou | fMRI analysis of cocaine addiction using k-support sparsity[END_REF][START_REF] Gkirtzou | fmri analysis with sparse weisfeiler-lehman graph statistics[END_REF]. Although our method works in an implicitly high dimensional space, we empirically observe that Elastic Net controls the complexity at each stage of the pipeline. The first learning step selects approximately 1100 voxels. Using our method, we generate a feature vector of length 6 × 10 5 , but with a sparsity of ∼ 2%. The second application of Elastic Net selects only 2K dimensions. In each step, the method retains complexity much lower than a "simple" linear function over tens of thousands of voxels as in previous works.

Several broad observations are apparent from our quantitative results. From Table 4, we note that subtree patterns up to depth two performs best, and that deeper subtree patterns begin to reduce average performance. This indicates that the big-O complexity of the graph representation is only slightly higher than using a simple linear function as we use very small values of h in practice (see Section 3.3). The proposed method performs significantly better than the baselines (see Table 4 and Figure 4). In our final experiment of combining the raw voxel values with the subtree pattern features, we found that performance decreased slightly from that of only considering subtree pattern features.

Experiments on 3D shape classification

In this section, we evaluate the pyramid quantized Weisfeiler-Lehman kernel on two 3D shape categorization tasks. In our approach, we denote that the 3D shape can be viewed as graphs G(V, E), where the finite set of points in the 3D space will represent the vertices V and the connection between two points in order to form triangles or curved surfaces will represent the edges E. This perspective specifies the topology of a graph, but does not explicitly encode relative vertex positions or other geometric properties. Therefore, we extend the notion of the graph to incorporate node labels that encode properties such as local curvature of the surface.

3D shapes datasets

The evaluation is performed on two 3D shape datasets. The first dataset is the neuromuscular dystrophy dataset, which consists of 41 subjects: 27 are affected by a neuromuscular dystrophy, while the remaining 14 subjects are healthy [START_REF] Neji | Diffusion tensor imaging of the human skeletal muscle: Contributions and applications[END_REF][START_REF] Gkirtzou | Sparse classification with MRI based markers for neuromuscular disease categorization[END_REF]. The subjects were imaged in the calf using a 1.5 T MRI scanner. An example of the T1-weighted MR images of the calf from a healthy and patient subject can be seen in Figure 5.

The T1 weighted MR images were manually segmented by an expert separating 7 important calf muscle groups and each segmented muscle is then transformed into a 3D surface mesh using the itk-snp program. 3 The SHREC 2013 dataset consists of 20 classes of generic objects. 4 Each class contains 18 different large-scale models, resulting in a total of 360 objects.

Node label description

We select as labels on the vertices of the 3D surface mesh a weighted average over the principal curvature features of the immediately adjacent triangulated faces [START_REF] Rusinkiewicz | Estimating curvatures and their derivatives on triangle meshes[END_REF]. In the SHREC 2013 dataset, we also used as labels multi-viewpoint rendering features, a successful method for 3D shape classification [START_REF] Ohbuchi | Distance metric learning and feature combination for shape-based 3d model retrieval[END_REF], but only for the 3D mesh enclosed within a fixed radius. As we cannot assume a canonical basis for specifying the 3D coordinates of the surface control points, we use PCA to determine one. For comparison, we also develop a multiviewpoint rendering baseline. Similarly, we use PCA to determine a basis. We then render images in these canonical bases and linear, polynomials of 2nd and 3rd degree and Gaussian kernels are computed.

Results

For both datasets we use the same experimental setup, a nested cross-validation procedure where the inner loop is used for automatic model selection. We only report results from the outer cross-validation procedure ensuring the statistical validity of our model selection. The performance is evaluated as the accuracy and the area under the ROC curve (AUC). We also compare our method with two other baselines on both datasets (a) a pyramid bag of words model and (b) a multi-viewpoint rendering procedure following the same experimental setup. We also present the results obtained from the combination of the multi-viewpoint rendering representation with the pyramid quantized Weisfeiler-Lehman kernel.

The performance on the neuromuscular dystrophy dataset is shown in Table 5. The pyramid quantized Weisfeiler-Lehman kernel outperforms both baseline methods.

The overall best performance is achieved when we combined our method with the multiviewpoint rendering approach with a accuracy of ∼ 83% and an AUC of 0.6648. The performance for the SHREC 2013 dataset is shown in Table 6. The overall best performance is achieved when we combined the pyramid quantized Weisfeiler-Lehman kernel with the multi-viewpoint rendering images with an AUC of ∼ 0.85 across all 20 classes. A Wilcoxon signed-rank test showed that the combined method performed better than all other methods with high statistical significance (p < 10 -3 ). We further show the learned weight of In the neuromuscular dystrophy dataset, the pyramid quantized Weisfeiler-Lehman kernel performs substantially better than both baselines. These techniques clearly contain complementary information as the combined method performs best. We have further confirmation from the SHREC 2013 dataset since the combined approach gives the best average performance with statistical significance.

Discussion

In this article, we have presented a fully automated, efficient and statistically sound framework for graph classification with continuous vector node labels using the pyramid quantized Weisfeiler-Lehman graph representation. It combines the Weisfeiler-Lehman graph kernel [START_REF] Shervashidze | Weisfeiler-Lehman graph kernels[END_REF] with the pyramid match kernel [START_REF] Grauman | The pyramid match kernel: Efficient learning with sets of features[END_REF], in order to learn from graphs with continuous labels through a pyramid quantization scheme. The method was evaluated on a wide variety of datasets and outperformed other machine learning techniques with statistical significance. We have demonstrated the utility of such graph representations in two main areas: fMRI analysis and 3D shape classification.

The fMRI application was evaluated on a real world cocaine addiction dataset and outperformed other methods with statistical significance, including kernel ridge regression and the Elastic Net. This validates our hypothesis for fMRI analysis: that the interconnections between voxels can contain additional information about brain structure that is not apparent in a linear function on the raw voxel values. The 3D shape classification The mean AUC on the SHREC 2013 dataset. The compared methods are (left to right) the pyramid quantized Weisfeiler-Lehman kernel (WLpyramid), the pyramid bag of words model (pyramid BoW), the multi-viewpoint rendering procedure (Rendering), and a combination of the multi-viewpoint rendering procedure with the pyramid quantizes Weisfeiler-Lehman kernel (Combined). In bold are the one-vs-rest classifier where the combined classifier outperforms the multi-viewpoint rendering procedure.

tasks showed consistent improvement over baselines using the proposed method. We have not directly incorporated any features in our node labels capturing surface reflectance, color, or texture. This is an interesting area for future research. Learned shape retrieval by discriminative training of a Mahalanobis metric [START_REF] Weinberger | Distance metric learning for large margin nearest neighbor classification[END_REF] is another interesting possible future direction that is directly applicable with the proposed representation.

We have made five main contributions in this work, (i) we applied a generalization of the Weisfeiler-Lehman graph kernel to continuous vector node labels, the pyramid quantized Weisfeiler-Lehman kernel, (ii) we showed that graph representations provide a rich family of functions for fMRI analysis, (iii) we developed a novel framework for shape classification based on the interpretation of shape meshes as annotated graphs, (iv) we performed experiments on fMRI analysis, medical imaging, and semantic shape classification tasks with improved classification accuracy across diverse applications and (v) we showed visualizations of the learned discriminant function, providing rich information about the discriminative power of the method.

A

  lg o r it h m C o m p le x it y U n la b e le d D is c r e t e C o n t in u o u s V e c t o r WalksGärtner et al.[START_REF] Gärtner | On graph kernels: Hardness results and efficient alternatives[END_REF]

  Given labeled graphs G and G .

  (a) Labeling schema for class 0. (b) Labeling schema for class 1.

Figure 2 :

 2 Figure 2: The spatial distribution of node labels for the synthetic experiments in Section 5.

Figure 3 :

 3 Figure 3: A visualization of the areas of the brain selected by Elastic Net as well as a visualization of the learned functions on the quantized Weisfeiler-Lehman representation.

Figure 4 :

 4 Figure 4: Mean accuracy and standard error on the cocaine addiction dataset. The compared methods are (left to right) Gaussian kernel ridge regression (GKRR), the Elastic Net on raw voxels, pyramid quantized Weisfeiler-Lehman (WLpyramid), and the Elastic Net with a concatenation of the raw voxels and the pyramid quantized Weisfeiler-Lehman features (Combined EN+WL). The horizontal red line indicates chance performance.

Figure 5 :

 5 Figure 5: On the left a slice of the MR image from a healthy subject, and on the right a slice of the MR image from a patient with a neuromuscular disease.

4 2 Figure 6 :

 26 Figure 6: An example of the learned weights of the pyramid quantized Weisfeiler-Lehman Kernel of the class bird from the SHREC 2013 dataset for three different subtree depths h ∈ {0, 1, 2}. (Figure best viewed in color).

Figure 7 :

 7 Figure 7: Visualization of the learned weights of the pyramid quantized Weisfeiler-Lehman kernel of subtree patterns with depth h = 1 of all classes for the SHREC2013 dataset. (Figure best viewed in color.)

Table 3 :

 3 A comparison of the Weisfeiler-Lehman pyramid kernel with a graph matching based approach. The Weisfeiler-Lehman kernel is adaptable to different graph topologies, while graph matching can fail when the topology has significant local random variations as in the k-NN graph construction. ± 1.[START_REF] Odone | Building kernels from binary strings for image matching[END_REF] 48.08 ± 2.12 48.08 ± 2.24 53.33 ± 2.51 64.00 ± 2.55 63.67 ± 2.15 64.33 ± 2.26 k-NN ± 1.57 44.80 ± 0.68 82.22 ± 0.51 84.15 ± 0.31 85.77 ± 0.57 84.88 ± 0.58 84.60 ± 0.6

	Experiment	Graph Matching	Bag of Words	k=1	k=2	WLpyramid k=3	k=4	k=5
	Grid	72.42

Table 4 :

 4 Mean accuracy on the cocaine addiction dataset of the pyramid quantized Weisfeiler-Lehman graph representation for four different subtree pattern depths, h ∈ {0, 1, 2, 3}. Maximum performance is achieved with subtree patterns up to depth two.

		Pyramid Quantized Weisfeiler-Lehman
	h	0	1	2	3
	Accuracy 54.00% 57.14% 64.28% 63.42%
			Accuracy for different methods
		0.75
		0.7
	Mean Accuracy	0.55 0.6 0.65
		0.5
		GKRR	Elastic Net	WLpyramid Combined WL+voxels

Table 5 :

 5 The mean accuracy and the mean AUC on the neuromuscular dystrophy dataset. The compared methods are (left to right) the pyramid quantized Weisfeiler-Lehman kernel (WLpyramid), the pyramid bag of words model (pyramid BoW), the multi-viewpoint rendering images procedure (Rendering) and a combination of the multi-viewpoint rendering procedure with the pyramid quantized Weisfeiler-Lehman kernel (Combined).

		WLpyramid pyramid BoW Rendering Combined
	Accuracy	78.00%	73.00%	75.50%	82.93%
	AUC	0.6410	0.6361	0.6300	0.6648
	Class	WLpyramid pyramid BoW Rendering Combined
	Bird	0.85	0.83	0.85	0.86
	Bicycle	0.84	0.87	0.90	0.90
	Biped	0.89	0.88	0.99	0.99
	Biplane	0.60	0.63	0.68	0.69
	Bird	0.73	0.73	0.80	0.80
	Bottle	0.76	0.76	0.79	0.80
	Car	0.78	0.79	0.80	0.80
	CellPhone	0.74	0.80	0.88	0.89
	Chair	0.69	0.68	0.70	0.72
	Cup	0.85	0.84	0.88	0.88
	Desklamp	0.80	0.80	0.88	0.89
	Fish	1.00	1.00	1.00	1.00
	Floorlamp	0.80	0.77	0.89	0.89
	Insect	0.64	0.60	0.62	0.66
	Monoplane	0.84	0.82	0.88	0.90
	Mug	0.82	0.82	0.85	0.87
	Phone	0.83	0.74	0.72	0.83
	Quadruped	0.89	0.86	0.97	0.98
	Sofa	0.76	0.75	0.74	0.75
	Wheelchair	0.81	0.79	0.88	0.90
	Average	0.80	0.79	0.84	0.85

Table 6 :

 6

In practice, we do not use the fixed weighting across quantization levels proposed by Grauman and Darrell[START_REF] Grauman | The pyramid match kernel: Efficient learning with sets of features[END_REF], but instead discriminatively optimize over the function space (see Section 4).

This section is based in part on results originally published in[START_REF] Gkirtzou | fmri analysis with sparse weisfeiler-lehman graph statistics[END_REF].

http://www.itksnap.org/pmwiki/pmwiki.php

Acknowledgments

This work was partially funded by ERC Grant 259112, FP7-MC-CIG 334380, NIA 1R21DA034954-01, and the AFM-Telethon foundation. We thank Jean-Franc ¸ois Deux, Guillaume Bassez, Alain Rahmouni, and Thibault Varacca at Henri Mondor University Hospital, Créteil, France for collecting the muscle dataset. We thank Rita Goldstein at Icahn School of Medicine, Mount Sinai, and Dimitris Samaras and Jean Honorio at 14 Stony Brook University for collecting and preprocessing the fMRI dataset.