Jie Chen
email: chenjieg@sina.com

Senior Member, IEEE Jose Carlos

M Bermudez
email: j.bermudez@ieee.org

Senior Member, IEEE Paul Honeine
email: paul.honeine@utt.fr

Non-negative least-mean-square algorithm

Dynamic system modeling plays a crucial role in the development of techniques for stationary and non-stationary signal processing. Due to the inherent physical characteristics of systems under investigation, non-negativity is a desired constraint that can usually be imposed on the parameters to estimate. In this paper, we propose a general method for system identification under non-negativity constraints. We derive the so-called non-negative least-mean-square algorithm based on stochastic gradient descent, and we analyze its convergence. Experiments are conducted to illustrate the performance of this approach and consistency with the analysis.

I. Introduction

In many real-life phenomena including biological and physiological ones, due to the inherent physical characteristics of systems under investigation, non-negativity is a desired constraint that can be imposed on the parameters to estimate in order to avoid physically absurd and uninterpretable results. For instance, in the study of a concentration field or a thermal radiation field, any observation is described with non-negative values (ppm, joule). Non-negativity as a physical constraint has received growing attention from the signal processing community during the last decade. For instance, consider the following non-negative least-square inverse problem min

x 1 2 Ax -b 2 subject to [x] i ≥ 0, ∀i, (1)
with A a real M ×N matrix of rank k ≤ min(M, N ), b an M -length real vector, and x an N -length real vector.

