Analysis and algorithms for $\ell_p$-based semi-supervised learning on graphs

TL;DR

This paper investigates $ abla_p$-based semi-supervised learning on graphs, providing convergence theory, developing scalable algorithms, and demonstrating the effectiveness of Lipschitz learning ($p= fty$) on real datasets with few labels.

Contribution

It introduces new convergence results for $p$-Laplacian on $k$-NN graphs, develops efficient algorithms for $p$-Laplace equations, and empirically validates the advantages of Lipschitz learning in semi-supervised tasks.

Findings

$p$-Laplacian retains data distribution information on $k$-NN graphs as $p o fty$.

Lipschitz learning ($p= fty$) performs well with few labels on $k$-NN graphs.

Algorithms demonstrate convergence and scalability on synthetic and real datasets.

Abstract

This paper addresses theory and applications of $\ell_p$-based Laplacian regularization in semi-supervised learning. The graph $p$-Laplacian for $p>2$ has been proposed recently as a replacement for the standard ($p=2$) graph Laplacian in semi-supervised learning problems with very few labels, where Laplacian learning is degenerate. In the first part of the paper we prove new discrete to continuum convergence results for $p$-Laplace problems on $k$-nearest neighbor ($k$-NN) graphs, which are more commonly used in practice than random geometric graphs. Our analysis shows that, on $k$-NN graphs, the $p$-Laplacian retains information about the data distribution as $p\to \infty$ and Lipschitz learning ($p=\infty$) is sensitive to the data distribution. This situation can be contrasted with random geometric graphs, where the $p$-Laplacian forgets the data distribution as $p\to \infty$. We also present a general framework for proving discrete to continuum convergence results in graph-based learning that only requires pointwise consistency and monotonicity. In the second part of the paper, we develop fast algorithms for solving the variational and game-theoretic $p$-Laplace equations on weighted graphs for $p>2$. We present several efficient and scalable algorithms for both formulations, and present numerical results on synthetic data indicating their convergence properties. Finally, we conduct extensive numerical experiments on the MNIST, FashionMNIST and EMNIST datasets that illustrate the effectiveness of the $p$-Laplacian formulation for semi-supervised learning with few labels. In particular, we find that Lipschitz learning ($p=\infty$) performs well with very few labels on $k$-NN graphs, which experimentally validates our theoretical findings that Lipschitz learning retains information about the data distribution (the unlabeled data) on $k$-NN graphs.

Figures

Tables

Table 1. Table 1. Newton’s method with homotopy solved problem S for p = 50 , ϵ < 10 − 12 formulae-sequence 𝑝 50 italic-ϵ superscript 10 12 p=50,\epsilon<10^{-12} in just 56 56 56 iterations. We take increasingly large steps in p 𝑝 p while remaining in the quadratic convergence regime, suggesting we could go much farther than p = 50 𝑝 50 p=50 .

	$p$
$N$	$3$	$4$	$6$	$8$	$10$	$15$	$20$	$25$	$30$	$40$	$50$
1	2e-01	3e-02	1e-02	8e-04	2e-04	8e-05	4e-06	3e-07	3e-08	1e-09	2e-11
2	5e-02	3e-03	4e-04	1e-05	7e-07	1e-06	2e-08	6e-10	3e-11	5e-12	3e-14
3	7e-03	2e-04	4e-05	7e-08	2e-10	3e-10	3e-13	4e-15	1e-16	9e-17	-
4	8e-04	1e-05	4e-06	2e-09	1e-13	2e-14	-	-	-	-	-
5	5e-05	4e-07	8e-07	5e-12	-	-	-	-	-	-	-
6	4e-07	1e-09	7e-08	2e-16	-	-	-	-	-	-	-
7	2e-11	2e-14	6e-10	-	-	-	-	-	-	-	-
8	8e-15	-	1e-12	-	-	-	-	-	-	-	-
9	-	-	2e-16	-	-	-	-	-	-	-	-

Table 2. Table 2. The Newton-like method with homotopy solved problem S for p = 50 𝑝 50 p=50 and ϵ < 10 − 12 italic-ϵ superscript 10 12 \epsilon<10^{-12} in 52 52 52 iterations. Similar to Newton’s method, increasingly large step sizes in p 𝑝 p are possible, suggesting the ability to go even higher in p 𝑝 p .

	$p$
$N$	$3$	$4$	$6$	$8$	$10$	$15$	$20$	$25$	$30$	$40$	$50$
1	4e-02	1e-02	9e-03	4e-03	2e-03	3e-03	1e-03	7e-04	5e-04	6e-04	3e-04
2	4e-03	3e-03	3e-03	2e-03	1e-03	3e-03	2e-03	7e-04	3e-04	9e-04	6e-04
3	1e-03	8e-04	2e-03	6e-04	5e-04	2e-03	9e-04	1e-04	9e-05	1e-04	3e-05
4	1e-04	1e-04	5e-04	5e-04	6e-05	3e-04	5e-04	3e-16	3e-16	3e-16	3e-16
5	6e-06	1e-05	9e-05	2e-04	2e-07	6e-05	3e-16	-	-	-	-
6	2e-16	2e-16	2e-16	8e-06	4e-16	3e-16	-	-	-	-	-
7	-	-	-	3e-16	-	-	-	-	-	-	-

Equations

$$J_2(u):=\frac{1}{4}\sum _{x,y\in \mathcal{X}}{w}_{xy}(u(x)-u(y){)}^2.$$

$${\Delta }_2^{G}u(x):=\sum _{y\in \mathcal{X}}{w}_{xy}(u(y)-u(x)).$$

$$J_{p}(u):=\frac{1}{2p}\sum_{x,y\in\mathcal{X}}w_{xy}\big{|}u(x)-u(y)\big{|}^{p}.$$

$${\Delta }_p^{G}u(x):=\sum _{y\in \mathcal{X}}{w}_{xy}|u(x)-u(y){|}^{p-2}(u(y)-u(x)).$$

$$|u(x)-u(y)|\le C{\left({\int }_{\Omega }|\nabla u{|}^{p}dx\right)}^{1/p}|x-y{|}^{1-d/p}$$

$$J_{\infty }(u)=\underset{x,y\in \mathcal{X}}{\max }{w}_{xy}|u(x)-u(y)|.$$

$${\Delta }_{\infty }^{G}u(x):=\underset{y\in \mathcal{X}}{\min }{w}_{xy}(u(y)-u(x))+\underset{y\in \mathcal{X}}{\max }{w}_{xy}(u(y)-u(x)).$$

$$\sum _{u(y)>u(x)}{w}_{xy}(u(y)-u(x){)}^{p-1}=\sum _{u(x)>u(y)}{w}_{xy}(u(x)-u(y){)}^{p-1}.$$

$$\underset{y\sim x}{\max }(u(y)-u(x))=\underset{y\sim x}{\max }(u(x)-u(y)),$$

$$\underset{y\sim x}{\max }(u(y)-u(x))+\underset{y\sim x}{\min }(u(y)-u(x))=0,$$

$${\mathcal{L}}_p^{G}u(x)=\frac{1}{d_{x}p}{\Delta }_2^{G}u(x)+\lambda \left(1-\frac{2}{p}\right){\Delta }_{\infty }^{G}u(x),$$

$${\Delta }_{p}u:=\text{div}(|\nabla u{|}^{p-2}\nabla u)=0.$$

$$\underset{u}{\min }{\int }_{\Omega }|\nabla u{|}^{p}dx.$$

$${\Delta }_{p}u=|\nabla u{|}^{p-2}(\Delta u+(p-2){\Delta }_{\infty }u),$$

$${\Delta }_{\infty }u:=\frac{1}{|\nabla u{|}^2}\sum _{i,j=1}^d{u}_{x_i{x}_j}{u}_{x_i}{u}_{x_j}.$$

$$\frac{1}{p}\Delta u+\left(1-\frac{2}{p}\right){\Delta }_{\infty }u=0.$$

$$\underset{u:\mathcal{X}\to \mathbb{R}}{\min }{J}_p(u)\text{ subject to }u(x)=g(x)\text{ for all }x\in \mathcal{O},$$

$$J_{p}(u):=\frac{1}{2p}\sum_{x,y\in\mathcal{X}}w_{xy}\big{|}u(x)-u(y)\big{|}^{p}+\sum_{x\in\mathcal{X}}f(x)u(x),$$

