An Analysis of State-Relevance Weights and Sampling Distributions
on $L_1$-Regularized Approximate Linear Programming Approximation Accuracy

Gavin<EMAIL_ADDRESS>Connor Geer
David Piekut
United States Naval Academy,
572M Holloway Rd., Stop 9F, Annapolis, MD  21402-5002

Reinforcement Learning, Feature Selection, L1 Regularization,
Value Function Approximation

Recent interest in the use of $L_1$ regularization in the use of value
function approximation includes Petrik et al.'s introduction of
$L_1$-Regularized Approximate Linear Programming (RALP).  RALP is unique among
$L_1$-regularized approaches in that it approximates the optimal value
function using off-policy samples.  Additionally, it produces policies which
outperform those of previous methods, such as LSPI.  RALP's value function
approximation quality is affected heavily by the choice of state-relevance
weights in the objective function of the linear program, and by the
distribution from which samples are drawn; however, there has been no
discussion of these considerations in the previous literature.  In this
paper, we discuss and explain the effects of choices in the state-relevance
weights and sampling distribution on approximation quality, using both
theoretical and experimental illustrations.  The results provide insight not
only onto these effects, but also provide intuition into the types of MDPs
which are especially well suited for approximation with RALP.

§ INTRODUCTION

In recent years, the Reinforcement Learning community has paid considerable
attention to creating value function approximation approaches which perform
automated feature selection while approximating a value
function [Kolter & Ng, 2009, Johns et al., 2010, Mahadevan & Liu, 2012, Liu et al., 2012].  This approach frees
researchers from hand-selecting and -tuning feature sets, while greatly
increasing approximation accuracy.  One of these approaches, $L_1$-Regularized
Approximate Linear Programming (RALP) [Petrik et al., 2010, Taylor & Parr, 2012] is unique in that
it results in an approximation of the optimal value function and makes use of
off-policy samples.

However, important aspects of RALP have not been fully explored or explained.
In particular, the objective function of the linear program offers an
opportunity to create a better approximation of the value function in some
regions of the state space at the expense of others.  Optimal policies of many
realistic reinforcement learning problems heavily traffic some parts of the
state space while avoiding others, making an understanding of this flexibility

Therefore, in Section <ref>, we derive a new error bound for the RALP
approximation, which is tighter than those presented by Petrik et
al. icml10 and provides new insight into the result of changing the
state-relevance weights in the objective function.  In addition, this bound
provides an insight into the types of MDPs particularly well suited to the
RALP approach.  Finally, this section provides evidence that rather than
weighting all states equally, as was done by Petrik et al. icml10 and
Taylor and Parr taylor12, states should instead be weighted in
proportion to the stationary distribution under the optimal policy.

Additionally, if the sampling distribution is not uniform across the state
space, the approximation can be heavily affected.  In realistic domains,
sampling is rarely uniform.  Therefore, for RALP to be appropriate for these
domains, it is important to address this lack of analysis and understand how
the approximation is likely to be altered due to the problem's available
sampling scheme.

Section <ref> then discusses the impact on the approximation if
learning is performed using samples drawn other than uniformly from the state
space.  We demonstrate that sampling from a distribution acts as a de
facto alteration of the objective function.  We also discuss the effects of
sampling distributions in light of the bounds from Section <ref>.

The intuition provided by Sections <ref> and <ref> are then
demonstrated experimentally in Section <ref>.  Using a simple,
easily visualized domain, we demonstrate the effect of the various parameters
on approximation quality.

§ NOTATION AND PROBLEM STATEMENT

