Title: CO-MANIFOLD LEARNING WITH MISSING DATA

Abstract: Representation learning is typically applied to only one mode of a data matrix, either its rows or columns. Yet in many applications, there is an underlying geometry to both the rows and the columns. We propose utilizing this coupled structure to perform co-manifold learning: uncovering the underlying geometry of both the rows and the columns of a given matrix, where we focus on a missing data setting. Our unsupervised approach consists of three components. We first solve a family of optimization problems to estimate a complete matrix at multiple scales of smoothness. We then use this collection of smooth matrix estimates to compute pairwise distances on the rows and columns based on a new multi-scale metric that implicitly introduces a coupling between the rows and the columns. Finally, we construct row and column representations from these multi-scale metrics. We demonstrate that our approach outperforms competing methods in both data visualization and clustering.

1. INTRODUCTION:

Dimension reduction plays a key role in exploratory data analysis, data visualization, clustering and classification. Techniques range from the classical PCA and nonlinear manifold learning to deep autoencoders (Tenenbaum et al., 2000; Roweis & Saul, 2000; Belkin & Niyogi, 2003; Coifman & Lafon, 2006; Vincent et al., 2008; Rifai et al., 2011; Kingma & Welling, 2014) . These techniques focus on only one mode of the data, often the observations (columns) which are measurements in a high-dimensional feature space (rows), and exploit correlations among the features to reduce the dimension of the feature vectors and obtain the underlying low-dimensional geometry of the observations. Yet for many data matrices, for example in gene expression studies, recommendation systems, sensor networks, and word-document analysis, correlations exist among both observations and features. In these cases, we seek a method that can exploit the correlations among both the rows and columns of a data matrix to better learn lower-dimensional representations of both. Biclustering methods, which extract distinct biclusters along both rows and columns, give a partial solution to performing simultaneous dimension reduction on the rows and columns of a data matrix. In certain settings, however, assuming a bi-clustering model is too restrictive and results in breaking up smooth geometries into artificial disjoint clusters that do not match the actual structure of the data. This can occur when the true geometry is one of overlapping rather than disjoint clusters, for example in word-document analysis (Ahn et al., 2010) , or when the underlying structure is not one of clusters at all but rather a smooth manifold (Gavish & Coifman, 2012) . Thus, we consider a more general viewpoint: data matrices possess geometric relationships between their rows (features) and columns (observations) such that both modes lie on low-dimensional manifolds. Furthermore, the relationships between the rows may be informed by the relationships between the columns, and vice versa. Several recent papers (Gavish & Coifman, 2012; Ankenman, 2014; Mishne et al., 2016; Shahid et al., 2016; Mishne et al., 2017; Yair et al., 2017) exploit this coupled relationship to coorganize matrices and infer underlying row and column embeddings.

Further complicating the story is that such matrices may suffer from missing values, due to measurement corruptions and limitations. Missing values can sabotage efforts to learn the low dimensional manifold underlying the data. Specifically, kernel-based methods rely on calculating a similarity matrix between observations, whose eigendecomposition yields a new embedding of the data. As the number of missing entries grows, the distances between points are increasingly distorted, resulting in poor representation of the data in the low-dimensional space (Gilbert & Sonthalia, 2018) . Matrix completion algorithms assume the data is low-rank and fill in the missing values by fitting a global linear subspace to the data. Yet, this may fail when the data lies on a nonlinear manifold. The three components of our approach: 1) smooth estimates of a matrix with missing entries at multiple scales via co-clustering, 2) a multi-scale metric using the smooth estimates across all scales, yielding an affinity kernel between rows/columns, and 3) nonlinear embeddings of the rows and columns. The multiscale metric between two columns (red and orange) is a weighted Euclidean distance between those columns at multiple scales, given by solving the co-clustering for increasing values of the cost parameters γ r and γ c .

Manifold learning in the missing data scenario has been addressed by a few recent papers. Nonlinear Principle Component Analysis (NLPCA) (Scholz et al., 2005) uses an autoencoder neural network, where the middle layer serves to learn a low-dimensional embedding of the data, and the trained autoencoder is used to fill in missing values. Missing Data Recovery through Unsupervised Regression (Carreira-Perpin & Lu, 2011) first fills in the missing data with linear matrix completion methods, then calculates a non-linear embedding of the data and incorporates this embedding in an optimization problem to fill in the missing values. Recently Gilbert & Sonthalia (2018) proposed MR-MISSING which first calculates an initial distance matrix using only non-missing entries and then uses the increase-only-metric-repair (IOMR) method to fix the distance matrix so that it is a metric from which they calculate an embedding. None of these methods consider the co-manifold setting, where the coupled structure of the rows and the columns can be used to fill in the data, and to calculate an embedding.

