Sarao Mannelli Stefano

Giulio Biroli

Chiara Cammarota

Florent Krzakala

Pierfrancesco Urbani

Lenka Zdeborová

Stefano Sarao Mannelli

Complex Dynamics in Simple Neural Networks: Understanding Gradient Flow in Phase Retrieval

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

Introduction

In many machine learning applications one optimizes a non-convex loss function; this is often achieved using simple descending algorithms such as gradient descent or its stochastic variations. The positive results obtained in practice are often hard to justify from the theoretical point of view, and this apparent contradiction between non-convex landscapes and good performance of simple algorithms is a recurrent problem in machine learning.

A successful line of research has studied the geometrical properties of the loss landscape, distinguishing between good minima -that lead to good generalization error -and spurious minima -associated with bad generalization error. The results showed that in some regimes, for several problems from matrix completion [START_REF] Ge | Matrix completion has no spurious local minimum[END_REF] to wide neural networks [START_REF] Kawaguchi | Deep learning without poor local minima[END_REF][START_REF] Soudry | No bad local minima: Data independent training error guarantees for multilayer neural networks[END_REF], spurious minima disappear and consequently under weak assumptions [START_REF] Lee | Gradient descent only converges to minimizers[END_REF] gradient descent will converge to good minima. However, these results do not justify numerous other results showing that good and spurious minima are present, but systematically gradient descent works [START_REF] Safran | Spurious local minima are common in two-layer relu neural networks[END_REF][START_REF] Liu | Bad global minima exist and sgd can reach them[END_REF]. In [START_REF] Sarao Mannelli | Passed & spurious: analysing descent algorithms and local minima in spiked matrix-tensor model[END_REF] it was theoretically shown that in a toy model -the spiked matrix-tensor model -it is possible to find good minima with high probability in a regime where exponentially many spurious minima are provably present. In [START_REF] Sarao Mannelli | Who is afraid of big bad minima? analysis of gradient-flow in spiked matrix-tensor models[END_REF] it was shown that this is due to the presence of the so-called threshold states in the landscape, that play a key role in the dynamics of the gradient flow [START_REF] Leticia | Analytical solution of the off-equilibrium dynamics of a long-range spin-glass model[END_REF][START_REF] Antenucci | Glassy nature of the hard phase in inference problems[END_REF]: at first attracting it, and successively triggering the converge towards lower minima under certain conditions [START_REF] Tlh Watkin | Optimal unsupervised learning[END_REF][START_REF] Sarao Mannelli | Marvels and pitfalls of the langevin algorithm in noisy highdimensional inference[END_REF]. However, the spiked matrix-tensor model is an unsupervised learning model and it remained open whether the picture put forward in [START_REF] Sarao Mannelli | Passed & spurious: analysing descent algorithms and local minima in spiked matrix-tensor model[END_REF][START_REF] Sarao Mannelli | Who is afraid of big bad minima? analysis of gradient-flow in spiked matrix-tensor models[END_REF] happens also in learning with neural networks.

In this work we thus study learning with a simple single-layer neural network on data stemming from the well-known phase retrieval problem -that consists of reconstructing a hidden vector having access only to the absolute value of its projection onto random directions. The problem emerges naturally in a variety of imaging applications where the intensity is easier to access than the phase [START_REF] Walther | The question of phase retrieval in optics[END_REF][START_REF] Millane | Phase retrieval in crystallography and optics[END_REF][START_REF] Robert W Harrison | Phase problem in crystallography[END_REF][START_REF] Bunk | Diffractive imaging for periodic samples: retrieving one-dimensional concentration profiles across microfluidic channels[END_REF] but it appears also in acoustics [START_REF] Balan | On signal reconstruction without phase[END_REF] and quantum mechanics [START_REF] Corbett | The pauli problem, state reconstruction and quantum-real numbers[END_REF]. The phase retrieval problem considered here leads to a high-dimensional and non-convex optimization problem defined as follows.

Phase retrieval. Consider αN N -dimensional sensing vectors x x x m with unitary norm and generated according to a centered Gaussian distribution, and the true labels y 2 m with y m = x x x m , W W W * and W W W * an unknown teacher-vector (the signal in the phase retrieval literature) from S N -1 ( √ N ). Given the dataset, the goal is to build an estimator W W W of W W W * by minimizing the loss function

L(W W W ) = 1 2 αN m=1 x x x m , W W W 2 -x x x m , W W W * 2 2 = 1 2 αN m=1 (ŷ m , y m ), (1)
with ŷm = x x x m , W W W and (ŷ, y) = (ŷ 2 -y 2 ) 2 is a modified square loss commonly used in the literature [START_REF] Emmanuel | Phase retrieval via wirtinger flow: Theory and algorithms[END_REF][START_REF] Chen | Gradient descent with random initialization: Fast global convergence for nonconvex phase retrieval[END_REF] that ensures a smoother landscape compared to the square loss with the absolute values. The loss function, Eq. ( 1), is minimized using gradient-descent flow on the sphere starting from random initialization

Ẇ W W (t) = -∇ W W W L(t) + µ(t)W W W (t), W i (t = 0) ∼ N (0, 1) ∀ i, (2)
with µ(t) the Lagrange multiplier that enforces the spherical constraint during the dynamics. The value of µ(t) can be readily evaluated by taking the scalar product of the gradient of the loss with W W W and dividing by N

µ(t) = 1 2N αN m=1 ∂ ∂ ŷm (ŷ m (t), y m )ŷ m (t). (3)
We can finally define the Hessian of the problem, denoting δ i,j the Kronecker delta, it reads

H i,j (W W W ) = 1 2 αN m=1 ∂ 2 ∂ ŷ2 m (ŷ m , y m )x m,i x m,j -µδ i,j . (4)
Related work and our main contributions.

Numerous algorithms can be applied to achieve a good reconstruction of the hidden signal (teacher vector) [START_REF] Emmanuel | Phase retrieval via wirtinger flow: Theory and algorithms[END_REF][START_REF] Schniter | Compressive phase retrieval via generalized approximate message passing[END_REF]. For random i.i.d. data and labels generated by a teacher the information-theoretically optimal generalization error has been analyzed in [START_REF] Barbier | Optimal errors and phase transitions in high-dimensional generalized linear models[END_REF], showing that for α < α WR = 0.5 in the limit of N → ∞ no estimator is able to obtain a generalization error better than a random guess. On the other hand for α > α IT = 1 algorithms (not necessarily polynomial ones) exist that are able to achieve zero generalization error. While the weak-recovery threshold is achievable with efficient algorithms [START_REF] Mondelli | Fundamental limits of weak recovery with applications to phase retrieval[END_REF], the information-theoretic threshold is conjectured not to be and an algorithmic gap has been conjectured to exist, with perfect recovery achievable with the approximate message passing algorithm only for α > α alg = 1.13 [START_REF] Barbier | Optimal errors and phase transitions in high-dimensional generalized linear models[END_REF].

In the present paper we are interested in the algorithmic performance of the gradient descent algorithm. The motivation of the present work is not to compare vanilla gradient descent to other algorithms but rather use this well studied phase-retrieval problem to get insights on the properties of the gradient descent algorithm in high-dimensional non-convex optimization. The spurious minima are known to disappear in phase retrieval if the number of samples in the dataset scales as O(N log 3 N ) with N the input dimension [START_REF] Sun | A geometric analysis of phase retrieval[END_REF]. Later [START_REF] Chen | Gradient descent with random initialization: Fast global convergence for nonconvex phase retrieval[END_REF] showed that randomly initialized plain gradient descent will solve the phase retrieval problem with O(N poly(log N ))

samples. An open theoretical question is whether randomly initialized gradient descent is able to solve the phase retrieval problem with O(N ) samples. We explore this question in the present work.

Our main contributions can be summarized in the following points: The left panel has labels created using a teacher, while the right panel has random labels constructed using a Gaussian distribution with variance matching the teacher-labels. The left figure shows that also the simulations with a teacher approach the threshold energy before transiting to the global minimum, if α is large enough.

• We show empirically that in the phase retrieval problem, randomly initialized gradient flow is attracted by the so-called threshold states. We show that as the number of samples increases, the geometry of the threshold states changes from minima to saddles with the slope pointing toward the teacher-vector. This transition is akin to the BBP phase transition known in random matrix theory [START_REF] Baik | Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices[END_REF]. This transition affects gradient descent that, following the slope, achieves zero generalization error.

• We obtain a close formula for the number of samples per dimension α c at which this transition happens. This depends on the joint probability distribution of the labels of the student and the teacher at the threshold.

• Using the replica theory from statistical physics as a non-rigorous proxy we characterize the approximate distribution of the labels of the student and the teacher leading to a approximate prediction for threshold α c = 13.8, notably suggesting that the additional logarithmic factors in the previous works might not be needed.

BBP on the threshold states

In Fig. 1 (left) we show the loss as a function of iteration time for the phase retrieval problem, defined above, with varying number of samples to dimension ratio α in many different runs (full lines). In the right hand side of the figure we show the loss but this time for labels that do not come from a teacher network, but that are randomly reshuffled. We see that in that case the loss converges after a very long time towards a value marked by the dotted line (reproduced also in the left part), that we defined to be the so-called threshold energy. We see that for small α, e.g. α = 6 the train loss on the phase retrieval problem does not decrease to zero (nor the test one), while for the larger values α = 18 and 24 it does very rapidly. The value α = 12 is close to a critical regime where some realization find perfect generalization and other do not, with the dynamics staying for a long time close to the threshold energy.

