Year: 2018

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.

heading: 1 INTRODUCTION
text: 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.
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.

heading: 2 CO-CLUSTERING AN INCOMPLETE DATA MATRIX
text: We seek a collection of complete matrix approximations of a partially observed data matrix X ∈ Rm×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 ofm×nmatrices onto an index set Θ, i.e. [PΘ(X)]ij is xij if (i, j) ∈ Θ and is 0 otherwise. We seek a minimizer U(γr, γc) of the following function.
f(U; γr, γc) = 1
2 ‖PΘ(X)− PΘ(U)‖2F + γrJr(U) + γcJc(U). (1)
The quadratic term quantifies how well U approximates X on the observed entries, while the two roughness penalties, Jr(U) and Jc(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 Jr(U) = ∑
(i,j)∈Er
Ω(‖Ui· −Uj·‖2) and Jc(U) = ∑
(i,j)∈Ec
Ω(‖U·i −U·j‖2), (2)
where Ui· (U·i) denotes the ith row (column) of the matrix U. The index sets Er and Ec 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 Er and Ec are connected, i.e. the row graph is connected if for any pair of rows, indexed by i and j with i 6= j, there exists a sequence of indices i → k → · · · → l → j such that (i, k), . . . , (l, j) ∈ Er. A column graph is connected under analogous conditions.
Assumption 2.2 The function Ω : [0,∞) 7→ [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 Ω
Ω(z) = 1
2 ∫ z 0 1√ ζ + dζ, (3)
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 Jr(U) and Jc(U) by quadratic row and column penalties
Jr(U) = 1
2 ∑ (i,j)∈Er wij‖Ui· −Uj·‖22 and Jc(U) = 1 2 ∑ (i,j)∈Ec w̃ij‖U·i −U·j‖22, (4)
gives a version of matrix completion on graphs (Kalofolias et al., 2014; Rao et al., 2015), where wij and w̃ij 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 U1· = (1 0) T and U2· = (1 + δ 0) T and take Ω to be as defined in (3). Then,
‖U1· −U2·‖2 = δ, Ω(‖U1· −U2·‖2) ≈ √ δ, and ‖U1· −U2·‖22 = δ2.
Suppose there is a small difference between the first and second rows of U, e.g. δ = 10−4. Then
‖U1· −U2·‖2 = 10−4, Ω(‖U1· −U2·‖2) ≈ 10−2, and ‖U1· −U2·‖22 = 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. δ = 104. Then
‖U1· −U2·‖2 = 104, Ω(‖U1· −U2·‖2) ≈ 102, and ‖U1· −U2·‖22 = 108.
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 Jr or Jc 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.
3. Our work generalizes both Shahid et al. (2016) and Chi 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 in Shahid 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).

heading: 2.1 CO-CLUSTERING ALGORITHM
text: 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 U0, w̃r,ij , and w̃c,ij 2: repeat 3: X̃← PΘ(X) + PΘc(Ut) 4: {Ut+1, nr, nc} ← CONVEX-BICLUSTER ( X̃, γr, γc, {w̃r,ij}, {w̃c,ij}
) 5: w̃r,ij ← Ω′(‖Ut+1,i· −Ut+1,j·‖2) for all (i, j) ∈ Er 6: w̃c,ij ← Ω′(‖Ut+1,·i −Ut+1,·j‖2) for all (i, j) ∈ Ec 7: until convergence 8: Return { U(γr, γc) = Ut, X̃, nr, nc }
function f(U) by a surrogate function g(U | Ũ) anchored at Ũ. Majorization is a combination of the tangency condition g(U | Ũ) = f(Ũ) and the domination condition g(U | Ũ) ≥ f(U) for all U ∈ Rm×n. The associated MM algorithm is defined by the iterates Ut+1 = arg min
U g(U | Ut).
It is straightforward to verify that the MM iterates generate a descent algorithm driving the objective function downhill, i.e. that f(Ut+1) ≤ f(Ut) for all t. The following function
g(U | Ũ) = 1 2 ‖X̃−U‖2F + γr ∑ (i,j)∈Er w̃r,ij‖Ui· −Uj·‖2 + γc ∑ (i,j)∈Ec w̃c,ij‖U·i −U·j‖2 + κ
majorizes our objective function (1) at Ũ, where κ is a constant that does not depend on U and w̃r,ij and w̃c,ij are weights that depend on Ũ, i.e.
w̃r,ij = Ω ′(‖Ũi· − Ũj·‖2) and w̃c,ij = Ω′(‖Ũ·i − Ũ·j‖2), (5)
where Ω′ denotes the first derivative of Ω. We give a detailed derivation of this majorization in Appendix A.
Minimizing g(U | 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̃ = Ut and the weights w̃r,ij and w̃c,ij have been computed based on Ũ = Ut 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 nr and nc, which are respectively the number of distinct rows and distinct columns in U(γr, γc). The CO-CLUSTERMISSING algorithm has the following convergence guarantee.
Proposition 1 Under Assumption 2.1 and Assumption 2.2, the sequence Ut 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.

heading: 2.2 CO-CLUSTERING AT MULTIPLE SCALES
text: 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 Er and Ec 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 1: Initialize Er, Ec 2: Set d(X·i,X·j) = 0 and d(X·i,X·j) = 0 3: Set nr = m,nc = n, k = k0, and l = l0 4: while nr > 1 do 5: while nc > 1 do 6: { U(l,k), X̃ (l,k) , nr, nc } ← CO-CLUSTER-MISSING ( PΘ(X), γr = 2l, γc = 2k
) 7: Update row distances: d (Xi·,Xj·) += 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 Ar(Xi·,Xj·) and Ac(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 nr and nc 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. nr = nc = 1.
We start with small values of γr = 2l0 and γc = 2k0 , where l0, k0 < 0. We apply the co-clustering (Algorithm 1) and obtain the smooth estimate U(l0,k0) = U(2l0 , 2k0) 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 (nc = 1). We then increase γr by power of 2 and reset γc = 2k0 . We repeat this procedure at increasing scales of γr = 2l, γc = 2k, 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) .

heading: 3 CO-MANIFOLD LEARNING
text: 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 cal-
culated 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.

heading: 3.1 MULTI-SCALE METRIC
text: 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(Xi·,Xj·) = ∑ 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.

heading: 3.2 DIFFUSION MAPS
text: 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(Xi·,Xj·), 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{−d2(Xi·,Xj·)/σ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−1A, 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 Rd:
Ψ : Xi· → ( λ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).

heading: 4 NUMERICAL EXPERIMENTS
text: 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 {zi}N1i=1 ∈ R3 be points along a helix and
let {yj}N2j=1 ∈ R3 be a two dimensional surface, where we set N1 = 190, N2 = 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,
X[i, j] = ‖zi − yj‖2. (6) 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 {xi}N1i=1 ∈ R3 be composed of points in 3 Gaussian clouds in 3D and let {yj}N2j=1 ∈ R3 be a two dimensional surface as before, where we set N1 = 200, N2 = 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.
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 linkage2 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.

heading: 5 CONCLUSIONS
text: 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.

heading: A DERIVATION OF MAJORIZATION
text: We first construct a majorization of the data-fidelity term. It is easy to verify that the following function of U
g1(U | Ũ) = 1
2 ‖X̃−U‖2F, (7)
where X̃ = PΘ(X) + PΘc(Ũ), majorizes the data-fidelity term 12‖PΘ(X)− PΘ(U)‖ 2 F at Ũ.
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 γrJr(U) + γcJc(U) with the function g2(U | Ũ) = γr ∑
(i,j)∈Er
w̃r,ij‖Ui· −Uj·‖2 + γc ∑
(i,j)∈Ec
w̃c,ij‖U·i −U·j‖2 + κ, (8)
where κ is a constant that does not depend on U and w̃r,ij and w̃c,ij (5) are weights that depend on Ũ. The sum of functions (7) and (8)
g(U | Ũ) = g1(U | Ũ) + g2(U | Ũ) (9)
= 1
2 ‖X̃−U‖2F + γr ∑ (i,j)∈Er w̃r,ij‖Ui· −Uj·‖2 + γc ∑ (i,j)∈Ec w̃c,ij‖U·i −U·j‖2 + κ
majorizes our objective function (1) at Ũ.

heading: B PROOF OF PROPOSITION 1
text: 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
f(U) = 1
2 ‖PΘ(X)− PΘ(U)‖2F + γrJr(U) + γcJc(U). (10)
To reduce notational clutter, we suppress the dependency of f on γr and γc since they are fixed during Algorithm 1. We prove Proposition 1 in three stages. First, we show that all limit points of the MM algorithm are fixed points of the MM algorithm map. Second, we show that fixed points of the MM algorithm are stationary points of f in (10). Finally, we show that the MM algorithm has at least one limit point.
B.1 LIMIT POINTS ARE FIXED POINTS
The convergence theory of MM algorithms relies on the properties of the algorithm map ψ(U) that returns the next iterate given the last iterate. For easy reference, we state a simple version of Meyer’s monotone convergence theorem Meyer (1976), which is instrumental in proving convergence in our setting.
Theorem 1 Let f(U) be a continuous function on a domain S and ψ(U) be a continuous algorithm map from S into S satisfying f(ψ(U)) < f(U) for all U ∈ S with ψ(U) 6= U. Then all limit points of the iterate sequence Uk = ψ(Uk−1) are fixed points of ψ(U).
In order to apply Theorem 1, we need to identify elements in the assumption with specific functions and sets corresponding to the problem of minimizing (10). Throughout the following proof, it will sometimes be convenient to work with the column major vectorization of a matrix. The vector b = vec(B) is obtained by stacking the columns of B on top of each other.
The function f : Take S = Rm×n and f : S 7→ R to be the objective function in (10) and majorize f with g(U | Ũ) given in (9). The function f is continuous. Let ψ(Ũ) = arg min
U g(U | Ũ) denote
the algorithm map for the MM algorithm. Since g(U | Ũ) is strongly convex in U, it has a unique global minimizer. Consequently, f(ψ(U)) < f(U) for all ψ(U) 6= U. Continuity of the algorithm map ψ: Continuity of ψ follows from the fact that the solution to the convex biclustering problem is jointly continuous in the weights and data matrix (Chi et al., 2017)[Proposition 2]. The weight w̃r,ij(Ũ) = Ω′(‖Ui· − Uj·‖2) is a continuous function of Ũ, since Ω′ is continuous according to Assumption 2.2. The weight w̃c,ij(Ũ) is likewise continuous in Ũ. The data matrix passed into the convex biclustering algorithm is X̃ = PΘ(X) +PΘc(Ũ), which is a continuous function of Ũ since the projection mapping PΘc is continuous.
B.2 FIXED POINTS ARE STATIONARY POINTS
Let Lij = (ei− ej)T⊗ I and L̃ij = I⊗ (ei− ej)T, where ⊗ denotes the Kronecker product. Then
vec(Ui· −Uj·) = Liju and vec(U·i −U·j) = L̃iju. The directional derivative of f in the direction v at a point u is given by
Ω′(‖Liju‖2; v) = { Ω′(‖Liju‖2)〈Lijv, Liju‖Liju‖2 〉 Liju 6= 0 Ω′(‖Liju‖2)‖Lijv‖2 otherwise.
A point u is a stationary point of f , if for all direction vectors v 0 ≤ 〈PΘ(u− x),v〉+ γr ∑
(i,j)∈Er
Ω′(‖Liju‖2; v) + γc ∑
(i,j)∈Ec
Ω′(‖L̃iju‖2; v),
where PΘ(u− x) = vec(PΘ(U)− PΘ(X)). A point u is a fixed point of ψ, if 0 is in the subdifferential of g(u | u), i.e.
0 ∈ {PΘ(u− x)}+ γr ∑
(i,j)∈Er
Ω′(‖Liju‖2)∂‖Liju‖2 + γc ∑
(i,j)∈Ec
Ω′(‖L̃iju‖2)∂‖L̃iju‖2, (11)
where the set on the right is the subdifferential ∂g(u | u). If Liju 6= 0, then ∂‖Liju‖2 = {
LTij Liju ‖Liju‖2
} . On the other hand, if Liju = 0, then ∂‖Liju‖2 =
∂‖0‖2 = {d : ‖d‖2 ≤ 1}.
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
Ω′(‖Liju‖2)〈dij ,v〉+ γc ∑
(i,j)∈Ec
Ω′(‖L̃iju‖2)〈dij ,v〉, (12)
where dij ∈ ∂‖Liju‖2 and d̃ij ∈ ∂‖L̃iju‖2.
Then the supremum of the right hand side of (12) over all dij ∈ ∂‖Liju‖2 and d̃ij ∈ ∂‖L̃iju‖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(Ut) → ∞ for any sequence ‖Ut‖F → ∞. Recall that according to Assumption 2.1 we assume that the row and column edge sets Er and Ec form connected graphs. Therefore, Jr(U) = Jc(U) = 0 if and only if U = a11T (Chi et al., 2017, Proposition 3). The edge-incidence matrix of the column graph Φc ∈ R|Ec|×n 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 Ut = at11T + Bt where 〈Bt,11T〉 = 0. Note that Jr(Ut) = Jr(Bt) and Jc(Ut) = Jc(Bt). Let Ut be a diverging sequence, i.e. ‖Ut‖F → ∞. There are two cases to consider.
Case I: Suppose that ‖Bt‖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,bt〉 = 0
‖Lbt‖2 ≥ σmin‖Bt‖F. (13)
Also note that
Lbt =
( vec(ΦrBt)
vec(BtΦ T c )
) .
Since the mapping x = ( xT1 x T 2 )T 7→ max{‖x1‖2, ‖x2‖2} is a norm, and all finite dimensional norms are equivalent, there exists some η > 0 such that
η‖Lbt‖2 ≤ max { ‖ΦrBt‖F, ‖BtΦTc ‖F } . (14)
By the triangle inequality
max { ‖ΦrBt‖F, ‖BtΦTc ‖F } ≤ max  ∑ (i,j)∈Er ‖Lijbt‖2, ∑ (i,j)∈Ec ‖L̃ijbt‖2  . (15) Let M = max{|Er|, |Ec|} then
max  ∑ (i,j)∈Er ‖Lijbt‖2, ∑ (i,j)∈Ec ‖L̃ijbt‖2  ≤M max { max (i,j)∈Er ‖Lijbt‖2, max (i,j)∈Ec ‖L̃ijbt‖2 } . (16)
Putting inequalities (13), (14), (15), and (16) together gives us
ησmin M ‖Bt‖F ≤ max
{ max
(i,j)∈Er ‖Lijbt‖2, max (i,j)∈Ec ‖L̃ijbt‖2
} . (17)
Since Ω is increasing according to Assumption 2.2, it follows that Ω (ησmin
M ‖Bt‖F
) ≤ max { Ω ( max
(i,j)∈Er ‖Lijbt‖2
) ,Ω ( max
(i,j)∈Ec ‖L̃ijbt‖2
)} . (18)
Inequality (18) implies that min{γr, γc}MΩ (ησmin
M ‖Bt‖F
) ≤ min{γr, γc}max {Jr(Ut), Jc(Ut)}
≤ γrJr(Ut) + γcJc(Ut).
Consequently, since Ω is increasing and ‖Bt‖F →∞ implies that f(Ut)→∞. Case II: Suppose ‖Bt‖F ≤ B for some B. Then |at| → ∞. Note that we have the following inequality
f(Ut) ≥ ∑
(i,j)∈Θ
(xij − bk,ij − at)2
≥ ∑
(i,j)∈Θ
a2t − 2at(xij − bk,ij)
= |Θ|a2t − 2at ∑
(i,j)∈Θ
(xij − bk,ij)
≥ |Θ|a2t − 2at sup ‖Bt‖F≤B ∑ (i,j)∈Θ (xij − bk,ij)
= |Θ| [ a2t − 2atC ] = |Θ| [ (at − C)2 − C2 ] ,
where C = |Θ|−1 sup ‖Bt‖F≤B
∑ (i,j)∈Θ(xij − bk,ij).
The function (at − C)2 diverges since |at| → ∞. Therefore, the function f is coercive.

heading: C FILLING IN MISSING DATA
text: We present the original underlying structure of 3D points used to generate the Euclidean distance matrix X for the datasets linkage and linkage2 in Figure 5 and Figure 8. In Figure 6 and Figure 9, on the left we plot the original complete matrix where the rows and columns have been ordered according to the geometry of the 3D points. On the right we plot the matrix we analyze whose rows and columns have been permuted and 50% of the entries have been removed. In Figure 7 and Figure 10 we display the matrix X̃ (l,k)
for three pairs of values l, k to demonstrate the smoothing that is occurring across the different scales of the rows and columns.

Review summaries:
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.