In this paper, we introduce a new method for performing joint dimension reduction on the rows and columns of a data matrix, which we term co-manifold learning, in the missing data setting. We build on two recent lines of work on co-organizing the rows and columns of a data matrix (Gavish & Coifman, 2012; Mishne et al., 2016; 2017) and convex optimization methods for performing coclustering (Chi et al., 2017; 2018) . The former provide a flexible framework for jointly organizing rows and columns but lacks algorithmic convergence guarantees. The latter provides convergence guarantees but does not take full advantage of the multiple scales of the data revealed in the solution.

In the first stage of our approach, rather than inferring biclusters at a single scale, we use a multiscale optimization framework to fill in the data at fine to coarse scales while imposing smoothness on both the rows and the columns. The scales of the solutions are encoded in a pair of joint cost parameters along the rows and columns. Next, we define a new multi-scale metric based on the filled-in matrix across all scales, which we then use to calculate nonlinear embeddings of the rows and columns. Thus our approach yields three results: a collection of smoothed estimates of the matrix, pairwise distances on the rows and columns that better estimate the geometry of the complete data matrix, and corresponding nonlinear embeddings (see Figure 1 ). We will demonstrate in experimental results that our method reveals meaningful representations in coupled datasets with missing entries, whereas other methods are capable of revealing a meaningful representation only along one of the modes.

The paper is organized as follows. We present the optimization framework in Section 2, the new multi-scale metric for co-manifold learning in Section 3 and experimental results in Section 4.

2. CO-CLUSTERING AN INCOMPLETE DATA MATRIX:

We seek a collection of complete matrix approximations of a partially observed data matrix X ∈ R m×n that have been smoothed along their row and columns to varying degrees. This collection will serve in computing row and column multi-scale distances to better estimate the pairwise distances of the complete data matrix. Let [m] denote the set of indices {1, . . . , m}, and let Θ ⊆ [m] × [n] be a subset of the indices that correspond to observed entries of X, and let P Θ denote the projection operator of m×n matrices onto an index set Θ, i.e. [P Θ (X)] ij is x ij if (i, j) ∈ Θ and is 0 otherwise. We seek a minimizer U(γ r , γ c ) of the following function.

EQUATION

The quadratic term quantifies how well U approximates X on the observed entries, while the two roughness penalties, J r (U) and J c (U), incentivize smoothness across the rows and columns of U.

The nonnegative parameters γ r and γ c tune the tradeoff between how well U agrees with X over Θ and how smooth U is with respect to its rows and columns. By tuning γ r and γ c , we obtain estimates of X at varying scales of row and column smoothness.

In this paper, we use roughness penalties of the following forms

EQUATION

where U i• (U •i ) denotes the ith row (column) of the matrix U. The index sets E r and E c denote the edge sets of row and column graphs that encode a preliminary data-driven assessment of the similarities between rows and columns of the data matrix. The function Ω, which maps [0, ∞) into [0, ∞), will be explained shortly. Variations on the optimization problem of minimizing (1) have been previously proposed in the literature. When there is no data missing, i.e.

Θ = [m] × [n]

and Ω is a linear mapping, minimizing the objective in (1) produces a convex biclustering problem (Chi et al., 2017) . Additionally, if either γ r or γ c is zero, then we obtain convex clustering (Pelckmans et al., 2005; Hocking et al., 2011; Lindsten et al., 2011; Chi & Lange, 2015) . If we take Ω to be a nonlinear concave function, problem (1) reduces to an instance of concave penalized regressionbased clustering (Pan et al., 2013; Marchetti & Zhou, 2014; Wu et al., 2016) . The convergence properties of our co-clustering procedure will rely on the following two assumptions.

Assumption 2.1 The row and column graphs E r and E c are connected, i.e. the row graph is connected if for any pair of rows, indexed by i and j with i = j, there exists a sequence of indices i → k → • • • → l → j such that (i, k), . . . , (l, j) ∈ E r . A column graph is connected under analogous conditions.

Assumption 2.2 The function Ω : [0, ∞) → [0, ∞) is (i) concave and continuously differentiable on (0, ∞), (ii) vanishes at the origin, i.e. Ω(0) = 0, (iii) is increasing on [0, ∞), and (iv) has finite directional derivative at the origin.

Figure 4 in Appendix A plots an example of a function Ω that satisfies Assumption 2.2. For concreteness, in the rest of this paper, we use the following function Ω

EQUATION

where is a small positive number, e.g. 10 -12 . The key feature of Ω in (3), which satisfies Assumption 2.2, is that when it is used in the roughness penalties (2) small differences between rows and columns are penalized significantly more than larger differences resulting in more aggressive smoothing of small noisy variations and leaving intact more significant systematic variations. Specifically, functions that satisfy Assumption 2.2 are not differentiable at the origin and therefore incentivize sparsity in the differences in the rows and columns. Consequently, small noisy variations between pairs of rows and columns are eliminated completely for sufficiently large γ r and γ c . This is in contrast with commonly used quadratic penalties. For example, replacing J r (U) and J c (U) by quadratic row and column penalties

