Machine Learning, 105(3):335-366, 2016. This is a preprint version. The final publication isavailable at http://link.springer.com/article/10.1007/s10994-016-5568-6.

Sparse and Low-Rank Matrix Regularization for Learning
Time-Varying Markov Networks
Jun-ichiro Hirayama1,* , Aapo Hyvärinen2,1 , and Shin Ishii3,1
1

Brain Information Communication Research Laboratory Group, Advanced
Telecommunications Research Institute International (ATR), Kyoto, Japan2
Department of Computer Science and HIIT, University of Helsinki, Helsinki, Finland
3
Graduate School of Informatics, Kyoto University, Kyoto, Japan*hirayama@atr.jp

Abstract
Statistical dependencies observed in real-world phenomena often change drastically withtime. Graphical dependency models, such as the Markov networks (MNs), must deal with this
temporal heterogeneity in order to draw meaningful conclusions about the transient natureof the target phenomena. However, in practice, the estimation of time-varying dependency
graphs can be inefficient due to the potentially large number of parameters of interest. To
overcome this problem, we propose such a novel approach to learning time-varying MNs
that effectively reduces the number of parameters by constraining the rank of the parameter
matrix. The underlying idea is that the effective dimensionality of the parameter spaceis relatively low in many realistic situations. Temporal smoothness and sparsity of the
network are also incorporated as in previous studies. The proposed method is formulated
as a convex minimization of a smoothed empirical loss with both the `1 - and the nuclearnorm regularization terms. This non-smooth optimization problem is numerically solvedby the alternating direction method of multipliers (ADMM). We take the Ising model as afundamental example of an MN, and we show in several simulation studies that the rankreducing effect from the nuclear norm can improve the estimation performance of timevarying dependency graphs. We also demonstrate the utility of the method for analyzing realworld datasets for improving the interpretability and predictability of the obtained networks.

1

Introduction

The sparse Markov Network (MN) [37, 52, 32, 18, 3, 23] is a powerful model for analyzing
statistical dependency in multivariate data. Like the ordinary MN, it describes conditional(in)dependences between random variables of interest in terms of an undirected graph, quantified by model parameters associated with the graph. A sparse MN, in particular, learns thegraph structure and parameters using sparse estimation techniques, such as the `1 -norm regularization [48, 11], formulated as a convex optimization that can be solved efficiently. This
computational efficiency of the sparse MN is a great practical advantage, when compared toclassical methods based on combinatorial optimization of the graph structure.
A challenging issue with the sparse MNs arises when attempting to deal with the dependency
graph that is potentially heterogeneous in a given dataset. In particular, statistical dependency
observed in dynamical phenomena often changes drastically with time. The dependency graphs
varying over time then must be considered in order to analyze any transient nature of the target1phenomena. The conventional “static” case of the sparse MNs is not useful for this purpose, as
it assumes that the dependency graph is consistent for all time points.
However, estimating time-varying dependency graphs is rather difficult due to the potentiallylarge number of parameters of interest: the number of edges or their weights grows quadraticallywith the number of variables and linearly with the sample size (i.e., the number of time points).
This high degree of freedom of the model typically degrades the parameter estimation due torelative insufficiency of data. Previous studies on the time-varying MNs [46, 28, 58, 29] haveaddressed this issue, by assuming that no substantial changes occur in the network over a
relatively short time period, which is likely to be reasonable in many cases.
In the present study 1 , we propose a new framework for learning time-varying MNs, focusing
on particular situations where the parameter matrix of the time-varying MN can be reasonably
assumed to be of relatively low rank according to some prior knowledge. As illustrated in Fig. 1,
every instance of the dependency graph of a time-varying MN is described by a single column ofthe parameter matrix Θ. The low-rank assumption of Θ equivalently means that the columns
and the rows of Θ are constrained to low-dimensional subspaces. Intuitively, the basis vectors
in the column and row subspaces correspond to the patterns in the dependency graph and to
the time courses of these patterns’ coefficients, respectively. The low rank then implies that the
number of these patterns is relatively small. Although, in reality, the true Θ may also contain
unstructured noise and thus may not be exactly rank-deficient, we can still expect that the true
Θ is well approximated by a low rank matrix in many situations.
We should emphasize here that we do not seek a method that works in every situation orfor every dataset. Rather, we focus on improving the estimation in specific situations where
the low-rank assumption is reasonable. However, we still expect the method to have a wide
range of applicability covering many practical situations. An example of this type of situation,suggested directly by the visualization in Fig. 1, is a network with a community (or cluster)
structure [38]. The low-rank property will emerge if the overall changes in the network are dominated by dynamic activities of major communities that produce specific dependency patterns
in the observed data. Another important situation is regime switching, a well-known concept
in time series analysis in which the network visits a limited number of regimes, possibly manytimes. Actually, the idea of the low rank parameter matrix is a generalization of regime switching because now each regime is not only selected from a finite set of dependency patterns but
each may also be any linear combination of them.
Here, we focus on the time-varying Ising model [46, 28], a fundamental instance of timevarying MNs, and extend its basic framework of sparse estimation [46, 28], formulated as suchto employ `1 -norm regularization and kernel smoothing. In particular, we introduce the nuclear(trace) norm regularization [15, 47] into this basic framework, which has been widely used
for learning low-rank matrices in recent years. The proposed sparse and low-rank estimation
framework of time-varying MNs can still be formulated as a convex optimization and thus solved
efficiently. We also introduce an optimization algorithm based on the alternating direction
method of multipliers (ADMM) [4, 14, 7], which solves the problem by effectively splitting thetwo non-smooth regularization terms.
The remainder of this paper is organized as follows. We first briefly introduce the Isingmodel and the sparse estimation method for this model, using the pseudolikelihood [5, 6] (Section 2). We then propose our new method for learning time-varying MNs based on the low-rank
assumption (Section 3), and we briefly summarize related studies (Section 4). We also present
simulation results using several artificially-generated datasets (Section 5), and demonstrate theapplicability of the proposed method to a real-world problem (Section 6). Finally, we discuss
1

A preliminary version of this paper appeared in a conference proceeding [22], which introduced the basic
algorithm and provided a simple simulation result (included as Section 5.1 in this extended version).

2

Figure 1: Illustration of the basic idea. In a time-varying Markov Network, the n-th instanceof (undirected) dependency graph (n = 1, 2, . . . , N ) is represented by the n-th column of theparameter matrix Θ, depicted as a block on the left-hand side of the equality. Each elementof these column vectors corresponds, for example, to an edge weight between two particular
nodes (out of D nodes), which describes a pairwise interaction between two random variables.
Here, in particular, we assume that the parameter matrix has a relatively low rank, implicitly
decomposed into the two matrices; i.e., the dependency patterns and the coefficient time-series,
depicted as the first and second blocks on the right-hand side, respectively.
the results and open issues (Section 7).

2

Preliminaries

In this section, we introduce the Ising model and its sparse estimation using pseudolikelihood.

2.1

Ising model

An MN represents conditional independences between random variables by means of an undirected graph, where the parameters associated with cliques (i.e., complete subgraphs) quantify
the full probability distribution of the random variables [40, 31]. The Ising model is one of thesimplest examples of an MN, having only parameters associated with nodes and edges so thatit represents, at most, pairwise interactions between binary variables.
Let y = (y1 , y2 , . . . , yD )> be a D-dimensional binary random vector that takes a value in{−1, 1}D . The probability distribution of y under the Ising model, parameterized by a realvalued vector θ, is then defined by
XX
1p(y; θ) =exp
θij yi yj +θi yi ,
(1)
Z(θ)
i<j

i

where the first summation in the exponent
is taken over allpairs (i, j) that satisfy i < j and the
PPP
partition function, Z(θ) = y exp
i<j θij yi yj +i θi yi , simply makes the total probability
equal to unity. The parameter vector θ contains the C (= D(D + 1)/2) network parameters,
θi ’s and θij ’s, in a certain fixed order. Here, we refer to θi ’s and θij ’s as node-wise and edgewise network parameters (or “edge weights”), respectively. For convenience, we also introduce
3auxiliary variables θji for any i < j, and use θij and θji interchangeably by implicitly assumingθji = θij .
The undirected graph associated with this model has D nodes, each of which corresponds toa single variable yi (i = 1, 2, . . . , D). The graph has an undirected edge (i, j) between the i-thand j-th nodes if θij is non-zero, and otherwise has no edge. If the two nodes are connected(disconnected), then the variables yi and yj are conditionally dependent (independent), given
their Markov blanket (i.e., all other variables that have direct connections to either or both ofthese variables). The dependence structure of the constituent binary variables is then specified
by the graph structure (i.e., which edge-wise parameters are zero and which are non-zero).

2.2

Estimation by maximum pseudolikelihood

In the Ising model, as in many other MNs, computing the likelihood or its gradient is often
difficult due to the intractability of the partition function Z(θ), unless D is sufficiently small.
One approach for overcoming the intractability of likelihood is to use the pseudolikelihood [5,6, 25, 30]; many variants are in popular use in the field of machine learning [52, 42, 23, 20].
The pseudolikelihood for a parametric model p(y;Q θ), given a single datum y, is generally
defined as a product of conditional likelihoods; i.e., D
i=1 p(yi | y \i ; θ). Here, y \i denotes all thevariables of interest other than yi . The corresponding loss function l(y, θ) is introduced as thenegative logarithm of the pseudolikelihood, which is given specifically for the Ising model byl(y, θ) = −DXlog p(yi | y \i ; θ) =

i=1

D
X

{log cosh(ξi ) − yi ξi + log 2} ,

(2)

i=1

P
x
−x
where ξi :=j:j6=i θij yj + θi and cosh(x) = (e + e )/2. Its partial derivatives needed fornumerical optimization are given by∂l(y, θ) = tanh(ξi )yj + yi tanh(ξj ) − 2yi yj ,∂θij∂l(y, θ) = tanh(ξi ) − yi , for any i,∂θifor any i < j,

(3a)
(3b)

P
n
where tanh(x) = (ex − e−x )/(ex + e−x ). The minimization of the empirical loss, N
n=1 l(y , θ),then yields the maximum pseudolikelihood estimator of θ, which has been shown to be consistent [25]. Note that taking the second derivatives readily reveals that this loss function l isconvex with respect to θ.

2.3

Sparse estimation using `1 -norm regularization

