* De Paolis et al. (2020) De Paolis, F., Nucita, A. A., Strafella, F., Licchelli, D., & Ingrosso, G. 2020, Mon. Not. Roy. Astron. Soc., 499, L87, doi: 10.1093/mnrasl/slaa140
  * Del Pozzo (2012) Del Pozzo, W. 2012, Phys. Rev. D, 86, 043011, doi: 10.1103/PhysRevD.86.043011
  * Diaz & Mukherjee (2021) Diaz, C. C., & Mukherjee, S. 2021.  https://arxiv.org/abs/2107.12787
  * Diego (2019) Diego, J. M. 2019, Astron. Astrophys., 625, A84, doi: 10.1051/0004-6361/201833670
  * Ding et al. (2019) Ding, X., Biesiada, M., Zheng, X., et al. 2019, JCAP, 04, 033, doi: 10.1088/1475-7516/2019/04/033
  * Drlica-Wagner et al. (2018) Drlica-Wagner, A., Sevilla-Noarbe, I., Rykoff, E. S., et al. 2018, ApJS, 235, 33, doi: 10.3847/1538-4365/aab4f5
  * Eisenstein & Hu (1997) Eisenstein, D. J., & Hu, W. 1997, Astrophys. J., 511, 5, doi: 10.1086/306640
  * Eisenstein & Hu (1998) —. 1998, Astrophys. J., 496, 605, doi: 10.1086/305424
  * Ezquiaga & Holz (2021) Ezquiaga, J. M., & Holz, D. E. 2021, Astrophys. J. Lett., 909, L23, doi: 10.3847/2041-8213/abe638
  * Farmer et al. (2019) Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, ApJ, 887, 53, doi: 10.3847/1538-4357/ab518b
  * Farr et al. (2019) Farr, W. M., Fishbach, M., Ye, J., & Holz, D. 2019, Astrophys. J. Lett., 883, L42, doi: 10.3847/2041-8213/ab4284
  * Feeney et al. (2019) Feeney, S. M., Peiris, H. V., Williamson, A. R., et al. 2019, Phys. Rev. Lett., 122, 061105, doi: 10.1103/PhysRevLett.122.061105
  * Finke et al. (2021) Finke, A., Foffa, S., Iacovelli, F., Maggiore, M., & Mancarella, M. 2021, JCAP, 08, 026, doi: 10.1088/1475-7516/2021/08/026
  * Fishbach & Holz (2017) Fishbach, M., & Holz, D. E. 2017, ApJ, 851, L25, doi: 10.3847/2041-8213/aa9bf6
  * Fishbach et al. (2018) Fishbach, M., Holz, D. E., & Farr, W. M. 2018, Astrophys. J. Lett., 863, L41, doi: 10.3847/2041-8213/aad800
  * Fishbach & Kalogera (2021) Fishbach, M., & Kalogera, V. 2021, ApJ, 914, L30, doi: 10.3847/2041-8213/ac05c4
  * Fishbach et al. (2019) Fishbach, M., Gray, R., Magaña Hernandez, I., et al. 2019, ApJ, 871, L13, doi: 10.3847/2041-8213/aaf96e
  * Fishbach et al. (2021) Fishbach, M., Doctor, Z., Callister, T., et al. 2021, ApJ, 912, 98, doi: 10.3847/1538-4357/abee11
  * Freedman (2017) Freedman, W. L. 2017, Nature Astron., 1, 0121, doi: 10.1038/s41550-017-0121
  * Fryer et al. (2001) Fryer, C. L., Woosley, S. E., & Heger, A. 2001, Astrophys. J., 550, 372, doi: 10.1086/319719
  * Gallegos-Garcia et al. (2021) Gallegos-Garcia, M., Berry, C. P. L., Marchant, P., & Kalogera, V. 2021, arXiv e-prints, arXiv:2107.05702.  https://arxiv.org/abs/2107.05702
  * Górski et al. (2005) Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759, doi: 10.1086/427976
  * Graham et al. (2020) Graham, M. J., Ford, K. E. S., McKernan, B., et al. 2020, Phys. Rev. Lett., 124, 251102, doi: 10.1103/PhysRevLett.124.251102
  * Gray et al. (2021) Gray, R., Messenger, C., & Veitch, J. 2021, arXiv e-prints, arXiv:2111.04629.  https://arxiv.org/abs/2111.04629
  * Gray et al. (2020) Gray, R., et al. 2020, Phys. Rev. D, 101, 122001, doi: 10.1103/PhysRevD.101.122001
  * Hannuksela et al. (2019) Hannuksela, O. A., Haris, K., Ng, K. K. Y., et al. 2019, ApJ, 874, L2, doi: 10.3847/2041-8213/ab0c0f
  * Heger & Woosley (2002) Heger, A., & Woosley, S. E. 2002, Astrophys. J., 567, 532, doi: 10.1086/338487
  * Hirata et al. (2010) Hirata, C. M., Holz, D. E., & Cutler, C. 2010, Phys. Rev. D, 81, 124046, doi: 10.1103/PhysRevD.81.124046
  * Holz & Hughes (2005) Holz, D. E., & Hughes, S. A. 2005, Astrophys. J., 629, 15, doi: 10.1086/431341
  * Holz & Wald (1998) Holz, D. E., & Wald, R. M. 1998, Phys. Rev. D, 58, 063501, doi: 10.1103/PhysRevD.58.063501
  * Jimenez et al. (2019) Jimenez, R., Cimatti, A., Verde, L., Moresco, M., & Wandelt, B. 2019, JCAP, 03, 043, doi: 10.1088/1475-7516/2019/03/043
  * Kochanek et al. (2001) Kochanek, C. S., Pahre, M. A., Falco, E. E., et al. 2001, ApJ, 560, 566, doi: 10.1086/322488
  * Komatsu et al. (2011) Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18, doi: 10.1088/0067-0049/192/2/18
  * Kudritzki & Puls (2000) Kudritzki, R.-P., & Puls, J. 2000, ARA&A, 38, 613, doi: 10.1146/annurev.astro.38.1.613
  * Kushnir et al. (2016) Kushnir, D., Zaldarriaga, M., Kollmeier, J. A., & Waldman, R. 2016, MNRAS, 462, 844, doi: 10.1093/mnras/stw1684
  * Linder (2003) Linder, E. V. 2003, Phys. Rev. Lett., 90, 091301, doi: 10.1103/PhysRevLett.90.091301
  * Lyke et al. (2020) Lyke, B. W., Higley, A. N., McLane, J. N., et al. 2020, ApJS, 250, 8, doi: 10.3847/1538-4365/aba623
  * MacLeod & Hogan (2008) MacLeod, C. L., & Hogan, C. J. 2008, Phys. Rev. D, 77, 043512, doi: 10.1103/PhysRevD.77.043512
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Makarov et al. (2014) Makarov, D., Prugniel, P., Terekhova, N., Courtois, H., & Vauglin, I. 2014, A&A, 570, A13, doi: 10.1051/0004-6361/201423496
  * Mandel et al. (2019) Mandel, I., Farr, W. M., & Gair, J. R. 2019, Mon. Not. Roy. Astron. Soc., 486, 1086, doi: 10.1093/mnras/stz896
  * Mapelli et al. (2021) Mapelli, M., Bouffanais, Y., Santoliquido, F., Arca Sedda, M., & Artale, M. C. 2021, arXiv e-prints, arXiv:2109.06222.  https://arxiv.org/abs/2109.06222
  * Mapelli et al. (2019) Mapelli, M., Giacobbo, N., Santoliquido, F., & Artale, M. C. 2019, MNRAS, 487, 2, doi: 10.1093/mnras/stz1150
  * Mastrogiovanni et al. (2021) Mastrogiovanni, S., Haegel, L., Karathanasis, C., Hernandez, I. M. n., & Steer, D. A. 2021, JCAP, 02, 043, doi: 10.1088/1475-7516/2021/02/043
  * Mastrogiovanni et al. (2021) Mastrogiovanni, S., Leyde, K., Karathanasis, C., et al. 2021, Phys. Rev. D, 104, 062009, doi: 10.1103/PhysRevD.104.062009
  * Matas et al. (2020) Matas, A., et al. 2020, Phys. Rev. D, 102, 043023, doi: 10.1103/PhysRevD.102.043023
  * Messenger & Read (2012) Messenger, C., & Read, J. 2012, Phys. Rev. Lett., 108, 091101, doi: 10.1103/PhysRevLett.108.091101
  * Mukherjee et al. (2021a) Mukherjee, S., Broadhurst, T., Diego, J. M., Silk, J., & Smoot, G. F. 2021a, Mon. Not. Roy. Astron. Soc., 501, 2451, doi: 10.1093/mnras/staa3813
  * Mukherjee et al. (2021b) —. 2021b, Mon. Not. Roy. Astron. Soc., 506, 3751, doi: 10.1093/mnras/stab1980
  * Mukherjee et al. (2021c) Mukherjee, S., Lavaux, G., Bouchet, F. R., et al. 2021c, Astron. Astrophys., 646, A65, doi: 10.1051/0004-6361/201936724
  * Mukherjee et al. (2021) Mukherjee, S., Wandelt, B. D., Nissanke, S. M., & Silvestri, A. 2021, Phys. Rev. D, 103, 043520, doi: 10.1103/PhysRevD.103.043520
  * Mukherjee et al. (2020) Mukherjee, S., Wandelt, B. D., & Silk, J. 2020, Mon. Not. Roy. Astron. Soc., 494, 1956, doi: 10.1093/mnras/staa827
  * Nair et al. (2018) Nair, R., Bose, S., & Saini, T. D. 2018, Phys. Rev. D, 98, 023502, doi: 10.1103/PhysRevD.98.023502
  * Nakamura (1998) Nakamura, T. T. 1998, Phys. Rev. Lett., 80, 1138, doi: 10.1103/PhysRevLett.80.1138
  * Namikawa et al. (2016) Namikawa, T., Nishizawa, A., & Taruya, A. 2016, Phys. Rev. Lett., 116, 121302, doi: 10.1103/PhysRevLett.116.121302
  * Ng et al. (2018a) Ng, K. K., Wong, K. W., Broadhurst, T., & Li, T. G. 2018a, Phys. Rev. D, 97, 023012, doi: 10.1103/PhysRevD.97.023012
  * Ng et al. (2018b) Ng, K. K. Y., Vitale, S., Zimmerman, A., et al. 2018b, Phys. Rev. D, 98, 083007, doi: 10.1103/PhysRevD.98.083007
  * Nishizawa (2017) Nishizawa, A. 2017, Phys. Rev. D, 96, 101303, doi: 10.1103/PhysRevD.96.101303
  * Nissanke et al. (2010) Nissanke, S., Holz, D. E., Hughes, S. A., Dalal, N., & Sievers, J. L. 2010, ApJ, 725, 496, doi: 10.1088/0004-637X/725/1/496
  * Norberg et al. (2002) Norberg, P., Cole, S., Baugh, C. M., et al. 2002, MNRAS, 336, 907, doi: 10.1046/j.1365-8711.2002.05831.x
  * Oguri (2016) Oguri, M. 2016, Phys. Rev. D, 93, 083511, doi: 10.1103/PhysRevD.93.083511
  * Oguri (2018) —. 2018, Mon. Not. Roy. Astron. Soc., 480, 3842, doi: 10.1093/mnras/sty2145
  * Ossokine et al. (2020) Ossokine, S., et al. 2020, Phys. Rev. D, 102, 044055, doi: 10.1103/PhysRevD.102.044055
  * Palmese et al. (2021) Palmese, A., Fishbach, M., Burke, C. J., Annis, J. T., & Liu, X. 2021, Astrophys. J. Lett., 914, L34, doi: 10.3847/2041-8213/ac0883
  * Palmese et al. (2020) Palmese, A., deVicente, J., Pereira, M. E. S., et al. 2020, ApJ, 900, L33, doi: 10.3847/2041-8213/abaeff
  * Perlmutter et al. (1999) Perlmutter, S., et al. 1999, Astrophys. J., 517, 565, doi: 10.1086/307221
  * Pratten et al. (2021) Pratten, G., et al. 2021, Phys. Rev. D, 103, 104056, doi: 10.1103/PhysRevD.103.104056
  * Renzo et al. (2020) Renzo, M., Farmer, R., Justham, S., et al. 2020, Astron. Astrophys., 640, A56, doi: 10.1051/0004-6361/202037710
  * Riess et al. (2019) Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D. 2019, Astrophys. J., 876, 85, doi: 10.3847/1538-4357/ab1422
  * Riess et al. (1996) Riess, A. G., Press, W. H., & Kirshner, R. P. 1996, Astrophys. J., 473, 88, doi: 10.1086/178129
  * Riess et al. (2016) Riess, A. G., et al. 2016, Astrophys. J., 826, 56, doi: 10.3847/0004-637X/826/1/56
  * Schneider et al. (1992) Schneider, P., Ehlers, J., & Falco, E. E. 1992, Gravitational Lenses (Springer), doi: 10.1007/978-3-662-03758-4
  * Schutz (1986) Schutz, B. F. 1986, Nature, 323, 310, doi: 10.1038/323310a0
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Soares-Santos et al. (2019) Soares-Santos, M., Palmese, A., Hartley, W., et al. 2019, ApJ, 876, L7, doi: 10.3847/2041-8213/ab14f1
  * Spergel et al. (2003) Spergel, D., et al. 2003, Astrophys. J. Suppl., 148, 175, doi: 10.1086/377226
  * Spergel et al. (2007) —. 2007, Astrophys. J. Suppl., 170, 377, doi: 10.1086/513700
  * Suyu et al. (2010) Suyu, S. H., Marshall, P. J., Auger, M. W., et al. 2010, ApJ, 711, 201, doi: 10.1088/0004-637X/711/1/201
  * Takahashi & Nakamura (2003) Takahashi, R., & Nakamura, T. 2003, Astrophys. J., 595, 1039, doi: 10.1086/377430
  * Talbot & Thrane (2018) Talbot, C., & Thrane, E. 2018, Astrophys. J., 856, 173, doi: 10.3847/1538-4357/aab34c
  * Taylor & Gair (2012) Taylor, S. R., & Gair, J. R. 2012, Phys. Rev. D, 86, 023502, doi: 10.1103/PhysRevD.86.023502
  * Taylor et al. (2012) Taylor, S. R., Gair, J. R., & Mandel, I. 2012, Phys. Rev. D, 85, 023535, doi: 10.1103/PhysRevD.85.023535
  * Thompson et al. (2020) Thompson, J. E., Fauchon-Jones, E., Khan, S., et al. 2020, Phys. Rev. D, 101, 124059, doi: 10.1103/PhysRevD.101.124059
  * Thrane & Talbot (2019) Thrane, E., & Talbot, C. 2019, PASA, 36, e010, doi: 10.1017/pasa.2019.2
  * Umeda et al. (2020) Umeda, H., Yoshida, T., Nagele, C., & Takahashi, K. 2020, Astrophys. J. Lett., 905, L21, doi: 10.3847/2041-8213/abcb96
  * van Son et al. (2021) van Son, L. A. C., de Mink, S. E., Callister, T., et al. 2021, arXiv e-prints, arXiv:2110.01634.  https://arxiv.org/abs/2110.01634
  * Vitale et al. (2020) Vitale, S., Gerosa, D., Farr, W. M., & Taylor, S. R. 2020.  https://arxiv.org/abs/2007.05579
  * Wang et al. (1996) Wang, Y., Stebbins, A., & Turner, E. L. 1996, Phys. Rev. Lett., 77, 2875, doi: 10.1103/PhysRevLett.77.2875
  * Weatherford et al. (2021) Weatherford, N. C., Fragione, G., Kremer, K., et al. 2021, ApJ, 907, L25, doi: 10.3847/2041-8213/abd79c
  * White et al. (2011) White, D. J., Daw, E. J., & Dhillon, V. S. 2011, Classical and Quantum Gravity, 28, 085016, doi: 10.1088/0264-9381/28/8/085016
  * Wong et al. (2020) Wong, K. C., et al. 2020, Mon. Not. Roy. Astron. Soc., 498, 1420, doi: 10.1093/mnras/stz3094
  * Ye & Fishbach (2021) Ye, C., & Fishbach, M. 2021, Phys. Rev. D, 104, 043507, doi: 10.1103/PhysRevD.104.043507
  * You et al. (2021) You, Z.-Q., Zhu, X.-J., Ashton, G., Thrane, E., & Zhu, Z.-H. 2021, Astrophys. J., 908, 215, doi: 10.3847/1538-4357/abd4d4
  * Yu et al. (2020) Yu, J., Wang, Y., Zhao, W., & Lu, Y. 2020, Mon. Not. Roy. Astron. Soc., 498, 1786, doi: 10.1093/mnras/staa2465
  * Zonca et al. (2019) Zonca, A., Singer, L., Lenz, D., et al. 2019, Journal of Open Source Software, 4, 1298, doi: 10.21105/joss.01298
  * Zou et al. (2019) Zou, H., Gao, J., Zhou, X., & Kong, X. 2019, ApJS, 242, 8, doi: 10.3847/1538-4365/ab1847