J r (U) = 1 2 (i,j)∈Er w ij U i• -U j• 2 2 and J c (U) = 1 2 (i,j)∈Ec wij U •i -U •j

gives a version of matrix completion on graphs (Kalofolias et al., 2014; Rao et al., 2015) , where w ij and wij are fixed row and column weights inferred from the data. Shahid et al. (2016) also use quadratic row and column penalties to perform joint linear dimension reduction on the rows and columns of a data matrix. Penalties like the ones given in (4), unlike the ones considered in Assumption 2.2, smooth out more significant systematic variations more aggressively and shrink, but do not completely eliminate, small noisy variations.

A simple concrete example illuminates the differences between roughness penalties considered in this paper and commonly used existing penalties. Let

U 1• = (1 0) T and U 2• = (1 + δ 0)

T and take Ω to be as defined in (3). Then,

U 1• -U 2• 2 = δ, Ω( U 1• -U 2• 2 ) ≈ √ δ, and U 1• -U 2• 2 2 = δ 2 .

Suppose there is a small difference between the first and second rows of U, e.g. δ = 10 -4 . Then

U 1• -U 2• 2 = 10 -4 , Ω( U 1• -U 2• 2 ) ≈ 10 -2 , and U 1• -U 2• 2 2 = 10 -8 .

Thus, small differences are penalized the most using the concave Ω and the least by the convex quadratic penalty. Suppose there is a large difference between the first and second rows of U, e.g. δ = 10 4 . Then

U 1• -U 2• 2 = 10 4 , Ω( U 1• -U 2• 2 ) ≈ 10 2 , and U 1• -U 2• 2 2 = 10 8 .

Thus, large differences are penalized the least using the concave Ω and the most by the convex quadratic penalty.

The practical consequence of the differences highlighted by the example above is that the convex penalties, either when Ω is linear or quadratic, do not introduce enough smoothing for small differences and too much smoothing for large differences. Indeed, (Pan et al., 2013; Marchetti & Zhou, 2014; Wu et al., 2016) showed that superior clustering results could be had, when only J r or J c is used, by using concave Ω. Chi et al. (2017) also showed that empirically that the solution to the convex biclustering problem tended to identify too many biclusters and consequently also introduced a one-step reweighted convex biclustering refinement, which recovered the true biclusters more accurately in simulation experiments. The reweighting refinement can be seen as taking a single step in an iterative algorithm for minimizing (1) inexactly when Ω is concave (Chi et al., 2018; Chi & Steinerberger, 2018) . As our method combines biclustering of incomplete data at different scales, or values of γ r and γ c , consequently we extend the reweighting refinement introduced by Chi et al.

(2017) in Section 2.1.

Our problem formulation (1) is distinct from related problem formulations in the following ways:

1. Rows and columns of U are simultaneously shrunk towards each other as the parameters γ r and γ c increase. Note that this shrinkage procedure is fundamentally different from methods like the clustered dendrogram, which independently cluster the rows and columns as well as alternating partition tree construction procedures (Gavish & Coifman, 2012; Mishne et al., 2016) .

2. Our ultimate goal is not to perform matrix completion (though this is a by-product of our approach) but rather to perform joint row and column dimension reduction. 

2.1. CO-CLUSTERING ALGORITHM:

We now introduce a majorization-minimization (MM) algorithm (Sun et al., 2017) for solving the minimization in (1). The basic strategy behind an MM algorithm is to convert a hard optimization problem into a sequence of simpler ones. The MM principle requires majorizing the objective Algorithm 1 CO-CLUSTER-MISSING(P Θ (X), γ r , γ c ) 1: Initialize U 0 , wr,ij , and wc,ij 2: repeat 3:

X ← P Θ (X) + P Θ c (U t ) 4: {U t+1 , n r , n c } ← CONVEX-BICLUSTER X, γ r , γ c , { wr,ij }, { wc,ij } 5: wr,ij ← Ω ( U t+1,i• -U t+1,j• 2 ) for all (i, j) ∈ E r 6: wc,ij ← Ω ( U t+1,•i -U t+1,•j 2 ) for all (i, j) ∈ E c 7: until convergence 8: Return U(γ r , γ c ) = U t , X, n r , n c function f (U) by a surrogate function g(U | Ũ) anchored at Ũ. Majorization is a combination of the tangency condition g(U | Ũ) = f ( Ũ) and the domination condition g(U | Ũ) ≥ f (U) for all U ∈ R m×n . The associated MM algorithm is defined by the iterates U t+1 = arg min U g(U | U t ).

It is straightforward to verify that the MM iterates generate a descent algorithm driving the objective function downhill, i.e. that f (U t+1 ) ≤ f (U t ) for all t.

