Semisupervised regression in latent structure networks on unknown manifolds

We shall denote a vector by a bold lower-case letter, $\mathbf{x}$ for instance. Bold upper-case letters such as $\mathbf{P}$ will be used to represent matrices. The Frobenius norm and the maximum row-norm of a matrix $\mathbf{B}$ will be denoted respectively by $\left\lVert\mathbf{B}\right\rVert_{F}$ and $\left\lVert\mathbf{B}\right\rVert_{2,\infty}$. The $i$-th row of a matrix $\mathbf{B}$ will be denoted by $\mathbf{B}_{i*}$, and the $j$-th column of $\mathbf{B}$ will be denoted by $\mathbf{B}_{*j}$. We will denote the $n\times n$ identity matrix by $\mathbf{I}_{n}$, and $\mathbf{J}_{n}$ will denote the $n\times n$ matrix whose each entry equals one. Also, $\mathbf{H}_{n}=\mathbf{I}_{n}-\frac{1}{n}\mathbf{J}_{n}$ will denote the $n\times n$ centering matrix. Unless otherwise mentioned, $\left\lVert\mathbf{x}\right\rVert$ will represent the Euclidean norm of a vector $\mathbf{x}$. The set of all orthogonal $d\times d$ matrices will be denoted by $\mathcal{O}(d)$. The set of all positive integers will be denoted by $\mathbb{N}$ and for any $n\in\mathbb{N}$, $[n]$ will denote the set $\left\{1,2,3...,n\right\}$.

A graph is an ordered pair $(V,E)$ where $V$ is the set of vertices (or nodes) and $E\subset V\times V$ is the set of edges connecting vertices. An adjacency matrix $\mathbf{A}$ of a graph is defined as $\mathbf{A}_{ij}=1$ if $(i,j)\in E$, and $\mathbf{A}_{ij}=0$ otherwise. Here, we deal with hollow and undirected graphs; hence $\mathbf{A}$ is symmetric and $\mathbf{A}_{ii}=0$ for all $i$. Latent position random graphs are those for which each node is associated with a vector that is called its latent position, denoted by $\mathbf{x}_{i}$, and the probability of formation of an edge between the $i$-th and $j$-th nodes is given by $\kappa(\mathbf{x}_{i},\mathbf{x}_{j})$ where $\kappa$ is a suitable kernel.

Random vectors drawn from any arbitrary probability distribution cannot be latent positions of a random dot product graph, as their magnitudes can be unbounded whereas probabilities must lie in the interval $[0,1]$. The following definition allows us to work with a restricted class of distributions more amenable to random dot product graphs.

Next, we define the random dot product graphs, the basis of the models considered in this paper.

The latent positions of a random dot product graph are typically unknown and need to be estimated in practice. The following definition puts forth an estimate of these latent positions via adjacency spectral embedding.

Now, we present two results from the literature which give us the consistency and asymptotic normality of suitably rotated adjacency spectral estimates of the latent positions of a random dot product graph.

Henceforth, we will denote this optimally rotated adjacency spectral embedding by $\tilde{\mathbf{X}}^{(n)}$, that is, $\tilde{\mathbf{X}}^{(n)}=\hat{\mathbf{X}}^{(n)}\mathbf{W}^{(n)}$. To simplify notations we will often omit the superscript $n$, and we will use $\tilde{\mathbf{x}}_{i}$ to denote the $i$-th row of $\tilde{\mathbf{X}}$. Observe that in the attempt to estimate the true latent positions from the adjacency matrix, we encounter an inherent non-identifiability issue: for any $\mathbf{W}\in\mathcal{O}(d)$, $\mathbb{E}(\mathbf{A}|\mathbf{X})=\mathbf{X}\mathbf{X}^{T}=(\mathbf{X}\mathbf{W})(\mathbf{X}\mathbf{W})^{T}$. This is the reason why the adjacency spectral embedding needs to be rotated suitably so that it can approximate the true latent position matrix.