##  Appendix A Population prior models

###  A.1 Models for background cosmologies

We use a flat $\Lambda$CDM cosmological model with dark energy density as a
function of redshift $z$ described by Linder (2003)

$\rho_{\Lambda}(z)=\rho_{\Lambda,0}(1+z)^{3(1+w_{0})},$ (A1)

where $\rho_{\Lambda}$ is the dark energy density and $w_{0}$ is a
phenomenological parameter. If the dark energy density is constant during the
cosmic expansion, as it is in the standard cosmological model, then
$w_{0}=-1$. The luminosity distance is calculated as

$D_{\mathrm{L}}(z)=\frac{c(1+z)}{H_{0}}\int_{0}^{z}\frac{{\rm
d}z^{\prime}}{\sqrt{\Omega_{\rm
m}(1+z^{\prime})^{3}+\Omega_{\Lambda,0}(1+z^{\prime})^{3(1+w_{0})}}},$ (A2)

where $\Omega_{\rm m}$ and $\Omega_{\Lambda,0}$ are the present-day
dimensionless matter and dark energy densities respectively and
$\Omega_{\Lambda,0}=1-\Omega_{\rm m}$.

We consider two sets of priors for the cosmological background model which are
indicated in Table 5. For the first set of priors, we adopt a flat
$\Lambda$CDM cosmology that restricts the Hubble constant to only the range
compatible with the $H_{0}$ tension, while for the second set we adopt more
general, wide priors.

|  Restricted priors ($H_{0}$-tension) |  |
---|---|---|---
Parameter |  Description |  |  Prior
$H_{0}$ |  Hubble constant expressed in $\rm{km\,s^{-1}}\,Mpc^{-1}$ in the $H_{0}$-tension region. |  |  $\mathcal{U}$($65$, $77$)
$\Omega_{\rm m}$ |  Present-day matter density of the Universe fixed to the mean value inferred from measurements of the CMB in Ade et al. (2016) |  |  0.3065
$w_{0}$ |  Dark energy equation of state parameter fixed to the value that corresponds to a constant density. |  |  -1
|  Wide priors |  |
Parameter |  Description |  |  Prior
$H_{0}$ |  Hubble constant expressed in $\rm{km\,s^{-1}}\,Mpc^{-1}$ |  |  $\mathcal{U}$($10.0$, $200.0$)
$\Omega_{\rm m}$ |  Present-day matter density of the Universe. |  |  $\mathcal{U}$($0.0$,$1.0$)
$w_{0}$ |  Dark energy equation of state parameter. |  |  $\mathcal{U}$($-3.0$,$0.0$)
Table 5:  Summary of the priors on the cosmological parameters, for the two
sets of priors considered.

###  A.2 Merger rate and redshift distribution priors

We model the binary merger rate using a phenomenological model introduced with
the form of Madau & Dickinson (2014), motivated by the fact that the binary
formation rate might follow the star formation rate. The parameterization that
we use Callister et al. (2020) for the merger rate in the detector frame is

$\frac{{\rm d}N}{{\rm d}t_{\rm d}{\rm d}z}=R_{0}[1+(1+z_{\rm
p})^{-\gamma-k}]\frac{\partial V_{\rm c}}{\partial z}(\Lambda_{\rm
c})\frac{1}{1+z}\,\frac{(1+z)^{\gamma}}{1+[(1+z)/(1+z_{\rm p)}]^{\gamma+k}},$
(A3)

where $R_{0}$ is the binary merger rate today, $V_{\rm c}$ is the comoving
volume, $\gamma$ and $k$ are the slopes of the two power-law regimes before
and after a turning point $z_{\rm p}$ and $\Lambda_{\rm c}$ are a set of
parameters describing the cosmological expansion. The extra $1/(1+z)$ factor
encodes the clock difference in the source and detector, while the factor
$(1+z_{\rm p})$ ensures that today the merger rate is $R(z=0)=R_{0}$. The
redshift prior, normalized over the redshift, can be expressed as

$\pi(z|\gamma,\kappa,z_{\rm p},\Lambda_{\rm c})=\frac{1}{C}\frac{{\rm
d}N}{{\rm d}t_{\rm d}{\rm d}z},$ (A4)

where $C$ is the normalization factor calculated from Eq. (A3).

The priors that we use on the merger rate hyper-parameters are indicated in
Table 6. The prior ranges that we model are wide enough to include the effect
of a possible time delay between the formation and the merger of the binary.

Parameter |  Description |  |  Prior
---|---|---|---
$R_{0}$ |  BBH merger rate today in Gpc-3 yr-1 |  |  $\mathcal{U}$($0.0$, $100.0$)
$\gamma$ |  Slope of the powerlaw regime for the rate evolution before the point $z_{p}$ |  |  $\mathcal{U}$($0.0$, $12.0$)
$k$ |  Slope of the powerlaw regime for the rate evolution after the point $z_{\rm p}$ |  |  $\mathcal{U}$($0.0$, $6.0$)
$z_{p}$ |  Redshift turning point between the powerlaw regimes with $\gamma$ and $k$ |  |  $\mathcal{U}$($0.0$, $4.0$)
Table 6:  Summary of the prior hyper-parameters used for the merger rate
evolution models adopted in this paper.

###  A.3 Phenomenological Mass priors

The three phenomenological mass models that we implement are a superposition
of two probability density distributions and are compatible with the
phenomenological priors used for BBHs in Abbott et al. (2021c, d), although
the prior ranges on the population parameters are different. The first is a
truncated power law $\mathcal{P}(x|x_{\rm min},x_{\rm max},\alpha)$ described
by slope $\alpha$, and lower and upper bounds $x_{\rm min},x_{\rm max}$ at
which there is a hard cut-off

$\mathcal{P}(x|x_{\rm min},x_{\rm
max},\alpha)\propto\begin{cases}x^{\alpha}&\left(x_{\rm min}\leqslant
x\leqslant x_{\rm max}\right)\\\ 0&\mathrm{Otherwise}.\end{cases}$ (A5)

The second is a Gaussian distribution with mean $\mu$ and standard deviation
$\sigma$,

$\mathcal{G}(x|\mu,\sigma,a,b)=\frac{1}{\sigma\sqrt{2\pi}}\exp{\left[-\frac{(x-\mu)^{2}}{2\sigma^{2}}\right]}.$
(A6)

The source mass priors for the BBHs population that we consider are factorized
as

$\pi(m_{1},m_{2}|\Phi_{m})=\pi(m_{1}|\Phi_{m})\pi(m_{2}|m_{1},\Phi_{m}),$ (A7)

where $\pi(m_{1}|\Phi_{m})$ is the distribution of the primary mass component
while $\pi(m_{2}|m_{1},\Phi_{m})$ is the distribution of the secondary mass
component given the primary. For all of the mass models, the secondary mass
component $m_{2}$ is described with a truncated power-law with slope $\beta$
between a minimum mass $m_{\rm min}$ and a maximum mass $m_{1}$

$\pi(m_{2}|m_{1},m_{\rm min},\alpha)=\mathcal{P}(m_{2}|m_{\rm
min},m_{1},\beta),$ (A8)

while the primary mass is described with several models discussed in the
following paragraphs.

For some of the phenomenological models, we also apply a smoothing factor to
the lower end of the mass distribution

$\pi(m_{1},m_{2}|\Phi_{m})=[\pi(m_{1}|\Phi_{m})\pi(m_{2}|m_{1},\Phi_{m})]S(m_{1}|\delta_{m},m_{\rm
min})S(m_{2}|\delta_{m},m_{\rm min}),$ (A9)

where $S$ is a sigmoid-like window function that adds a tapering of the lower
end of the mass distribution. See Eq. (B6) and Eq. (B7) of Abbott et al.
(2021c) for the explicit expression for the window function.

The three phenomenological mass models are highlighted in the following list.
In Table 7, we report the prior ranges used for the population hyper-
parameters.

  * •

Truncated model: It describes the distribution of the primary mass $m_{1}$
with a truncated power law with slope $-\alpha$ between a minimum mass $m_{\rm
min}$ and a maximum mass $m_{\rm max}$.