Determination of the graph structure can be facilitated by accommodating sparse estimation
based on the `1 -norm regularization [37, 52, 32, 18, 3, 23] with the maximum pseudolikelihood.
In a continuous optimization of the parameter vector θ, the additional regularization term
encourages the optimal θ to have many entries that are exactly zero (i.e., effectively determining
the graph structure).
The loss function introduced above now formulates the sparse estimation of the Ising model
as a convex minimization problem, such thatN

minimizeθ

1 Xl(y n , θ) + kλ
N
n=1

4

θk1 ,

(4)

where k · k1 denotes the `1 -norm (i.e., the sum of all absolute values of the elements), λ is avector of regularization coefficients, anddenotes element-wise multiplication. The objectivefunction is convex, since both the loss function l and the `1 -norm are convex with respect to θ
and any sum of convex functions is again convex [8].
The second term of Eq. (4) can also be written asXXkλ θk1 =λij |θij | +λi |θi |.
(5)
i<j

i

We typically use λi = 0 for any i, because only the edge-wise parameters are relevant for making
the graph structure sparse. In addition, for simplicity, we often use a common λ for any i < j.

3

Proposed method

The time-varying MN [46, 1, 28, 58] extends the sparse MN to incorporate potential heterogeneityof the dependency graph over time. It assumes that N successive observations or measurements
y 1 , y 2 , . . . , y N are generated independently from the MNs p(y n ; θ n ) with different parameter
vectors θ n at each time step. The corresponding time-varying dependency graphs are thenassociated with the parameter matrix Θ = (θ 1 , θ 2 , . . . , θ N ) ∈ RC×N , and each column specifies
an instance of the MN. The above sparse estimation framework now requires the incorporationof additional assumptions about Θ other than sparsity; otherwise, each instance of the MNp(y n ; θ n ) is estimated only with a single datum y n .
In this section, we extend a previous approach [46, 58, 28] by incorporating an additional assumption whereby the parameter matrix Θ has a relatively low rank. Based on this assumption,
we formulate the sparse and low-rank estimation of the time-varying MN as a convex optimization problem; in particular, we use joint regularization by both the `1 -norm for sparsity andthe nuclear norm for rank reduction. We then introduce an algorithm for solving this problem
based on the ADMM.

3.1

Previous kernel smoothing method

Previous studies assumed that temporally adjacent θ n ’s are similar to each other. Here, we
extend the existing approach based on kernel smoothing [46, 58, 28]. According to this approach,parameters are estimated using locally weighted averages of the loss function based on kernelsmoothing, which makes the estimates temporally smooth. In the context of the present study,
the smoothed version of empirical loss, averaged over all instances, is given by 2
N

N

1 XX
fw (Θ) :=ϕw (m − n)l(y m , θ n ),N(6)

n=1 m=1

