Non-approximate Inference for Collective Graphical Models on Path Graphs via Discrete Difference of Convex Algorithm
Authors: Yasunori Akagi
Naoki Marumo
Hideaki Kim
Takeshi Kurashima
Hiroyuki Toda
Venue: NeurIPS 2021 Poster (2021)
Keywords: Collective Graphical Model
Aggregated Data
DC Algorithm
Minimum Cost Flow

Abstract: The importance of aggregated count data, which is calculated from the data of multiple individuals, continues to increase. Collective Graphical Model (CGM) is a probabilistic approach to the analysis of aggregated data. One of the most important operations in CGM is maximum a posteriori (MAP) inference of unobserved variables under given observations. Because the MAP inference problem for general CGMs has been shown to be NP-hard, an approach that solves an approximate problem has been proposed. However, this approach has two major drawbacks. First, the quality of the solution deteriorates when the values in the count tables are small, because the approximation becomes inaccurate. Second, since continuous relaxation is applied, the integrality constraints of the output are violated. To resolve these problems, this paper proposes a new method for MAP inference for CGMs on path graphs. Our method is based on the Difference of Convex Algorithm (DCA), which is a general methodology to minimize a function represented as the sum of a convex function and a concave function. In our algorithm, important subroutines in DCA can be efficiently calculated by minimum convex cost flow algorithms. Experiments show that the proposed method outputs higher quality solutions than the conventional approach. 

Non-approximate Inference for Collective Graphical Models on Path Graphs via Discrete Difference of Convex Algorithm
Yasunori Akagi
Naoki Marumo NTT Communication Science Laboratories,
NTT Human Informaties Laboratories,
NTT Corporation
NTT Corporation
yasunori.aken www.aaaaaaaaaaa.a. cu@hco. co.
naoki marumo. ec@hco ntt.to........ co.
Hideaki Kim
Takeshi Kurashima
NTT Human Informatics Laboratories,
NTT Human Informatics Laboratories,
NTT Corporation
NTT Corporation
hideaki kin, cn@hco nt..to.t...
takeshi .krrshhiaa uf @fc .tt. co. jp
Hiroyuki Toda
NTT Human Informatics Laboratories,
NTT Corporation
hiroyuki..m toda. xb@hco ntt..c..j
Abstract
The importance of aggregated count data, which is calculated from the data multiple individuals, continues to increase. Collective Graphical Model (CGM) is a probabilistic approach to the analysis of aggregated data. One of the most important operations in CGM is maximum posteriori (MAP) inference of unob- served variables under given observations. Because the MAP inference problem for general CGMs has been shown be NP-hard, approach that solves approximate problem has been proposed. However, this approach has two major drawbacks. First, the quality of the solution deteriorates when the values in the count tables are small, because the approximation becomes inaccurate. Second, since continuous relaxation is applied, the integrality constraints of the output are violated. To resolve these problems, this paper proposes a new method for MAP inference for CGMs on path graphs. Our method based on the Difference Convex Algorithm (DCA), which is a general methodology to minimize function represented as the sum of convex function and concave function. In our algo- rithm, important subroutines in DCA can be efficiently calculated by minimum convex cost flow algorithms. Experiments show that the proposed method outputs higher quality solutions than the conventional approach.
Introduction
In recent years, the importance of aggregated count data, which is calculated from the data of multiple individuals, has been increasing 21 27 Although technologies for acquiring individual data such as sensors and GPS have greatly advanced, it is still very difficult to handle individual data due to privacy concerns and the difficulty of tracking individuals. However, there are many situations where data aggregated from multiple individuals can be obtained and utilized easily. For example, Mobile Spatial Statistics 222] which is the hourly population data of fixed-size square grids calculated from
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

