extfieldsofstudy:
Computer Science
Mathematics
Physics

text:
Landscape Complexity for the Empirical Risk of Generalized Linear Models

We present a method to obtain the average and the typical value of the number of critical points of the empirical risk landscape for generalized linear estimation problems and variants. This represents a substantial extension of previous applications of the Kac-Rice method since it allows to analyze the critical points of high dimensional non-Gaussian random functions. We obtain a rigorous explicit variational formula for the annealed complexity, which is the logarithm of the average number of critical points at fixed value of the empirical risk. This result is simplified, and extended, using the non-rigorous Kac-Rice replicated method from theoretical physics. In this way we find an explicit variational formula for the quenched complexity, which is generally different from its annealed counterpart, and allows to obtain the number of critical points for typical instances up to exponential accuracy.

1 Introduction and main results

Introduction
Characterizing the landscape of the empirical risk is a key issue in several contexts. Many current machine learning problems are both non-convex and high-dimensional. In these cases, the analysis of optimization algorithms, such as gradient descent and its stochastic variants, represents a very hard challenge. In recent years, there has been a series of works that developed a landscape-based approach to tackle this challenge.
The key idea is to study the statistical properties of the empirical risk landscape, and to use these findings to obtain results on the performance of algorithms. Without the aim of being exhaustive this research avenue includes analysis of the landscape of neural networks, matrix completion, tensor factorization and tensor principal component analysis [Fyodorov, 2004, Fyodorov and Nadal, 2012, Kawaguchi, 2016, Soudry and Carmon, 2016, Ge et al., 2016, Freeman and Bruna, 2019, Bhojanapalli et al., 2016, Park et al., 2017, Du et al., 2018, Ge and Ma, 2017, Lu and Kawaguchi, 2017, Ling et al., 2019, Ben Arous et al., 2019, Ros et al., 2019, Sarao Mannelli et al., 2019b,a, Ben Arous et al., 2020, Biroli et al., 2019. The majority of these works identifies the region of parameters where the landscape is "easy", i.e. it focuses on the regime where there shouldn't be any bad local minima and it proves that indeed there are none. However, gradient descent and other landscape-based algorithms are often observed to work even very far from the region described above where the landscape can be proved mathematically to be "easy". A possible reason is that the bounds obtained rigorously are not tight enough. Another, more interesting, is that the landscape is "hard", i.e. spurious minima are present, but their basins of attraction are small and the dynamics is able to avoid them [Sarao Mannelli et al., 2019b,a].
Here we develop a general method that allows to study and focus directly on the "hard" regime, where the empirical risk displays a huge number of bad minima. Our aim is to obtain explicit formulas for the

Main results
We consider two classes of high-dimensional random functions. The first one is a kind of energy that arises in a simple model of neural networks (the perceptron, cf. Engel and Van den Broeck [2001]) and in mean-field glass models [Franz and Parisi, 2016]: where φ : R → R is a smooth activation function (the hypotheses on φ are precised later), x ∈ S n−1 , the unit sphere in n dimensions, and ξ µ are i.i.d. random variables generated from the standard Gaussian distribution in R n . The second class of functions we will consider are related to the loss functions of generalized linear models (GLMs) Nelder and Wedderburn [1972], McCullagh [1984]. In generalized linear estimation an observer has to infer a hidden vector x ∈ S n−1 from the observation of the mdimensional output vector Y = {φ(ξ µ · x )} m µ=1 . In this sense, the GLMs generalize the usual linear regression by allowing the output function to be non-linear 1 . We consider here random GLMs, meaning that the data (or measurement) matrix ξ is taken random, with an i.i.d. standard Gaussian distribution, and we assume that the function φ and the data matrix ξ are given to the observer. This naturally leads to the mean square loss L 2 : GLMs (and their random versions) arise in many different areas of statistics, such as e.g. compressed sensing, phase retrieval, logistic regression, or in random artificial neural networks; we refer to Barbier et al. [2019] for a review of its numerous applications.
Here, we are interested in the statistics of the number of critical points, or more precisely the complexity of the associated empirical risk (2). For any open intervals B ⊆ R + and Q ⊆ (−1, 1), we consider the (random) number Crit n,L 2 (B, Q) of critical points of the function L 2 with loss value in B and overlap with the signal q ≡ x · x in Q: Here grad is the Riemannian gradient on S n−1 . For L 1 we define the similar quantity Crit n,L 1 (B), dropping the notion of overlap. The average E Crit n,L 2 (B, Q) is the quantity that can be analyzed rigorously. Its logarithm divided by n is called the annealed complexity. However, since the random variable Crit n,L 2 (B, Q) is in general strongly fluctuating and scales exponentially with n, its typical value is different from the mean and can be obtained by taking the exponential of E ln Crit n,L 2 (B, Q). This last quantity (divided by n) is called the quenched complexity. It is in general different from the annealed one, with very few exceptions [Subag, 2017, Crisanti andSommers, 1995].
Our main results consist in explicit formulas for the annealed and quenched complexities for L 1 and L 2 . The formula for the annealed case is obtained by a rigorous Kac-Rice method, whereas the one for the quenched complexity is obtained by theoretical physics methods combining the Kac-Rice method with replica theory. We consider the limit n, m → ∞ with m/n → α > 1, a setting called in the statistical physics literature the thermodynamic limit. The condition α > 1 is essential, as can be seen for instance in eq. (1) : if m < n, for each realization of {ξ µ }, the function L 1 has an infinite number of critical points in the set of unit-norm x orthogonal to all the {ξ µ }, and counting the critical points in this case is meaningless (or one has to quotient the space to take care of the degeneracy). Our results hold for many classical activation functions φ, such as e.g. the hyperbolic tangent, the arctangent, the sigmoid, or a smoothed and leaky version of the ReLU activation function 1 . Henceforth we shall denote C + the strict upper-half complex plane, and M(R k ) the set of probability measures on R k . For two probability measures µ and ν we define the relative entropy H(µ|ν) ≡ ln(dµ/dν)dµ if µ is absolutely continuous with respect to ν, and +∞ otherwise. Finally, µ G is a generic notation for the standard Gaussian measure on any R k .
We can now present our main results.
Theorem 1 (The annealed complexity of L 1 ). Let B ⊆ R a non-empty open interval and denote M φ (B) the set of probability measures ν on R such that ν(dt)φ(t) ∈ B. Given: • Let z ∈ R n×m an i.i.d. standard Gaussian matrix, and y ∈ R m a vector with components taken i.i.d. from a probability measure ν. Let D (ν) the diagonal matrix of size m with elements D (ν) µ = φ (y µ ). We define µ α,φ [ν] as the asymptotic spectral measure of zD (ν) z /m.
The precise hypotheses on the activation function φ are precised in Section 2. one has 1 : A note on free probability Interestingly, the measure µ α,φ [ν] can be interpreted as the free multiplicative convolution of the Marchenko-Pastur law (at ratio α) and the asymptotic spectral distribution of the matrix D (ν) , cf. e.g. Voiculescu [1987], Anderson et al.
[2010] 2 . We describe in Section 3 how to explicitly compute the density of µ α,φ [ν], or its linear spectral statistics (as e.g. κ α,φ (ν, C)), via the computation of its Stieljtes transform. We turn to our second annealed result: Theorem 2 (The annealed complexity of L 2 ). Let B ⊆ R + and Q ⊆ (−1, 1) two non-empty open intervals. For q ∈ (−1, 1) we denote M φ (B, q) the set of probability measures ν on R 2 such that Given: • Let z ∈ R n×m an i.i.d. standard Gaussian matrix, and Y ∈ R m×2 with components taken i.i.d. from a probability measure ν on R 2 . Let D (ν,q) the diagonal matrix of size m with elements D (ν,q) µ = f q (Y µ ). We define µ α,φ [q, ν] as the asymptotic spectral measure of zD (ν,q) z /m.
then one has: The proof of Theorem 1 is presented in Section 2. The proof of Theorem 2 is a straightforward generalization, and sketched in Appendix C.1. The variational problems in Theorems 1 and 2 are challenging, as they imply an optimization on a set of measures, and they involve transforms of this measure that are very hard to access numerically. In Section 3 we present a drastic simplification: a heuristic calculation that allows one to reduce the supremum over the probability measure ν to a much more straightforward optimization over a small number of parameters.
As we have already stressed, the annealed complexity, although interesting in itself, is generically not representative of the landscape corresponding to a given typical instance of the empirical risk. In 1 A fully rigorous statement would imply a lower and an upper bound given by a supremum over the adherence and the interior of M φ (B). For reasons of lightness and clarity of the presentation we write it in the simpler presented form.
2 Free multiplication is usually defined for positively-supported measures, however one can generalize it here by explicitly separating the positive and negative parts of φ (we can show freeness of the resulting two random matrices). order to obtain the value of the quenched complexity we use the replicated Kac-Rice method, which is an extension to non-Gaussian functions of the one developed in Ros et al. [2019]. Although the replica method is non-rigorous, it is considered an exact method in theoretical physics and it has been proven to give correct results for spin-glasses and inference problems [Talagrand, 2006, Barbier et al., 2019. We have obtained an explicit formula 1 for the quenched complexity of L 1 and L 2 at fixed values of the empirical risk, and overlap with the solution (in the L 2 case).
For L 1 , using the notations of Theorem 1 we have: Dξ ln I(ξ) .
Here, ξ = (ξ q , ξ a , ξ c , ξ c ) and Dξ is the standard Gaussian probability measure on R 4 . The extr denotes extremization with respect to all variables (A,Â, a,â, C,Ĉ,ĉ, {g(λ)}). We denoted In this formula, the notation extr denotes that one should set the partial derivatives with respect to the involved variables to zero. This notation arises from the replica calculation, which mixes saddle-point computations with Lagrange multipliers associated to certain constraints, and the precise meaning of this extremization (as a supremum or infimum) would have to be clarified by a more rigorous method. On a numerical point of view, one would have to solve the associated saddle-point equations, so that this precise meaning is not crucial for applications. We can state a very similar result for L 2 : Result 2 (Quenched complexity of L 2 ). Let B ⊆ R, Q ⊆ (−1, 1) two open intervals and: • For m ∈ (−1, 1), M φ (B, m) is the space of probability measures ν on R 2 that satisfy: Let z ∈ R n×m an i.i.d. standard Gaussian matrix, and Y ∈ R m×2 with components taken i.i.d. from a probability measure ν on R 2 . Let D (ν) the diagonal matrix of size m with elements D as the asymptotic spectral measure of zD (ν) z /m.
1 Here we used a replica symmetric structure, which is correct in many cases, and a very good approximation in others were replica symmetry has to be broken [Mézard et al., 1987].

Conclusion, outcomes and perspectives
We have obtained analytical results for the annealed and quenched complexities of statistical models with non-Gaussian loss functions arising in generalized linear estimation and simple models of glasses and neural networks. Our method is versatile and can be easily extended to other cases. We describe in Appendix C.2 three other inference models to which it applies: binary linear classification, a mixture of two Gaussians, and a simple model of unsupervised learning.
As a sanity check of our results, we have analytically verified by explicit solution that for a linear activation function, the annealed complexities of L 1 is null. It is again a tedious but straightforward computation to check that the annealed complexity of L 1 with a quadratic activation φ(x) = x 2 is also null, as the number of critical points in this case is linear with n. Note that for L 2 , even the case of a linear activation is non trivial, as shown in the very recent analysis of Fyodorov and Tublin [2020].
Our results allow for a complete characterization of the empirical loss landscapes of generalized linear models. The main issue ahead is determining for which class of functions φ and in which regimes (e.g. values of α), the annealed and quenched complexities become positive, i.e. when the associated landscape is rough. This will allow to study the connection between landscape properties and dynamics induced by local algorithms. In particular, it will shed light on the relationship between the roughness of the empirical loss landscape and the existence of "hard" phases in the learning of generalized linear models [Barbier et al., 2019]. It will also provide an interesting benchmark for obtaining the algorithmic thresholds of gradient descent (and variants) only through the knowledge of the landscape properties [Sarao Mannelli et al., 2019a. Based on ongoing works, we can for instance conjecture the existence of a rough landscape for small enough α in phase retrieval [Lucibello et al., 2019] and retarded learning [Engel and Van den Broeck, 2001].
Addressing these questions requires additional work, which is beyond the scope of this paper. We have explained in Section 3 how to make tractable the variational problem associated to the annealed and quenched complexity. Its analysis for specific models is an ongoing direction of research and will be presented elsewhere.
Another important extension of our results consists in counting the critical points of a fixed index (i.e. with a fixed number of negative directions in the spectrum of the Hessian). This would provide additional interesting information, in particular it would allow to differentiate local minima from the other critical points of the landscape, as was done for spin glass models in Auffinger et al. [2013]. Such a counting would require to understand the large deviations properties of the eigenvalues of the Hessian arising for generalized linear models, a random matrix problem that is hard but hopefully tractable by building on recent developments [Maida, 2007, Ben Arous et al., 2019. This is an ongoing work that we are pursuing, with already promising results.
As final note, it is an open problem to generalize our methods to neural network models with many nodes and hidden layers; the random matrix analysis of the Hessian in this case is a particularly exciting challenge.

Proof of Theorem 1 for the annealed complexity
In this section we prove Theorem 1. The technique leverages the Kac-Rice formula and Sanov's theorem on the large deviations of the empirical measure of i.i.d. variables. First we precise our hypotheses on φ, that we will take in the following set of "well-behaved" activation functions: Definition 1. φ : R → R is "well-behaved" if it is of class C 3 and if, for y ∼ N (0, 1), the random variable a = φ (y) admits a continuous probability density in a neighborhood of a = 0.

The Kac-Rice formula
The first step is to apply the Kac-Rice formula to the random function L 1 : Lemma 3 (Kac-Rice formula). For any x ∈ S n−1 , denote gradL 1 (x) and HessL 1 (x) the (Riemannian) gradient and Hessian of L 1 at the point x. Then gradL 1 (x) has a well defined density (on the tangent space T x S n−1 R n−1 ) in a neighborhood of zero, that we denote ϕ grad L 1 (x) . Denote µ S the usual surface measure on S n−1 . One has: The proof of this lemma uses necessary conditions for a random function to be a.s. Morse 1 , that are stated in Azaïs and Wschebor [2009]. The details are given in Appendix A.1.

The complexity at finite n
In this section we state the result of the Kac-Rice method. For y ∈ R m , let Λ(y) ∈ R m×m : in which we (abusively) denote φ (y) ≡ (φ (y µ )) m µ=1 , and D(y) ∈ R m×m the diagonal matrix with elements D(y) µ = D(y µ ) = n m φ (y µ ). The main result of this section is: Lemma 4 (Complexity at finite n).
in which C n is exponentially trivial, meaning lim n→∞ (1/n) ln C n = 0. The variable y ∈ R m follows N (0, I m ), and z ∈ R (n−1)×m has i.i.d. standard Gaussian matrix elements, independent of y. H Λ n (y) is a square matrix of size (n − 1) with the following distribution : The rest of the section is devoted to the proof of Lemma 4. We start from the result of Lemma 3. The following proposition specifies the joint distribution of (L 1 (x), gradL 1 (x), HessL 1 (x)): Proposition 5 (Distribution of the gradient and Hessian). Let x ∈ S n−1 . Then (L 1 (x), grad L 1 (x), Hess L 1 (x)) follows the following joint distribution: ∼ N (0, I n−1 ), and all {y µ , z ν } are independent. We identified in these equations the tangent spaces T x S n−1 with R n−1 .
Proof. Denote P ⊥ x the orthogonal projection on {x} ⊥ . For a smooth function f : S n−1 → R, ∇f and ∇ 2 f are its Euclidean gradient and Hessian. The Riemannian structure on S n−1 induces the gradient and Hessian of f as grad f (x) = P ⊥ x ∇f and Hess f ( Applying these formulas yields: Letting y µ ≡ ξ µ · x and z µ ≡ P ⊥ x ξ µ (identified to an element of R n−1 ) yields the result.
The joint distribution of eq. (10) is invariant with respect to x, thus we can chose x to be the North pole x = e n = (δ i,n ) n i=1 . With ω n ≡ 2π n/2 /Γ(n/2) the volume of S n−1 , we obtain from Lemma 3: Removing the e n indication and conditioning on the distribution of y in eq. (10), we reach: Once conditioned on y, eq. (10b) describes a Gaussian density so we can directly compute: in which ln C n = O n (n) (using Stirling's formula). The conditioning of the Hessian by grad L 1 = 0 at fixed y reduces to a linear conditioning on z. One thus obtains by classical Gaussian conditioning: in which H Λ n (y) is defined by eq. (9). This ends the proof of Lemma 4.

Large deviations
This section is devoted to the end of the proof of Theorem 1. Denote ν m y ≡ 1 m m µ=1 δ yµ the empirical distribution of y. We take the notations of Theorem 1 and Lemma 4. We first state an important lemma on the concentration of E z | det H Λ n (y)| 1 : Lemma 6. There exists η > 0 such that for all t > 0: The proof of Lemma 6 is detailed in Appendix A.2. Note that we expect this result to actually be valid up to η = 1, as the large deviations of the spectral distribution of random matrices is typically on the n 2 scale [Ben Arous andGuionnet, 1997, Hiai andPetz, 1998]. The following moment condition, proven in Appendix A.3, will be important: Lemma 7. For every γ ∈ (1, α) we have: Then we can conclude using Sanov's large deviation principle [Sanov, 1958, Dembo andZeitouni, 1998], and Varadhan's lemma. Using our previous lemmas, we obtain the statement of Theorem 1: The proof of eq. (18) is detailed in Appendix A.4.
3 Towards a numerical solution to the variational problem is given by the unique solution in C + to the implicit equation (as shown for instance in Silverstein and Bai [1995]): For any µ ∈ M(R) and t ∈ R we define the logarithmic potential as and an arbitrary t ∈ R, see Theorem 1. For clarity, we will write µ for µ α,φ [ν] for the remainder of this section. Let us define for any z ∈ C + , is well defined and holomorphic on C + . Moreover, from the Chapter II of Faraut [2014], we know that From this it is clear that in order to get the logarithmic potential in the bulk, we need to be able to evaluate G(z) for z ∈ C + . Define, for z, g ∈ C + 2 : At any fixed z, F (z, g) is an holomorphic function of g on C + . Its Wirtinger derivative is: Thus g(z) (the Stieltjes transform of µ, cf. eq. (19)) is the only g ∈ C + such that ∂F ∂g (z, g) = 0. Moreover, by definition g(z) is an holomorphic function on C + with values in C + . We can thus apply the usual composition of derivatives and obtain: Furthermore we know dG dz = −g(z). Computing the remaining constant by investigating the limit Re[z] → ∞, we reach that G(z) = F (z, g(z)) for every z ∈ C + . We thus have the crucial relation: This allows for an efficient numerical derivation of the logarithmic potential of µ α,φ [ν], as we will see in more details below.

Heuristic derivation of the simplified fixed point equations corresponding to Theorem 1
We present here an heuristic derivation of scalar fixed point equations for the numerical resolution of Theorem 1. This technique could be easily extended to Theorem 2 as well as the quenched calculations presented afterwards, but we restrict to this simpler case for the sake of the presentation.

Expressing
From eq. (23) we know that for every t ∈ R: For every t ∈ R, g(t + i ) is the only solution in C + to the partial derivative of the previous equation: So heuristically, we can write that for a small enough : with g = g r + ig i (in practice one iterates over g r and g i successively).

Heuristic solution to Theorem 1
We start from the result of Theorem 1. For a function f , we write E[f (X)] ≡ ν(dt)f (t). We introduce Lagrange multipliers to fix the conditions E[φ(X)] ∈ B, and we fix the values of E[φ (X) 2 ] and E[Xφ (X)]. We obtain: Note that now the supremum over ν is unconstrained over the set M(R) of probability distributions. We now make use of eq. (25) to write, with 2K(α) ≡ −1 + ln α − 2α ln α and a small > 0: For any scalar function F , the maximum sup ν [E[F (X)] − αH(ν|µ G )] is attained in ν * with density proportional to e −x 2 /2+F (x)/α , which is called the Gibbs measure in statistical physics. This gives (D is the standard Gaussian measure on R): Plugging this into our previous equation for the annealed complexity yields: This can be further simplified, as the extrema over A, t are trivially solved and give the value of λ 2 = Re[g] and λ 1 = (2A) −1 . Thus we obtain: Let us now denote then the fixed point equations of eq. (30) can be written as: These equations are to be iterated over λ 0 , λ 1 , g, and l (while enforcing the constraint l ∈ B). From experience, the best procedure is to start from the solution of the unconstrained problem (without any constraint on the loss value), before smoothly following the solution while adding the constraint. In the case of L 2 (x), one would follow a similar procedure.

The quenched complexity and the replica method
In this section we detail the principle of the quenched calculation that gives rise to Results 1-2. For the sake of the presentation we restrict to Result 1, while Result 2 will be discussed in Appendix C.1. We therefore focus on the function L 1 of eq. (1). As the very basis of this calculation is non-rigorous we present this calculation in a fashion closer to theoretical physics standards, differently from Section 2 in which we present rigorous results on the annealed complexity. Some technicalities will be postponed to Appendix B.

The replica trick and the p-th moment
The replica method is a heuristic tool of theoretical physics that allows to compute the quenched values of observables in the thermodynamic (i.e. n → ∞) limit from the knowledge of their integer moments, under some assumptions. It is based on the non-rigorous identity (note that it involves an inversion of limits), for a strictly positive function f of a n-dimensional random vector x: Mézard et al. [1987] gives a comprehensive introduction to the replica method and its (many) Here, ϕ {grad L(x a )} p a=1 (0) represents the joint density of the p gradients, taken at 0. Note that the nonlinearity L 1 (x) only depends on the parameters y a µ ≡ ξ µ · x a , so we will often write L 1 (y) ≡ L 1 (x). Proceeding as in the annealed case, we can rewrite the expectations by conditioning over {y a } p a=1 : The gradient and Hessian at x a live in the tangent plane to the sphere at x a , identified with R n−1 . Note that the {y a µ } are Gaussian variables with zero mean and covariance E[y a µ y b ν ] = δ µν q ab , with q ab ≡ x a · x b the "overlap" between replicas a and b. We introduce the variables {q ab } via delta functions in eq. (34): Since we fixed the {q ab }, the distribution of the {y a } is fixed, as well as the joint distribution of the loss, gradients and Hessians, as we will explicit in the following. As the number of overlap variables is p(p − 1)/2 = O n (1), we will perform a saddle-point over the variables {q ab } in the thermodynamic limit. The replica-symmetric assumption (see Mézard et al. [1987]) is a crucial hypothesis that can be made in the framework of the replica method. It amounts to assume that, once the saddle-point is performed, the extremizing {q ab } are "symmetric" over the different replicas of the system. Concretely, we assume that the variables {q ab } satisfy q aa = 1, q ab = q for a = b. Assuming this structure of the overlap matrix allows to extend the expression of the moments to arbitrary non-integer p, and then to take the p → 0 + limit as needed in eq. (32). We used a replica symmetric structure, which is correct in many cases, and a very good approximation in others were replica symmetry has to be broken.

The phase volume factor
Let us first compute the phase space factor in eq. (35). More precisely, the term: in which we denoted q aa = 1. As we detail in Appendix B.1 we reach, when p → 0 + and n → ∞:

The joint density of the gradients
We will now compute the joint density of the gradients at {x a }, conditioned on the values of {y a }. The calculation is an extension of Sections V.C and V.E of Ros et al. [2019]. We consider two vectors x a and x b of overlap q ab = q. It is easy to see that E[grad L(x a )|{y b } p b=1 ] = 0 from eq. (10b), so we will focus on the covariance matrix E[grad L(x a )grad L(x b ) |{y c } p c=1 ]. After some calculations detailed in Appendix B.2 we get the gradient density at leading exponential order: in which the auxiliary functions (z p (q), f p (q), f 0 p (q)) are defined in eq. (86).

Factorization of the mean product of determinants
The argument of this section is very close to Section V.F of Ros et al. [2019]. We consider the term: We make two important remarks, which are straightforward transpositions of the arguments of Ros et al.
[2019] to our problem, and we refer to this work for more extensive physical justifications.
• The conditioning over the gradients being zero, as in the annealed calculation, only gives a finiterank change to the Hessians Hess L 1 (x a ) and thus does not modify the limit at the scale e Θ(n) . At this scale, the statistics of the p matrices {Hess L 1 (x a )} p a=1 are identical.
• The spectral measure of Hess L 1 (x a ) concentrates at a rate at least n 1+ for a small enough > 0 (we expect that the actual rate is n 2 ), so that at the order e Θ(n) the expectation value factorizes and we can assume all the Hessians to be independent.
Before stating the consequences of such remarks, we give some definitions: • µ G,q is the Gaussian probability measure on R p with zero mean and covariance E[X a X b ] = (1 − q)δ ab + q. Note that {y µ } m µ=1 are i.i.d. variables distributed according to µ G,q .
• We define ν y as the empirical measure of (y 1 , · · · , y m ), that is ν y ≡ 1 m µ δ y µ . For every a, we denote ν a y its marginal distribution: ν a y (dλ a ) ≡ b( =a) ν y (dλ). Then ν a y is also the empirical distribution of (y a µ ) m µ=1 . Our remarks show that we can use the results of the annealed calculation, and we have here by factorization of the expectation of the determinants, at leading exponential order: 4.2 Decoupling the replicas and the p → 0 + limit We can then apply Sanov's theorem to the empirical measure ν y ∈ M(R p ). Recall that we have constraints on this measure by the density of the gradient and the fixation of the energy level. More precisely, we denote M (p) φ (q, B) the set of probability measures on R p that satisfy the following: Recall that the functions (z p (q), f 0 p (q), f p (q)) are defined in eq. (86). Leveraging from the results of eqs. (36), (37) and (39), we obtain from Sanov's theorem and Varadhan's lemma: Recall that ν a is the marginal distribution of ν for the variable λ a . We can then decouple the replicas under an assumption on the measure ν that amounts for replica symmetry. We stress that this replica symmetric assumption in the Kac-Rice calculation actually corresponds to a 1-step replica symmetry breaking (1RSB) structure of the zero-temperature Gibbs measure, that is an exponential number of single-point metastable states that all have the same two-point overlap. While possibly not exact, this assumption should already yield a good approximation to the landscape, and could be analytically checked by studying the stability of the replica-symmetric ansatz within replica theory. This allows to take subsequently the p → 0 + limit, and after some simplifications, we reach from eq. (41)  A Technical steps of the proof of Theorem 1

A.1 Proof of Lemma 3
We will apply the Kac-Rice machinery in the form of the remark made in Paragraph 6.1.4 of Azäis and Wschebor [2008]. We recall it as a theorem: Theorem 8 (Azais-Wschebor). Let k, d ∈ N . Let Z : U → R d be a random field, in which U is an open subset of R d . Assume that for every t ∈ U , we can write Z(t) = H [Y (t)], such that: (i) {Y (t), t ∈ U } is a Gaussian random field with values in R k , C 1 paths, and such that for every t ∈ U , the distribution of Y (t) is non-degenerate.
(ii) H : R k → R d is a C 1 function.
(iii) For all t ∈ U , Z(t) has a density ϕ Z(t) (x), which is a continuous function of (t, x) ∈ U × R d .
Define, for every compact set B ⊆ U , N (Z, B) to be the (finite) number of zeros of Z in B. Then: We wish to apply this theorem to the gradient gradL 1 (x). Verifying its hypotheses will end the proof of Lemma 3. We denote ξ ∈ R n×m the matrix {ξ iµ } = {(ξ µ ) i }, ∇L 1 the Euclidean gradient of L 1 , and P ⊥ x the orthogonal projection on T x S n−1 . Since gradL 1 (x) = P ⊥ x ∇L 1 (x) we have: We will apply Theorem 8 with d = n − 1 and k = m × n. The Gaussian random field Y (x) ∈ R n×m is defined as Y (x) ≡ P ⊥ x ξ 1 · · · P ⊥ x ξ m ξ 1 · x · · · ξ m · x . Since Y (x) is just ξ written in an orthonormal basis of R n whose first vector is x, its distributions is non-degenerate. H : R n×m → R n−1 is defined as: Since φ is C 2 , H is C 1 . This verifies (i) and (ii). We turn our attention to verifying (iii). One can write the distribution of the gradient of eq. (43) as grad L 1 (x) ∼ N (0, I n−1 ), and all {y µ , z ν } are independent. Since the distribution of grad L 1 (x) does not depend on x, it is enough to check that its density exists and is a continuous function. To do so, we will show that its characteristic functionφ gradL 1 (x) ∈ L 1 (R n−1 ). We denoteφ a the characteristic function of the random variable a ≡ φ (y), and one obtains: Since α > 1, if [qE zφa (qz)] = O q→∞ (1) we can conclude that φ gradL 1 (x) 1 < ∞. And: 2q 2 e iaz = 1 q E e − q 2 a 2 2 , by Fubini's theorem. Therefore q E zφa (qz) → q→∞ ϕ a (0) by continuity of ϕ a around a = 0 (Definition 1), so φ gradL 1 (x) 1 < ∞. Thus, gradL 1 (x) admits the following probability density: which is a continuous function of u, sinceφ gradL 1 (x) ∈ L 1 (R n−1 ). This shows (iii). In order to show (iv), we will use Proposition 6.5 of Azäis and Wschebor [2008], that we recall here: Lemma 9 (Azais-Wschebor). Let d ∈ N * , and U a compact subset of R d . Consider Z : U → R d a random field, such that (a): The paths of Z are of class C 2 , and (b): There exists C > 0 such that for all t ∈ U and all u in a neighborhood of 0, the density ϕ Z(t) of Z verifies ϕ Z(t) (u) ≤ C. Then P [∃t ∈ U s.t. Z(t) = 0 and det Z (t) = 0] = 0.
Since φ is assumed to be of class C 3 , hypothesis (a) is verified for Z = grad L 1 . Notice then that we can fix C > 0 such that E z∼N (0,1)φa (qz) ≤ C 1+q for all q ≥ 0. Starting from eq. (45): with C c , D n constants depending only on n, using that m ≥ n (α > 1). This shows (b), so by Lemma 9, hypothesis (iv) of Theorem 8 follows. This ends the proof of Lemma 3.

A.2 Proof of Lemma 6
The proof is done in several parts, of which some are inspired by arguments of Silverstein [1995], Silverstein and Bai [1995], Silverstein and Choi [1995], Bai and Silverstein [2010].

A.2.1 Technicalities on the Hessian
We begin by a quick lemma on Λ(y), defined in eq. (8).
Lemma 10 (Low-rank perturbation). Since the distributions of z and y are independent, by rotation invariance we can assume that Λ(y) is a diagonal matrix with elements Λ µ (y). There exists a constant, denoted ||D|| ∞ , such that for all n, y, |D(y)| ≤ ||D|| ∞ . Then we have: (ii) Let Z ∈ R (n−1)×m be i.i.d. variables with zero mean and unit variance. We denote µ Proof. Recall that |D(y)| = (n/m)|φ (y)|. Since m/n → α > 1 and φ is bounded, |D(y)| is bounded (uniformly over n, y) by a constant that we denote ||D|| ∞ . Note that sup 1≤µ≤m |Λ µ (y)| = sup ||u||=1 u Λ(y)u. Using eq. (8) and denoting v(y) ≡ φ (y)/||φ (y)||, we reach in which we used the uniform boundedness of |D(y µ )|, and the Cauchy-Schwarz inequality. This proves (i). We note that 1 n zΛ(y)z and 1 n zD(y)z differ by a rank-2 matrix. (ii) is thus an immediate application of the following lemma (from a course of C. Bordenave): Lemma 11. Let n ≥ 1, and A, B two symmetric matrices of size n, such that the rank of A − B is r. Denote F A (resp. F B ) the c.d.f. of the empirical spectral distribution of A (resp. B). Then This ends the proof of Lemma 10.
Lemma 11. We note λ 1 (A) ≥ · · · ≥ λ n (A) the eigenvalues of A (and similarly for B). Recall weak Weyl's interlacing inequalities [Weyl, 1912]: for every 1 ≤ i ≤ n, λ i+r (A) ≤ λ i (B) ≤ λ i−r (A) (we use the convention λ 1−i = +∞ and λ n+i = −∞ for i ≥ 1). Let t ∈ R, and i, j be the smallest indices such that λ i (A) ≤ t and λ j (B) < t. By the interlacing inequalities, |i − j| ≤ r. And We have some control of the boundedness of the Hessian, summarized in two subsequent lemmas: Lemma 12 (Moment bound). For all γ > 0, one has: Proof. We begin by bounding the extremal eigenvalues λ min , λ max of H Λ n (y) (denoted H for lightness): We used (i) of Lemma 10. Note that this bound is independent of y. In the same way we can bound λ min , and we reach: Using this identity, we have the bound: It is then a very classical result of random matrix theory that the largest eigenvalue of a Wishart matrix zz /n satisfies a large deviation principle in the scale n. This is stated e.g. in Theorem 2.4 of Biroli and Guionnet [2020],which gives moreover the behavior of the rate function I(x). We recall some of its properties: • I(x) : [s max , +∞) → R + is continuous and increasing.
Combining this inequality with eq. (46) ends the proof of Lemma 12.
(ii) The support of µ α,φ [ν m y ] is included in (−C, C) uniformly over y and n.
(iii) For all y ∈ R m , µ α,φ [ν m y ] has a well-defined and continuous density outside x = 0.
Points (ii) and (iii) of Lemma 13 are consequences of Theorem 1.1 of Silverstein and Choi [1995], while item (i) follows from the boundedness of Λ(y) by Lemma 10, and the one of xφ (x).

A.2.2 The cut-off and the logarithmic potential
For any > 0, define ln : x ∈ R + → ln (max(x, )), then x → ln |x| is a −1 -Lipschitz function on R. Let δ ∈ (0, 1). In this section, we show that a cut-off n = n −δ on the smallest eigenvalues does not perturb the logarithmic potential at the thermodynamical scale. We rely on the following result, proven in an ongoing work of Ben Arous, Bourgade and McKenna. For any δ ∈ (0, 1), there exists η > 0 such that for all t > 0: Physically, this makes explicit that, with large probability, there should not be enough eigenvalues of H Λ n (y) around zero so that they contribute macroscopically to the logarithmic potential. This is a consequence the natural fluctuations and repulsion of the eigenvalues of H N n (y). Denote {λ i } n−1 i=1 the (sorted) eigenvalues of H Λ n (y). We can now state: Lemma 14. There exists η > 0 such that for all K > 0: Proof. We consider η given by eq. (47a). Let t > 0. We denote A (n) We have for all y and t > 0 (Ā (n) t being the complementary event to A (n) t ): So that (using ln n (x) ≥ ln(x) for all x > 0): We know ln n |x| ≥ −δ ln(n). By Lemma 12, for all γ > 0: lim sup n→∞ sup y∈R m 1 n ln E z e γ n−1 i=1 ln n |λ i | < +∞.
Fixing γ > 1 and using Hölder's inequality, there exists therefore C > 0 such that for all K > 0 and t ∈ (0, K): in which we used the Markov inequality in (a) and eq. (47a) in (b).
In the following, for simplicity we will often abusively denote ln n | det H Λ n (y)| ≡ e n−1 i=1 ln n |λ i | and ln n E| det H Λ n (y)| ≡ ln E e n i=1 ln n |λ i | .

A.2.3 Concentration of the logarithmic potential with a cut-off
We show here that discarding the eigenvalues of the Hessian that are close to 0 using a cut-off n ≡ n −δ , we have concentration of the logarithmic potential.

Let us recall the Hoffman-Wielandt inequality Hoffman and Wielandt [2003]:
Lemma 16 (Hoffman-Wielandt inequality for the L 2 norm). Let k ∈ N , and A, B ∈ S k (R) be two symmetric matrices with respective eigenvalues λ 1 (A) ≤ · · · ≤ λ k (A) and λ 1 (B) ≤ · · · ≤ λ k (B). Then In particular if A and B are positive matrices one has Tr[AB] ≤ i λ i (A)λ i (B). We use this in eq. (50) along with the n δ -Lipschizity of f n : in which we used Lemma 10. We denote A the event It is a classical concentration result (cf e.g. Chapter 3.1 of Vershynin [2018]) that there exists c > 0 such that: In particular, this implies Let us now show that it suffices to prove the bound of Proposition 15 assuming that A does not occur. Indeed, n −2(1−δ) ln P z [A] ≤ −cn 2δ for a constant c > 0, and: We used Markov's inequality in (a). The inequality (b) can be obtained by very similar arguments than the one used to prove Lemma 12: by rescaling z by z , it is easy to see that the event A will not change the scaling of the large deviations of the largest eigenvalue of the Wishart matrix zz /n, so that the bound of Lemma 12 will also apply when conditioning by the event A. All in all, Using the law of total expectation, we can safely ignore the occurrence of A when showing Proposition 15. Assuming that A is not occurring yields: Recall the Lipschitz concentration of independent variables with laws satisfying the logarithmic Sobolev inequality with a uniform constant c (see for instance Anderson et al. [2010] for a proof and an introduction to the logarithmic Sobolev inequalities): Lemma 17 (Herbst). Let n ∈ N and P be a probability distribution on R n satisfying the Logarithmic Sobolev Inequality with constant c > 0. Let G be a Lipschitz function on R n with Lipschitz constant G L . Then for all t > 0, P[|G − EG| ≥ t] ≤ 2 exp[−t 2 /(2c G 2 L )]. It is easy to check that the Gaussian standard law of z, conditioned by the (extremely probable) event A, satisfies the Logarithmic Sobolev Inequality with constant c = 1 + O n (1). Applying Lemma 17 alongside eq. (54) finishes the proof.

A.2.4 The logarithmic potential of the asymptotic measure
In this part, we relate the expected logarithmic potential to the logarithmic potential of the measure µ α,φ [ν m y ], cf Theorem 1. Proposition 18 (Concentration on κ α,φ ). There exists η > 0 such that for all t > 0: Proof. The proof goes in two parts. First, we show that there exists η 1 > 0 such that 1 : We will then conclude by showing that there exists η 2 > 0 such that for all t > 0: We begin by eq. (57). We take η 2 given by eq. (47b). We have Since ln (x) ≥ ln(x) , we reach from this: and using eq. (47b), we reach eq. (57). Let us show eq. (56). Its proof is based on the following lemma, a consequence of the analysis of Silverstein and Bai [1995], Bai and Silverstein [2010]: Lemma 19. Denote g n (z) the Stieltjes transform of 1 n zΛ(y)z , and g α,φ [ν m y ](z) the one of µ α,φ [ν m y ], for z ∈ C + . Then there exists η ∈ (0, 1) such that for all z ∈ C + : The proof of this lemma is postponed to Appendix A.5. Let us fix η given by this lemma. As stated for instance in Theorem 2.4.4 of Anderson et al. [2010], a consequence of the Stieltjes inversion theorem is that for every Borel set E ⊆ R, we have in which µ n is the empirical spectral distribution of 1 n zΛ(y)z . Fix η 1 < η. We have, uniformly over y: Let us fix C > 0 given by item (i) of Lemma 13. Moreover, we can also bound t φ (ν m y ) by ||xφ || ∞ . This gives that for n large enough the quantity of eq. (60) is bounded (uniformly over y) by: Using eq. (59) and since η 1 < η, this shows eq. (56).

A.2.5 Conclusion of the proof
Let us conclude the proof of Lemma 6 from all our results of Section A.2. We fix δ > 0 such that δ < 1/2. Note that Proposition 15 is a uniform result on y, much stronger than what we required. Since it shows that the concentration (as a function of z) of the log-determinant is in a scale n 1+ with > 0, it implies that there exists η > 0 such that for all t > 0: Combining this identity with Lemma 14 and Proposition 18, we reach the conclusion of Lemma 6.
A.4 Proof of eq. (18) Let t > 0, and fix η > 0 given by Lemma 6. We define E A n is related to the complexity by Lemma 4, and by Varadhan's lemma (and Lemma 7), we have: The factor α in front of the relative entropy arises as we consider the empirical distribution of m = αn i.i.d. variables. For all t > 0, we have by definition of A n , B n : Using Hölder's inequality and Lemma 7, there exists γ > 1 and a constant C > 0 such that: • Assume that lim B n = −∞ and lim sup A n > −∞. Let us fix a lower-bounded sub-sequence A ϕ(n) of A n , so that lim[A ϕ(n) − B ϕ(n) ] = +∞. However, by eq. (67) and Lemma 6, we have n ] → −∞. So we showed that lim B n = −∞ ⇒ lim A n = −∞, which shows eq. (18) in this case.
• Let us now assume that lim B n > −∞. Using the left inequality of eq. (67) and Lemma 6, we reach in the same way that lim inf[A n − B n ] ≥ −t, which implies that lim inf A n > −∞. Thus we can use the right inequality of eq. (67) to show similarly that lim sup[A n − B n ] ≤ t. Taking the t → 0 limit finishes the proof of eq. (18).

B Technical aspects of the quenched calculation B.1 The phase volume factor
Introducing the Fourier transform of the delta functions, we reach at leading exponential order in n: The replica symmetric assumption can be made on the variablesq that achieve this supremum :q aa =q 0 andq ab = −q for a = b. This leads to detq = (q 0 +q) p−1 (q 0 − (p − 1)q), and after taking the p → 0 + limit, we reach: The diverging term − 1 2 log n will be canceled out by the joint density of the gradients as we will see later. The solution of the supremum is easy to carry out, and we finally reach eq. (36).

B.2 The joint density of the gradients
We denote S = Span ({x a } p a=1 ). Following Ros et al. [2019], for every 1 ≤ a ≤ p we can construct an orthonormal basis of S, denoted (e a b ) 1≤b≤p for which x a is the first vector, that is e a a = x a . This basis is convenient, since {x a } ⊥ ∩ S = Span ({e a b } b( =a) ). We can also chose an arbitrary orthonormal basis (e p+1 , · · · , e n ) of S ⊥ . With this choice of basis, we can see that the gradient grad L(x a ) is identified with the vector in R n−1 with components: Recall that ∇L(x a ) = 1 m µ ξ µ φ (y a µ ). Let us make a few remarks: • For every a, the basis (e a b ) p b=1 is only a function of the values of the overlaps {q ab }.
• We consider the joint density of the gradients conditioned by the value of {y a }. In particular, this means that for every a = b, ∇L(x a ) · e a b is fixed by the values of {y c } p c=1 and the overlaps q ab . In particular, the first (p − 1) components of eq. (84) are deterministic, thus their density will yield delta functions that are constraints on {y a } and {q ab }.
• The last n − p components of eq. (84) are (at fixed {y a }) zero mean Gaussian random variables with covariance given by E gradL(x a ) i gradL( . Their joint density taken at 0 is thus at leading exponential order in the n → ∞ limit: Given these remarks and eq. (85), in order to complete the calculation of the joint gradient density we need to compute the quantities (∇ L(x a ) · e a b ) for every a = b as a function of {y a µ } and {q ab }. In order to simplify the calculation, we will already make use of the replica-symmetric assumption on q, that is we assume q aa = 1 and q ab = q for a = b. Let us now describe a possible construction for the basis (e a b ) p b=1 . We introduce three auxiliary quantities that are functions of q and p: With these definitions, we can consider: It is straightforward to check from eq. (87) that we have for all a, b, c that e a b · e a c = δ bc . We can now see that the delta term of the joint density of the gradients taken at 0 is: The product of eq. (85) and eq. (88) gives eq. (37).
B.3 Decoupling the replicas and taking the p → 0 + limit B.3.1 Replica symmetry and decoupling In order to apply the replica method, we need to be able to take the p → 0 limit, by analytically extending eq. (41) to all p > 0. The main idea is that we expect replica symmetry to influence the measure ν that appears in eq. (41). More precisely, we expect that for all permutation π ∈ S p , we have ν(dλ 1 , · · · , dλ p ) = ν(dλ π(1) , · · · , dλ π(p) ). We make in eq. (41) the substitution:  (89) In this last expression, the replica symmetric assumption leads us in particular to assume that µ a = µ for all a. In order to make the remaining calculation tractable we will also need to fix some linear statistics of ν via Lagrange multipliers: • For every a ≤ b, we will fix the linear statistics ν(dλ) φ (λ a )φ (λ b ) = A ab , with Lagrange multipliersÂ ab . Note that by replica symmetry, we actually assume that A ab = a for a = b and A aa = A (and samely for the Lagrange multipliers).
• For all a, b we will fix the linear statistics ν(dλ) φ (λ a )λ b = B ab , with Lagrange multipliersB ab . By replica symmetry, we assume that B aa = B and B ab = b.
Combining these remarks, we reach that the ν-dependent term of eq. (41) is equal to: In the Kac-Rice formula, one has now to integrate over the overlap q ≡ x · x as well. Moreover, we condition over the joint values of a µ ≡ ξ µ · x and b µ ≡ (1 − q 2 ) −1/2 [(ξ µ · x ) − qa µ ], rather than just ξ µ · x (as we did for L 1 ). Note that (a µ , b µ ) follows a joint standard Gaussian distribution.Using these definitions, we can obtain the counterpart of Lemma 4 for L 2 : E Crit n,L2 (B, Q) = C n Q dq e n(1+ln α+ln(1−q 2 )) 2 E a,b δ(P n (a, b))1 L2(a,b)∈B e −nEn(a,b) E z |det H n (a, b)| , in which C n is exponentially trivial, and we defined: Here z ∈ R (n−2)×m is an i.i.d. standard Gaussian matrix. T The condition P q (a, b) = 0 arises from the conditioning on the nullity of the gradient in the linear subspace spanned by x , and E n (a, b) from the density of the gradient in the subspace orthogonal to {x, x }. A crucial feature of this equation is that the joint distribution of (L 2 (x), gradL 2 (x), HessL 2 (x)) only depends on x via the overlap q = x · x with the "true" solution. Once conditioned over the values of q, it thus becomes clear why the calculations made for L 1 will generalize here. As in Section 2.3, one can then show the concentration of the empirical logarithmic potential on the functional κ α,φ (q, ν m a,b ), in which ν m a,b ∈ M(R 2 ) is now the empirical measure of {a µ , b µ } m µ=1 . We obtain the counterpart of Lemma 6: there exists η > 0 such that for all t > 0: lim n→∞ 1 n 1+η ln P 1 n ln E z |det H n (a, b)| − κ α,φ (q, ν m a,b ) ≥ t = −∞.
Thanks to this result, we perform then a saddle point on the overlap q and the empirical measure ν ∈ M(R 2 ), using Sanov's theorem and Varadhan's lemma. This yields the result of Theorem 2. As a final note, there exists similar results to the one presented in Section 3 that allow to compute the density (and the logarithmic potential) of µ α,φ [q, ν], via the computation of its Stieltjes transform.

C.2 Generalizations to other models
Our calculations, both annealed and quenched, can be generalized straightforwardly to many other loss functions and models. As is clear for instance in the annealed computation of Section 2, the key features that must be present are a Gaussian distribution for the data, and a loss function L(x) that only depends on the data samples ξ µ via their projection over a few vectors (as x for L 1 (x) and x, x for L 2 (x)). We give thereafter three examples of models, that can be found in Engel and Van den Broeck [2001], Mei et al. [2018], and for which the calculations can be performed.
Model 2 (Mixture of two Gaussians). Consider n, m ≥ 1 such that m/n → α > 1. We are given m samples y µ ∈ R n , generated as y µ i.i.d.
∼ 2 a=1 p a N (θ 0 a , I n ). The proportions p 1 , p 2 are known, and we wish to recover θ 0 1 and θ 0 2 by minimizing the maximum-likelihood estimator: Model 3 (Simple unsupervised learning model). Consider n, m ≥ 1 such that m/n → α > 1. Let φ : R → R a smooth activation function, V : R → R + a "potential", and x 0 ∈ S n−1 a fixed vector. We assume that we are given i.i.d. data samples {ξ µ } m µ=1 ∈ R n distributed such that their projection on x 0 has a probability density P (ξ µ · x = h) ∝ e − 1 2 h 2 −V (h) , and the other coordinates of ξ µ are i.i.d. standard Gaussian variables. We wish to recover the vector x 0 by minimizing: For each of these three models one can replicate the annealed and quenched calculations of Sections 2 and 4, under suitable technical hypotheses.
A note on non-spherical priors Finally, it is clear from the calculation of Appendix B (particularly Section B.1) that we can also generalize these techniques (at least heuristically) to non-spherical prior distributions on the vectors x. The most natural hypothesis that allows the computation to be generalized is that the prior distribution takes the decoupled form P (dx) = i P (dx i ).


extfieldsofstudy:
Computer Science
Engineering
Mathematics
Environmental Science

text:
Machine learning approach to remove ion interference effect in agricultural nutrient solutions

High concentration agricultural facilities such as vertical farms or plant factories consider hydroponic techniques as optimal solutions. Although closed-system dramatically reduces water consumption and pollution issues, it has ion-ratio related problem. As the root absorbs individual ions with different rate, ion rate in a nutrient solution should be adjusted periodically. But traditional method only considers pH and electrical conductivity to adjust the nutrient solution, leading to ion imbalance and accumulation of excessive salts. To avoid those problems, some researchers have proposed ion-balancing methods which measure and control each ion concentration. However, those approaches do not overcome the innate limitations of ISEs, especially ion interference effect. An anion sensor is affected by other anions, and the error grows larger in higher concentration solution. A machine learning approach to modify ISE data distorted by ion interference effect is proposed in this paper. As measurement of TDS value is relatively robust than any other signals, we applied TDS as key parameter to build a readjustment function to remove the artifact. Once a readjustment model is established, application on ISE data can be done in real time. Readjusted data with proposed model showed about 91.6 ~ 98.3% accuracies. This method will enable the fields to apply recent methods in feasible status.

INTRODUCTION
Development of IoT applications enabled the propagation of smart farms around agricultural industries. Even small farm owners now adopts smart farm facilities to reduce human resource. Highly concentrated facilities such as vertical farm or plant factories are one of the greatest issue in smart farm area. Aqua cultural methods such as NFT(nutrient film technique), DWT(deep water technique), water dripping or aeroponics are now considered as optimal solutions for highconcentrated smart farms because it becomes significantly efficient to control nutritional environment of whole facility with hydroponic approaches.
Opened hydroponic system discards once-used nutrient solution because the nutritional ions are already uptaken by plant roots. This approach tries to supply stable and balanced nutrition to plants. However, the wastage of nutrition water is hazard to environment so closed hydroponic systems becomes popular [1]. Closed system recycles the nutrient solution so it significantly reduces the cost and possibility of pollution.
It is one of the most important issue in closed system to maintain nutritional balance. Traditional approaches measure electrical conductivity and pH of the nutrient solution to monitor the fertilization [2][3]. However, those methods cannot distinguish the imbalance of nutrient ions and accumulation of useless components, caused by the difference between absorption rate of individual ions at plant roots. It requires additional maintenance labors such as foliar spray to avoid nutrition disorder which may lead to poor yield or death of crops.
Researchers recently suggests application of the ion selective electrodes(ISEs) approaches to measure individual ion's concentration from the nutrient solution [1,[4][5][6]. Readjustment of low-concentrated ions with ion selective measurement seems to enable avoidance of nutrient balance collapse. However, applying ISEs on nutrient solution measurement have significant problem: ion interference effect.
ISEs acquires concentration of ion by measuring electrical potential. ISE has ion-selective permeable film and the ionexchange of specific ion fit for the film provides Nernst potential, described on equation (1). F is Faraday constant and z is the number of ions exchanged along the membrane.
However, because of the nature of ISE, ion-exchange along the membrane is interfered by other ions. This phenomenon is described with Nikolsky-Eisenman equation [7] given on equation (2) where is the activity, is the selectivity coefficient and is the exchange of interfering ion i along the membrane.
The activity parameter for target ion and the selectivity coefficients of all the other ions of solution are required to predict the magnitude of ion interference potential. It is not feasible. For example, Yamazaki's nutrient solution for lettuce [8] contains 23 different ions. Acquisition of those parameter with experiments to maintain agricultural fertilizer is not desirable. Liberti et al. suggested readjustment method using Gran's plots [9] but it also requires prior knowledge for some specific parameter.
We suggest a machine learning approach to train a readjustment function, which removes ion interference effect from ISE measured data without any prior knowledge chemicophysical parameter.

II. ION INTERFERENCE EFFECT
We measured ion interference effect with ISEs to figure out the exact magnitude of the artifact. Also the data was used for the machine learning algorithm.

A.
Equipments for experiment Ion selective electrodes were used to measure the concentration of ions. We used Vernier's Go Direct® sensor series. GDX-NO3 was used to measure NO3ion, and GDX-NH4 was used for NH4 + , GDX-CA was used for Ca 2+ , and GDX-K was used for K + ion concentration. The sensors were connected to Arduino to transform analog signal into digital value.
Yamazaki's nutrient solution for lettuce [8] was selected as experimental environment because it is one of the most popular fertilizer in hydroponics area. But we eliminated the low-concentrated ions in order to avoid ISE-unmeasurable ions. KNO3, Ca(NO3)2-4H2O and NH4H2PO4 were used to reproduce the simplified Yamazaki's solution. The experimental solution contains K + , NO3 -, Ca 2+ , NH4 + , which are measurable with ISEs, and HxPO4 3-x families which are the only unmeasurable value. The concentration of H + is measurable with pH meter but we did not measured it because pH of the solution is around 7 so we did not expect any significant interference effect on pH meter.

B.
Calibration with machine learning algorithm Vernier's manual recommended to measure 2 standard buffer solution to acquire 2 voltage value, and perform linear regression on exp(volt) term along the theoretical value of ion concentration. They presents voltage-concentration equation as equation (3), which fits for Nernst potential described on equation (1). Generally, α describes the sensitivity of sensor and β describes standard potential offset of the electrode.
However, we figured out that Vernier's ISE does not work as equation (3). 10 times of buffer solution measurement showed that the correlation between concentrations and exp(volt) values showed exponential behaviors. The distribution of measured buffer solution data is described on Figure 1.
It means that Vernier's ISE does not work robustly. The behavior of sensor in thin solution differs from the behavior in dense solution. The concentration of buffer solutions, measured voltage and coefficients for equation (4) of each ion is presented in Table 1. Coefficients are driven with linear regression between log scaled concentration and exp(volt) values. The R 2 for each regression lines are also listed on Table 1. Regression models with equation (4) fits for the experimental data.

C.
Ion interference effect measurement We prepared 100 times more concentrated version of simplified Yamazaki's solution. We added this ampoule 10 times to 1 liter water with 4 ISEs installed. The theoretical values of ionic concentration for each steps and the experimental values are presented on Table 2.
Comparison between the experimental values and the theoretical values are described on Figure 2. The interference effect becomes larger at higher-concentration status because the flux of undesirable ions through the membrane of ISE becomes greater magnitude of ion interference effect is understood as the distance between the blue and orange lines.

III. METHOD
A readjustment method to remove ion interference effect is proposed in this section. In order to provide a rapid and feasible solution to the industrial field, we avoided evaluation of chemicophysical parameters. We applied simple machine learning algorithm which doesn't even requires any GPU devices.

A.
Equipments The regression process was performed on Microsoft's Azure VM standard D4s v3, which has4 vcpus and 16 GiB memory. The operation system is ubuntu 18.04 and we used R and R Studio Server for data analysis and regression.

B.
Data preperation Experimental data on Table 2 from ion interference effect measurement was used. We randomly separated 3 columns for test data and the other 7 columns for training data. Column 5, 7 and 9 were selected as test data and the others were fed to the model for training.

C.
Readjustment method We applied simple readjustment method in order to application of the trained model on readjustment as simple as possible because real-time and on-site readjustment is one of the biggest issue in the agricultural industry [1].
As TDS(total dissolved solids) is calculated regardless of ion composition in the solution, we took TDS value as key independent variable. The readjustment method is described as equation (5), where Cr is readjusted value for concentration, CISE is experimentally driven value for concentration and μ is a function of TDS value.
We assumed the concentration of each ion and the magnitude of ion interference effect as functions of TDS in Yamazaki's nutrient solution system because the imbalance of ionic balance does not occur fast. Also, with continuous readjustment on the nutrient solution, the ionic balance is maintained continuously.  Regression of function TDS values were calculated from the theoretical value of ion concentrations from Table 2. TDS is calculated as equation (6). Mion is the molecular weight and [ion] is molarity.
And we expect the readjusted value Cr to converge toward theoretical value Ct. Therefore, we launched regression to fit equation (9).
The values of μ(TDS) driven from 7 points of training data are plotted on Figure 3 as blue dots.
As the behavior of voltage along ion concentration from equation (1) and the interference effect driven from equation (2) are both logarithm, we expected that the error correction coefficient μ to show logarithmic distribution. However, quadratic regression fit much better than log regression. Mean R 2 from log regression was 0.70615 while quadratic regression showed 0.96158. The regression lines are described as orange lines on Figure 3. Parameters for quadratic regression model equation (10) and R 2 values for the models are also presented in Table 3.
IV. EXPERIMENT

A. Evaluation on training data
To evaluate the trained model on training data, we calculated μ(TDS) from equation (10), rather than equation (9) because value from equation (9) will provide 100% correctness, which is meaningless. The driven μ are multiplied with CISE to show whether the ion interference effect were removed. The results are presented on Figure 4. The model removed 45.8752%, 63.2386%, 81.5443% and 92.434% of ion interference effect from K + , Ca 2+ , NO3and NH4 + . R 2 of Ca 2+ is greater than that of NO3but NO3showed better result.

B.
Evaluation on test data Evaluation on test data was performed by trained model. We apply untrained TDS value as input to μ function, and multiplied that to CISE of test data.

V. RESULTS
Test data was applied on the trained quadratic model. The result is described on Figure 5 and Table 3-4. The error was defined as equation (11), where C stands for both CISE and corrected concentration. Accuracy was defined as 100%error.
The errors on whole date are presented on Figure 3. Except only one case, column 5 from K + , proposed method successfully removed the ion interference effect. Especially, NH4 + model showed remarkable result. Average accuracy values are presented on Table 4. The model worked on both training data and test data. The accuracy between training data and test data from Ca + was significantly different. However, the proposed method successfully increased the accuracy of test data, 59.331%, into 93.209 %.

VI. CONCLUSION
Machine learning algorithm with low computational complexity to remove ion interference effect was presented. Ion interference effect exists and gets stronger in denser environment. Quadratic regression on equation (9) successfully removed ion interference effect from almost all the data point, except only one case. As the computational resource for using pre-trained quadratic model is small enough, the presented method is highly feasible to industries. Even an embedded device can run the inference process in milliseconds. We expect applying higher-variance machine learning algorithms such as deep neural network, rather than high-bias and low variance model such as quadratic regression, would enhance the ability of artifact removal.
ACKNOWLEDGMENT Special thanks to Janghun Lee for selecting ISE devices and establishing measurement pipeline with Arduino.  Table 4. Error rate of raw data and proposed method (in %) The column 5, 7 and 9 are training data and the other columns are used for training. Table 5. Accuracy of raw data and proposed method (in %)


Neural Machine Translation with Recurrent Highway Networks

Recurrent Neural Networks have lately gained a lot of popularity in language modelling tasks, especially in neural machine translation(NMT). Very recent NMT models are based on Encoder-Decoder, where a deep LSTM based encoder is used to project the source sentence to a fixed dimensional vector and then another deep LSTM decodes the target sentence from the vector. However there has been very little work on exploring architectures that have more than one layer in space(i.e. in each time step). This paper examines the effectiveness of the simple Recurrent Highway Networks(RHN) in NMT tasks. The model uses Recurrent Highway Neural Network in encoder and decoder, with attention .We also explore the reconstructor model to improve adequacy. We demonstrate the effectiveness of all three approaches on the IWSLT English-Vietnamese dataset. We see that RHN performs on par with LSTM based models and even better in some cases.We see that deep RHN models are easy to train compared to deep LSTM based models because of highway connections. The paper also investigates the effects of increasing recurrent depth in each time step.

Introduction
Neural Machine Translation(NMT) is a recent approach towards machine translation( [6], [2], [1] , [3]).In contradiction with the conventional Statistical Machine Translation(SMT) system ( [7]) which consists of many small sub-components that are tuned separately, in neural machine translation the whole neural network is jointly trained to maximize the conditional probability of a correct translation given a source sentence, using the bilingual corpus. The use of Neural Networks for Machine Translation leads to fluent translation.Popular neural machine translation models( [2], [1], [3]) use stacked LSTMs, which read through the source sentence one word at a time till it reaches the end of sequence symbol< eos > tag, then it starts to emit target words one by one till it generates the end of sequence < eos > tag. We move in the same direction, but instead of stacking layers of LSTM, we stack layers within each time step, as illustrated in Figure 1. Each time block is recurrent highway network following the work of [4] . With this architecture we demonstrate their effect on machine translation when layers within a time step is increased one by one.Parallelly with the neural machine translation(NMT), the concept of "attention" has gained a lot of popularity in recent years( [2], [3]). We follow the attention model by [3] as our base model. With IWSLT English-Vietnamese dataset we found that without attention increasing the recurrent depth does not affect the BLEU [8] score much, whereas with attention it shows a significant increase in the BLEU score with increasing recurrent depth. The problem with these vanilla attention models is that translation generated by these models often lack adequacy.There is a problem of over-translation and under-translation i.e some words are translated more than once ,while some words are not translated at all.To address this problem we use reconstructor model( [9]) which we will discuss in Section 4.

Related Work
Even prior to the recent end-to-end NMT models neural networks were used for SMT based translation with some success.The concept of end-to-end learning for machine translation has been attempted in the past with limited success. From then its a long journey to the end-to-end models( [1,2,3,6]) which have outperformed the previous SMT baselines( [7]). Our work being the first application of RHN in neural machine translation task, is very close to the work [1], [3] and [9]. At the same time we take the motivation to apply the RHN in this area from [5]; who has shown the success of recurrent highway network(RHN) in automatic speech recognition.

Recurrent Highway Network
Recurrent state transition in RNN is described by Similarly, a Recurrent Highway Network(RHN) layer has multiple highway layers in the recurrent state transition.RHN were first introduced in [4].
Let, W H,T,C ∈ R n×m and R H l ,T l ,C l ∈ R n×n represent the weights matrices of the H nonlinear transform and the T and C gates at layer l = 1, ..., L. The biases are denoted by b H l ,T l ,C l ∈ R n and let s l denote a intermediate output at layer l with s[t 0 ] = y [t −1] . Then an RHN layer with a recurrence depth of L is described by: where and It is to be observed that a RHN layer with L = 1 is essentially a basic variant of an LSTM layer. The computation graph for a RHN block is illustrated in Figure 2.

Model
Many sequential processing tasks require complex nonlinear transition functions from one step to the next. However, recurrent neural networks with "deep" transition functions remain difficult to train, even when using Long Short-Term Memory (LSTM) networks as shown by ([4], [11]).We use existing neural machine translation models with recurrent highway networks(which extend the LSTM architecture to allow step-to-step transition depths larger than one) as their recurrent block .

Model 1
We propose a neural machine translation model using the recurrent highway network architecture given by [4], and we use encoder-decoder system on this architecture with attention as given by [3] as our two base models.We use RHN instead of LSTM to model conditional probability of target sentence given source.

Model 2
The main problem with previous models and attention mechanisms was that sometimes translation produced by them lacked adequacy.There are two reasons for this: -Over-Translation-Some words are translated again and again.This is because there is no mechanism in the previous architecture to see that a word which is translated is not heavily attended.
-Under-Translation-Some words are not translated at all.This is because there is no mechanism to make sure that all words are attended or translated.
To tackle this problem we use model 1 with a reconstructor block as in [9] .The encoder-decoder framework is used to translate source sentence into target sentence.The decoder-reconstructor framework is used to convert target sentence into source sentence again.There is maximum information flow from encoder to decoder so that the reconstructor is able to reproduce source sentence from target sentence.By using reconstructor,model 2 tries to make sure that almost all source words are translated into target sentence so that reconstructor can reconstruct source words later.Thus ,translation produced is fluent and adequate.We use RHN instead of LSTM to model conditional probability of target sentence given source and to model conditional probability of source sentence given target hidden states (in reconstructor block). The thing to be noted is that we use the reconstructor block only during training phase.During inference,we use the regular encoder-decoder framework without reconstructor block.Unlike [9],we only Fig. 3. Example of an reconstructor-based NMT system as described in ( [9]).Green blocks are part of reconstructor. For clarity, the embedding and projection layers are not shown in Figure. use reconstructor in training phase and not during inference.The loss function used for training comprises of two terms,decoder loss and reconstructor loss. We train the neural network by minimizing two loss terms.The first term decoder loss L d ,is the log probability of a correct translation T given the source sentence S. So the training objective is The second term reconstructor loss L r ,is the log probability of a correct translation S given target hidden states H.
where S is the training set. The final loss function is given by: where β is a hyperparameter which reflects the importance of reconstruction compared to translation task and should be tuned based on the dataset used.

Experiments
We applied our method on the IWSLT English-Vietnamese dataset in both directions. Reference sentences are used to compute BLEU score of the translated sentences. We changed the number of hidden units per cell and recurrence depth(no. of layer in each time step) while keeping all other parameters constant. We report the performance of these translation models in Table 1 , and visualize the resulting translations. We report the performance of these translation models using LSTM in Table 2.

Dataset details
IWSLT English-Vietnamese data set is used for our experiment. The dataset has 133k sentence pairs of training sentences. We used 'tst2012'(en-vi) and 'tst2013'(en-vi) as our development and test dataset respectively. The dataset is available at Stanford NLP web page.Out of vocab words produced during translation, are replaced by an '< unk >' token.

Evaluation Metrics
-Perplexity :In natural language processing, perplexity is a way of evaluating language models. A language model is a probability distribution over entire sentences or texts.
-BLEU Score:BLEU [8], or the Bilingual Evaluation Understudy, is a score for comparing a candidate translation of text to one or more reference translations.

Training details
-We have used SGD with a learning rate of 0.1 .
-We have used a dropout of 20%.
-A batch size of 32 is used.
-Roughly all the experiments ran for 30 epochs.

Sample Translations
A few example translations produced by our model1 and model2 along with the reference sentences are shown in Table 3.From this randomly picked sample it is visible that, NMT2 avoids problems of over translation as seen in NMT1.Now we check the translation of the randomly picked sentences for recurrence depth 3 and 4 in Table 4. From this randomly picked sample it is visible that, increasing the recurrence depth helps to capture the semantics of the language better.  Depth-3 When I was a little bit , I thought , "Can see this country is the best country in the world and I often sing " what I can't envy".
Depth-4 When I was a little boy , I think < unk > is the best country in the world and I often sing "We doesn't have to envy".

Ref
Since my family couldn't understand Chinese , I thought my family was going to be arrested . Depth-3 Because my family didn't understand Mandarin , so I thought they would be arrested . Depth-4 Because my family didn't know Chinese ,so I thought they would be arrested .

Model analysis
To check how good our model is learning with the change of layers i.e. recurrence depth, we plot the training perplexity of the model 1 with different number of recurrence depth.
It can be seen that, there is a visible difference in the curve when we change depth from 2 to 3, and a bit less difference when we change from 3 to 4. Perplexity decreases with the increment of recurrence depth( Figure 4).As we can see in Table 5,RHNN models have less number of parameters as compared to LSTM and are able to perform as good as LSTM models or in some cases even better.

Conclusion
We use Recurrent Highway Networks, a powerful new model designed to take advantage of increased depth in the recurrent transition while retaining the ease of training of LSTMs.Our model works reasonably well with the IWSLT English-Vietnamese dataset in terms of BLEU score.The model shows significant amount of increase in score when we increase the recurrence depth one layer at a time. The increasing trend is largely visible when we change from 2 to 3 layers, a bit less from 4 to 5 layers. It is almost negligible when we change the recurrence depth from 5 to 6. Hence, we claim in machine translation task, we get the advantage of increasing the recurrence depth , which is similar to the results as mentioned in [1]; regarding the stack depth of LSTM's. Also to further produce adequate translation ,we use reconstructor model which also results in significant improvement.By using reconstructor block along with attention,we get improvement of almost 1 BLEU score compared to model without reconstruction.Also, we see that with beam search we were further able to improve translation and got improvement of around 1.5 BLEU score compared to benchmark result provided by [10]. [4] showed that deep LSTM are difficult to train compared to deep RHN .We can also see from Table 5 that the results are similar to work of [4] which shows that in language modelling task also with almost half the parameters RHNN are able to outperform LSTM based language model. The same is true with neural machine translation task.