for any Θ. In Eq. (6), l is the loss function of Eq. (2), and ϕw is the smoothing kernel (i.e.,a nonnegative symmetric function centered at zero [34]), depending on the kernel width w.
Examples of the smoothing kernels are the box kernel,(
1/(2w + 1) |z| ≤ w
ϕw (z) =,(7)0otherwise2Here, we assume that y n are sampled regularly in time. For non-regular sampling intervals, the smoothingkernel can be defined from the difference between the time stamps of two data points y m and y n , as in [28].

5

0.08

Box
Epa.

0.06

w=10

0.04
w=20

0.02

w=50

0
−50

0

50

Figure 2: Smoothing kernels: Box (Eq. (7)) and Epanechnikov (Eq. (8)) with the kernel widthsof w = 10, 20, and 50.
and the Epanechnikov kernel,(3(1 − (z/w)2 )/(4w) |z| < w
ϕw (z) =.
0
otherwise(8)

Figure 2 illustrates these smoothing kernels. The kernel width determines the smoothness andshould be carefully selected, whereas the precise form of the kernel function is usually not very
important [34].

3.2

Proposed sparse and low-rank estimation of time-varying MN

We now focus on situations where the parameter matrix Θ can be reasonably restricted to haverelatively low rank according to any prior knowledge. This further reduces the effective degreeof freedom associated with the large number of parameters in the time-varying MN, even when
the temporal smoothness may not be very effective by its own.
Every instance θ n of the parameter vectors is now implicitly represented in terms of K basis
vectors (dependency patterns) ak and their coefficients snk :
θn =

K
X

snk ak .

(9)

k=1This also implies that the matrix Θ can be decomposed as Θ = AS with A = (a1 , a2 , . . . , aK )and S = (s1 , s2 , . . . , sK )> , as depicted in Fig. 1. We assume K < min{C, N } so that thematrix Θ is rank-deficient. Note that, in reality, true Θ might not be strictly rank-deficient
but is instead approximately of low rank with the singular values approaching to zero relatively
quickly; Eq. (9) then can be seen as an approximation of the true system.
The underlying idea is that networks in the real world often exhibit a limited number ofcommunities or regimes, as already stated in Section 1. Intuitively, each ak in Eq. (9) may
represent a specific dependency graph that corresponds to a community or a regime in the
network. Equation (9) actually represents a generalized form of the regime switching because
each network θ n now can be any linear combination of the dependency patterns.
Here, we propose to extend the previous kernel smoothing method into that with the lowrank assumption. This could be achieved by minimizing the kernel-weighted empirical loss of6

Eq. (6) with respect to A and S under the low-rank decomposition Θ = AS described above.
However, joint minimization of the two decomposition terms is not a convex problem and thus
cannot avoid the issue of local minima, which may limit its applicability in practice. We thus
adopt an alternative approach which formulates the problem as a convex optimization directlyon Θ without explicitly decomposing it.
The proposed sparse and low-rank estimation of Θ is now formulated asminimizeΘ

fw (Θ) + kΛ

Θk1 + ηkΘk∗ ,(10)where the two regularization terms are introduced in the objective function, as well as thekernel-weighted empirical loss. Here, the `1 -norm k · k1 is of a long vector that concatenates thecolumns. The nuclear (or trace) norm [15, 47] k · k∗ is defined asXkΘk∗ =σk (Θ),(11)k

where σk (Θ) ≥ 0 is the k-th singular value of Θ. The nuclear norm therefore means the `1 norm of the vector of singular values, and its minimization introduces sparsity on the singularvalues, encouraging the matrix to have a low rank [15, 47]. The two regularization coefficients(i.e., the non-negative matrix Λ and the positive scalar η) control the strengths of the `1 - andthe nuclear-norm regularizations, respectively. For simplicity, we set Λ = (λ, λ, . . . , λ) with acommon λ in every column. Note that the problem (10) is convex because all three terms in
the objective function are convex [8].
In this formulation, the rank reduction by the nuclear norm introduces similarity between any
columns θ n or between any rows, as it is invariant under permutation of columns or rows. The
kernel smoothing also introduces a similarity, but only between temporally-adjacent columns.
The nuclear norm is thus expected to be particularly useful when some networks located at
distant time points (relative to the specified kernel width) are potentially similar, or when theweights of some edges potentially have similar time courses. Note that the nuclear norm by itself
never encourages any two networks θ n1 and θ n2 to differ. If the loss terms l(·, θ n1 ) and l(·, θ n2 )actually promote differences in the networks, the nuclear norm oppositely encourages them tobe similar, cooperatively with kernel smoothing.

3.3

Estimation algorithm by ADMM

The minimization problem (10) is convex and unconstrained, but the objective function is nonsmooth due to the regularization terms. Hence, we cannot directly apply standard unconstrained
optimization techniques. In this subsection, we derive a simple first-order algorithm to solve
this problem (10) based on the alternating direction method of multipliers (ADMM) [4, 14, 7],
which is a variant of the augmented Lagrangian (AL) method. The standard AL method has
been effectively used to solve problems with either the `1 -norm regularization (e.g., in [50]) or
the nuclear-norm regularization (e.g., in [51] 3 ), but the joint use of the `1 and nuclear norms
prevents an effective application of the standard AL. This motivated us to use the ADMM.
In the following sections, we first introduce the ADMM in a general form and then apply
the ADMM to an equality-constrained problem that is equivalent to the problem consideredherein (10).
3Note that these studies applied the standard AL to the Fenchel dual of the original `1 - or nuclear-regularized
problems, rather than directly to the original problems, which is also possible.

73.3.1

General framework

Consider the following optimization problem that includes two convex functions φ and γ:
minimize

φ(x) + γ(z)subject to Jx = z

,

(12)

with respect to real-valued vectors x and z, where J is a matrix of appropriate size that represents
linear constraints. To solve this problem, the ADMM iterates the following three steps from any
initial conditions x(0) , z (0) , r (0) until a given convergence criterion is satisfied:
o
n
α
2
(13a)
Jx − z (t−1) 2 ,x(t) = argmin φ(x) + hr (t−1) , Jxi +
2
x
n
o
α
2
z (t) = argmin γ(z) − hr (t−1) , zi +Jx(t) − z 2 ,
(13b)
2
z

r (t) = r (t−1) + α Jx(t) − z (t) ,(13c)where h·, ·i and k · k2 are the standard inner product and the `2 -norm, and α > 0 is any
positive constant. The superscript (t) denotes the iteration number. Each single iteration stepof the ADMM can be seen as a single cycle in an alternating minimization of the augmented
Lagrangian [4, 39]αkJx − zk22 ,2

Lα (x, z, r) = φ(x) + γ(z) + hr, Jx − zi +(14)

with respect to primal vectors x and z, followed by the update of the dual vector r. Note thatin the standard AL method, the primal vectors x and z are simultaneously updated to jointly
minimize Lα with the dual vector r fixed, instead of updating them by Eqs. (13a) and (13b)
only once.
In the ADMM, the constant α can be chosen rather freely, as in the standard AL method.
This is in contrast to the penalty method [39], in which the strength of the penalty should be
sufficiently large, which can cause numerical instability.
3.3.2

Derivation for the problem of the present study

We apply the ADMM by reformulating the problem (10) in a similar manner as in [4, 16] asfollows. This reformulation allows us to deal separately with the loss term and the two different
regularization terms, leading to a simple iterative algorithm summarized in Algorithm 1.
We first introduce auxiliary variables Z1 , Z2 , and Z3 in RC×N and defineγ(Z) = fw (Z1 ) + kΛZ2 k1 + ηkZ3 k∗ ,

(15)

>
> >
with Z = (Z>
1 , Z2 , Z3 ) . We can then rewrite the original problem (10) as

minimize
X,Z

φ(X) + γ(Z)subject to JX = Z

,

(16)

where φ(X) ≡ 0, J = (IC , IC , IC )> and IC is a C × C unit matrix. Note that the equality
constraint implies that X = Zq for any q. The two problems (10) and (16) are equivalent in thesense that, given the optimal solutions Θ? for (10) and (X? , Z? ) for (16), Θ? = X? = Z?q holds
for any q.

8

>
> >
Therefore, with the dual variables M = (M>
1 , M2 , M3 ) , the ADMM procedure can beapplied in a straightforward manner to this problem. The first step, Eq. (13a), is now a quadratic
minimization and thus has a closed-form solution:
1  (t−1)
1  (t−1)
(t−1)
(t−1)
(t−1)
(t−1)
X(t) =
Z1
+ Z2
+ Z3
M1
+ M2
+ M3
−
.
(17)
3
3α

The second step, Eq. (13b), can be separately written as
o
n
α
(t)
(t)
Z1 = argmin fw (Z1 ) + kZ1 − G1 k2F ,
2
Z1
o
n
α
(t)
(t)
Z2 = argmin kΛ Z2 k1 + kZ2 − G2 k2F ,
2
Z2
n
o
α
(t)
(t)
Z3 = argmin ηkZ3 k∗ + kZ3 − G3 k2F ,
2
Z3
(t)

(18a)
(18b)
(18c)

(t−1)

where Gq = X(t) + α−1 Mq(q = 1, 2, 3) and k · kF is the Frobenius norm. We can then(t)(t)
obtain Z1 numerically using any standard unconstrained optimization technique, and Z2 and(t)
Z3 can be obtained as follows:


(t)
(t−1)
Z2 = Soft X(t) + α−1 M2, Λ/α ,
(19a)


(t)
(t−1)
Z3 = Svt X(t) + α−1 M3, η/α .
(19b)
Here, Soft(·, ·) is an element-wise application of the soft-thresholding [19, 26, 12] operator,soft(a, b) = sign(a) max(|a| − b, 0), and Svt(·, ·) denotes the singular value thresholding [9, 49]
operator defined by
Svt(A, b) = Udiag(Soft(σ, b))V> ,(20)

where A = Udiag(σ)V> is the singular value decomposition of A, and σ is the vector of
singular values. Finally, the dual update step, Eq. (13c), is given separately for Mq (q = 1, 2, 3)
(t)
(t−1)
(t)
by Mq = Mq
+ α(X(t) − Zq ).
3.3.3

Stopping criterion and final estimates of Θ

The procedure described in Section 3.3.2 converges to the optimal solution under certain conditions on the accuracies of the two inner minimization problems denoted by Eqs. (13a) and(13b) [14]. Recent theoretical analysis on convergence behavior of ADMM can also be found,for example, in [21]. Although a detailed theoretical analysis of the proposed method is beyond
the scope of our current study, the basic characteristics can be understood from these general
results. In the simulation experiments in Section 5, there was no run that did not converge,
which has empirically proved that the proposed method converges in ordinary situations.
A criterion for stopping an ADMM algorithm was presented in [7]: The algorithm terminates
after the t-th iteration if the following two conditions are satisfied.
δp(t) := JX(t) − Z(t) F ≤ (t)
p ,
(t)

(t)

(t)

δd := α J> (Z(t) − Z(t−1) ) F ≤ d .

(21)

(t)

Here, the tolerances p and d are practically chosen by√(t)3CN abs + rel max{kJX(t) kF , kZ(t) kF },
p =
√
(t)
d = CN abs + rel kJ> M(t) kF ,
9

(22a)
(22b)

Algorithm 1: The ADMM for the problem (16) (equivalent to (10))
Given: w (kernel width), λ (`1 -norm coefficient) and η (nuclear-norm coefficient)X ← 0; Z1 ← 0; Z2 ← 0; Z3 ← 0; M1 ← 0; M2 ← 0; M3 ← 0;
while a specific convergence criterion is not satisfied, do
1
1(Z1 + Z2 + Z3 ) −(M1 + M2 + M3 )
3
3α
o
n
α
Z1 ← argmin fw (Z1 ) + kZ1 − (X + α−1 M1 )k2F
2
Z1

−1
Z2 ← Soft X + α M2 , Λ/αZ3 ← Svt X + α−1 M3 , η/αX←Mq ← Mq + α (X − Zq )for q = 1, 2, 3

end
where abs and rel are the absolute and relative tolerances, respectively, set at small positive
numbers, e.g., (abs , rel ) = (10−5 , 10−4 ) [7]. In addition, slow convergence is avoided by adjusting
(t)
(t)
(t)
(t)
α in each iteration, such that α(t+1) = 2α(t) , 0.5α(t) and α(t) , if δp > 10δd , δd > 10δp andotherwise, respectively.
b Z
b 1, Zb 2 , andAfter the ADMM algorithm has converged, we obtain the primal solutions, X,
b
Z3 . Although in theory, these solutions should be equal to each other and also to Θ? , inpractice, these solutions are not exactly the same because we stop the algorithm when theequality constraints hold only approximately. In the experiments in Section 5, we define theb by X,b whereas each element of Θb is explicitly set at zero if the corresponding
final estimate Θbelement of Z2 is zero. The estimated graph structure can finally be given by whether eachb is zero or non-zero.
element of Θ

4

Related studies

The proposed method particularly assumes η > 0 in order to explicitly include the nuclear normin the objective function. If we set η = 0, the optimization problem of Eq. (10) almost reduces tothe one solved by KELLER [46, 28] (kernel-reweighted logistic regression); KELLER minimizes
the same objective function but with a slightly relaxed constraint. To see this, suppose that
λi = 0 and λij = λ for any i < j, as well as η = 0. It then follows from Eq. (2) and Eq. (6) thatEq. (10) has an equivalent optimization problem, given byo
PN PD n in nPn| ,minimize
f
(θ
)
+
(λ/2)
|θ
w
i
n=1
i=1
j:j6=i ij(23)nn
subject to θij = θji for any i < j and any n.
Here, the optimization variables include both θij and θji for any i < j, with the equality
constraints explicitly introduced. The smoothed empirical loss fw in Eq. (6) is decomposed into
the N D terms, given byNfwin (θ ni ) :=1 X
nϕw (m − n){− log p(yim |y m\i ; θ i )},
N

(24)

m=1

n : j 6= i} is the parameter subset on which the conditional probability
where θ ni := {θin } ∪ {θijactually depends. The problem (23) now can be seen as jointly solving the N D different `1 regularized logistic regression problems based on the weighted loss, while every set of D problems

10

cannot be further decomposed due to the equality constraints. On the other hand, KELLER
relaxes the equality constraints, so that the N D problems can be solved separately. However,it requires a post-processing step to resolve the inconsistency between θbij and θbji .
A large body of literature describes regime switching or statistical change-point detection,with a motivation closely related to that of the present study. For example, Carvalho andWest [10] studied Gaussian covariances in their dynamic matrix-variate graphical modelingframework. Yoshida et al. [53] estimated time-dependent gene networks from microarray databy dynamic linear models with Markov switching. The low-rank assumption on Θ actually
generalizes the concept of regime switching, and this assumption is likely valid in change-point
detection when the number of change-points are relatively small. Note that the method proposed
here has an advantage in its computational efficiency due to the convexity of the optimization
problem.
We introduced a new idea of using low-rank regularization in the specific context of the timevarying MNs. A closely-related idea has recently been presented in [2], which used a low-rank
decomposition to improve the estimation of the time-varying Gaussian Graphical Model (GGM),
the fundamental MN for real-valued data. This explicitly decomposed the GGM parameters into
the fixed patterns and their time-varying coefficients, as in Eq. (9), and directly estimated bothof them. It solved a non-convex optimization problem which would suffer from the issue of
local optima; the convex nature of our method can be beneficial in practice, while a detailed
comparison between these two approaches is left for a future study.
Joint `1 - and nuclear-norm regularization has been studied recently in e.g., [13, 41, 35,57], independently from our preliminary work [22]. However, their target applications werequite different from ours: finding rank-one submatrices [13], multi-task learning [35], matrix
completion [41], and estimating block or community structures in a single, static network [57].
Algorithms similar to our ADMM algorithm were developed in [41, 57]. In particular, Zhou etal. [57] applied the joint `1 - and nuclear-norm regularization, using a slightly different form of
ADMM, to enforce a point-process network model (multi-dimensional Hawkes model) to havea static community structure. Although technically very similar, the underlying idea is quite
different from ours. That is, they considered learning only a single network at a time, assumingthat its (weighted) adjacency matrix is of low rank. In contrast, we here consider learning many
instances of time-varying networks simultaneously, assuming that the entire parameter matrix
Θ is of low rank but the adjacency matrix of each instance (obtained from the corresponding
column θ n appropriately) is not necessarily of low rank. The ideas are thus mutually orthogonal
and they could even be combined in principle.

5

Simulation experiments

We conducted several simulation studies in order to validate the effect of the rank reduction
by nuclear-norm regularization in estimating time-varying dependency graphs. In all theseexperiments, the regularization coefficient for the `1 -norm was set such that λi = 0 and λij = λfor all i < j (see Section 2.3). The minimization in the first step of the ADMM, Eq. (18a), wasnumerically solved by a quasi-Newton method 4 .

5.1

Simulation study I: Basic effects of rank reduction and kernel smoothing

First, we examined the effect of the nuclear-norm regularization in a simple simulation setting
corresponding to Fig. 1. The dataset consisted of a seven-dimensional (D = 7) binary time4

We used a Matlab implementation [43] of a limited-memory BFGS presented by Mark Schmidt, which isavailable at http://www.cs.ubc.ca/˜schmidtm/Software/minFunc.html

11

1
7
1
4 0.5
6
5
0
1
0.5
0
1
0.5
0
0
2

3

50

100 150
Time

200

Figure 3: Simulation study I: Three graphs corresponding to the basis vectors (left) and thetime-series of their coefficients (right).
series of length 200 (N = 200). We sampled y n ∈ {−1, 1} at each time step from an Ising modelp(y n ; θ n ) according to the joint probabilities directly calculated by Eq. (1); every θ n ∈ R28 wascreated byθ n = sn1 a1 + sn2 a2 + sn3 a3 ,(25)

with three basis vectors {a1 , a2 , a3 } and their coefficients {sn1 , sn2 , sn3 }, so that every θ n was
effectively constrained to be in a three-dimensional subspace, and the rank of Θ was at mostthree.
The dependency graphs corresponding to the three vectors ak (k ∈ {1, 2, 3}) are shown in
Fig. 3. Specifically, for each k, we set an edge-wise parameter akij at a positive constant of 0.5 if
the k-th graph has an edge between the nodes i and j, or set it at zero otherwise; we simply set
all the node-wise parameters aki (i = 1, 2, . . . , D) at zero. Figure 3 also shows the time-series of
their coefficients, sn1 , sn2 , and sn3 , which switched between 0 or 1 periodically with different cycle
lengths.
We ran the proposed algorithm with every combination of the following values of the tuning
parameters (w, λ, η): w ∈ {2, 4, 6, 8}, λ from zero and positive values such that log10 λ’s wereregularly spaced by 0.05 in [−4, −2], and η from positive values such that log10 η’s were regularly
spaced by 0.05 in [−3, −1]. We used the box kernel for kernel smoothing (see Fig. 2). We also
ran a baseline algorithm with η = 0 for comparison, which basically solves a problem equivalent
to that of KELLER [46, 28] as described in Section 4. In this experiment, the algorithm was(t) (t)(t)terminated simply when max{δz , δm } ≤  (= 10−5 ) was achieved, as in [55], where δz and(t)δm are the maximum absolute values of all of the elements in Z(t) − Z(t−1) and M(t) − M(t−1) ,respectively. The constant α was initially set at 10−3 , but in order to accelerate the final(t)convergence, it was multiplied by a factor of 1.5 at each iteration after reaching δm ≤ 0.8.
The evaluation was performed in a similar manner to that of binary classification. In other
12

0.95

AUC

0.9
0.85

w=2
w=4
w=6
w=8

0.8
0.75
0.7
-3

-2
log10 η-1

Figure 4: Simulation study I: The area under the ROC curve (AUC) versus log10 η (η, which is
the regularization coefficient for the nuclear norm). For each w, the AUC value with η = 0 is
depicted by a horizontal line.
n according to the final estimate Θb waswords, as described in Section 3.3.3, each edge weight θbij
examined to determine whether the weight was nonzero or exactly zero, corresponding to edge
existence or non-existence, respectively. Every such binary estimate was then compared withthe truth for all i < j and n. We quantified the performance by examining the area under the
ROC curve (AUC) for every choice of w and η, as shown in Fig. 4. This figure clearly showsthat for a wide range of η, the proposed method outperformed the previous method without thenuclear-norm regularization (η = 0).
Figure 5 also shows some examples of the ROC curves. For each w, the curves correspondto the best value of the AUC in Fig. 4 and the baseline algorithm of η = 0. The verticaland horizontal axes respectively denote the true positive and false positive rates (within all
correct and incorrect outputs, respectively), where positive refers to a nonzero weight (i.e., theexistence of an edge). Note that the ROC curve was plotted in a non-standard manner: rather
than varying the threshold for classification (which is lacking in the setting of the present study),
each single curve shows the variation of the classification performance according to the value ofthe regularization coefficient λ (i.e., the strength of sparsity). The point at the top-right cornercorresponds to λ = 0, and the other points to positive λ’s whose values increase as the bottom
left corner is approached.
These results showed that the rank reduction can improve the overall performance of edge
detection in terms of the AUC measure. We will use the same scheme also in the following
simulation studies (Sections 5.2, 5.3 and 5.4) to evaluate the average performance.
It should be noted that kernel smoothing may locally have a detrimental effect on edge
detection especially when the smoothing is strong (i.e., with large w). To illustrate this, we
examined a simple change-point problem in which the data were generated in the same manner
as above except that the true dependency structures (i.e., true Θn ) changed only at the 101-thtime point. The two graphs before and after the change point are illustrated at the top of Fig. 6.
Only the three edge weights changed from zero to nonzero (edges 1 and 3) or nonzero to zeron being set at 0.5. We set η as log η = −1.5, and λ as(edge 2), with every nonzero weight θij
10
log10 λ = −2.5 in Fig. 6 (a) and log10 λ = −2.3 in Fig. 6 (b). To clarify the effect of smoothing,13

True Positive Rate

1

0.5

0
0

True Positive Rate

ROC curve (w=4)

No nuclear norm
log10 η = -1.4
0.5
False Positive Rate
ROC curve (w=6)

1

0.5

0
0

No nuclear norm
log10 η = -1.65
0.5
False Positive Rate

1

0.5

0
0

1

True Positive Rate

True Positive Rate

ROC curve (w=2)

0.5
False Positive Rate
ROC curve (w=8)

1

1

0.5

0
0

1

No nuclear norm
log10 η = -1.6

No nuclear norm
log10 η = -1.55
0.5
False Positive Rate

1

Figure 5: Simulation study I: Estimation of graph structure improved by nuclear-norm regularization. The ROC curves were plotted by varying the strength of sparsity (regularization
coefficient of the `1 norm) λ, instead of varying the threshold for classification (which is lacking
in the proposed method). Each panel corresponds to a specific choice of the kernel width w.
The vertical and horizontal axes denote the true positive and false positive rates, respectively,
and the two curves show the results without the nuclear norm (η = 0) and with a positive η (η:the regularization coefficient for the nuclear norm), where only the result with the best η value(in terms of AUC shown in Fig. 4) is shown for clarity.

14

3
1

True edge

0.2

Estimated
0.1 changepoints
0
0.2
0.1
0
0.2
0.1
0
0

50

100

150

3

2

Edge 3 Edge 2 Edge 1

Edge 3 Edge 2 Edge 1

2

200

Time

1

0.1
0.05
0
0.1
0.05
0
0.1
0.05
0
0

50

100

150

200

Time
(b) Sparsity high: log10 λ = −2.3(a) Sparsity low: log10 λ = −2.5

Figure 6: Simulation study I: Strong kernel smoothing (w = 32, box kernel) possibly biases theedge detection around steep change points (while improving it in stationary regimes), dependingon the value of regularization coefficient: log10 λ = −2.5 in (a) and −2.3 in (b), and log10 η =−1.5 (both (a) and (b)). Top: True dependency graphs in the former and latter halves of thetotal duration N = 200, with positive edge weights fixed at 0.5. Bottom: Time series of theestimated edge weights (gray lines, for 20 randomly-generated datasets) corresponding to thethree edges indicated in the graphs. Distribution of the estimated change points are indicated byvertical solid (median) and dashed (10- and 90-percentiles) lines. Black horizontal bars indicate
the duration of the existence of true (nonzero) edges.
we set the kernel width relatively large, i.e., w = 32. In both figures, estimated time series ofn (for 20 randomly-generated datasets) are drawn by gray lines, andthe three edge weights θbij
the median and two percentiles (of 10% and 90%) of the detected change points (from zero tononzero in edges 1 and 3, and from nonzero to zero in edge 2) are indicated by the solid anddashed vertical lines, respectively. The black horizontal bars represent the time period during
which the true edge existed.
As is obvious in Fig. 6 (a), the relatively strong kernel smoothing enabled the accurate
identification of zero or nonzero within the stationary regimes, while it caused the estimatedweights (gray lines) to be wrongly nonzero around the true change points, so that the timing of
estimated change points (vertical solid/dashed lines) was biased. Strong kernel smoothing may
thus lead to locally inaccurate edge detection around a steep change point. As seen in Fig. 6(b), a stronger sparsity (or possibly another threshold on weights) could compensate for such abiased estimation of change points by shrinking the estimated weights toward zero. However,the combined effect should depend on different choices of tuning parameters, and is not easily
predictable.
However, such a detrimental effect of kernel smoothing would be minor on the whole due tothe relatively small fraction of change points among all the time points, so that the AUC can bea valid measure for evaluating overall edge detection performance. If the objective is to know
the precise timing of change points, results must be carefully interpreted. We discuss possible
solution to this issue in Section 7.

15

5.2

Simulation study II: Smooth transition with dependency patterns

Next, we evaluated the performance of the estimation of graph structures with a variety ofsynthetic time-varying networks. The objective of this experiment is to clarify the situations
where the rank reduction works efficiently.
Here, the number of variables was D = 12, and a total of N = 360 instances y n were
sampled from time-varying Ising models p(y n ; θ n ), where θ n ∈ R66 is the vector of edge-wiseparameters at the n-th time step; all the node-wise parameters were set to zero for simplicity.
Every parameter
vector θ n was given by a linear combination of K basis vectors accordingPKn k
to θ n =k=1 sk a . Each basis vector again defines a three-node clique (triangle) in thedependency graph, as in Section 5.1, where the K cliques were randomly chosen in each run
without any overlapped edges between cliques; every nonzero element of ak was fixed at one.
We created more realistic time courses of the network than those in the previous section by
generating the time-varying basis coefficients sni in a similar scheme to that used in [58]. The
network smoothly transitions between stationary regimes, and each dependency pattern appears
non-periodically with non-constant strength.
Figure 7 illustrates this scheme. First, we divided the 360 time points into 360/duration + 1
blocks, where duration means the length of the time period during which the structure ofa dependency graph was maintained. The first and the last blocks had half of the length ofthe others. In the first block, a fixed number num edges/3 of active (nonzero) coefficients was
randomly generated. As each clique has three edges, num edges means the number of total edges
in an instance graph. The nonzero values of these active coefficients were sampled uniformly
in the range of [2, 6]. To generate the next block, we randomly set half (rounded-up if not
an integer) of the active coefficients at zero, and activated the same number of other inactive
coefficients with their values sampled in the same way as in the first block. This procedure wasrepeated until the final block was generated. Finally, the transitions between active and inactive
coefficients were replaced by smooth sigmoidal curves.
Here, the values of duration and num edges respectively controlled the level of stationarity
and the level of sparsity of the generated time-varying networks. For each specification of thesevalues, we randomly created ten different time-series of time-varying networks. A single binary
time-series was generated from each of these by the corresponding Ising model. For each single
time-series, we ran the proposed algorithm with various combinations of the tuning parameters(w, λ, η) with the Epanechnikov kernel (see Fig. 2) and the stopping criterion described in
Section 3.3.3. We also ran the baseline algorithm with η = 0 for comparison.
Figure 8 shows how the AUC varied against η for different kernel widths w. Each panel
corresponds to a particular combination of duration and num edges. The number of underlying
dependency patterns (maximum rank) K was also varied as K = 7 (a) and K = 14 (b). For thelarger maximum rank of K = 14, the number of active coefficients (i.e., num edges/3) was also
set larger than that of K = 7. Note that when num edges = 3, only a single coefficient sni wasnonzero at each time, resembling the ordinary situation of regime switching. In every panel, the
AUC of the proposed method increased over that of the baseline method (log10 η → −∞) as
log10 η increased, until log10 η reached around −1.5; the increase was more obvious with smaller
kernel widths w. This tendency is more clearly seen in those panels corresponding to smaller
num edges and smaller maximum rank (i.e., K = 7). Note that sparser Θ tends to have a lower
rank. The rank reduction was thus more effective when the true rank was lower.
In this simulation, a particular dependency (regime) may repeatedly appear especially when
the network changes more frequently. Figure 8 thus confirmed our statement in Section 3.2
that the rank reduction works efficiently, especially when temporally-distant networks (relative
to the kernel width) can be similar to each other. This is most evident if we select the best
kernel widths (among the four specific ones) for duration = 5, 15 and 60 as w = 10, 30 and
16

60, respectively, with which the highest AUC value was achieved. The increase in AUC thenappeared to be large, moderate, and small when the duration was 5, 15, and 60, respectively.
The more significant improvement of the AUC for a smaller duration clearly indicates that the
proposed method effectively incorporated the similarity between temporally-distant networks
which could not be (re-)used by the local kernel smoothing. The same thing is also seen when
the kernel width w was set smaller than the best one (e.g., with w = 3 or 10) for the most casesof duration = 60.
The conjoint effect of the temporal smoothing and the rank reduction in our method naturally
raises a question: How does the proposed method work solely with the rank reduction but
without temporal smoothing? Figure 9 examines this issue. This again shows the AUC versus
log10 η in the same setting as that of Fig. 10 (b) (i.e., K = 14), but in particular with w = 1 (i.e.,
the kernel smoothing was not effective). Clearly, in every panel, even the highest AUC was notcomparable to the AUC obtained by the proposed or the baseline method with kernel smoothing(see Fig. 10 (b)). This implies that the temporal smoothing is essential for successful estimationof time-varying dependency graphs; the rank reduction then further improves its performancein the particular situations considered here.

5.3

Simulation study III: No explicit dependency patterns

To examine a limitation of the proposed method, we also studied how the proposed method works
when the network does not a priori have a limited number of communities or regimes. Here,n without explicitly using any dependency
we directly generated time-varying edge weights θijpatterns ak . In other words, we replaced the three-node cliques ak by all the edges (K = 66)
and generated their coefficients snk in the same manner as in Section 5.2. Note that the meaningof the two simulation parameters duration and num edges does not change at all.
Figure 10 shows the result in a similar manner to that in Figs. 8 and 12. Here, we examined
the four levels of sparsity, num edges = 3, 6, 9 and 12, and the three levels of stationarity,
duration = 5, 15, and 60. In this figure, the increase in AUC by rank reduction disappeared inmany cases. This is reasonable because the parameter matrix is now not necessarily of low rank.
However, a slight increase in AUC over that of the baseline method in some kernel widths w isstill seen especially in those panels of small num edges (say 3 or 6), where the true rank could beapproximately low. Most notably, the increase in AUC was rather large when duration was 60
with small kernel width w = 3 or 10 (Interestingly, the top-right panel of Fig. 10 is quite similar
to that of Fig. 8). This was probably because the number of stationary regimes was limited aposteriori in these cases, as a consequence of the small number of transitions during the entire
time-series, as shown intuitively by Fig. 11. The rank reduction was thus sometimes helpful
to improve the estimation, even when the network does not a priori have a limited number ofcommunities or regimes.

5.4

Simulation study IV: Small sample size (short time-series)

A difficult situation arises in practice because the total sample size, i.e., the length of the timeseries, is smaller than the total number of parameters. In order to see how the proposed method
works in this situation, we repeated the simulation study II (Section 5.2) with a reduced sample
size of N = 60, which was smaller than the number of edge-wise parameters.
The result is shown in Fig. 12 in the same format as in Fig. 8. Note that duration = 60 now
means that the transition of the network occurred only once during the whole period. Every plot
exhibited larger variance as expected, while the increase in AUC with increasing log10 η is still
evident in some cases. For example, in the upper three panels of duration = 5, the AUC clearly
increased over that of the baseline method within an appropriate range of log10 η, particularly
17with kernel width of w = 3 or 10. Similarly, in the upper three panels of duration = 60,
the AUC increased particularly with kernel width of w = 10 or 30 (but not of w = 3). Thus,
even if the time-series length is quite short, the rank reduction can be useful, at least when theestimation variance is not too large, as in the bottom-most panels of Fig. 12. The choice of
kernel width seems to be important for the method to be effective.

5.5

Computation time

Here, we compared the computation time between the ADMM algorithms and the regressionbased KELLER [46, 28]. All the algorithms were implemented on Matlab 7.14 and run on
a Linux computer with 4 cores, 3.30-GHz CPU and 126-GB RAM. We modified the original
implementation of KELLER 5 in order to use the same smoothing kernel and the same coefficienton the `1 -norm (according to Eqs. (23) and (24)). We measured the elapsed CPU time atevery iteration of the ADMM algorithm on the ten datasets above of N = 360 (Section 5.2),
with K = 14, num edges = 6 and duration = 5. Throughout this experiment and the othersimulation studies described so far, we observed no runs that did not converge.
Figure 13 shows the values of the objective function versus the CPU time for three different
λ’s. Note that the baseline ADMM (η = 0) and KELLER minimize the same objective function,
while the low-rank ADMM includes the additional term ηkΘk∗ , with η being specifically set as
0.01. Each point for KELLER also indicates the objective value versus the total CPU time of asingle run (i.e. to solve the N D logistic regression problems), obtained with varying termination
conditions in the logistic regression. Note that the ADMM does not necessarily decrease theobjective value monotonically (as in the case of log10 λ = −2.5). The result shows that KELLER
achieved almost the minimum objective value in about 20-50 seconds in any case, while our
ADMM algorithms took more than 100 seconds for convergence. Thus, if the rank reduction is
not necessary, the computational efficiency of KELLER is favorable in practice; whereas, if therank reduction is desirable, the ADMM (η > 0) would be a better option.
In this experiment, the mean CPU time (with standard deviation) per iteration of the ADMM
(η = 0.01) was 0.82 ± 0.13 [sec], which was slightly heavier than 0.78 ± 0.09 [sec] of ADMM
(η = 0) due to the additional computation of SVD for non-zero η. On the other hand, Fig. 13
shows that the ADMM (η = 0.01) converged earlier than the ADMM (η = 0) did. This implies
that the ADMM (η = 0.01) needed fewer iterations to converge than that by the ADMM (η = 0).
Figure 13 also suggests that with larger λ, the ADMM (η = 0) tended to converge earlier.

6

Demonstration with a real-world network

Here, we demonstrate the effect of rank reduction when estimating time-varying MNs for areal-world social network: the network of one hundred senators during the 109th Congress of
the U.S. Senate exhibited in their recorded votes during this congress 6 . This dataset has been
analyzed in previous studies [3, 1, 28]. We followed [3] in particular for preparing the data.
In this dataset, each binary instance y n means the record of a single roll-call vote by the 100
senators (D = 100) for the n-th bill. Each vote by the i-th senator was recorded as yin = −1 for“nay” and yin = 1 for “yea” with missing votes simply treated as “nay” [3]. The total numberof the bills, i.e., the length of the time-series, was N = 645. See [3] for more details about the
preprocessing.
5
6

Matlab code was available on http://cogito-b.ml.cmu.edu/keller/downloads.html
The summary of roll-call votes is available from the website (http://www.senate.gov/)

18

6.1

Comparison of predictive performanceFirst, we evaluated the time-varying networks obtained by the proposed and the baseline methods in terms of their abilities for predicting unobserved data, instead of their performance inrecovering graph structures. This was done because the true dependency graph between the
100 senators is not available. To this end, estimation (training) was performed in each run
based only on the training dataset consisting of 516 (4/5 of the total) randomly selected instances, and a score of predictive performance was computed for the test dataset consisting ofb The predictability score was specifithe other 129 instances, using the obtained estimate Θ.
callytest log-pseudolikelihood, computed by summing-up the log-pseudolikelihood
P givenn bynthebn ) over the test instances. This evaluates the performance of conditional prelog
p(y
|y
;
θ
i
i
\i
diction on the outcome of each node given those of the other nodes.
During training, the values of the smoothing kernel at the test instances were explicitly
replaced with zeros in order to effectively exclude them from the training dataset. We used theEpanechnikov kernel (see Fig. 2), and employed the stopping criterion described in Section 3.3.3.
For comparison, we also trained and tested the static sparse MN (Section 2.3) in the same manner
as above.
Figure 14 shows the test log-pseudolikelihood versus log10 λ for some kernel widths w. Eachpanel shows the results for four specific values of η at a particular w, and also shows the results
for the baseline method of η = 0. It also indicates the best performance achieved by the static
MN, which was obtained at the best value of λ that maximized the test log-pseudolikelihood(Fig. 15). The rank reduction with appropriate combinations of η and λ greatly improved theprediction performance especially when the kernel width was relatively small. The performanceby the proposed method was at a high level already at w = 30, especially with good choices ofη and λ (say, log10 η = −1.75 and log10 λ ∈ [−5, −4.5]). The larger kernel widths then further
improved the performance, at least until the width reached to w = 90. On the other hand, theperformance of the baseline method at w = 30, with optimized λ, was even lower than that of
the static network; this low performance was improved as the kernel width increased.
The social network considered here is greatly expected to have a community structure because of the two large political parties (i.e., Democratic and Republican), probably with several(transient) sub-communities within or across them. It is also reasonably expected that the
network does not change very frequently, because it likely reflects the political positions of thesenators. Thus, the result shown here is reasonable because in such a case, the estimation canbe greatly improved by rank reduction, especially when the kernel width is relatively small, as
actually seen in the simulation studies in Section 5.

6.2

Effect of rank reduction in estimated network

We obtained an intuitive understanding of how the rank reduction affects the quality of estimatednetwork by again analyzing the whole dataset of the 645 instances by the proposed method with
log10 η = −1.75 and by the baseline method with η = 0, both with the kernel width of w = 70.
The value of λ was commonly set at log10 λ = −4.5, so that the two methods had a comparable
level of predictability, as seen in Fig. 14.
Figure 16 shows the time courses of 50 randomly selected edge weights. The proposedn over time, which were
method clearly reduced undesirable fluctuations of estimated weights θij
unfavorably observed in the baseline method. This additional smoothing effect was particularly
due to the rank reduction, as the two methods used the same kernel width: The underlying
dependency patterns ak effectively divide different edges into groups and allow each group toshow a consistent time course described by its common coefficient snk .
Figure 16 also implies that most of the edge weights at each time point took non-zero values,
19

in particular in the result by the proposed method. Actually, we observed that maximizing thepredictability score in this dataset tended to produce rather dense dependency graphs. This
is not strange by itself, because the senators’ votes likely reflect some typical political stances,
depending especially on the political party to which they belong, so that their network should
not consist of sparse local connections seen in social networks of friendship.

6.3

Implication for subsequent analysis

In practice, determining useful knowledge from hundreds of large-scale dependency graphs is
not easy. This is especially the case in the context of data mining, where we often do nothave much prior knowledge or solid hypotheses about the data. Subsequent analysis of thesedependency graphs will then be important for summarizing relevant information in order todisclose embedded knowledge.
Here, we demonstrate a simple case of this type of a post-analysis, and show that the improved quality of the dependency graphs obtained by the proposed method is actually beneficial
at this stage. Specifically, we consider partitioning a graph into several clusters, each represented by an influential node called an exemplar. This is the problem considered in [17] forwhich an effective algorithm, the affinity propagation clustering (APC) [17], has been proposed.
We applied this APC to every instance of the dependency graph in order to examine the changes
in the exemplars and the clusters over time. The similarity matrix between any pair of nodes
to which the APC was applied was defined as a matrix where both (i, j) and (j, i) entries are
n . A nice property of APC is that it does not suffer fromgiven by the estimated edge weight θbij
the issue of local optima, and it automatically determines the number of clusters by specifying
a positive “preference” parameter [17]. This parameter was set at 0.01 based on a preliminary
analysis.
For comparison, the APC was performed on both the estimates by the proposed and theb η and Θ
b 0 , respectively. The kernel width was set at w = 70,
baseline methods, denoted as Θand the regularization coefficients were specifically set at log10 λ = −5 and log10 η = −1.75
for the proposed method, and log10 λ = −4.5 (and η = 0) for the baseline method, so as
to roughly maximize the test log-pseudolikelihood in each method (Fig. 14). Note that theb 0 was more
smaller predictability score of the baseline method indicates that the estimate Θlikely overfitted to the training dataset.
Figure 17 shows the temporal changes of cluster members when one of the six nodes (senators)
was chosen by APC as an exemplar. The six exemplars displayed here were the most frequentlyb η . The cluster members displayed herechosen ones (in all the 645 time points) based on Θ
for each exemplar were those belonging to the exemplar at least five times, which was alsob η . The same exemplars and cluster members are also displayed for thedetermined based on Θb 0 for the sake of comparison. As seen in this figure, the overall patterns ofresult based on Θb η and Θ
b 0 , indicated by “Proposed” and “No nuclear norm” at the top of
the clustering in Θeach block, respectively, were quite similar to each other. However, the APC recovered transient
clusters more stably in “Proposed” than in “No nuclear norm.” For example, the proposed
method indicated that Senator Sununu as the exemplar had a cluster with Senator Gregg and
Senator Kyl from the time points from 1 to about 300. However, this stable cluster is not clearlyseen in “No nuclear norm” as it is fragmented into many discontinuous segments. Thus, the rank
reduction in the estimation stage can improve the quality of the result of subsequent analysis,
from which one may draw clearer and more specific speculation into human relationships.

20

7

Discussion

We proposed a method for learning time-varying MNs based on the assumption that the rank ofparameter matrix is relatively low, in addition to the sparsity and temporal smoothness of the
network structure. The problem was formulated as a non-smooth convex optimization problem
where the objective function (i.e., kernel-weighted log-pseudolikelihood) was jointly regularizedwith both `1 and nuclear norms. We proposed to solve this problem by using the ADMM, which
led to a simple first-order optimization algorithm. The main contributions can be summarized
as follows: 1) we proposed a novel approach to learning time-varying MNs with a new type of
application of the joint `1 - and nuclear-norm regularization, and 2) we empirically showed thatour proposed method can outperform existing methods in terms of the estimation performanceas well as the predictability, as recapitulated below. Although we did not give any theoretical
performance guarantee, the empirical results successfully demonstrated the applicability of the
proposed approach to various situations.
In the simulation studies, we examined the effect of rank reduction on estimation of timevarying dependency graphs, in cases where 1) the network exhibits a relatively small number of
dependency patterns (Secs. 5.2 and 5.4), resembling community structures or regime switching,and 2) the network a priori has no limited number of dependency patterns (Section 5.3). In case
1, the rank reduction was particularly effective when the network changes frequently, so that
temporally-distant networks may have similar dependency, which has not been dealt with very
well by local kernel smoothing; the effect was rather limited if the time-series was very short,
but was still evident in some situations with appropriate kernel widths. In case 2, the effect ofrank reduction became weaker in most situations as expected, while it was surprisingly helpfulwhen the number of dependency patterns was limited a posteriori; namely, when the stationary
regimes in the entire time-series were relatively small due to the stability of the network. The
rank reduction is therefore useful for dealing with datasets with a small number of change points,
even when the network neither exhibits community structure nor repeats any previous regimes.
We also demonstrated the effect of rank reduction with a real-world dataset of US Senate voting records (Section 6). The rank reduction improved the predictability of the obtained networkin terms of the test pseudolikelihood, and also led to more interpretable results in subsequent
analysis of the collection of estimated dependency graphs. In particular, we demonstrated theuse of the APC for clustering the nodes (senators).
As seen in Section 6, the test pseudolikelihood can be used for choosing η jointly with λ,which evaluates the predictability of the model and, in principle, can avoid overfitting to thetraining dataset. The kernel width w may also be automatically selected in a similar manner,
while the choice should also reflect our prior knowledge about the time scale of the phenomenaof interest. Thus, in practice, we recommend interactive selection of w from a relatively small
set of candidates around the time scale of interest. A joint search for w, λ and η could be
conducted by an heuristic greedy scheme for ease of computation. For example, w, λ, and η
may be sequentially determined with λ and η initialized by zeros. Alternatively, the BIC-like
information criterion developed in the previous study [28] could also be modified to be suitable
for the method proposed here, which will further reduce the computational burden.
Our present study successfully demonstrated the feasibility and empirical performance of thejoint `1 - and nuclear-norm regularization method for learning time-varying networks throughextensive simulation studies as well as a real-world data analysis. To further clarify the applicability of our approach beyond the situations therein, we will need theoretical analyses of
its statistical performance, but such analyses are beyond the scope of the current study. Yet,
recent theoretical studies [41, 35] of joint `1 - and nuclear-norm regularization give insights into
our approach. In these studies, error bounds on the matrix estimation were theoretically given21

