Christophe Jaffrelot

India, the United States' New Ally in Asia

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

WTO to denounce Washington's payment of anti-dumping duty proceeds to United States industry and that Canada, co-complainant, had imposed sanctions. 5 -Although Indian Airlines purchased 43 aircraft from Airbus in 2004, the following year Air India announced an order for 50 Boeing 787 Dreamliners in circumstances that baffled its European competitor: While the A380 had just finished its first test flights, the Indian company chose an aircraft for which tests will not begin before 2007. 6 -The Indian-Israeli rapprochement has also followed its course. This is not unconnected to India-U.S. relations. The Indian initiative of reestablishing diplomatic relations with Israel in 1991 was in part aimed at Washington: it was a way of "getting at" the United States. For the coalition parties in power today the Indian-Israeli rapprochement, which culminated in September 2003 with Ariel Sharon's visit to New Delhi, was much more open to criticism than the Indian-U.S. rapprochement. The communists were appalled by it.

A year after Manmohan Signh came to power nothing, or nearly nothing, has changed: military cooperation between the two countries has not slackened 7 and the government has not really readjusted its relations with the Middle East in favor of its old Arab allies. It is significant that the only Indian official to have made the trip to attend the funeral of Yasser Arafat -a close friend of the Nehru/Gandhi family -was Natwar Singh, minister of foreign affairs.

The communists may have secured the suspension of joint military maneuvers between the Indian and Israeli armies, but this is the only concession that has been made to them so far.

American ambivalence and Indian expectations

At first glance, the pursuance of the Indian-United States rapprochement intrigues observers all the more since New Delhi does not seem able to obtain what it demands in exchange. 5 The Hindu, 04/09/2005 6 EADS officials were openly offended at not being given the same treatment as their competitors, in particular for not being allowed to present their product under the same conditions. The Indian authorities argued that since the decision involved renewing their fleet, the airline companies were free to make their own choices. However, given that these are state-owned companies, it is highly likely that political power influenced -or accompanied -the final decision. 7 The third meeting of the Joint Working Group on military cooperation was held in Jerusalem in December 2004. Joint development of a ballistic missile was allegedly on the agenda.

Up to now, the NSSP has not fulfilled its promises. This agreement pertains to four types of issues: closer collaboration in the civilian nuclear and space industries, more trade in high technologies, and an expanded dialogue on missile defense.

A year and a half after the NSSP was signed, only the first two items on this agenda have been broached in any significant manner -but without yielding any noteworthy results.

In the field of nuclear energy, India has not obtained the aid it hoped the United States would grant for the maintenance -and modernization -of its four nuclear energy plants, despite their coming under AIEA rules. First, the United States claims -not without reason -to be bound by the rules of the NPT, which India did not sign; second, it fears that Indian scientists will encourage proliferation. In October 2004, they imposed penalties on two retired Indian engineers suspected of having helped Iran develop its nuclear program.

There seems to have been some thaw in this regard in February 2005 when the US Nuclear Regulatory Commission invited Indian scientists to discuss the risks that natural disasters posed for the country's nuclear installations. Members of India's Atomic Energy Regulatory Board went to the United States for six months in this context -which did not escaped the Pakistani press!8 All the same, India just announced that it planned to work with Russia in building a new nuclear power plant -when a similar type of installation is already under construction in Tamil Nadu.9 This announcement could very well have been made to put pressure on the United States.

In space technology, the United States has only made two concessions to India to date: the headquarters of the Indian Space Research Organisation was struck from the list of entities banned from receiving any technology transfer, and the seven ISRO subsidiaries remaining on this list are now allowed to import materials that can be used for dual purposes figuring on the two lists of sensitive technologies with much flexibility.

However, these gestures seem like a fool's game because, first of all, ISRO headquarters does not import high technology and, second, the two lists of materials now accessible to its subsidiaries contain no sophisticated technology.

Moreover, India has not managed to loosen the ties created between the United States and Pakistan in the aftermath of September 11, 2001, as attest the granting of 3 billion dollars in civil and military aid secured by Islamabad and Washington's complacency in the case of A.Q. Khan, the father of the Pakistani bomb whose nuclear dissemination (even proliferation) activities can scarcely be doubted. 10 The symbol of this "full-scale strategic relationship" -to use the terms of the American ambassador in Pakistan, Ryan C. Crocker, was, however, the announcement of the sale of F-16 fighter jets to Islamabad in March 2005. The Congress Party, which while in the opposition had already objected to Pakistan's being named as one of its "Non NATO allies" by the United States in 2004, immediately objected that such a decision amounted to ignoring the use Pakistan could make of these weapons against India. 11 Manmohan Singh voiced his disappointment to G. Bush in no uncertain terms. 12 Despite this setback, the authorities claim to be as convinced as ever of the excellent relations their country enjoys with the United States.

India is persuaded that the NSSP and other aspects of the cooperation between New Delhi and Washington are tantamount to recognition of India's status as a nuclear power. In the eyes of many Indians, this is a remarkable evolution.

New Delhi also entertains the conviction that the United States would given it access to highly sophisticated military equipment -as had already been the case with the Phalcon system. Now the Americans have offered to sell the Indians F-16s or F-18 "Hornets" when they announced the sale of F-16s to Pakistan. But the Indian army has always been very reluctant to purchase American equipment due to the fickleness of the U.S. Congress, quick to accuse its clients of misconduct and impose sanctions on them -and hence stop delivery of spare parts. As soon as the American offer was made public, former Air Force chief Anil Tipnis and former Army chief Shankar Roychowdhury each separately made known to the government in general and the defense minister in particular their hostility to such a revolution in Indian military procurement. But the civil power may well disregard these objections, all the more so since the Americans may agree to technology transfers that would allow parts of these aircraft to be manufactured in India. 13 unstable areas abroad. This "New Framework for the Indo-U.S. Defense Relationship," signed by Mukherjee and Rumsfeld, heralded a new era in relations between the two countries reflecting "common principles and shared national interest."

Last, New Delhi is counting on United States support to take advantage, when the time comes, of the UN Security Council reform which will allow India to occupy a permanent seat. However, Washington does not seem prepared to back India's candidacy (or Japan's, Brazil's and Germany's for that matter) unless it forgoes the right to veto that the permanent member status allows today. After having been discussed at length internally, this possibility seems excluded by the Indian authorities. The U.S. attitude on this issue has clearly caused deep disappointment on the Indian side.

Beyond the diplomatic and strategic compensations that India expects from the United States, its good will toward the United States can be explained more structurally by the expansion of economic and social ties linking the two countries.

Structural factors of the Indian-United States rapprochement.

An trade partnership with a strong high tech content. The United States is India's largest trading partner today. In both trade and investments, the expansion dates from the economic reform undertaken in 1991. In 10 years, from 1992 to 2002, U.S. exports to India multiplied threefold, 14 whereas American imports of Indian goods multiplied by 3.5. 15 In 2004-2005, the United States remain India's more important trading partner with over 11.1 % of its foreign trade.

The change is even more spectacular in the field of investments. The United States is the second largest foreign investor -after Maurice, through which country many American NRI 16 investments are channeled for tax reasons -with 17.08% of the total in cumulated figures since 1991. The increase in American investments in India is linear and rapid, for the went from 319 million dollars in 2002-2003 to 522 million dollars in 2004. 17 The number of American corporations in India has multiplied by 14 since 1991, numbering now over a thousand. 18 Along with the mass consumer products sold by Mc Donald's, Dominos, Pizza Hut, Pepsi, Coca-Cola, Reebok, Nike and Avon, American investments are also made in 14 From 2.8 billion dollars in 1992 to 6.4 billion dollars in 2002. 15 From 4 to 14 billion dollars. 16 Non Resident Indians constitute the business elite of the Indian diaspora in the United States. 17 The magnitude of less sophisticated outsourcing will further increase the interlocking of Indian-U.S. business activities. Many banks, such as American Express, Citibank and Lehman Brothers and automobile manufacturers (Ford, General Motors) for instance have their accounting and much of their back office activity done in India.

And the number of call centers is constantly on the rise. The largest U.S. company in the field, Convergys Corps, employs over 5,000 people in India although it only set up its first offices there in 2001.

The influence of the diaspora and the student community. In addition to business activities, social ties play an extremely important role in the Indian-U.S. rapprochement. The primary factor in these is the Indian diaspora to the United States, a demographic mass that has doubled in 10 years to 2 million people. The flow of legal Indian immigrants is now greater than that of the Chinese. Many of them come to study computer engineering. If

Oracle established offices in India so early on it is because among its management there are many Indians. Similarly, Adobe decided to set up in India on an initiative of its Indian managers of the company in the United States. In May 2005, the US Citizenship and Immigration Services (USCIS) announced it would grant an additional 20,000 H1B visas 19 to foreign graduates of American universities with a Master's or higher level and specialized in state-of-the-art technology. This measure, aimed primarily to meet the demands of U.S.

corporations in terms of computer engineers, works to the great benefit of Indian students seeking a job in the United States.

Most Indian migrants, however, are students. Their numbers rose from 30,000 in 1996-1997 to 75,000 in 2003-2004. 20 Since 2002-2003, the contingent of Indian students is higher than any other group of foreign students. 21

