Jean-Baptiste Hiriart-Urruty

Hai Yen Le
email: hyle@math.univ-toulouse.fr

A VARIATIONAL APPROACH OF THE RANK FUNCTION

Keywords: Rank of matrix, Optimization, Nonsmooth analysis, Moreau-Yosida regularization, Generalized subdifferentials

In the same spirit as the one of the paper [24] on positive semidefinite matrices, we survey several useful properties of the rank function (of a matrix) and add some new ones.

Since the so-called rank minimization problems are the subject of intense studies, we adopt the viewpoint of variational analyis, that is the one considering all the properties useful for optimizing, approximating or regularizing the rank function.

Introduction

Associated with (square) matrices are some familiar notions like the trace, the determinant, etc. Their study from the variational point of view, via the usual differential calculus, is easy and well-known. We consider here another function of (not necessarily square) matrices, the rank function. The rank function has been studied for its properties in linear algebra (or matrix calculus), semi-algebraic geometry, etc. One task in the present paper is to survey some of the main properties of the rank in these classical areas of mathematics. But we are more interested in considering the rank function from the variational viewpoint. Actually, besides being integervalued, the rank function is lower-semicontinuous; this is the only valuable topological property it enjoys. But, since modern nonsmooth analysis allows us to deal with non-differentiable (even 1 discontinuous) functions, a second task we intend to carry out is to answer questions like: what is the generalized subdifferential (in whatever sense) of the rank function?

If we are interested in the rank function from the variational viewpoint, it is because the rank function appears as an objective (or constraint) function in various modern optimization problems. The archetype of the so-called rank minimization problems is as follows:

(P)    Minimize f (A) := rank of A subject to A ∈ C,
where C is a subset of M m,n (R) (the vector space of m by n real matrices). The constraint set is usually rather "simple" (expressed as linear equalities, for example), the main difficulty lies in the objective function. Of course, we could have an optimization problem

(P 1 )    Minimize g(A)
subject to A ∈ C and rank A ≤ k, with a rather simple objective function but a fairly complicated constraint set. Both problems (P) and (P 1 ) suffer from the same intrinsic difficulty: the occurence of the rank function.

A related (or cousin) problem to (P), actually equivalent in terms of difficulty, stated in R n this time, consists in minimizing the so-called counting (or cardinality) function x = (x 1 , . . . , x n ) ∈ R n -→ c(x) := number of nonzero components x i of x:

(Q)    Minimize c(x) subject to x ∈ S,
where S is a subset of R n . Often c(x) is denoted as x 0 , although it is not a norm. Problems (P) and (Q) share some bizarre and/or interesting properties, from the optimization or variational viewpoint. We shall review some of them like those related to relaxation, global optimization, Moreau-Yosida approximations, generalized subdifferentials.

Here is the plan of our survey. The properties of the rank function are well-known in the context of linear algebra or matricial calculus. Any book on these subjects presents them; there even are compilations of them (see [START_REF] Bernstein | Matrix mathematics: theory, facts, and formulas[END_REF] for example). We recall here the most basic ones: While properties like (b) are the ones used in solving systems of linear equations for example, properties and inequalities in (c)-(d) are those which are useful in optimization. The second property in (c) shows us that it is unnecessary to look at the rank of matrices lying outside some ball centered at 0 and of positive radius. Property (a) involves positive semidefinite (symmetric) matrices A T A and AA T ; therefore, a priori there should not be a loss of generality in considering positive semidefinite matrices, but it turns out that it is not the case in rank minimization problems (we have to treat of A ∈ M m,n (R) as they are given in the data of the problems).

Since rank is a complicated integer-valued function, are there ways of bounding it by using more familiar and easier to handle functions like the trace, the determinant, the eigenvalue functions? Indeed, there are some, not very convincing however. Here are two prosaic ones.

Theorem 1. (i) Let

A ∈ M n (R) be non-null diagonalizable, with all real eigenvalues (for example A ∈ S n (R)). Then

rank A ≥ (tr A) 2 tr (A 2 ) , ( 1
)
where tr B denotes the trace of the matrix B.

(ii) Let A ∈ S n (R) be non-null and positive semidefinite. Then

tr A λ max (A) ≤ rank A ≤ tr A λ min (A) , ( 2
)
where λ max (A) (resp. λ min (A)) stands for the maximal (resp. minimal) eigenvalue of A.

(If c > 0, c divided by 0 is put equal to +∞).

Proof. We content ourselves by sketching the proof of the inequality (1).

For given real numbers λ 1 , . . . , λ k , let us denote by λ :

= λ 1 +•••+λ k k their mean value. From the inequality k i=1 (λ i -λ) 2 ≥ 0 comes the inequality below k i=1 λ 2 i ≥ ( k i=1 λ i ) 2 k = k λ2 . ( 3
)
Let thus A = 0 of rank k. If λ 1 , . . . , λ n denote the eigenvalues of A, we may suppose, without loss of generality, that λ 1 = 0, . . . , λ k = 0, λ k+1 = 0, . . . , λ n = 0 (since k ≥ 1).

Now, tr A = k i=1 λ i , tr (A 2 ) = k i=1 λ 2 i . We then infer from (3) k ≥ ( k i=1 λ i ) 2 k i=1 λ 2 i = (tr A) 2 tr (A 2
) .

Comments and examples

• Having an equality in (1) requires too specific properties on A, like: there exists k such that A 2 = kA.

• The eigenvalues of A do not appear explicitly in the inequality (1). If the eigenvalues λ i of A are available, a more informative inequality is

rank A ≥ ( k i=1 |λ i |) 2 tr(A 2 ) .
This comes from the Cauchy-Schwarz inequality.

• Bounds ( 1) and ( 2), although rough, could be used to relax constraints on the rank of a parameterized matrix: rank [A(x)] ≤ k could be replaced by [tr

A(x)] 2 ≤ k tr [A(x) 2 ] or tr [A(x)] λmax[A(x)] ≤ k.
To take a simple example, consider

x ∈ R -→ A(x) =   1 x x 1   ∈ S 2 (R). Requiring that A(x) is positive semidefinite amounts to having x ∈ [-1; 1]; requiring futhermore, that rank [A(x)] ≤ 1 leads to x = -1 or x = 1. Now keeping the constraint "A(x) is positive semidefinite", one relaxes the constraint rank [A(x)] ≤ 1 with tr [A(x)] 2 ≤ tr [A(x) 2 ]. This leads to x 2 ≥ 1, that is x = -1 or x = 1 again.
Let us look at the second relaxation, the one coming from (2). Since tr [A(x)] = 2 and

λ max [A(x)] = 1 + |x|, substituting the relaxed form tr [A(x)] λmax[A(x)] ≤ 1 for rank [A(x)] ≤ 1 gives rives to 2 1+|x| ≤ 1, that is again x = -1 or x = 1. • Note that A → tr A is a linear function, A → λ max (A) a convex one, A → λ min (A) a
concave one. The bounding functions in (2) are nonsmooth but continuous on the domains where they are defined.

The rank from the topological and semi-algebraic viewpoint

The only (useful) topological property of the rank function is that it is

lower-semicontinuous: if A ν → A in M m,n (R) when ν → +∞, then lim inf ν→+∞ rank A ν ≥ rank A. ( 4
)
This is easy to see if one thinks of rank A characterized as the maximal integer r such that a (r, r)-submatrix extracted from A is invertible; the continuity of the determinant function makes the rest.

Since the rank function is integer-valued, a consequence of the inequality ( 4) is that the rank function does not decrease in a sufficiently small neighborhood of any matrix A.

For k ∈ {0, 1, . . . , p}, consider now the following two subsets of M m,n (R):

S k := {A ∈ M m,n (R)| rank A ≤ k}, Σ k := {A ∈ M m,n (R)| rank A = k}.
S k is the sublevel-set (at level k) of the lower semicontinuous function rank; it is therefore closed.

But, apart from the case k = 0 (where S 0 = Σ 0 ), what about the topological structure of Σ k ?

The answer is given in the following statement.

Theorem 2. (i) Σ

p is an open dense subset of S p = M m,n (R).

(ii) If k < p, the interior of Σ k is empty and its closure is S k .

The singular value decomposition of matrices (see below) will show how intricate the subsets S k and Σ k may be. For example, if rank A = k, in any neighborhood of A, there exist matrices

of rank k + 1, k + 2, . . . , p.
Now, what about S k from the algebraic geometry viewpoint? Consider m = n for the sake of simplicity. The geometric structure of the S k 's is well understood (cf. [START_REF] Parrilo | The convex algebraic geometry of rank minimization[END_REF]): since S k is defined by the vanishing of all (k +1, k +1)-minors of A, it is thus a solution set of polynomial equations, therefore a so-called semi-algebraic variety. It is non-singular, except on those points (matrices) of rank less than or equal to k -1. Its dimension is (2n -k)k. At a smooth point (matrix) A of S k , the tangent space T S k (A) can be made explicit, provided that one knows a singular value decomposition of A ( [45, p.15]).

Given A ∈ M m,n (R) of rank r and an integer k ≤ r, how are made the matrices in S k closest to A? This best approximation problem needs firstly that a distance (via a norm for instance)

be defined on M m,n (R). Secondly, even if the existence of best approximants does not offer any difficulty (remember that S k is closed), the question of uniqueness as well as that of an explicit form of best approximants remain posed. It turns out that there is a beautiful theorem answering these questions.

Before going further, we recall a technique of decomposition of matrices which is central in numerical matricial analysis and in statistics: the singular value decomposition (SVD). Here it is: Given A ∈ M m,n (R), there is an (m, m) orthogonal matrix U , an (n, n) orthogonal matrix V , a "pseudo-diagonal" matrix D1 of the same sizes as A, such that A = U DV2 .

Figure 1: The singular value decomposition

The picture above helps to understand the decomposition.

• D has the same number of columns and rows as A, it is a sort of skeleton of A: all the "non-diagonal" entries of D are null; on the "diagonal" of D are the singular values σ 1 , σ 2 , . . . , σ p of A, that are the square roots of the eigenvalues of A T A (or AA T ). By definition, all the σ i 's are nonnegative, and exactly r of them (if r = rank A) are positive.

By changing the ordering in columns or rows in U and V , and without loss of generality, we can suppose that

σ 1 ≥ σ 2 ≥ • • • ≥ σ k > σ k+1 = • • • = σ p = 0.
• U and V are orthogonal matrices of appropriate sizes (so that the product U DV of matrices can be performed). Of course, these U and V are not unique. We denote by O(m, n) A the set of (U, V ) such that U and V are orthogonal matrices and U diag m,n (σ 1 (A), . . . , σ p (A))V T is a singular value decomposition of A.

Clearly, if one has a SVD of A, A = U DV , one also has one of A T , namely

A T = V T D T U T .
The SVD helps to define some norms on M m,n (R). The three most important ones are:

-A -→ max i=1,...,p σ i (A), called the spectral norm of A, which we denote by A sp ( A 2 is another much used notation);

-A -→ tr (A T A) = p i=1 σ 2 i (A)
, called the Frobenius-Schur norm of A, which we denote by A F . This norm is derived from the inner product A, B := tr (A T B) on M m,n (R); thus (M m,n (R), . F ) is an Euclidean space.

-A -→ p i=1 σ i (A), called the nuclear norm (or trace norm, or 1-norm, etc.) which we denote by A * . These three norms are the "matricial cousins" of the usual norms . ∞ , . 2 , . 1 on R p .

The best approximation problem that we consider in this section is as follows: Given A ∈ M m,n (R) of rank r,

(A k )    Minimize A -M M ∈ S k
The problem of course depends on the choice of . ; it is solved when . is either the Frobenius-Schur norm or the spectral norm.

Theorem 3 (Eckart and Young, Mirsky [START_REF] Higham | Matrix nearness problems and applications[END_REF]). Let A ∈ M m,n (R) and let U DV be a SVD of A. Choose . as either . F or . sp . Then