The following function

g(U | Ũ) = 1 2 X -U 2 F + γ r (i,j)∈Er wr,ij U i• -U j• 2 + γ c (i,j)∈Ec wc,ij U •i -U •j 2 + κ

majorizes our objective function (1) at Ũ, where κ is a constant that does not depend on U and wr,ij and wc,ij are weights that depend on Ũ, i.e.

EQUATION

where Ω denotes the first derivative of Ω. We give a detailed derivation of this majorization in Appendix A.

Minimizing g(U | Ũ) is equivalent to minimizing the objective function of the convex biclustering problem for which efficient algorithms have been introduced (Chi et al., 2017) . Thus, in the t + 1th iteration, our MM algorithm solves a convex biclustering problem where the missing values in X have been replaced with the values of Ũ = U t and the weights wr,ij and wc,ij have been computed based on Ũ = U t according to (5). Note that the weights are continuously updated throughout the optimization as opposed to the fixed weights in (Chi et al., 2017) . This introduces a notion of the scale of the solution into the weights.

Algorithm 1 summarizes our MM algorithm, CO-CLUSTER-MISSING, which returns a smooth output matrix U(γ r , γ c ), a filled-in matrix X = P Θ (X) + P Θ c (U(γ r , γ c )) as well as n r and n c , which are respectively the number of distinct rows and distinct columns in U(γ r , γ c ). The CO-CLUSTER-MISSING algorithm has the following convergence guarantee.

Proposition 1 Under Assumption 2.1 and Assumption 2.2, the sequence U t generated by Algorithm 1 has at least one limit point, and all limit points are stationary points of (1).

The proof of Proposition 1 is presented in Appendix B.

2.2. CO-CLUSTERING AT MULTIPLE SCALES:

Initializing Algorithm 1 is very important as the objective function in (1) is not convex. The matrix U (0) is initialized to be the mean of all non-missing values. The connectivity graphs E r and E c are initialized at the beginning using k-nearest-neighbor graphs, and remain fixed throughout all considered scales. If we observed the complete matrix, employing a sparse Gaussian kernel is a natural way to quantify the local similarity between pairs of rows and pairs of columns. The challenge is that we do not have the complete data matrix X but only the partially observed one P Θ (X). Therefore, we rely only on the observed values to calculate the k-nearest-neighbor graph, based on the distance used by Ram et al. (2013) in an image inpainting problem.

Algorithm 2 Co-manifold learning on an Incomplete Data Matrix while n c > 1 do 6:

1: Initialize E r , E c 2: Set d(X •i , X •j ) = 0 and d(X •i , X •j ) = 0 3: Set n r = m, n c = n, k = k 0 ,

U (l,k) , X(l,k) , n r , n c ← CO-CLUSTER-MISSING P Θ (X), γ r = 2 l , γ c = 2 k 7: Update row distances: d (X i• , X j• ) += d X(l,k) i• , X(l,k) j• 8: Update column distances: d (X •i , X •j ) += d X(l,k) •i , X(l,k) •j 9: k ← k + 1 10:

end while 11:

l ← l + 1 12: end while 13: Calculate affinities A r (X i• , X j• ) and A c (X •i , X •j ) 14: Calculate embeddings Ψ r , Ψ c

Solving the co-clustering problem in Algorithm 1 at a single scale yields the smooth estimate U, the filled-in data matrix X, and n r and n c which are the number of distinct row and column clusters, respectively, identified at that scale through columns and rows merging in U. To obtain a collection of estimates at multiple scales, we solve the optimization problem for pairs of values for γ r , γ c set at logarithmic scale (Chi & Steinerberger, 2018) until we have converged to single global bi-cluster, i.e. n r = n c = 1.

We start with small values of γ r = 2 l0 and γ c = 2 k0 , where l 0 , k 0 < 0. We apply the co-clustering (Algorithm 1) and obtain the smooth estimate U (l0,k0) = U(2 l0 , 2 k0 ) used to fill in the data matrix X(l0,k0)

. Keeping γ r fixed, we continue increasing γ c by power of 2 and applying the biclustering until the algorithm converges to one cluster along the columns (n c = 1). We then increase γ r by power of 2 and reset γ c = 2 k0 . We repeat this procedure at increasing scales of γ r = 2 l , γ c = 2 k , until we have converged to a single global bicluster. This multiscale procedure yields a collection of filled-in matrices at all scales X(l,k)

l,k

. Note that the l and k denote the power of 2 taken for specific row and column cost parameters (γ c , γ r ) in the solution. This is intended as a compact notation that corresponds a pair of parameters (γ r , γ c ) to their solution U (l,k) and filled in estimate X(l,k) .

3. CO-MANIFOLD LEARNING:

Kernel-based manifold learning relies on constructing a good similarity measure between points, and a dimension reduction method based on this similarity. The eigenvectors of these kernels are typically used as the new low-dimensional coordinates for the data. Here we leverage having calculated an estimate of the filled-in matrix at multiple scales

