Laurent Bako
email: laurent.bako@ec-lyon.fr

Subspace clustering through parametric representation and sparse optimization

Keywords: Subspace Clustering, Sparse Optimization, Subspace Arrangement

We consider the problem of recovering a finite number of linear subspaces from a collection of unlabeled data points that lie in the union of the subspaces. The data are such that it is not known which data point originates from which subspace. To address this challenge, we show that the clustering problem is amenable to a sparse optimization problem. Considering a candidate subspace and the distances of the data points to that subspace, the foundation of the proposed method lies in the maximization of the number of zero distances. This can be relaxed into a convex optimization. Efficiency of the relaxation can be significantly increased by solving a sequence of reweighted convex optimization problems.

I. INTRODUCTION

We consider the problem of estimating a finite number s of linear subspaces {S k } s k=1 of the euclidean vector space R n , with equal dimension d. More precisely, assume we are given a finite set {y i } N i=1 of noise-free data points drawn from the subspace arrangement S 1 ∪ • • • ∪ S s . Then the problem of interest in this paper is to infer the individual subspaces S 1 , . . . , S s from the collected data. This can indeed be viewed as a clustering problem because the available data are mixed in the sense that we do not know which data point originates from which subspace.

The literature on subspace clustering/segmentation contains a variety of methods that can be roughly divided into four categories: algebraic methods [START_REF] Vidal | Generalized principal component analysis (GPCA)[END_REF], [START_REF] Ma | Estimation of subspace arrangements with applications in modeling and segmenting mixed data[END_REF], [START_REF] Kanatani | Motion segmentation by subspace separation and model selection[END_REF], statistical methods [START_REF] Fischler | RANSAC random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography[END_REF], matrix factorization-based methods [START_REF] Costeira | A multibody factorization method for independently moving objects[END_REF], [START_REF] Gruber | Multibody factorization with uncertainty and missing data using the EM algorithm[END_REF], sparse representation-based methods [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF], [START_REF] Liu | Exact subspace segmentation and outlier detection by low-rank representation[END_REF], [START_REF] Liu | Robust recovery of subspace structures by low-rank representation[END_REF]. For a detailed exposition of some of those methods we refer to the recent survey [START_REF] Vidal | Subspace clustering[END_REF]. The more recent trend in the treatment of the subspace segmentation problem builds on some ideas borrowed from the field of compressed sensing [START_REF] Donoho | Compressed sensing[END_REF], [START_REF] Candès | Reflections on compressed sensing[END_REF]. The subspace segmentation methods reported in [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF] and [START_REF] Liu | Exact subspace segmentation and outlier detection by low-rank representation[END_REF] allow for an elegant retrieval of the subspaces, their dimensions and their number by solving a convex optimization problem and then applying spectral clustering. However, these methods assume that the subspaces to be estimated are independent, a requirement that imposes a constraint on the number and the dimensions of the subspaces. Another relevant method is the one described in [START_REF] Adler | Probabilistic subspace clustering via sparse representations[END_REF] which relies on sparse representation over a fixed but optimized dictionary in a probabilistic framework.

In this note, we describe a technique for subspace clustering that uses sparse optimization. The method makes it possible to compute one subspace from the entire mixed dataset by minimizing the number of nonzero vectors obtained when the data are projected onto the orthogonal complement of a single subspace. In contrast to the methods mentioned above, the method presented here can estimate subspaces that are not independent, for example an arrangement of hyperplanes. The principle of the method is as follows. Considering a candidate subspace and the distances of the data points to that subspace, we propose to maximize the number of zero distances. Because this is a hard sparse optimization problem, we later consider a convex relaxation of it. While the initial sparse formulation requires only very mild conditions, the convex surrogate might necessitate a higher level of sparsity to yield the desired solution. This means that there must exist a subspace in the arrangement which contains the majority of the data points. To overcome this limitation, a sequence of reweighted convex optimization problems can be solved [START_REF] Candès | Enhancing sparsity by reweighted 1 minimization[END_REF].

II. DESCRIPTION OF THE PROPOSED METHOD

We will assume, for the sake of exposition clarity, that all subspaces have equal and known dimension d. Also, the subspaces S 1 , . . . , S s are assumed to be discernible in the sense that for any pair (k, ) of indices, S k ⊆ S . But the method to be presented can be extended to arrangements of subspaces with different dimensions.

For any i = 1, . . . , s, let B i ∈ R n×d , with rank(B i ) = d, be a basis of the subspace S i . Denote with S ⊥ i the orthogonal complement of S i in R n . If we assume noiseless data, then for any y j ∈ S i , the orthogonal projection of y j onto S ⊥ i is equal to zero, that is

