A Distributed Method for Fitting Laplacian Regularized Stratified Models
Abstract
Stratified models are models that depend in an arbitrary way on a set of selected categorical features, and depend linearly on the other features. In a basic and traditional formulation a separate model is fit for each value of the categorical feature, using only the data that has the specific categorical value. To this formulation we add Laplacian regularization, which encourages the model parameters for neighboring categorical values to be similar. Laplacian regularization allows us to specify one or more weighted graphs on the stratification feature values. For example, stratifying over the days of the week, we can specify that the Sunday model parameter should be close to the Saturday and Monday model parameters. The regularization improves the performance of the model over the traditional stratified model, since the model for each value of the categorical ‘borrows strength’ from its neighbors. In particular, it produces a model even for categorical values that did not appear in the training data set.
We propose an efficient distributed method for fitting stratified models, based on the alternating direction method of multipliers (ADMM). When the fitting loss functions are convex, the stratified model fitting problem is convex, and our method computes the global minimizer of the loss plus regularization; in other cases it computes a local minimizer. The method is very efficient, and naturally scales to large data sets or numbers of stratified feature values. We illustrate our method with a variety of examples.
1 Introduction
We consider the problem of fitting a model to some given data. One common and simple paradigm parametrizes the model by a parameter vector , and uses convex optimization to minimize an empirical loss on the (training) data set plus a regularization term that (hopefully) skews the model parameter toward one for which the model generalizes to new unseen data. This method usually includes one or more hyperparameters that scale terms in the regularization. For each of a number of values of the hyperparameters, a model parameter is found, and the resulting model is tested on previously unseen data. Among these models, we choose one that achieves a good fit on the test data. The requirement that the loss function and regularization be convex limits how the features enter into the model; generally, they are linear (or affine) in the parameters. One advantage of such models is that they are generally interpretable.
Such models are widely used, and often work well in practice. Least squares regression [1, 2], lasso [3], logistic regression [4], support vector classifiers [5] are common methods that fall into this category (see [6] for these and others). At the other extreme, we can use models that in principle fit arbitrary nonlinear dependence on the features. Popular examples include tree based models [7] and neural networks [8].
Stratification is a method to build a model of some data that depends in an arbitrary way on one or more of the categorical features. It identifies one or more categorical features and builds a separate model for the data that take on each of the possible values of these categorical features.
In this paper we propose augmenting basic stratification with an additional regularization term. We take into account the usual loss and regularization terms in the objective, plus an additional regularization term that encourages the parameters found for each value to be close to their neighbors on some specified weighted graph on the categorical values. We use the simplest possible term that encourages closeness of neighboring parameter values: a graph Laplacian on the stratification feature values. (As we will explain below, several recent papers use more sophisticated regularization on the parameters for different values of the categorical.) We refer to a startified model with this objective function as a Laplacian regularized stratified model. Laplacian regularization can also be interpreted as a prior that our model parameters vary smoothly across the graph on the stratification feature values.
Stratification (without the Laplacian regularization) is an old and simple idea: Simply fit a different model for each value of the categorical used to stratify. As a simple example where the data are people, we can stratify on sex, i.e., fit a separate model for females and males. (Whether or not this stratified model is better than a common model for both females and males is determined by outofsample or crossvalidation.) To this old idea, we add an additional term that penalizes deviations between the model parameters for females and males. (This Laplacian regularization would be scaled by a hyperparameter, whose value is chosen by validation on unseen or test data.) As a simple extension of the example mentioned above, we can stratify on sex and age, meaning we have a separate model for females and males, for each age from 1 to 100 (say). Thus we would create a total of 200 different models: one for each age/sex pair. The Laplacian regularization that we add encourages, for example, the model parameters for 37 year old females to be close to that for 37 year old males, as well as 36 and 38 year old females. We would have two hyperparameters, one that scales the deviation of the model coefficients across sex, and another that scales the deviation of the model coefficients across adjacent ages. One would choose suitable values of these hyperparameters by validation on a test set of data.
There are many other applications where a Laplacian regularized stratified model offers a simple and interpretable model. We can stratify over time, which gives timevarying models, with regularization that encourages our timedependent model parameters to vary slowly over time. We can capture multiple periodicities in timevarying models, for example by asserting that the 11PM model parameter should be close to the midnight model parameter, the Tuesday and Wednesday parameters should be close to each other, and so on. We can stratify over space, after discretizing location. The Laplacian regularization in this case encourages the model parameters to vary smoothly over space. This idea naturally falls in line with Waldo Tobler’s “first law of geography” [9]:
Everything is related to everything else, but near things are more related than distant things.
Compared to simple stratified models, our Laplacian regularized stratified models have several advantages. Without regularization, stratification is limited by the number of data points taking each value of the stratification feature, and of course fails (or does not perform well) when there are no data points with some value of the stratification feature. With regularization there is no such limit; we can even build a model when the training data has no points with some values of the stratification feature. In this case the model borrows strength from its neighbors. Continuing the example mentioned above, we can find a model for 37 year old females, even when our data set contains none. Unsurprisingly, the parameter for 37 year old females is a weighted average of the parameters for 36 and 38 year old females, and 37 year old males.
We can think of stratified models as a hybrid of complex and simple models. The stratified model has arbitrary dependence on the stratification feature, but simple (typically linear) dependence on the other features. It is essentially nonparametric in the stratified feature, and parametric with a simple form in the other features. The Laplacian regularization allows us to come up with a sensible model with respect to the stratified feature.
The simplest possible stratified model is a constant for each value of the stratification feature. In this case the stratified model is nothing more than a fully nonparametric model of the stratified feature. The Laplacian regularization in this case encourages smoothness of the predicted value relative to some graph on the possible stratification feature values. More sophisticated models predict a distribution of values of some features, given the stratified feature value, or build a stratified regression model, for example.
In this paper, we first describe stratified model fitting with Laplacian regularization as a convex optimization problem, with one model parameter vector for every value of the stratification feature. This problem is convex (when the loss and regularizers in the fitting method are), and so can be reliably solved [10]. We propose an efficient distributed solution method based on the alternating direction method of multipliers (ADMM) [11], which allows us to solve the optimization problem in a distributed fashion and at very large scale. In each iteration of our method, we fit the models for each value of the stratified feature (independently in parallel), and then solve a number of Laplacian systems (one for each component of the model parameters, also in parallel). We illustrate the ideas and method with several examples.
2 Related work
Stratification, i.e., the idea of separately fitting a different model for each value of some parameter, is routinely applied in various practical settings. For example, in clinical trials, participants are often divided into subgroups, and one fits a separate model for the data from each subgroup [12]. Stratification by itself has immediate benefits; it can be helpful for dealing with confounding categorical variables, can help one gain insight into the nature of the data, and can also play a large role in experiment design (see, e.g., [13, 14] for some applications to clinical research, and [15] for an application to biology).
The idea of adding regularization to fitting stratified models, however, is (unfortunately) not as well known in the data science, machine learning, and statistics communities. In this section, we outline some ideas related to stratified models, and also outline some prior work on Laplacian regularization.
Regularized lasso models.
The datashared lasso is a stratified model that encourages closeness of parameters by their difference as measured by the norm [16]. (Laplacian regularization penalizes their difference by the norm squared.) The pliable lasso is a generalization of the lasso that is also a stratified model, encouraging closeness of parameters as measured by a blend of the  and norms [17, 18]. The network lasso is a stratified model that encourages closeness of parameters by their difference as measured by the norm [19, 20]. The network lasso problem results in groups of models with the same parameters, in effect clustering the stratification features. (In contrast, Laplacian regularization leads to smooth parameter values.) The ideas in the original papers on the network lasso were combined and made more general, resulting in the software package SnapVX [21], a general solver for convex optimization problems defined on graphs that is based on SNAP [22] and cvxpy [23]. One could, in principle, implement Laplacian regularized stratified models in SnapVX; instead we develop specialized methods for the specific case of Laplacian regularization, which leads to faster and more robust methods. We also have observed that, in most practical settings, Laplacian regularization is all you need.
Varyingcoefficient models.
Varyingcoefficient models are a class of regression and classification models in which the model parameters are smooth functions of some features [24]. A stratified model can be thought of as a special case of varyingcoefficient models, where the features used to choose the coefficients are categorical. Estimating varyingcoefficient models is generally nonconvex and computationally expensive [25].
Geographically weighted regression.
The method of geographically weighted regression (GWR) is widely used in the field of geographic information systems (GIS). It consists of a spatiallyvarying linear model with a location attribute [26]. GWR fits a model for each location using every data point, where each data point is weighted by a kernel function of the pairwise locations. (See, e.g., [27] for a recent survey and applications.) GWR is an extremely special case of a stratified model, where the task is regression and the stratification feature is location.
Graphbased feature regularization.
One can apply graph regularization to the features themselves. For example, the group lasso penalty enforces sparsity of model coefficients among groups of features [28, 29], performing feature selection at the level of groups of features. The group lasso is often fit using a method based on ADMM, which is similar in vein to the one we describe. The group lasso penalty can be used as a regularization term to any linear model, e.g., logistic regression [30]. A related idea is to use a Laplacian regularization term on the features, leading to, e.g., Laplacian regularized leastsquares (LapRLS) and Laplacian regularized support vector machines (LapSVM) [31]. Laplacian regularized models have been applied to many disciplines and fields, e.g., to semisupervised learning [32, 33], communication [34, 35, 36], medical imaging [37], computer vision [38], natural language processing [39], and microbiology [40].
Graph interpolation.
In graph interpolation or regression on a graph, one is given vectors at some of the vertices of a graph, and is tasked with inferring the vectors at the other (unknown) vertices [41]. This can be viewed as a stratified model where the parameter at each vertex or node is a constant. This leads to, e.g., the total variation denoising problem in signal processing [42].
Multitask learning.
In multitask learning, the goal is to learn models for multiple related tasks by taking advantage of pairwise relationships between the tasks (see, e.g., [43, 44]). Multitask learning can be interpreted as a stratified model where the stratification feature is the task to be performed by that model. Laplacian regularization has also been investigated in the context of multitask learning [45].
Laplacian systems.
Laplacian matrices are very well studied in the field of spectral graph theory [46, 47]. A set of linear equations with Laplacian coefficient matrix is called a Laplacian system. Many problems give rise to Laplacian systems, including in electrical engineering (resistor networks), physics (vibrations and heat), computer science, and spectral graph theory (see, e.g., [48] for a survey). While solving a general linear system requires order flops, where is the number of vertices in the graph, Laplacian system solvers have been developed that are (nearly) linear in the number of edges contained in the associated graph [49, 50, 51]. In many cases, however, the graph can be large and dense enough that solving a Laplacian system is still a highly nontrivial computational issue; in these cases, graph sparsification can be used to yield results of moderate accuracy in much less time [52]. In practice, the conjugate gradient (CG) method [53], with a suitable preconditioner, can be extremely efficient for solving the (positive definite) Tikhonovregularized Laplacian systems that we encounter in this paper.
Laplacian regularized minimization problems.
Stratified model fitting with Laplacian regularization is simply a convex optimization problem with Laplacian regularization. There are many generalpurpose methods that solve this problem, indeed too many to name; in this paper we derive a problemspecific algorithm using ADMM. Another possibility would be majorizationminimization (MM), an algorithmic framework where one iteratively minimizes a majorizer of the original function at the current iterate [54]. One possible majorizer for the Laplacian regularization term is a diagonal quadratic, allowing block separable problems to be solved in parallel; this was done in [55]. The MM and ADMM algorithms have been shown to be closely connected [56], since the proximal operator of a function minimizes a quadratic majorizer. However, convergence with MM depends heavily on the choice of majorizer and the problem structure [57]. In contrast, the ADMMbased algorithm we describe in this paper typically converges (at least to reasonable practical tolerance) within 100–200 iterations.
3 Stratified models
In this section we define stratified models and our fitting method; we describe our distributed method for carrying out the fitting in §4.
We consider data fitting problems with items or records of the form . Here is the feature over which we stratify, is the other features, and is the outcome or dependent variable. For our purposes, it does not matter what or are; they can consist of numerical, categorical, or other data types. The stratified feature values , however, must consist of only possible values, which we denote as . What we call can include several original features, for example, sex and age. If all combinations of these original features are possible, then is the product of the numbers of values these features can take on. In some formulations, described below, we do not have ; in this case, the records have the simplified form .
Base model.
We build a stratified model on top of a base model, which models pairs (or, when is absent, just ). The base model is parametrized by a parameter vector . In a stratified model, we use a different value of the parameter for each value of . We denote these parameters as , where is the parameter value used when . We let denote the collection of parameter values, i.e., the parameter value for the stratified model.
Local loss and regularization.
Let , , denote a set of training data points or examples. We will use regularized empirical loss minimization to choose the parameters . Let be a loss function. We define the th (local) empirical loss as
(1) 
(without when is absent).
We may additionally add a regularizer to the local loss function. The idea of regularization can be traced back to as early as the 1940s, where it was introduced as a method to stabilize illposed problems [58, 59]. Let be a regularization function or regularizer. Choosing to minimize gives the regularized empirical risk minimization model parameters, based only on the data records that take the particular value of the stratification feature . This corresponds to the traditional stratified model, with no connection between the parameter values for different values of . (Infinite values of the regularizer encode constraints on allowable model parameters.)
Laplacian regularization.
Let be a symmetric matrix with nonnegative entries. The associated Laplacian regularization is the function given by
(2) 
(To distinguish the Laplacian regularization from in (1), we refer to as the local regularization.) We can associate the Laplacian regularization with a graph with vertices, which has an edge for each positive , with weight . We refer to this graph as the regularization graph. We can express the Laplacian regularization as the positive semidefinite quadratic form
where denotes the Kronecker product, and is the (weighted) Laplacian matrix associated with the weighted graph, given by
for .
The Laplacian regularization evidently measures the aggregate deviation of the parameter vectors from their graph neighbors, weighted by the edge weights. Roughly speaking, it is a metric of how rough or nonsmooth the mapping from to is, measured by the edge weights.
Fitting the stratified model.
To choose the parameters , we minimize
(3) 
The first term is the sum of the local objective functions, used in fitting a traditional stratified model; the second measures the nonsmoothness of the model parameters as measured by the Laplacian. It encourages parameter values for neighboring values of to be close to each other.
Convexity assumption.
We will assume that the local loss function and the local regularizer are convex functions of , which implies that that the local losses are convex functions of . The Laplacian is also convex, so the overall objective is convex, and minimizing it (i.e., fitting a stratified model) is a convex optimization problem. In §4 we will describe an effective distributed method for solving this problem. (Much of our development carries over to the case when or is not convex; in this case, of course, we can only expect to find a local minimizer of ; see §6.)
The two extremes.
Assuming the regularization graph is connected, taking the positive weights to forces all the parameters to be equal, i.e., . We refer to this as the common model (i.e., one which does not depend on ). At the other extreme we can take all to , which results in a traditional stratified model, with separate models independently fit for each value of .
Hyperparameters.
The local regularizer and the weight matrix typically contain some positive hyperparameters, for example that scale the local regularization or one or more edges weights in the graph. As usual, these are varied over a range of values, and for each value a stratified model is found, and tested on a separate validation data set using an appropriate true objective. We choose values of the hyperparameters that give good validation set performance; finally, we test this model on a new test data set.
Modeling with no data at nodes.
Suppose that at a particular node on the graph (i.e., value of stratification feature), there is no data to fit a model. Then, the model for that node will simply be a weighted average of the models of its neighboring nodes.
3.1 Data models
So far we have been vague about what exactly a data model is. In this section we list just a few of the many possibilities.
3.1.1 Point estimates
The first type of data model that we consider is one where we wish to predict a single likely value, or point estimate, of for a given and . In other words, we want to construct a function such that .
Regression.
In a regression data model, and , with one component of always having the value one. The loss has the form , where is a penalty function. Some common (convex) choices include the square penalty [1, 2], absolute value penalty [60], Huber penalty [61], and tilted absolute value (for quantile regression). Common generic regularizers include zero, sum of squares, and the norm [3]; these typically do not include the coefficient in associated with the constant feature value one, and are scaled by a positive hyperparameter. Various combinations of these choices lead to ordinary least squares regression, ridge regression, the lasso, Huber regression, and quantile regression. Multiple regularizers can also be employed, each with a different hyperparameter, as in the elastic net. Constraints on model parameters (for example, nonnegativity) can be imposed by infinite values of the regularizer.
The stratified regression model gives the predictor of , given and , given by
The stratified model uses a different set of regression coefficients for each value of . For example, we can have a stratified lasso model, or a stratified nonnegative least squares model.
An interesting extension of the regression model, sometimes called multivariate or multitask regression, takes with . In this case we take , so is a matrix, with as our predictor for . An interesting regularizer for multitask regression is the nuclear norm, which is the sum of singular values, scaled by a hyperparameter [62]. This encourages the matrices to have low rank.
Boolean classification.
In a boolean classification data model, and . The loss has the form , where is a penalty function. Common (convex) choices of are (square loss); (hinge loss), and (logistic loss). The predictor is
The same regularizers mentioned in regression can be used. Various choices of these loss functions and regularizers lead to least squares classification, support vector machines (SVM) [63], logistic regression [4, 64], and so on [6]. For example, a stratified SVM model develops an SVM model for each value of ; the Laplacian regularization encourages the parameters associated with neighboring values of to be close.
Multiclass classification.
Point estimates without .
As mentioned above, can be absent. In this case, we predict given only . We can also interpret this as a point estimate with an that is always constant, e.g., . In a regression problem without , the parameter is a scalar, and corresponds simply to a prediction of the value . The loss has the form , where is any of the penalty functions mentioned in the regression section. If, for example, were the square (absolute value) penalty, then would correspond to the average (median) of the that take on the particular stratification feature , or . In a boolean classification problem, , and the predictor is , the most likely class. The loss has the form , where is any of the penalty functions mentioned in the boolean classification section. A similar predictor can be derived for the case of multiclass classification.
3.1.2 Conditional distribution estimates
A more sophisticated data model predicts the conditional probability distribution of given (rather than a specific value of given , in a point estimate). We parameterize this conditional probability distribution by a vector , which we denote as . The data model could be any generalized linear model [68]; here we describe a few common ones.
Logistic regression.
Here and , with one component of always having the value one. The conditional probability has the form
(So .) The loss, often called the logistic loss, has the form , and corresponds to the negative loglikelihood of the data record , when the parameter vector is . The same regularizers mentioned in regression can be used. Validating a logistic regression data model using the average logistic loss on a test set is the same (up to a constant factor) as the average log probability of the observed values in the test data, under the predicted conditional distributions.
Multinomial logistic regression.
Here and is finite, say with one component of equal to one. The parameter vector is a matrix . The conditional probability has the form
The loss has the form
and corresponds to the negative loglikelihood of the data record .
Exponential regression.
Here and (i.e., the nonnegative reals). The conditional probability distribution in exponential regression has the form
where is the parameter vector. This is an exponential distribution over , with rate parameter given by . The loss has the form
which corresponds to the negative loglikelihood of the data record under the exponential distribution with parameter . Exponential regression is useful when the outcome is positive. A similar conditional probability distribution can be derived for the other distributions, e.g., the Poisson distribution.
3.1.3 Distribution estimates
A distribution estimate is an estimate of the probability distribution of when is absent, denoted . This data model can also be interpreted as a conditional distribution estimate when . The stratified distribution model consists of a separate distribution for , for each value of the stratification parameter . The Laplacian regularization in this case encourages closeness between distributions that have similar stratification parameters, measured by the distance between their model parameters.
Gaussian distribution.
Here , and we fit a density to the observed values of , given . For example can parametrize a Gaussian on , . The standard parametrization uses the parameters , with ( denotes the set of positive definite matrices.) The probability density function of is
The loss function (in the standard parametrization) is
where . This loss function is jointly convex in and ; the first three terms are evidently convex and the fourth is the matrix fractional function (see [10, p76]). Some common (convex) choices for the local regularization function include the trace of (encourages the harmonic mean of the eigenvalues of , and hence volume, to be large), sum of squares of (shrinks the overall conditional dependence between the variables), norm of (encourages conditional independence between the variables) [69], or a prior on or .
When , this model corresponds to the standard normal distribution, with mean and variance . A stratified Gaussian distribution model has a separate mean and covariance of for each value of the stratification parameter . This model can be validated by evaluating the average log density of the observed values , under the predicted distributions, over a test set of data.
Bernoulli distribution.
Here . The probability density function of a Bernoulli distribution has the form
where . The loss function has the form
which corresponds to the negative loglikelihood of the Bernoulli distribution.
Poisson distribution.
Here (i.e., the nonnegative integers). The probability density function of a Poisson distribution has the form
where . If , then . The loss function has the form
which corresponds to the negative loglikelihood of the Poisson distribution. A similar data model can be derived for the exponential distribution.
Nonparametric discrete distribution.
Here is finite, say where . We fit a nonparametric discrete distribution to , which has the form
where , i.e., the probability simplex. The loss function has the form
A common (convex) regularization function is the negative entropy, given by .
Exponential families.
A probability distribution that can be expressed as
where is a concave function of , is an exponential family [70, 71]. Some important special cases include the Bernoulli, multinomial, Gaussian, and Poisson distributions. A probability distribution that has this form naturally leads to the following (convex) loss function
3.2 Regularization graphs
In this section we outline some common regularization graphs, which are undirected weighted graphs where each vertex corresponds to a possible value of . We also detail some possible uses of each type of regularization graph in stratified model fitting.
Complete graph.
A complete graph is a (fully connected) graph that contains every possible edge. In this case, all of the models for each stratification feature are encouraged to be close to each other.
Star graph.
A star graph has one vertex with edges to every other vertex. The vertex that is connected to every other vertex is sometimes called the internal vertex. In a stratified model with a star regularization graph, the parameters of all of the noninternal vertices are encouraged to be similar to the parameter of the internal vertex. We refer to the internal vertex as the common model. A common use of a star graph is when the stratification feature relates many vertices only to the internal vertex.
It is possible to have no data associated with the common model or internal vertex. In this case the common model parameter is a weighted average of other model parameters; it serves to pull the other parameter values together.
It is a simple exercise to show that a stratified model with the complete graph can also be found using a star graph, with a central or internal vertex that is connected to all others, but has no data. (That is, it corresponds to a new fictitious value of the stratification feature, that does not occur in the real data.)
Path graph.
A path graph, or linear/chain graph, is a graph whose vertices can be listed in order, with edges between adjacent vertices in that order. The first and last vertices only have one edge, whereas the other vertices have two edges. Path graphs are a natural choice for when the stratification feature corresponds to time, or location in one dimension. When the stratification feature corresponds to time, a stratified model with a path regularization graph correspond to a timevarying model, where the model varies smoothly with time.
Cycle graph.
A cycle graph or circular graph is a graph where the vertices are connected in a closed chain. Every vertex in a cycle graph has two edges. A common use of a cycle graph is when the stratification feature corresponds to a periodic variables, e.g., the day of the week; in this case, we fit a separate model for each day of the week, and the model parameter for Sunday is close to both the model parameter for Monday and the model parameter for Saturday. Other examples include days of the year, the season, and nontime variables such as angle.
Tree graph.
A tree graph is a graph where any two vertices are connected by exactly one path (which may consist of multiple edges). Tree graphs can be useful for stratification features that naturally have hierarchical structure, e.g., the location in a corporate hierarchy, or in finance, hierarchical classification of individual stocks into sectors, industries, and subindustries. As in the star graph, which is a special case of a tree graph, internal vertices need not have any data associated with them, i.e., the data is only at the leaves of the tree.
Grid graph.
A grid graph is a graph where the vertices correspond to points with integer coordinates, and two vertices are connected if their maximum distance in each coordinate is less than or equal to one. For example, a path graph is a onedimensional grid graph.
Grid graphs are useful for when the stratification feature corresponds to locations. We discretize the locations into bins, and connect adjacent bins in each dimension to create a grid graph. A spacestratified model with grid graph regularization is encouraged to have model parameters that vary smoothly through physical space.
Entity graph.
There can be a vertex for every value of some entity, e.g., a person. The entity graph has an edge between two entities if they are perceived to be similar. For example, two people could be connected in a (friendship) entity graph if they were friends on a social networking site. In this case, we would fit a separate model for each person, and encourage friends to have similar models. We can have multiple relations among the entities, each one associated with its own graph; these would typically be weighted by separate hyperparameters.
Products of graphs.
As mentioned earlier, can include several original features, and , the number of possible values of , is the product of the number of values that these features take on. Each original stratification feature can have its own regularization graph, with the resulting regularization graph for given by the (weighted) Cartesian product of the graphs of each of the features. For example, the product of two path graphs is a grid graph, with the horizontal and vertical edges having different weights, associated with the two original stratification features.
4 Distributed method for stratified model fitting
In this section we describe a distributed algorithm for solving the fitting problem (3). The method alternates between computing (in parallel) the proximal operator of the local loss and local regularizer for each model, and computing the proximal operator of the Laplacian regularization term, which requires the solution of a number of regularized Laplacian linear systems (which can also be done in parallel).
To derive the algorithm, we first express (3) in the equivalent form
minimize  (4)  
subject to 
where we have introduced two additional optimization variables and . The augmented Lagrangian of (4) has the form
where and are the (scaled) dual variables associated with the two constraints in (4), respectively, and is the penalty parameter. The ADMM algorithm (in scaled dual form) for the splitting and consists of the iterations
Since the problem is convex, the iterates , , and are guaranteed to converge to each other and to a primal optimal point of (3) [11].
This algorithm can be greatly simplified (and parallelized) with a few simple observations. Our first observation is that the first step in ADMM can be expressed as
where is the proximal operator of the function [72] (see the Appendix for the definition and some examples of proximal operators). This means that we can compute and at the same time, since they do not depend on each other. Also, we can compute (and ) in parallel.
Our second observation is that the second step in ADMM can be expressed as the solution to the regularized Laplacian systems
(5) 
These systems can be solved in parallel. Many efficient methods for solving these systems have been developed; we find that the conjugate gradient (CG) method [73], with a diagonal preconditioner, can efficiently and reliably solve these systems. (We can also warmstart CG with .) Combining these observations leads to Algorithm 4.

