MIXED-FEATURES VECTORS & SUBSPACE SPLITTING
Year: 2021

Abstract:
Motivated by metagenomics, recommender systems, dictionary learning, and related problems, this paper introduces subspace splitting (SS): the task of clustering the entries of what we call a mixed-features vector, that is, a vector whose subsets of coordinates agree with a collection of subspaces. We derive precise identifiability conditions under which SS is well-posed, thus providing the first fundamental theory for this problem. We also propose the first three practical SS algorithms, each with advantages and disadvantages: a random sampling method , a projection-based greedy heuristic , and an alternating Lloyd-type algorithm ; all allow noise, outliers, and missing data. Our extensive experiments outline the performance of our algorithms, and in lack of other SS algorithms, for reference we compare against methods for tightly related problems, like robust matched subspace detection and maximum feasible subsystem, which are special simpler cases of SS.

text: Motivated by metagenomics, recommender systems, dictionary learning, and related problems, this paper introduces subspace splitting (SS): the task of clustering the entries of what we call a mixed-features vector, that is, a vector whose subsets of coordinates agree with a collection of subspaces. We derive precise identifiability conditions under which SS is well-posed, thus providing the first fundamental theory for this problem. We also propose the first three practical SS algorithms, each with advantages and disadvantages: a random sampling method , a projection-based greedy heuristic , and an alternating Lloyd-type algorithm ; all allow noise, outliers, and missing data. Our extensive experiments outline the performance of our algorithms, and in lack of other SS algorithms, for reference we compare against methods for tightly related problems, like robust matched subspace detection and maximum feasible subsystem, which are special simpler cases of SS.

heading: 1 INTRODUCTION
text: As the reach of data science expands, and as we continuously improve our sensing, storage and computing capabilities, data in virtually all fields of science keeps becoming increasingly highdimensional. For example, the CERN Large Hadron Collider currently “generates so much data that scientists must discard the overwhelming majority of it, hoping that they’ve not thrown away anything useful” [1], and the upcoming Square Kilometer Array is expected to produce 100 times that [2]. Fortunately, high-dimensional data often has an underlying low-dimensional structure. Inferring such structure not only cuts memory and computational burdens, but also reduces noise and improves learning and prediction. However, higher dimensionality not only increases computational requirements; it also augments data’s structure complexity. In light of this, several research lines have explored new low-dimensional models that best summarize data, going from principal component analysis (PCA) [3–11] and single subspaces [12–21] to unions of subspaces [22–40], and algebraic varieties [41].
This paper introduces mixed-features vectors (MFV’s): a new model that describes the underlying structure of data arising from several modern applications that is not captured by existing lowdimensional models. The main idea is that each entry of a MFV comes from one out of several classes, and that the entries of the same class lie in an underlying subspace. In particular, MFV’s are motivated by megatenomics [42–46] and recommender systems [47–59]: in metagenomics each gene segment comes from one of the several taxa present in a microbiome; in recommender systems each rating may come from one of several users sharing the same account. However, MFV’s also have applications in robust estimation (e.g., robust PCA [3–11] and robust dictionary learning [60–67]), matrix completion [48–59], subspace clustering [22–39], and more.
This paper also introduces subspace splitting (SS): the task of clustering the entries of a MFV according to its underlying subspaces. SS is tightly related to other machine learning problems. In particular, SS can be thought as a generalization of robust matched subspace detection (RMSD) [12–17], and maximum feasible subsystem (MAXFS) [68–74]. However, the added complexity of SS renders existing approaches for these problems inapplicable, which calls the attention for specialized SS theory and methods. In these regards, (i) we derive precise identifiability conditions under which SS is well-posed, and (ii) we propose the first three SS algorithms.