In a different model, the spiked matrix-tensor, Ref. [START_REF] Sarao Mannelli | Who is afraid of big bad minima? analysis of gradient-flow in spiked matrix-tensor models[END_REF] described exactly this phenomenology and showed that gradient flow starting from random initial conditions has a transition when the Hessian of the spurious minima that trap the dynamics, the so-called threshold states, display a BBP transition [START_REF] Baik | Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices[END_REF]. This leads to the emergence of a descending direction toward the informative minimum which is correlated with the ground truth signal W W W * . In Fig. 2 we argue that the same mechanism is at play in the phase retrieval problem and based on these insights we derive an analytic equation for the corresponding threshold in Sec. 2.1. On the left figure, we show the evolution of the density of the bulk of the eigenvalues, from zero density in white to high density in blue, the smallest eigenvalue in orange, and the second smallest in green. The picture shows that a BBP transition occurs when the training loss approaches the threshold energy. The right panel depicts the evolution of the training loss (purple) and the generalization loss (dark blue) in time. The training loss rapidly approaches a plateau at the level of the threshold states (black dashed line) and converges towards the teacher after the BBP transition. Namely, when the smallest eigenvalue detaches from the rest of the bulk. The bottom panels show the distribution of the eigenvalues at five different instants in the dynamics. At iteration ≈ 1210 we cross the threshold energy and we observe an isolated eigenvalue detaching from the bulk.

Theory for the BBP threshold

Based on the numerical results just presented, we aim to obtain an equation determining the value of threshold α c such that for α > α c the BBP transition occurs, whereas for α < α c it does not. For α < α c the system is at long times trapped in the threshold states and not able to recover, even weakly, the signal. We define P (ŷ, y) the long-time limit of the distribution of the estimated labels and the true labels; P (ŷ, y) allows us to study the Hessian of the threshold states, which is the random matrix defined in (4) with ŷµ and y µ distributed following the law P (ŷ, y).

A type of random matrix M i,j with similar structure as the contribution to the Hessian coming from the Loss function, H i,j = -M i,j -µδ i,j , has been studied recently in [START_REF] Yue | Phase transitions of spectral initialization for high-dimensional nonconvex estimation[END_REF] and the convergence in probability for large N of the largest and the second largest eigenvalues of such matrix was proven.

We applied the results of the Theorem 1 of [START_REF] Yue | Phase transitions of spectral initialization for high-dimensional nonconvex estimation[END_REF] to determine the behavior of the smallest and second smallest eigenvalue of the Hessian. Call

Ψ α (λ) = λ 1 α -E ŷ,y α∂ 2 ŷ (ŷ, y) 2λ + α∂ 2 ŷ (ŷ, y) , Φ(λ) = -λE ŷ,y α∂ 2 ŷ (ŷ, y)y 2 2λ + α∂ 2 ŷ (ŷ, y) ; (5)
let λ = arg min Ψ α (λ) and ξ α (λ) = Ψ α (max{λ, λ}) then the largest eigenvalue of M i,j is λ 1 → P ξ α (λ * α ) with λ * α being the solution of ξ α (λ) = Φ(λ). The second largest is λ 2 → P Ψ α ( λα ). A BBP transition occurs at α BBP , the largest α such that λ 1 = λ 2 . Following [START_REF] Yue | Phase transitions of spectral initialization for high-dimensional nonconvex estimation[END_REF] this leads to an equation on λ and α BBP which reads:

1 α BBP = E ŷ,y α BBP ∂ 2 ŷ (ŷ, y) 2 λ + α BBP ∂ 2 ŷ (ŷ, y) 2 . ( 6
)
We now use an additional assumption, which comes from studies of gradient descent dynamics of mean-field spin-glasses [START_REF] Cugliandolo | Dynamics of glassy systems[END_REF], that the threshold states are marginal, i.e. the smallest eigenvalue of their Hessian is null. As the smallest eigenvalue of the Hessian is determined by the largest eigenvalue of M i,j , this imposes the additional condition λ 1 = λ 2 = -µ. Using the second Eq. in [START_REF] Safran | Spurious local minima are common in two-layer relu neural networks[END_REF] to enforce this last condition and the definition of µ in Eq. ( 3) one finds (see SM for more details):

µ = λE ŷ,y α BBP ∂ 2 ŷ (ŷ, y)y 2 2 λ + α BBP ∂ 2 ŷ (ŷ, y) , µ = α BBP 2 E ŷ,y ∂ ∂ ŷ (ŷ, y)ŷ . (7)
Together the three eqs. ( 6), and the two in [START_REF] Sarao Mannelli | Passed & spurious: analysing descent algorithms and local minima in spiked matrix-tensor model[END_REF] allow to determine the three unknown α BBP , λ, µ.

Assuming that α c = α BBP , they therefore provide an equation for the algorithmic threshold once P (ŷ, y) at threshold states is known.

Further numerical justifications

We now illustrate and test the theory by numerical simulations. Fig. 2 shows that the training loss slowly tends to the threshold energy, before departing in direction of the global minimum. According to the theory described in the previous section, the phenomenon occurring is a BBP transition. In order to confirm this we characterize the spectrum of the Hessian and focus on the smallest eigenvalues. On left figure of Fig. 2 we show the density of eigenvalues in blue and we show separately the smallest eigenvalue in orange and the second smallest eigenvalue in green. During the dynamics the spectrum moves compactly forming a bulk -except for the largest eigenvalues that do not play role in the transition. Approaching the threshold states, for α > α c the dynamics feels the presence of a descending direction associated with the smallest eigenvalue. This becomes increasingly strong and when the negative eigenvalue pops out of the bulk, the dynamics will follow and will converge to the global minimum. In the lower panel of Fig. 2 we show 5 snapshots of the spectral density during the dynamics.

In Fig. 3 we study the spectral properties of the Hessian at the time when the training loss hits the threshold energy. In the first panel from the left, we show the difference between the second smallest eigenvalue and the smallest eigenvalue (solid line) and the difference between the third smallest and the second smallest eigenvalue (dotted). The two quantities are close, and very small, until a certain value value of α -that depends on the size and is reported in the caption -after which the former increases linearly whereas the latter remains small. The second property that we investigate, central figure, is the overlap of the first and second eigenvectors with the teacher-vector W W W * (respectively solid and dotted lines). This overlap is zero before the transition, then as λ 2 -λ 1 , it suddenly increases and finally saturates. The first two panels are provide a strong evidences that a BBP transition is taking place. To further corroborate this findings, in the right panel we consider the fluctuation of the Figure 4: On the left panel is the loss of the threshold states from the simulations and from the analysis with the 1RSB replica method, evaluated for different values of α. The errorbars given by mean and standard deviation of 1000 simulations with input dimension N = 2048 and shuffled labels.

As expected increasing the number of datapoints the number of violated constraints increases and so does the threshold energies. On the right we compare the moments of the label distribution Eq. ( 17) (solid lines) with the moments obtained in the simulations (circles).

overlap (Fig. 3 right panel). In statistical physics terms, the overlap plays the role of order parameter, and its fluctuations diverge at the BBP transition. We indeed find that at the value of α at which the BBP transition seems to take place fluctuations peak (the more so the larger is N as expected for a phase transition).

Finally, we also note the presence of strong finite size effects that shift the effective value of α at which the transition (cross-over for finite N ) takes place.

Characterization of threshold states

The pivotal point of our analysis is that the gradient-flow dynamics is that the spurious minima trapping the dynamics and hindering weak recovery are the threshold states. The study of threshold states has been started in statistical physics of disordered systems, where it has been conjectured and verified that they play a prominent role in the gradient descent dynamics of such random systems [START_REF] Leticia | Analytical solution of the off-equilibrium dynamics of a long-range spin-glass model[END_REF]. according to this theory, threshold states are the most numerous minima and the Hessian associated to the threshold states has a spectrum that is positive semi-definite and gap-less.

In the previous section we have obtained a closed set of equations that allow us to obtain α BBP once the distribution P(ŷ, y) is known. Such distribution can be derived using tools from statistical physics, in particular replica theory [START_REF] Antenucci | Glassy nature of the hard phase in inference problems[END_REF][START_REF] Monasson | Structural glass transition and the entropy of the metastable states[END_REF][START_REF] Franz | Universality of the sat-unsat (jamming) threshold in non-convex continuous constraint satisfaction problems[END_REF]. The main quantity of interest is the partition function, which is defined as the normalization constant of the Gibbs distribution associated with the loss Eq. ( 1),

Z = S α(N -1) ( √ N ) e -βL(W W W ) dW W W , ( 8
)
where β is a parameter associated to the inverse temperature in the physical analog of the problem. We consider the β → ∞ limit to study the minima of L that are relevant in gradient-flow dynamics.

The disordered systems approach focuses on the average properties of the systems. In order to ensure concentration properties as the input dimension goes to infinity, the quantity to average is the logarithm of the partition function. This quantity in general is hard to compute and we resort to replica theory. The idea is to apply replica method, described below and in more detail in the SM, to move from the moments of the partition function to average of its logarithmic value. The analysis leads to the free energy density, Eq. ( 11), that we use to compute the distribution of the labels.

We start writing the moments of Eq. ( 8)

Z n = S α(N -1) ( √ N ) E (x x x,y) e -βnL(W W W ) dW W W , (9)
where we represent the average over the possible datasets with the overbar. The high-dimensional integral in Eq. ( 9) can be written in term of the overlap matrix Q a,b = W W W a • W W W b /N that encodes the similarity between two configurations extracted with the Gibbs measure, of the labels ŷµ and y µ ,

Z n = e N nS(Q Q Q) n a≥b=0 dQ a,b with S(Q Q Q) = 1 2 log det Q Q Q + α log exp   1 2 n a,b=0 Q a,b ∂ 2 ∂h a ∂h b   exp - β 2 a (h a , h 0 ) {ha=0 ∀a} . ( 10
)
We perform a so-called 1-step replica symmetry breaking (1RSB) analysis [START_REF] Mézard | Spin glass theory and beyond[END_REF] that allows us to obtain an explicit expression for the distribution of (ŷ, y). Using the 1RSB ansatz on this problem we reduce the number of parameters in Eq. ( 10) to χ and z: where χ is related to amplitude of the fluctuations in minimum when this is perturbed, and z is a Lagrange multiplier that we use to enforce that the minima in consideration have a Hessian with gapless spectrum, i.e. they are threshold states. The general prescription for finding the global minimum of this kind of problem would require to optimize over χ and z, however, the same formalism can be used to characterize the threshold states by imposing a value for z and optimizing only over χ [START_REF] Monasson | Structural glass transition and the entropy of the metastable states[END_REF][START_REF] Zamponi | Mean field theory of spin glasses[END_REF]. We evaluate the integral of the partition function Eq. ( 8) at the threshold using Laplace's method on Eq. ( 10). Finally we consider the replica method, by taking the analytic continuation of the moments of the partition function and the limit n → 0 + . With this considerations we obtain the free energy density in low temperature limit