X(l,k) l,k

, to define a new metric between rows and columns. This metric will encompass all bi-scales as defined by joint pairs of optimization cost parameters γ r , γ c . Given a new metric we employ Diffusion maps (Coifman & Lafon, 2006) to obtain a new embedding of the rows and columns. The full algorithm is given in Algorithm 2.

3.1. MULTI-SCALE METRIC:

We define a new metric to estimate the geometry of the complete data matrix both locally and globally. For a given pair γ r , γ c , we calculate the Euclidean distance between rows for the filled-in matrix at that joint scale, weighted by the cost parameters:

d X(l,k) i• , X(l,k) j• = (γ r γ c ) α X(l,k) i• - X(l,k) j• 2

where X(l,k) = P Θ (X) + P Θ c (U (l,k) ), and we set α = -1/2 to favor local over global structure in our simulations. Our goal is to aggregate distances between a pair of rows (columns) across multiple scales of the solution, to calculate a metric that better recovers the local and global geometry of the data despite the missing values, thus "fixing" the missing data metric. Having solved for multiple pairs from the solution surface, we sum over all the distances to obtain a multi-scale distance on the data rows:

d(X i• , X j• ) = l,k d X(l,k) i• , X(l,k) j•

.

An analogous multi-scale distance is computed for pairs of columns. Note that if there are no missing values, this metric is just the Euclidean pairwise distance scaled by a scalar, so that we recover the original embedding of the complete matrix.

This metric takes advantage of solving the optimization for multiple pairs of cost parameters and filling in the missing values with increasingly smooth estimates (as γ r and γ c increase). It also alleviates the need to identify the ideal scale at which to fill in the points; it is not clear that a single "optimal" scale actually exists, but rather different points in the matrix may have different optimal scales. As opposed to the partition-tree based metric of Mishne et al. (2017) , this metric takes into account all joint scales of the data as the matrix U is smoothed across rows and columns simultaneously, thus fully taking advantage of the coupling between both modes.

3.2. DIFFUSION MAPS:

Having calculated a multi-scale metric on the rows and columns throughout the joint optimization procedure, we can now construct a pair of low-dimensional embeddings based on these distances. Specifically we use Diffusion maps (Coifman & Lafon, 2006) , but any dimension reduction technique relying on the construction of a distance kernel could be used instead. We briefly review the construction of the diffusion maps for the rows (features) of a matrix but the same can be applied to the columns (observations). Given a distance between two rows of the matrix d(X i• , X j• ), we construct an affinity kernel on the rows. We choose an exponential function, but other kernels can be considered depending on the application:

A[i, j] = exp{-d 2 (X i• , X j• )/σ 2 },

where σ is a scale parameter. The exponential function enhances locality, as pairs of samples whose distance exceed σ have negligible affinity. One possible choice for σ is to be the median of distances within the data.

We derive a row-stochastic matrix P by normalizing the rows of matrix A: P = D -1 A, where D is a diagonal matrix whose elements are given by D[i, i] = j A[i, j]. The eigendecomposition of P yields a sequence of positive decreasing eigenvalues: 1 = λ 0 ≥ λ 1 ≥ . . ., and right eigenvectors {ψ } . Retaining only the first d eigenvalues and eigenvectors, the mapping Ψ embeds the rows into the Euclidean space R d :

Ψ : X i• → λ 1 ψ 1 (i), λ 2 ψ 2 (i), . . . , λ d ψ d (i) T .

The embedding integrates the local connections found in the data into a global representation, which enables visualization of the data, organizes the data into meaningful clusters, and identifies outliers and singular samples. This embedding is also equipped with a noise-robust distance, the diffusion distance. For more details on diffusion maps, see Coifman & Lafon (2006) .

4. NUMERICAL EXPERIMENTS:

The model we consider in the paper is such that the data is not represented by a biclustering model but rather at least one of the modes (rows/columns) lies on a low-dimensional manifold. In our experiments we consider three such examples. In the first a manifold structure exists along both rows and columns, and for the second and third the columns belong to disjoint clusters while the rows lie on a manifold:

• linkage A synthetic dataset with a one-dimensional manifold along the rows and a twodimensional manifold along the columns. Let {z i } N1 i=1 ∈ R 3 be points along a helix and let {y j } N2 j=1 ∈ R 3 be a two dimensional surface, where we set N 1 = 190, N 2 = 300. We analyze the matrix of Euclidean distances between the two spatially distant sets of points to reveal the underlying geometry of both rows and columns,

EQUATION

)

