\section{Introduction}

\IEEEPARstart{P}{hotovoltaic} (PV) smart inverter technology introduced in recent years enables solar panels to act as distributed energy resources (DERs) that can provide bi-directional reactive power support to electric power grid operations \cite{Farivar2012, Seuss_2015, Kekatos2015}. This support can be used to regulate local and system-wide voltages in distributed grids, and the IEEE Standard 1547-2018 \cite{IEEE_1547_2018} provides requirements on the use of such support. Voltage regulation is critical for network safety, both at the transmission and distribution levels.

\inserted{In the distribution grid, voltage regulation is usually controlled either through discrete switching (e.g. tap transformers, capacitor banks) or devices with continuous set points (e.g. PV inverters)}. Broadly speaking there are two categories of control and information structure to address the voltage regulation in distribution grids. The first category of solutions assumes complete or partial knowledge of system parameters and topology (such as a sample of references \cite{Su_2014, Zhang2015, Zhu2016, Lin2018, LiNa_Optimal, LiNa_limited_comm, ElHelou2020}). The second category of solutions primarily relies on observation data and do not explicitly use the knowledge of a physical model of the distribution network (such as a subset of references \cite{ADG, Xu2012, TwoTimescaleDeep, Brandon, SafeRL, ChangfuLi, WeiWang_2019}).

In the first set of solutions, control schemes are adopted based on assumed system models and on knowledge of line parameters. Given the broad set of references, we provide representative ones to illustrate some key directions of research. In \cite{Su_2014}, a multi-objective OPF problem is solved over a given unbalanced distribution grid model using Sequential Quadratic Programming (SQP), where the set of objectives include minimizing line losses, voltage deviation from nominal, voltage phase unbalance and power generation and curtailment costs. In \cite{Zhang2015}, a model-based voltage regulation problem is shown to be solvable by an equivalent Semi-Definite Programming (SDP) problem, and the sufficient conditions under which it can be solved as a result of convexifying the problem are outlined. In \cite{Zhu2016, Lin2018, LiNa_Optimal}, the distribution grid is modelled using the widely adopted linearized flow model, known as LinDistFlow, that assumes a tree (radial) grid structure and negligible line losses. In these papers, the same voltage regulation problem is solved, and an extension to these works are \cite{LiNa_limited_comm, ElHelou2020} in which limited or no communication between buses is needed and the same LinDistFlow model is adopted to provide theoretical guarantees on convergence and stability of the proposed control schemes.