Under suitable regularity conditions, a combined use of Theorem 2 and Delta method gives us the asymptotic distribution of $\sqrt{n}(\gamma(\tilde{\mathbf{x}}_{i})-\gamma(\mathbf{x}_{i}))$ for a function $\gamma:\mathbb{R}^{d}\to\mathbb{R}$, which depends on the true distribution $F$ of the latent positions. Therefore, $\mathrm{var}(\gamma(\tilde{\mathbf{x}}_{i})-\gamma(\mathbf{x}_{i}))$ can be approximated from the optimally rotated adjacency spectral estimates $\tilde{\mathbf{x}}_{i}$ and their empirical distribution function. In a random dot product graph for which the latent positions lie on a known one-dimensional manifold in ambient space $\mathbb{R}^{d}$, and the nodal responses are linked to the scalar pre-images of the latent positions via a simple linear regression model, we can use the approximated variance of $(\gamma(\tilde{\mathbf{x}}_{i})-\gamma(\mathbf{x}_{i}))$ (motivated by works in Chapters $2$, $3$ of fuller1987measurement ) to improve the performance (in terms of mean squared errors) of the naive estimators of the regression parameters, which are obtained by replacing $\mathbf{x}_{i}$ with $\tilde{\mathbf{x}}_{i}$ in the least square estimates of the regression parameters. We demonstrate this method in detail in Section 7 (see Figure 7).

After stating the relevant definitions and results pertinent to random dot product graphs, in the following section we introduce the manifold learning technique we will use in this paper. Just for the sake of clarity, the topic of manifold learning in general has nothing to do with random dot product graph model; hence Section 2.2 and Section 2.3 can be read independently.

Our main model is based on a random dot product graph whose latent positions lie on a one-dimensional Riemannian manifold. Since one-dimensional Riemannian manifolds are isometric to one-dimensional Euclidean space, we wish to represent the latent positions as points on the real line. This is the motivation behind the use of the following manifold learning technique, which relies upon approximation of geodesics by shortest path distances on localization graphs (tenenbaum2000global , bernstein2000graph , trosset2021rehabilitating ). Given points $\mathbf{x}_{1},\dots\mathbf{x}_{n}\in\mathcal{M}$ where $\mathcal{M}$ is an unknown one-dimensional manifold in ambient space $\mathbb{R}^{d}$, the goal is to find $\hat{z}_{1},\dots\hat{z}_{n}\in\mathbb{R}$, such that the interpoint Euclidean distances between $\hat{z}_{i}$ approximately equal the interpoint geodesic distances between $\mathbf{x}_{i}$. However, the interpoint geodesic distances between $\mathbf{x}_{i}$ are unknown. The following result shows how to estimate these unknown geodesic distances under suitable regularity assumptions.

Given the dissimilarity matrix $\mathbf{D}=\left(d_{n,\lambda}(\mathbf{x}_{i},\mathbf{x}_{j})\right)_{i,j=1}^{n}$, the raw-stress function at $(z_{1},\dots z_{n})$ is defined as

where $w_{ij}\geq 0$ are weights. For the purpose of learning the manifold $\mathcal{M}$, we set $w_{ij}=1$ for all $i,j$, and compute

Since the scalars $\hat{z}_{i}$ are obtained by embedding $\mathbf{D}$ into one-dimension upon minimization of raw-stress, we shall henceforth refer to $\hat{z}_{i}$ as the one-dimensional raw-stress embeddings of $\mathbf{D}$.

Suppose $\psi:\mathbb{R}\to\mathbb{R}^{d}$ is a known bijective function and let $\mathcal{M}=\psi([0,L])$ be a one-dimensional compact Reimannian manifold. Consider a random dot product graph for which the nodes (with latent positions $\mathbf{x}_{i}$) lie on the known one-dimensional manifold $\mathcal{M}$ in $d$-dimensional ambient space. Let $t_{1},\dots t_{n}$ be the scalar pre-images of the latent positions such that $\mathbf{x}_{i}=\psi(t_{i})$ for all $i\in[n]$, where $n$ is the number of nodes of the graph. Suppose, for each $i$, the $i$-th node is coupled with a response $y_{i}$ which is linked to the latent position via the following regression model