$\pi(m_{1}|m_{\rm min},m_{\rm max},\alpha)=\mathcal{P}(m_{1}|m_{\rm
min},m_{\rm max},-\alpha).$ (A10)

  * •

Power Law + Peak model: It describes the primary mass component as a
superposition of a truncated PL, with slope $-\alpha$ between a minimum mass
$m_{\rm min}$ and a maximum mass $m_{\rm max}$, plus a Gaussian component with
mean $\mu_{\rm g}$ and standard deviation $\sigma_{\rm g}$,

$\pi(m_{1}|m_{\rm min},m_{\rm max},\alpha,\lambda_{\rm g},\mu_{\rm
g},\sigma_{\rm g})=[(1-\lambda_{\rm g})\mathcal{P}(m_{1}|m_{\rm min},m_{\rm
max},-\alpha)+\lambda_{\rm g}\mathcal{G}(m_{1}|\mu_{\rm g},\sigma_{\rm g})].$
(A11)

  * •

Broken Power Law model: It describes the distribution of $m_{1}$ as a PL
between a minimum mass $m_{\rm min}$ and a maximum mass $m_{\rm max}$. The
Broken Power Law model is characterized by two PL slopes $\alpha_{1}$ and
$\alpha_{2}$ and by a breaking point between the two regimes at $m_{\rm
break}=b(m_{\rm max}-m_{\rm min})$, where $b$ is a number $\in[0,1]$. The
broken PL model is

$\pi(m_{1}|m_{\rm min},m_{\rm
max},\alpha_{1},\alpha_{2})=\mathcal{P}(m_{1}|m_{\rm min},m_{\rm
break},-\alpha_{1})+\frac{\mathcal{P}(m_{\rm break}|m_{\rm min},m_{\rm
break},-\alpha_{1})}{\mathcal{P}(m_{\rm break}|m_{\rm break},m_{\rm
max},-\alpha_{2})}\mathcal{P}(m_{1}|b,m_{\rm max},-\alpha_{2}).$ (A12)

|  Truncated |  |
---|---|---|---
Parameter |  Description |  |  Prior
$\alpha$ |  Spectral index for the PL of the primary mass distribution. |  |  $\mathcal{U}$($1.5$, $12.0$)
$\beta$ |  Spectral index for the PL of the mass ratio distribution. |  |  $\mathcal{U}$($-4.0$, $12.0$)
$m_{\rm min}$ |  Minimum mass of the PL component of the primary mass distribution. |  |  $\mathcal{U}$($2.0\,\rm{M}_{\odot}$, $10.0\,\rm{M}_{\odot}$)
$m_{\rm max}$ |  Maximum mass of the PL component of the primary mass distribution. |  |  $\mathcal{U}$($50.0\,\rm{M}_{\odot}$, $200.0\,\rm{M}_{\odot}$)
|  Power Law + Peak |  |
Parameter |  Description |  |  Prior
$\alpha$ |  Spectral index for the PL of the primary mass distribution. |  |  $\mathcal{U}$($1.5$, $12.0$)
$\beta$ |  Spectral index for the PL of the mass ratio distribution. |  |  $\mathcal{U}$($-4.0$, $12.0$)
$m_{\rm min}$ |  Minimum mass of the PL component of the primary mass distribution. |  |  $\mathcal{U}$($2.0\,\rm{M}_{\odot}$, $10.0\,\rm{M}_{\odot}$)
$m_{\rm max}$ |  Maximum mass of the PL component of the primary mass distribution. |  |  $\mathcal{U}$($50.0\,\rm{M}_{\odot}$, $200.0\,\rm{M}_{\odot}$)
$\lambda_{\rm g}$ |  Fraction of the model in the Gaussian component. |  |  $\mathcal{U}$($0.0$, $1.0$)
$\mu_{\rm g}$ |  Mean of the Gaussian component in the primary mass distribution. |  |  $\mathcal{U}$($20.0\,\rm{M}_{\odot}$, $50.0\,\rm{M}_{\odot}$)
$\sigma_{\rm g}$ |  Width of the Gaussian component in the primary mass distribution. |  |  $\mathcal{U}$($0.4\,\rm{M}_{\odot}$, $10.0\,\rm{M}_{\odot}$)
$\delta_{m}$ |  Range of mass tapering at the lower end of the mass distribution. |  |  $\mathcal{U}$($0.0\,\rm{M}_{\odot}$, $10.0\,\rm{M}_{\odot}$)
|  Broken Power Law |  |
Parameter |  Description |  |  Prior
$\alpha_{1}$ |  PL slope of the primary mass distribution for masses below $m_{\mathrm{break}}$. |  |  $\mathcal{U}$($1.5$, $12.0$)
$\alpha_{2}$ |  PL slope for the primary mass distribution for masses above $m_{\mathrm{break}}$. |  |  $\mathcal{U}$($1.5$, $12.0$)
$\beta$ |  Spectral index for the PL of the mass ratio distribution. |  |  $\mathcal{U}$($-4.0$, $12.0$)
$m_{\rm min}$ |  Minimum mass of the PL component of the primary mass distribution. |  |  $\mathcal{U}$($2.0\,\rm{M}_{\odot}$, $50.0\,\rm{M}_{\odot}$)
$m_{\mathrm{max}}$ |  Maximum mass of the primary mass distribution. |  |  $\mathcal{U}$($50.0\,\rm{M}_{\odot}$, $200.0\,\rm{M}_{\odot}$)
$b$ |  The fraction of the way between $m_{\text{min}}$ and $m_{\text{max}}$ at which the primary mass distribution breaks. |  |  $\mathcal{U}$($0.0$,$1.0$)
$\delta_{m}$ |  Range of mass tapering on the lower end of the mass distribution. |  |  $\mathcal{U}$($0.0\,\rm{M}_{\odot}$, $10.0\,\rm{M}_{\odot}$)
Table 7:  Summary of the priors used for the population hyper-parameters for
the three phenomenological mass models.

##  Appendix B Full results from the population analysis and effect of
different SNR cuts

We provide extra details on the joint inference of cosmological and population
parameters using BBHs.

In Fig. 13 we show the corner plots for the posterior associated with the
Power Law + Peak model (the one implemented for the galaxy catalog analysis)
adopting wide priors on the cosmological parameters and for the population of
BBHs. As was also shown in Fig. 5, the main mass-related population parameters
correlating with the cosmological parameters, and in particular $H_{0}$, are
the position of the Gaussian component (BBH excess in the mass distribution)
and the higher end of the source mass distribution.

The other parameter that correlates with the estimation of $H_{0}$ is the rate
evolution parameter $\gamma$. This parameter models a power-law increasing
merger rate with the redshift. We find that higher values of $\gamma$ support
lower values of $H_{0}$, which is due to the fact that lowering $H_{0}$ will
place events at lower redshifts, which are incompatible with the observed mass
distribution; therefore $\gamma$ tries to correct this by supporting higher
redshifts. However, the posterior on the rate evolution is well within the
statistical uncertainties given in (Abbott et al., 2021d).

Figure 13: Corner plots for the preferred Power Law + Peak model parameters
and cosmological parameters, fitted to BBHs with SNR $>11$.

We run additional systematic studies using different SNR cuts for the
selection of the BBHs to include in our analysis. We explore a higher SNR cut
of 12 (more pure) that selects a sample of $35$ events. We also explore a
lower SNR cut of 10, allowing $60$ events with a IFAR$>0.5$ yr and no
plausible instrumental origin. For this set of events, GW190426_152155 and
GW190531_023648 are excluded as their secondary mass extends to lower masses
in the NS region.

The marginal $H_{0}$ posterior for all the mass models is shown in Fig. 14,
where we show that including more events always produces a posterior on the
$H_{0}$ within the statistical uncertainties of other selection criteria. In
all of these cases, the excess of BBHs around $35\rm{M}_{\odot}$ is present
for all the SNR cuts and it is responsible for the preference observed in the
$H_{0}$ posterior.

Figure 14: Marginal posterior probability distributions for $H_{0}$ with using
BBH events with three different SNR cuts and the three mass models.

##  Appendix C Schechter luminosity function studies

In this Appendix we show comparisons of the Schechter luminosity function (LF)
for the galaxies with the K and $B_{J}$ bands galaxies reported in GLADE+.
Wrong assumptions on the LF, or incorrect description of the selection biases,
could potentially introduce a bias in the inferred $H_{0}$. Indeed, one of the
key assumptions that we have made to construct our completeness corrections is
that the galaxy catalog is _magnitude limited_ , i.e. galaxies are not
detected only because they are too faint. This cannot be the case if another
selection bias (based on e.g., colors or spectral features) were present, or
if color and evolution corrections are not implemented properly.

In Fig. 15 we show a comparison of the assumed LF and the number density of
galaxies per comoving volume present in GLADE+. In the case that the galaxy
catalog can be correctly described as magnitude-limited, we expect that the
distribution of the GLADE+ galaxies distribution will match the assumed LF at
its bright end, and then will start to decrease when we reach (and exceed) the
corresponding absolute magnitude threshold. Galaxies in the K–band are well
described by this behaviour and missing galaxies can be explained by the
impact of the apparent magnitude threshold, while for the $B_{J}$–band there
seems to be some additional missing galaxies at low redshift. This observed
behaviour motivated our decision to present our main results using the K–band
magnitudes compiled in the GLADE+ catalog.

Figure 15: Left: Comparison of the assumed $B_{J}$–band galaxy luminosity
function (pink dashed line) and the differential number density of galaxies in
different redshift bins (solid colored lines). The vertical dashed line
indicates the median absolute magnitude threshold for galaxy detection, as
computed in our code. Right: Same comparison, but for the K–band luminosity
function. In each panel the solid lines are calculated based on the median
apparent magnitude threshold, computed over all sky directions, then
collecting and making a histogram of all galaxies in the GLADE+ catalog that
lie in the appropriate redshift bin and are brighter than the apparent
magnitude threshold.


Reconsidering Dependency Networks from an Information Geometry Perspective

Kazuya Takabatake<EMAIL_ADDRESS>
Human Informatics and Interaction Research Institute

National Institute of Advanced Industrial Science and Technology

Central 2: 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan.
Shotaro Akaho<EMAIL_ADDRESS>
Human Informatics and Interaction Research Institute

National Institute of Advanced Industrial Science and Technology

Central 2: 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan.

Dependency networks [Heckerman et al., 2000] are potential probabilistic graphical models for systems comprising a large number of variables.
Like Bayesian networks, the structure of a dependency network is represented by a directed graph, and each node has a conditional probability table.
Learning and inference are realized locally on individual nodes; therefore, computation remains tractable even with a large number of variables.
However, the dependency network's learned distribution is the stationary distribution of a Markov chain called pseudo-Gibbs sampling and has no closed-form expressions.
This technical disadvantage has impeded the development of dependency networks.
In this paper, we consider a certain manifold for each node.
Then, we can interpret pseudo-Gibbs sampling as iterative m-projections onto these manifolds.
This interpretation provides a theoretical bound for the location where the stationary distribution of pseudo-Gibbs sampling exists in distribution space.
Furthermore, this interpretation involves structure and parameter learning algorithms as optimization problems.
In addition, we compare dependency and Bayesian networks experimentally.
The results demonstrate that the dependency network and the Bayesian network have roughly the same performance in terms of the accuracy of their learned distributions.
The results also show that the dependency network can learn much faster than the Bayesian network.

Dependency Networks, Graphical models, Information Geometry, Pseudo-Gibbs sampling, Learning

§ INTRODUCTION

The primary purpose of this paper is to investigate dependency networks [Heckerman et al., 2000] from an information geometry perspective and offers new insights into dependency networks.

Bayesian and Markov probabilistic graphical models <cit.> are frequently employed to construct multivariate probabilistic models.
Both Bayesian and Markov graphical models have well-studied theoretical foundations, and their structure/parameter learning and inference algorithms are also well-established.
Here, let $X=\{X_0,...,X_{n-1}\}$ be variables in the entire network and $p^*(X)$ be the true distribution underlying training data.
Then, the basic mechanisms of Bayesian and Markov networks are as follows:

Learning. Given training data, the entire network learns $p^*(X)$, that is, the entire network builds a distribution $\pi(X)$ that is an estimation of $p^*(X)$.

Inference. Given inference condition $V=v$, the entire network computes the posterior distribution $\pi(U|V=v)$, where $U,V$ are certain variables in the network.