f (χ|α, z) = -lim β→∞ lim N →∞ 1 βN log Z/(βN ) that reads f (χ|α, z) = - 1 2z log χ + z χ - α z R dy 1 √ 2π e -y 2 2 log γ 1 e -z 2 V (ŷ,y|χ) ŷ=0 , (11)
where V (ŷ, y|χ) . = min ỹ (ŷ-ỹ) 2 χ + (ỹ, y), the symbol represent a convolution and we write γ 1 as a compact notation to represent a centered Gaussian distribution with unit variance. The arguments χ and z are given via implicit equations that respectively impose that the Laplace's approximation on f and that the free energy describes the threshold states (further details are in the SM):

1 χ(χ + z) = α 4 ∞ -∞ 1 √ 2π e -y 2 2 ∞ -∞ 1 √ 2π e - ŷ2 2 e -z 2 V (ŷ,y|χ) (∂ ŷ V (ŷ, y|χ)) 2 dŷ ∞ -∞ 1 √ 2π e -ŷ2 2 e -z 2 V (ŷ,y|χ) dŷ dy; (12)
1 = α 4 χ 2 ∞ -∞ 1 √ 2π e -y 2 2 ∞ -∞ 1 √ 2π e - ŷ2 2 e -z 2 V (ŷ,y|χ) ∂ 2 ŷ V (ŷ, y|χ) 2 dŷ ∞ -∞ 1 √ 2π e -ŷ2 2 e -z 2 V (ŷ,y|χ) dŷ dy . (13)
In order to obtain P(ŷ, y) we follow the strategy introduced in [START_REF] Franz | Universality of the sat-unsat (jamming) threshold in non-convex continuous constraint satisfaction problems[END_REF]. The partition function in Eq. ( 8) can be written as a functional integral over the empirical labels ρ quen (ŷ, y) = 1 αN αN m=1 δ(yy m )δ(ŷ -ŷm ):

Z = D[ρ quen (ŷ, y)] e -βαN 2 dŷdyρquen(ŷ,y) (ŷ,y)+N S[ρquen(ŷ,y)] , (14)
where S[ρ quen (ŷ, y)] is the entropic factor (evaluated on configurations corresponding to the threshold states). This makes clear that the free-energy can also be obtained in terms of a large-deviation principle:

f = Min ρ α 2 dŷdyρ quen (ŷ, y) (ŷ, y) -S[ρ quen (ŷ, y)] (15)
The minimizer of this variational problem corresponds by definition to P(ŷ, y). By taking the functional derivative of f with respect to (h, h 0 ) one obtains:

δf δ (ŷ, y) = α 2 P(ŷ, y). (16)
Since we have an explicit expression of f in terms of , see eq. ( 11), we can obtain this distribution directly. In taking the functional derivative of Eq. ( 11) we must be careful in considering the implicit dependence of V (ŷ, y|χ) from (ŷ, y). Finally inverting Eq. 9160 we obtain the label distribution

P 1RSB (ŷ, y) = 1 √ 2π e -y 2 2 exp - ŷ2 2 -z 2 V (ŷ, y|χ) R exp - ỹ2 2 -z 2 V (ỹ, y|χ) dỹ . ( 17
)
The index 1RSB denotes that the free-energy has been obtained within the 1RSB scheme.

We run numerical experiments to characterize the threshold states in terms of their energy and their moments with respect to the label probability distribution. In Fig. 4 (left) the 1RSB threshold energy is plotted together with the value of the plateau obtained from the simulations in Fig. 1. The 1RSB ansatz appears to be a good approximation of the empirically obtained energy. In the inset we highlight the discrepancy with the empirical line, which may be due to both finite size effects as well as to the fact that the 1RSB scheme is an approximation and more involved schemes must be employed to obtain the exact distribution [START_REF] Mézard | Spin glass theory and beyond[END_REF]. Finally in Fig. 4 (right) we show the first 4 moments of the label distribution from the 1RSB analysis and we compare them with the numerical results. The 1RSB moments give a nice agreement in relation to the empirical ones.

Encouraged by the reasonable, though not perfect, accuracy of the 1RSB approximation, we use the expression [START_REF] Balan | On signal reconstruction without phase[END_REF] as input for the three eqs. in [START_REF] Liu | Bad global minima exist and sgd can reach them[END_REF][START_REF] Sarao Mannelli | Passed & spurious: analysing descent algorithms and local minima in spiked matrix-tensor model[END_REF]. Their solutions leads to a finite threshold at α 1RSB BBP ≈ 13.8. This value could be compatible with the numerical results presented in Fig. 3 if the finite size shift saturates for yet larger sizes.

Discussion

In Fig. 5 we present the fraction of runs of gradient descent that converge to zero test error in the phase retrieval problem under consideration. We see that a large fraction of simulations achieves convergence before the α 1RSB BBP threshold, in general before α ≈ 11. The mechanism behind this difference is not clear to us, as well as it is not clear whether it will hold as the input dimension tends to infinity. The numerical results are affected by strong finite size effects and may be both consistent with a transition taking place at finite alpha in the infinite dimensional limit and with curves that shift in alpha as a logarithm in the input dimension. Empirically we do not observe any unsuccess above α 1RSB BBP which is consistent with our theoretical scenario for any N > 64. In [START_REF] Sarao Mannelli | Passed & spurious: analysing descent algorithms and local minima in spiked matrix-tensor model[END_REF] the authors showed that large finite size effects in the initialization affect the location of the BBP transition. Whether this happens also in this case is unclear and deserves further investigations.

Let us conclude by commenting on the limitations and possible generalizations of the results presented in this paper. A key element of the phase retrieval model is the existence of a phase at small α in which the best achievable generalization error is not better than random guess. This arises in models that present a symmetry, such as the ±1 symmetry due to the absolute value in the phase retrieval problem. The picture shown in Fig. 1 where the threshold states are characterized using the gradient-flow dynamics in a variant of the model with randomized labels is enabled by this symmetry and the rest of our analysis relies on the simplifications stemming from this. We expect that the BBP picture presented in this work generalizes to all learning where at small sample complexity error is not better than random guess. Of course working out the details, even for two layer neural networks, is an open problem of interest for future work. In most learning problems a naive linear regression often achieves a better than random guess performance, and thus an extension of our theory for problems of that type (lacking a symmetry) would be an interesting direction for future work. Figure 5: Fraction of simulations with gradient descent that achieve zero test error for various input dimensions N . The number of simulations increases with the smaller input dimension in order the account for the larger fluctuations. The number varies from 10000 simulations for N = 4192 to 100000 simulations for N = 64. The learning rate is 0.0002 and the simulation runs until either the loss goes below 10 -8 or the number of steps exceeds 1500 log 2 N iterations. The logarithmic term in the time accounts for the fact that the average initial overlap with the ground truth is ∼ 1/

√

N . In the inset we zoom in the region where all the simulations converge showing the indication of a crossing with large fluctuations due to finite size effects.

Broader Impact

Our work is theoretical in nature, and as such the potential societal consequence are difficult to foresee. We anticipate that deeper theoretical understanding of the functioning of machine learning systems will lead to their improvement in the long term.

A BBP transition and phase transition of spectra initialization

In non-convex estimation problems, such as Phase Retrieval, big advantages follow from the development of well tailored spectral methods to be used as initialization step. Recently the outcomes of one such spectral method widely used in Phase Retrieval and based on the construction of a data matrix

M i,j = 1 αN αN m=1 T (Y m )u m,i u m,j = 1 α αN m=1 T (Y m )x m,i x m,j (18)
from sensing vectors u u u m , with elements of order one (or sensing vectors x x x m on the unitary sphere, according to our definition) and measurements Y m has been exactly derived [START_REF] Yue | Phase transitions of spectral initialization for high-dimensional nonconvex estimation[END_REF]. The method involves a pre processing function T (Y m ) that can be optimized to further improve the results. Once the data matrix is constructed the eigenvector v 1 corresponding to the largest eigenvalue λ 1 can be used as an estimator of the signal W * .

To obtain the performances of this kind of spectral initialization it is assumed [START_REF] Yue | Phase transitions of spectral initialization for high-dimensional nonconvex estimation[END_REF] that the measurements Y m are independently drawn according to a density function conditional on y m = x m W * associated to the particular acquisition process, and it is recalled that y m are themselves Gaussian random variables, due to the definition of the problem. Finally in the large N limit the empirical average used to construct the data matrix can be replaced by the expected value E Y,y over these two distributions.

The result [START_REF] Yue | Phase transitions of spectral initialization for high-dimensional nonconvex estimation[END_REF] goes as follows. Given two functions defined as

Ψ α (λ) ≡ λ 1 α + E Y,y T (Y ) λ -T (Y ) (19)
and

Φ(λ) = λE Y,y T (Y )y 2 λ -T (Y ) (20)
with λ > T (Y ), and given

ξ α (λ) = Ψ α (max{λ, λ}) (21) with λ = arg min Ψ α (λ) , (22)
the two largest eigenvalues of M, λ 1 and λ 2 , are such that

λ 1 → P ξ α (λ * α ) , (23)
with λ * α the solution of ξ α (λ) = Φ(λ), and

λ 2 → P Ψ α ( λ). (24)
A phase transition occurs at the largest α such that λ 1 = λ 2 , which can be evaluated by imposing