heading: 2 PROBLEM STATEMENT AND FUNDAMENTAL THEORY
text: Let U1, . . . ,UK be subspaces of Rd, and let Ω0,Ω1, . . . ,ΩK denote a partition of [d] := {1, . . . ,d}. For any subspace, matrix or vector that is compatible with a set of indices Ω ⊂ [d], we will use the subscript Ω to denote its restriction to the coordinates/rows in Ω. For example, U1ΩK ⊂ R
|ΩK| denotes the restriction of U1 to the coordinates in ΩK. Define x ∈ Rd as the mixed-features vector (MFV) such that xΩk ∈ U k Ωk
for each k = 1, . . . ,K, and the entries of xΩ0 are outliers. Let ∈ Rd denote a noise vector with variance σ2. Given U1, . . . ,UK, and an incomplete observation yΩ = xΩ + Ω, the goal of subspace splitting (SS) is to determine the subsets Ω1 ∩ Ω, . . . ,ΩK ∩ Ω indicating the observed coordinates of y that match with each subspace.
Example 1. Consider the following setup, with 1-dimensional subspaces U1,U2 spanned by U1,U2:
U1 =  1 1 1 1 1 1  , U2 =  1 2 3 4 5 6  , x =  1/2 1/2 6 8 9 10  , =  0.1 −0.1 −0.1 0.1 −0.1 0.1  , yΩ =  0.51 0.49 5.9 8.1 8.9 ·  .
It is easy to see that Ω1 = {1, 2}, Ω2 = {3, 4}, Ω0 = {5}, because xΩ1 = 12U 1 Ω1 and xΩ2 = 2U 2 Ω2 .
The keen reader will immediately wonder: is there another partition {Ω′1, . . . ,Ω ′ K} different from {Ω1, . . . ,ΩK} such that xΩ′k ∈ U k Ω′k
for every k? In other words, is this problem well-posed, and if so, under what conditions? Our main theoretical result answers this question, showing that under the next assumptions, Ωk can be recovered if and only if it has more elements than the dimension of Uk.
A1 Each Uk is drawn independently with respect to the uniform measure over the Grassmannian.
A2 Each xΩk is drawn independently according to an absolutely continuous distribution with respect to the Lebesgue measure on UkΩk .
In words, A1 essentially requires that U1, . . . ,UK are in general position with no particular relation with one another. Similarly, A2 requires that each piece of x is in general position over its corresponding piece of subspace. This type of genericity assumptions are becoming increasingly common in compressed sensing, matrix completion, subspace clustering, tensor theory, and related problems [10, 21, 30, 31, 33–38, 41, 56–59]. All our statements hold with probability 1 with respect to the measures in A1 and A2. We point out that A1 and A2 do not imply coherence or affinity (other typical assumptions in related theory that quantify alignment with the canonical axes or between subspaces [3, 6, 11, 23, 26, 27, 48–50, 54, 55]) nor vice-versa. For example, bounded coherence and affinity assumptions indeed allow subspaces perfectly aligned on some coordinates. However, they rule-out cases that our assumptions allow, for example the non-zero measure set of highly coherent or affine subspaces that are somewhat aligned with the canonical axes or with one another. To sum up, these assumptions are different, not stronger nor weaker than the usual coherence and affinity assumptions. With this, we are ready to state our main theorem, showing that subspace splitting is possible if and only if x contains more than dim(Uk) entries of each subspace Uk.
Theorem 1. Suppose A1 and A2 hold. Given x and U1, . . . ,UK, one can identify Ω1, . . . ,ΩK if and only if |Ωk| > dim(Uk) for every k.
Example 1 shows a case where the conditions of Theorem 1 are met (|Ω1| = |Ω2| = 2 > dim(U1) = dim(U2) = 1), and consequently subspace splitting is well-posed (there exist no partition other than the true {Ω1,Ω2} that splits x into U1 and U2). Conversely, the following Example shows a case where the conditions of Theorem 1 are not satisfied, and at least some Ωk is unidentifiable.
Example 2. Consider the following setting:
U1 = 111 1  , U2 = 123 4  , U3 = 014 9  , x =  1/2 1/2 6 3  . Now |Ω2| = |Ω3| = 1 = dim(U2) = dim(U3). As Theorem 1 shows, there exist multiple ways to split x into U1, U2, and U3. Here the partitions could be {Ω1,Ω2,Ω3} = {{1, 2}, {3}, {4}} or {{1, 2}, {4}, {3}}, and there is no way of telling which is the true partition from U1, U2, U3, and x. In other words, Ω2 and Ω3 are unidentifiable. Remark 1. We point out that the constructions in our examples are not generic (i.e., they were not constructed according to A1 and A2). We chose them for their simplicity to build intuition and make a point. However, our results still apply to these constructions, showing that in addition to all generic cases, our theory also holds for some non-generic ones. An exact characterization of all the non-generic cases that our theory covers requires a careful study of notions of sketching and partial coordinate discrepancy [31], which are out of the scope of this paper.
The proof of Theorem 1 follows by the next two lemmas. Lemma 1 shows that any subset of r or fewer entries of any vector will always match with any r-dimensional subspace in general position. Lemma 2 shows that r + 1 entries of a vector won’t match with a random r-dimensional subspace by chance. Said in other words, r + 1 entries of a vector will match with an r-dimensional subspace if and only if such entries truly come from that r-dimensional subspace. Lemma 2 is effectively the key towards Theorem 1, as it allows us to try all combinations of r + 1 entries that fit in an r-dimensional subspace, knowing that we will never get a false match. All proofs are in Appendix A.
Lemma 1. Suppose A1 holds. Let Ω ⊂ [d]. If |Ω| ≤ dim(Uk), then xΩ ∈ UkΩ for every x ∈ Rd. Lemma 2. Suppose A1 and A2 hold. Let Ω be an arbitrary subset of [d] with exactly dim(Uk) + 1 elements. Then xΩ ∈ UkΩ if and only if Ω ⊂ Ωk.
Corollary 1 extends these results to account for noise and missing data, replacing A1 and A2 with:
A1’ The coherence of Uk is upper bounded by µ, and the geodesic distance over the Grassmannian between Uk and U` is lower bounded by ϕ, for every k, ` = 1, . . . ,K.
A2’ The coherence of xΩk is upper bounded by ν, and its norm is lower bounded by ψ, for every k = 1, . . . ,K.
Corollary 1. Suppose A1’ and A2’ hold with µ, ν < Cσ, and ϕ,ψ > cσ for some constants C and c. Let yΩ and U
1, . . . ,UK be given. Suppose |Ωk ∩ Ω| > dim(Uk) for every k. Then with probability decreasing in C and increasing in c, one can identify Ω1 ∩ Ω, . . . ,ΩK ∩ Ω.
To guarantee identifiability, Corollary 1 requires that subspaces and samples are sufficiently incoherent and separated to overcome the noise level σ. Notice that since now there is missing data (we observe yΩ, instead of x as in Theorem 1), Corollary 1 requires that there are enough observed entries per subspace, i.e., that |Ωk ∩ Ω| are sufficiently large (rather than Ωk). Similarly, only the observed entires can be classified, so only the intersections Ωk ∩ Ω are identifiable (rather than Ωk).

