論文 / 著書情報
Article / Book Information
題目(和文)

パラメータ制約付き特異モデルの統計的学習理論

Title(English)

Statistical Learning Theory of Parameter-Restricted Singular Models

著者(和文)

林 直輝

Author(English)

Naoki Hayashi

出典(和文)

学位:博士(理学),学位授与機関:東京工業大学,
報告番号:甲第12028号,
授与年月日:2021年6月30日,
学位の種別:課程博士,
審査員:渡邊 澄夫,高安 美佐子,金森 敬文,山下 真,澄田 範奈

Citation(English)

Degree:Doctor (Science),
Conferring organization: Tokyo Institute of Technology,
Report number:甲第12028号,
Conferred date:2021/6/30,
Degree Type:Course doctor,
Examiner:,,,,

学位種別(和文)

博士論文

Type(English)

Doctoral Thesis

Powered by T2R2 (Tokyo Institute Research Repository)

2021FY Ph.D. Dissertation（令和 3 年度 博士論文）

Statistical Learning Theory of Parameter-Restricted
Singular Models（パラメータ制約付き特異モデルの統計的学習理論）

First Submit
ReviseFinal Submit
Supervisor
Ch. Examiner
Examiners

Aﬃliation

2021/02/17
2021/03/21, 2021/04/16
2021/05/20
Prof. Sumio Watanabe
Prof. Sumio Watanabe
Prof. Takafumi Kanamori
Prof. Misako Takayasu
Prof. Makoto Yamashita
Prof. Hanna Sumita
Tokyo Institute of Technology
School of Computing
Department of Mathematical and Computing Science

Naoki Hayashi

ii

Abstract
Statistical models used in machine learning are called learning machines. It is well-known thatlearning machines are widely applied to predict unknown events and discover knowledge bycomputers in many fields. Indeed, machine learning has grown over the last several decades.
They are used for statistical learning/inference and usually have hierarchical structures. These
structures are eﬀective for generalizing to the real world. Statistical learning theory is a theory
to clarify the generalization performances of learning machines.
Singular learning theory is a mathematical foundation for statistical inference using singularmodels. Typical hierarchical models, such as neural networks, tree and forest model, mixture
model, matrix factorization, and topic model, are statistically singular since a map froma parameter to a probability density function is not one-to-one. Clarifying generalization
behaviors in singular models is an important problem to estimate suﬃcient sample sizes,
design models, and tune hyperparameters. However, conventional statistics theory cannot
be applied to these models because their likelihoods cannot be approximated by any normal
distribution. Singular learning theory provides a general view for this problem; birational
invariants of an analytic set (a.k.a. algebraic variety) determine the generalization error.
That is defined by zero of a Kullback-Leibler (KL) divergence between the data-generating
distribution and the model. Algebraic structures of statistical models are essential in singularlearning theory; thus, it can be interpreted as an intersection between algebraic statistics andstatistical learning theory.
One of such invariants is a real log canonical threshold (RLCT). An RLCT is a negativemaximum pole of a zeta function defined by an integral of a KL divergence. Determining
an RLCT of a concrete model is performed by resolution of singularities. In fact, algebraic
statisticians and machine learning researchers have derived the exact values or upper boundsof the RLCTs for several singular models. The theoretical value of the RLCT is eﬀective instatistical model selection such as sBIC proposed by Drton and Plummer. Besides, Nagata
proposed a tuning method using RLCTs for exchange Monte Carlo.
On the other hand, from the practical point of view, the parameter region of the model
is often restricted to improve interpretability. Non-negative matrix factorization (NMF) andlatent Dirichlet allocation (LDA) are well-known examples of parameter-restricted singularmodels. In general, such constraints make the generalization error changed. However, foreach singular model and condition, the quantitative eﬀect of those constraints has not yet beenclarified because the singularities in the above analytic set are also changed by the restriction
to the parameter region.
In this dissertation, as a foundation to establish a singular learning theory of parameterrestricted statistical models, we theoretically study the asymptotic behavior of the Bayesian
generalization error in NMF and LDA. NMF and LDA are two typical singular models whoseparameter regions are constrained. In NMF, we derive an upper bound of the RLCT and alower bound of the variational approximation error. In LDA, we prove that its RLCT is equal
to that of matrix factorization with simplex restriction and clarify the exact asymptotic formof the generalization error, i.e. we determine the exact value of the RLCT of LDA. These
results provide quantitative diﬀerences of generalization errors from matrix factorization
whose parameter space is not restricted.

iii

Contents
Chapter 1

1.1
1.2
1.3
1.4
Chapter 2

2.1
2.2
Chapter 3

3.1
3.2
3.3
3.4
Chapter 4

4.1
4.2
4.3
4.4
4.5
4.6
Chapter 5

5.1
5.2
5.3
5.4
5.5
5.6
Chapter 6

6.1
6.2
Appendix A

Introduction

1

Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Research Goal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Dissertation Structure . . . . . . . . . . . . . . . . . . . . . . . . . .
Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
5
5
6

Bayesian Inference
7
Basic Concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
Framework of Bayesian Inference . . . . . . . . . . . . . . . . . . . . 10
Singular Learning Theory
Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Singularity Resolution Theorem and Zeta Function . . . . . . . . . .
Key Results of Singular Learning Theory . . . . . . . . . . . . . . .
Information Criteria from Singular Learning Theory . . . . . . . . . .

15

Bayesian Generalization Error in Non-negative Matrix Factorization
Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Main Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Preparation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Proof of Main Theorem . . . . . . . . . . . . . . . . . . . . . . . . .
Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

Bayesian Generalization Error in Latent Dirichlet Allocation
Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Main Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Preparation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Proof of Main Theorem . . . . . . . . . . . . . . . . . . . . . . . . .
Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61

15
17
21
24

31
33
35
47
49
58

61
64
67
70
76
83

Conclusion
85
Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
Questions and Answers in Defense

89

Contents

iv
A.1
A.2

Singular Learning Theory . . . . . . . . . . . . . . . . . . . . . . . . 89
Main Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

Acknowledgement

95

Acknowledgement (in Japanese)

97

Bibliography

99

List of Publications

105

v

List of Figures
1.1

Horizontal bar plots of regression coeﬃcients. Left graph is the nonrestricted result and right one is the result with non-negative restriction. . . .

4

4.1This image represents the critical lines of the RLCT of NMF and the learning
coeﬃcient in VBNMF when M = N . The horizontal and the vertical axises
are corresponding to the hyperparameter ϕU and ϕV , respectively. The usedfour points {(0.25, 0.25), (0.5, 0.5), (1, 1), (2, 2)} are plotted as black disks.
The red-dashed-line is the critical line of the RLCT of NMF by Theorem
4.2. On the other hand, the blue-line is that of the learning coeﬃcient in
VBNMF by Theorem 4.1. Obviously, they are thoroughly diﬀerent: VB is
not equivalent to Bayesian inference. . . . . . . . . . . . . . . . . . . . . . 59

5.1This figure gives an overview of LDA. The categorical distributions Cat
that depend on the documents. Words in the uppercase such as NAME,
FOOD, and CODING are topics. There are categorical distributions that are
diﬀerent for each topic; the words (Ciel, curry, lambda, ...) are generated
from them. Hence, we can explain LDA as a mixture of categorical mixture
models. This model has a hierarchical structure. This figure is quoted andmodified from the author’s works [38, 34] . . . . . . . . . . . . . . . . . . 62(a) In this chapter, we give the exact value of the learning coeﬃcient of LDA
λ. The learning coeﬃcient is smaller than half of the parameter dimension
d/2, since LDA is a singular statistical model. The dotted blue line drawn
by the circles in this figure represents the learning coeﬃcients of LDA when
the number of topics H is increased. If LDA was a regular statistical model,
its learning coeﬃcient would be the dotted yellow line drawn by the squares.
The behavior of them are so diﬀerent. (b) This figure shows the theoreticallearning curve of LDA and that of a regular statistical model whose parameter
dimension d is same as LDA. The former is the solid blue line and the latter isthe dashed yellow line. The vertical axis means the expected generalization
error E[Gn ] and the horizontal one is the sample size n. This is based on
Theorem 3.4 and the exact value of λ which is clarified by our result. . . . . 78

5.2

List of Figuresvi
5.3

This figure is drawn based on Table 5.2 and Theorem 5.2. It comparesnumerically-calculated RLCTs λ̂ (Numerical λˆ, as the dashed yellow linewith the error bars) and theoretical ones λ (Theoretical λ, as the solid blueline) for H = 2, 3.4.5. The horizontal line means the number of topics H
and the vertical one is the numerically-calculated or theoretical value of the
RLCT. Each error bar of experimental results is the 1-standard deviation
range. The line of Numerical λˆ and that of Theoretical λ are very close andthe standard deviations are suﬃciently smaller than the scale of the RLCTs. . 82

vii

List of Tables
4.1

Numerically Calculated and Theoretical Values of the Learning Coeﬃcients

58

5.1
5.2

Description of Variables in LDA Terminology . . . . . . . . . . . . . . . . 64
Numerically-Calculated and Theoretical Values of RLCTs . . . . . . . . . . 82

1Chapter 1

Introduction
In this chapter, we introduce the research investigated by this dissertation. First,
in Sec. 1.1, we describe the background of this study. In Sec. 1.2, the aim of thisresearch is stated. In Sec. 1.3, the structure of this thesis is explained. And lastly,in Sec. 1.4, we show symbols commonly used in this dissertation.

1.1 Background1.1.1 Statistical Learning Theory
Machine learning is a ubiquitous technology that tackles real-world problems such as future
prediction and knowledge discovery from data, and it has been developing pattern recognition[14]. It is significant to make scientific or technical judgments based on the data; hence,machine learning has been widely applied as same as statistics. For example, non-negative
matrix factorization (NMF) [61, 18] has been used for signal processing [51], text mining[88], bioinformatics [46], and purchase analysis [47]. Another example is latent Dirichlet
allocation (LDA). LDA has been applied to text analysis [15, 30], image recognition [52],market research [72], and geology [98]. Statistical learning theory (learning theory) is atheory that aims to build a foundation for solving these problems faced in the real world andfor evaluating the results obtained by machine learning. When n-data are obtained from thetrue distribution (data-generating distribution) q(x) *1 , we cannot know it in reality. The
objectives of learning theory are the following three things. First, to build a theoreticalfoundation for estimating the true distribution q(x) with a predictive distribution p∗ (x) basedon the data X n = (X1 , . . . , Xn ), by learning it with a model p(x|θ), where x is a variablein the data space X ⊂ RN and θ is a parameter (θ ∈ W ⊂ Rd ). Second, to clarify
the behavior of the error of the learning results. Third, to devise an algorithm (learningalgorithm) for estimating the true distribution with high accuracy in practice. Such situations
specifically emerge when we create artificial intelligence that recognizes images and sounds
by learning from examples, or we conduct knowledge discovery such as reducing dimension
and clustering, or future prediction such as anomaly detection or time series prediction.
*1 Actually, q(x) is a probability density/mass function; however, we call it a probability distribution according
to the convention.

2

Chapter 1 Introduction
Typical learning algorithms include maximum likelihood estimation and Bayesian inference[85]. In maximum likelihood estimation, the predictive distribution becomes p(x|θ∗ ), whereθ∗ is a parameter that maximizes the likelihood (i.e. the probability density when the modelobserves the data). On the other hand, in Bayesian inference, the prior distribution φ(θ)
is defined with the probability model p(x|θ), i.e. the statistical model is denoted by thesimultaneous distribution p(x|θ)φ(θ). The posterior distribution ψ(θ|X n ) ∝ p(x|θ)φ(θ) is
calculated∫ and the predictive distribution becomes the expectation of p(x|θ) by ψ(θ|X n ):p∗ (x) = p(x|θ)ψ(θ|X n )dθ. This prediction distribution is called the Bayesian prediction
distribution. The predictive distribution is often also written as p(x|X n ), in the sense that itis the probability distribution of x given the data X n .
Recently-used statistical models and learning machines often have hierarchical structuresor latent variables. For example, neural networks and matrix factorization have hierarchicalstructures, and hidden Markov model and topic model have latent variables. In such models,a map from the parameter set to the probability distribution set is not one-to-one. Its loglikelihood function cannot be approximated by any quadratic form; its likelihood and posterior
cannot be approximated by any normal distribution. If the above map is injective, then thestatistical model is called regular (or the regular model). The statistical model which is not
regular is called singular (or the singular model). Hence, the learning machines mentioned
above are singular models. Learning theory for singular models is called singular learningtheory. From singular learning theory, neither maximum likelihood estimator nor maximum
a posteriori estimator has asymptotic normality and consistency. Besides, Bayesian inferenceis superior to the maximum likelihood method and maximum a posteriori method in the senseof generalization, i.e. Bayesian inference can make the error between the true distribution
and the predictive one (the generalization error) smaller than those point-estimation methods[79, 82, 85]. This fact has been proved not only numerically but also mathematically. Almost
all practical learning machines are singular; thus, the determination of the generalization
error in singular models is one of the most important issues in machine learning and statistics
community.
In Bayesian inference, the asymptotic behavior of the generalization error had been clarified.
When the statistical model can realize the true distribution and is regular, the following theorem
is proved [4, 85].
Theorem 1.1 Let n be the number of the data and d be the dimension of the parameter space.
The expected Bayesian generalization error E[Gn ] in realizable and regular statistical modelhas the following asymptotic behavior:( )d1
E[Gn ] =
+o
,
(1.1)
2n
n
∫ ∏n
where the operator E[·] means expectation on overall datasets: E[·] = [·] i=1 q(xi )dxi .
In the general case, with supposing some technical assumption, the following theorem holds[79, 82, 85]:
Theorem 1.2 Let n be the number of the data. For the expected Bayesian generalization
error E[Gn ] in a realizable statistical model, there exists the constant λ and the following

1.1 Background3

asymptotic behavior holds:
λE[Gn ] = + o
n

( )
1
.
n

(1.2)

The constant λ is in the leading term of the asymptotic behavior of the generalization error;
thus, it is called the learning coeﬃcient *2 in learning theory. This coeﬃcient is the sameas the real log canonical threshold (RLCT) in algebraic geometry. Conversely, algebraic
geometry is needed to prove Theorem 1.2. In the case that the model is regular, the RLCT isequal to d/2; thus, Theorem 1.2 includes Theorem 1.1 as a special case. Therefore, singularlearning theory can provide general knowledge for statistical modeling and machine learning,in particular when we use Bayesian inference. RLCTs are birational invariants in algebraic
geometry and they depend on the algebraic varieties characterized by statistical models. A
detailed review of these theories will be given in later chapters. It is important for model
selection and hyperparameter tuning to clarify the theoretical generalization error. Besides, if
a theoretical relation between the sample size n and the generalization error Gn is known, thenwe can estimate the suﬃcient sample size to realize the needed generalization performance.
We can also evaluate the correctness of numerical experiments. On the contrary, the absenceof such theoretical facts means that there is no method to verify the accuracy of the numericalexperiments, therefore the correctness of the empirical results cannot be guaranteed. As amore direct application, a precise model-selection method that uses RLCTs had been proposed[24]. RLCTs are also useful to tune the inverse temperature in exchange Monte Carlo [56].
RLCTs of several models are studied. For example, three-layered neural network [80, 9, 8],
reduced rank regression [77, 10], normal mixture model [94], Poisson mixture model [64],hidden Markov model [96], Markov model [100], naive Bayes method [62], Bayesian network[95], Boltzmann machine [97, 6, 7], latent Gaussian tree [23], and etc. RLCTs of them areanalyzed by using Atiyah’s form [11] of resolution of singularity theorem [39].
The above theorems assume that the statistical model can realize the true distribution.
However, since the asymptotic behavior of the generalization error has also been clarified
and they are characterized by RLCTs in the case the model is not realizable [82, 85], therealizability assumption cannot hurt the value of our research; determine of the RLCTs ofstatistical models. Moreover, the realizable case is essential for practical uses because weface a situation that the model is more redundant than the true distribution at the level known
from a given sample in real-data analysis situations. The given sample is finite in real-world
problems, and the complexity of the appropriate model is finite to the extent that it can be
known from that sample. In this case, some models can be redundant, i.e. they seem toinclude the true distribution. If the model is too simple and suﬀers from underfitting, we use
more complex models. It is the goal of model selection and hypothesis testing to compare
appropriate and redundant models in such cases.

1.1.2 Parameter Restriction
From the practical point of view, the parameter region of the statistical model is often restricted
in order to have consistency on domain knowledge and improve interpretability. As a toy*2 The learning coeﬃcient characterizes the learning curve when the sample size n increases, not the number of
epochs. The learning coeﬃcient is not the learning rate in gradient descent methods. In this dissertation, wecall it the real log canonical threshold to avoid confusion between learning coeﬃcient and rate.

Chapter 1 Introduction

4

example, we consider a problem to predict purchase existence for a product by using a logistic
regression model. An usual logistic regression treats linear classification. Suppose that there
are an objective variable y (purchase existence label) and some covariates like advertisement
placements aTVCM , aDM , . . . and ratings r in EC-sites, and the classifier learns a relation
between the labels (purchased or not) and these covariates from the sample. The logistic
regression model is formalized as the following. That the product is purchased is denotedby y = 1. If it is not purchased, then let y be 0. The label y ∈ {0, 1} is considered to begenerated by a Bernoulli distribution
Ber(y|u) = uy (1 − u)1−y ,(1.3)

u = σ(s),(1.4)

s = β0 + β1 aTVCM + β2 aDM + . . . + βd r,(1.5)

σ(s) = 1/(1 + exp(−s)),

(1.6)

where

β0 is an intercept term, and (βi ) are regression coeﬃcients. Intuitively, the contributionsof these covariates are expected to be non-negative *3. However, the coeﬃcients of thecovariates sometimes become negative. If there is a negative coeﬃcient, the contributions ofeach covariate cancel each other out, making it diﬃcult to interpret. Non-negative restriction
to the coeﬃcients (βk ≧ 0 for k = 1, . . . , d) is a solution to resolve such an issue (like Fig.
1.1).

Fig. 1.1: Horizontal bar plots of regression coeﬃcients. Left graph is the non-restricted result
and right one is the result with non-negative restriction.
There are several typical restrictions: non-negative and stochastic. Non-negative restriction,
as the name implies, constrains the parameter space to non-negative real numbers. Stochastic
one constrains the parameter vectors into simplexes, i.e. they are non-negative and the sumof those is equal to one. Practical examples of these restrictions are non-negative matrix
factorization (NMF) and latent Dirichlet allocation (LDA). NMF decomposes an observed*3 Excessive advertising may discourage purchases; however, this eﬀect is ignored in this toy problem.

1.2 Research Goal
matrix into a product of two matrices whose entries are non-negative. Non-negative restriction
in NMF makes discovered-latent patterns clearer to interpret than that of non-restricted cases.
Kohjima has reviewed that NMF is useful and interpretable for analyzing relations between
users and items in purchase logs [48]. LDA is one of the topic models. It represents aword-generating process and its parameters state the probability of words appearance and the
topic proportion; thus, the stochastic restriction is imposed on them [15]. LDA is intuitively
explained as stochastic matrix factorization (SMF): NMF with stochastic restriction [2, 98].
In fact, a mathematical equivalence between LDA and SMF can be proved in the sense ofparameter regions; these RLCTs are the same [38] (see Chap. 5). Moreover, in practice, it
seems to be a large demand for parameter restrictions. A probabilistic programming language
called Stan [28, 17, 71] provides types of variables in order to restrict the parameter region
for statistical models which users make. Thus, parameter-restricted models have been widely
applied. However, although theoretical elucidation of the behavior of the generalization error
is an important problem in learning theory as mentioned in the previous subsection, it has not
yet been clarified how the generalization performance changes with the addition of parameter
restrictions.

1.2 Research Goal
In this dissertation, we theoretically clarify the asymptotic behavior of the Bayesian generalization errors in NMF and LDA, as two of the most popular parameter-restricted models. We
use a resolution of singularity and mathematically analyze the RLCTs of these models.
Specifically, we theoretically consider the Bayesian inference of NMF and LDA based on
the framework described in later chapters and clarify the RLCTs. For NMF, we derive theexact value of the RLCT in some cases, and for the general case, we give a theoretical upperbound by using these exact values. In addition, we clarify the eﬀect of hyperparameters in
the case that the prior is gamma distribution, which is often used in NMF [18, 49]. For LDA,
we prove that the RLCT of LDA is equal to the one of SMF and determine the exact value of
the RLCT in all cases. We also clarify a relationship between LDA and non-restricted matrix
factorization (whose RLCT is equal to one of reduced rank regression).
In this study, we provide theorem proofs based on pure theory. Although it is unknown asa direct problem that the assumptions of the theorems are satisfied in real problems, the real
problem is solved as an inverse problem; thus, the theoretical results under certain assumptions
are essential for solving real issues. Once the theoretical results of certain assumptions areclarified, the structure of the real world can be elucidated by comparing them with thesituations we face in reality. NMF and LDA are widely used; however, their mathematicalcharacteristics have not been clarified. Therefore, this dissertation is fundamental research to
make statistical inferences with NMF and LDA.

1.3 Dissertation Structure
The structure of the rest of this dissertation is as follows. In Chap. 2, we describe a theoreticalframework of Bayesian inference. In Chap. 3, we briefly introduce statistical learning theory
for singular models (singular learning theory). In Chap. 4, the results and discussions of thetheoretical analysis of NMF are presented. Similarly, in Chap. 5, those of LDA are described.

5

Chapter 1 Introduction

6

Lastly, in Chap. 6, we conclude this research. In Appendix A, we summarize the questions
and the answers discussed in the defenses of this dissertation.

1.4

Symbols