Algorithm 4.1 Distributed method for fitting stratified models with Laplacian regularization. given Loss functions , local regularization function , graph Laplacian matrix , and penalty parameter . Initialize. . repeat in parallel 1. Evaluate proximal operator of . 2. Evaluate proximal operator of . 3. Solve the regularized Laplacian systems (5) in parallel. 4. Update the dual variables. until convergence
Evaluating the proximal operator of .
Evaluating the proximal operator of corresponds to solving a fitting problem with sumofsquares regularization. In some cases, this has a closedform expression, and in others, it requires solving a small convex optimization problem. When the loss is differentiable, we can use, for example, the LBFGS algorithm to solve it [74]. (We can also warm start with .) We can also reuse factorizations of matrices (e.g., Hessians) used in previous evaluations of the proximal operator.
Evaluating the proximal operator of .
The proximal operator of often has a closedform expression. For example, if , then the proximal operator of corresponds to the projection onto . If is the sum of squares function and , then the proximal operator of corresponds to the (linear) shrinkage operator. If is the norm and , then the proximal operator of corresponds to softthresholding, which can be performed in parallel.
Stopping criterion.
The primal and dual residuals
converge to zero [11]. This suggests the stopping criterion
where and are given by
for some absolute tolerance and relative tolerance .
Selecting the penalty parameter.
Algorithm 4 will converge regardless of the choice of the penalty parameter . However, the choice of the penalty parameter can affect the speed of convergence. We adopt the simple adaptive scheme [75, 76]
where , , and are parameters. (We recommend , , and .) We found that this simple scheme with worked very well across all of our experiments. When , we must rescale the dual variables, since we are using with scaled dual variables. The rescaling is given by
where .
Regularization path via warmstart.
Our algorithm supports warm starting by choosing the initial point as an estimate of the solution, for example, the solution of a closely related problem (e.g., a problem with slightly varying hyperparameters.)
Software implementation.
We provide an (easily extensible) implementation of the ideas described in the paper,
available at www.github.com/cvxgrp/cvxstrat.
We use numpy
for dense matrix representation and operations [77],
scipy
for sparse matrix operations and statistics functions [78],
networkx
for graph factory functions and
Laplacian matrix computation [79],
torch
for LBFGS and GPU computation [80],
and multiprocessing
for parallel processing.
We provide implementations of a number of common stratified models, all of which support the method
Here , , and are numpy
matrices with the data and stratification features
and is a networkx
graph with nodes corresponding to describing the Laplacian.
When model
is a distribution estimate, its signature is model.fit(Y,Z,G)
.
We also provide methods that automatically compute cross validation scores.
We note that our implementation is for the most part expository and for experimentation with stratified models; its speed could evidently be improved by using a compiled language (e.g., C or C++). Practictioners looking to use these methods at a large scale should develop specialized codes for their particular application using the ideas that we have described.
5 Examples
In this section we illustrate the effectiveness of stratified models by combining base fitting methods and regularization graphs to create stratified models. We note that these examples are all highly simplified in order to illustrate their usage; better models can be devised using the very same techniques illustrated. In each example, we fit three models: a stratified model without Laplacian regularization (which we refer to as a separate model), a common model without stratification, and a stratified model with handpicked edge weights. (In practice the edge weights should be selected using a validation procedure.) We find that the stratified model significantly outperforms the other two methods in each example.
The code is available online at https://github.com/cvxgrp/cvxstrat. All numerical experiments were performed on an unloaded Intel i78700K CPU.
5.1 Mesothelioma classification
We consider the problem of predicting whether a patient has mesothelioma, a form of cancer, given their sex, age, and other medical features that were gathered during a series of patient encounters and laboratory studies.
Dataset.
We obtained data describing patients from the Dicle University Faculty of Medicine [81, 82]. The dataset is comprised of males and females between (and including) the ages of and , with (29.6%) of the patients diagnosed with mesothelioma. The 32 medical features in this dataset include: city, asbestos exposure, type of MM, duration of asbestos exposure, diagnosis method, keep side, cytology, duration of symptoms, dyspnoea, ache on chest, weakness, habit of cigarette, performance status, white blood cell count (WBC), hemoglobin (HGB), platelet count (PLT), sedimentation, blood lactic dehydrogenise (LDH), alkaline phosphatise (ALP), total protein, albumin, glucose, pleural lactic dehydrogenise, pleural protein, pleural albumin, pleural glucose, dead or not, pleural effusion, pleural thickness on tomography, pleural level of acidity (pH), Creactive protein (CRP), class of diagnosis (whether or not the patient has mesothelioma). We randomly split the data into a training set containing 90% of the records, and a test set containing the remaining 10%.
Data records.
The ‘diagnosis method’ feature is perfectly correlated with the output, so we removed it (this was not done in many other studies on this dataset which led to nearperfect classifiers). We performed rudimentary feature engineering on the raw medical features to derive a feature vector . The outcomes denote whether or not the patient has mesothelioma, with meaning the patient has mesothelioma. The stratification feature is a tuple consisting of the patient’s sex and age; for example, corresponds to a 62 year old male. Here the number of stratification features .
Data model.
We model the conditional probability of contracting mesothelioma given the features using logistic regression, as described in §3.1.2. We employ sum of squares regularization on the model parameters, with weight .
Regularization graph.
We take the Cartesian product of two regularization graphs:

Sex. The regularization graph has one edge between male and female, with edge weight .

Age. The regularization graph is a path graph between ages, with edge weight .
Model  Test ANLL  Test error 

Separate  
Stratified  0.76  0.27 
Common 
Results.
We used for all experiments, and and for the stratified model. We compare the average negative log likelihood (ANLL) and prediction error on the test set between all three models in table 1. The stratified model performs slightly better at predicting the presence of mesothelioma than the common model, and much better than the separate model. Figure 1 displays the stratified model glucose parameter over age and sex.
5.2 House price prediction
We consider the problem of predicting the logarithm of a house’s sale price based on its geographical location (given as latitude and longitude), and various features describing the house.
Data records.
We gathered a dataset of sales data for homes in King County, WA from May 2014 to May 2015. Each record includes the latitude/longitude of the sale and features: the number of bedrooms, number of bathrooms, number of floors, waterfront (binary), condition (15), grade (113), year the house was built (19002015), square footage of living space, square footage of the lot, and a constant feature (1). We randomly split the dataset into 16197 training examples and 5399 test examples, and we standardized the features so that they have zero mean and unit variance.
Data model.
The data model here is ordinary regression with square loss. We use the sum of squares local regularization function (excluding the intercept) with regularization weight .
Regularization graph.
We binned latitude/longitude into equally sized bins. The regularization graph here is a grid graph with edge weight .
Model  No. parameters  Test RMSE 