• denotes the Euclidean 2-norm and [ • ] i the i-th entry of the vector. This problem has been addressed in various contexts, with applications ranging from image deblurring in astrophysics [START_REF] Benvenuto | Nonnegative least-squares image deblurring: improved gradient projection approaches[END_REF] to deconvolution of emission spectra in chemometrics [START_REF] Van Benthem | Fast algorithm for the solution of large-scale non-negativity-constrained least squares problems[END_REF]. Another similar problem is the non-negative matrix factorization (NMF), which is now a popular dimension reduction technique [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF], [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF], [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations[END_REF]. Given a matrix X with non-negative entries, the squared error version of this problem can be stated as follows min

W ,H X -W H 2 F subject to [W ] ij ≥ 0, [H] ij ≥ 0, ∀i, j (2)
where • F denotes the Frobenius norm. This problem is closely related to the blind deconvolution one, and has found direct application in hyperspectral imaging [START_REF] Berry | Algorithms and applications for approximate nonnegative matrix factorization[END_REF]. Separation of non-negative mixture of non-negative sources has also been considered in [START_REF] Plumbley | Algorithms for nonnegative independent component analysis[END_REF], [START_REF] Moussaoui | Separation of non-negative mixture of non-negative sources using a bayesian approach and MCMC sampling[END_REF].

Over the last fifteen years, a variety of methods have been developed to tackle non-negative least-square problems (NNLS). Active set techniques for NNLS use the fact that if the set of variables which activate constraints is known, then the solution of the constrained least-square problem can be obtained by solving an unconstrained one that only includes inactive variables. The active set algorithm of Lawson and Hanson [START_REF] Lawson | Solving Least Squares Problems[END_REF] is a batch resolution technique for NNLS problems. It has become a standard among the most frequently used methods. In [START_REF] Bro | A fast non-negativity-constrained least squares algorithm[END_REF], Bro and De Jong introduced a modification of the latter, called fast NNLS, which takes advantage of the special characteristics of iterative algorithms involving repeated use of non-negativity constraints. Another class of tools is the class of projected gradient algorithms [START_REF] Rosen | The gradient projection method for nonlinear programming. part 1: Linear constraints[END_REF], [START_REF] Calamai | Projected gradient methods for linearly constrained problems[END_REF], [START_REF] Barzilai | Two-point step size gradient methods[END_REF], [START_REF] Theodoridis | Adaptive learning in a world of projections: a unifying framework for linear and nonlinear classification and regression tasks[END_REF]. They are based on successive projections on the feasible region. In [START_REF] Lin | Projected gradient methods for nonnegative matrix factorization[END_REF], Lin used this kind of algorithms for NMF problems. Low memory requirements and simplicity make algorithms in this class attractive for large scale problems. Nevertheless, they are characterized by slow convergence rate if not combined with appropriate step size selection. The class of multiplicative algorithms is very popular for dealing with NMF problems [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF], [START_REF] Lin | On the convergence of multiplicative update algorithms for nonnegative matrix factorization[END_REF].

Particularly efficient updates were derived in this way for a large number of problems involving non-negativity constraints [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF]. These algorithms however require batch processing, which is not suitable for online system identification problems.

In this paper, we consider the problem of system identification under non-negativity constraints on the parameters to estimate. The Karush-Kuhn-Tucker (KKT) conditions are established for any convex cost function, and a fixed-point iteration strategy is then applied in order to derive a gradient descent algorithm.

Considering the square-error criterion as a particular case, a stochastic gradient scheme is presented. A convergence analysis of this algorithm is proposed. The resulting model accurately predicts the algorithm behavior for both transient and steady-state conditions. Finally, experiments are conducted to evaluate the algorithm performance and its consistency with the analysis.

II. System identification with non-negativity constraints

Consider an unknown system, only characterized by a set of real-valued discrete-time responses to known stationary inputs. The problem is to derive a transversal filter model

y(n) = α x(n) + z 1 (n), (3)
with α = [α 1 , α 2 , . . . , α N ] the vector of the model parameters, and

x(n) = [x(n), x(n -1), . . . , x(n -N + 1)]
the observed data vector. The input signal x(n) and the desired output signal y(n) are assumed stationary and zero-mean. The sequence z 1 (n) accounts for measurement noise and modeling errors.

Due to the inherent physical characteristics of systems under investigation, in this paper, non-negativity is a desired constraint that is imposed on the coefficient vector α. Therefore, the problem of identifying the optimum model can be formalized as follows

α o = arg min α J(α) subject to α i ≥ 0, ∀i, (4)
with J(α) a continuously differentiable and strictly convex cost function in IR N , and α o the optimal solution to the constrained optimization problem.

A. A fixed-point iteration scheme

In order to solve the problem (4), let us consider its Lagrangian function Q(α, λ) given by [START_REF] Boyd | Convex Optimization[END_REF] Q(α, λ) = J(α)λ α, where λ is the vector of non-negative Lagrange multipliers. The Karush-Kuhn-Tucker conditions must necessarily be satisfied at the optimum defined by α o , λ o , namely,

∇ α Q(α o , λ o ) = 0 α o i [λ o ] i = 0, ∀i
where the symbol ∇ α stands for the gradient operator with respect to α.

Using ∇ α Q(α, λ) = ∇ α J(α) -λ,
these equations can be combined into the following expression

α o i [-∇ α J(α o )] i = 0, (5)
where the extra minus sign is just used to make a gradient descent of J(α) apparent. To solve Equation [START_REF] Cichocki | Nonnegative Matrix and Tensor Factorizations[END_REF] iteratively, two important points have to be noticed. The first point is that D(-∇ α J(α)) is also a gradient descent of J(α) if D is a symmetric positive definite matrix. The second point is that equations of the form ϕ(u) = 0 can be solved with a fixed-point iteration algorithm, under some conditions on function ϕ, by considering the problem u = u + ϕ(u). Implementing this strategy with Equation (5) leads us to the component-wise gradient descent algorithm

α i (n + 1) = α i (n) + η i (n)f i (α(n)) α i (n)[-∇ α J(α(n))] i (6)
with η i (n) a positive step size required to get a contraction scheme and to control the convergence rate.

Function f i (α) > 0 in [START_REF] Berry | Algorithms and applications for approximate nonnegative matrix factorization[END_REF] is the i-th entry of a diagonal matrix D. It is an arbitrary positive function of α. Some criteria J(α) are defined only for inputs α with positive entries, e.g., Itakura-Saito distance, Kullback-Leibler divergence. If necessary, this condition can be managed by an appropriate choice of the step size parameter. Let us assume that α i (n) ≥ 0. Non-negativity of α i (n + 1) is guaranteed if

1 + η i (n)f i (α(n))[-∇ α J(α(n)] i ≥ 0. (7)
If [∇ α J(α(n))] i ≤ 0, condition [START_REF] Plumbley | Algorithms for nonnegative independent component analysis[END_REF] is clearly satisfied and non-negativity does not impose any restriction on

the step size. Conversely, if [∇ α J(α(n))] i > 0, non-negativity of α i (n + 1) holds if 0 ≤ η i (n) ≤ 1 f i (α(n)) [∇ α J(α(n))] i . ( 8
)
Using a single step size η(n) in [0, η max (n)] for all entries of α so that

η max (n) = min i 1 f i (α(n)) [∇ α J(α(n))] i , i = 1, . . . , N (9)
the update equation can be written in vector form as

α(n + 1) = α(n) + η(n) d(n), (10)
where the weight adjustment direction d(n), whose i-th entry is defined as follows

[d(n)] i = f i (α(n)) α i (n)[-∇ α J(α(n))] i (11)
is a gradient descent direction because

f i [α(n)] α i (n) ≥ 0.
It should be noted that condition (9) on the step size η(n) guarantees the non-negativity of α(n) for all n, but does not ensure the stability of the algorithm.

B. The non-negative least-mean-square algorithm

Let us now consider the mean-square error criterion J mse (α) to be minimized with respect to α, that is,

α o = arg min α E{[y(n) -α x(n)] 2 } subject to α o i ≥ 0, ∀i, (12)
where we have included the non-negativity constraint only on the optimum solution because J mse (α) is defined for all α, that is, for all positive and negative entries α i . The gradient of J mse (α) can be easily computed as

∇ α J(α) = 2 (R x α -r xy ) (13)
with R x the autocorrelation matrix of x(n) and r xy the correlation vector between x(n) and y(n). Using [START_REF] Bro | A fast non-negativity-constrained least squares algorithm[END_REF] and [START_REF] Rosen | The gradient projection method for nonlinear programming. part 1: Linear constraints[END_REF] with f i (α) =1 2 for all i, the update rule for minimizing the mean-square error under non-negativity constraints is given by

α(n + 1) = α(n) + η(n) D α (n) (r xy -R x α(n)) (14)
where D α (n) is the diagonal matrix with diagonal entries given by α(n). Following a stochastic gradient approach, the second-order moments R x and r xy are replaced in ( 14) by the instantaneous estimates x(n) x (n) and y(n) x(n), respectively. This leads to the stochastic approximation of ( 14) given by 1

α(n + 1) = α(n) + η(n) e(n) D x (n) α(n), η(n) > 0 ( 15
)
where D x (n) stands for the diagonal matrix with diagonal entries given by x(n), and e(n

) = y(n)-α (n) x(n).
It is interesting to notice how the term α(n) in the update term on the right-hand side (r.h.s.) of [START_REF] Lin | Projected gradient methods for nonnegative matrix factorization[END_REF] affects the dynamics of the coefficient update when compared with the well know LMS algorithm [START_REF] Sayed | Adaptive Filters[END_REF]. Note that the extra multiplying factor α i (n) in the update term of the i-th row of (15), which is not present in the LMS update, provides extra control of both the magnitude and the direction of the weight update, as compared to LMS. For a fixed step size η, the update term for the i-th component of α(n) is proportional to -α i (n) e(n) x i (n), the stochastic gradient component. Thus, compared to the LMS stochastic gradient component -e(n) x i (n), the constrained algorithm includes the multiplying factor α i (n). A negative α i (n)

will then change the sign of the LMS adjustment, which on average tends to avoid convergence to negative coefficients of the unconstrained solution. Thus, coefficients that would normally converge, on average, to negative values using unconstrained LMS will tend to converge to zero using the constrained algorithm. In addition, α i (n) close to zero will tend to slow its own convergence unless the magnitude of e(n) x i (n) is very large. Finally, α i (n) = 0 is clearly a stationary point of Equation ( 15).

In the following, the adaptive weight behavior of the adaptive algorithm ( 15), called non-negative LMS, is studied in the mean and mean-square senses for a time-invariant step size η.

III. Mean behavior analysis

We shall now propose a model to characterize the mean behavior of the non-negative LMS algorithm.

α * = arg min α E{[y(n) -α x(n)] 2 }. ( 16
)
whose solution α * satisfies the Wiener-Hopf equations

R x α * = r xy . (17)
The residual signal 1 accounts for measurement noise and modeling errors.

z(n) = y(n) -x (n) α * in Figure
It is assumed in the following that z(n) is stationary, zero-mean with variance σ 2 z and statistically independent of any other signal. Thus, E{z(n) D x (n} = 0.

The adaptive algorithm [START_REF] Lin | Projected gradient methods for nonnegative matrix factorization[END_REF] attempts to estimate the optimum α o for the constrained problem [START_REF] Calamai | Projected gradient methods for linearly constrained problems[END_REF].

The analytical determination of the optimal solution α o is not trivial in general. In the particular case of independent and identically distributed (i.i.d.) input samples, however, R x = σ 2

x I where I is the identity matrix. In this case, the Karush-Kuhn-Tucker conditions imply that α o is obtained by turning the negative entries of α * to zero

α o = {α * } + ( 18
)
where {u} + = max{0, u}. The minimum mean-square error produced by solution α o is then

J ms min = σ 2 y -2 r xy {α * } + + σ 2 x {α * } + {α * } + (19)
with σ2 y the variance of y(n).

A. Mean weight behavior model

Defining the weight-error vector

v(n) = α(n) -α * = [v 1 (n), v 2 (n), . . . , v N (n)]
, the update equation [START_REF] Lin | Projected gradient methods for nonnegative matrix factorization[END_REF] can be written as

v(n + 1) = v(n) + η e(n) D x (n) v(n) + α * . ( 20
)
Using

e(n) = y(n) -α (n) x(n) = z(n) -v (n) x(n) leads us to the following expression v(n + 1) = v(n) + η z(n) D x (n) v(n) + η z(n) D x (n) α * (21) -η D x (n) v(n) v (n) x(n) -η D x (n) α * x (n) v(n).
Taking the expectation of expression [START_REF] Alligood | Chaos: an introduction to dynamical systems[END_REF], neglecting the statistical dependence of x(n) and v(n), 2 and using

that E{z(n) D x (n)} = 0 yields E{v(n + 1)} ≈ I -η E{D x (n) α * x (n)} E{v(n)} -η E{D x (n) v(n) v (n) x(n)}. (22)
The first expectation on the r.h.s. of ( 22) is given by In order to evaluate the second expectation on the r.h.s. of ( 22), let us compute the i-th component of the

E{D x (n) α * x (n)} = E{D α * x(n) x (n)} = D α * R x . ( 23
)
z(n) y(n) x(n) α(n) Algo. + + + - e(n) α *
vector D x (n) v(n) v (n) x(n). We have [D x (n) v(n) v (n) x(n)] i = N j=1 x(n -i + 1) v i (n) v j (n) x(n -j + 1) (24)
Taking the expectation of this expression, defining K(n) = E{v(n) v (n)}, and neglecting the statistical dependence of x(n) and v(n), we obtain

[E{D x (n) v(n) v (n) x(n)}] i ≈ N j=1 r x (j -i) [K(n)] ij (25) = [R x K(n)] ii This implies that E{D x (n) v(n) v (n) x(n)} ≈ diag{R x K(n)}
, where diag{A} denotes the vector whose i-th entry is defined by [A] ii . Using these results with Equation ( 22) yields the following mean weight-error vector update equation

E{v(n + 1)} = I -η D α * R x E{v(n)} -η diag{R x K(n)}. ( 26
)
This equation requires second-order moments defined by K(n) in order to update the first-order one E{v(n)}.

A recursive model will be derived for K(n) in Section IV, see Equation (39). That model can be used along with [START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF] to predict the mean weight behavior of the algorithm. Nevertheless, we have found that a sufficiently accurate and more intuitive mean behavior model can be obtained using the following separation approximation

K(n) ≈ E{v(n)} E{v (n)}. ( 27
)
Using ( 27) in [START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF] we obtain the following result

E{v(n + 1)} = I -η D α * R x E{v(n)} -η diag{R x E{v(n)} E{v (n)}}. ( 28
) Approximation (27) assumes that Cov{v i (n), v j (n)} E{v i (n)}E{v j (n)} (29)
In general, (29) is valid when the adaptive weights are far from convergence, as the mean weight-error component tends to be much larger than the weight-error fluctuation about the mean. For correlated x(n), the level of the weight-error fluctuations at convergence tends to be much less than the values of the nonzero optimal weights. For white input signals E{v i (n)} tends to zero for those indexes corresponding to the positive coefficients of α o . In this case, approximation (29) will in fact tend to eliminate the weight estimation error at convergence. Extensive simulation results have shown that the simplified model in (28) yields a prediction of the mean weight behavior which is sufficient for design purposes. Furthermore, this simplification allows the more detailed analytical study of the mean weight behavior shown in the next section.

B. Special case of a white input signal

In general, the behavior of (28) can become very complex to be studied analytically [START_REF] May | Simple mathematical models with very complicated dynamics[END_REF]. In order to obtain analytical results that allow some understanding of the mean weight behavior, we study here the particular case of x(n) i.i.d. and drawn from a zero-mean distribution. Unit variance σ 2

x is also assumed without loss of generality. Using R x = I in (28) yields the component-wise expression

E{v i (n + 1)} = (1 -η α * i ) E{v i (n)} -η E{v i (n)} 2 . ( 30
)
Function E{v i (n + 1)} in Equation ( 30) is a parabola in the variable E{v i (n)} with roots at E{v i (n)} = 0

and E{v i (n)} = 1-η α * i η
. It reaches its maximum value

(1-η α * i ) 2 4η at E{v i (n)} = 1-η α * i 2η . Fixed points are found by solving E{v i (n + 1)} = E{v i (n)}, which yields E{v i (n)} = 0 or E{v i (n)} = -α * i .
This result is consistent with solution [START_REF] Boyd | Convex Optimization[END_REF] where

v o i =    0 if α * i ≥ 0 -α * i otherwise (31) with v o i the i-th entry of v o = α o -α * .
Let us now derive conditions ensuring convergence of (30) to 0 and -α

* i . Writing u(n) = η 1-η α * i E{v i (n)},
where the index i has been dropped to simplify the notation, we obtain the following difference equation known as the logistic map [START_REF] May | Simple mathematical models with very complicated dynamics[END_REF], [START_REF] Alligood | Chaos: an introduction to dynamical systems[END_REF], [START_REF] Perrin | La suite logistique et le chaos[END_REF] u(n

+ 1) = ρ u(n) (1 -u(n)) (32)
with ρ = 1η α * i , which is assumed nonzero. Fixed points defined in (31) now correspond to u = 0 and u = ρ-1 ρ , respectively. Convergence of the logistic map to these values depends on the parameter ρ and on the initial condition u(0) as follows. See [START_REF] May | Simple mathematical models with very complicated dynamics[END_REF], [START_REF] Alligood | Chaos: an introduction to dynamical systems[END_REF], [START_REF] Perrin | La suite logistique et le chaos[END_REF] for details and Figure 2 for illustration.

Case 1: 0 < ρ < 1 An illustration of this case is shown in Figure 2 (left). The fixed point u = 0 attracts all the trajectories originating in the interval ρ-1 ρ ; Fixed points become unstable. New fixed points appear between which the system alternates in stable cycles of period 2 k , with k tending to infinity as ρ increases. This case may even lead to a chaotic behavior, and falls out of the scope of our study.

1 ρ . The logistic map u(n) is identically equal to ρ-1 ρ for n ≥ 1 if u(0) = ρ-1 ρ or u(0) = 1 ρ . Outside the interval, u(n) diverges to -∞.
To derive conditions for convergence of the difference equation (30) to 0 or -α * i , we must consider separately components of E{v i (n)} associated with positive or negative unconstrained optimum α * i , respectively. On the one hand, based on the analysis of the logistic map (32), convergence of (30) to 0 corresponds to the conditions on ρ and u(0) satisfying Case 1 and Case 2 above. This yields 0 < η < 1

α * i -α * i < v i (0) < 1 η (33)
in the case where α * i > 0. If α * i = 0, these two conditions become η > 0 and 0 < v i (0) < 1 η . On the other hand, ρ and u(0) must obey the conditions presented in Case 3 for convergence of Equation (30) to -α * i . This leads to 0 < η ≤ -2

α * i 0 < v i (0) < 1 η -α * i ( 34
)
in the case where α * i < 0. Finally, combining these inequalities leads to the following theoretical conditions for convergence of E{v(n)}:

0 < η ≤ min i 1 |α * i | and 0 < v i (0) < 1 η for all i (35)
or, using also Equations ( 33) and ( 34), for convergence of E{α(n)}:

0 < η ≤ min i 1 |α * i | and 0 < α i (0) < 1 η for all i. (36)
Conditions ( 35) and (36) on v i (0) and α i (0) show that there is more freedom in choosing α i (0) for small values of η. They guarantee the convergence of the difference equation (30).

C. Simulation examples for the first-order moment analysis

This section presents simulation examples to verify the validity of the first-order moment analysis of the non-negative LMS algorithm. We illustrate the accuracy of the model (30) through a first example where inputs x(n) and noise z(n) are i.i.d. and drawn from a zero-mean Gaussian distribution with variance σ 2 x = 1 and σ 2 z = 10 -2 , respectively. The impulse response α * is given by

α * = [0.8 0.6 0.5 0.4 0.3 0.2 0.1 -0.1 -0.3 -0.6] ( 37
)
The initial impulse response α(0) is drawn from the uniform distribution U([0; 1]), and kept unchanged for all the simulations. The algorithm's stability limit was determined experimentally to be η max ≈ 5 × 10 -3 . As usually happens with adaptive algorithms, this limit is more restrictive than the mean weight convergence limit given by (36), as stability is determined by the weight fluctuations [START_REF] Sayed | Adaptive Filters[END_REF]. The mean value E{α i (n)} of each coefficient is shown in Figure 3 for η = 10 -3 = η max /5 and η = 5 × 10 -5 = η max /10. The simulation curves (blue line) were obtained from Monte Carlo simulation averaged over 100 realizations. The theoretical curves (red line) were obtained from model (30). One can notice that all the curves are perfectly superimposed and, as predicted by the result [START_REF] Boyd | Convex Optimization[END_REF], each coefficient α i (n) tends to {α * i } + as n goes to infinity. It is interesting to note that convergence towards the solution {α * } + agrees with the theoretically predicted behavior of (32). For each positive entry α * i of α * , the corresponding value of ρ = 1η α * i is in ]0; 1[. This corresponds to Case 1 in Section III-B, where the fixed point u = 0 attracts all the trajectories and v i (n) converges to zero. It can also be verified that each ρ associated with a negative entry α * i is in ]1; 3]. This corresponds to Case 3 where u = (ρ -1)/ρ attracts all the trajectories and lim n→∞ v i (n) = -α * i .

The second simulation example illustrates the accuracy of the model (30) for inputs x(n) correlated in time. We consider a first-order AR model given by x(n) = a x(n -1) + w(n), with a = 1 2 . The noise w(n) is i.i.d. and drawn from a zero-mean Gaussian distribution with variance σ 2 w = 1 -1 4 , so that σ 2 x = 1 as in the first example. The other parameters of the initial experimental setup remain unchanged, except for the step size values. In order to verify the model's accuracy also for large step sizes we performed the simulations for η = 2.5×10 -3 = η max /2 and η = 5×10 -5 = η max /10. The mean value E{α i (n)} of each coefficient is shown in Figure 4. As before, the simulation curves (blue line) and the theoretical curves (red line) are superimposed. It can be noticed that α(n) no longer converges to {α * } + since the input samples x(n) are now correlated.

We can easily verify that E{e 2 (n)} = 4.97 dB using {α * } + , and E{e 2 (n)} = 3.82 dB at convergence of the non-negative LMS algorithm.

IV. Second-order moment analysis

We now present a model for the behavior of the second-order moments of the adaptive weights. To allow further analysis progress, we assume in this section that the input x(n) is Gaussian.

The expected values in (40)-( 48) are calculated in the following. In order to keep the calculations mathematically tractable, the following statistical assumptions are employed: A1: The input vector x(n) is zero-mean Gaussian.

A2: The weight-error vector v(n) is statistically independent of x(n) x (n). The reasoning for this approximation has been discussed in detail in [START_REF] Minkoff | Comment: On the unnecessary assumption of statistical independence between reference signal and filter weights in feedforward adaptive systems[END_REF].

A3:

The statistical dependence of v(n) v (n) and v(n) is neglected. This assumption follows the same reasoning valid for assuption A2, see [START_REF] Minkoff | Comment: On the unnecessary assumption of statistical independence between reference signal and filter weights in feedforward adaptive systems[END_REF].

A4: v(n) and (x (n) v(n)) 2 are statistically independent given A2. This is because (x (n) v(n)) 2 is a linear combination of the entries of v(n) v (n). Thus, this approximation follows basically the same reasoning discussed in [START_REF] Minkoff | Comment: On the unnecessary assumption of statistical independence between reference signal and filter weights in feedforward adaptive systems[END_REF] to justify A2.

Expected value P 1

This expected value has been already calculated in [START_REF] Minkoff | Comment: On the unnecessary assumption of statistical independence between reference signal and filter weights in feedforward adaptive systems[END_REF]. Remember that

P 1 = E{D x (n) α * } = E{D α * x(n) x (n)} = D α * R x . ( 49
)
Expected value P 2

Basic linear algebra gives

P 2 = E{D x (n) α * α * D x (n)} = E{D α * x(n) x (n) D α * } = D α * R x D α * . (50)
Expected value P 3

Neglecting the statistical dependence of x(n) and v(n) yields

P 3 = E{D x (n) v(n) α * D x (n)} ≈ E{D v (n)} R x D α * . ( 51
)
Expected value P 4

The (i, j)-th entry of the matrix within the expectation in P 4 is given by

[D x (n) v(n) v (n) D x (n)] ij = x(n -i + 1) v i (n) v j (n) x(n -j + 1). ( 52
) Using A2, E{x(n -i + 1) v i (n) v j (n) x(n -j + 1)} ≈ E{x(n -i + 1) x(n -j + 1)} E{v i (n) v j (n)} and P 4 ≈ R x • K(n) (53)
where • denotes the so-called Hadamard entry-wise product.

Expected value P 5

Defining D v (n) as the diagonal matrix with diagonal entries given by v(n), we first note that

E{v(n) x (n) v(n) v (n) D x (n)} = E{v(n) v (n) x(n) x (n) D v (n)}. ( 54
)
Now, taking the expected value with respect to v(n),

E{(x (n) v(n)) 2 x(n) x (n)} ≈ trace{R x K(n)} R x + 2 R x K(n) R x . (64)
Then we obtain the final result

P 8 ≈ E{D v (n)} (trace{R x K(n)} R x + 2 R x K(n) R x ) D α * . ( 65
)
Expected value P 9

Computing the (i, j)-th entry of matrix P 9 within the expectation, and using A2,

[P 9 ] ij = k E{x(n -i + 1) [v(n) v (n)] ik [x(n) x (n)] k [v(n) v (n)] j x(n -j + 1)} = k E{x(n -i + 1) x(n -k + 1) x(n -+ 1) x(n -j + 1)} E{v i (n) v j (n) v k (n) v (n)}. (66)
For x(n) zero-mean Gaussian (A1), we know that [START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF] E{x

(n-i+1) x(n-k+1) x(n-+1) x(n-j+1)} = r x (k-i) r x (j-)+r x ( -i) r x (j-k)+r x (j-i) r x ( -k). ( 67
)
The expectation

E{v i (n) v j (n) v k (n) v (n)
} cannot be evaluated directly, as the statistics of v(n) are unknown.

Approximate expressions can be obtained using numerous different approaches. We have chosen to use the following approximation which preserves relevant information about the second moment behavior of the adaptive weights while keeping the mathematical problem tractable. We first note that

E{v i (n) v k (n) v (n) v j (n)} = E{v i (n) v j (n)} E{v k (n) v (n)} + Cov{v i (n) v j (n), v k (n) v (n)}. (68)
Now, writing v i (n + 1) v j (n + 1) = v i (n) + η ∆v i (n) v j (n) + η ∆v j (n) , we see that the fluctuations in v i (n + 1) v j (n + 1) are proportional to η. Using the same reasoning for v k (n) v (n) we finally note that the covariance in (68) is proportional to η 2 . The higher order moments of the entries of v(n) in (68) will then be proportional to η p with p ≥ 2. Thus, for sufficiently small values of η, neglecting these terms yields the approximation

E{v i (n) v k (n) v (n) v j (n)} ≈ E{v i (n) v j (n)} E{v k (n) v (n)}. ( 69
)
This approximation is supported by the simulation results presented in Section IV-B. Substituting the two equations above into the expression of [P 9 ] ij leads to

[P 9 ] ij ≈ r x (j -i) k r x ( -k) [K(n)] k [K(n)] ij + k r x (k -i) r x (j -) [K(n)] k [K(n)] ij + k r x ( -i) r x (j -k) [K(n)] k [K(n)] ij . (70)
The first right-hand term of Equation ( 70) can be expressed as follows

r x (j -i) k r x ( -k) [K(n)] k [K(n)] ij = [R x ] ij k ( [R x ] k [K(n)] k )[K(n)] ij = [trace{R x K(n)} R x ] ij [K(n)}] ij . (71)
The second and third right-hand terms write

k r x (k -i) r x (j -) [K(n)] k [K(n)] ij = k [R x ] ik [K(n)] k x ] j [K(n)] ij = [R x K(n) R x ] ij [K(n)}] ij . (72)
This leads to the following close-form expression