(I n -P i ) y j = 0 (1)
where P i is the projection matrix defined by P i = B i (B i B i ) -1 B i and I n stands for the identity matrix of order n. Note in passing that P i is uniquely determined by the subspace S i ; it does not depend on any specific basis B i of S i . For, if the matrix B i forms a basis for S i , any other basis has necessarily the form Bi = B i T with T ∈ R n×n being a nonsingular matrix. As a consequence, it can be checked that the corresponding projection matrix Pi is equal to P i . We can therefore describe the subspaces S i with the matrices P i , i = 1, . . . , s. Each P i and I n -P i are orthogonal projections matrices that project the ambiant space onto S i and S ⊥ i respectively. Let us formalize. Definition 1. A square matrix P ∈ R n×n is called an orthogonal projection matrix (or an orthogonal projector) if P is idempotent and symmetric, that is, P 2 = P and P = P .

At this point, it is worth recalling some useful properties of orthogonal projectors. Lemma 1.

1) If P i is an orthogonal projection matrix onto S i , then im(P i ) = ker(I n -P i ) = S i and ker(P i ) = im(I n -P i ) = S ⊥ i , 2) The eigenvalues of an orthogonal projector are 0 or 1.

Hence, any orthogonal projector P i satisfies I n P i 0, where A B means that A -B is a positive semidefinite matrix.

3) The trace of an orthogonal projector equals its rank.

Let now Y = y 1 • • • y N ∈ R n×N
be a matrix formed with all the available data samples and define E(P ) to be the projection of Y onto S ⊥ i , i.e., E(P ) = (I n -P ) Y.

(

) 2
The method to be presented is based on the fundamental observation that the matrix E(P ) contains a significant number of zero column vectors. We shall then say that E(P ) is sparse (to some degree) as a sequence of vectors. More precisely, all the columns of the form (I n -P i )y j with y j ∈ S i , are equal to zero. Let E 2 (P ) be the vector obtained by applying the euclidean norm • 2 to the columns of E(P ), i.e.,

E 2 (P ) = (I n -P ) y 1 2 . . . (I n -P ) y N 2 ∈ R N .
We can then search for one matrix P i by solving the sparse optimization problem

Here, • 0 refers to the 0 pseudo-norm which counts the number of nonzero entries in a vector.

Assumption 1. The data {y j } N j=1 are in general position in each subspace S i , that is, any group {y j1 , . . . , y j d } ⊂ S i of d data samples is formed of linearly independent vectors. Theorem 1. Let P denote any solution to (3). For any P ∈ R n×n , define I (P ) = {y j : (I -P ) y j = 0}. Under Assumption 1, if |I (P i )| ≥ sd for all i = 1, . . . , s, then

P ∈ {P 1 , . . . , P s } .
Proof. Because the matrices P i are feasible points of the optimization problem (3), it holds that |I (P )| ≥ |I (P i )| for i = 1, . . . , s. We claim that there is an index i * ∈ {1, . . . , s} such that |I (P ) ∩ I (P i * )| ≥ d. Assume, for contradiction, that this is not the case. Then by letting n i denote the cardinality of the set I (P ) ∩ I (P i ), we have n i < d, for i = 1, . . . , s. Since all the data points lie in It follows that P possesses 1 as an eigenvalue with geometric multiplicity larger or equal to d = dim(S i * ). Since P is symmetric, it is diagonalizable. Hence, there exists an orthonormal matrix U ∈ R n×n such that P can be factorized in the form

I (P 1 ) ∪ • • • ∪ I (P s ), it follows that |I (P )| ≤ n 1 + • • • + n s < sd ≤ |I (P i )| , i = 1, . . . ,
P = U ΣU , with Σ = diag(λ 1 , . . . , λ d , λ d+1 , . . . , λ n ), λ 1 = λ 2 = • • • = λ d = 1.
Recalling that rank(P ) = d, we have necessarily, λ d+1 = • • • = λ n = 0 so that P obeys (P ) 2 = P . It can also be noticed that rank(I n -P ) = rank(U (I n -Σ)U ) = n -d, which, by the rank theorem, implies that dim (ker(I n -P )) = d and S i * = ker(I n -P ).

From this, we deduce that P is in fact the orthogonal projection matrix onto S i * . Hence P = P i * concluding the proof.