Other functions of the distance can also be used such as the elastic or Coulomb potential operator (Coifman & Gavish, 2011) . Missing values correspond to having access to only some of the distances between pairs of points across the two sets. Note that this is unlike MDS as we do not have pairwise distances between all datapoints, but rather distances between two sets of points with different geometries. • linkage2 A synthetic dataset with a clustered structure along the rows and a twodimensional manifold along the columns. Let {x i } N1 i=1 ∈ R 3 be composed of points in 3 Gaussian clouds in 3D and let {y j } N2 j=1 ∈ R 3 be a two dimensional surface as before, where we set N 1 = 200, N 2 = 300.

• lung500 A real-world dataset composed of 56 lung cancer patients and their gene expression (Lee et al., 2010) . We selected the 500 genes with the greatest variance from the original collection of 12,625 genes. Subjects belong to one of four subgroups; they are either normal subjects (Normal) or have been diagnosed with one of three types of cancers: pulmonary carcinoid tumors (Carcinoid), colon metastases (Colon), and small cell carcinoma (Small Cell).

For all datasets, the rows and columns of the data matrix are randomly permuted so their natural order does not play a role in inferring the geometry. We evaluate results both qualitatively and quantitatively. of missing values for our multiscale approach (blue), diffusion maps of missing data matrix (red), NLPCA (yellow) and FRPCAG for two different cost parameter settings (green and purple). We evaluate using the adjusted Rand Index (ARI) compared to the ground-truth labels of (left) the 4 cancer types for the lung500 dataset, and (right) 3 Gaussian clusters of the linkage2 dataset.

We compare our embeddings to three approaches: NLPCA with missing data completion (Scholz et al., 2005) , Fast Robust PCA on Graphs (FRPCAG) (Shahid et al., 2016) and Diffusion maps (DM) (Coifman & Lafon, 2006) on the missing data. FRPCAG provides a linear embedding of a low-rank estimate of a data matrix using quadratic row and column penalties, while the other methods output nonlinear embeddings. NLPCA and DM are applied to each mode separately, while our method and FRPCAG take into account the coupled geometry. Comparing to Diffusion maps demonstrates how missing values corrupt the embedding. Note that this is also equivalent to applying our approach for only a single scale of the cost parameters (γ r , γ c → ∞), as for this choice of parameters the solution U converges to the grand mean of the data.

In Figure 2 , we display the embeddings of the different methods for each of the three datasets for both their rows (top) and their columns (bottom), where 50% of the entries have been removed. Both NLPCA and DM reveal the underlying 2D surface structure on the rows in only one of the linkage datasets, and err greatly on the other. DM correctly infers a 1D path for the linkage dataset but it is increasingly noisy. For NLPCA the 1D embedding is not as smooth and clean as the embedding inferred by the co-manifold approach. Our method reveals the 2D surface in both cases. FRPCAG is unsuccessful in uncovering the non-linear manifolds underlying either rows or columns in the linkage dataset. For linkage2, the rows do not separate cleanly into disjoint clusters, and for the columns only one of the parameters of the 2D plane is uncovered.

For the lung500 data, NLPCA and DM embed the cancer samples such that the normal subjects (yellow) are close to the Colon type (cyan), whereas both our method and FRPCAG separates the normal subjects from the cancer types. This is due to taking into account the coupled structure of the genes and the samples. As opposed to the clustered structure along the samples, the three nonlinear methods reveal a smooth manifold structure to the genes, which is different than the assumed clustered structure a biclustering method would infer. The representation yielded by FRPCAG is the least structured but also does not reveal disjoint gene clusters. For plots presenting the datasets and filled-in values at multiple scales see Appendix C.

Manifold learning is not used only for data visualization but also for calculating new data representations that can then be used for signal processing and machine learning tasks. Here we calculate clustering accuracy of clustering the low-dimensional representation of the data. We apply k-means to the column embeddings of each method, with k set to the correct number of clusters in the data, as we want to evaluate the ability of the methods to properly represent the data without being sensitive to empirical estimation of the number of clusters in the data. We use the Adjusted Rand Index (ARI) (Hubert & Arabie, 1985) , to measure the similarity between the k-means clustering of the embedding and the ground-truth labels.

We note that Shahid et al. (2016) do not provide guidelines by which to select the appropriate cost parameters γ r and γ c in the FRPCAG solution. Therefore we calculated the solution using a wide range of possible values and plot the results for worst and best performance to demonstrate the sensitivity to setting these parameters. In contrast, our approach alleviates the need to perform this parameter selection.

The left panel of Figure 3 compares clustering the embedding of the cancer patients in lung500 by each method for increasing percentage of missing values in the data, where we averaged over 30 realizations of missing entries. For low values of missing data, FRPCAG performs the best but its performance degrades as the percentage of missing values increases. For higher values of missing data (50% and above), our embedding (blue plot) gives the best clustering result and its performance is only slightly degraded by increasing the percentage of missing values, as opposed to Diffusion maps (red plot). This demonstrates that the metric we calculate is a good estimate of the metric of the complete data matrix. NLPCA (yellow plot) performs worst.