[P 9 ] = (trace{R x K(n)} R x + 2 R x K(n) R x ) • K(n). (73)
Using the expected values P 1 to P 9 in Equation (39), we finally obtain a recursive analytical model for the behavior of K(n). This result can be used to study the convergence properties of E{e 2 (n)}, and can be applied to design. 4 The next section illustrates the model accuracy in predicting the non-negative LMS algorithm behavior.

B. Simulation examples for the second-order moment analysis

This section presents simulation examples to verify the accuracy of the model (39). Figures 5 and6 show the behavior of the excess MSE J emse (n) = trace{R x K(n)} corresponding to the experiments conducted in Section III-C. The simulation curves (blue line) were obtained from Monte Carlo simulation averaged over 100 realizations. The theoretical curves (red line) were obtained from model (39). Note the model's accuracy even for step sizes as large as η max /2 (left side of Figure 6). Also note that the theoretical value of the minimum excess mean-square error J emse min is represented in Figure 5. 5 It can be observed that J emse (n) tends to J emse min as n goes to infinity. Figure 7 highlights the performance of the model for uncorrelated and correlated input signals x(n) through the same experimental setup as described before, except that the noise variance σ 2 z is now set to 1. All these experiments illustrate the accuracy of the model, which can provide important guidelines for the use of the non-negative LMS algorithm in practical applications.