In this section, we define the symbols often used in this thesis.
Symbols of sets are denoted as follows. Let N, R and C be the set of positive integers,
real numbers, and complex numbers, respectively. Put D ⊂ R. Let R≧0 := {x ∈ R|x ≧ 0}
and R>0 := {x ∈ R|x > 0}. M(M, N, D) is denoted by the set of M × N matrices whose
entries are in D. Put M, N ∈ N and E ⊂ [0, 1]. Let Onehot(N ) := {w = (wj ) ∈ {0, 1}N |∑N
j=1 wj = 1} = {(1, 0, . . . , 0), . . . , (0, . . . , 0, 1)} be an N -dimensional one-hot vector set∑N
and Sim(N, E) := {c = (cj ) ∈ E N | j=1 cj = 1} be an N -dimensional simplex. Let
S(M, N, E) = Sim(M, E)N be a set of M × N stochastic matrices whose elements are inE.
Besides, symbols of logical operation are defined as follows. Let ∧ and ∨ be logicalconjunction and disjunction, respectively. Also let ⇒ and ⇔ be implication and logicalequivalence, respectively. Iverson bracket is defined as{
1 proposition is true.
[(proposition)] =0 proposition is false.
Let ∼ be a binomial relation such that the functions K1 (w) and K2 (w) have same RLCT
if K1 (w) ∼ K2 (w).

7

Chapter 2

Bayesian Inference
In this chapter, we describe the framework of Bayesian inference. First, in Sec.
2.1, we explain basic concepts to define what Bayesian inference is. Second, inSec. 2.2, we state the framework of Bayesian inference and its definition.
The selected references in this chapter are [85] for Sec. 2.1 and 2.2.

2.1 Basic Concepts
In this section, we introduce a part of probability theory and Kullback-Leibler (KL) divergence
which we need to define Bayesian inference.

2.1.1 Probability Theory
Let X and Y be random variables in measurable spaces (RM , BM ) and (RN , BN ), respectively *1. A pair (X, Y ) is also a measurable function (i.e. a random variable) for a σ
algebra BM ⊗BN , where Bk is denoted by k-dimensional Borel algebra. First, we define thesimultaneous probability distribution p(x, y) of X and Y .
Definition 2.1 (Simultaneous Probability Distribution) A simultaneous probability density function (or distribution) is defined by the probability density function p(x, y) of (X, Y ),where∫p(x, y)dxdy = 1,
and a probability of an event C ∈ BM ⊗BN is∫Pr[(X, Y ) ∈ C] =p(x, y)dxdy.
CMoreover, a simultaneous probability is defined by the above probability Pr[(X, Y ) ∈ C].
*1 Let (Ω, F, P) be a (complete) probability space. Strictly speaking, a random variable X : Ω → R is
defined as a F/BM -measurable map. However, in this dissertation, we do not need to deal with the elementary
event F ∈ F; thus, we define random variables without explicitly specifying the probability space (Ω, F, P).
In fact, statistics and machine learning mainly aim at the data-generating distribution Q : BM → [0, 1](Q(A) := P(X −1 (A))), i.e. a probability measure induced by X and its density q(x) rather than P.
MChapter 2 Bayesian Inference

8

The marginal probability distributions of X and Y are denoted by pX (x) and pY (y),respectively.
Definition 2.2 (Marginal Probability Distribution) About the above simultaneous probability distribution, put∫∫
pX (x) := p(x, y)dy, pY (y) := p(x, y)dx.
Then, pX (x) and pY (y) are the marginal probability distribution of X and Y , respectively.
From the above probability distributions, the conditional probability distribution is defined.
Definition 2.3 (Conditional Probability Distribution) The conditional probability distribution (or density function) of Y given X is defined byp(y|x) :=p(x, y).
pX (x)

This can be interpreted as a probability density function of y with a parameter x. The
conditional probability distribution (or density function) of X given Y is defined in the sameway:p(x, y).
p(x|y) :=pY (y)
Note that the aboves are not defined when the denominator is zero. If the denominator is zero,
the conditional density is defined as zero.
From the definition of the conditional probability distribution, we immediately have thefollowing “theorem”.
Theorem 2.1 (Bayes’s Theorem)
p(x, y) = p(y|x)pX (x) = p(x|y)pY (y)andp(y|x) =p(x|y)pY (y)pX (x)

hold.
Here, we define independent-identically-distributed (i.i.d.) random variables.
Definition 2.4 (independent Random Variables) Let p1 (x) and p2 (y) be the probability
density functions of X and Y , respectively. Two random variables X and Y are called
independent if the simultaneous probability density function p(x, y) of X and Y is equal top1 (x)p2 (y).
Definition 2.5 (Identically-distributed Random Variables) Assume that random variables
X and Y have the same images: Im(X) = Im(Y ). Let p1 (x) and p2 (y) be the probability
density functions of X and Y , respectively. Two random variables X and Y are called
identically-distributed if these probability density functions are equal: p1 = p2 .

2.1 Basic Concepts

9

Definition 2.6 (I.i.d. Random Variables) Two random variables X and Y are called i.i.d.
if they are independent and identically-distributed.
In this thesis, we mainly assume that the random variables X1 , . . . , Xn which mean thesample (X n ∏= (X1 , . . . , Xn )) are i.i.d. Hence, the simultaneous probability of the samplennis q(x ) = i=1 q(xi ), where q(x) is the data-generating distribution and xi ∼ q(x) fori = 1, . . . , n *2.

2.1.2 Kullback-Leibler Divergence
Definition 2.7 (Kullback-Leibler Divergence) Let p(x) and q(x) be probability density
functions on a Euclidean space. Suppose that the supports of p(x) and q(x) are equal
to each other. The Kullback-Leibler (KL) divergence KL is defined by∫q(x)KL(q∥p) := q(x) log
dx.
p(x)If the support of p(x) and q(x) is a discrete set X = {x1 , . . . , xn , . . .}, i.e. p(x) and q(x)are probability mass functions, then the KL divergence is defined byKL(q∥p) :=∞∑q(xi ) log

i=1

q(xi ).
p(xi )The KL divergence is the same as the relative entropy in statistical mechanics. This is basedon a simulation error from an i.i.d. sample generated by the model to the empirical distribution
obtained by the true distribution (Sanov’s theorem) [63]. Moreover, some “distance-like”
properties hold.
Proposition 2.1 KL divergence satisfies the followings:• Non-negativity: KL(q∥p) ≧ 0,
• Non-degenerateness: KL(q∥p) = 0 if and only if q(x) = p(x) almost everywhere(a.e.).
Thus, the KL divergence has been widely used in machine learning and statistics for a metricof discrepancy from a probability distribution to another one. Note that it is not exchangeable:
KL(q∥p) ̸= KL(p∥q) by the definition.
Example 2.1 (Average Code Length) In information theory, an information loss from an
unknown information source q(x) to a receiver p∗ (x) is called an average code length. This
is KL(q∥p∗ ) and means redundancy of the receiver.
Example 2.2 (Generalization Error) In statistics and machine learning, an error in thestatistical inference of the true distribution q(x) by the predictive distribution p∗ (x) is called
the generalization error. This is KL(q∥p∗ ).
M
dis*2 Let Xi be a random variable in (R , BM ) for i = 1, . . . , n. X1 , . . . , Xn induce∩ the same probability∏ntribution Q and its density function is q(x). For arbitrary A1 , . . . , An ∈ BM , Q( n
A
)
=
Q(A
i)
i=1 i
i=1
holds.

Chapter 2 Bayesian Inference

10Here, we prove Proposition 2.1.
Proof of Proposition 2.1. Let f be a function
f : R → R,
f (t) = t + e−t − 1.

Because of f ′ (t) = 1 − e−t and f ′′ (t) = e−t > 0, f is a convex function and it becomes theminimum f (0) = 0 if and only t = 0. Thus, f (t) ≧ 0 and f (t) = 0 ⇔ t = 0. Here, we have∫q(x)KL(q∥p) = q(x) log
dx
p(x)
∫
q(x)
= q(x) log
dx + 1 − 1
p(x)
∫
∫
∫
q(x)
dx + p(x)dx − q(x)dx= q(x) log
p(x)
∫
∫
∫
q(x)
p(x)
= q(x) log
dx + q(x)dx − q(x)dx
p(x)
q(x)
)
∫ (
p(x)
q(x)
+ q(x)− q(x) dx=
q(x) log
p(x)
q(x)
(
)
∫
q(x) p(x)= q(x) log+− 1 dxp(x) q(x)
(
)
∫
q(x)
= q(x)f logdx.
p(x)Owing to the above and q(x) ≧ 0, we obtain KL(q∥p) ≧ 0. Besides, since q(x) = p(x), a.e. xq(x)= 0, a.e. x, we immediately have the second property.
if and only if log p(x)□2.2

Framework of Bayesian InferenceWith the above preparations, we define Bayesian inference. Let a sample X n = (X1 , . . . , Xn )be a collection of i.i.d. random variables subject to a fixed probability distribution q(x). The
sample X n is also a random variable subject toq(x1 , . . . , xn ) =n∏q(xi ),

i=1

and we call q(x) the true probability distribution or the true probability density function.
Simply, it is called the true distribution. We estimate it from the model and the priordistribution (prior), mentioned below. Let W ⊂ Rd be a parameter set. The conditional
probability distribution of x ∈ RN given θ ∈ W is written as p(x|θ) and it is called theprobability model (or the model, simply). Besides, we define the prior distribution (prior) as
a probability distribution of θ ∈ W.

2.2 Framework of Bayesian Inference

11

The structure of the framework of Bayesian inference is as follows. First, we define theposterior distribution (posterior) by using the model, the prior, and the sample. Second, theBayesian predictive distribution is defined by an expectation of the model by the posterior.
Bayesian inference is to infer that “the true distribution q(x) may be the Bayesian predictive
distribution p∗ (x)”; hence, Bayesian inference is a distributional estimation.
We define the posterior as follows.
Definition 2.8 (Posterior Distribution) The posterior distribution (posterior) of the parameter θ given the sample X n is defined byn∏1
ψ(θ|X ) :=φ(θ)
p(Xi |θ),Zni=1nwhere Zn is the normalizing constant which makes ψ(θ|X n ) satisfy
∫
Zn :=

φ(θ)
W

n
∏

∫

ψ(θ|X n )dθ = 1:

p(Xi |θ)dθ.

i=1

The normalizing constant Zn is also called the marginal likelihood. In statistical mechanics,
it is called the partition function.
The Bayesian predictive distribution is defined by the following.
Definition 2.9 (Bayesian Predictive Distribution) The expectation of the model on the parameter with regard to the posterior∫∗p (x) =
p(x|θ)ψ(θ|X n )dθW

is called the Bayesian predictive distribution (or the predictive, simply). To emphasize thatthe predictive is depend on the sample X n , it is also represented by p(x|X n ).
Bayesian inference is statistical inference, not propositional logic-based reasoning. Therefore, it is not correct *3. However, we can mathematically evaluate the error of the inferred
result: how diﬀerent the predictive is from the true distribution. The generalization loss andthe free energy are typical criteria for that.
Definition 2.10 (Generalization Loss and Free Energy) The generalization loss Gn andthe free energy F n are respectively defined by
(∫
)
∫
∫
n
Gn := − q(x) log
p(x|θ)ψ(θ|X )dθ dx = − q(x) log p(x|X n )dx,
W

F n := − log Zn = − log

∫
φ(θ)
W

n
∏

p(Xi |θ)dθ.

i=1

Obviously, minimizing F n is equivalent to maximizing Zn . On the other hand, minimizing
Gn is neither a necessary nor suﬃcient condition for minimizing F n . Two methods have been
*3 There is no “correct” statistical inference. The philosophical interpretation of probability (subjective v.s.
objective or Bayesian v.s. frequentist) never makes either one “correct”.

Chapter 2 Bayesian Inference

12

devised to select an appropriate model: minimizing generalization losses and minimizing
free energy and they are nowadays used as the accuracy of prediction and the certainty ofknowledge discovery, respectively.
Both of the generalization loss and free energy can be normalized *4 in the sense that they
remove the term that depends only on the true distribution. The normalized generalization
loss is called the generalization error.
Definition 2.11 (Normalized Generalization Loss and Free Energy) The normalized generalization loss, i.e. the generalization error Gn is defined by
{∫
(
) }
∫
q(x)
n
Gn := − q(x) log
ψ(θ|X ) exp − logdθ dx.
p(x|θ)W
The normalized free energy Fn is defined as(∫
Fn := − logWφ(θ) exp − logn∏q(Xi )

p(Xi |θ)
i=1

)
dθ.

The meaning of the above normalizing is formalized by the following proposition. All wehave to do is analyzing the generalization error Gn and the normalized free energy Fn since
this normalizing decomposes them the model-depend term and the term which only depends
on the true distribution q(x). The behavior of Gn is described in Theorem 1.2.
Proposition 2.2 Let S and Sn be the entropy and the empirical entropy, respectively:∫
S := − q(x) log q(x)dx,1∑
log q(Xi ).
n i=1nSn := −
The following equalities hold:

Gn = S + Gn ,F n = nSn + Fn .
Proof. Sincelog p(x|θ) = log q(x) − log
i.e.

(

q(x)
,
p(x|θ)q(x)p(x|θ) = q(x) exp − logp(x|θ))*4 Generally, even when dealing with models that are not realizable, the generalization losses and free energies
can be normalized similarly. See also the third chapter in [85].

2.2 Framework of Bayesian Inference

13

holds, we obtain(∫)∫
n
Gn = − q(x) log
p(x|θ)ψ(θ|X )dθ dx
W
(∫
(
)
)
∫
q(x)
n
= − q(x) log
q(x) exp − logψ(θ|X )dθ dxp(x|θ)
W
)
)
(
(
∫
∫
q(x)
n
ψ(θ|X )dθ dx= − q(x) log q(x)exp − logp(x|θ)
W
{
(∫
(
)
)}
∫
q(x)
n
= − q(x) log q(x) + log
exp − logψ(θ|X )dθdxp(x|θ)
W
(∫
(
)
)
∫
∫
q(x)
n
= − q(x) log q(x)dx − q(x) log
exp − logψ(θ|X )dθ dxp(x|θ)W= S + Gn .
Besides, we have
∫
Zn =

φ(θ)
W

n
∏

p(Xi |θ)dθ

i=1
n
∏

(

)
q(Xi )q(Xi ) exp − log
=
φ(θ)
dθ
p(Xi |θ)
W
i=1
)}
( n
){
( n
∫
∏
∑
q(Xi )dθ=φ(θ)
q(Xi )exp −logp(Xi |θ)
W
i=1
i=1
)}
( n
)∫
{
( n
∏
∑
q(Xi )dθ.
=q(Xi )
φ(θ) exp −
log
p(X
|θ)
i
W
i=1
i=1
∫

Thus, the following hold:F n = − log Zn
[( n
)∫
∏
= − logq(Xi ){()} ]q(Xi )
log
φ(θ) exp −dθp(Xi |θ)
W
i=1
i=1
)}
{
( n
∫
n
∑
∑
q(Xi )=−
log q(Xi ) − log
φ(θ) exp −logdθp(Xi |θ)
W
i=1
i=1
n
∑

= nSn + Fn .
□The generalization error in Definition 2.2 seems to be diﬀerent from that in Definition 2.11;
however, we can prove that they are equal if the model is realizable.

14

Chapter 2 Bayesian Inference
Example 2.2. From Definition 2.11, we have
{∫
(
) }
∫
q(x)
n
Gn = − q(x) log
ψ(θ|X ) exp − logdθ dxp(x|θ)
W(∫)∫
p(x|θ)n= − q(x) log
ψ(θ|X )dθ dx
W q(x)(∫)∫
p(x|θ)ψ(θ|X n )dθW= − q(x) log
dx
q(x)
∫
p∗ (x)= − q(x) log
dx
q(x)
∫
q(x)
= q(x) log ∗
dx.
p (x)
Besides, from Definition 2.7, the most-right-hand-side term is KL(q∥p∗ ). Therefore, we have
Example 2.2.
□

15

Chapter 3

Singular Learning Theory
In this chapter, we outline statistical learning theory for singular models: singular learning theory. This chapter is organized as follows. In Sec. 3.1, we
mention the motivations behind applying algebro-geometric methods to the theory of statistical learning. In Sec. 3.2, as a mathematical preparation, we depict
the framework of analyzing Fn and Gn through the use of algebraic geometry.
Thirdly, in Sec. 3.3, we introduce the key consequences of singular learning theory, which shows the relationship between RLCTs and Fn and Gn . Lastly, in Sec.
3.4, as applications of singular learning theory, we describe information criteria:
WAIC, WBIC, sBIC, and WsBIC.
The selected textbooks as the references of this chapter are [81] and [85]. For
introducing sBIC and WsBIC, the author refers the original papers [24, 43]. In
addition, note that much of this chapter is taken from [38], one of the author’s
papers.

3.1 Motivation
First, we explain motivation why we apply algebraic geometry to statistical learning theory.
As described in the above chapter, statistical learning encounters a situation that the true
distribution q(x) is not known although a plurality of data (a.k.a. sample) X n can be obtained,
where the number of data or the sample size is n. Researchers and practitioners design learning
machines or statistical models p(x|θ) to estimate q(x) by making the predictive distribution
p(x|X n ). There is a problem, “How diﬀerent are our model and the true distribution?”
This issue can be characterized as the model selection problem, “Which model is suitable?”This “suitableness” criteria are the normalized free energy Fn and the generalization error
Gn , as mentioned above *1. However, calculating F n is very high cost for computers andGn cannot be computed since q(x) is unknown. These normalized values Fn and Gn also

*1 In actual data analysis situations, the model should be evaluated according to the purpose of analysis andthe domain knowledge of the data. In singular learning theory, we consider Fn and Gn as typical evaluation
criteria used in generic situations [73, 20, 40, 50, 54]. Watanabe (2021) [86], “their proposal is widely acceptedin statistics, data science, and machine learning, on which many statistical systems and learning machines arebeing applied to scientific and practical problems” (p. 2).

Chapter 3 Singular Learning Theory16

depend on∏q(x). We should estimate them from the data. Assume that the likelihood functionnL(θ) = l=1 p(Xl |θ) and the posterior distribution ψ(θ|X n ) can be approximated by aGaussian function of θ. This case is called regular and the model in regular case is calleda regular model. For a regular model, we can estimate Fn and Gn , by using Bayesian
information criterion (BIC) [66] and Akaike information criterion (AIC) [3], respectively.
AIC and BIC are respectively defined by
1∑d
AIC = −
log p(Xi |θ̂) +
n i=1
n
n

and
BIC = −

n∑log p(Xi |θ̂) +i=1

dlog n,
2

where θ̂ is the maximum likelihood estimator or the maximum posterior estimator and d is theparameter dimension. AIC and BIC are derived by not using algebraic geometry; however,they are asymptotically equal to Gn and F n if L(θ) and ψ(θ|X n ) can be approximated by anormal distribution. In order to briefly describe the situation, we define the following error
functions for learning theory.
Definition 3.1 (Error Functions) Suppose p(x|θ) can realize q(x), i.e. there exists θ0 suchthat p(x|θ0 ) = q(x). The function∫q(x)dx = KL(q∥p)K(θ) := q(x) log
p(x|θ)
is called an average error function. In addition, the function
1∑
q(Xi )Kn (θ) :=log
n i=1
p(Xi |θ)n

is called an empirical error function.
In a regular case, the zero point θ0 of K(θ) *2 is unique and ∇2 K(θ0 ) is a strictly positive
definite matrix. Therefore, we can approximate K(θ) as1(θ − θ0 )T ∇2 K(θ∗ )(θ − θ0 )(3.1)
2
by Taylor’s theorem, where there exists t ∈ R (0 < t < 1) such that θ∗ = tθ +(1−t)θ0 *3. Byusing this approximation, Akaike had derived AIC and Schwarz had done BIC, respectively.
However, in general, we cannot estimate Gn and Fn by using AIC and BIC. This is because
the above approximation does not hold since ∇2 K(θ0 ) has eigenvalues which are zeros. θ0
is not uniquely determined. K −1 (0) includes singularities in the parameter space. Thus, weneed algebraic geometry to study the singularities.
K(θ) =−1*2 Because the model can realize the true distribution, argminK(θ) is equal to K (0). If p(x|θ) is notrealizable q(x), we can obtain similar results by considering argminK(θ) instead of K −1 (0). For example,
Takeuchi information criterion (TIC) [70, 85] and widely applicable information criterion (WAIC) [82] do notneed the realizableness assumption.
*3 From the definition, K(θ0 ) = 0 and ∇K(θ0 ) = 0 hold.

3.2 Singularity Resolution Theorem and Zeta Function

3.2 Singularity Resolution Theorem and Zeta Function
Second, the framework of analyzing Gn and Fn , which uses algebraic geometry, is explained.
We consider K(θ) and its zero set K −1 (0): this is an algebraic variety which can havesingularities. We use the following form by [11] of the singularities resolution theorem [39].
This form was originally derived by Atiyah for the analysis of distributions (hyperfunctions);
however, Watanabe proved that it is useful for constructing singular learning theory [79, 82,
81, 85].
Theorem 3.1 (Singularity Resolution Theorem) Let F be a non-negative analytic function
on the open set W ′ ⊂ Rd and assume that there exists θ ∈ W ′ such that F (θ) = 0. Then,
there are d-dimensional manifold M and an analytic map g : M → W ′ such that for each
local chart of M,
2kd
1
F (g(u)) = u2k
1 . . . ud ,|g ′ (u)| = b(u)|uh1 1 . . . uhd d |,where |g ′ (u)| is the Jacobian of g, kj and hj are non-negative integers and b : M → R isstrictly positive analytic: b(u) > 0.
A pair of the above manifold and map (M, g) is called a resolution of singularity.
This theorem does not remove singularities but makes them easier to handle. In general,a zero set of a non-negative analytic function F has singularities; however, (F ◦ g)−1 (0)2kd
1
has only singularities which are the following form: the zero set of u2k
1 . . . ud . Suchsingularities are called normal-crossing singularities.
Thanks to Theorem 3.1, the following analytic theorem is proved [11, 12, 65].
Theorem 3.2 Let F : Rd → R be an analytic function of a variable θ ∈ Rd . a : W → R isdenoted by a C ∞ -function with compact support W. Then∫ζ(z) =|F (θ)|z a(θ)dθW

is a holomorphic function on Re(z) > 0. Moreover, ζ(z) can be analytically continued toa unique meromorphic function on the entire complex plane C. The poles of the extended
function are all negative rational numbers.
The KL divergence is non-negative and analytic; thus, we can apply Theorem 3.1 to K(θ)
on K −1 (0) ∩ W ′ . Then, we obtain
2kd
1
K(g(u)) = u2k
1 . . . ud ,|g ′ (u)| = b(u)|uh1 1 . . . uhd d |.
Assuming the domain of prior φ(θ) is W and it satisfies W ⊂ W ′ , we can also apply Theorem
3.2 to (K(θ), φ(θ)) and obtain the following zeta function.

17

Chapter 3 Singular Learning Theory

18

Definition 3.2 (Zeta Function in Statistical Learning Theory) Let K(θ) be the average error function. Suppose that the support of the prior φ(θ) is compact. The zeta function in
statistical learning theory is defined by the following univariate complex function:∫
ζ(z) =
K(θ)z φ(θ)dθ.
(3.2)
W

If the prior can be zero or infinity on the zero set K −1 (0) of the average error function,then we must consider the zeta function above as defined by Definition 3.2. However, if theprior is positive and bounded on K −1 (0), then the prior has no eﬀect on the maximum pole.
With this in mind, we can consider the following zeta function:∫
ζ(z) =
K(θ)z dθ.
(3.3)
W

∫

Although K(θ)z φ(θ)dθ is a holomorphic function on Re(z) > 0, we can prove thatthe zeta function in statistical learning theory has an analytic continuation on C as a uniquemeromorphic function and its poles are negative rational numbers. Applying Theorem 3.1 toK(θ)φ(θ), there exists a manifold M and an analytic map g such that
2kd
1
K(g(u)) = u2k
1 . . . ud ,φ(g(u))|g ′ (u)| = b(u)|uh1 1 . . . uhd d |for each chart of M, where |g ′ (u)| is a Jacobian of g and b is a positive-analytic function. Let
2kd
1k = (k1 , . . . , kd ) and h = (h1 , . . . , hd ) be non-negative multi-indexes. Put u2k=:
1 . . . ud
∑
hd
h1
2k
h
u and u1 . . . ud =: u . Considering a partition of unity a ϕa (u) for M, the zetafunction becomes as follows:∑∫
ζ(z) =
K(g(u))z |g ′ (u)|ϕa (u)du
a

∑∫

=

u2kz |uh |b(u)ϕa (u)du

a

∑∫

=

a

u2kz uh b(u)du.

[0,1]dSince b(u) is positive, it does not eﬀect the maximum pole of ζ(z). Hence, we only have to
treat
∫
d ∫ 1
∏
2k z+hj2kz hu u du =uj j
duj .
[0,1]dj=1

0

Calculating the integral, we have
∫
2kz h

u

u du =[0,1]d

d ∫ 1
∏
j=1

=

d
∏

2k z+hj

uj j

duj

0

1
.
2kj z + hj + 1
j=1

3.2 Singularity Resolution Theorem and Zeta Function
Thus, allowing for duplication, the poles (−λj )dj=1 are as follows:z = −λ1 := −h1 + 12k1..
.
z = −λj := −hj + 12kj..
.
z = −λd := −hd + 1.
2kd

The zeta function in statistical learning theory is a meromorphic function which is the resultof analytic continuation. Hence, its poles are negative rational numbers. The definition of areal log canonical threshold is the negative maximum pole of it.
Definition 3.3 (Real Log Canonical Threshold and its Multiplicity) Let k = (k1 , . . . , kd )and h = (h1 , . . . , hd ) be non-negative d-dimensional multi-indexes. Assume that k =(k1 , . . . , kd ) is at least one of them positive. For a sequence
hj + 12kj(j = 1, . . . , d),if kj = 0, we define
hj + 1= ∞.
2kj
The minimum of the above sequence
hj + 1j=1
2kj
d

λ := minis called a real log canonical threshold (RLCT) and the number m of elements which areequal to λ is called its multiplicity:
m :=

d [
∑
hj + 1
j=1

2kj

]
=λ ,

where [(proposition)] constitutes Iverson bracket.
From the above discussion of the zeta function and the definition of an RLCT, an RLCT
λ is a sign reversal of the maximum pole (−λ) of the zeta function ζ(z), and its multiplicity
m is equal to the order of the maximum pole of the zeta function. Namely, the maximum
pole is corresponding to the deepest singularity in the analytic set K −1 (0). Furthermore, an
RLCT is independent of methods of singularity resolution; an RLCT is a birational invariant
in algebraic geometry. The RLCT is determined only for the statistical model p(x|θ), theprior φ(θ) and the true distribution q(x) although there are infinite∫ pairs (M, g) to resolve thesingularities in the zero set of the given analytic function K(θ) = q(x) log(q(x)/p(x|θ))dx.

19

Chapter 3 Singular Learning Theory

20

An RLCT can be interpreted as a volume dimension of the analytic set K −1 (0) in theparameter space W [81]. Let t > 0.
Theorem 3.3 Let V (t) be the volume of the set {θ ∈ W | K(θ) < t} in the sense of themeasure φ(θ)dθ, i.e.
∫
V (t) :=

φ(θ)dθ.
K(θ)<tThen, the following equality holds:
log V (t)
.
t→+0
log t

λ = limThis is a base of a method to numerically compute an RLCT by using Monte Carlosimulation [89]. Besides, the above limit is similar to a Minkowski dimension, one of fractal
dimensions.
Definition 3.4 (Minkowski Dimension) Let S be a subset of Rd . The Minkowski dimensionof S is defined as the following d∗ :log Vnbhd (S, t)
,
t→+0
log t

d∗ := d − lim

where Vnbhd (S, t) is the volume of the t-neighborhood of S:∫
Vnbhd (S, t) :=dθ, dist(S, θ) = inf{∥s − θ∥ | s ∈ S}.
dist(S,θ)<tλ and d∗ are similar but diﬀerent concepts because λ is a positive rational number whereasd∗ can be an irrational number. However, both of λ and d∗ can be interpreted as intrinsic
dimensions of subsets included by a Euclidean space. Indeed, they are useful for learningtheory. As described later, the RLCT of K(θ) = KL(q||p) dominates the Bayesian generalization error and the free energy (see Theorem 3.4). On the other hand, convergence
rates of the approximation and generalization errors by deep neural networks depend on the
Minkowski dimension of the data and the rates are optimal in the minimax sense [58]. Thus,
we can refer an RLCT to an intrinsic dimension of the model.
In order to compute the exact value of an RLCT for a particular statistical model, one needs
to find the manifold M and the analytic map g in the Singularity Resolution Theorem. By using
the case when it is easy to calculate the exact value of the RLCT, some studies have elucidated
the theoretical upper bound of the RLCT. As mentioned in Sec. 1.1, one of statistical models
whose RLCT is well-known is reduced rank regression [10]. In this case, the upper bound of
the RLCT was derived before its exact value was determined [77]. In almost all of the others,
upper bounds are only clarified, i.e. it is much more diﬃcult to find the exact value of the RLCT
than to derive a non-trivial upper bound of that. Even clarifying a non-conservative upperbound of the RLCT is challenging since there is no standard method to find a resolution of
singularities (M, g) for a collection of analytic functions like the KL divergence of statisticalmodels. Instead, researchers have been studying RLCTs by developing novel methods toanalyze RLCTs for each statistical model [80, 9, 8, 94, 64, 96, 100, 62, 95, 97, 6, 7, 23] (seealso Sec. 1.1).

3.3 Key Results of Singular Learning Theory

3.3 Key Results of Singular Learning Theory
We introduce the theorem which shows the relationship between RLCTs and Gn and Fn[79, 81, 85].
Theorem 3.4 (Watanabe) Let q(x), p(x|θ), and φ(θ) be the true distribution, the learning
machine, and the prior distribution, where x is a point of RN and θ is an element of thecompact subset W of Rd . Let K(θ) be the average error function and λ is denoted by the
RLCT of (K(θ), φ(θ)). If there exists at least one θ0 such that q(x) = p(x|θ0 ) (i.e. the model
can realize the true distribution), then the asymptotic behavior of the generalization error Gn
and the normalized free energy Fn is as belows:
(
)
m−1
1
λ
+o
,
E[Gn ] = −n n log nn log nFn = λ log n − (m − 1) log log n + Op (1).
Theorem 3.5 (Watanabe) If there exists at least one θ0 such that q(x) = p(x|θ0 ) andmaximum likelihood or posterior method is applied (i.e. the predictive distribution is p∗ (x) =p(x|θ̂), where θ̂ is the maximum likelihood or posterior estimator), then there is a constantµ > d/2 such that the asymptotic behavior of the generalization error Gn is as belows:
( )
1
µ
.
E[Gn ] = + onnK(θ) depends on q(x) and p(x|θ); thus, Theorem 3.4 can be understood as we can clarify Gn
and Fn if the RLCT is clarified, which is determined by (q(x), p(x|θ), φ(θ)). As introducedin Chap. 1, there are several researches to find the RLCT of a statistical model by analyzing
the maximum pole of the zeta function. Their studies are based on Theorem 3.4 and the zetafunction derived by Theorem 3.2. The researchers have found the singularity resolution map
g for the exact value or an upper bound of Φ(θ), and have obtained the one of the RLCT since
the RLCT is order isomorphic: if ∫Φ(θ) ≦ Ψ(θ), then λΦ ∫≦ λΨ , where (−λΦ ) and (−λΨ )are the maximum pole of ζ1 (z) = Φ(θ)z dθ and ζ2 (z) = Ψ(θ)z dθ, respectively [81].
Recently, RLCTs of statistical models have been used in other ways. Analyzing andtuning exchange probabilities in exchange Monte Carlo methods [56], deriving the asymptotic
behavior of the Bayesian estimation accuracy for latent variables [90, 91], devising singular
Bayesian information criterion [24] (see Sec. 3.4) and constructing deep learning theory[55] are such examples. Singular learning theory assumes that the optimal distribution isessentially unique [85]; however, in some of cases when that assumption is not satisfied,
RLCTs and its multiplicities also determine the behavior of the free energy [57]. Moreover,singular learning theory in non-i.i.d. cases is also studied such for conditional independent
samples [84, 85], in the exchangeable case [87], and for structured data [93].
From the practical point of view, Theorem 3.5 shows that Bayesian inference makes the freeenergy and the generalization error smaller than maximum likelihood or posterior method in
singular case since µ > d/2 ≧ λ [85]. Hence, if the RLCT is clarified, then we can draw theE[Gn ]-n learning curve and estimate the sample size which satisfies the required inference
performance.

21

Chapter 3 Singular Learning Theory

22

There is no standard method to find an RLCT for a statistical model (family of functions).
Here, we show a fundamental method to find the RLCT for a non-negative analytic function:
this is called blowing-up [39]. We explain blowing-up which is used to study learning
machines, based on a concrete example and [81, 85]. If a reader needs the rigorous definitionof blowing-up, then see [39]. Let
K(θ) = θ12 + . . . + θd2(3.4)

and θi , i = 1, . . . , d be independent parameters. Especially, we treat the case d = 2.
Blowing-up of K(θ) is a transformation of the coordinate that is defined{(2) (2)(1)
θ1 = θ1 = θ1 θ2 ,.
(2)(1) (1)θ2 = θ1 θ2 = θ2
Using this blowing-up,
(1)

(1)

(2)

(2)

K(θ) = (θ1 )2 {1 + (θ2 )2 } = (θ2 )2 {(θ1 )2 + 1},
and the absolute of Jacobian |J| of this transformation is
(1)

(2)

|J| = |θ1 | = |θ2 |.
(j)

From the applied mathematical point of view, 1 + (θi )2 is strictly positive thus the RLCT
can be calculated. The zeta function ζ(z) is∫∫(1) 2z+1(1) 2 z(2)(2)(1)
ζ(z) = (θ1 ){1 + (θ2 ) } dθ = (θ2 )2z+1 {1 + (θ1 )2 }z dθ(2)(j)

and it is immediately proved that 1 + (θi )2 does not eﬀect the RLCT λ. Then all we haveto consider is the function∫∫(1) 2z+1(2)(1)
ζ(z) = (θ1 )dθ = (θ2 )2z+1 dθ(2) ,which are analytically connected to C as a unique meromorphic function, respectively. Therefore, we getc1c2
ζ(z) ==,
2(z + 1)2(z + 1)andλ = min{1, 1} = 1,
where c1 and c2 are positive constants. By the same way, in general d, the RLCT λ is equal
to
d
λ= .
(3.5)
2
If the statistical model is regular, then its average error function can be approximated bya quadratic form as Eq. (3.1). By diagonalizing the positive-definite matrix ∇2 K(θ0 ) andcentralizing that form, we arrive at Eq. (3.4). Because of Eq. (3.5), we can reconstruct theasymptotic behavior of the Bayesian generalization error in regular case: Theorem 1.1.
Let ∼ be a binomial relation such that the functions K1 (w) and K2 (w) have same RLCT
if K1 (w) ∼ K2 (w). There are some propositions to evaluate RLCTs [81].

3.3 Key Results of Singular Learning Theory

23

Proposition 3.1 Let Φ : Rd → [0, ∞) be a non-negative analytic function. Let λ and m be
the RLCT and its multiplicity defined by the following zeta function∫
ζ(z) = Φ(θ)z φ(θ)dθ,
respectively. Let a(θ) > 0 and b(θ) > 0 be smooth functions of θ which are strictly positiveon the neighborhood of Φ−1 (0). Then, the following zeta function has the same maximum
pole as that of the above: the RLCT and its multiplicity are same∫
ζa,b (z) = (a(θ)Φ(θ))z b(θ)φ(θ)dθ.
Proposition 3.2 RLCTs save the orders: if Φ(θ) ≦ Ψ(θ) and φ∫1 (θ) = φ2 (θ), then λΦ ≦ λΨ ,where(−λΦ ) and (−λΨ ) are the maximum pole of ζ1 (z) = Φ(θ)z φ1 (θ)dθ and ζ2 (z) =∫Ψ(θ)z φ2 (θ)dθ, respectively.
From the aboves, the following is immediately derived.
Proposition 3.3 Let λΦ and λΨ be RLCTs in the above proposition. If φ1 (θ) = φ2 (θ) andthere are two positive constants c1 > 0 and c2 > 0 such thatc1 Φ(θ) ≦ Ψ(θ) ≦ c2 Φ(θ),
then λΦ = λΨ holds and their multiplicities are also same.
By using the relationship between ideals and analytic sets, the following property can beproved.
Proposition 3.4 Suppose s, t ∈ N, and let f1 (w), . . . , fs (w), g1 (w), . . . , gt (w) be real polynomials. Furthermore, let
I := ⟨f1 , . . . , fs ⟩,J := ⟨g1 , . . . , gt ⟩be the generated ideals of (f1 , . . . , fs ) and (g1 , . . . , gt ), respectively. We put
F (w) :=s∑2fi (w) ,G(w) :=

i=1

Then, F ∼ G if I = J.
Corollary 3.1 Assume that F (w) =t∑gj (w)2 .

j=1∑s

F (w) +

2
i=1 fi (w) . Then

( s
∑

)2
fi (w)∼ F (w).

i=1

The prior can aﬀect to the RLCT as the following.
Proposition 3.5 Let λ1 := λΦ and λ2 := λΨ be RLCTs in Proposition 3.2. If Φ(θ) = Ψ(θ)
and φ1 (θ) ≦ φ2 (θ), then λ1 ≧ λ2 holds.

Chapter 3 Singular Learning Theory

24

3.4

Information Criteria from Singular Learning Theory

From the practical point of view, singular learning theory provides novel information criteria:
WAIC [82] and WBIC [83]. They are useful even if the posterior cannot be approximated byany normal distribution, i.e. neither AIC nor BIC can be applied.
Definition 3.5 (Widely Applicable Information Criterion (WAIC)) Widely applicable information criterion (WAIC) is defined by the following random variable Wn :Wn := Tn + Vn /n,where Tn is the empirical loss and Vn is the functional variance:
n
n
1∑
1∑
∗
Tn = −
log p (Xi ) = −log Eθ [p(Xi |θ)],
n i=1
n i=1

Vn =

n [
∑

Eθ [(log p(Xi |θ)) ] − {Eθ [log p(Xi |θ)]}
2

2

]
=

n
∑

i=1

Vθ [log p(Xi |θ)].

i=1

Wn − Sn is called the WAIC error.
WAIC asymptotically approximates the Bayesian generalization loss in the general case *4 inthe sense of the expectation [82].
Theorem 3.6 The expected WAIC is asymptotically equal to the expected generalization loss
and these diﬀerence is O(1/n2 ):E[Wn ] = E[Gn ] + O(1/n2 ).
This means that WAIC is asymptotically equal to leave-one-out cross validation (LOOCV)
loss. LOOCV loss is collapsed in the case when the data has leverage points if LOOCV loss
is calculated by the importance sampling (ISCV) [25]; however, WAIC seems consistent inthat case [54]. Some outliers can become leverage points in practical uses. Thus, WAIC isconsidered to be numerically more robust than the ISCV.
The mathematical property of the variance of WAIC is as follows.
Theorem 3.7 The WAIC error has same variance as the Bayesian generalization error:
Wn − Sn + Gn = 2λ/n + op (1/n).
Note that Wn − Sn and Gn are random variables but the summation of them is asymptoticallydeterministic 2λ/n. This property means that the WAIC error has an inverse correlation to
the Bayesian generalization error.
Widely applicable Bayesian information criterion (WBIC) is defined as follows. It uses the
expectation by the tempered posterior distribution.
*4 Some technical assumptions (such the relatively finite variance of the model) are needed but they are practicallyconsistent [85]. As far as we know, widely used statistical models are considered to satisfy the assumptions.

3.4 Information Criteria from Singular Learning Theory
Definition 3.6 (Tempered Posterior) Let β > 0 be a positive constant. The tempered posterior ψ β (θ|X n ) is defined by
n
∏
1
φ(θ)
p(Xi |θ)β ,Zn (β)i=1

ψ β (θ|X n ) :=where the partition function Zn (β) is equal to∫Zn (β) =

φ(θ)
W

n
∏

p(Xi |θ)β dθ

i=1and the scale constant β is called the inverse temperature in statistical mechanics.
Definition 3.7 (Widely Applicable Bayesian Information Criterion (WBIC)) Let Eβθ bethe expectation operator by the tempered posterior:∫βEθ [·] :=[·]ψ β (θ|X n )dθ.
W

Widely applicable Bayesian information criterion (WBIC) is defined by the following randomvariable Wn′ :Wn′ := Eβθ [nLn (θ)],
where Ln (θ) is the negative log likelihood
1∑
log p(Xi |θ)Ln (θ) = −n i=1nand β = 1/ log n.
WBIC asymptotically approximates the free energy.
Theorem 3.8 WBIC Wn′ satisfies the following asymptotic equality:√
Wn′ = F n + Op ( log n).
The probability that the true model is selected by minimizing WBIC converges to one whenn → ∞ Moreover, if the model is regular, then the diﬀerence between WBIC and BIC
converges to zero in probability.
√
This is because Wn′ = nSn + λ log n + Op ( log n) and F n = nSn + λ log n − (m −1) log log n + Op (1) have been proved. In this way, WBIC matches the free energy only up
to the leading term. WBIC tends to underestimate the free energy and the accuracy of themodel selection may be low. The underestimating term has been theoretically analyzed andan adjustment term has been proposed in order to make WBIC asymptotically unbiased [44].
Also, more precise approximation methods have been proposed: sBIC [24] and WsBIC [43].
Singular Bayesian information criterion (sBIC) have been proposed by Drton and Plummer[24]. sBIC is derived as a novel expansion of BIC based on Theorem 3.4. Let θ̂ be the (local)

25

Chapter 3 Singular Learning Theory26

maximum likelihood estimator. Namely, sBIC is defined as the first and second leading termof the free energy:sBIC ≜ nLn (θ̂) + λ log n − (m − 1) log log n
and the sBIC seems to match the free energy higher-order terms than WBIC, where A ≜ B
means that A is namely (intuitively) equal to B but is not mathematically equal to it *5.
However, in fact, sBIC is calculated by solving fixed point simultaneous equation system andhas consistency, i.e. the probability that the true model is selected converges to one whenn → ∞ [24]. sBIC uses the theoretical value of λ and m. In general, they are depend on
the true distribution; however, in computing sBIC, the true distribution is referred to one ofsubmodels. For example, we theoretically consider a model selection problem: selecting thenumber of hidden units i in three-layered neural network model (more general, the considered
model is dominated by a control variable i ∈ N ∪ {0} and all submodels are monotonically
included: model j ⊂ model i if j < i). For the sake of simplicity, the model whose numberof hidden units is i is called the model i (written as pi (x|θ)) or that the model is i. In thesame way, let φi (θ) be the prior when the model is i. When the candidate models are set
to i = 0, . . . , H, for each i, the submodels are considered to j = 0, . . . , i. Let λij and mij
be the RLCT and its multiplicity when the model is i and the true distribution is the modelj, respectively. If the RLCT and its multiplicity are clarified, we can immediately have thefollowing two matrices:


λ00
0
...
0
 λ10 λ11 . . .
0 


Λ= .
(3.6)
..
..  ,
.
.
.
 .
.
.
. 
λH0

λH1

...

λHH

m00
 m10

M = .
 ..

0
m11
..
.

...
...
..
.

0
0
..
.

mH0

mH1

...

mHH






.


(3.7)

If an upper bound of the RLCT is only clarified, then the matrix Λ is constructed by thebounds and the positive entries of M are set to one. Let θ̂i be the (local) maximum likelihood
estimator (MLE) of the model i and L′ij satisfy− log L′ij = nL(i)(3.8)
n (θ̂i ) + λij log n − (mij − 1) log log n,∏n(i)where the likelihood of the model i is denoted by Pi (X n |θ) = l=1 pi (Xl |θ) and Ln (θ) =− n1 log Pi (X n |θ). Obviously, we have(log n)mij −1.
nλij
In this case, sBIC is calculated as Algorithm 1 *6.
L′ij := Pi (X n |θ̂i )(3.9)*5 In scare quotes sense, A is “equal” to B (but A ̸= B).
*6 The original paper defines sBIC as approximation of the log marginal likelihood −F n . Thus, in this thesis,
the sign is reversed in some parts compared to the original version.

3.4 Information Criteria from Singular Learning Theory

27

Algorithm 1 Calculation Algorithm of sBIC
Require: H:candidate model size, Λ:matrix of RLCTs, M :matrix of multiplicities
Pi :likelihood of the model i, X n :data, θ̂i : MLE for the model i
Ensure: Calculation of sBIC
Allocate L′ij := Pi (X n |θ̂i )(log n)mij −1 /nλij
Allocate arrays L, sBIC ← initialize
for i = 0 to H do
if i == 0 thenL[i] ← L′ii
sBIC[i] ← − log L[i]else∑φ (θ )b ← −L′ii + j<i L[j] φji (θji )∑φ (θ )c ← j<i L[j]L′ij φji (θji )√L[i] ← (−b + b2 + 4c)/2sBIC[i] ← − log L[i]end ifend forreturn sBIC

In this algorithm, if i ̸= 0, we solve the quadratic equation system
L[i]2 + bL[i] − 4c = 0, i = 1, . . . , H.

(3.10)

It is proved that this simultaneous equation√ has a unique positive solution [24]. Therefore,
we compute and plug-in L[i] ← (−b + b2 + 4c)/2. The value of sBIC for the model i
is sBIC[i] in the above pseudo-code. Thus, the model argmin{sBIC} is selected. ActuallysBIC requires theoretical RLCTs; however, if they are clarified, it is empirically known thatsBIC is more precise than WBIC in the sense of the accuracy of the model selection [24, 43].
For real computation, the value of the likelihood becomes small not to correctly evaluate theresult of floating-point arithmetics in a computer. However, if the purpose of the computation
is to solve the model selection problem by sBIC, we need the order of the maximal likelihoods{Pi (X n |θ̂i )} of the model i as well as the exact maximal value of the maximal likelihood
MLn := maxi {Pi (X n |θ̂i )}. Hence, we can normalize sBIC by using Pi (X n |θ̂i )/MLn
instead of Pi (X n |θ̂i ) for i = 0, . . . , H in order to make the computation stable (private
communication with Mathias Drton and Fumito Nakamura).
It is well-known that MLEs of singular models are neither stable nor unique even if thesample size is large enough to be treated as an asymptotic scale. Thus, there are some heuristics
to make sBIC more precise. One of them is using variational Bayes estimator (VBE) insteadof MLE [35]. For some of singular models like Gaussian mixture model (GMM), we canderive expectation-maximization (EM) algorithm to calculate its MLE; however, it is provedthat there is no MLE in such mixture models [32]. As far as we know, when EM algorithm
can be derived, deterministic variational Bayes method can also be constructed, for example
GMM; hence, this heuristics does little to reduce the applicability of sBIC. According to thetechnical report [35], the model selection with the original (using MLE calculated by EM)

Chapter 3 Singular Learning Theory

28

sBIC failed and that with the proposed (using VBE) method succeeded when the model is
GMM with non-diagonal variance-covariance matrices. Another heuristics is considering
the diﬀerence of the RLCTs and sBIC for i = 0, . . . , H [69]. This heuristic method returns
the first point when the slope of the sequence (i, sBIC[i]) becomes negative as the model
selection result.
Lastly, we introduce widely applicable singular Bayesian information criterion (WsBIC).
WsBIC is defined by a variant of sBIC which uses numerically-calculated RLCTs insteadof theoretical ones [43]. In order to obtain numerical RLCTs, Imai proposed a consistentestimator λ̂sV of an RLCT, where s ∈ N is the number of independent simulations.
Definition 3.8 (Imai’s estimator of an RLCT) Let Vβθ be the variance operator along withthe tempered posterior:
[
{
}2 ]
β
β
β
2
Vθ [·] = Eθ [(·) ] − Eθ [·].
λ̂1V (X n ) is defined as

λ̂1V (X n ) = β 2 Vβθ [log P (X n |θ)],where P (X n |θ) is the likelihood: log P (X n |θ) = −nLn (θ).
Assume that there are (simulated) independent datasets (D1 , . . . , Ds ) and these size is sameas nsimulated , i.e. |Dk | = nsimulated and Dk is generated by the submodel for k = 1, . . . , s.
Let λ̂sV be defined ass1∑ 1
λ̂V (Dk ).
λ̂sV =s
k=1

Theorem 3.9 (Imai) Let λ be the theoretical RLCT of the considered model. λ̂1V (X n ) is aconsistent estimator, i.e. we have√λ̂1V (X n ) = λ + Op (1/ log n).
Besides, λ̂sV is asymptotically normal. In other words, the asymptotic distribution of λ̂sV isdetermined as
λ̂sV →d N (λ, σλ ),
where
σλ2 =

vM
,
4s log nsimulated

vM is a positive constant and N (µ, σ) be a normal distribution whose mean and standard
deviation are µ and σ, respectively.
Imai’s estimator contains useful properties as the above. Moreover, it is also proved that
there is a unique inverse temperature β0 = 1/ log n + op (1/ log n) which is depend on(q(x), p(x|θ), φ(θ)) and makes λ̂1V (X n ) an unbiased estimator of λ. On the other hand, it
requires Markov chain Monte Carlo (MCMC) sampling from the tempered posterior; thus,
it cannot directly apply to Gibbs sampling which is diﬃcult to sample from the tempered
posterior in general.

3.4 Information Criteria from Singular Learning Theory
WsBIC is obtained by replacing the theoretical RLCT λ to Imai’s estimator λ̂sV . λ̂sV canbe calculated if the tempered posterior is realized; thus, the applicability of WsBIC is nearly
equal to that of WBIC. Moreover, WsBIC is based on sBIC; thus, it is also more precise than
WBIC nevertheless sBIC is less widely applicable than WBIC. Except for the computationalcost, WsBIC outperforms sBIC and WBIC in terms of accuracy and applicability.

29

31

Chapter 4

Bayesian Generalization Error in
Non-negative Matrix Factorization
In this chapter, we report the result of theoretical analysis for Bayesian generalization error in NMF. This chapter consists of five parts. First, in Sec. 4.1, we
describe motivation of this theoretical study. Second, in Sec. 4.2, we state themain theorem with regard to Bayesian generalization error in NMF. Third, in Sec.
4.3, we prepare lemmas for the proof of the theorem. Fourth, in Sec 4.4, we provethe main theorem. Lastly, in Sec. 4.5, we discuss the theoretical results.
In the followings, θ = (U, V ) is a parameter and x = X is an observed randomvariable. Note that this chapter is based on the author’s papers [37, 36, 33].

4.1 Motivation
Non-negative matrix factorization (NMF) [61, 18] has been applied to text mining [88], signal
processing [51], bioinformatics [46], consumer analysis [47], and recommender systems [16].
NMF experiments discover the knowledge and predict the future unknown structures in thereal world, however, the method suﬀers from many local minima and seldom reaches theglobal minimum. In addition, the results of numerical experiments strongly depend on theinitial values; a rigorous method has not yet been established.
In order to resolve this diﬃculty, Bayesian inference for NMF has been established [18]. It
uses, for numerical calculation of the Bayesian posterior distribution, Gibbs sampling method
which is a kind of Markov chain Monte Carlo method (MCMC). Bayesian NMF is known as
a more robust method than usual recursive methods of NMF since it numerically realizes theposterior distribution; the parameters are subject to a probability distribution and that makes
it possible to grasp the degree of fluctuation of the learning/inference result. As is describedlater, in general, Bayesian method has higher estimation accuracy than maximum likelihood
estimation and maximum posterior estimation if the model has hierarchical structures or
hidden variables, like NMF.
On the other hand, the variational Bayesian algorithm (VB) for NMF has also been established [18], with being inspired the mean field approximation. The variational Bayesian NMF
algorithm (VBNMF) also results more numerically stable than usual recursive algorithms as

Chapter 4 Bayesian Generalization Error in Non-negative Matrix Factorization

32

VB approximates the Bayesian posterior distribution. Moreover, VBNMF computes fasterthan usual Bayesian inference such an MCMC. However, its free energy (called the variational
free energy) is larger than the Bayesian free energy, since VB ascends the evidence lower
bound but it is not the true model evidence. Note that the marginal likelihood is also called
the model evidence and the negative logarithm value of the evidence lower bound is equal tothe variational free energy. From the above, it is important to clarify the approximation errorof variational inference for not only theoretical reasons but also practical points of view.
As mentioned Chap. 3, researchers have studied RLCTs of many statistical models toclarify theoretical behaviors of the generalization error and the free energy. Moreover, forseveral statistical models, the variational free energy was proved that it asymptotically equals
nSn + λvb log n + Op (1),where λvb is a learning coeﬃcient and it depends on the model.
Normal mixture models [78], hidden Markov models [42], and NMF [49] are such examples.
Kohjima proved the following theorem about the learning coeﬃcient of VBNMF [49].
Theorem 4.1 (Kohjima) Let the elements of the data matrices xij (i = 1, . . . , M ; j =1, . . . , N ) be independently generated from the Poisson distribution whose mean is equal tothe (i, j) element of U0 V0 , where the number of columns in U0 (= the number of rows in V0 )is H0 ; called the non-negative rank of U0 V0 [19].
Let the likelihood model and the prior be the following Poisson and gamma distributions,
respectively:
M ∏N∏((U V )ij )xij

e−(U V )ij ,x!
ij
i=1 j=1
)
(
)
(
H ∏
N
M ∏
H
ϕV
ϕU
∏
∏
θ
θU
ϕV −1 −θV vkj
VuϕikU −1 e−θU uikvkje,φ(U, V ) =
Γ(θ
)
Γ(θ
)
U
V
j=1
i=1
p(X|U, V ) =

k=1

k=1

where ϕU , θU , ϕV , θV > 0 are hyperparameters, the size of U and V are M × H and H × N ,and (U V )ij is the (i, j) entry of U V , respectively.
Then, the variational free energy F̄nvb satisfies the following asymptotic equality:F̄nvb = nSn + λvb log n + Op (1) (n → ∞),where{(H − H0 )(M ϕU + N ϕV ) + 12 H0 (M + N ), if M ϕU + N ϕV < M +N
2
λvb = 1H(M+N),
otherwise.
2
In general, the learning coeﬃcient of VB may not be equal to but becomes an upper boundof the RLCT: λvb ≧ λ, since the variational free energy is larger than the free energy even
if the sample size diverges infinity. Unfortunately, the variational generalization error is notequal to λvb /n, asymptotically. Besides, as described below, the variational inference seeks
the mode of the true posterior and the variational posterior distributions tend to concentrate ina fraction of true posterior distributions. One might say that a variational posterior distribution
is approximately equivalent to a posterior distribution; however, for these reasons, variational
inference is just an approximation of Bayesian inference and they are not equivalent. As far as
we know, there has been no direct theoretical comparison between them, except for [92]. The

4.2 Main Theorem

33

diﬀerence between VB and the local variational approximation, which is an approximation of
VB, has been studied theoretically [76]; however, for Bayesian inference and VB, there exists
few theoretical comparison.
VBNMF has been devised [18], and the exact learning coeﬃcient of VBNMF has been
derived [49]. Nevertheless, the variational approximation error has not been clarified since
the RLCT of NMF has been unknown. If the prior distribution is strictly and entirely positiveand bounded analytic function on the domain, then an upper bound of the RLCT of NMF has
been proved [37, 36]. If the non-negative restriction is not assumed for matrix factorization,then the exact value of the RLCT has been clarified as those of reduced rank regression
models [10]. However, the RLCT has been unknown in the case of that the prior is a gamma
distribution, which may be zero.
In this chapter, we give an upper bound of the RLCT of NMF λ when the prior is a gamma
distribution, which determines an upper bound of the Bayesian generalization error and thefree energy in that case. Moreover, by comparing λ with λvb , we also derive a lower boundof the variational approximation error in VBNMF.

4.2 Main Theorem
Let M(M, N, C) be a set of M × N matrices whose elements are in C, where C is a subsetof R. Let K be a compact subset of R≧0 = {x ∈ R|x ≧ 0} and let K0 be a compact
of subset R>0 = {x ∈ R|x > 0}. We denote that U ∈ M(M, H, K), V ∈ M(H, N, K)
and U0 ∈ M(M, H0 , K0 ), V0 ∈ M(H0 , N, K0 ) are the NMF result of U0 V0 such that H0
is the non-negative rank [19] of U0 V0 , i.e. they give the minimal H0 , where H ≧ H0 and{(x, y, a, b) ∈ K × K0 |xy = ab} ̸= ∅.
Definition 4.1 (RLCT of NMF) Assume that the largest pole of the function of one complex
variable z,∫∫()zζ(z) =dU
dV ∥U V − U0 V0 ∥2 φ(U, V )M(M,H,K)

M(H,N,K)

is equal to (−λ). Then λ is said to be the RLCT of NMF.
According to one of the first Bayesian NMF paper [18], we set the prior as gamma distributions:φ(U, V ) = Gam(U |ϕU , θU )Gam(V |ϕV , θU ),
where
M ∏
H
ϕU
∏
θU
Gam(U |ϕU , θU ) =uϕikU −1 e−θU uik ,Γ(θU )i=1
k=1

Gam(V |ϕV , θV ) =

H ∏
N
∏
k=1

θVϕV ϕV −1 −θV vik
vik e
,
Γ(θ
)
V
j=1

and ϕU , θU , ϕV , θV > 0.
In this chapter, we prove the following theorems.

Chapter 4 Bayesian Generalization Error in Non-negative Matrix Factorization

34

Theorem 4.2 If the prior is the above gamma distributions, then the RLCT of NMF λ satisfies
the following inequality:λ≦1[(H − H0 ) min{M ϕU , N ϕV } + H0 (M + N − 2) + δH0 ] ,2where δH0 = [H0 ≡ 1(mod 2)]. The equality holds when H = H0 = 1 or 2.
We prove this theorem in the next section. As two applications of this theorem, we obtainan upper bounds of the free energy and Bayesian generalization error of NMF in this case.
The following theorem shows a statistical bound of Bayesian estimation of NMF.
Theorem 4.3 Let the probability density functions of X ∈ M(M, N, K) be q(X) andp(X|U, V ), which represent a true distribution and a learning machine respectively definedbyq(X) = Poi(X|U0 V0 ),p(X|U, V ) = Poi(X|U V ),where
Poi(X|A) =N
M ∏∏(aij )xij
i=1 j=1

xij !

M,Ne−aij , X = (xij )M,N
i=1,j=1 , A = (aij )i=1,j=1 .

Also let φ(U, V ) = Gam(U |ϕU , θU )Gam(V |ϕV , θU ). Then, the normalized free energy Fn
and the expected generalization error E[Gn ] satisfies the following inequality:
1[(H−H0 ) min{M ϕU , N ϕV } + H0 (M +N −2) + δH0 ] log n + Op (1),
2
( )
1
1
E[Gn ] ≦[(H−H0 ) min{M ϕU , N ϕV } + H0 (M +N −2) + δH0 ] + o
.
2n
n
Fn ≦

In Theorem 4.3, we study a case when a set of random matrices X n = X1 , X2 , ..., Xn
are observed and the true decomposition U0 and V0 are statistically estimated. Actuallysometimes NMF has studied in the case when only one target matrix is decomposed, however,in general, decomposition of a set of independent matrices should be studied because target
matrices are often obtained daily, monthly, or diﬀerent places for purchase study [47]. Other
situation for observing multiple data matrices is traﬃc data analysis [21]. In such cases,
decomposition of a set of matrices results in statistical inference. We consider this situationand it is common to [49]. A statistical model p(X|U, V ) which has parameters (U, V ) are
employed for estimation. Then the free energy and generalization error of Bayesian estimation
is given by this theorem.
Moreover, we also derive a lower bound of the variational approximation error in VBNMF as
Theorem 4.4 in the following. In VB, the posterior ψ(θ|X n ) is approximated by the variational
posterior ψ vb (θ|X n ) which is a distribution of independent parameters like ψ vb (θ|X n ) :=
ψ1 (θ1 |X n )ψ2 (θ2 |X n ). The variational approximation error is defined by the KL divergence4.3 Preparation

35

from the variational posterior to the true one: KL(ψ vb ∥ψ). Developing it, we have∫
ψ vb (θ|X n )vb
KL(ψ ∥ψ) = dθψ vb (θ|X n ) log(4.1)
ψ(θ|X n )
∏n
(
)
∫
vb
n
vb
n
i=1 p(Xi |θ)φ(θ)(4.2)= dθψ (θ|X ) log ψ (θ|X ) − log
Zn
(
)
∫
n
∏
= dθψ vb (θ|X n ) log ψ vb (θ|X n ) − logp(Xi |θ)φ(θ) + log Zn
i=1

(4.3)

∫
=

ψ vb (θ|X n )− Fn .
i=1 p(Xi |θ)φ(θ)

dθψ vb (θ|X n ) log ∏n(4.4)

The first term in Eq. (4.4) is called the variational free energy F̄nvb . Because of KL(ψ vb ∥ψ) ≧0, F̄nvb ≧ F n . It is the objective function of VB. If θ1 and θ2 such that θ = (θ1 , θ2 ) wereindependent, then the variational posterior would be same as the true posterior and thevariational approximation error would become zero. However, in general, θ1 and θ2 such that
θ = (θ1 , θ2 ) are not independent. Thus, VB minimizes the variational free energy.
Theorem 4.4 Let the variational free energy of VBNMF be F̄nvb . Then, the following inequality is attained:F̄nvb − F n ≧ λ log n + Op (1),
where
{
λ=

1
2 [(H−H0 )(M ϕU +N ϕV + max{M ϕU , N ϕV })−δH0 ]+H0 ,1
2 [(H−H0 )(M +N − min{M ϕU , N ϕV })−δH0 ]+H0 ,if M ϕU +N ϕV < M +N2otherwise.

The framework of VB and the definition of the variational free energy are described inAppendix. Theorem 4.4 gives a lower bound of the diﬀerence of the free energy between thevariational approximation and the true.

4.3 Preparation
In order to prove Theorem 4.2, we use the following five lemmas.
Lemma 4.1 Let q(X) and p(X|U, V ) be the probability density functions of a non-negative
matrix X, which represent a true distribution and a learning machine respectively defined byq(X) = Poi(X|U0 V0 ),p(X|U, V ) = Poi(X|U V ),where Poi(X|W ) is a probability density function of the Poisson distribution with average∑q(X)W . Then, the RLCT of the KL divergence X q(X) log p(X|U,V) is equal to the RLCT of
NMF in Definition 4.1:
∥U V − U0 V0 ∥2 ∼∑X

q(X) logq(X).
p(X|U, V )Chapter 4 Bayesian Generalization Error in Non-negative Matrix Factorization

36

Lemma 4.2 Let λ be the absolute value of the maximum pole of∫∫)z(ζ(z) =dU dV ∥U V ∥2 Gam(U |ϕU , θU )Gam(V |ϕV , θV ).
Then,

H min{M ϕU , N ϕV }2holds; this is the equality of Theorem 4.2 in the case H0 = 0.
λ=Lemma 4.3 If H0 = H = 1, the equal sign of the Theorem4.2 holds:
M +N −1.
2
Lemma 4.4 If H0 = H = 2, the equal sign of the Theorem4.2 holds:
λ=λ = M + N − 2.
Lemma 4.5 If H = H0 , the Theorem4.2 is attained:
H0 (M + N − 2) + [H0 ≡ 1(mod 2)].
2
Let the entries of the matrices (U, V ) beλ≦U = (u1 , . . . , uH ), uk = (uik )M
i=1 ,
V = (v1 , . . . , vH )T , vk = (vkj )N
j=1 ,and the ones of (U0 , V0 ) be0
U0 = (u01 , . . . , uH), u0k = (u0ik )M
i=1 ,
0
0 T
0 N
V0 = (v10 , . . . , vH
) , vk0 = (vkj)j=1 ,
0respectively.
In the following, these lemmas are proved.
0
Proof of Lemma 4.1. Let (wij ) = W = U V and (wij
) = W0 = U0 V0 . Calculating the KL∑q(X)divergence X q(X) log p(X|U,V ) , we have∑q(X) log

X

=

=

=

q(X)
p(X|U, V )
∏M ∏N

0 xij
0
(wij
)
−wij
i=1
j=1
xij ! eq(X) log ∏M ∏N (w )xij
ij
−wij
X
i=1
j=1
xij ! e

∑
∑

q(X)

M ∑
N
∑

X

i=1 j=1

∑

M ∑
N
∑

X

q(X)

i=1 j=1

(4.5)

(4.6)

0
0
{xij log wij− log xij ! − wij− (xij log wij − log xij ! − wij )} (4.7)

{

}
0
wij
0
xij log− (wij− wij ) .
wij

(4.8)

4.3 Preparation
Because of

37
∞
0 xij
∑
(wij
)
xij =0

xij !

0
−wij

e

{
0
wij
x i′ j ′ =x i′ j ′(i, j) = (i′ , j ′ ),(i, j) ̸= (i′ , j ′ ),(4.9)

we have∑q(X)p(X|U, V )
X
}
{
M ∑
N
0
∑
∑
wij
0
=
q(X)
xij log− (wij− wij )w
ij
i=1 j=1
X

{
}
M
N
M
N
0
x
0
ij
∑ ∑ ∑ ∏ ∏ (wij )
w
0
ij
0

=
e−wij  xij log− (wij− wij )x!
w
ij
ij
i=1 j=1i=1 j=1
q(X) log(4.10)(4.11)(4.12)

X

}
0
wij
0
e
xij log=− (wij− wij )x!
w
ij
ij
X i=1 j=1
{
}
M ∑
N ∑
∞
0 xij
0
∑
(wij
)
wij
0
−wij
0
=
e
xij log− (wij − wij )
xij !
wij
i=1 j=1 xij =0
{
}
M ∑
N
0
∑
w
ij
0
0
=
wij
log
− (wij− wij ) .
w
ij
i=1 j=1
M ∑
N
0 xij
∑∑
(wij
)

{

0
−wij

(4.13)(4.14)(4.15)All we have to prove is00(wij− wij )2 ∼ wij
log

0
wij
0
− (wij− wij )
wij

for i = 1, . . . , M and j = 1, . . . , N , since ∥W − W0 ∥2 =Let f (a, b) = (a − b)2 andg(a, b) = a log(4.16)

∑M ∑N
i=1

a
− (a − b)b02
j=1 (wij − wij ) .

(4.17)on (a, b) ∈ (R≧0 )2 . Owing to∂a g(a, b) = loga − logb,∂b g(a, b) = 1 − a/b,
and that a log function is monotone increasing, we obtain∂a g(a, b) = ∂b g(a, b) = 0 ⇔ a = b.

(4.18)

38

Chapter 4 Bayesian Generalization Error in Non-negative Matrix Factorization
Signs of the above partial derivations are∂a g(a, b) > 0 ∧ ∂b g(a, b) < 0 in case of a > b,

(4.19)∂a g(a, b) < 0 ∧ ∂b g(a, b) > 0 in case of a < b.

(4.20)On account of the above and smoothness, the increase or decrease and convexity of g(a, b)
is the same as those of f (a, b). Hence there exists c1 , c2 > 0 such thatc1 f (a, b) ≦ g(a, b) ≦ c2 f (a, b).

(4.21)

i.e. K(a, b) ∼ (b − a)2 . Therefore, we have
∥U V − U0 V0 ∥2 ∼∑q(X) log

X

q(X).
p(X|U, V )Put Φ(U, V ) = ∥U V ∥ . The prior φ(U, V ) is gamma distributions, hence,(M H) H N∪ ∪∪ ∪φ−1 (0) ={(U, V ) | uik = 0} ∪ {(U, V ) | vkj = 0} .

(4.22)□

2

i=1 k=1

k=1 j=1

If Φ−1 (0) ̸= ∅, φ−1 (0) ̸= ∅, and Φ−1 (0) ∩ φ−1 (0) = ∅, then φ(U.V ) does not correspondto the maximum pole of the zeta function. However, in the case Lemma 4.2, the set Φ−1 (0)has intersections with φ−1 (0). Here, we prove Lemma 4.2.
Proof of Lemma 4.2. We consider simultaneous resolution ∥U V ∥2 = 0 and φ(U, V ) = 0
since {(U, V ) | ∥U V ∥2 = 0} has intersections with φ−1 (0).
The zeta function is equal to∫∫()zζ(z) =dU dV ∥U V ∥2 Gam(U |ϕU , θU )Gam(V |ϕV , θV ).
(4.23)
Since elements of matrices is nonnegative, the inequalityH∑2u2ik vkj
≦

(H
∑

)2
≦H

uik vkj

k=1

k=1

H
∑

2
u2ik vkj(4.24)

k=1

holds. Thus, we have
M ∑
N ∑
H
∑

2
u2ik vkj≦ ∥U V ∥2 ≦ H

i=1 j=1 k=1

N ∑
H
M ∑
∑

2
u2ik vkj.

(4.25)

i=1 j=1 k=1

As a log canonical threshold is not changed by any constant coeﬃcient and it is orderisomorphic; ∃(c1 , c2 ) ∈ R2 s.t. c1 F ≦ G ≦ c2 F ⇒ F ∼ G, all we have to do is calculating
the RLCT of
(M
)( N
)
M ∑
N ∑
H
H
∑
∑
∑
∑
2 2
2
2
uik vkj =uikvkj .
(4.26)
i=1 j=1 k=1

k=1

i=1

j=1

4.3 Preparation

39

The RLCT λ becomes a sum of each ones about k. For each k, we blow-up ofvariables (uik ), (vkj ) each other. Simultaneously, we also consider blowing-ups for∏H ∏M ϕU −1 ∏N ϕV −1) and the determinant |J| of their Jacobi matrix J. |J| is
k=1 ( i=1 uik
j=1 vkj
called Jacobian.
Let A be either U or V and brc [A] : (U, V ) 7→ (U ′ .V ′ ) be a blowing-up such that the′′(i, k)-entry of U ′ and the (k, j)-entry of V ′ are equal to Uik
and Vkj
:

{

uik

′


Uik = u u(i = r and k = c) or k ̸= c
′
, Vkj
= vkj
i
=
̸
r
and
k
=
c
rc ik
{

vkj
k = r and j = c

′
′


Uik = uik , Vkj = v votherwiserc kj

A = U,.

(4.27)A=V

For example, let M = N = H = 2 and apply b11 [U ] to
∥U V ∥2 ∼

(M
H
∑
∑
k=1

Then, we have

(

′

U =

u11
u11 u21

u2ik

)( N
∑

i=1

)
u12
,
u22

)
2
vkj
.

(4.28)

j=1

′

V =

(

v11
v21

v12
v22

)
(4.29)

and
( 2
2
∑
∑
k=1

i=1

u2ik

)( 2
∑

)
2
vkj

2
2
2
2
7→ u211 (1 + u221 )(v11+ v12) + (u212 + u222 )(v21+ v22), (4.30)

j=1

2
2 ∏
2
∏
∏
2(ϕ −1)
ϕV −1
ϕV −1 ϕV −1 ϕU −1 ϕU −1 ϕV −1 ϕV −1
ϕU −1
( uik
vkj
) 7→ u11 U uϕ21U −1 v11v12 u12 u22 v21 v22 .
k=1 i=1

j=1(4.31)The Jacobi matrix is an 8 × 8 matrix as the below:
1 u12 u21 u22
0 1
0
0

0 0 u110′′∂(U , V ) 0 0
0
1
=
J=
0
0
0
0
∂(U, V )

0 0
0
0

0 0
0
0
0 0
0
0
and its determinant is |J| = (u11 )1 .

v11
0
0
0
1
0
0
0

v12
0
0
0
0
1
0
0

v21
0
0
0
0
0
1
0


v22
0 

0 

0 

0 

0 

0 
1

(4.32)

40

Chapter 4 Bayesian Generalization Error in Non-negative Matrix Factorization
We arbitrarily take k ∈ {1, . . . , H} and fix it. In the general case, because of the symmetryof the variables, only the blowing-ups b1k [U ] and bk1 [V ] should be treated. Hence, we have

(M
)( N
)
(
)
M
N
∑
∑
∑
∑
2
2
2 
u2ik
vkj
7→ u21k vk1
1+
u2ik 1 +vkj,(4.33)
i=1
M
∏

j=1

uϕikU −1

i=1

N
∏

i=2

ϕV −1
ϕU −M N ϕV −N
vkj
7→ uM
vk1
1k

j=1

j=2
M
∏

uϕikU −1i=2

−1 N −1
The Jacobian |J|k is equal to |J|k = uM
vk1 . The term
1k

(

N
∏

ϕV −1
vkj
.

(4.34)

j=2∑M2
i=2 uik)(∑N2
j=2 vkj)

1+
1+
is strictly positive; thus, we should consider the maximum pole of the following meromorphic
function:


∫∫
M
N
∏
∏
ϕV −1 2 z  M ϕU −M N ϕV −Nζ̃(z) =dU dV (u21k vk1
)
u1k
vk1
uϕikU −1
vkj
|J|k (4.35)
∫∫
=C
∫∫
=C
=D

i=2

j=2

2z M ϕU −M N ϕV −N M −1 N −1
du1k dvk1 u2z
vk1
u1k vk1
1k vk1 u1k(4.36)

ϕU −1 2z+N ϕV −1
du1k dvk1 u2z+M
vk1
1k(4.37)

1
1
×,2z + M ϕU
2z + N ϕV(4.38)

where C and D are positive constants. The poles are z = −M ϕU /2, −N ϕV /2; therefore,
the RLCT λk ismin{M ϕU , N ϕV }
λk =
.
2
(
)(
)
∑
∑
M
N
N
2
2
Let Uk be (uik )Mfor
i=1 , Vk be (vkj )j=1 , and Φk (Uk , Vk ) be
i=1 uik
j=1 vkjeach k. The RLCT λ becomes a sum of each ones about k;
λ=

H
∑

λk ,

(4.39)

k=1

since Φ(U, V ) ∼∑H

k=1 Φk (Uk , Vk ) and the RLCT of Φk is λk . Thus, we haveλ=H min{M ϕU , N ϕV }.
2

(4.40)□

Remark 4.1 Lemma 4.2 means that the equality in Main Theorem holds if U0 V0 = O, i.e.
H0 = 0.
Next, the exact values of the RLCTs in the case H = H0 = 1 and H = H0 = 2 (Lemma
4.3 and 4.4) are derived. To prove Lemma 4.3 and 4.4, we show the following lemma.

4.3 Preparation

41

Lemma 4.6 The RLCT of reduced rank regression λMF ,which is equal to the RLCT of matrix
factorization, is a lower bound of the RLCT of NMF λ, i.e.
λMF ≦ λ.

(4.41)

Proof of Lemma 4.6. Let A, B, A0 and B0 be M × H, H × N , M × r and r × N matrices,
respectively. The RLCT of matrix factorization is defined by the absolute maximum pole ofthe zeta function∫∫
ζ(z) =∥AB − A0 B0 ∥2z φMF (A, B)dAdB,(4.42)
where φMF is a prior distribution of (A, B) whose support includes negative real numbers (seealso Definition 4.2 in the below). On the other hand, φ is a prior distribution of non-negative
matrices (U, V ); thus, there exists c > 0 such that
φ(A, B) ≦ cφMF (A, B).

(4.43)

Because of the property of the maximum pole of the zeta function [81], we obtain λMF ≦ λ.
□Remark 4.2 As compared to the RLCT of reduced rank regression, the domain of the priorof NMF is considered that includes negative real numbers as 0 density:{φ(U, V ) A = U ≧ 0 ∧ B = V ≧ 0,
φ(A, B) =(4.44)
0
A or B has negative entries.
The RLCT of reduced rank regression was clarified in the all case [10] as the follwoing
theorem.
Theorem 4.5 (Aoyagi) The RLCT of reduced rank regression λMF is as follows:(1) If N + r ≦ M + H ∧ M + r ≦ N + H ∧ H + r ≦ M + N ,(1−1) in the case M + H + N + r is even,λMF =−(H +r)2 −M 2 −N 2 + 2{(H +r)(M +N )+M N }.
8(1−2) in the case M + H + N + r is odd,
λMF =1−(H +r)2 −M 2 −N 2 + 2{(H +r)(M +N )+M N }.
8(2) Else if M + H < N + r,
λMF =HM − Hr + N r
.
2

λMF =

HN − Hr + M r.
2(3) Else if N + H < M + r,

42

Chapter 4 Bayesian Generalization Error in Non-negative Matrix Factorization(4) Or else, i.e. in the case M + N < H + r,
λMF =

MN
.
2

Especially, if H = r, thenλMF =H(M + N − H).
2

In the case (1−2), the multiplicity is two: mMF = 2. Otherwise, it equals one: mMF = 1.
By using Theorem 4.5 and Lemma 4.6, the following inequalities hold.
Corollary 4.1 If H0 = H = 1, thenM +N −1≦ λ.
2Besides, if H0 = H = 2, thenM + N − 2 ≦ λ.

In the case H = H0 , we can treat the entries of U and V as positive numbers since theentries of U0 and V0 are positive for any (U0 , V0 ) ∈ Φ−1 (0) and there exists a neighborhood
such that its all elements are positive matrices owing to continuity of real numbers. Thus,in this case, the prior φ(U, V ) is positive and bounded; it does not aﬀect the RLCT. Then,because of Corollary 4.1, we only have to prove inversed inequalities of the aboves in orderto derive Lemma 4.3 and 4.4.
Put Φ(U, V ) = ∥U V − U0 V0 ∥2 . In general, if Φ−1 (0) has an R-dimensional sub-manifold,then the rank of its Hesse matrix ∇2 Φ is at most H(M + N ) − R. According to [81] (byusing implicit function theorem), we haveλ≦H(M + N ) − R.
2(4.45)

Proof of Lemma 4.3. Let H = H0 = 1 and p > 0. Because we consider the situation
Φ(U, V ) = ∥U V − U0 V0 ∥2 = 0 and U and V are vectors, we can treat the elements in U and
V as strictly positive. Hence, the entries of pU and p−1 V are also strictly positive. We put
P := {p ∈ R|p > 0} and it is a 1-dimensional sub-manifold of Φ−1 (0); thus, by using theinequality (4.45), we have
M +N −1λ≦.
(4.46)
2
Owing to Corollary 4.1, we obtain
λ=M +N −1.
2(4.47)□

Remark 4.3 The author initially derived Lemma 4.3 in a diﬀerent way: by using degeneratingof ideal and mathematical induction about M and N . The initial proof of Lemma 4.3 isdescribed in Appendix.

4.3 Preparation

43

In a way similar to Proof of Lemma 4.3, we prove Lemma 4.4 in the following.
Proof of Lemma 4.4. Let H = H0 = 2 and arbitrarily take a parameter matrix pair (U, V ) ina suﬃcient small neighborhood N of Φ−1 (0). A set of 2 × 2 real regular matrices is denotedby GL(2, R). In N , we can treat the entries of U and V as positive numbers. Let
P := {P ∈ GL(2, R)|U P > 0, P −1 V > 0}.

(4.48)

P is a sub-manifold of Φ−1 (0); thus, all we have to do is prove that P is a 4-dimensionalvariety. P ∈ P is denoted by
(
)
p q
P :=
(4.49)
r s
and we put ∆ := det P . Then, we immediately get
(
)
1
s −q
−1
P =
.
−r p
∆

(4.50)

Hence, we have
u11 p + u12 r
u11 q + u12 s


..
..
UP = 
,.
.
uM 1 p + uM 2 r uM 1 q + uM 2 s()
1−v21 q + v11 s . . . −v2N q + v1N s
−1
P V =
.
v21 p − v11 r . . . v2N p − v1N r∆(4.51)(4.52)

Without loss of generality, we can consider the case ∆ > 0 because we can derive thelemma in the case ∆ < 0 by using inversed inequalities. We can represent the condition
U P > 0 and P −1 V > 0 as the following simultaneous inequalities:
u11 p + u12 r > 0..
(4.53).
u p + u r > 0,
M1
M2



u11 q + u12 s > 0..
(4.54).

u q + u s > 0,
M1
M2



−v21 q + v11 s > 0..
(4.55).

−v q + v s > 0,
2N
1N



v21 p − v11 r > 0..
(4.56).

v p − v r > 0,2N
1N

Chapter 4 Bayesian Generalization Error in Non-negative Matrix Factorization

44

for uik > 0 and vkj > 0 (i = 1, . . . , M ; k = 1, 2; j = 1, . . . , N ). Therefore, (p, q, r, s) must
be in upper domains of lines in rp-plain and qs-plain. The domains in rp-plain and qs-plain
are 2-dimensional manifold, respectively. The Cartesian product of these is also a manifold
and its dimension is 2 + 2 = 4. Hence, we haveλ≦2(M + N ) − 4= M + N − 2.
2

(4.57)

By applying Corollary 4.1, Lemma 4.4 is proved:λ=2(M + N ) − 4= M + N − 2.
2

(4.58)□As the last part of this section, we drive an upper bound of the RLCT of NMF when
H = H0 . Put
u1(2k−1) u1(2k)


..
..
Uk := ,.
.
uM (2k−1)
(

)

...
...

v(2k−1)NvH(2k)

u01(2k−1)
..
Uk0 := .
0uM (2k−1)

u01(2k)..
,.
0uM (2k)Vk :=v(2k−1)1vH(2k−1)

uM (2k)(
Vk0 :=0v(2k−1)10vH(2k−1)

...
...

0
v(2k−1)N0vH(2k),).

Proof of Lemma 4.5. We divide this proof into two cases: whether H = H0 is even.
Case (1): H=H0 ≡0 (mod 2), i.e. there exists H ′ ∈N such that H=H0 =2H ′ .
There is a positive constant C > 0 such that∥U V − U0 V0 ∥2
=

M ∑
N
∑

0
0
(ui1 v1j + ... + uiH vHj − u0i1 v1j− u0iH vHj)2

i=1 j=1

=

(H
M ∑
N
∑
∑
i=1 j=1(4.59)(4.60))20(uik vkj − u0ik vkj)(4.61)

k=1
′

≦C

M ∑
N ∑
H
∑

0
0
(ui(2k−1) v(2k−1)j +ui(2k) v(2k)j −u0i(2k−1) v(2k−1)j−u0i(2k) v(2k)j)2 =: Φ1 .

i=1 j=1 k=1(4.62)

4.3 Preparation

45

Thus, we have
′

Φ1 ∼

M ∑
N ∑
H
∑

0
0
(ui(2k−1) v(2k−1)j + ui(2k) v(2k)j − u0i(2k−1) v(2k−1)j− u0i(2k) v(2k)j)2i=1 j=1 k=1

(M N
H
∑
∑∑

(4.63)
)

′

=

k=1

0
0
)2
−u0i(2k) v(2k)j(ui(2k−1) v(2k−1)j +ui(2k) v(2k)j −u0i(2k−1) v(2k−1)j

i=1 j=1(4.64)′=H∑

Uk Uk − Uk0 Vk02=: Φ2 .

(4.65)

k=1
2

Let λk be the RLCT of Uk Uk − Uk0 Vk0 for k = 1, . . . , H ′ . The RLCT of Φ2 is denotedby λ. As all intersections between {(Uk , Vk )} and {(Uk′ , Vk′ )} are empty sets if k ̸= k ′ ,
′

λ=

H
∑

λk

(4.66)

k=1

holds. For each k, because of Lemma 4.4, we have λk = M + N − 2.
Thus, on account of H ′ = H0 /2, we obtain′H∑(M + N − 2)

(4.67)= H (M + N − 2)
M +N −2= H0.
2(4.68)

λ=k=1′(4.69)

Case (2): H=H0 ≡1 (mod 2), i.e. there exists H ′ ∈N such that H=H0 =2H ′ −1.
In the same way as Case (1), there exists a positive constant C ′ > 0 such that∥U V − U0 V0 ∥2
(H
)2
N
M ∑
∑
∑
0
(uik vkj − u0ik vkj)=
i=1 j=1≦ C′(4.70)(4.71)

k=1

M ∑
N
∑

(H ′ −1∑i=1 j=1k=100(ui(2k−1) v(2k−1)j + ui(2k) v(2k)j − u0i(2k−1) v(2k−1)j− u0i(2k) v(2k)j
)2

)
0
2
+ (ui(2H ′ −1) v(2H ′ −1)j − u0i(2H ′ −1) v(2H
′ −1)j )(4.72)(4.73)

Chapter 4 Bayesian Generalization Error in Non-negative Matrix Factorization

46
∼

M ∑
N
∑

(H ′ −1∑i=1 j=1k=100(ui(2k−1) v(2k−1)j + ui(2k) v(2k)j − u0i(2k−1) v(2k−1)j− u0i(2k) v(2k)j)2(4.74))02+ (ui(2H ′ −1) v(2H ′ −1)j − u0i(2H ′ −1) v(2H
′ −1)j )

=

′
H
−1
∑

(M N
∑∑

k=1

i=1 j=1(4.75))00(ui(2k−1) v(2k−1)j + ui(2k) v(2k)j − u0i(2k−1) v(2k−1)j− u0i(2k) v(2k)j)2(4.76)
+

M ∑
N
∑

0
2
(ui(2H ′ −1) v(2H ′ −1)j − u0i(2H ′ −1) v(2H
′ −1)j )(4.77)

i=1 j=1

=

′
H
−1
∑

Uk Uk − Uk0 Vk0

2

′

0
T 2
(4.78)+ ∥u2H ′ −1 (v2H ′ −1 )T − u02H ′ −1 (v2H
′ −1 ) ∥ =: Φ .

k=1∑H ′ −1′′2
Let λ1 and λ2 be the RLCT of k=1 Uk Uk − Uk0 Vk0 and ∥u2H ′ −1 (v2H ′ −1 )T −T 2
0
u02H ′ −1 (v2H
′ −1 ) ∥ , respectively. The former is calculated by Lemma 4.4 as same as in Case
(1) and the latter is derived by Lemma 4.3. Hence, we have
′
λ1 =

′
H
−1
∑

(M + N − 2) = (H ′ − 1)(M + N − 2),(4.79)

k=1
′

λ2 =

M +N −1.
2(4.80)′′The RLCT of Φ is denoted by λ . By the definition, H ′ = (H0 + 1)/2 holds. Since all
combinations of intersections in ({(Uk , Vk )}, {(Uk′ , Vk′ )}, {(u2H ′ −1 , v2H ′ −1 )}) are empty
sets if k ̸= k ′ , we have′

M +N −1
2(H0 − 1)(M + N − 2) + (M + N − 1)=2
(H0 − 1)(M + N − 2) + (M + N − 2) + 1=
2
H0 (M + N − 2) + 1=.
2

λ = (H ′ − 1)(M + N − 2) +(4.81)(4.82)(4.83)(4.84)′The RLCT is order isomorphic: λ ≦ λ[H0 ≡ 0 (mod 2)]+λ [H0 ≡ 1 (mod 2)]. Therefore,summarizing Case (1) and (2), we obtainλ≦1{H0 (M + N − 2) + [H0 ≡ 1 (mod 2)]} .
2(4.85)□

4.4 Proof of Main Theorem

47

4.4 Proof of Main Theorem
Based on the above preparation, Theorem 4.2, 4.3 and 4.4 are proved.
Proof of Theorem 4.2. Put δH0 := [H0 ≡ 1(mod 2)]. There exists a positive constant D > 0such that∥U V − U0 V0 ∥2
=

M ∑
N
∑

(4.86)

00(ui1 v1j + . . . + uiH vHj − u0i1 v1j− . . . − u0iH0 vH)20j

i=1 j=1

=

(H
M ∑
N
∑
∑
i=1 j=1

=

k=1

(H
M ∑
N
0
∑
∑
i=1 j=1

uik vkj −

H0
∑

)2
0
u0ik vkj

0(uik vkj − u0ik vkj
)+

H
∑

)2
uik vkj(4.89)0 T 2+δH0 D∥uH0 (vH0 )T −u0H0 (vH) ∥0(4.90)

k=1∑(4.88)

k=1

k=H0 +1(H0 −δH0 )/2≦D(4.87)Uk Vk −Uk0 Vk0

2

k=1




2
u1(H0 +1) . . . u1Hv(H0 +1)1 . . . v(H0 +1)N


..
.. 
..
..
..
..
+D 
 .
.
.
.
. 
.
.
uM (H0 +1) . . . uM H
vH1
...
vHN

(4.91)Let K 1 , K 2 and K 3 be
(H0 −δH0 )/2

K1 =

∑

Uk Vk −Uk0 Vk02,(4.92)k=1
0 T 2
K 2 = ∥uH0 (vH0 )T −u0H0 (vH
) ∥ ,
0


2
u1(H0 +1) . . . u1Hv(H0 +1)1 . . . v(H0 +1)N


..
.. 
..
..
..
..
K3 = 
 ,
.
.
.
. 
.
.

uM (H0 +1) . . . uM H

vH1

...

(4.93)(4.94)

vHN

respectively. The RLCT of K 1 , K 2 and K 3 are respectively denoted by λ1 , λ2 and λ3 . The
constant D does not aﬀect the RLCTs; thus, we only have to consider λ1 , λ2 and λ3 .
According to Lemma 4.5, we have
λ1 + λ2 =1{H0 (M + N − 2) + δH0 }.
2(4.95)

Moreover, as K 3 is equivalent to ∥U V ∥2 replaced from H to H − H0 in Lemma 4.2, we
obtain1(4.96)
λ3 = (H − H0 ) min{M ϕU , N ϕV }.
2Chapter 4 Bayesian Generalization Error in Non-negative Matrix Factorization48

Summarizing the aboves, the RLCT of NMF λ satisfies
λ ≦ λ1 + λ2 + λ3
1
1
= {H0 (M + N − 2) + δH0 } + (H − H0 ) min{M ϕU , N ϕV }2
2
1
= {(H − H0 ) min{M ϕU , N ϕV } + H0 (M + N − 2) + δH0 } .
2(4.97)(4.98)(4.99)

Therefore, Theorem 4.2 is proved:
λ≦1{(H − H0 ) min{M ϕU , N ϕV } + H0 (M + N − 2) + δH0 } .
2(4.100)□

Proof of Theorem 4.3. Let λ be the upper bound in Theorem 4.2:
λ=1{(H − H0 ) min{M ϕU , N ϕV } + H0 (M + N − 2) + δH0 } .
2

By using Theorem 3.4 and Theorem 4.2 λ ≦ λ, we obtain
(
)
λ
m−1
1
E[Gn ] = −+o,n n log nn log n
( )
1
λ
≦ +o
n
n

(4.101)(4.102)(4.103)andFn = λ log n − (m − 1) log log n + Op (1)≦ λ log n + Op (1).

(4.104)(4.105)Thus, Theorem 4.3 is derived.
□
Proof of Theorem 4.4. Let Sn be the empirical entropy. By using Proposition 2.2 and Theorem 4.3, we have
Fn = nSn + Fn

(4.106)≦ nSn + λ log n + Op (1).

(4.107)F̄nvb = nSn + λvb log n + Op (1)(4.108)

Also, because of Theorem 4.1,

holds, where{(H − H0 )(M ϕU + N ϕV ) + 21 H0 (M + N ), M ϕU + N ϕV < M +N
2
λvb = 1
(4.109)
M +N
H(M
+
N
),
M
ϕ
+
N
ϕ
≧
U
V
2
2 .

4.5 Discussion

49

Thus, we compute their diﬀerenceF̄nvb − F n = (λvb − λ) log n + Op (1)(4.110)≧ (λvb − λ) log n + Op (1).

(4.111)

Case(1): M ϕU + N ϕV < M +N
holds. We have2
1
λvb − λ = (H − H0 )(M ϕU + N ϕV ) + H0 (M + N )(4.112)
2
1− [(H − H0 ) min{M ϕU , N ϕV } + H0 (M + N − 2) + δH0 ](4.113)2[]
1= (H − H0 ) M ϕU + N ϕV − min{M ϕU , N ϕV }(4.114)2
1
δH(4.115)+ H0 (M + N − M − N + 2) − 0
221= [(H−H0 )(M ϕU +N ϕV +M ϕU +N ϕV − min{M ϕU , N ϕV })+δH0 ]+H02
(4.116)1= [(H − H0 )(M ϕU + N ϕV + max{M ϕU , N ϕV }) + δH0 ] + H0 .
(4.117)2
holds. We haveCase(2): M ϕU + N ϕV ≧ M +N
2
1
1
H(M +N )− [(H−H0 ) min{M ϕU , N ϕV }+H0 (M +N −2)+δH0 ] (4.118)
2
2
1
1
1
δH
= H(M +N )− (H−H0 ) min{M ϕU , N ϕV }− H0 (M +N )− 0 +H02
2
2
2(4.119)λvb −λ =11δH(H − H0 )(M + N ) − (H − H0 ) min{M ϕU , N ϕV } − 0 + H0 (4.120)
222
1= [(H − H0 )(M + N − min{M ϕU , N ϕV }) + δH0 ] + H0 .
(4.121)
2=

Therefore, we obtain Theorem 4.4.
□

4.5 Discussion
Here, we will discuss the results of this chapter from four points of view. After that, we will
describe the numerical behavior of the theoretical result by conducting numerical experiments.

4.5.1 Application to Model Selection
First, we will explain an application of the Main Theorems. In this paper, we theoreticallyclarify the diﬀerence between the variational free energy and the usual free energy in NMF.
From a practical point of view, the free energy F n can be calculated from the data; however,Chapter 4 Bayesian Generalization Error in Non-negative Matrix Factorization50

the entailed numerical integration is very hard and the sampling approximation, such an
exchange Monte Carlo method, spends a long time for finding F n . On the other hand, we cancompute the variational one F̄nvb more easily than F n . If the estimator of VBNMF is found,
all we have to do is to substitute it for the functional whose minimum value is equal to F̄nvb .
It has not been clarified how much the variational free energy diﬀers from the free energy;
however, the Main Theorems give the lower bound. We can use the lower bound to approximate the free energy from the variational one. Namely, when F̄nvb is known, we have theapproximation
F n ≈ F̄nvb − λ log n.
The usual VBNMF gives F n ≈ F̄nvb ; here though, we can obtain a more accurate value *1.
In this way we should be able to more accurately select the model in VBNMF by usingF̄nvb − λ log n.

4.5.2

Generalization Error

Second, we describe the generalization error in NMF. Theorem 4.3 also gives an upper boundof the generalization error Gn as well as the free energy Fn . Generally speaking, the learning
coeﬃcients that control the asymptotic behavior of the Fn and Gn are the same RLCTs [85];hence, we can clarify both behaviors at once. Since the situation in which the probability
model p(X|U, V ) is a Poisson distribution and the prior φ(U, V ) is a gamma distribution isa case where the Gibbs sampling [18] of NMF is performed, it can be regarded that not only
Fn but also Gn are theoretically clarified when Gibbs sampling is applied.
By contrast, in Theorem 4.1, only the learning coeﬃcient of the variational free energy F̄nvb
is determined. This is because the learning coeﬃcient of the variational generalization error
is not equal to the one of F̄nvb . Generally, in the case of VB, no zeta function is capable ofuniformly handling Fn and Gn , and the RLCT cannot obtain the learning coeﬃcient *2. For
example, in VB of three-layered linear neural networks, the asymptotic behaviors are clarified
not only with the variational free energy but also the variational generalization error [60], andtheir learning coeﬃcients are diﬀerent. A linear neural network is also known as a reduced
rank regression, a dimension reduction model, and the parameters are equivalent to a matrix
factorization model without a non-negative value constraint. In contrast in Bayesian inferencein matrix factorization and NMF, the RLCT of matrix factorization is a lower bound for the
RLCT of NMF, and it is known that the non-negative rank is dominant rather than the rank of
the matrix in NMF, as described in [37, 36]. Therefore, we cannot directly apply the resultsof linear neural networks to VBNMF.
In this way, theoretical generalization error in VB is rarely clarified, although that in
Bayesian inference has been clarified with the free energy. The Main Theorems show that
Gibbs sampling is more reliable than VB, in the sense that it gives a theoretical guarantee not
only about the free energy but also the generalization error. We can estimate the sample sizeto achieve the needed inference performance and tune the hyperparameters. Although various
factors determine whether Gibbs sampling or VB is appropriate, our research can answer thequestion of whether or not the theoretical generalization error is clarified.
*1 Actually λ has the true non-negative rank H0 ; however, in the same way as sBIC [24], we can avoid using
the true knowledge by considering H0 = 0, . . . , H.
2
* The learning coeﬃcient of VB is not equal to the RLCT.

4.5 Discussion

51

4.5.3 Effect of Non-negative Restriction
Third, the eﬀect of non-negative restriction to the parameter region is considered. In NMF,
the entries of the parameter matrices (U, V ) are non-negative. On the other hand, we canconsider non-restricted matrix factorization: the entries of the parameter matrices can be
negative. Usually, this is just called matrix factorization (MF) and its RLCT is defined asfollows.
Definition 4.2 (RLCT of MF) Let C be a compact set of R and let U MF ∈ M(M, H, C),
V MF ∈ M(H, N, C), U0MF ∈ M(M, r, C) and V0MF ∈ M(r, N, C). The largest pole of thefollowing univariate complex function
∫
∫
(
)z
MF
ζ(z)=
dU
dV MF ∥U MF V MF −U0MF V0MF ∥2 φMF (U MF , V MF )M(M,H,C)

M(H,N,C)is denoted by (−λMF ) and its order is denoted by mMF , where φMF is a prior distribution of(U MF , V MF ). Then λMF is called the RLCT of MF and mMF is called its multiplicity.
The RLCT of MF was clarified in [10] as that of reduced rank regression (see Theorem4.5). As proved above, it is a lower bound of the RLCT of NMF (see Theorem 4.6). We
discuss what the case λ = λMF holds or what the case λ > λMF holds. The control variableof MF is H = rankU V and it is compared with r = rankU0 V0 .
To discuss relationship between NMF and MF, there is an important concept; the smallest
inner matrix dimension of NMF is called a non-negative rank [19]. A useful property of a
non-negative rank is as follows.
Theorem 4.6 (Cohen) Suppose W ∈ M(M, N, K). Let rank+ W be the nonnegative rankof W . The following inequality holds:
rankW ≦ rank+ W ≦ min{M, N }.
A suﬃcient condition of rankW = rank+ W is rankW ≦ 2 ∨ M ≦ 3 ∨ N ≦ 3.
From the above theorem, if r ≦ H0 ≦ 2, then r = H0 . An example of non-negative
matrices whose non-negative rank is strictly larger than its rank is as follows:
1 1 0 0
 1 0 1 0 W = 0 1 0 1 .
0 0 1 1
The above 4 × 4 non-negative matrix W satisfies rankW = 3 and rank+ W = 4 [19].
In the case H = H0 = 1 or H = H0 = 2, H0 = r = 1 or H0 = r = 2 hold, respectively;thus, the RLCT of NMF is equal to that of MF in these case. On the other hand, the case that
λ > λMF holds is also included in Theorem 4.2. If H0 = 0, then λ can be strictly greaterthan λMF . We prove this property below.

52

Chapter 4 Bayesian Generalization Error in Non-negative Matrix Factorization
Proof. According to Theorem 4.6, r = 0 holds from H0 = 0. Hence, we consider Theorem
4.5 in the case r = 0 and Theorem 4.2 in the case H0 = 0. Then, we haveλ=H min{M ϕU , N ϕV }.
2

(4.122)

To keep the condition consistent with Theorem 4.5, let ϕU = ϕV = 1. Here, we calculate
λMF subject to Theorem 4.5 in the followings.
Case (1): N ≦ M + H ∧ M ≦ N + H ∧ H ≦ M + N ∧ M + H + N ≡ 0 (mod 2).
We assume N ≦ M , i.e. λ = HM/2. Owing to Theorem 4.5,{}
λMF = 2H(M + N ) − (M − N )2 − H 2 /8{}
= 2(M H + HN + N M ) − M 2 − N 2 − H 2 /8

(4.123)(4.124)holds; thus we have{}
λMF − λ = − M 2 + H 2 + N 2 − 2(M H + HN + N M ) − 4HM /8(4.125){ 2}(4.126)= − M + H 2 + N 2 − 2(−M H + HN + N M ) /8[]
= − (−M )2 + (−H)2 + N 2 + 2{(−M )(−H) + (−H)N + N (−M )} /8(4.127)= −(N − M − H)2 /8.

(4.128)Because of N ≦ M and H ≧ 1, N < M + H holds. Therefore, λMF < λ. If N > M , thatcan be derived in the same way as above.
Case (2): N ≦ M + H ∧ M ≦ N + H ∧ H ≦ M + N ∧ M + H + N ≡ 1 (mod 2).
We assume N ≦ M , i.e. λ = HM/2. In the same way as Case (1), we have
λMF − λ = 1/8 − (N − M − H)2 /8(4.129)= 12 /8 − (M + H − N )2 /8(4.130)= −(M + H − N + 1)(M + H − N − 1)/8.

(4.131)

We derive N + 1 ≦ M + H by reductio ad absurdum. We suppose M + H < N + 1.
Using the assumption of Case (2), N < M + H < N + 1 holds. M + H and N are
positive integers; thus the above inequality is inconsistent. That is why N + 1 ≦ M + H and−(M + H − N + 1)(M + H − N − 1) ≦ 0. Therefore, λMF ≦ λ. If N > M , that can bederived in the same way as above.
Case (3): N + H < M i.e. N < N + H < M .
On account of N < M ,λ = HM/2 = λMF .
Case (4): M + H < N i.e. M < M + H < N .
In the same way as Case (4), λ = HN/2 = λMF .
Case (5): M + N < H i.e. N < M + N < H ∧ M < M + N < H.
On account of M < H ∧N < H, M N < HN ∧M N < HM i.e. M N < H min{M, N }holds. Thus, we have λ = H min{M, N }/2 > M N/2 = λMF .
□4.5 Discussion

53

4.5.4 Robustness on Probability Distributions
Third, let us discuss the true distribution and the model of the data. In this study, we
consider the case in which the probability model p(X|U, V ) is a Poisson distribution andthe prior φ(U, V ) is a gamma distribution in the same way as in the derivation of the Gibbs
sampling algorithm of NMF by Cemgil [18]. These assumptions are necessary for Gibbs
sampling and the derivation of VB; however, other models can be considered when using
other MCMC methods. For example, we may want to set that the entries of the data matrix
X are non-negative real numbers. Is the main result applicable to these cases?
To tell the truth, several distributions satisfy the condition that the RLCT of NMF is equal
to the absolute value of the maximum pole of the following zeta function∫∫()zζ(z) =dU
dV ∥U V − U0 V0 ∥2 φ(U, V ).
M(M,H,K)

M(H,N,K)

Specifically, when the elements of the data matrix follow a normal distribution, a Poisson
distribution, an exponential distribution, or a Bernoulli distribution, the asymptotic behaviorof the free energy and the generalization error can be described using the same RLCT definedby the above zeta function: Theorem 4.2. Therefore, if the prior distribution is a gamma
distribution, Theorem 4.2 and Theorem 4.3 hold not only when the probability model andthe true distribution are Poisson distributions but also when they are normal distributions,
exponential distributions, or Bernoulli distributions.
We also consider the case in which X is generated by an exponential distribution. Then,the KL divergence has the same RLCT as the square error if elements of U V are positive.
Proposition 4.1 Let the probability density functions of X ∈ M(M, N, K) be q(X) andp(X|U, V ), which represent the true distribution and the model respectively defined byq(X) ∝ Exp(X|U0 V0 ),p(X|U, V ) ∝ Exp(X|U V ),where Exp(X|W ) is a probability density function of an exponential distribution with average
W:
M ∏
N
∏
e−xij /wij
Exp(X|W ) :=, X = (xij ), W = (wij ).
wij
i=1 j=1
Also let φ(X, Y ) be a probability density function such that it is positive on a compact subsetof M(M, H, K0 ) × M(H, N, K0 ). Then, the Kullback-Leibler divergence has same RLCT as
the square error.
Proof. Let x > 0, a > 0 and b > 0. We putp(x|a) :=g(a, b) :=

e−x/a
,
∫ a
p(x|a)log(4.132)p(x|a)dx.
p(x|b)(4.133)

Chapter 4 Bayesian Generalization Error in Non-negative Matrix Factorization

54

By the similar way to the proof of Lemma 4.1, all we have to do is prove g(a, b) ∼ (a − b)2 .
By using
log

p(x|a)b= log e−x/a+x/bp(x|b)
ax x= log b − log a − + ,ab(4.134)(4.135)

we have∫p(x|a)dx = 1,

(4.136)xp(x|a)dx = E[x] = a.

(4.137)∫Developing the terms, we obtain∫(x x)g(a, b) = p(x|a) log b − log a − +dxaba= log b − log a − 1 + .
b(4.138)Hence, we immediately get
1 1
− ,
b a
1
a
∂b g(a, b) = − 2
b b
b−a
= 2
b

∂a g(a, b) =and increase or decrease of g(a, b). Thus, this proposition can be proved in the same way as
Lemma 4.1.
□
Remark 4.4 The right sides of Eqs. (4.17) and (4.138) are respectively equal to ”Idivergence”[74] and ”Itakura-Saito-divergence”[45] which are used as criteria of diﬀerence
between observed matrix and reproduced matrix in NMF[51, 26, 27]. Moreover, in the same
way as the proof of Theorem 4.5 [10], we can prove that the square error ∥U V − U0 V0 ∥2 has
the same RLCT as that of the KL divergence when the entries of the data matrix are subject
to a normal distribution.
As a result, we can apply Theorem 4.2 to Bayesian inference if we use a Kullback-Leibler
divergence or a square error as a criterion of diﬀerence between U V and U0 V0 in cases where
elements of matrices are generated by normal, Poisson or exponential distributions. Thus,
the upper bound of the expected generalization error E[Gn ] and the free energy Fn can beclarified if the size and inner dimension of the observed matrix and reproduced matrix are
given.

4.5 Discussion

55

4.5.5 Experiment
Here, we run numerical experiments to check the numeric behavior of the theoretical results.
Theorem 4.1 gives the exact value of the learning coeﬃcient λvb of VBNMF and its validity
was confirmed in Kohjima’s previous research [49]. The core result of this chapter is Theorem 4.2. Therefore, we only have to run experiments for it; i.e., the RLCT λ of Bayesian NMF
is calculated by using Gibbs sampling.
∫
Let Eθ be an expectation operator of the posterior: Eθ [·] = dθψ(θ|D)[·]. Let λ̂ be thenumerically calculated RLCT. The widely applicable information criterion (WAIC) [82] isdefined by the following random variable Wn :Wn = Tn + Vn /n,where Tn is the empirical loss and Vn is the functional variance:
n
n
1∑
1∑
∗
log p (Xi ) = −log Eθ [p(Xi |θ)],Tn = −
n i=1
n i=1

Vn =

n [
∑
i=1

2

Eθ [(log p(Xi |θ)) ] − {Eθ [log p(Xi |θ)]}2(4.139)]=n∑Vθ [log p(Xi |θ)].

(4.140)i=1Even if the posterior distribution cannot be approximated by any normal distribution (i.e.,the model is singular), the expected WAIC E[Wn ] is asymptotically equal to the expected
generalization loss E[Gn + S] [82];E[Gn + S] = E[Wn ] + o(1/n2 ).
Moreover, the generalization error and the WAIC error Wn − Sn have the same variance [85]:Gn + Wn − Sn = 2λ/n + op (1/n).

(4.141)

Eq. (4.141) is useful for computing λ̂ because the leading term 2λ/n is deterministic. Nevertheless, the left hand side is probabilistic. This means that the needed number of simulation∑nTt)
D is smaller than that in the case using this approximation: Gn ≈ n1 t=1
log pq(X∗ (X ) , wheretnT is the number of the test data and (X ∗ )nT = (X1∗ , . . . , Xn∗T ) is the test data generated byq(X).
The method was as follows. First, the training data X n was generated from the true
distribution q(X). Second, the posterior distribution was calculated by using Gibbs sampling[18] (see also Algorithm 3). Third, Gn and Wn − Sn were computed by using the training
data X n and the artificial test data (X ∗ )nT generated from q(X). These three steps wererepeated and each value of n(Gn + Wn − Sn )/2 was saved. After all repetitions have been
completed, n(Gn + Wn − Sn )/2 was averaged over the simulations. This average was λ̂.
The pseudo-code is listed in Algorithm 2, where nK is the sample size of the parameter
subject to the posterior. We used the programing language named Julia 1.3.0 [13] for thisexperiment. The implementation is available at the following github page: https://github.
com/chijan-nh/LearningCoefficient-RLCT-ofNMF-usingGS.

56

Chapter 4 Bayesian Generalization Error in Non-negative Matrix Factorization
Algorithm 2 How to Compute λ̂Require: ϕ = (ϕU , θU , ϕV , θV ) > 0: the hyperparameters,
U0 : the true parameter matrix whose size is (M, H0 ),V0 : the true parameter matrix whose size is (H0 , N ),GS: the Gibbs sampling function whose return value consists of the samples from theposterior. See also Algorithm 3.
Ensure: The numerical computed RLCT λ̂.
Allocate an array Λ[D].
for d = 1 to D do
Generate X n ∼ q(X) = p(X|w0 ), where w0 = (U0 , V0 ).
Allocate arrays U[M, H, nK ] and V[H, N, nK ].
Get U, V ← GS(X n , ϕ).
Generate (X ∗ )nT ∼ q(X).
∑nT∑nq(X ∗ )
Calculate Gn ≈ n1T t=1
log Eθ [p(Xt∗ |θ)] , Sn = − n1 i=1 log q(Xi ),t ∑∑nn
and Wn ≈ − n1 i=1 Eθ [p(Xi |θ)] + n1 i=1 Vθ [log p(Xi |θ)],∑nK
where Eθ [f (θ)] ≈ n1K k=1
f (θk ) and θk = (U[:, :, k], V[:, :, k]).
Save Λ[d] ← n(Gn + Wn − Sn )/2.
end for∑D
1
Calculate λ̂ = D
d=1 Λ[d].

We set M = N = 4, H = 2, H0 = 1, and nT = 100n. To examine the behavior given
diﬀerent sample sizes, we set n = 500, 1000 and nK = 1000, 2000, respectively. To decrease
the probabilistic eﬀect of Eq. (4.141), we conducted the simulations twenty times: D = 20.
The hyperparameters were set to θU = θV = 1 and(ϕU , ϕV ) = (0.25, 0.25), (0.5, 0.5), (1, 1), (2, 2).
We chose four pairs of (ϕU , ϕV ) in view of a theoretical point: the critical lines of the
RLCT of NMF and the learning coeﬃcient in VBNMF. The RLCT or the learning coeﬃcient
drastically changes on these lines. Figure 4.1 shows this phenomenon called phase transition.
Under the condition M = N , ϕU + ϕV = 1 is the phase transition line (see Theorem 4.1).
Each point on the straight line ϕU = ϕV is characterized as follows. (0.25, 0.25) is before
the phase transition line, and (0.5, 0.5) is a phase transition point. (1, 1) is a case where theprior distribution φ(U, V ) is strictly positive and bounded. (2, 2) is one of points beyond thecritical line.
In the Gibbs sampling, we had to conduct a burn-in to decrease the eﬀect of the initialvalues and thin the samples in order to break the correlations. The sample size for the burn-in
was 20000, while the sample size for the thinning was 20; thus, the sample sizes of theparameter was 20000 + 20K = 40000 and 60000 (K = 1000 and 2000) and we used the(20000 + 20k)-th sample as the entry of U[:, :, k] and V[:, :, k] for k = 1 to K.
The experimental results are shown in Table 4.1. The symbol λ denotes the theoretical
upper bound of the RLCT λ in Theorem 4.2. There are columns for each sample size, andeach row contains the hyperparameter (Hyperparam.) and the learning coeﬃcient (Coeﬀ.).
The experimental values have three significant digits.

4.5 Discussion
Algorithm 3 Gibbs Sampling for NMF
GS(X n , ϕ, K = K, burnin = 20000, thin = 20)
Require: X n = (Xl )nl=1 : the data where Xl is an M × N non-negative integer matrix,
ϕ = (ϕU , θU , ϕV , θV ) ∈ R4>0 : the hyperparameter of the Gamma priors for the nonnegative matrices U and V ,Gam(W |Φ, Θ): a Gamma distribution of a matrix whose (i, j)-entry is generated by
Gam(w|Φ[i, j], Θ[i, j]),Multi(s|x, π): a Multinomial distribution whose trials and event probabilities are x ∈ N
and π ∈ Sim(H, [0, 1]), respectively.
Ensure: Sampling non-negative matrices from the numerical posterior.
Let iter = butnin + thin ∗ K.
Allocate arrays U[M, H, K], U[M, H, iter], V[H, N, K] and V[H, N, iter].
# Inital sampling for U and V from the prior:
Generate U, V ∼ Gam(U |ϕU , θU )Gam(V |ϕV , θV ).
# Sampling from the posterior:for k = 1 to iter do## Sampling the hidden variable s.
Allocate an array s[M, H, N, n].
for i = 1, j = 1, and l = 1 to M , N , and n dofor h = 1 to H do
∑H
Let π[i, h, j] = U [i, h]V [h, j]/ h=1 U [i, h]V [h, j].
end for
Generate s[i, :, j, l] ∼ Multi(s|Xl [i, j], π[i, :, j]).
end for## Sampling the non-negative matrix U .
for i = 1 and h =∑1 to Mand H don ∑NLet ϕ̂U [i, h] = l=1 j=1 s[i, h, j, l] + ϕU .
∑NLet θ̂U [i, h] = n j=1 V [h, j] + θU .
end for
Generate U ∼ Gam(U |ϕ̂U , θ̂U ).
Put U[M, H, k] ← U .
## Sampling the non-negative matrix V .
for h = 1 and j = 1 to H and N do∑n ∑MLet ϕ̂V [h, j] = l=1 i=1 s[i, h, j, l] + ϕV .
∑MLet θ̂V [h, j] = n i=1 U [i, h] + θV .
end for
Generate V ∼ Gam(V |ϕ̂V , θ̂V ).
Put V[H, N, k] ← V .
end for# Burn-in and thinning.
for k = 1 to K do
U[M, H, k] ← U[M, H, burnin + thin ∗ k].
V[H, N, k] ← V[H, N, burnin + thin ∗ k].
end for
Return U, V.

57

Chapter 4 Bayesian Generalization Error in Non-negative Matrix Factorization

58

Table 4.1: Numerically Calculated and Theoretical Values of the Learning Coeﬃcients
Hyperparam. / Coeﬀ.
λvb
ϕU =ϕV =0.25
λθU =θV =1λ̂
λ − λ̂
λvb
ϕU =ϕV =0.5
λ
θU =θV =1λ̂
λ − λ̂
λvb
ϕU =ϕV =1λθU =θV =1λ̂
λ − λ̂
λvb
ϕU =ϕV =2λθU =θV =1λ̂
λ − λ̂

n = 500
6
4
3.74 ± 0.0412
0.260 ± 0.0412
8
9/2
4.05 ± 0.0706
0.450 ± 0.0706
8
11/2
4.52 ± 0.0492
0.976 ± 0.0492
8
15/2
4.76 ± 0.0454
2.74 ± 0.0454

n = 1000643.73 ± 0.0508
0.268 ± 0.0508
8
9/2
4.15 ± 0.0842
0.346 ± 0.0842
8
11/2
4.54 ± 0.0628
0.965 ± 0.0628
8
15/2
4.79 ± 0.0477
2.71 ± 0.0477

As shown in Table 4.1, all numerically calculated values are smaller than the theoretical
upper bound, thus, the Theorem 4.2 is consistent with the experimental result. Since λ̂ islarger in the case ϕU = ϕV = 2 than that in the case ϕU = ϕV = 1, it is conjectured that thelarger ϕU and ϕV are, the larger λ will be, while λvb saturates to H(M + N )/2 = 8 (owing
to Theorem 4.1) and the upper bound of λ looks less tight than that in the other case. In thecase ϕU = ϕV = 0.5, when the hyperparameter is on the phase transition line, λ̂ fluctuates
more than in the other case. We can see that learning is unstable on the critical line. Hence,
the hyperparameters should be set to avoid the neighborhood of the phase transition line.

4.6

Conclusion

We give an upper bound of the RLCT for NMF whose priors are gamma distributions (Theorem
4.2) and describe theoretical applications to Bayesian and variational inference. According
to Theorem 4.3, a theoretical upper bounds of the generalization error and the free energy
are derived and hyperparameters make them change: there are a phase transition structure.
Moreover, owing to Theorem 4.4, the variational approximation error, i.e., the diﬀerence
between the variational free energy and the free energy, are quantitatively evaluated. These
diﬀerence depends on the true non-negative rank and the hyperparameters from the gamma
prior distributions. The numerical results are consistent with the theoretical results and theysuggest the exact values of the RLCT and the stability of learning. Future tasks include
conducting large-scale experiments and clarifying the exact value of the RLCT.

4.6 Conclusion

Fig. 4.1: This image represents the critical lines of the RLCT of NMF and the learning coeﬃcient in VBNMF when M = N . The horizontal and the vertical axises arecorresponding to the hyperparameter ϕU and ϕV , respectively. The used four points{(0.25, 0.25), (0.5, 0.5), (1, 1), (2, 2)} are plotted as black disks. The red-dashed-line isthe critical line of the RLCT of NMF by Theorem 4.2. On the other hand, the blue-line isthat of the learning coeﬃcient in VBNMF by Theorem 4.1. Obviously, they are thoroughly
diﬀerent: VB is not equivalent to Bayesian inference.

59

61

Chapter 5

Bayesian Generalization Error in
Latent Dirichlet Allocation
In this chapter, we report the result of theoretical analysis for Bayesian generalization error in LDA. This chapter consists of five parts. First, in Sec. 5.1, we
describe motivation of this theoretical research. Second, in Sec. 5.2, we state themain theorem with regard to Bayesian generalization error in LDA. Third, in Sec.
5.3, we prepare for the proof of the theorem. Fourth, in Sec 5.4, we prove themain theorem. Lastly, in Sec. 5.5, we discuss the theoretical results.
In the followings, θ = (A, B) is a parameter and x is an observed randomvariable. Note that this chapter is based on the author’s papers [38, 34].

5.1 Motivation
Topic model [29] is a ubiquitous learning machine used in many research areas, including
text mining [15, 30], computer vision [52], marketing research [72], and geology [98]. Latent Dirichlet allocation (LDA) [15] is one of the most popular Bayesian topic models. It
has been devised for text analysis, and it can extract essential information in documents bydefining the topics of the words. The topics are formulated as one-hot vectors subject tocategorical distributions which are diﬀerent for each document (Fig. 5.1). This formulationrefers documents to bags of words: a text is simply referred to a multiset of words [31]. For
example, a text “I think that that that that that boy wrote is wrong” is referred to a multiset{I, think, that, that, that, that, that, boy wrote, is, wrong}m . To consider the frequency of appearance of each word, a muliset allows duplication. By using one-hot-encoding, this multiset
is transformed to the vectors                     
0
0
0
0
0
0
0
0
0
0
1
0 1 0 0 0 0 0 0 0 0 0                     
0 0 1 1 1 1 1 0 0 0 0                     
0 , 0 , 0 , 0 , 0 , 0 , 0 , 1 , 0 , 0 , 0 .
(5.1)
                     
0 0 0 0 0 0 0 0 1 0 0                     
0 0 0 0 0 0 0 0 0 1 0
1
0
0
0
0
0
0
0
0
0
0

62

Chapter 5 Bayesian Generalization Error in Latent Dirichlet Allocation
If the text is the given document, LDA models that the above vectors are generated from the
mixture of categorical mixture distribution in Fig. 5.1 in the case N = 1. Practically, a single
document has more words and the number of documents is also larger *1. In fact, topic model(including LDA) can treat bags of anythings by formulating the raw data to the bags; thus,
it has many application described above. The standard inference algorithms, such as Gibbs
sampling [30] and the variational Bayesian method [15], require the appropriate number of the
topics to be set in advance. If the chosen number of topics is too small, then LDA suﬀers fromunderfitting. On the other hand, if the chosen number of topics is too large, it suﬀers fromoverfitting on the training data. In practical applications, neither the optimal number of topicsof the ground truth nor the true distribution is known; thus, researchers and practitioners facea situation in which the number of topics they set may be larger than the optimal one. Thus,
clarifying the behavior of the generalization error is necessary as a theoretical foundation toresolve model selection problems. However, the mathematical property of LDA has not yet
been clarified, because it has a hierarchical structure; the regular statistical theory can not beapplied. To rephrase, LDA is a singular model.

Fig. 5.1: This figure gives an overview of LDA. The categorical distributions Cat that dependon the documents. Words in the uppercase such as NAME, FOOD, and CODING are topics.
There are categorical distributions that are diﬀerent for each topic; the words (Ciel, curry,
lambda, ...) are generated from them. Hence, we can explain LDA as a mixture of categoricalmixture models. This model has a hierarchical structure. This figure is quoted and modified
from the author’s works [38, 34]
Matrix factorization (MF) is also used in machine learning frequently. MF decomposes
the data matrix into a product of two matrices and discovers hidden structures or patterns,*1 Besides, for practical uses, one may remove articles (e.g. “a” and “the”) and demonstratives (e.g. “this” and“that”) to make a language model [68]. Also, one may refer a text to a bag of n-grams and vectorize it.

5.1 Motivation

63

hence it has been experimentally used for knowledge discovery in many fields. However,MF has no guarantee of reaching the unique factorization, and it is sensitive to the initialvalue of the numerical calculation. This non-uniqueness interferes with data-driven inference
and interpretations of the results. Besides, the sensitivity to the initial value causes thefactorization result to have low reliability. From the viewpoint of data-based prediction, thisinstability may lead to incorrect predictions. To improve interpretability, non-negative matrix
factorization (NMF) [61, 51] has been devised; it is a restricted MF wherein the elements ofthe matrix are non-negative. Thanks to the non-negativity constraint, the extracted factors
are readily interpretable, therefore NMF is frequently used for extracting latent structures andpatterns in many fields (see also Sec. 4.1).
Stochastic matrix factorization (SMF) was devised by Adams [2]; it can be understood asa restriction of NMF in which at least one matrix factor is “stochastic”: the elements of thematrix factors are non-negative and the sum of the elements in a column is equal to one. In
other words, the columns of the matrix are in a simplex. By making two further assumptions,
Adams proved the uniqueness of the results of SMF [1, 2]. For a statement of these two
conditions, let us consider a data matrix X whose size is M × N and factor matrices A and B
which are “stochastic” and whose sizes are M × H and H × N , respectively. H might be therank of X but the “stochastic” condition makes this determination non-trivial. In other words,
SMF can be viewed as a method that finds a factor matrices pair (A, B) such that X = AB
for a given X and H. The non-uniqueness property has been paraphrased as the existence of
H × H regular matrix P ̸= IH such that
X = AP P −1 B,(5.2)

where IH is an H × H identity matrix. Thus, uniqueness means that Eq. (5.2) is attained if
and only if P = IH . Adams assumed thatAP ≧ 0 and P −1 B ≧ 0,

(5.3)

i.e., the elements of AP and P −1 B are non-negative, and P −1 B =: (b′kj ) satisfies
H∑
k=1

b′kj = 1 orN∑

b′kj = 1.

(5.4)

j=1

Adams claimed that these assumptions are “natural” and applied SMF to image recognition(the same problem analyzed in [51]) and text mining[2]. It is emphasized that, in this thesis,
we consider that it is not clear whether Adam’s assumptions (5.3) and (5.4) are mathematically“natural”. In the following, we do not assume Eqs. (5.3) and (5.4).
The MF methods described so far, including SMF, are understood as a deterministic
procedure. As will be shown later, for hierarchical learning machines such as MF, we study
probabilistic procedures, because Bayesian inference has higher predictive accuracy thandeterministic methods or maximum likelihood estimation. The same is also true regarding theaccuracy of the discovered knowledge. Moreover, the probabilistic or statistical view gives
a wider application. Indeed, Bayesian NMF [75, 18] has been applied to image recognition[18], audio signal processing [75], and recommender systems [16]. From a statistical point ofview, the data matrices are random variables subject to the true distribution. As mentioned inChapter 5 Bayesian Generalization Error in Latent Dirichlet Allocation

64

Chap. 4, the factorization of a set of independent matrices should be studied because the target
matrices are often obtained daily, monthly, or in diﬀerent places. More importantly, as proved
later, the SMF has the same learning curve as LDA; if the Bayesian generalization error in
SMF has been clarified, then that of LDA is also determined. That is why the decompositionof a set of matrices with stochastic restriction is considered to be a statistical inference inorder to clarify the Bayesian generalization error in LDA.
In this chapter, we derive the exact asymptotic form of the Bayesian generalization error
by determination of the exact real log canonical threshold λ (RLCT; see also Chap. 3 for thegeneral definition and property) in LDA. For finding λ, we prove that the RLCT of LDA isequal to that of SMF, i.e. LDA is equivalent to SMF in the sense of generalization.

5.2

Main Theorem

Now let us introduce the main result of this chapter.
Let K be a compact subset of [0, 1] = {x ∈ R|0 ≦ x ≦ 1} and let K0 be a compact
subset of (0, 1) = {x ∈ R|0 < x < 1}. Let Onehot(N ) := {w = (wj ) ∈ {0, 1}N |∑N
j=1 wj = 1} = {(1, 0, . . . , 0), . . . , (0, . . . , 0, 1)} be an N -dimensional one-hot vector set∑N
and Sim(N, K) := {c = (cj ) ∈ K N | j=1 cj = 1} be an N -dimensional simplex. Let
S(M, N, E) = Sim(M, E)N be a set of M × N stochastic matrices whose elements are in
E, where E is a subset of [0, 1], and M, N ∈ N. In addition, we set H, H0 ∈ N and H ≧ H0 .
In LDA terminology, the number of documents and the vocabulary size are denoted by N
and M , respectively. Let H0 be the true or optimal number of topics and H be the chosenone. In this situation, the sample size n is the number of words in all of the given documents.
See also Table 5.1. This table is quoted and modified from the author’s works [38, 34].
Table 5.1: Description of Variables in LDA Terminology
Variablebj =(bkj ) ∈ Sim(H, K)
ak =(aik ) ∈ Sim(M, K)
x=(xi ) ∈ Onehot(M )y=(yk ) ∈ Onehot(H)z=(zj ) ∈ Onehot(N )∗0 and ∗0

Description
probability that topic is k when document is j
probability that word is i when topic is kword i is defined by xi = 1
topic k is defined by yk = 1
document j is defined by zj = 1
true or optimal variable corresponding to ∗Index
k=1, ..., H
i=1, ..., M
i=1, ..., M
k=1, ..., H
j=1, ..., N-

Assume that M ≧ 2, N ≧ 2, and H ≧ H0 ≧ 1. We define A = (aik ) ∈ S(M, H, K) andB = (bkj ) ∈ S(H, N, K), and assume that A0 = (a0ik ) ∈ S(M, H0 , K0 ) and B0 = (b0kj ) ∈S(H0 , N, K0 ) are SMFs such that they give the minimal factorization of A0 B0 .
Let q(x|z) and p(x|z, A, B) be conditional probability density functions of x ∈ Onehot(N )given z ∈ Onehot(M ), which represent the true distribution and the model of LDA, respec-

5.2 Main Theorem

65

tively,q(x|z) =p(x|z, A, B) =

N
∏

(H
0
∑

j=1

k=1

(H
N
∏
∑

j=1

b0kj

M
∏

)zj
0 xi
(aik
)

i=1

bkj

M
∏

,

(5.5)

.

(5.6)

)zj
(aik )xi

i=1

k=1

These distributions are the marginalized ones of the following simultaneous ones with respect
to the topics y 0 ∈ Onehot(H0 ) and y ∈ Onehot(H):
N
∏



(

M
∏

)yk0 zj
 ,

(5.7)

[H (
)yk ]zjNM∏∏∏.
p(x, y|z, A, B) =bkj
(aik )xi(5.8)

q(x, y 0 |z) =

j=1

j=1



H0
∏

b0kj

k=1

k=1

(a0ik )xi

i=1

i=1

In practical cases, the topics are not observed; thus, we use Eqs. (5.5) and (5.6). Besides, letφ(A, B) > 0 be a probability density function such that it is positive on a compact subset ofS(M, H, K) × S(H, N, K) including Φ−1 (0) i.e. (A0 , B0 ).
Definition 5.1 (RLCT of LDA) Let KL(A, B) be the KL divergence from q(x|z) top(x|z, A, B) in the aboves:
KL(A, B) :=∑∑q(x|z)q ′ (z) logz∈Onehot(M ) x∈Onehot(N )q(x|z),
p(x|z, A, B)

where q ′ (z) is the true distribution of the document. In LDA, q ′ (z) is not observed andassumed that it is positive and bounded.
Then, the holomorphic function of one complex variable z (Re(z) > 0)∫∫
ζ(z) =dA
dB KL(A, B)z
S(M,H,K)

S(H,N,K)

can be analytically continued to a unique meromorphic function on the entire complex plane
C and all of its poles are rational and negative (see also Theorem 3.2). If the largest pole is(−λ), then λ is said to be the RLCT of LDA.
Definition 5.2 (RLCT of SMF) Set Φ(A, B) = ∥AB − A0 B0 ∥2 . Then the holomorphic
function of one complex variable z (Re(z) > 0)∫∫
ζ(z) =dA
dB Φ(A, B)z
S(M,H,K)

S(H,N,K)

can also be analytically continued to a unique meromorphic function on C and its all poles
are rational and negative. If the largest pole is z = −λ, then λ is the RLCT of SMF.

Chapter 5 Bayesian Generalization Error in Latent Dirichlet Allocation

66

In this paper, we prove the following two theorems.
Theorem 5.1 (Equivalence of the LDA and SMF) Let λSMF be the RLCT of SMF and λLDA
be the RLCT of LDA. Then, the following equality holds:
λSMF = λLDA .
In the same way,
mSMF = mLDA
holds where mSMF and mLDA are the multiplicities of SMF and LDA, respectively.
In order to prepare to state Main Theorem of this chapter, we define an intrinsic value r of\(M,H0 )\(H ,1)and B̃0 0
the true distribution defined by stochastic matrices A0 and B0 . Let Ã0
be matrices defined as\(M,H0 )−1,H0 −1
= (a0ik )M
(5.9)
Ã0
i=1,k=1
and\(H0 ,1)B̃0(\M ),H0

respectively. Also, Ã0
are same such that(\M ),H0Ã0
and(\H0 ),1and B̃0(5.10)are denoted by matrices whose column vectors−1,H0 −1
= (a0(1:M −1)H0 , . . . a0(1:M −1)H0 ) = (a0iH0 )M
i=1,k=1(5.11)

H0 −1,N −1= (b0(1:H0 −1)N , . . . b0(1:H0 −1)1 ) = (b0k1 )k=1,l=1,(5.12)(\H0 ),1B̃00 −1,N= (b0kj )H
k=1,j=2 ,−10
0 H0 −1
where a0(1:M −1)H0 = (a0iH0 )M
i=1 and b(1:H0 −1)1 = (bk1 )k=1 are an (M − 1)-dimensional
vector and an (H0 − 1)-dimensional vector, respectively. Then, let\(M,H0 )U0 = Ã0\(H0 ,1)V0 = B̃0(\M ),H0− Ã0(\H0 ),1

− B̃0

,

,

(5.13)(5.14)

and r = rank(U0 V0 ). Obviously, r depends on H0 .
The main result of this paper is the following theorem.
Theorem 5.2 (Main Theorem) Suppose M ≧ 2, N ≧ 2, and H ≧ H0 ≧ 1. Let r be therank of U0 V0 which is a product of two matrices (U0 , V0 ) defined in Eqs. (5.13) and (5.14).
The RLCT of LDA λ = λLDA and its multiplicity m = mLDA are as follows.
1. If N + r + 1 ≦ M + H and M + r + 1 ≦ N + H and H + r + 1 ≦ M + N ,（a）and if M + N + H + r is odd, then
λ=1
1{2(H + r + 1)(M + N ) − (M − N )2 − (H + r + 1)2 } − N, m = 1.
82（b）and if M + N + H + r is even, then

λ=1
1{2(H + r + 1)(M + N ) − (M − N )2 − (H + r + 1)2 + 1} − N, m = 2.
82

5.3 Preparation

67

2. Else if M + H < N + r + 1, thenλ=1
{M H + N (r + 1) − H(r + 1) − N }, m = 1.
2

3. Else if N + H < M + r + 1, thenλ=1
{N H + M (r + 1) − H(r + 1) − N }, m = 1.
2

4. Else (i.e. M + N < H + r + 1), thenλ=1
(M N − N ), m = 1.
2

We prove Main Theorem in the next section. As two applications of this theorem, we
obtain the exact form of the free energy and Bayesian generalization error of LDA by applying
Theorem 3.4.
Theorem 5.3 Let Fn be the normalized free energy and E[Gn ] be the expected generalization
error in LDA, respectively. Then, these asymptotic forms are as the followings:
Fn = λ log n − (m − 1) log log n + Op (1),
(
)
1
λ
m−1
E[Gn ] = −+o,n n log nn log nwhere λ and m are the RLCT of LDA and its multiplicity which are determined in Theorem
5.2.
Theorem 5.3 can be immediately derived if Theorem 5.2 is proved. Thus, we prove Theorem
5.2 in the following Secs. 5.3 and 5.4. To prove Theorem 5.2, we derive Theorem 5.1 in thenext section.

5.3 Preparation
In order to prove Theorem 5.1, we show the following facts.
0 M
Proposition 5.1 Let w = (wi )M
i=1 ∈ Sim(M, K) and w0 = (wi )i=1 ∈ Sim(M, K0 ). Then,
there exist positive constants c1 > 0 and c2 > 0 such that

c1

M
∑
i=1

I.e.

(wi − wi0 )2 ≦

M
∑

M
∑
wi0
wi log≦ c2(wi − wi0 )2 .
w
i
i=1
i=1

∑M

∑M
wi0
0 2
i=1 wi log wi has the same RLCT as
i=1 (wi − wi ) .

This proposition was derived in [53] in order to clarify the RLCT of categorical mixture
model. By using Proposition 5.1, we prove Theorem 5.1.

Chapter 5 Bayesian Generalization Error in Latent Dirichlet Allocation

68

Proof of Theorem 5.1. Without loss of generality, we can rewrite the notation of q(x|d) andp(x|d, A, B) as follows:q(x|zi = 1) =

H0
∑

b0kj

M
∏

(a0ik )xi ,p(x|zi = 1, A, B) =

i=1

k=1

H
∑

bkj

M
∏

(aik )xi .

(5.15)

aik bkj .

(5.16)

i=1

k=1

The word x is a one-hot vector; hence, we obtainq(xj = 1|zi = 1) =H0∑

a0ik b0kj ,p(xj = 1|zi = 1, A, B) =

k=1

H
∑
k=1

Then, the conditional Kullback-Leibler divergence between q(x|d) and p(x|d, A, B) isequal toKL(A, B) =∑∑q(x|z)q ′ (z) logz∈Onehot(M ) x∈Onehot(N )M ∑N∑q(x|z)
p(x|z, A, B)

q(xj = 1|zi = 1)p(xj = 1|zi = 1, A, B)
i=1 j=1)(H∑H0 0 0NM0∑∑∑
k=1 aik bkj
a0ik b0kj log ∑H=q ′ (zi = 1).
k=1 aik bkj
j=1
i=1 k=1=q(xj = 1|zi = 1)q ′ (zi = 1) log(5.17)(5.18)(5.19)

Owing to A = (aik ) ∈ S(M, H, K), B = (bkj ) ∈ S(H, N, K), A0 = (a0ik ) ∈S(M, H0 , K0 ), and B0 = (b0kj ) ∈ S(H0 , N, K0 ), the (i, j) entries of AB and A0 B0 are∑H∑H0 0 0(AB)ij := k=1 aik bkj and (A0 B0 )ij := k=1
aik bkj . We have
KL(A, B) =N∑′q (zi = 1)M∑

(A0 B0 )ij log

i=1

j=1

(A0 B0 )ij.
(AB)ij(5.20)∑M(A0 B0 )ij
in Eq. (5.20) has the same
According to Proposition 5.1, i=1 (A0 B0 )ij log (AB)
ij
∑M
2
′
RLCT of i=1 ((AB)ij − (A0 B0 )ij ) . In addition, q (zi = 1) is positive and bounded.
Accordingly, we have
KL(A, B) =N∑

(5.21)((AB)ij − (A0 B0 )ij )2(5.22)((AB)ij − (A0 B0 )ij )2 = ∥AB − A0 B0 ∥2 .

(5.23)q (zi = 1)

j=1

∼

N
∑

i=1
′

q (zi = 1)

j=1

∼

M
∑

(A0 B0 )ij(AB)ij

′

N ∑
M
∑

M
∑

(A0 B0 )ij log

i=1

j=1 i=1

Therefore, KL(A, B) ∼ ∥AB − A0 B0 ∥2 ; i.e., the RLCT of LDA equals the RLCT of SMF.

5.3 Preparation

69□
The RLCT of MF was clarified in the all case [10] as the follwoing theorem. Note that thistheorem is same as Theorem 4.5 in Chap. 4; however, for self-containedness of this chapter,
we repeat it.
Theorem 5.4 (Aoyagi) Let M , N , H be positive integers. Matrices U and V are denotedby an M × H and an H × N matrix whose entries are real numbers. Let U0 and V0 be
an M × ∗ matrix and an ∗ × N matrix, respectively. Suppose U0 and V0 are constants andput r = rank(U0 V0 ) and Φ(U, V ) = ∥U V − U0 V0 ∥2 . Let λMF and mMF be the RLCT of
MF and its multiplicity, respectively. I.e. (−λMF ) is the maximum pole of the following zetafunction and mMF is its order:∫∫ζMF (z) =dU dV Φ(U, V )z .
Then, λMF and mMF are as follows:(1) If N + r ≦ M + H ∧ M + r ≦ N + H ∧ H + r ≦ M + N ,(1−1) in the case M + H + N + r is even,λMF =−(H +r)2 −M 2 −N 2 + 2{(H +r)(M +N )+M N }.
8(1−2) in the case M + H + N + r is odd,
1−(H +r)2 −M 2 −N 2 + 2{(H +r)(M +N )+M N }
λMF =
.
8
(2) Else if M + H < N + r,
λMF =HM − Hr + N r
.
2

λMF =

HN − Hr + M r.
2(3) Else if N + H < M + r,

(4) Or else, i.e. in the case M + N < H + r,
λMF =

MN
.
2

In the case (1−2), the multiplicity is two: mMF = 2. Otherwise, it equals one: mMF = 1.
The RLCT of MF and its multiplicity depend on (M, N, H, r); thus, we write them
λMF = λMF (M, N, H, r) and mMF = mMF (M, N, H, r), respectively.
In the next section, by using changes of variables, we come down Theorem 5.2 (Main
Theorem) to Theorem 5.4. This is why we consider the rank of U0 V0 , where U0 and V0 are
defined in Eqs. (5.13) and (5.14).

Chapter 5 Bayesian Generalization Error in Latent Dirichlet Allocation

70

5.4

Proof of Main Theorem

Based on the above preparation, Theorem 5.2 and 5.3 are proved.
The structure of the proof of Main Theorem is as follows. First, we summarize terms in∥AB − A0 B0 ∥2 and consider degeneration of a polynomial ideal. Second, we resolve thenon-negative restriction by variable transformations which are isomorphic maps. Third, we
verify that the problem can come down to finding the RLCT of reduced rank regression.
Lastly, we calculate the concrete value of the RLCT in each case.
Proof of Main Theorem. Because of Theorem 5.1 we only have to consider the analytic set
defined by{(A, B) | ∥AB − A0 B0 ∥2 = 0, A and B are stochastic matrices.}to determine the RLCT of LDA λ and its multiplicity m.
The first part is same as the first half of the proof of Appendix A in our previous research[38]. For the sake of self-containedness, we write down the process of developing the terms
in the above paper. Let ∼ be a binomial relation such that the functions K1 (w) and K2 (w)have same RLCT if K1 (w) ∼ K2 (w). Summarizing the terms, we have∥AB − A0 B0 ∥2 =

N M
−1 ∑
H
∑
∑

aik bkj −j=1 i=1 k=1

H0∑

a0ik b0kj +

N ∑
H
∑

aM k bkj −j=1 k=1

k=1

H0∑

a0M k b0kj .

k=1(5.24)Put
Kij :=

H∑

aik bkj −Lj :=a0ik b0kj ,(5.25)

a0M k b0kj ,(5.26)

k=1

k=1
H
∑

H0
∑

aM k bkj −

H0
∑
k=1

k=1

then we get
−1
N M
∑
∑

∥AB − A0 B0 ∥2 =j=1 i=1

2
Kij
+

N
∑

L2j .

j=1

∑M −1
∑H−1
∑M −1
Using aM k = 1 − i=1 aik , bHj = 1 − k=1 bkj , a0M k = 1 − i=1 a0ik , and b0H0 j =∑H0 −1 0
1 − k=1
bkj , we have
Kij =

H−1
∑

(aik − aiH )bkj −

H∑
0 −1

k=1

Lj = −M−1 H−1∑∑i=1 k=1(a0ik − a0iH0 )b0kj + (aiH − a0iH0 ),(5.27)

k=1(aik − aiH )bkj +

M
−1 H∑
0 −1
∑
i=1

k=1

0
(aik
− a0iH0 )b0kj −M−1∑(aiH − a0iH0 ),i=1(5.28)

5.4 Proof of Main Theorem

71

thus
L2j =

(M −1
∑

)2
Kij

.

i=1

Therefore∥AB − A0 B0 ∥ =
2

=

N M
−1
∑
∑

2
Kij
+

∑M −1
i=1

L2j

(5.29)

j=1 i=1

j=1

N M
−1
∑
∑

N
∑

(M −1
∑

j=1

i=1

2
Kij
+

j=1 i=1

Since the polynomial
have

N
∑

)2
Kij

.

(5.30)−1,N
Kij is contained in the ideal generated from (Kij )M
i=1,j=1 , we∥AB − A0 B0 ∥ ∼
2

N M
−1
∑
∑

2
Kij
,

j=1 i=1i.e.
∥AB − A0 B0 ∥2
{
}2
H∑
N M
−1 H−10 −1∑∑∑∼(aik − aiH )bkj −(a0ik − a0iH0 )b0kj + (aiH − a0iH0 )j=1 i=1

=

k=1

[ −1
−1 H∑
N M
0
∑
∑
j=1 i=1(5.31)(5.32)

k=1{(aik −aiH )bkj −(a0ik −a0iH0 )b0kj }+

k=1

H−1
∑

]2
(aik −aiH )bkj +(aiH −a0iH0 )

.

k=H0

(5.33)
aik = aik − aiH , k < H
Let ci = aiH − a0iH0 ,
bkj = bkj

(5.34)and put a0ik = a0ik − a0iH0 . Then we have∥AB − A0 B0 ∥2 ∼

N M
−1
∑
∑

{H −1
0
∑

j=1 i=1

k=1(aik bkj − a0ik b0kj ) +

H−1
∑

}2
aik bkj + ci.

(5.35)

k=H0

We had derived an upper bound of λ by using some inequalities of Frobenius norm and theexact value of λ in special cases [38]. However, in this paper, we use changes of variables
which resolve non-negative restrictions and find the RLCT in the all cases.
The transformation (5.34) resolves the non-negative restrictions of aik (k < H) and ci for
i = 1, . . . , M − 1. The changed variables aik (k < H) and ci can be negative. We call theChapter 5 Bayesian Generalization Error in Latent Dirichlet Allocation

72

determinant of the Jacobian matrix Jacobian for the sake of simplicity. The Jacobian of thetransformation (5.34) equals one.
k<H
aik = aik ,∑
∑H−1
H0 −1
0
0
Let xi = ci + k=1 (aik bk1 − aik bk1 ) + k=H0 aik bk1 ,(5.36)

bkj = bkj .
It is immediately derived that the Jacobian of this map is equal to one. About the transform(5.36), for j = 2, . . . , N , we haveH∑
0 −1(aik bkj − a0ik b0kj ) +

k=1

H−1
∑

aik bkj + ci

k=H0

= xi −

H∑
0 −1

H−1
∑

(aik bk1 − a0ik b0k1 ) −k=1

aik bk1 +

k=H0

H∑
0 −1

(aik bkj − a0ik b0kj ) +

k=1

H−1
∑

aik bkj .

k=H0(5.37)
Substituting this for∑H0 −1k=1(aik bkj − a0ik b0kj ) +∑H−1

k=H0 aik bkj + ci in Eq. (5.35), we have∥AB − A0 B0 ∥2
{ −1
}2
H−1
M
−1 H∑
0
∑
∑
(aik bk1 − a0ik b0k1 ) +aik bk1 + ci
∼
i=1

+

{ −1
N M
−1 H∑
0
∑
∑
j=2 i=1

=

M
−1
∑

x2i +

i=1

+

k=H0

k=1

(aik bkj − a0ik b0kj ) +

k=1
N M
−1
∑
∑

=

M
−1
∑
i=1

xi −

H∑
0 −1

x2i +

−1
N M
∑
∑

aik bkj + ci

0 0(aik bk1 − aikbk1 ) −

k=1

(aik bkj − a0ik b0kj ) +

k=1

}2

k=H0

{

j=2 i=1

H∑
0 −1

H−1
∑

H−1
∑

}2

H−1
∑

aik bk1

k=H0

aik bkj

k=H0

[
xi +

j=2 i=1H∑
0 −1{aik (bkj −bk1 )−a0ik (b0kj −b0k1 )}+

k=1

H−1
∑

]2
aik (bkj −bk1 )

.

k=H0

(5.38)
Put
gij =

H∑
0 −1

{aik (bkj − bk1 ) − a0ik (b0kj − b0k1 )} +

k=1

H−1
∑

aik (bkj − bk1 ).

k=H0

From Eq. (5.38), we have∥AB − A0 B0 ∥ ∼
2

M
−1
∑
i=1

x2i +

N M
−1
∑
∑
j=2 i=1(xi + gij )2 .

(5.39)

5.4 Proof of Main Theorem73

M −1,N−1
Let J be a polynomial ideal ⟨(xi )M
i=1 , (gij )i=1,j=2 ⟩. On account of xi + gij ∈ J, we have
M
−1
∑

x2i +

i=1

N M
−1
∑
∑

(xi + gij ) ∼2

j=2 i=1

M
−1
∑
i=1

x2i +

N M
−1
∑
∑

2
(x2i + gij),

j=2 i=1i.e.
∥AB − A0 B0 ∥2
∼

N M
−1
∑
∑

2
(x2i + gij)

j=2 i=1

∼

=

M
−1
∑

x2i +

N M
−1
∑
∑

2
gij

i=1

j=2 i=1

M
−1
∑

N M
−1
∑
∑

[H −1
0
∑

j=2 i=1

k=1

i=1

x2i +

{aik (bkj − bk1 ) − a0ik (b0kj − b0k1 )} +

]2

H−1
∑

aik (bkj − bk1 )

.

k=H0

(5.40)
aik = aik ,x = x ,iiLetbk1 = bk1 ,
bkj = bkj − bk1

k<H(5.41)
j > 1.

For k = 1, . . . , H − 1 and j = 2, . . . , N , non-negative restrictions of bkj can be resolved.
The Jacobian of the transformation (5.41) is one. Apply this map to Eq. (5.40) and put
b0kj = b0kj − b0k1 . Then, we have∥AB − A0 B0 ∥ ∼
2

=

N M
−1
∑
∑

{H −1
0
∑

i=1

j=2 i=1

k=1

M
−1
∑

N M
−1
∑
∑

(H−1
∑

j=2 i=1

k=1

M
−1
∑

i=1

x2i +

x2i +

(aik bkj − a0ik b0kj ) +

H−1
∑

}2
aik bkj

k=H0

aik bkj −

H∑
0 −1

)2

a0ik b0kj.

(5.42)

k=1

There are not bk1 (k = 1, . . . , H − 1) in the right hand side; thus, we can regard thenon-negative restrictions of the all variables are resolved after applying the transformation(5.41).
M −1,H−1
H−1,N −1
Real matrices U , V , U0 , and V0 are denoted by U := (uik )i=1,k=1, V := (vkl )k=1,l=1
,
M
−1,H
−1
H
−1,N
−1
0
0
U0 := (u0ik )i=1,k=10 , and V0 := (vkl
)k=1,l=1, respectively. Here, we have uik = aik ,0000
vkl = vk(j−1) = bkj , uik = aik , and vkl = vk(j−1) = b0kj for i = 1, . . . , M − 1, k =1, . . . , H − 1 and j = 2, . . . , N . Note that these U0 and V0 are same as in Eqs. (5.13)
and (5.14) because of the above transformations for the entries of A0 and B0 . Therefore,r := rank(U0 V0 ) is also equal to that of Main Theorem.

74

Chapter 5 Bayesian Generalization Error in Latent Dirichlet Allocation
Now, let us start coming down the problem from LDA to reduced rank regression.
∥U V − U0 V0 ∥2 =

=

N
−1 M
−1
∑
∑

(H−1
∑

l=1 i=1

k=1(
N M−1 H−1∑∑∑j=2 i=1

uik vkl −

H∑
0 −1

)2
0
u0ik vkl

k=1

aik bkj −

k=1

H∑
0 −1

)2
a0ik b0kj(5.43)

k=1

holds; thus, from Eq. (5.42) and (5.43), we have∥AB − A0 B0 ∥ ∼
2

M
−1
∑

x2i + ∥U V − U0 V0 ∥2 .

(5.44)

i=1

Let (λ1 , m1 ) and (λ2 , m2 ) be pairs of the RLCT and its multiplicity of the first and the second
term, respectively. There is no common variable between {(xi )Mi=1 } and {(U, V )}; hence, wehave
λ = λ1 + λ2 ,(5.45)

m = m1 + m2 − 1.

(5.46)

By simple calculation, λ1 = (M −1)/2 and m1 = 1 hold. Besides, the entries of the matrices
U and V can be real as well as non-negative. Thus, λ2 is the RLCT of non-restricted MF, i.e.
that of reduced rank regression [10]. The same is true for the multiplicity m2 . Therefore, we
obtainλ=M −1+ λMF (M − 1, N − 1, H − 1, r),2
m = mMF (M − 1, N − 1, H − 1, r),(5.47)(5.48)

where r = rank(U0 V0 ).
Finally, we concretely calculate λ and m. According to [10], the RLCT and its multiplicityof MF are as follows.
(1) If N + r + 1 ≦ M + H and M + r + 1 ≦ N + H and H + r + 1 ≦ M + N and
M + N + H + r + 1 is even (M + N + H + r is odd), thenλMF (M −1, N −1, H−1, r) =1
{2(H+r−1)(M +N −2) − (M −N )2 − (H+r−1)2 },8(5.49)mMF (M − 1, N − 1, H − 1, r) = 1.

(5.50)(2) Else if N + r + 1 ≦ M + H and M + r + 1 ≦ N + H and H + r + 1 ≦ M + N and
M + N + H + r + 1 is odd (M + N + H + r is even), thenλMF (M −1, N −1, H−1, r) =1
{2(H+r−1)(M +N −2) − (M −N )2 − (H+r−1)2 + 1},
8(5.51)mMF (M − 1, N − 1, H − 1, r) = 2.

(5.52)

5.4 Proof of Main Theorem

75(3) Else if M + H < N + r + 1, then1{(M − 1)(H − 1) + (N − 1)r − (H − 1)r}, (5.53)
2
mMF (M − 1, N − 1, H − 1, r) = 1.
(5.54)λMF (M − 1, N − 1, H − 1, r) =(4) Else if N + H < M + r + 1, then1{(N − 1)(H − 1) + (M − 1)r − (H − 1)r}, (5.55)
2
mMF (M − 1, N − 1, H − 1, r) = 1.
(5.56)λMF (M − 1, N − 1, H − 1, r) =(5) Else (i.e. M + N < H + r + 1), then1(M − 1)(N − 1),
2
mMF (M − 1, N − 1, H − 1, r) = 1.

λMF (M − 1, N − 1, H − 1, r) =(5.57)(5.58)

Since the multiplicity is clear, we find the RLCT. We develop the terms in each case by using
Eq. (5.47).
In the case (1), we haveM −1 1+ {2(H + r − 1)(M + N − 2) − (M − N )2 − (H + r − 1)2 }(5.59)
2
8
M −1 1=+ {2(H + r)(M + N ) − 2(M + N ) − 4(H + r) + 4
2
8(5.60)− (M − N )2 − (H + r)2 + 2(H + r) − 1}1
= {4M − 4 + 2(H+r)(M +N ) − (M −N )2 − (H+r)2 − 2(M +N ) − 2(H+r)+3}8(5.61)
1= {2(H + r)(M + N ) − (M − N )2 − (H + r)2 + 2(M + N ) − 2(H + r) − 1 − 4N }
8
(5.62)
N1(5.63)= {2(H + r + 1)(M + N ) − (M − N )2 − (H + r + 1)2 } − .
82

λ=In the case (2), by the same way as the case (1), we have
λ=1
N{2(H + r + 1)(M + N ) − (M − N )2 − (H + r + 1)2 + 1} − .
82(5.64)

Chapter 5 Bayesian Generalization Error in Latent Dirichlet Allocation76

In the case (3), we haveM −1 1+ {(M − 1)(H − 1) + (N − 1)r − (H − 1)r}2
2
M −1 1+ {M H − (M + H) + 1 + N r − r − Hr + r}=
2
2
1
= (M H − M − H + 1 + N r − Hr + M − 1 + N − N )21= {M H + 1 + N (r + 1) − H(r + 1) − N }21N= {M H + 1 + N (r + 1) − H(r + 1)} − .
22

λ=(5.65)(5.66)(5.67)(5.68)(5.69)In the case (4), by the same way as the case (3), we have
λ=N
1
{N H + 1 + M (r + 1) − H(r + 1)} − .
22(5.70)In the case (5), we haveM −1 1+ (M − 1)(N − 1)221
= (M − 1)N
2
1
N= MN − .
22

λ=(5.71)(5.72)(5.73)

From the above, Main Theorem is proved. Comparing the RLCT of MF [10], we also obtainλ = λMF (M, N, H, r + 1) −N.
2Therefore, we obtain Theorem 5.2.
□

5.5

Discussion

Here, we will discuss the results of this chapter from three points of view. After that, we will
describe the numerical behavior of the theoretical result by conducting numerical experiments.

5.5.1

Parameter Restriction

The RLCT of LDA can be represented by using that of MF. Namely, Main Theorem can beinterpreted as that the learning coeﬃcient of LDA is that of the unconstrained MF minus thepenalty due to the simplex constraint. In fact, it can be proved thatλ = λMF (M, N, H, r + 1) −N
2

(5.74)

5.5 Discussion
holds (see also the rigorous proof of Main Theorem in Appendix). The dimension of thestochastic matrix AB with the degrees of freedom is (M − 1)N = M N − N . The subtracted
N is the dimension of the parameter that can be uniquely determined from the parametersof the other (M − 1)N dimensions in the matrix AB. This part can be regarded as an
N-dimensional regular statistical model, whose RLCT is N/2. This is the reason of theabove statement. Note that Main Theorem and its proof are not trivial. A hermeneutic
explanation cannot be a mathematical proof. In addition, the actual parameter dimension is(M − 1)H + (H − 1)N = (M + N − 1)H − N because we have to consider the matrices
A and B rather than AB. We cannot reach the result of this paper simply by maintaining
consistency of the degrees of freedom. Algebraic geometrical methods are used to solve thisproblem in learning theory: what the learning coeﬃcient of LDA is.

5.5.2 Theoretical Application
Since LDA is a knowledge discovery method, marginal-likelihood-based model selection often
tends to be preferred. However, BIC [66] cannot be used for LDA because it is a singularstatistical model. Although Gibbs sampling is usually used for full Bayesian inference of
LDA, it is diﬃcult to achieve a tempered posterior distribution; thus, we need other Markov
chain Monte Carlo method (MCMC) to calculate WBIC [83] and WsBIC [43]. The resultof this study allows us to perform a rigorous model selection of LDA with sBIC [24], which
is MCMC-free. Even when the marginal likelihoods are computed directly by the exchange
Monte Carlo method, our result is useful for the design of the exchange probability [56].
Furthermore, it may be possible to evaluate how precise MCMC approximates the posterior
with the exact values of that [85, 43].
One may use BIC for model selection of LDA; however, using it causes that too small
models are chosen. This is because there exists a large diﬀerence in values and behaviors
between d/2 and λ. In a regular statistical model, the learning coeﬃcient is half of theparameter dimension d/2. In LDA, d/2 = (M + N − 1)H/2 − N/2 holds; hence, it linearly
increases as the number of topics H does. On the other hand, the RLCT of LDA λ does not.
In addition, λ is much smaller than d/2. Fig. 5.2a shows how the RLCT of LDA λ behaves
when the number of topics H increases, with λ-value in the vertical axis and H-value in thehorizontal axis. If λ was equal to d/2, then it would linearly increase (the square markers
dotted plot in Fig. 5.2a). However, in fact, λ is given by Main Theorem and its curve is
obviously non-linear (the circles dotted plot in Fig. 5.2a). Hence, their values and behaviors
are very diﬀerent. BIC is based on d/2 from the asymptotics of regular statistical models.
In contrast, the foundation of sBIC is singular learning theory; thus, it uses λ instead of d/2.
That is why sBIC is theoretically recommended for LDA.

5.5.3 Behavior of Learning Curve
We can draw the theoretical learning curve like the solid line in Fig. 5.2b, with E[Gn ]-value
in the vertical axis and n-value in the horizontal axis. We also namely draw a curve likethe dashed line in Fig. 5.2b. This dashed curve is not only an upper bound of the learning
curve of LDA in Bayesian inference but also a lower bound of that in maximum likelihood or
posterior estimation methods. Let GMAP
and µ be the generalization error and the learningn
coeﬃcient of LDA in maximum posterior (MAP) methods, respectively. This is well-defined,

77

Chapter 5 Bayesian Generalization Error in Latent Dirichlet Allocation

78

(a)

(b)

Fig. 5.2: (a) In this chapter, we give the exact value of the learning coeﬃcient of LDA λ.
The learning coeﬃcient is smaller than half of the parameter dimension d/2, since LDA is asingular statistical model. The dotted blue line drawn by the circles in this figure represents
the learning coeﬃcients of LDA when the number of topics H is increased. If LDA was aregular statistical model, its learning coeﬃcient would be the dotted yellow line drawn by thesquares. The behavior of them are so diﬀerent.
(b) This figure shows the theoretical learning curve of LDA and that of a regular statisticalmodel whose parameter dimension d is same as LDA. The former is the solid blue line andthe latter is the dashed yellow line. The vertical axis means the expected generalization error
E[Gn ] and the horizontal one is the sample size n. This is based on Theorem 3.4 and theexact value of λ which is clarified by our result.

] = µ/n + o(1/n) holds (see also Theorem 3.5). On the basis of the same priori.e. E[GMAPndistribution, Watanabe proved the following inequality [85]:λ < d/2 < µ.

(5.75)

This means E[GMAP] > E[Gn ] + o(1/n) and the leading term of these diﬀerence is (µ −nλ)/n > (d − 2λ)/2n. Owing to Main Theorem, we immediately have the exact value ofd − 2λ. Therefore, our result shows at least how much Bayesian inference improves thegeneralization performance of LDA compared to MAP method. If the prior distribution is a
uniform one, then µ equals the learning coeﬃcient of LDA in maximum likelihood estimation.
Hence, the above consideration can be applied to maximum likelihood estimation.

5.5.4

ExperimentNow, we run numerical experiments to check the behavior of Main Theorem when the sample
size is finite. Theorem 5.2 gives the exact asymptotic form of Bayesian generalization error in
LDA by using Theorem 3.4. We calculate the Bayesian generalization error in LDA by using
Gibbs sampling and compare the numerically-calculated RLCT with the theoretical one. Our
experimental approach and its description in this section is based on [33] since it also treats thenumerical experiment to compute the RLCT by using Gibbs sampling to verify the numerical5.5 Discussion

79

behavior of its theoretical result.
∫
Let Eθ be an expectation operator of the posterior: Eθ [·] = dθψ(θ|D)[·]. Let λ̂ be thenumerically calculated RLCT. The widely applicable information criterion (WAIC) [82] isdefined by the following random variable Wn :Wn = Tn + Vn /n,(5.76)

where Tn is the empirical loss and Vn is the functional variance:
Tn = −
Vn =

n
n
1∑
1∑
log p∗ (Xi ) = −log Eθ [p(Xi |θ)],
n i=1
n i=1

(5.77)

n [
n
] ∑
∑
2
Eθ [(log p(Xi |θ))2 ] − {Eθ [log p(Xi |θ)]} =Vθ [log p(Xi |θ)].
i=1(5.78)

i=1

Even if the posterior distribution cannot be approximated by any normal distribution (i.e.,the model is singular), the expected WAIC E[Wn ] is asymptotically equal to the expected
generalization loss E[Gn + S] [82];E[Gn + S] = E[Wn ] + O(1/n2 ).

(5.79)

Moreover, the generalization error and the WAIC error Wn − Sn have the same variance[82, 85]:
Gn + Wn − Sn = 2λ/n + Op (1/n).

(5.80)

We need to repeat the simulations to compute λ̂ to decrease the random eﬀect caused by Gn ,Wn and Sn . Thus, Eq. (5.80) is useful for computing λ̂ because the leading term 2λ/n isdeterministic, nevertheless the left hand side is probabilistic. This means that the needed
number of simulations D can be less than that in the case using λ ≈ nE[Gn ] from Theorem(3.4).
The method was as follows. First, the training data D was generated from the true distribution q(X|Z). Second, the posterior distribution was calculated by using Gibbs sampling[30] (see also Algorithm 5). Third, Gn and Wn − Sn were computed by using the training
data D and the artificial test data D∗ = (Xt∗ ) generated from q(X|Z). These three steps wererepeated and each value of n(Gn + Wn − Sn )/2 was saved. After all repetitions have been
completed, n(Gn + Wn − Sn )/2 was averaged over the simulations. This average was λ̂.
The pseudo-code is listed in Algorithm 4, where K is the sample size of the parameter
subject to the posterior and nT is the sample size of the synthesized test data. We used theprograming language named Julia [13] for this experiment.
We set M = 10, N = 5, H0 = 2, r = 1, n = 1000, and nT = 200n = 200000. To
examine the numerical behavior when the number of topics H ≧ H0 in the model increases,
we set H = 2, 3, 4, 5 and carried out experiments in each case. To decrease the probabilistic
eﬀect of Eq. (5.80), we conducted the simulations one hundred times, i.e. D = 100.
In the Gibbs sampling, we had to conduct a burn-in to decrease the eﬀect of the initial values
and thin the samples in order to break the correlations. The length of the burn-in was 10000,
while the length of the thinning was 20; thus, the sample sizes of the parameter was 10000 +80

Chapter 5 Bayesian Generalization Error in Latent Dirichlet Allocation
Algorithm 4 How to Compute λ̂Require: D: the number of simulations,
A0 : the true parameter matrix whose size is (M, H0 ),B0 : the true parameter matrix whose size is (H0 , N ),GS: the Gibbs sampling function whose return value consists of the samples from theposterior. See also Algorithm 5.
Ensure: The numerical computed RLCT λ̂.
Allocate an array Λ[D].
for d = 1 to D do
Generate D from the true distribution.
Allocate arrays A[M, H, K] and B[H, N, K].
Get A, B ← GS(D).
Generate D∗ from the
distribution.
∑true∑nq(X ∗ )
nT
Calculate Gn ≈ n1T t=1
log Eθ [p(Xt∗ |θ)] , Sn = − n1 i=1 log q(Xi ),t ∑∑nn
and Wn ≈ − n1 i=1 Eθ [p(Xi |θ)] + n1 i=1 Vθ [log p(Xi |θ)],∑K
1
where Eθ [f (θ)] ≈ K
k=1 f (θk ) and θk = (A[:, :, k], B[:, :, k]).
Save Λ[d] ← n(Gn + Wn − Sn )/2.
end for∑D
1
Calculate λ̂ = D
d=1 Λ[d].

20K = 50000 (K = 2000) and we used the (10000 + 20k)-th sample as the entry of A[:, :, k]and B[:, :, k] for k = 1 to K. The implementation is available at the following github page:https://github.com/chijan-nh/LearningCoefficient-RLCT-ofLDA-usingGS.
The experimental results are shown in Table 5.2 and visualized in Fig. 5.3. The symbol
λ denotes the exact value of the RLCT λ in Theorem 5.2. There are three columns foreach H, and each row contains the model settings (Settings), symbols of calculated values(RLCTs), and the theoretical or numerical values (Values). The experimental values have four
significant digits. The numerically-calculated RLCT λ̂ is an average of n(Gn + Wn − Sn )/2
obtained from each simulations; hence, we also show the standard deviation of it as the rightnext to the plus-minus sign ± in each setting.
As shown in Table 5.2 and Fig. 5.3, all numerically calculated values are nearly equal to thetheoretical RLCTs, i.e. these diﬀerences are suﬃciently smaller than the standard deviation
overall simulations. Note that the parameter dimensions (M − 1)H + (H − 1)N are 23, 37,
51 and 65 for H = 2, 3, 4 and 5; thus, we consider that the sample size n = 1000 is not an
asymptotic scale. Moreover, we also consider that it is natural to fix the sample size while thenumber of topics increases because we compare some models for a dataset (the sample size
is fixed) in practical situations. Although the sample size is finite (not an asymptotic scale)and fixed, the theoretical values are included in the 1-standard deviation ranges for each case.
Besides, because of Fig. 5.3, the standard deviations are suﬃciently small for the scale of the
RLCTs. Therefore, Theorem 5.2 is consistent with the experimental result.

5.5 Discussion

81

Algorithm 5 Gibbs Sampling for LDA
GS(D, α = 1H , β = 1N , K = K, burnin = 10000, thin = 20)
L(=n/N )Require: D = {(zl , xl )}l=1: the data where zl ∈Onehot(N ) and xl ∈Onehot(M ) area document and a word, respectively. When zlj =1, xl is the l-th word in the document j.
α ∈ RH>0 : the hyperparameter of the Dirichlet prior for the stochastic matrix A,N
β ∈ R>0 : the hyperparameter of the Dirichlet prior for the stochastic matrix B,
Dir(C|Γ): a Dirichlet distribution of a stochastic matrix whose h-th column is generated
by Dir(c|Γ[:, h]).
Ensure: Sampling stochastic matrices from the numerical posterior.
Let iter = butnin + thin ∗ K.
Allocate arrays A[M, H, K], A[M, H, iter], B[H, N, K] and B[H, N, iter].
# Inital sampling for A and B from the prior:
Generate A, B ∼ Dir(A|α)Dir(B|β).
# Sampling from the posterior:for k = 1 to iter do## Sampling the topic y.
Allocate an array y[L, H].
for l = 1 to L dofor h = 1 to H do ∑∑MNLet η[l, h] = exp( j=1 zlj (log bhj + i=1 xli log aih )).
∑H
Put η[l, h] ← η[l, h]/ h=1 η[l, h].
end for
Generate y[l, :] ∼ Cat(y|η[l, :]).
end for## Sampling the stochastic matrix A.
for h = 1 and i =∑1L to H and M do
Let α̂[i, h] = l=1 y[l, h]x[l, i] + αh .
end for
Generate A ∼ Dir(A|α̂).
Put A[M, H, k] ← A.
## Sampling the stochastic matrix B.
for j = 1 and h = 1 to N and H do
∑L
Let β̂[h, j] = l=1 z[l, j]y[l, h] + βj .
end for
Generate B ∼ Dir(B|β̂).
Put B[H, N, k] ← B.
end for# Burn-in and thinning.
for k = 1 to K do
A[M, H, k] ← A[M, H, burnin + thin ∗ k].
B[H, N, k] ← B[H, N, burnin + thin ∗ k].
end for
Return A, B.

82

Chapter 5 Bayesian Generalization Error in Latent Dirichlet Allocation
Table 5.2: Numerically-Calculated and Theoretical Values of RLCTs
Settings
H=2(M = 10, N = 5
H0 = 2, r = 1)
H=3
(M = 10, N = 5
H0 = 2, r = 1)
H=4
(M = 10, N = 5
H0 = 2, r = 1)
H=5
(M = 10, N = 5
H0 = 2, r = 1)

RLCTs
Theoretical: λNumerical: λ̂
Diﬀerence: |λ − λ̂|
Theoretical: λNumerical: λ̂
Diﬀerence: |λ − λ̂|
Theoretical: λNumerical: λ̂
Diﬀerence: |λ − λ̂|
Theoretical: λNumerical: λ̂
Diﬀerence: |λ − λ̂|Values
21/2
10.79 ± 0.8591
0.2901 ± 0.8591
12
12.25 ± 0.9510
0.2534 ± 0.9510
27/2
13.57 ± 1.036
0.07114 ± 1.036
15
14.80 ± 1.143
0.2049 ± 1.143

Fig. 5.3: This figure is drawn based on Table 5.2 and Theorem 5.2. It compares numericallycalculated RLCTs λ̂ (Numerical λˆ, as the dashed yellow line with the error bars) and theoreticalones λ (Theoretical λ, as the solid blue line) for H = 2, 3.4.5. The horizontal line means
the number of topics H and the vertical one is the numerically-calculated or theoretical valueof the RLCT. Each error bar of experimental results is the 1-standard deviation range. Theline of Numerical λˆ and that of Theoretical λ are very close and the standard deviations are
suﬃciently smaller than the scale of the RLCTs.

5.6 Conclusion

5.6 Conclusion
We determine the RLCT of LDA in general cases (Theorem 5.2) and describe theoreticalapplications to Bayesian inference. According to Theorem 5.3, Theorem 5.2 means that theexact asymptotic form of the generalization error and the free energy are theoretically derived
for LDA.
The RLCT of LDA monotonically increases when the number of topic increases (Fig 5.2a),although its representation has complex branches as Theorem 5.2. Besides, the RLCT saturates
to a constant which is not depend on the number of the model’s topics and the rank determinedby the true distribution. Hence, its behavior is non-linear and bounded. On the other hand, theparameter dimension d monotonically and linearly increases and it is not bounded. If LDA
was regular, its RLCT would be d/2. However, LDA is singular and its RLCT is less thand/2. The theoretical result shows how they are diﬀerent. Moreover, these diﬀerence d/2 − λgives a lower bound of the diﬀerence between MAP and Bayesian generalization error. We
numerically compute the RLCT when the sample size is finite (not asymptotic scale) and theexperimental result is consistent to the main result.
One of future works is clarifying the eﬀect of the hyperparameter to the generalization error
and the free energy. This is formulated to finding simultaneous resolution of singularities
when the prior distribution is a Dirichlet distribution. A density function of a non-uniform
Dirichlet distribution has zero or diverged points; thus, the hyperparameter (i.e. the parameterof the Dirichlet distribution) may aﬀect the RLCT of LDA.

83

85

Chapter 6

Conclusion
6.1 Conclusion
In this dissertation, the author aims to establish statistical learning theory when the parameterof the model region is restricted. As a foundation of that, we studied the asymptotic behaviorsof the Bayesian generalization errors in non-negative matrix factorization (NMF) and latent
Dirichlet allocation (LDA). NMF and LDA are two typical singular models whose parameterregions are restricted.
Singular learning theory is a mathematical foundation for statistical inference using singularmodels. It describes the asymptotic behaviors of the Bayesian generalization error and thefree energy by using algebraic geometry; a real log canonical threshold (RLCT) rules them.
Thus, we analyzed those of NMF and LDA in the framework of singular learning theory. The
contributions of the author’s works are as follows:• An upper bound of the RLCT of NMF is theoretically derived. The numerical experiment was consistent with the theoretical result and provided knowledge about thestability of learning.
– The upper bound provides theoretical upper bounds of the Bayesian generalization
error and the free energy in NMF.
– The bound includes the exact value in some cases and depends on the hyperparameter.
– The Bayesian generalization error and the free energy in NMF become strictly
larger than that of non-restricted matrix factorization when the entries of the true
parameter matrices are zero.
– The diﬀerence between the learning coeﬃcient of VBNMF and the RLCT of NMF
gives the variational approximation error. Thus, using the upper bound, we also
derived a lower bound of the variational approximation error in NMF.
– This theoretical analysis for NMF shows a phase transition structure; there isa critical line of hyperparameters. The Bayesian generalization error and thevariational approximation error drastically changes beyond that line.
– The phase transition line we found is diﬀerent from that of variational Bayesian
NMF. That is because the variational posterior of NMF is diﬀerent from theBayesian posterior distribution of NMF.
• The exact value of the RLCT is mathematically determined in general cases for LDA.
The numerical experiment was consistent with the theoretical result.

Chapter 6 Conclusion

86– The RLCT provides the exact asymptotic form of the Bayesian generalization error
and the free energy in LDA.
– We clarify that the RLCT of LDA is equal to that of stochastic matrix factorization(SMF). SMF is one of the matrix factorization models whose parameter region isrestricted. By using this fact, we obtain the exact value of the RLCT.
– SMF is NMF whose matrices have columns in simplexes. Thus, the main resultshows the eﬀect of the simplex restriction in matrix factorization.
– When the number of topics increases, the RLCT monotonically and non-linearly
grows but is bounded, whereas the parameter dimension linearly does and is notbounded.
– This mathematical study for LDA shows the RLCT of LDA is much smaller thanthat of a regular model whose parameter dimension is the same as LDA. This factprovides at least how much the generalization performance improves when we use
Bayesian inference.

6.2

Future Work

The future research topics for establishing singular learning theory of parameter-restricted
models include the followings. There are two layers of topics: general theory and problems
with concrete models such as NMF and LDA.
First, we briefly mention the research topics for creating a general theory when there exists
a constraint for the parameter region. Even if the sample is i.i.d. from the data-generating
distribution, the eﬀect to the generalization error caused by parameter restriction has not yet
been clarified. In general, like Lemma 4.6, when the model is the same, the support of theprior distribution is larger, the RLCT is smaller. This fact causes that parameter restriction
which does not decrease the dimension, such as non-negative restriction in NMF, increases
the Bayesian generalization error. However, this is only a qualitative evaluation. Thus, there
has not been a quantitative evaluation method. Besides, the restriction that decreases thedimension, such as simplex restriction in LDA and SMF, does not generally increase the
RLCT. In fact, in SMF, there are simplex restrictions to the parameter region and the RLCTof SMF is smaller than that of non-restricted matrix factorization.
Because of the diﬃculty in the general case, a possible policy is to construct a theory
for each class of algebraic varieties that characterize RLCTs. Vandermonde matrix typesingularities, which give RLCTs of mixture distribution models and neural networks, are
known and these RLCTs have been studied [5, 9, 8]. However, the quantitative eﬀect of theconstraints to these parameter regions is still unknown except for the RLCT of the Poisson
mixture model [64]; the average parameter of each component is non-negative. In this thesis,
the strategy was to deal with singularities formed by the squared error of matrix factorization.
For these singularities, we considered two typical constraints (NMF and LDA). Future topics
include other restrictions such as semi non-negative [22] and convolutive non-negative [67].
Second, we describe future works for singular learning theory of NMF and LDA.
For NMF, clarifying the exact value of the RLCT is considered one of the most significant
but demanding open problems. The non-negative rank of a non-negative matrix is generally
diﬀerent from the usual rank of that; thus, we cannot apply the proof of Aoyagi’s result [10]to the RLCT of NMF. Moreover, analysis and comparison of non-parametric Bayesian NMF

6.2 Future Work[41] is also included in future tasks.
For LDA, determining the RLCT when the prior is a Dirichlet distribution is considered as
one of the most important future research directions. This situation is formulated to finding
simultaneous resolution of singularities. A density function of a non-uniform Dirichlet distribution has zero or diverged points; thus, the hyperparameter (the parameter of the Dirichlet
distribution) may aﬀect the RLCT of LDA. If we conducted the simultaneous resolution,
the eﬀect of the hyperparameter on the RLCT would be clarified. Moreover, variational approximation error could be quantitatively evaluated in the same way as Theorem 4.4 since thevariational free energy and its learning coeﬃcient in LDA has been determined [59]. Analysis
and comparison of non-parametric LDA [99] is also considered as one of the future research
topics.
Note that the extracted clusters strongly depend on the hyperparameter in both nonparametric NMF and LDA. Hence, we have to reveal the free energy and the generalization
error as functions of the hyperparameter. For singular models, the eﬀect of the hyperparameter is complex since they have phase transition structures. Besides, in the frameworkof non-parametric Bayesian inference, the model size such as the non-negative rank and thenumber of topics is referred to as infinity; thus, we may not be able to directly apply singularlearning theory to the analysis of non-parametric NMF and LDA. Hence, we consider that theanalysis of non-parametric cases is challenging.

87

89

Appendix A

Questions and Answers in Defense
In this chapter, the author writes the questions and the answers at the defense of this dissertation. The defenses were held on 2021/02/22 (pre-examination) and 2021/03/29 (publicpresentation).

A.1 Singular Learning Theory
Q. 1 What is the benefit of determining RLCTs of statistical models in the situation ofreal data analysis?
（a）One of the most important benefits is that we can use sBIC [24] (see also Sec. 3.4)
to select the appropriate model in knowledge discovery sense. The criterion sBIC
uses the theoretical value of the RLCT and its multiplicity. If they are unknown,
WsBIC [43] is one of the choices because it estimates the RLCT. However, thecalculation cost of estimating the RLCT is very high and this is just estimation:
we cannot obtain the exact value. Thus, clarifying the RLCT and its multiplicity
is useful for model selection.
Furthermore, a tuning method for the inversed temperature parameter in exchange
Monte Carlo has been proposed, which uses the theoretical RLCT [56]. It has
been also considered that the correctness of MCMC is verified by comparing theestimated RLCT and the theoretical one [85, 43]. We can estimate the suﬃcient
sample size n∗ such that E[Gn ] < ε by using λ/n∗ < ε ⇔ n∗ > λ/ε if thetheoretical value of the RLCT λ is clarified. This fact means that there is apotential application to propose how much data should be collected.
Q. 2 Let n be the sample size. Can singular learning theory treat the limit of the model
size D such that D, n → ∞ and D/n = γ > 0? For example, D is the size of thedata matrix in NMF or the number of topics in LDA.
（a）At present, no it cannot. In singular learning theory, the model and the prior are
fixed to the sample. However, we can consider the limit of D → ∞ in the RLCT.
For example, in LDA, when the number of topics H tends to the infinity, the RLCTof LDA λ converges to a constant:lim λ = (M N − N )/2,

H→∞

where M and N are the number of vocabulary and that of documents, respectively(see Chap. 5).

Appendix A Questions and Answers in Defense

90

Q. 3 Does singular learning theory treat only the case that the parameter region is not
restricted?
（a）Singular learning theory can treat not only the case that the parameter region is
not restricted but also that is constrained, even if the model is singular, a map froma parameter to a probability distribution is not injective and its likelihood andposterior cannot be approximated by any normal distribution. However, the eﬀectof the parameter restriction on the generalization error has yet been unknown.
Thus, this dissertation aims to construct statistical learning theory for parameterrestricted singular models. As a foundation of it, the author studies two typicalmodels: NMF and LDA.
Q. 4 Is it really that the parameter restriction increases the generalization error? Thequestioner guesses that the constraint seems to make the estimation result fit agood area.
（a）In Chap. 4, it is mathematically proved that the non-negative restriction strictlyincreases the RLCT in some cases. Intuitively, this is because the constraint makes
the parameter space narrow without decreasing in dimension, which eliminates
regions that are more amenable to generalization.
Q. 5 Some models are often trained with a penalty term to its loss function to relax
their non-identifiability. Can we consider a singular model as a regular model byapplying norm constraints with a penalty?
（a）No, we cannot. When the loss function corresponds to the likelihood of the model,
the penalty is referred to as the prior distribution. For example, the L2 penalty
to the negative log-likelihood function is equivalent to the standard normal prior:log N (θ|0, 1) ∝ ∥θ∥2 . Although the model p(x∥θ) and the prior φ(θ) characterize
the singularity of the model, their contributions are significantly diﬀerent. The
scale of the eﬀect from the model is larger than that from the prior. Typically, theaverage error function is determined by the model and its RLCT rules the behaviorof the free energy and the generalization error (see Chap. 3).
Q. 6 How do we derive a non-trivial upper bound of the RLCT?
（a）We find an upper bound U (θ) of the average error function K(θ). It is desiredthat the RLCT of the bounding function U (θ) can be easily calculated. If we candetermine the RLCT λ of U (θ), we obtain an upper bound of the RLCT λ ofK(θ) by using Proposition 3.2. This proposition means that the RLCT is orderisomorphic; thus, K(θ) ≦ U (θ) ⇒ λ ≦ λ is derived. The tighter U (θ) is, thetighter λ is. Let d be the parameter dimension. If the prior is strictly positive and∑dbounded, a trivial bounding function is U (θ) ∝ ∥θ∥2 = i=1 θi2 : λ = d/2. If
the model is regular, λ = d/2. Thus, if the found bound satisfies λ > d/2, theupper bound is vacuous.

A.2 Main Results
Q. 1 In the work of NMF, where does the non-negative restriction essentially aﬀect themain result?
（a）The non-negative restriction is used to prove Lemmas 4.2, 4.3 and 4.4. In partic-

A.2 Main Resultsular, when the true non-negative rank is zero (H0 = 0, Lemma 4.2), the RLCT of
NMF is strictly larger than that of non-restricted matrix factorization. By using
these lemmas, the upper bound is derived for the general case (Theorem 4.2), andthe bound immediately provides the upper bound of the Bayesian generalization
error and the free energy in NMF (Theorem 4.3).
Q. 2 In the work of NMF, what is the novel point of this research comparing with theprior study by Kohjima [49]?
（a）This work is a theoretical analysis for inference methods diﬀerent from [49]. Kohjima’s previous work [49] has determined the exact learning coeﬃcient of thevariational free energy in NMF (see also Theorem 4.1). In this dissertation, theauthor treats Bayesian NMF and evaluates the Bayesian generalization error andthe free energy in NMF by analyzing the RLCT. The RLCT rules the asymptotic behavior of the Bayesian generalization error and the free energy. On theother hand, the RLCT is not equal to (is smaller than) the learning coeﬃcient of
the variational free energy and the coeﬃcient does not dominate the variational
generalization error. Besides, with merging this study and Kohjima’s result, thevariational approximation error is also theoretically evaluated: a lower bound of
it is derived (Theorem 4.4).
Q. 3 In the work of NMF, how do we understand the fact that the variational posterior
is essentially diﬀerent from the true one, in particular when we carry out real
data analysis (like making a predictive model)?
（a）Let θ = (θ1 , θ2 ) be a parameter of the model. One might think that the variational
posterior ψ1 (θ1 )ψ2 (θ2 ) is equivalent to the true one ψ(θ1 , θ2 ), i.e. variational
inference is faster and equally accurate with comparing to MCMC (full-Bayesianinference). However, this is false if the parameters θ1 and θ2 are not independent.
Thus, if we want to realize the posterior, then the variational inference can beinappropriate (MCMC is recommended). If we want a fast algorithm rather thanan accurate one, then the variational inference can be helpful but it cannot realize
the true posterior.
Q. 4 In the work of NMF, when M = N , can it be immediately proved that thephase transition lines of Bayesian inference and variational Bayes method areorthogonal?
（a）Yes. Let l1 and l2 be the phase transition lines of Bayesian NMF and VBNMF in
Theorem 4.4 (see also Fig. 4.1), respectively. In the (ϕU , ϕV )-plane, l1 and l2 areparameterized as M ϕU = N ϕV and M ϕU + N ϕV = (M + N )/2, respectively.
Because of M = N , we have l1 : ϕV = ϕU and l2 : ϕV = −ϕU + 1. The slopeof l1 is 1 and that of l2 is −1; thus, they are orthogonal. Note that they are neither
orthogonal nor parallel if M and N (M ̸= N ) are finite but they are formallyparallel when either M → ∞ ∧ N < ∞ or N → ∞ ∧ M < ∞.
Q. 5 In the work of NMF, what order is the term “phase transition”?
（a）The asymptotic (variational) free energy is dominated by the learning coeﬃcient
who has a phase transition line of the hyperparameters (ϕU , ϕV ). On the phasetransition lines l1 and l2 described above, the learning coeﬃcient is continuous
but not diﬀerentiable. Hence, this is a second-order phase transition.
Q. 6 In the work of LDA, why does not the author analyze the variational approximation error by the same way of NMF?

91

92

Appendix A Questions and Answers in Defense（a）Actually, the variational free energy in LDA is theoretically analyzed and its

learning coeﬃcient is determined [59]. However, the eﬀect of the hyperparametersof the Dirichlet prior distribution has not been clarified. Therefore, at this point,
it is diﬃcult to compare the variational free energy and the true one in LDA andto evaluate the variational approximation error. The author considers that this is
an important future task (see Sec. 6.2).
Q. 7 In the work of LDA, the main result immediately derives at least how much
Bayesian inference improves the generalization performance compared to maximum a posteriori (MAP) method and maximum likelihood (ML) method. On theother hand, what about the computational cost?
（a）In general, MAP and ML methods suﬀer from non-identifiability and hardly reachthe global minima for singular models. They are numerically unstable and hard toconverge. In this sense, their computational cost is not low although they carry out
point estimations. Moreover, like Fig. 5.2a, the gap between half of the parameter
dimension d/2 and the RLCT λ is large, and d/2 − λ gives a lower bound of thediﬀerence of the generalization errors of MAP (or ML) and Bayesian inference.
Thus, even if the computational cost is reduced by avoiding MCMC and using
MAP method, the generalization error will be so large that it is not worth theadvantage of reducing the cost.
Q. 8 Both the work of NMF and LDA, to prove the main results, do the prior distributions have to be set to specific distributions (gamma and Dirichlet)?
（a）For the main theorem of NMF, the distribution form is important; however, it is
not limited to the gamma distribution. The gamma distribution Gam(w|ϕ, θ) =θϕ
ϕ−1 −θwe
can become zero or unbounded because of the power term wϕ−1
Γ(θ) wand this term aﬀects the upper bound (see Theorem 4.2). Therefore, if the prior
has the same power term and the other terms are positive and bounded, the main
theorem holds.
（b）For the main result of LDA, the prior must not aﬀect the RLCT; thus, it must be
positive and bounded. In other words, if the prior is positive and bounded, thedistribution form of it is not limited to the uniform Dirichlet distribution (like theGaussian distribution on the simplex).
Q. 9 What is the correspondence between the author’s publications and the theorems
in this dissertation?
（a）The author’s peer-reviewed publications which construct this dissertation are as
follows: [37], [36], [33], [38] and [34] (see also Bibliography and List of Publications).
（b）The two journal papers [37, 33] and one international conference paper [36]correspond to the study of NMF (Chap. 4). The paper [37] is the first study for the
RLCT of NMF. In [37], Lemmas 4.2 and 4.3 was proved in the case ϕU = ϕV = 1.
Besides, by using these lemmas, an upper bound of the RLCT was also derived.
In the paper [36], Lemma 4.4 and a tighter upper bound was obtained and thebound is equal to the upper bound in Theorem 4.2 in the case ϕU = ϕV = 1.
In [33], the eﬀect of the hyperparameter was clarified, i.e. the case ϕU > 0 and
ϕV > 0 was treated. Moreover, a lower bound of the variational approximation
error in NMF was also derived (Theorem 4.4. Rigorously speaking, [33] shows aA.2 Main Results
looser upper bound of the RLCT than that in this dissertation because it could not
be proved that Lemma 4.4 for general hyperparameter ϕU > 0 and ϕV > 0 at thattime. In fact, Lemma 4.4 is true in that case; thus, Theorem 4.2 is proved in thisdissertation.
（c）The two journal papers [38, 34] are corresponding to the work of LDA (Chap. 5).
The paper [38] is the first paper for the RLCT of LDA. In [38], it was proved that
the RLCT of stochastic matrix factorization is equal to that of LDA (Theorem 5.1.
By using this relation, an upper bound of the RLCT of LDA was also derived. In
the paper [34], the exact value of the RLCT of LDA was determined (Theorem
5.2).

93

95

Acknowledgement
I would like to express my great appreciation for the following help, supports, and encouragement. I had carried out this research under the laboratory of Prof. Sumio Watanabe,
Department of Mathematical and Computing Science, Tokyo Institute of Technology. Prof.
Watanabe taught me the basics of statistical learning theory and algebraic geometry and has
widely and profoundly supervised me studying since I was a master’s student in 2016. I
learned the illustriousness of investigating the nature and establishing mathematical theory
from Prof. Watanabe. From research to real life and career counseling, Prof. Watanabe encouraged me with gentle words. Some of the travels and papers for this research were partially
supported by the Ministry of Education, Science, Sports and Culture in Japan, Grant-in-Aid
for Scientific Research 15K00331 from Prof. Watanabe. Besides, this study was also partially supported by the expenses of NTT DATA Mathematical Systems Inc. where I haveworked as a data scientist. I thank them for their patience. Former and present members of
Watanabe Group, Ph.D. Fumito Nakamura, Ph.D. Masahiro Kohjima, Ph.D. Natsuki Kariya,
M.S. Natsuki Kohda, M.Eng. Toshihiko Kondo, M.S. Shoji Sugai, M.Eng. Rito Takeuchi,
M.S. Hikaru Fujimura, M.S. Kohya Ohashi, M.S. Takashi Okano, M.S. Satoshi Kataoka,
M.S. Kenichiro Sato, M.S. Raiki Tanaka, M.S. Atsuyoshi Muta, M.S. Shuya Nagayasu, M.S.
Toshiki Mikami, M.S. Kento Yoshimura, B.S. Takumi Watanabe, B.S. Tatsuyuki Baba, B.S.
Rui Hirose, B.S. Satoshi Nakagawa, B.S. Joe Hirose, B.S. Nozomi Maki, and B.S. Ko Takemura (degree, grade, and alphabetical order), lively discussed this research and other science
topics with me. Prof. Mathias Drton, Department of Mathematics, Technical Universityof Munich, kindly answered my questions about sBIC and encouraged this research. Prof.
Miki Aoyagi, Department of Mathematics, Nihon University, taught me in detail about the
Vandermonde matrix type singularity theory, which is one of the most advanced topics insingular learning theory. Prof. Kota Matsui, Department of Biostatistics, Nagoya University,
invited me to StatsML Symposium’20. I was blessed with the valuable opportunity having
discussed and introduced the singular learning theory and the contents of my research in that
invited talk. Prof. Takafumi Kanamori, Prof. Misako Takayasu, Prof. Makoto Yamashita,
and Prof. Hanna Sumita carefully reviewed this dissertation. Lastly, I greatly appreciate myparents for warmly raising and supporting me, my friends for kindly encouraging each other,and all the other people who took care of me though I cannot write in this chapter.

97

Acknowledgement (in Japanese)
本研究は様々な支援を受けました．指導教官である渡邊澄夫教授からは統計的学習理論及
び代数幾何学・特異点論の基礎から具体的かつ詳細な研究指導に至るまで広く深くご教授いただき，実世界を探索し数理科学の理論を構築することの尊さを学ぶことができました．修士課程から博士課程まで，途中就職により研究室に在籍していないときでも，渡邊澄夫教授は筆者の研究相談に乗ってくださいました．研究から実生活・進路相談に至るまで，渡邊澄夫教授は優しい言葉をかけて私を励ましてくださいました．一部の論文投稿や学会のための渡航などについても，渡邊澄夫教授から科学研究費補助金 15K00331 の援助を受けました．また，筆者が所属する株式会社 NTT データ数理システムの経費も，一部の論文投稿や学会参加の費用として用いました．研究室で時を同じくして過ごした中村文
士博士，幸島匡宏博士，仮屋夏樹博士，香田夏輝修士，近藤稔彦修士，須貝将士修士，竹内理人修士，大橋耕也修士，岡野高志修士，藤村光修士，片岡諭史修士，佐藤件一郎修士，
田中来輝修士，永安修也修士，牟田篤兄修士，三上敬生修士，吉村健斗修士，渡邉匠学士，
中川哲学士，馬場達之学士，廣瀬瑠伊学士，竹村航学士，広瀬青学士，槇望学士（学位・学年・あいうえお順）は様々なディスカッションに付き合ってくださいました．ミュンヘ
ン工科大学の Mathias Drton 教授は sBIC に関する筆者の拙い質問に丁寧に答えてくださり，本研究を応援してくださいました．日本大学の青柳美樹教授からは特異学習理論における最先端理論である Vandermonde 行列型特異点論について詳細にご教示いただきました．名古屋大学の松井孝太講師からは「第 5 回統計・機械学習若手シンポジウム」に招待いただき，特異学習理論と本研究内容を紹介してディスカッションを行う貴重な機会が得られました．金森敬文教授，高安美佐子教授，山下真教授，澄田範奈講師には本博士論文を審査いただきました．最後に，筆者を育ててくださった両親，博士課程期間中に励まし合うことができた友人たち，その他本章に書き記しきれないお世話になりましたすべての方々に感謝します．

99

Bibliography[1] Christopher P. Adams. Finite mixture models with one exclusion restriction. The
Econometrics Journal, Vol. 19, No. 2, pp. 150–165, 2016.
[2] Christopher P. Adams. Stochastic matrix factorization. SSRN Electronic Journal, pp.
1–24, 2016. Available at SSRN: https://ssrn.com/abstract=2840852.
[3] Hirotogu Akaike. Information theory and an extension of the maximum likelihood
principle. In Selected papers of hirotugu akaike, pp. 199–213. Springer, 1998.
[4] Hirotugu Akaike. A new look at the statistical model identification. IEEE transactions
on automatic control, Vol. 19, No. 6, pp. 716–723, 1974.
[5] Miki Aoyagi. A bayesian learning coeﬃcient of generalization error and vandermondematrix-type singularities. Communications in Statistics―Theory and Methods, Vol. 39,
No. 15, pp. 2667–2687, 2010.
[6] Miki Aoyagi. Stochastic complexity and generalization error of a restricted boltzmann
machine in bayesian estimation. Journal of Machine Learning Research, Vol. 11, No.
Apr, pp. 1243–1272, 2010.
[7] Miki Aoyagi. Learning coeﬃcient in bayesian estimation of restricted boltzmann
machine. Journal of Algebraic Statistics, Vol. 4, No. 1, pp. 30–57, 2013.
[8] Miki Aoyagi. Learning coeﬃcient of vandermonde matrix-type singularities in model
selection. Entropy, Vol. 21, No. 6, p. 561, 2019.
[9] Miki Aoyagi and Kenji Nagata. Learning coeﬃcient of generalization error in bayesianestimation and vandermonde matrix-type singularity. Neural Computation, Vol. 24,
No. 6, pp. 1569–1610, 2012.
[10] Miki Aoyagi and Sumio Watanabe. Stochastic complexities of reduced rank regression
in bayesian estimation. Neural Networks, Vol. 18, No. 7, pp. 924–933, 2005.
[11] Michael Francis Atiyah. Resolution of singularities and division of distributions.
Communications on pure and applied mathematics, Vol. 23, No. 2, pp. 145–150, 1970.
[12] Joseph Bernstein. The analytic continuation of generalized functions with respect to aparameter. Funktsional’nyi Analiz i ego Prilozheniya, Vol. 6, No. 4, pp. 26–40, 1972.
[13] Jeﬀ Bezanson, Alan Edelman, Stefan Karpinski, and Viral B Shah. Julia: A fresh
approach to numerical computing. SIAM review, Vol. 59, No. 1, pp. 65–98, 2017.
[14] Christopher M Bishop. Pattern recognition and machine learning. springer, 2006.
[15] David M Blei, Andrew Y Ng, and Michael I Jordan. Latent dirichlet allocation. Journalof machine Learning research, Vol. 3, No. Jan, pp. 993–1022, 2003.
[16] Jesús Bobadilla, Rodolfo Bojorque, Antonio Hernando Esteban, and Remigio Hurtado.
Recommender systems clustering using bayesian non negative matrix factorization.
IEEE Access, Vol. 6, pp. 3549–3564, 2018.
[17] Bob Carpenter, Andrew Gelman, Matthew D Hoﬀman, Daniel Lee, Ben Goodrich,Michael Betancourt, Marcus A Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell.

100

Bibliography
Stan: a probabilistic programming language. Grantee Submission, Vol. 76, No. 1, pp.
1–32, 2017.
[18] Ali T. Cemgil. Bayesian inference in non-negative matrix factorisation models. Computational Intelligence and Neuroscience, Vol. 2009, No. 4, p. 17, 2009. Article ID
785152.
[19] Joel E. Cohen and Uriel G. Rothblum. Nonnegative ranks, decompositions, and factorizations of nonnegative matrices. Linear Algebra and Its Applications, Vol. 190, pp.
149–168, 1993.
[20] Peter D Congdon. Bayesian hierarchical models: with applications using R. CRC
Press, 2019.
[21] Dingxiong Deng, Cyrus Shahabi, Ugur Demiryurek, Linhong Zhu, Rose Yu, and Yan
Liu. Latent space model for road networks to predict time-varying traﬃc. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery andData Mining, pp. 1525–1534, 2016.
[22] Chris HQ Ding, Tao Li, and Michael I Jordan. Convex and semi-nonnegative matrix
factorizations. IEEE transactions on pattern analysis and machine intelligence, Vol. 32,
No. 1, pp. 45–55, 2008.
[23] Mathias Drton, Shaowei Lin, Luca Weihs, Piotr Zwiernik, et al. Marginal likelihood
and model selection for gaussian latent tree and forest models. Bernoulli, Vol. 23,
No. 2, pp. 1202–1232, 2017.
[24] Mathias Drton and Martyn Plummer. A bayesian information criterion for singularmodels. Journal of the Royal Statistical Society Series B, Vol. 79, pp. 323–380, 2017.
with discussion.
[25] Ilenia Epifani, Steven N MacEachern, Mario Peruggia, et al. Case-deletion importancesampling estimators: Central limit theorems and related results. Electronic Journal of
Statistics, Vol. 2, pp. 774–806, 2008.
[26] Cédric Févotte, Nancy Bertin, and Jean-Louis Durrieu. Nonnegative matrix factorization with the itakura-saito divergence: With application to music analysis. Neuralcomputation, Vol. 21, No. 3, pp. 793–830, 2009.
[27] Lorenzo Finesso and Peter Spreij. Nonnegative matrix factorization and i-divergence
alternating minimization. Linear Algebra and its Applications, Vol. 416, No. 2-3, pp.
270–287, 2006.
[28] Andrew Gelman, Daniel Lee, and Jiqiang Guo. Stan: A probabilistic programming language for bayesian inference and optimization. Journal of Educational and Behavioral
Statistics, Vol. 40, No. 5, pp. 530–543, 2015.
[29] Daniel Gildea and Thomas Hofmann. Topic-based language models using em. In Sixth
European Conference on Speech Communication and Technology, 1999.
[30] Thomas L Griﬃths and Mark Steyvers. Finding scientific topics. Proceedings of the
National academy of Sciences, Vol. 101, No. suppl 1, pp. 5228–5235, 2004.
[31] Zellig S. Harris. Distributional structure. WORD, Vol. 10, No. 2-3, pp. 146–162, 1954.
[32] JA Hartigan. A failure of likelihood asymptotics for normal mixtures. In Proceedingsof the Berkeley Conference in Honor of J. Neyman and J. Kiefer, 1985, pp. 807–810,
1985.
[33] Naoki Hayashi. Variational approximation error in non-negative matrix factorization.
Neural Networks, Vol. 126, pp. 65–75, 2020.
[34] Naoki Hayashi. The exact asymptotic form of bayesian generalization error in latent

101
dirichlet allocation. Neural Networks, Vol. 137, pp. 127–137, 2021.
[35] Naoki Hayashi and Fumito Nakamura. Experimental analysis of variational bayesian
method in model selection of gaussian mixture model by singular bayesian information
criterion. In Information-Based Induction Sciences and Machine Learning (IBISML),
IEICE Technical Report, Vol. 117, pp. 19–26, 9 2017. in Japanese.
[36] Naoki Hayashi and Sumio Watanabe. Tighter upper bound of real log canonicalthreshold of non-negative matrix factorization and its application to bayesian inference.
In IEEE Symposium Series on Computational Intelligence (IEEE SSCI), pp. 718–725,
11 2017.
[37] Naoki Hayashi and Sumio Watanabe. Upper bound of bayesian generalization error in
non-negative matrix factorization. Neurocomputing, Vol. 266C, No. 29 November, pp.
21–28, 2017.
[38] Naoki Hayashi and Sumio Watanabe. Asymptotic bayesian generalization error in latent
dirichlet allocation and stochastic matrix factorization. SN Computer Science, Vol. 1,
No. 2, pp. 1–22, 2020.
[39] Heisuke Hironaka. Resolution of singularities of an algbraic variety over a field of
characteristic zero. Annals of Mathematics, Vol. 79, pp. 109–326, 1964.
[40] N Thompson Hobbs and Mevin B Hooten. Bayesian models: a statistical primer for
ecologists. Princeton University Press, 2015.
[41] Matthew D Hoﬀman, David M Blei, and Perry R Cook. Bayesian nonparametric matrix
factorization for recorded music. In ICML, 2010.
[42] Tikara Hoshino, Kazuho Watanabe, and Sumio Watanabe. Stochastic complexity ofvariational bayesian hidden markov models. In International Joint Conference on
Neural Networks, Vol. 2, pp. 1114–1119 vol. 2, July 2005.
[43] Toru Imai. Estimating real log canonical thresholds. arXiv preprint arXiv:1906.01341,
2019.
[44] Toru Imai. On the overestimation of widely applicable bayesian information criterion.
arXiv preprint arXiv:1908.10572, 2019.
[45] F. Itakura and S. Saito. Analysis synthesis telephony based on the maximum likelihood
method. In In Proc. 6th of the International Congress on Acoustics, 1968.
[46] Hyunsoo Kim and Haesun Park. Sparse non-negative matrix factorizations via alternating non-negativity-constrained least squares for microarray data analysis. Bioinformatics, Vol. 23, No. 12, pp. 1495–1502, 2007. doi:10.1093/bioinformatics/btm134.
PMID 17483501.
[47] Masahiro Kohjima, Tatsushi Matsubayashi, and Hiroshi Sawada. Probabilistic nonnegative inconsistent-resolution matrices factorization. In Proceeding of CIKM ’15
Proceedings of the 24th ACM International on Conference on Information and Knowledge Management, Vol. 1, pp. 1855–1858, 2015.
[48] Masahiro Kohjima, Tatsushi Matsubayashi, and Hiroshi Sawada. Multiple data analysis
and non-negative matrix/tensor factorization [i]: multiple data analysis and its advances.
The journal of the Institute of Electronics, Information and Communication Engineers(IEICE), Vol. 99, No. 6, pp. 543–550, 2016. in Japanese.
[49] Masahiro Kohjima and Sumio Watanabe. Phase transition structure of variational
bayesian nonnegative matrix factorization. In International Conference on Artificial
Neural Networks, pp. 146–154. Springer, 2017.
[50] Ben Lambert. A student s guide to Bayesian statistics. Sage, 2018.

102

Bibliography[51] Daniel D. Lee and H. Sebastian Seung. Learning the parts of objects with nonnegative
matrix factorization. Nature, Vol. 401, pp. 788–791, 1999.
[52] Fei-Fei Li and Pietro Perona. A bayesian hierarchical model for learning naturalscene categories. In Proceedings of the 2005 IEEE Computer Society Conference on
Computer Vision and Pattern Recognition (CVPR’05) - Volume 2 - Volume 02, CVPR
’05, pp. 524–531, Washington, DC, USA, 2005. IEEE Computer Society.
[53] Ken Matsuda and Sumio Watanabe. Weighted blowup and its application to a mixtureof multinomial distributions. IEICE Transactions, Vol. J86-A, No. 3, pp. 278–287,
2003. in Japanese.
[54] Richard McElreath. Statistical Rethinking: A Bayesian Course with Examples in R and
Stan. CRC Press, 2016.
[55] Daniel Murfet, Susan Wei, Mingming Gong, Hui Li, Jesse Gell-Redman, and Thomas
Quella. Deep learning is singular, and that’s good. arXiv preprint arXiv:2010.11560,
2020.
[56] Kenji Nagata and Sumio Watanabe. Asymptotic behavior of exchange ratio in exchangemonte carlo method. Neural Networks, Vol. 21, No. 7, pp. 980–988, 2008.
[57] Shuya Nagayasu and Sumio Watanabe. Asymptotic behavior of free energy when
optimal probability distribution is not unique. arXiv preprint arXiv:2012.08338, 2020.
[58] Ryumei Nakada and Masaaki Imaizumi. Adaptive approximation and generalizationof deep neural network with intrinsic dimensionality. Journal of Machine Learning
Research, Vol. 21, No. 174, pp. 1–38, 2020.
[59] Shinichi Nakajima, Issei Sato, Masashi Sugiyama, Kazuho Watanabe, and Hiroko
Kobayashi. Analysis of variational bayesian latent dirichlet allocation: Weaker sparsity than map. In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K. Q.
Weinberger, editors, Advances in Neural Information Processing Systems, Vol. 27, pp.
1224–1232. Curran Associates, Inc., 2014.
[60] Shinichi Nakajima and Sumio Watanabe. Variational bayes solution of linear neuralnetworks and its generalization performance. Neural Computation, Vol. 19, No. 4, pp.
1112–1153, 2007.
[61] Pentti Paatero and Unto Tapper. Positive matrix factorization: A non-negative factor
model with optimal utilization of error estimates of data values. Environmetrics, Vol. 5,
No. 2, pp. 111–126, 1994. doi:10.1002/env.3170050203.
[62] Dmitry Rusakov and Dan Geiger. Asymptotic model selection for naive bayesian
networks. Journal of Machine Learning Research, Vol. 6, No. Jan, pp. 1–35, 2005.
[63] Ivan N Sanov. On the probability of large deviations of random variables. UnitedStates Air Force, Oﬃce of Scientific Research, 1958.
[64] Kenichiro Sato and Sumio Watanabe. Bayesian generalization error of poisson mixture
and simplex vandermonde matrix type singularity. arXiv preprint arXiv:1912.13289,
2019.
[65] Mikio Sato and Takuro Shintani. On zeta functions associated with prehomogeneous
vector spaces. Annals of Mathematics, pp. 131–170, 1974.
[66] Gideon Schwarz. Estimating the dimension of a model. The annals of statistics, Vol. 6,
No. 2, pp. 461–464, 1978.
[67] Paris Smaragdis. Non-negative matrix factor deconvolution; extraction of multiple
sound sources from monophonic inputs. In International Conference on Independent
Component Analysis and Signal Separation, pp. 494–499. Springer, 2004.

103[68] Wee Meng Soon, Hwee Tou Ng, and Daniel Chung Yong Lim. A machine learning
approach to coreference resolution of noun phrases. Computational Linguistics, Vol. 27,
No. 4, p. 521âĂŞ544, December 2001.
[69] Sayaka Suzuki, Souta Shina, and Miki Aoyagi. The improving method of singularbayesian information criterion by analyzing learning coeﬃcients. In InformationBased Induction Sciences and Machine Learning (IBISML), IEICE Technical Report,
Vol. 117, pp. 71–76, 3 2018. in Japanese.
[70] Kei Takeuchi. Jouhou-toukeiryou no bunpu to model no tekisetsu-sa no kijun. Suurikagaku, Vol. 14, No. 3, pp. 12–18, 1976. in Japanese.
[71] Stan Development Team. Stan Modeling Language Users Guide and Reference Manual,
2.26 edition, 2019. https://mc-stan.org.
[72] Seshadri Tirunillai and Gerard J Tellis. Mining marketing meaning from online chatter: Strategic brand analysis of big data using latent dirichlet allocation. Journal ofMarketing Research, Vol. 51, No. 4, pp. 463–479, 2014.
[73] M Antónia Amaral Turkman, Carlos Daniel Paulino, and Peter Müller. Computational
Bayesian Statistics: An Introduction, Vol. 11. Cambridge University Press, 2019.
[74] G. Tusnady, I. Csiszar, and G. Tusnády. Information geometry and alternating minimization procedures. Statistics and Decisions: Supplement Issues, Vol. 1, pp. 205–237,
1984.
[75] Tuomas Virtanen, A Taylan Cemgil, and Simon Godsill. Bayesian extensions to nonnegative matrix factorisation for audio signal modelling. In Acoustics, Speech andSignal Processing, 2008. ICASSP 2008. IEEE International Conference on, pp. 1825–
1828. IEEE, 2008.
[76] Kazuho Watanabe. An alternative view of variational bayes and asymptotic approximations of free energy. Machine learning, Vol. 86, No. 2, pp. 273–293, 2012.
[77] Kazuho Watanabe and Sumio Watanabe. Upper bounds of bayesian generalization
errors in reduced rank regression. IEICE Transactions, Vol. J86-A, No. 3, pp. 278–
287, 2003. in Japanese.
[78] Kazuho Watanabe and Sumio Watanabe. Stochastic complexities of gaussian mixturesin variational bayesian approximation. Journal of Machine Learning Research, Vol. 7,
No. Apr, pp. 625–644, 2006.
[79] Sumio Watanabe. Algebraic analysis for non-regular learning machines. Advances in
Neural Information Processing Systems, Vol. 12, pp. 356–362, 2000. Denver, USA.
[80] Sumio Watanabe. Algebraic geometrical methods for hierarchical learning machines.
Neural Networks, Vol. 13, No. 4, pp. 1049–1060, 2001.
[81] Sumio Watanabe. Algebraix Geometry and Statistical Learning Theory. Cambridge
University Press, 2009.
[82] Sumio Watanabe. Asymptotic equivalence of bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning
Research, Vol. 11, No. Dec, pp. 3571–3594, 2010.
[83] Sumio Watanabe. A widely applicable bayesian information criterion. Journal of
Machine Learning Research, Vol. 14, No. Mar, pp. 867–897, 2013.
[84] Sumio Watanabe. Diﬀerence between bayes cross validation and waic for conditional
independent samples. In Tenth Workshop on Information Theoretic Methods in Scienceand Engineering (WITMSE2017), pp. 38–41, 2017.
[85] Sumio Watanabe. Mathematical theory of Bayesian statistics. CRC Press, 2018.

104

Bibliography[86] Sumio Watanabe. Information criteria and cross validation for bayesian inference inregular and singular cases. Japanese Journal of Statistics and Data Science, pp. 1–19,
2021.
[87] Sumio Watanabe. Waic and wbic for mixture models. Behaviormetrika, pp. 1–17,
2021.
[88] W. Xu, X. Liu, and Y. Gong. Document clustering based on non-negative matrix
factorization. In Proceedings of the 26th annual international ACM SIGIR conference
on Research and development in information retrieval. Association for Computing
Machinery, pp. 267–273, 2003.
[89] Keisuke Yamazaki. Asymptotic Expansion of Stochastic Complexities in Singular
Learning Machines. PhD thesis, Tokyo Institute of Technology, 2003.
[90] Keisuke Yamazaki. Asymptotic accuracy of bayes estimation for latent variables withredundancy. Machine Learning, Vol. 102, No. 1, pp. 1–28, 2016.
[91] Keisuke Yamazaki. Bayesian estimation of multidimensional latent variables and its
asymptotic accuracy. Neural Networks, Vol. 105, pp. 14–25, 2018.
[92] Keisuke Yamazaki and Daisuke Kaji. Comparing two bayes methods based on the freeenergy functions in bernoulli mixtures. Neural Networks, Vol. 44, pp. 36–43, 2013.
[93] Keisuke Yamazaki and Yoichi Motomura. Hidden node detection between observablenodes based on bayesian clustering. Entropy, Vol. 21, No. 1, p. 32, 2019.
[94] Keisuke Yamazaki and Sumio Watanabe. Singularities in mixture models and upperbounds of stochastic complexity. Neural Networks, Vol. 16, No. 7, pp. 1029–1038,
2003.
[95] Keisuke Yamazaki and Sumio Watanabe. Stochastic complexity of bayesian networks.
In Uncertainty in Artificial Intelligence (UAI’03), 2003.
[96] Keisuke Yamazaki and Sumio Watanabe. Algebraic geometry and stochastic complexityof hidden markov models. Neurocomputing, Vol. 69, pp. 62–84, 2005. issue 1-3.
[97] Keisuke Yamazaki and Sumio Watanabe. Singularities in complete bipartite graph-typeboltzmann machines and upper bounds of stochastic complexities. IEEE Transactions
on Neural Networks, Vol. 16, pp. 312–324, 2005. issue 2.
[98] Kenta Yoshida, Tatsu Kuwatani, Takao Hirajima, Hikaru Iwamori, and Shotaro Akaho.
Progressive evolution of whole-rock composition during metamorphism revealed bymultivariate statistical analyses. Journal of Metamorphic Geology, Vol. 36, No. 1, pp.
41–54, 2018.
[99] Elias Zavitsanos, Georgios Paliouras, and George A Vouros. Non-parametric estimationof topic hierarchies from texts with hierarchical dirichlet processes. Journal of Machine
Learning Research, Vol. 12, No. 10, 2011.
[100] Piotr Zwiernik. An asymptotic behaviour of the marginal likelihood for general markov
models. Journal of Machine Learning Research, Vol. 12, No. Nov, pp. 3283–3310,
2011.

105

List of Publications
Peer-reviewed Journal Papers
1. Naoki Hayashi, Sumio Watanabe. "Upper Bound of Bayesian Generalization Error in
Non-Negative Matrix Factorization", Neurocomputing, Volume 266C, 29 November
2017, pp.21-28. doi: 10.1016/j.neucom.2017.04.068.
2. Naoki Hayashi, Sumio Watanabe. "Asymptotic Bayesian Generalization Error in Latent Dirichlet Allocation and Stochastic Matrix Factorization", SN Computer Science,
Volume 1, 69 (2020), pp.1-22. doi: 10.1007/s42979-020-0071-3.
3. Naoki Hayashi.
"Variational Approximation Error in Non-negative Matrix
Factorization", Neural Networks, Volume 126, June 2020, pp.65-75.
doi:
10.1016/j.neunet.2020.03.009.
4. Keita Harada, Naoki Hayashi, Katsushi Kagaya. "Individual behavioral type captured
by a Bayesian model comparison of cap making by sponge crabs", PeerJ 8:e9036,
pp.1-26. doi: 10.7717/peerj.9036.
5. Naoki Hayashi. "The Exact Asymptotic Form of Bayesian Generalization Error in
Latent Dirichlet Allocation". Neural Networks, Volume 137, May 2021, pp.127-137.
doi: 10.1016/j.neunet.2021.01.024.

Peer-reviewed International Conferences
1. Naoki Hayashi, Sumio Watanabe."Tighter upper bound of real log canonical thresholdof non-negative matrix factorization and its application to Bayesian inference," 2017
IEEE Symposium Series on Computational Intelligence (SSCI), Honolulu, HI, 2017,
pp.1-8, doi: 10.1109/SSCI.2017.8280811.
2. Akira Ito, Masaru Okutsu, Masatoshi Yukishima, Ryosuke Matsushita, Naoki Hayashi,
Aiko Furukawa, Gaku Shoji, Takanobu Suzuki. "Trial of damage prediction fortelecommunication conduits using machine learning". 17 World Conference on Earthquake Engineering (17 WCEE), Sendai, Miyagi, 2020, pp.1-10, 6a-0011.

Not-peer-reviewed International Conferences
1. Naoki Hayashi. "Bayesian Generalization Error and Real Log Canonical Thresholdin Non-negative Matrix Factorization and Latent Dirichlet Allocation". Algebraic
Statistics 2020, Online (virtual mini conference). June. 22 - 26, 2020.

106

List of Publications

Not-peer-reviewed Domestic Conferences (in Japanese)
1. 林直輝, 渡邊澄夫. "非負値行列分解の実対数閾値と Bayes 学習への応用", 第 19回情報論的学習理論ワークショップ (IBIS2016), 京都, 信学技報, Vol.116, No.300,
pp.215-220. (2016/11/17 発表).
• (ENG) Naoki Hayashi, Sumio Watanabe."A Real Log Canonical Threshold ofNonnegative Matrix Factorization and Its Application to Bayesian Learning",
IBIS2016, IEICE Technical Report Vol.116, No.300, pp.215-220. (2016/11/17).
2. 林直輝, 渡邊澄夫. "非負値行列分解における実対数閾値の実験的考察", ニュー
ロコンピューティング研究会 (NC), 東京, 信学技報, Vol.116, No.521, pp.85-90.
(2017/3/13 発表).
• (ENG) Naoki Hayashi, Sumio Watanabe. "Experimental Analysis of Real Log
Canonical Threshold in Non-negative Matrix Factorization", Neurocomputing(NC), IEICE Technical Report, Vol.116, No.521, pp.85-90. (2017/3/13).
3. 林直輝, 中村文士. "特異 Bayes 情報量規準による混合正規分布のモデル選択における変分 Bayes 法の実験的考察", 情報論的学習理論と機械学習（IBISML）, 東京, 信
学技報, Vol.117, No.211, pp.19-26. (2017/9/15 発表).
• (ENG) Naoki Hayashi, Fumito Nakamura. "Experimental Analysis of Variational
Bayesian Method in Model Selection of Gaussian Mixture Model by Singular
Bayesian Information Criterion", Information-Based Induction Sciences and Machine Learning (IBISML), IEICE Technical Report, Vol.117, No.211, pp.19-26.
(2017/9/15).
4. 林直輝, 渡邊澄夫. "確率行列分解の実対数閾値と Bayes 学習への応用", 第 20 回
情報論的学習理論ワークショップ (IBIS2017), 東京, 信学技報, Vol.117, No.293,
pp.23-30. (2017/11/9 発表).
• (ENG) Naoki Hayashi, Sumio Watanabe. "Real Log Canonical Threshold of
Stochastic Matrix Factorization and its Application to Bayesian Learning",
IBIS2017, IEICE Technical Report, Vol.117, No.293, pp.23-30. (2017/11/9).
5. 林直輝, 渡邊澄夫. "ハミルトニアンモンテカルロ法を用いた確率行列分解における実対数閾値の実験的考察", ニューロコンピューティング研究会 (NC), 東京. 信学技報, vol. 117, no. 508, NC2017-89, pp. 127-131. (2018/3/14 発表).
• (ENG) Naoki Hayashi, Sumio Watanabe. "Experimental Analysis of Real Log
Canonical Threshold in Stochastic Matrix Factorization using Hamiltonian Monte
Carlo Method", Neurocomputing (NC), Tokyo. IEICE Technical Report, vol. 117,
no. 508, NC2017-89, pp. 127-131. (2018/3/14).
6. 林直輝．"非負値行列分解における変分近似精度の理論解析"，第 21 回情報論的学習
理論ワークショップ (IBIS2018), 札幌. 信学技報, vol. 118, no. 284, IBISML2018-51,
pp. 53-60. (2018/11/5 発表).
• (ENG) Naoki Hayashi. "Variational Approximation Accuracy in Non-negative
Matrix Factorization", IBIS2018, Sapporo. IEICE Technical Report, vol. 118,
no. 284, IBISML2018-51, pp. 53-60. (2018/11/5).
7. 西郷彰, 林直輝, 伊藤孝太朗. "YOLOv3 とドメイン知識を用いた CT 画像の病変部位検出", 第 33 回日本人工知能学会全国大会 (JSAI2019), 新潟. (2019/6/5 発表).
• (ENG) Akira Saigo, Naoki Hayashi, Kotaro Ito. "Lesion Detection in Computed

107
Tomography Images using YOLOv3 and Domain Knowledge", The 33rd Annual
Conference of the Japanese Society for Artificial Intelligence 2019(JSAI2019),
Niigata. (2019/6/5).
8. 林直輝．"LDA における汎化誤差の厳密な漸近形"，第 23 回情報論的学習理論
ワークショップ (IBIS2020), オンライン. (2020/11/26 発表). 動画リンク：https://www.youtube.com/watch?v=e8t3ZGluc6U. 優秀発表賞ファイナリスト.
• (ENG) Naoki Hayashi. "The Exact Asymptotic Form of Generalization Error
in LDA", IBIS2020, Online. (2020/11/26). Outstanding Presentation Award
Finalist.

Domestic Invited Talks (in Japanese)
1. 林直輝．"パラメータ制約付き行列分解のベイズ汎化誤差解析"，第 5 回統計・機械
学習若手シンポジウム (StatsML Symposium’20), オンライン. (2020/12/5 発表). 動
画リンク：https://www.youtube.com/watch?v=ZbVfah9pnb4.
• (ENG) Naoki Hayashi. "Theoretical Analysis of Bayesian Generalization Error
in Parameter-restricted Matrix Factorization", StatsML Symposium’20, Online.
(2020/12/5).

Awards
1. Outstanding Presentation Award Finalist at IBIS2020 (2020/11/26 published, 2021/1/8
gotten).