$$\{\begin{array}{rlrl}-{\Delta }_p^{G}u(x)& =f(x)& & \text{if }x\in \mathcal{X}\setminus \mathcal{O}\\ u(x)& =g(x)& & \text{if }x\in \mathcal{O},\end{array}$$

$$\{\begin{array}{rlrl}-{\mathcal{L}}_p^{G}u(x)& =f(x)& & \text{if }x\in \mathcal{X}\setminus \mathcal{O}\\ u(x)& =g(x)& & \text{if }x\in \mathcal{O},\end{array}$$

$$\beta \le \rho (x)\le {\beta }^{-1}\text{ for all }x\in \Omega .$$

$${\mathcal{X}}_n:=\{X_1,X_2,\dots ,X_n\}.$$

$${\int }_{B(0,1)}\eta (z)dz=1.$$

$$r\eta (r)+\theta (r-r_0{)}^2\le r_0\eta (r_0)\text{ for all }0\le r\le 1.$$

$$\eta (r)=\{\begin{array}{ll}2,& \text{if }0\le r\le \frac{1}{2}\\ \frac{1}{r},& \text{if }\frac{1}{2}\le r\le \frac{3}{4}\\ 0,& \text{if }r\ge \frac{3}{4}.\end{array}$$

$${\mathcal{L}}^{n,\epsilon }u(x)=\frac{2}{{\sigma }_{\eta }n{\epsilon }^{d+2}}\sum _{y\in {\mathcal{X}}_n}{\eta }_{\epsilon }(|x-y|)(u(y)-u(x)).$$

$$d^{n,\epsilon }(x)=\sum _{y\in {\mathcal{X}}_n}{\eta }_{\epsilon }(|x-y|).$$

$${\mathcal{L}}_{rw}^{n,\epsilon }u(x)=\frac{2}{{\sigma }_{\eta }{\epsilon }^2{d}^{n,\epsilon }(x)}\sum _{y\in {\mathcal{X}}_n}{\eta }_{\epsilon }(|x-y|)(u(y)-u(x)).$$

$${\mathcal{L}}_{\infty }^{n,\epsilon }u(x)=\frac{1}{r_0^2\eta (r_0){\epsilon }^2}\left(\underset{y\in {\mathcal{X}}_n}{\min }\left\{{\eta }_{\epsilon }(|x-y|)(u(y)-u(x))\right\}+\underset{y\in {\mathcal{X}}_n}{\max }\left\{{\eta }_{\epsilon }(|x-y|)(u(y)-u(x))\right\}\right),$$

$${\mathcal{L}}_p^{n,\epsilon }u(x)=\frac{1}{p}{\mathcal{L}}_{rw}^{n,\epsilon }u(x)+\left(1-\frac{2}{p}\right){\mathcal{L}}_{\infty }^{n,\epsilon }u(x).$$

Full text

In the first part of the paper we prove new discrete to continuum convergence results for $p$-Laplace problems on $k$-nearest neighbor ($k$-NN) graphs, which are more commonly used in practice than random geometric graphs. Our analysis shows that, on $k$-NN graphs, the $p$-Laplacian retains information about the data distribution as $p\to\infty$ and Lipschitz learning ($p=\infty$) is sensitive to the data distribution. This situation can be contrasted with random geometric graphs, where the $p$-Laplacian forgets the data distribution as $p\to\infty$. We also present a general framework for proving discrete to continuum convergence results in graph-based learning that only requires pointwise consistency and monotonicity.

In the second part of the paper, we develop fast algorithms for solving the variational and game-theoretic $p$-Laplace equations on weighted graphs for $p>2$. We present several efficient and scalable algorithms for both formulations, and present numerical results on synthetic data indicating their convergence properties. Finally, we conduct extensive numerical experiments on the MNIST, FashionMNIST and EMNIST datasets that illustrate the effectiveness of the $p$-Laplacian formulation for semi-supervised learning with few labels. In particular, we find that Lipschitz learning ($p=\infty$) performs well with very few labels on $k$-NN graphs, which experimentally validates our theoretical findings that Lipschitz learning retains information about the data distribution (the unlabeled data) on $k$-NN graphs.

Source Code: https://github.com/mauriciofloresML/Laplacian_Lp_Graph_SSL.git

Funding: The authors gratefully acknowledge National Science Foundation grants 1713691, 1821266, 1830418, and a University of Minnesota Grant in Aid Award.

1. Introduction

Data science problems, such as regression and classification, are pervasive in today’s world, and the size of datasets is growing rapidly. In the supervised setting, data needs to be labeled, requiring substantial effort (e.g. writing a transcript for speech recognition), or may require expert input (deciding whether a brain scan is healthy or not). In contrast, unlabeled data can often be acquired in large quantities with substantially less effort. Semi-supervised learning harnesses the additional information present in unlabeled data to improve learning tasks. This can include geometric or topological properties of unlabeled data, which can provide valuable information about where to place decisions boundaries, for instance. This can be contrasted with fully supervised algorithms, which only make use of labeled data. Fully supervised learning algorithms typically learn parameterized functions and require an abundant amount of labeled data.

A common setting within semi-supervised learning is graph-based semi-supervised learning, which is concerned with propagating label information on graphs. Here, we are given an undirected weighted graph $G=(\mathcal{X},\mathcal{W})$, where $\mathcal{X}$ are the vertices and $\mathcal{W}=\{w_{xy}\}_{x,y\in\mathcal{X}}$ are nonnegative edge weights, which are chosen so that $w_{xy}\approx 1$ when $x$ is similar to $y$, and $w_{xy}\approx 0$ when $x$ and $y$ are dissimilar. Each vertex $x$ in the observation set $\mathcal{O}\subset\mathcal{X}$ is assigned a label $g(x)$, where $g:\mathcal{O}\rightarrow\mathbb{R}^{k}$. In a classification problem with $k$ classes, the $i^{\rm th}$ class is usually assigned the label vector $g(x)=e_{i}$, where $e_{i}$ is the $i^{\rm th}$ standard basis vector in $\mathbb{R}^{k}$, that is, the vector with all zeros and a one in the $i^{\rm th}$ coordinate (called a “one-hot” vector in machine learning). The task of graph-based semi-supervised learning is to extend the labels from the observation set $\mathcal{O}$ to a label function $u:\mathcal{X}\to\mathbb{R}^{k}$ on the whole graph in some meaningful way. In practice, the equations that are solved for propagating labels are separable among the coordinates of $\mathbb{R}^{k}$, and the problem reduces to solving for $k$ functions $u_{i}:\mathcal{X}\to\mathbb{R}$, one for each class $i=1,\dots,k$, and setting $u=(u_{1},\dots,u_{k})$ (this is referred to as the “one-vs-rest” approach in machine learning). Thus, we can, without loss of generality, focus on the scalar case $k=1$ for algorithms and analysis in this paper.

Since the problem of extending labels is a priori an ill-posed problem (there are infinitely many solutions), one usually makes the semi-supervised smoothness assumption, which asks that the learned labeling function $u:\mathcal{X}\to\mathbb{R}$ should be smooth in dense regions of the graph [18]. The smoothness assumption is often enforced by defining a functional (or regularizer) $J(u)$ that measures the smoothness of a labeling $u:\mathcal{X}\to\mathbb{R}$, and then minimizing $J(u)$ subject to either hard label constraints $u(x)=g(x)$ for $x\in\mathcal{O}$, or a soft penalty constraint like the mean squared error in the labels. Soft penalties are useful when the labels are corrupted by noise. In this paper, we are concerned with learning problems with very few labels, in the range of one label per class, so a basic assumption we make is that the labels are clean and are not corrupted by noise. Thus, the hard constraint is natural to impose and nothing is gained by considering a soft constraint. In fact, the soft-constraint is more likely to be ignored, when there are very few labels, unless the penalty parameter is chosen sufficiently large so that the constraint is essentially a hard one. The soft-constrained problem would be interesting to consider in the context of more moderate label rates, where noisy labels can be better tolerated. All of the techniques we discuss in this paper extend directly, with minor modifications, to problems with soft constraints.

One of the most widely used methods in semi-supervised learning is Laplacian regularization [62], which uses the smoothness functional

[TABLE]

Minimizing $J_{2}$ attempts to force similar data points in dense regions of the graph to have similar labels. Minimizers of $J_{2}$ are graph harmonic, and solve the graph $2$-Laplace equation $\Delta^{G}_{2}u=0$, where

[TABLE]

In classification, the values of $u$ are rounded to the nearest label. Laplacian regularization, and ideas based upon it, are very widely used in machine learning [3, 58, 59, 61, 31, 32, 54, 56, 57], and have achieved great successes. However, it has been noted first in [43] and later in [21], that Laplacian regularization becomes ill-posed (degenerate) in problems with very few labels. We say a graph-based learning problem is ill-posed in the limit of infinite unlabeled data and finite labeled data if the sequence of learned functions does not continuously attain the labeled (e.g., boundary) data in the continuum limit. In this case, the learned function becomes nearly constant on the whole graph, with sharp spikes near the labeled data. Thus, even with a hard constraint the labels are almost entirely ignored. See Figure 1a for a depiction of this degeneracy. In the continuum, this is merely reflecting the fact that the capacity of a point is zero in dimension $d\geq 2$ [38].

To address this issue, El Alaoui et al. [21] proposed a class of $\ell_{p}$-based Laplacian regularizers, which use the smoothness functional

[TABLE]

Choosing $p>2$ places a heavier penalty on large gradients $|u(x)-u(y)|$, which discourages the solution from developing sharp spikes. We also note that choosing $1\leq p<2$, often $p=1$, encourages the gradient $|u(x)-u(y)|$ to be sparse, and can be viewed as a relaxation of graph-cut energies. This can yield good results for classification at moderate label rates (see, e.g., [35]), but at very low label rates the issue with spikes is even more pronounced when $p<2$, and the results are similar to Laplace learning.

We note that minimizers of $J_{p}$ satisfy the graph $p$-Laplace equation $\Delta_{p}^{G}u=0$, where

[TABLE]

We call (1.4) the variational graph $p$-Laplacian. Figure 1 depicts $\ell_{p}$-regularization for different values of $p$. As $p$ increases the learned function transitions more smoothly between labeled and unlabeled points. From a continuum perspective, the energy $J_{p}$ is related to the $p$-Dirichlet energy $\int_{\Omega}|\nabla u|^{p}\,dx$, and the Sobolev embedding $W^{1,p}(\Omega)\hookrightarrow C^{0,1-d/p}(\Omega)$ allows isolated boundary points when $p>d$, where $d$ is the dimension.111Here, $\Omega\subset\mathbb{R}^{d}$ is an open, bounded domain, and $W^{1,p}(\Omega)$ is the Sobolev space of functions $u:\Omega\to\mathbb{R}$ such that $\int_{\Omega}u^{p}+|\nabla u|^{p}\,dx<\infty$. Indeed, by Morrey’s inequality [25] we have

[TABLE]

whenever $p>d$ and $|x-y|\leq\tfrac{1}{2}\text{dist}(x,\partial\Omega)$. Morrey’s inequality implies that $u$ is Hölder continuous, and prevents spikes in Figure 1a from occurring. In particular, the continuum $p$-Dirichlet problem (see (1.10)) with constraints at isolated points is well-posed (e.g., admits a unique solution attaining the boundary data continuously) if and only if $p>d$.

The variational graph $p$-Laplacian (1.4) has appeared previously in machine learning [2, 6, 60], but was first suggested for problems with few labels in [21] with $p\geq d+1$. Recently, it was rigorously proven that $\ell_{p}$-based regularization is ill-posed (its minimizer is degenerate) for $p\leq d$, and well-posed for $p>d$ in the continuum limit of infinite unlabeled and finite labeled data [50]. This justifies the continuum heuristics described above.

Formally sending $p\to\infty$ in (1.3) one obtains Lipschitz learning [36, 41], which corresponds to the smoothness functional

[TABLE]

We note that minimizers of (1.6) are not unique. To see why this is the case, note that if $w_{xy}|u(x)-u(y)|<J_{\infty}(u)$, then we can change the values of $u(x)$ and $u(y)$ slightly, without changing the value of $J_{\infty}(u)$. So minimizers cannot be unique unless the maximum gradient is attained everywhere. Among the non-unique minimizers, one generally looks for one whose gradient cannot be locally improved (i.e., made smaller by adjusting some values of $u(x)$). More precisely, in [36] the authors show that there is a unique minimizer whose gradient is smallest in the lexicographical order, called the lex-minimizer. This turns out to be equivalent to the notion of absolutely minimal, which has been used in the partial differential equation (PDE) and analysis community to select the unique Lipschitz extension in the continuum [4].

Lex-minimizers of (1.6) satisfy the graph $\infty$-Laplace equation $\Delta^{G}_{\infty}u=0$ where

[TABLE]

To see why (1.7) is the correct form for the graph $\infty$-Laplacian, consider a graph $p$-harmonic function $u$, which satisfies $\Delta^{G}_{p}u(x)=0$, where $\Delta^{G}_{p}$ is defined in (1.4). We split the terms in the sum defining $\Delta^{G}_{p}u(x)$ by their sign to obtain

[TABLE]

Taking the $p^{\rm th}$ root of both sides and sending $p\to\infty$ yields

[TABLE]

where we write $y\sim x$ if $w_{xy}>0$. This can be simplified to

[TABLE]

which is the graph $\infty$-Laplacian defined in (1.7) for an unweighted graph. To obtain the weighted graph $\infty$-Laplacian (1.7) in the limit as $p\to\infty$, we simply replace $w_{xy}$ in (1.4) with $w_{xy}^{p}$.

It was proven in [10] that Lipschitz learning is well-posed with arbitrarily few labels, and the continuum limit on random geometric graphs is the continuum $\infty$-Laplace equation (see (1.11)). Several papers (see, e.g., [21, 10]) have noted that the continuum $\infty$-Laplacian (1.11) does not involve the data distribution, making it presumably unsuitable for semi-supervised learning, which is supposed to use properties of the unlabeled data, often through its distribution. In [10], it was shown how to re-weight the graph to introduce varying degrees of sensitivity to the data distribution in Lipschitz learning.

In order to combine the well-posedness of Lipschitz learning with the distributional sensitivity of Laplacian regularization, it is natural to augment the $2$-Laplacian with a small $\infty$-Laplace term and solve an equation of the form $\Delta^{G}_{2}u+\varepsilon\Delta^{G}_{\infty}u=0$. To this end, we define the game theoretic $p$-Laplacian on the graph (the name will be explained shortly)

[TABLE]

where $d_{x}=\sum_{y\in\mathcal{X}}w_{xy}$ is the degree of vertex $x$, and $\lambda>0$ is a constant. For semi-supervised learning with the game-theoretic $p$-Laplacian we solve $\mathscr{L}^{G}_{p}u=0$ subject to $u=g$ on $\mathcal{O}$ (see Section 1.1 for precise definitions). The second author proved in [8] that the game theoretic $p$-Laplacian is well-posed with very few labels for $p>d$, and argued for the use of this formulation as an alternative regularization for semi-supervised learning on graphs. In the context of these results, $\lambda$ is chosen (explicitly) depending on the kernel used to define the weights $w_{xy}$, in order to produce the consistency results described below. Since having two parameters, $p$ and $\lambda$, is redundant in practice, we take $\lambda=1$ in the numerical sections of the paper. Compared to the variational $p$-Laplacian, the game-theoretic $p$-Laplacian appears better conditioned numerically when $p$ is large, since it does not require computing large powers of $p$. Another main difference is that the game-theoretic $p$-Laplacian does not arise through an optimization problem, and so the methods for solving the equation are somewhat different. We note that the game-theoretic graph $p$-Laplacian (and similar models) have been used very recently for data clustering and learning problems [24, 22, 23, 30], though not in the context of very few labeled data points. A related definition of the game-theoretic $p$-Laplacian on graphs was also studied in [42]. We also mention recent work [49, 15] that approaches the semi-supervised learning problem with few labels by re-weighting the graph so that the weights $w_{xy}$ are large near labels.

Both the variational (1.4) and the game-theoretic (1.8) graph $p$-Laplace equations are consistent in the continuum, on random geometric graphs, with the $p$-Laplace equation

[TABLE]

It is important to point out that the weights $1/p$ and $\lambda(1-2/p)$ in (1.8) are chosen precisely so that $\mathscr{L}^{G}_{p}$ is consistent with $\Delta_{p}$. The operator $\Delta_{p}$ is called the $p$-Laplacian, and solutions of (1.9) are called $p$-harmonic functions [40]. The $p$-Laplace equation arises as the necessary conditions (Euler-Lagrange equation) for the $p$-Dirichlet problem

[TABLE]

Note that we can expand the $p$-Laplacian to obtain

[TABLE]

where $\Delta_{\infty}$ is the $\infty$-Laplacian, given by

[TABLE]

Thus, any solution of $\Delta_{p}u=0$ also satisfies

[TABLE]

The left hand side in (1.12) is often called the game-theoretic or homogeneous $p$-Laplacian since it arises in two player stochastic tug-of-war games [47, 39]. This justifies the definition (1.8) of the game-theoretic graph $p$-Laplacian. We note that while the $p$-Laplace equation (1.9) is equivalent to the game theoretic version (1.12) at the continuum level, these are different formulations at the discrete level.

Given the recent interest in graph $p$-Laplacian models in machine learning, it is important to have both strong theoretical results that are relevant in practice, and efficient and scalable algorithms for solving the equations in real-world settings. Most of the literature on discrete to continuum convergence in graph-based learning, such as the recent work on the $p$-Laplacian [50, 8], assumes the graph is a random geometric graph. However, such graphs have poor sparsity properties222As an example of the poor sparsity properties of random geometric graphs, consider the MNIST dataset, which has $70,000$ images of handwritten digits. A $k$-nearest neighbor graph on MNIST with $k=3$ neighbors is connected, while the sparsest random geometric graph that is connected requires a bandwidth of $\varepsilon=8.5$ and has on average $294$ neighbors per image, which is roughly $100$ times less sparse than the $k$-nearest neighbor graph. Graph connectivity is one of the basic properties required by most graph-based learning algorithms., and practitioners almost always use some form of a $k$-nearest neighbor ($k$-NN) graph instead. Thus, it is important to develop theory for the graph $p$-Laplacian on $k$-NN graphs in order to study the graphs that are used in practice. It turns out, looking forward to Section 2, that the discrete to continuum theory for the $p$-Laplacian on $k$-NN graphs is fundamentally different from the existing theory on random geometric graphs, and this discrepancy accounts for the results of our numerical experiments conducted later in the paper.

From a computational perspective, there are relatively few works on fast algorithms for graph $p$-Laplacians. Kyng et al., [36] developed an efficient algorithm for Lipschitz learning ($p=\infty$). Their algorithm has a poor worst case complexity analysis (roughly quadratic complexity in the number of data points), but seems to run very fast in practice. Oberman [45] considers the game-theoretic formulation on regular grids in dimensions $d=2,3$, and developed a fast semi-implicit solution method, as well as gradient-descent methods. One contribution of this paper is an adaptation of Oberman’s semi-implicit method to the graph setting. Other works [24, 22, 23] use slow iterative methods, such as Jacobi iteration or gradient descent. It was suggested in [21] to use Newton’s method for the variational $p$-Laplacian, but the method was not investigated in any depth. The energy $J_{p}$ is smooth and convex, but not strongly convex when $p>2$. Other works, such as [36], suggest to use convex programming to solve the variational $p$-Laplacian.

1.1. Main results and contributions

Before summarizing the main contributions of our paper, let us give precise definitions of the problems we study. We let $G=(\mathcal{X},\mathcal{W})$ be a connected graph with vertices $\mathcal{X}$ and edge weights $\mathcal{W}=\{w_{xy}\}_{x,y\in\mathcal{X}}$. The subset of labeled vertices is denoted $\mathcal{O}\subset\mathcal{X}$, and the label function is $g:\mathcal{O}\to\mathbb{R}$. We denote the number of data points in $\mathcal{O}$ by $m$, and the number of data points in $\mathcal{X}$ by $n$. The $\ell_{p}$-based Laplacian regularized learning problem [21] is given by

[TABLE]

where we extend the definition of $J_{p}$ to be

[TABLE]

where $f$ is a source function. The unique minimizer $u:\mathcal{X}\rightarrow\mathbb{R}$ of (1.13) satisfies the optimality conditions

[TABLE]

where $\Delta^{G}_{p}$ is defined in (1.4). In machine learning applications we always take $f\equiv 0$ in (1.14) and (1.15). However, it is useful to formulate the equation in more generality so that we can construct exact solutions by choosing $u(x)$ and computing $f(x)=-\Delta^{G}_{p}u(x)$ accordingly. This allows us to evaluate and compare the convergence rates of different methods. We will refer to (1.13) as the variational problem, and to (1.15) as the variational $p$-Laplace equation.

The game-theoretic graph $p$-Laplacian equation is given by

[TABLE]

where $\mathscr{L}^{G}_{p}$ is the game-theoretic graph $p$-Laplacian defined in (1.8). We will refer to (1.16) as the game-theoretic problem. As before, we are only concerned with $f\equiv 0$ in machine learning, but it is convenient to proceed in generality. The Lipschitz learning problem corresponds to the game-theoretic problem with $p=\infty$.

We now summarize the contributions of the paper.

1.1.1. Discrete to continuum on $k$NN graphs

In Section 2, we give a detailed analysis of discrete to continuum convergence of graph-based algorithms on $k$-nearest neighbor graphs. We prove that the continuum limit of the game-theoretic $p$-Laplacian on symmetrized $k$-NN graphs has a significantly different form compared to random geometric graphs. In particular, we show that the continuum operator has an additional drift term along the gradient of the data distribution that does not vanish as $p\to\infty$. The drift term arises through the symmetrization process in the $k$-NN graph construction. This means the conventional wisdom that Lipschitz learning is not sensitive to the data distribution (see, e.g., [21, 10, 50]) is a phenomenon specific to random geometric graphs, and does not hold true for $k$-NN graphs, which highlights the need to perform analysis on different graph constructions. This result is also borne out in our experimental results in Section 4, which show that the $p=\infty$ graph Laplacian gives good results for semi-supervised learning with very few labels, indicating sensitivity to the data distribution as predicted by our theoretical results.

We also prove a discrete to continuum convergence result for general elliptic learning algorithms on graphs. The result is very general and only requires the equation to have a certain monotonicity property (essentially ellipticity) that allows for the comparison principle to be used. The results show that all elliptic semi-supervised learning algorithms are well-posed at label rates as low as $O(\varepsilon)$, where $\varepsilon>0$ is the bandwidth of the graph (or average distance to the $k^{\rm th}$ nearest neighbor in a $k$-NN graph). This gives a baseline for comparing algorithms at low label rates, and indicates that algorithms should only be claimed to be superior at low label rates if they are well-posed at label rates significantly less than $O(\varepsilon)$. For example, recent work [16] by the second author shows that Laplacian regularization ($p=2$) is well-posed at label rates as low as $O(\varepsilon^{2})$, and previous work [50, 8] showed that the $p$-Laplacian is well-posed for arbitrarily low label rates when $p>d$. It is an open problem to determine the lowest label rates for which $p$-Laplacian regularization is well-posed for $2<p<d$.

1.1.2. Efficient algorithms

We develop and study a range of algorithms for solving both the variational and game-theoretic $p$-Laplace equations, and identify the algorithms that are efficient and scalable in each setting. In each case, we conduct numerical experiments on synthetic data to measure execution time and its dependence on the intrinsic dimensionality of the graph,

For the variational $p$-Laplacian, we propose to use Newton’s method with homotopy on $p$. Other solvers, such as iteratively reweighted least squares (IRLS) and a primal-dual method, were considered in [27]. The primal dual method is slower than Newton with homotopy, and IRLS converges only for $p<3$, though recent work has found ways to bypass this restriction [1].

For the game theoretic $p$-Laplacian, we study a gradient descent-like algorithm, a new Newton-like method, and a semi-implicit method, which is an extension of Oberman’s method [45] to graphs. The Newton-like method with homotopy on $p$, and the semi-implicit method converge faster than gradient descent, but do not have provable convergence guarantees. The gradient descent-like method is slower, though provably convergent (see Theorem 3.7).

All algorithms are presented and analyzed in Section 3, while the numerical experiments on synthetic data are presented in Section 4

1.1.3. Experimental study

We conduct a thorough experimental study of $p$-Laplacian semi-supervised learning on real data, including MNIST [37], Fashion MNIST [55], and Extended MNIST [19]. In particular, we study classification problems with very few labels, and show that graph $p$-Laplacian learning with $p>2$ is superior to Laplacian regularization. Our results show that $p$-Laplacian learning becomes more accurate when provided with more unlabeled data, which confirms the semi-supervised learning paradigm. The experiments on real data are presented in Section 5.

1.1.4. Source code

The code for all numerical experiments is available online333https://github.com/mauriciofloresML/Laplacian_Lp_Graph_SSL.git., and the $p$-Laplacian semi-supervised learning algorithm is implemented in the GraphLearning Python package [12].

2. Continuum limits on k-nearest neighbor graphs

Since random geometric graphs (also called $\varepsilon$-ball graphs) generally have poor sparsity, it is common in practice to use $k$-nearest neighbor ($k$-NN) graphs, were each point is connected to its $k$-nearest neighbors. However, there are very few discrete to continuum or consistency results for graph Laplacians on $k$-NN graphs. The only results on $k$-NN graphs that we are aware of are pointwise consistency results (without rates) [51], $\Gamma$-convergence results for variational problems [28], and spectral convergence rates [14].

We give here a detailed analysis of pointwise consistency and discrete to continuum convergence for graph $p$-Laplace equations on various $k$-NN and $\varepsilon$-ball graphs (Sections 2.2, 2.3, and 2.4). After a careful study, we find that previous observations about $p$-Laplacian regularization on $\varepsilon$-graphs do not hold on $k$-NN graphs. In particular, previous work [21, 8, 50] has shown that $p$-Laplacian regularization forgets the distribution of the unlabeled data as $p\to\infty$, which renders the algorithm unsuitable for semi-supervised learning. We show in the following sections (see Remarks 2.6, 2.11, and 2.15) that this observation is true only for $\varepsilon$-ball graphs, and the situation is completely different for symmetrized $k$-NN graphs. For symmetrized $k$-NN graphs, even Lipschitz learning ($p=\infty$) is sensitive to the distribution of unlabeled data.

We also give, in Section 2.5, a general framework for proving discrete to continuum convergence results for elliptic equations on graphs. The main result, Theorem 2.23, gives a general discrete to continuum convergence result for a wide class of elliptic equations on graphs, and shows that all suitable semi-supervised learning algorithms are well-posed at label rates of $O(\varepsilon)$ and higher. This result indicates that algorithms for semi-supervised learning at low label rates should be judged by their ability to operate below the $O(\varepsilon)$ threshold, and theoretical results should aim to establish this.

Let us remark that the pointwise consistency results in this section share some similarities, in terms of proof techniques, to the results in [8] and [14]. In [8], the second author established pointwise consistency, without any rate, for the game theoretic $p$-Laplacian on $\varepsilon$-ball graphs. In Theorem 2.3 we extend this result to an $O(\varepsilon)$ convergence rate. This is already well-known for $p=2$ (see [33]), and since the $p$-Laplacian is a convex combination of the $p=2$ and $p=\infty$ Laplacians, we quote existing results for $p=2$ (actually, we quote the $p=2$ result for $\varepsilon$-ball graphs from [14] in the proof of Theorem 2.3, since the form of the result from [14] is simpler to use).

For $k$-NN graphs, the only existing pointwise consistency results with rates are the $p=2$ results in the manifold setting established by the second author in [14]. The only other work we are aware of on pointwise consistency for $k$-NN graphs is [51], which considers $p=2$ and proves consistency without any rates. In this section, we prove pointwise consistency with linear convergence rates for the game-theoretic $p$-Laplacian on $k$-NN graphs, both symmetrized (Theorem 2.14) and non-symmetrized (Theorem 2.10). The analogous results for $p=2$ were established by the second author in [14], and these results can be viewed as an extension of those to the $p$-Laplacian with $p>2$. Naturally, we make use of the $p=2$ results from [16], since the game-theoretic $p$-Laplacian is a convex combination of the $p=2$ and $p=\infty$ graph Laplacians. For both $k$-NN and $\varepsilon$-ball graphs, we prove linear rates of $O(\varepsilon)$ (where $\varepsilon=\left(\frac{k}{n}\right)^{1/d}$ for $k$-NN graphs). For $\varepsilon$-ball graphs and $p=2$, it is well-known that sharper $O(\varepsilon^{2})$ pointwise consistency rates are available (see, e.g.,[11, Remark 5.26]). It would be interesting to investigate whether these sharper rates can be extended to $k$-NN graphs.

For reference, we include here a table of notation used in this section.

Notation

• $n$:

Number of vertices in the graph.

• $\varepsilon$:

Graph connectivity length scale for $\varepsilon$-ball graphs.

• $k$:

Number of neighbors in a $k$-NN graph.

• $\eta_{\varepsilon}$:

Rescaled weight kernel $\eta_{\varepsilon}(t)=\eta\left(\tfrac{t}{\varepsilon}\right)$.

• $\mathcal{X}_{n}$:

Vertices of our graphs—an i.i.d. sample of size $n$ with density $\rho$ on $\Omega\subset\mathbb{R}^{d}$.

• $\beta$:

Bound on the density $0<\beta\leq\rho\leq\beta^{-1}$.

• $r_{0}$:

Unique maximum of $r\mapsto r\eta(r)$.

• $d^{n,\varepsilon}$:

Degrees on an $\varepsilon$-ball graph.

• $N_{\varepsilon}(x)$:

The number of points in the set $\mathcal{X}_{n}\cap B(x,\varepsilon)$.

• $\varepsilon_{k}(x)$:

The distance from $x$ to its $k^{\rm th}$ nearest neighbor in $\mathcal{X}_{n}$.

• $s_{k}(x)$:

Typical distance from $x$ to its $k^{\rm th}$ nearest neighbor (satisfies $\alpha(d)ns_{k}(x)^{d}\rho(x)=k$).

• $\varepsilon_{k}(x,y)$:

Symmetrization of $\varepsilon_{k}(x)$, given by $\varepsilon_{k}(x,y)=\max\{\varepsilon_{k}(x),\varepsilon_{k}(y)\}$.

• $s_{k}(x,y)$:

Symmetrization of $s_{k}(x)$, given by $s_{k}(x,y)=\max\{s_{k}(x),s_{k}(y)\}$.

• $\mathscr{L}^{n,\varepsilon}$:

Unnormalized graph Laplacian on an $\varepsilon$-ball graph.

• $\mathscr{L}^{n,\varepsilon}_{rw}$:

Random walk graph Laplacian on an $\varepsilon$-ball graph.

• $\mathscr{L}^{n,k}_{a,rw}$:

Random walk graph Laplacian on a nonsymmetric $k$-NN graph.

• $\mathscr{L}^{n,k}_{s,rw}$:

Random walk graph Laplacian on a symmetric $k$-NN graph.

• $\mathscr{L}^{n,\varepsilon}_{p}$:

Game-theoretic graph $p$-Laplacian on an $\varepsilon$-ball graph.

• $\mathscr{L}^{n,k}_{a,p}$:

Game-theoretic graph $p$-Laplacian on a nonsymmetric $k$-NN graph.

• $\mathscr{L}^{n,k}_{s,p}$:

Game-theoretic graph $p$-Laplacian on a symmetric $k$-NN graph.

2.1. Graph construction

Some parts of the construction of the random graphs are common to $\varepsilon$-ball graphs and $k$-NN graphs, and we review this now. Let $X_{1},X_{2},\dots,X_{n}$ be a sequence of i.i.d. random variables on an open connected domain $\Omega\subset\mathbb{R}^{d}$ with a probability density $\rho:\Omega\to[0,\infty)$. We assume the boundary $\partial\Omega$ is smooth, $\rho\in C^{2}(\overline{\Omega})$ and there exists $\beta>0$ such that

[TABLE]

Let $\partial_{r}\Omega=\{x\in\Omega\,:\,\text{dist}(x,\partial\Omega)\leq r\}$ and $\Omega_{r}=\Omega\setminus\partial_{r}\Omega$. The vertices of our graph are

[TABLE]

Let $\eta:[0,\infty)\to[0,\infty)$ be smooth and nonincreasing such that $\eta(t)\geq 1$ for $0\leq t\leq\frac{1}{2}$, and $\eta(t)=0$ for $t>1$. For $\varepsilon>0$ define $\eta_{\varepsilon}(t)=\eta\left(\tfrac{t}{\varepsilon}\right)$ and set $\sigma_{\eta}=\int_{\mathbb{R}^{n}}|z_{1}|^{2}\eta(|z|)\,dz$. Since the graph is unchanged under scaling the weights by a constant, we may as well assume that

[TABLE]

Then, in particular, $\int_{B(0,\varepsilon)}\eta_{\varepsilon}(z)\,dz=1$ as well. As in [8], we assume there exists $r_{0}\in(0,1)$ and $\theta>0$ so that

[TABLE]

Remark 2.1*.*

The condition (2.4) says that $r\mapsto r\eta(r)$ has a unique maximum, at $r=r_{0}$, and is strongly concave near the maximum. Looking forward to the definition of the graph $\infty$-Laplacian (e.g., (2.8)) and Lemma 2.4, the $\max$ and $\min$ over $y\in\mathcal{X}_{n}$ in (2.8) turn out to occur at a distance of $|x-y|\sim\varepsilon r_{0}$ from $x$, at least asymptotically as $n\to\infty$ and $\varepsilon\to 0$, and the quantitative assumption (2.4) allows us to control this approximation error. If the maximum of $r\mapsto r\eta(r)$ is not unique, then the graph $\infty$-Laplacian may not be asymptotically consistent with the $\infty$-Laplace operator $\Delta_{\infty}$.

It is possible to construct nonincreasing $\eta$ where (2.4) is not satisfied. For example, consider the kernel

[TABLE]

The maximum of $r\mapsto r\eta(r)$ is attained on the interval $[\tfrac{1}{2},\tfrac{3}{4}]$ where $r\eta(r)=1$. We can mollify $\eta$ to obtain a smooth kernel wth the same property on the slightly smaller interval $[\tfrac{1}{2}+\delta,\tfrac{3}{4}-\delta]$ for any $\delta>0$. We remark that the commonly used Gaussian kernel $\eta(r)=e^{-\frac{r^{2}}{2\sigma^{2}}}$ satisfies (2.4) (i.e., the maximum of $f(r):=r\eta(r)$ occurs at $r=\sigma$ and $f^{\prime\prime}(\sigma)=-\frac{2f(r)}{\sigma}<0$.

Remark 2.2*.*

We note it is also common to make the manifold assumption, where $X_{1},X_{2},\dots,X_{n}$ are a sequence of i.i.d. random variables sampled from an $m$-dimensional manifold $\mathcal{M}$ embedded in $\mathbb{R}^{d}$, where $m\ll d$. We expect that the analysis here will carry over to the manifold setting with additional technical details.

2.2. $\varepsilon$-ball graphs

The graph constructed with vertices $\mathcal{X}_{n}$ and edge weights $w_{xy}=\eta_{\varepsilon}(|x-y|)$ is called a random geometric graph, or sometimes an $\varepsilon$-ball graph, and is the most widely used graph construction in theoretical analysis of graph-based learning algorithms. Here, we review pointwise consistency for the $p$-Laplacian on $\varepsilon$-ball graphs. Much of the consistency theory was established previously in [8], so this section is mostly review, with the additional observation that the arguments in [8] establish pointwise consistency rates for the $\infty$-Laplacian.

We define the graph Laplacian on a random geometric graph by

[TABLE]

We also define the degree

[TABLE]

The random walk graph Laplacian is defined by

[TABLE]

The random walk Laplacian is the generator for a random walk on the graph with transition probabilities $\eta_{\varepsilon}(|x-y|)/d^{n,\varepsilon}(x)$ of transitioning from $x$ to $y$. The graph $\infty$-Laplacian is defined by

[TABLE]

and the game-theoretic graph $p$-Laplacian is defined by

[TABLE]

In the continuum, $\mathscr{L}^{n,\varepsilon}_{p}$ is consistent with the weighted $p$-Laplace operator

[TABLE]

as is shown in the following theorem.

Theorem 2.3** (Consistency on $\varepsilon$-graphs).**

There exists $C_{1},C_{2},C_{3}>0$ such that for any $\varepsilon>0$ with $n\varepsilon^{d}\geq 1$ and any $0<\lambda\leq 1$, the event that

[TABLE]

holds for all $u\in C^{3}(\overline{\Omega})$ has probability at least $1-C_{2}n\exp\left(-C_{3}n\varepsilon^{3d/2}\lambda^{d/2}\right)-C_{2}n\exp\left(-C_{3}n\varepsilon^{d+2}\lambda^{2}\right)$.

Proof.

It was shown in [8, Theorem 5] that the event that

[TABLE]

holds for all $u\in C^{3}(\overline{\Omega})$ has probability at least $1-Cn\exp\left(-cn\varepsilon^{d+2}\lambda^{2}\right)$ for constants $C,c>0$ and any $0<\lambda\leq 1$. This is a uniform (over functions $u$) version of pointwise consistency for the graph Laplacian. The weaker nonuniform version dates back to results in [34, 33]. The proof follows by combining (2.12) with Lemma 2.4 below. ∎

Lemma 2.4**.**

There exists $C_{1},C_{2},C_{3}>0$ such that for any $\varepsilon>0$ with $n\varepsilon^{d}\geq 1$, the event that

[TABLE]

holds for all $u\in C^{3}(\overline{\Omega})$ has probability at least $1-C_{2}n\exp\left(-C_{3}n\varepsilon^{3d/2}\lambda^{d/2}\right)$.

Before proving Lemma 2.4, we make a few remarks.

Remark 2.5*.*

In fact, pointwise consistency of the graph $\infty$-Laplacian (Lemma 2.4) requires $\nabla u(x)\neq 0$, since the $\infty$-Laplacian is discontinuous (as a function of $\nabla u$) at $\nabla u(x)=0$. We interpret the right hand side of (2.13) to be $\infty$ if $\nabla u(x)=0$. The viscosity solution framework for proving discrete to continuum convergence does not require consistency at points where $\nabla u(x)=0$, since the viscosity sub- and supersolution conditions are not required to hold at such points (see, e.g., [8]).

Remark 2.6*.*

We note that the continuum operator $\Delta_{p}$ forgets about the distribution $\rho$ as $p\to\infty$. This was first observed in [21], and proved rigorously in [8, 50]. On the other hand, for semi-supervised learning problems with very few labels, the theory in [21, 8, 50] suggests one should use $d<p<p+\delta$ for some small $\delta>0$, so that the algorithm is well-posed with very few labels, and still has maximal dependence on the data distribution $\rho$, which is essential for semi-supervised learning. We show in Section 2.4 that this observation is completely dependent on the $\varepsilon$-ball graph construction, and in particular, the story is very much different on $k$-NN graphs.

We now turn to the proof of Lemma 2.4. By [8, Theroem 6], we know that

[TABLE]

holds for all $u\in C^{3}(\overline{\Omega})$ and $x\in\Omega$ with $\nabla u(x)\neq 0$ with probability one, provided $\varepsilon_{n}\to 0$ so that

[TABLE]

Lemma 2.4 strengthens this to pointwise consistency with a convergence rate. For this, we need some additional notation. We define

[TABLE]

We need a couple of preliminary results before proving Lemma 2.4. The first concerns the approximation of the maximum value of a function over $\Omega$ with the maximum over the point cloud $\mathcal{X}_{n}$.

Proposition 2.7**.**

Let $\psi\in C^{2}(\overline{\Omega})$ and assume $\psi$ attains its maximum value over $\overline{\Omega}$ at some $x_{0}\in\Omega_{\delta_{n}}$. Then we have

[TABLE]

Proof.

Since $\nabla\psi(x_{0})=0$ and $\psi$ is $C^{2}$, for any $x\in\mathcal{X}_{n}$ we have

[TABLE]

and therefore

[TABLE]

The second preliminary result we quote directly from [10, Proposition 4.2].

Proposition 2.8**.**

For any $p\in\mathbb{R}^{d}$ with $p\neq 0$ and $A\in\mathbb{R}^{d\times d}$, we have

[TABLE]

We note that we prove a generalization of Proposition 2.8 in Section 2.4 when studying symmetrized $k$-NN graphs (see Proposition 2.20).

We now give the proof of Lemma 2.4.

Proof of Lemma 2.4.

First, we claim that

[TABLE]

where

[TABLE]

To see this, consider the $\max$ term in $\mathscr{L}^{n,\varepsilon}_{\infty}u(x)$,

[TABLE]

We use Proposition 2.7 to obtain

[TABLE]

Applying Proposition 2.8 we have

[TABLE]

since $\Delta_{\infty}u(x)=|\nabla u(x)|^{-2}\nabla u(x)^{T}\nabla^{2}u(x)\nabla u(x)$. We can apply a similar argument to the $\min$ term from $\mathscr{L}^{n,\varepsilon}_{\infty}u(x)$, and this establishes the claim (2.18).

Now, let $r^{\pm}\in[0,1]$ such that

[TABLE]

By (2.4) we have

[TABLE]

Therefore

[TABLE]

Since $r\mapsto r^{2}\eta(r)$ is Lipschitz continuous we have

[TABLE]

and we deduce that

[TABLE]

Since

[TABLE]

we have

[TABLE]

Similarly, since

[TABLE]

we have

[TABLE]

Combining this with (2.18) and Proposition 2.9 below completes the proof. ∎

Proposition 2.9**.**

There exists $C_{1},C_{2}>0$ such that for every $t>0$ with $nt^{d}\geq 1$ we have $\mathbb{P}(\delta_{n}>t)\leq C_{1}n\exp\left(-C_{2}\beta nt^{d}\right)$.

Proof.

Let $t>0$ with $nt^{d}\geq 1$. Since $\Omega$ has a smooth (and hence Lipschitz) boundary, for any $h>0$ we can find a covering of $\Omega$ by balls $B(x_{1},h),B(x_{2},h),\dots,B(x_{M},h)$ so that $|B(x_{i},h)\cap\Omega|\geq C_{1}h^{d}$ and $M\leq C_{2}h^{-d}$ for some universal constants $C_{1},C_{2}>0$. Since $\rho\geq\beta>0$, the probability that $B(x_{i},h)\cap\Omega\cap\mathcal{X}_{n}$ is empty is bounded by $(1-C_{1}\beta h^{d})^{n}\leq\exp(-C_{1}\beta nh^{d})$. Hence, the event that $B(x_{i},h)\cap\Omega\cap\mathcal{X}_{n}$ has at least one point for all $i=1,\dots,M$ has probability at least $1-C_{2}h^{-d}\exp(-C_{1}\beta nh^{d})$. Since the balls $B(x_{i},h)$ cover $\Omega$, for each $x\in\Omega$ there exists $x_{i}$ such that $|x-x_{i}|\leq h$. If $B(x_{i},h)\cap\Omega\cap\mathcal{X}_{n}$ is nonempty, then there exists $y\in\mathcal{X}_{n}$ such that $|x-y|\leq 2h$. Thus, if $\delta_{n}>2h$ then at least one of $B(x_{i},h)\cap\Omega\cap\mathcal{X}_{n}$ is empty. Therefore

[TABLE]

The proof is completed by setting $h=t/2$ and recalling $h^{-d}=2^{d}t^{-d}\leq 2^{d}n$. ∎

2.3. Nonsymmetric k-nearest neighbor graphs

We now consider the simplest $k$-nearest neighbor graph. Let

[TABLE]

be the number of samples in an $\varepsilon$-neighborhood of $x\in\Omega$. Here, $B(x,\varepsilon)$ is the closed ball of radius $\varepsilon>0$ centered about $x$. For $k\leq n$ we define

[TABLE]

The value of $\varepsilon_{k}(x)$ for $x\in\mathcal{X}_{n}$ is the distance from $x$ to its $k^{\rm th}$ nearest neighbor in $\mathcal{X}_{n}$. The nonsymmetric $k$-nearest neighbor random walk graph Laplacian is then given by

[TABLE]

We note that since $\varepsilon_{k}(x)$ satisfies, on average, $\alpha(d)n\varepsilon_{k}(x)^{d}\rho(x)\approx k$, the normalization in (2.21) is equivalent, up to the factor $\rho(x)$, with the normalization in the case of the $\varepsilon$-ball graph given in (2.7).

The graph $\infty$-Laplacian is defined by

[TABLE]

We now define the game-theoretic $k$-nearest neighbor graph $p$-Laplacian to be

[TABLE]

The main result in this section is consistency for the nonsymmetric $k$-nearest neighbor graph Laplacian.

Theorem 2.10** (Consistency on nonsymmetic $k$-NN graphs).**

There exists $C_{1},C_{2},C_{3},c_{1},c_{2}>0$ such that for $k\leq c_{1}n$ and $0<\lambda\leq 1/4$, the event that

[TABLE]

holds for all $u\in C^{3}(\overline{\Omega})$ and $x\in\Omega_{c_{2}(k/n)^{1/d}}\cap\mathcal{X}_{n}$, where $B_{k}=B(x,c_{2}(k/n)^{1/d})$, has probability at least

[TABLE]

Remark 2.11*.*

Notice that, up to the factor $\rho^{-2/d}$, the nonsymmetric $k$-NN graph $p$-Laplacian is consistent with the same operator $\Delta_{p}$ as the $\varepsilon$-ball graph Laplacian. This is due to the ability to treat the $k$-NN Laplacian as an $\varepsilon$-ball graph Laplacian with spatially varying $\varepsilon$, since there is no symmetrization to be concerned with. The same observation as in Remark 2.6 holds here; that is, the continuum operator $\Delta_{p}$ forgets the distribution $\rho$ as $p\to\infty$. We also remark that the rates in Theorems 2.10 and 2.3 are essentially the same, since the quantity $\left(\frac{k}{n}\right)^{1/d}$ is the average distance to the $k^{\rm th}$ nearest neighbor (obtained by setting $n\varepsilon^{d}\sim k$), and plays the same role as $\varepsilon$ in Theorem 2.3.

We now turn to the proof of Theorem 2.10. To connect $\mathscr{L}^{n,k}_{a,p}$ to $\mathscr{L}^{n,\varepsilon}_{p}$, we define

[TABLE]

The radius $s_{k}(x)$ satisfies $\alpha(d)ns_{k}(x)^{d}\rho(x)=k$, so that the ball $B(x,s_{k}(x))$ contains on average $k$ points from $\mathcal{X}_{n}$. Then by algebraic manipulations we have

[TABLE]

This identity allows us to view the nonsymmetric $k$-nearest neighbor graph Laplacian $\mathscr{L}^{n,k}_{a,p}$ as an $\varepsilon$-ball graph Laplacian for spatially varying $\varepsilon=\varepsilon(x)$, and consistency will follow by applying Theorem 2.3 and a covering argument.

Before proving consistency, we recall some facts about $N_{\varepsilon}(x)$ and $\varepsilon_{k}(x)$. These facts can be found in [14], in the manifold setting, but we include the proofs here for the reader’s convenience, as they are simpler in the Euclidean setting.

Proposition 2.12**.**

For any $0<t\leq 1$ and $x\in\Omega_{\varepsilon}$ we have

[TABLE]

where $p(x,\varepsilon)=\int_{B(x,\varepsilon)}\rho(y)\,dy$.

Proof.

Note that $N_{\varepsilon}(x)$ is the sum of Bernoulli zero/one random variables with parameter $p(x,\varepsilon)=\int_{B(x,\varepsilon)}\rho(x)\,dx$. By the Chernoff bounds we have

[TABLE]

for any $0<t\leq 1$. The proof is completed by noting that $p(x,\varepsilon)\geq\alpha(d)\beta\varepsilon^{d}$. ∎

Lemma 2.13**.**

There exists $C_{1},C_{2},c_{1},c_{2}>0$ such that for any $0<t\leq 1/4$, $k\leq c_{1}n$, and $x\in\Omega_{c_{2}(k/n)^{1/d}}$, we have

[TABLE]

Proof.

For any $\varepsilon>0$ we have

[TABLE]

Choose $\varepsilon>0$ so that $np(x,\varepsilon)=k(1+t)$ for $t>0$. Then we have

[TABLE]

provided $t\leq 1$. By Proposition 2.12 we have

[TABLE]

for all $0<t\leq 1$. We now note that

[TABLE]

and

[TABLE]

Therefore

[TABLE]

For a more refined estimate, we have

[TABLE]

Therefore

[TABLE]

and so

[TABLE]

provided $C\varepsilon^{2}\leq\frac{1}{2}$. Hence, due to (2.28), there exists $c>0$ such that $k\leq cn$ implies $C\varepsilon^{2}\leq\frac{1}{2}$ and so

[TABLE]

It follows that there exists $C_{1},C_{2}>0$ such that

[TABLE]

for any $0<t\leq 1$. The proof for the estimate in the other direction is similar, but requires a restriction $t\leq 1/4$. ∎

Proof of Theorem 2.10.

Fix $0<a<b\leq 1$ with $na^{d}\geq 1$. Let $C_{1},C_{2},C_{3}>0$ so that Theorem 2.3 holds and let $x\in\Omega_{b}$. For $0<\delta\leq a/2$ and $0\leq\lambda\leq 1$, let $A_{\delta,\lambda}$ denote the event that

[TABLE]

holds for all $u\in C^{3}(\overline{\Omega})$ and all $\varepsilon\in[a-\delta,b+\delta]\cap\mathbb{Z}_{\delta}$, where $\mathbb{Z}_{\delta}=\delta\mathbb{Z}$. By Theorem 2.3 and a union bound we have

[TABLE]

For $\varepsilon>0$, let us denote by $\lfloor\varepsilon\rfloor_{\delta}$ the largest number belonging to the set $\mathbb{Z}_{\delta}:=\delta\mathbb{Z}$ that is less than or equal to $\varepsilon$, and write $\lceil\varepsilon\rceil_{\delta}=\lfloor\varepsilon\rfloor_{\delta}+\delta$. Noting that

[TABLE]

provided $|x-y|\leq\varepsilon$, we can compute that

[TABLE]

holds for all $u\in C^{3}(\overline{\Omega})$ and $\varepsilon\in[a-\delta,b+\delta]$ with probability at least $1-2\exp\left(-cna^{d}\right)$ for $c>0$, due to Proposition 2.12. Therefore, when $A_{\delta,\lambda}$ occurs and $\varepsilon_{k}(x)\in[a,b]$ we have

[TABLE]

provided we choose $\delta=a^{3}$. Using (2.25) we can write

[TABLE]

We now invoke Lemma 2.13, set $a=c_{1}(k/n)^{1/d}$ and $b=c_{2}(k/n)^{1/d}$ for constants $0<c_{1}<c_{2}$, and union bound over $x\in\Omega_{c_{2}(k/n)^{1/d}}\cap\mathcal{X}_{n}$. The proof is completed by noting that $\delta^{-1}\leq Cn^{3/d}\leq Cn^{2}$. ∎

2.4. Symmetric k-nearest neighbor graphs

We now consider symmetrized $k$-nearest neighbor graphs. Define

[TABLE]

The symmetric $k$-nearest neighbor random walk graph Laplacian is then defined by

[TABLE]

where the degree $d^{n,k}_{s}(x)$ is defined by

[TABLE]

We can also define a symmetrized graph using $\varepsilon_{k}(x,y)=\min\{\varepsilon_{k}(x),\varepsilon_{k}(y)\}$, which would produce the mutual $k$-NN graph, or any other suitable symmetric combination of $\varepsilon_{k}(x)$ and $\varepsilon_{k}(y)$. We expect the consistency results in this section to hold, with minor modifications, for many other types of symmetrization.

The graph $\infty$-Laplacian on the symmetric $k$-nearest neighbor graph is defined by

[TABLE]

and the graph $p$-Laplacian is defined by

[TABLE]

The graph $p$-Laplacian on a symmetric $k$-NN graph is consistent in the continuum with the operator

[TABLE]

as is shown in the following theorem.

Theorem 2.14** (Consistency on symmetric $k$-NN graphs).**

There exists $C_{1},C_{2},C_{3},c_{1},c_{2}>0$ such that for $0<t\leq 1/4$ and $k\leq c_{1}n$, the event that

[TABLE]

holds for all $x\in\mathcal{X}_{n}\cap\Omega_{\varepsilon}$, where $\varepsilon=c_{2}(k/n)^{1/d}$, and for all $u\in C^{3}(\overline{\Omega})$, has probability at least $1-C_{3}n\exp\left(-C_{2}k\left(\tfrac{k}{n}\right)^{2/d}t^{2}\right)$.

Proof.

The proof simply combines Lemma 2.18 and Lemma 2.21, both of which are proved below. ∎

Remark 2.15*.*

It is important to point out the additional drift term $-\frac{1}{d}\nabla\log\rho\cdot\nabla u$ in $\Delta_{p}^{s}$, compared to $\Delta_{p}$. This term is due to the symmetrization of the $k$-nearest neighbor relation, and is a result of the Taylor expansion

[TABLE]

which makes an appearance in the proof of Lemma 2.16 below. Since this term persists in the limit as $p\to\infty$, this shows that $p$-Laplacian learning for large $p$, and in particular, Lipschitz learning ($p=\infty$) are sensitive to the distribution $\rho$ of the unlabeled data provided the graph is constructed as a symmetrized $k$-nearest neighbor graph, as is often done in practice. This should be contrasted with the case of $\varepsilon$-ball graphs and nonsymmetric $k$-NN graphs, where the continuum operator $\Delta_{p}$ forgets the distribution $\rho$ as $p\to\infty$.

We first record a result about $\varepsilon_{k}(x,y)$, which was proved in [14], but not explicitly stated as a lemma. We define $a_{-}=\min\{a,0\}$.

Lemma 2.16**.**

There exists $C_{1},C_{2},c_{1},c_{2}>0$ such that for any $0<t\leq 1/4$, $k\leq c_{1}n$ and $x\in\Omega_{\varepsilon}$, the event that

[TABLE]

holds for all $y\in\mathcal{X}_{n}\cap B(x,\varepsilon)$, with $\varepsilon=c_{2}(k/n)^{1/d}$, has probability at least $1-2n\exp\left(-C_{2}kt^{2}\right)$.

Proof.

Let $x\in\Omega_{\varepsilon}$. By Lemma 2.13 and a union bound, there exists $C_{1},C_{2},c>0$ such that for any $0<t\leq 1/4$ and $k\leq cn$ we have

[TABLE]

holds with probability at least $1-2n\exp\left(-C_{2}kt^{2}\right)$, where $\varepsilon=c_{2}(k/n)^{1/d}$. We assume for the rest of the proof that (2.36) holds. For $y\in\mathcal{X}_{n}\cap B(x,\varepsilon)$ it follows from (2.36) that

[TABLE]

where $s_{k}(x,y)=\max\{s_{k}(x),s_{k}(y)\}$. Since $|x-y|\leq\varepsilon$, we compute

[TABLE]

Therefore

[TABLE]

which completes the proof. ∎

Lemma 2.17**.**

There exists $C_{1},C_{2},C_{3},c>0$ such that for each $x\in\Omega_{\varepsilon}$ with $\varepsilon=c(k/n)^{1/d}$

[TABLE]

holds with probability at least $1-C_{2}\exp\left(-C_{3}kt^{2}\right)$ for any $0<t\leq 1/4$.

Proof.

By Lemma 2.16 there exists $c>0$ so that for every $x\in\Omega_{\varepsilon}$, where $\varepsilon=c(k/n)^{1/d}$, we have

[TABLE]

for all $y\in\mathcal{X}_{n}\cap B(x,\varepsilon)$ with probability at least $1-2n\exp\left(-C_{2}kt^{2}\right)$. Notice the drift term involving $\nabla\log\rho(x)\cdot(y-x)_{-}$ form Lemma 2.16 is absorbed into the $O(\varepsilon)$ error term on the right hand side, since $|x-y|\leq\varepsilon$. Assuming (2.38) holds we have

[TABLE]

for $C>0$, provided $t\leq 1$ and $\varepsilon\leq 1$. By Proposition 2.12 we have $N_{Cs_{k}(x)}(x)\leq\overline{C}k$ with probability at least $1-\exp\left(-ck\right)$. Thus,

[TABLE]

By the Bernstein inequality (see, e.g., [8, Theorem 5]) we have that

[TABLE]

with probability at least $1-2\exp\left(-ckt^{2}\right)$. This completes the proof. ∎

We now prove pointwise consistency for the random walk graph Laplacian $\mathscr{L}^{n,k}_{s,rw}$ on a symmetric $k$-nearest neighbor graph. The analogous result for the unnormalized graph Laplacian was established in [14].

Lemma 2.18**.**

There exists $C_{1},C_{2},C_{3},c_{1},c_{2}>0$ such that for $0<t\leq 1/4$ and $k\leq c_{1}n$, the event that

[TABLE]

holds for all $u\in C^{3}(\overline{\Omega})$ and $x\in\Omega_{c_{2}(k/n)^{1/d}}\cap\mathcal{X}_{n}$, where $B_{k}=B(x,c_{2}(k/n)^{1/d})$, has probability at least $1-C_{2}n\exp\left(-C_{3}\left(\tfrac{k}{n}\right)^{2/d}kt^{2}\right)$.

Proof.

Let $u\in C^{3}(\overline{\Omega})$ and $x\in\Omega$. Let us define the unnormalized $k$-nearest neighbor graph Laplacian

[TABLE]

Then by [14, Theorem 3.6] we have that

[TABLE]

holds for all $x\in\Omega_{\varepsilon}\cap\mathcal{X}_{n}$ with probability at least $1-Cn\exp\left(-c\left(\tfrac{k}{n}\right)^{2/d}kt^{2}\right)$ for $0\leq t\leq 1/4$, $k\leq c_{2}n$, and $\varepsilon=c_{2}(k/n)^{1/d}$ with $c_{1},c_{2}>0$. Noting that

[TABLE]

and invoking Lemma 2.17 we have

[TABLE]

holds for all $x\in\Omega_{\varepsilon}\cap\mathcal{X}_{n}$ with probability at least $1-Cn\exp\left(-c\left(\tfrac{k}{n}\right)^{2/d}kt^{2}\right)$.

To prove uniformity over $u\in C^{3}(\overline{\Omega})$, we fix $x\in\Omega_{\varepsilon}$ and $u\in C^{3}(\overline{\Omega})$ and we first Taylor expand $u$ to find

[TABLE]

for $|x-y|\leq\varepsilon$. Therefore

[TABLE]

where $p_{i}(y)=y_{i}-x_{i}$ and the error term is controlled with probability at least $1-C\exp\left(-ck\right)$, due to Lemma 2.17. By the argument above

[TABLE]

and

[TABLE]

hold for all $x\in\Omega_{\varepsilon}\cap\mathcal{X}_{n}$ with probability at least $1-Cn\exp\left(-c\left(\tfrac{k}{n}\right)^{2/d}kt^{2}\right)$, where $\delta_{ij}=1$ if $i=j$ and $\delta_{ij}=0$ otherwise. Substituting this into (2.39) we have

[TABLE]

which completes the proof. ∎

We require the following simple proposition.

Proposition 2.19**.**

Let $\psi\in C^{1}(\overline{\Omega})$ and assume $\psi$ attains its maximum value over $\overline{\Omega}$ at some $x_{0}\in\Omega$. Then we have

[TABLE]

Proof.

For any $x\in\mathcal{X}_{n}$ we have

[TABLE]

and therefore

[TABLE]

where we recall $\delta_{n}$ was defined in (2.16). ∎

We define

[TABLE]

The following proposition can be viewed as an extension of Proposition 2.8, which applies when $v=0$. It is used to evaluate the $\max$ and $\min$ in the graph $\infty$-Laplacian, asymptotically after Taylor expansions, in Lemma 2.21 below.

Proposition 2.20**.**

For any $\varepsilon,r>0$, $p\in\mathbb{R}^{d}$ with $p\neq 0$ and $X\in\mathbb{R}^{d\times d}$ we have

[TABLE]

Proof.

Let $z_{r}$ such that $|z_{r}|=r$ and

[TABLE]

Set $w_{r}=rp/|p|$. Choosing $z=w_{r}$ yields

[TABLE]

Note that for any unit vectors $a,b$ we have

[TABLE]

Therefore

[TABLE]

Combining this with the observations above we have

[TABLE]

Therefore

[TABLE]

and so

[TABLE]

which completes the proof. ∎

Lemma 2.21**.**

There exists constants $C_{1},C_{2},c_{1},c_{2}>0$ such that for $0<t\leq 1$ and $k\leq c_{1}n$, the event that

[TABLE]

holds for all $x\in\mathcal{X}_{n}\cap\Omega_{\varepsilon}$, where $\varepsilon=c_{2}(k/n)^{1/d}$, and for all $u\in C^{3}(\overline{\Omega})$, has probability at least $1-2n\exp\left(-C_{2}k\left(\tfrac{k}{n}\right)^{2/d}t^{2}\right)$.

Proof.

By Lemma 2.16 we have

[TABLE]

for $|x-y|\leq\varepsilon_{k}(x,y)$, where $v=-\tfrac{1}{d}\nabla\log\rho(x)$, with probability at least $1-2n\exp\left(-C_{2}k\varepsilon^{2}t^{2}\right)$, where $0<t\leq 1$ and $\varepsilon=c(k/n)^{1/d}$. Therefore, by Proposition 2.19 we have

[TABLE]

We now make a change of variables, setting

[TABLE]

For $y$ sufficiently close to $x$, depending only on $|v|$, the mapping $y\mapsto\Phi(y)$ is invertible and

[TABLE]

where we write $s_{k}=s_{k}(x)$ for simplicity, and note that $s_{k}=O(\varepsilon)$. Therefore, using Proposition 2.20 we have

[TABLE]

Let us set

[TABLE]

and let $r_{1}\in[0,1]$ so that

[TABLE]

By (2.4) we have

[TABLE]

It follows that

[TABLE]

and so

[TABLE]

This yields

[TABLE]

Inserting this above we have

[TABLE]

We can apply the same argument as above to $-u$ to obtain

[TABLE]

Therefore

[TABLE]

∎

2.5. Discrete to continuum convergence

Many types of discrete to continuum convergence results have been proven recently for various learning problems on graphs, using tools like $\Gamma$-convergence [52, 28, 29], the maximum principle and viscosity solutions [8, 10], and even Martingale techniques [16]. The $\Gamma$-convergence results are variational in nature and do not apply to the game theoretic $p$-Laplacian. The viscosity solution approach in [8] does not require a variational structure, but used uniform equicontinuity of the sequence of learned functions to establish compactness.

We present here a very general technique for proving discrete to continuum convergence results in a general setting that applies to all the symmetric and nonsymmetric operators discussed in Sections 2.2, 2.3 and 2.4, and any others that satisfy certain monotonicity properties. The framework only requires the graph problem to have a maximum principle, and to be pointwise consistent with a well-posed equation in the continuum. It is not necessary to prove that the sequence of functions is equicontinuous, or that the solution of the continuum PDE is smooth. The framework is essentially an adaptation of the Barles-Souganidis framework [5] to convergence of discrete problems on graphs.

Let $\Omega\subset\mathbb{R}^{d}$ be an open and bounded domain. For each $n\geq 1$, let $\mathcal{X}_{n}\subset\Omega$ be a collection of $n$ points in $\Omega$, and let $L^{2}(\mathcal{X}_{n})$ denote the space of functions $u:\mathcal{X}_{n}\to\mathbb{R}$. Let

[TABLE]

denote our graph operator, generalizing a graph Laplacian, and suppose $u_{n}\in L^{2}(\mathcal{X}_{n})$ is a solution of the boundary value problem

[TABLE]

where $g:\Omega\to\mathbb{R}$ is a given continuous function, and $\varepsilon_{n}>0$ represents the length scale on which the graph is connected. For convenience, we recall that $\partial_{r}\Omega=\{x\in\Omega\,:\,\text{dist}(x,\partial\Omega)\leq r\}$ and $\Omega_{r}=\Omega\setminus\partial_{r}\Omega$. The graph equation (2.42) represents our generalization of a semi-supervised learning problem on a graph with labels $g$ on the set $\mathcal{X}_{n}\cap\partial_{\varepsilon_{n}}\Omega$. This is just one model for labeled data, and others are possible (see [16] for other models).

We now lay out simple and general conditions on $\mathcal{X}_{n}$, $u_{n}$ and $F_{n}$ that ensure $u_{n}$ has a well-posed continuum limit; that is, $u_{n}$ converges uniformly to the solution of a continuum PDE. We say the points $\mathcal{X}_{n}$ are space filling if

[TABLE]

We say $F_{n}$ is monotone if for all $u,v\in L^{2}(\mathcal{X}_{n})$, $t\in\mathbb{R}$ and $x\in\Omega$

[TABLE]

We say $F_{n}$ is proper if for all $u\in L^{2}(\mathcal{X}_{n})$, $s,t\in\mathbb{R}$, and $x\in\Omega$

[TABLE]

We say the operators are consistent as $n\to\infty$ with the differential operator444$\mathbb{S}(d)$ denotes the space of symmetric real-valued matrices of size $d\times d$.

[TABLE]

if for every $\varphi\in C^{\infty}(\mathbb{R}^{d})$ and every sequence of real numbers $\xi_{n}\to 0$ we have

[TABLE]

We note that the connectivity length scale $\varepsilon_{n}$ of the graph is encoded into the consistency statement, since we do not assume consistency near the boundary, where the ball $B(x,\varepsilon_{n})$ overlaps with $\partial\Omega$.

Remark 2.22*.*

In the context of Sections 2.2, 2.3 and 2.4, we would set $F_{n}(u,u(x),x)=-\mathscr{L}u(x)$, where $\mathscr{L}$ is any of the graph Laplacians in those sections. For example, if $\mathscr{L}$ is the unnormalized graph Laplacian (1.2), then

[TABLE]

The first argument of $F_{n}$ encodes the dependence of $F_{n}$ on the neighboring values $u(y)$ for $y\in\mathcal{X}_{n}$ with $y\neq x$, and the second argument $t$ encodes the dependence on $t=u(x)$. Since all graph Laplacians are increasing functions of the difference $u(y)-t$, they are both monotone and proper. The consistency results established in Sections 2.2, 2.3 and 2.4 show that (2.46) holds with probability one for any choice of $\mathscr{L}$ from those sections, provided $\varepsilon=\varepsilon_{n}\to 0$ sufficiently slowly. For example, in Theorem 2.3, in order to make sure the results hold with high probability we require

[TABLE]

This can be rewritten as lower bounds on the length scale $\varepsilon$ as follows

[TABLE]

We now follow the Barles-Souganidis framework [5] to prove that the sequence $u_{n}$ converges uniformly to the solution of the boundary value problem

[TABLE]

provided the equation (2.48) is well-posed in the viscosity sense with generalized Dirichlet condition $u=g$ on $\partial\Omega$. We review the definition of viscosity solution and the generalized Dirichlet problem in the appendix in Section A.

Theorem 2.23**.**

Assume $g$ is continuous, (2.48) enjoys strong uniqueness, $F_{n}$ is monotone, proper, and consistent with $F$, $\mathcal{X}_{n}$ is space filling, and $\varepsilon_{n}\to 0$ as $n\to\infty$. Let $u_{n}$ be a sequence of solutions of (2.42) that are uniformly bounded. Then

[TABLE]

where $u\in C(\overline{\Omega})$ is the unique viscosity solution of (2.48).

Remark 2.24*.*

In Theorem 2.23, by uniformly bounded we mean that there exists $C>0$ such that $\max_{x\in\mathcal{X}_{n}}|u_{n}(x)|\leq C$ for all $n\geq 1$. Bounds of this nature are generally proved using the discrete maximum principle on the graph, and often rely on graph-connectivity. For a simple example, consider

[TABLE]

For an appropriate choice of normalization parameter $\alpha_{n}$, the equation $F_{n}=0$ would be consistent with the elliptic equation $u-\Delta u=1$ in the continuum limit. To establish the uniform bound on $u_{n}$, we use a maximum principle argument. Let $x\in\mathcal{X}_{n}$ be a point where $u_{n}$ attains its maximum value over $\mathcal{X}_{n}$. If $x\in\partial_{\varepsilon_{n}}\Omega$, then $u_{n}(x)\leq g(x)$. If $x\in\Omega_{\varepsilon_{n}}$, then we use that $u_{n}(x)-u_{n}(y)\geq 0$ for any $y\in\mathcal{X}_{n}$ to obtain

[TABLE]

Therefore $u_{n}(x)\leq 1$ and we obtain the bound

[TABLE]

A bound on the minimum of $u_{n}$ is obtained similarly.

Proof of Theorem 2.23.

We define the weak upper and lower limits

[TABLE]

Due to (2.43) and the assumption that $u_{n}$ are uniformly bounded, we have $-\infty<\underline{u}(x)\leq\overline{u}(x)<\infty$ for all $x\in\overline{\Omega}$. Furthermore, it is straightforward to check (see, e.g., [9]) that $\overline{u}\in\text{USC}(\overline{\Omega})$ and $\underline{u}\in\text{LSC}(\overline{\Omega})$, where $\text{USC}(\overline{\Omega})$ and $\text{LSC}(\overline{\Omega})$ denote the spaces of upper semicontinuous and lower semicontinuous functions on $\overline{\Omega}$, respectively. We claim that $\overline{u}$ is a viscosity subsolution of (2.48) and $\underline{u}$ is a viscosity supersolution. Once we establish this, it follows from strong uniqueness that $\overline{u}\leq\underline{u}$. Therefore $\overline{u}=\underline{u}$ and (2.49) immediately follows.

We will show that $\overline{u}$ is a viscosity subsolution; the proof that $\underline{u}$ is a supersolution is similar. Let $x_{0}\in\overline{\Omega}$ and $\varphi\in C^{\infty}(\mathbb{R}^{d})$ such $\overline{u}-\varphi$ has a local maximum at $x_{0}$ with respect to $\overline{\Omega}$. Define

[TABLE]

Then $\psi(x_{0})=\overline{u}(x_{0})$, and we can choose $C>0$ large enough so that $\overline{u}-\psi$ has a strict global maximum at $x_{0}$ relative to $\overline{\Omega}$. It follows that there exists $n_{k}\to\infty$, $\mathcal{X}_{n_{k}}\ni x_{n_{k}}\to x_{0}$ with $u_{n_{k}}(x_{n_{k}})\to\overline{u}(x_{0})$ such that $u_{n_{k}}-\psi$ attains its maximum value over $\mathcal{X}_{n_{k}}$ at the point $x_{n_{k}}$ for each $k$. Set $\xi_{k}=u_{n_{k}}(x_{n_{k}})-\psi(x_{n_{k}})$ so that $u_{n_{k}}\leq\psi+\xi_{k}$. Then since $F_{n}$ is monotone (i.e., (2.44) holds) we have

[TABLE]

We now have two cases.

1. If $x_{0}\in\Omega$ then $x_{n_{k}}\in\mathcal{X}_{n_{k}}\cap\Omega_{\varepsilon_{n_{k}}}$ for $k$ sufficiently large, and since $\xi_{k}\to 0$ as $k\to\infty$ we can combine consistency (i.e., (2.46)) with (2.50) to find that

[TABLE]

Since $\nabla^{2}\psi(x_{0})=\nabla^{2}\varphi(x_{0})$, $\nabla\psi(x_{0})=\nabla\varphi(x_{0})$ and $\psi(x_{0})=\overline{u}(x_{0})$ we have, as desired, the viscosity subsolution condition

[TABLE]

2. If $x_{0}\in\partial\Omega$, then we can pass to a further subsequence, if necessary, so that either $x_{n_{k}}\in\mathcal{X}_{n_{k}}\cap\Omega_{\varepsilon_{n_{k}}}$ or $x_{n_{k}}\in\mathcal{X}_{n_{k}}\cap\partial_{\varepsilon_{n_{k}}}\Omega$ for all $k$. In the former case we again find that (2.51) holds. In the latter case, we have $u_{n_{k}}(x_{n_{k}})=g(x_{n_{k}})$ for all $k$, and since $g$ is continuous we have that $\overline{u}(x_{0})=g(x_{0})$. Thus, when $x\in\partial\Omega$ we have

[TABLE]

which is the viscosity subsolution condition on the boundary. This completes the proof. ∎

Remark 2.25*.*

In the context of semi-supervised learning, Theorem 2.23 shows that the discrete graph problems are well-posed in the continuum with $O(n\varepsilon_{n})$ labeled data points, which is the number of points in $\mathcal{X}_{n}\cap\partial_{\varepsilon_{n}}\Omega$, provided $\mathcal{X}_{n}$ are roughly evenly spread (e.g., an i.i.d. sequence with Lebesgue density). This is a labeling rate of $O(\varepsilon_{n})$, which vanishes as $n\to\infty$. However, it may vanish very slowly, since pointwise consistency requires lower bounds on $\varepsilon_{n}$, as explained in Remark 2.22 (see Eq. (2.47)). We emphasize that this is a general result that is independent of the structure of the learning algorithm or of the continuum equation. For the game-theoretic $p$-Laplacian with $p>d$, it was shown in [8] that the $\varepsilon$-ball game-theoretic $p$-Laplacian is well-posed in the continuum with $O(1)$ labels, which is a labeling rate of $O(1/n)$. This is far smaller than $O(\varepsilon_{n})$ due to (2.47). A similar result was proved for the variational $p$-Laplacian in [50], though here there is a restriction on $\varepsilon_{n}$ for well-posedness, even when $p>d$. In [16], it is shown that when $p=2$, Laplacian regularization is well-posed at the lower label rate of $O(\varepsilon_{n}^{2})$ using random walk techniques. For the $p$-Laplacian with $2\leq p\leq d$, it is an open problem to determine the lowest labeling rate at which the algorithm has a well-posed continuum limit. One would conjecture the lowest rate is $O(\varepsilon_{n}^{p})$.

3. Algorithms for p-Laplacian learning

We now present algorithms for $p$-Laplacian learning with both the game-theoretic and variational $p$-Laplacians. Section 3.1 reviews how to apply Newton’s method to the variational $p$-Laplace equation, and discuss how to apply homotopy on $p$ to accelerate convergence. Section 3.2 presents three algorithms for solving the game-theoretic $p$-Laplacian on a graph: a gradient descent approach, a Newton-like method, and a semi-implicit algorithm.

3.1. Newton’s Method for variational p-Laplacian

Since $J_{p}$ is smooth and convex, it is natural to use Newton’s method to minimize $J_{p}$. We give here the explicit details of the Newton iteration for minimizing $J_{p}$. It is useful to first rewrite the function $J_{p}(u)$ using vector notation. Let $\mathcal{X}=\{x^{1},\dots,x^{n+m}\}\subset\mathbb{R}^{d}$, where $\mathcal{O}=\{{x}^{n+1},\dots,{x}^{n+m}\}$ is the observation set. We define $u_{i}=u(x^{i})$ and set ${\bf u}=(u_{1},\dots,u_{n})\in\mathbb{R}^{n}$. Similarly, set $w_{ij}=w_{x_{i}x_{j}}$, $f_{i}=f(x^{i})$, $g_{i}=g(x^{i+n})$, ${\bf f}=(f_{1},\dots,f_{n})\in\mathbb{R}^{n}$, and ${\bf g}=(g_{1},\dots,g_{m})\in\mathbb{R}^{m}$. Then, subject to the constraints in (1.13), we can write

[TABLE]

Newton’s method corresponds to the iteration

[TABLE]

For notational convenience, define $a_{ij}({\bf u})=w_{ij}|u_{i}-u_{j}|^{p-2}$ and $b_{ij}({\bf u})=w_{i,j+n}|u_{i}-g_{j}|^{p-2}$, and

[TABLE]

Then, we write A({\bf u})=\big{(}a_{ij}({\bf u})\big{)}_{ij}\in\mathbb{R}^{n\times n},B({\bf u})=\big{(}b_{ij}({\bf u})\big{)}_{ij}\in\mathbb{R}^{n\times m} and $D({\bf u})=\text{diag}(d_{i}({\bf u}))\in\mathbb{R}^{n\times n}$. In this notation, we compute

[TABLE]

where $L({\bf u}):=D({\bf u})-A({\bf u})$, and thus the Newton update is given by

[TABLE]

The inversion of $L({\bf u}^{k})$ is performed with an iterative method, such as the preconditioned conjugate gradient. The matrix $L({\bf u}^{k})$, being a graph-Laplacian, is always positive semi-definite. If the graph is connected and ${\bf u}$ is nondegenerate, then it is also non-singular.

Proposition 3.1**.**

If the graph on $n$ nodes with weights $(a_{ij}({\bf u}))_{i,j=1}^{n}$ is connected, and $b_{ij}({\bf u})>0$ for some $i,j$, then $L({\bf u})$ is positive definite.

Proof.

Since $L({\bf u})$ is positive semi-definite, we simply have to prove that $L({\bf u})$ is non-singular. The proof follows a maximum principle argument. If $L({\bf u})\mathbf{x}=0$, $\mathbf{x}=(x_{1},\dots,x_{n})$, then we have

[TABLE]

for all $i=1,\dots,n$.

Let $i$ be an index for which $x_{i}\geq x_{j}$ for all $j$; that is, the node where $\mathbf{x}$ attains its maximum over the graph. We have two cases now.

Case 1. If $\sum_{j=1}^{m}b_{ij}({\bf u})=0$, then it follows from (3.5), and the fact that $x_{i}-x_{j}\geq 0$ for all $j$, that

[TABLE]

Thus, at any neighbor in the graph where $a_{ij}({\bf u})>0$ we have $x_{i}=x_{j}$.

Case 2. If $\sum_{j=1}^{m}b_{ij}({\bf u})>0$, then it follows from (3.5) that

[TABLE]

and so $x_{i}\leq 0$.

Note the observations above hold for any node $i$ where $\mathbf{x}$ attains its maximum value. We claim that these observations imply that $\max_{1\leq j\leq n}x_{j}\leq 0$. To see this, choose an index $i$ for which $\mathbf{x}$ attains it maximum value, and let $k\in\{1,\dots,n\}$ be a node for which

[TABLE]

which is guaranteed to exist by assumption. Since the graph with weights $(a_{ij}({\bf u}))_{i,j=1}^{n}$ is assumed connected, we can construct a path $i=i_{0},i_{1},\dots,i_{q}=k$ between node $i$ and $k$ for which $a_{i_{\ell},i_{\ell+1}}({\bf u})>0$ for $\ell=0,\dots,q-1$. We can assume that

[TABLE]

or else we can redefine $k$ as the earliest node along the path for which (3.6) holds. Thus, we can apply Case 1 above along the path to show that $x_{i_{\ell}}=x_{i_{\ell+1}}$ for $\ell=0,\dots,q-1$. Thus, $x_{i}=x_{k}$. Then we apply Case 2 above, since (3.6) holds, to show that $x_{i}=x_{k}\leq 0$, which proves the claim.

We have proved that if $L({\bf u})\mathbf{x}=0$ then $x_{i}\leq 0$ for all $i$. Since the equation is linear, we also have $L({\bf u})(-\mathbf{x})=0$ and so $-x_{i}\leq 0$ for all $i$. Thus, $L({\bf u})\mathbf{x}=0$ implies $\mathbf{x}\equiv 0$, and so $L({\bf u})$ is nonsingular. ∎

Remark 3.2*.*

The conditions in Proposition 3.1 hold when the original graph $G=(\mathcal{X},\mathcal{W})$ is connected, and $u_{i}\neq u_{j}$ for all $i\neq j$, and $u_{i}\neq g_{j}$ for some $i,j$ with $w_{i,j+n}>0$. In other words, if all the values $\{u_{1},\dots,u_{n},g_{1},\dots,g_{m}\}$ are unique then $L({\bf u})$. We can get into trouble when ${\bf u}$ is constant, or locally constant, since the conditions in the Proposition 3.1 fail to hold. In practice we find that $L({\bf u})$ remains non-singular throughout the Newton iterations.

Remark 3.3*.*

By the Newton-Kantorovich Theorem [46], Newton’s method is guaranteed to converge provided the initial guess is sufficiently close to the true solution. Since $J_{p}$ is convex, but not strongly convex for $p>2$, convergence may not be quadratic. In fact, according to [27], the number of iterations required for Newton’s method to converge may scale at least linearly in some cases.

Remark 3.4* (Homotopy on $p$).*

Let us mention that one can significantly speed up Newton’s method with a good starting point ${\bf u^{0}}$, since most of the computational cost of Newton’s method comes from approaching the solution, and once we reach the quadratic convergence region, only a few steps are required for convergence. As a result, we can solve (1.13) efficiently with a homotopy method in $p$. That is, we compute the solution for an increasing sequence of values of $p$, starting at $p=2$, and initializing Newton’s method each time from the solution from the previous value of $p$. If the steps in $p$ are small enough, this initialization falls in the quadratic convergence region and only a handful of Newton iterations are required for each increment in $p$. This approach is known in the literature as homotopy, and has been applied to a variety of problems in order to improve performance on related problems. For example, see [53, 26]. We illustrate the effectiveness of Newton’s method with homotopy in Section 4.2.

3.2. Algorithms for the game theoretic formulation

We now consider algorithms for solving the game theoretic problem (1.16). We shall discuss a gradient descent-type method in Section 3.2.1, a Newton-like algorithm in Section 3.2.2, and a semi-implicit method in Section 3.2.3.

3.2.1. Gradient Descent Approach

We first consider a gradient descent-type approach to solving (1.16), which is based on iterating

[TABLE]

where $\alpha>0$ is the time step. We call this a gradient descent-type iteration, where $\mathscr{L}^{G}_{p}u^{k}(x)+f(x)$ plays the role of a gradient, but it is important to point out it is not gradient descent, since the game theoretic $p$-Laplacian does not arise from a variational principle. The most straightforward stopping condition is to fix $\varepsilon>0$ and iterate until

[TABLE]

for all $x\in\mathcal{X}\setminus\mathcal{O}$. However, this stopping condition does not guarantee that $u^{k}$ is within $\varepsilon$ of the solution $u$ of the $p$-Laplacian learning problem (1.16), since the stability of the operator $\mathscr{L}^{G}_{p}$ depends in complicated ways on the graph and choice of boundary nodes $\mathcal{O}$. We present here a modification of the gradient descent method that uses the comparison principle to inform the stopping condition.

Lemma 3.5** (Comparison principle).**

Assume $w_{xy}\leq 1$ for all $x,y\in\mathcal{X}$. Let $u^{k},v^{k}:\mathcal{X}\to\mathbb{R}$ satisfy

[TABLE]

and

[TABLE]

for all $x\in\mathcal{X}\setminus\mathcal{O}$ and $0\leq k\leq T-1$, where $T\in\mathbb{N}$, and assume that $u^{0}\leq v^{0}$ and $u^{k}(x)\leq v^{k}(x)$ for all $x\in\mathcal{O}$ and $1\leq k\leq T$. If $\alpha\leq p/(2p-3)$ then $u^{k}\leq v^{k}$ on $\mathcal{X}$ for all $0\leq k\leq T$.

Proof.

Fix $u,v:\mathcal{X}\to\mathbb{R}$ and assume $u(x)\leq v(x)$ for all $x\in\mathcal{X}$. Fix $x\in\mathcal{X}$ and let $y_{1},y_{2}\in\mathcal{X}$ such that

[TABLE]

and

[TABLE]

Since $w_{xy}\leq 1$ and $\alpha\leq p/(2p-3)$ we have that

[TABLE]

We can now compute

[TABLE]

Thus, we have shown that when $\alpha\leq p/(2p-3)$ we have

[TABLE]

The proof is completed by using (3.8) and (3.9) to write

[TABLE]

for $x\in\mathcal{X}\setminus\mathcal{O}$ and using induction. ∎

We present our method in the case that $f(x)=0$. We define

[TABLE]

and

[TABLE]

and define $\overline{u}^{k}(x)$ and $\underline{u}_{k}(x)$ for $k\geq 1$ by

[TABLE]

for $x\in\mathcal{X}\setminus\mathcal{O}$, and $\overline{u}^{k+1}(x)=\underline{u}^{k+1}(x)=g(x)$ for $x\in\mathcal{O}$.

Remark 3.6*.*

The extension of this method to $f\neq 0$ is not immediately obvious. It is important for the convergence analysis in Theorem 3.7 below that $\overline{u}^{0}$ and $\underline{u}^{0}$ are super- and subsolutions of (1.16), respectively. When $f$ is nonzero, it is not clear how to construct such super- and subsolutions to initialize the method. We leave the extension of the method to nonzero $f$ to future work. We recall that in semi-supervised learning we always take $f=0$, so this extension is not needed for machine learning applications (though it may be of interest in numerical analysis).

We prove below that this iteration scheme converges to the solution $u_{*}:\mathcal{X}\to\mathbb{R}$ of (1.16) with $f=0$. We remark that the solution $u_{*}$ of (1.16) is unique when the graph is connected [8, 10].

Theorem 3.7** (Convergence).**

Assume $w_{xy}\leq 1$ for all $x,y\in\mathcal{X}$, and assume the graph is connected. Let $u_{*}$ be the solution of (1.16) with $f=0$, and define $u^{k}(x)=\tfrac{1}{2}(\overline{u}^{k}(x)+\underline{u}_{k}(x))$. If $\alpha\leq p/(2p-3)$ then for all $k\geq 1$ and all $x\in\mathcal{X}$ we have

[TABLE]

and

[TABLE]

Remark 3.8*.*

Theorem 3.7 proves convergence of the gradient descent-type scheme, and gives us a simple way to set the stopping condition. If we fix $\varepsilon>0$ and iterate until

[TABLE]

then Theorem 3.7 guarantees that $u^{k}(x)=\tfrac{1}{2}(\overline{u}^{k}(x)+\underline{u}_{k}(x))$ satisfies $|u^{k}(x)-u_{*}(x)|\leq\varepsilon$.

Proof.

By Lemma 3.5 we have

[TABLE]

for all $k\geq 0$ and $x\in\mathcal{X}$. It follows that

[TABLE]

and

[TABLE]

This establishes (3.13).

We now prove (3.14). We claim that $\mathscr{L}^{G}_{p}\overline{u}^{k}(x)\leq 0$ for all $k\geq 0$ and $x\in\mathcal{X}\setminus\mathcal{O}$. The proof is by induction. Since $\overline{u}^{0}(y)-\overline{u}^{0}(x)\leq 0$ for all $x\in\mathcal{X}\setminus\mathcal{O}$ and all $y\in\mathcal{X}$, the case $k=0$ is trivial. Now assume $\mathscr{L}^{G}_{p}\overline{u}^{k}(x)\leq 0$ for all $x\in\mathcal{X}\setminus\mathcal{O}$. Fix $x_{0}\in\mathcal{X}\setminus\mathcal{O}$ and define

[TABLE]

Then $w\leq\overline{u}^{k}$ and so

[TABLE]

as in the proof of Lemma 3.5. Since $w(x_{0})=\overline{u}^{k+1}(x_{0})$ and $w\geq\overline{u}^{k+1}$, we have

[TABLE]

which establishes the claim.

By a similar argument, we have that $\mathscr{L}^{G}_{p}\underline{u}^{k}(x)\geq 0$ for all $k\geq 0$ and $x\in\mathcal{X}\setminus\mathcal{O}$. It follows that $\overline{u}^{k+1}\leq\overline{u}^{k}$ and $\underline{u}^{k+1}\geq\underline{u}^{k}$ for all $k\geq 0$. Therefore, there exists $\overline{u},\underline{u}:\mathcal{X}\to\mathbb{R}$ such that

[TABLE]

By continuity of $\mathscr{L}^{G}_{p}$ we have that $\overline{u}$ and $\underline{u}$ both satisfy (1.16). Since the graph is connected, solutions of (1.16) are unique, and hence $\overline{u}=\underline{u}=u_{*}$, which completes the proof. ∎

If we make a minor modification to the equation, then we can obtain a linear convergence rate. Let $\varepsilon>0$ and consider the iteration

[TABLE]

for $x\in\mathcal{X}\setminus\mathcal{O}$, and $u^{k+1}(x)=g(x)$ for $x\in\mathcal{O}$. We prove below that this iteration scheme converges at a linear rate to the solution $u_{\varepsilon}:\mathcal{X}\to\mathbb{R}$ of

[TABLE]

When $\varepsilon>0$, the solution $u_{\varepsilon}$ of (3.17) is unique even when the graph is disconnected. We have the following convergence theorem, which is an adaptation of the contraction argument from [44].

Theorem 3.9** (Convergence rate).**

Assume $w_{xy}\leq 1$ for all $x,y\in\mathcal{X}$ and $\alpha\leq p/((2+\varepsilon)p-3)$. Let $u^{0}:\mathcal{X}\to\mathbb{R}$ and define $u_{\varepsilon}^{k}$ by (3.16), and let $u_{\varepsilon}$ be the solution of (3.17). Then for every $k\geq 0$ we have

[TABLE]

Remark 3.10*.*

Letting $u_{*}$ solve (1.16) and $u_{\varepsilon}$ solve (3.17), it follows from Theorem 3.9 that

[TABLE]

If we can quantify the error $|u_{\varepsilon}-u_{*}|$, then this would prove a convergence rate for the original problem (1.16) with $\varepsilon=0$. However, it seems that proving error estimates between (1.16) and (3.17) would require some additional strong assumptions on properties of the graph.

Proof.

Define

[TABLE]

Then we have $u^{k}_{\varepsilon}=\Phi^{k}_{\varepsilon}[u^{0}]$. As in the proof of Lemma 3.5 we have

[TABLE]

We also have

[TABLE]

for any constant $C>0$. It follows that

[TABLE]

Since

[TABLE]

we have

[TABLE]

provided

[TABLE]

which is equivalent to

[TABLE]

It follows that $\Phi_{\varepsilon}$ is a contraction in the norm $\|u\|_{\infty}:=\max_{x\in\mathcal{X}}|u(x)|$, and so there exists a unique fixed point $u_{\varepsilon}:\mathcal{X}\to\mathbb{R}$ such that $\Phi_{\varepsilon}[u_{\varepsilon}]=u_{\varepsilon}$ and

[TABLE]

for any $u:\mathcal{X}\to\mathbb{R}$. This completes the proof. ∎

3.2.2. A Newton-like Algorithm

We now consider a Newton-like method for solving (1.16). Newton’s method is based on iteratively solving a linearized version of the problem. In order to linearize (1.16), we define $y^{k}_{+}$ and $y^{k}_{-}$ by

[TABLE]

and define $\beta^{k}_{xy}$ by

[TABLE]

We also define

[TABLE]

The Newton-like iteration computes $u^{k+1}$ as the solution of

[TABLE]

We note that if $y^{k+1}_{\pm}=y^{k}_{\pm}(x)$ for all $x$, then $f=-\mathscr{L}^{G}_{p,k}u^{k+1}=-\mathscr{L}^{G}_{p}u^{k+1}$. Hence, we obtain convergence to the exact solution as soon as the locations of the min and max in (3.20) are correct. Thus, the algorithm can converge to the exact solution in a finite number of iterations. It seems rather difficult to construct a convergence proof for the Newton-like iterations, and we leave this to future work. We note that $\mathscr{L}^{G}_{p,k}$ is a linear operator, but may not be symmetric, which is one drawback of the method.

3.2.3. Semi-implicit Approach

Here, we extend the semi-implicit method of Oberman [45] to the graph setting. Given $\theta(x)\geq 1$, we add $-\theta(x)\Delta^{G}_{2}u(x)/(2d_{x})$ to both sides of (1.16), to obtain

[TABLE]

Solving for $\Delta^{G}_{2}u(x)$ on the left hand side reduces this equation to

[TABLE]

which suggests the iterative scheme

[TABLE]

where we have defined

[TABLE]

One fundamental advantage of the semi-implicit iteration (3.24) is that, unlike (3.23), the iteration (3.24) requires the solution of the same symmetric positive definite system at each iteration. This means we can dramatically speed up the solver, for large scale problems, by pre-computing a Cholesky factorization and using it at each iteration, or pre-computing an Incomplete Cholesky factorization, to be used as a preconditioner for CG. Another favorable feature of this scheme is the fact that the conditioning of the system is independent of $p$, allowing the scheme to reliably solve problems for any $p$.

The choice of $\theta(x)$ affects stability and convergence of the scheme. We give a heuristic argument here suggesting the iteration (3.24) is a contraction when

[TABLE]

Define $y^{k}_{\pm}(x)$ as in (3.20) and write $y_{\pm}$ in place of $y^{k}_{\pm}(x)$. Noting that (3.25) implies $-1\leq 2\gamma(x)-1\leq 1$, we have

[TABLE]

Therefore, for $f\equiv 0$ it follows from the definition of the iteration (3.24) that

[TABLE]

where $0<\beta(x)<1$. While this is not a contraction, it is suggestive of what we observe in practice, namely that the semi-implicit iteration is a contraction when (3.25) is satisfied. It seems that a convergence proof for the semi-implicit method not straightforward. Indeed, even in the case of a uniform grid in $2$-dimensions, a proof of convergence for the semi-implicit method is not available [45].

4. Algorithm Comparisons

We now give a numerical study of the performance of each algorithm on synthetic problems. This allows us, in particular, to control the intrinsic dimension of the graph—the dimension of the ambient Euclidean space or manifold from which the graph is sampled—and study how the intrinsic dimension affects convergence rates.

4.1. Experiment Design

This section specifies design choices for our synthetic data experiments.

4.1.1. Problem S (synthetic data)

We sample $\mathcal{X}$ from a uniform distribution on $[0,1]^{d}$, and label $10$ of these points with labels sampled from a uniform distribution on $[0,1]$. The experiments we report on use $d=2$, $d=5$ and $d=10$, as the computational cost did not change substantially for $d>10$.

4.1.2. Graph Construction

Define the relation $\sim_{K}$ on $\mathcal{X}$ by $x\sim_{K}y$ if $x$ is among the $k$ nearest neighbors of $y$ in Euclidean distance. We construct a symmetric $K$-nearest neighbor graph as follows

[TABLE]

The constant $\sigma$ plays the role of being a typical length scale for the problem. In our experiments, we compute it as \sigma=\frac{1}{2}\max\big{\{}|x-y|\,:\,w_{xy}>0\big{\}}. We use $K=10$ in all experiments. Other choices of $K$ would produce similar accuracy, but larger $K$ would increase the cost of each iteration. The only requirement is that $K$ needs to be large enough to ensure the graph is connected. It should be noted that, even though we prescribe $K=10$, after we apply (4.1), most vertices will have a higher number of neighbors (so that symmetry can be enforced).

4.1.3. Error Reporting

For the variational formulation, we report

[TABLE]

which amounts to scaling the largest residual in equation (1.13), by the number of points, and the typical length scale $\sigma$. This scaling ensures a fair comparison across different $n,p,d$. For the game theoretic formulation, the scaling of the residual in equation (1.16) is different, and we report

[TABLE]

4.2. Homotopy Results for Newton and Newton-like Methods

Homotopy is very effective at improving the performance and reliability of Newton’s method for higher $p$. We illustrate this by solving problem S with $d=10$, $p=50$, $m=10$ and $n=10^{4}$. We use the solution to the $p=2$ problem as ${\bf u}^{0}$ for the $p=3$ problem. Then, we apply homotopy until $p=50$ is reached. For each $p$, we require $\epsilon<10^{-12}$ according to (4.2). Table 1 reports $\epsilon$ at each iteration.

We repeat the experiment using the Newton-like algorithm for the game-theoretic formulation (now $\epsilon$ is given by (4.3)). As we can see from Table 2, homotopy is equally effective in this case.

4.3. Computational Cost vs. Dimension $d$ and Size $n$

Our timing results were performed on a laptop with 16GB of RAM memory. For the gradient descent method, we used a C-language implementation. Other methods were implemented in MATLAB 2018, as the computational cost was dominated by matrix operations. For these methods, we observed the cost of direct solvers increase by as much as $100\times$, when going from $d=2$ up to $d=100$, hence we used iterative methods instead. For Newton’s method, we used a conjugate gradient’s method (CG), preconditioned with an Incomplete Cholesky factorization. Given that the matrix system for Newton-like’s method is no longer symmetric, we used a Generalized Minimal Residual Algorithm (GMRES) [48], preconditioned with an Incomplete LU factorization instead. The semi-implicit method performed very well with a preconditioned CG method. Since the matrix to be inverted does not change, we compute an Incomplete Cholesky factorization, with a drop tolerance of $10^{-1}$, at the very beginning, and use it as our preconditioner in each iteration.

We explore how the Newton, Newton-like, Semi-Implicit, and Gradient Descent algorithms scale with $n,d$. To measure this relationship, we solved problem S for $p=11$, until $\epsilon<10^{-7}$, according to (4.2) and (4.3) respectively. These results are reported in Figure 2. In the figure, the complexity of Newton’s method scales roughly like $n^{0.8}$, the Newton-like method scales roughly like $n^{1.3}$, the semi-implicit method scales like $n^{1}$ for $d=2$ and $d^{0.75}$ for $d=5,10$, and the gradient descent-type method scales like $n^{2}$ for $d=2$, and $n^{1.5}$ for $d=5,10$. We note that Figure 2a and Figures 2(b-d) cannot be directly compared, since the parameters $p,\epsilon$ do not have the same meaning in both problems, but give a rough idea of the relative performance of the variational and game-theoretic solvers. We also illustrate our ability to solve large scale problems by solving problem S, with $m=10$, $K=10$, $n=5\cdot 10^{5}$ and $d=10$. For the variational formulation, we use Newton’s method starting from the solution of the $p=2$ problem ($1$-step homotopy). For the game theoretic formulation, we choose the semi-implicit method. In both cases, iterative solvers were used. We plot the results, averaged over $5$ trials in Figure 3. In the figure, the CPU time for Newton’s method scales like $\varepsilon^{-0.05}$ with tolerance $\varepsilon>0$, while the CPU time for the semi-implicit method scales like $\varepsilon^{-0.2}$.

Remark 4.1*.*

It is interesting to point out that the results in Figure 2 show that the cost of all four methods decreases (sometimes significantly) as the dimension $d$ increases. This is due to the graph Laplacian matrices having better condition numbers in higher dimensions. To see why this is the case, we recall that on a uniform grid, the condition number for a discrete Laplacian is on the order of $\Delta x^{-2}$, where $\Delta x$ is the grid resolution. When $\Delta x$ is small, the Laplacian is poorly conditioned and iterative methods are stiff without some form of preconditioning. The same scaling for the condition number holds for graph Laplacians constructed from random geometric graphs, due to the spectral convergence results in, for example, [14], but now $\Delta x$ should be replaced by the length scale $\varepsilon>0$ on which the graph is constructed. In dimension $d$, the length scale is at least $\varepsilon\geq n^{-1/d}$, to ensure graph connectivity. This implies that as a function of $n$ and $d$, the condition number of the graph Laplacian should scale like $n^{-2/d}$. Thus, the curse of dimensionality improves the condition number of graph Laplacians in higher dimensions, and explains why iterative methods perform better for larger $d$. The stiffness of Laplace equations is essentially a phenomenon of high resolution meshes, which are only possible to construct in low dimensional spaces (e.g., $d=2,3$). By the manifold assumption in machine learning [18], we expect our graphs to have intrinsic dimensionality significantly higher than $d=3$, and so the graph Laplacians we encounter in practice are better conditioned than what we may expect from experience with solving PDEs in $d=2$ or $d=3$ dimensions.

4.4. Synthetic Experiment Conclusions

The experiments shown demonstrate that both variational and game-theoretic formulations can be solved with the algorithms we propose for large scale graphs with high intrinsic dimensionality, similar to what will be observed in practice. For the game-theoretic problem, the semi-implicit formulation is likely to perform the best, although it lacks the convergence guarantees enjoyed by the gradient descent method. The Newton and Newton-like algorithms show strong resemblance to one another, and both methods improve their performance and reliability substantially when homotopy is utilized.

5. Experiments with Real Data

We now give an experimental study of $p$-Laplacian semi-supervised learning classification on real datasets, including MNIST, and the more complex datasets Fashion MNIST and Extended MNIST.

5.1. Description of the Experiments and Datasets

On each dataset, we shall perform two experiments. For the first, we explore the performance as the number of labels varies, with the goal of understanding how different models perform with increasing, but very small, amounts of labeled data. For the second experiment, we fix the number of labels to $1$ per class, and explore the performance as a function of the amount of unlabeled data. This experiment helps to illustrate the degeneracy of the $2$-Laplacian, while also illustrating how other models can profit from an increasing amount of unlabeled data, which is the premise upon which semi-supervised learning is built. Before proceeding to the results, we provide some detail on the datasets we work with.

In all three problems, we first preprocess each image via a Scattering Transform [7] in order to extract features upon which to build the graph. The graphs are constructed as $K$-nearest neighbor graphs in the feature space using the weights (4.1). In all experiments we chose $K=10$, to ensure graph connectivity and sufficient sparsity. We also experimented with $K=25$ and $K=50$ and found that increasing the number of neighbors in the graph does not provide any advantages, and results in significantly higher computational cost.

We solve the $M$-class classification problem with the standard one-vs-rest approach, whereby we solve $M$ binary classification problems classifying each digit against the rest. This gives $M$ scores for each unlabeled image, and the final classification is chosen as the class with the highest score.

5.1.1. MNIST

The MNIST dataset is a standard benchmark that consists of $70,000$ images of handwritten digits [math] through $9$ [37]. Each image is a $28\times 28$ pixel grayscale image, meaning it can be represented as a vector in $\mathbb{R}^{784}$ dimensions. The classes are well balanced, with roughly $7,000$ examples for each digit.

5.1.2. Fashion MNIST

This dataset, recently introduced by [55], was designed as a drop-in replacement for MNIST, in order to test classification accuracy for a significantly harder problem. This dataset also consists of $70,000$ grayscale images of size $28\times 28$ pixels, except the $10$ handwritten digits are replaced by $10$ classes of clothing items, such as sandals, and dresses. For an example, see Figure 4.

5.1.3. Extended MNIST

The Extended MNIST (EMNIST) dataset [19] contains images of handwritten letters, of $28\times 28$ pixels. Given that letters may be uppercase or lowercase, and there are now $26$ classes, the semisupervised task on this dataset is substantially more challenging. To simplify the task slightly, we downsample the dataset for the first experiment (performance vs. number of labels) so that all classes have the same number of labels. We end up with $3419$ samples for each of the $26$ classes, for a total of $n=88,894$ images. For the second task, we measure performance with subsampled datasets of $2^{5},2^{6},\dots,2^{11}$ samples per class, always having $1$ labeled sample.

5.2. Experimental Results

We report the results of our experiments in Figure 5. In all experiments, we report results using $5$ formulations: the standard $2$-Laplacian, the Weighted Non-local Laplacian (WNLL) [49], and the game-theoretic $p$-Laplacian for $p=5$, $p=9$ and $p=\infty$ (e.g. Lipschitz learning). All models were solved to a tolerance of $\epsilon<10^{-2}$, which provides consistent results. The results we report paint a similar picture across all three datasets. For the first type of experiment (left column), all formulations improve their accuracy as the number of labels increases, and they gradually approach a similar level as $m$ grows. This means that the choice of a good model is of particular importance in the regime when labeled data is extremely limited. Both the $p$-Laplace and WNLL methods clearly outperform the $2$-Laplacian at low label rates, while the WNLL and $p$-Laplacian give fairly similar results.

The second experiment (right column) provides further evidence of the superiority of the $p$-Laplacian model over the $2$-Laplacian. In this case we fix $1$ label per class, which means $m=10$ for MNIST and Fashion MNIST, while $m=26$ for EMNIST. The premise of semi-supervised is that we may achieve superior performance by including both labeled and unlabeled data. If we keep $m$ fixed, we should expect the accuracy to improve as the number of unlabeled images increases. In all three datasets, the $2$-Laplacian classification performance decreases substantially as the amount of unlabeled data $n$ grows. This illustrates the need for alternative models when $n\gg m$. On the other hand, the $p$-Laplacian models (for sufficiently large $p$) do not degenerate as $n$ grows, and in fact, their performance increases slightly as $n$ grows. The Weighted Non-local Laplacian exhibits mixed results. For the MNIST and Fashion MNIST datasets, performance decreases slightly when $n=70,000$. Meanwhile, this model outperforms others on the EMNIST dataset.

Let us make some final remarks comparing and contrasting WNLL and $p$-Laplace learning. The experimental results in Figure 5 show that neither method is strictly better than the other. For the most part, the methods are comparable. The largest difference is seen at 1 label per class on MNIST where $p$-Laplace is roughly 10% better than WNLL. The difference is much more pronounced in terms of the theoretical guarantees for each method. It was shown in [8] (see [50] for the variational $p$-Laplacian) that $p$-Laplace learning is well-posed with arbitrarily few labeled examples; in fact, we can take the number of labeled examples to be finite while sending the number of unlabeled examples to infinity, and still obtain a well-posed continuous extension of the label values. On the other hand, the WNLL was shown in [15] to be ill-posed in the same setting of finite labeled data and infinite unlabeled data. In particular, [15, Corollary 3.8] shows that the WNLL with finite labeled data converges to a constant labeling function as the amount of unlabeled data tends to infinity. To the best of our knowledge, the only semi-supervised learning algorithm with theoretical guarantees at arbitrarily low label rates, such as the ones given in [8, 50], is the graph $p$-Laplacian.

6. Conclusions

This paper is focused on theory and applications of $\ell_{p}$-Laplacian regularized semi-supervised learning. We preformed a detailed analysis of discrete to continuum theory for the $p$-Laplacian on $k$-NN graphs, which are more commonly used in practice compared to random geometric graphs, and made the surprising discovery that the $p$-Laplacian models retain information about the data distribution as $p\to\infty$ on $k$-NN graphs, contrary to conventional wisdom from the existing $p$-Laplace theory on random geometric graphs. We also presented a simple and very general framework for proving discrete to continuum convergence results that only requires pointwise consistency and a monotonicity property. We expect this framework to be useful in future work.

We also studied and developed algorithms for solving the variational and game-theoretic formulations of the $p$-Laplacian on a weighted graph. The variational formulation may be solved efficiently using Newton’s method with homotopy. The semi-implicit method is the fastest method for solving the game-theoretic formulation, while the gradient-descent approach enjoys rigorous convergence guarantees. Our experiments with real data show that $p$-Laplacian learning is superior to Laplace learning ($p=2$) at very low label rates on common image classification datasets including MNIST, FashionMNIST and EMNIST.

Appendix A Review of viscosity solutions

Viscosity solutions are a notion of weak solution for partial differential equations that obeys the maximum principle and enjoys strong stability and uniqueness theorems. The theory is especially useful for passing from discrete to continuum limits (see, e.g., [9, 17, 13]). We review here the basic definitions. Let $\text{USC}(\overline{\Omega})$ (resp. $\text{LSC}(\overline{\Omega})$) denote the collection of functions that are upper (resp. lower) semicontinuous at all points in $\overline{\Omega}$. We make the following definitions.

Definition A.1**.**

We say $u\in\text{USC}(\overline{\Omega})$ is a viscosity subsolution of (2.48) if for all $x\in\overline{\Omega}$ and every $\varphi\in C^{\infty}(\mathbb{R}^{n})$ such that $u-\varphi$ has a local maximum at $x$ with respect to $\overline{\Omega}$

[TABLE]

Likewise, we say that $u\in\text{LSC}(\overline{\Omega})$ is a viscosity supersolution of (2.48) if for all $x\in\overline{\Omega}$ and every $\varphi\in C^{\infty}(\mathbb{R}^{n})$ such that $u-\varphi$ has a local minimum at $x$ with respect to $\overline{\Omega}$

[TABLE]

Finally, we say that $u$ is a viscosity solution of (2.48) if $u$ is both a viscosity sub- and supersolution. In this case, we say that the boundary conditions in (2.48) hold in the viscosity sense, which is also known as the generalized Dirichlet sense.

Definition A.2**.**

We say that (2.48) enjoys strong uniqueness if whenever $u\in\text{USC}(\overline{\Omega})$ is a subsolution of (2.48) and $v\in\text{LSC}(\overline{\Omega})$ is a supersolution, we have $u\leq v$ on $\overline{\Omega}$.

We refer the reader to [20] for the proof of the comparison principle for viscosity solutions with generalized Dirichlet boundary conditions, which implies strong uniqueness.