cell phone network data in Japan, are available for purchase; such data is being used for disaster prevention and urban planning [20]. In traffic networks, traffic volume data at each point can be obtained more easily by sensors or cameras than the trajectories of individual cars, and the data is useful for managing traffic congestion 99 (6266.
Collective Graphical Model (CGM) 17 is a probabilistic model to describe aggregated statistics of a sample drawn from a graphical model. CGM makes it possible to conduct various practical tasks on aggregated count data, such as estimating movements from population snapshots, parameter learning of the underlying graphical model, and interpolation and denoising of count tables. Particularly, the case where the underlying graph is a path graph is important because CGMs on path graphs can treat time series data which the states of interest follow Markov chains. In fact, most of the real- applications of CGMs utilize CGMs on path graphs to represent the collective movement of humans and animals 2 19. Detailed analyses of time series of collective people movements from limited observations would be useful for controlling people flow to avoid congestion and to maintain social distancing in urban spaces.
One of the most important operations in CGM is maximum a posteriori (MAP) inference. MAP inference is the discrete (combinatorial) optimization problem of finding an assignment of unobserved variables that maximizes the posterior probability under given observations. MAP inference makes it possible to interpolate missing values of aggregated data and to estimate more detailed information that lies behind the observations. Unfortunately, MAP inference for general CGMs has been shown to be NP-hard [16] and thus is difficult to solve exactly and efficiently. Therefore, an alternative approach that solves an approximate problem, which is derived by applying Stirling's approximation and continuous relaxation, has been proposed [6]. Subsequent studies have focused on solving this
approximate problem efficiently [12 18 19
23
However, there are inherent problems with this approach of solving the approximate problem. First, this approach tends to output a solution with a low posterior probability when the values in the count tables are small, because Stirling's approximation, log x!! log x - x, is inaccurate when is small. Such situation frequently occurs when the number of values that each variable in the graphical model takes is large, or when the sample size is small. Second, since continuous relaxation is applied, the integrality constraints of count table values are violated in the output. As a result, values that should be integers (e.g., the number of people) are no longer integers, which not only reduces interpretability, but also makes the output less sparse, resulting in high memory consumption to maintain the output. is possible to obtain integer-valued results by rounding the output, but this rounding process destroys the sum constraints among the estimated counts, e.g., the sum of the count table values at each node may not match the sample size.
To resolve these issues, in this paper, we propose a new method for MAP inference for CGMs on path graphs. We first show that the objective function of the problem can be expressed as the sum of univariate discrete convex functions and discrete concave functions. Based on this expression, we utilize the idea of the Difference of Convex Algorithm (DCA) 6 DCA is a framework to minimize a function expressed as the sum of convex function and concave function. In DCA, solution is obtained by repeatedly minimizing a surrogate function that upper-bounds the objective function, and the objective function value decreases monotonically in each iteration. In addition, the algorithm terminates in a finite number of iterations in our case since the variables are discrete not continuous.
The key to make the DCA-based algorithm efficient is a fast minimization algorithm for the surrogate function. Because the feasible region of our problem is limited to integer lattice points, continuous optimization methods such as the gradient descent, which are usually used in DCAs, cannot be applied to minimize our surrogate function. Instead, we utilize the special structure of path graphs; it enables us to formulate the minimization problem of the surrogate function as combinatorial optimization problem called the minimum convex cost flow problem. Fast algorithms for the minimum convex cost flow problem are known and we can minimize the surrogate function efficiently by using these algorithms.
The proposed method has several practical advantages. First, since the proposed method does not use Stirling's approximation, it offers accurate inference even when the values in the count tables are small. This makes it possible to output solutions with much higher posterior probability than the approximation based approach. Second, because the proposed method does not apply continuous relaxation, the obtained solution is guaranteed to be integer-valued, which results in sparse and interpretable outputs. In Section [ we show results gained from synthetic and real- wolldwordd datasets;
2

---

they indicate that the proposed method outputs higher quality solutions than the existing approach. We show that the superiority of the proposed method is much greater when the sample size is not very large the number of states on nodes in the graphical model is large.
Collective Graphical Models (CGMs)
Collective Graphical Model (CGM) is a probabilistic generative model that describes the distributions of aggregated statistics of a sample drawn from certain graphical model [17]. Let G = (V, E) be an undirected tree graph (i.e., a connected graph with no cycles). We consider pairwise graphical model over discrete random variable X = (X=)uEv defined by Pr(X x) Z [((,)) 9hn (Ius To), where Puo (Tu, Fy) is local potential function on edge (u, v) and Z 2 II ((( Puo (Tas Te) the partition function for normalization. In this paper, we assume that xi takes values the set [R] for all u E V. where [k] denotes the set (1,2..... } for positive integer k.
We draw an ordered sample (x(1), ,..., X (N)) independently from the graphical model, where M the sample size. Let nu (n (1))semm and nu == (t, j.) where n.((i) = Km x-(m) {}} and (i. = Km x(() i, x-(m) = Each entry of ny and nae is the number of occurrences of particular variable setting. We call (Na) node contingency table and (Niw)) (u.))E edge contingency table, and denote := ((n.) (eew) (U,))E) We assume that observations y == (yu) are generated by adding noise to the node contingency table (Nu)uev, and the distribution y is given by Pr(gy(n) I The (Vu (m,)) (i)), where Pui the noise distribution. additional assumption described below.
Assumption 1. For u E V and E [R), log (y)n)) a concave function n.
Assumption D is quite common assumption in CGM studies
[6 19
Commonly used noise
distributions such as Gaussian distribution Pui (v. (i))
V
exp
(v. and
20
Poisson distribution (Mu (z. (i)) = ne (1)...)))) /yu (i)! exp(-n.. satisfy Assumption
The MAP inference problem for CGM is to find that maximizes the posterior probability Pr(nly). Since Pr(nly) y) = Pr(n. y) Pr(y) from Bayes' rule, it suffices to maximize the joint probability
Pr(n, y) Pr(n) Pr(yy]n) Pr(n) is called CGM distribution and calculated as follows [I9]:
Pr(n)= = F(n) I(n E L2)),
(1)
Iuev ILERR (n (1))
M!
F(n) ZM
II
6u (()
(2)
[[ II.JE(R) (i,().
(u,))EE E.jEe
L = {n Z |
M = Σ (i) (u EN (
n, Σ (i. ((u,)) E F, i [R])] jE[R]
Here, II() is the indicator function, Vú is the degree of node u in G, and L is the set of possible contingency tables. Using the above notations, the MAP inference problem can be written
min NEL MM
log F(n) - log (r(y)n).
3 CGMs on Path Graphs
Hereafter, we focus on CGMs on path graphs, which is the main topic of this paper. Path graph P,, is an undirected graph whose vertex set V = [N] and edge set E = + | graphical model (not CGM) on path graph is the most basic graphical model that represents a time series generated by Markov model; that is, the current state depends only on the previous state. A CGM on path graph represents the distribution of aggregated statistics when there are many individuals whose state transition is determined by Markov model. In the rest of this paper, we use the notation Bu == (i))), Ray nl.t+ (i. and ouz 0r,t+ (i for simplicity, From I 4
[N- 1]}.
3

---

the MAP inference problem for CGMs on path graphs can be written as follows:
N
N
min
Σ
Sus(nes))
+
Σ
Σ
g(m))
Σ
Σ
he((nee)
2m2
(E(R)
FE(R)
s.t.
Σ
Tei
M (NE)
E
Daig
Ths
(LE - 11, 1,, iE [R]),
(5)
iE
JE[R]
Σ
Rej = Re-1d
(N - 11, jE [R]), Rujj Res Zzo.
iE|m
where (2) log z! 2 log Øesj. ((2) -a=--lop log 2!, (11(2) log Pts (yul2). For the details of the derivation, please see Appendix.
We give an example of a CGM on a path graph which models human mobility. Consider that space is divided into R distinct areas and that M people are moving around in the space. The random variable x-() represents the area which person m belongs at time step t, and the time series is determined by the graphical model p(x) 3 II 0tzeme+1 Here, t + 1. nh: represents the number of people in area i at time step i, and nęi represents the number of people who moved from area to at time step
x(m) (x(m)) ,...,NNN' (r))
Dtsy is the affinity between two areas and j at time step
+ 1. We have noisy observations Yua for t E [N] and E [R], which are generated by adding noise to Yas The MAP inference problem we want to solve is to find the true number of people of each area each time step, (Mr)) and the true number of people moving between each two areas, (N je) with the highest posterior probability given the observation (Mes)e N),
4 Proposed Method
4.1 Application of DCA
To solve problem (5). we propose utilizing the idea of the Difference of Convex Algorithm (DCA).
Before describing our method, we review the core idea of DCA based on the description in (I1].
DCA is a general framework to solve the minimization problem min P(n) = Q(n) + R(n). where ((n) is convex function and Rnn) is concave function. DCA does this by using the following procedure generate feasible solution sequence n² ..... (*) that satisfies P(n()) 2 P(n(2)) 2
((n( First, we choose arbitrary feasible solution (1) E D. When already have the sequence r1)), ..... n2 (5), find a function ((()()) that satisfies the following
three conditions: (i) ((() (n()) (+) R(n(()), (ii) (n) 2 R(n) (Vn E D), (ii) p(s) (n) = R())+R())(m) can be minimized efficiently D. Because R(n) is concave, by setting (n)
R(n!*)) (R(n(*)) (n ()))), which is linear approximation of R(n) at (*(), conditions (i)- hold. Using this function, we get new feasible solution by n2( (=-1) arg min, ((() This can be done easily because condition (iii) holds. Then, because ((n()+1)) < (()))(
p(+) (n (+)) (((n)))))())))))))) we get P(n(")) 2 ((n(2)) 2 ... ((n(+)) 2 (+n(++ by induction.
((n())),
To apply the framework of DCA, the objective function must be expressed the sum of convex and concave functions. The following proposition shows that our MAP inference problem in (5) has such structure.
Definition 1. A function f: Z>o + RU RU(+00) is called discrete convex function when f(z + 2) + f(z) 2 f(z for all Zzo -f a discrete convex function, is called a discrete concave function.
Propusition 1. fig is a discrete convex function. Under Assumption 7 his a discrete convex function. is a discrete concave function.
The proof is given in the Appendix. Hereafter, we set Q(n) = EN Σ Eije[R fos (mis)) + MN Exe[R] he (nx) and R(n) Σm EiE(R) g(re)). Thanks to Proposition we can ap- ply DCA to our problem. The following proposition provides a function R(() (m) that satisfies conditions (i) and (ii) required for DCA.

---

U1
W
U2
W
U2
W2
21.1
W1
1211
W2.1
23.11
W3,1
b.
M
d
by
M1.2
W22
22.2
W2.2
M3.2
W3,2
Cost function
(2)
(2)
2. + (0)
(2)
(z
Figure An example the MCFP instance defined in Proposition when N and R =
Proposition 2. Let Dri (z) = log(n/; 1) at (z-- where dis is real number which sat- isfies - log(n (x) 1) S ag
(x) S - log i Then, the function R())(n) = => 2 Evenn Did (na)
N-11
satisfies (n()) = ((n( and (n)) 2 R(n).
Please the Appendix for the proof. Intuitively, Ies is a tangent of g at Alri.
4.2 Minimum Cost Flow Algorithm for the Subroutine
The most important and difficult part to derive an efficient OCA-based algorithm is designing efficient algorithms for the problem min p(s) (n) (condition (iii). To achieve this, we show that minneD (D=)(n) can be formulated as the Minimum Convex Cost Flow Problem (C-MCFP), which the efficiently solvable subelass of the Minimum Cost Flow Problem (MCFP). The (non-linear) MCFP is a combinatorial optimization problem on directed graph G = (V, Each node E V has (+oo), MCFP is the problem of finding a minimum cost flow on G that satisfies the supply constraints at all nodes. MCFP can be described follows:
supply value 6; E Z. and each edge (i,j) e & has cost function Cjj Zzo
min #E 20
Σ cog(21g)
S.t.
Σ ZSj
Σ J-(J))EE
23 = b, E
(KJ)EE
j((jj)E
Note that z takes only integer values (i.e., E Z A subclass of MCFP in which all cost
functions are discrete convex functions (see Definition efficiently solvable M The following proposition shows that the subproblem minnen (())()) can be formulated as C-MCFP.
I is called the C-MCFP; is known to be
Proposition Define the MCFP instance as follows:
the node V is defined by V = {o, d3 U (ULEIN)(L WW)),, where Us (ue,)sE(R)) W, = (E)
the edge set E consists of four rypes of edges,
edges (o, u1,6, 0) and (wN, d, for [R].
edges (u1,3, Wt,i hes (z)) r = 1, N and i [R).
edges (ur,s, We,s Hia (z) + hex(2)) for = 2,...,
and E [R].
edges (We,ss U441.6 f(is(z) (2)))foorrr) for E [N - 1] and j E [R].
where (u, v, c(2)) represents a directed edge from node node with cost function c(2). the supply values (bi)tev are defined by b, M, ba M, and b 0 V {o, d}. Let z" is an optimal solution of this MCFP instance, and define n by n == 20,, and niis
to Wer..." Then, n² is an optimal solution of the problem minneD (()).)). Furthermore, the MCFP instance belongs C-MCFP.
The proof is given in the Appendix. Figure illlssraates example of the MCFP instance defined in Proposition 3] This MCFP can be interpreted as the problem of finding way to push M flows from node node d with minimum cost. The above proposition enables to solve the subproblem minnED p(a) (n) efficiently by applying existing algorithms for C-MCFP.

---

4.3 Overall View of the Proposed Method and Time Complexity Analysis
From the above arguments, we can construct an efficient optimization algorithm, described in Al- gorithm [ for the MAP inference problem (5)
Algorithm DCA for problem (
The algorithm is guaranteed to terminate after a finite number of iterations because ((n(*)) mono-
1: n t- 0
tonically decreases and D is finite set.
2: for = 1.2, ....... do
...
3:
n (++1) <- (n" defined in Proposition 3
We analyze the time complexity of one iteration of the proposed method (Lines 3-5 in Algorithm T The computational bottleneck is solving C-MCFP
4:
P(n(")) = (n((+)) then
5:
return n (s)
in Line 3. There are several algorithms to solve C-MCFP and time complexity varies depending on which one is adopted. In this paper, we consider two typical methods, the Successive Shortest Path algorithm (SSP) and the Capacity Scaling algorithm (CS) (..
SSP is an algorithm that successively augments unit flow along the shortest path from a supply node (i.e. 6; > 0) to demand node (i.e. calculated from the current flow. Given a C-MCFP instance with graph G = (V,E), the shortest path in the residual graph can be computed in O(E log [I) time by Dijkstra's ' algorithm with binary heap, and the augmentation of the flow can be done in O(E)) time. The augmentation is performed the total time complexity is O log [V])- CS resembles
0) in the residual graph, which is auxiliary graph
B = (Eie V 122 times totally, SSP, but it differs in that it tries to push a large amount of flow, rather than a unit amount of flow, in a single augmentation. In CS, the number of shortest path calculations and flow augmentations can be bounded O(E log U') times, where U == maxiey [b:]. so the total computational complexity is O(EE log |V log U).
Because | = O R). [] =()=)())). B (M) and = O(M) in our problem, the time com- plexity in one iteration is O(MNR² log(NR)) when SSP is applied and O Æ log(NN R) log M) when CS is applied. This result implies that each method has its own advantages and disadvantages: SSP has small time complexity for N and R, while CS has small time complexity for M. This difference is confirmed empirically in Section 5.1
5 Experiments
We perform experiments to evaluate the effectiveness of the proposed method using synthetic and real-world instances. All experiments are conducted on 64-bit macOS machine with Intel Core i7 CPUs and 16 GB of RAM. All algorithms are implemented in C++ (gee 9.1.0 with -O3 option).
5.1 Synthetic Instances
Settings. We solve randomly generated synthetic instances the MAP inference problem (5) and compare the attained objective values. We fix N to 5 and vary the values of R and M. The input observation yıs is independently drawn from uniform distribution on the set of integers {1, 2, ...,....22.
LA []. As the noise distribution, we use Gaussian distribution Pta Pua(Vs/na x exp (- 0.01 (yus- ). We use two types of potential functions as follows. (1) uniform. Pzgj is independently drawn from uniform distribution on the set of integers (1, (1,2,-. ........ 10). (2) distance. We set daij T-
This potential models the movement of individuals in one-dimensional space: the state indices and j represent coordinates in the space, and the closer the two points are, the more likely are movements between them occur.
Proposed Method.
To construct surrogate functions in the proposed method, we can choose
arbitrary at which satisfies the condition - log(n (x) 1) S a < log n (see Proposition
(s)
(s)
2). To investigate the influence of the choice of a we try three strategies to decide o (s), (1)
(s)
P log(n(e) (+) (2) (a) = -1 (log(n(.) (n(( log(n(4) + 1)), (3) (*) - log(n (a)
1). We call
them Proposed (L), Proposed (M), Proposed (R), respectively.

---

Table 1: Attained objective functions in synthetic instances. For each setting, generated 10 instances and average values are shown. The smallest value is highlighted for each setting. and D mean the "uniform" and "distance" potential settings, respectively.
M
10"
10"
10
20
30
10
20
30
Proposed (L)
-9,97e+01
-8,90++01
-8,74c+01
-I11e+003
-1.19++03
-1.22e+03
-1.07+004
-1.3e+04
-1.40er44
U
-9.81c+01
-8.90c+01
-8.74c+01
-1.11c005
-1.19c+08
-1.22c+03
-1.07c404
-1.31e404
-1.40c+04
Proposed (R)
-9.44c+01
-8.76e401
-8.74e+01
-1.1
-1.18c103
-1.21e103
-1.07e404
-1.31e004
-1.40e404
NLBPP
7.19c+01
-7.0e++1
-7.0e==01
1.08e+03
9.77+02
9.02e402
1.07e404
1.300+04
-1.37e004
Proposed (.)
3.35c-01
5.00c-01
5.00c41
-548ee01
-303er01
-1.18e01
3.83400
-.06ee02
-LO1er03
D Proposed (M)
3.35e-01
5.00c.01
5.00c.01
S.43++01
291e+01
-1.44401
5.82+400
9.06++02
-1.0e+03
Proposed (R)
3.35e-01
5.00e-01
5.00c401
-5.39e001
2.89x001
-1.06400
-58804000
-9.0e+022
-1.0le033
NLBP
3.20+401
4.56e+01
5. 5.28c+01
-1.38e+01
1.77e402
3.25e+02
1.20e+00
8.02e+02
5.31e+02
Compared Method. As the compared method, we use Non-Linear Belief Propagation (NLBP) 19 which is message-passing style algorithm to the solve approximate MAP inference problem derived by applying Stirling's approximation and continuous relaxation. Because the output of NLBP is not integer-valued and log(2!) is defined only if 2 is an integer, we cannot calculate the objective function of (5) directly. To address this, we calculate it by replacing the term log(z!) by linear interpolation of
[=]] and log([=]] which given by [[х] - 2) log[[=]!)) + (z [=]) log([ > !!). Note that although there various algorithms to solve the approximate MAP inference problem (see Section
6.I). the objective function values attained by these algorithms are the same. This is because the
approximate problem is convex optimization problem [16].
Comparison of attained objective values. The results shown in Table [ We generated 10 instances for each parameter setting and determined the average of attained objective function values. Because the objectíve function P(n) is equal to log Pr(n/y) + const.., P(n) takes both positive and negative values, and the difference of the objective function values is essential; when ((m))-P(nn) = 8, Pr(ni y) = exp(-8) Pr(nely) holds.
All the proposed methods consistently have smaller objective function values than the compared method. The difference tends to be large when R is large and M is small. This would be because small values appear in the contingency table more frequently when R is large and M is small, and the effect of the inaccuracy of Stirling's 's approximation becomes larger. Among the three proposed methods, there was not much difference in obtained objective function values although Proposed (L) was found to consistently achieve slightly smaller objective function values than others. This indicates that the proposed method is robust with respect to the choice of hyperparameters Oe
Characteristics of the output solution. To compare the characteristics of solutions obtained by proposed (L) and NLBP, we solve an instance with R 20, M 102. and uniform potential by each method. Obtained edge contingency tables 14 are shown in Figure 2 as heat maps. We also show the edge contingency table obtained by rounding the NLBP solution to integers. We observe that the proposed method outputs sparse solutions while the solutions by NLBP are blurred and contain lot of non-zero elements. This difference quantified by "sparsity", which is calculated by 1.0- (# of non-zero (> 10 10-2) clements)/(# of elements): sparsity of the output of proposed (L) 77%, while the sparsity of the output of NLBP is 0%. This is caused by its application of continuous relaxation and the inaccuracy of Stirling's approximation around O. In the solution of NLBP (rounded), many near-zero values are rounded to O, and constraints of the problem (5) are totally violated; for example, the sum of the edge contingency table values does not match the sample size, additional experiments, we observed that the outputs of the three methods become closer as M increases. For more details, please see the Appendix.
Comparison of computation time. We compare the computation time of each algorithm. As explained in Section 4.3) we can choose an arbitrary C-MCFP algorithm as the subroutine in the proposed method and the time complexity varies depending the choice. We compare proposed (L)
with SSP, proposed (L) with CS, and NLBP.
Figure 3 shows the relationship between input size and computation time, and Figure 4 shows the relationship between running time and objective function value. These results are consistent with the complexity analysis results in Section 4.3 SSP is efficient when R is large but becomes inefficient
7

---

Proposed (L), Sparsity: 77%
NLBP, Sparsity: 0%
NLBP (rounded), Sparsity: 99%
Figure 2: Comparison of solutions yielded by proposed method (L), NLBP, and NLBP (rounded), We solve an instance with R = 20, M 102 and uniform potential. The obtained edge contingency table ma is presented matrix heatmap with the maximum value of color map 3.
Table Attained objective function values and NAEs for real-world instances, For each setting, we generated 10 instances and averages are shown. The smallest value is highlighted.
56
208
56
208
208
Proposed (L)
-3002e+12
2,97++02
-6.61e+43
-2.62e-003
-1.70e+04
-9,76eH3
Obj.
1.30c+003
4.70c=03
NLBP trounded)
Proposed (1.)
0.690
0.441
1.002
0.829
L,073
0.956
NAE
NLBP NLBP (rounded)
1.38
1.544
1.241
1.438
1,209
1.381
0,877
0,870
1.079
0,907
1.128
0.991
when M is large, and the converse is true for CS. The results also suggest that it is important to choose the algorithm depending on the size parameter the input. The proposed method is not much worse than the existing method in terms of computation time by choosing appropriate C-MCFP algorithm according to the size the input. Appropriately chosen proposed methods attain the minimum of the existing method more quickly; proposed methods take lot of time to achieve smaller objective function value than the minimum of the existing method.
5.2 Real-world Instances
We conduct experiments using real-world population datasets. The datasets are generated from 8694 car trajectories collected by navigation application in the Greater Tokyo area, Japan We randomly sample M (M = 100, 500, 1000) trajectories from this data and create aggregated population data of each geospatial space into fixed-size grid cells. The grid size is set to 10km 10km (R = x 7 = 56) and 5km 5km (R 16 13 208), and time interval is 60 minutes (N 24). As the noise distribution, we use Gaussian distribution Peo (yt (as)) X exp (-( ) We construct the potential duj = exp (- dist(i, j)), where dist(i, j) is the Euclidean distance between the centers of cell and cell in the grid space. We create 10 instances by random sampling and averaged the attained objective function values for each setting. We also evaluate the estimation accuracy of We normalized absolute error (NAE) as the
fixed time intervals. The areas decided by dividing the targeted
the edge contingency table (NN
evaluation metric, which is defined ΣIE(N Eigrry Σje(r) Pe55 -1445 where n is the true
EieNN-1 Ene(R LFEER teue
value and miy is the estimated value of the edge contingency table. In addition to Proposed (L) and NLBP, we also evaluate the estimation accuracy of NLBP (rounded), which a method that rounds the output of NLBP to integer values. Note that do evaluate the objective function values of NLBP (rounded). This because the output of NLBP (rounded) completely violates the constraints summation in the MAP inference problem (5) and the objective function value cannot evaluate whether the optimization problem is successfully solved or not.
We use the data collected by the smartphone car navigation application of NAVITIME JAPAN Co., Ltd. (http://corporate navitime.iocoooo wwatttimm.c...jp.ee... jp/en//) The data collected with consent and appropriately anonymized.

---

NLBP
102
100
CS
SSP
time (sec)
10
10-1
100
NLBP
10
CS
10-1
10'
10a
10a
10
10 20 30
40 50 60 70 80 90100
M
R
Figure The computation time of each algorithm. The values are averages of synthetic instances when R is fixed to 20 (left) an M is fixed to 10 (right). N is set and the uniform potential is used.
-17800
13400
NLBP
NLBP
SSP
SSP
13600
tg0 Val.
-117850
13800
-117900
14000
0.0
0.5
1.0
1.5
2.0
0.0
0.5
1.0
1.5
time (sec)
time (sec)
Figure 4: The relationship between running time and objective function value. The left figure shows the result of instance of R 20, .M = 10a and uniform potential, and right figure shows that of R = 30, M 10 and uniform potential.
Table Sssows the results. We observe that Proposed (L) consistently attain smaller objective values and NAEs than the existing method. The superiority of the proposed method increase when R is large and M is small, and this is the same trend as the results with synthetic data. The NAE of NLBP relatively large; this is due the fact that small values that should be 0, which is the same phenomenon seen in Figure 2] The NAE values are improved to some extent by rounding, but the proposed method still superior.
assigned the elements of the output
Related Work
6.1
MAP inference for CGMs
Several methods have been proposed the MAP inference of CGMs, but most of them take the approach of solving the approximate problem [16. which is derived by applying Stirling's
approximation and continuous relaxation. For example, the interior point method [16], projected gradient descent (23). message passing [19 and Bankhorm-Knopp algorithm (18) have been used
to solve the approximate problem. In particular, [12] proposes a method to use DCA to solve this approximate problem. Although this approach is similar to proposal in that it uses DCA, the purpose of applying DCA is totally different: our focus is to solve the MAP inference problem without using any approximation continuous relaxation.
One of the few exceptions the method proposed 3, which solves the original MAP inference problem directly without using approximation. Our method follows this line of research, but there are two major differences. First, their method can only be applied to CGM on graph with two vertices, and thus applicability is very limited. Since our method is consistent with this method when applied

---

to CGM on graph with two vertices, our method can be regarded as generalization of their method. Second, their work assumes accurate observations and does not handle observation noise.
13 solves related collective MAP inference problems on path graphs. The problems addressed in this paper are different from ours; their purpose is finding the most likely assignments of the entire variables for each individual, while our purpose is finding the most likely node and edge contingency tables. In their settings, non-linear terms in the log posterior probability vanish, and the MAP inference problem can be solved easily by linear optimization approaches.
6.2 Difference of Convex Algorithm (DCA)
DCA, which is sometimes called Convex Concave Procedure [25]. is framework to minimize function expressed as the sum of a convex function and a concave function [6]. DCA was originally proposed as method for optimization in continuous domains. DCA has been used in various machine learning fields, such as feature selection
[7. reinforcement learning [4]. support vector machines
24 and Boltzmann machines
[3]
Several studies have applied DCA to discrete optimization problems. This line research sometimes called discrete DCA [8. 11 propose algorithms to minimize the sum of a submodular function and supermodular function. This algorithm generalized to yield the minimization of the sum of an M/L.-convex function and an M/L--onnaae function [8] where M-convex function
and L-convex function are classes of discrete convex functions [[O]. Although our work is closely related these studies, it is not part of them. This is because problem can be regarded as the minimization the sum two M-convex functions and separable concave function, and this is not included in the class of functions dealt with in g
7 Conclusion
In this paper, we propose a non-approximate method to solve the MAP inference problem for CGMs on path graphs. Our algorithm is based application of DCA. In the algorithm, surrogate functions can be constructed closed-form and minimized efficiently by C-MCFP algorithms. Experimental results show that our algorithm outperforms approximation-base based methods in terms of quality of solutions.
References
[1] Ravindra K. Ahuja, Thomas Magnanti, and James Orlin. Network Flows: Theory, Algorithms, and Applications. Prentice-Hall, Inc., 1993.
[2] Yasunori Akagi, Takuya Nishimura, Takeshi Kurashima, and Hiroyuki Toda. A fast and accurate method for estimating people flow from spatiotemporal population data. IJCAI, pages 3293-3300, 2018.
[3] Yasunori Akagi, Takuya Nishimura, Yusuke Tanaka, Takeshi Kurashima, and Hiroyuki Toda. Exact and efficient inference for collective flow diffusion model via minimum convex cost flow algorithm. In AAAI, pages 3163-3170, 2020.
[4] Jiali Du, Akshat Kumar, and Pradeep Varakantham. On understanding diffusion dynamics of patrons at theme park. In AAMAS, pages 1501-1502, 2014.
[5] Rishabh Iyer and Jeff Bilmes. Algorithms for approximate minimization of the difference between submodular functions, with applications. In UAI, pages 407-417, 2012.
[6] Hoai An Thi and Tao Pham Dinh. DC programming and DCA: Thirty years of developments. Mathematical Programming, 169(1):5-68, 2018
[7] Hoai An Le Thi, Hoai Minh Le, and Tao Pham Dinh. Feature selection in machine learning: an exact penalty approach using difference convex function algorithm. Machine Learning, 101(1):163-186 2015.
AA separable convex function is both L-convex and M-convex, but sum of two M-convex functions is neither M-convex nor L-convex.
10

---

[8] Takanori Maehara and Kazuo Murota. A framework of discrete DC programming by discrete convex analysis. Mathematical Programmiing. 152(1-2):435 2015.
[9] Tetsuro Morimura, Takayuki Osogami, and Tsuyoshi Idé. Solving inverse problem of markov chain with partial observations. In NIPS, pages 1655-1663, 2013.
[10] Kazuo Murota. Discrete convex analysis. Mathematical Programming, 83(1):313-371, 1998.
[11] Mukund Narasimhan and Jeff Bilmes. submoodula--superrmduu procedure with applications
to discriminative structure learning. In UAI, page 404-412, 2005.
[12] Thien Nguyen, Akshat Kumar, Hoong Chuin Lau, and Daniel Sheldon. Approximate inference using DC programming for collective graphical models. In AISTATS, pages 685-693, 2016.
[13] Atsushi Nitanda and Taiji Suzuki. Stochastic difference of convex algorithm and its application to training deep Boltzmann machines. In AISTATS, pages 470-478. PMLR, 2017.
[14] Bilal Piot, Matthicu Geist, and Olivier Pietquin. Difference of convex functions programming for reinforcement learning. In NIPS, pages 2519-2527, 2014.
[15] Daniel Sheldon, MA Saleh Elmohamed, and Dexter Kozen. Collective inference on markov models for modeling bird migration. In NIPS, pages 1321-1328, 2007.
[16] Daniel Sheldon, Tao Sun, Akshat Kumar, and Tom Dietterich. Approximate inference in collective graphical models. In ICML, pages 1004-1012, 2013.
[17] Daniel R. Sheldon and Thomas G. Dietterich. Collective graphical models. In NIPS, pages 1161-1169, 2011.
[18] Rahul Singh, Isabel Haasler, Qinsheng Zhang, Johan Karlsson, and Yongxin Chen. Inference with aggregate data: An optimal transport approach. 2020. 21V1:2003.119333
[19] Tao Sun, Daniel Sheldon, and Akshat Kumar. Message passing for collective graphical models.
In ICML, pages 853-861, 2015.
[20] Toshihiro Suzuki, Masashi Yamashita, and Masayuki Terada. Using mobile spatial statistics in field of disaster prevention planning. NTT DOCÓMO Tech. J, 14(3):37-4,, 2013.
[21] Yusuke Tanaka, Toshiyuki Tanaka, Tomoharu Iwata, Takeshi Kurashima, Maya Okawa, Yasunori Akagi, and Hiroyuki Toda. Spatially aggregated gaussian processes with multivariate areal outputs. In NeurlPS, pages 3000-3031, 2019.
[22] Masayuki Terada, Tomohiro Nagata, and Motonari Kobayashi. Population estimation technology for mobile spatial statistics. NTT DOCOMO Technical Journal, 14(3):10-15, 2013.
[23] Luke Vilnis, David Belanger, Daniel Sheldon, and Andrew McCallum. Bethe projections for non-local inference. In UAI, pages 892-901, 2015.
[24] Hai-Ming Xu, Hui Xue, Hong-Hong Chen, and Yun-You Yun Wang. Solving indefinite kernel support vector machine with difference of convex functions programming. In AAAI, pages 2782-2788. 2017.
[25] Alan L Yuille and Anand Rangarajan. The concave-convex procedure (CCCP). In NIPS, pages 1040-1040, 2001.
[26] Shanghang Zhang, Guanhang Wu, Joao P Costeira, and Jose MF Moura. Understanding traffic density from large-scale web camera data. In CVPR, pages 5898-5907, 2017.
[27] Yivan Zhang, Nontawat Charoenphakdee, Zhenguo Wu, and Masashi Sugiyama. Learning from aggregate observations. In NeurlPS, pages 470-478, 2020.

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] mentioned that the proposed
method may be inferior to the existing methods terms of computation time depending
on the parameters settings of inputs in Section [
(c) Did you discuss any potential negative societal impacts of your work? [N/A] Our study unlikely to have negative impact on society.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] All assumptions required for the theoretical results are clarified.
(b) Did you include complete proofs of all theoretical results? [Yes] See Appendix.
3. If you ran experiments..
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or as a URL)? [No] Because they the properties of the company to which I belong, I cannot disclose the details.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] See Section
(c) Did you report error bars (e.g., with respect to the random seed after running experi- ments multiple times)? [Yes] See Section
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, cloud provider)? [Yes] See Section 5
4. If you using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] See the footnote in
Section 5
(b) Did you mention the license of the assets? [No] We have purchased the data and have a license use it.
(c) Did you include any new assets either in the supplemental material or URL? [No] The data is the company's 's property and cannot be disclosed.
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [Yes] See the footnote in Section 5
(c) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [Yes] See the footnote in Section
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A] We do not use crowdsourcing and conducted research with human subjects. The answers to the following questions are the same.
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent participant compensation? [N/A]
12


Reviews:
[{"review_id": "wQQO6e_Drf", "reviewer": "Reviewer_6fXy", "summary": "The paper proposes a new method of estimating MAP solutions from a collective graphical model. ", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The proposed method is very solid. Solving the MAP problem in the discrete space does not need approximating the log joint function as in continuous approximations.\n\nThe performance of the proposed method outperforms the baseline. \n\nThe writing is generally clear. \n\nThe work also has a few limitations. \n1. The proposed method only applies to a small scope of problems. With the few characterizations (CGM, path-graph, and small counts), the proposed method can only show strength in a small number of problems. \n\n2. The experiment on the real data uses true counts as the ground-truth, but this may not be correct. With the model specification in the experiment, true counts may not be the MAP solution. \n\nSmall issues: \nIs it really useful to differentiate the three algorithms L, M, and R given their similarity? If so, can the authors put more analysis of how differences arise from the three configurations? \n\nI am not sure the two propositions make the writing clear. They are only analyses of the problem in focus. The notations are from the previous text, so the propositions do not stand by themselves. \n\nA little more insight into C-MCFP can assist the understanding. I have a rough idea that one can get a good solution with (even greedy) descending steps. But can you indicate the complexity here? \n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "No potential negative societal impact is detected. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "tpo7NhKgq5O", "reviewer": "Reviewer_2Z5W", "summary": "The paper proposes a non-approximate method for solving MAP inference problem in collective graphical model (CGM). The paper uses difference of convex algorithm (DCA) strategy as a solver for the optimization problem since the objective function can be expressed as a sum of convex and concave function. The solution is obtained via minimization of a surrogate function that upper bounds the objective function. The technical details appear to be correct and the simulation results show the superiority of the proposed algorithm over its competitors. Further, a real world application involving car trajectories is performed to validate the performance of the proposed methodology.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The contribution appears to be technically correct and is a valid addition to the graphical model literature. Although the DCA has been used in other contexts its application to solve the MAP inference in CGM is novel. The algorithm also not too much worse compared to its competitors in terms of scalability. The paper is overall well-written although at some places (especially in the section where minimum cost flow algorithm was described) a bit more clarity would have been ideal. I have a couple of minor comments for the author(s)\n\n- How does the choice of the surrogate function in DCA affect the MAP inference?\n\n-How does one choose the value of M in the c-MCFP algorithm?\n\n-Is there any practical guideline when to use the successive shortest path (SSP) or the capacity scaling (CS)?\n\n\n******* Change after rebuttal **********\n\nI am satisfied with the response from the author(s). The rebuttal has answered most of the queries related to the technical details of the paper. It also has also clarified some of the questions regarding the experiments with some additional results. I am keeping my previous score \"Good paper, accept\" unchanged.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Nothing as such.", "review_text": "", "extra_scores": {"ethical_concerns": "None.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3 hours"}, "extra_text": {}}, {"review_id": "Lre_pLXBG_", "reviewer": "Reviewer_fJvy", "summary": "The authors consider the problem of MAP estimate of count data given noisy count data. For example $X_1(t), X_2(t), \\ldots, X_M(t)$ describe the location of an agent at time $t$, $N_i(t)$ is the number of agents at location $i$ at time $t$, and $Y_i(t)$ is a noisy estimate of $N_i(t)$. The goal is to estimate $N$ given $Y$. The noise distribution is assumed to be log-convex.\n\nThe authors develop an approach based on the difference of convex algorithm, which is an extension of the convex-concave procedure. They show that this procedure can monotonically improve a solution by repeatedly solving the minimum convex flow problem, which is known to be efficiently solvable.\n\nThe authors demonstrate that their approach improves over non-linear belief propagation with a continuous approximation of the discrete time problem.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Pros:\n\nThe paper uses several interesting techniques to solve the problem of MAP estimate of count data on path graphs. I found the application of DCA and the use of min-cost convex flow to be interesting. They are able to work with the discrete problem, rather than previous work which needed to do a continuous relaxation.\n\nThe paper is mostly clear, however the notation is somewhat dense. This can be improved by adding names to the key variables. For example, instead of just $M$ the authors can write \"the number of samples $M$\", $\\mathcal{Q}$ can be \"the convex term $\\mathcal{Q}$\", $\\bar{\\mathcal{R}}$ can be \"the proxy concave term $\\bar{\\mathcal{R}}$\", etc. This will give the reader a reminder of what the variable represents.\n\nCons:\n\nThe algorithm presented is a heuristic without much analysis. It can be seen as somewhat incremental over [12], which also uses DCA.\n\nThe authors could compare to prior work on the bird migration benchmark, as used in [12]. Additionally, the experiments considered are all small scale and can be solved in less than a second. Can larger experiments be completed? This should be possible, at least for the synthetic experiments.\n\nQuestions:\n\nWhere would the noise on the count data come from in practice? To me, it seems more likely that there is noise on the positions $X_j(t)$. The authors should ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "This work can be used for population tracking. However I am not too concerned about this paper increasing the negative aspects of this, as there are already methods that can solve the problem considered.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors present a new algorithm for MAP inference in chain-structured Collective Graphical Models that uses flow techniques to solve the problem in discrete space, as opposed to previous works which use continuous relaxations. Three reviewers found the techniques to be correct and represent an advance over prior methods, with sufficient evidence presented to show the advantages of working in discrete space. The scores were 6, 6, 7. One reviewer raised their score after the author response. The meta-reviewer recommends accept. 
Decision: Accept (Poster)