Ψ α (λ * α ) = 0 or equivalently 1 α = E Y,y T (Y ) 2 ( λ -T (Y )) 2 , ( 25
)
which corresponds to Ψ α ( λ) = 0 as at that point λ * α = λ. At larger α the largest eigenvalue pops out from the spectrum bulk (λ 1 = λ 2 ) and the corresponding eigenvector develops a finite correlation with the signal, in a phenomenon called BBP transition [START_REF] Baik | Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices[END_REF], hence the definition of α BBP when this occurs.

It is interesting to note that the structure of the data matrix is closely reminiscent of the structure of the first term in the Hessian of our problem (see Eq. ( 4)). Indeed incidentally the original idea of this spectral method for initialization can be traced back to the study of Hessian's principal directions [START_REF] Li | On principal hessian directions for data visualization and dimension reduction: Another application of stein's lemma[END_REF]. In particular we observe that

-H i,j (W W W ) = - 1 2 αN m=1 ∂ 2 ∂ ŷ2 m (ŷ m , y m )x m,i x m,j + µδ i,j = M i,j + µδ i,j , (26)
provided that the pre processing function is

T (Y ) = - α 2 ∂ 2 ∂ ŷ2 (ŷ, y) = - α 2 ∂ 2 ŷ (ŷ, y) . ( 27
)
Note that we consider a case for which measurements Y have a one to one correspondence with y (i.e. Y = y 2 ). Moreover in the problem discussed empirical averages involve not only measurements Y but also estimated ŷ = xW in correspondence of some W of interest. Therefore the expected value, E ŷ,y , should be taken over the relative joint probability distribution function P (ŷ, y). In conclusion, the results just mentioned tell immediately what is the largest α (i.e. α BBP ):

1 α BBP = E ŷ,y α BBP ∂ 2 ŷ (ŷ, y) 2 λ + α BBP ∂ 2 ŷ (ŷ, y) 2 , ( 28
)
before the smallest eigenvalue of the Hessian pops out from the spectrum bulk, being also associated to an eigenvector with finite projection on the signal to be detected.

In this case we are focusing on performance of a gradient descent dynamics, which for mean-field spin-glasses naturally gets stuck on what are called threshold states [START_REF] Cugliandolo | Dynamics of glassy systems[END_REF]. We argue that the gradient descent dynamics applied to Phase Retrieval, when retrieval fails, will also approach threshold states which are mainly characterized by their property of being marginal, i.e. the smallest eigenvalue of their Hessian is null. This qualifies the relevant W as a typical configuration belonging to threshold states and P (ŷ, y) as the joint probability distribution at threshold states. Moreover it allows to introduce the marginal condition λ 2 = λ 1 = -µ, which can be re-expressed by equating λ 2 = Ψ α ( λ) to -µ:

µ = λE ŷ,y α BBP ∂ 2 ŷ (ŷ, y)y 2 2 λ + α BBP ∂ 2 ŷ (ŷ, y) . ( 29
)
Finally the definition of the spherical parameter (3) in the main text must be considered to close the system of equations

µ = α BBP 2 E ŷ,y (∂ ŷ (ŷ, y)ŷ) (30)
to be used to determine α BBP , λ, µ in correspondence of P (ŷ, y) from threshold states.

The resulting picture is as follows. In the small α regime, gradient descent will systematically approach threshold states and remain stuck there. However starting from α BBP in the Hessian of these states, that is otherwise marginal, an isolated eigenvalue pops out from the bulk immediately becoming negative. Moreover the eigenvector associated to such negative direction, naturally followed by gradient descent, has a finite overlap with the signal. Therefore, we argue, it is from that point on that the signal should be easily retrieved.

B Replica Analysis

The field of physics of disordered systems has developed numerous tools to deal with random systems [START_REF] Mézard | Spin glass theory and beyond[END_REF][START_REF] Nishimori | Statistical physics of spin glasses and information processing: an introduction[END_REF]. At an abstract level of thinking, using those tools means that we identify the inference problem with a physical systems that subject to a certain potential. The randomness comes from the having a dataset made of random projections. The estimator is mapped into a spherical spin and the loss function becomes the energy -or Hamiltonian. Finally the system the temperature in which the system lives is sent to zero and the system tend to the lowest energy, herefore minimizing the loss. The ground truth in the inference problem becomes equivalent to a minimum planted in the energetic landscape [START_REF] Krzakala | Hiding quiet solutions in random constraint satisfaction problems[END_REF]. For example this formulation is equivalent to a physical system that before the experiment is liquid, but as we cool it down it can become either a crystal -an ordered solid -or a glass -an amorphous solid. Finding the crystal means reconstructing the signal.

B.1 Partition function

Moving to the mathematical aspects of the problem. We define a Gibbs distribution associated with the problem and evaluate its normalization constant -the partition function Z. The partition function, and in particular its logarithm divided by the temperature -the free energy -, contains all the information we aim to understand in the problem. By taking the proper derivatives, and possibly add an external field, we can compute relevant macroscopic properties such as: average overlap with the ground truth, average loss achieved. In disordered system we have to consider the additional complication given by the randomness. Therefor, we need to consider the average free energy that is the average of the logarithm of a high-dimensional integral, that can be done by the simple observation:

log Z = lim n→0 Z n -1 n . (31)
Which for arbitrary n is not simpler than computing the logarithm, but it is much simpler if n ∈ N and we perform an analytic continuation to n ∈ R. Under this -replica trick -the average of the logarithm is equivalent to compute the average of the nth moment of the partition function and take the limit. Formally the nth moment correspond to the partition function of n identical -replicatedsystem that do not interact but with the same realization of the disorder.

The problem is now computing the moments of the partition function which in general is prohibitive and we have to use an ansatz on the specific form of the solution, in particular we use the so called replica symmetry breaking ansatz [START_REF] Mézard | Spin glass theory and beyond[END_REF]. This largely reduces the number of parameters. Finally the average free energy can be evaluated by set of saddle point equations.

We can now move to the analysis. The partition function already defined in the main text is

Z = S N -1 ( √ N ) DW W W S exp - β 2 µ ( x x x µ , W W W S , x x x µ , W W W T ) (32)
and consider its nth moment -the replicated partition function -

Z n = E {x x xµ} S n(N -1) ( √ N ) n a=1 DW W W a exp - β 2 µ (| x x x µ , W W W S |, | x x x µ , W W W T |) . (33)
This is formally equivalent to have n independent systems. Introduce the overlaps with the projector r (a) µ = x x x µ , W W W a with indices a = 0, . . . , n where a = 0 is the overlap with the ground truth and the others are the overlaps with the estimators of the n systems. Introduce those quantities in the replicated partition function via Dirac's deltas using their Fourier transform

Z n = E {x x xµ} a DW W W a D(r, r) exp - β 2 a,µ (r (a) µ , r (0) µ ) + i a,µ r(a) µ r (a) µ -i a,µ x x x µ , r(a) µ W W W a = = a DW W W a D(r, r) exp    - β 2 a,µ (r (a) µ , r (0) µ ) + i a,µ r(a) µ r (a) µ - 1 2N µ n a,b=0 r(a) µ r(b) µ W W W a , W W W b    . ( 34
)
Where D contains the normalization factor of the Fourier transform. We introduce the matrix of overlaps between estimators and ground ground truth

Q Q Q, Q ab = 1 N W W W a , W W W b
. This is done using the same idea of introducing delta function and gives a contribution N 2 log det Q Q Q to the action [START_REF] Zamponi | Mean field theory of spin glasses[END_REF][START_REF] Kurchan | Exact theory of dense amorphous hard spheres in high dimension. ii. the high density regime and the gardner transition[END_REF]. The equation can now be factorized in µ, so we can drop the µs from r and r. Observing that

e -1 2 ab r(a) r(b) Q ab = e 1 2 ab Q ab ∂ 2 ∂ha ∂h b e -i ha r(a) {ha=0}a , (35)
we can integrate over r and r, and write a simplified replicated partition function

Z n = n a≥b=0 dQ ab e N S(Q Q Q) with action S(Q Q Q) = 1 2 log det Q Q Q+α log exp   1 2 n a,b=0 Q a,b ∂ 2 ∂h a ∂h b   exp - β 2 a (h 0 , h a ) {ha=0}a . (36)
where the first term is an entropic term that accounts for the degeneracy of the matrix Q Q Q in the space of symmetric matrices. And the second term is an energetic term that accounts for the potential acting on the system.

Observe that so far we did not make any ansatz on the structure of the overlap matrix Q Q Q. In the next subsections we will consider the 1-step replica symmetry breaking ansatz (1RSB).

B.2 1 step replica symmetric breaking

The 1RSB scheme consists in making an ansatz on the structure of the overlap matrix Q Q Q [START_REF] Mézard | Spin glass theory and beyond[END_REF]. The assumption is that not all the replicated systems will have the same overlap -which correspond to the replica symmetric ansatz -but the systems are clustered. Systems inside the same cluster will have a larger overlap, systems outside the cluster will have a smaller overlap. This translates into the following parameters: q 1 overlaps inside the same cluster, q 0 overlaps in different clusters, x dimension of the clusters, finally m the overlap with the signal. Schematically we have

Q Q Q =         1 m m m 1 q 1 q 1 q 1 1 q 1 q 1 q 1 1 q 0 m q 0 1 q 1 q 1 q 1 1 q 1 q 1 q 1 1         (37)
with Q Q Q of dimension (n + 1) × (n + 1), and the inner matrices of dimension x × x.

The analysis proceed as in a standard way: we derive the 1RSB free energy with the associated saddle point equations and move to the zero temperature limit [START_REF] Mézard | Spin glass theory and beyond[END_REF][START_REF] Nishimori | Statistical physics of spin glasses and information processing: an introduction[END_REF], we impose the marginal stabilitycorresponding to the threshold states - [START_REF] Monasson | Structural glass transition and the entropy of the metastable states[END_REF], finally we derive the label distribution [START_REF] Franz | Universality of the sat-unsat (jamming) threshold in non-convex continuous constraint satisfaction problems[END_REF].