Problem ( 3) is in general very hard to solve exactly because it involves a combinatorial search. The difficulty arises essentially from the presence of the 0 -norm objective and the rank constraint. To get a more tractable optimization problem, a classical idea intensively used in the field of compressed sensing [START_REF] Candès | Reflections on compressed sensing[END_REF], [START_REF] Donoho | Compressed sensing[END_REF], is to approximate the 0 pseudo-norm with the 1 -norm1 , the main advantage being that the last one is convex and therefore efficiently solvable. This leads to the following formulation min P W E (P ) 1 s.t. P = P , I n P 0, rank(P ) = d [START_REF] Costeira | A multibody factorization method for independently moving objects[END_REF] where W = diag(w 1 , • • • , w N ) is a user-defined weighting diagonal matrix with positive entries. Prior knowledge, when available, can be incorporated in W ; by default W is set to the identity matrix I N . The next step is to relax the rank constraint. Knowing from Lemma 1 that the rank of a projection matrix equals its trace, we can try to enforce the rank constraint by replacing it with a constraint on the trace. As a result, we obtain a completely convex relaxation of (3) in the following form

min P N i=1 w i (I n -P ) y i 2 s.t. P = P , I n P 0, tr(P ) = d. (6)
We first look at well-posedness by showing that if a single subspace was to be estimated, then ( 6) is a convenient formulation for recovering the orthogonal projection matrix.

Lemma 2. Assume s = 1, i.e., the data are drawn from a single subspace S 1 . If rank(Y ) = dim(S 1 ) = d, then problem (6) has a unique solution P which corresponds exactly to the orthogonal projection matrix onto S 1 .

Proof. Let P 1 be the orthogonal projection matrix onto S 1 . It is clear that P 1 is a feasible point for problem [START_REF] Gruber | Multibody factorization with uncertainty and missing data using the EM algorithm[END_REF] and minimizes the associated cost function (whose optimal value is actually equal to zero). Therefore P 1 is obviously one solution to problem [START_REF] Gruber | Multibody factorization with uncertainty and missing data using the EM algorithm[END_REF]. We just need to prove uniqueness. We do so by showing that any solution P to (6) is necessarily equal to P 1 . In effect, if P solves (6), then (I n -P ) Y = 0. Since rank(Y ) = d, we have (I n -P ) y = 0 for any y ∈ S 1 , the span of Y . In other words, 1 is an eigenvalue of P with geometric multiplicity at least d. Let λ 1 , . . . , λ n be the eigenvalues of P enumerated in descending order. We have

λ 1 = • • • = λ d = 1, so that λ d+1 +• • •+λ n = tr(P )-d = 0.
With P 0, this leads to λ d+1 = • • • = λ n = 0. P is the orthogonal projection matrix onto S 1 and therefore P = P 1 .

Consider now the general case where s ≥ 2.

Next we state a condition under which a subspace can be exactly recovered by solving [START_REF] Gruber | Multibody factorization with uncertainty and missing data using the EM algorithm[END_REF].

Let H = ∆ ∈ R n×n : ∆ = ∆ , tr(∆) = 0, ∆ 2 ≤ 1 . Theorem 2. Assume that w 1 + • • • + w N = 1. Let r(Y ) = max ∆ 2 s.t. ∆ ∈ H, ∆Y 2,col = 1 (7)
where A 2,col is, for any matrix A, the sum of the 2-norms of the column vectors of A. If there is a P i satisfying j∈I (Pi)

w j > 1 - 1 2ηr(Y W ) (8)
with η = max j y j 2 , then P i is the unique solution to problem (6).

Proof. P i uniquely solves problem (6) if, for any nonzero matrix ∆ in R n×n such that P i + ∆ is feasible, it holds that N j=1 w j (I n -P i ) y j 2 < N j=1 w j (I n -P i -∆) y j 2 .

Feasibility of P i + ∆ translates into tr(∆) = 0, ∆ = ∆ and I n -P i ∆ -P i . These, together with the fact that P i is a projection matrix implies that ∆ lies H. An equivalent form of the inequality above is

j / ∈I (Pi) w j (I n -P i ) y j 2 -(I n -P i -∆) y j 2 < j∈I (Pi)
w j ∆y j 2 .

From the triangle inequality property of vector 2-norm, the following identity a 2 -a + b 2 ≤ b 2 holds for any two vectors a and b of compatible dimensions. By applying this, it can be observed that for P i to be a solution to problem [START_REF] Gruber | Multibody factorization with uncertainty and missing data using the EM algorithm[END_REF], it is sufficient that j / ∈I (Pi)

w j ∆y j 2 < j∈I (Pi)
w j ∆y j 2

for any ∆ ∈ R n×n obeying the conditions above. Adding the term j / ∈I (Pi) w j ∆y j 2 on each side of the previous inequality yields, after simple algebraic manipulations,