\inserted{In the second set of solutions, reinforcement learning (RL) approaches are used to bypass the need to \textit{explicitly} model the system. {Nonetheless,} some form of a system model and a power flow solver is still needed to simulate the effects of actions on states, as it may not be feasible to learn  by directly interacting with a live power grid. However, knowledge of such a model or its parameters need not be \textit{explicitly} known by the solver of the RL problem, since optimal policies can be inferred solely by interacting with the simulation environment and receiving reward signals. This flexibility to system models is in part what makes RL approaches attractive, in contrast with the aforementioned approaches that rely heavily on a very specific class of system models. Furthermore, the set of objective functions commonly used in conventional control methods tend to be restrictive (e.g. quadratic), whereas in RL-based approaches, reward functions can be arbitrary and are usually designed to more directly reflect the user's underlying objective.}

In \cite{ADG}, for example, Batch RL is adopted to solve the optimal setting of voltage regulation transformers, where a virtual transitions generator is used to allow the RL agent to collect close-to-real samples, for learning, without jeopardizing real-time operation.
In \cite{Xu2012}, the Optimal Reactive Power Dispatch (ORPD) problem is solved using tabular Q-learning, where the objective is to minimize line losses, and actions include various forms of discrete reactive power re-dispatch, including switching control of tap transformers and of capacitor banks. Tabular Q-learning works well when there is a relatively small number of discrete states and actions, but suffers from the \textit{curse of dimensionality} and does not scale up well for larger systems.
In \cite{TwoTimescaleDeep}, \textit{deep} RL is used to optimize reactive power support over two timescales: one for discrete capacitor configuration and the other for continuous inverter setpoints. Here \textit{deep} refers to the representation of Q-value functions (used in Q-learning) in deep neural network form, as opposed to tabular form. This approach is well known as DQN \inserted{(Deep Q-learning)}. This is also applied in \cite{Brandon} to control voltages on the transmission level. While DQN uses neural networks to approximate value functions over continuous state spaces, it still relies on the assumption that the action space is discrete.
Policy gradients, which we review later in this paper, are an alternative set of RL approaches that enable continuous action spaces, and in \cite{SafeRL, ChangfuLi, WeiWang_2019}, policy gradients are used to solve Volt-VAR control problems where reactive power support is selected from a continuous set using a deep neural network that maps states directly to actions.

Such methods are inherently limited by physical constraints on reactive power support which are conventionally assumed to be uncontrollable. Here, we discuss the flexibility of such constraints and the value of relaxing them. \inserted{For example, conventional practices of maximum power-point tracking (MPPT) have been state-of-the-art, wherein each PV inverter is designed to extract the maximum real/active power from the solar panel. However, with a growing number of PV panels in the distribution grid, it becomes important to fully investigate the benefits and costs of always absorbing the maximum real power from the sun into the grid in real-time. By absorbing less real power, for instance, there is more room for reactive power support. In practice, real power curtailment at generators is performed either by solving model-based OPF problems, such as in \cite{Su_2014}, or by using a so-called VoltWatt (VW) droop controller, which requires carefully tuning the droop curves at each generator.}

In this paper, an RL approach is proposed for directly learning decentralized joint active and reactive power support. We illustrate a set of scenarios where instead of injecting all of the solar power into the network, it might rather be better to save or store the power and to inject it at a later time. Even in the absence of a storage system, under deep enough photovoltaic penetration, we
\inserted{discover that the RL agent surprisingly learns by itself that it might} be better to draw only parts of the available power, in order to avoid over-voltage, \inserted{especially} if there is an insufficient amount of reactive power resources available. This is accomplished by designing a reward function, used in the RL approach, that strikes a balance between voltage regulation and real power absorption.

\inserted{Despite modern advancements in deep RL, and in artificial intelligence and machine learning more broadly, the adoption rate of such tools in the power sector is still substantially low compared to other sectors. This is primarily due to the difficulty in interpreting deep neural networks, and due to the complexity of setting up processes that are required for the use of such tools in practice. These issues make it quite risky for a power grid operator to both trust those tools and to prefer them over conventional ones for operating over critical infrastructure. It is thus necessary for the adoption of RL by grid operators that the proposed approaches are heavily simplified and made much more interpretable, which is what we aim to contribute in this study.}

The key contributions of this paper are suggested as follows:

\begin{enumerate}

    \item \inserted{Joint optimization of real and reactive power injection from the PV generation is formulated, with a parameterized reward function tailored for a multi-agent RL approach. Such parameterization is shown to facilitate more interpretable and easier to train deep RL policies, for a more user-friendly experience.}

    \item \inserted{A variant of PPO (Proximal Policy Optimization), a popular RL approach for handling continuous action spaces, is proposed for \textit{decentralized} settings, to dramatically simplify the search for optimal policies for relatively large systems. Not only can it be shown to train as well as or better than a centralized policy architecture, but it is much easier to implement than RL-based decentralized approaches in the existing body of literature.}

\end{enumerate}

\inserted{In our framework, an RL agent per bus observes voltages locally and incrementally updates real and reactive power setpoints at photovoltaic inverters,} similar to an integral droop controller (e.g. \cite{Katiraei2006}). However, it does not rely on any explicit knowledge of network topology or line parameters, and is fully decentralized, requiring minimal communication infrastructures for practical implementation. \inserted{This is illustrated in Fig. \ref{fig:agents_in_grid}, wherein agents are shown to be able to communicate with other agents or a central controller. The word \textit{fully} in the title of this paper refers to the fact that under our proposed framework, each agent in real-time (online) executes its policy locally without communicating with any other agent, not even any of its neighbors. During training though (while exploring and updating policies), observations are aggregated centrally so that all agents can adapt to one another. Since this is done on a much longer-term basis, the communication technologies used do not need to be as sophisticated for this purpose as they might be for an alternative online communication-based approach.

\begin{figure}[!ht]
\centering
\includegraphics[width=0.45\textwidth]{agent_in_grid.pdf}
\caption{\inserted{Every agent in the distribution grid controls a local photovoltaic system, with possible communication between agents. Note: each colored dot in the ``Distribution Grid" box refers to one agent.}}
\label{fig:agents_in_grid}
\end{figure}

Finally, the term ``load" in Fig. \ref{fig:agents_in_grid} is not restrictive. For example, it can include individual residential units or larger and more aggregate neighborhoods. In the former case, examples include houses with rooftop solar and the agent is controlling a few panels, and in the latter case, one bus could correspond to a solar farm and the agent would control many units simultaneously.}

The remainder of the paper is organized as follows. In Section \ref{Preliminaries}, the voltage regulation objective with joint real and reactive power compensation is formulated. In Section \ref{RL_section}, we provide a general review of Markov Decision Processes, and one specific to our problem in Section \ref{Voltage_Regulation_MDP}, with modifications to simplify the task. In Section \ref{policy_optimization}, centralized and decentralized policy architectures are proposed, and they are evaluated in Section \ref{numerical_simulation} with numerical simulations.

\section{Preliminaries}
\label{Preliminaries}

We consider a  three-phase balanced distribution network which consists of a set  $\mathcal{N} =\{0,1,\hdots,N\}$ of buses and a set $\mathcal{L} \subset\mathcal{N}\times\mathcal{N}$ of distribution lines connecting the buses. Bus $0$ represents the substation that acts as the single point  of connection to a bulk power grid. For each distribution line $(i, j ) \in  \mathcal{L}$,   $y_{i,j}$ denotes its admittance.

The set of algebraic power flow equations that govern this three-phase balanced (single-phase equivalent) network are:
\begin{align}
    \label{eq:PF_balanced}
    P_i - \mathbf{i}Q_i = \widetilde{V}^*_i\sum_{j\in\mathcal{N}}y_{ij}\widetilde{V}_j,
    \quad \forall i\in\mathcal{N},
\end{align}
where $P_i$ and $Q_i$ are the net injection of real and reactive power,  $\widetilde{V}_i$ is the complex phasor voltage at bus $i$. $\widetilde{V}^{*}_i$  denotes the complex conjugate of $\widetilde{V}_i$ and $\mathbf{i}:=\sqrt{-1}$.

To model a distribution grid that is not three-phase balanced, or \textit{unbalanced} for short, you may simply replace bus indices with phase indices in Eq. (\ref{eq:PF_balanced}), and $\mathcal{N}$ with the set of all phases. With this, you can generalize over two-phase and single-phase buses, which are common in real distribution grids.

Let  $V_i$ denotes  the \textit{positive sequence} voltage magnitude at bus $i$.  At the substation bus, $V_0$ is fixed at 1.0 p.u. as it is modeled as an ideal voltage source.

Let $\mathcal{C}\subset\mathcal{N}$ be the set of buses that are equipped with solar panels and controllable smart inverters, and $n:=|\mathcal{C}|$. Let $P^c_i$ and $Q^c_i$ be the total real and reactive power, respectively, injected by the PV inverter at bus $i \in \mathcal{C}$. Each inverter has an apparent power capacity $S_i$, which limits $P^c_i$ and $Q^c_i$ as
\begin{align}
    \label{eq:inverter_constraints}
    (P^c_i)^2 + (Q^c_i)^2 \leq S_i^2, \hspace{0.4cm}  P^c_i \leq p^{\text{env}} \leq 0.9S_i,
\end{align}
where $p^{\text{env}}_i$ is the maximum amount of real power that can be drawn from the solar panel at a given moment in time. It depends on exogenous \textit{environmental} factors (irradiance, temperature, etc.), hence the superscript. The upper bound on this quantity is $0.9S_i$ since each inverter in the network is assumed to obey standard IEEE 1547-2018 \cite{IEEE_1547_2018}. We let this injected power be evenly distributed across all phases per bus.

Strictly speaking, if $P^c_i(t)$ is the actual real power injected by the inverter at time $t$, and $\overline{P^c_i}(t)$ is the setpoint, then those two cannot be equal at the same time. There is a small time delay ($\sim$ 10 ms, or less than one 60 Hz cycle) between when the setpoint is assigned and when the actual quantity tracks it. We let both the discrete time step and the tracking time be 10 ms.  This allows us to treat the system as a quasi-steady state system, as illustrated in Fig. \ref{fig:quasi}.

\begin{figure}[!ht]
\centering
\includegraphics[width=0.45\textwidth]{quasi_clear.png}
\caption{Illustration of quasi-steady state behavior.}
\label{fig:quasi}
\end{figure}

More precisely, net real and reactive power injections at bus $i$ at time $t$ can be expressed as $P_i(t) = P^c_i(t) - P^l_i(t)$ and $Q_i(t) = Q^c_i(t) - Q^l_i(t)$,  where $P^l_i(t)$ and  $Q^l_i(t)$ are the (uncontrollable) load consumption   at bus $i$ at time $t$. Due to the one time step delay between setpoint and actual injections, $P^c(t+1)=\overline{P^c}(t)$ and $Q^c(t+1)=\overline{Q^c}(t)$. Then, due to the quasi-steady state nature, $V(t+1) = f(\overline{P^c}(t), \overline{Q^c}(t))$, where $f(\cdot)$ represents the solution to the algebraic power flow equation given in \eqref{eq:PF_balanced}.

We address the problem of voltage regulation in the distribution grid through joint real and reactive power control of PV inverter setpoints. The control objective is to track desired voltage levels while not wasting solar power in the process. The voltage regulation problem  can be formulated as follows:
\begin{subequations}
    \begin{align}
        \begin{split}
            \label{eq:objective}
            \underset{\overline{P^c}, \overline{Q^c}}{\text{maximize}} &\quad \mathbb{E}\left[\sum_{t=0}^{T}\sum_{i\in\mathcal{C}}R_{V_i(t)} + \mu_i R_{P^c_i(t)}\right],
        \end{split}\\
        \begin{split}
            \label{eq:net_P}
            \text{s.t. } &\quad P_i(t) = P^c_i(t) - P^l_i(t),
        \end{split}\\
        \begin{split}
            \label{eq:net_Q}
            &\quad Q_i(t) = Q^c_i(t) - Q^l_i(t),
        \end{split}\\
        \begin{split}
            &\quad \text{Eq. (\ref{eq:PF_balanced}), i.e. \textit{power flow}, } \forall t,
        \end{split}\\
        \begin{split}
            &\quad \text{Eq. (\ref{eq:inverter_constraints}), i.e. \textit{inverter constraints}, } \forall t,
        \end{split}\\
        \begin{split}
            \label{eq:R_V}
            & R_{V_i(t)} = \cfrac{1}{0.05}\min\left\{\delta - \left|1-V_i(t)\right|,0\right\},
        \end{split}\\
        \begin{split}
            \label{eq:R_P}
            & R_{P^c_i(t)} = \cfrac{P^c_i(t)}{0.9S_i},
        \end{split}
    \end{align}
\end{subequations}
where $\mu_{i}$'s and $\delta$ are positive scalars. Unlike with $P^c$ and $Q^c$ (controllable inverter), $P^l$ and $Q^l$ (uncontrollable load) are generally not evenly distributed across all individual phases. For clarity, subscript $i$ refers to $i^\text{th}$ bus index, and terms where the subscript is dropped refer to all buses collectively.

\begin{figure}[!ht]
\centering
\includegraphics[width=0.49\textwidth]{R_V.png}
\caption{Voltage deviation reward function $R_V$ for different $\delta$ (see Eq. (\ref{eq:R_V})).}
\label{fig:R_V}
\end{figure}

Voltage deviation (from nominal 1.0 p.u.) at each bus is considered acceptable if it is kept within some user-defined $\delta$. Deviations greater than this are assigned negative rewards, as depicted in Fig. \ref{fig:R_V}, to signify an undesirable voltage profile. The reward term $R_V$ in \eqref{eq:R_V} quantifies this criteria.

In conventional voltage regulation schemes, reactive power injection or consumption is used to alleviate voltage deviation problems. Reactive power support is limited by physical constraints. For example in the case of PV inverters, as expressed in Eq. (\ref{eq:inverter_constraints}), real power injection by the solar panel directly limits available reactive power. In our proposed objective, we not only consider reactive power but real power provision as well as a decision variable to add a degree of freedom and to relax constraints on reactive power support if needed. Since we seek to extract as much real power as possible from the solar panel, physically bounded by $0.9S_i$ as expressed in Eq. (\ref{eq:inverter_constraints}), we assign a positive reward  to more power drawn from the solar panels, as quantified in \eqref{eq:R_P} and \eqref{eq:objective}.

Under this framework with joint provision of real and reactive power, the user (e.g. utility) selects parameter $\mu$. This parameter acts as a balancing term, between voltage deviation minimization and solar production maximization, considering the fact that over-injection of power leads directly both to over-voltage and to tighter constraints on reactive power. A value for $\mu$ chosen too high yields a control scheme equivalent to conventional maximum power point tracking (MPPT) control with limited reactive power support. On the other hand, as $\mu$ approaches zero, much less real power is likely to be drawn from solar panels by the optimal controller. In section \ref{numerical_simulation}, we demonstrate a balanced choice of $\mu$.

The following assumptions are made about variables that are not explicitly controlled:
\begin{itemize}
    \item While provided to the simulator, neither network topology nor line parameters are used by the controller at any time during training or execution. That is, no a priori knowledge of such values is needed.
    \item No load or solar forecasting is made available to the controller, neither upon training nor during execution.

    \item Net load at a control bus is measured by the controller before supplying a setpoint to the solar panel inverter.
\end{itemize}
The voltage regulation problem formulated in this section is posed as a Markov Decision Process (MDP) in Section \ref{Voltage_Regulation_MDP}, but first, in Section \ref{RL_section}, MDP terminology and formalism is introduced for a more general class of control problems.

\section{Markov Decision Process and Reinforcement Learning}
\label{RL_section}
In this section, we give a brief review of the Markov Decision Processes (MDP) terminology. We will also review a specific RL algorithm, called Proximal Policy Optimization (PPO) \cite{PPO}, which we use to solve the voltage regulation problem.

MDP is a standard  formalism for modeling and solving  control problems. The goal is to solve sequential decision making (control) problems  where the control actions can influence the evolution of the state of the system. An MDP can be defined as a four-tuple $(\mathcal{S},\mathcal{A},\Bar{P},R)$, where $\mathcal{S}$ is the state space and $\mathcal{A}$ is the action space. $\Bar{P}(s'|s,a)$ is the probability of transitioning from state $s$ to $s'$ upon taking action $a$, and $R(s)$ is the reward collected at this transition. We consider a finite horizon MDP setting with horizon (episode) length of $T$. We assume that the rewards depend only on the state, not on the actions. A control policy $\pi: \mathcal{S} \rightarrow \mathcal{A}$ specifies the control action to take in each possible state. The performance of a policy is measured using the metric of value of a policy, $\Bar{V}_{\pi}$, defined as,
\begin{align}
    \label{eq:value_fcn}
    \Bar{V}_{\pi}(s) = \mathbb{E}_{\pi} \left[\sum^{T-1}_{t=0} R(s(t)) | s(0) = s \right],
\end{align}
where, $s(t+1) \sim \Bar{P}(\cdot|s(t),a(t)), a(t) = \pi(s(t)),$ and $s(t)$ is the state of the system at time $t$ and $a(t)$  is the action taken at time $t$. The goal is to find the optimal policy $\pi^{*}$ that achieves the maximum value, i.e, $\pi^{*} = \arg \max_{\pi} \Bar{V}_{\pi}$. The corresponding value function, $\bar{V}^{*} = \bar{V}_{\pi^{*}}$, is called the optimal value function. $\pi^{*}$ and $\bar{V}^{*}$ satisfy the Bellman equation,
\begin{align}
    \pi^*(s) = \underset{a\in\mathcal{A}}{\text{argmax}}\left[R(s) + \sum_{s'\in\mathcal{S}} \bar{P}(s'|s,a) \bar{V}^{*}(s')\right]
\end{align}

When the system model $\bar{P}$ is known, the optimal policy and value function can be computed using dynamic programming \cite{sutton1998}. However, in most real-world applications, the  system model is either unknown or extremely difficult to model. In the voltage regulation problem, the line parameters and/or the topology of the network may not be known a priori. Even if they were known, it would still be difficult to model the effect of actions on states in a feed-forward fashion, due to the algebraic nature of their relationship.  In such scenarios, the optimal policy has to be learned from sequential state/reward observations by interacting with an environment, which in this paper is a simulation environment as described in Section \ref{Voltage_Regulation_MDP}. Reinforcement learning is the approach for computing the optimal policy for an MDP when the model is unknown.

Policy gradient algorithms are  a popular class of RL algorithms. In a policy gradient algorithm, we represent the policy as $\pi_{\theta}$, where $\theta$ denotes parameters of the neural network used to represent the policy. Let $J(\theta) = \mathbb{E}_{s}[\bar{V}_{\pi_{\theta}}(s)],$ where the expectation is w.r.t. to a given initial state distribution.  The goal is to find the optimal parameter $\theta^{*} = \arg \max_{\theta} J(\theta)$. This is achieved by implementing a gradient descent update, $\theta_{k+1} = \theta_{k} + \alpha_{k} \nabla J(\theta_{k})$, where $\alpha_{k}$ is the learning rate. The gradient $\nabla J(\theta)$ is given by the celebrated policy gradient theorem \cite{sutton1998} as $\nabla J(\theta) = \mathbb{E}_{\pi_{\theta}}[\Bar{Q}_{\pi_{\theta}}(s,a) \nabla \log \pi_{\theta}(s,a)]$, where expectation is w.r.t. the state and action distribution realized by following the policy $\pi_{\theta}$. Here $\Bar{Q}_{\pi_{\theta}}$ is the Q-value function corresponding to the policy $\pi_{\theta}$. Often, the Q-value function is represented using a neural network of its own (different from the one used for policy representation). The neural network which represents the policy is called the actor network and that which represents the value function is called the critic network. The terminology is due to the fact that the policy network determines actions given observations and the value network provides the `critic' feedback to update the policy parameter, as clear from the expression for $\nabla J(\theta)$. This class of algorithms is also called actor-critic algorithms. The goal is to incrementally update the parameters of both networks in such a way that they converge to yield parameters corresponding to optimal policy and Q-value functions.

Trust region policy optimization (TRPO) \cite{TRPO} is a recent variant of policy gradient algorithms. For improving the sample efficiency and ensuring reliable convergence, TRPO modifies the policy update as
\begin{align}
    \theta_{k+1} &= \underset{\theta}{\text{argmax}}~ \mathbb{E}_{\theta_{k}} \left[\frac{\pi_{\theta}(s,a)}{\pi_{\theta_{k}}(s,a)} \Bar{Q}_{\pi_{\theta_{k}}}(s, a)\right] \\
    &\text{s.t.}~~ \mathbb{E}\left[D_{\text{KL}}(\pi_{\theta_{k}}(s, \cdot), \pi_{\theta}(s, \cdot)) \right] \leq d,
\end{align}
where, $D_{\text{KL}}(\cdot, \cdot)$ is the Kullback-Leibler divergence between two policies, and $\pi_{\theta}(s,a)$ denotes probability of selecting action $a$ given state $s$. Constant $d$ is a user-defined threshold.

Proximal policy optimization (PPO) algorithm \cite{PPO} builds upon the TRPO framework by modifying the objective function and optimization update which enables improved data efficiency and easier implementation since the KL divergence constraint is dropped and the objective is clipped (or saturated) as described in \cite{PPO}. Clipping of the objective is performed to discourage the optimizer from over-updating $\theta$. We adapt this state-of-the-art algorithm to a decentralized setting to solve the voltage regulation problem we consider.

\section{Voltage Regulation  as an RL Problem}
\label{Voltage_Regulation_MDP}
We use the MDP formalism to model the voltage regulation problem presented in Section \ref{Preliminaries}. As a reminder, $n:=|\mathcal{C}|$.

\subsubsection{State space}
The state space $\mathcal{S}\subset\mathbb{R}^{2n}$ is the set of real power injections and voltage measurement at all controllable buses. For convenience, each state $s\in\mathcal{S}$ is defined as an affine transformation of those measurements. More precisely, the state of the system at time $t$, $s(t)$, is given by
\begin{align}
    \label{eq:s_PV}
    s(t) := (&s^P_1(t),~s^V_1(t),~\cdots,~s^P_n(t),~s^V_n(t)),  \\
    \text{where,}~ &s^P_i(t) \gets \cfrac{P^c_i(t)}{0.9S_i}-1,~~ s^V_i(t) \gets \cfrac{1-V_i(t)}{0.05}. \nonumber
\end{align}
So, when $P^c_i(t)$ is its maximum allowable value, $s^P_i(t)$ is zero. Also, when the voltage is equal to the nominal value,  $s^V_i(t)$ is zero. Thus, ideal scenarios correspond to the state value zero, and critical scenarios correspond to magnitudes of order one or less, assuming critical voltages exceed $1\pm0.05$. This scaling helps to initialize and train the RL algorithm.

\subsubsection{Action space}
Given voltage measurements at every bus in $\mathcal{C}$, we can choose two approaches to determine $\overline{P^c}$ and $\overline{Q^c}$: 1) Directly determine optimal real and reactive power setpoints, i.e. $\left(\overline{P^c},\overline{Q^c}\right)$, by algebraically tying to voltage, or 2) change setpoints incrementally, i.e. $\left(\Delta\overline{P^c},\Delta\overline{Q^c}\right)$, similar to an integral controller. The first approach requires the design and memorization of a highly non-linear function that is likely dependant on system operating conditions. Due to the lack of tracking in this approach, forecasting would be required to respond to different operating conditions. The second approach, on the other hand, enables tracking a desired state in a simple and incremental way \cite{ElHelou2020}. We use the second approach in this paper.

The action space $\mathcal{A} \in [-1,1]^{2n}$ is the set of possible \textit{scaled increments} in real and reactive power setpoints. The action at time $t$, $a(t)$, is given by $a(t) = ( a^P_1,~a^Q_1, \cdots,~a^P_n,~a^Q_n)$  and the increment to those setpoints are defined respectively as $\Delta\overline{P^c_i} = a^P_i\cdot\Delta_{\max}^P$ and $\Delta\overline{Q^c_i} = a^Q_i\cdot\Delta_{\max}^Q$. Here $\Delta_{\max}$ explicitly limit the size of actual (as opposed to scaled) increments $\left(\Delta\overline{P^c},\Delta\overline{Q^c}\right)$.

\subsubsection{Transition Model}
We assume that \textit{next states} are obtained by interaction either with a real-world distribution grid or with a simulator, such as OpenDSS \cite{OpenDSS}. In the case of a simulator, provided changes in loads $(P^l,~Q^l)$ and current state and action based on $\mathcal{S}$ and $\mathcal{A}$, next states can be computed directly. For example, actions are mapped to states using OpenDSS as follows:
\begin{align}
    V(t+1) = \text{OpenDSS}(\overline{P^c}(t),\overline{Q^c}(t))
\end{align}
We choose to use OpenDSS simulator for two main reasons: (a) It can solve power flow for unbalanced distribution grids, and (b) one can directly interact with it using Python, where RL methods are easier to implement.

\subsubsection{Reward function}
System-wide reward at every time step is obtained as follows:
\begin{align}
    R_t = \cfrac{1}{n}\sum_{i\in\mathcal{C}}R_{V_i(t)} + \mu_i R_{P^c_i(t)}
\end{align}
where $R_{V_i(t)}$ and $R_{P^c_i(t)}$ are defined in Eq. (\ref{eq:R_P},\ref{eq:R_V}).

Note that the state and action spaces have been defined in such a way that each element ranges from $-1$ to $1$, with an exception where the voltage-related state may exceed $\pm 1$ if the p.u. voltage exceeds $1\pm 0.05$ under abnormal conditions. This is a suitable choice for training an RL policy as it allows for initialization and adjustment of policy parameters $\theta$ in a standard way by exploiting the  state of the art algorithms (most of  which requires that state and action spaces be a \textit{box} inside $\pm 1$ along all dimensions).

\label{nested_RL_integral}
Based on the definition of action space, the RL agent seeks to learn the \textit{magnitude and direction} in which to incrementally change the setpoints, for every starting state. This raises the question: what information does the agent need to guide this action? The state defined in Eq. (\ref{eq:s_PV}) has the following advantage: If both the voltage term and the power term are zero (i.e. maximum power drawn and nominal voltage), then the scenario is ideal and no \textit{extra} injection is needed. If the load changes in the system, though, a simple amendment to the RL controller is needed:
\begin{equation}
    \label{eq:nested}
    \begin{aligned}
        \Delta\overline{P^c} &\gets \Delta_{max}^P\cdot a_P +\Delta P^l \\
        \Delta\overline{Q^c} &\gets \Delta_{max}^Q\cdot a_Q +\Delta Q^l \\
    \end{aligned}
\end{equation}
where $a_P$ and $a_Q$ (both in $[-1,+1]$) are determined by the RL agent's zero-centered policy $\pi$, and $\left(\Delta P^l, \Delta Q^l\right)$ is the observed change in load at the controllable buses.

The strategy adopted in Eq. (\ref{eq:nested}) is termed an integral controller since setpoint $(\overline{P^c},\overline{Q^c})$ behaves as a discrete-time integrator of changes in operating conditions. Moreover, this controller tracks the state to zero in steady state, within resource limits, since all terms in Eq. (\ref{eq:nested}) go to zero if $s=0$. Under scarcity of resources, one or more of the state terms in Eq. (\ref{eq:s_PV}) will be non-zero, which calls for a balance between \textit{maximum power point tracking} and \textit{voltage regulation}.

Note that state-tracking incremental setpoint changes are bounded by $\Delta_{max}^P$ and $\Delta_{max}^Q$ to limit fluctuations. These values are chosen heuristically as $0.09S_i$ and $0.2S_i$ respectively since those are one-tenth of the maximum possible jumps in setpoints $\overline{P^c}$ and $\overline{Q^c}$.

\section{Control Policy Architecture and Optimization}
\label{policy_optimization}
In this section, we present the design and architecture of our RL algorithm for voltage regulation. We build on the PPO \cite{PPO}  algorithm and extend it to a decentralized setting. Fig. \ref{fig:framework} summarizes the RL-based control policy architecture we propose. \textit{State Observer} refers to Eq. (\ref{eq:s_PV}) and \textit{Integral Controller} refers to Eq. (\ref{eq:nested}). $\Delta\theta_\pi$ refers to changes in policy $\pi$ determined by the PPO algorithm. The policy and critic network architectures used to implement PPO are detailed in the remainder of this section.
\begin{figure}[!ht]
\centering
\includegraphics[width=0.45\textwidth]{Loop.png}
\caption{RL-based Control Architecture.}
\label{fig:framework}
\end{figure}

The standard PPO algorithm assumes that there is a single centralized agent which fully observes the state of the system and can take any control action. However, this centralized control policy may not be feasible for the voltage regulation problem we consider. Firstly, the RL algorithm may not be able to scale to a large network with many nodes. Secondly, even if a centralized training of an RL algorithm is feasible,  the implementation of such a centralized control policy in the real-world system may not be possible due to the communication infrastructure needed. We propose a decentralized RL algorithm that overcomes these challenges.

Recall that $s_i$ and $a_i$ refer local state and action (at bus $i$) respectively. As defined in Eq. (\ref{eq:s_PV}), each state $s_i$ contains two terms per bus, relating to real power and voltage measurements. Similarly, action $a_i$ was defined in such a way that it also contains two terms per bus, relating to changes in real and reactive power setpoints. Our goal is to find the optimal policy parameter $\theta^*_i$ for each bus $i$ that maps local state $s_{i}(t)$ to local control action $a_{i}(t)$, i.e. $a_{i}(t) = \pi_{\theta^{*}_{i}}(s_i(t))$, in an optimal way. The objective is to maximize the cumulative global (system-wide) reward.

To replace the centralized policy, used in standard PPO, with a decentralized policy, we propose a neural network architecture for $\pi$ that connects input to output only at the same bus, rendering it equivalent to a decentralized controller, to compete with conventional methods. That is, there are $n$ neural networks in parallel, each with just 2 inputs and 2 outputs. Each of those smaller networks is parametrized by a group of weights and biases, denoted collectively as $\theta_i$, and notation $\pi_{\theta_i}$ is shortened to $\pi_i$. This time, the PPO algorithm optimizes over $\theta_1,\hdots,\theta_n$, in search of optimal policies $\pi_1,\hdots,\pi_n$, where
\begin{align}
    a_i \gets \pi_i(s_i) \quad \forall i \in \mathcal{C}
\end{align}

Note that the only difference between this case and the centralized case (optimizing over $\theta$), is that here we enforce the strict rule that all neural network weights connecting states at bus $i$ to actions at bus $j$ are fixed at zero iff $i\neq j$. One can also modify this architecture by replacing the condition $i\neq j$ with $(i,j)\notin\mathcal{L}$, if the desired setup involves neighboring buses communicating with one another. For training and implementation purposes, we perform orthogonal initialization on neural network weights for all policies and assign very small initial values to those in the last layer to prevent instability in the feedback controller.

In PPO, value function $\Bar{V}$ (see Eq. (\ref{eq:value_fcn})) is also trained at every iteration when the policy is trained. Even though we desire a decentralized control setting, actor policies are not trained to optimize over local reward maximization, rather they are trained to maximize global (system-wide) reward. For this reason, instead of $n$ value functions ($\Bar{V}_i$ for each $i\in\mathcal{C}$), there is a single value function (critic) for all actors combined, and this function's argument is the system-wide state. In contrast, Multi-Agent Deep Deterministic Policy Gradients (MADDPG) \cite{MADDPG}, a method that also employs multiple policy networks in a policy gradient approach, expects each `agent' not only to observe and act locally (with its own policy network), but to also have its own critic network (value function estimate) and to estimate the policy network of others. In that paper, each agent is assumed to have their own reward function or objective. However, since there's always a global objective (single reward signal), then only one critic function suffices, and the policy networks train simultaneously to adapt to each other.

Due to the way the policy architecture was chosen, there is still a single policy function $\pi$ from a standard PPO algorithm perspective, with $2n$ inputs and $2n$ outputs in this formulation. Since there is a single value function, the PPO agent is trained as usual, without any explicit changes to the algorithm, other than the aforementioned restriction that the policy network connects input to output only at the same bus. To implement this, the optimizer (e.g. Adam optimizer in PyTorch) is told to ignore the weights (initialized and left at zero) in the policy network which link actions at one bus to states at another.

\begin{algorithm*}
\caption{\inserted{One Episode of Interacting With the Distribution Grid}}\label{alg:algo_episode}
\inserted{
\begin{algorithmic}[1]
\renewcommand{\algorithmicloop}{\textbf{function }}
\Loop \texttt{RunEpisode}(\texttt{in\_training}, $\mathcal{D}$) \Comment{\texttt{in\_training} is True or False, $\mathcal{D}$ is experience buffer}
\State Specify load \& solar scenario (random \& time of day).
\If{\texttt{in\_training}}
    \State Set agent policies to sample actions stochastically
\Else
    \State Set agent policies to sample actions deterministically
\EndIf
\While{episode not done}
    \For{each controllable bus $i$ in grid}
        \State Observe local voltage and active power measurements.
        \State Use policy $\pi_i$ to sample changes to active \& reactive power setpoints at local PV panels.
    \EndFor
    \State Collect rewards and append actions, observations, \& rewards to $\mathcal{D}$. \Comment{Do this only if \texttt{in\_training}}
\EndWhile
\EndLoop
\end{algorithmic}
}
\end{algorithm*}

\begin{algorithm*}
\caption{\inserted{Training Decentralized PPO for Distribution Grids}}\label{alg:algo_train}
\inserted{
\begin{algorithmic}[1]
\State Initialize the following:
\begin{itemize}
    \item policy neural network, $\pi_i$ for each bus $i$
    \item value neural network, $\bar{V}$
    \item PPO hyper-parameters
    \item experience buffer, $\mathcal{D}$ \Comment{memory of actions, observations and rewards used in PPO}
\end{itemize}
\State \texttt{in\_training} $\gets$ True
\While{number of episodes less than threshold}
    \For{each episode}
        \State \texttt{RunEpisode}(\texttt{in\_training}, $\mathcal{D}$) \Comment{Defined in Algorithm \ref{alg:algo_episode}}
    \EndFor
    \renewcommand{\algorithmicwhile}{\textbf{every}}
    \While{pre-specified number of episodes}
        \State Use $\mathcal{D}$ to update policy and value networks according to PPO.
    \EndWhile
\EndWhile
\end{algorithmic}
}
\end{algorithm*}

\section{Numerical Simulation}
\label{numerical_simulation}
In this section, we apply the proposed policy architecture and use PPO to solve the MDP. Numerical simulations are conducted on a 240-node distribution grid (see Fig. \ref{fig:grid}) using OpenDSS to solve unbalanced power flow, as described in Section \ref{Voltage_Regulation_MDP}. All parameters associated with this network are obtained from real line parameters and real load data, based on an anonymous distribution grid in Midwest U.S. \cite{IowaState}. Experiment details (e.g. software and hardware details) are found in the Appendix \ref{Appendix_Experiment_Details}.

Based on numerical simulations, as shown in Fig. \ref{fig:rewards}, we have found that the decentralized agent is more sample efficient and trains with fewer fluctuations and variance in episodic rewards over the learning process. On the other hand, the centralized agent takes a bit less computation time (about 20\% less) per iteration, yet takes more iterations to converge.

\subsection{Simulation Setup}
The RL agent interacts with the distribution grid every 10 ms (the time step), and each episode contains 100 time steps, for a total of one second per episode. $\mu$ is set to $0.1$ to favor voltage regulation over solar production maximization. We use the distribution grid shown in Fig. \ref{fig:grid}, where $N=240$, and we select $n=16$ and $n=194$ for the case studies that follow. 194 is the number of controllable nodes provided originally with the OpenDSS model of this grid. For each of these 194 nodes, we have one year \inserted{(2017)} of real historical load data $(P^l,~Q^l)$, which we take advantage of to generate random samples for our simulation at the start of every episode. \inserted{A short summary of the load data is shown in Fig. \ref{fig:load_stat}.}

\begin{figure}[!ht]
\centering
\includegraphics[width=0.45\textwidth]{load_stat.png}
\caption{\inserted{Example of real load data from grid located in the Midwest U.S. \cite{IowaState}. (Top) about 200 load profiles in the first week of January 2017. (Bottom) load distribution over the entire year of 2017.}}
\label{fig:load_stat}
\end{figure}

Since each episode is 1 second, it is fair to assume that fluctuations in $p^{env}$ and $(P^l,~Q^l)$ are negligible within one episode. For this reason,
\inserted{at the beginning of each episode, as indicated in Algorithm \ref{alg:algo_episode}, we} randomly generate and fix $p^{env}$ and $(P^l,~Q^l)$ for the remainder of the episode. \inserted{Maximum solar power output, $p^{env}$, for training purposes is randomly selected in each episode as a multiple of the net load at the same bus. For example, at a given bus, if the net load is $x$ kW, then $p^{env}$ at the same bus is chosen uniformly between 0 and $2x$ kW. This helps generate a diverse set of scenarios, ranging from net over-consumption to net over-production at each bus.}

\inserted{The implementation of each episode is detailed in Algorithm \ref{alg:algo_episode}. If the episode is being run for the purpose of training, experiences are collected and stored in the experience buffer (also known as the replay buffer). Furthermore, each agent acts deterministically during online interaction with the environment. However, during training, actions are sampled from a Gaussian, whose mean and variance are determined by the policy's neural network. The \texttt{in\_training} boolean (True or False) is used in the algorithm to distinguish between online execution only, or training. Finally, Algorithm \ref{alg:algo_train} describes how interaction with the environment is used to train the policy and value networks ($\pi_i~\forall i$ and $\bar{V}$).}

\begin{figure}[!ht]
\centering
\includegraphics[width=0.45\textwidth]{diagram_small.png}
\caption{Real distribution grid located in the Midwest U.S. \cite{IowaState}, with 240 nodes excluding the substation node.}
\label{fig:grid}
\end{figure}

\subsection{Case Study on a smaller (16-bus) subsystem}
In this case study, we compare the use of a centralized policy to that of a decentralized policy, presented in Section \ref{policy_optimization}.

Since $n=16$, neural networks of both centralized and decentralized policies have 32 inputs and 32 outputs. The standard choice of 2 hidden layers with 64 neurons per layer is made for the centralized policy, with $tanh(\cdot)$ activation functions, whereas the decentralized policy splits into 16 sub-policies each with 2 inputs, 2 outputs, and two hidden layers each with 4 neurons. This gives both the centralized and decentralized policies a `height' of 64 neurons in the hidden layer ($4\times16 = 64$), but a total of 8352 parameters to tune for the former and 672 for the latter. In fact, in the decentralized case, we assign 16 different Adam optimizers, one to tune each sub-policy, so it's not so much 672 parameters to optimize per PPO iteration, rather 42 per optimizer, compared to 8352 per (single) optimizer in the centralized setting.

The training curve for each is shown in Fig. \ref{fig:rewards}, where each `PPO iteration' on the $x$-axis refers to 2048 steps of interacting with the environment (or 20s, considering 10ms time step). It is evident that the centralized agent does not out-perform the decentralized agent, and is clearly less interpretable, and requires a wide communication infrastructure to implement in practice. Note: in both centralized and decentralized cases, \textit{value} function $V$ is centralized (fully connected neural network). That is, the RL agent is centralized during training (computer simulation), but decentralized during execution (real-world).

\begin{figure}[!ht]
\centering
\includegraphics[width=0.45\textwidth]{rewards_mixed.png}
\caption{RL training curve for centralized vs. decentralized policies. Decentralized agent trains more monotonically and with less variance in episodic rewards. For both, rewards below 0 indicate voltages are not within user-defined safety boundary.}
\label{fig:rewards}
\end{figure}

In classic RL benchmarks, a threshold on average cumulative rewards is chosen to determine when the learning problem is \textit{solved}. This helps the person who monitors and debugs the learning process to get a sense of progress to save time and effort. In our context, the threshold is set to $0$, as shown in Fig. \ref{fig:rewards}, for the following reasons. We know that $R_V\leq0$ and $R_P\geq0$, from Eq. (\ref{eq:R_P},\ref{eq:R_V}). Both reward terms have been designed in such a way that the magnitudes of the rewards are of order 1 or less during normal operating conditions. Moreover, say that the user desires to keep voltages within $1\pm\delta$. It is then a fact that $(R_V+\mu R_P)\geq0$ at every bus only if voltages are kept within the desired region at all buses. It logically follows that if the inequality does not hold, then the voltage at least at one bus must be outside the desired region. Thus, we can state that by simply monitoring the training curve, one may claim that not all voltages are inside $1\pm\delta$ if the curve is still below the threshold, and that certainly more time needs to be given before ending the training progress. This necessary condition on voltage serves as a useful tool for users who seek to implement this approach. \textit{Note:} in that figure, the term 'voltage safety' merely refers to voltages being inside $1\pm\delta$. Fig. \ref{fig:rewards} demonstrates that the decentralized agent permanently crosses this threshold after 2 iterations, while the centralized takes 4 iterations to do so.

By these results, we claim that one can obtain results for a decentralized agent that are similar to, or even better than, those for a centralized agent, simply by manipulating the neural network's architecture.

\subsection{Case Study on a larger (194-bus) subsystem}
In the previous subsection, we compared centralized and decentralized policy architectures. In this subsection, we dig deeper to examine our proposed framework from purely a power systems perspective. We ask the following question: what is the impact of joint real and reactive power control (as opposed to just the latter) on system-wide voltage profile in the midst of deep photovoltaic penetration?

Consider $n=194$ buses, and the same grid as before, with controllable real and reactive power inverter setpoints. As shown in Fig. \ref{fig:compare_PV}, when maximum real power is drawn from the solar panels, leaving less reactive power support, deep photovoltaic penetration causes over-voltage.
With joint provision of real and reactive power, the RL agent manages to keep voltages within the user-defined desired region ($1\pm\delta$). Surprisingly, a small reduction in real power injection was needed to achieve this effect. Fig. \ref{fig:Pp_distribution} shows the steady-state distribution of real power consumption per bus, as a ratio to maximum possible injection ($p^{env}$). It is worth noting how well the voltage was improved system-wide, even though most solar panels produced near maximum output (note the 0.85 on the y-axis of both figures \ref{fig:compare_PV} and \ref{fig:Pp_distribution}). This justifies the value in considering joint provision of real and reactive power support.

\begin{figure}[!ht]
\centering
\includegraphics[width=0.42\textwidth]{PV_BusIndex.pdf}
\caption{\inserted{Comparison between control policies under a scenario with deep photovoltaic penetration. Joint provision of real and reactive power yields improved voltage profile with minor power sacrifice.}}
\label{fig:compare_PV}
\end{figure}

\inserted{Indeed, there exist plenty of studies in the body of literature that implement some form of power curtailment (to provide reactive power support) for the purposes of voltage regulation. We propose the novelty here is twofold. Firstly, by using RL to automatically learn to balance between voltage regulation and maximum power utilization, the user does not need to design heuristics to reach similar results or to explicitly rely on any understanding of how the system works. Secondly, as shown in Fig. \ref{fig:compare_PV}, very little power curtailment is required to significantly improve the voltage profile and keep it well within the desired thresholds. This is made possible by the structure of the proposed decentralized PPO policy which allows all the agents to simultaneously train to adapt to one another during exploration, despite the lack of communication between them during real-time operation.}

\begin{figure}[!ht]
\centering
\includegraphics[width=0.42\textwidth]{Pp_distribution.png}
\caption{Histogram of $P^c/p^{env}$ under RL policy, using results of Fig. \ref{fig:compare_PV}. Most buses inject near-maximum real power.}
\label{fig:Pp_distribution}
\end{figure}

\section{Concluding Remarks}
This paper introduces a reinforcement learning-based voltage control strategy with joint provision of real and reactive power
for distribution grids with  deep photovoltaic penetration. The
\inserted{joint real power and voltage support problem is formulated} as a Markov Decision Process with rewards parametrized to balance between voltage deviation minimization and solar production maximization.

\inserted{Compared with conventional multi-agent RL algorithms, we develop a tailor-designed \textit{decentralized PPO (Proximal Policy Optimization)} algorithm that would both work well with large continuous action spaces and simplify the optimization process by updating multiple policies simultaneously. We demonstrate that by reducing the search space, for the simper 16-bus subsystem, from over 8000 parameters for a centralized setting to under 700 parameters for a decentralized setting, we still achieve similar, and in some cases better, average rewards. This size reduction is on the order of 10 or so, but for the 194-bus subsystem, it is on the order of 100.} Numerical simulations on a 240-node distribution grid based on real parameters show that it is not always the best strategy to absorb all the solar power available. This observation implies that it would benefit the distribution grid the most if some fraction of the real-time power produced from the PV panels can be locally absorbed.

This paper also proposes and verifies a fully decentralized (communication-free) approach for this type of control, which can be implemented on existing physical infrastructure, helping alleviate problems related to communication failure or cyber-attacks. In future work, competition between agents is considered, whereby the inverter at each bus seeks to maximize local, not system-wide, rewards. Further research could also investigate the optimal combination of local energy storage together with PV panels for real-time operation.

\appendices
\section{Experiment Details}
\label{Appendix_Experiment_Details}

\inserted{Software:} Simulations are conducted in Python 3.7, interfacing with OpenDSS \cite{OpenDSS} (unbalanced distribution grid simulator) and using PyTorch \cite{PyTorch} (Python-based library to train neural networks) to model, build and train actor and critic neural networks.
\inserted{Hardware:} Lenovo, 64-bit Windows 10, Intel\textregistered Core\texttrademark i7-6700HQ CPU @ 2.60Ghz, 16.0 GB RAM.

\inserted{Proximal policy optimization (PPO) parameters: Number of steps per training update of 512 to 2048 time steps. Considering each episode/scenario involved 100 time steps, this implies on the order of 10 episodes every actor and critic network update. Batch size of 16 time steps with 10 epochs. This means that for each training update, the aforementioned neural networks are updated by performing nonlinear regression 10 times, each using 16 time steps randomly selected from the experience buffer (refer to Algorithm \ref{alg:algo_train}).}

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}


\section{Introduction}
Latent Dirichlet allocation (LDA) \cite{David2003LDA} is one of topic models \cite{Gildea1999topic} which is a ubiquitous statistical model used in many research areas. Text mining \cite{David2003LDA, Griffiths2004LDAGS}, computer vision \cite{Li2005LDA4CV}, marketing research \cite{Tirunillai2014LDA4MR}, and geology \cite{Yoshida2018LDA4Geo} are such examples. LDA had been originally proposed for natural language processing and it can extract essential information from documents by defining the topics of the words.
The topics are formulated as one-hot vectors subject to categorical distributions which depend on each document. The parameters of those categorical distributions express the topic proportion and they are the object of inference. In addition, the words are also formulated as one-hot vectors generated by other categorical distributions whose parameters represent appearance probability of the words in each document. This appearance probability is also estimated.

In the standard inference algorithms such as Gibbs sampling \cite{Griffiths2004LDAGS} and variational Bayesian method \cite{David2003LDA}, LDA requires setting the number of topics (the dimension of the topic one-hot vector) in advance. The optimal number of topics in the ground truth is unknown, thus researchers and practitioners face to typical model selection problem; the chosen number of topics may be larger than the optimal one. In this situation, the estimated parameter cannot be uniquely determined. From the theoretical point of view, LDA is non-identifiable, i.e. a map from a parameter set to a probability density function set is not injective. Besides, LDA has a degenerated Fisher information matrix and its likelihood and posterior distribution cannot be approximated by any normal distribution. Such models are called singular statistical models and LDA is one of them.

If a map from parameters to probability density functions is injection in a statistical model, then the model is called a regular statistical model. In a regular statistical model, its expected generalization error is asymptotically equal to $d/2n + o(1/n)$, where $d$ is the dimension of the parameter and $n$ is the sample size \cite{Akaike1974AIC}. Moreover, its negative log marginal likelihood (a.k.a. free energy) has asymptotic expansion represented by $nS_n + (d/2) \log n + O_p(1)$, where $S_n$ is the empirical entropy \cite{Schwarz1978BIC}.

\begin{figure}[t]
\begin{minipage}{0.5\hsize}
\begin{flushleft}
\includegraphics[width=7.5cm]{./RLCTandDover2_rvd.eps}
\subcaption{}
\label{fig:RLCT}
\end{flushleft}
\end{minipage}
\begin{minipage}{0.5\hsize}
\begin{flushright}
\includegraphics[width=7.5cm]{./learningcurve_rvd.eps}
\subcaption{}
\label{fig:LearningCurve}
\end{flushright}
\end{minipage}
\caption{(a) In this paper, we give the exact value of the learning coefficient of LDA $\lambda$.
The learning coefficient 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 coefficients of LDA when the number of topics $H$ is increased.
If LDA was a regular statistical model, its learning coefficient would be the dotted yellow line drawn by the squares.
The behavior of them are so different.
\\ (b) This figure shows the theoretical learning 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 is the dashed yellow line.
The vertical axis means the expected generalization error $\mathbb{E}[G_n]$ and the horizontal one is the sample size $n$.
This is based on Eq. (\ref{thm-watanabeG}) and the exact value of $\lambda$ which is clarified by our result.}
\end{figure}
On the other hand, in the general case, by using resolution of singularity \cite{Hironaka}, Watanabe had proved that the asymptotic forms of its generalization error $G_n$ and marginal likelihood $Z_n$ are the followings \cite{Watanabe1,Watanabe2,SWatanabeBookE}:
\begin{align}
\label{thm-watanabeG}
\mathbb{E}[G_n] &= \frac{\lambda}{n} - \frac{m-1}{n \log n} + o\left(\frac{1}{n \log n} \right), \\
\label{thm-watanabeF}
-\log Z_n &= nS_n + \lambda \log n - (m-1) \log\log n + O_p(1),
\end{align}
where $\lambda$ is a positive rational number, $m$ is a positive integer, and $\mathbb{E}[\cdot]$ is an expectation operator on the overall datasets. The constant $\lambda$ is called a learning coefficient since it is dominant in the leading term of the above forms which represent learning curves (Fig. \ref{fig:LearningCurve}). The above forms hold not only in the case the model is regular but also in the case that is singular. For example, neural networks are singular statistical models and these generalization errors can be represented by Eq. (\ref{thm-watanabeG}) \cite{Watanabe2}. LDA is also singular; thus, this study is one of researches which clarify the generalization errors of singular statistical models.
In the regular case, $\lambda=d/2$ and $m=1$ hold. However, in the singular case, they are depend on the model.
Let $K(w)$ be the Kullback-Leibler (KL) divergence between the true distribution to the statistical model, where $w$ is a parameter of the model. The function $K(w)$ is non-negative and analytic. The constants $\lambda$ and $m$ are characterized by a set of the zero points of the KL divergence: $K^{-1}(0)$. $K^{-1}(0)$ is an analytic set (a.k.a. algebraic variety). $\lambda$ is called a real log canonical threshold (RLCT) and $m$ is called a multiplicity in algebraic geometry. In LDA, if the number of topics changes, then $K^{-1}(0)$ does (Fig. \ref{fig:RLCT}).
A model selection method, called {\it sBIC} which uses RLCTs of statistical models, has been proposed by Drton and Plummer \cite{Drton}.
Drton and Imai have empirically verified that sBIC is precise to select the optimal and minimal model if the exact values or tight bounds of RLCTs are clarified \cite{Drton, Drton2017forest, Imai2019estimating}.
Other application of RLCTs is a design procedure for exchange probability in exchange Monte Carlo method by Nagata \cite{Nagata2008asymptotic}. To determine $\lambda$ and $m$, we should consider resolution of singularity for those concrete varieties. In general, we should find RLCTs to a family of functions to clarify a learning coefficient of a singular statistical model. There is no standard method to find RLCTs to a given collection of functions; thus, researchers study RLCTs with considering different procedures for each statistical model. In fact, RLCTs of several models has been analyzed in both statistics and machine learning. For instance, the RLCTs had been studied in Gaussian mixture model \cite{Yamazaki1}, Poisson mixture model \cite{SatoK2019PMM}, reduced rank regression \cite{Aoyagi1}, three-layered neural networks \cite{Watanabe2}, naive Bayesian networks \cite{Rusakov2005asymptotic}, Bayesian networks \cite{Yamazaki3}, Boltzmann machines \cite{Yamazaki4,Aoyagi2,Aoyagi3}, Markov models \cite{Zwiernik2011asymptotic}, hidden Markov models \cite{Yamazaki2}, Gaussian latent tree and forest models \cite{Drton2017forest}, and non-negative matrix factorization \cite{nhayashi2,nhayashi5,nhayashi8}. Note that clarifying the exact value of the RLCT in the all case is challenging problem. Whereas we would like to emphasize that this is not to deny the value and novelty of these researches, in deed, they cannot have clarified the exact value except for Aoyagi's result in 2005 \cite{Aoyagi1}.

In this paper, we derive the exact asymptotic form of the Bayesian generalization error by determination of the exact RLCT in LDA.
This article is divided to four parts. First, we introduced background of this research in the above. Second, we describe the framework of Bayesian inference and relationship between its theory and algebraic geometry. Third, we state the Main Theorem. Fourth, we report result of numerical experiment which is carried out in order to verify the behavior of Main Theorem if the sample size is finite. Lastly, we discuss about this theoretical result and we conclude this paper. We prove the Main Theorem in Appendix.

\section{Framework of Bayesian Inference and its Theory}

In this section, we briefly explain the framework of Bayesian inference and its mathematical theory.
Let $\mathcal{D} = (X_1, \ldots, X_n)$ be a sample (a.k.a. dataset: a collection of random variables) of $n$ independent and identically distributed from a data generating distribution (a.k.a. true distribution). The densities of the true distribution and a statistical model is denoted by $q(x)$ and $p(x|w)$, respectively. These domain $\mathcal{X}$ is a subset of a finite-dimensional real Euclidean or discrete space. Let $\varphi(w)$ be a probability density of a prior distribution.
We define a posterior distribution as the following density function on $W$:
\begin{equation}
\label{def-posterior}
\varphi^*(w|\mathcal{D}) = \frac{1}{Z_n}\varphi(w) \prod_{i=1}^n p(X_i|w),
\end{equation}
where $Z_n$ is a normalizing constant to satisfy the condition $\int \varphi^*(w|\mathcal{D})dw=1$:
\begin{equation}
\label{def-marginallikelihood}
Z_n =  \int \mathrm{d}w \varphi(w) \prod_{i=1}^n p(X_i|w).
\end{equation}

This is called a marginal likelihood or a partition function. Its negative log value is called a free energy $F_n = -\log Z_n$.
Note that the marginal likelihood is a probability density function of a dataset.

The free energy appears in a leading term of the difference between the true distribution and the model in the sense of dataset generating process.
An entropy of the true distribution and an empirical one are denoted by
\begin{align}
\label{def-entropy}
S &= -\int \mathrm{d}x q(x) \log q(x), \\
S_n &= -\frac{1}{n} \sum_{i=1}^n \log q(X_i).
\end{align}
The KL divergence between the true distribution to the statistical model is denoted by
\begin{equation}
\label{def-aveerror}
K(w) = \int \mathrm{d}x q(x) \log \frac{q(x)}{p(x|w)}.
\end{equation}
By definition, $X_i \sim q(x)$ for $i=1,\ldots,n$ and $\mathcal{D} \sim \prod_{i=1}^n q(x_i)$ hold; thus let $\mathbb{E}_n[\cdot]$ be an expectation operator on overall dataset defined by
\begin{equation}
\label{def-expectation}
\mathbb{E}_n[\cdot] = \int \prod_{i=1}^n \mathrm{d}x_i q(x_i) [\cdot].
\end{equation}
The symbol $\mathbb{E}_n$ is referred to as $\mathbb{E}$ so long as there is no risk of confusion.
Obviously, $\mathbb{E}_n[S_n]$ is $S$.
Then, we have the following KL divergence
\begin{align}
\int \prod_{i=1}^n \mathrm{d}x_i q(x_i) \log \frac{\prod_{i=1}^n q(x_i)}{Z_n}
&= -\mathbb{E}_n\left[nS_n \right] - \mathbb{E}_n[\log Z_n] \\
&= -nS + \mathbb{E}_n[F_n].
\end{align}
The expected free energy is an only term depend on the model and the prior.
For this reason, the free energy is used as a criterion to select the model.

A predictive distribution is defined by the following density function on $\mathcal{X}$:
\begin{equation}
\label{def-predictive}
p^*(x | \mathcal{D} ) = \int \mathrm{d}w \varphi^*(w|\mathcal{D}) p(x|w).
\end{equation}
This is a probability distribution of a new data. It is also important for statistics and machine learning to evaluate the dissimilarity between the true and the model in the sense of a new data generating process.
A generalization error $G_n$ is defined by a KL divergence between the true distribution and the predictive one:
\begin{equation}
\label{def-gerror}
G_n = \int \mathrm{d}x q(x) \log \frac{q(x)}{p^*(x|\mathcal{D})}.
\end{equation}

Bayesian inference is defined by inferring that the true distribution may be the predictive one.
For an arbitrary finite $n$, by the definition of the marginal likelihood (\ref{def-marginallikelihood})
and the predictive distribution (\ref{def-predictive}), we have
\begin{align}
p^*(X_{n+1}|\mathcal{D}) &= \frac{1}{Z_n}\int \mathrm{d}w \varphi(w) \prod_{i=1}^n p(X_i|w) p(X_{n+1}|w) \\
&= \frac{1}{Z_n} \int \mathrm{d}w \varphi(w) \prod_{i=1}^{n+1} p(X_i|w) \\
&= \frac{Z_{n+1}}{Z_n}.
\end{align}
Thus, the following equation holds:
\begin{align}\label{Eq-GandFrelation}
-\log p^*(X_{n+1}|\mathcal{D}) &= -\log Z_{n+1} - (-\log Z_n) \\
&= F_{n+1} - F_n.
\end{align}
We can treat the expectation by $X_{n+1} \sim q(x_{n+1})$ as $\int \mathrm{d}x q(x)[\cdot]$ since the sample $(\mathcal{D},X_{n+1})$ is independent. Hence, we have
\begin{equation}
\int \mathrm{d}x_{n+1} q(x_{n+1}) [\cdot]
= \int \mathrm{d}x q(x) [\cdot].
\end{equation}
By definition,
\begin{align}
&\quad -\int \mathrm{d}x_{n+1} q(x_{n+1}) \log p^*(x_{n+1}|\mathcal{D}) \\
&= \int \mathrm{d}x_{n+1} q(x_{n+1}) \log \frac{q(x_{n+1})}{p^*(x_{n+1}|\mathcal{D})} - \int \mathrm{d}x_{n+1} q(x_{n+1}) \log  q(x_{n+1}) \\
&=G_n + S
\end{align}
holds; thus, with integrating the both sides in Eq. (\ref{Eq-GandFrelation}) by $\mathrm{d}x_{n+1} q(x_{n+1})$, we obtain
\begin{equation}
G_n + S = \int \mathrm{d}x_{n+1} q(x_{n+1}) F_{n+1} - F_n.
\end{equation}
Considering expectation on overall dataset, we get
\begin{align}
\mathbb{E}_{n}[G_n] + S &= \mathbb{E}_{n+1}[F_{n+1}] - \mathbb{E}_{n}[F_n].
\end{align}
Hence, $G_n$ and $F_n$ are important random variables in Bayesian inference.
In mathematical theory of Bayesian statistics (a.k.a. singular learning theory),
we consider how they asymptotically behave in the general case \cite{SWatanabeBookMath}.
To establish this theory, resolution of singularity in algebraic geometry has been needed.

Now, we briefly explain the relationship between singular learning theory and algebraic geometry.
As technical assumptions, we suppose the parameter set $W \subset \mathbb{R}^d$ is sufficiently wide compact and the prior is positive and bounded on $K^{-1}(0)$: $0<\varphi(w)<\infty$ for any $w \in K^{-1}(0)$. Moreover, $\varphi(w)$ is a $C^{\infty}$-function on with the compact support $W$.
Considering $K(w)$ in Eq. (\ref{def-aveerror}) and its zero points $K^{-1}(0)$,
we use the following analytic form by \cite{Atiyah1970resolution} of the singularities resolution theorem \cite{Hironaka}.
Atiyah has derived this form of the singularities resolution theorem in order to analyze the relationship between a division of distributions (hyperfunctions) and local type zeta functions \cite{Atiyah1970resolution}.
Watanabe has proved that it is useful for constructing singular learning theory \cite{Watanabe1, Watanabe2, SWatanabeBookE}.
\begin{thm}[Singularities Resolution Theorem]\label{thm-singular}
Let $H$ be a non-negative analytic function on $W \subset \mathbb{R}^d$.
Assume that $H^{-1}(0)$ is not an empty set.
Then, there are an open set $W'$, a $d$-dimensional smooth manifold $\mathcal{M}$, and an analytic map $g:\mathcal{M} \rightarrow W'$ such that $g: \mathcal{M}\setminus g^{-1}(H^{-1}(0)) \rightarrow W'\setminus H^{-1}(0)$ is isomorphic and
\begin{gather*}
H(g(u))=u_1^{2k_1} \ldots u_d^{2k_d}, \\
|\det g'(u)| =b(u)|u_1^{h_1} \ldots u_d^{h_d}|
\end{gather*}
hold for each local chart $U \ni u$ of $\mathcal{M}$,
where $k_j$ and $h_j$ are non-negative integer for $j=1,\ldots,d$, $\det g'(u)$ is the Jacobian of $g$ and $b:\mathcal{M} \rightarrow \mathbb{R}$ is strictly positive analytic: $b(u)>0$.
\end{thm}
Let $\mathrm{Re}(z)$ be the real part of a complex number $z$: $\mathrm{Re}(\alpha+\beta\sqrt{-1})=\alpha$, where $\alpha$ and $\beta$ are real numbers and $\sqrt{-1}$ is the imaginary unit. By using Theorem \ref{thm-singular},
the following theorem is proved \cite{Atiyah1970resolution, Bernstein1972, Sato1974zeta}.
\begin{thm}\label{thm-zeta}
Let $F: \mathbb{R}^d \rightarrow \mathbb{R}$ be an analytic function of a variable $w \in \mathbb{R}^d$.
$a: W \rightarrow \mathbb{R}$ is denoted by a $C^{\infty}$-function with compact support $W$. Then
$$\zeta(z) = \int_W |F(w)|^z a(w) \mathrm{d}w$$
is a holomorphic function in $\mathrm{Re}(z)>0$. Moreover, $\zeta(z)$ can be analytically continued to a unique meromorphic
function on the entire complex plane $\mathbb{C}$. Its all poles are negative rational numbers.
\end{thm}
\begin{example}
As a simple example of Theorem \ref{thm-zeta}, put $F(\theta)=\theta^2$ and $a(\theta)=1$, where $\theta \in [0,1]$.
Consider the univariate function of a complex number $z$
\begin{equation*
\zeta(z) = \int_0^1 \theta^{2z} \mathrm{d} \theta.
\end{equation*}
If $\mathrm{Re}(z)$>0, then this integral is immediately calculated as
\begin{equation*}
\int_0^1 \theta^{2z} \mathrm{d} \theta = \frac{1}{2z+1}.
\end{equation*}
Hence,  $\zeta(z)$ is a holomorphic function in $\mathrm{Re}(z)>0$.
However, $z \mapsto 1/(2z+1)$ is a meromorphic function on the entire complex plane and its pole is $z=-1/2$.
Therefore, we can analytically continue $\zeta(z)$ as the meromorphic function whose pole is $z=-1/2$.
\end{example}
Applying Theorem \ref{thm-singular} to the KL divergence in Eq. (\ref{def-aveerror}), we have
\begin{gather}
K(g(u))=u_1^{2k_1} \ldots u_d^{2k_d}, \\
|\det g'(u)| =b(u)|u_1^{h_1} \ldots u_d^{h_d}|.
\end{gather}
Suppose the prior density $\varphi(w)$ has the compact support $W$ and the open set $W'$ satisfies $W \subset W'$.
By using Theorem \ref{thm-zeta}, we can define a zeta function of learning theory.

\begin{defi}[Zeta Function of Learning Theory]
\label{def-learnzeta}
Let $K(w) \geqq 0$ be the KL divergence mentioned in Eq. (\ref{def-aveerror}) and $\varphi(w)\geqq 0$ be a prior density function which satisfies the above assumption. A zeta function of learning theory is defined by the following univariate complex function
$$\zeta(z) = \int_W K(w)^z \varphi(w) \mathrm{d}w.$$
\end{defi}

\begin{defi}[Real Log Canonical Threshold]
Let $\zeta(z)$ be a zeta function of learning theory represented in Definition \ref{def-learnzeta}.
Consider an analytic continuation of $\zeta(z)$ from Theorem \ref{thm-zeta}.
A real log canonical threshold (RLCT) is defined by the negative maximum pole of $\zeta(z)$ and its multiplicity is defined by the order of the maximum pole.
\end{defi}

Here, we describe how to determine the RLCT $\lambda>0$ of the model corresponding to $K(w)$.
We apply Theorem \ref{thm-singular} to the zeta function of learning theory.
For each local coordinate $U$, we have
\begin{align}
\zeta(z) &= \int_U K(g(u))^z \varphi(g(u)) |\det g'(u)| \mathrm{d}u \\
&= \int_U u_1^{2k_1z+h_1} \ldots u_d^{2k_dz+h_d}\varphi(g(u))b(u) \mathrm{d}u.
\end{align}
The functions $\varphi(g(u))$ and $b(u)$ are strictly positive in $U$;
thus, we should consider the maximum pole of
\begin{equation}
\int_U u_1^{2k_1z+h_1} \ldots u_d^{2k_dz+h_d} \mathrm{d}u = \frac{C_1(z)}{2k_1 z + h_1}\ldots\frac{C_d(z)}{2k_d z + h_d},
\end{equation}
where $(C_j(z))_{j=1}^d$ are non-zero functions of $z \in \mathbb{C}$. Hence, we give the RLCT in the local chart $U$ as follows
\begin{equation}
\lambda_U = \min_{j=1}^d \left\{ \frac{h_j + 1}{2k_j} \right\}.
\end{equation}
By considering the duplication of indices, we can also find the multiplicity $m$.
Therefore, we can determine the RLCT as $\lambda = \min_U \lambda_U$.
The RLCT is equal to the learning coefficient because of Eq. (\ref{thm-watanabeG}) and (\ref{thm-watanabeF}).
That is why we need resolution of singularity to clarify the behavior of the Bayesian generalization error $G_n$ and the free energy $F_n$ via determination of the learning coefficient.

\section{Main Theorem}
\label{sec-main}
In this section, we state the Main Theorem: the exact value of the RLCT of LDA.
\begin{defi}[Stochastic Matrix]
A stochastic matrix is defined by a matrix whose columns are in simplices,
i.e. the sum of the entries in each column is equal to one and each entries are non-negative.

An $M \times N$ matrix $C=(c_{ij})_{i=1,j=1}^{M,N}$ is stochastic if and only if $c_{ij} \geqq 0$ and $\sum_{j=1}^N c_{ij} = 1$ hold.
\end{defi}

In the following, the parameter $w=(A,B)$ is a pair of stochastic matrices and the data $x$ is a one-hot vector.
By definition, a set of stochastic matrices is compact. Let $\mathrm{Onehot}(N):=\{w = (w_j) \in \{0,1\}^N \mid \sum_{j=1}^N w_j =1\}$ be a set of $N$-dimensional one-hot vectors and
$\mathrm{\Delta}_N:=\{c=(c_j)_{j=1}^N \mid \sum_{j=1}^N c_j=1\}$ be an $N$-dimensional simplex.
In LDA terminology, the number of documents and the vocabulary size is denoted by $N$ and $M$, respectively.
Let $H_0$ be the optimal (or true) number of topics and $H$ be the chosen one. In this situation, the sample size $n$ is the number of words in all of the given documents.
Note that $x_i$ is the $i$-th entry of the one-hot vector $x$ in this section (it is not the realization of the $i$-th data $X_i$ of the dataset like the previous section).
See also Table \ref{params} (this table is quoted and modified from our previous study \cite{nhayashi7}).

\begin{table}[htb]
  \centering
  \caption{Description of Variables in LDA Terminology}
  \begin{tabular}{|c|c|c|} \hline
    Variable & Description & Index \\
    \hline\hline
    $b_j=(b_{kj}) \in \mathrm{\Delta}_H$ & topic proportion of topic $k$ in document $j$ & for $k=1,\ldots,H$ \\
    $a_k =(a_{ik}) \in \mathrm{\Delta}_M$ & appearance probability of word $i$ in topic $k$ & for $i=1,\ldots,M$ \\
    \hline
    $x=(x_i) \in \mathrm{Onehot}(M)$ & word $i$ is defined by $x_i=1$ & for $i=1,\ldots,M$  \\
    $y=(y_k) \in \mathrm{Onehot}(H)$ & topic $k$ is defined by $y_k=1$ & for $k=1,\ldots,H$  \\
    $z=(z_j) \in \mathrm{Onehot}(N)$ & document $j$ is defined by $z_j=1$ & for $j=1,\ldots,N$  \\
    \hline
    $*_0$ and $*^0$ & optimal or true variable corresponding to $*$ & - \\
  \hline
  \end{tabular}
\label{params}
\end{table}

Let $A=(a_{ik})_{i=1,k=1}^{M,H}$ and $B=(b_{kj})_{k=1,j=1}^{H,N}$ be $M \times H$ and $H \times N$ stochastic matrix, respectively.
Assume that the pair of stochastic matrices $(A_0, B_0)$ is one of optimal parameters, where $A_0=(a^0_{ik})_{i=1,k=1}^{M,H_0}$ and $B_0=(b^0_{kj})_{k=1,j=1}^{H_0,N}$. Suppose $A_0$ and $B_0$ are full rank.

Here, we define the RLCT of LDA in the below. This definition is also quoted and modified from the statement in \cite{nhayashi7}.

\begin{defi}[LDA]
\label{def-LDAmodel}
Assume that $M\geqq 2$, $N \geqq 2$, and $H \geqq H_0 \geqq 1$.
Let $q(x|z)$ and $p(x|z,A,B)$ be conditional probability mass functions of $x \in \mathrm{Onehot}(M)$ given $z \in \mathrm{Onehot}(N)$ as the following:
\begin{align}
\label{topic-true}
q(x|z) & = \prod_{j=1}^N \left(\sum_{k=1}^{H_0} b^0_{kj} \prod_{i=1}^M (a^0_{ik})^{x_{i}}\right)^{z_j}, \\
\label{topic-model}
p(x|z,A,B) & = \prod_{j=1}^N \left(\sum_{k=1}^{H} b_{kj} \prod_{i=1}^M (a_{ik})^{x_{i}}\right)^{z_j}.
\end{align}
The prior density function is denoted by $\varphi(A,B)$.
The conditional mass $q(x|z)$ and $p(x|z,A,B)$ represent the true distribution of LDA and the statistical model of that, respectively.
\end{defi}

These distributions are the marginalized ones of the followings which contain the true topics $y^0 \in \mathrm{Onehot}(H_0)$ and the one of the model $y \in \mathrm{Onehot}(H)$:
\begin{align*}
q(x,y^0|z) & = \prod_{j=1}^N \left[\prod_{k=1}^{H_0} \left(b^0_{kj} \prod_{i=1}^M (a^0_{ik})^{x_{i}}\right)^{y^0_k}\right]^{z_j}, \\
p(x,y|z,A,B) & = \prod_{j=1}^N \left[\prod_{k=1}^{H} \left(b_{kj} \prod_{i=1}^M (a_{ik})^{x_{i}}\right)^{y_k}\right]^{z_j}.
\end{align*}
In fact, $q(x|z) = \sum_{y^0 \in \mathrm{Onehot}(H_0)} q(x, y^0 | z)$ and $p(x|z, A, B) = \sum_{y \in \mathrm{Onehot}(H_0)} p(x,y|z,A,B)$ hold because of $z \in \mathrm{Onehot}(N)$.
In practical cases, the topics are not observed; thus, we use Eqs. (\ref{topic-true}) and (\ref{topic-model}) as the definition of LDA.

\begin{defi}[RLCT of LDA]
\label{def-LDARLCT}
Let $K(A,B)$ be the KL divergence between $q(x|z)$ and $p(x|z,A,B)$:
$$K(A,B)=\sum_{z \in \mathrm{Onehot}(M)} \sum_{x \in \mathrm{Onehot}(N)} q(x|z)q'(z) \log \frac{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 and assumed that it is positive and bounded.
Assume that $\varphi(A,B) >0$ is positive and bounded on $K^{-1}(0) \ni (A_0,B_0)$.
Then, the zeta function of learning theory in LDA is
the holomorphic function of univariate complex variable $z \ (\mathrm{Re} (z) >0)$
\[
\zeta(z)=\iint
 K(A,B)^z \mathrm{d}A\mathrm{d}B.
\]
Because of Theorem \ref{thm-zeta}, $\zeta(z)$ can be analytically continued to a unique meromorphic function on the entire complex plane $(z \in \mathbb{C})$ and all of its poles are negative rational numbers.
The RLCT of LDA is defined by $\lambda$ if the largest pole of $\zeta(z)$ is $(-\lambda)$.
Its multiplicity $m$ is defined as the order of the maximum pole.
\end{defi}

In order to prepare to state Main Theorem,
we define an intrinsic value $r$ of the true distribution defined by stochastic matrices $A_0$ and $B_0$.
Let $\tilde{A}_0^{\setminus (M,H_0)}$ and $\tilde{B}_0^{\setminus (H_0, 1)}$ be matrices defined as
\begin{equation}
\tilde{A}_0^{\setminus (M,H_0)} = (a_{ik}^0)_{i=1, k=1}^{M-1, H_0-1}
\end{equation}
and
\begin{equation}
\tilde{B}_0^{\setminus (H_0, 1)} = (b_{kj}^0)_{k=1, j=2}^{H_0-1, N},
\end{equation}
respectively. Also, $\tilde{A}_0^{(\setminus M), H_0}$ and $\tilde{B}_0^{(\setminus H_0), 1}$ are denoted by matrices
whose column vectors are same such that
\begin{equation}
\tilde{A}_0^{(\setminus M), H_0} = (a^0_{(1:M-1)H_0}, \ldots a^0_{(1:M-1)H_0}) = (a^0_{iH_0})_{i=1,k=1}^{M-1, H_0-1}
\end{equation}
and
\begin{equation}
\tilde{B}_0^{(\setminus H_0), 1} = (b^0_{(1:H_0-1)N}, \ldots b^0_{(1:H_0-1)1}) = (b^0_{k1})_{k=1,l=1}^{H_0-1, N-1},
\end{equation}
where $a^0_{(1:M-1)H_0}=(a^0_{iH_0})_{i=1}^{M-1}$ and $b^0_{(1:H_0-1)1}=(b^0_{k1})_{k=1}^{H_0-1}$ are an ($M-1$)-dimensional vector and an ($H_0-1$)-dimensional vector, respectively.
Then, let
\begin{align}
U_0 &=\tilde{A}_0^{\setminus (M,H_0)} - \tilde{A}_0^{(\setminus M), H_0}, \label{defU0} \\
V_0 &= \tilde{B}_0^{\setminus (H_0, 1)} - \tilde{B}_0^{(\setminus H_0), 1}, \label{defV0}
\end{align}
and $r=\mathrm{rank} (U_0V_0)$. Obviously, $r$ depends on $H_0$.

The main result of this paper is the following theorem.

\begin{thm}[Main Theorem]
\label{thm-main}
Suppose $M\geqq 2$, $N \geqq 2$, and $H \geqq H_0 \geqq 1$.
Let $r$ be the rank of $U_0V_0$ which is a product of two matrices $(U_0, V_0)$ defined in Eqs. (\ref{defU0}) and (\ref{defV0}).
The RLCT of LDA $\lambda$ and its multiplicity $m$ are as follows.

\begin{enumerate}
\item If $N+r+1 \leqq M+H$ and $M+r+1\leqq N+H$ and $H+r+1\leqq M+N$,
    \begin{enumerate}
        \item and if $M+N+H+r$ is odd, then
        $$\lambda = \frac{1}{8}\{ 2(H+r+1)(M+N)-(M-N)^2-(H+r+1)^2 \} - \frac{1}{2}N, \ m = 1.$$
        \item and if $M+N+H+r$ is even, then
        $$\lambda = \frac{1}{8}\{ 2(H+r+1)(M+N)-(M-N)^2-(H+r+1)^2 +1 \} - \frac{1}{2}N, \ m = 2.$$
    \end{enumerate}
\item Else if $M+H<N+r+1$, then
$$\lambda = \frac{1}{2}\{ MH+N(r+1)-H(r+1)-N \}, \ m = 1.$$
\item Else if $N+H<M+r+1$, then
$$\lambda = \frac{1}{2}\{ NH+M(r+1)-H(r+1)-N \}, \ m = 1.$$
\item Else (i.e. $M+N<H+r+1$), then
$$\lambda = \frac{1}{2}(MN-N), \ m = 1.$$
\end{enumerate}
\end{thm}

To prove Main Theorem, we use the RLCT of matrix factorization (MF).

\begin{defi}[RLCT of MF]
Let $U$, $V$, $U_0$ and $V_0$ be $M \times H$, $H \times N$, $M \times H_0$ and $H_0 \times N$ real matrix, respectively.
Set $M \geqq 1$, $N \geqq 1$, $H \geqq H_0 \geqq 0$. Assume that they in a compact subset $W$ of $(M+N)H$-dimensional Euclidean space and put $r = \mathrm{rank}(U_0V_0)$.
The RLCT of MF $\lambda_{\mathrm{MF}} = \lambda_{\mathrm{MF}}(M,N,H,r)$ is defined by the negative maximum pole of the following zeta function
$$\zeta_{\mathrm{MF}}(z) = \iint_W \lVert UV -U_0V_0 \rVert^{2z} \mathrm{d}U\mathrm{d}V.$$
Its multiplicity $m_{\mathrm{MF}}$ is defined as the order of the maximum pole.
\end{defi}

The exact value of $\lambda_{\mathrm{MF}}$ had been clarified as that of reduced rank regression (a.k.a. three-layered linear neural network) by Aoyagi in \cite{Aoyagi1}. By making the RLCT of LDA come down to that of MF, we prove Main Theorem.
The rigorous proof of Main Theorem is described in Appendix.

\begin{proof}[Sketch of Proof]
Let $\zeta_{\mathrm{SMF}}(z)$ be a zeta function of learning theory in stochastic MF:
\begin{equation}
\zeta_{\mathrm{SMF}}(z) = \iint \lVert AB - A_0B_0 \rVert^{2z} \mathrm{d}A\mathrm{d}B.
\end{equation}
Let $\lambda_{\mathrm{SMF}}$ and $m_{\mathrm{SMF}}$ be the negative maximum pole and its order of
$\zeta_{\mathrm{SMF}}(z)$, respectively.
According to \cite{Matsuda1-e} and \cite{nhayashi7}, there exists two positive constants $C_1>0$ and $C_2>0$ such that
\begin{equation}
C_1 K(A,B) \leqq \lVert AB - A_0B_0 \rVert^2 \leqq C_2 K(A,B),
\end{equation}
i.e.
\begin{gather}
\lambda = \lambda_{\mathrm{SMF}}, \
m = m_{\mathrm{SMF}}
\end{gather}
hold; thus, we only have to consider $\lambda_{\mathrm{SMF}}$ and $m_{\mathrm{SMF}}$.

Developing $\lVert AB - A_0B_0 \rVert^2$,
performing several changes of variables
and considering the integral range of transformed variables in the zeta function,
we have
\begin{gather}
\lambda_{\mathrm{SMF}} = \frac{M-1}{2} + \lambda_{\mathrm{MF}}(M-1,N-1,H-1,r), \
m_{\mathrm{SMF}} = m_{\mathrm{MF}}.
\end{gather}

We calculate $\lambda_{\mathrm{MF}}$ and $m_{\mathrm{MF}}$ by following \cite{Aoyagi1},
then we obtain Main Theorem.
\end{proof}

\section{Experiment}

Now, we run numerical experiments to check the behavior of Main Theorem when the sample size is finite.
Theorem \ref{thm-main} gives the exact asymptotic form of Bayesian generalization error in LDA by using Eq. (\ref{thm-watanabeG}).
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 \cite{nhayashi8} since it also treats the numerical experiment to compute the RLCT by using Gibbs sampling to verify the numerical behavior of its theoretical result.

Let $\mathbb{E}_w$ be an expectation operator of the posterior: $\mathbb{E}_w[\cdot]=\int \mathrm{d}w \varphi^*(w|\mathcal{D})[\cdot]$. Let $\hat{\lambda}$ be the numerically calculated RLCT.
The widely applicable information criterion (WAIC) \cite{WatanabeAIC} is defined by the following random variable $W_n$:
\begin{equation}
W_n = T_n + V_n / n,
\end{equation}
where $T_n$ is the empirical loss and $V_n$ is the functional variance:
\begin{align}
T_n &= -\frac{1}{n} \sum_{i=1}^n p^*(X_i) = -\frac{1}{n} \sum_{i=1}^n \mathbb{E}_w [p(X_i | w)], \\
V_n &= \sum_{i=1}^n \left[ \mathbb{E}_w [(\log p(X_i|w))^2] - \left\{ \mathbb{E}_w [ \log p(X_i|w)] \right\}^2 \right]
= \sum_{i=1}^n \mathbb{V}_w[\log p(X_i | w)].
\end{align}
Even if the posterior distribution cannot be approximated by any normal distribution (i.e., the model is singular),
the expected WAIC $\mathbb{E}[W_n]$ is asymptotically equal to the expected generalization loss $\mathbb{E}[G_n + S]$ \cite{WatanabeAIC};
\begin{equation}
\mathbb{E}[G_n + S] = \mathbb{E}[W_n] + o(1/n^2).
\end{equation}
Moreover, the generalization error and the WAIC error $W_n - S_n$ have the same variance \cite{WatanabeAIC, SWatanabeBookMath}:
\begin{align}\label{var_RLCT}
G_n + W_n - S_n = 2\lambda /n + o_p(1/n).
\end{align}
We need to repeat the simulations to compute $\hat{\lambda}$ to decrease the random effect caused by $G_n$, $W_n$ and $S_n$.
Thus, Eq. (\ref{var_RLCT}) is useful for computing $\hat{\lambda}$ because the leading term $2\lambda /n$ is deterministic, 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 Eq. (\ref{thm-watanabeG}).

The method was as follows.
First, the training data $\mathcal{D}$ was generated from the true distribution $q(X|Z)$.
Second, the posterior distribution was calculated by using Gibbs sampling \cite{Griffiths2004LDAGS}.
Third, $G_n$ and $W_n - S_n$ were computed by using the training data $\mathcal{D}$ and the artificial test data $\mathcal{D^*}=(X_t^*)$ generated from $q(X|Z)$.
These three steps were repeated and each value of  $n(G_n + W_n - S_n)/2$ was saved.
After all repetitions have been completed,
$n(G_n + W_n - S_n)/2$ was averaged over the simulations. The average value was $\hat{\lambda}$.

The pseudo-code is listed in Algorithm \ref{exp_alg}, where $K$ is the sample size of the parameter subject to the posterior and $n_T$ is the sample size of the synthesized test data.
We used the programing language named Julia \cite{Bezanson2017julia} for this experiment.

\begin{algorithm}
\caption{How to Compute $\hat{\lambda}$}
\label{exp_alg}
\begin{algorithmic}
\REQUIRE $D$: the number of simulations, \\
$A_0$: the true parameter matrix whose size is $(M, H_0)$,\\
$B_0$: the true parameter matrix whose size is $(H_0, N)$,\\
$\mathrm{GS}$: the Gibbs sampling function whose return value consists of the samples from the posterior.
\ENSURE The numerical computed RLCT $\hat{\lambda}$.
\STATE Allocate an array $\Lambda[D]$.
\FOR{$d=1$ to $D$}
\STATE Generate $\mathcal{D}$ from the true distribution.
\STATE Allocate arrays $\mathcal{A}[M, H, K]$ and $\mathcal{B}[H, N, K]$.
\STATE Get $\mathcal{A}, \mathcal{B} \leftarrow \mathrm{GS}(\mathcal{D})$.
\STATE Generate $\mathcal{D^*}$ from the true distribution.
\STATE Calculate $G_n \approx \frac{1}{n_T} \sum_{t=1}^{n_T} \log \frac{q(X^*_t)}{\mathbb{E}_w [p(X^*_t | w)]}$,
$S_n = -\frac{1}{n} \sum_{i=1}^n \log q(X_i)$,
\STATE and $W_n \approx -\frac{1}{n} \sum_{i=1}^n \mathbb{E}_w [p(X_i | w)] + \frac{1}{n} \sum_{i=1}^n \mathbb{V}_w[\log p(X_i | w)]$,
\STATE $\quad$ where $\mathbb{E}_w[ f(w) ] \approx \frac{1}{K} \sum_{k=1}^K f(w_k)$ and $w_k = (\mathcal{A}[:,:,k], \mathcal{B}[:,:,k])$.
\STATE Save $\Lambda[d] \leftarrow n(G_n + W_n - S_n)/2$.
\ENDFOR
\STATE Calculate $\hat{\lambda} = \frac{1}{D} \sum_{d=1}^D \Lambda[d]$.
\end{algorithmic}
\end{algorithm}

We set
$M=10$, $N=5$, $H_0=2$, $r=1$, $n=1000$, and $n_T = 200n = 200000$.
To examine the numerical behavior when the number of topics $H \geqq H_0$ in the model increases, we set $H=2,3,4,5$ and carried out experiments in each case.
To decrease the probabilistic effect of Eq.~(\ref{var_RLCT}), 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 effect 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 $20000 + 20K = 50000$ ($K=1000$) and we used the $(10000 + 20k)$-th sample as the entry of $\mathcal{A}[:, :, k]$ and $\mathcal{B}[:, :, k]$ for $k=1$ to $K$.

The experimental results are shown in Table~\ref{exp-result} and visualized in Fig.~\ref{fig:exp-visual}.
The symbol $\lambda$ denotes the exact value of the RLCT $\lambda$ in Theorem~\ref{thm-main}.
There are three columns for each $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 $\hat{\lambda}$ is an average of $n(G_n + W_n - S_n)/2$ obtained from each simulations;
hence, we also show the standard deviation of it as the right next to the plus-minus sign $\pm$ in each setting.

\begin{table}[h]\centering
\caption{Numerically-Calculated and Theoretical Values of RLCTs}
  \begin{tabular}{|c|c|c|c|} \hline
   \multicolumn{1}{|l|}{Settings}
      & \multicolumn{1}{l|}{RLCTs}
      & \multicolumn{1}{l|}{Values} \\ \hline

    $H=2$ &  Theoretical: $\lambda$ & 21/2 \\
    $(M=10,N=5$ & Numerical: $\hat{\lambda}$ & $10.79 \pm 0.8591$ \\
    $H_0=2,r=1)$ & Difference: $|\lambda-\hat{\lambda}|$ & $0.2901 \pm 0.8591$ \\ \hline
    $H=3$ &  Theoretical: $\lambda$ & 12 \\
    $(M=10,N=5$ & Numerical: $\hat{\lambda}$ & $12.25 \pm 0.9510$ \\
    $H_0=2,r=1)$ & Difference: $|\lambda-\hat{\lambda}|$ & $0.2534 \pm 0.9510$ \\ \hline
    $H=4$ &  Theoretical: $\lambda$ & 27/2 \\
    $(M=10,N=5$ & Numerical: $\hat{\lambda}$ & $13.57 \pm 1.036$ \\
    $H_0=2,r=1)$ & Difference: $|\lambda-\hat{\lambda}|$ & $0.07114 \pm 1.036$ \\ \hline
    $H=5$ &  Theoretical: $\lambda$ & 15 \\
    $(M=10,N=5$ & Numerical: $\hat{\lambda}$ & $14.80 \pm 1.143$ \\
    $H_0=2,r=1)$ & Difference: $|\lambda-\hat{\lambda}|$ & $0.2049 \pm 1.143$ \\ \hline
  \end{tabular}\label{exp-result}
\end{table}

As shown in Table~\ref{exp-result} and Fig.~\ref{fig:exp-visual}, all numerically calculated values are nearly equal to the theoretical RLCTs, i.e. these differences are sufficiently 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 the number 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.~\ref{fig:exp-visual}, the standard deviations are sufficiently small for the scale of the RLCTs.
Therefore, Theorem~\ref{thm-main} is consistent with the experimental result.

\begin{figure}[h]
\centering
\includegraphics[width=11cm]{./ExpResult-NumTheoRLCT2.eps}
\caption{This figure is drawn based on Table \ref{exp-result} and Theorem~\ref{thm-main}. It compares numerically-calculated RLCTs $\hat{\lambda}$ (Numerical $\lambda^{\hat{}}$, as the dashed yellow line with the error bars) and theoretical ones $\lambda$ (Theoretical $\lambda$, 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 value of the RLCT. Each error bar of experimental results is the 1-standard deviation range. The line of Numerical $\lambda^{\hat{}}$ and that of Theoretical $\lambda$ are very close and the standard deviations are sufficiently smaller than the scale of the RLCTs.}
\label{fig:exp-visual}
\end{figure}

\section{Discussion and Conclusion}
In this paper, we described how the exact RLCT (i.e. learning coefficient) of LDA is determined in the general case.
Main Theorem
Using this result, we also clarified the exact asymptotic forms of the Bayesian generalization error and the marginal likelihood in LDA.

The RLCT of LDA can be represented by using that of MF.
Namely, Main Theorem can be interpreted as that the learning coefficient of LDA is that of the unconstrained MF minus the penalty due to the simplex constraint. In fact, it can be proved that
\begin{equation}
\lambda = \lambda_{\mathrm{MF}}(M,N,H,r+1)-\frac{N}{2}
\end{equation}
holds (see also the rigorous proof of Main Theorem in Appendix).
The dimension of the stochastic matrix $AB$ with the degrees of freedom is $(M-1)N=MN-N$.
The subtracted $N$ is the dimension of the parameter that can be uniquely determined from the parameters of 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 the above 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 this problem in learning theory: what the learning coefficient of LDA is.

Since LDA is a knowledge discovery method, marginal-likelihood-based model selection often tends to be preferred.
However, BIC \cite{Schwarz1978BIC} cannot be used for LDA because it is a singular statistical model.
Although Gibbs sampling is usually used for full Bayesian inference of LDA, it is difficult to achieve a tempered posterior distribution; thus, we need other Markov chain Monte Carlo method (MCMC) to calculate WBIC \cite{WatanabeBIC} and WsBIC \cite{Imai2019estimating}.
The result of this study allows us to perform a rigorous model selection of LDA with sBIC \cite{Drton}, 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 \cite{Nagata2008asymptotic}.
Furthermore, it may be possible to evaluate how precise MCMC approximates the posterior with the exact values of that \cite{SWatanabeBookMath, Imai2019estimating}.

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 difference in values and behaviors between $d/2$ and $\lambda$.
In a regular statistical model, the learning coefficient is half of the parameter 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 $\lambda$ does not.
In addition, $\lambda$ is much smaller than $d/2$.
Fig. \ref{fig:RLCT} shows how the RLCT of LDA $\lambda$ behaves when the number of topics $H$ increases,
with $\lambda$-value in the vertical axis and $H$-value in the horizontal axis.
If $\lambda$ was equal to $d/2$, then it would linearly increase (the square markers dotted plot in Fig. \ref{fig:RLCT}).
However, in fact, $\lambda$ is given by Main Theorem and its curve is obviously non-linear (the circles dotted plot in Fig. \ref{fig:RLCT}). Hence, their values and behaviors are very different.
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 $\lambda$ instead of $d/2$.
That is why sBIC is theoretically recommended for LDA.

We can draw the theoretical learning curve like the solid line in Fig. \ref{fig:LearningCurve},
with $\mathbb{E}[G_n]$-value in the vertical axis and $n$-value in the horizontal axis.
We also namely draw a curve like the dashed line in Fig. \ref{fig:LearningCurve}.
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 $G_n^{\mathrm{MAP}}$ and $\mu$ be the generalization error and the learning coefficient of LDA in maximum posterior methods, respectively. This is well-defined, i.e. $\mathbb{E}[G_n^{\mathrm{MAP}}] = \mu /n + o(1/n)$ holds.
On the basis of the same prior distribution, Watanabe proved the following inequality \cite{SWatanabeBookMath}:
\begin{equation}
\lambda < d/2 < \mu.
\end{equation}
This means $\mathbb{E}[G_n^{\mathrm{MAP}}] > \mathbb{E}[G_n]  + o(1/n)$
and the leading term of these difference is $(\mu - \lambda)/n > (d-2\lambda)/2n$.
Owing to Main Theorem, we immediately have the exact value of $d-2\lambda$.
Therefore, our result shows at least how much Bayesian inference improves the generalization performance of LDA compared to maximum posterior method.
If the prior distribution is a uniform one, then $\mu$ equals the learning coefficient of LDA in maximum likelihood estimation.
Hence, the above consideration can be applied to maximum likelihood estimation.

One of future works is finding simultaneous resolution of singularities when the prior is a Dirichlet distribution.
A density function of a non-uniform Dirichlet distribution has zero or diverged points; thus, it may affect the learning coefficient.