In Bayesian network structure learning <cit.>, the network searches the graphical structure $G$ that minimizes $scost(G)$ for a structure cost function $scost$.
The set of possible $G$ becomes extremely large as the number of nodes increases.
Even using greedy algorithms, the network must evaluate vast $scost(G)$ values, which leads to very slow structure learning.

In Markov network parameter learning <cit.>, the network searches parameters $\theta=\{\theta_j\}(j\in J)$ that minimize $pcost(\theta)$ for a parameter cost function $pcost$.
$|J|$(=number of parameters) becomes large as the number of nodes increases.
In such cases, parameter learning becomes a minimization problem in a very high-dimensional space that is intractable even for high-dimensional minimization algorithms such as L-BFGS <cit.>.

Dependency network graphical models [Heckerman et al., 2000] hold potential for large systems.
Like Bayesian networks, the structure of a dependency network is represented by a directed graph, and each node has a conditional probability table (CPT).
However, differing from Bayesian networks, dependency network graphs can be cyclic.

In dependency network structure learning, the structure cost function $scost$ can be divided into local costs as follows:

$$scost(G)=\sum_i scost_i(Y_i),$$

where $Y_i$ is the local structure of the node $i$ (Section <ref>).
In this case, we can minimize $scost$ by independently minimizing $scost_i$.
Since the search space for $Y_i$ is much smaller than the search space for $G$,
structure learning remains tractable even in the cases where the number of nodes is large.

In dependency network parameter learning, we can also divide the parameter cost function $pcost$ into local costs as follows:

$$pcost(\theta)=\sum_i pcost_i(\theta_i).$$

Therefore, we can perform parameter learning by independent “local" learning of $\theta_i$ (Section <ref>).
This independence of local learning keeps parameter learning tractable even with a large number of nodes [Bayesian network parameter learning also has this property.].

Let $X_{-i}$ be variables in the entire network except for $X_i$, and let $Y_i$ be certain variables included in $X_{-i}$.
For a distribution $p(X)$, $p(X_i|X_{-i})$ is called the full conditional distribution <cit.>.
Now, the framework of the dependency networks can be described as follows:

Learning. Given training data, each node learns its full conditional distribution independently, that is, each node builds CPT $\theta_i(X_i|Y_i)$ that is an estimation of $p^*(X_i|X_{-i})$.

Inference. Given CPTs and an inference condition $V=v$, a Markov chain Monte Carlo method called pseudo-Gibbs sampling [Heckerman et al., 2000] draws samples from its stationary distribution[We refer to these samples as output data (or shortly outputs) of pseudo-Gibbs sampling.], and this stationary distribution becomes the estimation of $p^*(U|V=v)$.
Users know the estimation by counting samples in the output data.

Learning the full conditional distribution is a local task of each node; therefore, it is much more tractable than learning the joint distribution.
Various regression techniques can be applied to learn full-conditional distributions.
If a user use a regression technique with input variable selection, the user can ignore most variables in $X_{-i}$ and can approximate $p^*(X_i|X_{-i})$ by $\theta_i(X_i|Y_i)$, where $Y_i$ consists of only several variables.

The inference task is also tractable even in cases where the number of nodes is large because pseudo-Gibbs sampling consists of simple local tasks (Section <ref>).

Although the dependency network mechanism consists of such simple local tasks, the distribution represented by the dependency network is complicated.
Let $\pi(X)$ be the joint distribution represented by a graphical model.
For Bayesian networks, $\pi(X)$ is expressed as follows:

\begin{equation}\label{eq:BN represent}
\pi(X)=\prod_i \pi(X_i|Y_i),
\end{equation}

where $Y_i$ denotes the parents of $X_i$ <cit.>.
For Markov networks, $\pi(X)$ is expressed as follows:

\begin{equation}\label{eq:MN represent}
\pi(X)=\frac1Z\prod_c\phi_c(X_c),
\end{equation}

where $c$ denotes clique in the network, $X_c$ denotes variables included in $c$, $\phi_c$ denotes clique potential, and $Z$ denotes partition function <cit.>.
In contrast, the joint distribution represented by a dependency network is the stationary distribution of a Markov chain called pseudo-Gibbs sampling [Heckerman et al., 2000], and this stationary distribution has no closed-form expressions such as Eqs.(<ref>) and (<ref>).
This disadvantage is significant; for example, the conventional maximum likelihood estimation method is not applicable because the likelihood itself is unknown.
The progress of dependency networks has been disturbed by this problematic property.

Heckerman et al., 2000 investigated a specific class of dependency networks, called consistent dependency networks, which are equivalent to Markov networks, and explained general dependency networks by perturbation from consistent dependency networks.
They showed a theoretical bound for the difference between consistent dependency networks and general dependency networks.
However, they also found that their theoretical bound was too loose in practice.

Takabatake and Akaho, 2012 independently studied graphical models, which they referred to as firing process networks, which are equivalent to general dependency networks.
They considered a certain manifold for each node and interpreted pseudo-Gibbs sampling as iterative m-projections onto these manifolds.
This interpretation offers the following consequences:

* Providing a theoretical bound for the location where the stationary distribution exists in distribution space;

* Providing how to construct structure/parameter learning algorithms as optimization problems.

The remainder of this paper is organized as follows.
In Section 2, we review the structure and mechanism of dependency networks.
In Section 3, from a perspective of information geometry, we discuss pseudo-Gibbs sampling, which is a Markov chain Monte Carlo method to synthesize a joint distribution from CPTs owned by the nodes.
In Section 4, we derive structure and parameter learning algorithms based on the interpretation obtained in Section 2.
Probabilistic inference in dependency networks is discussed in Section 5.
In Section 6, we generalize the input variables referred by CPT to enhance the expressive power of dependency networks.
In Section 7, we compare dependency and Bayesian networks and examines the behavior of dependency network learning experimentally.
Conclusions and suggestions for future work are provided in Section 8.
Appendix section provides proofs of theorems.

§ DEPENDENCY NETWORKS

We explain some notations used in this paper.
A random variable is denoted by a capital letter, such as $X$, and a value that $X$ takes is denoted by a lower case letter, such as $x$.
Concatenated variables such as $(X,Y,Z)$ are shortened and denoted $XYZ$.
Note that concatenated variables are often treated as a single variable.
When $X$ is a concatenation of certain variables, the set of the variables included in $X$ are also denoted by $X$, for example, if $X=X_0X_1X_2X_3$ and $Y=X_1X_2$ then $Y\subseteq X$ and $X\setminus Y=X_0X_3$.
$\left<f(X)\right>_{p(X)}$ denotes the expectation of $f(X)$ with distribution $p(X)$, that is, $\left<f(X)\right>_{p(X)}=\sum_x p(x)f(x)$.
For a distribution $p$, an empirical distribution of data drawn from $p$ is denoted as $\tilde p$.

§.§ Dependency Network Structure

Node $i$

The structure of dependency networks resembles the structure of Bayesian networks.
A dependency network is represented by a directed graph with $n$ nodes indexed by the numbers $0,...,n-1$.
As shown in Figure <ref>, each node has a variable $X_i$ and a conditional probability table (CPT) $\theta_i(X_i|Y_i)$, where $Y_i$ denotes certain variables in $X_{-i}$.
We call $Y_i$ input variables (or simply inputs) of node $i$.

It is possible to consider that $X_i$ takes continuous values.
In this case, $\theta_i$ is a density function, and CPT is implemented using a regression technique for continuous variables.
However, for simplicity, we focus on the cases where $X_i$ takes discrete finite values in this paper.

The dependency network graph consists of edges directed from node $i$'s input variables to node $i$.
Unlike the Bayesian network, the graph of a dependency network can be cyclic.

§.§ Pseudo-Gibbs Sampling

Here, we explain how to synthesize a joint distribution from CPTs owned by the nodes.

Pseudo-Gibbs sampling [Heckerman et al., 2000] is a Markov chain Monte Carlo method that is applied to synthesize a joint distribution $\pi(X)$ from CPTs owned by the nodes.
The following algorithm is an example of pseudo-Gibbs sampling.
Pseudo-Gibbs sampling
$X=X_0...X_{n-1}$: variables in the entire network

$o$: output data

$N$: number of samples to be drawn
[1]
$o\leftarrow\emptyset$, $X\leftarrow x^0$ arbitrary initial value
$t=0$ $N-1$
Append $X$ into $o$
Select a node $i$
Referring to inputs $Y_i$, draw a sample $x_i$ from distribution $\theta_i(X_i|y_i)$ and let $X_i\leftarrow x_i$

There are two ways to select a node in line <ref>.
One is to select a node in a sequential and cyclic manner such as $0,1,...,n-1,0,1,...$.
We refer to this type of pseudo-Gibbs sampling as ordered-pseudo-Gibbs sampling.
The second option is to randomly select a node $i$ with probability $c_i$ such that

\begin{equation}\label{eq:c_i}
\sum_i c_i=1.
\end{equation}

Typically, $c_i=1/n$ is used unless there is a particular reason (Section <ref>).
We refer to this type of pseudo-Gibbs sampling as random-pseudo-Gibbs sampling.
The operation in line <ref> is mentioned frequently in this paper; therefore, we will refer to this operation as firing node $i$.

Let $O^N=\{X^0,...,X^{N-1}\}$ be $N$ samples of output data—as a sequence of random variables—of pseudo-Gibbs sampling algorithm.
Here, we assume the ergodicity of random-pseudo-Gibbs sampling.
Then, the following consequences hold <cit.>.

* There uniquely exists a distribution $\pi(X)$, called stationary distribution, such that for any initial value $x^0$,

\begin{equation}\label{eq:stationary}
\pi(X)=\lim_{t\to\infty}p(X^t|x^0).
\end{equation}

* For any real valued function $f(X)$ that satisfies $\left<|f(X)|\right>_{\pi(X)}<\infty$,

\begin{equation}\label{eq:ergodicity}
\lim_{N\to\infty}\frac{f(X^0)+...+f(X^{N-1})}N=\left<f(X)\right>_{\pi(X)}\quad(a.s.).
\end{equation}

Substituting $1(X=x)$ (indicator function) for $f$ in Eq.(<ref>), we obtain

\begin{equation}\label{eq:represented}
\lim_{N\to\infty}\frac{N_x}N=\pi(x)\quad(a.s.),
\end{equation}

where $N_x$ denotes the number of occurrences of $X^t=x(t=0,...,N-1)$ in the output dataset.
Equation (<ref>) implies that, even if it is not possible to know $\pi(x)$ directly, we can estimate $\pi(x)$ by a Monte Carlo method.

Ordered-pseudo-Gibbs sampling is an inhomogeneous Markov chain <cit.> and does not have a stationary distribution.
Even in this case, the subsequence $X^iX^{i+n}X^{i+2n}...$ becomes a homogeneous Markov chain and has a stationary distribution $\pi_i(X)$ under the assumption of ergodicity.
Therefore, the limit in the left side of Eq.(<ref>) still exists and

\begin{equation}\label{eq:ordered}
\lim_{N\to\infty}\frac{N_x}N=\frac{\pi_0(x)+...+\pi_{n-1}(x)}n.
\end{equation}

Contrary to pseudo-Gibbs sampling, we refer to usual Gibbs sampling as actual-Gibbs sampling in this paper.
Actual-Gibbs sampling is a special case of the pseudo-Gibbs sampling such that there exists a joint distribution $p(X)$ that satisfies

\begin{equation}\label{eq:Gibbs}
\theta_i(X_i|Y_i)=p(X_i|X_{-i})
\end{equation}

for all $i$.
In this case, $p$ becomes the stationary distribution because firing any node $i$ does not move $p$.

Here, we build the following rough hypothesis.

Hypothesis 1. If every CPT is a good estimation of $p^*(X_i|X_{-i})$, then the stationary distribution of pseudo-Gibbs sampling becomes a good estimation of $p^*(X)$.

Hypothesis 1 plays the central role in dependency network learning.
We will justify this hypothesis in Section <ref>.

§.§ Local Learning

The learning task involves building a distribution $\pi(X)$ that is a good estimation of the true distribution $p^*(X)$.
In dependency networks, controlling the stationary distribution $\pi(X)$ by adjusting CPTs is difficult.
Therefore, rather than learning $p^*(X)$ directly, each node in the dependency network learns the full-conditional distribution $p^*(X_i|X_{-i})$ independently from training data.
A favorable property here is that node $i$ can learn $p^*(X_i|X_{-i})$ without knowing other nodes' behavior.

Local learning is realized using regression algorithms with input variable selection.
For example, [Heckerman et al., 2000] used a decision tree algorithm.
In the process of constructing the tree, variables $Y_i$ are automatically selected from $X_{-i}$.

