Fast Algorithm for the Lasso based L1-Graph Construction
Yasuhiro Fujiwara†‡∗ , Yasutoshi Ida† , Junya Arai† , Mai Nishimura† , Sotetsu Iwamura††NTT Software Innovation Center, 3-9-11 Midori-cho Musashino-shi, Tokyo, Japan‡NTT Communication Science Laboratories, 2-4 Seika-Cho Soraku-gun, Kyoto, Japan∗Osaka University, 1-5 Yamadaoka, Suita-shi, Osaka, Japan{fujiwara.yasuhiro, ida.yasutoshi, arai.junya, nishimura.mai,
iwamura.sotetsu}@lab.ntt.co.jp

ABSTRACTinteraction, medical analysis, scientific data exploration, andinformation retrieval [10].
Graph is a fundamental data model in finding the hidden structures of high-dimensional data. The nodes andedges of a graph represent data points and the relationships
among them, respectively. Since graphs are useful, various
approaches to exploit graphs have been proposed in many
communities [8, 9, 11, 19]. However, the database community has paid relatively little attention to developing graphconstruction approaches although it is a quite important research problem in real-world applications. Conceptually, agraph should capture the intrinsic complex structures induced by high dimensionality of the data. A primitive approach for graph construction is k-nearest neighbor (k-NN).
In a k-NN graph, each node is connected to its k-nearest
neighbor nodes by exploiting Euclidean distance. However,Euclidean distance is sensitive to data noise, a problem thatis inherent in real-world applications [4].
In order to address this problem, Meinshausen et al. proposed to construct the L1 -graph by solving the L1 -optimization
problem of the lasso [16]. The lasso is a popular least squares
regression approach that represents each data point as thelinear combination of the remaining data points along withcoefficients [13]. By adding an L1 -norm regularization term,
the lasso achieves the sparsity of the solution; it assigns exactly zero coefficients to the most irrelevant data points in
the regression. The idea of Meinshausen et al. is to exploit the lasso to construct graphs; they applied the lasso toeach node, where the regression coefficients are used as edge
weights in the graph. In the approach, the neighborhood relationships and edge weights are simultaneously obtained bysolving the L1 -optimization problem of the lasso. Thus, theresulting graph is referred as the L1 -graph. The L1 -graph isfundamentally different from traditional k-NN graphs. Since
the L1 -graph utilizes the higher order relationships among
data points, it is robust to data noise unlike k-NN graphs[4]. In addition, the L1 -graph has small memory consumption since the lasso can sparsely represent each node by using
the small number of remaining nodes [16].
Even though the lasso-based L1 -graph is used in many
applications, its construction incurs high computation cost
[4, 14, 15]. The graph construction approach proposed by
Meinshausen et al. picks up nodes one by one and computes
the lasso for each node. In terms of processing speed, thecoordinate descent [5] is an efficient approach to the lasso
that iteratively updates edge weights one at a time to compute the regression result. However, it is infeasible to apply
the coordinate descent to large-scale data sets. This is be-

The lasso-based L1 -graph is used in many applications since
it can effectively model a set of data points as a graph. The
lasso is a popular regression approach and the L1 -graph represents data points as nodes by using the regression result.
More specifically, by solving the L1 -optimization problem of
the lasso, the sparse regression coefficients are used to obtain the weights of the edges in the graph. Conventionalgraph structures such as k-NN graph use two steps, adjacency searching and weight selection, for constructing thegraph whereas the lasso-based L1 -graph derives the adjacency structure as well as the edge weights simultaneously byusing a coordinate descent. However, the construction costof the lasso-based L1 -graph is impractical for large data sets
since the coordinate descent iteratively updates the weights
of all edges until convergence. Our proposal, Castnet, canefficiently construct the lasso-based L1 -graph. In order toavoid updating the weights of all edges, we prune edges thatcannot have nonzero weights before entering the iterations.
In addition, we update edge weights only if they are nonzero
in the iterations. Experiments show that Castnet is significantly faster than existing approaches.

1.

INTRODUCTION

We are now living in the big data era. With the rapid development of database technologies, high-dimensional data
has become prevalent in many application domains which
demand techniques to process data in an effective manner[21, 17]. From the data mining perspective, it is important
to detect the hidden structures of high-dimensional datathat can be collectively revealed only by processing large
amounts of data. This is because the hidden structures allow
us to associate, and thus utilize, the semantic information of
data; data belonging to the same cluster share the same semantic information [9]. From a practical perspective, finding
the hidden structures of high-dimensional data is a fundamental process in computational biology, human-computerThis work is licensed under the Creative Commons AttributionNonCommercial-NoDerivatives 4.0 International License. To view a copyof this license, visit http://creativecommons.org/licenses/by-nc-nd/4.0/. For
any use beyond those covered by this license, obtain permission by emailing
info@vldb.org.
Proceedings of the VLDB Endowment, Vol. 10, No. 3
Copyright 2016 VLDB Endowment 2150-8097/16/11.

229

opinions for decision making. In computer-aided diagnosis, analyzing CT images is an essential step in generating
3D models. Although it is conceptually simple that lesions
are just regions whose features differ from those of normal
anatomical structures, it is technically challenging to detect
lesions in CT images due to the poor contrast between thelesion and surrounding tissues, and the high variability of lesion shape. The lasso-based L1 -graph is a robust approachto lesion detection where a voxel in a CT image corresponds
to a node in the graph [15]. The approach identifies lesion
regions in CT images by applying a semi-supervised learning
approach that propagates labels of manually-labeled nodes
to unlabeled nodes in the graph. The approach is motivatedby the superior discriminant power of the sparse representations yielded by the lasso. Due to the sparsity constraint
imposed by the L1 -norm regularization term, the L1 -graphhave a sparse structure where only nodes of similar features
are connected to each other. As a result, the approach caneffectively determine the lesion regions from noisy data. The
approach experimentally shows that it can distinctly identify prostate cancer; one of the leading causes of male cancer
death in the United States.

cause it iteratively computes the weights of all edges untilconvergence even if weights of almost all edges are zero dueto the L1 -norm regularization term of the lasso.

1.1 Problem Statement
If N and M are the numbers of nodes and dimensions,
respectively, a set of high-dimensional data is represented as
matrix X of N × M size. We address the following problem:

Problem (L1 -graph construction).
Given: Matrix X of the high-dimensional data and tuning parameter λ.
Construct: the L1 -graph of N nodes with edge weights
computed by the lasso, efficiently.

In this problem, tuning parameter λ is a positive regularization parameter that balances the trade-off between theloss function and the L1 -norm regularization of the lasso.
Our approach can be used in various applications as shown
below due to the generality of the problem.
Brain Analysis. Recent advances in computational neuroscience allow researchers to use large-scale gene expression
databases in analyzing genes and anatomies. It is crucial
to detect gene-anatomy association for understanding brain
function from molecular and genetic information. The lassobased L1 -graph has been used to model the anatomical organization of brain structure [14]. In the graph, each node
represents the spatial location of a gene, and the edges between nodes encode the correlations between locations in
the brain. Since spatially adjacent regions tend to exhibit
correlated expression patterns, most of the edges connect
to spatially adjacent regions in the graph. However, thegraph revealed an apparent exception; the nodes annotatedas dentate gyrus (DG) were found to be highly correlated tomany other nodes in distant regions of the brain. DG is an
important area for learning and memory to process spatial
information. According to neuroanatomy, DG receives multiple sensory inputs including visual and auditory. Thus,
DG plays an important role in blocking and filtering excitatory activity from the inputs.
Motion Segmentation. In the field of computer vision, motion segmentation is a fundamental process in decomposing
a video into moving objects and a background. Motion segmentation is a pre-processing step in surveillance, tracking,and action recognition. It is important to develop a motion segmentation approach that is robust against noise sincedata is never clean in real world applications. For example,
in human image recognition, the data exhibits different degrees of noise depending on rotation, lighting, scale, and soon. Although many approaches have been proposed in theliterature, the approach based on the affine camera modelis arguably the most popular for motion segmentation. An
important intuition behind the approach is to use motion;
we can easily discern independently moving objects by simply seeing their motion. In the approach, all the trajectories
associated with a single rigid object lie in a linear subspace.
Therefore, motion segmentation can be achieved by clustering trajectories into different motion subspaces [4]. Morespecifically, the approach constructs the L1 -graph from the
trajectory data to apply spectral clustering since the L1 graph is robust against noise unlike other types of graphs.
The approach is more effective than existing approaches.
Lesion Detection. Automatic lesion detection is an important process in early medical screening or providing second