For notational convenience we call γ a (x) the probability density function of a Gaussian with zero mean and (co-)variance a, and we use the symbol to indicate the convolutions, i.e. f g(x) = R f (x -t)g(t)dt.

We can plug Eq. (37) into Eq. ( 36) and obtain 1RSB formulation of the action.

1 n S 1RSB (Q Q Q) = 1 2 log(1 -q 1 ) + 1 2x log 1 -q 1 + x(q 1 -q 0 ) 1 -q 1 + 1 2 q 0 -m 2 1 -q 1 + x(q 1 -q 0 ) + + α R 2 dŷdyγ Σ Σ Σ (ŷ, y) 1 x log γ q1-q0 γ 1-q1 e -β 2 (ŷ,y) x . (38)
and using Eq. ( 31) the free energy is -1 nβ S 1RSB . We can now take derivatives to find the saddle point equations.

1 x 1 1 -q 1 - 1 1 -q 1 + x (q 0 -q 1 ) + q 0 -m 2 [1 -q 1 + x (q 0 -q 1 )] 2 = = α R 2
dŷdyγ Σ (ŷ, y)e -x f (x,ŷ,y) γ q1-q0 e x f (1,ŷ,y) d dŷ f (1, ŷ, y)

2 ;

(39)

q 0 -m 2 [1 -q 1 + x (q 0 -q 1 )] 2 = α R 2 dŷdyγ Σ (ŷ, y) d dh f (x, ŷ , y) 2 ;
(40)

m 1 -q 1 + x (q 0 -q 1 ) = α R 2 dŷdyγ Σ (ŷ, y) d 2 dŷdy f (x, ŷ, ŷ); (41)
where we define

f (1, ŷ, y) = log γ 1-q1 e -β 2 (ŷ,y) ; f (x, ŷ, y) = 1 x log γ q1-q0 γ 1-q1 e -β 2 (ŷ,y) x = 1 x log γ q1-q0 e xf (1,ŷ,y) ; (ŷ, y) = (ŷ 2 -y 2 ) 2 ;
and Σ is a 2 × 2 covariance matrix with entries Σ 11 = 1, Σ 12 = m and Σ 22 = q 0 .

Zero temperature and parameter ansatz As the zero temperature goes to zero the order parameters q 1 and x needs to be rescaled with the temperature as q 1 ≈ 1 -χT and x ≈ zT with χ, z ∼ O(1). Instead m and q 0 will not be affected by the limit. The reason why some parameters need to be rescaled is that as long as there is a positive temperature the replicated systems exploit this thermal energy to fluctuate in the basin of the minimum, therefor their overlap q 1 is given by average overlap in the basin of attraction. When the temperature drops to zero the thermal goes to zero and all the system shrink to a point. χ represent the fluctuation that systems can have when they receive an infinitesimal amount of thermal energy. With the same physical reason, also the cluster itself shrinks to a point. The rate of convergence to the point is given by z.

Under those observation we obtain the 1RSB free energy

-f (Q Q Q) ≈ β 1 1 2z log χ + z(1 -q 0 ) χ + 1 2 q 0 χ + z(1 -q 0 ) + + α y R dyγ Σ Σ Σ (ŷ, y) log γ 1-q0 e -y 2 V (ŷ,y|χ) ŷ=0 . (42)
and the corresponding saddle point equations from Eqs. (39-41).

However, the solution of those equation will lead to the global minimum of the loss, while we are interested in the threshold states. We follow the idea of [START_REF] Monasson | Structural glass transition and the entropy of the metastable states[END_REF] where z is used as a Lagrange multiplier that selects the threshold states. z is fixed by imposing the marginal stability condition, i.e. that the spectrum of the Hessian of the minima in consideration touches zero, following [START_REF] Franz | Universality of the sat-unsat (jamming) threshold in non-convex continuous constraint satisfaction problems[END_REF] this is given by

1 = α 4 χ 2 R dy 1 √ 2π e -y 2 2 γ 1 V (ŷ, y|χ) 2 e -z 2 V (ŷ,y|χ) γ 1 e -z 2 V (ŷ,y|χ) ŷ=0 . ( 43
)
At the threshold the overlap with the signal is zero, m = 0, and for symmetries of the problem also q 0 = 0. In fact m = 0 is always a solution of Eq. ( 41), and if m = 0 then q 0 = 0 is also solution of the second saddle point equation, Eq. ( 40), as in the RHS the Gaussian becomes degenerate in h and d h f (x, h, h 0 )| h=0 being an odd function in h. Therefore we can restrict the equations to just the one for χ, Eq. (39) that becomes

β 2 z 1 χ - 1 χ + z = α R dyγ 1 (y)e -x f (x,0,y) R dŷγ 1 (ŷ) e x f (1,ŷ,y) d dŷ f (1, ŷ, y) 2 . (44)
Rewriting these equation expliciting the convolutions we obtain the equations presented in the main text. With those element we can obtain the label distribution presented in the main text.

C Additional details on the numerical experiments C.1 Dynamics with shuffled labels

In the main text we showed that the dynamics is always attracted by the threshold states before (possibly) finding the good direction that leads to the optimal solution. To visually complement the explanation, in Fig. 6 we consider the dynamics in the same setting of Fig. 2, same dataset and same initialization, but using shuffled labels. The upper figure shows the density of eigenvalues during the evolution while the lower figure considers the histogram for five cuts at iteration 1, 453, 906, 1208 and 1510. We remark that in the initial stages of their evolution the two simulations follow the same dynamics (iteration 1, 453, and 906). Only after a transient the simulation with the correct labels shows the presence of the BBP transition (iteration 1208) and converge to the solution (iteration 1510). Instead, the simulation with the shuffled labels keeps approaching the threshold states for the rest of the simulation time tending to the marginal condition. We remark that in the upper figure the smallest (dashed orange line) and second smallest eigenvalues (solid green line) that are always very close, as it is expect since in this case the distribution of the eigenvalues always forms a bulk.

C.2 Numerical threshold

On the choice of the learning rate. The dynamics considered in the main text is gradient flow, while numerically we rely on the discretized version of the algorithm, namely gradient descent. In order to have agreement between the theoretical analysis and the simulations we must consider a learning rate sufficiently small to reduce the discrepancy between the two algorithms. We consider different input dimensions and tested different learning rates, in Fig. 7 we report this process for N = 1024. In the On the left we show the overlap W W W • W W W * /N in time for 10 different simulations at α = 7. The different line styles refer to different learning rates. On the right we show the fraction of success, as reported in Fig. 5, as we change the learning rate. In the simulations we adapted the stopping criterion to the learning rate so that the curve get to the same simulation time, more specifically for the left panel the stopping times are : 1000 log 2 N for η = 2 × 10 -4 , 2000 log 2 N for η = 1 × 10 -4 , and 4000 log 2 N for η = 5 × 10 -5 .

On the nature of the transition. Our analysis assumes a first order transition between: the easy phase, where gradient flow with high probability solves the problem; and the hard phase, where gradient flow with high probability does not find the optimal solution and this solution is information theoretically achievable. According to this hypothesis threshold identifies a jump in the overlap between estimator and ground truth, which defines the first order transition. We bring as additional support on the matter Fig. 8, where we show the final overlaps reached by the simulations for different size of the dataset α. As α increases, the figures show that densities become bimodal and finally concentrate on 1 meaning that the problem is easy. This is an indicator of a first order transition. Contrary to the second order (continuous) transition, where the distribution of the overlaps should be unimodal.

Figure 1 :

 1 Figure 1: Evolution of the training loss for systems of size N = 16384 with number of samples αN .The left panel has labels created using a teacher, while the right panel has random labels constructed using a Gaussian distribution with variance matching the teacher-labels. The left figure shows that also the simulations with a teacher approach the threshold energy before transiting to the global minimum, if α is large enough.

Figure 2 :

 2 Figure2: Properties of the Hessian for phase retrieval of a system of size N = 2048 at α = 10. On the left figure, we show the evolution of the density of the bulk of the eigenvalues, from zero density in white to high density in blue, the smallest eigenvalue in orange, and the second smallest in green. The picture shows that a BBP transition occurs when the training loss approaches the threshold energy. The right panel depicts the evolution of the training loss (purple) and the generalization loss (dark blue) in time. The training loss rapidly approaches a plateau at the level of the threshold states (black dashed line) and converges towards the teacher after the BBP transition. Namely, when the smallest eigenvalue detaches from the rest of the bulk. The bottom panels show the distribution of the eigenvalues at five different instants in the dynamics. At iteration ≈ 1210 we cross the threshold energy and we observe an isolated eigenvalue detaching from the bulk.

Figure 3 :

 3 Figure 3: The three images show the occurring of a BBP transition at the moment when the training loss crosses the threshold energy. The images investigate the BBP as the input dimension is increased form N = 512 to N = 8192. At the BBP transition the smallest eigenvalue pops out of the bulk of eigenvalues and the associated eigenvector contains information on the signal. The left figure is the difference of the smallest eigenvalues (the second with the first in solid line, the third with the second in dotted line). The central image shows the overlap of the first eigenvector (full) and second eigenvector (dotted) with the signal. The transition appears to shift as in the left image. Finally on the right the fluctuations of the overlap first eigenvector with the teacher are shown, the peak corresponds to the transition.

Figure 6 :

 6 Figure6: Comparison of the properties of the Hessian for phase retrieval between the problem with the correct labels and the problem with shuffled labels. The system simulated has size N = 2048 and α = 10. The database considered in the simulation in the upper figure is the one used for Fig.2with the shuffled labels. The figures represent the same quantities plotted in Fig.2and are intended to make the comparison between the evolution in the two settings. The figure below represents the distribution of the eigenvalues during the dynamics for the problem with correct and shuffled labels, respectively in blue and orange.