where $\epsilon_{i}\sim^{iid}N(0,\sigma^{2}_{\epsilon})$ for all $i\in[n]$. Our goal is to estimate $\alpha$ and $\beta$.
If the true regressors $t_{i}$ were known, we could estimate $\alpha$ and $\beta$ by their ordinary least square estimates given by

Since the true latent positions $\mathbf{x}_{i}$ are unknown, we estimate the true regressors $t_{i}$ by
$\hat{t}_{i}=\arg\min_{t}\left\lVert\tilde{\mathbf{x}}_{i}-\psi(t)\right\rVert$ where $\tilde{\mathbf{x}}_{i}$ is the optimally rotated adjacency spectral estimate for the $i$-th latent position $\mathbf{x}_{i}$. The existence of $\hat{t}_{i}$ is guaranteed by the compactness of the manifold $\mathcal{M}=\psi([0,L])$. We then substitute $t_{i}$ by $\hat{t}_{i}$ in $\hat{\alpha}_{true}$ and $\hat{\beta}_{true}$ to obtain the substitute (or the plug-in estimators) given by

The steps to compute $\hat{\alpha}_{sub}$ and $\hat{\beta}_{sub}$ are formally stated in Algorithm 1a.

Here, we assume that $\psi:[0,L]\to\mathbb{R}^{d}$ is unknown and arclength parameterized, that is, $\left\lVert\dot{\psi}(t)\right\rVert=1$ for all $t$. Additionally, assume that $\mathcal{M}=\psi([0,L])$ is a compact Reimannian manifold. Consider a $n$-node random dot product graph whose nodes (with latent positions $\mathbf{x}_{i}$) lie on the unknown manifold $\mathcal{M}=\psi([0,L])$ in ambient space $\mathbb{R}^{d}$. Assume that the first $s$ nodes of the graph are coupled with a response covariate and the response $y_{i}$ at the $i$-th node is linked to the latent position via a linear regression model

where $\epsilon_{i}\sim^{iid}N(0,\sigma^{2}_{\epsilon})$ for all $i\in[s]$. Our goal is to predict the response for the $r$-th node, where $r>s$. First, we compute the adjacency spectral estimates $\hat{\mathbf{x}}_{i}$ of the latent positions of all $n$ nodes. We then construct a localization graph on the adjacency spectral estimates $\hat{\mathbf{x}}_{i}$ under the following rule: join two nodes $\hat{\mathbf{x}}_{i}$ and $\hat{\mathbf{x}}_{j}$ if and only if $\left\lVert\hat{\mathbf{x}}_{i}-\hat{\mathbf{x}}_{j}\right\rVert<\lambda$, for some pre-determined $\lambda>0$ known as the neighbourhood parameter. Denoting the shortest path distance between $\hat{\mathbf{x}}_{i}$ and $\hat{\mathbf{x}}_{j}$ by $d_{n,\lambda}(\hat{\mathbf{x}}_{i},\hat{\mathbf{x}}_{j})$, we embed the dissimilarity matrix $\mathbf{D}=\left(d_{n,\lambda}(\hat{\mathbf{x}}_{i},\hat{\mathbf{x}}_{j})\right)_{i,j=1}^{l}$ into one-dimension by minimizing the raw-stress criterion, thus obtaining

where $l$ is such that $s<r<l\leq n$. We then use a simple linear regression model on the bivariate data $(y_{i},\hat{z}_{i})_{i=1}^{s}$ to predict the response for $\hat{z}_{r}$ corresponding to the $r$-th node. The abovementioned procedure to predict the response at the $r$-th node from given observations is formally described in Algorithm 1b.

Recall that we observe a random dot product graph with $n$ nodes for which the latent positions $\mathbf{x}_{i}$ lie on a one-dimensional manifold. Our following result shows that we can consistently estimate $(\alpha,\beta)$ by $(\hat{\alpha}_{sub},\hat{\beta}_{sub})$.