V. Conclusion

In many real-life phenomena, due to the inherent physical characteristics of systems under investigation, non-negativity is a desired constraint that can be imposed on the parameters to estimate in order to avoid physically absurd and uninterpretable results. In this paper, we proposed a general method for system identification under non-negativity constraints, and we derived the so-called non-negative LMS based on stochastic gradient descent. This algorithm switches automatically between a gradient descent mechanism and a gradient ascent one depending whether the non-negativity constraint is violated or not. Finally, we analyzed the algorithm convergence in the mean sense and in the mean-square sense. In future research efforts, we intend to explore these models in practical applications since they provide important guidelines to algorithm designers. We also plan to derive variants of this approach, e.g., in the spirit of the normalized-LMS and the sign-LMS algorithms.

Fig- ure 1

 1 shows a block diagram of the problem studied in this paper. The input signal x(n) and the desired output signal y(n) are assumed stationary and zero-mean. Let us denote by α * the solution of the unconstrained least-mean-square problem

Fig. 1 .

 1 Fig. 1. Adaptive system under study

Case 2 :Fig. 2 .

 22 Fig. 2. Convergence of the logistic map, in the Case 1 (left) and in the Case 3 (right). The dashed line is the line of equation u(n + 1) = u(n).

Fig. 3 .Fig. 4 .

 34 Fig. 3. Convergence of the coefficients α i (n) in the case where input x(n) and noise z(n) are i.i.d. Two different step sizes are considered: η = 10 -3 on the left figure, and η = 5 × 10 -4 on the right figure. The theoretical curves (red line) obtained from the model (30) and simulation curves (blue line) are perfectly superimposed.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Convergence of J emse (n) in the case where input x(n) and noise z(n) are i.i.d. Two different step sizes are considered: η = 10 -3 on the left figure, and η = 5 × 10 -4 on the right figure. The theoretical curves (red line) obtained from (38) and the model (39) and simulation curves (blue line) are perfectly superimposed.