In this section, we formally define Markov decision processes and linear value
function approximation.  A Markov decision process (MDP) is a tuple
$(\sS,\sA,P,R,\gamma)$, where $\sS$ is the measurable, possibly infinite set
of states, and $\sA$ is the finite set of actions. $P: \sS \times \sS \times
\sA \mapsto [0,1]$ is the transition function, where $P(s'|s,a)$ represents
the probability of transitioning from state $s$ to state $s'$, given action
$a$. The function $R: \sS \mapsto \Re$ is the reward function, and $\gamma,$ a
number between 0 and 1, is the discount factor, representing the comparative
desire for reward at the current time step to the desire for reward at the
next time step.

We are concerned with finding a value function $V$ that maps each state
$s\in\sS$ to the expected total $\gamma$-discounted reward for the
process.  Value functions can be useful in creating or analyzing a policy
$\pi: \sS \times \sA \rightarrow [0,1]$ such that for all
$s\in\sS$, $\sum_{a\in\sA} \pi(s,a) = 1$.  The transition and reward
functions for a given policy are denoted by $P_\pi$ and $R_\pi$.  We denote
the Bellman operator for a given policy as $T_\pi,$ and the max
Bellman operator simply as $T$.  That is, for some state $s\in\sS$:

\[
T_\pi V(s)=R(s)+\gamma\int_{\sS}P(ds'|s,\pi(s))V(s')
\]

\[
TV(s)=\max_{\pi\in\Pi}T_\pi V(s).
\]

We additionally denote the Bellman operator for selecting a particular action
$a$ as

\[
T_a V(s)=R(s)+\int_{\sS}P(ds'|s,a)V(s').
\]

The optimal value function $V^*$ satisfies $T V^*(s) = V^*(s)$ for all

For simplicity, in this paper we will assume no noise exists in the MDP;
results can be easily extended to noisy domains using Taylor and
Parr's taylor12 approach of local smoothing.

Sets of samples, therefore, are defined as $\Sigma \subseteq \{(s,a,r,s' |
s,s'\in\sS, a\in\sA\}$, where $s'$ is the state the agent arrived at given
that it started in state $s$ and took action $a$, and $r=R(s)$.  An individual
sample in the set $\Sigma$ will be denoted $\sigma$, and an element of a
sample $\sigma$ will be denoted with superscripts; that is, the $s$ component
of a $\sigma=(s,a,r,s')$ sample will be denoted $\sigma^s$.

We focus on linear value function approximation for discounted
infinite-horizon problems, in which the value function is represented as a
linear combination of possibly nonlinear basis functions (vectors). For each
state $s$, we define a vector $\Phi(s)$ of features. The rows of the basis
matrix $\Phi$ correspond to $\Phi(s)$, and the approximation space is
generated by the columns of the matrix. That is, the basis matrix $\Phi$, and
the approximate value function $\hat V$ are represented as:

\[
\Phi = \begin{pmatrix} - & \Phi(s_1) & - \\  &
\vdots & \end{pmatrix} \qquad \hat V = \Phi w .
\]

This form of linear representation allows for the calculation of an
approximate value function in a lower-dimensional space, which provides
significant computational benefits over using a complete basis; if the number
of features is small and the environment is noisy, this framework can also
guard against overfitting any noise in the samples.

If we define $\Phi$ to be overcomplete, with potentially far more features
than sampled states, then to receive the above benefits we must perform
feature selection. In this process, a few features are chosen from the set, the
span of which will represent the available linear approximation space.  We can
use $L_1$ regularization to calculate a sparse $w$, in which nearly all
features receive a weight of 0, thereby performing automated feature selection.

§ PREVIOUS WORK

RALP was introduced by Petrik et al. icml10 to extend the
capabilities of the linear programming approach to value function
approximation [d'Epenoux, 1963, Schweitzer & Seidmann, 1985, de Farias & Van Roy, 2003].  Given a set of
samples $\Sigma$, the linear program is defined as follows:

\begin{equation}
\label{eqn:ralp}
\begin{array}{rl}
\displaystyle\min_{w} &\rho^T \Phi w\\
\mbox{s.t.} &T_{\sigma^a}\Phi(\sigma^s) w\leq
\Phi(\sigma^s) w ~~\forall \sigma\in\Sigma \\ &\norm{w}_{-1}
\leq \psi,
\end{array}
\end{equation}

where $\rho$ is a distribution, which we call the state-relevance weights, in
keeping with the (unregularized) Approximate Linear Programming (ALP)
terminology of de Farias and Van Roy deFarias03.
$\norm{w}_{-1}$ is the $L_1$ norm of the vector consisting of all weights
excepting the one corresponding to the constant feature.

This final constraint, which contributes $L_1$ regularization, provides
several benefits.  First, regularization in general ensures the linear program
is bounded, and produces a smoother value function.  Second, $L_1$
regularization in particular produces a sparse solution, producing automated
feature selection from an overcomplete feature set.  Finally, the sparsity
results in few of the constraints being active, speeding the search for a
solution by a linear program solver, particularly if constraint generation is

Other techniques have used $L_1$ regularization in similar ways.
LARS-TD [Kolter & Ng, 2009] and LC-MPI [Johns et al., 2010] both approximate the fixed
point of the $L_1$-regularized LSTD problem.  Mahadevan and
Liu mahadevan12 introduced the use of mirror descent, which has a
computation complexity which allows it to be better suited than many other
approaches for online reinforcement learning problems.  These above approaches
are most reliable when samples are collected on-policy.  Liu et
al. liu12 introduced RO-TD, which converges to an approximation of
the value function of a given policy, even when trained on off-policy samples.

In contrast to the above approaches, RALP provides an approximation to the
value function of the optimal policy, even when samples are drawn from
non-optimal or random policies.  Approximations produced by RALP have bounded
error, and have performed well experimentally in comparison to other
approaches.  Finally, in noisy domains, the well-known weakness of linear
programming approaches to value function approximation can be mitigated or
eliminated using local smoothing [Taylor & Parr, 2012].

This previous work on RALP has largely ignored the state-relevance weights
$\rho$ in the objective function, setting $\rho=\bone$ without discussion.
However, a change in the objective function would obviously affect the
solution of the linear program.  In certain practical situations this would be
useful to understand.
For example, consider the task of calculating a value function for an aircraft
in flight.  The space of possible flight attitudes, velocities, etc. is very
large.  However, the percentage of this space trafficked in non-catastrophic
flight is small; it is likely worthwhile to improve the approximation quality
in this relevant portion, at the expense of accuracy in the remainder of the

Some previous results exist regarding the effects of changing the
state-relevance weights in the closely-related ALP [de Farias & Van Roy, 2003].
However, the assumptions and types of appropriate problems are very different
between ALP and RALP, making these previous results insufficient.  First, ALP
assumes a sample is drawn from every state-action pair, an assumption which is
not required for RALP.  This means it was not necessary with ALP to consider
the behavior of the approximation between samples or in a continuous space.
Furthermore, it was not necessary to consider the effects of sampling
distributions at all.  This assumption was later weakened by a followup
paper [de Farias & Van Roy, 2004], but not to a degree necessary for large or continuous
problems, particularly when large numbers of samples are not available.  The
second difference is ALP is unregularized, simplifying the definition of the
feasible space of the linear program.  Despite these differences, these
previous results will serve as a useful guide.

Additionally, previous work does not cover the effects of sampling schemes on
RALP.  If states are sampled heavily in one portion of the state space, the
linear program will choose a solution which tightens constraints in that
portion of the space over others.  In realistic settings, it can be difficult
to sample uniformly, making it especially important to understand the effect
of other sampling distributions on the resulting approximate value function.

The remainder of this document fills in these gaps in the previous work.

§ STATE-RELEVANCE WEIGHTS

The theoretical results presented on RALP in the literature thus far offer no
insights into the behavior of the approximation as the state-relevance weights
are altered.  Therefore, it is necessary to derive new bounds for RALP which
contain $\rho$ to understand its effects.

The approach we take follows the example of a proof introduced by de Farias
and Van Roy deFarias03 to bound the ALP approximation, but we extend
it match the weaker assumptions of RALP, along with the requirement that the
weights be $\lone$ regularized.

We begin by defining the relevant notation.  In the following definitions, we
will use $\Re_+$ to refer to the set of non-negative real numbers.

We introduce an operator $H$, defined by

\[
\]

for all $L: \sS\rightarrow\Re_+$.

Therefore, $(HL)(s)$ represents the expected value of $L$ of the next
state if actions are chosen to maximize $L$.

A non-negative function $L:\sS\rightarrow \Re_+$ is a Lyapunov
function if there exists a subset of states $\sB$ and a $\beta_L<1$
such that for all $s\in \sS\setminus\sB$, $\gamma(HL)(s)\leq \beta_L

For an example of a Lyapunov function defined over a MDP, consider the simple
case of a random walk along the non-negative number line, so that
$s\in\ints_+$.  Assume a single action, in which the probability

\[
\]

\[

$L(s)=s$ and $\sB=\{0\}$, then $L$ is a valid Lyapunov function, because
$L(s)$ is expected to decrease for all $s\neq 0$.

Lyapunov functions are often used to prove stability of Markov processes.
Definition <ref> differs from the definition commonly used for
stability analysis in a few ways.  First, in stability analysis, it is
required that $\sS$ be countable, and that $\sB$ be finite.  We have made
neither of these assumptions, though our bounds will be tightest when $\sB$ is
small.  The second difference is we have added a multiplicative term of
$\gamma$.  Because of these differences, a Lyapunov function as defined in
Definition <ref> may not strictly evidence stability.

Besides stability analysis, Lyapunov functions have also previously appeared
in Reinforcement Learning literature, though in different contexts from our
application [Perkins & Barto, 2003, Rohanimanesh et al., 2004].

In the remainder of this section, we will occasionally refer to the weighted
max-norm, where for a vector $U$ and a function $F$,
$\norm{U}_{\infty,F}=\max_i \abs{U(s_i)\cdot F(s_i)}$, and the weighted $L_1$
norm, where $\norm{U}_{1,F}=\sum_i \abs{U(s_i)\cdot F(s_i)}$.

We will start with the following Lemma.  To conserve space, and because the
proof is similar to one presented by de Farias and Van
Roy deFarias03, we reserve the proof for Appendix <ref>.

Assume samples have been drawn from every possible state-action pair.  Let
$\sW=\{w:\norm{w}_{-1}\leq\psi\},$ and let $w^*=\min_{w\in\sW}\norm{V^*-\Phi
w}_\infty$.  Additionally, for a given Lyapunov function $\Phi w_L$, let

\[
\bar w = w^*+\lnorm{V^*-\Phi w^*}\betafrac w_L.
\]

If a Lyapunov function $\Phi w_L$ is constructed such that $\bar w\in\sW$,

\[
\norm{V^*-\Phi \tilde w}_{1,\rho}\leq \frac{2\rho^T\Phi
w_L}{1-\beta_{\Phi w_L}} \min_{w\in\sW} \lnorm{V^*-\Phi
\]

We note that proving the existence of a Lyapunov function as required in the
above lemma is trivial.  First we construct a weight vector $w_L$ with all
zeros but for a positive weight corresponding to the bias feature; this
results in $\Phi w_L$ being a valid Lyapunov function.  Second, we note that
in this case $\norm{\bar w}_{-1}=\norm{w^*}_{-1}$, meeting the requirement
that $\bar w\in\sW$.

We must now remove the assumption that a sample exists for every state-action
pair.  To enable us to bound the behavior of the value function between
samples, we make the following assumption, similar to the sufficient sampling
assumption made by Petrik et al. icml10:

Assume sufficient sampling, that is, for all $s\in\sS$ and $a\in\sA$, there
exists a $\sigma\in\Sigma$ such that $\sigma^a=a$ and:

\begin{align*}
\norm{\phi(\sigma^s)-\phi(s)}_\infty\leq&\delta_\phi\\
\norm{R(\sigma^s)-R(s)}_\infty\leq&\delta_R\\
\norm{p(s'|\sigma^s,a)-p(s'|s,a)}_\infty\leq&\delta_P~\forall s'\in\sS
\end{align*}

This assumption is not unrealistic.  For example, if the reward function,
basis functions, and transition functions are Lipschitz continuous, then
appropriate values of $\delta_\phi, \delta_R,$ and $\delta_P$ are easily
calculated given the greatest distance between any point in the state space
and a sampled point.

We describe the maximum difference between the RALP solution and the true
solution by using the limits from Assumption <ref> to
demonstrate the following Lemma:

Let $M_1$ be an MDP with optimal value function $V^*_1$, and let $\Sigma$ be
an incomplete set of samples drawn from $M_1$ such that not all state-action
pairs are sampled, but Assumption <ref> is fulfilled.
Therefore, the RALP for $M_1$ has the constraint
$T_{\sigma^a}\Phi(\sigma^s)w\leq\Phi(\sigma^s)w$ for all $\sigma\in\Sigma$,
and the bounded $L_1$ constraint, but is missing all other possible RALP

There exists an MDP $M_2$ with an optimal value function $V^*_2$, identical
in every way to $M_1$ but for the reward function, such that the RALP
solution with no missing constraints is equal to the RALP solution
constructed on $\Sigma$, and $\norm{V^*_1-V^*_2}_{\infty}\leq
\frac{2(\delta_\phi\psi+\delta_R+\delta_P\psi)}{1-\gamma}$.

Proof Sketch: We first show that if $R_1$ and $R_2$ are the
respective reward functions of $M_1$ and $M_2$,

\[
\norm{R_1-R_2}_\infty\leq 2(\delta_\phi\psi + \delta_R+\delta_P\psi).
\]

We then show that if
$\norm{R_1-R_2}_\infty\leq \delta$,

\[
\norm{V^*_1-V^*_2}_{\infty}\leq\frac{\delta}{1-\gamma}.
\]

We leave the details of the proof for Appendix <ref>.

We are now prepared to present the first result of this paper.

Let $\Phi w_L$ be a Lyapunov function as required by Lemma
<ref>.  Define $\Sigma$, MDPs $M_1$ and $M_2$ and their
respective optimal value functions $V^*_1$ and $V^*_2$ as in Lemma
<ref>. Define $\epsilon_p=\delta_\phi\psi+\delta_R+\delta_P\psi$.
Let $\tilde w$ be the RALP solution to $M_1$.

\[
\norm{V_1^*-\Phi \tilde w}_{1,\rho}\leq \frac{2\rho^T\Phi
w_L}{1-\beta_{\Phi w_L}} \min_{w\in\sW} \lnorm{V^*_2-\Phi
\]

Because $\tilde w$ is an optimal solution given all samples from $M_2$, Lemmas
<ref> and <ref> allow us

\[
\norm{V_2^*-\Phi \tilde w}_{1,\rho}\leq \frac{2\rho^T\Phi
w_L}{1-\beta_{\Phi w_L}} \min_{w\in\sW} \lnorm{V_2^*-\Phi
\]

Additionally, Lemma <ref> gave us

Because $\rho$ is a probability distribution,

\[
\norm{V_1^*-V_2^*}_{1,\rho}\leq\norm{V_1^*-V_2^*}_\infty
\leq\frac{2\epsilon_p}{1-\gamma}.
\]

Due to the triangle inequality, Theorem <ref> follows.

§.§ Discussion

This bound is not only tighter than those presented in previous literature,
but also allows us to analyze the RALP approximation quality in new ways.
First, we can observe which Lyapunov functions would result in a better
approximation, and discuss the characteristics of MDPs which allow for those
Lyapunov functions, and therefore lend themselves particularly well to value
function approximation by RALP.  Second, as $\rho$ now appears in our bound,
the bound provides a way of relating our choice of $\rho$ to approximation
quality, allowing for more intuitive and successful parameter assignments.  We
address these in turn.

The Lyapunov function $\Phi w_L$ appears in the first term of our bound in
three places, namely the dot product with $\rho$, the definition of
$\beta_{\Phi w_L}$, and in the norm defining the “optimal" $w$ to which we
compare our approximation.  We first note that the bound becomes smaller as
$\beta_{\Phi w_L}$ decreases.  This suggests that the more stable the MDP, the
better RALP can approximate the value function.

This interpretation of Theorem <ref> leads to other intuitive
explanations.  Consider an MDP with only $L(s)=1~\forall s\in\sS$ as a
Lyapunov function.  Now assume two nearby samples, one where a “good" action
is taken, in the direction of positive reward, and another where a “bad"
action is taken, in the direction of negative reward.  When the linear program
is solved, due to the proximity of the two samples, the constraint
corresponding to the “bad" sample is nearly certain to be loose, and may as
well be removed.  However, if the MDP is highly stable, then these two
extremely different samples would be unlikely, and both constraints are
candidates for being tight.  Therefore, more samples from an MDP with a small
$\beta_L$ are likely to be involved in defining the feasible space,
potentially resulting in an improved approximation.

The appearance of the Lyapunov function in the norm of the bound indicates the
bound is tighter when the feature space $\Phi$ allows for a close
approximation in areas where the Lyapunov function is small.  Bertsimas et
al. bertsimas98 demonstrated that a small Lyapunov function value
correlates in expectation with a higher probability in the stationary
distribution of a Markov chain.  This is particularly interesting when
considering the appearance of the dot product between the Lyapunov function
and $\rho$.  This dot product makes it apparent that the approximation
improves when $\rho$ is large only where $\Phi w_L$ is small.  This provides
evidence that the stationary distribution of the MDP under the optimal policy
may be an advantageous setting for $\rho.$  This evidence meshes well with the
intuition that greater accuracy is most useful in frequently-visited

§ SAMPLING DISTRIBUTION

Imagine an MDP with a small finite state space and a single action.  Ideal
sampling would provide a single sample from each state, giving us an objective
function of $\sum_{s\in\sS}\rho(s)\Phi(s)w$.  However, if sampling from a
distribution across the state space, this ideal situation would be unlikely;
some states would go unsampled, while others would be sampled multiple times.
Because the objective function is defined on samples, this means states that
were sampled multiple times would appear in the objective function multiple
times, causing the linear program to tighten constraints at those states at
the expense of accuracy in other states.

Of course, a similar scenario occurs in infinite state spaces as well.
Multiple states near to each other may be sampled, while other regions have
very few samples; this encourages the linear program to choose features and
an approximate value function which tightens constraints in heavily-sampled
regions at the expense of sparsely-sampled regions.  In this section we
discuss the effects of sampling from an arbitrary distribution over the state
space $\mu$ and ways this can help the researcher understand how to design
sampling methods.

Let $\Sigma_1$ and $\Sigma_\mu$ be sample sets of equal cardinality $N$
drawn from the state space from the uniform distribution and an arbitrary
distribution $\mu$, respectively.  Let $\Phi_1$ and $\Phi_\mu$ be the
feature matrices defined over the states of sets $\Sigma_1$ and
$\Sigma_\mu$.  For any weight vector $w$,

\[
\expect{\mu\Phi_1 w}=\expect{\bone\Phi_\mu w}.
\]

This observation is easy to support; both expectations equal $\int_\sS
\mu(s)\Phi(s) w~ds$.

This means that for any $w$, sampling from a non-uniform distribution $\mu$
provides an equivalent objective function in expectation to sampling from a
uniform distribution but setting the state-relevance weights equal to $\mu$.
We note this is not equivalent to expecting the same approximate value
function as the constraints remain different; this makes the effect of
altering the sampling distribution greater than that of altering the objective
function alone.

Additionally, the bound presented in Theorem <ref> offers an
interpretation of results from sampling from a distribution.  As sampling
becomes less uniform, $\epsilon_p$ likely increases, due to the existence of
larger unsampled regions.  For example, this would happen in the
Lipschitz-continuous case for fulfilling Assumption <ref>.
However, for a sampling distribution which is dense where the Lyapunov value
is small, then the total Lyapunov value in the numerator of the first addend
is small, as well.  Therefore, it may be that the ideal sampling distribution
is one which is dense where the Lyapunov value is low, but still provides
sufficient coverage for $\epsilon_p$ to be small.

§ EXPERIMENTAL RESULTS

In this section, we demonstrate experimentally the conclusions drawn in the
previous sections.  Previous literature has already demonstrated RALP's
effectiveness in common benchmark domains; the purpose of this section
therefore, is to clearly illustrate the conclusions of the previous sections.
In a simple, easily visualized domain, we make a series of comparisons.
First, we compare the approximation accuracy of sampling from a domain with a
stable Lyapunov function to the accuracy resulting from sampling from a domain
without such a function.  Next, we compare the accuracy of the approximation
resulting from sampling uniformly to the accuracy of the approximation
resulting from sampling from two different nonuniform distributions.  Finally,
we compare the approximation accuracy of calculating an approximation with
$\rho=\bone$ to the approximation accuracy of calculating an approximation
when $\rho$ is nonuniform.  This is demonstrated using two different,
nonuniform distributions.

The results were obtained by drawing samples and calculating an approximation
500 times for each compared approach; the error $\abs{V^*-\Phi w}$
was then calculated, and averaged across all 500 trials.  Finally, we
calculate and display the difference between the average errors from the two
approaches.  So, if $\hat V_i^A$ is the approximation from the $i$-th run on
approach A, then when comparing two approaches, A and B, a point on the graphs
of Figure <ref> equals

\[
\sum_{i=1}^{500}\frac{\abs{V^*(s)-\hat
V_i^A(s)}}{500}-\sum_{i=1}^{500}\frac{\abs{V^*(s)-\hat V_i^B(s)}}{500}.
\]

§.§ Domain

Room Domain

All of the experiments were run on a domain defined by a 25 by 25 grid world.
This world included four reward regions in the corners of the grid, each of
which consisted of 9 states, as can be seen in Figure <ref>.
Two reward regions (colored gold) had a reward of 1, while the others had a
reward of -1 (colored red).  The remaining states had a reward of 0. Actions
were to move one square in any of the four directions, unless constrained by a
wall, in which case that action would result in no movement.  The discount
factor $\gamma$ was set to 0.95.  For all trials, the feature set consisted of
symmetric Gaussian features centered around each $\sigma^s$ with variances of
2, 5, 10, 15, 25, 50, and 75, plus the bias feature, resulting in $9n+1$
features for $n$ samples.  The optimal value function can be seen in Figure

This value function is an easy one for RALP to approximate with proper
settings of the regularization parameter and sufficient sampling.  The number
of samples and choices of the regularization parameter $\psi$ were therefore
chosen to illustrate the differences between the results of the compared
methods, not to optimize performance.

[Average error from sampling from stable domain - Average error
from sampling from unstable domain]

[Average error from sampling from $\bone$ - Average error from
sampling from $\zeta$]

[Average error from $\rho=\bone$ - Average error from

[Average error from sampling from $\bone$ - Average error from
sampling from $(\bone-\zeta)$]

[Average error from $\rho=\bone$ - Average error from

Difference in Average Error

§.§ Lyapunov Stable Domain

First, we demonstrate the improvement in RALP's approximation when the domain
has a stable Lyapunov function.  A stable Lyapunov function can be created by
forcing the actor into a defined area in the state space.  In order to keep
the representational difficulty of the optimal value functions the same, we
created a Lyapunov function by eliminating actions which move the actor
further from the nearest positive reward.  This preserves the optimal policy
of the unaltered domain, keeping the optimal value functions identical, making
approximation accuracy a fair comparison.

In the domain without a stable Lyapunov function, the actor was free to move
in the state space based on a random choice among the four actions. However,
in the domain with a stable Lyapunov function, the actor was only allowed to
move in the two directions which would not move it further from the nearest
goal.  We note that not all remaining actions are optimal, so sampling still
includes off-policy samples.

This creates a Lyapunov function where $L(s)$ equals the Manhattan distance
from $s$ to the nearest of state (1,1) or (25,25), and
$\sB=\{(1,1),(25,25)\}$.  In each trial, we uniformly sampled 20 samples.  The
regularization parameter $\psi$ was set to 0.2 and $\rho$ was $\bone$.

The result of subtracting the average errors of the approximation from the
domain with a Lyapunov function from the domain without a Lyapunov function
can be seen in Figure <ref>.  Therefore, positive values
indicate higher error from the domain without a Lyapunov function.

The results show that an approximation learned from samples drawn from a
stable domain is more accurate than an approximation learned from a less
stable domain everywhere except for where $s$ is near equal distance from the
two goal states.  This reinforces the intuition from Subsection
<ref> that samples drawn from a
stable domain are more effective than those that are not, particularly near
the most heavily-visited regions.

§.§ Sampling from a Nonuniform Distribution

Next we illustrate the change in approximation accuracy when sampling from a
nonuniform distribution $\mu$.  Section <ref> presents evidence that
a distribution which is most dense where $L(s)$ is smallest may be
advantageous.  To create such a distribution, an agent was started at a random
state, and was allowed to take the optimal policy for 25 steps.  This was done
10,000 times, and the number of visits to each state was tabulated and
normalized.  This defined our distribution, which was heaviest on the edges
and reward corners, and otherwise slightly increasing with increasing
proximity to the positive reward regions.  We will refer to this distribution
as $\zeta$.

For these trials, 20 samples were drawn per run from the domain with the
stable Lyapunov function as discussed in Subsection <ref>.
The regularization parameter $\psi$ was set to 1.5, and $\rho=\bone$.  The
average error from the 500 uniformly sampled runs was subtracted from the
average error from the 500 runs with $\mu=\zeta$; the result can be seen
in Figure <ref>.  Because the error from the nonuniform sampling
was subtracted from the error from uniform sampling, the positive difference
indicates the results from sampling from $\zeta$ were superior.  The results show the distribution met the goal from Section <ref>;
sampling was varied enough to keep $\epsilon_p$ low, while dense enough in the
areas where $L(s)$ was small.

We then subtracted $\zeta$ from $\bone$ and normalized, making a distribution
we will refer to as $\bone-\zeta$, which was largest where an agent was least
likely to traverse in the stable domain.  We subtracted the average error from
sampling from $\mu=\bone-\zeta$ from the average error from sampling
uniformly, to produce Figure <ref>.  A positive value would
indicate larger error from the approximations on uniformly sampled states.
However, through the entirety of the state space, and particularly in the most
trafficked areas, sampling from $\bone-\zeta$ gave us an inferior result.

This provides evidence for the conclusions of Section <ref> that a
sampling distribution which is densest in the areas where the Lyapunov
function is smallest would produce the best approximations.

§.§ Changing the State-Relevance Weights

Lastly, we illustrate the effect on the approximation of changing the
state-relevance weights.  200 samples were drawn uniformly from the state
space; whereas the effects from the previous two experiments are most
pronounced when samples are sparse, the effects from altering $\rho$ are most
pronounced when a number of constraints can be tightened in a given region.

We would prefer to set $\rho$ to the stationary distribution; however, because
the domain is not recurrent, this is not an option.  However, the distribution
$\zeta$ created for sampling in Subsection <ref> is large
where the Lyapunov function is small, making it a reasonable replacement.  In
one set of trials, $\rho(s)$ was set to the value of this distribution at that
state; in the other, $\rho=\bone$.  $\psi$ was set to 4.  Average error from
the approximation resulting from a nonuniform $\rho$ was subtracted from
average error from the approximation resulting from a uniform $\rho$.
Therefore, a positive value indicates a better approximation from the
nonuniform $\rho.$ Sampling was done uniformly from the stable domain.

Figure <ref> shows that nearly the entire state space was more
accurate with a nonuniform $\rho$, except for where $\zeta(s)\approx 0$.  The
difference is small because with a large number of samples, both
approximations were quite accurate.

In addition, we compared the uniform $\rho$ approximation to an approximation
using $\rho=\bone-\zeta$, which is large where the Lyapunov value is large,
resulting in an increased dot product $\rho\Phi w$ in Theorem <ref>.
Again, we subtracted the error of the approximation resulting from
$\rho=\zeta-1$ from the error of the approximation resulting from
$\rho=\bone$, producing Figure <ref>.  A positive value
indicates the nonuniform $\rho$ approximated that state better than did using
a uniform $\rho$.  However, there are few positive values as the use of
$\rho=\bone-\zeta$ resulted in a dramatically inferior approximation, particularly
in areas where $\rho$ was small.

From both figures, it is clear that a higher $\rho$ value in a given portion
of the state space resulted in an improved approximation, particularly if
designed with Theorem <ref> in mind.

§ CONCLUSION

The experimental success of RALP in previous literature, along with its
easily-fulfilled assumptions, suggests promise for its application to
real-life, complicated, and complex domains.  Despite this promise, and
despite the evidence that the effects are dramatic, no theory had been
produced to analyze changes in the approximation quality given changes to the
objective function parameter $\rho,$ or due to differences in sampling
strategies.  These considerations are essential to the use of RALP in the real
world; it is rarely possible to sample uniformly, and the importance of
accuracy across the state space is rarely consistent.

In this paper, we demonstrate the importance of understanding these ideas, and
produce a bound on the approximation error of RALP which is tighter and more
informative than previous bounds.  This bound provides intuition into the
quality of the RALP approximation as a function of state-relevance weights and
sampling distributions.  In addition, we demonstrated that the quality of a
RALP approximation is particularly good when the domain is stable and has a
Lyapunov function with a small $\beta_L$.

Future work remains, particularly in the area of solving the linear program
quickly in the presence of large amount of data.  Though convex optimization
solvers are considered “fast," with large amounts of data, the memory and
time requirements may be too large for realistic use.  Fortunately, it may be
that the structure of the problem, the small percentage of tight constraints,
and small percentage of active features will avail itself to faster, but
equivalent, approaches.

§ ACKNOWLEDGEMENTS

Thank you to the anonymous reviewers for their help in improving this paper.
Additionally, we are grateful for support from the Naval Research Laboratory
Information Management & Decision Architecture Branch (Code 5580), as well as
financial support by the Office of Naval Research, grant numbers
N001613WX20992 and N0001414WX20507.

[Bertsimas et al., 1998]
Bertsimas, Dimitris, Gamarnik, David, and Tsitsiklis, John N.
Geometric bounds for stationary distributions of infinite markov
chains via lyapunov functions.
Technical report, Massachusetts Institute of Technology, 1998.

[de Farias & Van Roy, 2003]
de Farias, Daniela Pucci and Van Roy, Benjamin.
The Linear Programming Approach to Approximate Dynamic Programming.
Operations Research, 2003.

[de Farias & Van Roy, 2004]
de Farias, Daniela Pucci and Van Roy, Benjamin.
On Constraint Sampling for the Linear Programming Approach to
Approximate Dynamic Programming.
Mathematics of Operations Research, pp.  462–478, 2004.

[d'Epenoux, 1963]
d'Epenoux, F.
A probabilistic production and inventory problem.
Management Science, 1963.

[Johns et al., 2010]
Johns, Jeffrey, Painter-Wakefield, Christopher, and Parr, Ronald.
Linear complementarity for regularized policy evaluation and
In Lafferty, J., Williams, C. K. I., Shawe-Taylor, J., Zemel, R.S.,
and Culotta, A. (eds.), Advances in Neural Information Processing
Systems 23, pp.  1009–1017, 2010.

[Kolter & Ng, 2009]
Kolter, J. Zico and Ng, Andrew.
Regularization and Feature Selection in Least-Squares Temporal
Difference Learning.
In Bottou, Léon and Littman, Michael (eds.), Proceedings
of the 26th International Conference on Machine Learning, pp.  521–528,
Montreal, Canada, June 2009. Omnipress.

[Liu et al., 2012]
Liu, Bo, Mahadevan, Sridhar, and Liu, Ji.
Regularized Off-Policy TD-Learning.
In Proceedings of the Conference on Neural Information
Processing Systems (NIPS), 2012.

[Mahadevan & Liu, 2012]
Mahadevan, Sridhar and Liu, Bo.
Sparse Q-Learning With Mirror Descent.
In Conference on Uncertainty in Artificial Intelligence, 2012.

[Perkins & Barto, 2003]
Perkins, Theodore J. and Barto, Andrew G.
Lyapunov Design for Safe Reinforcement Learning.
The Journal of Machine Learning Research, 3:0
803–832, 2003.

[Petrik et al., 2010]
Petrik, Marek, Taylor, Gavin, Parr, Ronald, and Zilberstein, Shlomo.
Feature selection using regularization in approximate linear programs
for markov decision processes.
In Proceedings of the 27th International Conference on Machine
Learning, 2010.

[Rohanimanesh et al., 2004]
Rohanimanesh, Khashayar, Platt, Robert, Mahadevan, Sridhar, and Grupen,
Coarticulation in Markov Decision Processes.
In Advances in Neural Information Processing Systems, 2004.

[Schweitzer & Seidmann, 1985]
Schweitzer, Paul J. and Seidmann, Abraham.
Generalized Polynomial Approximations in Markovian Decision
Journal of mathematical analysis and applications,
1100 (6):0 568–582, 1985.

[Taylor & Parr, 2012]
Taylor, Gavin and Parr, Ronald.
Value function approximation in noisy environments using locally
smoothed regularized approximate linear programs.
In de Freitas, Nando and Murphy, Kevin (eds.), Conference on
Uncertainty in Artificial Intelligence, pp.  835–842, Catalina island,
California, June 2012.

§ PROOF OF LEMMA 1

As stated in Section 4, this proof is very similar to Theorem 3 by de Farias
and van Roy (2003), but we include it for clarity nonetheless.  The Lemma
constructs a point in the feasible space of the linear program which provides
an approximation with bounded error, and then shows the point chosen by RALP
must be no further than that constructed point.  The proof first requires a
series of additional Lemmas.

For any functions $V$ and $\bar V$,

\[
\abs{T\bar V-TV}\leq\gamma \max_\pi P_\pi \abs{\bar V-V}.
\]

For any $V$ and $\bar V$,

\begin{align*}
T\bar V - TV =& \max_\pi(R+\gamma P_\pi \bar V)-\max_\pi(R+\gamma P_\pi V)\\
=&R+\gamma P_{\pi_{\bar V}} \bar V-R-\gamma P_{\pi_V} V \\
\leq&\gamma \max_\pi P_\pi (\bar V - V)\\
\leq&\gamma \max_\pi P_\pi \abs{\bar V - V},
\end{align*}

where $\pi_V$ and $\pi_{\bar V}$ represent the greedy policies with respect to
value functions $V$ and $\bar V$.  By reversing the terms, we can show
$TV-T\bar V\leq\gamma \max_\pi P_\pi \abs{\bar V - V}$, leading to our result.

For any vector $L$ with positive components and any vector $V$,

\[
TV\leq V+(\gamma HL+L)\|V-V^*\|_{\infty,\frac{1}{L}}.
\]

Note that

\[
\abs{V^*(s)-V(s)}\leq\|V-V^*\|_{\infty,\frac{1}{L}}V(s).
\]

Because of Lemma <ref>,

\begin{align*}
\leq&\gamma\max_\pi\sum_{s'\in\sS}P_\pi(s,s')\abs{V(s')-V^*(s')}\\
\leq&\gamma\|V-V^*\|_{\infty,\frac{1}{L}}\max_{a\in\sA}\sum_{s'\in\sS}P_a(s,s')L(s')\\
\end{align*}

Define $\epsilon=\|V-V^*\|_{\infty,\frac{1}{L}}$.

\begin{align*}
TV(s)\leq& V^*(s)+\gamma\epsilon(HL)(s)\\
\leq& V(s)+\epsilon L(s)+\gamma\epsilon(HL)(s).
\end{align*}

Let $w_L$ be a weight vector such that $\Phi w_L$ is a Lyapunov function, $w$
be an arbitrary weight vector, and

\[
\bar w = w+\lnorm{V^*-\Phi w}\betafrac w_L.
\]

Then, $T\Phi\bar w\leq \Phi \bar w$.

Let $\epsilon=\lnorm{V^*-\Phi w}$.  For any state

\begin{align*}
|(T\Phi&\bar w)(s)-(T\Phi w)(s)|\\
=&\abs{\left(T\left[(\Phi w + \epsilon\betafrac \Phi
w_L\right]\right)(s)-(T\Phi w)(s)}\\
\leq&\gamma\max_\pi \sum_{s'\in\sS}P_\pi(s,s')\\
&\cdot\abs{(\Phi w(s') +
\epsilon\betafrac \Phi w_L(s'))-\Phi w(s')}\\
\leq&\gamma\max_\pi \sum_{s'\in\sS}P_\pi(s,s')\epsilon\betafrac(\Phi
=&\gamma\epsilon\betafrac(H\Phi w_L)(s).
\end{align*}

The first line is a replacement of $\bar w$ with its definition, the second is
due to Lemma <ref>, the third is due to the cancellation of the
two $\Phi w$ terms and the fact that because $\Phi w_L$ is a Lyapunov
function, $\betafrac>0$ (note this is true for states in sets $\sB$ and
$\sS\setminus\sB$ as defined in Definition 1), and the final line is due to
Definition 2.

From this, we can conclude

\[
T\Phi\bar w\leq T\Phi w + \gamma \epsilon\betafrac H\Phi w_L.
\]

We can apply Lemma <ref> to get

\[
T\Phi \bar w \leq \Phi w + \epsilon(\gamma H\Phi w_L+\Phi w_L),
\]

and therefore,

\begin{align*}
T\Phi\bar w \leq& \Phi w + \epsilon(\gamma H\Phi w_L+\Phi w_L) +
\gamma\epsilon\betafrac H\Phi w_L\\
=&\Phi w+\epsilon\betafrac \Phi w_L-\epsilon\betafrac \Phi w_L\\
&+ \epsilon(\gamma H\Phi w_L+\Phi w_L)+ \gamma\epsilon\betafrac H\Phi w_L\\
=&\Phi \bar w-\epsilon\betafrac \Phi w_L+\epsilon(\gamma H\Phi w_L+\Phi
&+ \gamma\epsilon\betafrac H\Phi w_L\\
=&\Phi \bar w+ \epsilon(\gamma H\Phi w_L+\Phi w_L)- \epsilon\betafrac(\Phi w_L-\gamma H\Phi w_L)\\
\leq&\Phi \bar w+ \epsilon(\gamma H\Phi w_L+\Phi w_L)-\epsilon(\Phi w_L+\gamma H\Phi w_L)\\
=&\Phi \bar w.
\end{align*}

The penultimate line can be shown given that $\Phi w_L-\gamma H\Phi w_L >0$

\begin{align*}
\frac{2}{1-\beta_{\Phi
w_L}}-1=&\frac{2}{1-\max_{s\in\sS\setminus\sB}((\gamma(H\Phi w_L)(s))/((\Phi
=&\max_{s\in\sS\setminus\sB}\frac{(\Phi w_L)(s)+\gamma(H\Phi w_L)(s)}{(\Phi w_L)(s)-\gamma(H\Phi w_L)(s)}
\end{align*}

Lemma <ref> demonstrates that all constraints in RALP will be
satisfied by $\bar w$, with the exception of the constraint enforcing the
$L_1$ regularization.  However, we have required even this constraint to be
satisfied by requiring $\bar w\in\sW$.  Therefore, $\bar w$ lies in the
feasible region for RALP.

If every state-action pair is represented with a constraint in the RALP, a
vector $\tilde w$ solves the RALP if and only if it solves

\begin{align*}
\argmin_w \hspace{.3cm}&\norm{V^*-\Phi w}_{1,\rho}\\
s.t. \hspace{.5cm}&T_a\Phi(s) w\leq\Phi(s) w
~~\forall s\in\sS, a\in\sA\\
\hspace{.5cm}&\norm{w_{-1}}_1\leq \psi
\end{align*}

For any policy $\pi$, the Bellman operator $T\pi$ is a contraction in max
norm.  If the Bellman error is one-sided, $T$ is also monotonic.  Therefore,
for any $V$ such that $V\geq TV$,

\[
V\geq TV \geq T^2V \geq V^*.
\]

Therefore, any $w$ that is a feasible solution to a RALP satisfies
$\Phi w\geq V^*$.  From this, we can conclude

\begin{align*}
\norm{V^*-\Phi w}_{1,\rho}=&\sum_{x\in S}\rho(x)\abs{V^*(x)-\Phi(x) w}\\
=&\rho^T\Phi w-\rho^T V^*.
\end{align*}

Because $V^*$ is constant, minimizing $\rho^T\Phi w$ with RALP
constraints is equivalent to minimizing $\norm{V^*-\Phi w}_{1,\rho}$
with RALP constraints.

Given Lemmas <ref> and <ref>, we can finally prove
Lemma 1.

\begin{align*}
\norm{V^*-\Phi \tilde w}_{1,\rho}\leq& \norm{V^*-\Phi\bar w}_{1,\rho}\\
=&\sum_{s\in\sS}\rho(s)\abs{V^*-(\Phi\bar w)(s)}\\
=&\sum_{s\in\sS}\rho(s)(\Phi w_L)(s)\frac{\abs{V^*-(\Phi\bar
w)(s)}}{(\Phi w_L)(s)}\\
\leq&\left(\sum_{s\in\sS}\rho(s)(\Phi w_L)(s)\right)
\max_{s'\in\sS}\frac{\abs{V^*-(\Phi\bar w)(s')}}{(\Phi w_L)(s')}\\
=&\rho^T\Phi w_L\norm{V^*-\Phi\bar w}_{\infty,1/\Phi w_L}\\
\leq&\rho^T\Phi w_L\left(\norm{V^*-\Phi w^*}_{\infty,1/\Phi
w_L}+\norm{\Phi\bar w-\Phi w^*}_{\infty,1/\Phi w_L}\right)\\
\leq&\rho^T\Phi w_L (\norm{V^*-\Phi w^*}_{\infty,1/\Phi
&+\norm{V^*-\Phi w^*}_{\infty,1/\Phi w_L} \left(\frac{2}{1-\beta_{\Phi
w_L}}-1\right)\norm{\Phi w_L}_{\infty,1/\Phi w_L})\\
=&\frac{2\rho^T\Phi w_L}{1-\beta_{\Phi w_L}}\norm{V^*-\Phi
w^*}_{\infty,1/\Phi w_L}.
\end{align*}

The penultimate line is due to the definition of $\bar w$, and the final line
occurs because $\norm{\Phi w_L}_{\infty,1/\Phi w_L}=1$.$\blacksquare$

§ PROOF OF LEMMA 2

When a constraint does not exist in RALP for some state, this does not mean
the value at that state is completely unconstrained; because we bounded the
rate of change of all components of the approximate and true value functions
in Assumption 1, the existence of a constraint constructed on a nearby state
means the existence of what we will call an implied constraint.

This lemma explicitly constructs these implied constraints, and quantifies the
maximum distance from the true constraint which would have existed had that
state been sampled.  It does this by building two MDPs, identical in every
way, but for the reward function.  $M_1$ has been incompletely sampled, with
sample set $\Sigma$.  Every state-action pair therefore has either an explicit
or implied constraint in the corresponding RALP.  $M_2$, however, has been
completely sampled.  Every state-action pair in the set $\Sigma$ is identical
in $M_2$, but all state-action pairs not in $\Sigma$ have a sample producing a
constraint identical to the implied constraints of $M_1$.  Because the
constraints are the same, the RALP solution is the same.  We demonstrate the
difference in the reward functions $R_1$ and $R_2$ is bounded, and thus, the
difference in the optimal value functions $V^*_1$ and $V^*_2$ is bounded.

Given an MDP $M_1$ such that Assumption 1 is true and incomplete sample set
$\Sigma$, an MDP $M_2$ exists such that constructing the RALP with constraints
for all state-action pairs results in an identical RALP solution to that of
the RALP constructed from $\Sigma$, and

Consider an arbitrary state-action pair $s,a$, which is not represented by a
sample in $\Sigma$.  This means we are missing the constraint

\begin{equation}
\label{eqn:wish}
R_1(s)+\gamma \sum_{x\in\sS}\left[p(x|s,a)\Phi(x)\right] w \leq
\Phi(s) w.
\end{equation}

Let us refer to the sample in $\Sigma$ which fulfills the sampling assumption
with $s$ and $a$ as $\sigma$.  We can now construct a bound for how incorrect
each component of this constraint can be if we use the constraint at $\sigma$
and our sampling assumption to replace the missing constraint.  For instance,
the reward function $R(s)$ is easily bounded.

\begin{equation*}
\label{eqn:rewbound}
R_1(\sigma^s)-\delta_R\leq R_1(s) \leq R_1(\sigma^s)+\delta_R
\end{equation*}

We now bound $\Phi(s) w$.  Because the sampling assumption allows each basis
function to change only a finite amount, and because $\norm{ w_{-1}}_1\leq
\psi$, and $\bone(s)=\bone(\sigma^s)$,

\begin{equation*}
\label{eqn:phibound}
\Phi(\sigma^s) w-\delta_\Phi\psi
\leq\Phi(s) w\leq\Phi(\sigma^s) w+\delta_\Phi \psi
\end{equation*}

The final component is $\gamma \sum_{x\in \sS}p(x|s,a)\Phi(x) w$,
which expresses our expected value at the next state.  It will be convenient
to separate the bias feature $\bone$ from the rest of $\Phi$.  We will denote
the remainder of the design matrix as $\Phi_{-1}$, and the weights that
correspond to $\Phi_{-1}$ as $w_{-1}$.  Similarly, we will denote the weight
corresponding to $\bone$ as $w_1$.

\begin{align*}
\sum_{x\in \sS}p(x|s,a)\Phi(x) w=&
\sum_{x\in\sS}p(x|s,a) w_1
+ \sum_{x\in\sS}p(x|s,a)\Phi_{-1}(x) w_{-1}\\
=&  w_1+\sum_{x\in\sS}p(x|s,a)
\Phi_{-1}(x) w_{-1}
\end{align*}

Again, we have bounded the allowable change in our expression of probability.

\begin{align*}
w_1+&\sum_{x\in\sS}p(x|s,a) \Phi_{-1}(x) w_{-1}\\
\leq& w_1+\sum_{x\in\sS}\left[p(x|\sigma^s,\sigma^a)+
\delta_P\right]\Phi_{-1}(x) w_{-1}\\ =&
w_1+\sum_{x\in\sS}p(x|\sigma^s,\sigma^a)\Phi_{-1}(x) w_{-1}+
\delta_P\sum_{x\in\sS}\Phi_{-1}(x) w_{-1}
\end{align*}

Because each basis function $\Phi$ is can be standardized such that
$\norm{\Phi}_1=1$, and because $\norm{w_{-1}}_1\leq \psi$, the second
summation can be at most $\psi$.  So,

\begin{align*}
\label{eqn:pbound}
\sum_{x\in\sS}p(x|\sigma^s,\sigma^a)\Phi(x) w-\delta_P \psi&\leq
\sum_{x\in\sS}\left[p(x|s,a)\Phi(x)\right] w\\
&\leq \sum_{x\in\sS}p(x|\sigma^s,\sigma^a)\Phi(x) w+\delta_P\psi.
\end{align*}

We now combine these results, and construct our implied
constraint to take the place of the missing constraint expressed by Equation
<ref>.  We see that the maximum possible change by the
approximate value function is $\delta_\Phi \psi + \delta_R + \delta_P \psi$.
So, the total cumulative error in the constraint is at most $2(\delta_\Phi
\psi + \delta_R + \delta_P \psi)$.  So, we effectively have the following

\begin{equation*}
\label{eqn:qer}
R_1(s)+q- \gamma \sum_{x\in\sS}\left[p(x|s,a) \Phi(x)\right] w \geq
\Phi(s) w,
\end{equation*}

where $\abs{q}\leq 2(\delta_\Phi \psi + \delta_R + \delta_P \psi)$.

Let $M_2$ be an MDP which is identical in every way to $M_1$, except
$R_2(s)=R_1(s)+q$.  The RALP solution for $M_1$ will be equivalent to the RALP
solution for $M_2$, and $\norm{R_1-R_2}_\infty\leq

An illustration of Lemma 2.  The blue bars are constraints at sampled
points $s_1,s_2\in\sS$.  The red and purple lines indicate the maximum rate
of change of the value function, given our settings of $\psi, \delta_\phi,
\delta_R,$ and $\delta_P$.  The center diamond is therefore the feasible
area for the approximate value function, and the red bar is the
implied constraint at some novel point $s'\in\sS$.  Because
$\epsilon_p=\delta_\Phi\psi+\delta_R+\delta_P\psi$ is the maximum change, we
see that the difference between the best possible setting of $\Phi(s') w$
and the worst possible setting of
$\Phi(s') w$ is at most $2\epsilon_p$.

Let $M_1$ and $M_2$ be MDPs that differ only in their reward vectors $R_1$ and
$R_2$.  Let $V_1^*$ and $V_2^*$ be their optimal value functions.  Then, for

Let $s$ be an arbitrary point in the sets $\sS_1$ and $\sS_2$,
and define $r_{1i}(s)$ and $r_{2i}(s)$ to be the $i$-th reward received in
exploring $M_1$ and $M_2$ from state $s$, respectively.  Note that

\begin{equation*}
\end{equation*}

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

\begin{align*}
\left(\sum_{i=0}^\infty\gamma^i\expect{r_{1i}(s)}+
\sum_{i=0}^\infty\gamma^i\delta\right)\\
\end{align*}

Because this is true for an arbitrary $s$,

Lemma 2 is trivially proven by combining Lemmas <ref> and
<ref>, and is illustrated by Figure <ref>.


In the same manner as in the previous subsection, using (4.4.9)–(4.4.14) and
(4.5.3)–(4.5.5), we can derive

(4.5.6)

[This is a generalization of (4.4.15). Please remind of the proof just after
(4.4.15). The characteristics 2) mentioned there holds also for the vertex,
which can emitts an additonal gauge boson with color $c$, Lorentz index $\nu$,
and the zero momentum.] Next, to generalize (4.5.6) to (4.5.2), we use the
induction with respect to the number of loops; in this case it is only to
repeat the method used in Subsection 4.4 from (4.4.17) to (4.4.20), starting
from the following graph,

[Therefore, we can skip to write the proof here.] Then, the proof of (4.5.3)
has been finished. Here, we have used the fact that the second term in (4.5.2)
has no surface term in the N-dimensional method, which is a generalization of
(4.4.23) in Subsection 4.4. This is depicted in this case as

(4.5.7)

[Please refer to (2.2.8) and (2.2.11)–(2.2.13).] Accordingly, using the
cancellation mechanism by the surface term (4.4.23), (4.5.7), we have the
following formula from (4.5.2),

(4.5.8)

Let us apply this to our process of the fermion scattering by the colorless
external source. That is, we consider the difference of two processes in the
soft momentum limit of two gauge bosons; one process is the emission of two
soft gauge bosons, having the same color $c$ with momenta $k_{\mu}$ and
$l_{\nu}$, while the other process is the emission of colorless two soft
photons with the same momenta. Pictorially, it is

(4.5.9)

Here, (4.5.8) has been used in the deformation of the first line. The meaning
of the limit $k_{\mu},l_{\mu}\to 0$ is to take the limit under the
regularization by the N-dimensional method, as was discussed in the last
subsection. In (4.5.9), two cases of soft gauge boson emissions, one from the
external lines and the other from the internal lines, are not specified. For
example, it includes the following two soft gauge bosons’ emission, which
implies

(4.5.10)

[We should examine the Thomson limit in QCD, using this formula. A difficulty
appears, however, that a soft gauge boson can be easilly resolved into two or
three soft gauge bosons via the self-interaction between the soft gauge
bosons. Therefore, it is difficult to restrict the number of soft gauge bosons
to be one in the initial or in the final state.]

By this (4.5.9), soft gauge boson emission in QCD can be reduced to the soft
photon emission as follows:

(4.5.11)

Apply the Eikonal identities (2.2.14), (2.2.15) to the right hand side of this
equation, we obtain

(4.5.12)

[when using (2.2.14), (2.2.15), the particle with momenta
$l_{1},l_{2},\cdots,l_{n}$, should be considered as a gauge boson with a
color]. In this way (4.4.3) is shown to hold in case of two soft gauge bosons’
emission. Including the external corrections to the external gauge bosons, and
multiplying physical polarization vectors, it becomes (4.4.2). [47]

Here, we will notice by comparison that (4.4.3) obtaiend here and (2.2.7) or
(2.2.17) in Subsection 2.2 have the same form up to the color factor
$[t^{c}t^{c}]_{ji}$. In Section 2.2, this (2.2.7) plays an important role,
which contributes to derive the differential equation (2.2.3) [or (2.2.30)]
controling the heighest infrared divergences in all orders; this differential
equation is easily solved to give (2.2.1) [or (2.2.31)]. Therefore, using
(4.4.3) in QCD, corresponding to (2.2.7) in QED, it is easily expected to be
able to derive the differential equation which controlls the highest inftared
divergences in QCD. [This derivation was given by Cornwall and Tiktopoulos
[20], but the proof of Low’s theorem in (4.4.3) was not given at that time. On
the other hand, Kinoshita and Ukawa [48] develpped a general theory to
estimate the highest logarithmic infrared divergences, and derived the above
differentail equation. Their derivation is, however, too complicated to
understood easilly. The merit of the following derivation of the differential
equation is its simpleness. The new point here is that the Low’s theorem
(4.4.3) has been proved in all orders. Why the proof becomes successful
largely owes the usage of the axial gauge.]

Now, let us repeat the derivation of the differential equaion in QCD, as was
done in QED in Subsection 2.2. There, we have introduced the gauge boson mass
$\lambda$ as the infrared regularization. Then, applying
$\lambda\frac{\partial}{\partial\lambda}$ to the fermionic scattering
amplitude $S$ by the colorless external source, [even if it has the same
notation in Subsection 2.2, but it is here the amplitude in QCD], we obtain,
in the similar way to (2.2.6), the following

$\lambda\frac{\partial}{\partial\lambda}S=\frac{1}{2}\int\frac{d^{N}k}{(2\pi)^{N}}\sum_{c}K^{\mu\nu}_{cc}(p,p^{\prime};k,-k)\lambda\frac{\partial}{\partial\lambda}\tilde{D}_{\mu\nu}(k;\lambda).$
(4.5.13)

[Here, $N$ was used as the reularization of the ultraviolet divergences.] This
$K^{\mu\nu}_{cc}(p,p^{\prime};k,-k)$ in QCD corresponds to Figure (4.5.1) in
QED.

(Figure 4.5.1)

Two leggs of gauge bosons are connected by the dressed propagator, and have
the same color $c$. In the following we proceed by referring to Subsection 2.2
in QED. To examine the infrared divergences of (4.5.13), we have to know the
behavior of $K^{\mu\nu}_{cc}(p,p^{\prime};k,-k)$ in the limit of $k_{\mu}\to
0$ and the behavior of
$\lambda\frac{\partial}{\partial\lambda}\tilde{D}_{\mu\nu}(k;\lambda)$. First,
for the $K^{\mu\nu}_{cc}(p,p^{\prime};k,-k)$, we can use (4.4.3), the Low’s
theorem for two gauge bosons’ emission with the same color, proved in this
subsection. In the present notation it reads

$\displaystyle\sum_{c}K^{\mu\nu}_{cc}(p,p^{\prime};k,l)\underset{k_{\mu},l_{\nu}\to
0}{\approx}+(g_{B})^{2}\sum_{c}[t^{c}t^{c}]$
$\displaystyle\times\left[\frac{p^{{}^{\prime}\mu}}{p^{\prime}\cdot
k}-\frac{p^{\mu}}{p\cdot
k}\right]\left[\frac{p^{{}^{\prime}\mu}}{p^{\prime}\cdot
l}-\frac{p^{\mu}}{p\cdot l}\right]S$ (4.5.14)

Here, the color index $c$ is summed over, so
$\sum_{c}t^{c}t^{c}=C_{2}(R)\bm{1}$ holds, [$\bm{1}$ is a unit matrix].
[Please refer to (3.4.9) and the subsequent explanation.] This corresponds to
(2.2.7). When using (4.5.14) in (4.4.13), it is necessary to have an
attention. It is because in the proof of (4.5.14) in this subsection, we used
the N-dimensional method for both ultraviolet and infrared regularizations. In
deriving the differential equation to control the infrared divergences, we
will keep the N-dimensional method for ultraviolet regularization, but the
different regularization method will be used for the infrared divergences;
introduction of non-zero mass $\lambda$ to the gauge boson for the infrared
regularization. Therefore, it is better to try again the proof of (4.5.14) [or
Low’s theorem (4.4.13)]. However, the difference appears only in (4.4.11)
among the relations (4.4.9)–(4.4.14) for simplest graphs; for which the
$O(\lambda^{2})$ breaking appears in its right hand side. Namely, (4.4.11) is
replaced by the following expression:

(4.5.15)
$\displaystyle\text{where}~{}B^{cab}_{\mu\nu\lambda}=gf^{cab}\frac{1}{(q^{2}-\lambda^{2}+i\varepsilon)^{2}}\frac{\partial}{\partial
q_{\mu}}K_{\nu\lambda}(q),$ (4.5.16) $\displaystyle K_{\nu\lambda}(q)\equiv
g_{\nu\lambda}-\frac{q_{\nu}n_{\lambda}+n_{\nu}q_{\lambda}}{q\cdot n}+(\alpha
q^{2}+n^{2})\frac{q_{\nu}q_{\lambda}}{(q\cdot n)^{2}}~{}[\text{numerator of
the propagator}]$ (4.5.17)

Accordingly, each step of proof and the final result of Low’s theorem have
always $O(\lambda^{2})$ breakings. These breqakings do not finally remain
after taking the limit of the cutoff $\lambda\to 0$, so that we can ignore
these breakings. This proof of Low’s theorem in this subsection corresponds to
(2.2.8)–(2.2.17) in Section 2.2.

Now, we can use (4.5.14) in (4.5.13), and obtain the following differential
equation,

$\displaystyle\lambda\frac{\partial}{\partial\lambda}\ln[S]_{\text{IR}}=$
$\displaystyle-\frac{1}{2}g_{B}^{2}C_{2}(R)\int_{0}^{\delta}\frac{d^{N}k}{(2\pi)^{N}}\left[\frac{p^{{}^{\prime}\mu}}{p^{\prime}\cdot
k}-\frac{p^{\mu}}{p\cdot
k}\right]\left[\frac{p^{{}^{\prime}\nu}}{p^{\prime}\cdot
k}-\frac{p^{\nu}}{p\cdot k}\right]$
$\displaystyle\times\lambda\frac{\partial}{\partial\lambda}\tilde{D}_{\mu\nu}(k;\lambda).$
(4.5.18)

[Here, the part of $k_{\mu}\to 0$ contributes to the highest infrared
singularities.]

This differential equation can be solved easily, and we have

$[S]_{\text{IR}}=e^{-\frac{1}{2}g_{B}^{2}C_{2}(R)\int_{0}^{\delta}\frac{d^{N}k}{(2\pi)^{N}}\left[\frac{p^{{}^{\prime}\mu}}{p^{\prime}\cdot
k}-\frac{p^{\mu}}{p\cdot
k}\right]\left[\frac{p^{{}^{\prime}\nu}}{p^{\prime}\cdot
k}-\frac{p^{\nu}}{p\cdot k}\right]\times\tilde{D}_{\mu\nu}(k;\lambda)}\times
S^{(0)},$ (4.5.19)

where
$S^{(0)}=\bar{u}(p^{\prime}){\ooalign{\hfil/\hfil\crcr$J$}}_{\text{colorless}}u(p)$.
There is the arbitrariness in the initial condition. [(4.5.19) implies that
only the highest divergences are equal in both hand sides and the
arbitrariness for the finite parts is allowed.]

In summary, the highest infrared divergences included in the fermion
scattering amplitude by a colorless source is determined, as in (4.5.18), by
the behavior of the dressed gauge boson propagator,
$\tilde{D}_{\mu\nu}(k;\lambda)$. [Of cource this result holds under the axial
gauge condition.]

Comparing with QED, the results are formally equal up to the group factor
$C_{2}(R)$ [$\frac{N^{2}-1}{2N}$ for $SU(N)$]. The essential difference
exists, however, in $\tilde{D}_{\mu\nu}(k;\lambda)$. That is, in QED $Z_{3}$
is infrared finite [as in (2.2.27) and (2.2.28)], while in QCD the infrared
divergence exists in $Z_{3}$. As was shown by Frenkel and Meuldermans [29] in
the one loop calculation, the infrared divergence included in $\Pi^{(1)}$ near
the mass shell, is

$\Pi^{(1)}(k^{2}=0;N)_{\text{IR}}=-\frac{g^{2}}{2\pi^{2}}c_{2}(G)\frac{1}{(N-4)^{2}},$
(4.5.20)

while $\Pi^{(2)}(k^{2}=0;N)_{\text{IR}}$ has no double pole at $N=4$, so that
it can be ignored, in comparison with $\Pi^{(1)}$. In our regularization
adopted here, this gives

$\Pi^{(1)}(k^{2}=0;\lambda)_{\text{IR}}=-\frac{g^{2}}{2\pi^{2}}c_{2}(G)\left(\ln\frac{m}{\lambda}\right)^{2}$
(4.5.21)

[Please refer to the definition of $\Pi^{(1)},\;\Pi^{(2)}$ in (4.3.27).] Thus,
there is the infrared divergence in $Z_{3}=1/1+\Pi^{(1)}$ (on-shell), not
appearing in QED. Of course, if we perform the renormalization on the mass
shell, we have

$\displaystyle g_{B}^{2}\tilde{D}_{\mu\nu}(k;\lambda)=$
$\displaystyle(g_{B}Z_{3}^{1/2})^{2}\frac{\tilde{D}_{\mu\nu}(k;\lambda)}{Z_{3}}$
$\displaystyle=$ $\displaystyle(g_{R})^{2}\tilde{D}_{\mu\nu}(k;\lambda)_{R},$
(4.5.22)

so that divergence of $Z_{3}$ in $\tilde{D}_{\mu\nu}(k;\lambda)$ is
renormalized into $g_{R}$. If we perform the off-mass-shell renormalization
[as in (4.4.40)–(4.4.43)], however, we have

$g_{B}^{2}\tilde{D}_{\mu\nu}(k;\lambda)=g(\mu)^{2}\tilde{D}_{\mu\nu}(k;\lambda,\mu)\left[\frac{Z_{1}(\mu)}{Z_{2}(\mu)}\right]^{2},$
(4.5.23)

so that the infrared divergence remains in
$\tilde{D}_{\mu\nu}(k;\lambda,\mu)$, not being renormalized into the off-shell
coupling $g(\mu)$. [Here, $\tilde{D}_{\mu\nu}(k;\lambda,\mu)$ is the dressed
gauge boson propagator, renormalizd at a renormalization point $\mu$.]
Because, for $\mu\neq 0$, $Z_{1}(\mu),Z_{2}(\mu),Z_{3}(\mu)$ are infrared
finite. Here again, the characteristic feature of QCD, that the infrared
divergence appears in the charge renormalization, plays the important role. In
short, the problem is reduced to investigate the behavior at $k_{\mu}\sim 0$
of the dressed gauge boson propagator, $\tilde{D}_{\mu\nu}(k;\lambda,\mu)$.
This is the same kind of problem as that of investigating how does $g(\mu)$
behave for $\mu\to 0$ [with referring (3.4.23) and (3.4.24)], so that this
will be the important target from now on.

## Acknowledgements

At the end of the thesis, the author gives his sincere thanks to Dr. Ukawa for
introducing the author to this field, and to Mr. Nakagawa and Mr. Yamamoto for
their collaboration in Reference paper II. He is grateful to thank Dr.
Shimizu, Dr. Higashijima, Dr. Sato, Dr. Arafune, Dr. Fukugita, Mr. Nishimura
(A.), Dr. Nakazawa, Dr. Okamura, and Dr. Okabayasi for their valuable
discussions and advices as well as encouragements to him, since his beginning
to work in this field. He gives his gratitude to Professor Yamaguchi and the
members of the particle physics group at University of Tokyo for their warmly
watching the author since his entrance to the university.

##  5 Appendix: Application of M. Sato’s microfunction to the cancellation
problem of infrared divergences

This Appendix includes the latter half of the author’s report, which was
submitted to Prof. Mikio Sato in August (1976), to get the credit of the
graduate course in physics at U. of Tokyo. The title of the course was
“Microlocal Analysis and Quantum Field Theory” [50].

The first half of the report (Section 1 and 2) was an introductory explanation
of the cancellation of infrared divergences in QED; $d\sigma_{V}+d\sigma_{B}=$
infrared finite, where $d\sigma_{V}$ is the cross section with the virtual
corrections by photon of fermion scattering by an external source, while
$d\sigma_{B}$ is the cross section associated with the emission of an
additonal photons (“B” means Bremmsstrahlung), and hence this part can be
omitted, except for important sentences and equations.

###  5.1 A claim stated in Section 1 of the report

The singular spectrum (S.S.) of a micro function gives the place where the
function becomes not microanalytic. Since the infrared divergence appearing in
the scattering amplitude belongs to the S.S. of a Feynman integral (which is a
micro function), the analysis by S.S. is especially effective, when the
divergent part and its cancellation are dicussed without referring to the
finite part.

###  5.2 Definition of $J_{a}$ and $J_{b}$ given in Section 2 of the report

The cancellation of infrared divergences can be understood by the Cutkosky
rule. For example, consider the forward amplitude, $I$, of the scattering of
an electron by an external source,

$\begin{minipage}{341.43306pt}
\centering\includegraphics[width=284.52756pt,clip]{eqapp.1.eps}
\@add@centering\end{minipage},$ (5.2.1)

where the two diffferent Cutkosky cuts at $a$ and $b$ are taken. The numerator
is common to both cases, so that we consider $J$, which is $I$ up to the
numerator:

$\displaystyle J=\int
dp_{0}^{\prime}\;dq_{0}\;dk_{0}\;\frac{1}{p^{{}^{\prime}2}-m^{2}+i0}\cdot\frac{1}{q^{2}-m^{2}+i0}\cdot\frac{1}{k^{2}+i0}\;\delta(p^{\prime}-q-k).$
(5.2.2)

Application of Cutkosky cuts at $a$ and $b$ gives,

$\displaystyle J_{a}=\int dp_{0}^{\prime}\;dq_{0}\;dk_{0}\;(-2\pi
i)\theta(p^{\prime}_{0})\delta(p^{{}^{\prime}2}-m^{2})\cdot\frac{1}{q^{2}-m^{2}+i0}\cdot\frac{1}{k^{2}+i0}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}\hskip
56.9055pt\times\delta(p^{\prime}-q-k),$ (5.2.3) $\displaystyle J_{b}=\int
dp_{0}^{\prime}\;dq_{0}\;dk_{0}\;\frac{1}{p^{{}^{\prime}2}-m^{2}+i0}\cdot(-2\pi
i)\theta(q_{0})\delta(q^{2}-m^{2})\cdot(-2\pi i)\theta(k_{0})\delta(k^{2})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}\hskip
56.9055pt\times\delta(p^{\prime}-q-k).$ (5.2.4)

If the sum of $J_{a}$ and $J_{a}$ reduces the degree of singularity by one,
then the infrared divergences are cancelled between $a$ (virtual process) and
$b$ (emission process of photon).

###  5.3 Section 3– Application of S.S. of micro function to infrared
divergences

We will estimate the singular spectra (S.S.) of Eq.(5.2.3) and Eq.(5.2.4).
Before that we summarize the formulae for S.S [50]:222(Footnote added in 2022)
The singular spectrum (S.S.) gives the “singularity structure” of a
microfunction, to which Feynman amplitudes belong. Feynman amplitude is given
by a product of propagators with $i\epsilon(=i0)$ rule, which means the
singularities are avoided by shifting the pole positions a little in the
imaginary direction. Thus, the microfunction is defined as a boundary value on
the real axis of an analytic function given in the complex plane. To manifest
the shift in the complex direction of the pole, the differential (cotangent
vector) giving the pole shift in the imaginary direction, is specified
additionally in the definition of S.S.. An example of the differential can be
seen for the propagator as
$\displaystyle\sqrt{-1}\frac{\alpha}{2}d(p^{2}-m^{2})=\sqrt{-1}\left(\alpha\sqrt{\vec{p}^{2}+m^{2}}\right)\left\\{-d(p^{0}+\sqrt{\vec{p}^{2}+m^{2}})+d(p^{0}-\sqrt{\vec{p}^{2}+m^{2}})\right\\},$
(5.3.1) which shows the directions of the shifts for the positive and negative
energy poles are opposite.

$\displaystyle\widehat{S.S.}\left(\frac{1}{p^{2}-m^{2}+i0}\right)=\left\\{\left(p,\sqrt{-1}\frac{\alpha}{2}d(p^{2}-m^{2})\right)\Big{|}\;\alpha\geq
0,\;\alpha(p^{2}-m^{2})=0\right\\},$ (5.3.2)
$\displaystyle\widehat{S.S.}\left(\theta(p_{0})\delta(p^{2}-m^{2})\right)=\left\\{\left(p,\sqrt{-1}\frac{\alpha}{2}d(p^{2}-m^{2})\right)\Big{|}\;\alpha\in\mathbb{R},\;p^{2}-m^{2}=0,\;p_{0}>0\right\\}.~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(5.3.3)

The product and integral formulae are given as follows [50]:

$\displaystyle\widehat{S.S.}\left(f(x)g(x)\right)=\left\\{(x,\sqrt{-1}\xi
dx)|\;\xi=\xi_{1}+\xi_{2}~{}s.t.\right.$ $\displaystyle\hskip
99.58464pt\left.(x,\sqrt{-1}\xi_{1}dx)=\widehat{S.S.}f(x),~{}(x,\sqrt{-1}\xi_{2}dx)=\widehat{S.S.}g(x)\right\\},$
(5.3.4) $\displaystyle\widehat{S.S.}\left(\int
f(x,y)dy\right)=\left\\{(x,\sqrt{-1}\xi dx)\Big{|}~{}\exists
y~{}s.t.~{}(x,y;\sqrt{-1}\xi
dx)\in\widehat{S.S.}f(x,y)\right\\}.~{}~{}~{}~{}~{}$ (5.3.5)

####  5.3.1 $\widehat{S.S.}(J_{a})$

If we choose the coefficients of cotangent vectors in the directions of
$p^{{}^{\prime}2}-m^{2}$, $q^{2}-m^{2}$, $k^{2}$, and $p^{\prime}-q-k$, be
$\frac{\alpha}{2}$, $\frac{\beta}{2}$, $\frac{\gamma}{2}$, and $y^{\mu}$,
respectively, then we have

$\displaystyle\widehat{S.S.}(J_{a})=\left\\{\left(\vec{p^{\prime}},\vec{q},\vec{k};\;\sqrt{-1}[\vec{u}d\vec{p^{\prime}}+\vec{v}d\vec{q}+\vec{w}d\vec{k}]\right)\right\\},$
(5.3.6)

where there exist $\alpha\in\mathbb{R}$, $\beta\geq 0$, $\gamma\geq 0$,
$y^{\mu}\in\mathbb{R},\;p^{{}^{\prime}0},\;q^{0},\;k^{0}$, such that the
following condition holds,

$\displaystyle\vec{u}d\vec{p^{\prime}}+\vec{v}d\vec{q}+\vec{w}d\vec{k}=\frac{\alpha}{2}d(p^{{}^{\prime}2}-m^{2})+\frac{\beta}{2}d(q^{2}-m^{2})+\frac{\gamma}{2}d(k^{2})+\vec{y}d(\vec{p^{\prime}}-\vec{q}-\vec{k}),$
(5.3.7)

which is satisfied, when the following relations are satisfied,

$\displaystyle\vec{u}=\alpha\vec{p^{\prime}}+\vec{y},\;\vec{v}=\beta\vec{q}-\vec{y},\;\vec{w}=\gamma\vec{k}-\vec{y},$
(5.3.8) $\displaystyle 0=\alpha p^{{}^{\prime}0}+y^{0},\;0=\beta
q^{0}-y^{0},\;0=\gamma k^{0}-y^{0},$ (5.3.9) $\displaystyle
p^{{}^{\prime}0}>0,\;p^{{}^{\prime}2}-m^{2}=0;\;\beta(q^{2}-m^{2})=0,\;\gamma(k^{2})=0,$
(5.3.10) $\displaystyle(p^{\prime}-q-k)^{\mu}=0.$ (5.3.11)

Let us examine the conditions Eq.(5.3.8)-Eq.(5.3.11).

Eq.(5.3.9) appears due to the integrations over $p^{{}^{\prime}0},\;q^{0}$,
and $k^{0}$. Eqs.(5.3.10)(5.3.11) give the on shell condition, or means the
singularity arises when the particle flowing the propagator becomes on the
mass shell. The kinematics shows that three momenta, $p^{\prime},q,k$, can not
be on the mass shell at the same time. Therefore, the following two cases
remain,

$\displaystyle(Case1):\beta=0,\;k^{2}=0,~{}~{}(Case2):q^{2}-m^{2}=0,\;\gamma=0.$
(5.3.12)

First we consider (Case2); from Eq.(5.3.9), $y^{0}=0$ and $\alpha
p^{{}^{\prime}0}=\beta q^{0}=0$. Since on the mass shell conditions tell
$p^{{}^{\prime}0},\;q^{0}\neq 0$, we have

$\displaystyle\alpha=\beta=\gamma=y^{0}=0.$ (5.3.13)

This implies the corresponding cotangent vector (the direction to avoid the
singularity) vanishes,

$\displaystyle\vec{u}d\vec{p^{\prime}}+\vec{v}d\vec{q}+\vec{w}d\vec{k}=\vec{0},$
(5.3.14)

so that $J_{a}$ is microanalytic in this case.

Thus we consider (Case1); from Eq.(5.3.9), $y^{0}=0$ and since
$p^{{}^{\prime}0}\neq 0$, we have

$\displaystyle\alpha=0,\;k^{0}=0,$ (5.3.15)

where $\gamma=0$ holds as before.

Now we have

$\displaystyle\alpha=\beta=y^{0}=0;\;|\vec{k}|=0,\;\vec{u}=\vec{v}=\vec{0},\;\vec{w}=\gamma\vec{k}.$
(5.3.16)

Of course, if $|\vec{k}|=0$, then we have $\vec{w}=\vec{0}$, and $J_{a}$
becomes microanalytic. This happens because we have not do the regularization.
In reality we have to introduce the regularization (photon mass) and consider

$\displaystyle|\vec{k}|\sim\lambda,\;\vec{w}=\gamma\vec{k}.$ (5.3.17)

Summarizing the above consideration, we have arrived at

$\displaystyle\widehat{S.S.}(J_{a})=\left\\{\left(\vec{p^{\prime}},\vec{q},\vec{k};\sqrt{-1}\gamma(\vec{k}d\vec{k})\right)\big{|}\;\gamma>0,\;|\vec{k}|\sim\lambda,\;k^{2}=0,\;\vec{p^{\prime}}=\vec{q}+\vec{k}\right\\}.$
(5.3.18)

Taking into account the supplementary condition of

$\displaystyle\alpha=\beta=y^{0}=0,$ (5.3.19)

the result of $\widehat{S.S.}(J_{a})$ becomes quite reasonable from the
physical view point.

We will examine this in the following: $S.S(J_{a})$ arises from the region of
$\alpha=\beta=0,\;\gamma\neq 0,\;|\vec{k}|\sim\lambda,\;k^{2}=0$. If we
consider $\alpha,\;\beta,\;\gamma$ be “resistance” in the circuit analogy,
$\alpha=\beta=0$ implies that the large electric currents $p^{\prime},\;q$
flow along the corresponding paths, while $\gamma\neq 0$ implies the small
current $|\vec{k}|\sim\lambda$ flows along the path, which gives the region of
infrared divergence in the limit of soft photon. That is,

$\displaystyle\widehat{S.S}(J_{a})\Leftrightarrow\mathrm{Infrared~{}divergence~{}region~{}of}~{}J_{a}.$
(5.3.20)

####  5.3.2 $\widehat{S.S.}(J_{b})$

In the same manner, $\widehat{S.S}(J_{b})$ can be obtained. Namely we have

$\displaystyle\widehat{S.S.}(J_{b})=\left\\{\left(\vec{p^{\prime}},\vec{q},\vec{k};\;\sqrt{-1}[\vec{u}d\vec{p^{\prime}}+\vec{v}d\vec{q}+\vec{w}d\vec{k}]\right)\right\\},$
(5.3.21)

where there exist $\alpha\geq 0$, $(\beta,\gamma,y^{0})\in\mathbb{R}$,
$p^{{}^{\prime}0},\;q^{0},\;k^{0}$, such that the following conditions hold,

$\displaystyle\vec{u}=\alpha\vec{p^{\prime}}+\vec{y},\;\vec{v}=\beta\vec{q}-\vec{y},\;\vec{w}=\gamma\vec{k}-\vec{y},$
(5.3.22) $\displaystyle 0=\alpha p^{{}^{\prime}0}+y^{0},\;0=\beta
q^{0}-y^{0},\;0=\gamma k^{0}-y^{0},$ (5.3.23)
$\displaystyle\alpha(p^{{}^{\prime}2}-m^{2})=0;\;q^{0}>0,\;q^{2}-m^{2}=0;\;k^{0}>0,\;k^{2}=0,$
(5.3.24) $\displaystyle(p^{\prime}-q-k)^{\mu}=0.$ (5.3.25)

Let us examine these conditions. From Eq.(5.3.24) we have
$p^{{}^{\prime}2}\neq m^{2}$, so $\alpha=0$, which implies $y^{0}=0$ and
$\beta=0$ and $\gamma k^{0}=0$.

In the same manner as before, if $\gamma=0$, $J_{b}$ is microanalytic, so that
we consider $\gamma\neq 0$. Since $k^{0}\sim\lambda$ and $k^{2}=0$, we have
$|\vec{k}|\sim\lambda$. Summarizing this, we have $\alpha=\beta=\gamma=0$,
yielding

$\displaystyle\widehat{S.S.}(J_{b})=\left\\{\left(\vec{p^{\prime}},\vec{q},\vec{k};\sqrt{-1}\gamma(\vec{k}d\vec{k})\right)\big{|}\;\gamma\neq
0,\;|\vec{k}|\sim\lambda,\;k^{2}=0,\;\vec{p^{\prime}}=\vec{q}+\vec{k}\right\\}.$
(5.3.26)

This result implies $p^{\prime}$ and $q$ are hard and $k$ is soft, giving the
infrared divergence of $J_{a}$. Schematically, we have

$\displaystyle\widehat{S.S}(J_{b})\Leftrightarrow\mathrm{Infrared~{}divergence~{}region~{}of}~{}J_{b}.$
(5.3.27)

It is important to note that $\widehat{S.S.}(J_{a})$ and
$\widehat{S.S.}(J_{b})$ overlap with each other, so that the sum of these two
can reduce the degree of infrared divergence.

####  5.3.3 Sum of $J_{a}$ and $J_{b}$

Instead of making a sum of $J_{a}$ and $J_{b}$, consider the forward
scattering amplitude $J$ without cut at $a$ or $b$. Tnen, we have

$\displaystyle\widehat{S.S.}(J)=\left\\{\left(\vec{p^{\prime}},\vec{q},\vec{k};\;\sqrt{-1}[\vec{u}d\vec{p^{\prime}}+\vec{v}d\vec{q}+\vec{w}d\vec{k}]\right)\right\\},$
(5.3.28)

where there exist $(\alpha,\beta,\gamma)\geq
0,\;y^{0}\in\mathbb{R};\;p^{{}^{\prime}0},\;q^{0},\;k^{0}$, such that the
following conditions hold,

$\displaystyle\vec{u}=\alpha\vec{p^{\prime}}+\vec{y},\;\vec{v}=\beta\vec{q}-\vec{y},\;\vec{w}=\gamma\vec{k}-\vec{y},$
(5.3.29) $\displaystyle 0=\alpha p^{{}^{\prime}0}+y^{0},\;0=\beta
q^{0}-y^{0},\;0=\gamma k^{0}-y^{0},$ (5.3.30)
$\displaystyle\alpha(p^{{}^{\prime}2}-m^{2})=0;\;\beta(q^{2}-m^{2})=0;\;\gamma
k^{2}=0,$ (5.3.31) $\displaystyle(p^{\prime}-q-k)^{\mu}=0.$ (5.3.32)

Let us examine these conditions. We can consider four cases.

(Case1): $\alpha=\beta=\gamma=0$ gives the zero cotangent vector, so that this
case can be dropped.

(Case2): $p^{\prime}$ is on the mass shell, $p^{{}^{\prime}2}-m^{2}=0$, and
others are off shell, $q^{2}-m^{2}\neq 0$, $k^{2}\neq 0$. This case also gives
the zero cotangent vector and is dropped.

(Case3): $q$ is on the mass shell, but others are off shell,
$p^{{}^{\prime}2}-m^{2}\neq 0$, and $k^{2}\neq 0$. This is also irrelevant.

(Case4): $k$ is on the mass shell $k^{2}=0$, and the others are off shell,
$p^{{}^{\prime}2}-m^{2}\neq 0$ and $q^{2}-m^{2}\neq 0$. This implies
$\alpha,\beta=0$, $y^{0}=0,\;\gamma k^{0}=0$. Choose $\gamma\neq 0$ and
$k^{0}\sim\lambda$, we have $|\vec{k}|\sim\lambda$. This can give the infrared
divergence.

(Case5): $p^{{}^{\prime}2}=m^{2},\;q^{2}=m^{2}$, but $k^{0}\neq 0$. This
implies $k^{0}\neq 0,\;\gamma=0,\;y^{0}=0,\;\alpha=\beta=0$, so that this is
an irrelevant case.

(Case6): The conditions $p^{{}^{\prime}2}=m^{2},\;k^{2}=0$, but $q^{2}\neq
m^{2}$ implies $\beta=0,\;y^{0}=0,\;\alpha=0$, and $\gamma k^{0}=0$, so that
$\gamma\neq 0$ and $k^{0}\sim\lambda$, and $|\vec{k}|\sim\lambda$. This is a
relevant case.

(Case7): We have $q^{2}=m^{2},\;k^{2}=0$, and $p^{{}^{\prime}2}\neq m^{2}$.
These are the same conditions to (Case6).

Since three vectors can not be on the mass shell at the same time, and hence
we have exhausted all the cases. Now the conditions lead to

$\displaystyle\alpha=\beta=y^{0}=0,\;\gamma>0,$ (5.3.33) $\displaystyle
k^{2}=0,\;k^{0}\sim|\vec{k}|\sim\lambda,$ (5.3.34)

and (Case4) $p^{\prime}$ and $q$ are both off the mass shell, and (Case6) or
(Case7) One of $p^{\prime}$ and $q$ is on the mass shell, remain as the
possible cases to generate infrared divergences. (Both $p^{\prime}$ and $q$
are off the mass shell is prohibited.) Thus, we have

$\displaystyle\widehat{S.S.}(J)=\left\\{\left(\vec{p^{\prime}},\vec{q},\vec{k};\sqrt{-1}\gamma(\vec{k}d\vec{k})\right)\big{|}\;\gamma>0,\;k^{2}=0,\;|\vec{k}|\sim\lambda,\;\vec{p^{\prime}}=\vec{q}+\vec{k}\right\\}.$
(5.3.35)

Schematically,

$\displaystyle\widehat{S.S}(J)\Leftrightarrow\mathrm{Candidate~{}of~{}infrared~{}contribution~{}to}~{}J.$
(5.3.36)

By the Cutkosky rule, we have

$\displaystyle J+J^{*}\propto J_{a}+J_{b},$ (5.3.37)

and know that the infrared singularities cancel between $J_{a}$ and $J_{b}$,
and that $J$ has no infrared singularities.

The microlocal analysis clarifies the position of singularities, but we need
to know furthermore “Order of Singularities” at the singular position.

Physically, we know the “Ord.” of singularities by power counting with respect
to the photon momentum $k$, that is,

$\displaystyle Ord.[J_{a}]=Ord.[J_{b}]=-2,$ (5.3.38)
$\displaystyle\mathrm{and}~{}~{}~{}Ord.[J]=-1.$ (5.3.39)

These equations show the appearance of infrared divergences and the
cancellation of them.

Next, we will give an “primitive way of power connoting”. First we introduce
the following criterions:

(C1) Assign -2 to the photon propagator,

(C2) Assign -1 to the internal fermion propagator which connects to the photon
on the mass shell and to the outgoing electron on the mass shell,

(C3) Assign +1 to the integration by photon momentum.

With these criterions we have

$\displaystyle Ord.[J_{a}]=Ord.[J_{b}]=-2-1+1=-2,$ (5.3.40) $\displaystyle
Ord.[J_{a}+J_{b}]\leq Ord.[J]=-2+1=-1,$ (5.3.41)

where the Cutkosky rule is taken into account. This reproduces the result
obtained by the explicit integration of the amplitudes.

We hope to have the technique in the micro functions to estimate the Order of
Singularities (Ord.).

(Note added in 2022) Ten years later, a similar work to this Appendix appeared
[51].

—————————————————–

## Acknowledgements and epilogue in 2022

The author deeply thanks Shiro Komata who kindly offered to translate a good
part of this thesis into English, written originally in Japanese in 1978.
Without his help, the English translation would not have been achieved.

He also thanks So Katagiri, Yoshiki Matsuoka, and Kimiko Yamashita for reading
and refining of the draft.

### Epilogue (2022)

Even now, the interplay between the infrared (IR) behavior and the ultraviolet
(UV) behavior in some theories (especially in quantum gravity), seems not to
be fully understood. During the past 45 years, Wilson’s lattice gauge theory
[52] was so successful in doing the strong coupling expansion, which gives not
only hadron masses but also scattering amplitudes consistent to the
experimental results. In this lattice gauge theory, a linearly rising
potential between quarks is given,

$V_{q\bar{q}}(k)\propto\frac{g(|k|)^{2}}{|k|^{2}}\propto\frac{1}{|k|^{4}}+\cdots,~{}\text{or}~{}V_{q\bar{q}}(r)\propto
r+\cdots,$

which was a mere guess in this thesis. Also the evolution equation, similar to
the differential equation derived in this thesis, was developed by Dokshitzer,
Gribov-Lipatov, and Altarelli-Parisi (DGLAP) [53], describing the evolution of
the parton distribution functions, consistently with experiments in QCD.

This thesis claims the importance of Ward-Takahashi identities in the
cancellation of infrared as well as of mass singularities to occur, while
guaranteeing the unitarity. For the cancellation to occur the gauge invariant
set of graphs is essential.

To represent the set, however, Wilson operators, such as W[Path=open] and
W[Path=closed] for non-Abelian gauge theory and W[Path=closed] for gravity,
seem not completely satisfactory. The more suitable definition of variables
surely exists, which is made up of massless quarks, gluons and gravitons,
moving almost in parallel with each other. Such variables include the soft
momentum cutoff $\lambda^{\prime}$ and the angle cutoff $\eta$ from the
parallel direction, or the fictitious mass $\lambda$ introduced to quarks,
gluons and gravitons. The new variables, if denoted as

$X[\text{open};\lambda]~{}\text{and}~{}X[\text{closed};\lambda],$

are expected to follow the string dynamics, with a tension given by the cutoff
$\lambda$,

$\alpha^{\prime}\propto\frac{1}{\lambda^{2}}.$

The theories are IR finite from the beginning, due to the cancellation of IR
divergences for the gauge invariant sets, studied in this thesis, where the
Ward-Takahashi identities were fully used.

Furthermore, if the modular symmetry ($\tau\to-1/\tau$) is imposed on the
partition function, then the UV finiteness can also be achieved. Here, $\tau$
is an evolution parameter in the partition function. The parameter is the
ratio per unit length $L$ of the object, of a time difference at finite
temperature $T$, i.e.

$\tau=\frac{c\Delta t}{L}+i\frac{\hbar c}{k_{B}TL}.$

It is quite plausible that the partition function for an infrared free set or
a linearly extended object can be consistently given when the modular
invariance is satisfied.

Recently, Ward-Takahashi identities associated with the Virasoro-like
symmetry, arising in the soft graviton limit of gravity, has been understood
by A. Strominger and others [54], as the soft graviton theorem. Some 60 years
ago, the symmetry was known by Bondi, van der Burg, Metzner and Sacks (BMS)
[55], and the soft graviton theorem was also known in [2], but the
interrelaton between them is clarified only recently. This new topics called
“Celestial Holography” uses the similar terminology as in this thesis, and
seems to pursuit the similar issue mentioned in the above.

Frankly speaking, many researchers of the infrared slavery at that time, were
surely aimimg at the issue mentioned above. A reason why they (including the
author) would not finish the issue is the appearance of the numerically
calculable lattice gauge theories. History repeats itself.

This motivates the author for translating his old thesis recently. He hopes
this English translation contributes to the study of the issue.

## References

  * [1] S. Weinberg, Rev. Mod. Phys. 46, 255 (1974).
  * [2] F. Bloch and A. Nordsieck, Phys. Rev. 52, 54 (1937);
D. R. Yennie, S. C. Frautschi, and H. Suura, Ann. Phys. (N.Y.) 13, 379 (1961).
The most easy-to-understand reference is S. Weinberg, Phys. Rev. 140, B516
(1965); J. M. Jauch and F. Rohrlich, The Theory of Photons and Electrons, 2nd
edition (Springer, Berlin, 1976).

  * [3] Works in Japan; T. Kinoshita, Prog. Theor. Phys. 5, 1045 (1950); On the cancellation in all orders up to Coulomb phase, N. Nakanishi, ibid. 19, 159 (1958). [This method can be extended to QCD.]
  * [4] Coherent state theory, T. Murota, Prog. Theor. Phys. 24, 1109 (1960);
V. Chung, Phys. Rev. 140, B1110 (1965);
P. P. Kulish and L. D. Faddeev, Theor. Math. Phys. 4, 745 (1970). [We wish to
construct the similar theory in QCD.]

  * [5] T. Appelquist et al., Phys. Rev. Lett. 36, 768 and 1161 (E) (1976);
T. Appelquist and J. Carrazone, Nucl. Phys. B120, 77 (1977). (This relates to
the infrared divergence cancellation of $e^{+}e^{-}\to\text{quark}+\text{X}$
at two loops.)

  * [6] Y. P. Yao, Phys. Rev. Lett. 36, 653 (1976);
Katsumi Matsubara, Soryushiron Kenkyu 53, 134 (1976);
L. Tyburski, Nucl. Phys. B116, 269 (1976). (These relate to quark-quark
scattering at one loop.)

  * [7] A. Sugamoto, Phys. Rev. D16, 1065 (1977). [Reference paper I]
  * [8] F. G. Krausz, Columbia Univ. preprint CU-TP-101 (1977).
  * [9] J. M. Cornwall and G. Tiktopoulos, Phys. Rev. D13, 3370 (1976); UCLA Report No. 76/TEP/22 (1976);
J. M. Cornwall, UCLA Report No. 76/TEP/10 (1976) (unpublished).

  * [10] C. P. Korthals Altes and E. de Rafael, Phys. Rev. Lett. 62B, 320 (1976); Marseille preprint 76/p. 877, December (1976).
  * [11] P. Cvitanović, Phys. Lett. 65B, 272 (1976); Univ. of Oxford preprint, Ref. 8/77.
  * [12] J. Frenkel et al., Phys. Lett. 64B, 211 (1976); Nucl. Phys. B121, 58 (1977).
  * [13] T. Kinoshita and A. Ukawa, Phys. Rev. D15, 1596 (1977); ibid. D16, 332 (1977).
E. Tomboulis, Phys. Lett. 67B, 414 (1977). (This assumes the power counting
theorem of Kinoshita and Ukawa [see [48] below]).

  * [14] H. Miyazawa, Japan-France joint seminar (1977) July, Univ. of Tokyo preprint UT-289;
J. M. Cornwall, UCLA Report No. 77/TEP/11 (1977).

  * [15] A. Sugamoto, H. Yamamoto and N. Nakagawa, Univ. of Tokyo preprint UT-291 (1977). [Reference paper II]
  * [16] In preparation.
  * [17] For non-leading log, J. Frenkel, M. L. Frenkel and J. C. Taylor, Nucl. Phys. B124, 268 (1977); E. C. Poggio, Brandeis Univ. preprint (1977).
  * [18] Use of the Bogoliubov’s renormalization group: V. P. Gusynin and Yu. A. Sitenko, Kiev preprint ITP-77-47E (1977).
  * [19] C. P. Korthals Altes and E. de Rafael, Nucl. Phys. B106, 237 (1976).
  * [20] The first paper of [9]: Phys. Rev. D13, 3370 (1976).
  * [21] F. E. Low, Phys. Rev. 110, 974 (1958).
  * [22] R. L. Arnowitt and S. I. Fickler, Phys. Rev. 127, 1821 (1962);
J. Schwinger, Phys. Rev. 130, 402 (1932);
E. S. Fradkin and I. V. Tyutin, Phys. Rev. D2, 2841 (1970).

  * [23] W. Konetschny and W. Kummer, Nucl. Phys. B100, 106 (1975);
W. Kummer, Acta Physica Austriaca, Suppl. XV, 423 (1976).

  * [24] J. Frenkel, Phys. Rev. D13, 2325 (1976).
  * [25] G. ’t Hooft, Nucl. Phys. B33, 173 (1971).
  * [26] W. Konetschny and W. Kummer, Nucl. Phys. B108, 397 (1976).
  * [27] K. Nishijima, Phys. Rev. 119, 485 (1960); M. Veltman, Physica 29, 186 (1963); G. ’t Hooft and M. Veltman, in Particle Interactions at Very High Energies, Part B, edited by F. Halzen et al. (Plenum, New York, 1974). (The name “diagrammar” appeared.)
  * [28] A similar remark was also made independently by K. J. Kim, Phys. Lett. 69B, 347 (1977).
  * [29] J. Frenkel and R. Meuldermans, Phys. Lett. 65B, 64 (1976). (It seems the calculation needs another check.)
  * [30] Among the divergences produced by the parameter integral, it is the ultraviolet divergence (UV) that is independent on the momentum configuration of the external line, and on the mass of the particle on the internal line. This is because the external momentum configurations and the mass of the internal particles, are irrelevant, when the internal loop momentum becomes infinitely large in the UV limit.
  * [31] Other than the ultraviolet divergence mentioned in [30], another divergence occurs, when the two conditions hold simultaneously; one is the external line being on the mass shell, and the other is a particular particle being massless. This divergence is called infrared divergence (IR), an example of which is the divergence occurring at $\alpha=\beta=0$. By the circuit theory, Feynman parameter acts as a resistance against the current (momentum flow). Therefore $\alpha=\beta=0$, $\gamma=1$ corresponds to the momentum configuration, in which $p^{\prime}+k$ and $p+k$ associated with $\alpha$ and $\beta$, respectively, being hard, and the momentum $k$ with $\gamma$ being soft.
  * [32] This symbol will work as a powerful tool for QCD in Section 4.
  * [33] However, in electron-positron pair creation and pair annihilation processes, two poles at $k^{0}=\frac{\bm{p^{\prime}}\cdot\bm{k}}{p^{\prime 0}}$ and $\frac{\bm{p}\cdot\bm{k}}{p^{0}}$ are split into the upper and lower half planes. Thus one of the two poles contributes to infrared divergence. This contribution acts as pure imaginary for (2.3.8) and produces a pure imaginary phase $e^{i\infty}$ in the scattering amplitude. This is called the “Coulomb phase”. In the cross section the divergence vanishes since $|e^{i\infty}|=1$. For the details, see reference [2] (in particular, the Appendix of S. Weinberg’s paper).
  * [34] Papers particularly discussing the forward scattering are C. T. Sachrajda, Phys. Rev. D14, 1072 (1976); L. Matsson and R. Meuldermans, CERN prepring Ref. TH. 2334-CERN (1977).
  * [35] In particular, for the specific form, see reference [7].
  * [36] The Feynman rule is shown in, for instance, Reference paper I.
  * [37] About the definition of $g(\mu)$, please refer for example to 3.4.23 and 4.4.40–4.4.42.
  * [38] From the experience from QED, $\frac{1}{N-4}\leftrightarrow\frac{1}{4}\ln\frac{1}{\mu^{2}}$ 2.1.25, but the coefficients in the relation are not clear.
  * [39] L. D. Fadeev and V. N. Popov, Phys. Lett. 25B, 29 (1967); Understandable and useful reference is Kiev Report No. ITP 67-36 [NAL-THY-57 (1972) english translation]; Standard trextbook is E. S. Abers and B. W. Lee, Phys. Rep. 9C, 1 (1973).
  * [40] R. P. Feynman, Acta Physica Polonica 24, 697 (1963).
  * [41] This is the so-called Kinoshita Theorem: TK̇inoshita, J. Math. Phys. 3, 650 (1962); TK̇inoshita and A. Ukawa, Phys. Rev. D13, 1573 (1976); The discussions in these papers are done in the Feynman parameter space. The discussion in the momentum space is E. C. Poggio and H. R. Quinn, Phys. Rev. D14, 578 (1976).
  * [42] It is noted that varous choices are possible for the independent momenta for the vertex.
  * [43] It is necessary to pay attention to the position where the arrow symbol is attached. First, fix the flow of independent external momenta and independent internal loop momenta. Then, the momentum flowing a given internal line is a linearly indepemdent sum of these external momenta and internal loop momenta. The arrow symbol associated with a given independent momentum is attached to the lines along which the momentum flow. It is recommended to make some examples for the better understanding.
  * [44] It is noted that, since the external field does not have a color, the first and the second terms in (4.4.36) have the same color factor $[t^{c}]_{ji}$.
  * [45] It is not clear in the off-mass-shell renormalization, however, whether a relation $Z_{1}(\mu)=Z_{2}(\mu)$ holds or not; which depends on how $\mu$ is chosen.
  * [46] The $-[A]_{\text{IR}}$ in Subsection 3.4 corresponds to (4.4.46) of the covariant gauge. However, there is no proper explanation of why $1/(N-4)$ in the covariant gauge becomes $1/(N-4)^{2}$ in the axial gauge. Therefore, it is necessary to recheck the calculation of this $[c^{(1)}(\mu)]_{\text{IR}}$.
  * [47] The application of off-mass-shell or on-mass-shell renormalization for (4.4.2) is the same as the renormalization of (4.4.1) in Subsection 4.4. Express formally (4.4.3) as  $\sum_{n=0}^{\infty}(g_{B})^{2n+2}D^{(n)}=(g_{B}Z_{3}^{1/2})^{2}\left(\frac{p^{\prime}\cdot\varepsilon}{p^{\prime}\cdot k}-\frac{p\cdot\varepsilon}{p\cdot k}\right)\left(\frac{p^{\prime}\cdot\varepsilon}{p^{\prime}\cdot l}-\frac{p\cdot\varepsilon}{p\cdot l}\right)\sum_{n=0}^{\infty}(g_{B})^{2n}B^{(n)},$ and define the renormalized amplitudes of $D^{(n)}$ as $D^{(n)}(\mu)$,  $\sum_{n=0}^{\infty}D^{(n)}=\sum_{n=0}^{\infty}[g(\mu)]^{2n+2}D^{(n)}(\mu),$ then the pertuabative expansion of the renormalized (4.4.3) with $\mu\neq 0$, is given for the off-mass-shell renoramalization by  $D^{(n)(\mu)}=\left(\frac{p^{\prime}\cdot\varepsilon}{p^{\prime}\cdot k}-\frac{p\cdot\varepsilon}{p\cdot k}\right)\left(\frac{p^{\prime}\cdot\varepsilon}{p^{\prime}\cdot l}-\frac{p\cdot\varepsilon}{p\cdot l}\right)\sum_{i,j\geq 0}^{\infty}c^{(i)}(\mu)c^{(j)}(\mu)B^{(n-i-j)},$ while for the on-mass-shell renormalization with $\mu=0$, by  $D^{(n)(0)}=\left(\frac{p^{\prime}\cdot\varepsilon}{p^{\prime}\cdot k}-\frac{p\cdot\varepsilon}{p\cdot k}\right)\left(\frac{p^{\prime}\cdot\varepsilon}{p^{\prime}\cdot l}-\frac{p\cdot\varepsilon}{p\cdot l}\right)B^{(n)}(0).$
  * [48] See the reference [13]. Let us recapitulate the formula by Kinoshita and Ukawa of counting the highest logarithmic powers of the infrared divergences. The power of the logarithms is less than the number of loops. Then, the graph for which the number of the logarithmic powers is equal to the number $L$ of loops [the highest logarithmic power graph], is a planar graph, from each loop of which a single logarithmic power arises. That is, from a graph $G$ with $L$ loops, choose a proper series of subgraphs, $S_{1},S_{2},\cdots,S_{L}=G$, and make reduced graphs $S_{i}/S_{i-1}$ [$i=1,2,\cdots,L$; where $S_{0}$ is a null set]. If each $S_{i}/S_{i-1}$ gives a logarithmic infrared divergence, then the infrared divergence of the original graph $G$ is $\left(\ln\frac{1}{\lambda}\right)^{L}$. [Here, $S_{i}/S_{i-1}$ is obtained from $S_{i}$ by reducing $S_{i-1}$ to a point.] By simply stating, the following factorization holds in the infrared limit:  $[G]_{\text{IR}}\approx[G/S_{L-1}]_{\text{IR}}\times[S_{L-1}/S_{L-2}]_{\text{IR}}\times\cdots\times[S_{2}/S_{1}]_{\text{IR}}$ [Please refer also to P. Cvitanović and T. Kinoshita Phys. Rev. D10, 3991 (1974).] The complexity of their proof comes from the usage of the Feynman parameters. It is interesting to perform the same thing in the momentum space; in this case the proof seems to be extremely simple.
  * [49] More precisely, the integral of (4.5.16) using $\tilde{D}_{\mu\nu}(k;\lambda;\mu)$, remains divergent.
  * [50] Mikio Sato, “Microlocal Analysis and Quantum Field Theory”, the lectures given at University of Tokyo in 1976.
  * [51] T. Kawai and H. P. Stapp, “Microlocal Analysis of Infrared Singularities”, LBL-24163 (1987), published in Algebraic Analysis, I, Academic Press (1988) pp.309-330.
  * [52] K. Wilson, Phys. Rev. D10, 2445 (1975);
About the QCD potential, see G. S. Bali, Phys. Rept. 343, 1 (2000).

  * [53] See for example, A. D. Martin, Acta Physica Polonica B39, 2225 (2008).
  * [54] T. He, V. Lysov, P. Mitra, and A. Strominger, JHEP 05, 151 (2015);
See also a good review, A.-M. Raclorin, arXiv: 2107.02075.

  * [55] H. Bondi, M. G. J. van der Burg, and A. W. K. Metzner, Proc. Roy. Soc. Lond. bf A269, 21 (1962);
R. K. Sachs, Proc. Roy. Soc. Lond. bf A270, 103 (1962); Phys. Rev. 128 2851
(1962).