max ∆Y W 2,col =1 ∆∈H j / ∈I (Pi) w j ∆y j 2 < 1/2.
Note that the maximum here could have been taken over the set

H ∩ ∆ : ∆Y W 2,col = 1 ∩ ∆ : I n P i + ∆ 0 .
Thanks to the identity max ∆∈A ∆ 2 ≤ max ∆∈B ∆ 2 whenever A ⊂ B, we can restrict the maximum to the first two members of the above chain of intersections and get a more conservative condition. Finally, from the inequality ∆y j 2 ≤ η ∆ 2 , it can be seen that the maximum above is bounded from above by η 1 -j∈I (Pi) w j r(Y W ). Therefore P i is the unique solution to (6) if

η 1 - j∈I (Pi) w j r(Y W ) < 1/2 as claimed.

A. Enhancing sparsity

As shown by Theorem 1, the formulation (3) can obtain a solution to the subspace clustering problem under mild conditions. However, for the convex relaxation [START_REF] Gruber | Multibody factorization with uncertainty and missing data using the EM algorithm[END_REF] to yield the desired solution, a relatively high level of sparsity might be needed. That is, a large majority of the data must pertain to the same subspace. This may not be the case in general. A possible solution is to solve a sequence of problems of the type [START_REF] Gruber | Multibody factorization with uncertainty and missing data using the EM algorithm[END_REF] with different weights computed iteratively [START_REF] Candès | Enhancing sparsity by reweighted 1 minimization[END_REF], [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF]. The iterative scheme can be defined for a fixed number r max of iterations as follows. At iteration r = 0, . . . , r max , compute P (r) = arg min with weights defined, for all j by w (0) j = 1/N , and w

(r) j = v (r) j / N j=1 v (r) j for r ≥ 1, where v (r) j = 1 (I n -P (r-1) )y j ν 2 + ζ , ζ >
0 is a small number preventing division by zero, r is the iteration number and ν ≥ 1. Since we are dealing here with a sequence of convex optimization problems, they can be numerically implemented using any convex solver, see e.g., [START_REF] Grant | CVX: Matlab software for disciplined convex programming[END_REF]. Many simulations results show that the iterative scheme increases significantly the recoverability of the P i 's in situations where the degree of sparsity is low. After a first parameter matrix P 1 is identified, we still need to identify the other subspaces. For this purpose, we can remove the data samples pertaining to the already identified subspace and repeat the iterative algorithm over the remaining set of data. The different steps are described in Algorithm II-A.

III. SPECIAL CASE OF INDEPENDENT SUBSPACES

In this section we analyze the particular case where the subspaces S 1 , . . . , S s are independent, that is when they satisfy

S i ∩ (S 1 + • • • + S i-1 + S i+1 + • • • + S s ) = {0}
for any i. If this is the case, then for all i, Si ⊕ j =i S j ⊂ ker(P i ) while im(P i ) = S i . Here ⊕ refers to direct sum of Algorithm 1 Subspace clustering algorithm Inputs: {y j } N j=1 1) Initialization: P ← ∅, J ← {1, . . . , N } 2) Repeat

• Compute, up to convergence, the sequence P (r) defined in [START_REF] Liu | Robust recovery of subspace structures by low-rank representation[END_REF] based on the data whose indices are contained in J. • Record the identified matrix P : P ← P ∪ P .

• Remove from J the indices of data y j satisfying (I n -P )y j 2 ≤ δ, where δ is a threshold. Until |J| = 0 3) Return P and set s = |P|. subspaces. For s > 2, the number of samples pertaining to subspace Si is generally larger than those pertaining to S i . As a result, solving problem ( 6) is likely to yield the orthogonal projector onto Si . More precisely, the following corollary of Theorem 2 holds. Proof. Denote with P the projection matrix onto Si o ⊕ j =i o S i . Since the subspaces are independent, tr( P ) = rank( P ) = dim( Si o ) = d(s -1). Hence P is feasible for the underlying optimization problem in [START_REF] Vidal | Subspace clustering[END_REF]. The condition min i=1,...,s |I (P i )| < 1 2ηr(Y ) is equivalent to I ( P ) > N -1 2ηr(Y ) , which, by Theorem 2, implies that P is the unique solution to [START_REF] Vidal | Subspace clustering[END_REF]. It follows that P = P . Data which are not in I (P ) are necessarily in S i o = span(I (P i o )).

The corollary suggests that it is presumably easier to segment the subspaces when they are independent.

IV. ROBUST FORMULATION

The method derived above can still work in the presence of a moderate amount of noise. However, when the noise level is high, we might need to resort to a more robust implementation. The convex problem (6) can, similarly as in [START_REF] Bako | Identification of MIMO switched state-space models[END_REF], be reformulated in the form