We will present our proposed learning algorithm in Section <ref>.

§.§ Probabilistic inference by Pseudo-Gibbs sampling

Let $U$ and $V$ be certain variables in the network.
Probabilistic inference involves estimating the posterior distribution $p^*(U|v)$.
Estimating $p^*(X)$ described in Section <ref> is a special case where $U=X,V=\emptyset$.

We can use pseudo-Gibbs sampling again to estimate $p^*(U|v)$ here.
After learning CPTs, the value of $V$ is clamped to $v$, and only the nodes out of $V$ are fired in pseudo-Gibbs sampling.
We refer to this type of pseudo-Gibbs sampling as clamped-pseudo-Gibbs sampling[In contrast, we refer to the pseudo-Gibbs sampling described in Section <ref> as free-pseudo-Gibbs sampling.].
Then, the stationary distribution of clamped-pseudo-Gibb sampling is used as the estimation of $p^*(U|v)$.

The following algorithm is an example of clamped-pseudo-Gibbs sampling.
Clamped-pseudo-Gibbs sampling
$X=X_0...X_{n-1}$: variables in the entire network

$V$: variables to be clamped to $v$

$\overline V$: all variables out of $V$

$o$: output dataset

$N$: number of samples to be drawn
[1]
$o\leftarrow\emptyset$, $\overline V\leftarrow \overline v^0$ arbitrary initial value, $V\leftarrow v$
$t=0$ $N-1$
Append $X$ into $o$
Select a node $i\not\in V$
Fire node $i$

In line <ref>, if we select a node in a random manner (random-pseudo-Gibbs sampling) then the probability to select node $i$ becomes as follows:

\frac1{n-|V|} &i\in V\\
0 &i\not\in V
\end{cases},$$

where $|V|$ is the number of nodes in $V$.

Let define as $\overline V:=X\setminus V$ and $Z_i:=X_{-i}\setminus V$.
Note that $\theta_i(X_i|Y_i)$, which is a function of $X_iY_i$, is a special case of $\theta_i(X_i|X_{-i})$, which is a function of $X$, because $X_iY_i\subseteq X$, that is, we can always rewrite $\theta_i(X_i|Y_i)$ as $\theta_i(X_i|X_{-i})$.
Then, the clamped-pseudo-Gibbs sampling is equivalent to the free-pseudo-Gibbs sampling with CPTs $\{\theta_i(X_i|Z_iv)\}(i\in\overline V)$.
Let $\pi_v(\overline V)$ be	the stationary distribution of clamped-pseudo-Gibbs sampling.
If $\theta_i(X_i|X_{-i})$ is a good estimation of $p^*(X_i|X_{-i})$, then $\theta_i(X_i|Z_iv)$ is also considered a good estimation of $p^*(X_i|Z_iv)$.
Therefore, $\pi_v(\overline V)$ is considered  a good estimation of $p^*(\overline  V|v)$ by Hypothesis 1.
Let $N_u$ be the number of samples such that $U=u$ in the $N$ samples of output data of clamped-pseudo-Gibbs sampling.
Then, $N_u/N$ becomes a good estimation of $p^*(u|v)$ for sufficiently large $N$.

Output data obtained by Algorithm <ref> are far from i.i.d. data.
To make output data closer to i.i.d. data, we use burn-in and thinning techniques <cit.>, that is, we discard the first $b$ samples and use only every $k$-th samples.
Then, Algorithm <ref> is modified as follows:
Clamped-pseudo-Gibbs sampling with burn-in and thinning
$X=X_0...X_{n-1}$: variables in the entire network

$V$: variables to be clamped to $v$

$\overline V$: all variables out of $V$

$o$: output dataset

$N$: number of samples to be drawn

$b$: The first $b$ samples are discarded from the output.

$k$: Every $k$-th samples are used for outputs.
[1]
$o\leftarrow\emptyset$, $X\leftarrow x^0$ arbitrary initial value, $V\leftarrow v$
$s=1$ $b$
Select a node $i\not\in V$
Fire node $i$
$t=0$ $N-1$
Append $X$ into $o$
$s=1$ $k$
Select a node $i\not\in V$
Fire node $i$

§.§.§ Lack of Joint Distribution in Inference

An ordinary learning machine constructs a joint distribution $\pi(X)$ as the estimation of the true distribution $p^*(X)$ in the learning phase and use $\pi(U|v)$ as the estimation of the posterior distribution $p^*(U|v)$ in the inference phase.
However, this approach often becomes intractable as the number of variables increases.

In contrast, the dependency network constructs $\pi_v(U)$—the stationary distribution of clamped-pseudo-Gibbs sampling— every time an inference condition $V=v$ is given.
Therefore, any joint distribution $p(X)$ such that satisfies $\forall v,p(U|v)=\pi_v(U)$ does not exist.
Inference by pseudo-Gibbs sampling involves a trade-off.
Computation tractability is achieved; however, the existence of a joint distribution consistent for all inference queries is not achieved.
However, in many practical applications, users want to estimate the distribution $p^*(U|v)$ for each condition $V=v$, and the lack of joint distribution consistent for all inference queries does not become a significant disadvantage.

§ INFORMATION GEOMETRY OF PSEUDO-GIBBS SAMPLING

In Section <ref>, we reviewed dependency networks and showed that pseudo-Gibbs sampling could be used for a tool to synthesize a joint distribution from CPTs.
In this section, we examine pseudo-Gibbs sampling from the perspective of information geometry.
As a result, we can interpret pseudo-Gibbs sampling as iterative m-projection onto certain manifolds.
Given this interpretation, we were motivated to consider making the stationary distribution close to any desired distribution.

We here prepare some definitions and consequences in information geometry.
In information geometry, a distribution is represented by a point in a distribution space $\mathcal P:=\{p|p(X)\}$.
For two distributions $p_0,p_1\in\mathcal P$, the one-dimensional manifold

$$\{p_\lambda|p_\lambda(X)=(1-\lambda)p_0(X)+\lambda p_1(X),\lambda\in[0,1]\}$$

is called m-geodesic between $p_0,p_1$, and the one-dimensional manifold

\begin{align*}
\left\{p_\lambda\left|\ln p_\lambda(X)=(1-\lambda)\ln p_0(X)+\lambda\ln p_1(X)-\ln Z,\lambda\in[0,1]\right.\right\}\\
Z=\sum_x e^{(1-\lambda)\ln p_0(x)+\lambda\ln p_1(x)}
\end{align*}

is called e-geodesic between $p_0,p_1$ <cit.>.
For a manifold $\mathcal M$, if $\mathcal M$ includes any e-geodesic/m-geodesic between any two distributions $p_0,p_1\in\mathcal M$ then $\mathcal M$ is said to be e-flat/m-flat.
An m-geodesic $pq$ and an e-geodesic $qr$ are said to be orthogonal if and only if

$$\sum_x(p(x)-q(x))(\ln q(x)-\ln r(x))=0.$$

Let $KL(*\|*)$ denote the Kullback-Leibler divergence defined by

For a distribution $p\in\mathcal P$ and a manifold $\mathcal M$, m-projection of $p$ onto $\mathcal M$ is defined by

$$m_p(p,\mathcal M):=\arg\min_{q\in\mathcal M}KL(p\|q).$$

If $\mathcal M$ is e-flat then the minimizer $q$ is uniquely determined and m-geodesic $pq$ is orthogonal to any e-geodesic $qr(r\in\mathcal M)$ as illustrated in Figure <ref> <cit.>.[e-flat manifolds are denoted by dotted lines in figures in this paper.]

m-projection onto e-flat manifold

Now, we are ready to discuss the information geometry of actual/pseudo-Gibbs sampling.
We define the following manifold:

\begin{equation}\label{eq:E}
E(\theta_i):=\{p\in\mathcal P|p(X_i|X_{-i})=\theta_i(X_i|Y_i)\}
\end{equation}

and refer to it as full-conditional-manifold of node $i$.
Any distribution $p\in\mathcal P$ is expressed as $p(X_{-i})p(X_i|X_{-i})$, and the manifold in Eq.(<ref>) determines the second factor that is the full-conditional distribution.
Then, the following theorem holds.

Any full-conditional-manifold is e-flat and m-flat.

The proof is given in Appendix.

Let $X^0,X^1...$ be an output sequence of pseudo-Gibbs sampling.
Here, we consider pseudo-Gibbs sampling as movements of $p^t(=p(X^t))$ in $\mathcal P$.

Suppose that node $i$ fires at a certain time $t$.
Then, $p^t$ moves to the following distribution:

\begin{equation}\label{eq:move to}
\end{equation}

that is, the full-conditional-distribution $p^t(X_i|X_{-i})$ is replaced by $\theta_i(X_i|Y_i)$.
As the following theorem shows, this $p^{t+1}$ is the m-projection of $p^t$ onto $E(\theta_i)$.

The m-projection of $p\in\mathcal P$ onto $E(\theta_i)$ is obtained by:

\begin{equation}\label{eq:m_p(p,E(theta_i))}
\end{equation}

Therefore, firing node $i$ moves the distribution $p$ into the m-projection of $p$ onto $E(\theta_i)$.

The proof is given in Appendix.

§.§ Ordered-(Actual/Pseudo)-Gibbs Sampling

As described in Section <ref>, in pseudo-Gibbs sampling, there are two options to select a node to fire.
The theoretical analysis for ordered-pseudo-Gibbs sampling is complicated because ordered-pseudo-Gibbs sampling is an inhomogeneous Markov chain.
However, ordered-pseudo-Gibbs sampling behavior is easily understood when represented graphically.

Movement of $p^t$ in ordered-actual-Gibbs sampling

Movement of $p^t$ in ordered-pseudo-Gibbs sampling

Figure <ref> illustrates the movements of $p^t$ in ordered-actual-Gibbs sampling.
In this case, the all full-conditional-manifolds intersect at a unique point $\pi$ because $\theta_i(X_i|Y_i)=\pi(X_i|X_{-i})$ for all $i$.
At first, node 0 fires, then the manifold $E(\theta_0)$ attracts $p^0$, and $p^0$ vertically falls down onto $E(\theta_0)$.
Next, node 1 fires, and $p^1$ falls down onto $E(\theta_1)$, and so on.
Given the depiction of the movement of $p^t$ in ordered-actual-Gibbs sampling shown in Figure <ref>, we can understand the convergence $p^t\to\pi$ intuitively.

On the other hand, Figure <ref> illustrates the case of ordered-pseudo-Gibbs sampling.
In this case, the full-conditional-manifolds do not have a common intersection.
However, if we consider a situation where every full-conditional-manifold is close to a certain point $\pi'$, as in Figure <ref>, $p^t$ is trapped in a small area around $\pi'$ after sufficient times of transitions.
To be more rigorous, $p^t$ moves along the cyclic orbit $\pi_0\to\pi_1\to...\to\pi_{n-1}\to\pi_0$, where $\pi_i$ is denoted in Eq.(<ref>).
As shown in Eq.(<ref>), $\pi$ is the centroid of $\{\pi_i\}(i=0,...,n-1)$.
Figure <ref> gives us an intuition that $\pi$ is close to $\pi'$.

§.§ Random-Pseudo-Gibbs Sampling

Comparing with ordered-pseudo-Gibbs sampling, we can perform more formal analysis for random-pseudo-Gibbs sampling because it is a homogeneous Markov chain and has a stationary distribution.
This formal analysis supports Hypothesis 1 “If every CPT $\theta_i(X_i|Y_i)$ is a good estimation of $p^*(X_i|X_{-i})$ then the stationary distribution of pseudo-Gibbs sampling becomes a good estimation of $p^*(X)$."

Under the condition that node $i$ fires at time $t$, $p^t$ moves to $p^{t+1}(X|i)=m_p(p^t,E(\theta_i))$ as shown in Theorem <ref>.
Therefore, under the condition that $i$ is unknown,

$$p^{t+1}(X)=\sum_i c_i p^{t+1}(X|i)=\sum_i c_i m_p(p^t,E(\theta_i)),$$

that is, a single step of random-pseudo-Gibbs sampling moves $p^t$ to $\sum_i c_i m_p(p^t,E(\theta_i))$.
Here, since $\pi$ is the stationary distribution,

\begin{equation}\label{eq:pi dont move}
\pi(X)=\sum_i c_im_p(\pi,E(\theta_i))(X)=\sum_i c_i\pi(X_{-i})\theta_i(X_i|Y_i)
\end{equation}

We here define KL-divergence between a distribution $p$ and a manifold $\mathcal M$ as