The right panel of Figure 3 compares clustering the embedding of the three Gaussian clusters in link-age2 for increasing percentage of missing values in the data, where we averaged over 30 realizations of the data itself and the missing entries. Our embedding (blue plot) gives the best clustering result (for 20% missing values and above) and its performance is unaffected by increasing the percentage of missing values up to 80%, as opposed to Diffusion maps (red plot) which is greatly degraded by the missing values. NLPCA (yellow plot) does not perform as well as our approach, with performance decreasing as the percentage of missing values increases. For a good parameter selection FRPCAG (purple plot) performs well for low percentage of missing values, and then performance is drastically impaired above 50%. For a poor parameter selection FRPCAG (green plot) has the worst performance for practically all percentages of missing values. Note that the overall poor performance of FRPCAG in this setting is due to the underlying manifold being non-linear and the data being high-rank, demonstrating the need for nonlinear embedding approaches.

5. CONCLUSIONS:

In this paper we presented a new method for learning nonlinear manifold representations of both the rows and columns of a matrix with missing data. We proposed a new optimization problem to obtain a smooth estimate of the missing data matrix, and solved this problem for different values of the cost parameters, which encode the smoothness scale of the estimate along the rows and columns. We leverage calculating these multi-scale estimates into a new metric that aims to capture the geometry of the complete data matrix. This metric is then used in a kernel-based manifold learning technique to obtain new representations of both the rows and the columns. In future work we will investigate additional metrics in a general co-manifold setting and relate them to optimal transport problem and Earth Mover's Distance (Coifman & Leeb, 2013) . We also intend to develop efficient solutions to accelerate the optimization in order to address large-scale datasets, in addition to the small-scale regime we demonstrate here. We note that the datasets considered here, while being small-scale in the observation domain are high-dimensional in the feature domain, which is a non-trivial setting, and indeed a challenge for supervised methods such as deep learning due to limited training data. 

A DERIVATION OF MAJORIZATION:

We first construct a majorization of the data-fidelity term. It is easy to verify that the following function of U

EQUATION

where X = P Θ (X) + P Θ c ( Ũ), majorizes the data-fidelity term 1 2 P Θ (X) -P Θ (U) 2 F at Ũ. We next construct a majorization of the penalty term. Recall that the first-order Taylor approximation of a differentiable concave function provides a global upper bound on the function. Therefore, under Assumption 2.2, we have the following inequality Ω(z) ≤ Ω(z) + Ω (z)(z -z), for all z, z ∈ [0, ∞).

Figure 4 shows the relationship between Ω given in (3) with = 10 -12 and its first-order Taylor approximation at z = 1.

Thus, we can majorize the penalty term γ r J r (U) + γ c J c (U) with the function

EQUATION

where κ is a constant that does not depend on U and wr,ij and wc,ij (5) are weights that depend on Ũ. The sum of functions ( 7) and ( 8)

g(U | Ũ) = g 1 (U | Ũ) + g 2 (U | Ũ) (9) = 1 2 X -U 2 F + γ r (i,j)∈Er wr,ij U i• -U j• 2 + γ c (i,j)∈Ec wc,ij U •i -U •j 2 + κ

majorizes our objective function (1) at Ũ.

B PROOF OF PROPOSITION 1:

The MM algorithm generates a sequence of iterates that has at least one limit point, and the limit points are stationary points of the objective function

EQUATION

Fix an arbitrary direction vector v. The inner product of v with an element in the set on right hand side of ( 11) is given by P Θ (u -x), v + γ r (i,j)∈Er

EQUATION

where

d ij ∈ ∂ L ij u 2 and dij ∈ ∂ Lij u 2 .

Then the supremum of the right hand side of ( 12) over all d ij ∈ ∂ L ij u 2 and dij ∈ ∂ Lij u 2 is nonnegative, because 0 ∈ ∂g(u | u). Consequently, all fixed points of ψ are also stationary points of f .

B.3 THE MM ITERATE SEQUENCE HAS A LIMIT POINT:

To ensure the existence of a limit point, we show that the function f is coercive, i.e. f (U t ) → ∞ for any sequence U t F → ∞. Recall that according to Assumption 2.1 we assume that the row and column edge sets E r and E c form connected graphs. Therefore, J r (U) = J c (U) = 0 if and only if U = a11 T (Chi et al., 2017, Proposition 3) . The edge-incidence matrix of the column graph Φ c ∈ R encodes its connectivity and is defined as

φ c,li =    1

If node i is the head of edge l, -1 If node i is the tail of edge l, 0 otherwise.

The row edge-incidence matrix Φ r ∈ R |Er|×m is defined similarly. Assume that Θ non-empty, i.e. at least one entry of the matrix has been observed. Finally, assume that Ω is also coercive.