Separate  22500  1.514 
Stratified  22500  0.181 
Common  10  0.316 
Random forest  985888  0.184 
Results.
In all experiments, we used . (In practice this should be determined using crossvalidation.) We compared a stratified model () to a separate model, a common model, and a random forest regressor with 50 trees (widely considered to be the best “out of the box” model). We gave the latitude and longitude as raw features to the random forest model. In table 2, we show the size of these models (for the random forest this is the total number of nodes) and their test RMSE. The Laplacian regularized stratified model performs better than the other models, and in particular, outperforms the random forest, which has almost two orders of magnitude more parameters and is much less interpretable.
In figure 2, we show the model parameters for each feature across location. In this visualization, the brightness of a particular location can be interpreted as the influence that increasing that feature would have on log sales price. (A dark spot means that an increase in that feature in fact decreases log sales price.) For example, the waterfront feature seems to have a positive effect on sales price only in particular coastal areas.
5.3 Senate elections
We model the probability that a United States Senate election in a particular state and election year is won by the Democratic party.
Dataset.
We obtained data describing the outcome of every United States Senate election from 1976 to 2016 (every two years, 21 time periods) for all 50 states [83]. At most one Senator is elected per state per election, and approximately 2/3 of the states elect a new Senator each election. We created a training dataset consisting of the outcomes of every Senate election from 1976 to 2012, and a test dataset using 2014 and 2016.
Data records.
In this problem, there is no . The outcome is , with meaning the candidate from the Democratic party won. The stratification feature is a tuple consisting of the state (there are 50) and election year (there are 21); for example, corresponds to the 1994 election in Georgia. Here the number of stratification features . There are 639 training records and 68 test records.
Data model.
Our model is a simple Bernoulli model of the probability that the election goes to a Democrat (see §3.1.3). For each state and election year, we have a single Bernoulli parameter that can be interpreted as the probability that state will elect a candidate from the Democratic party. To be sure that our model never assigns zero likelihood, we let and , where is some small positive constant; we used . Only approximately 2/3 of states hold Senate elections each election year, so our model makes a prediction of how a Senate race might unfold in years when a particular states’ seats are not up for reelection.
Regularization graph.
We take the Cartesian product of two regularization graphs:

State location. The regularization graph is defined by a graph that has an edge with edge weight between two states if they share a border. (To keep the graph connected, we assume that Alaska borders Washington and that Hawaii borders California. There are 109 edges in this graph.)

Year. The regularization graph is a path graph, with edge weight . (There are 20 edges in this graph.)
Results.
We used and . Table 3 shows the train and test ANLL of the three models on the training and test sets. We see that the stratified model outperforms the other two models. Figure 3 shows a heatmap of the estimated Bernoulli parameter in the stratified model for each state and election year. The states are sorted according to the Fiedler eigenvector (the eigenvector corresponding to the smallest nonzero eigenvalue) of the Laplacian matrix of the state location regularization graph, which groups states with similar parameters near each other. High model parameters correspond to blue (Democrat) and low model parameters correspond to red (Republican and other parties). We additionally note that the model parameters for the entire test years 2014 and 2016 were estimated using no data.
Model  Train ANLL  Test ANLL 

Separate  0.69  1.00 
Stratified  0.48  0.61 
Common  0.69  0.70 
5.4 Chicago crime prediction
We consider the problem of predicting the number of crimes that will occur at a given location and time in the greater Chicago area in the United States.
Dataset.
We downloaded a dataset of crime records from the greater Chicago area, collected by the Chicago police department, which include the time, location, and type of crime [84]. There are 35 types of recorded crimes, ranging from theft and battery to arson. (For our purposes, we ignore the type of the crime.) From the dataset, we created a training set, composed of the recorded crimes in 2017, and a test set, composed of the recorded crimes in 2018. The training set has records and the test set has records.
Data records.
We binned latitude/longitude into equally sized bins. Here there are no features , the outcome is the number of crimes, and the stratification feature is a tuple consisting of location bin, week of the year, day of the week, and hour of day. For example, could correspond to a latitude between 41.796 and 41.811, a longitude between 87.68 and 87.667, on the first week of January on a Tuesday between 1AM and 2AM. Here the number of stratification features .
Creating the dataset.
Recall that there are recorded crimes in our training data set. However, the fact that there were no recorded crimes in a time period is itself a data point. We count the number of recorded crimes for each location bin, day, and hour in 2017 (which could be zero), and assign that number as the single data point for that stratification feature. Therefore, we have training data points, one for each value of the stratification feature. We do the same for the test set.
Data model.
We model the distribution of the number of crimes using a Poisson distribution, as described in §3.1.3. To be sure that our model never assigns zero log likelihood, we let and , where is some small constant; we used .
Regularization graph.
We take the Cartesian product of three regularization graphs:

Latitude/longitude bin. The regularization graph is a twodimensional grid graph, with edge weight .

Week of the year. The regularization graph is a cycle graph, with edge weight .

Day of the week. The regularization graph is a cycle graph, with edge weight .

Hour of day. The regularization graph is a cycle graph, with edge weight .
The Laplacian matrix has over 37 million nonzero entries and the hyperparameters are , , , and .
Model  Train ANLL  Test ANLL 

Separate  0.068  0.740 
Stratified  0.221  0.234 
Common  0.279  0.278 
Results.
We ran the fitting method with
The method converged to tolerances of and in 464 iterations, which took about 434 seconds. For each of the three models, we calculated the average negative loglikelihood (ANLL) of the training and test set (see table 4). We also provide visualizations of the parameters of the fitted stratified model. Figure 4 shows the rate of crime, according to each model, in each location (averaged over time), figure 5 shows the rate of crime, according to each model, for each week of the year (averaged over location and day), and figure 6 shows the rate of crime, according to each model, for each hour of the day (averaged over location and week). (This figure shows oscillations in crime rate with a period around 8 hours. We are not sure what this is, but suspect it may have to do with work shifts.)
6 Extensions and variations
In §3 and §4 we introduced the general idea of stratified model fitting and a distributed fitting method. In this section, we discuss possible extensions and variations to the fitting and the solution method.
Varied loss and regularizers.
In our formulation of stratified model fitting, the local loss and regularizer in the objective is a sum of the local losses and regularizers, which are all the same (save for problem data). With slight alterations, the ideas described in this paper can extend to the case where the model parameter for each stratification feature has different local loss and regularization functions.
Nonconvex base fitting methods.
We now explore the idea of fitting stratified models with nonconvex local loss and regularization functions. The algorithm is exactly the same, except that we replace the and update with approximate minimization. ADMM no longer has the guarantees of converging to an optimal point (or even converging at all), so this method must be viewed as a local optimization method for stratified model fitting, and its performance will depend heavily on the parameter initialization. Despite the lack of convergence guarantees, it can be very effective in practice, given a good approximate minimization method.
Different sized parameter vectors.
It is possible to handle that have different sizes with minor modifications to the ideas that we have described. This can happen, e.g., when a feature is not applicable to one category of the stratification. A simple example is when one is stratifying on sex and one of the features corresponds to number of pregnancies. (It would not make sense to assign a number of pregnancies to men.) In this case, the Laplacian regularization can be modified so that only a subset of the entries of the model parameters are compared.
Laplacian eigenvector expansions.
If the Laplacian regularization term is small, meaning that if varies smoothly over the graph, then it can be well approximated by a linear combination of a modest number, say , of the Laplacian eigenvectors. The smallest eigenvalue is always zero, with associated eigenvector ; this corresponds to the common model, i.e., all are the same. The other eigenvectors can be computed from the graph, and we can change variables so that the coefficients of the expansion are the variables. This technique can potentially drastically reduce the number of variables in the problem.
Coordinate descent Laplacian solver.
It is possible to solve a regularized Laplacian system where by applying randomized coordinate descent [85] to the optimization problem
which has solution . The algorithm picks a random row , then minimizes the objective over , to give
(6) 
and repeats this process until convergence. We observe that to compute the update for , a vertex only needs the values of at its neighbors (since is only nonzero for connected vertices). This algorithm is guaranteed to converge under very general conditions [86].
Decentralized implementation.
Algorithm 4 can be implemented in a decentralized manner, i.e., where each vertex only has access to its own parameters, dual variables, and edge weights, and can only communicate these quantities to adjacent vertices in the graph. Line 1, line 2, and line 4 can be obviously be performed in parallel at each vertex of the graph, using only local information. Line 3 requires coordination across the graph, and can be implemented using a decentralized version of the coordinate descent Laplacian solver described above.
The full decentralized implementation of the coordinate descent method goes as follows. Each vertex broadcasts its value to adjacent vertices, and starts a random length timer. Once a vertex’s timer ends, it stops broadcasting and collects the values of the adjacent vertices (if one of them is not broadcasting, then it restarts the timer), then computes the update (6). It then starts another random length timer and begins broadcasting again. The vertices agree beforehand on a periodic interval (e.g., every 15 minutes) to complete line 3 and continue with the algorithm.
Acknowledgments
Shane Barratt is supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE1656518. The authors thank Trevor Hastie and Robert Tibshirani for helpful comments on an early draft of this paper.
References
 [1] A. Legendre. Nouvelles méthodes pour la détermination des orbites des comètes. 1805.
 [2] C. Gauss. Theoria motus corporum coelestium in sectionibus conicis solem ambientium, volume 7. Perthes et Besser, 1809.
 [3] R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society, 58(1):267–288, 1996.
 [4] D. R. Cox. The regression analysis of binary sequences. Journal of the Royal Statistical Society, 20(2):215–242, 1958.
 [5] B. Boser, I. Guyon, and V. Vapnik. A training algorithm for optimal margin classifiers. In Proceedings of the workshop on computational learning theory, pages 144–152. ACM, 1992.
 [6] T. Hastie, R. Tibshirani, and J. Friedman. Elements of statistical learning, 2009.
 [7] L. Breiman, J. Friedman, C.J. Stone, and R.A. Olshen. Classification and Regression Trees. The Wadsworth and BrooksCole statisticsprobability series. Taylor & Francis, 1984.
 [8] I. Goodfellow, Y. Bengio, and A. Courville. Deep Learning. MIT Press, 2016.
 [9] W. Tobler. A computer movie simulating urban growth in the Detroit region. Economic geography, 46(sup1):234–240, 1970.
 [10] S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, 2004.
 [11] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundation and Trends in Machine Learning, 3(1):1–122, 2011.
 [12] W. Kernan, C. Viscoli, R. Makuch, L. Brass, and R. Horwitz. Stratified randomization for clinical trials. Journal of clinical epidemiology, 52(1):19–26, 1999.
 [13] T. L. Lash K. J. Rothman, S. Greenland. Modern Epidemiology. Lippincott Williams & Wilkins, 3 edition, 1986.
 [14] B. Kestenbaum. Epidemiology and biostatistics: an introduction to clinical research. Springer Science & Business Media, 2009.
 [15] L. Jacob and J. Vert. Efficient peptide–MHCI binding prediction for alleles with few known binders. Bioinformatics, 24(3):358–366, 2007.
 [16] S. M. Gross and R. Tibshirani. Data shared lasso: A novel tool to discover uplift. Computational Statistics & Data Analysis, 101, 03 2016.
 [17] R. Tibshirani and J. Friedman. A Pliable Lasso. arXiv eprints, Jan 2018.
 [18] W. Du and R. Tibshirani. A pliable lasso for the Cox model. arXiv eprints, Jul 2018.
 [19] D. Hallac, J. Leskovec, and S. Boyd. Network lasso: Clustering and optimization in large graphs. In Proceedings of the ACM International Conference on Knowledge Discovery and Data Mining, pages 387–396. ACM, 2015.
 [20] D. Hallac, Y. Park, S. Boyd, and J. Leskovec. Network inference via the timevarying graphical lasso. In Proceedings of the ACM International Conference on Knowledge Discovery and Data Mining, pages 205–213. ACM, 2017.
 [21] D. Hallac, C. Wong, S. Diamond, A. Sharang, R. Sosic, S. Boyd, and J. Leskovec. SnapVX: A networkbased convex optimization solver. Journal of Machine Learning Research, 18(1):110–114, 2017.
 [22] J. Leskovec and R. Sosič. Snap: A generalpurpose network analysis and graphmining library. ACM Transactions on Intelligent Systems and Technology (TIST), 8(1):1, 2016.
 [23] S. Diamond and S. Boyd. CVXPY: A Pythonembedded modeling language for convex optimization. Journal of Machine Learning Research, 17(1):2909–2913, 2016.
 [24] T. Hastie and R. Tibshirani. Varyingcoefficient models. Journal of the Royal Statistical Society. Series B (Methodological), 55(4):757–796, 1993.
 [25] J. Fan and W. Zhang. Statistical methods with varying coefficient models. Statistics and Its Interface, 1:179–195, 02 2008.
 [26] C. Brunsdon, A. Fotheringham, and M. Charlton. Geographically weighted regression: a method for exploring spatial nonstationarity. Geographical analysis, 28(4):281–298, 1996.
 [27] D. McMillen. Geographically weighted regression: the analysis of spatially varying relationships, 2004.
 [28] M. Yuan and Y. Lin. Model selection and estimation in regression with grouped variables. Journal of the Royal Statistical Society, 68:49–67, 2006.
 [29] J. Friedman, T. Hastie, and R. Tibshirani. A note on the group lasso and a sparse group lasso. arXiv eprints, 2010.
 [30] L. Meier, S. Van De Geer, and P. Bühlmann. The group lasso for logistic regression. Journal of the Royal Statistical Society, 70(1):53–71, 2008.
 [31] M. Belkin, P. Niyogi, and V. Sindhwani. Manifold regularization: A geometric framework for learning from labeled and unlabeled examples. Journal of Machine Learning Research, 7(Nov):2399–2434, 2006.
 [32] X. Zhu, Z. Ghahramani, and J. Lafferty. Semisupervised learning using Gaussian fields and harmonic functions. In Proceedings of the International Conference on Machine Learning, pages 912–919, 2003.
 [33] B. Nadler, N. Srebro, and X. Zhou. Statistical analysis of semisupervised learning: The limit of infinite unlabelled data. In Y. Bengio, D. Schuurmans, J. D. Lafferty, C. K. I. Williams, and A. Culotta, editors, Advances in Neural Information Processing Systems, pages 1330–1338. Curran Associates, Inc., 2009.
 [34] S. Boyd. Convex optimization of graph Laplacian eigenvalues. In Proceedings International Congress of Mathematicians, pages 1311–1319, 2006.
 [35] J. Chen, C. Wang, Y. Sun, and X. Shen. Semisupervised Laplacian regularized least squares algorithm for localization in wireless sensor networks. Computer Networks, 55(10):2481 – 2491, 2011.
 [36] A. Zouzias and N. M. Freris. Randomized gossip algorithms for solving Laplacian systems. In Proceedings of the European Control Conference, pages 1920–1925, 2015.
 [37] D. Zhang and D. Shen. Semisupervised multimodal classification of alzheimer’s disease. In 2011 IEEE International Symposium on Biomedical Imaging: From Nano to Macro, pages 1628–1631, 2011.
 [38] Z. Wang, Z. Zhou, X. Sun, X. Qian, and L. Sun. Enhanced lapsvm algorithm for face recognition. International Journal of Advancements in Computing Technology, 4:343–351, 2012.
 [39] Z. Wang, X. Sun, L. Zhang, and X. Qian. Document classification based on optimal laprls. Journal of Software, 8(4):1011–1018, 2013.
 [40] F. Wang, Z.A. Huang, X. Chen, Z. Zhu, Z. Wen, J. Zhao, and G.Y. Yan. LRLSHMDA: Laplacian regularized least squares for human microbe–disease association prediction. Scientific Reports, 7(1):7601, 2017.
 [41] A. Kovac and A. D. A. C. Smith. Nonparametric regression on a graph. Journal of Computational and Graphical Statistics, 20(2):432–447, 2011.
 [42] L. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation based noise removal algorithms. Physica D: nonlinear phenomena, 60(14):259–268, 1992.
 [43] R. Caruana. Multitask learning. Machine Learning, 28(1):41–75, 1997.
 [44] Y. Zhang and Q. Yang. A survey on multitask learning. Computing Research Repository, abs/1707.08114, 2017.
 [45] D. Sheldon. Graphical multitask learning. Preprint, 2008.
 [46] D. Spielman. Spectral graph theory and its applications. In Foundations of Computer Science, pages 29–38. IEEE, 2007.
 [47] S. Teng. The Laplacian paradigm: emerging algorithms for massive graphs. In International Conference on Theory and Applications of Models of Computation, pages 2–14. Springer, 2010.
 [48] D. Spielman. Algorithms, graph theory, and linear equations in Laplacian matrices. In Proceedings of the International Congress of Mathematicians, pages 2698–2722. World Scientific, 2010.
 [49] N. Vishnoi. : Laplacian solvers and their algorithmic applications. Foundations and Trends in Theoretical Computer Science, 8(1–2):1–141, 2013.
 [50] M. B. Cohen, R. Kyng, G. L. Miller, J. W. Pachocki, R. Peng, A. B. Rao, and S. C. Xu. Solving sdd linear systems in nearly time. In Proceedings of the Fortysixth Annual ACM Symposium on Theory of Computing, STOC ’14, pages 343–352, New York, NY, USA, 2014. ACM.
 [51] M. Schaub, M. Trefois, P. van Dooren, and J. Delvenne. Sparse matrix factorizations for fast linear solvers with application to Laplacian systems. SIAM Journal on Matrix Analysis and Applications, 38(2):505–529, 2017.
 [52] V. Sadhanala, Y.X. Wang, and R. Tibshirani. Graph sparsification approaches for Laplacian smoothing. In Arthur Gretton and Christian C. Robert, editors, Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, volume 51 of Proceedings of Machine Learning Research, pages 1250–1259, Cadiz, Spain, 09–11 May 2016. PMLR.
 [53] M. Hestenes and E. Stiefel. Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards, 49(6), 1952.
 [54] Y. Sun, P. Babu, and D. Palomar. Majorizationminimization algorithms in signal processing, communications, and machine learning. IEEE Transactions in Signal Processing, 65(3):794–816, 2017.
 [55] J. Tuck, D. Hallac, and S. Boyd. Distributed majorizationminimization for Laplacian regularized problems. IEEE/CAA Journal of Automatica Sinica, 6(1):45–52, January 2019.
 [56] C. Lu, J. Feng, S. Yan, and Z. Lin. A unified alternating direction method of multipliers by majorization minimization. IEEE Transactions on Pattern Analysis and Machine Intelligence, 40(3):527–541, 2018.
 [57] D. R. Hunter and K. Lange. A tutorial on mm algorithms. The American Statistician, 58(1):30–37, 2004.
 [58] A. N. Tikhonov. On the stability of inverse problems. Doklady Akademii Nauk SSSR, 39(5):195–198, 1943.
 [59] A. N. Tikhonov. Solution of incorrectly formulated problems and the regularization method. Soviet Mathematics Doklady, 4:1035–1038, 1963.
 [60] R. Boscovich. De litteraria expeditione per pontificiam ditionem, et synopsis amplioris operis, ac habentur plura ejus ex exemplaria etiam sensorum impessa. Bononiensi Scientiarum et Artum Instuto Atque Academia Commentarii, 4:353–396, 1757.
 [61] P. Huber. Robust estimation of a location parameter. The Annals of Mathematical Statistics, 35(1):73–101, 1964.
 [62] L. Vandenberghe and S. Boyd. Semidefinite programming. SIAM review, 38(1):49–95, 1996.
 [63] C. Cortes and V. Vapnik. Supportvector networks. Machine Learning, 20(3):273–297, 1995.
 [64] D. W. Hosmer and S. Lemeshow. Applied Logistic Regression. John Wiley and Sons, Ltd, 2005.
 [65] J. Engel. Polytomous logistic regression. Statistica Neerlandica, 42(4):233–252, 1988.
 [66] X.Y. Yang, J. Liu, M.Q. Zhang, and K. Niu. A new multiclass svm algorithm based on oneclass svm. In Computational Science – ICCS 2007, pages 677–684, Berlin, Heidelberg, 2007. Springer Berlin Heidelberg.
 [67] C. D. Manning, P. Raghavan, and H. Schütze. Introduction to Information Retrieval. Cambridge University Press, 2008.
 [68] J. Nelder and R. Wedderburn. Generalized linear models. Journal of the Royal Statistical Society, 135(3):370–384, 1972.
 [69] J. Friedman, T. Hastie, and R. Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9(3):432–441, 2008.
 [70] B. O. Koopman. On distributions admitting a sufficient statistic. Transactions of the American Mathematical Society, 39(3):399–409, 1936.
 [71] E. J. G. Pitman. Sufficient statistics and intrinsic accuracy. Mathematical Proceedings of the Cambridge Philosophical Society, 32(4):567â579, 1936.
 [72] N. Parikh and S. Boyd. Proximal algorithms. Foundations and Trends in Optimization, 1(3):127–239, 2014.
 [73] M. R. Hestenes and E. Stiefel. Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards, 49:409–436, 1952.
 [74] D. Liu and J. Nocedal. On the limited memory bfgs method for large scale optimization. Mathematical programming, 45(13):503–528, 1989.
 [75] B. He, H. Yang, and S. Wang. Alternating direction method with selfadaptive penalty parameters for monotone variational inequalities. Journal of Optimization Theory and applications, 106(2):337–356, 2000.
 [76] S. Wang and L. Liao. Decomposition method with a variable parameter for a class of monotone variational inequality problems. Journal of optimization theory and applications, 109(2):415–429, 2001.
 [77] S. Van Der Walt, C. Colbert, and G. Varoquaux. The numpy array: a structure for efficient numerical computation. Computing in Science & Engineering, 13(2):22, 2011.
 [78] E. Jones, T. Oliphant, P. Peterson, et al. SciPy: Open source scientific tools for Python, 2001.
 [79] A. Hagberg, P. Swart, and D. Chult. Exploring network structure, dynamics, and function using NetworkX. Technical report, Los Alamos National Lab.(LANL), Los Alamos, NM (United States), 2008.
 [80] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer. Automatic differentiation in pytorch. In Advances in Neural Information Processing Systems, 2017.
 [81] O. Er, A. C. Tanrikulu, A. Abakay, and F. Temurtas. An approach based on probabilistic neural network for diagnosis of mesothelioma’s disease. Computers and Electrical Engineering, 38(1):75 – 81, 2012. Special issue on New Trends in Signal Processing and Biomedical Engineering.
 [82] D. Dua and C. Graff. UCI machine learning repository. http://archive.ics.uci.edu/ml, 2019.
 [83] MIT Election Data and Science Lab. U.S. House 1976â2016, 2017.
 [84] Chicago Police Department. Crimes  2001 to present. https://data.cityofchicago.org/PublicSafety/Crimes2001topresent/ijzpq8t2, 2019. Accessed: 20190409.
 [85] Y. Nesterov. Efficiency of coordinate descent methods on hugescale optimization problems. SIAM Journal on Optimization, 22(2):341–362, 2012.
 [86] P. Richtárik and M. Takáč. Iteration complexity of randomized blockcoordinate descent methods for minimizing a composite function. Mathematical Programming, 144(1):1–38, 2014.
 [87] D. M. Witten and R. Tibshirani. Covarianceregularized regression and classification for high dimensional problems. Journal of the Royal Statistical Society, 71(3):615–636, 2009.