1.2 Contributions
In this paper, we propose Castnet, a novel and efficientapproach to constructing the lasso-based L1 -graph. Before entering the iterations, we prune those edges whose
weights cannot be nonzero. In addition, we efficiently update nonzero weights in each iteration by pruning edges
whose weights are zero. Although our approach skips unnecessary weight updates similarly to the state-of-the-art
approach of the lasso [7], we substantially improve the efficiency compared to the previous approach by effectively
pruning the inner product computations used for graph construction. As a result, our approach is 130 times faster than
the state-of-the-art approach. Even though the lasso-based
L1 -graph was designed to overcome the sensitivity of k-NN
graphs to data noise, it is seldom applied to large-scale data
sets due to the high computation cost. However, our approach can efficiently construct graphs, which will improvethe usefulness of many applications such as brain analysis,
motion segmentation, and lesion detection.
The remainder of this paper is organized as follows: Section 2 describes related work. Section 3 gives an overview
of the background. Section 4 introduces our approach. Section 5 reviews our experiments and the results. Section 6
details case-studies. Section 7 provides our conclusions.

2. RELATED WORK
Meinshausen et al. proposed the lasso-based L1 -graphas more robust approach against noise than traditional kNN graphs [16]. Their approach picks up nodes one by
one and uses the lasso to compute weights for the node.
A popular approach to solving the L1 -optimization problemof the lasso is the coordinate descent [5]. To resolve the L1 optimization problem, it iteratively updates edge weights
one at a time. In terms of the efficiency, Sling is a state-ofthe-art approach based on the coordinate descent although
it is not so useful to construct the L1 -graph [7]. Sling uses
the sequential strong rule to discard unnecessary edges. It
is based on the assumption that a grid of tuning parameters (λ1 > λ2 > λ3 > . . .) is used in solving the lasso the230Table 1: Definitions of main symbols.
SymbolDefinition

N
M
m
λ
p
ep [u]
wp [u]
dp
K[u|wp ]K[u|wp ]K[u|wp ]xp
wp
XV
M[wp ]
C[wp ]

Number of nodes in the graph
Number of dimensions of each data point
Target rank of SVD
Tuning parameter
Node for weight computation
Edge from node u to p
Weight of edge ep [u]Number of edges to node p
KKT condition score of node u for weight vector wp
Upper bound of KKT condition score K[u|wp ]Lower bound of KKT condition score K[u|wp ]p-th row vector of matrix X
1 × N weight vector of node p
N × M matrix of high-dimensional data
Set of nodes in the graphSet of edges that must have nonzero weights for wp
Set of edges that can have nonzero weights for wp

where $ · $1 and $ · $2 are the L1 -norm and L2 -norm of avector, respectively, and λ > 0 is a tuning parameter. Notethat the graph does not have self-loop edges since wp [p] = 0
for each picked up node. In Equation (1), the first term corresponds to the squared loss by regression; the second term
corresponds to an L1 -norm constraint for the coefficients.
Therefore, Equation (1) implements the regularization ofweight vector wp by trading off the accuracy of the regression for a reduction in the sparsity of coefficients; it has theeffect of giving a solution to the objective function with few
edges of nonzero weights [23]. It is clear from Equation (1)that the graph has improved sparsity as we raise tuning parameter λ. As described in the previous paper [6], if somenodes are strongly correlated, a single node would be used torepresent node p; only an edge from the single node is likely
to have nonzero weight in that case. If matrix X has full row
rank, we have a unique solution for the optimization problem of Equation (1) [13]. Otherwise, which is necessarily thecase when N > M , there may not be a unique solution.
In order to efficiently solve the lasso, Friedman et al. proposed the coordinate descent that updates weights of edges
one at a time until they reach convergence [5]. The coordinate descent partially conducts optimization with respect
to weight wp [u] by supposing that it has already estimatedother weights for all nodes u such that u ∈ V\p. Specifically,
the coordinate descent updates weights as follows:same as the other screening techniques [24]; the tuning parameter controls the sparsity of the solution. In each iteration, Sling estimates parameters of the soft-thresholding
operator used in updating weights. In addition, Sling employs the covariance-based approach to efficiently updateweights by computing the inner products of all nodes before
entering each iteration [6]. After convergence, Sling checks
the Karush-Kuhn-Tucker (KKT) condition of all discarded
edges since the sequential strong rule may erroneously discard edges of nonzero weights [24]. Although Sling is fasterthan the original coordinate descent approach for the lasso,
it is not practical to construct the L1 -graph by directly using Sling. This is because the tuning parameter is constantwhen constructing a graph; the sequential strong rule is noteffective in pruning unnecessary edges since λ does not dynamically change when constructing a graph. In addition,
Sling computes the inner products of all nodes even if almost
all edges are zero as a result of the lasso. To the best of ourknowledge, Castnet is the first approach that can efficiently
construct the lasso-based L1 -graph.

3.

wp [u] ← S[zp [u|wp ], λ]where S[·, ·] is the following soft-thresholding operator [5]:S[zp [u|wp ], λ]!
zp [u|wp ] − λ= zp [u|wp ] + λ
0(zp [u|wp ] > 0 and |zp [u|wp ]| > λ)(zp [u|wp ] < 0 and |zp [u|wp ]| > λ)(|zp [u|wp ]| ≤ λ)(3)

In addition, zp [u|wp ], a parameter of node u for weight vector wp , is computed as follows:"M(u)1
zp [u|wp ] = M(4)
i=1 xu [i](xp [i] − x̃p [i])PRELIMINARY(u)where x̃p [i] is the regression result for element xp [i] without
using node u computed as follows:
"
(u)
(5)
x̃p [i] = v∈V\{p,u} wp [v]xv [i]We formally define the notations and introduce the background of this paper. Table 1 lists the main symbols andtheir definitions. In the lasso-based L1 -graph, nodes andedges correspond to data points and their relationships, respectively. Let p be a node in the graph and V be theset of nodes in the graph, the lasso is used to compute theweights of edges by picking up each node such that p ∈ V.
Since the lasso sets zero coefficients to almost all edges, thegraph has a sparse structure. Let X ∈ RN×M be the matrix of N data points that have M dimensional features and
xp = (xp [1], xp [2], . . . , xp [M ]) be the p-th row vector in matrix X, vector xp corresponds to the p-th data point or nodep. We assume the each row vector is centered and normalized [6]; each vector has average and variance of 0 and 1,
respectively. Let wp be a 1 × N weight vector whose uth element wp [u] is the weight of the edge from node u top. In the lasso-based L1 -graph, each node is represented
as the sparse linear superposition of other nodes by solving
the lasso optimization problem. Specifically, we computeweights of edges to node p by minimizing the following objective function of the regression where wp [p] = 0 [16]:1$xp − wp X$22 + λ$wp $1
minwp ∈RN 2M

(2)

(u)

As shown by the above equations, zp [u|wp ] and x̃p [i] are(u)required to update a weight by Equation (2). Since x̃p [i] isdifferent from each weight, above equations indicate that theoriginal coordinate descent approach needs different regression results for each weight. It takes O(M ) time to compute
zp [u|wp ] from Equation (4). In addition, Equation (4) ad(u)ditionally requires O(dp ) time to compute x̃p [i] for eachdimension by Equation (5) if dp is the number of nonzero
weight edges to node p; dp is degree of node p. As a result,if T is the number of updates until convergences, it takesO(dp M T ) time to compute edge weights to node p by using
Equation (2). As a result, the computation cost of construct
the L1 -graph can be prohibitive for large datasets.

4. PROPOSED METHODThis section presents Castnet, which can efficiently construct the lasso-based L1 -graph. Section 4.1 overviews theideas that underlie our approach. Section 4.2 describes ourapproach to determine the edge sets to update weights by

(1)

231

approach of skipping update computations for zero weightedges in Section 4.2.3.

using upper and lower bounds of KKT condition score. In
Section 4.3, we introduce the approach that computes thebounds of KKT condition score. Section 4.4 shows the graphconstruction algorithm along with its theoretical properties.

4.2.1 Definitions
Before entering the iterations, we compute the two edgesets of M[wp ] and C[wp ]. M[wp ] is the set of edges whose
weights must be nonzero when updated with weight vector wp . Similarly, C[wp ] is the set of edges that can havenonzero weights when updated with weight vector wp . Theoretically, edge set M[wp ] and C[wp ] are based on the KKT
condition [24]. We obtain these edge sets by computing theupper and lower bounds of KKT condition scores that are
given in Section 4.3. Before defining the edge sets, we show
how the KKT condition can check whether weight wp [u] hasnonzero score or not for weight vector wp as follows:

4.1 Main Ideas
In constructing the lasso-based L1 -graph, the original coordinate descent approach incurs high computation cost.
This is because (1) it updates the weights of all edges until
the convergence although almost all edges have zero weights
due to the L1 -norm constraint of the lasso and (2) it updates each weight by computing a different regression result
for each weight as described in Section 3.
For efficient graph construction, our approach is to prune
unnecessary edges that cannot have nonzero weights so as todetermine a set of edges for weight updating. More specifically, we first limit weight updates to the set of edges thatmust have nonzero weight until convergence. After that,
we update weights of the set of edges that can have nonzero
weight. As a result, we can efficiently obtain weights for eachnode. Theoretically, we efficiently determine edge sets thatmust/can have nonzero weight by computing the upper andlower bounds of KKT condition scores; KKT condition was
originally used to find edges that are erroneously discarded
by the sequential screening rule as described in Section 2[24]. In addition, we exploit two types of efficient update
computations for weights. The first type of update computation is based on the residual of the lasso that is the samefor each weight. Since we do not need to compute different regression results in each weight update, we can greatly
enhance computation speed comparing to the original coordinate descent approach. We use this type of computation for edges that have zero weight. Otherwise, we update
the weight using the covariance-based approach [6]; this is
the second type of update computation. Our approach canmore efficiently update weights than the previous approacheven though we exploit the same covariance-based approach[7]. This is because we use the residual-based update computation if the weight of an edge is zero. As described in
Section 2, the covariance-based approach requires the computations of inner products for all edges to update weights
before entering the iterations. Since we do not update alledge weights by the covariance-based approach, we can effectively prune the inner product computations unlike theprevious approach. Furthermore, we skip update computations for edges if they do not have nonzero weight. Thisis because the edge sets to update may include zero weightedges. By skipping the update computations for such edges,
we improve the efficiency with which nonzero weights are
updated in the iterations.
Even though we adopt a different strategy from the previous approach in constructing the L1 -graph, we can provably guarantee to output the same graph as the previous
approach if matrix X has full row rank. This is because
our approach converges on the unique lasso solution in thatcase. Therefore, our approach is effective and efficient inconstructing the lasso-based L1 -graph.

Definition 1 (KKT condition [24]). If K[u|wp ] is the
KKT score of node u for weight vector wp , we have1
K[u|wp ] = λ1 wp [u] + λM
(xp − wp X)x$u

(6)

where x$u is the transpose of xu . For all nodes u such thatu ∈ V/p, if the weight of the edge from node u to p is
updated with weight vector wp , K[u|wp ] has the following
scores based on weight wp [u] after the update:(iff wp [u] > 0)
 K[u|wp ] > 1
K[u|wp ] < −1
(iff wp [u] < 0)(7)

−1 ≤ K[u|wp ] ≤ 1 (iff wp [u] = 0)

It is clear that we can obtain a set of edges that havenonzero weights by directly computing the KKT condition
scores from Equation (6). This, however, incurs high computation cost; it takes O(N M ) time for determining each
KKT condition score in the worst case since the sizes of xp ,wp , X, and x$
u in Equation (6) are 1×M , 1×N , N ×M , andM ×1, respectively. To overcome this problem, we efficiently
compute the upper and lower bounds of the KKT condition
scores introduced in Section 4.3. If K[u|wp ] and K[u|wp ]are the upper and lower bounds of KKT condition scoreof K[u|wp ], respectively, we have K[u|wp ] ≤ K[u|wp ] ≤K[u|wp ]. By exploiting the bounds, the two edge sets are
given as follows:

Definition 2 (Edge Set M). For node u such that u (=p, the definition of edge set M[wp ] is given as follows:
M[wp ] = {ep [u] : K[u|wp ] < −1 or K[u|wp ] > 1}

(8)

Definition 3 (Edge Set C). Letting u (= p, the following equation defines edge set C[wp ]:C[wp ] = {ep [u] : K[u|wp ] > 1 or K[u|wp ] < −1}(9)In Equations (8) and (9), ep [u] is the edge from node u top. Therefore, an edge is included in M[wp ] if K[u|wp ] < −1or K[u|wp ] > 1 holds. Similarly, when we have K[u|wp ] > 1
or K[u|wp ] < −1, an edge is placed in edge set C[wp ]. The
KKT condition yields the following properties for edge sets:
Lemma 1 (Edge Set M). Each edge ep [u] such thatep [u] ∈ M[wp ] must have nonzero weight if its weight isupdated with weight vector wp .

4.2 Weight Updates
This section describes the approach to updating weights.
We first introduce the definitions of edge sets as well as theirproperties in Section 4.2.1. We then describe the weight
update computations in Section 4.2.2. We finally show our

Proof
For each edge included in edge set M[wp ], wehave K[u|wp ] < −1 or K[u|wp ] > 1 if its weight is updated for weight vector wp from Definition 2. In the case232

of K[u|wp ] < −1, we have K[u|wp ] ≤ K[u|wp ] < −1 sinceK[u|wp ] is the upper bound of K[u|wp ]. Therefore, weightwp [u] of edge ep [u] must have negative score from the KKT
condition. If K[u|wp ] > 1 holds for edge ep [u], we similarly
have K[u|wp ] ≥ K[u|wp ] > 1 since K[u|wp ] ≥ K[u|wp ]holds. As a result, edge ep [u] must have positive weight inthis case from the KKT condition.
!

Since the residual used in Equation (11) is the same foreach weight, we can efficiently compute parameter zp [u|wp ]if the regression result does not change after the update.
While Equation (4) needs O(dp M ) time, Equation (10) requires O(M ) time to compute parameter zp [u|wp ]. For
Equation (10), we have the following property:
Lemma 4 (Parameter z). Equation (4) and (10) yield
the same result in computing parameter zp [u|wp ].

Lemma 2 (Edge Set C). For each edge ep [u] such thatep [u] ∈ C[wp ], its weight wp [u] can have nonzero score if we
perform the update computation with weight vector wp .
Proof
To prove Lemma 2, we show that the weightof an edge must be zero after being updated with weightvector wp if K[u|wp ] ≤ 1 and K[u|wp ] ≥ −1 hold. In thecase of K[u|wp ] ≤ 1 and K[u|wp ] ≥ −1, we have −1 ≤K[u|wp ] ≤ K[u|wp ] ≤ K[u|wp ] ≤ 1. Therefore, such anedge cannot have nonzero weight from the KKT condition.
In addition, it is clear that such an edge cannot be included
in edge set C[wp ] from Definition 3. As a result, if we haveK[u|wp ] > 1 or K[u|wp ] < −1 for an edge, the edge can!
have nonzero weight and so is included in C[wp ].

Proof From Equation (5) and (11), we have
"
"
x̃(u)
p [i] =v∈V\{p,u} wp [v]xv [i] =v∈V\p wp [v]xv [i]−wp [u]xu [i]= xp [i]−rp [i]−wp [u]xu [i]Therefore, from Equation (4)"M(u)
1zp [u|wp ] = M
i=1 xu [i](xp [i] − x̃p [i])"M1= M
i=1 xu [i](rp [i] + wp [u]xu [i])
"M
"M
2
1
1
= M
i=1 xu [i]rp [i] + M wp [u]i=1 (xu [i])In terms of the relationship between edge set M[wp ] andC[wp ], we have the following property:
Lemma 3 (Edge Set M and C). If an edge is includedin M[wp ], the edge must be included in C[wp ].
ProofFrom Definition 2, we have K[u|wp ] < −1 orK[u|wp ] > 1 for an edge if the edge is included in M[wp ]. If
K[u|wp ] < −1 holds, we have K[u|wp ] ≤ K[u|wp ] < −1 forthe edge. In addition, if we have K[u|wp ] > 1, K[u|wp ] ≥K[u|wp ] > 1 holds for the edge. Therefore, it is clear thatsuch the edge is also included in C[wp ] from Definition 3. !
Lemma 3 indicates that, if an edge is determined to havenonzero weight by M[wp ], the edge is also determined tohave nonzero weight by C[wp ]. Similarly, Lemma 3 indicates
that, if an edge is determined not to have nonzero weightby C[wp ], the edge is also determined not to have nonzero
weight by M[wp ]; an edge cannot have nonzero weight if it
is not included in C[wp ]. We determine M[wp ] and C[wp ]before entering the iterations to efficiently perform update
computations for edges that must/can have nonzero weights.

Since the"variance of vector xu is 1 as described in Section 3,
2
we have M
i=1 (xu [i]) = M . As a result,"M1
zp [u|wp ] = M
i=1 xu [i]rp [i] + wp [u]which completes the proof from Equation (10).

!

We exploit Equation (10) if an edge has zero weight before
the update. Since such an edge is expected to have zeroweight again after the update [5], the residual is likely tohave the same score after the update. Therefore, we caneffectively utilize residual rp [i] again in the next iteration in
updating other weights. Even if the edge has nonzero weightafter the update, we can incrementally update residual rp [i]in O(1) time by exploiting Equation (11).
If an edge has nonzero weight before the update, we use
the covariance-based approach to update weights since it canefficiently compute parameter zp [u|wp ] in O(dp ) time [6]:'&"1
zp [u|wp ] = wp [u]+ M*xp ,xu +− v:|wp [v]|>0 wp [v]*xv ,xu + (12)where *xp , xu + is the inner product of vector xp and xu ,"i.e., *xp , xu + = M
i=1 xp [i]xu [i]. Note that Equation (4) and(12) give the same result. Although we exploit the same
covariance-based approach, our approach has less computation cost than the previous approach [7]. Before updating
weights, the covariance-based approach requires inner product computations of nonzero weights to all edges undergoing update computations. Since the previous approach updates the weights of all the edges by the covariance-based
approach, it needs O(dp N M ) time to compute inner products. On the other hand, we exploit the covariance-based
approach only if the edge has nonzero weight before theupdate. As a result, our approach takes O(d2p M ) time tocompute the inner products. Note that we have dp , N
since the L1 -graph has sparse structure.

4.2.2 Efficient Update Computation
The edge sets introduced in the previous section can, before entering the iterations, reduce the number of edges to
update. This section describes our approach to improving
the speed of a single update computation in the iterations.
As shown in Equation (2), the original coordinate descentapproach updates weights by exploiting soft-thresholding
operator S[zp [u|wp ], λ]. Since parameter zp [u|wp ] requires
different regression results for each weight as described in
Section 3, the original coordinate descent approach incurs
high computation cost. In our approach, we use two typesof efficient update computations. The first type of computation updates weight by computing parameter zp [u|wp ]without demanding different regression results as follows:"M1
zp [u|wp ] = M(10)i=1 xu [i]rp [i] + wp [u]4.2.3 Gradual Edge Addition
As shown in Section 4.2.1, we determine edge set M[wp ]and C[wp ] before entering the iterations. If U is an edge set
to update weights, we can naively reduce the number of updates by setting U = M[wp ] and U = C[wp ] to directly apply
the coordinate descent. However, this approach may iteratively perform update computations for zero weight edges
since we use the upper and lower bounds to obtain M[wp ]where rp [i] is the residual for element xp [i] in the regression
result by the lasso computed as follows:"(11)
rp [i] = xp [i] − v∈V\p wp [v]xv [i]233and C[wp ]; we exploit the approximations to determine M[wp ]and C[wp ]. This section describes the approach that skips
unnecessary update computations to improve the efficiencyof the L1 -graph construction1 .
For efficient graph construction, we add edges one by one
from edge set A to edge set U. In our approach, we firstobtain edge set U as edges of nonzero weight and perform
update computations until convergence for edge set U. After convergence, we compute M[wp ] and obtain set of addingedges as A = M[wp ]\U. Then, we perform an update computation for each edge in A one by one. If an edge hasnonzero weight after the update computation, we add theedge to U and update weights of U until convergence. Otherwise, we do not iteratively perform update computations forthe edge. Similarly, we compute A = C[wp ]\U and updateweights. Since we avoid iterative computations if an edge
does not have nonzero weight, we can improve the efficiency
of weight updates. Theoretically, this approach is based on
the following property of converged weights:
Lemma 5 (Converged Weights). Let wp% be a converged weight vector for edge set U[wp ] + u obtained by thecoordinate descent based on weight vector wp . If (1) weightvector wp reaches convergence by the coordinate descent withedge set U[wp ] and (2) wp [u] = 0 and wp% [u] = 0 hold fornode u, we have wp% = wp .

matrix and X̃ be the transformed matrix of X, matrix X
is represented as X = X̃V by SVD. Similarly, let vector x̃ube the transformed vector of xu , vector xu is represented as
xu = x̃u V. If r̃p = x̃p − wp X̃, the definitions of the upperand lower bounds are given as follows:
Definition 4 (Upper and Lower Bounds). If K[u|wp ]and K[u|wp ] are the upper and lower bounds of K[u|wp ], respectively, K[u|wp ] and K[u|wp ] are given as follows:
(
)
"
2
1
K[u|wp ] = λ1 wp [u]+2λM$r̃p $22 +M − m
i=1 (r̃p [i]− x̃u [i]) (13)
(
)
"
2
1
$r̃p $22 +M − m
K[u|wp ] = λ1 wp [u]−2λM
i=1 (r̃p [i]+ x̃u [i]) (14)In Definition 4, $r̃p $2 is the L2 -norm of vector r̃p , r̃p [i]is the i-th element of r̃p , and m is the target rank of SVD.
The following lemma shows that K[u|wp ] and K[u|wp ] givethe upper and lower bounds, respectively:

Lemma 6 (Upper and Lower Bounds). If we update
the weight of the edge from node u to p for weight vector wp ,we have K[u|wp ] ≥ K[u|wp ] and K[u|wp ] ≤ K[u|wp ] for theupper and lower bounds given in Definition 4.
Proof Since VV$ = I where I is the identity matrix,"M$ $$(xp−wp X)x$u = (x̃p−wp X̃)VV x̃u = r̃p x̃u =i=1 r̃p [i]x̃u [i]Since r̃p [i]x̃u [i] = 12 {(r̃p [i])2 +(x̃u [i])2 −(r̃p [i]−x̃u [i])2 } holds,
"M
2
2
2
1
(xp −wp X)x$u =2i=1 {(r̃p [i]) + (x̃u [i]) − (r̃p [i] − x̃u [i]) }"2
Note that M
i=1 (x̃u [i]) = M holds since SVD is orthogonaltransformation and each vector is normalized as describedin Section 3. Therefore, we have from Equation (6)
(
)
"
2
1
K[u|wp ] = λ1 wp [u] + 2λM$r̃p $22 + M − M
i=1 (r̃p [i] − x̃u [i])Proof Since wp [u] = 0 and wp% [u] = 0, the u-th elementof wp and wp% are obviously the same. If wp% [u] = 0 holds, it
is clear that weight wp% [u] does not have any impact on the
update results from Equation (5). In addition, weight vector
wp has already reached convergence by exploiting edge set
U[wp ]. Therefore, we have wp% [u] = wp [u] for all elements inU[wp ] after the update computations to obtain weight vector
wp% if we use weight vector wp in the coordinate descent. As!
a result, we have wp% = wp .

Since (r̃p [i] − x̃u [i])2 ≥ 0 holds, we have
(
)
"
2
1
K[u|wp ] ≤ λ1 wp [u] + 2λM$r̃p $22 + M − m
i=1 (r̃p [i] − x̃u [i])Lemma 5 indicates that we can skip the update computations for the edge from node u to p as its weight is zero in theiterations without changing the regression results yielded by
the coordinate descent. This skipping of unnecessary computations allows us to efficiently update nonzero weights.
We show the detail algorithm of this approach in Section 4.4
as a part of the graph construction algorithm.

Therefore, we have K[u|wp ] ≥ K[u|wp ] from Equation (13).
Similarly, we have K[u|wp ] ≤ K[u|wp ] since r̃p [i]x̃u [i] =1{−(r̃p [i])2 − (x̃u [i])2 + (r̃p [i] + x̃u [i])2 } holds.
!
2

From the above proof, it is clear that
" errors of the upper2
and lower bounds are proportional to M
i=m+1 (r̃p [i]− x̃u [i])"M2andi=m+1 (r̃p [i] + x̃u [i]) , respectively. Although SVD
gives the smallest squared error in approximating matrix
X, it is not guaranteed to give the tightest upper and lower
bounds since vector r̃p is distorted by vector wp in the formof r̃p = x̃p − wp X̃. However, it is expected that we can
improve the quality of the upper and lower bounds by using
SVD. This is because r̃p [i] and x̃u [i] are likely to take small
absolute values as i increases.
It takes O(m) time to compute the bounds from Equation (13) and (14) if we have norm $r̃p $2 . Note that norm$r̃p $2 is computed in O(dp m) time and used for all edges tonode p once we compute the norm. In addition, we can efficiently compute the SVD of matrix X in O(N M log m) timeby using an existing technique [12]. Note that we compute
SVD only once in constructing the graph.

4.3 Upper and Lower Bounds
As described in Section 4.2, we can prune edges that
do not have nonzero weight by computing KKT condition
scores. However, it needs the high cost of O(N M ) time tocompute the score for each edge because the size of data
matrix X is N × M . This section introduces our approaches
to efficiently compute the upper and lower bounds of KKT
conditions scores. Section 4.3.1 describes our approach tocompute the bounds by using SVD (singular value decomposition). Section 4.3.2 shows the incremental approach forbound computations.

4.3.1 SVD-Based Bound Computations
In this approach, we approximate matrix X by computing SVD. We adopt SVD since it gives high approximation
quality in terms of squared loss [20]. Let V be a unitary

4.3.2 Incremental Computation
As described in Section 4.2.3, we update weights by addingedges one by one to edge set U to update from edge set A
by computing A = M[wp ]\U and A = C[wp ]\U. We can reduce the computation time to obtain edge set M[wp ] and1

Although all the edges included in M[wp ] have nonzero weights after
the update computation with weight vector wp , the weight vector
changes after the updates. Therefore, M[wp ] can include edges whose
weights are zero as a result of update computations.

234

"In Equation (18) and (19), ∆p = κi=1 δpi . If the edge is notincluded in edge set A, the bounds are computed as follows:

C[wp ] by exploiting the SVD-based approach introducedin the previous section. However, the efficiency of the approach would only be moderate since it requires O(m) time
to compute the bounds for each edge every time we computethe edge sets. How can we cut down the computation costfor determining the bounds? This is the motivation behind
our approach of incremental computation; we incrementally
compute the bounds at O(1) time. This approach is basedon the approach that adds edges one by one as introducedin Section 4.2.3; we can effectively update the bounds since
the weight vector is not so different before and after theedge additions. This section first describes the relation between KKT condition score and parameter zp [u|wp ]. Then,it shows the incremental approach to computing the bounds.
The following lemma confirms that we can compute theKKT condition score from parameter zp [u|wp ]:

Proof

(15)

Lemma 8 (Incremental computation). After the νth edge addition, we have K[u|wpν ] ≥ K[u|wpν ] and K[u|wpν ] ≤K[u|wpν ] for the bounds given by Definition 5 if the weights
of edges are updated for weight vector wpν .
Proof From Equation (6) and (16), we have the following equation from the Cauchy-Schwarz inequality [22] forthe i-th edge addition (1 ≤ i ≤ ν):

$
i
i−1
1
(rip −ri−1[u])||K[u|wpi ]−K[u|wpi−1 ]| = λ1| Mp )xu +(wp [u]−wp

2
$ 2
i
i−1
i
i−1 2
i
i−1
1
1
≤ λ1{ M$rip−ri−1
p $2 $xu$2+$wp−wp $1} = λ{$rp−rp $2+$wp−wp $1}

1=M(wp [u]M +*xu , xp +−*xu , wp X+)

where *xp , xu + is the inner product of vector xp and xu .
Therefore, we have
')
(
" &
"
1
λK[u|wp ] = M
wp [u]M+ M
i=1 xu [i]xp [i]−xu [i]v∈V wp [v]xv [i]
&
')
(
"
"
1
=M
wp [u]M+ M
i=1xu [i] xp [i]− v∈V wp [v]xv [i]"2Since M
i=1 (xu [i]) = M and wp [p] = 0 as described in Section 3, we have the following equation from Equation (5):
&
')
("M
"
1
λK[u|wp ] = M
i=1 xu [i] wp [u]xu [i]+xp [i]− v∈V wp [v]xv [i]
&
')
("M
"
1
=M
i=1 xu [i] xp [i]− v∈V\{p,u} wp [v]xv [i]
&
"M
(u) '
1
=M
i=1 xu [i] xp [i]− x̃p [i]Therefore, λK[u|wp ] = zp [u|wp ] holds from Equation (4). !

Lemma 7 indicates that we can compute KKT condition
score in O(1) time if (1) an edge is included in the edge
set and (2) we have parameter zp [u|wp ] of the edge. We
incrementally compute the upper and lower bounds afteredge addition as follows:
Definition 5 (Incremental computation). Let ν be
the number of edges in edge set A that are added to edge
set U and wpi = (wpi [1], wpi [2], . . . , wpi [n]) be the converged
weight vector after the i-th edge addition (1 ≤ i ≤ ν). In
addition, let rip = (rpi [1], rpi [2], . . . , rpi [n]) be a residual vector
that corresponds to weight vector wpi given as follows:(16)

Let δpi be the score that corresponds to the difference of KKT
condition score before and after the i-th edge addition:2ii−1
δpi = λ1 {$rip − ri−1
p $2 + $wp − wp $1 }(21)

respond to the bounds after the edge additions. Therefore,
we can incrementally update the bounds by using zp [u|wpκ ],K[u|wp0 ], and K[u|wp0 ] from Definition 5. We have the following property for the bounds:

1
λK[u|wp ] = wp [u]+ M(xp −wp X)x$u

rip = xp − wpi X(20)

K[u|wpν ] = K[u|wp0 ] − ∆p

In Definition 5, K[u|wp0 ] and K[u|wp0 ] correspond to thebounds before the edge additions; K[u|wpν ] and K[u|wpν ] cor-

Lemma 7 (Equivalence of KKT condition score).
For node u and weight vector wp , we can compute KKT condition score K[u|wp ] from parameter zp [u|wp ] as follows:
K[u|wp ] = λ1 zp [u|wp ]From Equation (6), we haveK[u|wpν ] = K[u|wp0 ] + ∆p

(18)
(19)

1
z [u|wpκ ] − ∆p ≤ K[u|wpν ] ≤ λ1 zp [u|wpκ ] + ∆p
λ p

Similarly, for an edge that is not included in the edge set,K[u|wp0 ] − ∆p ≤ K[u|wpν ] ≤ K[u|wp0 ] + ∆p
which completes the proof.

!

As described in Section 4.2.2, we compute the residuals forthe added edges every time the edges change the regression
result to efficiently update the weights. Since the coordinate
descent updates weights one at a time, we can incrementally
compute δpi and ∆p at O(1) time from Equation (17). As
a result, we can incrementally compute the bounds in O(1)time by exploiting Definition 5.
In terms of computing the upper and lower bounds, theincremental approach of Definition 5 is the same as the SVDbased approach of Definition 4. However, they yield different
bounds since they adopt totally different processes to compute the bounds. In addition, the incremental approach cancompute the bounds more efficiently. Therefore, we first exploit the incremental approach to obtain the bounds for eachedge. Then, we use the SVD-based approach only for edges
that are determined to update weights by the incremental
approach as described in the next section.
Algorithm 1 gives a full description of Castnet, which efficiently constructs the lasso-based L1 -graph. In Algorithm 1,
W is an N ×N matrix whose p-th row corresponds to weightvector wp , and P is the set of nodes selected for weight computation. Castnet initializes weights based on the observation that the weight of an edge from node u to v is similar
to that from node v to u; wv [u] ≈ wu [v] [4].

We compute the upper and lower bounds after the iterations
based on Lemma 7 as follows if an edge is included in edge set
A and added to edge set U in the κ-th addition (1 ≤ κ ≤ ν):
K[u|wpν ] = λ1 zp [u|wpκ ] − ∆pThus, from Lemma 7, we have

4.4 Graph Construction Algorithm(17)

K[u|wpν ] = λ1 zp [u|wpκ ] + ∆p2This is because we have $x$
u $2 = M . Therefore, from Equation (17), |K[u|wpi ] − K[u|wpi−1 ]| ≤ δpi . As a result, for anedge included in the edge set,"
|K[u|wpν ] − K[u|wpκ ]| = | νi=κ+1 (K[u|wpi ] − K[u|wpi−1 ])|""≤ | νi=κ+1 δpi | ≤ | νi=1 δpi | = ∆p235

Algorithm 1 Castnet

have no edge to add. Therefore, it proceeds the step of constructing the graph if we have A = ∅ (lines 34-35). After theiterations, it sets the weights of edges to unselected nodes
based on the observation (lines 37-38).
In Algorithm 1, we implicitly assumed a single threadwhere we pick up nodes one by one from the entire set ofnodes. However, Algorithm 1 is easily parallelized; we divide the set of nodes into several sets and parallelly perform
Algorithm 1 for each obtained set. Since we compute edges
from all the nodes for each obtained set, we can obtain thesame graph as the single-thread approach even if we use theparallelization approach. We use the k-means clustering toeffectively divide the set according to the hidden structureof data points. Although the k-means clustering performed
in a single thread, we can efficiently divide the entire set ofnodes since (1) we use SVD to reduce the number of dimensions and (2) we apply k-means++ to increase the processing
speed of the k-means clustering [1].
For Algorithm 1, we have the following properties:

Input: matrix X, tuning parameter λ, target rank of SVD m
Output: matrix W
1: P = ∅;2: compute rank-m SVD of matrix X;
3: for i = 1 to N do
4:
p = argmax(#wu #1 |u ∈ V\P);5:
U = {ep [u] | wp [u] %= 0};
6:
update weights for U by Equation (12);7:for j = 1 to N do
8:K[uj |wp ] = −∞, K[uj |wp ] = ∞;9:step = 1;
10:while step ≤ 2 do
11:
for each u ∈ V\p do
12:compute the bounds by Definition 5;
13:
if step = 1 then14:M[wp ] = ∅;15:
for each u ∈ V\p do
16:
if K[u|wp ] < −1 or K[u|wp ] > 1 then17:
compute the bounds of node u by Definition 4;18:
if K[u|wp ] < −1 or K[u|wp ] > 1 then19:
add edge ep [u] to M[wp ];20:
A = M[wp ]\U;
21:
else
22:
C[wp ] = ∅;23:
for each u ∈ V\p do
24:
if K[u|wp ] > 1 or K[u|wp ] < −1 then25:
compute the bounds of node u by Definition 4;26:
if K[u|wp ] > 1 or K[u|wp ] < −1 then27:
add edge ep [u] to C[wp ];28:
A = C[wp ]\U;29:
for each u ∈ A do
30:
compute wp [u] for weight vector wp by Equation (10);31:
if wp [u] %= 0 then32:
add edge ep [u] to edge set U;
33:
update weights for U by Equation (12);34:
if A = ∅ then35:
step + +;36:
add node p to P;
37:
for each u ∈ V\P do
38:
wu [p] = wp [u];return weights of the graph;

Theorem 1 (Computation Cost). Let t be the number of update computations in our approach, Castnet takesO(dp t+M (N +log m+t+(dp )2 )) time to compute weights ofa selected node in the lasso-based L1 -graph.
Proof
Castnet first computes the SVD of matrix X
used in computing the bounds at O(M log m) time amortized on the number of nodes [12]. In the SVD-based boundcomputations, it takes O(dp m) time to compute norm $r̃p $2
and O(N m) time to compute the bounds. Similarly, weneed O(N ) time to compute the bounds by the incremental
computation approach. In addition, we can compute edgesets M[wp ], C[wp ], and A in O(N ) time. In performing update computations for edges whose weights are zero, it takes
O(M t) time to compute the residuals and O(M (N − dp ))time to update computations. In updating nonzero weights,
it needs O(M (dp )2 ) time to compute the inner products andO(dp t) time to update weights. Therefore, Castnet needs!
O(dp t+M (N +log m+t+(dp )2 )) time.

Algorithm 1 starts by initializing node set P and computing SVD of data matrix X to allow the upper and lower
bounds to be determined in the iterations (lines 1-2). It
picks up the node that has the maximum L1 norm for theweight vector since we can accurately perform the regression
for the node from the observation (line 4). It then performs
update computations for edges that have nonzero weights
by using Equation (12) and initializes the bounds (lines 58). Our approach has two steps in computing the weights of
selected nodes; the first step is based on edge set M[wp ] andthe second step exploits edge set C[wp ]. In order to prune
unnecessary updates, it computes the bounds determinedby the incremental approach as described in Section 4.3.2(lines 11-12) and checks the conditions of edge set M[wp ]and C[wp ] (line 16 and line 24). If an edge can meet theconditions of the edge sets (Definition 2 and 3), it computes the bounds again by using the SVD-based approachto determine the edge sets (lines 17-19 and lines 25-27). As
described in Section 4.2.3, it adds edges one by one to edge
set U so as to increase the efficiency. Therefore, it computes
edge set A after determining edge sets M[wp ] and C[wp ](line 20 and line 28). If the weight of an edge is zero, we
can efficiently compute its weight by using Equation (10)(line 30); such edges have no impact on the regression result
as shown in Lemma 5. Therefore, it updates weights untilconvergence only if the added edge has nonzero weight (lines
31-33). The weight vector clearly reaches convergence if we

In practice, the number of update computations t increases with the number of nodes N . However, we canreduce t if a node has small degree dp . We show these relationships in Section 5.
Theorem 2 (Construction result). If matrix X has
full row rank, the proposed approach converges on the same
regression result as the original coordinate descent approachin computing edge weights for a selected node.
Proof As shown in Algorithm 1, Castnet first updateweights if K[u|wp ] < −1 or K[u|wp ] > 1 holds. It thenupdate weights if K[u|wp ] > 1 or K[u|wp ] < −1 hold. As
shown in Lemma 3, we have K[u|wp ] > 1 or K[u|wp ] < −1
if K[u|wp ] < −1 or K[u|wp ] > 1 holds. Therefore, thesetwo processes of the proposed approach indicate that we do
not perform update computations for weights of edges suchthat K[u|wp ] ≤ 1 and K[u|wp ] ≥ −1. It is clear that suchedges cannot have nonzero weights from the KKT condition(Definition 1). As a result, our approach cannot prune edges
that have nonzero weights. Since there is a unique regression
result if matrix X has full row rank, the coordinate descentapproach converges to the unique result [6, 24]. Since ourapproach is based on the coordinate descent, it is clear that
Castnet converges the same result as the original coordinate
descent approach if matrix X has full row rank.
!

236

10

6

10

5

CastnetSling
Original

10

104
103
102
101
10

0

USPS Reuters Protein SensIT

MSD

INRIA

6

10

5

10

4

106

CastnetSling
Original

10
Wall clock time [s]7Wall clock time [s]Wall clock time [s]

10

103
102
10

1

10

0

USPS Reuters Protein SensIT(1) λ = 0.1

MSD

INRIA

5

CastnetSling
Original

104
10

3

10

2

10

1

10

0

USPS Reuters Protein SensIT(2) λ = 0.2

MSD

INRIA

(3) λ = 0.3

Figure 1: Graph construction time of each approach.
If matrix X does not have full row rank, which is necessarily the case when N > M , there may not be a uniqueregression result. In this case, our approach yields almost
the same result as the original coordinate descent approach.
In the next section, we detail experiments that show theefficiency and effectiveness of the proposed approach.

Table 2: Degree and updates in INRIA.

DegreeUpdates

5.

EXPERIMENTAL EVALUATION

We performed experiments to demonstrate the effectiveness of our approach. In the experiments, we used the six
datasets taken from various domains: USPS, Reuters, Protein, SensIT, MSD, and INRIA. USPS is a popular datasetof handwritten digits captured from envelopes by the U.S.
Postal Service2 . This dataset contains 2, 007 grayscale images, each 16 × 16 pixels; the number of features is 256.
Reuters is a corpus of newswire stories that contains 8, 293
documents3 . In this dataset, tf-idf is used as the document
feature; it has 18, 933 dimensions. Protein is a set of aminoacid sequences where the number of data and features are17, 766 and 357, respectively4 . SensIT is a dataset obtained
from a distributed wireless sensor network for vehicle typeclassification5 . In this dataset, the number of data and features are 78, 823 and 100, respectively. MSD is a music
dataset that contain 515, 345 songs of 90 features6 . In thisdataset, songs are mostly western and commercial tracks.
INRIA consists of 1, 000, 000 image features extracted fromboth personal and Flickr photos, which of each is represented
by a 128-D SIFT descriptor7 . Since we have N > M , matrix
X is not full row rank except for Reuters while matrix X isfull row rank for the dataset. The numbers of nodes in thedatasets increase in the order of USPS, Reuters, Protein,
SensIT, MSD, and INRIA.
In what follows, “Castnet”, “Sling”, and “Original” represents our approach, the state-of-the-art approach for thelasso [7], and the original coordinate descent approach [5],respectively. We set the target rank of SVD to 10 for the
proposed approach. All the experiments were conducted ona Linux 2.70 GHz Intel Xeon server with 1 TB of main memory. We implemented all the approaches using GCC.

λN = 250, 000

N = 500, 000

N = 1, 000, 000

0.1
0.2
0.3
0.1
0.2
0.3

16.77
9.91
6.76
5.44 × 105
4.82 × 105
4.22 × 105

15.02
9.34
6.82
1.03 × 106
9.27 × 105
7.56 × 105

15.32
9.59
7.17
2.04 × 106
1.82 × 106
1.53 × 106

parameter λ to 0.1, 0.2, and 0.3. In Figure 1, we omit theresult of Sling for MSD and INRIA datasets since it could
not construct the L1 -graph within a month. Similarly, we
omit the results of the original approach for SensIT, MSD,
and INRIA datasets for the same reason. Table 2 shows theaverage degree and update computations of each node; INRIA was used as the dataset and we changed the numberof nodes in the graph. To evaluate the scalability offeredby our parallelization approach, we plot speedup over single
thread vs. the number of threads for each graph in Figure 2. Figure 3 shows the computation overhead of k-meansclustering in the parallelization approach.
Figure 1 shows that our approach is much faster than theprevious approaches. Specifically, Castnet is up to 1, 300
times faster than the original approach. In order to obtainedge weights, the original approach requires a different regression result for each weight as described in Section 3.
Therefore, the original approach must compute the regression result every time it updates each weight. In addition,
the original approach iteratively updates weights until convergence even though almost all edges have zero weights dueto the sparse structure of the L1 -graph, which explains its
high computation cost. To improve the computation speed
for the lasso, Sling updates weights by using the covariancebased approach as described in Section 2. Since the covariancebased approach does not need the different regression results, Sling can more efficiently update weights in each iteration than the original approach. However, the covariancebased approach needs inner product computations for alledges to update the weights before entering the iterations.
On the other hand, our approach can avoid the heavy computations of inner products by using the two types of update
computations as described in Section 4.2.2. In addition, ourapproach selectively updates weights by computing edge setsfrom upper and lower bounds of KKT condition scores as
described in Section 4.2.1. As a result, our approach is up
to 130 times faster than the state-of-the-art approach. As
shown in Figure 1, computation time of our approach increases with dataset size. This is because we need more
update computations as the number of nodes increases as

5.1 Graph Construction TimeWe assessed the graph construction time needed for eachapproach. Figure 1 shows the results where we set tuning2

https://www.otexts.org/1577
http://www.cad.zju.edu.cn/home/dengcai/Data/data.html
4
http://www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets/5
http://mldata.org/repository/data/6
http://labrosa.ee.columbia.edu/millionsong/7
http://corpus-texmex.irisa.fr/
3

237

20
15
10

10

10

1

0

10

4-thread8-thread
16-thread
32-thread

-1

10

-2

10

-3

5
0

10

INRIA

Figure 2: Effect of parallelization
approach.

USPS Reuters Protein SensIT MSD

Figure 3:
lelization.

Overhead with paral-

10

USPS
Reuters
Protein
SensITMSD
INRIA

m = 10m = 2072.362 × 10
1.471 × 108
1.400 × 109
1.373 × 1010
4.258 × 1011
2.038 × 1012

Datasetm = 307

2.296 × 10
1.470 × 108
1.380 × 109
1.341 × 1010
1.975 × 1011
2.029 × 1012

3

102
1

USPS Reuters Protein SensIT

Figure 4:
number.

MSD

INRIA

Effect of target rank

Table 4: Regression results of each approach.

Table 3: Number of update computations.
Dataset

5

100

INRIA

Castnet(10)Castnet(20)Castnet(30)104

10USPS Reuters Protein SensIT MSD

7

106Wall clock time [s]25

10

4-thread8-thread
16-thread
32-threadWall clock time [s]Speedup over single thread30

Approach

Objectivefunction

Squaredloss

L1 -norm
constraint

CastnetSling
Original
CastnetSling
Original
CastnetSling
Original

0.13550.13550.13550.39850.39850.39850.33230.33230.33230.03480.03480.03480.33600.33600.33600.16510.16510.16510.10070.10070.10070.06250.06250.06250.16720.16720.16727

1.833 × 10
1.469 × 108
1.372 × 109
1.092 × 1010
7.196 × 1010
1.992 × 1012

USPS

Reuters

shown in Table 2. On the other hand, Figure 1 indicatesthat we can increase the efficiency by increasing tuning parameter λ. This is because, as shown in Table 2, the numberof update computations decreases if we reduce the degree of
each node by assigning large values to the tuning parameter.
Since our approach can effectively prune unnecessary computations, it can more efficiently construct the lasso-based
L1 -graph than the previous approaches.
As described in Section 4.4, Castnet can be parallelized
by using multiple threads. As shown in Figure 2, our parallelization approach is up to 25 times faster than the singlethread approach with 32-thread execution. Although ourparallelization approach needs the additional process of kmeans clustering in a single thread, its computation overhead has negligible impact on the graph construction time.
This is because we can efficiently perform k-means clustering by using SVD and k-means++ as shown in Figure 3; ittakes at most 5 seconds to perform the clustering for thelargest dataset. Therefore, we can effectively increase theprocessing speed by exploiting the parallelization approach.

Proteinwe can increase the efficiency of computing the bounds if
rank m has a small score. On the other hand, as shown in
Table 3, we can reduce the number of update computations
if rank m has a large score. This is because we can improve
the accuracy of SVD in approximating the data matrix byincreasing the target rank. In conclusion, we can efficiently
compute the bounds by reducing the target rank at the sacrifice of the number of update computations whereas we canefficiently update weights with the large target rank that
incurs high computation cost for the bounds. Due to thetrade-off derived from rank m, our approach is not sensitive
to the target rank of SVD.

5.3 Effectiveness of the Graph Structure
As described in Section 4.4, one major advantage of ourapproach is that it yields the same L1 -graph as the original approach if matrix X has full row rank. However, ourapproach can output different results from the original approach if matrix X does not have full row rank since there
may not be a unique lasso solution in terms of the objectivefunction introduced in Section 3 (Equation (1)). In this section, we evaluate the effectiveness of the regression result inconstructing the L1 -graph of each approach. Table 4 shows
the averages of the objective function in computing edge
weights for each node in the L1 -graph when λ = 0.1. In thetable, “Squared loss” and “L1 -norm constraint” correspondto the first and second terms of the objective function, respectively. Note that the datasets used in this experiment
do not have full row rank due to N > M except for Reuters
dataset. In addition, the same as our approach, Sling can
output different regression results from the original approachif matrix X does not have full row rank [7].
Table 4 indicates that our approach yields the same object function scores as the original approach for the Reuters
dataset, which does have full row rank. This is because ourapproach provably guarantees to output the same solution

5.2 Efficiency vs. Target Rank Number
As described in Section 4.3, our approach computes SVD
of rank m for the data matrix to increase the computation
speed of the upper and lower bounds. Since the bounds
are used in computing the edge sets for updating weights,
the target rank of SVD is expected to impact the graphconstruction time. In this section, we show the graph construction times recorded for different target rank values of
SVD. Figure 4 shows the results where we set λ = 0.1 as thetuning parameter of the lasso. In this figure, the results of
our approach are indicated by “Castnet(m)” where m is thetarget rank of SVD. In addition, Table 3 shows the number of update computations for all the nodes needed by ourapproach to construct the L1 -graph when λ = 0.1.
As shown in Figure 4, the graph construction time varies
only slightly against the target rank of SVD for the data
matrix; our approach is not sensitive to rank m in terms ofefficiency. As described in Section 4.3, it takes O(m) time
to compute the upper and lower bounds by SVD. Therefore,

238

Table 5: Label propagation.
Precision
Graph Construction [s]Number of edges

Table 6: Connected nodes by each approach.

Castnet

L1

k-NN

LN

0.581
19.19
34, 804

0.576
78.17
39, 3290.433
3.35
36, 9260.554
7.19
30, 426

Connected nodes
Approachfor the L1 -optimization problem as described in Theorem 2.
Although our approach can yield different regression results
if matrix X does not have full row rank, the graph structures
obtained by our approach were almost the same as those by
the original approach. This is because we perform update
computations for edges that must/can have nonzero weights
by computing the upper and lower bounds of KKT condition scores as shown in Algorithm 1. Since KKT condition
scores correspond to gradients in coordinate descent [5], ourapproach can effectively identify edges that effectively improve the objective functions by obtaining the edge set thatmust/can have nonzero weights. Therefore, our approachcan effectively compute the solution for the optimization
problem of the lasso by using the upper and lower bounds.
As a result, scores of the objective function of our approachand the original approach are the same in practice for USPS
and Protein datasets as shown in Table 4 even though thesedatasets do not have full row rank. Along with Figure 1,
Table 4 indicates that our approach significantly improves
the efficiency to construct the L1 -graph while the quality ofgraph structure is as high as the previous approaches.

6.

Castnet

L1

k-NN

CASE STUDIESLN

In order to effectively find the hidden structure of data
points, it is important to construct a graph that is robust
against the noise. Our approach allows users to efficiently
construct the lasso-based L1 -graph so as to enhance thequality of applications more than is possible with traditional
k-NN graphs. This section compares our approach to a kNN graph as well as the L1 -graph proposed by Cheng et al.
[2] and the linear neighborhood graph proposed by Wanget al. [25]. The L1 -graph by Cheng et al. uses only theL1 -norm constraint as the objective function while our approach uses the sum of squared loss and L1 -norm constraint
as described in Section 3. The linear neighborhood graphby Wang et al. computes the edge weight of each node as alinear combination of its neighbors.
We used the COIL-100 dataset that contains images of
100 objects8 . The objects were placed on a turntable thatwas rotated through 360 degrees to vary object pose withrespect to a fixed camera. Images of the object were taken
at pose intervals of 5 degrees; 72 poses per object, resulting
in 7, 200 images. We used bag of keypoints of 128 dimensions as image features where we exploited SIFT descriptor
to computer feature vector [3]. Since we have N > M , data
matrix X does not have full row rank for this dataset. In order to evaluate the effectiveness, we used label propagation[9]; it estimates labels of unlabeled data points from labeleddata points. To this end, it propagates labels to the unlabeled nodes in the given graph from a labeled seed node.
We randomly selected a data point from each object as theseed node. Table 5 shows labeling accuracy and graph construction time of each approach; labeling accuracy is ratioof labeled nodes that correspond to the same object as theseed node. Table 5 also shows the number of edges obtained
8

0.610

0.512

0.653

0.161

0.609

0.353

0.394

0.008

0.014

0.002

0.061

0.060

0.422

0.350

0.330

0.273

0.418

0.279

0.313

0.002

0.270

0.058

0.255

0.070

0.295

0.282

0.252

0.222

0.221

0.206

0.254

0.184

0.220

0.219

0.199

0.190

0.350

0.339

0.374

0.277

0.353

0.299

0.305

0.005

0.193

0.194

0.167

by each approach. In addition, Table 6 shows examples ofconnected nodes and their edge weights obtained by eachapproach. In Table 5 and 6, “L1 ” and “LN ” correspondto the results of the L1 -graph by Cheng et al. [2] and thelinear neighborhood graph by Wang et al. [25]. In constructing k-NN graph, we used FLANN that approximately
finds the nearest neighbors; it is a popular approach used in
OpenCV [18]. We set parameters of FLANN so that it yields
accuracy of 90% in computing the nearest neighbors since
this setting is well used in OpenCV. Note that FLANN outperforms LSH-based approaches in terms of efficiency. We
set λ = 0.4 to construct the L1 -graph and computed the
top four nodes to obtain edges in the k-NN graph. As aresult, each approach yielded similar number of edges asshown in Table 5. For the k-NN graph, we computed edge
weights from normalized graph Laplacians [9]. Note that
our approach yields the same labeling results as the original
coordinate descent approach.
Table 5 indicates that our approach can more effectively
identify semantically equivalent nodes than other graphs.
In addition, as shown in Table 6, our approach successfully
connected semantically equivalent nodes to the picked up
nodes. This is because, by effectively using the regression
approach, the loss-based L1 -graph can connect semanticallyequivalent nodes. Therefore, semantically equivalent nodes
comprise a cluster structure in the high dimensional feature
spaces. In terms of constructing the L1 -graph, the approachby Cheng et al. is similar to our approach. However, since
their approach does not use the squared loss as the objective function, it is not robust against noise as described in

http://www1.cs.columbia.edu/CAVE/software/239

the previous study [4]. This indicates that our approach can
improve the robustness against noise by using the squared
loss in the objective function. For example, in the third caseof Table 6 where a white car is a node, both our and theirapproaches connected a green car to the node. However, ourapproach assigns much lower edge weight to the semanticallydifferent node than semantically equivalent nodes. On theother hand, their approach assigns high edge weight to thesemantically different node. In addition, Table 5 indicatesthat our approach is more efficient than their approach inconstructing the L1 -graph. This is because our approachskips unnecessary weight updates. In terms of efficiency,
k-NN graph can be constructed more efficiently than thelasso-based L1 -graph as shown in Table 5. However, k-NN
graph as well as linear neighborhood graph offers lower labeling accuracy than our approach. This is because, as shownin Table 6, they are not so effective as our approach in connecting semantically equivalent nodes. For example, in thefirst case where a white cat is a node, k-NN graph connects
red doll to the node. In addition, they give high weights tosemantically different nodes as shown in Table 6.
Although the traditional k-NN graph is the most popular
approach, it is sensitive to data noise as shown in Table 5.
Effectiveness is quite important in many applications suchas lesion detection as described in Section 1.1; we can reduce
deaths caused by prostate cancer by effectively performing
lesion detection. The lasso base L1 -graph was proposed to
improve the sensitive to data noise. Although it overcomes
the drawback, the original coordinate descent approach does
not scale well to handle large data sets as shown in Figure 1.
The contribution of our approach is to significantly increase
the efficiency of constructing the L1 -graphs without sacrificing the usefulness of the graph structure. As a result, theproposed approach is an attractive option for the researchcommunity in constructing the lasso-based L1 -graphs.

7.

[4] E. Elhamifar and R. Vidal. Sparse Subspace Clustering. In
CVPR, pages 2790–2797, 2009.
[5] J. Friedman, T. Hastie, H. Höfling, and R. Tibshirani.
Pathwise Coordinate Optimization. Annals of Applied
Statistics, 1(2):302–332, 2007.
[6] J. H. Friedman, T. Hastie, and R. Tibshirani.
Regularization Paths for Generalized Linear Models via
Coordinate Descent. Journal of Statistical Software,
33(1):1–22, 2 2010.
[7] Y. Fujiwara, Y. Ida, H. Shiokawa, and S. Iwamura. Fast
Lasso Algorithm via Selective Coordinate Descent. In
AAAI, pages 1561–1567, 2016.
[8] Y. Fujiwara and G. Irie. Efficient Label Propagation. In
ICML, pages 784–792, 2014.
[9] Y. Fujiwara, G. Irie, S. Kuroyama, and M. Onizuka.
Scaling Manifold Ranking Based Image Retrieval. PVLDB,
8(4):341–352, 2014.
[10] Y. Fujiwara, M. Nakatsuji, H. Shiokawa, Y. Ida, andM. Toyoda. Adaptive Message Update for Fast Affinity
Propagation. In SIGKDD, pages 309–318, 2015.
[11] Y. Fujiwara, M. Nakatsuji, H. Shiokawa, T. Mishima, and
M. Onizuka. Fast and Exact Top-k Algorithm for
PageRank. In AAAI, 2013.
[12] N. Halko, P. Martinsson, and J. A. Tropp. Finding
Structure with Randomness: Probabilistic Algorithms forConstructing Approximate Matrix Decompositions. SIAM
Review, 53(2):217–288, 2011.
[13] T. Hastie, R. Tibshirani, and J. Friedman. The Elements ofStatistical Learning. Springer, 2011.
[14] S. Ji. Computational Network Analysis of the Anatomicaland Genetic Organizations in the Mouse Brain.
Bioinformatics, 27(23):3293–3299, 2011.
[15] S. Liao, Y. Gao, and D. Shen. Sparse Patch Based Prostate
Segmentation in CT Images. In MICCAI, pages 385–392,
2012.
[16] N. Meinshausen and P. Bühlmann. High-Dimensional
Graphs and Variable Selection with the Lasso. The Annalsof Statistics, 34(3):1436–1462, June 2006.
[17] T. Mishima and Y. Fujiwara. Madeus: Database Live
Migration Middleware under Heavy Workloads for Cloud
Environment. In SIGMOD, pages 315–329, 2015.
[18] M. Muja and D. G. Lowe. Scalable Nearest Neighbor
Algorithms for High Dimensional Data. IEEE Trans.
Pattern Anal. Mach. Intell., 36(11):2227–2240, 2014.
[19] M. Nakatsuji, Y. Fujiwara, A. Tanaka, T. Uchiyama, and
T. Ishida. Recommendations over Domain Specific User
Graphs. In ECAI, pages 607–612, 2010.
[20] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P.
Flannery. Numerical Recipes 3rd Edition. Cambridge
University Press, 2007.
[21] H. Shiokawa, Y. Fujiwara, and M. Onizuka. SCAN++:
Efficient Algorithm for Finding Clusters, Hubs and Outlierson Large-scale Graphs. PVLDB, 8(11):1178–1189, 2015.
[22] J. M. Steele. The Cauchy-Schwarz Master Class: An
Introduction to the Art of Mathematical Inequalities.
Cambridge University Press, 2004.
[23] R. Tibshirani. Regression Shrinkage and Selection via the
Lasso. Journal of the Royal Statistical Society, Series B,
58:267–288, 1996.
[24] R. Tibshirani, J. Bien, J. Friedman, T. Hastie, N. Simon,
J. Taylor, and R. J. Tibshirani. Strong Rules for Discarding
Predictors in Lasso-type Problems. Journal of the Royal
Statistical Society: Series B (Statistical Methodology),74(2):245–266, 2012.
[25] F. Wang and C. Zhang. Label Propagation through Linear
Neighborhoods. IEEE Trans. Knowl. Data Eng.,
20(1):55–67, 2008.

CONCLUSIONS

We addressed the problem of efficiently constructing thelasso-based L1 -graph. Our approach, Castnet, limits update computations to the edges that must/can have nonzero
weights before entering the iterations by computing the upper and lower bounds of KKT condition scores. In addition,
our approach efficiently updates nonzero weights in each iteration by pruning edges whose weights are zero. Experiments show that our approach can achieve higher efficiency
than existing approaches. Since the L1 -graph can effectively
capture the cluster structures that share useful semantic information, it plays a fundamental role in many applications
such as brain analysis, motion segmentation, and lesion detection. Our approach allows many applications to be implemented more efficiently, and will help to improve the effectiveness of future data mining applications.

8.

REFERENCES[1] D. Arthur and S. Vassilvitskii. k-means++: The
Advantages of Careful Seeding. In SODA, pages 1027–1035,
2007.
[2] B. Cheng, J. Yang, S. Yan, Y. Fu, and T. S. Huang.
Learning with l1 -graph for Image Analysis. IEEE Trans.
Image Processing, 19(4):858–866, 2010.
[3] C. Dance, J. Willamowski, L. Fan, C. Bray, and G. Csurka.
Visual Categorization with Bags of Keypoints. In ECCVInternational Workshop on Statistical Learning in
Computer Vision, 2004.

240