A rough sketch of proof for Theorem 4 is as follows. Note that $\hat{t}_{i}=\arg\min_{t}\left\lVert\tilde{\mathbf{x}}_{i}-\psi(t)\right\rVert$ where $\tilde{\mathbf{x}}_{i}$ is the optimally rotated adjacency spectral estimate of $\mathbf{x}_{i}$, and $t_{i}=\arg\min_{t}\left\lVert\mathbf{x}_{i}-\psi(t)\right\rVert$. Recall that Theorem 1 tells us that $\tilde{\mathbf{x}}_{i}$ is consistent for $\mathbf{x}_{i}$. This enables us to prove that $\hat{t}_{i}$ is consistent for $t_{i}$ which ultimately leads us to the consistency of $\hat{\alpha}_{sub}$ and $\hat{\beta}_{sub}$, because $\hat{\alpha}_{sub}$ and $\hat{\beta}_{sub}$ are computed by replacing the true regressors $t_{i}$ with $\hat{t}_{i}$ in the expressions of $\hat{\alpha}_{true}$ and $\hat{\beta}_{true}$ which are consistent for $\alpha$ and $\beta$ respectively. Theorem 4 gives us the consistency of the substitute estimators of the regression parameters under the assumption of boundedness of the gradient of the inverse function. Since continuously differentiable functions have bounded gradients on compact subsets of their domain, we can apply Theorem 4 whenever $\gamma=\psi^{-1}$ can be expressed as a restriction to a function with continuous partial derivatives. As a direct consequence of Theorem 4, our next result demonstrates that the substitute estimators have optimal asymptotic variance amongst all linear unbiased estimators.

Recall that our goal is to predict responses in a semisupervised setting on a random dot product graph on an unknown one-dimensional manifold in ambient space $\mathbb{R}^{d}$. We provide justification for the use of Algorithm 1b for this purpose, by showing that the predicted response $\tilde{y}_{r}$ at the $r$-th node based on the raw-stress embeddings approaches the predicted response based on the true regressors $t_{i}$ as $n\to\infty$.

Intuition suggests that in order to carry out inference on the regression model, we must learn the unknown manifold $\mathcal{M}$. We exploit the availability of large number of unlabeled nodes whose latent positions lie on the one-dimensional manifold, to learn the manifold. Observe that since the underlying manifold $\psi([0,L])$ is arclength parameterized, the geodesic distance between any two points on it is the same as the interpoint distance between their corresponding pre-images. Results from the literature (bernstein2000graph ) show that if an appropriate localization graph is constructed on sufficiently large number of points on a manifold, then the shortest path distance between two points approximates the geodesic distance between those two points. Therefore, on a localization graph of an appropriately chosen neighbourhood parameter $\lambda$, constructed on the adjacency spectral estimates $\hat{\mathbf{x}}_{1},\hat{\mathbf{x}}_{2},...,\hat{\mathbf{x}}_{n}$, the shortest path distance $d_{n,\lambda}(\hat{\mathbf{x}}_{i},\hat{\mathbf{x}}_{j})$ between $\hat{\mathbf{x}}_{i}$ and $\hat{\mathbf{x}}_{j}$ is expected to approximate the geodesic distance $d_{M}(\mathbf{x}_{i},\mathbf{x}_{j})=|t_{i}-t_{j}|$. Note that here is no need to rotate the adjacency spectral estimates $\hat{\mathbf{x}}_{i}$, because the shortest path distance is sum of Euclidean distances and Euclidean distances are invariant to orthogonal transformations. Thus, when the dissimilarity matrix $\mathbf{D}=\left(d_{n,\lambda}(\hat{\mathbf{x}}_{i},\hat{\mathbf{x}}_{j})\right)_{i,j=1}^{l}$ is embedded into one-dimension by minimization of raw-stress, we obtain embeddings $\hat{z}_{1},\dots\hat{z}_{l}$ for which interpoint distance $|\hat{z}_{i}-\hat{z}_{j}|$ approximates the interpoint distance $|t_{i}-t_{j}|$. In other words, the estimated distances from the raw-stress embeddings applied to the adjacency spectral estimates of the latent positions approximate the true geodesic distances, which is demonstrated by the following result. This argument is the basis for construction of a sequence of predicted responses based on the raw-stress embeddings, which approach the predicted responses based on the true regressors as the number of auxiliary nodes go to infinity.