in specific settings of matrix recovery [41] and multi-task regressions [35] both with quadratic
loss functions. As seen in Section 4, our loss function is closely related to a collection of manylogistic regressions, similar to the multi-task setting in [35]. Hence, although kernel smoothingmight need a special treatment, the approach in these studies could be extended to the problem
considered here by incorporating more general loss functions [41]. These error bounds basically
imply that true low-rank matrices can be correctly estimated when the true matrix is actuallysparse and of low rank and the regularization coefficients are set large enough compared to thetrue cardinality (i.e., number of nonzeros) and rank. We expect that similar results will qualitatively hold for our problem; detailed theoretical analyses will clarify the general conditions forthe proposed method to success or to fail. We leave this topic as an open issue for future study.
We finally discuss possible future extensions and other promising applications of our proposed framework. First, we focused on the fundamental application to the Ising model, but the
proposed approach is also applicable to other time-varying MNs possibly with slight modifications. For example, the GGM is a real-valued counterpart of the Ising model, where each edge is
given a scalar-valued (edge-wise) parameter defining the analog graph structure. The use of the`1 -norm, which induces element-wise sparsity, is then directly applicable. However, for general
MNs, each edge may be associated with a vector-valued parameter, or higher-order cliques may
even be involved. In these cases, the simple `1 -norm regularization can be replaced by otheradvanced techniques of sparse estimation, such as the group `1 [54, 36] or the hierarchical sparsity [56, 27, 44]. The ADMM is then slightly modified so that Eq. (18b) is replaced with anappropriate operation corresponding to the alternative regularization term.
Second, kernel smoothing may be replaced by other techniques if one needs to precisely
know the timing of abrupt network changes (see also Section 5.1). Some previous studies oftime-varying MNs [1, 28, 29] actually used a technique called total variation smoothing to allow thePnetwork parameters to be piecewise constants. It additionally introduces regularizationkθ n − θ n−1 k1 into the objective function with a non-smoothed version of empiricalterms Nn=2Pn nloss, (1/N ) N
n=1 l(y , θ ). The `1 -norm of parameter differences allows each element of the parameter to have jumps at a few time points, and thus it may be more suitable for situations likeregime switching. The additional low-rank regularization in our method then further exploits
the potential similarity between temporally distant pairs of networks, as well as the similarity
between pairs of the time-series of network parameters, which cannot be achieved by solely the
total variation smoothing. Note that ADMM then can be applied by including an additional
updating step (e.g.,
the additional regularization terms. A more conventional smoothP of Z4n) for n−1
ing regularizer, Nkθ−θk22 , is also available without complicating the algorithm (sincen=2
it is differentiable), but it may not be conceptually very different from kernel smoothing as it
introduces a similar smoothness into the network parameters.
Third, an interesting extension of the proposed time-varying MNs would be to use of thesame technique for auto-regressive models, relaxing the basic assumption in the Ising modelthat y n is independent of all past observations y 1 , y 2 , . . . , y n−1 given θ n . For example, we could
assume that y n depends on the past τ observations (i.e., y n ∼ p(y n |y n−1 , . . . , y n−τ ; θ n )). The
definition of θ n in the present study could then be simply modified to incorporate additionalterms of the lagged effects, such as yin yjn−m (1 ≤ m ≤ τ ). This only changes the number of rows
in Θ, and hence the ADMM algorithm will be directly applicable.
Finally, other than analyzing networks of social interactions as successfully demonstrated
with the US Senate data, the proposed method can also be useful in the context of analyzing
biological phenomena (e.g., gene-regulatory networks), as considered in the previous study [46].
Another promising application of the time-varying MNs would be the analysis of multi-channel
measurements of the brain activity, such as neuronal spike trains [45] or other brain imaging modalities like functional magnetic resonance imaging (fMRI) or magnetoencephalography22