These students -like their elders -are successful in the United States (when they don't return home). The U.S. year 2000 census in fact shows that Indian-Americans have an average per capita income of $60,093 compared to the national average of $38,885 (only 6 % of them reportedly live below poverty level), which can be explained by the fact that threequarters of this community have attended university. 22 The Indian diaspora is thus in a position to act as a lobby -particularly by financing the "Indian Caucus" in the House of Representatives and the Senate -, but also as a showcase for India, whose image has changed considerably in the United States: the age-old cliché of destitute poverty has been done away with.

19 A quota of 65,000 visas in this category had already been decided earlier in the year. 20 It doubled between 1998-1999 and 2002-2003 -in five years.

21 They already made up 13 % of total, 74,603 out of 586,323, in 2002-2003. 22 To temper the enthusiasm these figures incite, it is worth reading S. Khagram, M. Desai and J. Varughese, "Seen, rich but unheard? The politics of Asian Indians in the United States" in Gordon H. Chang (ed.) Asian Americans and Politics, Washington, Woodrow Wilson Center / Stanford University Press, 2001, pp. 258-284

Obstacles and contradictions

Is the weight of structural factors -as much in economic terms as in social termsenough to consolidate the Indian-U.S. rapprochement observable today in the diplomatic and strategic spheres?

India is probably justified in expecting more constructive gestures from the United States, given New Delhi's growing importance for Washington.

Although Pakistan remains a tactical ally in the "war on terror" and especially in the manhunt for Ben Laden and al-Qaida lieutenants, India is the strategic partner the Americans value today. This country -which Washington constantly touts as "the biggest democracy in the world" -is in fact more stable and, incidentally, better reflects the idea the neo-conservatives have of freedom. India weighs more and more heavily on the globe: the United States is probably counting on New Delhi to help it police the Indian Ocean, even counteract China's growing influence (whether this goes through Myanmar or Gwadar in Pakistan). And Washington can probably count on New Delhi's support on several sensitive issues such as the NMD (National Missile Defense).

Two shadows may, however, cloud the picture in the very near future. At a time precisely when China is manifesting its hostility toward the idea of India having a permanent seat on the UN Security Council -on the pretext that any expansion of this council would jeopardize the necessary reform of the current institution -, the United States may disappoint India's hopes. Indeed India may be deluding itself about American support.

A still more delicate issue has to do with Iran and the energy crunch India is faced with. With an economic growth of 7-8 % annually, at the end of the current decade India will replace South Korea as fourth largest consumer of energy behind the United States, China and Japan. This growth exacerbates the shortage in energy/raw materials the country already suffers from. Repeated power cuts attest to this and are detrimental to the economy.

To alleviate this obstacle, India does not hesitate to do business with countries more or less ostracized from the international community, not only Sudan, 23 Myanmar, country with which India will soon be linked via a gas pipeline, but Iran as well. In January 2005, the stateowned Indian Oil Corp. contracted with the Iranian company Petropars to exploit natural gas resources from Pars, an area which has the largest natural gas reserves in the world known to date. India is determined to purchase natural gas from Iran, even if it means putting the required pipeline through Pakistan. Such an energy supply contract would strengthen the ties between India and Iran, the two countries already having entered into a strategic partnership two years ago. Not only is India helping Iran to build a deepwater port in Chabahar to compete with Gwadar in Pakistan, but military cooperation between India and

Iran is allegedly about to take on a rather sophisticated technological dimension, with the development of submarine batteries that could allow Teheran to forego the Russian technology which is ill-suited to warm waters.

Such a rapprochement bothers Washington, which is doing its best to isolate Teheran so as to dissuade Iran from pursuing its nuclear program. On her visit to India last March, Condoleezza Rice again issued warnings that the American ambassador in New Delhi already proffered in attempt to divert India from this path. She claimed to be aware of the scope of the energy problem facing India, proposing, for one, to seek sources of supply in Turkmenistan -which would route the pipeline through Afghanistan -and also, to study together the possibility of cooperation in the area of civil nuclear power. 24

The Indian authorities reacted vehemently to what appeared to them as an outside pressure bordering on interference. Natwar Singh retorted that they had no problem whatsoever with Iran. The Communist Party of India (Marxist), a key element in the parliamentary coalition that backs Manmohan Singh's government, protested that the United States should not interfere with a decision made by three sovereign countries. 25 The United

States may well nevertheless try to sabotage the pipeline project between Iran and India by intensifying its pressure on Pakistan, an obligatory point of passage for the structure and its weak link: Islamabad cannot disregard the United States' wishes with as much pluck as India does.

The outstanding questions in India-U.S. relations are all of major importance, be it Iran, civil nuclear power, arms contracts or a seat for India at the UN. But for the moment, continuity in policy between the government borne out of the 2004 elections in India and the preceding administration is particularly striking with regard to its relations with the United

States, such that it is not an exaggeration to speak of a deepening.

24 The latter approach meets with strong support among several think tanks, including the Carnegie Endowment for International Peace, as can be seen in the recent report by Ashley J. Tellis that it published, entitled "South Asian Seesaw: A New US Policy on the Subcontinent." 25 It should be noted that the communists are also very hostile to India's purchasing weapons from the United States, not only for ideological reasons but also out of fear of making the country's defense excessively dependent on the United States.

  In the long run, India can also hope to acquire the Patriot anti-missile missile. Preliminary discussions between Indian Air Force officials and Defense Security Cooperation Agency officials in the U.S. took place in April 2005. They continued during defense minister P. Mukherjee's visit to the United States last June 28. The outcome demonstrated spectacular headway: a ten-year cooperation agreement providing for Indian-U.S. collaboration in joint production of military equipment, ballistic defense and intervention for peacekeeping operations in potentially 10 Regarding U.S.-Pakistani relations, see S. Shafqat, "Pakistan and the United States: a future unlike the past ?," Centre for the Study of Pakistan (Columbia University), Occasional Paper Series, March 2005. 11 The Hindu, 29/03/2005 12 Ibid, 26/03/2005 13 Besides, in addition to the F-18s, the United States is considering selling India the E-2C HawkeyeAirborne Early Warning and Control (AEW&C) aircraft, with which India would equip its Gorshkov aircraft carrier, as well as the P-3 Orion maritime reconnaissance plane.

  These figures are all taken from the Economic Survey 2004-2005, New Delhi, Government of India, 2005. 18 Most American companies established in India have grown at a very satisfactory rate. Mc Donald's, which for instance didn't open in Delhi until 1996, had 40 sales points in India seven years later, in 2003, and had sold 45 million burgers in the country! strategic areas of high technology, including in research and development. The computer engineers that India steadily produces and at meager salaries compared to labor costs in the United States, are a specific asset. Driven out of India by a wave of hostility toward U.S. multinationals in the 1970s, IBM is back again; with offices in 12 cities, the company employs 9,000 people and has developed a partnership with three of the best Indian Institutes of Technology, those in Delhi, Kanpur and Chennai, which supply it with internationally-qualified engineers. Hewlett-Packard, since it bought out Compaq in India, is still ahead in terms of labor with 10,000 employees and a network staking its presence in 120 cities. And Microsoft, whose largest foreign investments -other than industrial -have been made in India, has established its software development offices in Hyderabad. In general American computer companies all have set up major R & D offices in India. For instance, 80% of Oracle's 5,000 Indian employees work in its Indian Development Center -the first center of this type outside of the United States for the firm, which set up in Bangalore in 1994. By the same token, Adobe's development center was set up in Noida (near Delhi) in 1998. This is where Acrobat Reader was developed. As for Apple, its famous iPod was designed in the center that the company owns in Hyderabadbefore being mass-produced in factories in Taiwan. The Texas Instruments R & D center in Bangalore -approximately 900 engineers -has already filed for 225 patents, mainly for computer chips. Another industry giant, Intel, is also based in Bangalore. The largest American company R & D center in India, however, is General Electrics, which employs 22,000 people, including 1,800 computer engineers in Bangalore, which has enabled 95 patents to be filed for in the United States. In the Spring of 2005, the Indo-US Science and Technology Forum decided to set up R & D centers in India in various sectors, including design engineering, energy and ecoinformatics. A center devoted to this specialized field has already opened in Bangalore.

Director of the CERI (Sciences-Po/CNRS); recently published, Inde: la démocratie par la caste. Histoire d'une mutationsocio-politique (1885-2005), Fayard,

Dawn 13/02/2005.

The Russians are in fact building a 2,000-megawatt plant in Kudankulam, which should come into service in 2007.


Victorin Martin

Cyril Furtlehner

Yufei Han

Jean-Marc Lasgouttes

GMRF Estimation under Topological and Spectral Constraints

Keywords: Iterative proportional scaling, Gaussian belief propagation, walk-summability, Gaussian Markov Random Field

We investigate the problem of Gaussian Markov random field selection under a non-analytic constraint: the estimated models must be compatible with a fast inference algorithm, namely the Gaussian belief propagation algorithm. To address this question, we introduce the ⋆-IPS framework, based on iterative proportional scaling, which incrementally selects candidate links in a greedy manner. Besides its intrinsic sparsity-inducing ability, this algorithm is flexible enough to incorporate various spectral constraints, like e.g. walk summability, and topological constraints, like short loops avoidance. Experimental tests on various datasets, including traffic data from San Francisco Bay Area, indicate that this approach can deliver, with reasonable computational cost, a broad range of efficient inference models, which are not accessible through penalization with traditional sparsity-inducing norms.

Introduction

The Gaussian belief propagation algorithm [START_REF] Bickson | Gaussian Belief Propagation: Theory and Application[END_REF] (GaBP) is an efficient distributed inference algorithm, well adapted to online inference on large scale Gaussian Markov random fields (GMRF). However, since it may encounter convergence problems, especially with non-sparse structures, it can be of practical interest to construct off-line a GMRF which is compatible with GaBP. When selecting such a model from observations, we potentially face a difficult constrained problem. In the present work, we propose to solve it in an approximate but satisfactory manner, with good accuracy and limited computational cost. To achieve this, we combine various methods, which have been discussed in the context of sparse inverse covariance matrix estimation [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Scheinberg | Learning sparse Gaussian Markov networks using a greedy coordinate ascent approach[END_REF].

The GMRF distribution is naturally characterized by a mean vector µ ∈ R and a positive definite precision (or concentration) matrix A, which is simply the inverse of the covariance matrix C. Zero entries in the precision matrix A indicate conditionally independent pairs of variables. This gives a graphical representation of dependencies: two random variables are conditionally independent if, and only if, there is no direct edge between them. Observations are summarized in an empirical covariance matrix Ĉ ∈ R N ×N of a random vector X = (X i ) i∈{1,...,N } , and we look for a GMRF model with sparse precision matrix A. The model estimation problem can be expressed as the maximization of the log-likelihood: where S bp ++ formally represents the set of positive definite matrices corresponding to some GaBP-compatible GMRF.

A = argmax
Without any constraint on M, the maximum likelihood estimate is trivially A = Ĉ-1 . However, enforcing sparsity with simple thresholding of small magnitude entries may easily ruin the positive definiteness of the estimated precision matrix. In the context of structure learning, where meaningful interactions have to be determined, for instance among genes in genetic networks, the maximization is classically performed on the set of positive definite matrices, after adding to the log-likelihood a continuous penalty function P that imitates the L 0 norm. The Lasso penalty, a convex relaxation of the problem, uses the L 1 norm, measuring the amplitudes of off-diagonal entries in A [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF][START_REF] Hsieh | Sparse inverse covariance matrix estimation using quadratic approximation[END_REF]. Various optimization schemes have been proposed to solve it efficiently [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]. However, the L 1 norm penalty suffers from a modeling bias, due to excessive penalization of truly large magnitudes entries of A. To overcome this issue, concave functions, that perform constant penalization to the large magnitudes, have been proposed. Experimental results indicate promising improvements compared to Lasso penalty by reducing bias, while conserving the sparsity-introducing capability [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF][START_REF] Dicker | Variable selection and estimation with the seamless-L0 penalty[END_REF].

In our context, where compatibility with GaBP has to be imposed, sparsity is a desirable feature, albeit without much guarantee: specific topological properties, like the presence of short loops, are likely to damage the GaBP compatibility, even on a sparse graph. Some spectral properties, e.g. walk-summability [START_REF] Malioutov | Walk-sums and Belief Propagation in Gaussian graphical models[END_REF], which guarantee the compatibility with GaBP based inference, might be relevant too. In order to incorporate these explicitly, we propose an efficient constrained model selection framework called ⋆-IPS, where ⋆ stands for the imposed constraints. Approaches based on the iterative proportional scaling (IPS) procedure [START_REF] Speed | Gaussian Markov distributions over finite graphs[END_REF] have already been discussed for tackling the original sparse inverse covariance matrix problem [START_REF] Jalali | On learning discrete graphical models using greedy methods[END_REF][START_REF] Scheinberg | Learning sparse Gaussian Markov networks using a greedy coordinate ascent approach[END_REF]. A first contribution of this paper is to improve its performance by combining it with block update techniques used in [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF], along with providing some precision guarantee based on duality. Our second and main contribution is to exploit the incremental nature of the method to impose, for a reasonable cost, both topological and/or spectral constraints, to generate GMRF models compatible with GaBP, achieving a very good trade-off between computational cost and precision in inference tasks, as shown experimentally.

The paper is organized as follows. The principles of IPS are described in Section 2. Our method includes a likelihood maximization step at fixed graph structure, for which we give a stopping criterion based on duality. In Section 3, we propose several constraints improving GaBP compatibility of the estimated models, and show how to introduce them in our framework. In Section 4, we describe ⋆-IPS as a whole, discuss its complexity and provide some implementation details. Finally Section 5 is devoted to numerical experiments, both on synthetic data and on real traffic data coming from ≈ 10 3 fixed sensors in the San Francisco Bay Area, to illustrate the use of the method for traffic applications.

IPS-Based GMRF Selection

Iterative proportional scaling has been proposed for contingency table estimation [START_REF] Darroch | Generalized iterative scaling for log-linear models[END_REF] and extended further to MRF maximum likelihood estimation [START_REF] Speed | Gaussian Markov distributions over finite graphs[END_REF]. Assuming the structure of the graph is known, it appears to be less efficient than other gradient based methods [START_REF] Malouf | A comparison of algorithms for maximum entropy parameter estimation[END_REF]. Conversely, local changes based on single row-column updates have been shown to be very efficient, even in the first order setting [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]. In our work, we combine the benefits of the incremental characteristics of IPS to identify links (Section 2.1), with the efficiency of row-column update to optimize their parameters at fixed structure (Section 2.2).

Optimal 1-Link Perturbation

Suppose that we are given a set of single and pairwise empirical marginals pi and pij from a real-valued random vector X = (X i ) i∈{1...N } , and a candidate distribution P (n) , based on the dependency graph G (n) . Let us first describe optimal link addition in terms of likelihood. Let P (n) be the reference distribution, to which we want to add one factor ψ ij to produce the distribution

P (n+1) (x) = P (n) (x) × ψ ij (x i , x j ).
This is a special case of IPS and the optimal perturbation is

ψ ij (x i , x j ) = pij (x i , x j ) p (n) ij (x i , x j ) , (1)
where p (n) ij is the (i, j) pairwise marginal of P (n) . The correction to the loglikelihood can then be written as a Kullback-Leibler divergence:

∆L = D KL (p ij p (n) ij ) = pij (u, v) log pij (u, v) p (n) ij (u, v) dudv.
Sorting all the links w.r.t. this quantity yields the optimal 1-link correction to be made. Hence, the best candidate is the one for which the current model yields the joint marginal p (n) ij that is most divergent from pij . Note that the update mechanism can in fact also be applied if the link is already present.

In the general case, computing the pairwise marginals {p ij , (ij) / ∈ G (n) } is expensive. However, in the GMRF family, these marginals depend only on the covariance matrix associated to P (n) . The correction factor (1) reads in that case

ψ ij (x i , x j ) = exp - 1 2 (x i , x j ) T Ĉ-1 {ij} -C -1 {ij} (x i , x j ) ,
where C {ij} (resp. Ĉ{ij} ) represents the restricted 2 × 2 covariance matrix corresponding to the pair (X i , X j ) of the current model P (n) (resp. of the empirical distribution P) specified by precision matrix A = C -1 (resp. Â = Ĉ-1 ). Let [C {ij} ] denote the N × N matrix formed by completing C {ij} with zeros. The new model obtained after adding or changing link (i, j) reads

A ′ = A + [ Ĉ-1 {ij} ] -[C -1 {ij} ] def = A + [V], (2)
with a log-likelihood variation given by:

∆L = C ii Ĉjj + C jj Ĉii -2C ij Ĉij det(C {ij} ) -2 -log det( Ĉ{ij} ) det(C {ij} ) . ( 3
)
For a 2 × 2 perturbation matrix V = V {ij} , the Sherman-Morrison-Woodbury (SMW) formula allows us to efficiently compute the new covariance matrix as

C ′ = A ′-1 = A -1 -A -1 [C -1 {ij} ] I -[ Ĉ{ij} ][C -1 {ij} ] A -1 . (4)
The number of operations needed to maintain the covariance matrix -and to keep track of all pairwise marginals -after each addition is therefore O(N 2 ). This technical point is determinant to the usefulness of our approach. The identity det(A ′ ) = det(A) × det( Ĉ{ij} )/ det(C {ij} ) ensures that the new precision matrix remains definite positive when both C {ij} and Ĉ{ij} are non-degenerate.

It is also possible to remove links, so that, with help of a penalty coefficient per link, the model can be optimized with a desired connectivity level. For a GMRF with precision matrix A, removing the link (i, j) amounts to setting the entry A ij to zero, and thus ψ ij (x i , x j ) = exp A ij x i x j . The corresponding change of log-likelihood is then

∆L = log 1 -2A ij C ij -A 2 ij det(C {ij} ) + 2A ij Ĉij ,
and, using again the SMW formula, we get the new covariance matrix

C ′ = C - A ij 1 -2A ij C ij -A 2 ij det(C {ij} ) C[B {ij} ]C, with B {ij} def = ñ A ij C jj 1 -A ij C ij 1 -A ij C ij A ij C ii ô .
In this case, the positive-definiteness of A ′ needs to be checked and we have

det(A ′ ) = det(A) 1 -α C ij -C ii C jj 1 -α C ij + C ii C jj ,
so that A ′ is definite positive if the following condition is verified:

1 C ij -C ii C jj < A ij < 1 C ii C jj + C ij .

Block Updates

When a new link is added, existing links become detuned by a slight amount. As pointed out, the optimal update given in Section 2.1 is actually indifferent to whether the considered link exists or not. This means that, after a while, detuned links may be automatically updated if the likelihood gain exceeds the one obtained by adding a new link. We observe in practice that, when many links have been added, all the existing links are slightly detuned, which eventually causes suboptimal or bad decisions for the next links, resulting in a significant departure of the learning curve from the optimal one (see Fig. 2-left in Section 5). However, correcting existing links can become very time consuming, the update of one single link having the same computational cost O(N 2 ) as the addition of one link. There are various options to address this problem. To keep the algorithm fast, robust and simple, we choose to stay with the logic of coordinate descent, by remarking that local updates are still possible via a single rowcolumn update of the precision matrix, as originally proposed in [START_REF] Banerjee | Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data[END_REF] and refined in [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]. In our context, the method is based on the following expression of the log determinant of the precision matrix A:

log det(A) = log det(A \i\i ) + log(A ii -A T i A -1 \i\i A i ),
where A \i\i is the block matrix obtained after taking aside the i th row and column and A i is the i th column vector of A without A ii . The direct optimization of the log-likelihood w.r.t. A i and A ii yields the following updated values:

A ′ ii = 1 Ĉii + A T i A -1 \i\i A i and A ′ i = I V (i) A -1 \i\i I V (i) -1 I V (i) Ĉi Ĉii , (5)
where Ĉi represents the i th column vector of Ĉ, V (i) the set of neighbors of i in the current graph, and I V (i) the identity restricted to entries j ∈ {i} ∪ V (i).

Note that this solution involve the inverse A -1 \i\i of a matrix of size N -1. It is related to C = A -1 as follows:

A -1 \i\i = C \i\i - C i C T i C ii .
The overall cost for updating column (and row) i is thus O(|V (i)| 3 ) for the inversion of I V (i) A \i\i I V (i) and O(N 2 ) to update the covariance matrix C after this change. The log-likelihood gain ∆L reads

-log î Ĉii (A ii -A T i A -1 \i\i A i ) ó -2(A ′ i -A i ) T Ĉi - ï 1 Ĉii + A ′ i T A -1 \i\i A ′ i -A ii ò Ĉii .

Stopping Criterion

If the set of links to be optimized is given by some graph G, the likelihood optimization is a convex problem. Let us investigate its dual properties. Let A denote the precision matrix, and Π a Lagrange matrix multiplier, that imposes the structure given by G. The support of Π is the complementary graph of G: ∀(i, j) ∈ G, Π ij = 0, ∀i, Π ii = 0 and Π is symmetric. Then, given Π, we want to optimize

A Π = argmin M Tr(MΠ) + f (M), with f (M) def = Tr(M Ĉ) -log det(M) being convex for given support G. The explicit solution is A Π = ( Ĉ + Π) -1 . (6)
We assume that Π is such that Ĉ+Π is positive definite, so the dual optimization problem reads Y = argmax Π g(Π), with g(Π)

def = N + log det( Ĉ + Π).
The problem is now concave and, because of the barrier resulting from the log term, we are certain to have a positive definite solution. Thus, for any matrix Π, such that Ĉ + Π is definite positive, g(Π) is a lower bound of the log-likelihood. The support of Π represents the set of links to be removed from the precision matrix. Once optimality is reached for Π, all non-zero entries Π ij correspond to vanishing coefficients A ij in [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF]. We may proceed as before, by computing the potential log-likelihood gain ∆L for such local transformations of the covariance matrix. Local moves in the dual formulation deal with the covariance matrix instead of the precision matrix in the primal one. Let C and C ′ be two covariance matrices differing by a single modification on Π with A = C -1 and

A ′ = C ′-1 . We have det(C ′ ) = det(C) 1 + 2Π ij A ij -Π 2 ij det(A {ij} ) ,
with det(A {ij} ) > 0, since A is definite positive. Maximizing the log-likelihood variation yields the optimal values

Π ij = A ij det(A {ij} )
and ∆L = log

Ç 1 + A 2 ij det(A {ij} ) å .
In practice, we will not use this backward scheme: its computational cost is prohibitive, since the complementary graph, composed of links to be removed, is dense. However, this dual formulation will help us to build a confidence interval.

During the greedy procedure, we always have to maintain C = A -1 but C cannot be used directly to get a dual cost because, except at convergence, it does not fulfill the dual constraints

C ij = Ĉij , ∀(i, j) ∈ G.
Let Π be the correction matrix with coefficients

Π ij def = ( Ĉij -C ij )1 1 {(i,j)∈G} .
Provided that C def = A -1 + Π is definite positive, which happens when A is close enough to the optimum A ⋆ , it satisfies the dual constraints yielding the confidence bound

log det( C) + N ≤ -L(A ⋆ ) ≤ Tr(A Ĉ) -log det(A). We have log det( C) = -log det(A) + log det Ä I + AΠ ä ,
with both A and Π sparse matrices, so the determinant can be estimated in O(N 2 K) operations by expanding the logarithm at order 2 in AΠ . It leads to the following bound

∆L ≤ 1 2
Tr AΠ AΠ , which will be used in practice as a stopping criterion for the link updates.

Introducing Constraints for GaBP Compatibility

Usually, GMRF estimation intends to describe a dependency structure. We pursue here another aim: finding a model suitable for fast inference. While inference in GMRF models can always be performed exactly in O(N 3 ) through matrix inversion, this may not be fast enough for some "real-time" applications on large networks. The GaBP algorithm [START_REF] Bickson | Gaussian Belief Propagation: Theory and Application[END_REF] is a fast alternative to matrix inversion for sparse GMRF, which uses message passing along links in the graph G, and thus, assuming it converges, will perform the inference in O(mKN ), where K is the mean connectivity of G and m the maximum number of iterations before convergence. An important property of the GaBP algorithm is that, whenever it converges, it provides the exact mean values for all variables [START_REF] Weiss | Correctness of belief propagation in Gaussian graphical models of arbitrary topology[END_REF]. Variances are however generally incorrect [START_REF] Malioutov | Walk-sums and Belief Propagation in Gaussian graphical models[END_REF]. Having a sparse GMRF gives no guarantee about its compatibility with GaBP, so we need to impose more precise constraints on the precision matrix and to the graph structure. In this section, we make such constraints explicit and show how to impose them in the framework of Section 2.

Spectral Constraints

The most precise condition known for convergence of GaBP is walk-summability (WS) [START_REF] Malioutov | Walk-sums and Belief Propagation in Gaussian graphical models[END_REF]. Let R(A) def = Adiag(A) contain the off-diagonal terms of A, and let ρ(•) denote the spectral radius of a matrix, that is, the maximal modulus of its eigenvalues. The two equivalent necessary and sufficient conditions for WS that we will use are:

(i) The matrix W(A) def = diag(A) -|R(A)| is definite positive; (ii) ρ(|R ′ (A)|) < 1, with R ′ (A) ij def = R(A)ij √ AiiAjj .
Let us consider a GMRF with WS precision matrix A and investigate under which conditions the model remains WS after a perturbation of a link (i, j). The following proposition gives a sufficient condition: Proposition 1. Let A be a WS precision matrix and denote

W def = W(A). The matrix A ′ = A + [V {ij} ] is WS if Θ(α) > 0, ∀α ∈ [0, 1], ( 7
)
where Θ is the following function

Θ(α) def = det(W -1 {ij} ) α 2 V ii V jj -(|A ij | -|αV ij + A ij |) 2 + α Ä W -1 ii V ii + W -1 jj V jj ä + 2 (|A ij | -|αV ij + A ij |) W ij . (8)
In order to check condition [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF], it is necessary to solve two quadratic equations. Note however that knowledge of the matrix W(A) -1 is also mandatory. We will discuss this point at the end of this section.

Proof. The sufficient condition is obtained as follows: for α ∈ [0, 1] we have

W(A + αV) = W(A) + [φ(αV, A)],
with

φ(V, A) def = ï V ii |A ij | -|V ij + A ij | |A ij | -|V ji + A ji | V jj ò . W(A) being invertible, the determinant of W(A + αV) is expressed as det(W(A + αV)) = det (W(A)) det I + W(A) -1 [φ(αV, A)] ,
and we can check that Θ(α) = det(I + W(A) -1 [φ(αV, A)]), with Θ defined in [START_REF] Furtlehner | Spatial and temporal analysis of traffic states on large scale networks[END_REF]. The spectrum of W(A + αV) being a continuous function of α, as the roots of a polynomial, the condition (7) follows.

Note that the special case of removing one link of the graph always preserves the WS property. Indeed, it will change the matrix A in A ′ such as |R ′ (A ′ )| ≤ |R ′ (A)| where ≤ denotes element-wise comparison. Then, using elementary results on positive matrices [16, p. 22

], ρ(|R ′ (A ′ )|) ≤ ρ(|R ′ (A)|) and thus A ′ is WS whenever A is WS.
As we shall see in the numerical experiments, imposing WS is generally too restrictive. It is easy to find non WS models which are still GaBP compatible. The above principle allows us however to impose a weaker spectral constraint: imposing that the matrix diag(A) -R(A) remains definite positive. This is equivalent to constrain the spectral radius ρ(R ′ (A)) to be strictly lower than 1 and it is a necessary condition for GaBP convergence [START_REF] Malioutov | Walk-sums and Belief Propagation in Gaussian graphical models[END_REF]. We call that condition "weak walk-summability" (WWS) as a relaxation of the WS condition. We obtain the following condition Proposition 2. Let A be a WWS precision matrix, i.e. such as ρ(R ′ (A)) < 1, and S(A)

def = diag(A) -R(A). The matrix A ′ = A + [V ij ] is WWS if Γ (α) > 0, ∀α ∈ [0, 1], (9)
with Γ the following degree 2 polynomial and S def = S(A)

Γ (α) def = α 2 det(VS -1 {i,j} ) + α(V ii S -1 ii + V jj S -1 jj -2V ij S -1 ij ) + 1.
In order to check condition (9), we have to solve a quadratic equation. As in Proposition 1, we need to keep track of an inverse matrix, in this case S(A) -1 .

Proof. Mimicking the proof of Proposition 2, we define, for α ∈ [0, 1],

M(α) def = diag(A + α[V]) -R(A + α[V]) = S(A) + α diag([V]) -R([V]),
and we have

det(M(α)) = det(S(A)) det I + αS -1 [diag(V) -R(V)] = det(S(A))Γ (α).
The spectrum of M(α) being a continuous function of α, condition (9) follows.

Both ( 7) and ( 9) are only sufficient conditions for spectral constraint conservation after a pairwise perturbation. However, there are only a few cases where they lead to rejection of a valid perturbation. Indeed, it means that at least one eigenvalue goes to zero for some α ∈]0, 1[ and is positive again for α = 1.

We have pointed out that checking sufficient condition (7) (resp. ( 9)) imposes to keep track of the inverse matrix W(A) -1 (resp. S(A) -1 ). This will not impact the overall complexity of the algorithm since, using the SMW formula, it can be done in O(N 2 ) operations, like for keeping track of the covariance matrix.

Note that, if we want to maintain these spectral constraints, we will not be able in practice to use the column updates described in Section 2.2. Indeed, computing the optimal perturbation is in this case costly, and we have no easy way to check whether it leads to a new admissible model, since our method would involve higher order characteristic polynomials.

Topological Constraints

We present in this section another approach, based mainly on empirical knowledge about the belief propagation algorithm. Belief propagation has been designed as an exact procedure on trees [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Network of Plausible Inference[END_REF] and short loops are usually believed to cause convergence troubles. In the extreme case, where we forbid the addition of any loop, the best precision matrix estimate based on likelihood is known to be the max-spanning tree w.r.t. mutual information [START_REF] Chow | Approximating discrete probability distributions with dependence trees[END_REF]. Since this is usually not enough, we propose here that the estimated precision matrix contains no loops of size smaller than ℓ. This is quite easy to impose: when adding a link (i, j), we have to search if i is in the neighborhood of j of depth ℓ -1. The computational cost is O(K ℓ ), with K the connectivity of G.

We can impose a more precise condition using the fact that, in the absence of frustrated loops, the GaBP algorithm is always convergent [START_REF] Malioutov | Walk-sums and Belief Propagation in Gaussian graphical models[END_REF]. A frustrated loop is a loop along which the product of partial correlations is negative. Preventing the formation of frustrated loops is very similar to the previous loop constraint; the search cost is the same, the only difference is that we will avoid only this kind of loops. This last constraint cannot be imposed with guarantees during the block updates since the sign of partial correlations along edges may change. Prohibiting frustrated loops would require to store all the loops in the graph, which is by far too costly. However, experimental results show that a change of sign usually corresponds to small partial correlations, which are less likely to cause convergence issues.

Algorithm Description and Complexity

In this section, we give an overview of ⋆-IPS, leaving aside the backtracking option. A formal implementation4 is given in Algorithm 1. Note that we suppose that the initial point of the algorithm corresponds to an empty graph. We may as well start from any precision matrix A, provided that we have computed C = A -1 and W(A) -1 or S(A) -1 if we want to impose spectral constraints.

Algorithm 1 ⋆-constrained Iterative Proportional Scaling (⋆-IPS). The function Check constraint returns true if the perturbation of Û

ı leads to a model compatible with the given constraint ⋆.

Inputs: ⋆: constraint type (∅, WS, WWS, ℓ-LOOP, ℓ-FLOOP); Ĉ: empirical covariance matrix. Parameters: ǫ and ǫu: stopping criteria on the log-likelihood and for the update step; δK: connectivity increment after which an update step is performed.

1: A = C = W -1 = S -1 = diag(1). 2: while ∆Lmax > ǫ and niter < nmax do 3: ∆Lmax ← 0, niter ← niter + 1 4:
for all pairs of nodes Û ı do 5:

compute ∆L ij using (3) 6:

if ∆L ij > ∆Lmax then 7:

if Check constraint(⋆, Û ı, A) then 8:

∆Lmax ← ∆L ij and V ← [V {i,j} ] defined in (2) 9:

A ← A + V, update C using (4). 10:

if ⋆ = WS or WWS then 11: update W -1 or S -1 using the SMW formula. 12:

if connectivity has increased by δK then 13:

while Tr AΠ AΠ > ǫu do 14:

for Û ı | Aij = 0 do 15: ∆Lmax ← 0, compute ∆L ij using (3) 16: if ∆L ij > ∆Lmax then 17: if Check constraint(⋆, Û ı, A) then 18: ∆Lmax ← ∆L ij and V ← [V {i,j} ] defined in (2) 19: if ⋆ = WS

or WWS then 20:

A ← A + V, update C and W -1 or S -1 using the SMW formula. 21: else 22:

Block updates (5) of i and j corresponding to ∆Lmax. Update C using the SMW formula.

Let us clarify the complexity of this algorithm. Each link addition or update has a cost O(N 2 ) to update the covariance matrix. If spectral constraints are imposed, it is necessary to keep track of another inverse matrix, which requires as well O(N 2 ) operations and does not change the complexity. Adding M links will therefore require at least O(M N 2 ) operations. This means that our algorithm complexity is in O(N 3 ) in the sparse regions, whereas it becomes O(N 4 ) in the dense ones. Note that this complexity does not take into account the time spent in link updates. As pointed out in Section 2.2, this update step is useful to avoid departing from the optimal learning curve. For a given bound ǫ u , we observe numerically that the number of updates is O(N ), regardless of the mean connectivity in the sparse regime, so this adds up another O(N 3 ) computational cost. Note that the critical parts of the algorithm, which are the update of the covariance matrix and the search for the perturbation that maximizes likelihood increase, can easily be parallelized.

Let us anticipate on the application to emphasize the usefulness of our algorithm, which complexity is comparable to a direct covariance matrix inversion. Suppose that the workflow is:

1. select off-line a GMRF model based on an empirical covariance matrix; 2. use the above model to perform inference for a "real-time" application (which here means at most a few minutes).

We may allow the first task to take a few hours, and thus matrix inversion is acceptable for quite large networks. However, the resultant model will not be suited to GaBP in the inference task and we will have to resort to exact inference, through a matrix inversion of complexity O((N -n o ) 3 ), where n o is the number of observed variables. Using our sparse GaBP-compatible model instead, with a mean connectivity K, the approximate inference complexity O(mK(N -n o )) will then typically scale down from a few hours to a few seconds or a few minutes depending on the needed precision.

Experimental Results

To have some elements of comparison, let us first quickly describe traditional ways to tackle the maximum likelihood estimation problem with penalty-induced sparsity constraints. It involves a maximization of the form

A = argmax M∈S++ log det(M) -Tr(M Ĉ) + λP (M),
where S ++ is the set of positive definite matrices. A classical penalization function P is a continuous approximation to the discrete L 0 norm like the "seamless L 0 penalty" (SELO) proposed in [START_REF] Dicker | Variable selection and estimation with the seamless-L0 penalty[END_REF]

P (x) = log Å 2|x| + τ |x| + τ ã .
In the following, we set τ = 5.10 -3 , which is empirically good enough. We propose to use the Doubly Augmented Lagrange method [START_REF] Dong | An efficient algorithm for l0 minimization in wavelet frame based image restoration[END_REF] to solve this penalized log-determinant programming.

The second method used for comparison is QUIC [START_REF] Hsieh | Sparse inverse covariance matrix estimation using quadratic approximation[END_REF], which uses the L 1 norm as penalty. This is a second order optimization method, leading to superlinear convergence rates. We perform it directly on the empirical covariance matrix with different values for the regularization coefficient λ. Once the structure has been found, it is necessary to maximize the likelihood. According to our experiments, L 1 norm penalty leads to poor precision matrices in terms of likelihood, even if it may be very efficient to find an existing sparse graph structure.

We can now compare the performance of ⋆-IPS and sparsity penalized likelihood optimization. For generating one single GMRF with a given designed sparsity level, both methods are comparable in terms of computational cost, while ⋆-IPS is faster in very sparse regime. Due to its incremental nature, it also has the advantage of generating a full Pareto set of approximate solutions. To assess the quality of the ⋆-IPS model selection, we first look into data fitting accuracy through log-likelihood, and then investigate its compatibility with GaBP inference. Likelihood and GaBP Compatibility Trade-off. The first test is performed on a randomly generated GMRF of 100 variables. The structure of its precision matrix is an Erdős-Rényi random graph, where each link is assigned a value with random sign and magnitude (between 0.1 and 0.8). A diagonal term is added to make it definite positive. The results of the different algorithms are shown in Fig. 1. Both SELO and IPS algorithm are able to find the true graph with the exact covariance matrix. As expected, WS is a very strict constraint and yields low connectivity models. Using WWS instead yields better GaBP compatibility, but provides no guarantee about it. However, this constraint can be enough to get a GaBP-compatible model with (almost) maximal likelihood (Fig. 1-right). For both models, QUIC is clearly sub-optimal regarding sparsity/likelihood trade-off. Fig. 1-right illustrates ⋆-IPS performance in terms of overfitting. This overfitting starts after the Kullback-Leibler divergence reaches a minimum. This can happen as well as before as after the end of GaBP compatibility. Detecting this point is a classical but difficult statistical problem and further investigations are needed to find the best criterion adapted to our case. The second test (Fig. 2) is performed on traffic data from the San Francisco Bay Area5 . Each sample data is a Ndimensional vector of observed speeds { Vi , i = 1 . . . N }, giving a snapshot of the network at a given time of the day, as measured by a set of fixed sensors.

After filtering out inactive sensors, we finally kept 1020 variables, for which we had data from January to June 2013. The travel time distribution at each link, being bounded with heavy tail, is far from being Gaussian. In order to work with normal variables, we make the following transformation

Y i = Φ -1 1 + Fi (V i ) 2 , ∀i = 1 . . . N, (10)
which maps the speed V i to the positive domain of a standard Gaussian variable Y i , where Fi and Φ are respectively the cumulative distribution functions of the speeds and of a normal distribution. The input of the algorithms we are comparing is the covariance matrix of the vector Y. This mapping is important to use the selected GMRF for the inference tasks in the next section.

Fig. 2-left compares the performance of some ⋆-IPS variants with methods based on penalized norms. IPS with update but no constraint is comparable to SELO optimization, albeit much faster to generate the Pareto set, while QUIC is by far the weakest contender. In fact, QUIC is not performing well because, in this case, there is no true underlying sparse dependency graph. Both IPS and SELO loose the compatibility with GaBP at low likelihood and mean connectivity (respectively at K < 6.5 and K < 4). By contrast, imposing the "no frustrated loops of size 3" constraint (3-FLOOP-IPS) yields a nearly optimal L(K) path, up to some flat regime, which endpoint is still GaBP compatible. This is the best trade-off which can be found among all constraints that we have tested. While the WS constraint is again too restrictive, we notice that WWS yields models which are all the way compatible, but with a suboptimal L(K) path. This is partly due to the absence of block updates, replaced by less efficient local updates. Actually, we see that the WWS L(K) roughly follows the one obtained with ∅-IPS (no update at all), which also delivers only GaBP-compatible models. Our interpretation is that updating the links has the effect of reducing some uncorrelated noise otherwise present in the approximate model. At some point, it may spoil the GaBP compatibility because of stronger correlations being taken into account.

Inverse Models for GaBP Inference. Our original motivation for this work is to provide models for travel time inference for large scale traffic network in real time [START_REF] Furtlehner | Spatial and temporal analysis of traffic states on large scale networks[END_REF]. In this application, from an historical data set, we have to build a GMRF reflecting the mutual information between traffic levels among different segments of the traffic network. Then, in real time, GaBP runs on this GMRF to propagate the information given by observed segments to the other ones.

In our experiment with PeMS data, both the historical and test data sets contain samples of a N -dimensional vector of 5-minutes averaged speeds, obtained from fixed sensors. Given a sample, for which a proportion ρ of the variables is observed, we want to infer the states of the (1 -ρ)N unobserved variables. In practice, we proceed gradually on each test sample, by revealing the variables in a random order, and plot the relative L 1 error |v -v|/v made by the inference model on the unobserved variables as a function of ρ, aggregated over 100 different test samples. The error is measured on the speed, after inference has been done in the space defined by [START_REF] Jalali | On learning discrete graphical models using greedy methods[END_REF].

In principle, ⋆-IPS does not require complete samples and even knowledge of the whole covariance matrix is actually not mandatory. But, for this highway dataset, the samples have no missing values, which allows us to use a brute force k-NN predictor for the sake of comparison. The setting for k-NN is as follows: k samples out of the whole training set are selected according to their mean L 1 distance on the observed variables. Then, for each unobserved variable, the median value is extracted from the k selected samples as a predictor. In the experiments, the value k = 70 has been determined to yield the best k-NN predictions. and k-NN behave very differently. k-NN seems to capture rapidly (ρ ≤ 0.1) the global network behavior, but remains flat after that, despite the additional information. By contrast, GMRF models performance always improves with new information, because of their local nature. Moreover, while constraints applied to IPS offer us more precise models, the role of L as proxy is not completely respected. This is due to overfitting problems: as we can see, the full inverse covariance matrix model behaves worse than the K = 6.5 model obtained by simple IPS. This appears clearly in Fig. 3-left, where the prediction error reaches a minimum before increasing again with L. Finally, Fig. 3-right shows various trade-offs between precision and efficiency of the models. We clearly see that ⋆-IPS extracts more precise and less time-consuming models for traffic states reconstruction. For instance, the highest precision is obtained with WWS-IPS for K ≈ 50, leading to a 10-fold time reduction of the computational time w.r.t. the full inverse covariance model, with a gain of 5% in precision.

Conclusion

In this paper, we revisit IPS and show that it provides an efficient framework to find GMRF models with constraints more specific than sparsity. Comparisons show the merits of the proposed ⋆-IPS in terms of flexibility, likelihood values reached and diversity of solutions, since a Pareto set can be delivered for the computational cost of one estimation.

In terms of trade-off between sparsity and likelihood, ⋆-IPS is comparable to the SELO approach, with less computational cost. In contrast, L 1 based methods do not provide satisfying results in our problem setting.

In addition, the flexibility of ⋆-IPS allows one to embed additional and rather exotic but useful constraints for GaBP compatibility, which is not simple to do with traditional penalized likelihood approaches. Experiments show that both topological and spectral constraints are useful. While the walk-summability constraint seems too strict to be useful in practice, relaxing it to weak walksummability leads to very good models. At the same time, avoiding only frustrated triangles give very satisfactory results in our experiments.

Still in this context, the overfitting problem seems completely open to us. Classical information-theoretic criteria failed in our tests to locate properly where to stop in the incremental link addition process. In fact, we observe that both the link-addition and the link-update procedure can lead to overfitting, and the design of a specific criterion able to avoid this deserves further investigation.

  det(M) -Tr(M Ĉ),

Fig. 1 .

 1 Fig.1. Left: Log-likelihood L as a function of mean connectivity K for ⋆-IPS with different constraints, SELO and QUIC, all computed from the exact covariance matrix C. Right: Kullback-Leibler divergence to the actual distribution as a function of mean connectivity (estimations based on an empirical covariance matrix Ĉ, generated with 1000 samples). The end of GaBP compatibility for each algorithm is indicated by ×'s.

Fig. 2 -Fig. 3 .

 23 Fig. 3. Mean relative L1 reconstruction error as a function of log-likelihood (left) and of computational time cumulated over 100 experiments (right). The different results for k-NN are obtained with historical datasets of sizes 10 3 , 2.10 3 , 5.10 3 , 10 4 and 3.10 4 .

  Fig. 2. Left: log-likelihood L as a function of mean connectivity K. Right: mean relative L1 reconstruction error as a function of the fraction of observed stations on the San Francisco Bay Area network for various methods; results are averaged over 100 sample test experiments and normalized by the score obtained with daytime moving average predictor.

						0.8				IPS
										WS-IPS
										WWS-IPS
										3-FLOOP-IPS
	1600								k-NN Inverse covariance
	Log-likelihood L		L( Ĉ-1 ) IPS WS-IPS WWS-IPS 3-LOOP-IPS 3-FLOOP-IPS SELO (L0) QUIC (L1) IPS with no updates	Mean relative L1 error	0.6				SELO (L0)
	1000
						0.4
	0	50	100	150	200	0	0.2	0.4	0.6	0.8
			Mean connectivity K					Proportion of observed variables

The source code is available at http://www.rocq.inria.fr/~lasgoutt/star-ips/.

Available at the Caltrans PeMS website: http://pems.dot.ca.gov/.


Sami Karkar

Christophe Vergez

Bruno Cochelin

Toward the systematic investigation of periodic solutions in single reed woodwind instruments

Keywords: 43, 75, Ef 43, 75, Pq

Single reed woodwind instruments rely on the basic principle of a linear acoustic resonator -the air column inside the cylindrical or conical bore, coupled with a nonlinear exciter -namely the reed and the air jet entering the mouthpiece. The first one is described by its input impedance, which binds the acoustical pressure and flow at the entry of the bore through a linear relation, whereas the second one has a non-smooth, nonlinear characteristic which combines the pressure on both sides of the reed channel, the jet flow, and the reed motion. To find possible playing frequencies, one often analyses the input impedance spectrum in terms of central frequency, height and width of peaks -a method used in various recent publications on bore geometry optimisation. The exciter influence has rarely been taken into account, and in a few restrictive cases only : for precise, fixed value of control parameters ; through time domain simulations, which cannot give all information on the dynamics ; through simplifications of the equations, allowing analytical calculations of some parts of the bifurcation diagram. A more systematic investigation of a given instrument behaviors depending on control parameters requires the framework of dynamical systems and bifurcation theory, as well as specific numerical tools. In the present work, two continuation methods were used to obtain the bifurcation diagram of a clarinet, as comprehensive as possible. Stable and unstable, periodic and static solution branches are shown, revealing instrument characteristics such as oscillation, saturation, and extinction thresholds, as well as dynamic range.

INTRODUCTION

In this paper, a method for systematic investigation of static and periodic solutions of single reed woodwind instruments is presented. In the first part, the physical model used for this type of musical instruments is described through a set of nonlinear, non-smooth, ordinary differential equations. Then two numerical methods for continuing fixed points and periodic orbits of a parameter-dependent ODE system are described. In a third part, both tools are applied to the physical model of a clarinet revealing its bifurcation diagram. In the last part, the stable periodic solutions are reviewed in detail and characteristics of the instrument such as oscillation threshold, dynamic range, and playing frequency are extracted.

PHYSICAL MODEL OF SINGLE REED WOOD-WIND INSTRUMENTS

In this section, the physical model used in this study is presented. The most simple instrument that fits this model is a clarinet with all its holes closed, thus a cylindrical resonator with a single inward-striking reed at its input.

The versatility of this model makes one able to use a measured input impedance of any instrument for realistic modelling. Transposition to the case of lip-reed brass-like instruments as well as double-reed woodwind instruments will demand very little work, as the models are very similar.

Resonator acoustics

Let us consider the pressure P and the flow U at the input of the resonator. We use the decomposition of the input impedance into complex modes proposed in [START_REF] Silva | Emergence des auto-oscillations dans un instrument de musique à anche simple[END_REF] which leads to the following general form :

Z in (ω) = N m ∑ n=1 C n jω -s n + C * n jω -s * n ( 1
)
where N m is the number of modes taken into account, s n are the poles of impedance and C n their corresponding residues.

Remarks :

• whether using an analytical or a measured input impedance, either poles and residues numerical determination or curve fitting is needed to obtain the values of the coefficients s n and C n ; • I m(s n ) gives the angular frequency of the n th mode, Re(s n ) its damping, and |C n |/Re(s n ) the amplitude of the corresponding peak of the impedance spectrum ; • N m is typically between 10 and 20, because either the cut-off frequency of the tone-hole lattice or the reed low-pass filter behaviour makes high frequency modes useless ;

25-31 August 2010, Sydney and Katoomba, Australia

Proceedings of ISMA 2010

• the input impedance written as above is not equal to zero at null frequency, as it is usually the case, but has a modulus of the same order of magnitude as those of the lowest minima of the impedance spectrum. Even though the resistance of the bore to a constant flow should be derived from fluid dynamics, the impedance induced by a stationary flow is much lower than usual acoustical impedances, and the approximation should be valid.

Reed motion

The motion of the reed tip h(t) is driven by the pressure in the mouth of the player P m (a constant or slowly varying parameter) on one side and the pressure at the input end of the bore P(t) on the other side (see schematics on fig. 1 and2) through a single mode mechanical equation :

1 ω 2 r h ′′ (t) + q r ω r h ′ (t) + (h(t) -h 0 ) = P(t) -P m K (2)
where ω r and q r are the angular frequency and damping parameters of the reed, which may vary depending on the player's control of the embouchure, and K is the equivalent stiffness of the reed.

Introducing the characteristic pressure P c = Kh 0 (which represents the minimum static pressure difference needed to close the reed channel), and the dimensionless variables x = h/h 0 , y = x ′ /ω r , p = P/P c and γ = P m /P c , a dimensionless form of the previous equation is :

y ′ /ω r =1 -x + p -γ -q r y. (3)
Moreover, the reed tip position is limited by a unilateral contact with the mouthpiece, which means x(t) ≥ 0 at all times. As we will see later, it is quite difficult to solve numerically this kind of non-smooth contact law and the adopted regularisation sometimes let x be slightly negative.

Flow through the reed channel

As carefully explained in [START_REF] Hirschberg | Mechanics of Musical Instruments[END_REF][START_REF] Silva | Emergence des auto-oscillations dans un instrument de musique à anche simple[END_REF], the application of Bernoulli theorem between a remote point in the mouth of the player with pressure P m and velocity v m and a point M of the same field line in the reed channel (see fig. 2) with pressure P rc and velocity v rc leads to the following equation :

ρ v 2 m 2 + P m = ρ v 2 rc 2 + P rc . (4)
Assuming a constant velocity profile in the reed channel (this is obviously not true, but will lead to a mean velocity), the flow U is simply v rc multiplied by the reed tip opening h(t) and the reed channel width W . Let us now neglect v m because |v m | << |v rc |. Assuming that the jet formed at the output of the reed channel is totally dissipated by turbulence (without pressure recovery, so that P = P rc ) and using conservation of the flow, one gets to :

U(t) = W h(t) 2∆P/ρ (5)
where ∆P = P m -P(t) is the pressure difference between the mouth of the player and inside the mouthpiece.

In the case where ∆P < 0, the problem is almost symmetrical and the same considerations leads to :

U(t) = -W h(t) -2∆P/ρ (6)
Let us point out that this case has to be included because ∆P might become negative with h > 0.

Using the characteristic flow U 0 = W h 0 2P M /ρ, we define the dimensionless flow u = U/U 0 and get the following unique dimensionless form :

u(t) = H(x(t))sign(γ -p(t))x(t) |γ -p(t)| ( 7
)
where the Heaviside function has been added to deal with any possible negative x value (which, let us recall, might happen when solved numerically because of the regularisation of the non-smooth unilateral contact condition, but is not physical).

Pressure-flow equation

Back to the time domain, the input impedance presented above leads to the following relationship between the flow U, the pressure P and its modal components P n :

P ′ n (t) =C n U(t) + s n P n (t) (for n = 1..N m ) (8)
P(t) =2 N m ∑ n=1 Re P n (t) . (9)
As previously, we define the dimensionless variables p n = P n /P M and get to the following dimensionless equation :

p ′ n (t) =U 0 /P M C n u(t) + s n p n (t) (for n = 1..N m ) (10) p(t) =2 N m ∑ n=1 Re p n (t) . (11)

FRAMEWORK : SOLUTIONS OF A DYNAMICAL SYSTEM AND THEIR CONTINUATION

In this section, we recall a few basic definitions and results concerning dynamical systems, which will be useful when dealing with the numerical methods.

We consider a physical system that is described by an (a set of) ordinary differential equation(s) of the general form :

x ′ = f (x, λ ) (12)
where x(t) is the state vector of the system, x ′ (t) its time derivative, and λ some parameter that might have influence on the trajectories.

Continuation of fixed-points and periodic solutions

A fixed point of such a system is a particular solution of (12) that does not vary in time. Thus it is solution of the algebraic equation :

f (x, λ ) = 0 ( 13
)
Knowing a fixed point (x 0 , λ 0 ), there exists (under some conditions on f ) a continuum of similar solutions (x(λ ), λ ) in a neighbourhood around λ 0 . Investigating the influence of the parameter on the given solution is called the continuation (or A periodic orbit of the considered system is a periodic function of the time x : t → x(t), of shortest period T , that is solution of ( 12) at all times. In case of a periodic orbit, the same kind of results exists and allows the continuation a given periodic solution along a branch when the parameter is varied.

METHODS : NUMERICAL TOOLS FOR CONTIN-UATION

Whereas linear dynamical systems may be investigated analytically, non-linear and moreover non-smooth systems are often difficult to study without numerical tools. Time integration or simulation tools are one of the possible ways of dealing with the difficulty. However, it usually does not uncover all solutions of a given system, especially unstable solutions. Thus, numerical continuation tools are useful to investigate the general behaviour of a system.

Numerical tools are usually designed to solve algebraic equations and continue their solutions. It is thus straightforward to compute branches of fixed-points, whereas for continuing periodic orbits, a discretization is necessary to come down to an algebraic system. Two approaches are then possible :

• Time-domain discretization consists of sampling the solution x(t) over one period into a set of discrete values {x(t 0 ), x(t 1 ), ..., x(t N = t 0 + T ) = x(t 0 )}. The problem is then solved for each time sample using a dedicated method like collocation or differentiation scheme, leading to an algebraic problem. • Frequency-domain discretization is made possible by the Fourier representation of a periodic function. The harmonic balance method then leads to a set of algebraic equations.

Prediction-Correction Method

A classical Prediction-Correction Method (or PCM) consists of two steps. Starting from a given solution x i for the value λ i of the parameter and the quantity ẋi = dx dλ (x i ,λ i )

(which is the tangent of the curve representing the solution branch in a λ -x plane), the prediction step is the computation of a rough approximation of the solution x 0 i+1 = x i + ∆λ ẋi for the value of the parameter λ i+1 = λ i + ∆λ (see figure 3). Then an iterative correction algorithm is used to get a better approximation of the solution (at numerical precision). In the software AUTO [START_REF] Doedel | AUTO-07P : Continuation and Bifurcation Software for Ordinary Differential Equations[END_REF] that was used in this study, a specific parametrisation of the curve called pseudo arc-length is used (see [START_REF] Keller | Numerical Solution of Bifurcation and Nonlinear Eignevalue Problems[END_REF]), allowing to proceed with the continuation even through points where the "tangent", as defined above, tends to infinity (limit points). As this introduces a new variable into the system, namely the pseudo arc-length parameter, a new equation is needed. It is brought by the parametrisation equation, which locally defines this new variable : s = x(s)x i .x t i + λ (s)λ i .λ t i where the superscript t refers to the (newly defined) tangent. In this case, defining the length of the step as ∆s, the tangent predictor is :

x 0 i+1 =x i + ∆sx t i (14) λ 0 i+1 =λ i + ∆sλ t i . (15)
The correction step then uses a Newton-Chord algorithm to converge to the next solution (x i+1 , λ i+1 ).

Asymptotic Numerical Method

The MANLAB software (available online, see [START_REF] Manlab | an interactive continuation software[END_REF]), also used in this study, relies on the Asymptotic Numerical Method (or ANM). Its basic principle is a local high-order Taylor series expansion of all quantities as functions of a path parameter s, which is also the pseudo arc-length parameter in this case. Letting (x i , λ i ) be a known solution of ( 13) for (x(s), λ (s)) to lie on the solution branch for all s in some interval [0, s max ], we write :

x(s) = x i + +∞ ∑ n=1 s n x (n) i (16) λ (s) = λ i + +∞ ∑ n=1 s n λ (n) i . (17)
Then the x

(n) i and λ (n) i
are solutions of a new system of recursive linear equations with invariant left-hand side matrix. Truncating the series to a high order N (typically 15 to 30) gives an accurate, smooth, and continuous description of the branch for values of s below the convergence radius of the series. This way, one can compute a continuous portion of the branch of length s max , which is determined a posteriori as the value of s for which | f (x(s), λ (s))| reaches a user-defined tolerance threshold. A piecewise-continuous representation of a whole branch of solutions is thus obtained by connecting several portions of the branch.

Notes :

• fixing N = 1 is equivalent to the tangent predictor of the previous method with x 1 i = x t i , thus the ANM can be considered as a high-order prediction method ; • a classical Newton-Raphson algorithm has been implemented in case residues would accumulate through steps, but experience shows that usually no correction is needed.

Periodic solutions continuation

Concerning the first method (PCM), the AUTO software uses orthogonal collocation for the continuation of periodic orbits (see [START_REF] Doedel | AUTO-07P : Continuation and Bifurcation Software for Ordinary Differential Equations[END_REF][START_REF] Boor | Collocation at gaussian points[END_REF] for details). The periodic solution is sampled over a period using piecewise polynomial interpolation on each time interval. Writing [START_REF] Silva | Interaction of reed and acoustic resonator in clarinet-like systems[END_REF] for each collocation point leads to a system of algebraic equations. The solution of this system are then continued with the PCM. The number of time intervals for the discretization of the period is adjustable, as well as the number of collocation points per time interval, which defines the order of the interpolation.

In the MANLAB software, the investigation of periodic solutions is carried out combining the Harmonic Balance Method and the ANM. As a periodic solution is described with its truncated Fourier series, the HBM links the coefficients of this series and the angular frequency ω through a system of algebraic equations derived from the primary system (12) using trigonometric identities. The ANM is then used for continuing the solutions of this new system. The reader is refered to [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF] for a detailed description and practical application.

Remark : fixing ω = 0 as the phase equation of the HBM allows to compute static solutions (fixed-points).

RESULTS AND DISCUSSION

In this section, both methods are applied to the physical model of clarinet described in the first section. As we are interested in what happens when the player's blowing pressure varies, γ is chosen as the continuation parameter.

The static solution branch as well as several periodic solution branches, rising from direct Andronov-Hopf bifurcations of the static branch, are obtained. Stability analysis is performed and the only stable periodic regime is reviewed in details.

The values used for the parameters of the model are shown in table 1. The modal coefficients C n and s n are the poles and residues of the analytical input impedance for a cylinder of length L and radius r, determined using the MOREESC software (see [START_REF] Silva | Emergence des auto-oscillations dans un instrument de musique à anche simple[END_REF][START_REF] Moreesc | Modal Resonator-Reed Interaction Simulation Code[END_REF]). Visco-thermal losses are taken into account, but not radiation.

Static regime of the clarinet

In the static case (i.e. fixed-point solutions), the branch can be computed analytically as the solution of a third degree polynomial equation. The solution branch computed by both softwares are perfectly superimposed on each other and on the exact solution (up to numerical precision). However, the use of the AUTO software is preferred here since it performs a stability analysis and detects Andronov-Hopf bifurcations. The result is shown figure 4. The trivial solution where all quanti-Figure 4:

The static regime of the clarinet.

-: stable parts --: unstable parts H : Andronov-Hopf bifurcation.

ties are null for γ = 0 being a singular point (the jacobian of the system is not invertible), it cannot be used as a starting point for continuation. An analytical solution for 0 < γ < 1 has been used instead.

Several bifurcation points appear. All of them are of the Andronov-Hopf type, and a periodic solution branch rises from each one of them. As it will be shown further, the position of the first bifurcation gives the value of γ that is commonly referred to as the oscillation threshold. The computation of these periodic regimes has been carried out. Figure 5 shows the bifurcation diagram of our model : the root-mean-square value of p (computed over one period) is plotted versus γ for each periodic branch computed with AUTO.

Bifurcation diagram

Some branches are quite difficult to compute and a refined mesh of the period is needed to achieve convergence. The fifth branch is problematic, despite a very fine mesh and very small steps. Thus, the bifurcation diagram displayed here is slightly incomplete. However, it must be noticed that only the first periodic regime is stable, and other branches have less physical interest. Another important remark is that the first Andronov-Hopf bifurcation (which is the oscillation threshold) is direct in our case. An extended exploration of the parameter space (especially ω r and q r that can vary through the control of the embouchure by the player) within realistic ranges would be very interesting, to see if an inverse bifurcation is possible at the oscillation threshold with this model (as shown in [START_REF] Silva | Interaction of reed and acoustic resonator in clarinet-like systems[END_REF] and [START_REF] Silva | Emergence des auto-oscillations dans un instrument de musique à anche simple[END_REF]).

Remark : the extinction is a discontinuity-induced, degenerate bifurcation as all periodic branches converge towards the same point in γ = 1, with non-vertical tangents.

Stable periodic regime : details

As it has been pointed out, only the first periodic regime is stable in our case, and thus of great interest. Figure 6 shows this first periodic regime, computed with MANLAB. Root-meansquare as well as maximum value of p over one period (which represents the envelope of the periodic orbits) are plotted versus γ.

The branch differs slightly from the first one computed with AUTO in two regions : near saturation the branch seems to "oscillate" ; near the inverse bifurcation at extinction, the branch is not straight and does not converge to γ = 1. The first phenomenon is mainly caused by the Fourier truncation : in this 6 is computed with 25 harmonics and only exhibits little oscillations. The second phenomenon is due to the regularisation of the unilateral contact force in the mechanical equation. In MANLAB, the force is written in the form : F reg = ε/x 2 with ε << 1. Thus, in the region where x tends to 0 because of a high pressure difference, the regularisation introduces significant relative errors. A smaller regularisation parameter ε makes the branch to vanish closer to γ = 1, but makes the problem stiffer, thus leading to smaller steps and demanding a higher H to obtain convergence.

The lowest value for γ of this branch gives the oscillation threshold : γ osc = 0, 376. By monitoring the harmonics of the x part of the solution, one can deduce the threshold of beating reed : γ br = 0, 498. The highest point on the r.-m.-s. value curve gives the saturation threshold : γ sat = 1, 63, which is related to the dynamic range as it is the loudest possible sound : p max = 1, 237. The rightmost point of the curve gives the extinction threshold (for increasing blowing pressure) : γ ext = 1, 805, beyond which the reed is stuck against the mouthpiece and thus p = u = x = 0.

Remarks :

• all parameters other than γ have fixed values, so it does not represent the general dynamic range of the clarinet, but rather a theoretical one, for this set of parameters ; • experimental data are usually displayed in a way that shows only the envelope of the periodic orbits, which would here correspond to the red curve. However, the loudest sound that can be produced corresponds to the maximum of the r.-m.-s. value of p, which differs in our case from the maximum of its envelope because the signals are neither sinusoidal nor square and their harmonic content varies along the branch. Thus the dynamic range should be deduced from the blue curve ; • the extinction is an inverse bifurcation as there exists two stable regimes (one periodic, one static) for 1 ≤ γ ≤ γ ext , exhibiting a hysteretic behaviour : the extinction threshold for increasing γ and the oscillation threshold for decreasing γ are different. Variation of the fundamental frequency f 0 along the stable part of the periodic regime.

Frequency-domain features

A lot of frequency-domain related quantities are directly accessible along the branch. For instance, we plotted the playing frequency as a function of γ on the figure 7. It shows that the fundamental frequency variation amplitude is up to 18 cents, which is more than noticeable for any normal listener. Thus, the player will be forced to modify other parameters such as the reed's initial opening, or natural frequency and damping, using his lips, to keep the pitch as correct as possible while playing louder. Also, despite the oscillations visible at the end of the branch (causes have been discussed previously), there is a noticeable decrease of the playing frequency when γ is increased from 1, 1 to 1, 7. The playing frequency must be compared with the first resonance frequency of the pipe : I m(s 1 )/2π = 130, 44Hz.

This frequency is actually not representative of a real Bb clarinet : 65cm represents the total length, including the bell, whereas the effective length to consider turns out to be closer to 57cm. However, let us recall that it is only a parameter of the model that has been arbitrarily chosen. As the general behaviour is not sensibly affected, and as the main purpose of this paper is to present a new tool of investigation, all the results presented here are still valid. As for the results concerning a real clarinet, a measured impedance of the instrument would be better, and a precise identification of the reed's parameters would be necessary (which appears not to be simple, according to recent experiments on artificial mouth related in [START_REF] Silva | Emergence des auto-oscillations dans un instrument de musique à anche simple[END_REF]).

Another interesting point is to investigate Worman's rule [START_REF] Elliott | Self-sustained nonlinear oscillations of medium amplitude in clarinet-like systems[END_REF].

Using the representation adopted in previous studies (see [START_REF] Arthur | Fundamentals of Musical Acoustics[END_REF][START_REF] Dalmont | Some Aspects of Tuning and Clean Intonation in Reed Instruments[END_REF], the figure 8 shows the odd harmonics (1, 3, 5, 7, and 9) as functions of the first harmonic in logarithmic scale. Straight lines with respective slopes 1, 3, 5, 7 and 9 are plotted for comparison.

The result is in good agreement with Worman's modified rule (given and demonstrated in Ricaud [START_REF] Ricaud | Behavior of reed woodwind instruments around the oscillation threshold[END_REF])

P n = α n (γ -γ osc ) n/2
where the constant α n is different for each harmonic. However, it is very important to link this figure to the whole branch of solution : the agreement is only good until P 1 = -12dB, which correspond to γ = 0, 379. Going back to figure 6, one can see how narrow is the range of γ, from the oscillation threshold (0, 376) to this value. This results indicate that, despite the argument given by Benade in [START_REF] Arthur | Fundamentals of Musical Acoustics[END_REF], the so-called "change of feel", characterised by a change of slope of the curves P n = f (P 1 ), is not due to the reed beating against the mouthpiece : the beating reed threshold (γ br = 0, 498) lies far beyond the limit of this range.

The advantage of using the HBM with MANLAB is that one have a direct access to the amplitude of the harmonics. It is then very easy to plot the amplitude evolution of each (odd) harmonic with γ as shown on figure 9. It reveals that for the main part of the branch, the relative amplitude of each odd harmonic with respect to the first one is almost constant : P 1 -P 3 ≃ 10dB, P 1 -P 5 ≃ 16dB, P 1 -P 7 ≃ 21dB, and P 1 -P 9 ≃ 25dB.

Time-domain point of view

Reconstructing the time series for variables of interest (the internal pressure p, the flow u, and the reed tip position x) is also possible, demands little post-processing, and can be displayed for any point on the branch. Figure 10 shows such time-domain views for γ = 1, 80, just before extinction. As it can be seen, the reed channel is closed during 70 percents of the period, resulting in a null flow in the mean time. Though, little oscillations of x and u around 0 are visible. It is the result of the Fourier truncation (here H = 25) which is not high enough to render correctly tangent discontinuities.

CONCLUSIONS AND PROSPECTS

As analytical expressions for input impedance of various resonator are often available, allowing to compute easily resonance frequency, analytical work becomes very difficult to apply when one considers the whole dynamical system, especially in the case of highly non-smooth interaction like unilateral contact or dry friction. Then, numerical methods are useful for accessing all characteristics of a given model without simplification or restriction on the parameter values.

In this paper, a physical model for single-reed woodwind instruments has been presented. The static as well as periodic solutions of a clarinet-like instrument based on this model have been investigated. A classical numerical tool for bifurcation analysis, the AUTO software, has been used for the continuation of static branches and bifurcation detection, as well as periodic orbits continuation. A new tool based on the Asymptotic Numerical Method and a high-order harmonic balance formulation, the MANLAB software, has been presented and used for periodic solution continuation.

Lots of characteristics are accessible through the computation of the bifurcation diagram, which reveals the general behaviour of the studied instrument. This tool can also be used as a powerful method for comparison between models, and quantifying the influence of approximations.

In future works, the method will also be applied to the saxophone. Measured input impedance spectrum of real instruments will be used, allowing comparison with artificial mouth experimental data. An adaptation to the (very similar) physical model of brass-like instruments is in progress. Comparison with experiments is necessary and parameter identification would be very interesting. Also, specific path following methods will be applied for the continuation of special points (bifurcations, extrema, ...) with respect to a second parameter.

New developments allowing the investigation of quasi-periodic solutions are also considered.

THANKS

The authors would like to deeply thank Fabrice Silva, Philippe Guillemain, Didier Ferrand and Jean Kergomard for their fruitful discussions, remarks and contributions.

Figure 1 :

 1 Figure 1: Schematics of the clarinet and variables used in the physical model.

Figure 2 :

 2 Figure 2: Schematics of the mouthpiece showing how the flow expression is derived from the Bernoulli theorem.

  August 2010, Sydney and Katoomba, Australia path following) of a fixed-point solution and consists in calculating the branch formed by these continua.

Figure 3 :

 3 Figure 3: Illustration of the tangent predictor.

Figure 5 :

 5 Figure 5: Bifurcation diagram of the clarinet computed with AUTO. R.-m.-s. value of p vs γ. Plain line : stable parts ; dotted lines : unstable parts. The static branch is not visible at this scale.

  Figure 7:Variation of the fundamental frequency f 0 along the stable part of the periodic regime.

9 Figure 8 :

 98 Figure 8: Amplitude (log. scale) of odd harmonics P n as functions of the first one P 1 . Straight lines correspond to slope 3, 5, 7, and 9.

Figure 9 :

 9 Figure 9: Amplitudes (log. scale) of odd harmonics as functions of γ.

Figure 10 :

 10 Figure 10:Reconstructed time series for p, u and x just before exctinction.

Table 1

 1

	:	Model parameters for a clarinet-like in-
		strument
	L	r		N m	W	h 0
	650 mm 7 mm	12	12 mm .3 mm
	ω r /2π	q r		K	ρ
	1500 Hz	1	8.10 6 Pa/m 1, 185 kg/m 3