Fig. 7 .

 7 Fig. 7. Convergence of J emse (n) with step size η = 5 × 10 -4 , in the case where input x(n) is i.i.d. on the left figure, and generated by a first-order AR process on the right figure. Compared to Figure 5 (right) and 6 (right), the variance of the noise z(n) has been increased from 10 -2 to 1.

Note that Dα(n) x(n) = Dx(n) α(n).

This assumption is less restrictive than the well-known independence assumption[19, p. 247], as it does not require x(n) be Gaussian.

The two important properties of z(n) used in evaluating (39) are its independence of any other signal and its zero-mean.

This model can also be used in[START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF] for the mean weight behavior if needed. However, our experience has been that the simplified model (28) suffices for predicting the mean weight behavior for most practical needs. It also makes the analytical study presented in Section III-B tractable.

It can be easily shown, from Equations (17)-[START_REF] Sayed | Adaptive Filters[END_REF], that Jemse min = α * -(α * )+ 2 in the case where Rx = I.

This work has been partly supported by CNPq grant No. 305377/2009-4.

A. Second moment behavior model

Using e(n) = z(n)v (n) x(n), neglecting the statistical dependence of x(n) and v(n), and using the properties assumed for z(n) yields an expression for the mean-square estimation error (MSE):

Equation ( 26) clearly shows that the mean behavior of each coefficient is a function of a single diagonal entry of matrix R x K(n). In this case, approximation (28) could be used without compromising the accuracy of the resulting mean behavior model. This accuracy has been verified through Monte Carlo simulations in Section III-C. The MSE in (38), however, is a function of the trace of R x K(n). Thus, the effect of the second order moments of the weight-error vector entries on the MSE behavior becomes more significant than in [START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF], and in general cannot be neglected. Thus, we determine a recursion for K(n) starting again from the weight error update equation [START_REF] Alligood | Chaos: an introduction to dynamical systems[END_REF].

Premultiplying Equation ( 21) by its transpose, taking the expected value, and using the statistical properties of z(n), 3 yields

where matrices P 1 to P 9 are defined by

Now, using A2 and A3, the expectation can be approximated as

Finally, using again A2 we obtain

Expected value P 6

Basic linear algebra gives

Under A1 and applying the same methodology used to derive [24, Equation ( 29)],

Expected value P 7

Using basic algebra, A2 and A3 as done to obtain (55), we have

Finally, under A1 and applying the same methodology as in [START_REF] Hubscher | An improved statistical analysis of the least mean fourth (LMF) adaptive algorithm[END_REF]Equation (29)], yields

Expected value P 8

Using basic algebra we obtain

Using A4, P 8 becomes

The expected value E{(x (n) v(n)) 2 x(n) x (n)} for zero-mean Gaussian signal x(n) has already been evaluated in [24, equations ( 7)-( 9)], using results from [START_REF] Bershad | Stochastic analysis of gradient adaptive identification of nonlinear systems with memory for gaussian data and noisy input and output measurements[END_REF]. Following the same procedure as in [START_REF] Hubscher | An improved statistical analysis of the least mean fourth (LMF) adaptive algorithm[END_REF] yields

(63)


Paul Honeine
email: paul.honeine@utt.fr

Cédric Richard
email: cedric.richard@unice.fr

Keywords: kernel-based machines, pre-image problem, linear algebra, kernel-PCA, nonlinear denoising

technique in kernel-based machines since, on one hand it involves only linear algebra operations, and on the other it can be written directly using the kernel values, without the need to evaluate distances as with the MDS approach. The performance of the proposed approach is illustrated for denoising with kernel-PCA, and compared to state-of-the-art methods on both synthetic datasets and realdata handwritten digits.

A closed-form solution for the pre-image problem in kernel-based machines

Paul Honeine • Cédric Richard

Received: date / Accepted: date Abstract The pre-image problem is a challenging research subject pursued by many researchers in machine learning. Kernel-based machines seek some relevant feature in a reproducing kernel Hilbert space (RKHS), optimized in a given sense, such as kernel-PCA algorithms. Operating the latter for denoising requires solving the pre-image problem, i.e. estimating a pattern in the input space whose image in the RKHS is approximately a given feature. Solving the pre-image problem is pioneered by Mika's fixed-point iterative optimization technique. Recent approaches take advantage of prior knowledge provided by the training data, whose coordinates are known in the input space and implicitly in the RKHS, a first step in this direction made by Kwok's algorithm based on multidimensional scaling (MDS). Using such prior knowledge, we propose in this paper a new technique to learn the pre-image, with the elegance that only linear algebra is involved. This is achieved by establishing a coordinate system in the RKHS with an isometry with the input space, i.e. the inner products of training data are preserved using both representations. We suggest representing any feature in this coordinate system, which gives us information regarding its pre-image in the input space. We show that this approach provides a natural pre-image 1 Introduction

In the last decade or so, kernel-based machines have enjoyed increasing popularity, providing a breakthrough in both statistical learning theory and low computational complexity of nonlinear algorithms. Pioneered by Vapnik's Support Vector Machines (SVM) [START_REF] Vapnik | Statistical Learning Theory[END_REF], this concept attracted significant attention due to the everexpanding challenges in machine learning. Since then, many nonlinear algorithms have been developed, for supervised learning (or classification) such as kernel Fisher discriminant analysis [START_REF] Mika | Kernel fisher discriminants[END_REF] and least-squares SVM [START_REF] Suykens | Least Squares Support Vector Machines[END_REF], and for unsupervised learning (with unlabelled data) with kernel principal component analysis (kernel-PCA) [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF] and support vector domain description [START_REF] Tax | One-class classification; concept-learning in the absence of counter-examples[END_REF]. The main idea behind nonlinear algorithms in kernelbased machines is the kernel trick [START_REF] Aizerman | Theoretical foundations of the potential function method in pattern recognition learning[END_REF]. This concept gives rise to nonlinear algorithms based on classical linear ones, under the only requirement that the algorithm can be expressed only in terms of inner products between data. Then, data from the input space are (nonlinearly) mapped into a feature space. This mapping is achieved implicitly by substituting the inner product operator by a positive definite kernel, thus without much additional computational cost. This is the essence of the kernel trick. In order to provide the unified functional framework, common to many communities, this kernel is called the reproducing kernel while the induced feature space is the reproducing kernel Hilbert space (RKHS).

With the ever-increasing demands in machine learning, new challenges require computing the inverse map. For instance, while the kernel trick provides an elegant approach to apply denoising or compression techniques in the RKHS, we need to go back into the input space for the final result. This is the case in denoising an image (or a signal), the reconstructed image belongs to the input space of training images. However, getting back to the input space from the RKHS is not obvious in general, as most features of the latter may not have an exact pre-image in the former. This is the preimage problem in kernel-based machines, as one seeks an approximate solution. Solving this problem has received a growing amount of attention, with the most breakthrough given in [START_REF] Mika | Kernel pca and de-noising in feature spaces[END_REF] and [START_REF] Kwok | The pre-image problem in kernel methods[END_REF]. In the former work, Mika et al. present the problem and its ill-posedness, and derive a fixed-point iterative scheme to find an approximate solution. Hence, there is no guarantee that this leads to a global optimum, and may be unstable. In the latter work, Kwok et al. determine a relationship between the distances in the RKHS and the distances in the input data, based on a set of training data. Applying a multidimensional scaling technique (MDS) leads to an inverse map estimate and thus to the preimage. This approach opens the door to a range of other techniques taking prior knowledge from training data in both spaces, such as manifold learning [START_REF] Etyngier | Shape priors using manifold learning techniques[END_REF] and out-ofsample methods [START_REF] Bengio | Out-of-sample extensions for lle, isomap, mds, eigenmaps, and spectral clustering[END_REF][START_REF] Arias | Connecting the out-of-sample and pre-image problems in kernel methods[END_REF].

In this paper, we propose a novel approach to solve the pre-image problem. To achieve this, we learn a coordinate system, not necessarily orthogonal, in the RKHS having an isometry with the input space. In other words, the inner products of the training data are (approximately) equal in both representations. Thus, by representing any feature of the RKHS in this coordinate system, we get an estimate of the inner products between the training data and its counterpart in the input space. It turns out that this approach is natural to kernel-based machines, and essentially requires only linear algebra, with any off-the-shelf linear solver. The proposed method is universal in the sense of being independent, in its formulation, of both the type of the adopted kernel and of the feature under investigation. Moreover, it extends naturally to get the pre-images of a set of features, since the coordinate system is computed only once.