(MEG). Interestingly, time-varying functional connectivity (dependency) between cortical regions has recently been a very active research target in the neuroscience and brain imaging
communities [24, 33]. We will apply the method proposed here to these emerging new fields inour future studies.
Acknowledgements This research was supported by a contract entitled “Novel and innovative R&D making use of brain structures” with the Ministry of Internal Affairs and Communications, Japan, a contract entitled “Brain mapping by integrated neurotechnologies for diseasestudies” with the Japan Agency for Medical Research and Development, and JSPS KAKENHI
Grant Numbers 24300114, 25730155. A.H. was supported by the Academy of Finland, Centreof-Excellence in Inverse Problems Research.

References[1] A. Ahmed and E. P. Xing. Recovering time-varying networks of dependencies in social andbiological studies. Proceedings of the National Academy of Sciences of the United States ofAmerica, 106(29):11878–11883, 2009.
[2] P. Bachmann and D. Precup. Improved estimation in time varying models. In Proceedingsof the 29th International Conference on Machine Learning (ICML’12), 2012.
[3] O. Banerjee, L. E. Ghaoui, and A. d’Aspremont. Model selection through sparse maximum
likelihood estimation for multivariate Gaussian or binary data. Journal of Machine Learning
Research, 9:485–516, 2008.
[4] D. P. Bertsekas and J. N. Tsitsiklis. Parallel and distributed computation: Numericalmethods. Prentice-Hall, Inc., 1989.
[5] J. Besag. Statistical analysis of non-lattice data. The Statistician, 24(3):179–195, 1975.
[6] J. Besag. Efficiency of pseudo-likelihood estimation for simple Gaussian fields. Biometrika,
64:616–618, 1977.
[7] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed optimization andstatistical learning via the alternating direction method of multipliers. Foundations andTrends in Machine Learning, 3(11):1–122, 2011.
[8] S. Boyd and L. Vandenberghe. Convex Optimization. Cambridge University Press, 2004.
[9] J.-F. Cai, E. J. Candès, and Z. Shen. A singular value thresholding algorithm for matrix
completion. SIAM Journal on Optimization, 20(4):1956–1982, 2010.
[10] C. M. Carvalho and M. West. Dynamic matrix-variate graphical models. Bayesian Analysis,2(1):69–97, 2007.
[11] S. S. Chen, D. L. Donoho, and M. A. Saunders. Atomic decomposition by basis pursuit.
SIAM J. Sci. Comput., 20:33–61, 1998.
[12] I. Daubechies, M. Defrise, and C. De Mol. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Communications on Pure and Applied
Mathematics, 57(11):1413–1457, 2004.
[13] X. V. Doan and S. Vavasis. Finding approximately rank-one submatrices with the nuclear
norm and `1 -norm. SIAM Journal on Optimization, 23(4):2502–2540, 2013.
23

[14] J. Eckstein and D. Bertsekas. On the Douglas-Rachford splitting method and the proximal
point algorithm for maximal monotone operators. Mathematical Programming, 5:293–318,
1992.
[15] M. Fazel, H. Hindi, and S. Boyd. Rank minimization and applications in system theory. In
Proceedings American Control Conference, pages 3273–3278, 2004.
[16] M. Figueiredo and J. Bioucas-Dias. Restoration of Poissonian images using alternating
direction optimization. IEEE Transactions on Image Processing, 19(12):3133–3145, 2010.
[17] B. J. Frey and D. Dueck. Clustering by passing messages between data points. Science,
315(5814):972–976, 2007.
[18] J. Friedman, T. Hastie, and R. Tibshirani. Sparse inverse covariance estimation with thegraphical lasso. Biostat, 9(3):432–441, 2008.
[19] W. J. Fu. Penalized regressions: the Bridge versus the Lasso. Journal of Computational
and Graphical Statistics, 7(3):397–416, 1998.
[20] J. Guo, E. Levina, G. Michailidis, and J. Zhu. Joint structure estimation of markov networks. Submitted, 2010. Available from http://www.stat.lsa.umich.edu/˜guojian/.
[21] B. He and X. Yuan. On the o(1/n) convergence rate of the douglasrachford alternating
direction method. SIAM Journal on Numerical Analysis, 50(2):700–709, 2012.
[22] J. Hirayama, A. Hyvärinen, and S. Ishii. Sparse and low-rank estimation of time-varying
markov networks with alternating direction method of multipliers. In International Conference on Neural Information Processing (ICONIP’10), Lecture Notes in Computer Science,
volume 6443, pages 371–379, 2010.
[23] H. Höfling and R. Tibshirani. Estimation of sparse binary pairwise Markov networks using
pseudo-likelihoods. Journal of Machine Learning Research, 10:883–906, 2009.
[24] R. M. Hutchison et al. Dynamic functional connectivity: promise, issues, and interpretations. NeuroImage, 80:360–378, 2013.
[25] A. Hyvärinen. Consistency of pseudolikelihood estimation of fully visible Boltzmann machines. Neural Computation, 18(10):2283–2292, 2006.
[26] Aapo Hyvärinen. Sparse code shrinkage: Denoising of nongaussian data by maximum
likelihood estimation. Neural Computation, 11(7):1739–1768, 1999.
[27] R. Jenatton, J. Mairal, G. Obozinski, and F. Bach. Proximal methods for sparse hierarchicaldictionary learning. In Proceedings of the International Conference on Machine Learning(ICML), 2010.
[28] M. Kolar, L. Song, A. Ahmed, and E. P. Xing. Estimating time-varying networks. Annalsof Applied Statistics, 4(1):94–123, 2010.
[29] M. Kolar, L. Song, and E. P. Xing. Sparsistent learning of varying-coefficient models
with structural changes. In Proceedings of the 23rd Neural Information Processing Systems(NIPS’09), 2009.
[30] D. Koller and N. Friedman. Probabilistic Graphical Models: Principles and Techniques.
The MIT Press, Cambridge, MA, 2009.
24

[31] S. L. Lauritzen. Graphical Models. Clarendon Press, Oxford, 1996.
[32] S.-I. Lee, V. Ganapathi, and D. Koller. Efficient structure learning of Markov networks
using L1-regularization. In Advances in Neural Information Processing Systems (NIPS
2006), 2007.
[33] N. Leonardi et al. Principal components of functional connectivity: A new approach to
study dynamic brain connectivity during rest. NeuroImage, 83:937–950, 2013.
[34] C. Loader. Local Regression and Likelihood. Springer, 1999.
[35] S. Mei, B. Cao, and J. Sun. Encoding low-rank and sparse structures simultaneously inmulti-task learning. In Advances in Neural Information Processing Systems (NIPS), 2012.
[36] L. Meier, S. van de Geer, and P. Bühlmann. The group lasso for logistic regression. Journalof the Royal Statistical Society B, 70:53–71, 2008.
[37] N. Meinshausen, P. Bühlmann, and E. Zürich. High dimensional graphs and variable selection with the Lasso. Annals of Statistics, 34:1436–1462, 2006.
[38] M. E. J. Newman. Modularity and community structure in networks. Proceedings of the
National Academy of Sciences of the United States of America, 103(23):8577–8582, 2006.
[39] J. Nocedal and S. J. Wright. Numerical Optimization. Springer Series in Operations Research. Springer Verlag, 1999.
[40] J. Pearl. Probabilistic reasoning in intelligent systems: networks of plausible inference.
Morgan Kaufmann Publishers Inc., 1988.
[41] E. Richard, P.-A. Savalle, and N. Vayatis. Estimation of simultaneously sparse and low
rank matrices. In Proceedings of the 29th International Conference on Machine Learning(ICML’12), 2012.
[42] G. Rocha, P. Zhao, and B. Yu. A path following algorithm for sparse pseudo-likelihood
inverse covariance estimation (SPLICE). Technical Report 759, Statistics Department, UC
Berkeley, 2008.
[43] M. Schmidt. minFunc: unconstrained differentiable multivariate optimization in Matlab.
http://www.cs.ubc.ca/˜schmidtm/Software/minFunc.html, 2005.
[44] M. Schmidt and K. Murphy. Convex structure learning in log-linear models: Beyond pairwise potentials. In Proceedings of the Thirteenth International Conference on Artificial
Intelligence and Statistics, volume 9, pages 709–716, 2010.
[45] H. Shimazaki, S. Amari, E. N. Brown, and S. Grün. State-space analysis of time-varying
higher-order spike correlation for multiple neural spike train data. PLoS Comput Biol,
8(3):e1002385, 2012.
[46] L. Song, M. Kolar, and E. P. Xing. KELLER: estimating time-varying interactions between
genes. Bioinformatics, 25(12):i128–i136, 2009.
[47] N. Srebro, J. Rennie, and T. Jaakkola. Maximum margin matrix factorization. In Advancesin Neural Information Processing Systems (NIPS), volume 17, 2005.
[48] R. Tibshirani. Regression shrinkage and selection via the lasso. J. Royal. Statist. Soc B.,
58(1):267–288, 1996.
25

[49] K. C. Toh and S. W. Yun. An accelerated proximal gradient algorithm for nuclear norm
regularized least squares problems. Pacific J. Optimization, 6:615–640, 2010.
[50] R. Tomioka and M. Sugiyama. Dual augmented lagrangian method for efficient sparsereconstruction. IEEE Signal Proccesing Letters, 16(12):1067–1070, 2009.
[51] R. Tomioka, T. Suzuki, M. Sugiyama, and H. Kashima. A fast augmented lagrangian algorithm for learning low-rank matrice. In Proc. of the 27 th Annual International Conference
on Machine Learning (ICML 2010), 2010.
[52] M. J. Wainwright, P. Ravikumar, and J. D. Lafferty. High-dimensional graphical model
selection using `1 -regularized logistic regression. In Advances in Neural Information Processing Systems, volume 19, pages 1465–1472, Cambridge, MA., 2007. MIT Press.
[53] R. Yoshida, S. Imoto, and T. Higuchi. Estimating time-dependent gene networks from timeseries microarray data by dynamic linear models with markov switching. In Proc IEEE
Comput Syst Bioinform Conf., pages 289–298, 2005.
[54] M. Yuan and Y. Lin. Model selection and estimation in regression with grouped variables.
Journal of the Royal Statistical Society, Series B, 68(1):49–67, 2006.
[55] X. M. Yuan. Alternating direction methods for covariance selection models. Journal of
Scientific Computing, 51(2):261–273, 2012.
[56] P. Zhao, G. Rocha, and B. Yu. The composite absolute penalties family for grouped andhierarchical variable selection. Ann. Statist., 37(6A):3468–3497, 2009.
[57] K. Zhou, H. Zha, and L. Song. Learning social infectivity in sparse low-rank networks using
multi-dimensional Hawkes processes. In Proceedings of the 16th International Conferenceon Artificial Intelligence and Statistics (AISTATS’13), pages 641–649, 2013.
[58] S. Zhou, J. Lafferty, and L. Wasserman. Time varying undirected graphs. Machine Learning,80(2–3):295–319, 2010.

