Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
59,599 characters · 29 sections · 31 citation commands
Learning with latent group sparsity via heat flow dynamics on networks
Keywords. Latent group sparsity $|$ Networks $|$ Laplacian geometry $|$ Heat flow dynamics $|$ Gaussian Free Field $|$ Stochastic Block Model
The understanding and analysis of data with complex structure is a leitmotif of modern science and technology. The spectacular growth in the capacity and computational means to process gigantic volumes of data has motivated the development of novel analytical paradigms in recent years. A common theme that characterises many of these approaches is that they seek to incorporate the growing complexity that is inherent in such massive data sets national2013frontiers,marx2013big.
The intrinsic structure in data can manifest itself in various forms. These range from the almost ubiquitous scenario of sparsity in an appropriate basis, such as in compressive sensing and low-rank estimation problems achlioptas2007fast,grasedyck2013literature,chang2000adaptive,foucart2013invitation,berthet2013optimal, to algebraic constraints imposed by physical considerations, such as symmetries under rigid motions that are inherent in problems of cryo-electron microscopy cheng2015primer,singer2018mathematics,hadani2011representation,perry2019sample,fan2020likelihood,ghosh2021multi. In yet other instances, constraints may be stochastic in nature, pertaining to the statistical dependency structures that characterise the model mezard2009information,ros2019complex,ghosh2020fractal,lauritzen2019maximum,li2016fast,lavancier2015determinantal,ghosh2020gaussian,bardenet2021determinantal.
A significant structural feature that arises in many application scenarios is the clustering, or grouping, of some of the explanatory variables into a relatively limited number of categories, with the understanding that the quantities in the same category are strongly dependent on each other. A typical scenario is the existence of deterministic relationships governing the values of the variables in the same group. Further, it is often the case that only a few of these groups contribute meaningfully to the experimental observations, with the remaining variables being redundant or uninformative for predictive purposes. An important use-case of such structure is that of high dimensional supervised learning, where the explanatory variables are often clustered via natural constraints. For instance, meteorological measurements in spatially adjacent locations are likely to be highly correlated. Similarly, frequency of occurrence of certain words or phrases in spam emails are likely to be highly correlated.
A different setting in which a group structure on variables plays a significant role is that of community detection problems, where groups or clusters pertain to connectivity patterns in an underlying network. A typical example is that of a social network, where connectivity corresponds to friendship or acquaintance; yet another is that of collaboration network among scientists fortunato2010community,fortunato2007resolution,reichardt2006statistical. Due to obvious practical ramifications, this area has witnessed intense research activity in recent years, a significant achievement of which is the extensive theory of Stochastic Block Models (abbrv. SBM) abbe2017community,abbe2015exact,bandeira2016low,goldenberg2010survey,holland1983stochastic,karrer2011stochastic. In general, incorporating network geometry into standard statistical learning problems has been an area of recent interest hallac2015network, li2020graph, li2019prediction, li2020high.
The network structure brings into focus the geometric perspective on the clustering phenomenon, that is underpinned by the metric induced by the weighted graph distance in the network von2007tutorial. Intertwined with such geometry is the canonical dynamics associated to it -- in a very general Riemannian geometric setting, the metric structure gives rise to a Laplacian operator, which in turn serves as the generator of the so-called heat flow dynamics on the underlying space. These correspondences are classical in metric geometry and harmonic analysis rosenberg1997laplacian,jost2008riemannian.
Associated with the Laplacian geometry and heat flow dynamics is the canonical model in statistical physics referred to as the Gaussian Free Field (abbrv. GFF) (c.f. sheffield2007gaussian,berestycki2015introduction,friedli2017statistical), which has also emerged to be of independent interest as an important instance of Gaussian graphical models zhu2003semi,zhu2003combining,ma2013sigma,kelner2019learning,rasmussen2003gaussian. GFF-s complete the above picture from a statistical point of view, by embedding the stochastic dependency structure of Gaussian random variables in the setup of the geometric structure and dynamical properties of a weighted networks.
In the present work, we bring together these disparate strands into synergy -- clustering phenomena on variables or predictors in one direction, those in SBM-type network models in another, and the geometry of and dynamics on weighted networks in yet a third direction, along with their statistical physical implications. Leveraging their interplay, we obtain an algorithm to perform effective regression analysis in both high and low dimensional setup for variables with a latent group structure, using limited and local access to underlying network dependencies. In a more general setup, when an underlying graph may not be explicit in the problem description, we demonstrate a procedure to construct such a network based on the available data. Substantiated by rigorous mathematical analysis and robust experimental performance, our approach may be seen to outperform more classical approaches for group structured data. From an algorithmic perspective, our methodology alludes to interesting connections the with so-called diffusion mapping techniques, which have been effective as dimension reduction tools coifman2005geometric, coifman2006diffusion.
More generally, our approach opens the avenue to applications of similar diffusion-based techniques to classical statistical and data analytical problems, that are generally static in nature. The inherently local nature of the heat flow and related diffusion dynamics enables us to solve the relevant constrained optimization problems while being oblivious to the global geometry of the graph. In addition to economies of computational resources, such locality is of significance with regard to questions of privacy in data analysis, a problem that is gaining increasing salience in today's hyper-networked world.
Variable selection is a classic problem in statistics, which has become all the more important in the present age with the routine availability of large scientific datasets with measurements on tens of thousands of variables. Sparsity has become a key methodological instrument for meaningful inference from such “high-dimensional” datasets. Parsimonious models are easier to interpret and more resistant to overfitting. The Least Absolute Shrinkage and Selection Operator (abbrv. lasso) tibshirani1996regression, which employs a sparsity inducing $\ell_1$-penalty, is perhaps the most prominent method of variable selection.
A major problem with vanilla lasso is that it treats all variables equally. Thus, when there are natural groups in the variables, some variables in a group can get kicked off the model with other members still included. The group lasso penalty yuan2006model aims to solve this issue. Consider a supervised learning problem with $p$ predictors and corresponding parameter $\beta \in \mathbb{R}^p$. Denoting the groups by $\mathcal{C}_1$, \ldots, $\mathcal{C}_k$, the group lasso penalty uses a weighted $\ell_1$ norm of the groupwise $\ell_2$ norms: \[ \mathrm{GL}(\beta) = \sum_{\ell = 1}^k \sqrt{|\mathcal{C}_\ell|} \,\, \|\beta_{\mathcal{C}_\ell}\|_2. \] As a consequence, variables in a group exit the model together.
The group lasso penalty requires the groups to be known in advance. However, in many practical scenarios, the group information is a priori unknown to the statistician. There has been some work to address this issue. The cluster representative lasso (CRL) and cluster group lasso (CGL) algorithms of buhlmann2013correlated perform an initial clustering of variables into groups and then use the estimated clusters for grouped variable selection.
In this paper, we directly construct a penalty that automatically selects variables in groups without any prior group information. Notably, we dispense with the elaborate pre-processing step involving clustering of the variables, spectral or otherwise, which can be computationally intensive.
Suppose that we have explanatory variables $X_1, \ldots, X_p$, whose group structure is captured by a graph $G=(V,E)$ on $|V|=p$ vertices, each vertex corresponding to an $X_i$. In the simplest scenario, the groups of variables would correspond to the connected components $\{\mathcal{C}_i\}_{i=1}^k$ of $G$. More generally, we consider the nodes of $G$ to be a union of subsets $\{\mathcal{C}_i\}_{i=1}^k$ (corresponding to the variable groups), such that the $\mathcal{C}_i$-s are relatively densely intra-connected, but the inter-connections across the $\mathcal{C}_i$-s are relatively sparse. This is similar in flavour to the problem of multi-way partitioning of graphs, which has attracted considerable interest over the years lee2014multiway.
Let $L$ denote the Laplacian of $G$. We introduce the penalty
where $\odot$ denotes Hadamard/elementwise product of vectors, $\Phi(\beta) = (\sqrt{|\beta_1|}, \ldots, \sqrt{|\beta_p|})^\top$, and $\mathbf{1}$ is the all ones vector in $\mathbb{R}^p$. As we show in Lemma 1 and Corollary 2 in the appendix, in the limit of $t \to \infty$, the quantity $\Lambda_t(\beta)$ approaches the classical group lasso penalty $\mathrm{GL}(\beta)$ in the setting where the components $\{\mathcal{C}_i\}_{i=1}^k$ are fully disconnected from each other in the graph $G$. See Figure (ref) for a demonstration of this convergence.
Note that $\Lambda_t$ directly incorporates the Laplacian in terms of the heat flow operator $e^{-tL}$ on the underlying graph $G$. Unlike the group lasso penalty, $\Lambda_t$ is non-convex which is a potential problem vis-a-vis optimization. However, using the connection with heat flow, we can calculate $\Lambda_t$ in an efficient manner via random walks on the graph $G$.
We consider the penalised supervised learning problem
where $(X,y)$ denotes the training data, and $\mathcal{L}(\beta;X,y)$ is a suitable loss function defined with respect to the problem and computed using the training data at the parameter $\beta$. We expect that for suitably large $t$, the solution of this optimization problem will be close to that of the classical group lasso problem.
A fundamental example of this set-up is the penalised regression problem $\arg \min_{\beta} F_{t, \lambda}(\beta)$, where \[ F_{t, \lambda}^{\text{reg}}(\beta) = \frac{1}{2n} \|y - X\beta\|_2^2 + \lambda \Lambda_t(\beta). \] More generally, an important family of examples is accorded by the problem of penalised likelihood maximization, where $\mathcal{L}(\beta;X,y)$ is the negative log-likelihood of the observed data $(X,y)$ at the parameter value $\beta$. Yet another instance is that of the so-called Huber's loss, where the quadratic function of the $L_2$ norm in the regression set-up is replaced by a different convex function huber1992robust.
\algrenewcommand\algorithmicrequire{Input:} \algrenewcommand\algorithmicensure{Output:}
Our general algorithmic framework is underpinned by a heat-flow-based computation of vectors of the form $e^{-tL} f$ (or a subset of its co-ordinates) for a given $f \in \mathbb{R}^p$, where we view the latter as a function on the nodes of the graph $G$. In this setting, we observe that \[ (e^{-tL} f)_i = \mathbb{E}[f_{X(t)} \mid X(0) = i], \] where $X(t)$ is a continuous-time simple random walk on $G$, which is the canonical analogue of heat flow dynamics on in this setup. This indicates that if we start $B$ random walks $X^{(1)}, \ldots, X^{(B)}$ from the vertex $i$, then an estimate of $(e^{-tL} f)_i$ would be \[ \widehat{(e^{-tL} f)}_i = \frac{1}{B}\sum_{j = 1}^B f_{X^{(j)}(t)}. \]
In Algorithm (ref), we describe the pseudo-code for simulating $B$ heat flows from each vertex, run till time $t$. Note that the for-loops can be easily paralellised, rendering this algorithm highly efficient. The running time (i.e. computational complexity) for Algorithm (ref), even without factoring in parallelisation, is only $O(pBN_{\mathrm{step}})$, where $N_{\mathrm{step}}$ is the step count required for a single run of the heat flow. We demonstrate in (ref) that $N_{\mathrm{step}}$ is typically $O(\max\{\log p, \log n\})$; as such the computational complexity for Algorithm (ref) is only $O(B \cdot p \cdot \max\{\log p, \log n\})$.
The pseudo-code for approximating $e^{-tL} f$ (or a subset of its co-ordinates) using heat flow is presented in Algorithm (ref). Finally, the complete subgradient-based optimisation procedures are presented in Algorithms (ref) and (ref).
We use subgradient descent methods to minimise the penalised loss $F_{t,\lambda}(\beta;X, y)$ in (ref), which we denote henceforth as $F_{t,\lambda}(\beta)$ for reasons of brevity. The intermediary computation of the subgradients and losses is performed via a heat flow based local dynamics on the network. The precise algorithmic implementation of the latter is given in Algorithms (ref) and (ref), whereas the complete optimisation procedures (using these as subroutines) are encapsulated in Algorithms (ref) and (ref).
We now describe how the heat flow operator appears in the subgradient computations. Let us compute the subgradient of $F_{t, \lambda}(\beta)$. Set $h = e^{-tL} (\beta \odot \beta)$. Then $\Lambda_t(\beta) = \sum_{j = 1}^p \sqrt{|h_j|}$. Thus \[ \frac{\partial \Lambda_t(\beta)}{\partial \beta_{\ell}} = \sum_{j = 1}^p \partial s(h_j) \frac{\partial h_j}{\partial \beta_{\ell}} = \sum_{j = 1}^p \partial s(h_j) 2 (e^{-tL})_{j \ell} \beta_{\ell}, \] where $s(x) = \sqrt{|x|}$ so that $\partial s(x) = \frac{\mathrm{sgn}(x)}{2 s(x)}$. Set $\zeta = \zeta(\beta) = 2 (\partial s(h_j))_{1 \le j \le p}.$ Then, using the fact that $e^{-tL}$ is symmetric, we can write \[ \partial \Lambda_t(\beta) = (e^{-tL} \zeta) \odot \beta. \] Therefore \[ \partial F_{t, \lambda}(\beta) = \frac{1}{n}X^\top X \beta - \frac{1}{n} X^\top y + \lambda (e^{-tL} \zeta) \odot \beta. \] Thus the general subgradient step would be \[ \beta^{(m + 1)} = \beta^{(m)} - \alpha^{(m)} \partial F_{t, \lambda}(\beta^{(m)}), \] for some learning rare $\alpha^{(m)}$.
We may also use stochastic block co-ordinate descent updates. For this we need to compute $(e^{-tL} \zeta)_S$ for some $S \subset [p]$. This may be efficiently done using an on-demand computation of only the required co-ordinates of $h$, as described in Algorihm (ref).
We apply a hard-thresholding step to the output $\hat{\beta}$ of the optimisation to further enhance its support recovery properties. We apply $K$-means clustering with $K = 2$ on $|\hat{\beta}|$. Let $\mathcal{C}$ be the cluster whose mean is the closest to $0$. We zero out the values of $\hat{\beta}_{\mathcal{C}}$. The pseudo-code of this procedure appears in Algorithm (ref).
If the Laplacian $L$ (equivalently, the graph $G$) is not a priori available, we may estimate it using the covariance structure of the explanatory variables as follows (the pseudo-code appears in Algorithm (ref)). Let $\hat{\Sigma}$ be some estimate of the true covariance matrix $\Sigma$. Let $\hat{R}$ denote the estimate of the true correlation matrix, obtained by rescaling $\hat{\Sigma}$. We construct an adjacency matrix by thresholding $\hat{R}$: \[ A_{ij} = \mathbf{1}_{\{|\hat{R}_{ij}| \ge \tau(\hat{R})\}}, \] where $\tau(\hat{R})$ is an appropriate quantile of the $|\hat{R}_{ij}|$'s. (E.g., in our simulation experiments, we use the $0.75$-th quantile.) We then take the unnormalised Laplacian corresponding to $A$ as an estimate of $L$.
As for estimating $\Sigma$, we may use the sample covariance matrix in low to moderate dimensions, and appropriately thresholded versions of it (e.g., bickel2008covariance, bickel2008regularized) or other shrinkage estimators in high dimensions (e.g., friedman2008sparse, chen2010shrinkage).
A cornerstone of our algorithms is that we do not require any direct knowledge of the group structure. Clustering algorithms for estimating the group structure typically require the knowledge of the number of groups. Thus in order to use group Lasso with an estimated group structure, one needs to know the number of groups. Our algorithms, on the other hand, do not have such a requirement.
Further, we do not even require access to the full network. Rather, at each step of the optimisation algorithm, we find an approximation to the subgradient of the penalty (at the current value of the parameter $\beta$) by using an Monte Carlo based approach. Additional economy in computational resources is accorded by the fact that we are able to use the same heat flow throughout the optimization procedure, the initial and terminal nodes of which generated and stored for successive calls. The effectiveness of this economizing protocol is demonstrated in simulation experiments as well as applications to real world data.
In particular, running a continuous time random walk only entails exploring a local neighbourhood of the current state, which relieves us of the necessity to work with the entire graph at one go. The latter would be necessary, e.g., in an approach where we wanted to do an initial spectral clustering in order to arrive at a detailed understanding of the group structure buhlmann2013correlated. This may be prohibitively expensive in the setting of real-world massive networks, such as the large-scale social networks or the world-wide web.
In fact, we only require oracle access to a black box that returns the terminal state of the heat flow starting from a prescribed initial node. Such limited and local access to the network data can have significant implications with regard to considerations of privacy and security, which play an increasingly important role in modern statistical research.
The duration of the heat flow is determined by the goal to make the difference between our penalty and the classical group lasso penalty small, a bound on which is accorded by Theorem (ref). It follows therefrom that, under reasonable conditions, it suffices to have $t_{\mathrm{flow}} \gg \frac{1}{\lambda_g}\max \{ \log n, \log p\}$. In the most important setting of the clusters being fully disconnected from each other but densely connected within each other, the ground state $\lambda_g$ is the minimum of the ground states for the individual clusters. For the components being generic densely connected graphs of size $\Theta(p)$, we typically have $\lambda_g=\Theta(p)$, whence $t_{\mathrm{flow}} \gg \frac{1}{p}\max \{ \log n, \log p\}$; we refer to Section 3 in the appendix for details. We note that although the optimal choice of $t_{\mathrm{flow}}$ depends on $\lambda_g$, in practice we can perform cross-validation to choose the optimal value of $t_{\mathrm{flow}}$. It is precisely for this reason that our algorithms do not need the knowledge of the number of groups.
The number of steps of the heat flow dynamics is, roughly speaking, the heat flow time $t_{\mathrm{flow}}$ times the average number of steps per unit time. If the current state of the dynamics is a node $v \in G$ with degree $\mathrm{deg}(v)$, then the time until the next step in the random walk is an $\mathrm{Exponential}(\mathrm{deg}(v))$) random variable. Thus, step time distribution is stochastically dominated by an $\mathrm{Exponential}(d_{\max})$ random variable, where $d_{\max}$ is the maximum degree of $G$. Thus, on average, the total step count for the heat flow dynamics is given by $N_{\mathrm{step}}=O(d_{\max} \cdot t_{\mathrm{flow}})=O(d_{\max} \cdot \frac{1}{\lambda_g} \cdot \max \{ \log n, \log p\})$. For typical strongly intra-connected components as above, we have $d_{\max}=\Theta(p)$ (for details we refer to Section 3 in the appendix), whence we have $N_{\mathrm{step}} = O(\max \{ \log n, \log p\})$.
The logarithmic dependence on the problem dimensions as indicated above entails a light computational load for the heat flow based dynamic algorithm.
We note in passing that the spectral data of a heat flow operator (equivalently, an appropriate random walk transition matrix) is known to be useful for learning tasks, especially in the context of diffusion mapping for dimensionality reduction problems coifman2005geometric, coifman2006diffusion. However, the random walk in diffusion mapping takes place in a different space -- namely, on the space of the actual data points (therefore, having $n$ nodes); whereas in our setting, the random walk takes place on the co-ordinate indices of the data points (thereby entailing $p$ nodes). Nonetheless, it would be of interest to explore the possible interaction of these two diffusion-based approaches, in particular regarding the possibility of incorporating a dimension reduction step in our paradigm to achieve further economy of computational resources.
It may be noted that the map $\Psi$ that is embedded in our penalty is mildly non-convex, with an algebraically explicit square-root structure. This may be compared to the setting of classical group lasso when the groups are a priori not known, where non-convexity enters via its role in estimating the groups. This step involves solving a clustering problem, which is generally non-convex. In fact, the non-convex structure of the group recovery problem is known for its notoriety as a challenging optimization issue, and its rather intractable combinatorial nature makes it arguably a more complicated endeavour than the simple, algebraically explicit non-convexity posed by the map $\Psi$ in our approach. Exploiting this simple algebraic structure, perhaps by carrying out the optimization in a different co-ordinate system, would be an interesting direction for future research.
We first set some notations in order to lay out our main theoretical results. For definiteness, we focus on the setting of regression with a group structure on the parameters in our theoretical analysis. However, we note in passing that similar analysis would apply to a wide range of applications with our method, including logistic regression, generalised linear models and other use cases.
Consider the situation where the graph $G$ has exactly $k$ connected components. Then it is a well-known fact that the spectrum $0 = \lambda_0 \le \lambda_1 \le \cdots \le \lambda_p$ of the Laplacian matrix $L$ of $G$ then has exactly $k$ many zero eigenvalues. Let $\lambda_g := \min_{i > k} \lambda_i$ denote the spectral gap of $L$. Let \[ A(\beta) = \{ i \mid 1 \le i \le k, \|\beta_{\mathcal{C}_i}\|_2 \ne 0\}, \] and \[ I(\beta) = \cup_{i \in A(\beta)} \mathcal{C}_i. \]
Note that $\Lambda_t(\beta)$ is a non-convex penalty. As such, the objective function $F_{t, \lambda}$ can have multiple local optima. However, we will show that in an appropriate neighbourhood of the true parameter value $\beta^*$, the minimum of $F_{t, \lambda}(\beta)$ approximately minimises the group Lasso penalty for sufficiently large $t$.
As we shall see below, the presence of an additional restricted eigenvalue property leads to improved guarantees on the accuracy of our procedure, so we state this property below.
Property RE($s$): We say that the restricted eigenvalue property RE($s$) holds for $X$ with parameter $\kappa = \kappa(s)$ if
We denote $B(\beta^*,\varepsilon)$ to be the $L^2$ Euclidean ball in $\mathbb{R}^p$ with center $\beta^*$ and radius $\varepsilon$, and set \[ \Lambda(X; \eta) = \max_{1 \le j \le K} \sqrt{\|\frac{1}{n}X_j^\top X_j\|_{\mathrm{op}}}\bigg(1 + \sqrt{\frac{4\log \eta^{-1}}{|\mathcal{C}_j|}}\bigg). \] We may then state:
We further demonstrate that the RE($s$) property above for a random design matrix $X$ can be related, in the crucial setting of Gaussian covariates, to a similar property for the deterministic covariance matrix $\Sigma$ of the rows of $X$. The latter can often be much easier to verify; for instance, it may be easily seen to hold as soon as $\Sigma$ is well-conditioned.
Gaussian Free Fields (abbrv. GFF) have emerged as important models of correlated Gaussian fields, that are naturally commensurate with the geometry of their ambient space. In the case of graphs, the ambient geometry is spawned by the graph Laplacian. These Gaussian processes also have important applications in physics, where they are of interest in the context of Euclidean quantum field theories friedli2017statistical.
GFFs are in fact Gaussian Graphical Models (abbrv. GGM), where the precision matrix of the Gaussian random field aligns with the Laplacian of the graph, thereby leading to a rich interaction between the statistical properties of the GGM and the Laplacian geometry of the underlying graph zhu2003semi,zhu2003combining,ma2013sigma,kelner2019learning,rasmussen2003gaussian. GGMs have emerged as popular tools to model dependency relationships in data via a latent graph structure, the choice of Gaussian randomness being often motivated by the fact that a Gaussian distribution maximises entropy within the constraints of a given covariance structure. Applications of GGMs are ubiquitous, with use cases in diverse domains such as structural inference in biological networks, causal inference problems, speech recognition, and so on whittaker2009graphical,lauritzen1996graphical,edwards2012introduction,uhler2019gaussian. Since our approach exploits in an essential manner the Laplacian geometry of the graph, the GFF is a natural GGM to examine within its ambit.
For an in-depth introduction to the technical aspects of GFFs, we refer the reader to the excellent surveys sheffield2007gaussian and berestycki2015introduction; herein we will content ourselves with a brief description of its relevant features. Broadly speaking, a GFF is essentially a natural generalization of Brownian motion to general spaces, with time replaced by, e.g., the nodes of a network. We define the massive GFF $\mathbb{X}_\theta=(\mathbb{X}_\theta(v))_{v \in V}$ on a graph $G=(V,E)$ with mass parameter $\theta>0$ to be a mean-zero Gaussian field indexed by $V$ characterised by its precision matrix (i.e., the inverse of the covariance matrix) given by $(L + \theta I_{|V|})$, where $L$ is the (unnormalised) graph Laplacian on $G$ and $I_{|V|}$ is the $|V| \times |V|$ identity matrix (c.f. berestycki2015introduction). The covariance matrix of $\mathbb{X}_\theta$ is therefore given by $\Sigma=(L + \theta I_{|V|})^{-1}$. Note that $L$ itself is singular due to the all ones vector $\mathbbm{1}_V$ being in its kernel; therefore we are aided by the strong convexity accorded by the mass parameter $\theta>0$.
In this section, we consider a very different model of graph clustering that is motivated by Stochastic Block Models (abbrv., SBM) that have attracted intense focus in statistics and machine learning applications in recent years abbe2017community,goldenberg2010survey,holland1983stochastic,karrer2011stochastic. SBMs are underpinned by a block matrix structure, indexed by the vertices of a graph, with entries in $[0,1]$ that are constant on the blocks. These entries are the connection probabilities between the respective vertices. In our setting, it would be natural to consider the same block matrix pattern to generate the graph and the covariance structure.
For definiteness, we divide the vertex set $V$ in $k$ groups $\{\mathcal{C}_i\}_{i=1}^k$. We consider the $|V| \times |V|$ matrix $\mathcal{P}$ whose rows and columns are indexed by the nodes of the graph $G$. For $1 \le i,j \le k$, let $\mathcal{P}^{i,j}$ denote the submatrix indexed by the vertices in the groups $\mathcal{C}_i$ (along the rows) and $\mathcal{C}_j$ along the columns. For a vector $v \in \mathbb{R}^{|V|}$, we denote by $D(v)$ the diagonal matrix whose diagonal equals $v$. For brevity, we will denote by $\mathbbm{1}_i$ the indicator vector (in $\mathbb{R}^{|V|}$) of the nodes in $\mathcal{C}_i$. Let $a>b$ be numbers in $[0,1]$, such that $\mathcal{P}^{i,i} = a (\mathbbm{1}_i \mathbbm{1}_i^\top - D(\mathbbm{1}_i))$, and for $1 \le i \ne j \le k$, we have $\mathcal{P}^{i,j} = b \mathbbm{1}_i \mathbbm{1}_j^\top$. Typically, $k$ is small compared $|V|$ and $a-b$ is taken to be large enough.
The classical SBM is a random graph that is sampled with independent edges according to probabilities given by the matrix $\mathcal{P}$. In our setting, the graph $G$ would simply be a realisation of this SBM, with the parameter $b$ being small so as to ensure the graph to be sparsely connected across blocks. An alternative paradigm of graph generation from the matrix $\mathcal{P}$ would be to consider a deterministic, albeit weighted, graph $G$, whose edge weight between the vertices $u,v \in V$ equals $\mathcal{P}_{uv}$. The correlation structure of the covariates for either model is given by the covariance matrix $\Sigma=I_{|V|} + \mathcal{P}$. In practice, if a graph is constructed from such a block-structured sample covariance matrix via thresholding, then the resulting graph would be well-modelled by the above SBM.
In this section, we provide ballpark estimates on the quantities of interest in our random design models, that would be enough to guarantee the effectiveness of our approach. We record here the main conclusions of our analysis, postponing the details to the appendix.
For definiteness, we fix a polynomial decay of probability with which our recovery guarantees are to hold, which implies that the quantity $\eta=O(n^{-\alpha})$ in Theorem (ref) for a fixed $\alpha>0$. Under this error tolerance, the appropriate choice of $\lambda$ may be shown to be $\lambda \gtrsim \sigma_{\max}(\Sigma) \sqrt{\frac{\log n}{n}}$, which in turn leads to a prediction guarantee of
Our goal here is to understand the order of the flow time $t_{\mathrm{flow}}$ and the step count $N_{\mathrm{step}}$ (in terms of the other parameters of the problem) required to achieve a desired accuracy. We will bifurcate our analysis into two related parts.
Given the data size $n$, we investigate the order of $t=t_{\mathrm{flow}}$ at which the approximation error due to our heat flow based approach (roughly, the second term in (ref)) becomes comparable to the contribution to the prediction error bound for the classical group lasso methods that assume complete knowledge of the group structure (roughly, the first term in (ref)). It may be shown that under very general circumstances, we have
whereas for most models of interest, including the GFF and SBM based models in our purview, we may deduce the much simpler prescription $t_{\mathrm{flow}} ~\gtrsim \frac{1}{\lambda_g} \max \{ \log p , \log n \}$. The details are provided in Section 3 of the appendix.
We fix a threshold $\varepsilon$, and make explicit prescriptions for the order of $n$ and $t_{\mathrm{flow}}$ that will allow us to obtain a prediction error of order $O(\varepsilon)$. To this end, we posit that the two terms on the right hand side of (ref) are separately $O(\varepsilon^2)$. For $t_{\mathrm{flow}}$, this entails that $p^{3/2}e^{-\lambda_g t_{\mathrm{flow}}} \lesssim \varepsilon^2$, which translates into $t_{\mathrm{flow}} \gtrsim \frac{1}{\lambda_g}\left ( \log p + \log \frac{1}{\varepsilon} \right)$.
For a prescription for $n$ given a target prediction error $\varepsilon$, we have the bound
The details of the analysis are available in Section 3 of the appendix.
For typical clustered networks (for a concrete probabilistic model, see Section 3 in the appendix), we can further simplify above prescriptions on the heat flow time $t_{\mathrm{flow}}$ to the thresholds
resp. for the settings where data size $n$ is given and where the target prediction error $\varepsilon$ is given. These translate into the step count bounds
We refer the reader to Sec. 3 in the appendix for details. It may be noted that the flow time $t_{\mathrm{flow}}$ and the step count $N_{\mathrm{step}}$ do not depend on the statistical properties of the covariates $X_i$, but depend only on the geometric properties of the network $G$.
In the setting of GFF and SBM based random designs, we may obtain further explicit guarantees on the sample complexity $n$, via application of techniques from spectral graph theory; the detailed analysis is provided in Section 3 of the appendix.
For both SBM based and GFF based random designs on well clustered networks, we have the following prescriptions that suffice with high probability:
and
We observe that the sample complexity bound (ref) for an SBM based design matches, upto logarithmic factors, the analogous bound for classical sparse reconstruction problems. For GFF based random designs, it may be noted that the $\frac{1}{\varepsilon^2} \cdot \|\beta^*\|_0 \cdot p$ term in (ref) comes from the first term in (ref), which, roughly speaking, reflects the error incurred by classical group lasso. Thus, the linear dependence of the sample complexity on $p$ appears to be a fundamental characteristic of the problem for GFF based random designs, and is inherent to both classical group lasso and the present heat flow based methods. Empirical investigations also appear to corroborate this effect; for details we refer the reader to Section 4 in the appendix. Although our method is still effective in a high dimensional setup (up to logarithmic factors), devising methodologies with greater efficiency for GFF based random designs, as well as theoretical investigations of information theoretic lower bounds in this setting, would be an interesting direction for future research.
We take $n = 200, p = 100$ and $k = 4$ groups of relative sizes $(p_1, p_2, p_3, p_4)^\top / p = (0.16, 0.24, 0.40, 0.20)^\top$. We denote these groups by $\mathcal{C}_i, i = 1, \ldots, k$.
We consider two models for $X$:
The true parameter $\beta$ is generated in the following way \[ \beta_i
\]
Finally, we generate the response from the linear model \[ Y = X\beta + \varepsilon, \] where $\varepsilon$ is a noise vector, independent of $X$, with isotropic covariance matrix $\sigma^2 I_n$.
We compare the subgradient and the stochastic block co-ordinate descent versions of our procedure (referred to as “Heat flow (SD)” and “Heat flow (CD)”, respectively) against group lasso with group structure learned from spectral clustering von2007tutorial on $\hat{L}$, an estimate of the Laplacian matrix obtained using Algorithm (ref) with oracle knowledge of $k = 4$. For each algorithm, the tuning parameter $\lambda$ is chosen by cross-validation. In Tables (ref) and (ref), we report the prediction error, estimation error, and two measures of support recovery, namely, sensitivity and specificity:
The estimated $\beta$ for different methods in the experiment with block diagonal covariate structure is shown in Figure (ref).
From both the experiments we see that our proposed method based on the heat flow penalty has comparable performance to group lasso in terms of prediction/estimation error, without requiring an explicit knowledge of the group structure, or even of the number of groups. In terms of support recovery, the proposed method appears to be much superior.
In this section, we compare our methods against group lasso in terms of test-set performance in four real-world data sets. For group lasso, the group structure is learned by applying spectral clustering on $\hat{L}$, the estimated Laplacian matrix obtained by applying Algorithm (ref), with $k$ set to be the number of eigenvalues of $\hat{L}$ less than $0.01$. For each data set, we use an 80:20 split into training and test sets. The test-set errors are reported in Table (ref). We observe comparable performance in all four data sets.
We consider the well-known spambase data set\footnote{https://archive.ics.uci.edu/ml/datasets/spambase} containing 4601 emails classified as spam/non-spam. There are 57 explanatory variables. We fit a logistic regression model with group lasso and heat flow penalties. We report the test-set misclassification error in the second column of Table (ref).
We consider gene expression data from the microarray experiments of mammalian eye tissue samples of scheetz2006regulation. The response variable is the expression level of the TRIM32 gene. There are 200 predictor variables corresponding to different gene probes. The sample size is 120. We report the test-set mean-squared error (MSE) in the third column of Table (ref).
From the NCEP/NCAR reanalysis data set, we took the monthly average temperature of the Delhi-NCR region as the response variable. We took monthly average temperature, pressure, precipitation, wind-speed, etc. of $2.5^\circ \times 2.5^\circ$ blocks on the Bay of Bengal and the Arabian Sea as covariates. In total there were 101 such blocks, giving us $p = 606$ explanatory variables. We have these measurements for $n = 886$ months, starting January, 1947 till October 2021. We first removed seasonal variations and fitted a linear trend afterwards as preprocessing steps (as described in chatterjee2012sparse). We report the test-set MSE in the fourth column of Table (ref).
We have data on daily highs of the NIFTY 50 index from the National Stock Exchange (NSE) of India for 49 companies for $n = 2598$ days staring from November 4, 2010 till April 30, 2021. We use the mean index of $9$ companies in the financial sector as our response variable, and use the indices of the rest of the $p = 40$ companies as covariates. We report the test-set MSE in the fifth column of Table (ref).
In this work, we contribute an approach to learning under a group structure on explanatory variables that does not require prior information on the group identities. Our paradigm is motivated by the Laplacian geometry of an underlying network with a commensurate community structure, and proceeds by directly incorporating this into the penalty. In a more general setup, when an underlying graph may not be explicit in the problem description, we demonstrate a procedure to construct such a network based on the available data. Notably, we dispense with the elaborate pre-processing step involving clustering of the variables, spectral or otherwise, which can be computationally resource-intensive. Our paradigm is underpinned by rigorous theorems that guarantee effective performance and provide bounds on its sample complexity. In particular, we demonstrate that in a very wide range of settings, we need to run the heat flow dynamics for a time that is only logarithmic in the problem dimensions. We investigate in detail the interplay of our approach with key statistical physics paradigms such as the GFF and the SBM. We validate our approach by successful application to real-world data from diverse fields including computer science, genetics, climatology and economics.
Our approach opens the avenue to applications of similar dynamical techniques to classical statistical and data analytical problems, that are normally defined as static problems. The inherently local nature of the heat flow and related diffusion dynamics enables us to resolve the relevant constrained optimization problems while being oblivious to the global geometry of the graph (such as a complete understanding of the clustering structure of the variables). In addition to economies of computational resource, such locality is of significance in the context of questions of privacy in data analytical methodologies, a problem that is gaining increasing salience in today's hyper-networked world. On a related note, it would be of interest to enhance our approach to obtain similarly local algorithms that address additional structural features of the explanatory variables, such as smoothness or intra-group sparsity. Yet another intriguing direction would be to explore the interface of our approach and diffusion-mapping based techniques that have been effective for dimension reduction problems, and exploit their interplay to achieve further economy of scale and computational resources. In general, the interplay between the geometric structure provided by the Laplacian, the stochastic structure accorded by models such as the GFF and SBM and the inherent clustering structure of real world datasets raises the possibility of a rich mathematical theory and a suite of associated techniques to evolve.
S.G. is supported in part by the MOE grants R-146-000-250-133 and R-146-000-312-114. S.S.M. is supported by an INSPIRE Faculty Fellowship from the Department of Science and Technology, Government of India. The authors thank Snigdhansu Chatterjee for pointers to the NCEP/NCAR reanalysis data set.