$$KL(p\|\mathcal M):=\min_{q\in\mathcal M}KL(p\|q)=KL(p\|m_p(p,\mathcal M))$$

and introduce the following pseudo-distance:

\begin{align}\label{eq:FC}
FC(p\|q)&:=\sum_i c_i KL(p\|E_i(q)),\quad E_i(q):=\{p\in\mathcal P|p(X_i|X_{-i})=q(X_i|X_{-i})\}\nonumber\\
&=\sum_i c_i\left<\ln\frac{p(X_i|X_{-i})}{q(X_i|X_{-i})}\right>_{p(X)}.
\end{align}

We refer to this pseudo-distance as full-conditional-divergence.

Figure <ref> illustrates the geometrical interpretation of full-conditional-divergence, which is the average of KL-divergence between $p$ and the manifolds $E_i(q)$.
$FC(p||q)$ takes a finite value if it satisfies the condition $\forall x, p(x)>0\Rightarrow q(x)>0$.

The relationship between full-conditional-divergence and Besag's pseudo-likelihood <cit.> resembles the relationship between KL-divergence and likelihood.
We can rewrite Eq.(<ref>) as

\begin{equation}\label{eq:FC rewrite}
FC(p\|q)=\left<\sum_ic_i\ln p(X_i|X_{-i})\right>_{p(X)}
-\underbrace{\left<\sum_ic_i\ln q(X_i|X_{-i})\right>_{p(X)}}
\end{equation}

and the KL-divergence is given as:

\begin{equation}\label{eq:KL}
KL(p\|q)=\left<\ln p(X)\right>_{p(X)}
-\underbrace{\left<\ln q(X)\right>_{p(X)}}_{\text{log-likelihood}}.
\end{equation}

In Eqs.(<ref>) and (<ref>), the first terms cancel the second terms when $p=q$.

Pseudo-(log-)likelihood was introduced by Besag as a convenient alternative to  (log-)likelihood to avoid the expensive computation of the partition function [Besag, 1975].
Note that minimizing KL-divergence is equivalent to maximizing likelihood; thus, minimizing full conditional divergence is equivalent to maximizing pseudo-likelihood.

Furthermore, KL-divergence and full-conditional-divergence are special cases of Bregman divergence <cit.>.
Here, let $\Lambda$ be a convex subset of $\mathbb R^k$.
For a strictly convex differentiable function $f:\Lambda\to\mathbb R$, which is called the Bregman function, the Bregman divergence is defined as:

$$B_f(p\|q):=f(p)-f(q)-\nabla f(q)\cdot(p-q),$$

where $\nabla f(q)$ denotes the gradient vector of $f$ at $q$, and $\cdot$ denotes the Euclidean inner product.
We can treat a distribution $p\in\mathcal P$ as a vector in $\mathbb R^{|X|}$ whose $x$-th component is $p(x)$.
The KL-divergence is a Bregman divergence whose Bregman function is the negative entropy:

$$f(p)=-H_p(X)=\left<\ln p(X)\right>_{p(X)},$$

and the full-conditional-divergence is a Bregman divergence whose Bregman function is the average of negative conditional entropies:

$$f(p)=-\sum_i c_i H_p(X_i|X_{-i})=\sum_i c_i \left<\ln p(X_i|X_{-i})\right>_{p(X)}.$$

A noticeable feature of full-conditional-divergence is that it depends on the dependency among variables.
KL-divergence treats variables as a single concatenated variable; therefore, it does not reflect any dependency among variables.

The following theorem that indicates where the stationary distribution of random-pseudo-Gibbs sampling exists.

Let $p\in\mathcal P$ be an arbitrary distribution and $\pi$ be the stationary distribution of random-pseudo-Gibbs sampling using CPTs $\{\theta_i\}$.
Then, the following inequality holds.

The proof is given in Appendix.
This theorem implies that if every full-conditional-manifold is close to a distribution $p$, then the stationary distribution $\pi$ is also close to $p$.
Thus, Hypothesis 1 is justified.

In the following, we focus on only random-pseudo-Gibbs sampling to simplify the theoretical analysis.

§ LEARNING ALGORITHM DERIVED BY INFORMATION GEOMETRY

Assume a dependency network user seeks to approximate a given target distribution $p\in\mathcal P$ by the stationary distribution of pseudo-Gibbs sampling.
According to Theorem <ref>, the objective is to locate every full-conditional-manifold close to $p$, that is, the goal is to minimize $KL(p\|E(\theta_i))$ for each $i$.
This objective yields a structure/parameter learning algorithm.

Ideally, the target distribution should be the true distribution $p^*$; however, it is unknown.
Therefore, we employ $\tilde p^*$—empirical distribution of training data—as the target distribution and then consider certain regularization techniques to avoid overfitting to $\tilde p^*$.
Then, learning becomes a task to find inputs $Y_i$ and parameters $\theta_i(X_i|Y_i)$ that minimizes $KL(\tilde p^*\|E(\theta_i))+R_i$, where $R_i$ denotes a regularization term.
The minimization problem to be solved becomes as follows:

$$\min_{Y_i\subseteq X_{-i}}\min_{\theta_i}KL(\tilde p^*\|E(\theta_i))+R_i.$$

Here, we can divide this minimization problem into the following two stages of minimization.

Parameter learning.Given $Y_i$, find $\hat\theta_i$ such that

\begin{equation}\label{eq:parameter learning}
\hat\theta_i=\arg\min_{\theta_i}KL(\tilde p^*\|E(\theta_i))+R_i.
\end{equation}

Structure learning.Find $\hat Y_i$ such that

\begin{equation}\label{eq:structure learning}
\hat Y_i=\arg\min_{Y_i\subseteq X_{-i}}KL(\tilde p^*\|E(\hat\theta_i))+R_i.
\end{equation}

We firstly determine inputs $Y_i$ by structure learning and then determine CPT by parameter learning.
For the convenience of explanation, we firstly discuss parameter learning and then discuss structure learning in this section.

§.§ Parameter Learning

We here prepare some quantities used in information theory.
For a distribution $p(X)$, the entropy $H_p(X)$ is defined by

$$H_p(X)=-\left<\ln p(X)\right>_{p(X)},$$

and for a distribution $p(XY)$, the conditional entropy $H_p(X|Y)$ is defined by

$$H_p(X|Y)=H_p(XY)-H_p(Y)=-\left<\ln p(X|Y)\right>_{p(XY)}.$$

Both quantities are always non-negative.

In Eq.(<ref>), we use a regularization term $R_i$ that does not depend on the parameters; therefore, we can ignore $R_i$ in parameter learning.
Here, the following theorem holds.

Let $\hat\theta_i(X_i|Y_i)$ be the CPT that minimizes $KL(p\|E(\theta_i(X_i|Y_i)))$ for given inputs $Y_i$.
Then, the following equations hold.

\begin{align}
&\hat\theta_i(X_i|Y_i)=p(X_i|Y_i)\label{eq:hat theta_i}\\
\end{align}

The proof is given in Appendix.

Equation (<ref>) directly provides a parameter learning algorithm.
Since $\tilde p^*$ is the empirical distribution of training data,

\begin{equation}\label{eq:counting data}
\hat\theta_i(x_i|y_i)=\frac{N_{x_iy_i}}{N_{y_i}},\quad
\end{equation}

where $N_{x_iy_i}$ is the number of occurrence of $X_i=x_i\wedge Y_i=y_i$ in the training data.
Then, the parameter learning algorithm becomes as follows:
Parameter learning of node $i$
$d$: training dataset

$Y_i$: inputs of node $i$

$\theta_i(X_i|Y_i)$: CPT

$N_{x_iy_i}\leftarrow\text{number of occurrences $X_i=x_i\wedge Y_i=y_i$ in $d$}$
$\theta_i(x_i|y_i)\leftarrow N_{x_iy_i}/N_{y_i}$

§.§.§ A Trick to Keep CPT Positive

If $p^*$ is known, the objective function to be minimized in learning $p^*(X_i|X_{-i})$ is

However, if CPT contains entries such that $\theta_i(x_i|y_i)=0$ then $KL(p^*\|E(\theta_i))$ diverges to infinity.
A simple trick to avoid this phenomenon is adding one to $N_{x_iy_i}$ such that $N_{x_iy_i}=0$.
Then, the parameter learning algorithm is modified as follows:
Parameter learning of node $i$ (ensuring CPT positive)
$d$: training dataset

$Y_i$: inputs of node $i$

$\theta_i$: CPT
[1]
$N_{x_iy_i}\leftarrow\text{number of occurrences $X_i=x_i\wedge Y_i=y_i$ in $d$}$

$\theta_i(x_i|y_i)\leftarrow N_{x_iy_i}/N_{y_i}$

The for-loop between line <ref> and line <ref> is the trick to ensure $N_{x_iy_i}>0$.

This trick provides another benefit.
It guarantees the ergodicity of pseudo-Gibbs sampling because, for any two values $x^0,x^1$, pseudo-Gibbs sampling is able to move from the state $X=x^0$ to the state $X=x^1$ within $n$ transitions.

§.§ Structure Learning

As shown in Eq.(<ref>), structure learning is the task to find $\hat Y_i\subseteq X_{-i}$ such that

\begin{align*}
\hat Y_i&=\arg\min_{Y_i\subseteq X_{-i}}KL(\tilde p^*\|E(\hat\theta_i))+R_i\\
&=\arg\min_{Y_i\subseteq X_{-i}}H_{\tilde p^*}(X_i|Y_i)-H_{\tilde p^*}(X_i|X_{-i})+R_i\quad\text{(by Eq.\eqref{eq:KL(p||E)})}\\
&=\arg\min_{Y_i\subseteq X_{-i}}H_{\tilde p^*}(X_i|Y_i)+R_i.
\end{align*}

Therefore, the cost function to be minimized in structure learning becomes as follows:

\begin{equation}\label{eq:scost_i}
scost_i(Y_i)=H_{\tilde p^*}(X_i|Y_i)+R_i.
\end{equation}

Here, $H_{\tilde p^*}(X_i|Y_i)$ is obtained by

$$H_{\tilde p^*}(X_i|Y_i)=-\sum_{x_iy_i}\frac{N_{x_iy_i}}N\ln\frac{N_{x_iy_i}}{N_{y_i}},
\quad N_{y_i}=\sum_{x_i}N_{x_iy_i}.$$

Various information criteria such as AIC [Akaike, 1974] and MDL [Rissanen, 1989] are applicable for the regularization term:

\begin{align}
R_i&=\frac{k_i}{2N}\ln N\quad\text{(MDL)},\label{eq:MDL}
\end{align}

where $k_i=(|X_i|-1)|Y_i|$ is the degrees of freedom that CPT $\theta_i(X_i|Y_i)$ has, and $|X_i|,|Y_i|$ are the number of values $X_i,Y_i$ possibly takes.

It is practical to use a greedy search to minimize $scost_i$.
The following algorithm is an example of greedy search for structure learning.
Structure learning of node $i$
$Y_i$: inputs

$scost_i$: structure cost function (Eq.(<ref>))
[1]
Evaluate $scost_i(Y'_i)$ for all possible candidates $Y'_i$, and find
$scost_i(Y^{min})\ge scost_i(Y_i)$$Y_i$
$Y_i\leftarrow Y^{min}$

“Candidate" in line <ref> is a set of input variables obtained by adding a single node into $Y_i$ or removing a single node from $Y_i$.

Other than this greedy algorithm, tabu search <cit.> is also practical to search $\hat Y_i$.

§.§ Convergence to True Distribution

A significant interest is determining whether stationary distribution $\pi$ converges to the true distribution $p^*$ at the limit of $N\to\infty$, where $N$ is the number of samples in the training data.

Here, let $d^\infty=\{x^0,x^1,...\}$ be a set of infinite training data, $d^N$ be the first $N$ samples of $d^\infty$ and $\tilde p^*(X)$ be the empirical distribution of $d^N$.
Given $d^N$, the structure/parameter learning algorithm of the node $i$ determines its inputs $Y_i$ and the parameters $\theta_i$; therefore $Y_i$ and $\theta_i$ are the functions of $N$.
Here, we consider the regularization term $R_i$ as a function $R_i(Y_i,N)$.
For simplicity, we assume that Algorithm <ref> actually finds the minimizer of Eq.(<ref>) here.
The following theorem demonstrates that the full-conditional-manifolds approach $\tilde p^*$ as $N$ increases.

If $\lim_{N\to\infty}R_i(X_{-i},N)=0$ then

\begin{equation}\label{eq:convergence KL(p||E)}
\lim_{N\to\infty}KL(\tilde p^*\|E(\hat\theta_i))=0.
\end{equation}

The proof is given in Appendix.
It is clear that the regularization terms in Eqs.(<ref>) and (<ref>) satisfy the condition $\lim_{N\to\infty}R_i(X_{-i},N)=0$.

Note that the greedy search algorithm, such as Algorithm <ref>, may not find the global minimum of $scost_i$.
However, we can guarantee the convergence in Eq.(<ref>) by employing a trick to such cases.
If we add the following step as the last step of any search algorithm then the convergence is guaranteed.

Last step) If $scost_i(Y_i,N)>scost_i(X_{-i},N)$ then $Y_i\leftarrow X_{-i}$ and output $Y_i$.