Figure 7 :

 7 Figure 7: Choice of the learning rate for the phase retrieval problem with input dimension N = 1024.On the left we show the overlap W W W • W W W * /N in time for 10 different simulations at α = 7. The different line styles refer to different learning rates. On the right we show the fraction of success, as reported in Fig.5, as we change the learning rate. In the simulations we adapted the stopping criterion to the learning rate so that the curve get to the same simulation time, more specifically for the left panel the stopping times are : 1000 log 2 N for η = 2 × 10 -4 , 2000 log 2 N for η = 1 × 10 -4 , and 4000 log 2 N for η = 5 × 10 -5 .

Figure 8 :

 8 Figure 8: The five figures show the final overlaps between estimator and teacher for a large number of simulations at different values of α. The figure refer to different input dimensions, from left to right : 256, 512, 1024, 2048, and 4096.

Acknowledgments and Disclosure of Funding

We thank Carlo Lucibello for precious discussions and Federico Ricci-Tersenghi for interesting comments on the manuscript. We acknowledge funding from the ERC under theEuropean Union's Horizon 2020 Research and Innovation Programme Grant Agreement 714608-SMiLe;from the French government under management of Agence Nationale de la Recherche (ANR) grant PAIL,and grant "Investissements d'Avenir" LabEx PALM (ANR-10-LABX-0039-PALM) (StatPhysDisSys) and(ANR-19-P3IA-0001) (PRAIRIE 3IA Institute); and from the Simons Foundation collaboration Crackingthe Glass Problem (#454935, Giulio Biroli).


MAMA ABDELLI Akram Ben Aissa

Mama Abdelli
email: abdelli.mama@gmail.com

AND Alessandro Duca
email: alessandro.duca@univ-grenoble-alpes.fr

WELL-POSEDNESS AND EXPONENTIAL DECAY FOR THE EULER-BERNOULLI BEAM CONVEYING FLUID EQUATION WITH NON-CONSTANT VELOCITY AND DYNAMICAL BOUNDARY CONDITIONS

Keywords: 2010 Mathematics Subject Classification. 35L05, 35L15, 35L70, 93D15 Beam, Lyapunov's direct method, Uniform stability, Exponential stability, Internal fluid

Introduction

In this paper, we study the dynamics of a tubular structure, a flexible pipe, which conveys a moving fluid. We assume that the pipe diameter is negligible compared to its length and we model the system by a one dimensional problem. In this context, the motion of the structure is driven by the following Euler-Bernoulli beam conveying fluid equation (m p + 2m f ) w tt (t) + EIw xxxx (t) -(T -2m f V (t) 2 ) w xx (t) +cw t (t) + 2m f V t (t)w x (t) + 4m f V (t)w xt (t) = 0, 0 < x < L, t > 0, (1.1) endowed with the boundary conditions

w xx (0, t) = w xx (L, t) = w(0, t) = 0, t > 0, EIw xxx (L, t) -(T -2m f V (t) 2 )w x (L, t) + 2m f V (t)w t (L, t) = 0, t > 0. (1.2)
The parameters L and m p are the length and the mass per unit length of the flexible pipe, while EI and T are its bending stiffness and its tension. The terms m f and V (t) are the mass per unit length and the velocity of the internal fluid. We assume that the velocity V (t) is strictly positive or strictly negative (we refer to Section 4 for further details on the general situation). The solution of the system (1.1) with boundary conditions (1.2) represents the displacement of the flexible pipe at the position x and time t (see Figure 1). We study the well-posedness of (1.1)-(1.2) and the exponential stability of its solutions.

The equation (1.1)-(1.2) has been deduced when the velocity V (t) is constant by Liu et al in [START_REF] Liu | Modeling and boundary control of a flexible marine riser coupled with internal fluid dynamics[END_REF] in presence of the environmental disturbances and a boundary control authors consider the problem of a flexible marine riser and they compute the equation via a suitable energy functional. Finally, they prove the exponential decay of the solutions by Lyapunov method in presence of an additional boundary condition w x (0, t) = 0.

In our work, we assume that the tension T is larger than a specific value T * depending on the parameters of the problems. By exploiting such hypotheses, we firstly ensure the well-posedness of the (1.1)-(1.2). Secondly, we show the exponential stability of the energy of the solutions of (1.1)-(1.2) without considering additional boundary conditions.

Our stability result is obtained by exploiting a different energy functional from the one considered in [START_REF] Liu | Modeling and boundary control of a flexible marine riser coupled with internal fluid dynamics[END_REF] which validity is due to the assumption imposed on T . Such condition is not surprising from a practical point of view. It is reasonable to assume that the tension of the pipe has to be sufficiently strong, compared to the density and the velocity of the conveyed fluid, in order to have the stability. Otherwise, the fluid inside the pipe may dynamically interacts with its motion, possibly causing the flexible pipe to vibrate.

Euler-Bernoulli beam conveying fluid equations are found in many practical applications. They are used to model for instance risers of offshore platforms, pipes carrying chemical fluids, exhaust pipes in the engines, flue-gases stacks, air-conditioning ducts, tubes in heat exchangers and power plants, etc. A similar dynamics to (1.1) is studied in [START_REF] Khemmoudj | Stabilization of a viscoelastic beam conveying fluid[END_REF] by Khemmoudj where the internal damping cw t (t, x) is replaced by a viscoelastic term. There, the author considers suitable boundary conditions and he proves the exponential stability by Lyapunov method. In [START_REF] Conrad | On the stablization of a flexible beam with a tip mass[END_REF], Conrad et al. consider the equation w tt + w xxxx = 0 and they prove exponential stability in presence of specific dynamical boundary conditions. For other similar results, we refer to [START_REF] Canbolat | Boundary control of a cantilevered flexible beam with point-mass dynamics at the free end[END_REF][START_REF] Seghour | Control of a riser through the dynamic of the vessel[END_REF].

The paper is organized as follows. In Section 2, we present the well-posedness of the equation (1.1)-(1.2) in Theorem 2.2. In Section 3, we ensure our stability result in Theorem 3.1 by introducing a suitable Lyapunov functional.

Well-posedness of the problem

The aim of this section is to prove the existence and uniqueness of solutions for (1.1)-(1.2). To this purpose, we define the space

V = w ∈ H 2 (0, L), R w(0) = 0
and the Hilbert space

H = (w, v) w ∈ V, v ∈ L 2 (0, L), R
equipped with the norm • H induced by the scalar product

f 1 , f 2 H = L 0 ∂ 2 x w 1 ∂ 2 x w 2 + ∂ x w 1 ∂ x w 2 + v 1 v 2 dx, f 1 = w 1 , v 1 , f 2 = w 2 , v 2 ∈ H.
We define the family {A 0 (t)} t≥0 of operators in H such that

A 0 (t)f = A 0 (t) w v = v -EI mp+2m f w xxxx + T -2m f V 2 (t) mp+2m f w xx - 4m f V (t) mp+2m f v x ,
for every f = (w, v) in the domain D(A 0 (t)) defined by

D(A 0 (t)) = (w, v) ∈ H | w ∈ V ∩ H 4 , v ∈ V, w xx (0) = w xx (L) = 0, 2m f V (t)v(L) = -EIw xxx (L) + (T -2m f V 2 (t))w x (L) .
Let {B(t)} t≥0 be the family of bounded operators in H such that

B(t)f = 0 c mp+2m f v + 2m f Vt(t) mp+2m f w x , ∀f ∈ H, By imposing v = w t ,
)f = A(t)U (t, s)f.
Theorem 2.2 ensures the well-posedness of (2.1) endowed with the boundary conditions (1.2), and then of (1.1)-(1.2). The result is guaranteed when the tension T is sufficiently large with respect to the velocity V (t) and to the mass m f .

The following proposition is the core of the proof of Theorem 2.2. It ensures the existence of a unique evolution system on H generated by the family of operators (A 0 (t), D(A 0 (t))). Such result leads to Theorem 2.2 thanks to the properties of the family of operators B(t) . Proposition 2.3. Let V ∈ C 1 [0, +∞), R be a strictly positive or strictly negative function. Let T > 0 be such that T > 2m f sup t≥0 V 2 (t). There exists a unique evolution system

U (t, s) with 0 ≤ s ≤ t in H satisfying • U (t, s)D(A 0 (s)) = D(A 0 (t)) for every 0 ≤ s ≤ t; • for every s ≥ 0 and f ∈ D(A 0 (s)), the function t → U (t, s)f in H is continuously derivable with respect to the norm of H for t ≥ s and d dt U (t, s)f = A 0 (t)U (t, s)f.
Proof. Thanks to the assumptions on the velocity V (t), we can define the norm • t for every t ≥ 0 of H induced by the scalar product

f 1 , f 2 t = L 0 α(t)∂ 2 x w 1 ∂ 2 x w 2 + β(t)∂ x w 1 ∂ x w 2 + γ(t)v 1 v 2 dx, for every f 1 = (w 1 , v 1 ) and f 2 = (w 2 , v 2 ) ∈ H with α := EI 2m f , β(t) = T -2m f V 2 (t) 2m f , γ = m p + 2m f 2m f .
The assumptions on the potential V and on the tension T yield the existence of C > 1 such that

C -1 • H ≤ • t ≤ C • H , ∀t ≥ 0. (2.2)
The domain of A 0 (t) can be rewritten in terms of the parameters α and β(t) as follows

D(A 0 (t)) = (w, v) ∈ H | w ∈ V ∩ H 4 , v ∈ V, w xx (0) = w xx (L) = 0, V (t)v(L) = -αw xxx (L) + β(t)w x (L) .
1) Dissipative property. First, we prove that A 0 (t) is dissipative for every t ≥ 0 in H with respect to the scalar product •, • t . Let us denote

a 1 := - EI m p + 2m f , a 2 (t) := T -2m f V 2 (t) m p + 2m f , a 3 (t) := - 4m f V (t) m p + 2m f ,
For every f ∈ D(A 0 (t)), thanks to the relations

γa 1 = -α, γa 2 (t) = β(t), γ(t)a 3 (t) = -2V (t),
there holds