A k := U D k V,
(where D k is obtained from D by keeping σ 1 , . . . , σ k and putting 0 in the place of σ k+1 , . . . , σ r ) is a solution of the best approximation problem (A k ). For the Frobenius-Schur norm case, A k is the unique solution in (A k ) when σ k > σ k+1 .

The optimal values in (A k ) are as following:

min M ∈S k A -M F = r i=k+1 σ 2 i ; min M ∈S k A -M sp = σ k+1 .
This theorem is a classical result in numerical matricial analysis, proved for the Frobenius- ) is convex and smooth. But, due to the non-convexity of the constraint set S k , the optimization problem (A k ) is non-convex. It is therefore astonishing that one could provide (via Theorem 3) an explicit form of a solution of this problem. In short, since the S k 's are the sublevel-sets of the rank function, one always has at one's disposal a "projection" of (an arbitrary) matrix A on the sublevel-sets S k . Remark 2. From the differential geometry viewpoint, Σ k is a smooth manifold around any matrix A ∈ Σ k , while S k is not. However, S k enjoys a very specific property called prox-regularity 3 . Indeed, it is proved in various places like [START_REF] Lewis | Alternating projections on manifolds[END_REF], [START_REF] Absil | Projection-like retractions on matrix manifolds[END_REF], [START_REF] Luke | Prox-regularity of rank contraint sets and implications for algorithms[END_REF] that S k is prox-regularity at any A ∈ S k with rank A = k. Moreover, the projection on S k is also the projection on Σ k , at least locally; more precisely, according to [START_REF] Absil | Projection-like retractions on matrix manifolds[END_REF]: Let A ∈ Σ k , then for every A such that A -A ≤ σ k (A)/2, the projection of A onto Σ k exists, is unique, and is the projection of A onto S k (given in the Eckart and Young theorem). (

) 5
The MATLAB software uses this (SVD decomposition) to compute the rank of a matrix.

Among the many ways to compute the rank of a matrix, this one is maybe the most time consuming, but also the most reliable.

Global minimization of the rank

Consider the rank minimization problem (P) that we alluded to in the beginning of the paper. The constraint set C does not contain the zero matrix; if it were so, the zero matrix would be trivially a global minimizer in (P). If only local minimizers are sought, the problem does not make sense. As explained in [START_REF] Hiriart-Urruty | When only global optimization matters[END_REF], every point (matrix) A ∈ C is a local minimizer in (P) 4 . So, here, instead of tackling directly the hard problem (P), a much used approach in optimization is to consider a relaxed form of (P) by minimizing a more tractable objective function.

Relaxed forms of the rank function

Given the rank function rank : M m,n (R) -→ R, what is its closest "convex relative", that is to say its convex hull 5 ? Posed as it is, this question has no interest; indeed, due to the properties (c) in Section 1, it is easy to understand that co(rank), the convex hull of the rank function, is the zero function (on M m,n (R)). So, we restrict the rank function to an appropriate ball, by assigning the value +∞ out of this ball:

A ∈ M m,n (R) -→ rank R (A) :=    rank of A if A sp ≤ R; +∞ otherwise. (6)
We call this function the "restricted rank function". Now, an interesting question is: what is the convex hull of rank R ? To the best of our knowledge, this question was firstly answered in Fazel's PhD thesis ( [START_REF] Fazel | Matrix rank minimization with applications[END_REF]).

4 Such a phenomenon happens because the objective function is not continuous. Indeed, if C is a connected set, if f : C -→ R is a continuous function, and if every point in C is a local minimizer or a local maximizer of f , then f is necessarily constant on C ( [START_REF] Behrends | Functions for which all points are local extrema[END_REF]). 5 The convex hull or convex envelope of a function f is the largest possible convex function below f . Its epigraph (i.e., what is above the graph) is exactly the convex hull of the epigraph of f .

Theorem 4 (Fazel [18]). We have:

A ∈ M m,n (R) -→ co(rank R )(A) :=    1 R A * if A sp ≤ R; +∞ otherwise.
There are various ways of proving this theorem; we explain and comment three of them.

Proof 1. It consists in calculating the Legendre-Fenchel conjugate φ * of φ := rank R , and then the biconjugate φ * * (:= (φ * ) * ) of φ. We know from convex analysis that, since φ is bounded from below, φ * * is the closed convex hull of rank R , that is: φ * * = co(rank R ). But here, there is no difference between co(rank R ) and its closure co(rank R ). This is the way followed by Fazel ([18], Section 5.1.4).

Proof 2.

Here, one begins by convexifying the restricted counting (or cardinality) function (on R p ), and then one applies fine results by Lewis on how to get propertiers on rank R from those on c R .

The so-called restricted counting function c R on R p is defined as follows:

x ∈ R p -→ c R (x) :=    c(x) if x ∞ ≤ R; +∞ otherwise.
The next statement gives the convex hull of c R .

Theorem 5. We have

x ∈ R p -→ co(c R )(x) :=    1 R x 1 if x ∞ ≤ R; +∞ otherwise.
This result is part of the "folklore" in the areas where minimization problems involving counting functions appear (there are dozens of papers in signal recovery, compressed sensing, statistics, etc.). A direct proof is presented in [START_REF] Le | Confexifying the counting function on R p for convexifying the rank function on M m,n (R)[END_REF].

A. Lewis ([36], [START_REF] Lewis | The convex analysis of unitarily invariant matrix functions[END_REF]) showed that the Legendre-Fenchel conjugate of a function of matrices (satisfying some specific properties) could be obtained by just conjugating some associated function of the singular values of A. Using his results twice, one is able to calculate the Legendre-Fenchel biconjugate of the (restricted) rank function by calling on the biconjugate of the (restricted) counting function. Indeed, when dealing with matrices, we have: for

A ∈ M m,n (R), rank R (A) = c R [σ(A)], (7)
where σ(A) = (σ 1 (A), . . . , σ p (A)) is the vector of R p made up with the singular values σ i (A) of A. It happens that the assumptions of Lewis' "transfer results" from R p to M m,n (R) are satisfied in our context (expressed in [START_REF] Barioli | The maximal cp-rank of rank k completely positive matrices[END_REF]). In doing so, ones retrieves Theorem 4.

Proof 3.

Here we look at the problem with "geometrical glasses": instead of convexifying functions directly, we firstly convexify (appropriate) sets of matrices and, then, get convex hulls (or quasiconvex hulls) of functions as by-products. This is the path followed in ( [START_REF] Hiriart-Urruty | Convexifying the set of matrices of bounded rank: applications to the quasiconvexification and convexification of the rank function[END_REF]).

For that purpose, let

S R k := S k ∩ {A ∈ M m,n (R)| A sp ≤ R} = {A ∈ M m,n (R)| rank A ≤ k and A sp ≤ R}.
In this definition, k ∈ {0, 1, . . . , p} and R ≥ 0; the restriction " A sp ≤ R" plays the role of a "moving wall". Although S R k is a very complicated set of matrices (due to the definition of S k ), its convex hull has a fairly simple expression.

Theorem 6 (Hiriart-Urruty and Le, [START_REF] Hiriart-Urruty | Convexifying the set of matrices of bounded rank: applications to the quasiconvexification and convexification of the rank function[END_REF]). We have

co S R k = {A ∈ M m,n (R)| A sp ≤ R and A * ≤ Rk}. (8)
So, according to (8), the convex hull of S R k is the intersection of two balls, one for the nuclear norm . * , the other one for the spectral norm (the bound on . sp is the same in co S R k as in S R k ). Getting at such an explicit form of co S R k is due to the happy combination of these specific norms ( . sp and . * ). If . were any norm on M m,n (R) and

S R k := {A ∈ M m,n (R)| rank A ≤ k and A ≤ R},
due to the equivalence between the two norms . and . sp , we would get with (8) an inner estimate set and an outer estimate set of co

S R k . A function f : X -→ R ∪ {+∞} is called quasi-convex if all its sublevel sets [f ≤ α] := {x ∈ X| f (x) ≤ α} (sublevel set of f at the level α ∈ R)
are convex. Since the supremum of an arbitrary collection of quasiconvex functions is quasiconvex, one can define the largest quasiconvex function lower bounding f ; it is called the quasiconvex hull of f and denoted by f q . The way to construct f q from the sublevel sets of f is easy

x ∈ X -→ f q (x) = inf{α| x ∈ co[f ≤ α]}.
Following Theorem 6, we therefore get the explicit form of the quasiconvex hull of the (restricted) rank function.

Theorem 7 (Hiriart-Urruty and Le, [START_REF] Hiriart-Urruty | Convexifying the set of matrices of bounded rank: applications to the quasiconvexification and convexification of the rank function[END_REF]). We have:

A ∈ M m,n (R) -→ (rank R ) q (A) =    1 R A * if A sp ≤ R, +∞ otherwise,
where a stands for the smallest integer which is larger than (or equal to) a.

A final step remains to be carried out: to go from the quasiconvex hull to the convex hull of rank R . This is not difficult and has been done in [START_REF] Hiriart-Urruty | Convexifying the set of matrices of bounded rank: applications to the quasiconvexification and convexification of the rank function[END_REF].

Remark 4. If, in Theorem 4, the restriction " A sp ≤ R" is replaced by " A * ≤ R", the result (on the relaxed form) will be the same. But this does not hold true for a restriction " A ≤ R" written for any matricial norm . . Remark 6. The question of the convex relaxation of the restricted rank function, solved in Theorem 4, could also be posed for tensors (or hypermatrices). Actually, provided that the notion of rank of a tensor and various norms of a tensor have been defined in an appropriate way, it has been proved in [START_REF] Lim | Multiarray signal processing: Tensor decomposition meets compressed sensing[END_REF] that the nuclear norm of a tensor is a convex underestimator of the tensor rank on the unit ball for the Schatten ∞-norm (corresponding to the spectral norm for matrices). The question whether this nuclear norm is the largest convex underestimator of the rank (for tensors) was raised by the authors in [START_REF] Lim | Multiarray signal processing: Tensor decomposition meets compressed sensing[END_REF]; they conjectured that was the case. In this context, contrary to the world of matrices, there is no SVD decomposition (for tensors) as a tool at our disposal.

Surrogates of the rank function

When A is positive semidefinite, some surrogates for rank A can be defined with the help of trace and determinant functions; these give rise to heuristics in rank minimization problems ( [START_REF] Fazel | Matrix rank minimization with applications[END_REF],

Chapter 5). For general matrices, there is no hope to give reliable and good approximations of the (discontinuous) rank function by smooth or even continuous functions involing only norms or inner product of matrices. We nevertheless mention some of them, the first one due to Flores

([19], Chapter II.3).
A proposal of a continuous surrogate of the rank is as follows:

A ∈ M m,n (R) -→ g sf (A) :=      A 2 * A 2 F if A = 0, 0 if A = 0. ( 9
)
This function has a close relationship with the rank function, and also shares with it some similar properties.