heading: 3 RELATED WORK
text: Arguably, the closest link that subspace splitting has with other popular machine learning problems, is through robust matched subspace detection (RMSD) [12–17] and maximal feasible subsystem (MAXFS) [68–74]. In RMSD one observes a vector x ∈ Rd of the form x = Uθ+ z, where z ∈ Rd is a sparse vector of outlier entries, and θ ∈ Rr is the coefficient of the entries of x that match with the subspace spanned by U ∈ Rd×r. Given x and U, the goal is to equivalently find θ, z, or the support of z, which in turn determine the coordinates Ω where x matches U. Subspace splitting can be thought as the generalization of RMSD to the case where there are multiple subspaces, and we want to find the entries that match each one. Similarly, in MAXFS one has an inconsistent system of equations of the form x = Uθ, and wants to determine the solution θ that produces the largest
subset Ω of consistent equations. Subspace splitting can be thought as the generalization of MAXFS to the case where there are multiple consistent subsystems (one per subspace), and we want to find each of them. The three prevalent venues to solve these problems can be classified into three groups: (i) `1 minimization, (ii) mixed integer programs, and (iii) random sampling.
The first group essentially encompasses variants of minimizing ‖x−Uθ‖1 over θ ∈ Rr, which is tightly related to the Lasso [91]. The intuition is that the `1-norm will favor solutions that produce a sparse vector z = x−Uθ whose zero entries reveal Ω. While formulations like this work great for just one subspace (or subsystem), if there are two or more, then for each subspace, the entries of all other subspaces are outliers, in which case zk := x−Ukθk will no longer be sparse, and the solution to the `1-norm minimization will no longer reveal Ωk. The second group aims to directly recover Ω with variants of the mixed integer formulation that minimizes ‖z‖ subject to |x−Uθ| ≤ Mz over z ∈ {0, 1}d and θ ∈ Rr, where M is a tuning parameter. Here the main idea is to maximize the number of zero entries in z (which reveal Ω), where we are forcing x to match span{U}, because the constraint guarantees that |xi −Uθ| ≤ zi = 0. Unfortunately, because of the combinatorial nature of this approach, this method becomes exponentially slow as the fraction of outliers/inconsistencies increases, rendering it prohibitive in practice for even small values of K.The third group essentially uses the same principle as random sampling consensus (RANSAC) [92–98]. That is, one iteratively selects a candidate subset of entries at random, until it finds one where x agrees with span{U}. This is in fact the approach that we used in the proof of Theorem 1, and the main idea behind our first SS algorithm, which we present in Section 5.1, and study in Section 6.
Remark 2. We point out that SS is also related to subspace clustering [22–40]. The main difference is that in subspace clustering all entries of each column come from the same subspace.

heading: 4 MOTIVATING APPLICATIONS
text: This section details the main motivations behind SS, namely, metagenomics, recommender systems, and robust estimation, relevant for dictionary learning, matrix completion, and related problems.
Metagenomics Microbial communities affect the balance of entire ecosystems, and play a crucial role in the planet’s health [42]. Consequently, understanding microbiome compositions, interactions, and evolutions holds the key to the knowledge of the deep interconnectivity factors that play a role on Earth’s biodiversity, our agricultural sustainability, and adaptation to global threats like climate change. One cornerstone tool towards microbiome understanding lies in genomic analysis. In practice, to obtain an organism’s genome (e.g., a person’s genome), biologists feed a DNA sample (e.g., blood or hair) to a sequencer machine that produces a series of reads, which are short genomic sequences that can later be assembled and aligned to recover the entire genome. The challenge arises when the sequencer is provided a sample with DNA from multiple organisms, as is the case in any microbiome (e.g., the human gut microbiome), where any sample will contain a mixture of DNA from multiple taxa that cannot be trivially classified [43]. In this case, each read produced by the sequencer may correspond to a different taxa, resulting in a DNA sample with a mixture of genes. Existing approaches depend on the correct classification of taxonomic units, which in turn relies on the existence of reference tables, and often require human intervention [44]. However, reference sequences remain unknown for a vast majority of microbial biodiversity, which in turn precludes holistic metagenomic analyses, as scientists often have to discard unidentified data that cannot be currently categorized with existing reference tables [45]. In contrast, our work pioneers the theoretical foundations of a novel model tailored to automatically classify taxonomic units without any human intervention. Moreover, our model will naturally allow for missing data, and hence it will be robust to taxa with length-varying sequences, and fast-occurring mutations, which are in fact quite common in several microorganisms such as certain types of bacteria [46].
Recommender Systems and Image Inpainting In recommender systems like Amazon or Netflix, one aims to infer users’ preferences in order to make good recommendations [47]. Arguably the most renowned model for this task is low-rank matrix completion [48–58], which assumes that each user’s preferences vector can be explained as a linear combination of a few others. Said in other words, preferences lie in a linear subspace. However, often multiple types of users share the same account. In this case the vector of preferences will contain a mixture of entries from multiple subspaces. It is easy to see (details in the proof of Theorem 1 in Appendix A) that the coefficient corresponding to the entries in the kth subspace are given by θk = (UkΩ′k) −1xΩ′k , where Ω ′ k is any subset of Ωk with
more than dim(Uk) elements. This insight can be used to estimate the values corresponding to the kth subspace that would appear in the locations where there are observations from other subspaces, or where data is missing. In For example, in recommender systems, this would allow to infer the preferences of all users sharing an account. In image inpainting this would allow to estimate the missing, corrupted, or occluded pixels in an image [60–67].
Robust Learning In Section 3 we showed that SS is a generalization of RMSD and MAXFS, both of which are crucial subtasks in many importantproblems, including robust PCA [3–11], dictionary learning [60–67], matrix completion [48–58], and more. For example, in subspace tracking one needs to identify outliers in each new vector (meaning performing RMSD or MAXFS) before updating the underlying subspace [18–20]. Or in subspace clustering one needs to identify outliers in each vector before finding sparse representations [26, 28, 29, 32, 39]. These in turn are core techniques in target localization [76], medical imaging [79], communications [80], anomaly detection [81, 82], hyperspectral imaging [16], and networks estimation [84? –89], among many others [78]. Being a generalization of RMSD and MAXFS, subspace splitting is also applicable in these broader problems, bringing possibilities for improved performance.Moreover, as these problems evolve into more sophisticated models that allow mixed-features vectors (such as mixture matrix completion [59]), subspace splitting will gain importance and become a crucial subroutine of these emerging problems.

heading: 5 SUBSPACE SPLITTING ALGORITHMS
text: We now present our three subspace splitting algorithms: a random sampling method, an iterative projection-based greedy heuristic, and an alternating approach inspired by Lloyd’s algorithm [75].

heading: 5.1 RANDOM SAMPLING SPLITTING
text: Like Theorem 1 suggests, for each k we can iteratively try sets Ω′k with dim(U k) + 1 entries selected randomly until we find one such set where yΩ′k lies close to U k Ω′k (within the noise level σ). At that point we can estimate θk = (UkΩ′) −1yΩ′ , and subsequently Ωk as the entries in yΩ −U k Ωθ
k that are within σ. We call this approach random sampling splitting (RANSAS). The main advantage of RANSAS is that as consequence of Theorem 1, it is guaranteed to work after enough iterations. To see this, observe that if entries in x are uniformly distributed across subspaces and outliers, then the probability of randomly selecting r + 1 entries from the kth subspace is 1/(K+1)r+1. Consequently, the probability that r + 1 random entries come from the same subspace is K/(K+1)r+1. A simple Chernoff bound shows that the expected number of iterations before this happens is O(Kr). Since we want this to happen K times, we conclude that the expected number of iterations is upper bounded by O(Kr+1). We summarize this in the next theorem.
Theorem 2. Let A1 and A2 hold. Suppose |Ωk| > dim(Uk) for every k. Suppose P(i ∈ Ωk) = 1/(K+1) independently for every i, k. Then the expected number of iterations before RANSAS identifies Ω1, . . . ,ΩK is lower and upper bounded by O(Kmink dim(U k)+1) and O(Kmaxk dim(U k)+1).
We point out that a scenario where all entries are equally distributed is in fact the most difficult setting. Otherwise, identifying any subspace with probability pk > 1/(K+1) will take O(1/pkr+1) < O(1/(K+1)r+1) iterations to find. The main downside of this approach is that its computational complexity scales exponentially. So if K and r are small, or there is one dominant subspace, RANSAS may be a good idea. However, as we see in our experiments, this approach quickly becomes computationally prohibitive as K and r grow.

heading: 5.2 GREEDY SPLITTING
text: The main idea of our greedy splitting (GREEDYS) heuristic is to iteratively remove the most disruptive entry from each subspace. More precisely, for each k we will start with Ω0k = [d], and for each t > 0 we will define Ωtk := Ω t−1 k \it−1, where it−1 is the entry whose removal from Ω t−1 k minimizes the difference between yΩtk and its projection onto U k Ωtk
. The intuition is that removing it−1 gets yΩtk closer to UkΩtk . We repeat this procedure for each k until yΩtk is close to U k Ωtk
(within the noise level σ). We know from Lemma 1 that this procedure will terminate at some point, because any subset of r
or fewer entries of any vector will always match with any r-dimensional subspace in general position. If for some k this procedure terminates with a set Ωtk with more than rk = dim(U
k) entries, Theorem 1 suggests that yΩtk truly lies in U k Ωtk , that Ωtk ⊂ Ω, that we can estimate θ k = (UkΩtk )−1yΩtk , and identify Ωk as the entries in yΩ −U k Ωθ
k that are within σ. At this point we can prune all the entries corresponding to this subspace, and repeat the greedy procedure from the start without these entries. If at some point none of the remaining sets Ωtk has more than rk entries, we know by Lemma 1 that we cannot determine whether yΩtk truly lies in U k Ωtk
with these entries alone, then we stop and consider this a failure. Notice that for each k, GREEDYS removes the most disruptive entry among all d entries, then the next most disruptive among the remaining d− 1, until either all the remaining entries lie in the same subspace, or there remain no more than rk entries. This requires no more than d + (d−1) + (d−2) + · · ·+ (rk + 1) = O(d2) iterations. After identifying the entries of a subspace, if pruning is necessary, we potentially need to repeat the iterative removal process for the remaining K − 1 subspaces, and so on. In total, we potentially need to repeat the iterative removal process K + (K− 1) + (K− 2) + · · ·+ 1 = O(K2) times. Putting these two observations together, we obtain the following theorem, showing that GREEDYS has a polynomial computational complexity. Theorem 3. GREEDYS will terminate after no more than O(Kd)2 iterations.
The main caveat of this greedy approach, as the next Example shows, is that even in a noiseless setting, the most disruptive entry for Uk (that is, the entry whose removal minimizes the distance between yΩtk and its projection onto U k Ωtk ) may in fact correspond to Uk. Example 3. Consider the following construction where the first two entries of x correspond to U1, the last two correspond to U2, there are no outliers, i.e., Ω0 = ∅, there is no noise, i.e., = 0, and there is no missing data, i.e., Ω = [d], so that we observe all entries of y = x:
U1 = 223 4  , U2 = 131 1  , y = 442 2  . Pay special attention to U1 and y. Let Ω−i := [d]\i denote the set of all except the ith entry, and let ŷkΩ−i := U k Ω−i(U kT Ω−iU k Ω−i) −1UkTΩ−i ŷΩ−i denote the projection of yΩ−i onto U k Ω−i . For example, here:
U1Ω−1 = [ 2 3 4 ] , yΩ−1 = [ 4 2 2 ] , ŷ1Ω−1 = [ 1.572 2.2759 3.0345 ] .
Proceeding according to our greedy heuristic, we can compute the projection residuals when removing each entry:
‖yΩ−1 − x̂ 1 Ω−1‖ = 2.7038, ‖yΩ−2 − x̂ 1 Ω−2‖ = 2.7038,
‖yΩ−3 − x̂ 1 Ω−3‖ = 3.4641, ‖yΩ−4 − x̂ 1 Ω−4‖ = 2.7440.
Here we see that removing any of the first two entries makes y appear closer to U1, even though the first two entries in fact correspond to U1. We thus conclude that our greedy heuristic, while it works well in practice, it has no theoretical guarantee of working.

heading: 5.3 K-SPLITS
text: As the name suggests, our last algorithm, K-SPLITS, is inspired by the K-means algorithm. The main idea is to initialize partitioning randomly the entries of yΩ into Ω̂0, Ω̂1, . . . , Ω̂K, and then alternate until convergence between (i) estimating coefficients θ̂1, . . . , θ̂K in light of the assignment of entries, and (ii) reassigning entries according to these coefficients. More precisely, at each step after pruning when possible we:
(i) Estimate coefficients θ̂k = (Uk Ω̂k )−1yΩ̂k . (ii) Compute ŷk := Ukθ̂k, and assign entries as Ω̂k = { i ∈ [d] : |yi − ŷ k i | ≤ τ(σ)}, where
τ(σ) is a tuning thresholding parameter that depends on the noise level.
The assignment step is done on a cyclic manner, that is, in every iteration, each cluster is assigned the entry that agrees the most with it, and it receives no new entries until the next iteration. This changes the isotropic assumption of K-means with the assumption that each cluster has the same number of entries, as mentioned in Section 5.1, this is the most difficult setting for this problem, and whenever the assumption is not fulfilled, the algorithm will still be able to complete at least the most dominant subspace even if entries of that subspace contaminate other subspaces; for the next iteration that won’t be the case any more. Thanks to this, the assumption of every subspace having the same number of entries does not affect the effectiveness of the algorithm.
As we will see in our experiments, this method outperforms the rest. Its main downside, however, is that like most alternating methods, it suffers from local minima, and depends heavily on initialization. To mitigate the former, in our experiments we run several independent random initializations. We observe that in practice this has excellent performance. We conjecture that it is because this strategy will at some point spread the centers θk enough, in light of the subspaces. This sort of spread is crucial in other related clustering algorithms, such as K-means [75]. Our future work will further investigate this, convergence, complexity, and initialization strategies, such as parallels of K-means++ [99], which require careful analysis and are out of the scope of this paper.

heading: 6 EXPERIMENTS
text: In this section we study the performance of our three methods above, as a function of the number of subspaces K, the number of entries per subspace |Ωk| (which in turn is a proxy of the ambient dimension d) and the subspaces dimensions r. We measure error as the fraction of misclassified entries. In lack of other existing subspace splitting methods, for reference we compare them against the `1 minimization described in Section 3. Our results show that our SS methods perform as well as existing methods in the single subspace case, and succeed in the presence of multiple subspaces, where existing methods fail dramatically. In the interest of reproducibility, all our code is available here [100].
In our experiments we first generate matrices U1, . . . ,UK ∈ Rd×r with i.i.d. N(0, 1) entries, to use as bases of U1, . . . ,UK, with d = (K + 1)|Ωk|. Similarly, we generate θ1, . . . ,θK ∈ Rr, also with i.i.d. N(0, 1) entries, to use as coefficient vectors. Next we create the partition of entries {Ω0,Ω1, . . . ,ΩK} where each Ωk has the exact same number of elements. As discussed in Section 5.1, this is in fact the most difficult setting; otherwise the dominant subspace is easier to find, and one can simply prune its corresponding entries to simplify the problem. Then we create x as the vector such that xΩk = U k Ωk
θk for each k, and xΩ0 are i.i.d. N(0, 1) entries representing outliers. Finally, we generate ∈ Rd with i.i.d. N(0, σ2) entries and Ω ⊂ d with a fraction p of observed entries, to construct our observed data yΩ = xΩ + Ω.
In our first experiment we study the behavior of our algorithms as a function of K, with r = 5, |Ωk| = 40, p = 1/2 (so that on average 20 observed entries correspond to each subspace), and σ = 0, i.e., no noise. Figure 1 (a) shows that even for K = 2 (the smallest value of K), existing methodology is unable to split y. This is not surprising, as existing methods are meant for settings where most of the entries in y agree with a single subspace, and there are only a few outliers. In our setting, for each subspace, the entries of all other subspaces are outliers. Consequently, a fraction 1 − 1/K of the entries (at least half with K = 2, and a dominant majority as K increases) are outliers for each subspace, which is well-documented to make existing methods fail (see the discussion in Section 3 for details). In contrast, GREEDYS is better with K = 6 (which amounts to 83.3% outliers) than `1 is with K = 2, while K-SPLITS and RANSAS maintain a 100% success rate for every K we tried. Nonetheless, increasing in K comes at a price (time), as seen in Figure 1 (a’), which shows how the number of iterations of each algorithm grows with K, relative to the simplest case (K = 2). Notice that the speed of K-SPLITS decays nicely with K. For instance, with K = 6, K-SPLITS degrades an order of magnitude slower than RANSAS, while still maintaining the same 100% success rate.
In our second experiment we study the behavior of our algorithms as a function of |Ωk|, with r = 5, p = 1/2, σ = 0, and K = 2 fixed. Figures 1 (b) and (b’) show that as |Ωk| grows, the performance of all methods improves: `1 and GREEDYS achieve higher success rates, while RANSAS and K-SPLITS increase their speed and maintain the same 100% success rate.
In our third experiment we study performance as a function of r, with p = 1/2, σ = 0, and K = 2 fixed. Recall from Lemma 1 that as r increases, so does the number of entries per subspace required for identifiability (r + 1). Hence to allow a wide range of r, we fixed |Ωk| = 40. Figures 1 (c) and (c’), show that (conversely to our previous experiment) the performance of all methods declines as r grows: `1 and GREEDYS achieve lower success rates, while RANSAS and K-SPLITS increase their time (though they maintain the same 100% success rate).
Our next two experiments are consistent with the intuition that cases with fewer data (small |Ωk|) relative to the problem’s complexity (the dimension r) are harder. We verify this in our next experiments, where we study performance as a function of r, fixing p = 1/2, σ = 0, and K = 2 and the gap γ := p(|Ωk| − (r + 1)) = 20. This setting represents the case where we always have a few more (20) entries per subspace than the strictly required (r + 1) for splitting. This is arguably a better experiment to test the effect of the subspaces dimensions isolated from the sampling rate. Figures 1 (d) and (d’) now show that r slightly affects GREEDYS and K-SPLITS(GREEDYS in terms of accuracy, and K-SPLITS in terms of time), but dramatically affects RANSAS, causing an increase of iterations four orders of magnitude higher than the rest. This is consistent with our discussion in Section 5.1, showing that RANSAC will always succeed given enough time, which grows exponentially in r. Interestingly, K-SPLITS achieves the same 100% performance, but orders of magnitude faster, suggesting that K-SPLITS is in general the most promising SS algorithm.
In our last experiment we study the accuracy of our algorithms as a function of the noise level σ, with r = 5, p = 1/2, |Ωk| = 40, and K = 2 fixed. Consistent with Corollary 1, Figure 1 (e) shows that the performance of our algorithms decays nicely with noise, with K-SPLITS allowing up to σ = 10−3. The price, however, can be seen in execution time (Figure 1 (e’)), which grows up to the point where algorithms start to decay, in which case their execution time also starts decreasing. We point out that the fraction of missing data and outliers only affect results and performance in the sense that they reduce the number of available observed entries per subspace, so their effects are also captured in the experiments on |Ωk|. For example, a case where there are |Ωk| = 40 samples per subspace, but only half of the entries are observed (p = 1/2) is equivalent to observing |Ωk| = 20 samples per subspace with no missing data (p = 1). Appendix B shows additional results with K = 4, where our methods show even more dramatic improvements.
BROADER IMPACT
This paper introduces a new dimensionality reduction model, and a series of methods to infer such model, with broad applicability ranging from metagenomics to recommender systems. From a pragmatic standpoint, practitioners can use this research to better predict drug-target interactions, analyze the composition and dynamics soil or human microbiomes, and more. Ultimately, this can impact the well-being of society at many levels, improving drugs, medical diagnoses, treatments, agricultural sustainability, etc. From a learning perspective, we hope this work spurs discussions and insights, and motivates the data science community to explore new directions that stem from this initial work (for example near-optimal initialization strategies on our K-SPLITS algorithm, similar to K-means++), and that leads to better methods and understanding of subspace splitting.

heading: A PROOFS
text: In this section we present the proofs of all statements. We begin with the proof of Lemmas 1 and 2, which are required for the proof of Theorem 1.
Proof. (Lemma 1) Since Uk is in general position, if |Ω| ≤ dim(Uk), then UkΩ = R|Ω|. In other words, UkΩ is the entire |Ω|-dimensional space.
Proof. (Lemma 2) If Ω ⊂ Ωk then xΩ trivially lies in UkΩ by A2. To see the reverse implication, let r := dim(Uk). We know that xΩ ∈ UkΩ if and only if there exists a coefficient vector θ ∈ Rr such that
xΩ = U k Ωθ, (1)
where Uk ∈ Rd×r denotes a basis of Uk. Since |Ω| = r + 1 by assumption, equation 1 defines a system of r + 1 equations (the rows of xΩ) and r unknowns (the entries of θ). Assume without loss of generality that Ω = {1, 2, . . . , r + 1} (otherwise simply permute the rows accordingly), and let Ω′ = {1, 2, . . . , r} be the subset of Ω with the first r elements. Then we can rewrite equation 1 as
r { 1 {  xΩ′ xr+1  =  UkΩ′ Ukr+1 θ. Notice that UkΩ′ ∈ Rr×r is invertible because A1 guarantees that U
k is in general position with probability 1. Hence we can use the top block to obtain θ = (UkΩ′)
−1xΩ′ , and we can plug this in the last row to obtain:
xr+1 = U k r+1(U k Ω′) −1xΩ′ .
Let ki indicate the subspace corresponding to the ith entry of x. Then by A2 there exist coefficients θ1, . . . ,θK such that xi = Ukii θ ki . In light of this we can rewrite the last equation as
U kr+1 r+1 θ kr+1 = Ukr+1(U k Ω′) −1  Uk11 θ k1 Uk22 θ k2
... Ukrr θ kr  . (2) At this point equation 2 has no variables left (all Uk’s and θk’s are fixed), which means that equation 2 is a consistency condition. Notice that if Ω ⊂ Ωk, then ki = k, and equation 2 transforms into
Ukr+1θ k = Ukr+1(U k Ω′) −1UkΩ′θ k,
which further simplifies to Ukr+1θ k = Ukr+1θ k. In other words, if Ω ⊂ Ωk, the consistency condition equation 2 is met. On the other hand, if Ω 6⊂ Ωk, equation 2 does not simplify any further. Recall from A1 that U1, . . . ,UK are drawn independently. This implies that the entries in U1, . . . ,UK keep no relation with one another. Similarly, from A2 we know that xΩk is drawn independently. Equivalently, θ1, . . . ,θK keep no relation with one another. Intuitively, we can think of the entries in Uk’s and θk’s as independent random numbers drawn from a distribution with full measure. We thus conclude that if Ω 6⊂ Ωk, with probability 1 the consistency condition equation 2 would not hold.
To summarize, we know that xΩ ∈ UkΩ if and only if there is a vector θ that satisfies equation 1. We showed that this will be the case if and only if equation 2 holds, which in turn will happen if and only if Ω ⊂ Ωk, as claimed.
Notice that if Ω has dim(Uk) + ` entries, equation 1 becomes a system with dim(Uk) + ` equations, which result in ` consistency conditions similar to equation 2. By the same arguments as in the proof of Lemma 2, xΩ will agree with UkΩ if and only if these consistency conditions hold, which will happen if and only if Ω ⊂ Ωk. We thus obtain the following corollary
Corollary 2. Suppose A1 and A2 hold. Let Ω be an arbitrary subset of [d] with more than dim(Uk) elements. Then xΩ ∈ UkΩ if and only if Ω ⊂ Ωk.
With this, we are ready to give a proof of Theorem 1
Proof. (Theorem 1) Let rk := dim(Uk). To identify Ω1, . . . ,Ωk we will exhaustively search for combinations Ω of rk + 1 entries in x that match with Uk. By Lemma 2 we know that xΩ will match with Uk if and only if Ω ⊂ Ωk. By assumption there is at least one k for which |Ωk| > rk, so we know that at some point we will find such a subset Ω. Once we do, we can recover θk = (UkΩ′)
−1xΩ′ , where Ω′ can be any subset of Ω with exactly rk entries (recall that A1 guarantees that UkΩ′ is invertible). Finally, we can compute xk := Ukθk, and recover Ωk by inspection (the set of zero entries in xk − x), where Corollary 2 guarantees that there will be no additional false matching entries. Repeating this procedure for every k we can recover Ω1, . . . ,ΩK, as claimed.
To obtain the converse, suppose that rk or fewer entries of x correspond to Uk for some k. By Lemma 1, any subset Ω′ of Ωk ∪ Ω0 with |Ω′| ≤ rk satisfies xΩ′ ∈ UkΩ′ . Hence we can split the entries of x corresponding to Ωk and Ω0 in a non-unique manner, and still have them agree with Uk, which makes Ωk unidentifiable.
We now use the same strategy to prove Corollary 1, which extend these ideas to noisy settings.
Proof. (Corollary 1) Let PkΩ′ denote the projection operator onto U k Ω′ . Under the assumptions of the corollary, if Ω′ 6⊂ Ωk, then with high probability (decreasing in C and increasing in c), ‖xΩ′ −PkΩ′xΩ′‖2 > |Ω
′|σ2 (see Theorem 1 in [17]). In other words, the residual of projecting xΩ′ onto UkΩ′ will be too large if the entries of Ω
′ do not come from the kth subspace. Hence, instead of searching for a set Ω′ where xΩ′ perfectly fits in UkΩ′ (as in the proof of Theorem 1), we can search for subsets Ω′ ⊂ Ωk where the residual is within the noise level, i.e., ‖xΩ′ − PkΩ′xΩ′‖2 ≤ |Ω|σ2. The rest of the proof follows by the same arguments as the proof of Theorem 1.

Review summaries:
The paper introduces the problem of subspace spitting, in which an observed mixed-features vector is to be partitioned such that the identified partitions match with given subspaces. The main results of the paper lie in deriving sufficient and necessary conditions for identifiability of these partitions when the subspaces and the entries of the features are randomly positioned in the ambient dimension and the subspaces, respectively. The conditions simply require that there are more entries associated with each subspace than the dimension of the subspace. The paper also presents algorithms to perform the splitting. 

The authors propose a method to perform subspace splitting. That is, the task of clustering the entries of an input vector into sets of coherent subspaces. The contribution of the work is two-fold: (1) the theoretical characterization of the problem, and its well-posedness, and (2) the presentation of three algorithms for tackling the problem of subspace splitting. Quantitative analysis of the performance of the three algorithms is provided by means of synthetic experiments.