f, A 0 (t)f t = L 0 αw xx v xx + β(t)w x v x + γ a 1 w xxxx + a 2 (t)w xx + a 3 (t)v x v dx = -αw xxx (L)v(L) + β(t)w x (L)v(L) + L 0 α + γa 1 w xxxx + β(t) -γa 2 (t) w xx + γa 3 (t)v x v dx = -αw xxx (L) + β(t)w x (L) v(L) - L 0 2V (t)v x vdx = V (t)v(L) 2 -V (t)v(L) 2 = 0.
2) Surjectivity conditions. Second, fixed t ≥ 0, we ensure the surjectivity of the map

(λI -A 0 (t)) : D(A 0 (t)) ⊂ H -→ H
for every λ > 0. The property is equivalent to prove that, for every f * = (w * , v * ) ∈ H, there exists a unique solution f = (w, v) ∈ D(A 0 (t)) of the equation (λI

-A 0 (t))f = f * .
In other words, we need to study the existence of a unique solution of the following system of equations

         λw(x) -v(x) = w * (x), λv(x) -a 1 w xxxx (x) -a 2 (t)w xx (x) -a 3 (t)v x (x) = v * (x), x ∈ [0, L], V (t)v(L) = -αw xxx (L) + β(t)w x (L), w xx (0) = w xx (L) = w(0) = 0. (2.3) Solving the previous system is equivalent to study      λ 2 w -a 1 w xxxx -a 2 (t)w xx -a 3 (t)λw x = v * -a 3 (t)w * x + λw * , λV (t)w(L) -V (t)w * (L) = -αw xxx (L) + β(t)w x (L), w xx (0) = w xx (L) = w(0) = 0. Let g * = v * -a 3 (t)w *
x + λw * . The weak formulation of the previous system is given by

L 0 λ 2 wφ -a 1 w xx φ xx + a 2 (t)w x φ x -a 3 (t)λw x (x)φ dx -a 1 w xxx (L)φ(L) -a 2 (t)w x (L)φ(L) = L 0 g * φdx (2.4)
with φ ∈ V. The identity (2.4) can be rewritten as

L 0 λ 2 wφ -a 1 w xx φ xx + a 2 (t)w x φ x -a 3 (t)λw x (x)φ dx + γ -1 αw xxx (L) -β(t)w x (L) φ(L) = L 0 g * φdx,
Thanks to the validity of the boundary conditions in L, we have

L 0 λ 2 wφ -a 1 w xx φ xx + a 2 (t)w x φ x -a 3 (t)λw x φ dx -γ -1 λV (t)w(L)φ(L) = -γ -1 V (t)w * (L)φ(L) + L 0 g * φdx.
We consider the bilinear form •, • A 0 (t) in V such that, for every w, φ ∈ V,

w, φ A 0 (t) = L 0 λ 2 wφ -a 1 w xx φ xx + a 2 (t)w x φ x -a 3 (t)λw x φ dx -γ -1 λV (t)w(L)φ(L).
We notice that, for every w ∈ V,

w, w A 0 (t) = L 0 λ 2 w 2 -a 1 w 2 xx + a 2 (t)w 2 x dx - a 3 (t) 2 + γ -1 V (t) λw(L) 2 .
Thanks to the identity

a 3 (t) 2 + γ -1 V (t) = 0,
•, • A 0 (t) is a coercive bilinear form as -a 1 > 0 and a 2 (t) > 0. From the Lax-Milgram theorem, the weak formulation admits an unique solution. Thanks to the regularity of w ∈ H 4 ((0, L), R) and by using particular φ, it is possible to recover the boundary conditions in w and v is defined by (2.3) which is unique. Now, f = (w, v) ∈ D(A 0 (t)). This shows the surjectivity, for every λ > 0, of the map

(λI -A 0 (t)) : D(A 0 (t)) ⊂ H -→ H.
3) Conclusion. Thanks to the point 1), the operator A(t) is dissipative for t ≥ 0 with respect to the scalar product •, • t and then, for every λ > 0,

(λI -A 0 (t))ψ t ≥ λ ψ t , ψ ∈ D(A 0 (t)).

Now, we define a new norm of H

• t := m(t) • t , ∀t ≥ 0 by choosing a suitable strictly positive function m : R + -→ R + such that there exists C > 1 such that

C -1 • H ≤ • t ≤ • s ≤ C • H , ∀0 ≤ s ≤ t
(by keeping in mind the validity of the inequality (2.2)). Thus, we have

(λI -A 0 (t))ψ t ≥ λ ψ t , ∀λ > 0, ψ ∈ D(A 0 (t)).
Let us denote, for every t ≥ 0 and x ∈ [0, L],

f (t, x) = 1 + x L |V (t)| α -1 , g(t, x) = sign(V (t)) β(t) α x L .
We introduce the family of operators Γ t with t ≥ 0 such that

Γ t : (w, v) → w, f (t, •)v -g(t, •)w x .
We notice that ). The property is due to the the smoothness of the family of bounded operators B(t) which is due to hypotheses imposed on the velocity V (t).

