Edouard Ribes
email: edouard.augustin.ribes@gmail.com

Growth patterns of US professional services firms

Keywords: Firm growth Professional Services Firms Firm performance Mean Field Games Economic Classification [JEL]. E17, L16, L25, L84, D23

Growth patterns of US professional services firms.

Introduction

This work investigates the growth of US professional services firms. This is addressed both theoretically and empirically using both longitudinal and cross sectional datasets made available by the US Census bureau. The analysis of the PSFs sector shows that firm age and size impact PSFs growth in a fashion similar to other firms in the US economy, notably in the US manufacturing. However it shows that PSFs, to the difference of US manufacturing firms, develops themselves by creating new establishment, which reveals a question of market access and business development. If this may hold for the main body of PSFs, those patterns can be challenged for large PSFs. I therefore build upon the Hopenhayn framework to show that assuming that PSFs invests in programs to increase their labor productivity can explain the observed growth patterns. This works contributes to two main strands of the related literature. First this works builds upon the studies that stress the effects of size and age on firms growth. Since the pioneering work of Gibrat [START_REF] Gibrat | Les inégalités économiques[END_REF] that stated that the relationship between firm growth and its initial size can be described by a constant growth rate and a random element that is normally distributed, it has been empirically shown that firm growth is industry specific and exhibits power law patterns ( [START_REF] Mansfield | Entry, gibrat's law, innovation, and the growth of firms[END_REF], [START_REF] Michael Hr Stanley | Scaling behaviour in the growth of companies[END_REF]). However size is not the only relevant variable to describe firms growth dynamics. Starting with [START_REF] David | Tests of alternative theories of firm growth[END_REF], authors have started to complement the classical firms growth dependency in size by observing a dependency in age. The relationship between firm development and its age has spurred an increased interest over the last years ( [START_REF] Coad | Firm age and performance[END_REF], [START_REF] Haltiwanger | Who creates jobs? small versus large versus young[END_REF]). This has motivated numerous applications on firms entry and exit rates [START_REF] Anyadike | All grown up? the fate after 15 years of a quarter of a million uk firms born in 1998[END_REF], [START_REF] André Van Stel | The relationship between start-up motive and earnings over the course of the entrepreneurs business tenure[END_REF], firms internationalization strategies [START_REF] Grazzi | Small, young, and exporters: New evidence on the determinants of firm growth[END_REF] and resilience to extreme random events [START_REF] Cowling | Did firm age, experience, and access to finance count? sme performance after the global financial crisis[END_REF].However most of the growth study focuses manufacturing and little has been done on the service industry [START_REF] Islam A Nassar | The validity of gibrat's law in developed and developing countries (2008-2013): Comparison based assessment[END_REF]. To my knowledge there has been no study on professional services firm industry (see [START_REF] Empson | Researching professional service firms: An introduction and overview[END_REF] for a definition), which therefore constitutes one of the main contribution of this paper. Second this study goes beyond the descriptive mechanisms related to firm growth to explore its drivers. The initial explanations in terms productivity and technology shocks of [START_REF] Hugo A Hopenhayn | Entry, exit, and firm dynamics in long run equilibrium[END_REF] gave birth to numerous discussion between product and firm life cycle [START_REF] Cooley | Aggregate consequences of limited contract enforceability[END_REF], [START_REF] Florin O Bilbiie | Endogenous entry, product variety, and business cycles[END_REF]. This has caught up with the evolutionary literature focusing on age as an important determinant by providing some evidence of the learning behavior [START_REF] Arkolakis | Firm learning and growth[END_REF] of firms towards their customer preferences [START_REF] Haltiwanger | Who creates jobs? small versus large versus young[END_REF]. Growing beyond considerations on demand, the field has had a special interest in financial considerations associated to growth and transformation. The main areas of interest have lately revolved around convex entry costs [START_REF] Pablo | Investment and firm dynamics[END_REF] and capital management [START_REF] Thomas | Financial markets and firm dynamics[END_REF], [START_REF] Luca | Entry, exit, firm dynamics, and aggregate fluctuations[END_REF]. With respect to the firm growth driver literature, this paper provides a discussion around the Hopenhayn framework by transforming the assumption of random exogenous shocks in productivity in an assumption of deterministic endogenous productivity investments. This is shown to depict properly the growth patterns of the largest US PSFs and can be used to infer that large US PSFs face an arms race with respect to their investments. From a technical standpoint, this works builds upon the mean field game literature [START_REF] Lasry | Mean field games[END_REF] (referred to as MFG thereafter). This type of structure indeed mixes population dynamics through Kolmogorov forward equations (see [START_REF] Gabaix | Power laws in economics and finance[END_REF] for a review in economy and [START_REF] Malevergne | Zipf's law and maximum sustainable growth[END_REF] for an example on firm growth) with profit optimization expressed as Bellman equations (see [START_REF] Perthame | Career plans and wage structures: a mean field game approach[END_REF] for an example). If this serves well mechanistic description, it also leads to interesting results associated to costs optimization (see [START_REF] Doumic | Toward an integrated workforce planning framework using structured equations[END_REF] for an example). If the MFG related literature is growing, there has been numerous discussion about its inherent general complexity. However when paired with economic concepts that are technically related to power laws [START_REF] Alexander I Saichev | Theory of Zipf 's law and beyond[END_REF], the MFG frameworks can be simplified to yield closed formulas that can easily be interpreted. This is applied in this paper to demonstrate how sectoral firm considerations can be used to estimate productivity investments. Note that one of the key features of the MFG apparatus lies in its capacity to indirectly tackle competition questions that appear to have yet eluded the evolutionary field in economic with respect to firm growth. MFGs indeed address actors behavior in a given environment. The MFG application developed in this paper therefore opens an interesting doorway for further research. The remaining of this paper is organized as follows. Section 2 provides a descriptive analysis of the growth patterns of US PSFs sector and shows that PSFs growth is notably spurred by the creation of new establishments in a dynamic of market expansion. Section 3 further explore the growth of top 25% largest US PSFs, where market expansion may no longer be relevant and proposes an explanation of the observed pattern through the notion of labor productivity improvement investments. Section 4 concludes.

Important Legal Remarks. The findings and opinions expressed in this paper are those of the authors and do not reflect any positions from any company or institution.

US PSFs growth patterns.

To describe US PSFs growth patterns, a longitudinal database made available by the US Census bureau is used. This database, which features are depicted in 2.1, covers 25 years of history and the overall US economy. It is used to discuss the application of Gibrat's law to PSFs in 2.2 at establishment level, prior to a discussion of the growth dynamics in the field in 2.3.

1975 -2000 Longitudinal data description.

To analyze the growth of PSFs, this study leverages the US synthetic longitudinal database ( [37]) (synLDB). The synLDB covers the entire US economy over a 25 years period starting 1975 and provides establishment level records of employment. If the database leverages the old standard industrial classification (SIC) system, mapping tables are provided by the US bureau of labor statistics to link the data to the recent NAICS system. Therefore, describing PSFs through the NAICS code 54 is similar to describing PSFs through the SIC code 731 (Advertising), 733 (Computer Programming, Data Processing, And Other Computer Related Services), 811 (Legal services) and 871 (Engineering, Architectural, And Surveying), 872 (Accounting, Auditing, And Bookkeeping Services), 873 (Research, Development, And Testing Services) and 874 (Management And Public Relations Services). Between 1975 and 2000, the PSFs establishments landscape has stayed relatively unchanged from a distribution standpoint (see figure [START_REF] Anyadike | All grown up? the fate after 15 years of a quarter of a million uk firms born in 1998[END_REF]) with an average of 8 to 9 employees per PSFs establishment. However its has grown steadily over time displaying an average 5.4% growth in the number of establishments (see figure [START_REF] Arkolakis | Firm learning and growth[END_REF]). Interestingly the growth in the number of establishments appears to decrease over time.

2.2 PSFs age and size growth patterns.

PSFs establishment growth.

From a growth perspective, the US establishments PSFs growth changes with size. This means that for small firms with one or two establishments (i.e firms below 100 employees as per table (9), Gibrat's law won't apply. Establishments indeed have their growth rate decreasing and stabilizing as they become bigger. Interestingly, large PSFs growth do appear to follow Gibrat's law. This seems to stem from the fact that they are composed of numerous establishments which growth rate has stabilized, so that the impact of new establishment launch is mitigated.As seen on tables (5) and (5), growth rates yet highly fluctuates on a year on a year basis. Establishments have different dynamics though depending on whether they belong to a firm with multiple establishments or they are the firms themselves. The main difference is that when a new establishment is not the first one in a firm, its initial growth speed is much faster (i.e the sourcing of the first 20 employees). From an age standpoint, professional services establishments growth rates decrease with age as seen on tables [START_REF] Chesher | Testing the law of proportionate effect[END_REF] and [START_REF] Bureau | Bls handbook of methods[END_REF]. Additionally establishments which are newly created within an existing firm have a higher growth rate than the one of a new established firm. Those findings are similar to the one on the US manufacturing sector, where it was found that size does not effect growth rate for large firms (i.e. Gibrat's law holds) but that growth rates slightly decrease with size for small firms [START_REF] Bronwyn | The relationship between firm size and firm growth in the us manufacturing sector[END_REF]. Additionally, the growth rate in manufacturing decreases with age [START_REF] David | Tests of alternative theories of firm growth[END_REF].

PSFs growth persistence.

Over the 25 years period covered in the US SynLDB, PSFs mono (resp. multi) establishments have grown by 16.1% (resp. 23.9%) on average (see figure [START_REF] Florin O Bilbiie | Endogenous entry, product variety, and business cycles[END_REF]). When testing the PSFs growth rate time series for auto and cross correlation, it came that the growth rate of mono establishments doesn't present autocorrelation nor cross correlation effects with the multi establishment time series. However the multi establishments time series presents strong autocorrelation patterns as seen in figure [START_REF] Blundell | Gmm estimation with persistent panel data: an application to production functions[END_REF]. As seen on the previous subsection, Gibrat's law doesn't hold for professional services establishments as the growth rate changes with the establishment size and age. Additionally, the analysis of the aggregated sector growth rate hints at a persistence in growth. To further investigate this characteristic, Gibrat's law was tested for growth persistence over time in a manner similar to [START_REF] Chesher | Testing the law of proportionate effect[END_REF] by running the following regression at establishment level from 1979 to 2000.

z u i,t = z u i,t-1 .β u + γ u .z u i,t-2 + u i,t
where t is an index for time, i is an index for the establishments, u an index differentiating mono versus multiple firm establishments, and z t,i is the deviation of the logarithm of the size of company i at time t from the mean of the logarithms of the sizes of companies at time t. With this test, if (β; γ) = (1, 0) then Gibrat's law is considered valid, otherwise, if β < 1 (resp. β > 1) large establishments are expect to grow more slowly (resp.rapidly) than their smaller counterparts. The regression was performed on a sample of establishments who were still operating to avoid sampling bias as suggested by [START_REF] Mansfield | Entry, gibrat's law, innovation, and the growth of firms[END_REF]. The estimation results are displayed in table [START_REF] Thomas | Financial markets and firm dynamics[END_REF].

There are three main results on this analysis. First, whatever the their type, large establishments grow slower than smaller ones as β 0 < 1 and β 1 < 1. As PSFs are a collection of establishments, this means that smaller firms grow faster than their larger counterparts. Second, as suggested in the aggregated analysis, establishments that belong to a firm that has multiple location display a positive persistence in growth, which means they grow faster in the second period if they have growth faster in the past. However this is not something that holds for firms with only one establishment.Third Gibrat's law validity appears to fluctuate with the year, which could suggest that PSFs establishments are subject to a business cycle. From a comparison standpoint, the observed US PSFs persistence pattern in growth slightly differ from the current findings on US manufacturing. First of all, as discussed in [START_REF] Bottazzi | Common properties and sectoral specificities in the dynamics of us manufacturing companies[END_REF], US manufacturing establishments do not overall exhibit persistence in their grow patterns. However, when focusing on the largest firms with heavy investments in R&D, persistence in growth has been empirically demonstrated for US manufacturing firms [START_REF] Blundell | Gmm estimation with persistent panel data: an application to production functions[END_REF]. This could mean that the persistence in growth between the two industries is triggered by different mechanisms. While in the PSFs case, persistence could be linked to geographical expansion, persistence in manufacturing growth could be linked to heavy capital investments. This is consistent with macro level findings that state that at an aggregated U.S economy level, business growth is the sum of a cyclic component and a deterministic trend [START_REF] Hashem Pesaran | Persistence, cointegration, and aggregation: A disaggregated analysis of output fluctuations in the us economy[END_REF].

2.3 PSFs growth through geographical expansion.

US PSFs growth dynamics at an aggregated level.

When looking at the 2007 Statistics of US businesses (SUBS), it has to be noted that large firms have many establishments and that those establishments are larger than the ones of smaller firms. As seen in table [START_REF] Coad | Firm age and performance[END_REF] PSFs have only one establishment until they reach 20 employees, then additional locations appear as the firms become bigger. On average, it seems that establishments grow until about 50 to 60 people at which point other branches are appearing. Note that on average in the US, PSFs do have 1 establishment and employ about 9 to 10 persons. The professional services sector appears to have an establishment composition with firm size that differs from the overall economy and the manufacturing sector (see table [START_REF] Coad | Firm age and performance[END_REF]). Compared to the manufacturing sector, PSFs have much more establishments with suggest a geographical expansion associated to firm growth. However this growth pattern is not as strong in the PSFs than in the overall US economy. This would suggest that PSFs have to trade off between consolidating knowledge in hubs to enjoy returns to scale and winning local markets. Large PSFs indeed have a maximum number of establishments of 12 to 48 when their size grows above 2500 employees, while manufacturing have a maximum of 5 to 10 establishments. On the other hand, large firms in the overall US economy have between 60 to 400 locations. PSFs do not concentrate employment in a specific location to the difference of manufacturing. As seen in table [START_REF] Cooley | Aggregate consequences of limited contract enforceability[END_REF], PSFs have an employment footprint per establishment that is similar to the overall US economy. To this extent, large firms with more than 2000 people do employ an average of 65 persons per establishment, while the manufacturing sector concentrates between 180 and 280 persons at a same location. Therefore if PSFs growth is generated by the growth of its existing establishments, growth pattern at firms level may also be impacted by the development of a firm's establishment network. This may especially be true for small and medium firms (i.e < 250 employees). An additional analysis would be required here that is out of scope of this paper.

PSFs establishments: creation and destruction patterns.

To dig deeper in the geographical expansion of PSFs, it is important to look at the entry and exiting dynamics of PSFs at establishment level. From an entry standpoint, there a two main cases: either PSFs are created by single individuals from scratch or establishments are created as part of an existing firm with an average of 12 employees. The first option is more likely as seen in table [START_REF] Cowling | Did firm age, experience, and access to finance count? sme performance after the global financial crisis[END_REF] and over the 1976 to 2000 period there as been an average there is a 12% of newly created PSFs with on single employee. On the other hand, there is chance for a firm to open a new establishment 9%. Note that while the entry of single employee firm has been quite stable proportion wise, the market expansion through the opening of new establishment appears to fluctuate. Note that on average, the entry rate of new establishments has lowered over time as market expansion has slowed for existing firms.

From an exit pattern standpoint, PSFs establishments that disappear have 2 to 3 employees when belonging to single establishment firm or about 12 when belonging to an existing firm with multiple establishment. This suggests that PSFs cease their activity is they fail to launch. Overall about 7% of the PSFs establishments close on a year on year basis. Establishments belonging a a multi establishment firms are more robust as their exiting probability has been of 4.5% on average for the period, versus a 7.4% chance for the mono establishment counterparts. Interestingly, most of variance in the exit patterns is associated to the multi establishment firms. It therefore seems that the risk of a successful market expansion even though limited compared to the creation of a firm from scratch may be fluctuating because of markets evolutions. Finally, from a dynamic point of view, both size and age have a clear effect on the PSFs landscape establishment exiting probability, as seen on table [START_REF] Luca | Entry, exit, firm dynamics, and aggregate fluctuations[END_REF].The observed entry and exit patterns with respect to age, size and the difference between mono and multi establishments are consistent with the overall findings on US manufacturing [START_REF] Dunne | Patterns of firm entry and exit in us manufacturing industries[END_REF] and the rest of the US [START_REF] Rossi-Hansberg | Establishment size dynamics in the aggregate economy[END_REF] and developed economy in general [START_REF] Paul | What do we know about entry[END_REF]. The interesting element at stake is that when looking at entry and exit patterns of US manufacturing, one of the main elements is to consider whether or not the plants are being invested in and are diversifying, which spurs a reduction in exit rates. For PSFs, the driver of exit rate reduction appears slightly different as PSFs that expand geographically are more likely to survive. The firm growth patterns investigated in 2 depict most of the PSFs in the US, which are small to medium firms (i.e¡ 250 employees) composed of 1 to 4 establishments on average that expand geographically. However as those findings can be challenged for large firms as they do exhibit as they are composed of multiple established locations. Therefore geographical expansion and the addition of another 10 to 20 employees, may not be the main drivers at stake.

3 Large US PSFs growth.

To discuss what could be happening for large firms, an option is to leverage the traditionally notion of firm productivity, which has been used in the literature to explain firms growth patterns. In section (3.1), the Hopenhayn's framework [START_REF] Hugo A Hopenhayn | Entry, exit, and firm dynamics in long run equilibrium[END_REF] is revisited by replacing the hypothesis of exogenous random productivity shocks by endogenous productivity investment. The US QWI and US BLS dataset are then used to calibrate US PSFs productivity improvement investments in (3.2). Finally, in (3.3), productivity investment are tied to the distribution of PSFs in the case of large and high productive firms (i.e. > 250 employees) to discuss the implications of the framework in terms of growth speed with age and size, with the help of a firm level cross sectional dataset made available by the US census.

A toy model of firm growth

Similarly to [START_REF] Hugo A Hopenhayn | Entry, exit, and firm dynamics in long run equilibrium[END_REF], assume that PSFs supply an output q(φ, p, ω) (i.e. charged hours) to their clients, where φ is the PSF productivity (i.e utilization), ω the price of labor and p the price of the output (i.e the hourly rate). On the other hand, PSFs demand labor up to a quantity n(φ, p, ω) (i.e. number of employees), so that PSFs profit π can be expressed as:

π(φ, p, ω) = p.q(φ, p, ω) -ω.n(φ, p, ω)

Assume that PSFs are price taking, so that p t , ω t are known and driven by the overall market demand. The problem for the firm is to change its productivity φ t at a speed ρ(t)

(dφ t = ρ(φ)dt) v(t, φ t ) = max ρ(.) t+T t e -β.x (π(φ x , p x , ω x ) -c(ρ(x))dx (1)
where c(ρ(x)) represents the investment costs to change the firm productivity. According to ( [START_REF] Inklaar | Market services productivity across europe and the us[END_REF]), those investments in the case of professional services are mainly of two nature: investments in information and communication technology (ICT) and human capital. Assume that the price for the output is driven by the overall market demand p = D(Q) where Q = q(φ, p, ω).f (φ)dφ and that the price for the overall labor supply is driven by the firms demand ω = W (N ) where N = n(φ, p, ω)f (φ)dφ. Following the US bureau of labor statistics (BLS) methodology ( [START_REF] Bureau | Bls handbook of methods[END_REF]), productivity in professional services is defined as the amount of output divided by the amount of labor required to deliver this output.

In the case of services, the BLS recommendation is to measure output is as the dollar amount of sold of services adjusted for price changes in these over time. In the case of PSFs, the productivity can therefore be interpreted as number of charged hours (i.e. output) sold at a given price p that have been generated by a hour of work (i.e labor input). With respect to the previously developed notation, this means that:

φ t = q(φ t , p t , ω t ) n(φ t , p t , ω t )
Under this framework, competition between firms is in essence a competition in terms of productivity.

Intuitively at equal level of supply (i.e equal number of employees), PSFs with the highest productivity will be able to generate more profit and to increase their size faster than their competitors. This will be shown in the case of highly productive firms in (3.3). Note that under this view of productivity, comparing the productivity from one firm to another in a given year is not that straight forward. Sales and labor are easily obtained at firm level, but the output in terms of charged hours or the firm levels differences in terms of hourly rates are not accessible. Therefore in this particular study, the output in terms of charged hours is defined at a constant market price p which is shared across all firms. So that if firm A and firm B are selling the same amount of charged hours and have the same amount of workers, if firm A is selling services at a price above the price of firm B, its productivity will be higher. From the overall US level yearly time series provided by the US census, demand and supply curves can be modeled in a linear fashion for the sake of simplicity (see table [START_REF] Anyadike | All grown up? the fate after 15 years of a quarter of a million uk firms born in 1998[END_REF]): Assuming that q(φ, 0, ω) = 0 and n(φ, p, 0) = 0 and leveraging the previous definition of productivity as q/n = φ, the linear models imply that:

Insights from US macro level data

ω = W (N ) = a.N + b p = D(Q) = h.Q + k Slope Intercept
n(φ, p, ω) = A(φ).ω q(φ, p, ω) = A(φ).φ.p A(φ).(φ -1)f (φ)dφ = 0 Firm (resp. PSFs sector) level profit π (resp. Π) therefore can be simplified as:

π = A(φ).(φ.p 2 -ω 2 ) Π = ( A(φ)f (φ)dφ).(p 2 -ω 2 )
The optimization program (see eq. ( 1)) of the firm can therefore be expressed as:

∂ t v + max ρ (ρ.∂ φ v -β.v + π(φ, p, ω) -c(ρ)) = 0 (2)
This leads to ∂ φ v = ∂ ρ c. As we assume that ρ is dependent in φ, it is possible to look for a form on productivity improvement investment ĉ(φ) = c(ρ(φ)). Further differentiating (2) towards φ yields:

∂ φφ ĉ - β ρ ∂ φ ĉ = ∂ φ (π) ρ.∂ ρ .φ (3)
3.3 Link to the overall firm size -age growth patterns of high productive large PSFs?

Leveraging cross sectional data for the 2007 SBO public micro sample from the US Census, it is easy to access the distribution f (φ) (resp. the exiting firm distribution T (φ)). Interestingly figure [START_REF] Chesher | Testing the law of proportionate effect[END_REF] shows that new entrants have a lower productivity than existing firms. However this also shows that firms that cease operations have the same productivity than the existing ones. Although this may not be a consistent pattern over the year, this could suggest that the exit criteria based on productivity shocks detailed by [START_REF] Hugo A Hopenhayn | Entry, exit, and firm dynamics in long run equilibrium[END_REF] may not be the drivers for PSFs. When zooming on the top 25% most productive PSFs, extreme value theory yields that those firm are distributed in productivity according to a power law which exponent can be estimated according to the log rank technique of [START_REF] Gabaix | Power laws in economics and finance[END_REF]. Results and graphical representations are displayed in table [START_REF] Arkolakis | Firm learning and growth[END_REF] and figures ( 8) and ( 9). The main characteristic of those top firms is that they generate more than 107k$ per year and per employee, so that with an overall market price of 81.34$/h as per the 2007 US QSS, those firms are actually able to charge 1320 hours per employee per year on average. Knowing that an average year of work represents 260 days of 8 hours each, this means that top PSFs firms have on average a productivity/utilization above 50%. Taking a population dynamic point of view, and assuming that the PSFs distribution is time invariant, it comes that:

∂ φ (ρ.f ) = T (φ) f (s)
T (s) Pareto Exponent -1.86 (0.0004) -1.81 (0.003) R 2 99.9% 99.6%

Table 2: Highly productive PSFs -Pareto exponent Estimation.

As the top productive firms are distributed and exiting according to a power law, it comes that the productivity growth speed follows ρ(φ) = ρ 0 .φ γ , where γ ≈ 0.95 according on the estimates displayed in table [START_REF] Arkolakis | Firm learning and growth[END_REF]. When looking at the PSFs landscape dynamics and further restricting the scope to the high productive firm that employs more than 250 employees and investigating the problem under a firm size lens n, it comes that under the assumptions of time invariance:

∂ n (µ.m) = θ(n)
Where m (resp. θ) represents the distribution of existing (resp. entry )firm according to their size and µ is the growth speed of a firm (i.e dn t = µ(n t )dt). Restricting the analysis to the top productive large firms (i.e the distribution tail) leads to a power law form for m, θ and therefore µ = µ 0 .n Ψ . The power exponents can be estimated from the SBO PUMS sample, which features 760 large high productive high firms (i.e about 0.7% of the overall SBO sample),using [START_REF] Gabaix | Power laws in economics and finance[END_REF]. The estimation results shown in table (3) yield a value for Ψ ≈ 1. Interestingly this means, that for the subset of large PSFs with high productivity, Gibrat's law could hold.

m(s) θ(s) Pareto Exponent -2.35 (0.015) -2.43 (0.25) R 2 96.9% 90.8% Table 3: Highly productive large PSFs -Pareto exponent Estimation.

For a given firm, assuming the market is at the equilibrium (ω t = Cstte), the framework developed in the previous section entails that:

dn t = dφ t A (φ).ω.dt = µ(n t )dt ↔ ρ(φ).A (φ).ω = µ(A(φ).ω)
On the high end (i.e. large φ, large n), this means that:

dA(φ) A(φ) = ( µ 0 .ω ρ 0 .φ γ )dφ ↔ A(φ) = A(φ min ).e ([ µ 0 .ω ρ 0 .(1-γ) φ 1-γ ] φ φ min
From a firm growth standpoint; for high productive large PSFs, as dφ = ρ 0 .φ γ .dt and dn = µ 0 .n.dt, it comes that:

φ(t) = ((1 -γ).ρ.t + φ 1-γ 0 ) 1 1-γ n(t) = n 0 .e µ 0 .t
This complements the empirical findings of the synLBD, where establishment growth speed decrease with age and size. The rate of change in growth speed observed in the tables ( 4) and ( 6) indeed show that for old and large PSFs, the growth rate may not change much. It is then possible to conclude the analysis by pushing the framework with respect to productivity investments. As ρ = ρ 0 .φ γ ↔ φ = ( ρ ρ 0 ) 1 γ , the equation (3) becomes:

∂ φφ ĉ - β ρ 0 .φ γ ∂ φ ĉ = γ. A (φ)(φ.p 2 -ω 2 ) + A(φ)p 2 φ
Assume for the sake of simplicity that β = 0 (i.e. large firms do not discount their future cash flows). Asymptotically, for highly productive large firms A(φ) = Âe φ 1-γ . Under this approximation, it comes that:

∂ φφ ĉ = γ.( Â.(1 -γ).φ -γ .e φ 1-γ (φ.p 2 -ω 2 ) + Â.e φ 1-γ .p 2 ) > 0
As the higher the productivity, the higher the value of the firm ∂ φ v = ∂ ρ c > 0, it comes that ∂ φ c has the same sign as γ > 0.This means that for large highly productive PSFs, productivity investments are increasing with productivity to stay competitive and that their followers tend to invest relatively less and mimic the leader, so that there is a form of arms race competition.

Conclusion.

In this paper, both a longitudinal empirical dataset ranging from 1977 to 2000 and a cross-sectional dataset from 2007 from the US were used to show that within the professional services sector Gibrat's law doesn't hold. PSFs growth is indeed slowing both at equal size with age and with size. Directionally those findings are similar to the US manufacturing sector, which constitutes the main reference in the study of US firms growth. If the size and age patterns disappear for large PSFs, growth becomes persistent and exhibits 3 to 4 years cycles, which appears different from the manufacturing field. Additionally, while manufacturing firms concentrate their growth on a handful of localization, PSFs growth, similarly to the average US firm growth appears driven by geographical expansion.

As for large PSFs, the pattern of geographical expansion can be challenged, this paper proposes a variant to the classical Hopenhayn framework. It notably shows that, in the context of large and highly productive US PSFs, growth patterns can be explained by investment policies in labor productivity improvement programs. Those investments come as a form of arms race between the top competitors. Finally it could be interesting to further detail the developed considerations in terms of productivity. It could be for example worth to investigate what is happening for firms with a low productivity and perhaps to challenge their goals as profit maximizing organization by leveraging some of considerations of [START_REF] Herbert | Organizations and markets[END_REF]. Additionally, productivity investments from firms do not yield deterministic results. The framework proposed in section (3.1) could be easily randomized and numerical simulations performed. Eventually, the main next step of this analysis would be to drill down at a firm level to investigate the main productivity levers a PSF can pull (e.g. developing a multi local model and/or investing in automation, increasing its market access capabilities through its quality perception/reputation (see [START_REF] Ribes | What drives professional services firms prices?[END_REF])).

Compliance with Ethical Standards: This study has received no fundings.

5 Appendix.

Figure 1 :

 1 Figure 1: US PSFs establishment size distribution.Figure 2: US PSFs establishment number growth.

Figure 2 :

 2 Figure 1: US PSFs establishment size distribution.Figure 2: US PSFs establishment number growth.

Figure 3 :

 3 Figure 3: US Establishments growth rate evolution.

Figure 4 :

 4 Figure 4: US multi establishments growth rates autocorrelation.

  To estimate the PSFs demand function W in workers, the US quarterly workforce indicators (QWI) provided by the US Census were used over the last 10 years (i.e 2007 to 2017). N was taken as the total number of US PSFs employees (where PSFs are defined according to the NAICS code 54), while labor prices ω have been taken as the average monthly wage in the sector. On the other hand, the estimation of the PSFs supply function D was achieved by leveraging data from the quarterly service survey (QSS) provided by the US census. Price p where estimated as the total PSFs revenue divided by the total number of worked hours, while the output Q was associated to the number of worked hours. Results are displayed in figures (5) and[START_REF] Bureau | Bls handbook of methods[END_REF].

Figure 5 :

 5 Figure 5: US PSFs demand curve.Figure6: US PSFs supply curve.

Figure 6 :

 6 Figure 5: US PSFs demand curve.Figure6: US PSFs supply curve.

Figure 7 :

 7 Figure 7: US PSFs productivity repartitions.

Figure 8 :

 8 Figure 8: Top productive PSFsdensity estimation.

Figure 9 :

 9 Figure 9: Top productive exiting PSFs -density estimation.

Table 1 :

 1 US PSFs Demand and Supply Curve estimates

	R 2

Table 4 :

 4 US PS multi establishments Growth rate with size between 1977 & 2000 [%].

	Year/ Size 0	1	2	[3;5] [6;10] [11;20] 20+ Total
	1977	-9.42 10.0 2.89 2.18	1.03	0.77
	1982	-4.69 2.66 2.41 2.05	1.79	1.28
	1987	-4.94 3.67 2.94 2.44	1.77	1.60
	1992	-5.02 2.95 1.92 1.89	1.33	1.03
	1997	-1.65 1.66 1.32 1.20	1.15	0.96
	Year/ Size 0	1	2	[3;5] [6;10] [11;20] 20+ Total
	1977	-3.03 2.66 2.12 0.46	0.71	0.12
	1982	-2.32 1.70 1.06 0.97	1.07	1.05
	1987	-2.58 1.36 1.04 1.09	0.77	0.94
	1992	-2.43 1.24 1.16 0.94	0.89	0.66
	1997	-1.13 1.02 1.05 1.05	1.09	1.07

Table 5 :

 5 US PS mono establishments Growth rate with size between 1977 & 2000 [%].

	Size / Age	0	1	2	3	[4;5] [6;7] [8;11] [12;14] [15;20] 21+ Total
	0	-	-	-	-	-	-	-	-	-	-	-
	1	1.78 1.89 1.88 1.84 1.82 1.92 2.04	2.19	1.80	2.04	1.89
	2	1.60 1.62 1.61 1.57 1.54 1.55 1.68	1.77	1.61	1.77	1.62
	[3;5]	1.50 1.49 1.44 1.39 1.40 1.49 1.47	1.56	1.48	1.40	1.46
	[6;10]	1.34 1.32 1.31 1.30 1.25 1.29 1.32	1.34	1.30	1.27	1.31
	[11;20]	1.24 1.25 1.25 1.25 1.22 1.21 1.22	1.20	1.14	1.16	1.21
	20+	1.02 1.01 1.03 1.03 1.01 1.02 1.03	1.03	1.03	1.00	1.02
	Total	1.34 1.31 1.30 1.28 1.24 1.24 1.23	1.24	1.17	1.14	-

Table 6 :

 6 Average US PSFs multi establishments 1998 to 1999 Growth rate with size and age [%].

	Size / Age	0	1	2	3	[4;5] [6;7] [8;11] [12;14] [15;20] 21+ Total
	0	-	-	-	-	-	-	-	-	-	-	-
	1	1.29 1.26 1.26 1.25 1.25 1.24 1.24	1.23	1.23	1.22	1.25
	2	1.11 1.10 1.11 1.10 1.09 1.09 1.09	1.09	1.07	1.06	1.09
	[3;5]	1.06 1.06 1.06 1.06 1.05 1.05 1.05	1.04	1.03	1.04	1.05
	[6;10]	1.03 1.02 1.03 1.05 1.04 1.03 1.03	1.01	1.01	1.01	1.03
	[11;20]	1.09 1.07 1.06 1.06 1.06 1.05 1.04	1.02	1.01	1.00	1.03
	20+	1.16 1.21 1.22 1.17 1.16 1.10 1.06	1.06	1.03	1.02	1.08
	Total	1.20 1.18 1.18 1.18 1.17 1.15 1.14	1.13	1.12	1.11	-

Table 7 :

 7 Average US PSFs mono establishments 1998 to 1999 Growth rate with size and age [%].

	Size / Age	0	1	2	3	[4;5]	[6;7] [8;11] [12;14] [15;20] 21+ Total
	0	9.43 26.36 20.86 20.32 19.32 16.44 15.22 13.13	11.98 6.57	9.15
	1	4.11 12.41 18.43 17.64 16.50 13.95 12.80 11.01	10.04 5.52	9.71
	2	3.92 11.49 16.31 14.96 13.55 11.71 10.22	8.84	8.04	4.37	7.49
	[3;5]	3.48 10.89 15.28 12.87 12.01 10.38 8.86	7.46	6.97	4.05	6.42
	[6;10]	3.07 9.76 12.94 10.98 10.17 8.16	7.46	6.82	6.10	3.51	4.91
	[11;20]	2.74 7.79	9.76	8.60	8.19	6.30	5.90	5.35	4.44	2.87	3.18
	20+	2.20 5.66	7.57	8.33	8.12	5.99	5.17	4.87	4.28	2.14	2.50
	Total	6.04 14.54 14.40 12.63 12.87 10.76 9.54	8.27	7.52	4.04

Table 8 :

 8 Average US PSFs establishments exit probability [%] between 1975 and 2000.

	Firm Size Establishment(s) Establishment(s) Establishments(s)
		per PSFs	per Manufacturing	per firm
			firm	in the US
	0-4	1	1	1
	5-9	1	1	1,01
	10-14	1,03	1	1,08
	15-19	1,06	1	1,04
	20-24	1,11	1	1,13
	25-29	1,20	1,03	1,18
	30-34	1,23	1,04	1,25
	35-39	1,28	1,05	1,31
	40-44	1,31	1,08	1,35
	45-49	1,44	1,07	1,45
	50-74	1,56	1,11	1,61
	75-99	1,82	1,19	2,03
	100-149	2,19	1,27	2,59
	150-199	2,86	1,44	3,57
	200-299	3,28	1,67	4,94
	300-399	4,03	1,91	6,81
	400-499	4,73	1,98	9,26
	500-749	5,25	2,42	11,76
	750-999	6	3,26	15,48
	1,000-1,499	6,03	3,58	20,33
	1,500-1,999	6,78	3,62	29,97
	2,000-2,499	9,09	5,41	38,62
	2,500-4,999	12,21	5,59	61,61
	5,000 +	48,07	9,39	405,24
	Total	1,1	1,21	1,27

Table 9 :

 9 Number of US establishments with firm size according to 2007 SUBS in professional services, manufacturing and the overall economy.

	Firm Size	Employment per	Employment per	Employment per
		PSFs establishment Manufacturing establishment US establishment
	0-4	1,50	1,85	1,65
	5-9	6,49	6,68	6,49
	10-14	11,36	11,68	11,20
	15-19	15,67	16,56	15,45
	20-24	19,44	21,35	19,19
	25-29	22,05	25,88	22,56
	30-34	25,34	30,13	25,45
	35-39	28,20	34,51	27,90
	40-44	31,06	38,46	30,87
	45-49	31,78	42,75	32,12
	50-74	36,81	52,45	37,14
	75-99	43,78	67,4	42,04
	100-149	49,11	85,64	46,49
	150-199	49,80	97,78	48,12
	200-299	57,51	115,88	49,01
	300-399	62,90	127,54	50,58
	400-499	63,23	141,25	48,12
	500-749	66,45	143,95	51,54
	750-999	71,64	157,28	55,7
	1,000-1,499	93,63	162,92	59,89
	1,500-1,999	85,78	171,48	57,64
	2,000-2,499	66,85	167,52	57,56
	2,500-4,999	63,67	176,52	56,53
	5,000 +	53,24	269,88	50,55
	Total	9,48	40,19	15,65

Table 10 :

 10 Employment per establishment with firm size according to 2007 SUBS in professional services, manufacturing and the overall economy.

	Period	β 0	γ 0	β 1	γ 1
	1977 -1979 0.76 (0.05)	0.14 (0.05)	0.90 (0.04)	0.03 (0.05)
	1978 -1980 0.94 (0.05)* 0.04 (0.05)* 0.85 (0.05)	0.11 (0.05)
	1979 -1981 0.93 (0.05) -0.03 (0.06) 0.83 (0.06)	0.10 (0.06)
	1980 -1982 0.80 (0.05)	0.03 (0.05)	0.92 (0.04)	0.06 (0.04)
	1981 -1983 0.94 (0.05)	0.06 (0.04)	0.88 (0.05)	0.11 (0.06)
	1982 -1984 0.81 (0.05)	0.08 (0.06) 0.95 (0.05)* 0.05 (0.05)*
	1983 -1985 0.87 (0.05)	0.05 (0.05)	0.89 (0.04)	0.07 (0.05)
	1984 -1986 0.93 (0.05) -0.01 (0.05) 0.83 (0.05)	0.16 (0.05)
	1985 -1987 0.80 (0.05)	0.09 (0.05)	0.76 (0.05)	0.18 (0.06)
	1986 -1988 0.77 (0.05)	0.09 (0.05)	0.90 (0.05)	0.11 (0.05)
	1987 -1989 0.80 (0.05)	0.13 (0.05)	0.88 (0.05)	0.12 (0.05)
	1988 -1990 0.87 (0.05)	0.00 (0.05) 0.95 (0.05)* 0.02 (0.05)*
	1989 -1991 0.87 (0.05)	0.00 (0.05)	0.91 (0.05)	0.04 (0.005)
	1990 -1992 0.86 (0.05)	0.09 (0.05)	0.80 (0.05)	0.15 (0.05)
	1991 -1993 0.87 (0.05) -0.02 (0.05) 0.96 (0.04)* 0.00 (0.04)*
	1992 -1994 0.82 (0.05)	0.08 (0.05) 0.94 (0.05)* 0.04 (0.05)*
	1993 -1995 0.80 (0.05)	0.16 (0.05) 0.95 (0.05)* -0.01 (0.05)*
	1994 -1996 0.88 (0.05)	0.02 (0.05)	0.89 (0.04)	0.08 (0.04)
	1995 -1997 0.78 (0.05)	0.12 (0.05)	0.89 (0.06)	0.07 (0.06)
	1996 -1998 0.87 (0.05)	0.12 (0.05)	0.91 (0.04)	0.04 (0.04)
	1997 -1999 0.93 (0.05) -0.01 (0.05) 0.98 (0.05)* -0.02 (0.05)*
	1998 -2000 0.89 (0.05)	0.07 (0.05)	0.89 (0.05)	0.07(0.05)

Table 11 :

 11 Gibrat Law -Persistence test at establishment level.

	Year	N 0 t	N 1 t	E 0 t	E 1 t	N 0 e,t	N 1 e,t	E 0 e,t	E 1 e,t	N 0 c,t	N 1 c,t	E 0 c,t	E 1 c,t

Acknowledgments. I would like to thank Pierre Fleckinger for his continuous support throughout this work. I also would like to thank Angelo Secchi for his time and insightful remarks as well as two anonymous referees.

Conflict of Interest:

The author declare that he has no conflict of interest.


Ulrike Von Luxburg
email: ulrike.luxburg@tuebingen.mpg.de

Sébastien Bubeck
email: sebastien.bubeck@inria.fr

Stefanie Jegelka
email: stefanie.jegelka@tuebingen.mpg.de

Michael Kaufmann

Consistent Minimization of Clustering Objective Functions

Clustering is often formulated as a discrete optimization problem. The objective is to find, among all partitions of the data set, the best one according to some quality measure. However, in the statistical setting where we assume that the finite data set has been sampled from some underlying space, the goal is not to find the best partition of the given sample, but to approximate the true partition of the underlying space. We argue that the discrete optimization approach usually does not achieve this goal. As an alternative, we suggest the paradigm of "nearest neighbor clustering". Instead of selecting the best out of all partitions of the sample, it only considers partitions in some restricted function class. Using tools from statistical learning theory we prove that nearest neighbor clustering is statistically consistent. Moreover, its worst case complexity is polynomial by construction, and it can be implemented with small average case complexity using branch and bound.

Introduction

Clustering is the problem of discovering "meaningful" groups in given data. Many algorithms try to achieve this by minimizing a certain quality function Q n , for example graph cut objective functions such as ratio cut or normalized cut, or various criteria based on some function of the within-and between-cluster similarities. The objective of clustering is then stated as a discrete optimization problem. Given a data set X n = {X 1 , . . . , X n } and a clustering quality function Q n , the ideal clustering algorithm should take into account all possible partitions of the data set and output the one that minimizes Q n . The implicit understanding is that the "best" clustering can be any partition out of the set of all possible partitions of the data set. The algorithmic challenge is to construct an algorithm which is able to find this clustering. We will call this approach the "discrete optimization approach to clustering".

If we look at clustering from the perspective of statistical learning theory we assume that the finite data set has been sampled from an underlying data space X according to some probability measure. The ultimate goal in this setting is not to discover the best possible partition of the data set X n , but to learn the "true clustering" of the underlying space. In an approach based on quality functions, this "true clustering" can be defined easily. We choose a clustering quality function Q on the set of partitions of the entire data space X , and define the true clustering f * to be the partition minimizing Q. In this setting, a very important property of a clustering algorithm is consistency. Denoting the clustering constructed on the finite sample by f n , we require that Q(f n ) converges to Q(f * ) when n→∞. The most important insight of statistical learning theory is that in order to be consistent, learning algorithms have to choose their functions from some "small" function space only. To measure the size of a function space F one uses the quantity N F (x 1 , .., x n ) which denotes the number of ways in which the points x 1 , . . . , x n can be partitioned by functions in F. One can prove that in the standard setting of statistical learning theory, a necessary condition for consistency is that ❊ log N F (x 1 , .., x n )/n → 0 (cf. Theorem 2.3 in Vapnik, 1995, Section 12.4 of [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF].

Stated like this, it becomes apparent that the two viewpoints described above are not compatible with each other. While the discrete optimization approach on any given sample attempts to find the best of all (exponentially many) partitions, the statistical learning theory approach restricts the set of candidate partitions to have sub-exponential size. Hence, from the statistical learning theory perspective, an algorithm which is considered ideal in the discrete optimization setting is likely to overfit. One can construct simple examples (cf. [START_REF] Bubeck | Overfitting of clustering and how to avoid it[END_REF] which show that this indeed can happen: here the partitions constructed on the finite sample do not converge to the true clustering of the data space. In practice, for most cases the discrete optimization approach cannot be performed perfectly as the corresponding optimization problem is NP hard. Instead, people resort to heuristics. One approach is to use local optimization procedures potentially ending in local minima only (this is what happens in the k-means algorithm). Another approach is to construct a relaxation of the original problem which can be solved efficiently (spectral clustering is an example for this). In both cases, one usually cannot guarantee how close the heuristic solution is to the global finite sample optimum. This situation is clearly unsatisfactory: for most clustering algorithms, we neither have guarantees on the finite sample behavior of the algorithm, nor on its statistical consistency in the limit.

The following alternative approach looks much more promising. Instead of attempting to solve the discrete optimization problem over the set of all partitions, and then resorting to relaxations due to the NP-hardness of this problem, we turn the tables. Directly from the outset, we only consider candidate partitions in some restricted class F n containing only polynomially many functions. Then the discrete optimization problem of minimizing Q n over F n is no longer NP hard -it can trivially be solved in polynomially many steps by trying all candidates in F n . From a theoretical point of view this approach has the advantage that the resulting clustering algorithm has the potential of being consistent. In addition, it also leads to practical benefits: rather than dealing with uncontrolled relaxations of the original problem, we restrict the function class to some small enough subset F n of "reasonable" partitions. Within this subset, we then have complete control over the solution of the optimization problem and can find the global optimum. Put another way, one can also interpret this approach as some controlled way of sparsifying the NP hard optimization problem, with the positive side effect of obeying the rules of statistical learning theory.

Nearest neighbor clustering

In the following we assume that we are given a set of data points X n = {X 1 , . . . , X n } and pairwise distances d ij = d(X i , X j ) or pairwise similarities s ij = s(X i , X j ). Let Q n be the finite sample quality function to optimize on the sample. To follow the approach outlined above we have to optimize Q n over a "small" set F n of partitions of X n . Essentially, we have three requirements on F n : First, the number of functions in F n should be at most polynomial in n. Second, in the limit of n → ∞ the class F n should be rich enough to approximate any measurable partition of the underlying space. Third, in order to perform the optimization we need to be able to enumerate all members of this class, that is the function class F n should be "constructive" in some sense. A convenient choice satisfying all those properties is the class of "nearest neighbor partitions". This class contains all functions which can be generated as follows. Fix a subset of m ≪ n "seed points" X s1 , . . . , X sm among the given data points. Assign all other data points to their closest seed points, that is for all j = 1, . . . , m define the set Z j as the subset of data points whose nearest seed point is X sj . Then consider all partitions of X n which are constant on the sets Z j . More formally, for given seeds we define the set F n as the set of all functions f : X → {1, . . . , K} which are constant on the cells of the Voronoi partition induced by the seeds. Here K denotes the number of clusters we want to construct. The function class F n contains K m functions, which is polynomial in n if the number m of seeds satisfies m = O(log n). Given F n , the simplest polynomial-time optimization algorithm is then to evaluate Q n (f ) for all f ∈ F n and choose the solution f n = argmin f ∈Fn Q n (f ). We call the resulting clustering the nearest neighbor clustering and denote it by NNC(Q n ). In practice, the seeds will be chosen randomly among the given data points.

Consistency of nearest neighbor clustering

In this section we prove that nearest neighbor clustering is statistically consistent for many clustering quality functions. Due to the complexity of the proofs and the page restriction we can only present sketches of the proofs. All details can be found in von [START_REF] Luxburg | Supplementary material to "Consistent minimization of clustering objective functions[END_REF]. Let us start with some notation. For any clustering function f : | d → {1, . . . , K} we denote by the predicate A(f ) a property of the function which can either be true or false. As an example, define A(f ) to be true if all clusters have at least a certain minimal size. Moreover, we need to introduce a predicate A n (f ) which will be an "estimator" of A(f ) based on the finite sample only. Let m := m(n) ≤ n be the number of seeds used in nearest neighbor clustering. To simplify notation we assume in this section that the seeds are the first m data points; all results remain valid for any other (even random) choice of seeds. As data space we use X = | d . We define:

NN m (x) := NN m(n) (x) := argmin y∈{X1,...,Xm} x -y ( for x ∈ | d ) F := {f : | d → {1, . . . , K} | f continuous P-a.e. and A(f ) true} F n := F X1,...,Xn := {f : | d → {1, . . . , K} | f satisfies f (x) = f (NN m (x)), and A n (f ) is true} F n := X1,...,Xn∈| d F X1,...,Xn
Furthermore, let Q : F → | be the quality function we aim to minimize, and Q n : F n → | an estimator of this quality function on a finite sample. With this notation, the true clustering f * on the underlying space and the nearest neighbor clustering f n introduced in the last section are given by

f * ∈ argmin f ∈F Q(f ) and f n ∈ argmin f ∈Fn Q n (f ).
Later on we will also need to work with the functions

f * n ∈ argmin f ∈Fn Q(f ) and f * (x) := f * (NN m (x)).
As distance function between different clusterings f, g we will use

L n (f, g) := P(f(X) = g(X) | X 1 , . . . , X n )
(we need the conditioning in case f or g depend on the data, it has no effect otherwise).

Theorem 1 (Consistency of nearest neighbor clustering) Let (X i ) i∈◆ be a sequence of points drawn i.i.d. according to some probability measure P on | d , and m := m(n) the number of seed points used in nearest neighbor clustering. Let Q : F → | be a clustering quality function, Q n : F n → | its estimator, and A(f ) and A n (f ) some predicates. Assume that:

1. Q n (f ) is a consistent estimator of Q(f ) which converges sufficiently fast: ∀ε > 0, K m (2n) (d+1)m 2 sup f ∈ e Fn P(|Q n (f ) -Q(f )| > ε) → 0. 2. A n (f ) is an estimator of A(f )
which is "consistent" in the following way:

P(A n ( f * ) true) → 1 and P(A(f n ) true) → 1.
3. Q is uniformly continuous with respect to the distance L n between F and F n :

∀ε > 0 ∃δ(ε) > 0 ∀f ∈ F ∀g ∈ F n : L n (f, g) ≤ δ(ε) =⇒ |Q(f ) -Q(g)| ≤ ε. 4. lim n→∞ m(n) = +∞.
Then nearest neighbor clustering as introduced in Section 2 is weakly consistent, that is

Q(f n ) → Q(f * ) in probability.
Proof. (Sketch, for details see von [START_REF] Luxburg | Supplementary material to "Consistent minimization of clustering objective functions[END_REF]). We split the term 2) that the first term converges to 0. The main work consists in bounding the second term. As usual we consider the estimation and approximation errors

P(|Q(f n )-Q(f * )| ≥ ε) into its two sides P(Q(f n ) -Q(f * ) ≤ -ε) and P(Q(f n ) -Q(f * ) ≥ ε). It is a straightforward consequence of Condition (
P Q(f n ) -Q(f * ) ≥ ε ≤ P Q(f n ) -Q(f * n ) ≥ ε/2 + P Q(f *
First we bound the estimation error. In a few lines one can show that

P(Q(f n ) -Q(f * n ) ≥ ε/2) ≤ P(sup f ∈Fn |Q n (f ) -Q(f )| ≥ ε/4).
Note that even though the right hand side resembles the standard quantities often considered in statistical learning theory, it is not straightforward to bound as we do not assume that Q(f ) =

❊Q n (f ). Moreover, note that the function class F n is data dependent as the seed points used in the Voronoi partition are data points. To circumvent this problem, we replace the function class F n by the larger class F n , which is not data dependent. Using symmetrization by a ghost sample (cf. Section 12.3 of [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF], we then move the supremum out of the probability:

P sup f ∈Fn |Q n (f ) -Q(f )| ≥ ε/4 ≤ 2S K ( F n , 2n) sup f ∈ e Fn P |Q n (f ) -Q(f )| ≥ ε/16 inf f ∈ e Fn P |Q n (f ) -Q(f )| ≤ ε/8
(1)

Note that the unusual denominator in Eq. ( 1) emerges in the symmetrization step as we do not assume Q(f ) = ❊Q n (f ). The quantity S K ( F n , 2n) denotes the shattering coefficient, that is the maximum number of ways that 2n points can be partitioned into K sets using the functions in F n .

It is well known (e.g., Section 21.5 of [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF] that the number of Voronoi partitions of n points using m cells in | d is bounded by n (d+1)m 2 , hence the number of nearest neighbor

clusterings into K classes is bounded by S K ( F n , n) ≤ K m n (d+1)m 2 .
Under Condition (1) of the Theorem we now see that for fixed ε and n → ∞ the right hand side of (1) converges to 0. Thus the same holds for the estimation error. To deal with the approximation error, observe that if A n ( f * ) is true, then f * ∈ F n , and by the definition of f * n we have

Q(f * n ) -Q(f * ) ≤ Q( f * ) -Q(f * ) and thus P Q(f * n ) -Q(f * ) ≥ ε ≤ P(A n ( f * ) false) + P f * ∈ F n and Q( f * ) -Q(f * ) ≥ ε . (2)
The first expression on the right hand side converges to 0 by Condition (2) in the theorem. Using Condition (3), we can bound the second expression in terms of the distance L n to obtain

P f * ∈ F n , Q( f * ) -Q(f * ) ≥ ε ≤ P Q( f * ) -Q(f * ) ≥ ε ≤ P L n (f * , f * ) ≥ δ(ε) .
Now we use techniques from [START_REF] Fritz | Distribution-free exponential error bound for nearest neighbor pattern classification[END_REF] to show that if n is large enough, then the distance between a function f ∈ F evaluated at x and the same function evaluated at NN m (x) is small. Namely, for any f ∈ F and any ε > 0 there exists some b(δ(ε)) > 0 which does not depend on n and f such that

P(L n (f, f (NN m (•))) > δ(ε)) ≤ (2/δ(ε)) e -mb(δ(ε)) .
The quantity δ(ε) has been introduced in Condition (3). For every fixed ε, this term converges to 0 due to Condition (4), thus the approximation error vanishes.

Now we want to apply our general theorem to particular objective functions. We start with the normalized cut.

cut n (f k ) := 1 n(n-1) n i,j=1 f k (X i )(1 -f k (X j ))s(X i , X j ) cut(f k ) := ❊ X,Y f k (X)(1 -f k (Y ))s(X, Y ) vol n (f k ) := 1 n(n-1) n i,j=1 f k (X i )s(X i , X j ) vol(f k ) := ❊ X,Y f k (X)s(X, Y ) Ncut n (f ) := K k=1 cutn(f k ) voln(f k ) Ncut(f ) := K k=1 cut(f k ) vol(f k ) Note that ❊ Ncut n (f ) = Ncut(f ), but ❊ cut n (f ) = cut(f ) and ❊ vol n (f ) = vol(f ).
We fix a constant a > 0, a sequence (a n ) n∈◆ with a n ≥ a n+1 and a n → a and define the predicates Proof. We will check that all conditions of Theorem 1 are satisfied. First we establish that

A(f ) is true : ⇐⇒ vol(f k ) > a ∀k = 1, . . . , K A n (f ) is true : ⇐⇒ vol n (f k ) > a n ∀k = 1, . . . , K (3)
{| cut n (f k ) -cut(f k )| ≤ aε} ∩ {| vol n (f k ) -vol(f k )| ≤ aε} ⊂ {| cut n (f k ) vol n (f k ) - cut(f k ) vol(f k ) | ≤ 2ε}.
Applying the McDiarmid inequality to cut n and vol n , respectively, we obtain that for all f ∈ F n

P(| Ncut(f ) -Ncut n (f )| > ε) ≤ 4K exp - na 2 ε 2 8C 2 K 2 .
Together with m 2 log n/(n(aa n ) 2 ) → 0 this shows Condition (1) of Theorem 1. The proof of Condition ( 2) is rather technical, but in the end also follows by applying the McDiarmid inequality to vol n (f ). Condition (3) follows by establishing that for f ∈ F and g ∈ F n we have

| Ncut(f ) -Ncut(g)| ≤ 4CK a L n (f, g).
In fact, Theorem 1 can be applied to a large variety of clustering objective functions. As examples, consider ratio cut, within-sum of squares, and the ratio of between-and within-cluster similarity:

RatioCut n (f ) := K k=1 cutn(f k ) n k RatioCut(f ) := K k=1 cut(f k ) ❊f k (X) WSS n (f ) := 1 n n i=1 K k=1 f k (X i ) X i -c k,n 2 WSS(f ) := ❊ K k=1 f k (X) X -c k 2 BW n := K k=1 cutn(f k ) voln(f k )-cutn(f k ) BW := K k=1 cut(f k ) vol(f k )-cut(f k )
Here n k := i f k (X i )/n is the fraction of points in the k-th cluster, and c k,n :=

i f k (X i )X i /(nn k ) and c k := ❊f k (X)X/❊f k (X)

Implementation using branch and bound

It is an obvious question how nearest neighbor clustering can be implemented in a more efficient way than simply trying all functions in F n . Promising candidates are branch and bound methods. They are guaranteed to achieve an optimal solution, but in most cases are much more efficient than a naive implementation. As an example we introduce a branch and bound algorithm for solving NNC(Ncut) for K = 2 clusters. For background reading see [START_REF] Brusco | Branch-and-Bound Applications in Combinatorial Data Analysis[END_REF]. First of all, observe that minimizing Ncut n over the nearest neighbor function set F n is the same as minimizing Ncut m over all partitions of a contracted data set consisting of m "super-points" Z 1 , . . . , Z m (super-point Z i contains all data points assigned to the i-th seed point), endowed with the "super-similarity" function s(Z s , Z t ) := Xi∈Zs,Xj ∈Zt s(X i , X j ). Hence nearest neighbor clustering on the original data set with n points can be performed by directly optimizing Ncut on the contracted data set consisting of only m super-points. Assume we already determined the labels l 1 , . . . , l i-1 ∈ {±1} of the first i-1 super-points. For those points we introduce the sets A = {Z 1 , . . . , Z i-1 }, A -:= {Z j | j < i, l j = -1}, A + := {Z j | j < i, l j = +1}, for the remaining points the set B = {Z i , . . . , Z m }, and the set V := A ∪ B of all points. By default we label all points in B with -1 and, in recursion level i, decide about moving Z i to cluster +1. Analogously to the notation f k of the previous section, in case K = 2 we can decompose

Ncut(f ) = cut(f +1 ) • (1/ vol(f +1 ) + 1/ vol(f -1
)); we call the first term the "cut term" and the second term the "volume term". As it is standard in branch and bound we have to investigate whether the "branch" of clusterings with the specific fixed labels on A could contain a solution which is better than all the previously considered solutions. We use two criteria for this purpose. The first one is very simple: assigning at least one point in B to +1 can only lead to an improvement if this either decreases the cut term or the volume term of Ncut. Necessary conditions for this are max j≥i s(Z j , A + )s(Z j , A -) ≥ 0 or vol(A + ) ≤ vol(V )/2, respectively. If neither is satisfied, we retract. The second criterion involves a lower bound θ l on the Ncut value of -Set gi = +1.

-If Ncut(g) < Ncut(f ) return g, else return f . 3. // Pruning:

• If vol(A + ) > vol(A ∪ B)/2 and max j≥i (s(j, A + )s(j, A -)) ≤ 0 return f . • Compute lower bound θ l as described in the text.

• If θ l ≥ θu then return f . 4. // If no pruning possible, recursively call bbncut:

• all solutions in the current branch. It compares θ l to an upper bound θ u on the optimal Ncut value, namely to the Ncut value of the best function we have seen so far. If θ l ≥ θ u then no improvement is possible by any clustering in the current branch of the tree, and we retract. To compute θ l , assume we assign a non-empty set B + ⊂ B to label +1 and the remaining set B -= B \ B + to label -1. Using the conventions s(A, B) = Zi∈A,Zj ∈B sij and s(A, ∅) = 0, the cut term is bounded by

Set gi = +1, θ ′ u := min{Ncut(g), θu}, call g ′ := bbncut( S, g, i + 1, θ ′ u ) • Set gi = -1, θ ′′ u := min{Ncut(g ′ ), θ ′ u }, call g ′′ := bbncut( S, g, i + 1, θ ′′ u ) • If Ncut(g ′ ) ≤ Ncut(g ′′ ) then return g ′ , else return g ′′ . }
cut(A + ∪ B + , A -∪ B -) ≥ min j≥i s(Z j , A + ) if A -= ∅ s(A + , A -) + min j≥i s(Z j , A -) otherwise. (4)
The volume term can be maximally decreased in case vol(A + ) < vol(V )/2, when choosing B + such that vol(A

+ ∪ B + ) = vol(A -∪ B -) = vol(V )/2. If vol(A + ) > vol(V )/2
, then an increase of the volume term is unavoidable; this increase is minimal when we move one vertex only to A + :

1 vol(A + ∪B + ) + 1 vol(A -∪B -) ≥ 4/ vol(V ) if vol(A + ) ≤ vol(V )/2 vol(V )/max j≥i (vol(A + ∪ Z j ) vol(A -∪ B \ Z j )) otherw.
(5)

Combining both bounds we can now define the lower bound θ l as the product of Eq. ( 4) and ( 5). The entire algorithm is presented in Fig. 1. On top of the basic algorithm one can apply various heuristics to improve the retraction behavior and thus the average running time of the algorithm. For example, in our experience it is of advantage to sort the super-points by decreasing degree, and from one recursion level to the next one alternate between first visiting branch g i = 1 and g i = -1.

Experiments

The main point about nearest neighbor clustering is its statistical consistency: for large n it reveals an approximately correct clustering. In this section we want to show that it also behaves reasonably on smaller samples. Given an objective function Q n (such as WSS or Ncut) we compare the NNC results to heuristics designed to optimize Q n directly (such as k-means or spectral clustering).

As numeric data sets we used classification benchmark data sets from different repositories (UCI repository, repository by G. Rätsch) and microarray data from [START_REF] Spellman | Comprehensive identification of cell cycle-regulated genes of the yeast saccharomyces cerevisiae by microarray hybridization[END_REF]. Moreover, we use graph data sets of the internet graph and of biological, social, and political networks: COSIN collection, collection by M. Newman, email data by [START_REF] Guimerà | Self-similar community structure in a network of human interactions[END_REF], electrical power network by Watts and Strogatz (1998), and protein interaction networks of [START_REF] Jeong | Centrality and lethality of protein networks[END_REF] and [START_REF] Tsuda | Fast protein classification with multiple networks[END_REF]. Due to space constraints we focus on the case of constructing K = 2 clusters using the objective functions WSS and Ncut. We always set the number m of seed points for NNC to m = log n. In case of WSS, we compare the result of the k-means algorithm to the result of NNC using the WSS objective function and the Euclidean distance to assign data points to seed points. Note that one cannot run K-means on pure network data, which does not provide coordinates. In case of Ncut, we use the Gaussian kernel as similarity function on the numeric data sets. The kernel width σ is set to the mean distance of a data point to its k-th nearest neighbor. We then build the k-nearest neighbor graph (both times using k = ln n). On the network data, we directly use the given graph. For both types of data, we use the commute distance on the graph (e.g., [START_REF] Gutman | Generalized inverse of the Laplacian matrix and some applications[END_REF] as distance function to determine the nearest seed points for NNC.

In the first experiment we compare the values obtained by the different algorithms on the training sets. From the numeric data sets we generated z = 40 training sets by subsampling n/2 points. On each training set, we repeated all algorithms r = 50 times with different random initializations (the seeds in NNC; the centers in K-means; the centers in the K-means post-processing step in spectral clustering). Denoting the quality of an individual run of the algorithm by q, we then report the values mean z (min r q) ± standarddev z (min r q). For the network data sets we ran spectral clustering and NNC on the whole graph. Again we use r = 50 different initializations, and we report min r q. All results can be found in Table 1. For both the numeric data sets (left table, top lines) and the network data sets (right table) we see that the training performance of NNC is comparable to the other algorithms. This is what we had hoped, and we find it remarkable as NNC is in fact a very simple clustering algorithm.

In the second experiment we try to measure the amount of overfitting induced by the different algorithms. For each of the numeric data sets we cluster n/2 points, extend the clustering to the other n/2 points, and then compute the objective function on the test set. For the extensions we proceed in a greedy way: for each test point, we add this test point to the training set and then give it the label +1 or -1 that leads to the smaller quality value on the augmented training set. We also tried several other extensions suggested in the literature, but the results did not differ much. To compute the test error, we then evaluate the quality function on the test set labeled according to the extension. For Ncut, we do this based on the k-nearest neighbor graph on the test set only. Note that this experiment does not make sense on the network data, as there is no default procedure to construct the subgraphs for training and testing. The results on the numeric data sets are reported in Table 1 (left table, bottom lines). We see that NNC performs roughly comparably to the other algorithms. This is not really what we wanted to obtain, our hope was that NNC obtains better test values as it is less prone to overfitting. The most likely explanation is that both K-means and spectral clustering have already reasonably good extension properties. This can be due to the fact that as NNC, both algorithms consider only a certain subclass of all partitions: Voronoi partitions for K-means, and partitions induced by eigenvectors for spectral clustering. See below for more discussion.

Discussion

In this paper we investigate clustering algorithms which minimize quality functions. Our main point is that, as soon as we require statistical consistency, we have to work with "small" function classes F n . If we even choose F n to be polynomial, then all problems due to NP hardness of discrete optimization problems formally disappear as the remaining optimization problems become inherently polynomial. From a practical point of view, the approach of using a restricted function class F n can be seen as a more controlled way of simplifying NP hard optimization problems than the standard approaches of local optimization or relaxation. Carefully choosing the function class F n such that overly complex target functions are excluded, we can guarantee to pick the best out of all remaining target functions. This strategy circumvents the problem that solutions of local optimization or relaxation heuristics can be arbitrarily far away from the optimal solution.

The generic clustering algorithm we studied in this article is nearest neighbor clustering, which produces clusterings that are constant on small local neighborhoods. We have proved that this algorithm is statistically consistent for a large variety of popular clustering objective functions. Thus, as opposed to other clustering algorithms such as the K-means algorithm or spectral clustering, nearest neighbor clustering is guaranteed to converge to a minimizer of the true global optimum on the underlying space. This statement is much stronger than the results already known for K-means or spectral clustering. For K-means it has been proved that the global minimizer of the WSS objective function on the sample converges to a global minimizer on the underlying space (e.g., [START_REF] Pollard | Strong consistency of k-means clustering[END_REF]. However, as the standard K-means algorithm only discovers a local optimum on the discrete sample, this result does not apply to the algorithm used in practice. A related effect happens for spectral clustering, which is a relaxation attempting to minimize Ncut (see von [START_REF] Luxburg | A tutorial on spectral clustering[END_REF] for a tutorial). It has been shown that under certain conditions the solution of the relaxed problem on the finite sample converges to some limit clustering (e.g., von Luxburg et al., to appear). However, it has been conjectured that this limit clustering is not necessarily the optimizer of the Ncut objective function. So for both cases, our consistency results represent an improvement: our algorithm provably converges to the true limit minimizer of K-means or Ncut, respectively. The same result also holds for a large number of alternative objective functions used for clustering.

  Let s : | d ×| d → | + be a similarity function which is upper bounded by a constant C. For a clustering f : | d → {1, . . . , K} denote by f k (x) := ✶ f (x)=k the indicator function of the k-th cluster. Define the empirical and true cut, volume, and normalized cut as follows:

Theorem 2 (

 2 Consistency of NNC(Ncut n )) Let (X i ) i∈◆ be a sequence of points drawn i.i.d. according to some probability measure P on | d and s : | d × | d → | + be a similarity function which is upper bounded by a constant C. Let m := m(n) be the number of seed points used in nearest neighbor clustering, a > 0 an arbitrary constant, and (a n ) n∈◆ a monotonically decreasing sequence with a n → a. Then nearest neighbor clustering using Q := Ncut, Q n := Ncut n , and A and A n as defined in (3) is weakly consistent if m(n) → ∞ and m 2 log n/(n(aa n ) 2 ) → 0.

  Branch and bound algorithm for Ncut: f * = bbncut( S, i, f, θu){ 1. Set g := f ; set A -, A + , and B as described in the text 2. // Deal with special cases:• If i = m and A -= ∅ then return f .• If i = m and A -= ∅:

Figure 1 :

 1 Figure 1: Branch and bound algorithm for NNC(Ncut) for K = 2. The algorithm is initially called with the super-similarity matrix S, i = 2, f = (+1, -1, . . . , -1), and θ u the Ncut value of f .

Table 1 :

 1 Left: Numeric data. Results for K-means algorithm, NNC(WSS) with Euclidean distance; spectral clustering (SC); NNC(Ncut) with commute distance. The top line always shows the results on the training set, the second line the extended results on the test set. Right: Network data. NNC(Ncut) with commute distance and spectral clustering, both trained on the entire graph.

	Numeric data sets breast-c. diabetis	WSS K-means 6.95 ± 0.19 7.12 ± 0.20 6.62 ± 0.22 6.72 ± 0.22	Ncut SC 7.04 ± 0.21 0.11 ± 0.02 0.09 ± 0.02 NNC NNC 7.12 ± 0.22 0.22 ± 0.07 0.21 ± 0.07 6.71 ± 0.22 0.03 ± 0.02 0.03 ± 0.02 6.72 ± 0.22 0.04 ± 0.03 0.05 ± 0.05	Network data ecoli.interact ecoli.metabol helico beta3s AS-19971108	NNC 0.06 0.06 SC 0.03 0.04 0.16 0.16 0.00 0.00 0.02 0.02
	german	18.26 ± 0.27 18.56 ± 0.28 0.02 ± 0.02 0.02 ± 0.02	AS-19980402	0.01 1.00
		18.35 ± 0.30 18.45 ± 0.32 0.04 ± 0.08 0.03 ± 0.03	AS-19980703	0.02 0.02
	heart	10.65 ± 0.46 10.77 ± 0.47 0.18 ± 0.03 0.17 ± 0.02	AS-19981002	0.04 0.04
		10.75 ± 0.46 10.74 ± 0.46 0.28 ± 0.03 0.30 ± 0.07	AS-19990114	0.08 0.05
	splice	68.99 ± 0.24 69.89 ± 0.24 0.36 ± 0.10 0.44 ± 0.16	AS-19990402	0.11 0.10
		69.03 ± 0.24 69.18 ± 0.25 0.58 ± 0.09 0.66 ± 0.18	netscience	0.01 0.01
	bcw	3.97 ± 0.26	3.98 ± 0.26 0.02 ± 0.01 0.02 ± 0.01	polblogs	0.11 0.11
		3.98 ± 0.26	3.98 ± 0.26 0.04 ± 0.01 0.08 ± 0.07	power	0.00 0.00
	ionosph. 25.72 ± 1.63 25.77 ± 1.63 0.06 ± 0.03 0.04 ± 0.01	email	0.27 0.27
		25.76 ± 1.63 25.77 ± 1.63 0.12 ± 0.11 0.14 ± 0.12	yeastProtInt	0.04 0.06
	pima	6.62 ± 0.22	6.73 ± 0.23 0.03 ± 0.03 0.03 ± 0.03	protNW1	0.00 0.00
		6.73 ± 0.23	6.73 ± 0.23 0.05 ± 0.04 0.09 ± 0.13	protNW2	0.08 1.00
	cellcycle	0.78 ± 0.03	0.78 ± 0.03 0.12 ± 0.02 0.10 ± 0.01	protNW3	0.01 0.80
		0.78 ± 0.03	0.78 ± 0.02 0.16 ± 0.02 0.15 ± 0.03	protNW4	0.03 0.76

n ) -Q(f * ) ≥ ε/2 .