Theorem 5 shows that the one-dimensional raw-stress embeddings $\hat{z}_{1},\dots\hat{z}_{l}$ satisfy $(|\hat{z}_{i}-\hat{z}_{j}|-|t_{i}-t_{j}|)\to 0$ as $K\to\infty$, for all $i,j\in[l]$. This means that for every $i\in[l]$, $\hat{z}_{i}$ approximates $t_{i}$ up to scale and location transformations. Since in simple linear regression the accuracy of the predicted response is independent of scale and location transformations, we can expect the predicted response at a particular node based on $\hat{z}_{i}$ to approach the predicted response based on the true regressors $t_{i}$. The following theorem, our key result in this paper, demonstrates that this is in fact the case.

Recall that the validity of a simple linear regression model can be tested by an $F$-test, whose test statistic is dependent on the predicted responses based on the true regressors. Since we have a way to approximate the predicted responses based on the true regressors by predicted responses based on the raw-stress embeddings, we can also devise a test whose power approximates the power of the $F$-test based on the true regressors, as shown by our following result.

Thus, if one wants to perform a test for model validity in the absence of the true regressors $t_{i}$, then a test of approximately equal power, based on the raw-stress embeddings $\hat{z}_{i}$ for a graph of sufficiently large number of auxiliary nodes, can be used instead.

We take the manifold to be $\psi([0,1])$, where $\psi:[0,1]\to\mathbb{R}^{3}$ is the Hardy Weinberg curve, given by $\psi(t)=(t^{2},2t(1-t),(1-t)^{2})$. The number of nodes, $n$, varies from $600$ to $2500$ in steps of $100$. For each $n$, we repeat the following procedure over $100$ Monte Carlo samples. A sample $t_{1},.....t_{n}\sim^{iid}U[0,1]$ is generated, and responses $y_{i}$ are sampled via the regression model $y_{i}=\alpha+\beta t_{i}+\epsilon_{i}$, $\epsilon_{i}\sim^{iid}N(0,\sigma^{2}_{\epsilon})$, $i\in[n]$, where $\alpha=2.0$, $\beta=5.0$, $\sigma_{\epsilon}=0.1$. An undirected hollow random dot product graph with latent positions $\mathbf{x}_{i}=\psi(t_{i})$, $i\in[n]$ is generated. More specifically, the $(i,j)$-th element of the adjacency matrix $\mathbf{A}$ satisfies $\mathbf{A}_{ij}\sim\mathrm{Bernoulli}(\mathbf{x}_{i}^{T}\mathbf{x}_{j})$ for all $i<j$, and $\mathbf{A}_{ij}=\mathbf{A}_{ji}$ for all $i,j\in[n]$, and $\mathbf{A}_{ii}=0$ for all $i\in[n]$. We denote the true latent position matrix by $\mathbf{X}=[\mathbf{x}_{1}|\mathbf{x}_{2}|...|\mathbf{x}_{n}]^{T}$, and the adjacency spectral estimate of it by $\hat{\mathbf{X}}$. We compute