Γ t (D(A 0 (t))) = (w, v) ∈ H | w ∈ V ∩ H 4 , v ∈ V, w xx (0) = w xx (L) = 0, v(L) = -sign(V (

Exponential stability of the problem

Main result

In this section, we ensure the stability result for the solution of the problem (1.1)-(1.2). To the purpose, we introduce the following assumptions, Assumptions I. Let V ∈ C 2 [0, +∞), R be a strictly positive or strictly negative function. Let T > 0 be such that

T > 2m f sup t≥0 V (t) 2 + max{T 1 , T 2 },
where we denote

T 1 = L 2 4 (m p + 2m f ) + 2 √ 2Lm f sup t≥0 |V (t)|
and

T 2 = c 2 L 2 8(c -m p -2m f ) + 2m f sup t≥0 |V t (t)V (t)|.
Assumptions I ensure that the tension T of the beam is sufficiently strong with respect to the velocity V and to the parameters of the problem. In this framework, the quantity T 2 -m f V 2 (t) > 0 for every t ≥ 0 and it is valid the well-posedness result provided in the previous section by Theorem 2.2. In addition, the assumption on T allows us to consider the following energy functional associated to the solutions of (1.1)-(1.2)

E(t) = 1 2 (m p + 2m f ) L 0 w 2 t (t) dx + EI 2 L 0 w 2 xx (t) dx + T 2 -m f V 2 (t) L 0 w 2 x (t) dx. (3.1)
We are finally ready to state the stability result of the problem 1.1-(1.2).

(3.2) E(t) ≤ k 0 E(s)e -k 1 (t-s) ∀ 0 ≤ s ≤ t.
The proof of Theorem 3.1 is established in the final part of the section by gathering different results. We firstly compute the time derivative of the energy functional d dt E(t). After, we introduce a suitable Lyapunov functional L(t) for which there exist

C 1 , C 2 > 0 so that C 1 E(t) ≤ L(t) ≤ C 2 E(t). Finally, we show the existence of C 3 > 0 such that d dt L(t) ≤ -C 3 L(t)
. This identity implies L(t) ≤ L(s)e -C 3 (t-s) and Theorem 3.1 is proved by gathering the previous results.

Some preliminaries

We start by recalling the following Sobolev-Poincaré inequality (see [START_REF] Adams | Sobolev spaces[END_REF] for further details).

d dt E(t) = -c L 0 w 2 t (t) dx -2m f V t (t) L 0 w t (t)w x (t) dx -2m f V t (t)V (t) L 0 w 2 x (t) dx.
Proof. By multiplying the first equation in (1.1) with w t and by integrating over (0, L), we obtain

L 0 w t (t) (m p + 2m f )w tt (t) + EIw xxxx (t) -(T -2m f V 2 (t))w xx (t) dx + 4m f V (t) L 0 w t (t)w xt (t) dx = -c L 0 w 2 t (t) dx -2m f V t (t) L 0 w t (t)w x (t) dx. (3.3)
We study each term appearing in first integral of (3.3)

(3.4) (m p + 2m f ) L 0 w t (t)w tt (t) dx = 1 2 (m p + 2m f ) d dt L 0 w 2 t (t) dx, (3.5) EI L 0 w t (t)w xxxx (t) dx = EIw t (L, t)w xxx (L, t) + EI 2 d dt L 0 w 2 xx (t) dx
and

-(T -2m f V 2 (t)) L 0 w t (t)w xx (t) dx = 2m f V 2 (t) -T w t (L, t)w x (L, t) + T 2 d dt L 0 w 2 x (t) dx -m f V 2 (t) d dt L 0 w 2 x (t) dx. (3.6)
The right-hand side of (3.6) can be rewritten as follows

(2m f V 2 (t) -T )w t (L, t)w x (L, t) + T 2 d dt L 0 w 2 x (t) dx -m f d dt V 2 (t) L 0 w 2 x dx + 2m f V t (t)V (t) L 0 w 2 x dx. (3.7)
Now, we investigate the remaining term in the left-hand side of (3.3)

(3.8) 4m f V (t) L 0 w t (t)w xt (t) dx = 2m f V (t)w 2 t (L, t).
By using (3.4)-(3.8) into (3.3) yields the following expression

d dt E(t) + EIw xxx (L, t) -(T -2m f V 2 (t))w x (L, t) + 2m f V (t)w t (L, t) w t (L, t) = -c L 0 w 2 t (t) dx -2m f V t (t) L 0 w t (t)w x (t) dx -2m f V t (t)V (t) L 0 w 2 x (t) dx.
The boundary conditions (1.2) complete the proof.

A suitable Lyapunov functional

Let us introduce the functionals

G 1 (t) = (m p + 2m f ) L 0 w(t)w t (t) dx, G 2 (t) = 2m f V (t) L 0 w(t)w x (t) dx, and
G(t) = G 1 (t) + G 2 (t). (3.9)
We define the Lyapunov functional L such that (3.10)

L(t) = E(t) + G(t).
Lemma 3.4. Let the problem (1.1)-(1.2) satisfy Assumptions I. There exist two positive constants ξ 1 > 0 and ξ 2 > 0 depending on m p , m f , T and V (t) such that, for all t ≥ 0,

(3.11) ξ 1 E(t) ≤ L(t) ≤ ξ 2 E(t).
Proof. Let P be the Poincaré's constant from Lemma 3.2 and α 1 > 0. By using the Young's and the Poincaré inequalities, we obtain

|G 1 (t)| ≤ α 1 2 m p + 2m f L 0 w 2 t (t) dx + P 2α 1 m p + 2m f L 0 w 2 x (t) dx. (3.12)
By setting α 2 = 1 P , we have

|G 2 (t)| ≤2m f |V (t)| α 2 P 2 + 1 2α 2 L 0 w 2 x (t) dx =2m f |V (t)| √ P L 0 w 2 x (t) dx (3.13)
In the last relation, we impose α 2 = 1 P in order to minimize the function α 2 → { P 2 α 2 + 1 2α 2 } which attains its minimum in √ P exactly when α 2 = 1 P . By combining (3.1), (3.12) and (3.13), we have

L(t) ≤ T 2 -m f V 2 (t) + P 2α 1 m p + 2m f + 2m f |V (t)| √ P L 0 w 2 x (t) dx + (m p + 2m f ) 1 + α 1 2 L 0 w 2 t (t) dx + EI L 0 w 2 xx (t) dx. (3.14)
For the lower bound, we can see that

L(t) ≥ T 2 -m f V 2 (t) - P 2α 1 m p + 2m f -2m f |V (t)| √ P L 0 w 2 x (t) dx + (m p + 2m f ) 1 - α 1 2 L 0 w 2 t (t) dx + EI L 0 w 2 xx (t) dx. (3.15)
We recall that P = L 2 2 and, thanks to Assumptions I,

T > 2m f sup t≥0 V (t) 2 + L 2 4 (m p + 2m f ) + 2 √ 2Lm f sup t≥0 |V (
+ 2m f )(1 -α 1 2 ) > 0 and T 2 -m f inf t≥0 V 2 (t) ≥ T 2 -m f sup t≥0 V 2 (t) > 0
thanks to (3.16). By combining (3.14) and (3.15), we obtain (3.11). The constants ξ 1 and ξ 2 are explicitly given by

ξ 1 = max 1 T 2 -m f sup t≥0 V (t) 2 T 2 -m f V 2 (t) + P 2α 1 m p + 2m f + 2m f sup t≥0 |V (t)| √ P , 2 1 + α 1 2 , 2 , ξ 2 = min 1 T 2 -m f inf t≥0 V 2 (t) T 2 -m f V 2 (t) - P 2α 1 m p + 2m f -2m f sup t≥0 |V (t)| √ P , 2 1 - α 1 2 , 2 .
Lemma 3.5. Let w be the solution of (1.1)-(1.2). The functional G defined by (3.9) satisfies

d dt G(t) = -EI L 0 w 2 xx (t) dx -(T -2m f V 2 (t)) L 0 w 2 x (t) dx -c L 0 w(t)w t (t) dx + 2m f V (t) L 0 w x (t)w t (t) dx + (m p + 2m f ) L 0 w 2 t (t) dx. Proof. We know that d dt G(t) = d dt G 1 (t) + d dt G 2 (t). Now, d dt G 1 (t) = (m p + 2m f ) L 0 w(t)w tt (t) dx + (m p + 2m f ) L 0 w 2 t (t) dx
We use (1.1) and we obtain

d dt G 1 (t) = -EI L 0 w xxxx (t)w(t) + (T -2m f V 2 (t)) L 0 w xx (t)w(t) dx -2m f V t (t) L 0 w x (t)w(t) dx -4m f V (t) L 0 w xt (t)w(t) dx -c L 0 w(t)w t (t) dx + (m p + 2m f ) L 0 w 2 t (t) dx.
We integrate by parts and d dt G 1 (t) becomes

-EIw xxx (L, t) -T -2m f V 2 (t) w x (L, t) + 2m f V (t)w t (L, t) w(L, t) -EI L 0 w 2 xx (t) dx -T -2m f V 2 (t) L 0 w 2 x (t) dx -2m f V t (t) L 0 w x (t)w(t) dx -2m f V (t) L 0 w xt (t)w(t) dx -c L 0 w(t)w t (t) dx + (m p + 2m f ) L 0 w 2 t (t) dx.
By using the boundary conditions, we obtain

d dt G 1 (t) = -EI L 0 w 2 xx (t) dx -T -2m f V 2 (t) L 0 w 2 x (t) dx -2m f V t (t) L 0 w x (t)w(t) dx -2m f V (t) L 0 w xt (t)w(t) dx -c L 0 w(t)w t (t) dx + (m p + 2m f ) L 0 w 2 t (t) dx. (3.17) Now, we compute d dt G 2 (t) d dt G 2 (t) = 2m f V (t) L 0 w xt (t)w(t) dx + 2m f V (t) L 0 w x (t)w t (t) dx + 2m f V t (t) L 0 w(t)w x (t) dx. ( 3
d dt L(t) = -c L 0 w 2 t (t) dx -2m f V t (t)V (t) L 0 w 2 x (t) dx -EI L 0 w 2 xx (t) dx -(T -2m f V 2 (t)) L 0 w 2 x (t) dx -c L 0 w(t)w t (t) dx + (m p + 2m f ) L 0 w 2 t (t) dx. (3.19)
By using Young's inequality and Poincaré inequality, for every α 1 > 0, we have

-c L 0 w(t)w t (t) dx ≤ c α 1 2 L 0 w t (t) 2 dx + 1 2α 1 c L 0 w(t) 2 dx ≤ c α 1 2 L 0 w t (t) 2 dx + P 2α 1 c L 0 w x (t) 2 dx. (3.20)
Now, thanks to Assumptions I, there exists δ ∈ (0, c -m p -2m f ) such that

T > 2m f sup t≥0 V (t) 2 + c 2 P 4(c -m p -2m f -δ) + 2m f sup t≥0 |V t (t)V (t)|. (3.21) We set α 1 = 2 c-mp-2m f -δ c in (3.20) that we use in (3.19) in order to obtain d dt L(t) ≤ -γ 0 L 0 w 2 t (t) dx -γ 1 L 0 w 2 x (t) dx -EI L 0 w 2 xx (t) dx where γ 0 = c -c α 1 2 -m p -2m f , γ 1 = T -2m f sup t≥0 |V t (t)V (t)| -2m f sup t≥0 V (t) 2 -c P 2α 1 .
We notice that γ 0 > 0 thanks to the choice of α 1 and γ 1 > 0 thanks to (3.21)

Conclusions

In Section 2 and Section 3, we studied the well-posedness and the stability of the solutions of (1.1)-(1.2). We assumed that the velocity V (t) is a sufficiently smooth function with constant sign, while the tension T is larger than a specific value T * depending on the parameters of the problems. By exploiting such hypotheses, we proved the well-posedness of (1.1)-(1.2). Secondly, we ensured the exponential stability by introducing a suitable Lyapunov functional.

The choice of considering V with constant sign is due to the following reason. When V vanishes for some time, the dynamical boundary condition in (1.2) becomes a statical boundary condition. In this case, the problem lacks of a boundary condition on w t (t, x) and it is not clear which "natural" boundary condition appears in such a context. Nevertheless, the stability result from Section 3 could still be valid, at least from a formal point of view.

Finally, the assumption on T is not so surprising when we think to the nature of the problem modeled by (1.1)-(1.2). It is reasonable to assume that the tension of the pipe has to be sufficiently strong, compared to the density and the velocity of the conveyed internal fluid, in order to have the stability. From this perspective, it could be interesting to explore this phenomenon further, at least from a numerical point of view. One could seek for evidences of instability phenomena when the tension T is too low.

Figure 1 .

 1 Figure 1. The figure represents the displacement of a flexible pipe which conveys a moving fluid.

Theorem 3 . 1 .

 31 Let the problem (1.1)-(1.2) satisfy Assumptions I. There exist two positive constants k 0 and k 1 , such that for any solution w of (1.1)-(1.2), the energy functional (3.1) decays as follows

Lemma 3 . 2 . 2 . 3 . 3 .

 32233 (Sobolev-Poincaré inequality) For every v ∈ {w ∈ H 1 (0, L) | w(0) = 0}, there holds L 0 v 2 dx ≤ P L 0 v 2 x dx, with P := L 2Lemma Let w be a solution of the problem (1.1)-(1.2). The energy functional defined by (3.1) satisfies

Finally, by combining ( 3 .

 3 [START_REF] Lazzari | On the exponential decay of the Euler-Bernoulli beam with boundary energy dissipation[END_REF]) and (3.24), we obtain (3.2) with k 0 = ξ 2 ξ 1 , which completes the proof.

  Now, we are ready to prove our main result.

	.18)
	By gathering (3.17) and (3.18), we achieve the claim.
	Proof of Theorem 3.1
	Proof of Theorem 3.1. By the definition (3.10), we have d dt L(t) = d dt E(t) + d dt G(t). Thanks
	to Lemma 3.3 and Lemma 3.5, we obtain

  T -2m f inf t≥0 V (t) 2 , 2 .

					. Finally, we
	see that
	(3.22)		d dt	L(t) ≤ -ϑE(t)
	with
	ϑ := min	2γ 0 m p + 2m f	,	2γ 1
	A combination of (3.11) and (3.22) gives
	(3.23)	d dt	L(t) ≤ -k 1 L(t), ∀t ≥ 0,
	where k 1 = ϑ ξ 2 . We integrate (3.23) over (s, t) and
	(3.24)

L(t) ≤ L(s)e -k 1 (t-s) , ∀t ≥ s.

Acknowledgments. The authors are grateful to Kaïs Ammari for suggesting them the problem. The third author also thanks the colleagues Romain Joly and Nabile Boussaïd for the fruitful discussions on the well-posedness of the problem. The third author has been funded by the project ISDEEC ANR-16-CE40-0013.