The rest of the paper is organized as follows. In the next Section, we briefly present the framework behind kernel-based machines, with an illustration on kernel-PCA for denoising. In Section 3, the pre-image problem is described, and previous work on solving the problem are examined. The proposed method is described in Section 4, with connections to other methods. Experiments on synthetic and real datasets are presented in Section 5. Section 6 ends this paper with a brief conclusion.

2 Kernel-based machines, with application to nonlinear denoising using kernel-PCA

Kernel-based machines

Let X be a compact of IR p , equipped with the natural Euclidean inner product defined for any x i , x j ∈ X by x ⊤ i x j = p ℓ=1 x i,ℓ x j,ℓ , with x •,ℓ the ℓ-th entry of vector x • . Let κ(•, •) be a positive (semi-)definite kernel on X × X , where the positive (semi-)definitness is defined by the property i,j α i α j κ(x i , x j ) ≥ 0 for all α i , α j ∈ IR and x i , x j ∈ X . The Moore-Aronszajn theorem [START_REF] Aronszajn | Theory of reproducing kernels[END_REF] states that for every positive definite kernel, there exists a unique reproducing kernel Hilbert space (RKHS), and viceversa. With this one-toone correspondence between RKHS and positive definite kernels, the latter will be called reproducing kernels hereafter. Let H be the RKHS associated with κ, and let • , • H be the endowed inner product. This means that any arbitrary function ψ(•) in H can be evaluated at any x j ∈ X with

ψ(x j ) = ψ(•), κ(•, x j ) H . (1)
This expression shows that the kernel is the representer of evaluation. Moreover, replacing in this expression ψ(•) by κ(•, x i ) yields to the popular property

κ(x i , x j ) = κ(•, x i ), κ(•, x j ) H , (2)
for all x i , x j ∈ X . This is the reproducing property from which the name of reproducing kernel is derived.

Denoting by φ(•) the map that assigns to each input x ∈ X the kernel function κ(•, x), the reproducing property (2) implies that κ(x i , x j ) = φ(x i ), φ(x j ) H . The kernel then evaluates the inner product of any pair of elements of X mapped into H, without any explicit knowledge of either the mapping function φ(•) or the RKHS H. This is the well-known kernel trick. Examples of commonly used reproducing kernels are given in Table 1.

Table 1 Commonly used reproducing kernels in machine learning, with parameters β > 0, q ∈ IN, and σ > 0.