26

5

Coefficients

0
5
0
5
0
5
0
5
0
5
0
5
0
1

61

121

181

241

301

361

Time
Figure 7: Simulation study II: Illustration of the scheme generatingsmoothly-varying basisPKkn kcoefficients sni , from which the parameter vector is given by θ n =k=1 sk a where a is arandomly generated pattern vector corresponding to a 3-node clique (triangle) as in Fig. 3.
This illustration corresponds to K = 7 (number of underlying patterns ak ), num edges = 6(number of nonzero edges at each time point, i.e., pairs of coefficients sni being nonzero at atime), duration = 15 (length of near-stationary blocks), and N = 360 (length of time series; the361-th time point is included just for illustration). In each block segmented by vertical dottedlines, two selected coefficients are randomly given positive values which have been sampled
uniformly in the range of [2, 6] (white circle), while the rest are set to zero (black circle). The
piecewise constant values (dashed line) are then smoothed by replacing the transition periods
with sigmoidal curves (solid line). See the text for more details.

27duration

5

15

60

0.9
0.9
0.8

0.8
0.8

0.7

0.7

0.9
0.9

AUC

3

num edges
6

AUC

0.9

0.9

0.8
0.8
−3

−2

−1

0.8
−3

log10 η

−2

−1

−3

log10 η−2