Appendix A Proximal operators
Recall that for , the proximal operator of a function and parameter is defined as
In this section, we list some functions encountered in this paper and their associated proximal operators.
Indicator function of a convex set.
If is a closed nonempty convex set, and we define as 0 if and otherwise, then
i.e., the projection of onto . (Note: If we allow to be nonconvex, then is still of this form, but the projection is not guaranteed to be unique.) For example, if , then .
Sum of squares.
When , the proximal operator is given by
norm.
When , the proximal operator is given by softthresholding, or
norm.
The proximal operator of is
Negative logarithm.
The proximal operator of is
Sum of and squared.
The proximal operator of for , is
Quadratic.
If , with , , and , then
If, e.g., , then .
Poisson negative log likelihood.
Here , where and . If we let , then
Covariance maximum likelihood.
Here
We want to find . The firstorder optimality condition (along with the implicit constraint that ) is that the gradient should vanish, i.e.,
Rewriting, this is
implying that and share the same eigenvectors [87], stored in . Let denote the th eigenvalue of , and the th eigenvalue of . Then we have
To solve for , we solve these 1dimensional quadratic equations, with the th equation solved by
The proximal operator is then
Bernoulli negative loglikelihood.
Here
with is problem data and . We want to find
By setting the derivative of the above argument with respect to to zero, we find that
Solving this cubic equation (with respect to ) for the solution that lies in yields the result of the proximal operator.