Note that any sequence U t = a t 11 T + B t where B t , 11 T = 0. Note that J r (U t ) = J r (B t ) and J c (U t ) = J c (B t ). Let U t be a diverging sequence, i.e. U t F → ∞. There are two cases to consider.

Case I: Suppose that B t F → ∞. Let

L = I ⊗ Φ r Φ c ⊗ I ∈ R |Er|m+|Ec|n×mn ,

and let σ min denote the smallest singular value of L. Note that the null space of L is the span of 1. Therefore, since 1, b t = 0

EQUATION

Also note that

Lb t = vec(Φ r B t ) vec(B t Φ T c )

.

Since the mapping x = x T 1

x T 2 T → max{ x 1 2 , x 2 2 } is a norm, and all finite dimensional norms are equivalent, there exists some η > 0 such that X(-3,-2) , X(1,0) , X(5,2) . Rows and columns have been reordered based on the manifold embedding following (Ankenman, 2014) . 

EQUATION

:

, (4)

Figures and Tables:

Figure1: The three components of our approach: 1) smooth estimates of a matrix with missing entries at multiple scales via co-clustering, 2) a multi-scale metric using the smooth estimates across all scales, yielding an affinity kernel between rows/columns, and 3) nonlinear embeddings of the rows and columns. The multiscale metric between two columns (red and orange) is a weighted Euclidean distance between those columns at multiple scales, given by solving the co-clustering for increasing values of the cost parameters γ r and γ c .

3. Our work generalizes bothShahid et al. (2016) andChi et al. (2017) in that we seek the flexibility of performing non-linear dimension reduction on the rows and columns of the data matrix, i.e. a more general manifold organization than a co-clustered structure.4. Instead of determining an optimal single scale of the solution as inShahid et al. (2016);Chi et al. (2017), we recognize that the multiple scales of the different solutions can be aggregated to better estimate the underlying geometry, similar to the tree-based Earth mover's distance proposed in Ankenman (2014);Mishne et al. (2017).

and l = l 0 4: while n r > 1 do 5:

Figure2: Comparing row and column embeddings of NLPCA, FRPCAG, DM, Ours, for three datasets with 50% missing entries. For each dataset, top / bottom row is embedding of rows / columns of X. For the lung500 dataset, the color of the clusters are as follows: yellow -normal subjects, dark blue -carcinoid, cyan -colon , red -small cell carcinoma.

Figure3: Comparing k-means clustering applied to embedding of data for increasing percentages of missing values for our multiscale approach (blue), diffusion maps of missing data matrix (red), NLPCA (yellow) and FRPCAG for two different cost parameter settings (green and purple). We evaluate using the adjusted Rand Index (ARI) compared to the ground-truth labels of (left) the 4 cancer types for the lung500 dataset, and (right) 3 Gaussian clusters of the linkage2 dataset.

Figure 4: Majorization of the Ω function (black) given in (3) by its first-order Taylor approximation at 1 (blue).

Figure 5: Points in 3D used to generate the Euclidean distance matrix X in the linkage dataset. Rows correspond to the helix, columns to the 2D surface. Points are colored corresponding to the embedding of rows and columns in Figure 2.

Figure 6: linkage dataset: (Left) Complete matrix X. (Right) Matrix whose rows and columns and columns have been permuted and 50% of the values have been removed.

Figure 7: linkage dataset: Filled-in matrices X at multiple scales:X(-3,-2) , X(1,0) , X(5,2) . Rows and columns have been reordered based on the manifold embedding following(Ankenman, 2014).

Figure8: Points in 3D used to generate the Euclidean distance X in the linkage2 dataset. Rows correspond to the three 3D Gaussians, columns to the 2D surface. Points are colored corresponding to the embedding of rows and columns in Figure2

Peer review (iclr2019):
 This paper presents a joint learning method for filling missing value and bi-clustering. The method extends (Chi et al. 2017), using a penalized matrix approximation. The proposed method is tested on three data sets, where two are synthetic and one small real-world data matrix. The presented method is claimed to be better than two classical approaches Nonlinear PCA and Diffusion Maps. 1) Filling missing values is not new. Even co-clustering with missing values also exists. It is insufficient to defeat two methods which are older than ten years. More extensive comparison is needed but lacking here. Why not first use a dedicated method such as MICE or collaborative filtering, and then run embedding method on rows and columns? 2) The purpose of the learning is unclear. The title does not give any hint about the learning goal. The objective function reads like filling missing values. The subsequent text claims that minimizing such a objective can achieve biclustering. However, in the experiment, the comparison is done via visualization and normal clustering (k-means). 3) The empirical results are not convincing. Two data sets are synthetic. The only real-world data set is very small. Why k-means was used? How to choose k in k-means? 4) The choice Omega function after Proposition 1 needs to be elaborated. A function curve plot could also help. 5) What is Omega* in Eq. 2?