and finally, set $\tilde{\mathbf{X}}=\hat{\mathbf{X}}\hat{\mathbf{W}}$ to be the optimally rotated adjacency spectral estimate of the latent position matrix $\mathbf{X}$. Then we obtain $\hat{t}_{i}=\arg\min_{t}\left\lVert\tilde{\mathbf{x}}_{i}-\psi(t)\right\rVert$ for $i\in[n]$, and get $\hat{\alpha}_{sub}$ and $\hat{\beta}_{sub}$. Setting $\boldsymbol{\theta}=(\alpha,\beta)$, $\hat{\boldsymbol{\theta}}_{true}=(\hat{\alpha}_{true},\hat{\beta}_{true})$ and $\hat{\boldsymbol{\theta}}_{sub}=(\hat{\alpha}_{sub},\hat{\beta}_{sub})$, we compute the sample mean squared errors (MSE) of $\hat{\boldsymbol{\theta}}_{true}$ and $\hat{\boldsymbol{\theta}}_{sub}$ over the $100$ Monte Carlo samples and plot them against $n$. The plot is given in Figure 2.

We assume that the underlying arclength parameterized manifold is $\psi([0,1])$ where $\psi:[0,1]\to\mathbb{R}^{4}$ is given by $\psi(t)=(t/2,t/2,t/2,t/2)$. We take the number of nodes at which responses are recorded to be $s=20$. Here, $m$ denotes the number of auxiliary nodes, and $n=m+s$ denotes the total number of nodes. We vary $n$ over the set $\left\{500,750,1000,.....3000\right\}$. For each $n$, we repeat the following procedure over $100$ Monte Carlo samples. A sample $t_{1},....t_{s}\sim^{iid}U[0,1]$ is generated, and responses $y_{i}$ are generated from the regression model $y_{i}=\alpha+\beta t_{i}+\epsilon_{i}$, $\epsilon_{i}\sim^{iid}N(0,\sigma^{2}_{\epsilon})$ for all $i\in[s]$, where $\alpha=2.0,\beta=5.0,\sigma_{\epsilon}=0.01$. Additionally, we sample the pre-images of the auxiliary nodes, $t_{s+1},....t_{n}\sim^{iid}U[0,1]$. Thus, a random dot product graph with latent positions $\mathbf{x}_{i}=\psi(t_{i})$, $i\in[n]$ is generated and the adjacency spectral estimates $\hat{\mathbf{x}}_{i}$ of its latent positions are computed. A localization graph is constructed with the first $n/10$ of the adjacency spectral estimates as nodes under the following rule: two nodes $\hat{\mathbf{x}}_{i}$ and $\hat{\mathbf{x}}_{j}$ of the localization graph are to be joined by an edge if and only if $\left\lVert\hat{\mathbf{x}}_{i}-\hat{\mathbf{x}}_{j}\right\rVert<\lambda$, where $\lambda=0.8\times 0.99^{K-1}$ when $n$ is the $K$-th term in the vector $(500,750,1000,\dots 3000)$. Then, the shortest path distance matrix $\mathbf{D}=\left(d_{n,\lambda}(\hat{\mathbf{x}}_{i},\hat{\mathbf{x}}_{j})\right)_{i,j=1}^{l}$ of the first $l=s+1$ estimated latent positions is embedded into one-dimension by raw-stress minimization using mds function of smacof package in R and one-dimensional embeddings $\hat{z}_{i}$ are obtained. The responses $y_{i}$ are regressed upon the $1$-dimensional raw-stress embeddings $\hat{z}_{i}$ for $i\in[s]$, and the predicted response $\tilde{y}_{s+1}$ for the $(s+1)$-th node is computed. The predicted response $\hat{y}_{s+1}$ for the $(s+1)$-th node based on the true regressors $t_{i}$ is also obtained. Due to identicality of distributions of the true regressors $t_{i}$, the distribution of each of the predicted responses is invariant to the label of the node; hence we drop the subscript and use $\hat{y}_{true}$ to denote the predicted response based on the true regressors $t_{i}$, and use $\hat{y}_{sub}$ to denote the predicted response based on the raw-stress embeddings $\hat{z}_{i}$ (since raw-stress embeddings are used as substitutes for the true regressors in predicting the response). Finally, the sample mean of the squared distances $(\hat{y}_{sub}-\hat{y}_{true})^{2}$ over all the Monte Carlo samples is computed and plotted against $n$. The plot is given in Figure 3.