Theorem 8 (Flores,[START_REF] Flores | Problèmes d'optimisation globale en statistique robuste[END_REF]). We have:

(G1) g sf (cA) = g sf (A) for c = 0. (G2) 1 ≤ g sf (A) ≤ rank A for A = 0.
(G3) If all the non-zero singular values of A are equal, then g sf (A) = rank A.

(G4) rank A = 1 if and only if g sf (A) = 1.

Properties (G2) and (G4) above allow us to say that, for A = 0, the rank-one constraint could be replaced by a difference-of-convex (dc) inequality constraint; indeed, if A = 0,

rank A = 1 ⇐⇒ A * -A F ≤ 0. (10)
A property like (G4) was also used by Malick ([44]) in SDP problems with rank-one constraints.

Let

C := {A ∈ S n (R)| A positive semidefinite, diag A = (1, . . . , 1)}
the convex set of the so-called correlation matrices. In ( [START_REF] Malick | The spherical constraint in Boolean quadratic programs[END_REF], Theorem 1), the following was proved for A ∈ C:

rank A = 1 ⇐⇒ A F = n. ( 11
)
This property was at the origin of the designation "spherical constraint" as a subtitute for the rank-one constraint. The equality constraint A 2 F = n 2 , although nonlinear, can therefore be dualized; see [START_REF] Malick | The spherical constraint in Boolean quadratic programs[END_REF] for more on that.

The underestimator g sf of the rank funtion turns out to look similar to a large class of underestimators parameterized by matrices. The Theorem 9 below, due to Lokam ( [START_REF] Lokam | Spectral methods for matrix rigidity with applications to size-depth trade-offs and communication complexity[END_REF]), deserves to be better known, in our opinion. We provide here a proof for the convenience of the reader.

Theorem 9 (Lokam,[START_REF] Lokam | Spectral methods for matrix rigidity with applications to size-depth trade-offs and communication complexity[END_REF]). For any non-null B ∈ M m,n (R), define

A ∈ M m,n (R -→ g B (A) :=      | A, B | A sp B sp if A = 0, 0 if A = 0. Then g B (A) ≤ rank A for all A ∈ M m,n (R). ( 12
)
Proof.

(
(σ i -τ i ) 2 ≤ A -B 2 F . ( 13
)
Developing the left-hand side of (13) gives

p i=1 σ 2 i + p i=1 τ 2 i -2 p i=1 σ i τ i = A 2 F + B 2 F -2 p i=1 σ i τ i .
Now, if r denotes the rank of A (assumed non-null),

p i=1 σ i τ i = r i=1 σ i τ i ≤ r A sp B sp .
Since the norm . F is derived from the inner product ., . , developing the right-hand side of [START_REF] Drusvyatskiy | The dimension of semialgebraic subdifferential graphs[END_REF] gives

A -B 2 F = A 2 F + B 2 F -2 A, B .
Thus, combining all, a consequence of the inequality ( 13) is

A, B ≤ r A sp B sp .
Changing B into -B does not change B sp . Whence the desired inequality

| A, B | ≤ r A sp B sp is proved. Remark 7.
If B is chosen to be equal to A, Theorem 9 gives us an underestimator of the rank function,

A ∈ M m,n (R) -→ g(A) = A 2 F A 2 sp ,
which is of the same family as the g sf function.

Remark 8. If m = n and B is chosen to be the identity matrix, the inequality (12) reads as

|tr A| A sp ≤ rank A. ( 14
)

Regularization-approximation of the rank function

Besides the underestimates of the rank function proposed in the previous section, it is conceivable -although difficult, as we already said it-to design smooth or just continuous approximations R ε of the rank fucntion, depending on some parameter ε > 0. What we expect for R ε is a general comparison result with the rank, for example: R ε (A) ≤ rank A for all ε > 0, as also a convergence result: R ε → rank A when ε → 0. We propose here two classes of regularizationapproximation of the rank function: the first one consists of smoothed versions of the rank, the second one relies on the so-called Moreau-Yosida technique, widely used in the context of variational analysis and optimization.

Smoothed versions of the rank function

The underlying idea is as follows: Since

rank A = c[σ 1 (A), . . . , σ p (A)] (recall that c is the counting function on R p ) = p i=1 θ[σ i (A)], (15)
where θ(x) = 1 if x = 0, θ(0) = 0 (θ is the counting function on R), we have to design some smooth approximation of the θ function. Because all the σ i (A) are nonnegative, θ acts only on the nonnegative part of the real line, so we can choose even approximations of the θ function (like θ which is itself even).

A first example was proposed in [START_REF] Hiriart-Urruty | Deux questions de rang[END_REF], it is as following: For ε > 0, let θ ε be defined as

x ∈ R -→ θ ε (x) := 1 -e -x 2 /ε . ( 16
)
A ∈ M m,n (R) -→ R ε (A) := p i=1 [1 -e -σ 2 i (A)/ε ]. ( 17
)
An alternate expression of the R ε function is

R ε (A) = p -tr(e -A T A/ε ). ( 18
)
Then, R ε is a C ∞ (even analytic) function of A. The properties of R ε as an approximation of the rank function are summarized in the statement below.

Theorem 10. We have

(i) R ε (A) ≤ rank A for all ε > 0.
(ii) The sequence of functions (R ε ) ε>0 increases when ε decreases, and

R ε (A) → rank A for all A when ε → 0. (iii) If A = 0 and r = rank A, rank A -R ε (A) ≤ ε r i=1 1 σ 2 i (A) , ( 19
)
as also

rank A -R ε (A) ≤ ε 2 r i=1 1 σ 4 i (A) . ( 20
)
Sketch of the proof. The results in (i) and (ii) easily follow from the definitions ( 16) and [START_REF] Elad | Sparse and redundant representations: From theory to applications in signal and image processing[END_REF].

The upper bounds [START_REF] Flores | Problèmes d'optimisation globale en statistique robuste[END_REF] and [START_REF] Higham | Matrix nearness problems and applications[END_REF], which are not comparable, come from the study of the function x → e -x x and e -x x 2 on the half line R + .

We do not know if such an approximation scheme is of any interest for solving rank minimization problems.

Another proposal for approximating the rank function, a quite recent one, is due to Zhao

([51]
). It consists of using, for all ε > 0, the following even approximation of the θ function:

x ∈ R -→ τ ε (x) := x 2 x 2 + ε . ( 21
)
A ∈ M m,n (R) -→ Z ε (A) := p i=1 σ 2 i (A) σ 2 i (A) + ε . ( 22
)
Alternate expressions of the Z ε function are:

Z ε (A) = tr[A(A T A + εI n ) -1 A T ] = n -εtr(A T A + εI n ) -1 .
Here also, Z ε is a C ∞ (even analytic) function of A. The properties of Z ε as an approximation of the rank function are summarized in the next statement.

Theorem 11 (Zhao,[START_REF] Zhao | An approximation theory of matrix rank minimization and its application to quadratic equations[END_REF]). We have

(i) Z ε (A) ≤ rank A for all ε > 0.
(ii) The sequence of functions (Z ε ) >0 increases when ε decreases, and Z ε (A) → rank A for all A when ε → 0.

(iii) If A = 0 and r = rank A,

rank A -Z ε (A) = r i=1 ε σ 2 i (A) + ε ≤ ε r i=1 1 σ 2 i (A) . ( 23
)
The use of this function Z ε (instead of the rank function) in rank minimization problems as well as an application to solving a system of quadratic functions are discussed in ( [START_REF] Zhao | An approximation theory of matrix rank minimization and its application to quadratic equations[END_REF], Sections 3 and 4).

Moreau-Yosida approximation-regularization of the rank function

Although the rank function is a bumpy one, it is lower-semicontinuous and bounded from below; it therefore can be approximated-regularized in the so-called Moreau-Yosida way. Surprisingly enough, the rank function is amenable to such an approximation-regularization process, and we get at the end explicit forms of the approximated-regularized versions in terms of singular values (like in the previous subsection). Let us firstly recall what is known, as a general rule, for the Moreau-Yosida approximation-regularization technique in a nonconvex context (see [47, Section 1.G] for details, for example).

Let (E, . ) be an Euclidean space, let f : E -→ R ∪ {+∞} be a lower-semicontinuous function, bounded from below on E. For a parameter value λ > 0, the Moreau-Yosida approximate (or Moreau envelope6 ) function f λ and proximal set-valued mapping Prox λ f are defined by

f λ (x) := inf u∈E {f (u) + 1 2λ x -u 2 }, (24)
Prox λ f (x) := {u ∈ E| f (u) + 1 2λ x -u 2 = f λ (x)}. ( 25
)
Then:

(i) f λ is a finite-valued continuous function on E;

(ii) The sequence of functions (f λ ) λ>0 increases when λ decreases, and f λ (x) → f (x) for all x when λ → 0;

(iii) The set Prox λ f (x) is nonempty and compact.

(iv) The lower bounds of f and f λ on E are equal:

inf x∈E f (x) = inf x∈E f λ (x).
We now apply this process to the rank function (or restricted rank function). The context is therefore as following: E = M m,n (R), . F is the Frobenius-Schur norm and f : E -→ R∪{+∞} is the rank (or restricted rank) function. By definition,

(rank R ) λ (A) = inf B ∈ Mm,n(R) B sp≤R rank B + 1 2λ A -B 2 F . ( 26
)
The limiting case, i.e. with R = +∞, is

(rank) λ (A) = inf B∈Mm,n(R) rank B + 1 2λ A -B 2 F . ( 27
)
The next two theorems are taken from [START_REF] Hiriart-Urruty | From Eckart & Young approximation to Moreau envelopes and vice versa[END_REF] or [START_REF] Le | [END_REF]Chapter 3].

Theorem 12. We have, for all A ∈ M m,n (R) with rank r ≥ 1:

(i) (rank) λ (A) = 1 2λ A 2 F - 1 2λ r i=1 [σ 2 i (A) -2λ] + . ( 28
)
(ii) One minimizer in [START_REF] Hiriart-Urruty | Convex analysis and minimization algorithms[END_REF], i.e. one element in Prox λ (rank)(A), is provided by

B := U Σ B V ,
where

• (U, V ) ∈ O(m, n) A , i.e. U and V are orthogonal matrices such that A = U Σ A V , with Σ A = diag m,n [σ 1 (A), .
. . , σ r (A), 0, . . . , 0] (a singular value decomposition of A with

σ 1 (A) ≥ • • • ≥ σ r (A) > 0); • Σ B =              0 if σ 1 ≤ √ 2λ, Σ A if σ r (A) ≥ √ 2λ, diag m,n [σ 1 (A), . . . , σ k (A), 0, . . . , 0] if there is an integer k such that σ k (A) ≥ √ 2λ > σ k+1 (A). (29)
We may complete the result (ii) in the theorem above by determining explicitly the whole set Prox λ (rank)(A). Indeed, we have four cases to consider:

• If σ 1 (A) < √ 2λ, then Prox λ (rank)(A) = {0}. • If σ r (A) > √ 2λ, then Prox λ (rank)(A) = {A}. • If there is k such that σ k (A) > √ 2λ > σ k+1 (A), then for any (U 1 , V 1 ) and (U 2 , V 2 ) in O(m, n) A , U 1 diag m,n [σ 1 (A), . . . , σ k (A), 0, . . . , 0]V 1 = U 2 diag m,n [σ 1 (A), . . . , σ k (A), 0, . . . , 0]V 2 So, the set Prox λ (rank)(A) is a singleton and Prox λ (rank)(A) = {U diag m,n [σ 1 (A), . . . , σ k (A), 0, . . . , 0]V } with (U, V ) ∈ O(m, n) A .
• Suppose there is k such that σ k (A) = √ 2λ. We then define

k 0 := min{k| σ k (A) = √ 2λ}, k 1 := max{k| σ k (A) = √ 2λ}.
Then, Prox λ (rank)(A) is the set of matrices of the form U diag m,n (τ 1 , . . . , τ p )V , where

(U, V ) ∈ O(m, n) A and τ i = σ i (A) if i < k 0 , τ i = 0 if i > k 1 ; τ i = 0 or σ i (A) if k 0 ≤ i ≤ k 1 .

Comments

1. One could wish to express (rank) λ (A) in terms of traces of matrices as this was done for the smoothed versions of the rank function in Section 7.1. Indeed, A T A -2λI n is a symmetric matrix whose eigenvalues are σ 2 1 (A) -2λ, . . . , σ 2 r (A) -2λ, -2λ, . . . , -2λ. Its projection on the cone S + n (R) of positive semidefinite matrices has eigenvalues [σ 2 1 (A) -2λ] + , . . . , [σ 2 r (A) -2λ] + , 0, . . . , 0 ( [START_REF] Hiriart-Urruty | A fresh variational analysis look at the world of the positive semidefinite matrices[END_REF]). Thus, an alternate expression for (rank) λ (A) is:

(rank) λ (A) = 1 2λ tr(A T A) - 1 2λ tr[P S + n (R) (A T A -2λI n )]. (30)
2. If λ is small enough, say if

√ 2λ ≤ σ r (A), then (rank) λ (A) = rank A.
This easily comes from (28) since σ 2 i (A) -2λ ≥ 0 for all i and A 2 F = r i=1 σ 2 i (A). Therefore, the general convergence result that is known for the Moreau-Yosida approximates f λ of f is made much stronger here.

The counterpart of Theorem 12 for the counting function

x = (x 1 , . . . , x p ) ∈ R p -→ c(x) = p i=1 θ(x i ) (see Section 7.1)
is easier to obtain since c has a "separate" structure; it therefore suffices to calculate the Moreau envelope of the θ function. Actually,

• c λ (x) = 1 2λ x 2 -1 2λ p i=1 [x 2 i -2λ] + , • One element in Prox λ c(x) is p λ (x) ∈ R p , where [p λ (x)] i =    x i if |x i | ≥ √ 2λ, 0 otherwise.
These results were already observed in Example 5.4 of [START_REF] Attouch | Convergence of descent methods for semialgebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel method[END_REF].

As we saw it when considering the relaxed forms of the rank function (in Section 5), what is more useful and interesting for applications is the restricted rank function rank R . The calculations for its Moreau-Yosida approximates and proximal set-valued mappings are a bit more complicate than for the rank itself, of the same vein however. Here is the final and complete result.

Theorem 13. We have, for all A ∈ M m,n (R) of rank r ≥ 1:

(i) (rank R ) λ (A) = 1 2λ A 2 F - 1 2λ r i=1 {σ 2 i (A) -[(σ i (A) -R) + ] 2 -2λ} + . ( 31
)
(ii) One minimizer in [START_REF]Three-way arrays: Rank and uniqueness of trilinear decomposition, with application to arithmetic complexity and statistics[END_REF], i.e. one element in Prox λ (rank R )(A), is provided by

B := U Σ B V with Σ B = diag m,n [σ 1 (B), . . . , σ p (B)],
where

• (U, V ) ∈ O(m, n) A , i.e. U and V are orthogonal matrices such that A = U Σ A V , with Σ A = diag m,n [σ 1 (A), . . . , σ r (A), 0, . . . , 0] (a singular value decomposition of A with σ 1 (A) ≥ • • • ≥ σ r (A) > 0); • If √ 2λ ≥ R σ i (B) :=          R if σ i (A) > 2λ+λ 2 2λ , 0 or R if σ i (A) = 2λ+λ 2 2λ , 0 if σ i (A) < 2λ+λ 2 2λ . • If √ 2λ < R σ i (B) :=              R if σ i (A) > R, σ i (A) if √ 2λ < σ i (A) ≤ R, 0 or σ i (A) if √ 2λ = σ i (A), 0 if σ i (A) < √ 2λ.

Comments

1. As the positive parameter λ is supposed to approach 0 in the proximal approximation process, the second case of (ii) in the theorem above is more important than the first one. The pictures in Figure 2 and Figure 3 show, in the one dimensional case, the behaviour of (rank) λ := c λ and (rank R ) λ = (c R ) λ as λ → 0.

We know from Section 5 that the convex relaxed form of the restricted rank function rank

R is co(rank R ) = ψ R (A) =    1 R A * if A sp ≤ R, +∞ otherwise.
It is interesting to calculate explicitly the Moreau-Yosida approximations (ψ R ) λ of ψ R , and to compare them with those of rank R in Theorem 13. Here we are in a more familiar convex framework, so that calculations are easier to carry out. Since the proximal set-valued mapping Prox λ ψ R is actually single-valued on M m,n (R), we adopt the notation

Prox λ ψ R (A) = {prox λ ψ R (A)}.
Theorem 14. Let U and V be orthogonal matrices such that

A = U Σ A V , with Σ A = diag m,n [σ 1 (A), . . . , σ r (A), 0, . . . , 0] (a singular value decomposition of A with σ 1 (A) ≥ σ 2 (A) ≥ • • • ≥ σ r (A) > 0). We set p R λ (A) = (y 1 , . . . , y p ),
with

y i =          R if σ i (A) ≥ λ R + R, σ i (A) -λ R if λ R ≤ σ i (A) < λ R + R, 0 if σ i (A) < λ R . ( 32
)
Then, the proximal mapping and the Moreau envelope of ψ R are described as following.

(i) [ Proximal mapping] prox λ ψ R (A) = U diag m,n (y 1 , . . . , y p )V ; ( 33
) (ii) [ Moreau envelope]
We define, for t ∈ R,

f i λ R (t) := t 2 -2σ i (A)t + 2 λ R |t|,
and for x = (x 1 , . . . , x p )

f λ R (x) := p i=1 f i (x i ).
Then,

(ψ R ) λ (A) = 1 2λ A 2 F - 1 2λ f λ R [p R λ (A)]. ( 34
)
Moreover, for A such that A sp ≤ R,

(ψ R ) λ (A) = 1 2λ A 2 F - 1 2λ p R λ (A) 2 . ( 35
)
as well as how it compares with (rank R ) λ . It also illustrates the following fact: the convex hull (or closed convex hull) of (rank R ) λ is exactly (ψ R ) λ .

Figure 4: The Moreau-Yosida approximations of the restricted rank and nuclear norm.

The case where R = 1 deserves some additional comments. Recalling that

ψ 1 (A) =    A * if A sp ≤ 1 +∞ otherwise,
we have

(ψ 1 ) λ (A) = 1 2 A 2 F - 1 2 f λ [p 1 λ (A)], (36)
where p 1 λ (A) = (y 1 , . . . , y p ), with

y i =          1 if σ i (A) ≥ λ + 1, σ i (A) -λ if λ ≤ σ i (A) < λ + 1, 0 if σ i (A) < λ. ( 37
)
In short,

y i = [σ i (A) -λ] + -[σ i (A) -(λ + 1)] + for all i = 1, . . . , p, so that prox λ ψ 1 (A) = U diag m,n {([σ i (A) -λ] + -[σ i (A) -(λ + 1)] + ) i }V, (38)
(ψ 1 ) λ (A) = 1 2λ A 2 F - 1 2λ p i=1 f i λ (y i ) = 1 2λ r i=1 σ 2 i (A) - 1 2λ r i=1 f i λ ([σ i (A) -λ] + -[σ i (A) -(λ + 1)] + ). ( 39
)
These formulas ( 38) and ( 39) should be put side by side with the expressions of ( . * ) λ and prox λ ( . * ), such as given in [START_REF] Liu | An implementable proximal point algorithmic framework for nuclear norm minimization[END_REF] for example:

prox λ ( . * )(A) = U diag m,n [σ i (A) -λ] + V, ( 40
) ( . * ) λ (A) = 1 2 A 2 F - 1 2 r i=1 {([σ i (A) -λ] + ) i } 2 . = 1 2 r i=1 σ 2 i (A) - 1 2 r i=1 {[σ i (A) -λ] + } 2 . ( 41
)
As expected, since . * ≤ ψ 1 , one has ( . * ) λ ≤ (ψ 1 ) λ for all λ > 0. Also, for λ small enough, namely for λ ≤ σ r (A),

( . * ) λ (A) = (ψ 1 ) λ (A).
Note however that the convex relaxed form of (rank) λ is not ( . * ) λ ; as said before, to compare the relaxed form of the rank function with . * , as well as their corresponding Moreau-Yosida regularized forms, one has to consider their restricted versions on balls {A| A sp≤R }.

The formulas [START_REF] Lewis | Nonsmooth analysis of singular values[END_REF] and ( 41) are used for designing a proximal point algorithm scheme for nuclear norm minimization ( [START_REF] Liu | An implementable proximal point algorithmic framework for nuclear norm minimization[END_REF]).

The generalized subdifferentials of the rank function

Although the rank function is not differentiable, not even continuous, it is possible to calculate the so-called generalized subdifferentials of it (kinds of generalizations of gradients for differentiable functions). These mathematical objects are now much used in nonsmooth analysis and optimization, for theoretical purposes as well as for algorithmic ones. Dealing with the calculation of the generalized subdifferentials of the rank function directly is very difficult because:

• the behaviour of difference quotients like rank(A + tD) -rank A t when A → A and t > 0 goes to 0 is hard to control, and their limits impossible to determine;

• the singular values involved in the regularized versions of the rank function are themselves nonsmooth.

So, we adopted another strategy: to use the relationship between the counting and the rank function, as in [START_REF] Le | Confexifying the counting function on R p for convexifying the rank function on M m,n (R)[END_REF]. We first calculate the generalized subdifferential(s) of the counting function; then, by using the results of Lewis and Sendov ( [START_REF] Lewis | Nonsmooth analysis of singular values[END_REF], [START_REF] Lewis | Nonsmooth analysis of singular values[END_REF]), we obtain the generalized subdifferential(s) of the rank function. Good news: All the generalized subdifferentials turn out to be equal. Moreover, this (common) generalized subdifferential turns out to be a vector space.

Certainly 0 always belongs to it. This was foreseen by the fact that "Every point is a local minimizer of the rank function" (see Section 4). Finally, with the help of an alternate expression of the common subdifferential as a tensor product of vector spaces, we provide its dimension.

All the material of this section is excerpted from [START_REF] Le | The generalized subdifferentials of the rank function[END_REF].

Definitions and preliminary properties

The definitions and properties of generalized subdifferentials have been developed in several works, beginning with the case of locally Lipschitz functions ( [START_REF] Clarke | Optimisation and nonsmooth analysis[END_REF]). Then, they have been generalized for lower-semicontinuous functions. In this section, we only consider the lower semicontinuous case, since it is that of the rank function. Our context is therefore the following:

• f : E → R ∪ {+∞} is proper, lower-semicontinuous (l.s.c
) and x is a point at which f is finite-valued.

• (E, ., . ) is an Euclidean space, for example R p equipped with the usual inner product or M m,n (R) equipped with ., . .

Definition 1. A vector x * ∈ E is a F-subderivative of f at x if lim inf y→0 f (x + y) -f (x) -x * , y y ≥ 0. ( 42
)
The set of all F-subderivatives of f at x is called the Fréchet subdifferential of f at x, and

denoted as ∂ F f (x). Definition 2. A vector x * ∈ E is a viscosity subderivative of f at x if there exists a C 1 -function g : E → R such that g(x) = x * and f -g attains a local minimum at x.
If, in particular,

g(x) = x * , x -x -σ x -x 2
with some positive constant σ, then x * is called a proximal subgradient of f at x.

The set of all viscosity subderivatives and proximal subgradients of f at x are called the viscosity subdifferential and the proximal subdifferential of f at x and denoted as ∂ V f (x) and

∂ P f (x), respectively.
In a finite dimensional context, the Fréchet and the vicosity subdifferentials coincide. And this common subdifferential is also called "regular subdifferential" in some other works (see [START_REF] Lewis | Nonsmooth analysis of singular values[END_REF], [START_REF] Lewis | Nonsmooth analysis of singular values[END_REF], [START_REF] Rockafellar | Variational analysis[END_REF]). Definition 3. A vector x * ∈ E is a limiting subgradient of f at x if there is a sequence of points

x ν in E approaching x with values f (x ν ) approaching the finite value f (x), and a sequence of

y ν in ∂ F f (x ν ) approaching x * .
The set of all limiting subgradients is called the limiting subdifferential and denoted as

∂ L f (x). Definition 4. The Clarke subdifferential ∂ C f (x) of f at x is the set of all x * ∈ E such that ∀v ∈ E, x * , v ≤ f 0 (x, v) := lim ε↓0 lim sup y ↓ f x t ↓ 0 inf w∈v+εB f (y + tw) -f (y) t , ( 43
)
where B is the unit ball in E and y ↓ f x signifies that y and f (y) converge to x and f (x), respectively.

Remark 9. We can also define the Clarke subdifferential of f at x as the set of x * for which (x * , -1) lie in some appropriate normal cone to epif at (x, f (x)) (see [START_REF] Clarke | Optimisation and nonsmooth analysis[END_REF]). Thus, ∂ C f (x) is closed and convex for every x in E.

Moreover, since we are in a finite dimensional context, we have the next string of inclusions

∂ P f (x) ⊂ ∂ V f (x) = ∂ F f (x) ⊂ ∂ L f (x) ⊂ ∂ C f (x). ( 44
)

The main results

Let us fix some notations:

• As previously, for x ∈ R p , let diag m,n (x) denote an m×n matrix with entries diag m,n (x) i,i =

x i for all i, and diag m,n (x) i,j = 0 for i = j.

• O(n) is the set of orthogonal matrices in M n (R).

Theorem 15. All the generalized subdifferentials (proximal, Fréchet, viscosity, limiting, Clarke) of the rank function coincide. We denote ∂(rank) the common subdifferential set-valued mapping.

For A ∈ M m,n (R), ∂(rank)(A) is constructed as follows:

• Consider the matrices U ∈ O(m) and V ∈ O(n) such that U diag m,n (σ(A))V = A
(in other words, we collect all the orthogonal matrices U and V which give a singular value decomposition of A; elsewhere in this paper, the set of all such (U, V ) is denoted as

O(m, n) A ).
• Consider the "pseudo-diagonal" matrices diag m,n (x * ), where x * ∈ R p is such that x * i = 0 for all i = 1, . . . , r (recall that r = rank A).

• Then, collect all the matrices of the form U diag m,n (x * )V .

In a single formula,

∂(rank)(A) = {U diag m,n (x * )V | U ∈ O(m), V ∈ O(n) such that U diag m,n (σ(A))V = A,
x * i = 0 for all i = 1, . . . , r}.

A further interesting property of ∂(rank)(A) is that it is not only a closed convex set but a vector space.

Proposition 1.

For A ∈ M m,n (R), ∂(rank)(A) is a vector space.

An alternate expression of ∂(rank)(A), a more palpable one, is possible. The subdifferential of the rank function can be also represented as the tensor product of two vector spaces in R m and R n , as indicated in the following proposition.

Proposition 2. Let N (A) and N (A T ) denote the null spaces of matrices A and A T , respectively.

Then

∂(rank)(A) = N (A T ) ⊗ N (A),
where ⊗ denotes the tensor product. In a more detailed form,

N (A T ) ⊗ N (A) = i,j a ij α i β T j | (α i ) is a basis of N (A T ) , (β j ) is a basis of N (A)} .
Consequentely, the dimension of ∂(rank)(A) is (m -r)(n -r), where r is the rank of A.

We now have a good knowledge of ∂(rank); it is also possible to calculate ∂(rank λ ) and link Sections 7 and 8. The formulas obtained serve as an illustration or example of general results linking the generalized subdifferentials of f and of its Moreau envelope, like those displayed in

• [12, Section 5] on the so-called quadratic inf-convolutions (for example, properties at a point M where ∂ P (rank λ )(M ) is nonempty).

• [47, Section 10.32] on the "subsmoothness" of Moreau envelopes.

• [30, Section 5] on the limits (of various kinds) of the generalized subdifferentials of f λ when λ > 0 goes to 0.

Remark 10. In a paper that we came aware of after a first writing of the present paper, Luke Remark 11. The rank function is an example of semi-algebraic function for which every point is a critical point (in the sense that 0 belongs to the (here common) generalized subdifferential of the function at any point). For semi-algebraic extended-real-valued functions defined on R d , the graphs of the generalized subdifferentials are semi-algebraic subsets of R d × R d whose appropriately defined dimensions (local and global) can be defined and studied ( [START_REF] Drusvyatskiy | The dimension of semialgebraic subdifferential graphs[END_REF], [START_REF] Drusvyatskiy | Semi-algebraic functions have small subdifferentials[END_REF]). In 9 Further notions related to the rank: the spark, the rigidity, the cp-rank of a matrix.

The notion of rank of a matrix dates back to the mathematician Sylvester (1814-1897);

there have been since several further concepts whose definitions resemble that of rank, or whose definitions use that of rank. We mention here three of them: the spark, the rigidity and the cp-rank.

The spark of a matrix

Given A ∈ M m,n (R), the spark of A is the smallest positive integer k such that there exists a set of k columns of A which are linearly dependent. Remember that the rank of A is the largest number of columns of A which are linearly independent. The term spark seems to have been coined by Donoho and Elad in 2003 (see [START_REF] Bruckstein | From sparse solutions of systems of equations to sparse modeling of signals and images[END_REF] and [START_REF] Elad | Sparse and redundant representations: From theory to applications in signal and image processing[END_REF] for comments on it and some properties); it is also known as the girth in matroid theory, or even (within one unity) the Kruskal rank (see below).

Actually, the given definition of spark is a bit uncomplete: If A is of full column rank, i.e.

if rank A = n, there is no set of k columns of A which are linearly dependent. In that case, we should adopt +∞ as for the spark of A (the infimum over the empty set).

The other extreme case is when one column of A is a zero-column: then spark A = 1. In short, if A does not contain any zero-column and is not of full column rank, we have

2 ≤ spark A ≤ rank A + 1. ( 45
)
Just for an example, consider

A =      1 1 0 0 0 1 1 1 1      , so that A T =      1 0 1 1 0 1 0 1 1     
.

Here, rankA = rank A T = 2 while spark A = 2. Note that spark A T = 3, so that spark A = spark A T as a general rule.

The spark of A is much more difficult to calculate than the rank of A. There however is a variational formulation of spark A (i.e., as the result of some minimization problem). Indeed, to have columns A i 1 , A i 2 , . . . , A i k to be linearly dependent means that there exists a non-zero

λ = (λ i 1 , λ i 2 , . . . , λ i k ) ∈ R k such that k l=1 λ i l A i l = 0; in other words Aλ = 0 for λ = (0, . . . , λ i 1 , . . . , 0, . . . , λ i k , . . . , 0) ∈ R n .
To summarize in a single formula, spark A = min

Aλ = 0, λ = 0 c(λ). ( 46
)
Problem ( 46) is an example of rank minimization problem (see the Introduction), where the constraint set C := {λ ∈ R n , Aλ = 0, λ = 0} is nonempty when A is not of full column rank. A related notion is due to Kruskal (1977) who introduced it in the context of fitting trilinear arrays by simple models ( [START_REF]Three-way arrays: Rank and uniqueness of trilinear decomposition, with application to arithmetic complexity and statistics[END_REF]); this notion therefore has been called later the Kruskal rank, and denoted as K-rank. By definition, the Kruskal rank of A, K-rank A, is the maximal positive integer k such that any subset of k columns of A are linearly independent. For example, K-rank A = n if A is of full column rank, and K-rank A ≥ 1 except if A = 0 (in which case it is not defined). The relation between K-rank and spark is now clear

: if A = 0 is not of full column rank, spark A = K-rank A + 1. ( 47
)

The matrix rigidity

Suppose m = n. Then, the rigidity of A ∈ M n (R), denoted by R A (k), is the smallest number of entries that need to be changed in order to reduce the rank of A below k. This notion take roots in computer science ( [START_REF] Valiant | Graph theoretic arguments in low-level complexity[END_REF]), and obtaining bounds on rigidity has a number of implications in complexity theory ( [START_REF] Lokam | Spectral methods for matrix rigidity with applications to size-depth trade-offs and communication complexity[END_REF]). However, computing the rigidity of a matrix is in general a NP-hard problem. For A ∈ M n (R) and any integer k, one can check that

R A (k) ≤ (n -k) 2 . ( 48
)
A variational formulation of R A (k) can easily be proposed as a rank minimization problem:

R A (k) = min M ∈A+S k c(M ), (49)
where c(M ) denotes the number of nonzero entries in the matrix M (or the counting number of (m 11 , . . . , m nn )) and

S k = {S ∈ M n (R)| rank S ≤ k}.

The cp-rank of a matrix

Definition and first properties

The cp-rank is defined for a specific class of symmetric matrices called completely positive matrices; the prefix "cp" is precisely there for recalling the wording "completely positive". Before defining completely positive matrices, we recall what copositive matrices are.

A matrix M ∈ S n (R) is said to be copositive if

M x, x ≥ 0 for all x ∈ R n + .
The set of all copositive matrices is a closed convex cone of S n (R) which we denote by C n (R).

A matrix A ∈ S n (R) is then called completely positive if A, M ≥ 0 for all copositive matrices M. ( 50
)
Thus, the set of all completely positive matrices is again a closed convex cone of S n (R); we denote it by CP n (R). Hence, using the terminology of convex analysis, CP n (R) is the (positive) polar or dual cone of C n (R). Expressed in other ways,

CP n (R) = co{xx T | x ∈ R n + } (51) = {BB T | B is a nonnegative matrix}. (52)
A completely positive matrix A necessarily has nonnegative entries and is positive semidefinite.

Definition 5. The smallest number of columns of a nonnegative matrix B in a factorization A = BB T of a completely positive matrix A is called the cp-rank of A and denoted as cp-rank A.

The cp-rank of A is also the minimal number k of summands in a rank 1-decomposition of A as follows:

A = k i=1 b i b T i , with b i ∈ R n + . ( 53
)
For the zero matrix (which indeed lies in the cone CP n (R), we adopt by convention that its cp-rank is 0.

For matrices A ∈ S n (R) which are not in CP n (R), a decomposition like (53) is not possible;

as it is usual in variational analysis, we therefore pose cp-rank A = +∞.

More information about copositive matrices and completely positive matrices can be found in survey papers like [START_REF] Hiriart-Urruty | A variational approach to copositive matrices[END_REF], [START_REF] Dür | Copositive Programming -a survey[END_REF], [START_REF] Bomze | Copositive optimization -recent developments and applications[END_REF].

If one wished to compare the rank and the cp-rank of a matrix, one gets at the following inequality:

rank A ≤ cp-rank A, (54)
for all completely positive matrices, hence for all symmetric matrices. Indeed, in some cases, the cp-rank of A is equal to the rank of A; for example this holds true when n ≤ 3 or when rank A ≤ 2. To give an example where the inequality in (54) is strict, consider

A =        6 3 3 0 3 5 1 3 3 1 5 3 0 3 3 6        ∈ S 4 (R).
Then, rank A = 3, while cp-rank A = 4 ([4, page 140]).

The cp-rank however enjoys some properties similar to those of the rank (see Section 1 and Section 2). Loewy proved that cp-rank A ≤ n 2 4 for every completely positive matrix of order n ≥ 4 whose graph is triangle free ( [START_REF] Drew | Completely positive matrices associated with M -matrices[END_REF]). The fact that the bound n 2 4 was also valid in all other known cases led the authors to wonder whether this holds for every completely positive matrix of order n ≥ 4. Before going into the details of the proof, some observations are in order:

Theorem 16. (i) If

Conjecture (The DJL

• Like for Theorem 4, due to the homogeneity properties of the functions involved, it suffices to prove Theorem 18 for R = 1.

• When a matrix A is positive semidefinite, which is the case for A completely positive, To get at the converse inequality, decompose A as follows:

A * =
A = k i=1 b i b T i , b i ∈ R n + , b i = 0, k ≥ cp-rank A.
We "normalize" the decomposition above as: We complete the decomposition (64) with the zero matrix: Like for the spark of A (cf. 9.1), the cp-rank of A is difficult to determine, much more than the rank of A.

A = k i=1 α i c i c T i ,
A = k i=1 α i c i c T i + (1 - k i=1 α i )0.
Since the cp-rank is a lower-semicontinuous integer-valued function, like the rank function, it is tempting to ask ourselves: what are the generalized subdifferentials of the cp-rank function?

The question was answered in Section 8 for the rank function, we did not succeed in answering it for the cp-rank function.

By way of conclusion

In this survey, we have painted a broad panorama of the main properties of the rank function (and of some related notions) in the context of variational analysis and optimization. The algorithmic aspects of tackling the rank minimization problems concern the convex relaxed forms of the rank function (such as seen in Section 5); there is indeed a huge specific literature devoted to this way of doing. It remains to see if, even if we deal with a very bumpy (but fascinating) function, the variational analysis properties of the rank function itself could be used for considering directly the rank minimization problems and not their relaxed forms.

1 . 2 . 4 . 5 . 6 . 8 .1

 124568 The rank in linear algebra or matricial calculus. The rank from the topological and semi-algebraic viewpoint. Global minimization of the rank. Relaxed forms of the rank function. Surrogates of the rank function. 7. Regularization-approximation of the rank function. Generalized subdifferentials of the rank function. 9. Further notions related to the rank: the spark, the rigidity, the cp-rank of a matrix. The rank in linear algebra or matricial calculus Let p := min(m, n), let rank : M m,n (R) -→ {0, 1, . . . , p} A -→ rank A.

  (a) rank A = rank A T ; rank A = rank (AA T ) = rank (A T A). (b) If the product AB can be done, rank (AB) ≤ min(rank A, rank B) (Sylvester inequality).As a general rule, when the proposed products of matrices can be done,rank (A 1 A 2 . . . A k ) ≤ min i=1,...,k (rank A 1 , rank A 2 , . . . , rank A k ). When m = n, rank (A k ) ≤ rank A.(c) rank A = 0 if and only if A = 0; rank (cA) = rank A for c = 0. (d) |rank A -rank B| ≤ rank (A + B) ≤ rank A + rank B.

Remark 3 .

 3 As stated earlier, if σ(A) = (σ 1 (A), . . . , σ p (A)), the counting function c at σ(A) is exactly the rank of A, c[σ(A)] = rank A.

Remark 5 .

 5 Since Theorem 3 (of Eckart and Young) is a precise and powerful theorem solving the best approximation problem for any matrix to any sublevel-set of the rank function, we thought some time that it would be possible to derive Theorem 4 (on the convex hull of the restricted rank function) from Theorem 3; we did not succeed.

2 .

 2 When A sp = max i=1,...,p σ i (A) ≤ R, both Moreau-Yosida approximates (rank) λ and (rank R ) λ coincide at A.

Figure 2 :

 2 Figure 2: The Moreau-Yosida approximations of the rank.

Figure 3 :

 3 Figure 3: The Moreau-Yosida approximations of the restricted rank. R = 4 here.

(in [ 43 ,

 43 Proposition 3.6]) gives the common expression of various generalized normal cones to S k at A ∈ S k . Due to the existing relationship between generalized subdifferentials and generalized normal cones, via the indicator function, his results are consistent with our Theorem 15.

  the case of the rank function defined on M m,n (R) ≡ R m×n , the dimension of the subdifferential graph (both in the local sense and in the global one) is equal to the dimension of the underlying space, namely d = m × n 7 .

) 9 . 3 . 2 2 . 2 - 1 by 1 . 2 - 1 .

 932221121 A and B are completely positive, cp-rank (A + B) ≤ cp-rank A + cp-rank B. (55) (ii) If A is completely positive and c > 0, cp-rank (cA) = cp-rank A. (56) (iii) The cp-rank is a lower-semicontinuous function: If (A ν ) ν is a sequence of completely positive matrices converging to A, then A is completely positive and lim inf ν→+∞ cp-rank A ν ≥ cp-rank A. (57Upper bounds for the cp-rank. One of the most interesting challenges about the cp-rank is to provide an upper bound of it in terms of rank. In 1983, Hannay and Laffey ([21]) showed that the cp-rank of a completely positive matrix of rank k was necessarily less than or equal to k(k+1) Later, in 2003, this upper bound was sharpened to k(k+1) Barioli and Berman ([7]); they also proved that this upper bound could be achieved by a completely positive matrix of rank k. In short, Theorem 17. (i) For every completely positive matrix of rank k ≥ 2, cp-rank A ≤ k(k + 1) 2 -(ii) For every integer k ≥ 2, there exists a completely positive matrix A whose rank is k and whose cp-rank is k(k+1) By Theorem 17, if A is a n × n completely positive matrix, then cp-rank A ≤ n(na better upper bound on the cp-rank of n × n matrices? Drew, Johnson and

4 . ( 59 ) 9 . 3 . 3 Theorem 18 . 1 R

 459933181 Conjecture) If A is an n × n completely positive matrix, with n ≥ 4, then cp-rank A ≤ n 2 The conjecture was proved by Berman and Shaked-Monderer ([3]) for matrices whose comparison matrices are M-matrices and by Loewy and Tam ([41]) for 5 × 5 matrices whose graph is not complete. But then, for the first time, Barioli announced (in [6]) an example of 7 × 7 completely positive matrix of rank 5 and cp-rank 14. Such matrix is a counter-example to the DJL Conjecture. The convex relaxed form of the cp-rank.Due to properties like (56), it is easy to see that the convex hull of the cp-rank function onCP n (R) is the zero function. So, like for the rank function (in Section 5), we restrict it to some appropriate ball. For R > 0, considerA ∈ S n (R) → cp-rank R (A) :=    cp-rank A if A iscompletely positive and A * ≤ R; the following result, to be put in parallel with Theorem 4. We have A ∈ S n (R) → co(cp-rank R )(A) := A * if A is completely positive and A * ≤ R;

  sum of the eigenvalues of A = tr A. (62) Proof. (of Theorem 18)Let φ denote the function defined in the right-hand side of (61) for R = 1, proposed as the convex hull of the function cp-rank 1 . The domain of the function cp-rank 1 , i.e. the set of matrices where it is finite-valued, and that of φ are equal: it is the compact convex setCP n (R) ∩ {A ∈ S n (R)| A * ≤ 1}.Let therefore now A be chosen completely positive, with A * ≤ 1. Firstly, if A is the zero matrix, it is clear that co(cp-rank 1 )(0) = 0 = φ(0).Secondly, let us assume that A = 0. We have to prove that(cp-rank 1 )(A) = A * . Because cp-rank ≥ rank ≥ . * on {A ∈ S n (R)| A * ≤ 1},we get a first inequality co(cp-rank 1 )(A) ≥ A * .(63)

with c i := b i b i and α i := b i 2 .α 2 =α i ≤ 1 (

 221 i tr(c i c T i ) = k i=1 α i (because tr(c i c T i ) = c i because tr A = A * ≤ 1).

37 Thus,

 37 due to the convexity of co(cp-rank 1 ),co(cp-rank 1 )(A) ≤ k i=1 α i cp-rank 1 (c i c T i ) + (1 -k i=1 α i )cp-rank 1 (0) ≤ k i=1 α i (because cp-rank 1 (c i c T i ) = 1 and cp-rank 1 (0) = 0) ≤ tr A = A * . (65)Putting together the inequalities (63) and (65) yield the desired result.Remark 12. As for the rank function, there is no difference between co(rank R ) and its closure co(cp-rank R ).Remark 13. On the set where both co(rank R ) and co(cp-rank R ) are finite, i.e. on CP n (R)∩{A ∈ S n (R)| A * ≤ R}, they are equal and linear.

D "pseudo-diagonal" means that dij = 0 for i = j. One also uses the notation diag m,n [σ1, . . . , σp] for D.

The notation for the SVD is a bit nonstandard. Usually, a singular value decomposition takes the form U ΣV T (i.e., V T instead of V ). But there is no difference from the mathematical point of view.

A nonempty closed set S is an Euclidean space E is called prox-regular at a point x ∈ S if the projection set PS(x) of x on S is single-valued around x.

J.-J. Moreau presented this approximation process as an example of inf-convolution or epigraphic addition of two functions; it was in 1962, exactly 50 years ago. This date marks the birth of modern convex analysis and optimization.

D.Drusvyatskiy, personal communication (August 2012).

Acknowledgments. We thank the various readers or referees who, after the first circulation of this survey paper (July 2012), provided comments, improvements and additional references.


Hung Quoc

Lu

Soheib Fergani
email: sfergani@laas.fr

Carine Jauberthie

Fault detection combining adaptive degrees of freedom χ 2 -statistics and interval approach for nonlinear systems

Keywords: Adaptive Innovation-based method, Fault diagnosis, Adaptive amplifier coefficient, Automotive quarter model

An enlargement of the Adaptive degrees of freedom χ 2 -statistics (ADFC) method to fault detection for nonlinear systems with mixed uncertainties (stochastic and bounded uncertainties) is presented in this paper. The ADFC approach, primarily developped for fault detection in case of linear systems, is then combined with the Reinforced likelihood box particle filter (RLBPF). A residual generator is used, followed by the adaptive amplifier coefficient (a.a.c.) concept in the decision making stage. Then, the proposed approach is applied to a nonlinear Magneto-Rheological damper model to illustrate the efficiency of the method.

INTRODUCTION

In modern control system, many important issues are required to be solved such as the availability, cost efficiency, reliability, operating safety, environmental protection,... Therefore, fault diagnosis becomes indispensable and one of the main solutions for these requirements. A fault must be diagnosed as early as possible, even if it may be tolerable, to prevent any serious consequence. The fault diagnosis consists in three tasks: fault detection (FD), fault isolation and fault estimation (or fault identification), in which FD is the first stage and building block of any further task. For all the above reasons, this paper is focusing on the developement of a reliable FD strategy.

Model-based FD approach applied to dynamical systems using analytical redundancy is widely used. It consists in two stages: residual generation and residual evaluation or decision making. The residual is normally a difference between a measured signal and an estimate of the later. It should be zero-valued in the fault-free case and diverge from zero in the faulty case. In the second stage, the residual (or its transformation via a function) is compared to a constant threshold or an adaptive threshold depending on the fault detection method. Especially, using observerbased concept, a filter or an observer is the core of the residual generator which provides the estimate quantity of interest to deduce the residual. This is the structure of the FD method proposed in the paper (Fig. 1). The most advantage of this structure is that it can easily combine filter/observer and decision making method (pre-existent or to be developed) for adapting to various contexts and purposes.

The challenge of the present work is twofold: the non linearity of the system and the mixed uncertainties context (bounded and stochastic uncertainties). Indeed, in the literature, most of researches about FD are devoted to linear system while fewer articles provide solution for non linear counterpart. Nonlinear systems are usually assumed to be smooth "enough" with mild conditions for a nice linearization, then techniques for the linear case can be applied. For instance, the Standard Kalman Filter (SKF) [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF] and its extensions (included Extended Kalman Filter (EKF), Unscented Kalman Filter (UKF), Robust Kalman fitering, Interval Kalman Filter (IKF), Minimum Upper Bound of Variance Interval Kalman Filter (UBIKF), Optimal upper bounded Interval Kalman Filter (OUBIKF)) have been applied in residual generators of different FD schemes/methods. We refer to [START_REF] Zhong | A survey on model-based fault diagnosis for linear discrete time-varying systems[END_REF] which provides a survey on model-based fault diagnosis for linear discrete time-varying (LTV) systems with useful references concerning FD using SKF, EKF and UKF. Otherwise, [START_REF] Xiong | Fault detection using interval Kalman filtering enhanced by constraint propagation[END_REF] proposed FD method using IKF, [START_REF] Tran | Cadre unifié pour la modélisation des incertitudes statistiques et bornées -Application à la détection et isolation de défauts dans les systèmes dynamiques incertains par estimation[END_REF] with FD using UBIKF, [START_REF] Lu | A new scheme for fault detection based on Optimal Upper Bounded Interval Kalman Filter[END_REF] with FD based on OUBIKF, while output observer (non interval case) for FD in linear systems can be found in [START_REF] Mohamadi | Output observer for fault detection in linear systems[END_REF]. Furthermore, also for linear system, setmembership methods in fault diagnosis are concerned in [START_REF] Puig | Fault diagnosis and fault tolerant control using set-membership approaches: Application to real case studies[END_REF] while interval observer approach for FD is investigated, for instance, in [START_REF] Raka | Fault detection based on robust adaptive thresholds: A dynamic interval approach[END_REF]; [START_REF] Chevet | Robust sensor fault detection for linear parameter-varying systems using interval observer[END_REF]. Some of these works [START_REF] Xiong | Fault detection using interval Kalman filtering enhanced by constraint propagation[END_REF][START_REF] Tran | Cadre unifié pour la modélisation des incertitudes statistiques et bornées -Application à la détection et isolation de défauts dans les systèmes dynamiques incertains par estimation[END_REF][START_REF] Lu | A new scheme for fault detection based on Optimal Upper Bounded Interval Kalman Filter[END_REF] deal with mixed un-certainties context while the remaining deal either with bounded or stochastic uncertainties only.

In the present paper, an enlargement of the Adaptive degrees of freedom χ 2 -statistics (ADFC) method [START_REF] Lu | A new scheme for fault detection based on Optimal Upper Bounded Interval Kalman Filter[END_REF] to FD for non linear system is proposed. The method combines the Reinforced likelihood box particle filter (RLBPF) developed in [START_REF] Lu | Reinforced likelihood box particle filter[END_REF] (versus OUBIKF for linear system) together with the adaptive amplifier coefficient (a.a.c.) concept in the decision making stage. In doing so, the twofold challenge of the issue is met. Interval analysis is a primary tool for computation of the method, no strict condition is required for (non linear) dynamic functions of the system.

The paper is organized as follows. Preliminaries are provided in Section 2, including used notations and essential of interval analysis. Section 3 presents main results with the extended FD method. Simulations based on the nonlinear Magneto-Rheological damper model are provided in Section 4. Finally, Section 5 provides conlusions and perspectives.

PRELIMINARY

Essential of interval analysis

A real interval [x] = [x, x] = {t ∈ R | x ≤ t ≤ x} is a closed connected subset of R and characterized by two extreme values x ≤ x. A real interval matrix [X] = ([x ij ]) of dimension p × q, also called an interval in R p×q , is a matrix with real interval components [x ij ], i ∈ {1, ..., p}, j ∈ {1, ..., q}. Write X ∈ [X] to indicate a point matrix X = (x ij ) belonging element-wise to [X]. Define: • sup([X]) △ = (sup([x ij ])) ≡ X, • inf([X]) △ = (inf([x ij ])) ≡ X, • mid([X]) △ = (X + X)/2 , • rad([X]) △ = (X -X)/2 , • width([X]) △ = X -X ,
which are called respectively (resp.) the largest matrix, the smallest matrix, the midpoint matrix, the radius matrix and the width matrix of [X]. Define the hull of a closed set S and the hull of two interval [X 1 ], [X 2 ] of the same dimension as follows: hull{∅}

△ = ∅, hull{S} △ = [inf(S), sup(S)] and hull{[X 1 ], [X 2 ]} △ = [inf{X 1 , X 2 }, sup{X 1 , X 2 }].
We also write an interval as:

[X] = [X, X] = mid([X]) ± rad([X]).
Basic interval operators ⋄ ∈ {+, -, ×, ÷} defined in [START_REF] Jaulin | Applied Interval Analysis, with Examples in Parameter and State Estimation, Robust Control and Robotics[END_REF] can be used to compute directly all operations [u] ⋄ [v] and α ⋄ [u], for real intervals [u], [v] and α ∈ R, without any further approximation algorithm. Then, interval matrix computations are defined similarly to matrix computations using the basic operators while more general operators are constructed by means of inclusion function [f ] defined in Definition 1 below. In practice, the package Intlab [START_REF] Rump | INTLAB -INTerval LABoratory[END_REF] developed for Matlab (also existing in Octave and C/C++) is used for computations. Definition 1. Let f be any function from

D ⊆ R m to R n . An inclusion function of f is a function [f ] that maps any interval [x] in D to an interval [f ]([x]) in R n containing the image set f ([x]). Symbolically, [f ] : [x] → [f ]([x]) ⊃ f ([x]), ∀[x] ⊂ D, with [f ]([x]) is an interval in R n . The minimal inclusion function of f is a function [f ] * so that for every [x] ⊂ D, [f ] * ([x]) is the smallest interval containing f ([x]). □
If f is composed by a finite number of operators (+, -, ×, ÷) and elementary functions (exp, sin, cos, (.), ...) then a natural inclusion function of f is a function [f ] having the same expression of the former in which point variables are replaced by corresponding interval variables while operators and functions are replaced by their interval counterpart. There is other commonly used inclusion function such as centred inclusion function, mixed centred inclusion function, Taylor inclusion function,... (see [START_REF] Jaulin | Applied Interval Analysis, with Examples in Parameter and State Estimation, Robust Control and Robotics[END_REF] for more details).

Other notations

A positive semidefinite matrix X is denoted by

X ⪰ 0. Let M, N ∈ R n×n , define N ⪯ M if M -N ⪰ 0, M is called an upper bound of N and N a lower bound of M . Let ∅ ̸ = Ω ⊂ R n×n , M is an upper bound (resp. lower bound) of Ω, denoted Ω ⪯ M (resp. M ⪯ Ω) , if X ⪯ M (resp. M ⪯ X), ∀X ∈ Ω. Denote also: S + ([X]) △ = X ∈ [X] : X = X T and X ⪰ 0 , I(x) = 1 if the conditions x are true and I(x) = 0 otherwise, mean(x) = n i=1 x i /n, ∀x ∈ R n .
The notation p : l : q (p ≤ q) is used for a range from p to q with step l. For l = 1, we write p : q. A sequence of variables can be noted interchangeably as w 1 , ..., w k or w 1 : w k or w 1:k .

MAIN RESULTS

Model and problem formulation

Consider the following system (1) with additive measurement noises : (Σ) :

x k = f k (x k-1 , u k , w k ) , y k = h k (x k , u k ) + v k + f s k , k ∈ N * , (1)
where x k ∈ R nx and y k ∈ R ny represent state variables and measures respectively,

u k ∈ R nu inputs, w k ∈ R nx state disturbances/noises, v k ∈ R ny measurement noises, f s k ∈ R ny
additive sensor fault vectors. Each of its components corresponds to a sensor fault. Thus, the fault vector f s k can be of the multiple or single error type.

Assumptions (A) (A1) Given [x 0 ], [u k ], [w k ], [y k ], [µ k ], f k , h k , k ∈ N * . (A2) v k ∼ N (µ k , Σ k ) with unknown µ k and Σ k . (A3) x 0 ∈ [x 0 ], u k ∈ [u k ], w k ∈ [w k ], y k ∈ [y k ], µ k ∈ [µ k ] and Σ k ∈ [Σ k ]. (A4) [Σ k ] is unknown.
Remark 2. The measurements are given as intervals [y k ] due to the sensor precision. In many theoretical problems and applications, it is assumed that µ k ≡ [µ k ] ≡ 0 for simplicity. It is natural that [Σ k ] is unknown and the mention of this fact is necessary for the subsequent development. Assumptions (A3) are implicitly deduced that the belonging of unknown real (true) terms in their corresponding intervals are certain. This certainty can be understood as "with probability 1 ". System (1) with assumptions (A) can be adapted to a wide range of applications. The system is time varying, the uncertainties come from system disturbances or measurement noises with different kinds (stochastic/bounded uncertainties). □

To evaluate the fault detection performance, four indicators are introduced. Assume that system (1) is applied for N iterations among which faults occur in a time inteinte R with length l (0 ≤ l ≤ N ). Actually, R may be an interval or union of intervals but is called hereafter an error range for simplicity. Then, a detection signal is denoted by

π k which has values in {0; 1}. A signal π k = 1 is called a right detected signal if k ∈ R and a false detected signal if k ̸ ∈ R.
Then, the four indicators are defined as follows:

+ Detection Rate: DR = k∈R

I(π k =1) l × 100%, + False Alarm Rate: FAR = k̸ ∈R I(π k =1)
N -l × 100%, + Efficiency: EFF = DR -FAR, + No detection rate: NDR = 100% -DR.

Adaptive approach combining RLBPF to fault detection

In this section, the ADFC method is extended to a more general framework, including nonlinear dynamical systems. This approach has also the ability to combine with several set-membership estimation filters (w.r.t to some conditions in see Theorem 3). In the following, the proposed approach is applied with RLBPF and the extension is based on Theorem 3. Theorem 3. Consider system (Σ) with assumptions (A) and fault free case (f s k = 0, ∀k ∈ N * ). Let r k = y k -h k (x k , u k ) -µ k be the residual, [x k|k ] be interval estimates of an interval filter and

[x k+1|k ] △ = [f k+1 ] [x k|k ], [u k+1 ], [w k+1 ] , [r k+1 ] △ = [y k+1 ] -[h k+1 ] [x k+1|k ], [u k+1 ] -[µ k+1 ].
Then, for all k ≥ 1:

(a) y k |x k ∼ N (m k , Σ k ) where m k = h k (x k , u ) + µ k and
the notation y k |x k means for the stochastic variables y k given (or knowing) x k . When w k is a bounded disturbance, then x k is deterministic and

y k |x k ≡ y k . (b) r k ∼ N (0, Σ k ) where Σ k ∈ [Σ k ].
Consequently, a confidence interval of (more than) 99.7% of all admissible r k can be deduced as

CI(r k ) = ± 3. diag(Σ k ), with diag(A) = (a 11 , ..., a nn ) T , ∀A = (a ij ) ∈ R n×n .
Assuming further that x k ∈ [x k|k ] at every k ≥ 1, then:

(c) x k+1 ∈ [x k+1|k
] with certainty (probability 1), (d) r k+1 ∈ [r k+1 ] with certainty (probability 1). □

Proof. Knowing x k in either cases of w k (disturbance or stochastic noise), the term

h k (x k , u k ) is known (or deterministic). So, since v k ∼ N (µ k , Σ k ) then y k |x k = h k (x k , u k ) + v k is normally distributed with mean m k = h k (x k , u k ) + µ k and covariance Σ k .
Also, by the property of normal distribution, it concludes that r k ∼ N (0, Σ k ). So, statements (a) and (b) are proved.

To prove the statements (c) and (d), the inclusion function property is used, that is:

f (x) ∈ [f ]([x]), ∀x ∈ [x]
with certainty. By assumptions (A3) and x k ∈ [x k|k ], then

• f k+1 (x k , u k+1 , w k+1 ) = x k+1 ∈ [x k+1|k ], • h k+1 (x k+1 , u k+1 ) ∈ [h k+1 ]([x k+1|k ], [u k+1 ]), • µ k+1 ∈ [µ
(%) = N k=1 I(x k ∈ [x k|k ])/N × 100% measuring the percentage of x k ∈ [x k|k
] greater than some level (e.g. 80% or 90%) can be suited to the FD procedure of the ADFC method, although it may cause certain false alarms. □

Principle of the method

Since r k ∼ N (0, Σ k ) then ξ k △ = r T k Σ -1 k r k is χ 2 -distributed with n y degrees of freedom (d.f.) (ξ k ∼ χ 2 (n y )). Thus, for any a k > 0 so that S + ([Σ k ]) ⪯ a k I ny , ξ k ≥ r T k a -1 k r k = r T k r k /a k ≡ ξk . However, r k is unknown and r k ∈ [r k ] with a computable [r k ]
. So, we aim at using the statistic

U k = sup{abs([r k ] T [r k ]/a k )} = sup{abs([ ξk ])},
(2) as proposed in [START_REF] Lu | A new scheme for fault detection based on Optimal Upper Bounded Interval Kalman Filter[END_REF] with the a k aforementioned. The statistic U k is approximated by a χ 2 -distributed random variable with an appropriate d.f. κ k n y by determining a k and using an a.a.c. κ k .

Explication.

Theoretically, for ensuring S + ([Σ k ]) ⪯ a k I ny , it is required that max{diag(Σ k )} < a k [START_REF] Lu | Optimally bounded interval kalman filter[END_REF]. The more a k is large, the more ξk is small and vice versa. Besides, note that r k ∈ [r k ] is with probability 1 (Theorem 3), where

r k ∼ N (0, Σ k ), Σ k ∈ [Σ k ] . This implies [r k ] ⊃ CI(r k )
and [r k ] must be greater considerably than the confidence interval in order to contain r k with such a certainty. Therefore, an appropriate (optimal) a k must ensure a compromise between [r k ] and CI(r k ), says a function of [r k ] and [Σ k ], provided that:

max{diag(Σ k )} < a k = ϕ k ([r k ], [Σ k ]) < max{r k }. (3)
In the development of this section, a simple proposition for a k is that

a k = λ 1 • mean{width([r k ])} so that S + ([Σ k ]) ⪯ a k I, (4
) where λ 1 > 0 is a scaling factor changing depending on the considered application.

The a.a.c. κ k is involved with the threshold δ k to which U k is compared: P χ 2 (κ k n y ) > δ k = α, where α is a given significant level. So, U k ≈ χ 2 (κ k n y ) ≤ δ k with confident level 1 -α in the fault free case. κ k is chosen [START_REF] Lu | A new scheme for fault detection based on Optimal Upper Bounded Interval Kalman Filter[END_REF] as

κ k = mean{width([r k ])}. ( 5
)
so that it is sensitive to the fault, large enough to obtain a small FAR (e.g. ≤ 5%) in the fault free case, sensitively affected by width([r k ]) as well as U k but it does not increase as fast as the later whenever a fault occurs and affects on the width([r k ]). This choice also provided a good performance fault detection with FAR < 3% in several scenarios of simulations presented in the referenced paper. Thus, in the present work, it is raised naturally that 6) is chosen for adapting to different applications and purposes by using different level parameters λ 2 > 0.

κ k = λ 2 • mean{width([r k ])} (
Thanks to the use of RLBPF with interval analysis properties, the enlarged method is provided without assumptions about dynamic functions f k and h k nor any linearization is required. Remark 6. In linear systems, residuals are determined by r k = y k -h k (x k , u k )-µ k where xk is an estimate of x k and h k is a linear function. Thus, r k can be computed explicitly at every time step as well as its distribution under the SKF assumptions. Consequently, with additional bounded uncertainties, the covariance of [r k ] is well determined as some computable matrix [Σ k ] and a k has a more accurate choice using Theorem 2 in [START_REF] Lu | A new scheme for fault detection based on Optimal Upper Bounded Interval Kalman Filter[END_REF]

so that S + ([Σ k ]) ⪯ a k I ny . Evidently, if [Σ k
] is assumed to be known in assumption (A4), a k can be computed in the same manner. In this work, we consider

r k = y k - h k (x k , u k ) -µ k , although it is unknown, and compute [r k ] ensuring that r k ∈ [r k ] with certainty. □
The ADFC method is eventually enlarged to a more general framework of non linear system, as presented above. It can be summarized by the following algorithm.

Algorithm 1 ADFC method to nonlinear system

1: Initialization: α, λ 1 , λ 2 , [x 0|0 ] ≡ [x 0 ], f k , h k , [u k ], [w k ], [y k ], [µ k ], [Σ k ], k = 1, 2, ..., N . 2: for k = 1, 2, 3, ...N do 3: Use RLBPF to get: [x k-1|k-1 ] 4: [x k|k-1 ] = [f k ]([x k-1|k-1 ], [u k ], [w k ]) 5: [r k ] = [y k ] -[h k ]([x k|k-1 ], [u k ]) -[µ k ] 6: a k = λ 1 • mean{width([r k ])} 7: U k = sup{abs([r k ] T [r k ]/a k )} 8: κ k = λ 2 • mean{width([r k ])} 9: Find δ k s.t.: P(χ 2 (κ k n y ) > δ k ) = α. 10: Detection signal : π k = I(U k > δ k ). 11: end for 4. APPLICATION

Quarter vehicle model

The vertical quarter car model is often used to study the vertical behavior of a vehicle according to the suspension characteristic (passive or controlled) (Fig. 2). When controlled suspension is considered, the passive damper F c is removed and replaced by an actuator that provides a force u either active or semi-active depending on the chosen actuator (Fig. 3). In figures 2, the sprung mass m s and unsprung mass m us represent respectively the vehicle chassis and the vehicle wheel. z s and z us are respectively the relative vertical displacement of the vehicle chassis and the vehicle wheel with respect to the road. z r is considered as the road disturbance. Vertical efforts generated by the suspension and tire elements are nonlinear. Let recall that:

F tz = k t (z us -z r ) + c t ( żus -żr ) F sz = F k (z s -z us ) + F c ( żs -żus ) (passive suspension) F sz = F k (z s -z us ) + u (controlled suspension) (7)
where k t and c t are the linear tire stiffness and damping factors, F tz the tire force usually assumed to be linear and F sz the suspension force. The vertical quarter car model is given by the following dynamical equations,

m s zs = -(F sz + F dz ) m us zus = F sz -F tz (8)
where

• z def = (z s -z us ) is the suspension deflection,
• z s and z us are the chassis and unsprung masses bounce,

• m s and m us are sprung and unsprung masses,

• F k (.) is a nonlinear function of z def ,

• F c (.) is a nonlinear function of żdef ,

• F dz describes a vertical disturbance force (that can be caused by a load transfer, e.g. steering situation).

Then, according to the suspension model chosen, different kinds of quarter car models may be obtained:

• If u = F c ( żdef ), the suspension is passive. • If u = F c ( żdef , Ω), the suspension is semi-active,
where Ω is input parameter of the controlled damper that modifies the damping factor. • If u is an independent function, the quarter car is said to be active. Remark 7. In the vertical quarter vehicle model, the nonlinear phenomena come from the force description of the suspension elements and not from the equation structure. Therefore, the model can be set as a LPV system. The unsprung mass m us corresponds to the set of elements that compose the wheel, the suspension system and multiple links from the chassis to the"road". Without loss of generality, it is often referred to as the wheel since z us is the center of the wheel. □

Simulation

Consider the following nonlinear system modeling the MR (Magneto-Rheological) damper:

m s zs = -k s z def -F damper m us zus = k s z def + F damper -k t (z us -z r ) , (9)
F damper = c 0 żdef + k 0 z def + f I tanh (c 1 żdef + k 1 z def ) ,
where c 0 , k 0 , c 1 , k 1 are constant chosen according to [START_REF] Nino-Juarez | Minimizing the frequency effect in a black box model of a magneto-rheological damper[END_REF] such that c 0 = 1500 (Nsm -1 ), c 1 = 129 (sm -1 ),

k 0 = 989 (Nm -1 ), k 1 = 85 (m -1
), and f I is a controllable force depending on the input current I and satisfying the dissipativity constraint 0 < f min ≤ f I ≤ f max . In this simulation, we consider f min = 1000 N/m and f max = 1500 N/m. Other parameter values used in the simulation are presented in Table 1 issued from [START_REF] Fergani | Robust multivariable control for vehicle dynamics[END_REF]. Comparing to the general system (8), in the MR damper model ( 9), assume that F dz = 0 and F tz = k t (z us -z r ).

Putting

• x = [z s , żs , z us , żus ] T as state variable under consideration and x(i), i ∈ {1, ..., 4}, are its components, • u = f I as controllable input, • w = z r , then x, u, w are functions of time t and the state-space representation of ( 9) is expressed in the form

ẋt = f (t, x t ) =    f 1 (t, x t ) f 2 (t, x t ) f 3 (t, x t ) f 4 (t, x t )    , (10)
where

f 1 (t, x t ) = e T 2 x t , f 3 (t, x t ) = e T 4 x t , f 2 (t, x t ) = a T x t -u t tanh(b T x t ) /m s , f 4 (t, x t ) = c T x t + u t tanh(b T x t ) + k t w t /m us
, with e i 's are i-th column vectors of the corresponding identity matrix and

a =    -k s -k 0 -c 0 k s + k 0 c 0   , b =    k 1 c 1 -k 1 -c 1   , c =    k s + k 0 c 0 -k s -k 0 -k t -c 0   .
The system (10) will be discretized using the Fourth order Runge-Kutta method [START_REF] Kincaid | Numerical Analysis[END_REF] with a chosen sampling time T = 10 -4 (s). The resulted discrete time state dynamical system is denoted by:

x k = f (x k-1 , u k , w k ) + η k , k ∈ N * , (11)
where η k is assumed to be Gaussian noise with zero mean and covariances 10 -8 I nx . The corresponding observed measurements are assumed to be z def at every time step, thus the measurement dynamical equation can be expressed in the form

y k = h(x k ) + v k = Cx k + v k , C = [1, 0, -1, 0] , (12)
where v k is assumed to be Gaussian with mean

µ k ∈ [µ k ] = [-0.005, 0.005] and variance σ 2 k ∈ [σ 2 k ] = [1, 4] * 10 -6
. The precision of the sensors is assumed to be ±0.005 (m).

State and measurement simulation: Assume that the initial state is x 0 = (0, 0, 0, 0) T , the control force input is set to get its maximum value constantly (u = 1500) for all time instants and the road disturbance is set as w = 0.05 max{0, sin(πt)}. {x k , y k } k=1:N are then generated using ( 11) and ( 12) for N = 5000 steps. The measurements obtained will be intervals [y k ] = y k ± 0.005 because of sensor errors.

For the fault detection purpose, a fault of an amplitude of b = 0.02m is added to the simulated observed measurements in a range R with length l = 200 and the following choices are used: λ 1 = 0.02, λ 2 = 10, α = 0.03.

A simulation result is figured out in Figures 4 -6. In the error range, the residual deviate from 0 (Fig. 4) and most of the statistics U k (blue line) passe over the adaptive thresholds δ k (red line) (Fig. 5). 2). For a fault value b = 0.02m, the efficiency index (EFF) is about 66%. Table 2. ADFC method -Fault detection to a Quarter vehicle model.

CONCLUSION AND PERSPECTIVE

An adaptive approach of sensor fault detection applied to nonlinear discrete time dynamical system is proposed. The approach combines RLBPF with a hypothesis testing method using χ 2 -statistics with adaptive degrees of freedom. Theoretical framework is developed thanks to interval analysis. A great flexibility of adjusting several factors (parameters) makes the approach highly fitted to multiple applications.

Simulations are performed using the nonlinear Magneto-Rheological damper model. The results show the efficiency of the proposed method providing good performances w.r.t the fault magnitudes.

The method is developed however in the framework of (additive) sensor fault systems. Extend this method to deal with other kinds of fault (e.g. actuator faults) and with fault identification is a perspective of our future research.

Fig. 1 .

 1 Fig. 1. General structure of observer-based Fault detection.

  k+1 ] with certainty. Thus, it implies that r k+1 ∈ [r k+1 ] with certainty ensured by the inclusion function property. □ Remark 4. The notation [x k|k ], k ≥ 1, denotes the interval estimate provided by the corresponding interval filter at the end of the time instant k. [x k+1|k ] denotes the propagation box at the next iteration that allows to obtain [x k+1|k+1 ]. □ Remark 5. The assumption x k ∈ [x k|k ] at every k ≥ 1 is strong and related to the performance and convergence of the filter. Practically, an interval filter that has the O

Fig. 2 .

 2 Fig. 2. Quarter vehicle model

Fig. 4 .

 4 Fig. 4. ADFC method -Residual [r k ] for a Quarter vehicle model with sensor fault.

Fig. 6 .

 6 Fig. 6. ADFC method -Detection signal for a Quarter vehicle model. Then the fault detection procedure is replicated for L = 100 times where the error range R is chosen randomly and indicators DR, NDR, FAR, EFF are yielded as their corresponding means after L times of simulations (Table2). For a fault value b = 0.02m, the efficiency index (EFF) is about 66%.

Table 1 .

 1 Linearized Renault Mégane Coupé parameters of the quarter vertical model (front suspension).

	Symbol	Value	Unit	Signification
	ms	315	kg	sprung mass
	mus	37.5	kg	unsprung mass
	ks	29500	N/m suspension linearized stiffness
	kt	208000	N/m tire stiffness
	z def	[-0.09; 0.05] m	suspension bound (stroke limit)