−1

log10 η
w=60

w=30

w=10

w=3

(a) K = 7

duration

5

15

60

0.9

AUC

0.9
0.8

0.7
0.7

0.8

0.6
0.8

AUC

12

num edges

6

0.9
0.8

0.8

0.7
0.7
0.7

0.6

−3

−2

−1

0.6

−3

log10 η

−2

−1

−3

log10 η−2

−1

log10 η
w=60

w=30

w=10

w=3

(b) K = 14

Figure 8: Simulation study II: AUC versus the regularization coefficient η of the nuclear norm.
Each panel corresponds to a specific combination of num edges (level of sparsity) and duration(level of stationarity) to generate the true network, indicated respectively on the left and on the
top of the panels. The maximum rank of the true parameter matrix Θ was set at K = 7 (a) or
K = 14 (b). In each panel, the isolated markers on the left of the vertical solid line show theaverage AUC over ten runs by the baseline method (η = 0) for different w’s, with the error bars
indicating the standard deviation. On the right of the vertical solid line, the average AUC by
the proposed method (η > 0) versus log10 η is plotted, with the standard deviation indicated by
the shaded regions. The different markers correspond to different kernel widths shown in thelegend at the bottom.
28

5

AUC
AUC

num edges

6

0.7

12

duration
15

0.7