However, this trick has less practical meaning because $scost_i(X_{-i},N)\gg scost_i(Y_i,N)$ in most practical cases.

According to Theorem <ref> and Theorem <ref>, it is clear that $\lim_{N\to\infty}FC(\tilde p^*\|\pi)=0$, and this result directly deduces the following theorem.

If $\lim_{N\to\infty}R_i(X_{-i},N)=0$ and $\lim_{N\to\infty}\tilde p^*=p^*$ then

§ INFERENCE BY PSEUDO-GIBBS SAMPLING

According to Theorem <ref>, if $KL(p\|E(\theta_i))$ is small for each $i$, then the stationary distribution of pseudo-Gibbs sampling becomes a good estimation of $p(X)$.
As described in Section <ref>, clamped-pseudo-Gibbs sampling is equivalent to the free-pseudo-Gibbs sampling with CPTs ${\theta_i(X_i|Z_iv)}$ (see Section <ref> for the meaning of the symbols).
Then if $KL(p(\tilde V|v)\|\theta_i(X_i|Z_iv))$ is small for each $i$, the stationary distribution of the clamped-pseudo-Gibbs sampling becomes a good estimation of $p(\tilde V|v)$.

Consider the following equation:

\begin{align*}
&=\sum_{x_iz_iv} p(x_iz_iv)\ln\frac{p(x_i|z_iv)}{\theta_i(x_i|z_iv)}\\
&=\sum_v p(v)\sum_{x_iz_i}p(x_iz_i|v)\ln\frac{p(x_i|z_iv)}{\theta_i(x_i|z_iv)}\\
&=\sum_v p(v)\left<\ln\frac{p(X_i|Z_iv)}{\theta_i(X_i|Z_iv)}\right>_{p(\tilde V|v)}\\
&=\sum_v p(v)KL(p(\tilde V|v)\|E(\theta_i(X_i|Z_iv)).
\end{align*}

This equation means that $KL(p\|E(\theta_i))$ is the average of $KL(p(\tilde V|v)\|E(\theta_i(X_i|Z_iv))$.
In other words, if $KL(p\|E(\theta_i))$ is small then $KL(p(\tilde V|v)\|E(\theta_i(X_i|Z_iv)))$ is small on average.
Note that $KL(p(\tilde V|v)\|E(\theta_i(X_i|Z_iv))$ may become large if $p(v)$ is small, which implies that the inference for rare conditions becomes difficult.

§ GENERALIZED INPUT FOR CPT

In the sections so far, we let the inputs $Y_i$ be certain variables in $X_{-i}$.
We refer to this type of inputs as the direct input.
Here, we notice it is unnecessary to let $Y_i$ be variables in $X_{-i}$ because $Y_i$ is just a (concatenated) variable to be referred by CPT.
Even if we consider $Y_i$ is a variable determined by $Y_i=f_i(X_{-i})$, where $f_i$ is a function of $X_{-i}$, the most of discussion made in the sections so far remains valid.
Let us refer to this type of input as the generalized input.
The direct input is a special case of the generalized input where $f_i$ is a function to select variables in $X_{-i}$.

In the case we use the generalized input, the role of structure learning is to find a function $f_i$ that minimizes the structure cost in Eq.(<ref>).
Note that $|f_i|$, which is the size of $f_i$'s range, is suppressed to suppress the CPT's degrees of freedom and the regularization term.

Providing a function $f_i$ is equivalent to clustering all possible values $x_{-i}$ into $|f_i|$ (size of $f_i$'s range) groups.
A practical method to find $f_i$ in structure learning is using a decision tree learning algorithm [Heckerman et al., 2000].
In this case, a leaf of the decision tree corresponds to a group that shares a value of $f_i$.

§ EXPERIMENTS

In this section, we demonstrate a dependency network's performance by comparing it to a Bayesian network.

§.§ Compared Learning Machines

We compared two learning machines: a dependency network (DN) and a Bayesian network (BN).

§.§.§ Dependency Network DN

We used Algorithm <ref> for structure learning and Algorithm <ref> for parameter learning.
In structure learning, we used MDL in Eq.(<ref>) for the regularization term.
For inference, we used Algorithm <ref> and used the setting $b=k=n$(number of nodes) for burn-in and thinning.

§.§.§ Bayesian Network BN

Similarly to the learning discussed in Section <ref>, learning in BN comes down to the following minimization problems.

Parameter learning.Given graphical structure $G$, find $\hat\pi$ such that

$$\hat\pi=\arg\min_\pi KL(\tilde p^*\|\pi)+R.$$

Structure learning.Find $\hat G$ such that

$$\hat G=\arg\min_G KL(\tilde p^*|\hat\pi)+R.$$

We here skip the derivation and only show the following consequences.

\begin{equation}\label{eq:BN parameter learning}
\hat\pi(X_i|Y_i)=\tilde p^*(X_i|Y_i)
\end{equation}

$$KL(\tilde p^*\|\hat\pi)=-H_{\tilde p^*}(X)+\sum_i H_{\tilde p^*}(X_i|Y_i)$$

\begin{equation}\label{eq:BN scost}
scost(G)=\sum_i H_{\tilde p^*}(X_i|Y_i)+R
\end{equation}

In the experiments, we used the following MDL based regularization term.

$$R=\sum_i\frac{k_i}{2N}\ln N,\quad k_i=(|X_i|-1)|Y_i|,$$

where $k_i$ is the degrees of freedom that CPT of node $i$ has (see Eq.(<ref>)).

We used the following algorithm for structure and parameter learning.
Structure and parameter learning in BN
$G$: graph of BN

$N$: number of samples in training data

$scost$: structure cost function (Eq.(<ref>))
[1]
Evaluate $scost(G')$ for all possible candidates $G'$, and find

$scost(G^{min})\ge scost(G)$
Determine CPT by Eq.(<ref>)
$G$ and CPT
$G\leftarrow G^{min}$

“Candidate" in line <ref> is a graph obtained by one of the following modifications to $G$.

* Adding a single edge to $G$ unless the candidate has any cycles.

* Removing a single edge from $G$.

* Reversing the direction of a single edge in $G$ unless the candidate has any cycles.

We drew output data by the following algorithm.
Drawing output data from BN
$G$: graph of BN

$H$: a graph

$o$: output data

$a_*$: sequence of node number

$n$: number of nodes

$N$: number of samples in output data
[1]
$H\leftarrow G$, $o\leftarrow\emptyset$
$i=0$ to $n-1$
Select a node $j$ such that has no parents in $H$
$a_i\leftarrow j$;
Remove node $j$ and its outgoing edges from $H$

$t=0$ to $N-1$
$i=0$ to $n-1$
Fire node $a_i$
Append $X$ into $o$

The for-loop between line <ref> and line <ref> is a topological sorting <cit.>.

The output data drawn by this algorithm are i.i.d. data; thus, neither burn-in nor thinning techniques are needed.

§.§ Training Datasets

We used the following four training datasets for the experiments.
Each dataset is an i.i.d. dataset drawn from a known true distribution $p^*$.

§.§.§ BN20-37S, BN20-37L

BN20-37S and BN20-37L are datasets drawn from a Bayesian network with 20 nodes and 37 edges (Figure <ref>).
Every $X_i$ takes value 0 or 1.
BN20-37S includes 1000 samples, and BN20-37L has 100000 samples.

§.§.§ Ising5x5S, Ising5x5L

Ising5x5S and Ising5x5L are datasets drawn from a Ising spin model that has $5\times5$ nodes (Figure <ref>).
Every $X_i$ takes the value 0 or 1.
Ising5x5S has 1000 samples, and Ising5x5L has 100000 samples.

§.§ Platform for Experiments

All experiments were performed on a laptop PC(CPU: Intel Core i7-6700K @4GHz; Memory: 64GB; OS: Windows 10 Pro).
All programs were written in and executed on Java 8.

§.§ Results

Performance comparison between DN and BN
Dataset    $n$    $N$    $KL(\tilde p^*\|p^*)$    $KL(\tilde\pi\|p^*)$(DN)   $KL(\tilde\pi\|p^*)$(BN)
1.43	   0.78

BN20-37S	   20	   1000		   4.14	   1.43	   0.78

1.43	   0.78
0.69	   0.69

BN20-37L	   20	   100000	   0.67	   0.69	   0.69

0,69	   0.70
1.37	   1.37

Ising5x5S	   25	   1000		   3.87	   1.36	   1,37

1.39	   1.37
1.17	   1.16

Ising5x5L	   25	   100000	   1.16	   1.18	   1.15

1.19	   1.16

$n$: number of variables.
$N$: number of samples.
$\tilde p^*$: empirical distribution of training data.
$p^*$: true distribution. $\tilde\pi$: empirical distribution of output data.
The unit for KL-divergence is “nat".

Table <ref> demonstrates the performance comparison between DN and BN.
After learning, we drew 1000000 samples of output data respectively from DN and BN.
We measured the accuracy of learning by $KL(\tilde\pi\|p^*)$, where $\tilde\pi$ is the empirical distribution of output data, and $p^*$ is the true distribution of the training dataset.
The outputs depend on the random seeds used to draw outputs; therefore, we made three trials with different random seeds for each dataset.

For BN20-27S, BN has better performance than DN has.
One possible reason for this result is that BN20-37S is a dataset drawn from a  Bayesian network, so BN easily learns $p^*$ while DN can not imitate $p^*$ by a simple structure.
For the larger samples (BN20-37L), DN is allowed to take a more complex structure and learn $p^*$ correctly.

For the datasets BN20-37L, Ising5x5S, and Ising5x5L, $KL(\tilde\pi\|p^*)(DN)\approx KL(\tilde\pi\|p^*)(BN)$, that is, DN and BN have roughly the same performance.

For the datasets BN20-37L and Ising5x5L, both the training data and the output data have the same number of samples (100000); therefore, we can compare $KL(\tilde\pi\|p^*)$ and $KL(\tilde p^*\|p^*)$ directly.
The fact $KL(\tilde\pi\|p^*)\approx KL(\tilde p^*\|p^*)$ implies that the output data drawn from $\pi$ and the training data—i.i.d. data drawn from $p^*$—are indistinguishable; thus, the learning machine correctly learns $p^*$.

Learning behavior of node 0 in DN
Dataset	   $Y_0$	   $H(X_0)$	   $H(X_0|Y_0)$	   $KL(\tilde p^*\|E(\theta_0))$    $KL(p^*\|E(\theta_0))$
BN20-37S     $X_1X_2X_4X_{19}$								   0.68	   0.41	   0.41	   0.17

BN20-37L     $X_1X_2X_3X_4X_6X_7X_{11}X_{19}$	   0.68	   0.22	   0.09	   1.3E-3

Ising5x5S	   $X_1X_5$													   0.69	   0.44	   0.36	   2.9E-3

Ising5x5L    $X_1X_5$													   0.69	   0.44	   0.19	   1.4E-5

The unit for entropy and KL-divergence is “nat".

Here, we examine more details of the learning behavior of DN.
Table <ref> focuses on node 0.
We explain the results for BN20-37S as an example here.
The learning algorithm starts with no inputs ($Y_0=\emptyset$), and $X_0$ has 0.68nat of entropy.
The structure learning algorithm takes $X_1X_2X_4X_{19}$ as the inputs, and the entropy of $X_0$ is reduced to 0.41nat.
For the full conditional manifold $E(\theta_0)$ constructed by the learning algorithm, $KL(\tilde p^*\|E(\theta_0))=0.41$nat and $KL(p^*\|E(\theta_0))=0.17$nat.
The fact $KL(\tilde p^*\|E(\theta_i))\gg KL(p^*\|E(\theta_i))$ implies that overfitting to $\tilde p^*$ is successfully avoided, and the learning algorithm has high generalization ability.
We confirmed that the inequality $KL(\tilde p^*\|E(\theta_i))\gg KL(p^*\|E(\theta_i)))$ holds for every node.

Average of $KL(\tilde p^*\|E(\theta_i))$ and $KL(p^*\|E(\theta_i))$
Dataset    $\overline{KL(\tilde p^*\|E(\theta_0))}$    $\overline{KL(p^*\|E(\theta_0))}$
BN20-37S	   0.47	   4.5E-2

BN20-37L	   0.17	   1.1E-3

Ising5x5S    0.28	   7.0E-3

Ising5x5L    0.16	   5.3E-5

The unit for entropy and KL-divergence is “nat".

Table <ref> demonstrates the average of $KL(\tilde p^*\|E(\theta_i))$ and $KL(p^*\|E(\theta_i))$ defined by:

$$\overline{KL(\tilde p^*\|E(\theta_i))}=\sum_i c_i KL(\tilde p^*\|E(\theta_i)),\quad
\overline{KL(p^*\|E(\theta_i))}=\sum_i c_i KL(p^*\|E(\theta_i)).$$

As shown in the table, $\overline{KL(p^*\|E(\theta_i))}$ takes small values, and this fact implies that the stationary distribution $\pi$ is close to the true distribution $p^*$ because the inequality $FC(p^*\|\pi)\le\overline{KL(p^*\|E(\theta_i))}$ holds by Theorem <ref>.

Number of structure cost evaluations and CPU time for learning

Dataset		   $n$	   $N$		   #eval(DN)	   #eval(BN)	   time(DN)	   time(BN)
BN20-37S	   20		   1000		   1616				   16945 			   17				   73

BN20-37L	   20		   100000    2129				   23797			   2466			   13127

Ising5x5S	   25		   1000		   2497				   32298			   25				   144

Ising5x5L	   25		   100000    2545				   42365			   2419			   22676

$n$: number of nodes. $N$: number of samples in training data.

#eval: number of structure cost evaluations. time: CPU time to learn (ms).

Table <ref> demonstrates the number of structure cost evaluations and CPU time to learn the datasets.
Here, #eval(DN) is the number of evaluations in Algorithm <ref> line <ref>, and #eval(BN) is the number of evaluations in Algorithm <ref> line <ref>.
To measure CPU time, we made three trials for each dataset and showed the median in the table.

In the learning algorithms used for DN and BN, the structure cost evaluations dominate the computational costs.
The computational cost to evaluate one structure obeys $O(N)$; thus, the time to learn obeys $O(N\times\text{\#eval})$.
As shown in the table, DN learns the datasets with much smaller number of evaluations than BN and learns faster than BN.

§ CONCLUSION AND FUTURE WORK

In dependency network learning, the network does not directly learn the true distribution $p^*(X)$ underlying training data.
Rather than learning the joint distribution $p^*(X)$, each node independently learns the full conditional distribution $p^*(X_i|X_{-i})$.
This property—independence of learning tasks—keeps the dependency network learning tractable even where the number of nodes is large.

Pseudo-Gibbs sampling gathers the CPTs owned by the nodes and constructs a joint distribution as its stationary distribution.
If every node learns $p^*(X_i|X_{-i})$ successfully, that is, every CPT is a good estimation of $p^*(X_i|X_{-i})$, then the stationary distribution of pseudo-Gibbs sampling becomes a good estimation of $p^*(X)$.
It is challenging to know the stationary distribution numerically, however, we can draw as many samples as we want from the stationary distribution using pseudo-Gibbs sampling; thus, we can estimate the probabilities and expectation of variables using Monte Carlo methods.

From the information geometry perspective, we consider a manifold $E(\theta_i)$ for each node determined by its CPT.
Then, pseudo-Gibbs sampling is interpreted as iterative m-projections onto these manifolds.
This interpretation tells us how we should make the stationary distribution of pseudo-Gibbs sampling close to the desired distribution and simultaneously provides the structure/parameter learning algorithms as optimization problems.

We introduced a pseudo-distance $FC(p\|q)$ that is a Bregman divergence whose Bregman function is Besag's pseudo-log-likelihood.
We built a hypothesis “If every CPT is a good estimation of $p^*(X_i|X_{-i})$ then the stationary distribution of pseudo-Gibbs sampling is close to $p^*(X)$."
Using the pseudo-distance $FC$, we justified this hypothesis (Theorem <ref>).

We compared a dependency network (DN) and a Bayesian network (BN) experimentally.
Regarding the accuracy of the learned distribution, there was no significant difference between DN and BN.
The experiments also demonstrated that DN was much faster than BN in learning.

One possible future work is performance evaluation by real problems.
Suppose we perform a hold-out test for performance evaluation.
For real problems, the true distribution $p^*$ is unknown, and only $p^{test}$ that is an empirical distribution of data drawn from $p^*$ is known.
In the case of a dependency network, the learned distribution $\pi$ is also unknown, and only $\tilde\pi$, which is the empirical distribution of data drawn from $\pi$, is given.
Furthermore, both $p^{test}$ and $\tilde\pi$ are very “sparse", that is, $p^{test}(x)=0$ and $\tilde\pi(x)=0$ for most values $x$.
Therefore, we have to compare two sparse distributions $p^{test}$ and $\tilde\pi$.
However, evaluating the closeness between such two sparse distribution is a challenging problem.

Another possible future work is introducing hidden variables into dependency networks.
We expect that introducing hidden variables will provide much higher performance into dependency networks.

This work was supported by JSPS KAKENHI 17H01793.

§ APPENDIX

§.§ Proof of Theorem <ref>

Let $p_0,p_1$ be two distributions on the manifold $E(\theta_i)$, that is,

\begin{align*}
\end{align*}

and $p_\lambda$ be a distribution on the e-geodesic between $p_0$ and $p_1$ expressed by

\begin{align*}
\ln p_\lambda(X)&=(1-\lambda)\ln p_0(X)+\lambda p_1(X)-\ln Z\\
&=\ln(1-\lambda)p_0(X_{-i})+\lambda p_1(X_{-i})+\ln\theta_i(X_i|Y_i)-\ln Z.
\end{align*}

Then, we obtain:

$$p_\lambda(X)=\frac{\theta_i(X_i|Y_i)}Z e^{(1-\lambda)p_0(X_{-i})
+\lambda p_1(X_{-i})}$$

\begin{align*}
&=\frac1Z e^{(1-\lambda)p_0(X_{-i})+\lambda p_1(X_{-i})}\\
\end{align*}

Therefore, $E(\theta_i(X_i|Y_i))$ is e-flat.

Similarly, let $q_\lambda$ be a distribution on the m-geodesic between $p_0$ and $p_1$ expressed by

\begin{align*}
&=(1-\lambda)p_0(X)+\lambda p_1(X)\\
&=((1-\lambda)p_0(X_{-i})+\lambda p_1(X_{-i}))\theta_i(X_i|Y_i).
\end{align*}

Then, we obtain:

$$q_\lambda(X_{-i})=(1-\lambda)p_0(X_{-i})+\lambda p_1(X_{-i})$$

Therefore, $E(\theta_i)$ is m-flat.

§.§ Proof of Theorem <ref>

Let $q(X)=q(X_{-i})\theta_i(X_i|Y_i)$ be a distribution on $E(\theta_i)$.
Then, we obtain:

\begin{align*}
\end{align*}

It is clear that $KL(p\|q)$ is minimized when $q(X_{-i})=p(X_{-i})$; thus, we obtain:

$$m_p(p,E(\theta_i))=\arg\min_{q\in E(\theta_i)}KL(p\|q)=p(X_{-i})\theta_i(X_i|Y_i).$$

This result coincides with Eq.(<ref>).

§.§ Proof of Theorem <ref>

Let $\pi_i$ denote the m-projection of $\pi$ onto $E(\theta_i)$, that is,
Then, we obtain:

\begin{align*}
\sum_i c_i KL(p\|E(\theta_i))-FC(p\|\pi)
&=\sum_i c_i\left<\ln\frac{p(X_i|X_{-i})}{\theta_i(X_i|Y_i)}\right>_{p(X)}
-\sum_i c_i\left<\ln\frac{p(X_i|X_{-i})}{\pi(X_i|X_{-i})}\right>_{p(X)}\\
&=\sum_i c_i\left<\ln\frac{\pi(X_i|X_{-i})}{\theta_i(X_i|Y_i)}\right>_{p(X)}\\
&=\sum_i c_i\left<\ln\frac{\pi(X_{-i})\pi(X_i|X_{-i})}{\pi(X_{-i})\theta_i(X_i|Y_i)}\right>_{p(X)}\\
&=\sum_i c_i\left<\ln\frac{\pi(X)}{\pi_i(X)}\right>_{p(X)}\\
&=\left<-\sum_i c_i\ln\frac{\pi_i(X)}{\pi(X)}\right>_{p(X)}\\
&\ge\left<-\ln\sum_i c_i\frac{\pi_i(X)}{\pi(X)}\right>_{p(X)}\quad
\text{$\because-\ln$ is a convex function}\\
&=\left<-\ln\frac{\pi(X)}{\pi(X)}\right>_{p(X)}\quad\text{by Eq.\eqref{eq:pi dont move}}\\
\end{align*}

§.§ Proof of Theorem <ref>

\begin{align*}
\end{align*}

It is clear that this equation is minimized when $\theta_i(X_i|Y_i)=p(X_i|Y_i)$, and then

§.§ Proof of Theorem <ref>

Since $\hat Y_i$ is the minimizer of the cost function $scost$ in Eq.(<ref>), $scost_i(\hat Y_i,N)\le scost_i(X_{-i},N)$, that is,

$$H_{\tilde p^*}(X_i|\hat Y_i)+R_i(\hat Y_i,N)\le H_{\tilde p^*}(X_i|X_{-i})+R_i(X_{-i},N).$$

Therefore, we obtain:

\begin{equation}\label{eq:bound KL}
\begin{split}
KL(\tilde p^*\|E(\hat\theta_i))&=H_{\tilde p^*}(X_i|\hat Y_i)-H_{\tilde p^*}(X_i|X_{-i})\quad\text{by Eq.\eqref{eq:KL(p||E)}}\\
&\le R_i(X_{-i},N)-R_i(\hat Y_i,N)\\
&\le R_i(X_{-i},N).
\end{split}
\end{equation}

Considering Eq.(<ref>) at the limit of $N\to\infty$, we obtain:

$$0\le\lim_{N\to\infty}KL(\tilde p^*\|E(\hat\theta_i))\le\lim_{N\to\infty}R_i(X_{-i},N)=0.$$

[Akaike, 1974]
Hirotugu Akaike.
A new look at the statistical model identification.
IEEE transactions on automatic control, 190
(6):0 716–723, 1974.

[Amari, 2016]
Shun-ichi Amari.
Information geometry and its applications.
Springer, 2016.

[Besag, 1975]
Julian Besag.
Statistical analysis of non-lattice data.
Journal of the Royal Statistical Society: Series D (The
Statistician), 240 (3):0 179–195, 1975.

[Censor and Zenious, 1997]
Yair Censor and Stravos A. Zenious.
Parallel Optimization: theory, algorithms and applications.
Oxford University Press, 1997.

[Gilks et al., 1995]
Walter R Gilks, Sylvia Richardson, and David Spiegelhalter.
Markov chain Monte Carlo in practice.
Chapman and Hall/CRC, 1995.

[Glover and Taillard, 1993]
Fred Glover and Eric Taillard.
A user's guide to tabu search.
Annals of operations research, 410 (1):0
1–28, 1993.

[Heckerman et al., 2000]
David Heckerman, David Maxwell Chickering, Christopher Meek, Robert
Rounthwaite, and Carl Kadie.
Dependency networks for inference, collaborative filtering, and data
Journal of Machine Learning Research, 10
(Oct):0 49–75, 2000.

[Knuth, 1997]
Donald Ervin Knuth.
The art of computer programming, volume 1.
Addison-Wesley, 1997.

[Koller and Friedman, 2009]
Daphne Koller and Nir Friedman.
Probabilistic graphical models: principles and techniques.
MIT press, 2009.

[Liu and Nocedal, 1989]
Dong C Liu and Jorge Nocedal.
On the limited memory bfgs method for large scale optimization.
Mathematical programming, 450 (1-3):0
503–528, 1989.

[Rissanen, 1989]
Jorma Rissanen.
Stochastic complexity in statistical inquiry.
World Scientific, 1989.

[Seneta, 2006]
Eugene Seneta.
Non-negative matrices and Markov chains.
Springer Science & Business Media, 2006.

[Takabatake and Akaho, 2012]
Kazuya Takabatake and Shotaro Akaho.
Constraint-free graphical model with fast learning algorithm.
arXiv preprint arXiv:1206.3721, 2012.