where ε is a user-defined parameter to be tuned in function of the noise level. The rationale behind the formulation [START_REF] Donoho | Compressed sensing[END_REF] is that, for an appropriately chosen ε, the distances (I n -P i )y j 2 should be less than ε for all the y j close to the subspace spanned by P i .

V. EXPERIMENTS

This section provides some numerical results to illustrate the previous theory. Two experiments are considered: one with noise-free data and the other with noisy data. Both experiments are conducted in the same way as follows. We generate 100 datasets in R 6 . Each dataset contains 400 data points normally sampled from s = 4 different subspaces of the same dimension d. To challenge the proposed algorithm, each dataset is designed such that the number of samples originating from all subspaces is the same (this is indeed the most difficult scenario according to the discussions above). The algorithm is then run on each of the 100 datasets. The experiment is repeated for different values of the subspaces' dimension d = 1, . . . , 5. Results are reported in Tables I andII. For more comments on the effects of the user-defined parameters ζ and δ, see [START_REF] Bako | Identification of switched linear systems via sparse optimization[END_REF].

A number of observations can be made from these results: (a) The reweighted scheme (9) endows the proposed framework with a capability to accommodate subspace clustering problems where none of the subspaces contains the absolute majority of the data; (b) The proposed method can segment arrangements of subspaces with non trivial intersections. However, such arrangements seem more difficult to segment (more iterations are needed) than those formed of independent subspaces. Comparison with the SSC algorithm [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF] reveals the benefit of the proposed approach. dimension d 1 2 3 4 5 # iterations 2 5 5 10 10

VI. CONCLUSION

We have described a technique for estimating multiple subspaces from unlabeled data which lie in the union of the subspaces. The main idea hinges on solving a sequence of sparsity-inducing optimization problems. In contrast to comparable methods using sparse representation, the subspaces to be segmented need not be independent here. Theoretical correctness of the method is analyzed under some sufficient condition related to the proportion of data drawn from each subspace.

  . P = P , I n P 0, rank(P ) = d.

  -P )y j 2 , s.t. P = P , I n P 0, tr(P ) = d[START_REF] Liu | Robust recovery of subspace structures by low-rank representation[END_REF]

Corollary 1 .

 1 Assume that the subspaces S 1 , . . . , S s are independent. Define P as P = arg min P N j=1 (I n -P )y j 2 s.t. P = P , I n P 0, tr(P ) ≤ d(s -1).

( 10 )

 10 If min i=1,...,s |I (P i )| < 1 2ηr(Y ) , then {y j : j ∈ I c (P )} ⊂ S i o , where i o = arg min i |I (P i )| and I c (P ) = {1, . . . , N } \ I (P ).

w

  j ξ j s.t. (I n -P )y j 2 ≤ ε + ξ j , j = 1, . . . , N P = P , I n P 0, tr(P ) = d ξ j ≥ 0, j = 1, . . . , N

  s. However this contradicts the definition of P , hence proving the existence of an index i

* such that I (P ) ∩ I (P i * ) contains at least d different vectors. Note that span [I (P ) ∩ I (P i * )] ⊂ ker(I n -P i * ) = S i * (4) By invoking Assumption 1 with the fact that |I (P ) ∩ I (P i * )| ≥ d, the dimension of the span of I (P ) ∩ I (P i * ) is equal to d. Since dim(S i * ) = d, it can be concluded from (4) that span [I (P ) ∩ I (P i * )] = S i * ⊂ ker(I n -P ).

TABLE I NUMBER

 I OF ITERATIONS OF THE REWEIGHTED SCHEME (9) (ON NOISE-FREE DATA) NECESSARY FOR ACHIEVING 100% SUCCESSES IN THE CLUSTERING. PARAMETERS ARE ζ = 1e -3, δ = 0.1, ν = 5.

	dimension d	1	2	3	4	5
	this paper	0.022 0.027 0.030 0.029 0.035
	SSC [7]	0.074 0.158 0.350 0.519 1.028

TABLE II AVERAGE

 II ESTIMATION ERROR FOR 10 ITERATIONS OF THE REWEIGHTED SCHEME (9) ON NOISY DATA. THE NOISE IS GENERATED FROM A GAUSSIAN DISTRIBUTION OF MEAN ZERO AND VARIANCE 0.05.

The 1 -norm of a vector z = [z 1 . . . z N ] ∈ R N is given by z 1 = |z 1 | + . . . + |z N |.

ACKNOWLEDGMENT

The author thanks Prof. René Vidal for stimulating discussions.