0.6

0.6

0.6

0.5

0.5

0.5

0.4

0.4

0.4

0.7

0.7

0.7

0.6

0.6

0.6

0.5

0.5

0.5

0.4

0.4

0.4

−3

−2

−1

−3

−2

log10 η

60

0.7

−1

log10 η

−3

−2

−1

log10 ηw=1

Figure 9: Simulation study II: AUC versus the regularization coefficient η of the nuclear norm,
when the kernel smoothing was disabled, i.e., w = 1. The maximum rank was K = 14. See also
the caption of Fig. 8.

29

duration
5

0.9

15
0.9

AUC

0.8

3

60
0.9

0.8

0.7

0.7

0.6

0.8

0.6

0.5

0.9

AUC

0.9
0.8

0.7
0.8

0.7

0.6

0.6
0.8
0.7

AUC

9

num edges

6

0.8

0.7
0.7
0.6
0.6
0.6
0.5
0.7

AUC

12

0.6
0.6
0.6

0.5

−3

−2

−1

−3

−2

log10 η−1

−3

log10 η−2

−1

log10 η
w=60

w=30

w=10

w=3

Figure 10: Simulation study III: AUC versus the regularization coefficient η of the nuclear
norm. In this simulation, no explicit dependency patterns were used for generating data, withthe weight of each edge separately generated by the scheme of Fig. 7. See also the caption ofFig. 8.

30

Edge index (sorted)

5
60
40

40

0
1

60
20

30
30
20

20

30

15

15

20

10

10

10

0
360 1

0
360 1

5
0
360 1

120
15
10
5
0
360 1

360

Time (sample) index
Figure 11: Simulation study III: Examples of time-varying edge weights in different levels of
stationarity, duration, indicated by the number on the top of each panel. In each panel, thehorizontal axis denotes the time steps and the vertical axis denotes the edge indices (sorted, sothat an edge activated earlier with a nonzero weight had a smaller index). The varying thicknessof each horizontal line indicates the relative magnitude of the corresponding edge weight at eachtime point. Only the edges whose weights became nonzero at least once are displayed.

31duration
5

15
0.9

0.9

0.8
0.7

0.8

0.7

0.7

0.6

0.6

0.9

AUC

3

num edges
6

AUC

0.8

60

1

0.9

0.8

0.7

−3

−2

0.8

0.8

0.7

0.7

−1

−3

log10 η

−2

−1

−3

log10 η−2

−1

log10 η
w=60

w=30

w=10

w=3

(a) K = 7

duration
5

15
0.9

AUC

0.9

0.8

0.8

0.7
0.7

0.7
0.8

0.7
0.7

AUC

num edges

6

0.8

12

60

0.7
0.6

0.6
0.6
−3

−2

−1

−3

log10 η

−2

−1

−3

log10 η−2

−1

log10 η
w=60

w=30

w=10

w=3

(b) K = 14

Figure 12: Simulation study IV: AUC versus the regularization coefficient η of the nuclear norm,
with a reduced sample size of N = 60. The simulation setting except for the sample size was
the same as that of simulation study II. See also the caption of Fig. 8.

32

log10 λ = -2.5log10 λ = -3

log10 λ = -4

9

Objective

ADMM (η=0.01)
ADMM (η=0)
KELLER

7.5

7

8.5

6

7

5

8

6.5
4

7.5

6

3
10

100

CPU time [sec]

10

100

CPU time [sec]

10

100

CPU time [sec]Figure 13: Comparison of computation time between the ADMM algorithms (η = 0 and 0.01)
and the KELLER [46, 28] based on `1 -regularized (weighted) logistic regression. Solid (η = 0)
and dashed (η = 0.01) lines without markers indicate the objective value versus the elapsed CPU
time (log scale) at every iteration of the ADMM algorithm for ten different runs (for different
datasets of N = 360, K = 14, num edges = 6 and duration = 5; see Section 5.2 for details of
the data). Note that the baseline (η = 0) and KELLER minimize the same objective function,
while the low-rank ADMM includes the additional term ηkΘk∗ , with η being specifically set as
0.01. The three panels correspond to different settings of λ indicated on the top. The solid linewith markers indicates the objective value versus the total CPU time of a single run of KELLER(i.e., to solve the N D logistic regression problems). Each line is plotted by varying termination
conditions in the logistic regression on the same dataset.

33

−2800

w=30

w=50

w=70

w=90

Test log−pseudolikelihood

−3000

−3200
−2800

−3000

−3200
−2800

−6
w=110

−3000

−3200

−6

−5
log10 λ−4−5log10 λ−4

No nuclear norm
log10 η = -2.25
log10 η = -2
log10 η = -1.75
log10 η = -1.5

Figure 14: US Senate data: Test log-pseudolikelihood was improved by nuclear-norm regularization. Each panel shows the results obtained using the proposed method with four η’s (the four
thick lines for which markers are only shown on the vertical axis, where η is the regularization
coefficient of the nuclear norm) and that by the method without the nuclear-norm regularization(the dashed line with the marker ‘∗’), for a specific kernel width w. In each panel, the verticaland horizontal axes denote the test log-pseudolikelihood and log10 λ (λ is the regularization
coefficient of the `1 norm), respectively. The markers on the vertical axis indicate the results
obtained with λ = 0, and the horizontal dashed line indicates the best result by the static MN,
which is the maximum test log-pseudolikelihood shown in Fig. 15.

34

Test log−pseudolikelihood

−2000
−4000
−6000
−8000
−4

−3

−2
−1
log10 λ0

1Figure 15: US Senate data: Test log-pseudolikelihood of static MNs in which the parameter
vector θ was common for all the time steps. The horizontal axis denotes log10 λ, where λ is theregularization coefficient for the `1 norm.

0

0

0

0

0

1

500 1

500 1

500 1

500 1

500 1

500 1

500 1

500 1

500 1

500

Figure 16: US Senate data: Time courses of 50 randomly-selected edge weights with (black line)
and without (gray line) the nuclear-norm regularization. Vertical and horizontal axes in eachpanel show the edge weights (rescaled in each panel) and the time (from 1 to 645), respectively.
The ordering of these panels (going downward from the top left) was determined so that similar
time courses (solid line) are close to each other.

35

Sununu(R−NH)

Proposed

Collins(R−ME)

No nuclear norm

Feingold(D−WI)
Lieberman(D−CT)
Wyden(D−OR)
Allard(R−CO)
DeMint(R−SC)
Ensign(R−NV)
Frist(R−TN)
Gregg(R−NH)
Hagel(R−NE)
Kyl(R−AZ)
McCain(R−AZ)
Roberts(R−KS)
Santorum(R−PA)
Sessions(R−AL)
Smith(R−OR)
Specter(R−PA)
Thune(R−SD)
Voinovich(R−OH)
0

Isakson(R−GA)

300

Time index

600 0

Proposed

300

Time index

6000Lincoln(D−AR)

No nuclear norm

Cochran(R−MS)
Alexander(R−TN)
Allard(R−CO)
Bennett(R−UT)
Bond(R−MO)
Bunning(R−KY)
Burns(R−MT)
Domenici(R−NM)
Grassley(R−IA)
Hutchison(R−TX)
Lott(R−MS)
Lugar(R−IN)
McConnell(R−KY)
Murkowski(R−AK)
Roberts(R−KS)
Shelby(R−AL)
Specter(R−PA)
Stevens(R−AK)
Voinovich(R−OH)
Warner(R−VA)
0

No nuclear norm

300

Time index

Proposed

600 0

300

Time index

600

No nuclear norm

Baucus(D−MT)
Cantwell(D−WA)
Carper(D−DE)
Feinstein(D−CA)
Inouye(D−HI)
Kohl(D−WI)
Landrieu(D−LA)
Nelson(D−FL)
Nelson(D−NE)
Pryor(D−AR)
Rockefeller(D−WV)
Salazar(D−CO)
Stabenow(D−MI)
Burr(R−NC)
Chafee(R−RI)
Grassley(R−IA)

Allard(R−CO)
Allen(R−VA)
Bond(R−MO)
Bunning(R−KY)
Burr(R−NC)
Chambliss(R−GA)
Coburn(R−OK)
Cornyn(R−TX)
DeMint(R−SC)
Ensign(R−NV)
Enzi(R−WY)
Inhofe(R−OK)
Martinez(R−FL)
Murkowski(R−AK)
Santorum(R−PA)
Vitter(R−LA)
Warner(R−VA)
0

Proposed

Harkin(D−IA)
Mikulski(D−MD)
Nelson(D−FL)
Salazar(D−CO)
Bond(R−MO)
Chafee(R−RI)
Coleman(R−MN)
DeWine(R−OH)
Dole(R−NC)
Graham(R−SC)
Lugar(R−IN)
McCain(R−AZ)
Murkowski(R−AK)
Smith(R−OR)
Snowe(R−ME)
Specter(R−PA)
Talent(R−MO)
Thomas(R−WY)

300

Time index

600 0

Proposed

300

Time index

6000Crapo(R−ID)

No nuclear norm

300

Time index

Proposed

600 0

300

Time index

600

No nuclear norm

Allen(R−VA)
Burns(R−MT)
Craig(R−ID)
Domenici(R−NM)
Enzi(R−WY)
Graham(R−SC)
Grassley(R−IA)
Hatch(R−UT)Martinez(R−FL)
Shelby(R−AL)
Thomas(R−WY)
Thune(R−SD)
300

Time index

600 0

300

Time index

600

0

300

Time index

600 0

300

Time index

600

Figure 17: US Senate data: Temporal changes in clustering of nodes (senators) over the 645
time points found by the affinity propagation clustering (APC) [17]. The APC partitions every
instance of the dependency graph into several clusters each represented by an influential nodecalled an exemplar. Each of the six blocks shows the changes of the assignment of cluster
members listed at the left side to a particular exemplar node indicated at the top left in abox. The annotation after the name of each senator refers to the senator’s political party(D: Democrat, R: Republican, I: Independent) and state (as a two-letter abbreviation). Thebη
“Proposed” and “No nuclear norm” indicate that the APC was applied to the estimate Θbby the proposed method and on Θ0 by the baseline method (η = 0), respectively, both at
w = 70; the regularization coefficients η and λ were chosen so as to roughly maximize the test
log-pseudolikelihood of Fig. 15 (see text). Each black horizontal bar indicates a time period
during which the cluster member was actually assigned to the exemplar by the APC.
36