Polynomial κ(x i , x j ) = ( x i , x j + β) q Laplace κ(x i , x j ) = exp (-x i -x j /σ) Gaussian κ(x i , x j ) = exp `- x i -x j 2 /2σ 2 Ín
combination with the kernel trick, the representer theorem provides a powerful theoretical foundation for kernel-based machines. Initially derived in [START_REF] Kimeldorf | Some results on tchebycheffian spline functions[END_REF] and recently generalized in [START_REF] Schölkopf | A generalized representer theorem[END_REF], results of this theorem include SVM and kernel-PCA, where one seeks to maximize the separating margin between classes or the variance of projected data, respectively. This theorem states that any function ϕ * (•) of a RKHS H minimizing a regularized cost functional of the form

n i=1 J(ϕ(x i ), y i ) + g( ϕ 2 H ),
with predicted output ϕ(x i ) for input x i , and eventually the desired output y i , and g(•) a strictly monotonically increasing function on IR + , can be written as a kernel expansion in terms of available data

ϕ * (•) = n i=1 γ i κ(•, x i ). (3)
This theorem shows that even in an infinite dimensional RKHS, as with the Gaussian kernel, we only need to work in the subspace spanned by the n kernel functions of the training data, κ(•, x 1 ), . . . , κ(•, x n ).

Kernel-PCA for denoising

An elegant kernel-based machine is the kernel-PCA [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF], a nonlinear extension of one of the most used dimension reduction and denoising technique, the principal component analysis (PCA).

With PCA, one seeks principal axes that capture the highest variance in the data, that is, useful information as opposed to noise, and thus projecting data onto the space spanned by these relevant axes yields a denoising scheme. These principal axes are the eigenvectors ϕ k associated with the largest eigenvalues λ k of the covariance matrix R of data, i.e. solving the eigen-problem Rϕ k = λ k ϕ k . There exists another formulation of the PCA algorithm, using only inner products of the training data. By substituting the inner product operator with any valid reproducing kernel, we get an implicit nonlinear mapping of the data into a RKHS. This is the kernel-PCA, where each of the resulting principal functions takes the representer form (3), with

ϕ * k (•) = n i=1 γ i,k κ(•, x i ).
The weighting coefficients γ i,k are obtained from the eigen-decomposition of the so-called Gram matrix K, whose entries are κ(x i , x j ), for i, j = 1, . . . , n, by solving Kγ = nλ γ.

In order to have a PCA interpretation in feature space, two issues should be carried out. First, data is implicitly centered in feature space by substituting in this expression K with (1 -

1 n )K(1 -1 n )
, with 1 n the n-by-n matrix of entries 1/n and 1 the identity matrix; second, principal functions are normalized to 1, by scaling expansion coefficients such that n i=1 γ 2 i,k = 1/λ k . In the same spirit of the conventional PCA, one constructs a subspace of H spanned by the most relevant principal functions. Using kernel-PCA for denoising any given x ∈ X , we project the associated kernel function κ(•, x) onto that subspace. Since this subspace is spanned by most relevant principal functions, each of the form (3), any function from this subspace takes the same form, i.e. can be written as a kernel expansion in terms of available data. Let ϕ * (•) be this projection, with

ϕ * (•) = n i=1 γ i κ(•, x i ),
which is assumed to be noise-free by virtue of the PCA interpretation. From this denoised feature, we need to get its counterpart in the input space, e.g. a denoised image in the image space. As illustrated in Figure 1, this requires the estimation of the pattern x * from κ(•, x), by solving the pre-image problem.

A brief review of the pre-image problem

For supervised learning, one seeks a prediction value associated to any input such as in regression problems, while in classification this value is compared to a threshold, which yields a decision rule. While every optimal function ϕ * (•) takes the form (3), we obtain its evaluation at any x with n i=1 γ i κ(x i , x), thus requiring only computing values of the kernel. For pattern recognition with unsupervised learning, one is often interested in the feature in the feature space, or more precisely in its counterpart in the input space.

Estimating the input whose map is an arbitrary function in the RKHS is an ill-posed problem. To show

x x 1 x 2 x 3 xn x * φ(•) κ(•, x) κ(•, x 1 ) κ(•, x 3 ) κ(•, x 2 ) ϕ * (•) ? X H
Fig. 1 Illustration of the pre-image problem in kernel-based machines.

this, recall that the dimensionality of the feature space can be very high, and even infinite with some kernels such as the Gaussian kernel. Thus, (most) features ϕ * (•) ∈ H might not have an existing pre-image in X , i.e. a x * such that κ(•, x * ) = ϕ * (•). In order to circumvent this difficulty, one seeks an approximate solution, i.e. x * ∈ X whose map κ(•, x * ) is as close as possible to ϕ * (•). This is the pre-image problem in kernel-based machines. Methods for solving the pre-image problem are roughly classified into two categories: Fixed-point iterative methods and methods based on learning the inverse map.

The pre-image problem was initially studied by Mika et al. in [START_REF] Mika | Kernel pca and de-noising in feature spaces[END_REF]. They proposed to solve the optimization problem

x * = arg min x∈X ϕ * (•) -κ(•, x) 2 H , (4)
where • H denotes the norm in the RKHS. For this purpose, a fixed-point iterative scheme is used to solve the pre-image problem. However, since this optimization problem is highly non-convex, such iterative technique suffers from numerical instabilities and local minima. The pre-image will highly depend on the initial guess and is likely to get stuck in a local minimum. A further improvement of the fixed-point iterative scheme is presented in [START_REF] Rathi | Statistical shape analysis using kernel pca[END_REF], where the authors operate additional approximations by, roughly speaking, substituting the mapping κ(•, x) in (4) with its projection onto the subspace. It is worth noting that as an alternative to Mika's distance minimization, one may consider a collinearity maximization problem [START_REF] Arias | Connecting the out-of-sample and pre-image problems in kernel methods[END_REF], with

x * = arg max x∈X ϕ * (•) ϕ * (•) H , κ(•, x) κ(•, x) H H .
A key ingredient of these methods is a high dependence on the kernel type, since the fixed-point can only be applied to some specific kernels such as the Gaussian kernel, and only more recently extended to polynomial kernels in [START_REF] Kwok | The pre-image problem in kernel methods[END_REF].

Recent approaches take advantage of prior knowledge provided by some available training data, whose coordinates are available in both the input and the feature spaces. This approach is initiated by an algorithm based on multidimensional scaling (MDS), presented by Kwok et al. [START_REF] Kwok | The pre-image problem in kernel methods[END_REF]. This is achieved by computing distances between every pair of training data, in both spaces. For each pair, the classical Euclidean distance is used in the input space, as well as the distance in the RKHS which can be computed using kernel values. With these pairs of distances, a MDS technique is considered by performing a singular-value-decomposition1 . This yields an inverse map, in the same spirit of the outof-sample extension [START_REF] Arias | Connecting the out-of-sample and pre-image problems in kernel methods[END_REF]. In order to make this method tractable in practice, only the neighboring data affect the pre-image estimation. Learning the inverse map is studied in [START_REF] Bakir | Learning to find preimages[END_REF] by solving a regression problem, while alternative approaches can be based on the manifold learning [START_REF] Etyngier | Shape priors using manifold learning techniques[END_REF]. All these methods take advantage of prior knowledge, i.e. training data with information available in both input and feature spaces.

Using such prior knowledge, we propose in this paper to learn the inverse map without the need to compute distances, and does not require sophisticated optimization schemes. Only conventional linear algebra are needed. Furthermore, it is universe, in the sense that it does not depend on the kernel type, as opposed to fixed-point iterative techniques.

The proposed pre-image method

Given a set of training data {x 1 , x 2 , . . . , x n }, we seek the pre-image in X of some arbitrary ϕ * (•) of the RKHS H, denoted x * . The proposed method can be organized into two stages: learning the inverse map and operating a pre-image. To learn the inverse map, a coordinate system is constructed in the RKHS, having an isometry with the input space coordinates, where the isometry is defined with respect to the training data. In order to operate a pre-image, we represent ϕ * (•) in this coordinate system which, by virtue of the isometry, gives the values of the inner products of its pre-image with the training data in the input space. From these values we obtain the pre-image x * .

Stage 1: Learn the inverse map

In this stage, we provide a coordinate system in the RKHS that is isometric with the input space. In order to achieve such isometry, we consider a set of n training data {x 1 , x 2 , . . . , x n } ∈ X . By virtue of the representer theorem, we only need to consider the subspace spanned by their kernel functions {κ(•, x 1 ), κ(•, x 2 ), . . . , κ(•, x n )}. Within this subspace, we define the set of ℓ coordinate functions, denoted {ψ 1 (•), ψ 2 (•), . . . , ψ ℓ (•)} with ℓ ≤ n, and write

ψ k (•) = n i=1 α k,i κ(•, x i ),
for k = 1, 2, . . . , ℓ. For any kernel function κ(•, x), its coordinate on ψ k (•) is given by

ψ k (•), κ(•, x) H = ψ k (x) = n i=1 α k,i κ(x i , x),
where (1) is used. Therefore, its representation in this coordinate system is obtained by the ℓ coordinates, written vector-wise as

Ψ x = [ψ 1 (x) ψ 2 (x) • • • ψ ℓ (x)] ⊤ ,
where the k-th entry depends on the α k,i , for i = 1, . . . , n.

In order to estimate the coordinate functions, we propose an equivalence, between the inner products in this coordinate system and their counterparts in the canonic input space, using the model

Ψ ⊤ xi Ψ x j = x ⊤ i x j + ǫ ij , (5)
for all the training set, i.e. i, j = 1, 2, . . . , n, and where ǫ ij corresponds to the lack-of-fit of the model. We insist on the fact that this model is not coupled with any constraint on the coordinate functions, as opposed to the orthogonality between the functions resulting from the kernel-PCA. The only requirement we impose is the isometry defined in [START_REF] Bengio | Out-of-sample extensions for lle, isomap, mds, eigenmaps, and spectral clustering[END_REF]. The minimization of the variance of ǫ ij , a lack-of-fit criterion, consists of solving the optimization problem min ψ1,...,ψ ℓ 1 2 n i,j=1

x ⊤ i x j -Ψ ⊤ xi Ψ x j 2 + λ R(ψ 1 , . . . , ψ ℓ ).
As suggested in machine learning literature, we include in this expression a regularization term, where λ a tunable parameter controlling the tradeoff between the fitness to the model ( 5) and the smoothness of the solution. In order to penalize high norm functions, the regularization

R(ψ 1 , . . . , ψ ℓ ) = ℓ k=1 ψ k 2 H is used in this paper.
This optimization problem can be written in matrix form. This is done by a factorization of Ψ x into a matrix of unknowns and a vector of available information, with

Ψ x = A κ x , where κ x = [κ(x 1 , x) κ(x 2 , x) • • • κ(x n , x)] ⊤ and A is a ℓ × n matrix of unknowns whose (k, i)-th entry is α k,i .

This leads to the optimization problem

A = arg min A 1 2 n i,j=1 x ⊤ i x j -κ ⊤ xi A ⊤ A κ x j 2 +λ ℓ k=1 n i,j=1
α k,i α k,j κ(x i , x j ).

By denoting • F the Frobenius norm2 of a matrix and tr(•) its trace, this yields

A = arg min A 1 2 P -K A ⊤ A K 2 F + λ tr(A ⊤ AK),
where P and K are the Gram matrices with entries x ⊤ i x j and κ(x i , x j ), respectively. Taking the derivative of the above cost function with respect to A ⊤ A, rather than A, and setting it to zero, we get

A ⊤ A = K -1 P -λ K -1 K -1 . ( 6
)
In what follows, we show that only A ⊤ A is required to find the pre-image, rather than A. Fortunately, we do not need to compute the coefficients α k,j to generate the coordinate system in the RKHS; only their inner products are required.

Stage 2: Operate a pre-image Since the model ( 5) is valid for all the training data, we apply it to do the pre-image, as discussed in this stage. Let ϕ * (•) be any optimal function resulting from a kernel-based machine, with ϕ * (•) = n i=1 γ i κ(•, x i ) as given in [START_REF] Aronszajn | Theory of reproducing kernels[END_REF]. By virtue of the representer theorem, it belongs to the subspace spanned by the training kernel functions, and therefore can be expressed in terms of the provided coordinate system. The coordinate of ϕ * (•) associated to the coordinate function

ψ k (•) is ϕ * (•), ψ k (•) H = n i,j=1 α k,i γ j κ(x i , x j ).
Each of these ℓ coordinates are computed and collected into one vector, denoted Ψ ϕ * with some abuse of notation. Thus, we extend the model ( 5), and write

Ψ ⊤ xi Ψ ϕ * = x ⊤ i x * ,
for i = 1, 2, . . . , n, where x * is the pre-image to be estimated. This identity can be expressed matrix-wise with

K A ⊤ A Kγ = X ⊤ x * where γ = [γ 1 γ 2 • • • γ n ] ⊤ and X = [x 1 x 2 • • • x n ]
. By injecting the provided system (6) into this expression, we get

X ⊤ x * = (P -λ K -1 ) γ. (7)
This is a classical system of linear equations. Thus, the pre-image can be estimated by applying any offthe-shelf solver. For instance, one can solve the linear least-squares optimization problem

x * = arg min x X ⊤ x -(P -λ K -1 ) γ 2 , (8)
where any iterative or non-iterative technique can be used, such as the pseudo-inverse or the eigendecomposition3 , in the spirit of the Nyström method. It is worth noting that the optimization scheme is applied here to the input space, as opposed to high dimensional RKHS with the fixed-point iteration schemes. Moreover, one needs only to consider solution from the span of the training data, in coherence with previous work on the pre-image problem [START_REF] Mika | Kernel pca and de-noising in feature spaces[END_REF][START_REF] Kwok | The pre-image problem in kernel methods[END_REF]. The proposed method is universal in the sense of being independent, in its formulation, of both the type of the adopted kernel and of the feature under investigation.

In order to better understand this result, consider the potential theoretical setting of linear independent training data. In this case, the minimization problem (8) has a unique solution, given by solving the normal equations XX ⊤ x * = X P -λ K -1 γ. By using the pseudo-inverse matrix algebra with the identity

(XX ⊤ ) -1 X = X(X ⊤ X) -1 , we get x * = XP -1 P -λ K -1 γ. (9)
Extension to a set of features These expressions can be applied readily to a set of features in the RKHS to get their pre-images in the input space. This can be done straightforwardly by writing [START_REF] Hoffmann | Kernel PCA for novelty detection[END_REF] as

X ⊤ X * = (P -λ K -1 ) Γ ,
where each column of matrix Γ represents the coefficient vector γ, and each column of X * the corresponding pre-image. From the solution (9), we see that the matrix

M = XP -1 P -λ K -1
needs to be computed only once, and then applied with

X * = M Γ .
This corresponds to a matrix completion scheme, or more specifically the kernel matrix regression approach, as given in [START_REF] Yamanishi | Kernel matrix regression[END_REF][START_REF] Honeine | Localization in sensor networks -a matrix regression approach[END_REF].

Experiments

In this section, we compare the proposed method with two state-of-the art methods 4 : the fixed-point iterative technique [START_REF] Mika | Kernel pca and de-noising in feature spaces[END_REF] and the MDS-based approach [START_REF] Kwok | The pre-image problem in kernel methods[END_REF]. For this purpose, the kernel-PCA for denoising is applied on synthetic and real datasets. The Gaussian kernel κ(x i , x j ) = exp( x ix j 2 /2σ 2 ) is used, its bandwidth σ fixed to the same value for the three methods.

Synthetic datasets

We consider a family of four datasets in 2-D, each having a geometric form corrupted by a noise of bandwidth parameter ν. The data are sampled uniformly randomly within this area. We generate n train data to train the n eigen eigenfunctions and to construct the coordinate system. Then, we apply these results on another set of n pre-image generated data, in order to denoise using the pre-image techniques. For each dataset, the parameters values are summarized in Table 2.

The frame dataset consists of a square with sides of length 2. Data are generated uniformly randomly on each side and corrupted by a noise uniformly distributed on the interval [-ν, ν] normal to the side. The banana dataset is given by the parabola defined by the coordinates (x, x 2 + ξ), with x on the x-axis uniformly distributed on the interval [-1, 1], and ξ normally distributed with a standard deviation of ν. The spiral is defined by the coordinates (A(ϕ) cos(ϕ), A(ϕ) sin(ϕ)), with A(ϕ) = 0.07ϕ + ξ, where ϕ and ξ are generated uniformly on the intervals [0, 6π] and [0, ν], respectively. The sine dataset is defined by the coordinates (ϕ, 0.8 sin(2ϕ)), where ϕ is generated uniformly on the interval [0, 2π], and corrupted with an additive uniformly distributed noise in the range [0, ν] 2 . See [START_REF] Hoffmann | Kernel PCA for novelty detection[END_REF] for more information.

The fixed-point iterative algorithm is set with a stopping criterion of maximum 100 iterations, reaching the limit of reasonable cpu time. The initial estimate is chosen from the valid model x * = i γ i x i , with the weighting coefficients γ i generated uniformly on the interval [-1, 1]. The MDS-based algorithm operates using a global optimization scheme, which gives better results than the neighborhood setting. Since this algorithm is based on an eigen-decomposition technique, it results in a new coordinate system in the input space. Hence, we consider a procrustes technique to align it with the initial canonical one, by minimizing the mean-squares error.

In Figure 2, we show the four datasets, with on the one hand the training data (blue dots), and on the other the denoised estimates (red dots) obtained from another set of noisy data (not shown here, yet given by the unmarked ends of green lines). Green lines show the distance between the denoised and the initial noisy data.

The fixed-point iterative method suffers on one side from numerical instabilities, illustrated through many estimates falling outside the bounds of the images (following the long green lines), and on the other from local minima, illustrated with improper denoising (for in-stance, the upper border of the frame dataset (y-axis close to 1) are not denoised to the same area). It is obvious that the MDS-based approach is clearly inappropriate to any of the given datasets. The method presented in this paper gives good results with the four proposed datasets, with the smallest reconstruction error of all algorithms. It seems less sharper in denoising than the fixed-point iterative algorithm, without suffering from the drawbacks of the latter. However, it causes the estimates to fold over itself, in the same sense of manifold learning. This is illustrated for instance with the banana data, yet much less pronounced than the MDSbased results.

Real datasets

We illustrated the efficiency of our method on denoising real datasets. We consider the handwritten digit "2", obtained from the MNIST database of handwritten digits [START_REF] Lecun | The mnist database of handwritten digits[END_REF]. The images are (almost) binary images of 28-by-28 pixels. Hence, from a machine learning point of view, each image is simply a point in the 784-dimensional space. The original images were corrupted by adding a zero-mean white Gaussian noise with variance ν = 0.1. A set of n train = 1 000 images are used to train the kernel-PCA with the n eigen = 100 leading principal functions retained. We apply the Gaussian kernel to all three algorithms, with bandwidth set to σ = 105 . The parameter settings are summarized in Table 3.

To illustrate the denoised ability of each algorithm, another set of n pre-image = 10 images is considered under the same noise conditions. These images are illustrated in Figure 3 (first row), with results from the fixed-point iterative (second row), the MDS-based (third row) and the proposed (fourth row) algorithms. It is obvious that fixed-point iterative algorithm is inappropriate for such application, even with the number of maximum iterations set to 10 000 corresponding to an average total CPU time of up to 1 hour and a half. To take advantage of prior knowledge, the same training set is used for learning the inverse map. Realistic results can be obtained using the MDS-based algorithm, with 5 minutes and a half. The algorithm proposed in this paper achieves better denoised results, as illustrated in Figure 3. For this simulation, the regularization parameter was set to λ = 10 -9 , and the resulting average total CPU time is 1.3 seconds 5 . Fixed-point iterative method [START_REF] Mika | Kernel pca and de-noising in feature spaces[END_REF] MDS-based method [START_REF] Kwok | The pre-image problem in kernel methods[END_REF] Method proposed in this paper Fig. 2 Experimental results for the frame (first row), the banana (second row), the spiral (third row), and the sine (fourth row) datasets, using the fixed-point iterative (left), the MDS-based (middle), and the proposed (right) algorithms. Training data are represented by blue dots, estimated pre-images by red dots, and green lines illustrate the distance between these estimates and the initial noisy data (not shown).

10 digits corrupted by noise Fixed-point iterative method [START_REF] Mika | Kernel pca and de-noising in feature spaces[END_REF] MDS-based method [START_REF] Kwok | The pre-image problem in kernel methods[END_REF] Method proposed in this paper Fig. 3 Comparative analysis for denoising a set of ten "2" digits (first row), with the denoised images from the fixed-point iterative (second row), the MDS-based (third row), and the proposed (fourth row) algorithms.

Table 3 Values of the parameters for the real digit dataset.

n train n pre-image ν n eigen σ 1 000 10 0.1 100 10 5

Conclusion

In this paper, we presented a new method to solve the pre-image problem. As opposed to previous work, the proposed method neither suffers from numerical instability, nor requires computing the distances in the input and the RKHS spaces. We showed that using the inner product information in both spaces, we can provide a coordinate system in the RKHS to learn the inverse map. The efficiency of the proposed method were studied with experiments on both synthetic data and real handwritten digits, and compared to state-of-the-art methods. The major advantage of the proposed method resides on its simplicity in dealing with the optimization issue, thanks to conventional linear algebra.

Table 2

 2 Values of the parameters for the synthetic datasets.

	4 Matlab	codes	for	these	algorithms	are	avail-
	able from the Statistical Pattern Recognition Toolbox
	http://cmp.felk.cvut.cz/cmp/software/stprtool/

This is done by operating on the distances, transforming them into inner products, and then apply eigen-decomposition into the resulting Gram matrix to get the coordinates. This nicely captures our guiding intuition of the problem in contrast with the MDS : we propose to work exclusively on the inner products, without the need to compute distances.

The Frobenius norm of a matrix is the root of sum of squared (absolute) values of all its elements, or equivalently M 2 F = tr(M ⊤ M).

Doing eigen-decomposition gives the pre-image relative to the eigen-basis in the input space. A post-processing is required to set the pre-image relative to the training data; this is called the procrustes problem.

CPU times are given only as an indication of the computations required for the various algorithms.