{'fraction_non_alphanumeric': 0.036771474250278574, 'fraction_numerical': 0.009980136621287728, 'mean_word_length': 4.489893617021276, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Beyond kNN: Adaptive, Sparse Neighborhood Graphs via Optimal Transport

Nearest neighbour graphs are widely used to capture the geometry or topology of a dataset. One of the most common strategies to construct such a graph is based on selecting a fixed number k of nearest neighbours (kNN) for each point. However, the kNN heuristic may become inappropriate when sampling density or noise level varies across datasets. Strategies that try to get around this typically introduce additional parameters that need to be tuned. We propose a simple approach to construct an adaptive neighbourhood graph from a single parameter, based on quadratically regularised optimal transport. Our numerical experiments show that graphs constructed in this manner perform favourably in unsupervised and semi-supervised learning applications.


Introduction
The k-nearest neighbour (k-NN) graph is a staple of modern data science, enjoying near ubiquitous use as an initial step for extracting low-dimensional geometric or topological information from point cloud data. Given a collection of points {x i } N i=1 and a distance function d(x, y) ≥ 0, a graph is constructed in which each point x i is connected to the k points that are closest to it as measured by d. In many applications, data are generated by sampling from some underlying low-dimensional smooth manifold M ⊂ R d . In this scenario, the set of k-nearest neighbours of a point x i can be thought of as its local neighbourhood along the manifold, and the distances correspond to local distances along the manifold. Thus, neighbourhood graph construction comprises a common initial step that underpins most manifold learning methods, including Laplacian eigenmaps [1], diffusion maps [2], ISOMAP [3] and UMAP [4]. These tools are used in diverse areas of data science from image processing [5], to demography [6], to gene expression data [7]. (Right) Comparison of graphs (connections shown in black) constructed with quadratic OT and k-nearest neighbours, from samples (red) from a closed 1D curve embedded in 2D. Note that k-NN forms a graph that has more spurious edges in low-density regions.
In many applications, data are sampled non-uniformly from M and may be subject to noise. Nearest-neighbour graphs are known to struggle in these settings, since each point x i has a fixed neighbourhood size of k neighbours. In order to account for non-uniform sampling, the neighbourhood size k should instead be allowed to vary across the dataset -in regions of low sampling density a smaller choice of k may prevent spurious long-range connections. Various density-aware approaches have been developed to address this shortcoming, however these methods typically introduce additional parameters parameters that need tuning [8,7] or some form of density estimation [9].
In this paper, we present optimal transport (OT) as a flexible approach for constructing sparse local neighbourhood graphs that automatically adapts to variable sampling density and noise. Each point is assigned a unit of mass which may then be distributed to other points, forming the neighbourhood. This approach is related to recent work connecting entropy-regularised optimal transport to the Laplace-Beltrami operator in manifold learning: [10] establish that a bistochastic normalisation of the graph Laplacian can be constructed using entropy-regularized OT, and [11] show that this normalization improves robustness to heteroskedastic noise. However, entropic regularization has the drawback that the resulting kernel matrix is dense, representing a fully connected graph. We demonstrate that using instead a quadratic regularizer for the optimal transport problem [12] produces a sparse neighbourhood graph that can be interpreted similarly to a k-NN graph.
The rest of this paper is organized as follows. In Section 2 we motivate the use of optimal transport, and its quadratic regularization (quadratic OT) in particular, for constructing neighbourhood graphs. We investigate in Section 3 applications of this approach in several settings including manifold learning, semi-supervised learning, and denoising single-cell expression profiles.

Method
The approach we adopt follows [11] in defining a mass transport problem on the set {x i }. We assign to each point a unit of mass and allow for redistribution of mass to other points, subject to the constraint that the mass outgoing must equal the mass received. We prescribe the penalty incurred by sending a unit mass from x i to x j to be c ij = d(x i , x j ) p . Let π ij ≥ 0 denote the amount of mass sent from x i to x j : we will assume that π ij = π ji (symmetry) and that i π ij = 1 (mass balance). Thus, the matrix π can be treated as the adjacency matrix of a undirected weighted graph. The total penalty incurred by the transport plan π is i,j π ij c ij . Mass balance and symmetry assumptions can be written as the constraint set Π = {π : π ≥ 0, π = π , π1 = 1}. One may then consider the minimisation problem This problem admits a trivial solution π = I where each point sends its unit mass to itself (in which case the "neighbourhood" shrinks to a point), and thus reveals no information about local geometry. In order to encourage mass to be distributed away from a point and onto its "neighbours", we consider regularising the problem with a convex functional Ω and parameter ε > 0.
We refer the reader to [13,14] for an introduction to optimal transport with general convex regularisation. Note that something similar could also be achieved by modifying the cost function so that points cannot send mass to themselves (i.e. c ii = ∞). Previous work [11,10] considered the negative entropy Ω(π) = i,j π ij log π ij , which results in the symmetric entropy-regularised optimal transport problem. The nature of the entropy term ensures that the optimal π is strictly positive, i.e. one obtains a complete graph with local neighbourhood information encoded in the edge weights rather than presence of edges. Indeed, the optimal solution is characterised as a doubly-stochastic normalisation of the Gaussian kernel: π ij = e ui/ε e uj /ε e −cij /ε . Instead, we consider the quadratic regulariser Ω(π) = i,j π 2 ij , which yields a symmetric quadratically-regularised optimal transport problem [12]. As in the entropic case, the quadratic regulariser pushes the optimal π away from the identity. However, it is a well-known empirical result that the resulting matrices π that are sparse in the sense that most entries are identically zero. This can be seen from the characterisation of the optimal solution: π ij = 1 ε max{0, u i + u j − c ij }. Since the optimal π is sparse, the local neighbourhood of each point x i can be regarded as {x j : π ij > 0}. The matrix π therefore encodes a graph that can be interpreted similarly as k-NN graphs. The key difference here is that the parameter ε for the quadratic OT problem controls only the overall level of sparsity, without a fixed number of edges per node being prescribed as in the case of k-NN. [12] present a numerical method for solving the quadratically regularised problem (2), which can be specialised to the symmetric setting as we show in Algorithm 1 in the appendix.
The sparse nature of the quadratically regularized optimal transport plan can be intuited from the characterisation provided by [14], who show that the solution of (2) is equivalent to When Ω(π) = 1 2 π 2 , one has B Ω (π|ξ) = 1 2 π − ξ 2 . A straightforward computation shows that ξ = −c/ε, and so quadratic OT can be seen as finding a finding a Euclidean projection of −c/ε (a matrix with negative entries) onto Π (a set of transport plans constrained to be non-negative). Thus, one may expect that many entries of π will be identically zero.

Manifold learning with heteroskedastic noise
We investigate first the application of quadratic OT to the setting of a lowdimensional manifold embedded in a high-dimensional ambient space, with noise of varying intensity. This is the setting of 'heteroskedastic noise'where [11] proposed entropic OT. This setting occurs in single-cell RNA-sequencing data, where the noise level of each cell's expression profile is determined by the number of reads sequenced for that cell.
We generated the simulated dataset in Figure 2 with 10 3 points {x i } evenly spaced with respect to arc-length along the spiral parameterized by We then generated a random matrix R ∈ R 100×3 such that R R = I, i.e. the columns are orthogonal. Each point x i ∈ R 3 was then linearly embedded in the 100-dimensional ambient space with additive noise viax i = Rx i + η i . The noise model we use is where Z i follows the standard Gaussian distribution. In order to eliminate the effect of scaling on the effect of the spread parameter ε, we normalised all input data so that N −2 N i,j=1 x i −x j 2 = 1. Given the noisy high-dimensional dataset {x i }, we constructed local neighbourhood graphs using a variety of methods and a range of parameter values: • Quadratic OT: ε ∈ [10 −1.5 , 10 1.5 ] (logarithmic spacing) • Entropic OT [11]: ε ∈ [10 −2 , 10 1 ] (logarithmic spacing) • k-NN with Gaussian edge weights [1]: k-NN graph was constructed with k ∈ {5, 10, 15, 20, 25} and edge weights exp(− x i − x j 2 /ε), with ε in same range as entropic OT.
In Figure 3(a), we show the noisy data projected onto the first two principal components alongside representative graphs constructed using each method. This 2D embedding has a star-like shape with lobes that correspond to regions of large magnitudes of noise interspersed by regions of low noise which appear to be pinched inwards. Qualitative differences in the graphs are immediately visible -entropic OT forms edges between "lobes" of the spiral, corresponding to the dense nature of the graph. As visible in the zoomed region of high noise level, MAGIC (which relies on k-NN) forms few edges between points in the high noise region, with most edges branching away into the low-noise region. On the other hand, quadratic OT forms appreciably stronger connections between points in the high-noise region.
In Figures 3(b) and 4, we investigate the spectral properties of neighbourhood graphs are closely related to underlying geometry, and provide a convenient means of quantitative evaluation. To investigate this, we computed eigenmap embeddings from each of the constructed weighted graphs following [1]. Given a weighted adjacency matrix W computed by a method of choice, W was first random-walk normalized to form W : The spectral decomposition of W was then computed, yielding eigenvalues λ N ≤ λ N −1 ≤ . . . ≤ λ 1 = 1, and corresponding eigenvectors v N , . . . , v 1 . The top eigenvector is trivial and so the -dimensional eigenmap embedding for = 1, . . . , N − 1 is then the mapping As a reference, a graph was constructed from the clean 3-dimensional data {x i } (as opposed to {x i }) using a k-NN with k = 10 and a Gaussian kernel with ε = 0.025. From this, we computed reference eigenvectors {v ref i } i . To measure of how well the eigenvectors from the noisy high-dimensional data correspond to those computed from the clean low-dimensional data, we computed the mean principal angle [15] between the subspaces spanned by the top 10 non-trivial eigenvectors, i.e. span{v 2 , . . . , v 11 } and span{v ref 2 , . . . , v ref 11 }. We summarise these results in Figure 3(b). From this it is clear that the graph constructions obtained using entropic OT and a simple Gaussian kernel perform the worst overall -there are no parameter settings for which these methods perform well. Intuitively, this can be ascribed to the fact that the graphs are dense and try to model neighbourhoods using decaying edge weights. This effect can be seen also in Figure 3(a), where the graph constructed by entropic OT has visible edges between "lobes" of the spiral.
Methods that respect sparsity of the neighbourhood graph do significantly better, on the other hand. We note that k-NN and MAGIC perform comparably, with the key distinction that MAGIC has the capability of automatically adapting the effective bandwidth across regions of the dataset, while k-NN uses a universal bandwidth ε that needs to be tuned. We observe that quadratic OT outperforms both these methods. We believe that this is a result of the sparsity level and edge weights being not directly prescribed but instead controlled indirectly by the regularisation parameter ε.
Finally, we show in Figure 4 the best embeddings (in terms of the eigenspace error) in the top two non-trivial eigenvectors. As one would expect for eigenfunctions of a closed 1D curve with uniform density [9], the reference embedding is a perfect circle. Both entropic OT and the simple Gaussian kernel produce embeddings that are still star-shaped, making no improvement beyond the PCA embeddings displayed in 3(a). k-NN with Gaussian weights and MAGIC produce circular embeddings, however points are irregularly distributed, concentrating at several points. Only quadratic OT produces a near-perfect circle.

Application to semi-supervised learning
Next we turn to investigating an application of our graph construction method in which the downstream steps are more involved, namely semi-supervised learning (SSL). The problem setting here is that we are provided with data {x i } N i=1 , of which a small subset L ⊂ {1, . . . , N } are annotated with labels { i } i∈L . The semi-supervised learning problem is then to infer a set of labels {ˆ i } that (approximately) preserves the known labelsˆ i ≈ i , i ∈ L and also estimates the unknown labels {ˆ i } i / ∈L . These inference procedures typically operate under the heuristic that similarly labelled points should be in close proximity, necessitating the construction of a neighbourhood graph.
While there are many schemes for solving the SSL problem, we will focus on the classical method Learning Local and Global Consistency (LLGC) of [16], for which we now provide a brief description. LLGC aims to optimise label assignments to balance between matching the known labels (fitting) and consistency of inferred labels with the data geometry (consistency).
Let P be a matrix of dimensions N × C, with P ij = 1{x i has label j}. Thus, the rows of P are one-hot encodings for labelling known points, and unlabelled points correspond to all-zeros. The aim is to infer a matrix Q ∈ R N ×C ≥0 describing the labels of all points (both labelled and unlabelled), such that For a neighbourhood graph with row-stochastic weighted adjacency matrix W ij and regularisation parameter µ > 0, the LLGC optimisation problem is written as (3) We generated data from a multi-armed spiral similar to that used in [17]. We chose 10 arms started at angles θ = π/5, 2π/5, . . . , 2π, and for each initial angle θ the spiral arm was parameterised as We set t 0 = 1, t 1 = 5 and sampled 150 points along each arm to be uniformly spaced with respect to the arc-length of the branch. Following the approach of Section 3.1, the clean data x i were then embedded into 100 dimensions and noise was added For each branch, 2.5% of points were selected at random to be labeled. We then solved (3) with W taken to be the adjacency matrix constructed either by quadratic OT, or using a k-NN graph with (row-normalized) Gaussian weights. Each method was re-run with different parameter settings: ε ∈ [10 −2 , 10 1 ] (logarithmically spaced) for quadratic OT, and ε ∈ [10 −3 , 10 0 ] (logarithmically spaced), k ∈ [1, 50] (linearly spaced).
Results are summarised in Figure 5: from this it is apparent that in the presence of high-dimensional noise, quadratic OT significantly outperforms k-NN with Gaussian weights. In Figure 5(b), we show both inferred labels and neighbourhood graphs for the best results achieved by each method. Quadratic OT appears to have mostly correctly assigned labels and the neighbourhood graph is well-adapted to the underlying geometry, with very few edges between spiral arms. On the other hand, k-NN with Gaussian weights produces many incorrect labels. We note that this coincides with the presence of many edges connecting adjacent spiral arms.
In Figure 6, we explore the performance of both quadratic OT and kNN as the number of points along the spiral arms increases, and therefore the density increases. Here we generated the 10-armed spiral data following the same parametrization (4) with randomly sampled points along t ∈ [0.5, 5.0]. We labeled 1% of the data. The heat map summarizes accuracy for varying parameter ranges as the number of points per arm increases. From this it is clear that quadratic OT achieves comparable performance across problem scales for a wide range of ε, while the optimal value of k for kNN changes as we vary the number of points per arm. This can be attributed to the fact that quadratic OT balances the mass going into as well as the mass coming out of each point, instead of simply connecting each point to a fixed number k of neighbours.

Application to Single Cell RNA-sequencing
Single cell measurement technologies are revolutionizing the biological sciences by providing high-dimensional measurements of cell state. By sequencing the RNA molecules within individual cells, single-cell RNA-sequencing (scRNA-seq) quantifies the gene expression profiles of these cells (i.e. which genes are 'off' and which are 'on' and at what level). Recent droplet-based technologies have tremendously increased the throughput of scRNA-seq, so that large numbers of cells can be profiled, and have brought 'Big Data' to biology.
Sequencing more cells comes at the cost of sequencing fewer molecules of RNA per cell, and so scRNA-seq data is notoriously noisy. Lowly expressed genes may not be detected (i.e. they are 'dropped out'), which makes the analysis of gene-gene relationships difficult. This has motivated several groups to develop tools for denoising single-cell expression profiles. One of the most popular tools for this is MAGIC (Markov affinity-based graph imputation of cells) developed by [7].
MAGIC constructs a local neighbourhood graph from an adaptive Gaussian kernel to impute the data. For each point x i in the data, the Gaussian parameter σ i is chosen by σ i = d(x i , x j i ) where x j i is the j-th nearest neighbour of x i . This allows the kernel to have larger bandwidth in sparse areas and smaller bandwidth in dense areas, and was shown to be crucial in the original paper [7]. MAGIC then denoises the data by smoothing it according to a diffusion along this graph.
We tested whether quadratic OT could produce a neighbourhood graph without explicitly adjusting for density. Figure 7 shows the results of denoising a single cell dataset from [7] for neighbourhood graphs constructed by quadratic OT, entropic OT, and the original formulation of MAGIC with the adaptive Gaussian kernel. We find that quadratic OT performs nearly identical to MAGIC, which confirms that quadratic OT is adaptive to variable density. On the other hand, entropic OT fails to smooth points near the centre and loses geometric information that is present on the original.

Discussion
Quadratically regularised optimal transport provides a simple and robust method for constructing weighted neighbourhood graphs with only a single parameter ε.
Our numerical experiments provide evidence that, in many geometric machine learning applications, this graph construction can replace the role of k-nearest neighbour graphs to achieve similar or superior performance. This can be attributed to the fact that, instead of connecting each point to a fixed number k of neighbours, OT balances the mass coming into and going out of each point. The constructed graph along with its weight is adaptive to its domain which is a major advantage over classical k nearest neighbour graph based method especially on a high dimensional dataset. In future work, it will be of interest to study theoretical properties of the linear operator derived from the quadratically regularised optimal transport problem. Analysis of the behavior of this operator in the limit of infinite data may plausibly arrive at a convergence result to the Laplace-Beltrami operator, a result that is established in the literature for many other common graph constructions [18]. We hypothesise that such theoretical properties of the graph operator may involve a connection between quadratically regularised optimal transport and subdiffusions [19].

Quadratic OT
Entropic OT kNN+Gauss, k=10 Gaussian MAGIC, k=10 True  Note that for kNN the the optimal value of k depends on the number of points per arm, but for quadratic OT a single value (ε = 1) is optimal.