# An Expansion of Well Tempered Gravity

Eric V. Linder1,2, Stephen Appleby3 1Berkeley Center for Cosmological Physics
& Berkeley Lab, University of California, Berkeley, CA 94720, USA
2Energetic Cosmos Laboratory, Nazarbayev University, Nur-Sultan, 010000,
Kazakhstan
3Asia Pacific Center for Theoretical Physics, Pohang, 37673, Korea

###### Abstract

When faced with two nigh intractable problems in cosmology – how to remove the
original cosmological constant problem and how to parametrize modified gravity
to explain current cosmic acceleration – we can make progress by counterposing
them. The well tempered solution to the cosmological constant through
degenerate scalar field dynamics also relates disparate Horndeski gravity
terms, making them contrapuntal. We derive the connection between the kinetic
term $K$ and braiding term $G_{3}$ for shift symmetric theories (including the
running Planck mass $G_{4}$), extending previous work on monomial or binomial
dependence to polynomials of arbitrary finite degree. We also exhibit an
example for an infinite series expansion. This contrapuntal condition greatly
reduces the number of parameters needed to test modified gravity against
cosmological observations, for these “golden” theories of gravity.

##  I Introduction

Current cosmological acceleration is an overwhelming characteristic of our
universe, driving the expansion rate and shutting down the growth of large
scale structure. Yet in seeking its origin, physics explanations almost
invariably sweep under the rug the elephant in the room: the original
cosmological constant problem that a much higher energy scale vacuum energy
should have dominated the history of the universe, calling into question our
understanding of how gravity reacts to vacuum energy Weinberg:1988cp ;
Carroll:2000fy ; Martin:2012bt ; Nobbenhuis:2004wn ; Padilla:2015aaa .
Exploring a low energy cosmic acceleration by traversing a physics terrain
where there is an elephant under the rug is an uncomfortable position.

Well tempering us2018 ; nott2018 ; us2020a ; us2020b  aims to solve the
original cosmological constant problem by employing a dynamical scalar field
with certain degeneracy conditions in the equations of motion. To do so
requires modifications of general relativity, which is also one of the favored
approaches to explaining current cosmic acceleration. However, modified
gravity is quite difficult to test against cosmological observations in a
general manner, without a large number of model dependent assumptions and
multiple free functions. By contraposing the two problems, well tempering
shows a path toward resolving them both.

Modified gravity often works within Horndeski theory, the most general scalar-
tensor theory giving second order equations of motion. While generally this
involves four free functions, we take the simplest approach to setting the
speed of gravitational waves to be equal to the speed of light so that
$G_{5}=0$ and $G_{4}=G_{4}(\phi)$. Then Horndeski gravity involves three free
functions $K(\phi,X)$, $G_{3}(\phi,X)$, $G_{4}(\phi)$ of two variables: the
scalar field $\phi$ and its canonical kinetic form
$X=-g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi/2=\dot{\phi}^{2}/2$ for a
homogeneous scalar field in a Robertson-Walker spacetime.

To test against cosmological data, one could adopt the effective field theory
approach gubitosi ; eft1 ; glpv ; bellsaw ; eft2  that reduces three functions
of two variables to four functions of time, at the linear perturbation level.
This has the drawback that one loses information from the nonlinear regime
(including details of screening to satisfy solar system tests), where
cosmological data can give great insight. If we wish to keep the full
leverage, one generally has to assume a specific gravity theory, a specific
functional dependence within that theory, and specific parameters within those
functions. For example, one does not “test gravity” but says the data is or is
not consistent with, say, “the class of $f(R)$ gravity with a specific
function $f$ and specific ranges of parameters appearing within that function
$f$.” The conclusions tend to be highly model dependent. To put the general
case of three functions of two variables in perspective, figuring out how to
modify gravity yet obtain consistency with cosmological observations is like
trying to figure out how to kick a football to score a goal when 1) the pull
of gravity depends on the ball’s location and velocity, 2) the wind speed
varies with the ball’s location and velocity, and 3) the ball’s mass depends
on its location (e.g. $K(\phi,X)$, $G_{3}(\phi,X)$, $G_{4}(\phi)$
respectively).

In us2020a  we demonstrated that well tempering could not only address the
original cosmological constant problem, but reduce the modified gravity
parameter space from multiple functions to a small number ($\sim$ four) of
parameters. These give well defined, valuable theories that are predictive
across the full range of cosmological observations, including on nonlinear
scales, and could be regarded as the “golden” theories to use in cosmology,
without climbing over the elephant under the rug.

The solutions in us2020a  were predominantly within shift symmetric theories,
which have important protection against quantum corrections Finelli:2018upr ;
Heisenberg:2018vsk , and gave particular solutions to the well tempered
degeneracy equation by assuming various Ansätze for either $K(X)$ or
$G_{3}(X)$, with the other then determined by the equations. The degeneracy
equation is a nonlinear differential equation and so other solutions can exist
as well. Here we carry out a systematic analysis of all solutions having
certain conditions discussed below.

In Section II we restate the key degeneracy equation giving well tempering and
introduce a series expansion method. The remarkably simple solutions for a
finite series are presented in Sec. III and we discuss some special cases as
well, and how previous solutions are unified by our general expression. We
present an example of an infinite series solution in Sec. IV. We summarize and
conclude in Sec. V.

##  II Series Expansion Approach

Well tempering works by reducing the equations of motion involving $\ddot{a}$
and $\ddot{\phi}$, i.e. the expansion and scalar field evolutions, to have
degenerate solutions for de Sitter spacetime, referred to as “on shell”. This
empowers the scalar field dynamics to cancel the cosmological constant both on
and off shell, i.e. for the full cosmological history. See us2018  for
details.

The equation guaranteeing the degeneracy is

$\displaystyle(M-2g)\left\\{3h\dot{\phi}K_{X}+18h^{2}g-6h^{2}M+\lambda^{3}\right\\}$
$\displaystyle\qquad=(K_{X}+2XK_{XX}+6h\dot{\phi}g_{X})\left[-h\dot{\phi}(M-6g)+2XK_{X}\right]\,,$
(1)

where we have imposed shift symmetry so that $G_{4}=(M^{2}_{\rm pl}+M\phi)/2$,
and $G_{3}=G_{3}(X)$, $K(\phi,X)=F(X)-\lambda^{3}\phi$. For notational
convenience we write $g=XG_{3X}$, where a subscript $X$ denotes a derivative
with respect to $X$, and $h\equiv H_{\rm dS}$ is the de Sitter value of the
expansion rate.

This nonlinear differential equation relates the Horndeski terms $K$, $G_{3}$,
and $G_{4}$. A variety of solutions were presented and analyzed in us2020a .
Here we look for more general solutions using a series expansion method like
that introduced and used to good effect in nott2018 . Note that in their work
they assumed $K=X-V(\phi)$, which turns the nonlinear equation into a simpler
one (but their nonassumption of shift symmetry makes it a partial differential
equation). Here we impose shift symmetry, justified for its protection against
quantum loop corrections, but allow general $K(X)$ with the permitted
$\lambda^{3}\phi$ tadpole term.

By resummation and judicious insight, a series expansion can provide the
general functional dependence $K(X)$ or $G_{3}(X)$. This method can deliver an
algorithm for generating new solutions, and some general principles.

Since the quantities entering the degeneracy equation are $K_{X}$ and $g$, we
will carry out series expansion of these; they can be rewritten in terms of
$K$ and $G_{3}(X)$ as desired. Thus,

$\displaystyle g$ $\displaystyle=$ $\displaystyle\sum a_{n}X^{n/2}$ (2)
$\displaystyle K_{X}$ $\displaystyle=$ $\displaystyle\sum b_{p}X^{p/2}\,.$ (3)

The sums go from $n_{\rm min}$ to $n_{\rm max}$ and $p_{\rm min}$ to $p_{\rm
max}$ respectively, and we will see there are relations between these values.
The degeneracy equation becomes

$\displaystyle 0$ $\displaystyle=$ $\displaystyle 6h^{2}M\sum
a_{n}\left(5+n-\frac{\lambda^{3}}{3h^{2}M}\right)\,X^{n/2}-36h^{2}\left(\sum
a_{n}X^{n/2}\right)\,\left(\sum a_{n}(1+n)X^{n/2}\right)$ (4) $\displaystyle+$
$\displaystyle Mh\sqrt{2X}\sum b_{p}(4+p)X^{p/2}-2X\left(\sum
b_{p}X^{p/2}\right)\,\left(\sum b_{p}(1+p)X^{p/2}\right)$ $\displaystyle-$
$\displaystyle 6h\sqrt{2X}\left[\left(\sum a_{n}X^{n/2}\right)\,\left(\sum
b_{p}(1+p)X^{p/2}\right)+\left(\sum a_{n}(1+n)X^{n/2}\right)\,\left(\sum
b_{p}X^{p/2}\right)\right]$ $\displaystyle-$ $\displaystyle
6h^{2}M^{2}+M\lambda^{3}\,.$

##  III Finite Series Solution

One solves Eq. (4) by equating terms order by order in $X$ to obtain relations
between $a_{n}$ and $b_{p}$, and then summing the series to derive functional
relations between $K$ and $G_{3}$ (or $g$), with or without $M$ from the
$G_{4}$ term. We begin by considering a finite series.

###  III.1 General Results

The results are

$\displaystyle K_{X}$ $\displaystyle=$
$\displaystyle-3h\sqrt{2}X^{-1/2}g+h\sqrt{2}X^{-1/2}\left(M-\frac{\lambda^{3}}{6h^{2}}\right)\qquad{\rm[Branch\
A]}$ (5) $\displaystyle K_{X}$ $\displaystyle=$
$\displaystyle-3h\sqrt{2}X^{-1/2}\left(g-a_{0}\right)-\frac{3h\sqrt{2}}{2}\,X^{-1/2}\int
dX\,\frac{g-a_{0}}{X}+b_{-1}X^{-1/2}\qquad{\rm[Branch\ B]}\,.$ (6)

Appendix A shows the steps in the derivation but one can verify the solutions
by direct substitution into Eq. (1).

The coefficient $b_{-1}$ entering as $K_{X}^{[-1]}=b_{-1}X^{-1/2}$ plays a
special role (since $K_{X}^{[-1]}+2XK_{XX}^{[-1]}=0$ in the degeneracy
equation) and helps distinguish the two branches. For Branch A,

$b_{-1}=h\sqrt{2}M-\frac{\lambda^{3}}{3h\sqrt{2}}-3h\sqrt{2}\,a_{0}\qquad[{\rm
Branch\ A}]\,,$ (7)

i.e. $b_{-1}$ is specifically connected to $a_{0}$, the constant term in $g$
or the contribution $G_{3}^{[0]}=a_{0}\ln X$ (recall from us2020a  that
constant $g$ is a hallmark of Brans-Dicke type scalar-tensor theories as well
as $f(R)$ gravity and No Slip Gravity). The quantity $b_{-1}$ is also
connected to the $G_{4}$ mass $M$ and the tadpole scale $\lambda$ for Branch
A, so all the functions are woven together.

Branch B has different contrapuntal conditions between the functions, with

$b_{-1}=\\{{-h\sqrt{2}M,\ \rm arbitrary}\\},\quad
a_{0}=\frac{M}{2},\quad\lambda^{3}=3h^{2}M\qquad[{\rm Branch\ B}]\,.$ (8)

Here the connections are more concrete between $G_{3}$, $G_{4}$, and $K$, i.e.
$a_{0}$, $M$, and $\lambda$, but $b_{-1}$ also has the possibility of being
arbitrary, as we discuss in the next section.

Equations (5) and (6) are general solutions. They cover many of the particular
solutions given in us2020a  (hereafter called HV). For example, HV Eq. 3.36 is
a Branch A solution; Eqs. 3.34 and 3.38 are Branch B solutions. We have
essentially succeeded in deriving a unified solution for all the individual
solutions obtained in HV.

Some particular solutions need or merit special treatment and we deal with
these in the next subsections.

###  III.2 Special Case: Truncation Below $b_{-1}$

Under the standard solutions of the previous section, for some $a_{n}$, even
$a_{n_{\rm min}}$, there must be a $b_{n-1}$, as shown in Eq. (26). However,
we see that the $a^{2}_{n_{\rm min}}$ contribution (the second term in Eq. 4)
vanishes when $n=-1$ and so this term is an exception: if we set $n_{\rm
min}=-1$ we are free to set $b_{-2}=0$ (and all $b_{p<-2}$ are zero likewise).
This may be considered attractive in that it means that $K$ does not have any
negative powers of $X$, so when the field rolls slowly there is no blow up of
the kinetic term. Thus we consider this special case $n_{\rm min}=-1$, i.e.
$a_{n<-1}=0$.

With $b_{-2}=0$, for the next higher order equation, which involves $X^{-1/2}$
terms, $b_{-1}$ actually cancels out of the degeneracy equation. We continue
going through the intermediate powers and find that for Branch A, $b_{-1}$
takes the form in Eq. (7), unless $a_{n>0}=0$ and $a_{0}=M/2$, in which case
$b_{-1}$ remains arbitrary. The value of $a_{-1}$ does not affect $b_{-1}$
under these conditions. For Branch B, Eq. (7) reduces to $b_{-1}=-h\sqrt{2}M$;
however, if the above conditions hold and we also require $a_{-1}=0$, then
$b_{-1}$ is arbitrary. This explains the “$b_{-1}$ arbitrary” possibility
mentioned in the previous subsection.

Cases with arbitrary $b_{-1}$ can be seen in HV Eqs. 3.31, 3.35, and 3.39, all
Branch B solutions, and Eq. 3.37, a Branch A solution but where the arbitrary
$b_{-1}$ arises from an arbitrary term in $a_{0}$, i.e. it still follows Eq.
(7).

###  III.3 Special Case: $n_{\rm min}=0$, $n_{\rm max}=1$

One other special case of note is truncation where $n_{\rm min}=0$ and in
addition $n_{\rm max}=1$, i.e. $g=sX^{1/2}+a_{0}$. This gives only three
equations, so in addition to determining $b_{-1}$ and $b_{0}$ they must fix
$a_{0}$. The result shows that despite the tight restrictions this case
nevertheless follows the general Eqs. (5) and (6). That is, Branch A gives

$\displaystyle g$ $\displaystyle=$ $\displaystyle
sX^{1/2}+\frac{M}{3}-\frac{b_{-1}}{3h\sqrt{2}}-\frac{\lambda^{3}}{18h^{2}}$
(9) $\displaystyle K_{X}$ $\displaystyle=$
$\displaystyle-3h\sqrt{2}s+b_{-1}X^{-1/2}\,,$ (10)

with $b_{-1}$ arbitrary and $\lambda^{3}\neq 3h^{2}M$, following from Eq. (5),
and Branch B, i.e. Eq. (6), yields

$\displaystyle g$ $\displaystyle=$ $\displaystyle
sX^{1/2}+\frac{M}{3}+\frac{\lambda^{3}}{18h^{2}}$ (11) $\displaystyle K_{X}$
$\displaystyle=$
$\displaystyle-6h\sqrt{2}s-\frac{\lambda^{3}\sqrt{2}}{3h}X^{-1/2}\,.$ (12)

We can show this as follows. The $X^{1}$ equation gives
$b_{0}=\\{-3,-6\\}h\sqrt{2}\,a_{1}$, as expected from Eq. (25); for the choice
$b_{0}=-3h\sqrt{2}\,a_{1}$, i.e. following Branch A, the $X^{1/2}$ equation
gives Eq. (7) for $b_{-1}$, but for the choice $b_{0}=-6h\sqrt{2}\,a_{1}$ then
$b_{-1}$ cancels out and instead the equation imposes

$a_{0}=\frac{M}{3}+\frac{\lambda^{3}}{18h^{2}}\,,$ (13)

so this is a hidden version of Branch B. The hidden aspect arises because of
the lack of extra equations that would impose the usual additional consistency
conditions such that $\lambda^{3}=3h^{2}M$ and hence $a_{0}=M/2$. Finally, the
$X^{0}$ equation determines $b_{-1}$ following Eq. (7), which would only lead
to Eq. (8) upon setting the additional consistency conditions $a_{0}=M/2$,
$\lambda^{3}=3h^{2}M$, which do not apply here. Thus we have
$g=sX^{1/2}+a_{0}$ leading to HV Eq. 3.37 for Branch A and HV Eq. 3.33 for the
“hidden” Branch B.

###  III.4 Special Case: $\ln X$ Terms in $K_{X}$

A finite expansion in powers will not always work. Before we move on to an
infinite series in Sec. IV, let us consider the special case where $K_{X}$
contains terms involving $\ln X$. This is particularly of note since us2020a
did find solutions involving $\ln X$. Equations (21) and (22) below provide
the final solutions. To derive them we begin by expanding $g$ and studying the
form of the resulting degeneracy equation, writing

$\displaystyle g$ $\displaystyle=$ $\displaystyle\sum a_{n}X^{n/2}$ (14)
$\displaystyle K_{X}$ $\displaystyle=$ $\displaystyle X^{-1/2}\,B(X)\,,$ (15)

where $K$ still has a tadpole term $-\lambda^{3}\phi$.

The degeneracy equation becomes

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
XB_{X}\left(4B+12h\sqrt{2}\sum
a_{n}X^{n/2}-2h\sqrt{2}M\right)+h\sqrt{2}B\left[6\sum
a_{n}(1+n)X^{n/2}-3M\right]$ $\displaystyle-6h^{2}M\sum
a_{n}\left(5+n-\frac{\lambda^{3}}{3h^{2}M}\right)+36h^{2}\left(\sum
a_{n}X^{n/2}\right)\left(\sum
a_{n}(1+n)X^{n/2}\right)+6h^{2}M^{2}-M\lambda^{3}\,.$

A term involving $B\sim X^{m}\ln X$ has nothing to cancel against if $m\neq
0$, so we allow only a term like $B\sim\ln X$, besides standard powers.
Writing

$B=b_{r}\,\ln X+\sum b_{p}X^{p/2}\,,$ (17)

the degeneracy equation for the terms involving $\ln X$ becomes

$0=4b_{r}\ln X\,\left(b_{r}+\frac{1}{2}\sum
b_{p}pX^{p/2}\right)+h\sqrt{2}b_{r}\ln X\,\left[6\sum
a_{n}(1+n)X^{n/2}-3M\right]\,.$ (18)

If we set $b_{r}=0$, there is no $\ln X$ term and we return to the power
series of the previous sections. The solution for general $p\neq 0$, $-1$ is

$a_{p}=b_{p}\,\frac{p}{3h\sqrt{2}(1+p)}\,.$ (19)

However, looking at the terms in the degeneracy equation involving only powers
and not $\ln X$, i.e. $X^{p}$, we find the only consistent solution is
$b_{p}=0$ for all $p\neq 0$. For $p=0$, we find $b_{0}$ is arbitrary and

$a_{0}=\frac{M}{2}-\frac{b_{r}\sqrt{2}}{3h}\,.$ (20)

For the $X^{0}$ order (without $\ln X$), the solution requires either
$a_{0}=M/2$ (implying $b_{r}=0$ and hence reducing to the pure power expansion
without $\ln X$ as in the previous sections) or $\lambda^{3}=3h^{2}M$.

Thus the final solution is $B=b_{0}+b_{r}\ln X$, or

$K_{X}=b_{0}X^{-1/2}+b_{r}X^{-1/2}\ln X\,,$ (21)

with a tadpole term $-3h^{2}M\phi$ in $K$, and

$g=\frac{M}{2}-\frac{b_{r}\sqrt{2}}{3h}+a_{-1}X^{-1/2}\,,$ (22)

where $a_{-1}$ is arbitrary. This is the unique solution where $K_{X}$
involves a $\ln X$ term, and is equivalent to HV Eq. 3.20. When $g=0$, so
$b_{r}=3hM/(2\sqrt{2})$, we reproduce HV Eq. 3.11, and when $g=rM$ (as for
$f(R)$ and No Slip Gravity), we obtain HV Eq. 3.19.

##  IV Infinite Series Solution

When the series expansion is infinite then one must use the Cauchy product to
determine terms at a certain order. Since there is no finite $n_{\rm max}=N$,
we are no longer constrained by the lack of a counter term to
$a_{N}^{2}X^{N}$. This breaks the relation between $a_{n}$ and $b_{n-1}$.
There is very little one can say in general for such a situation. However, if
one restricts the series in some way then some progress can be made. For
example, consider the case where $b_{n\neq-1}=0$. Then the only series product
comes from the second term in Eq. (4), which we evaluate using the Cauchy
product, except for separating out where one term has $a_{-1}X^{-1/2}$. If we
choose $a_{n}$ as a semi-infinite series, with $a_{n<-1}=0$, then we can solve
the equation.

Starting with the $X^{-1/2}$ power equation and working up, we can obtain all
$a_{n}$ and hence $g$. A particularly compact form obtains for
$b_{-1}=-\lambda^{3}\sqrt{2}/(3h)$, in that then the (semi)infinite series
sums to the solution

$g=\frac{M}{2}+\frac{c}{2}\,X^{-1/2}\left[1\pm\sqrt{1+\frac{2(3h^{2}M-\lambda^{3})}{9h^{2}c}X^{1/2}}\
\right]\,.$ (23)

Note that the lower root gives $a_{-1}=0$, and hence $g$ involves only
nonnegative powers of $X$. (For the upper root $c=a_{-1}$.) When
$\lambda^{3}=0$, the entire kinetic term $K=2b_{-1}X^{1/2}-\lambda^{3}\phi=0$
and we have HV Eq. 3.45. When $\lambda^{3}>3h^{2}M$ we have a “speed limit” on
$X$, i.e. the scalar field motion $\dot{\phi}$, to keep the function real (see
for example eva1 ; eva2 ; eva3 ).

##  V Conclusions

The reduction of the description of modified gravity from three functions of
two variables to a single function of one variable, or a handful of constant
parameters, could open up powerful leverage on scanning theory space to
compare to observational data. Remarkably, we have shown this can be done with
the bonus of solving the original cosmological constant problem – through well
tempering – and protecting from at least some quantum corrections – through
shift symmetry. By contraposing two highly challenging problems we solve both.

Equations (5) and (6) give general solutions relating the Lagrangian terms
$K$, $G_{3}$, and $G_{4}$ for a wide range of gravity theories, derived using
a power series expansion under well tempering. We show that they unify the
disparate solutions found piecemeal in us2020a , while going well beyond them,
extending monomial or binomial cases to arbitrary finite polynomials.
Logarithmic terms are included as well. Branch A solutions are what us2020a
referred to as the $(\checkmark)$ class that gives a scalar field equation
that becomes trivial on shell, while Branch B solutions are fully well
tempering.

Allowing the series expansion to be infinite gives formal solutions but ones
difficult to sum to a compact functional form. We exhibit one example where
this can be done, generalizing a case from us2020a . For this result, the full
gravity theory can be described by four constant parameters: $M$, $\lambda$,
$h$, and $c$, which can readily be sampled for likelihood estimation compared
to data.

Such theories that possess highly desirable characteristics for fundamental
physics – solving rather than neglecting the cosmological constant problem,
and protecting against quantum corrections – soundness, and robust ability for
full comparison with observations (including nonlinear scales, in principle)
could be regarded as the favored “golden” gravity theories to work with. While
exciting work remains to investigate further their detailed properties, they
represent a significant step away from arbitrary functions toward true
benchmarks.

###### Acknowledgements.

SAA is supported by an appointment to the JRG Program at the APCTP through the
Science and Technology Promotion Fund and Lottery Fund of the Korean
Government, and was also supported by the Korean Local Governments in
Gyeongsangbuk-do Province and Pohang City. EL is supported in part by the
Energetic Cosmos Laboratory and by the U.S. Department of Energy, Office of
Science, Office of High Energy Physics under contract no. DE-AC02-05CH11231.

##  Appendix A Derivation for Finite Series

To begin the analysis of the finite series case order by order, let us look at
the maximum powers for each series. For $N\equiv n_{\rm max}$ we have one term
going as $X^{N}$. To match this we must have $P\equiv p_{\rm max}=N-1$, unless
$N=0$. That is, the $K_{X}$ series must cut off at one less power than the $g$
series. We find

$b_{N-1}=-3h\sqrt{2}\,a_{N}\left(1+\frac{1\mp
1}{2N}\right)=\left\\{-3h\sqrt{2}\,a_{N},\,-3h\sqrt{2}\,a_{N}\frac{N+1}{N}\right\\}\,.$
(24)

Thus there are two branches of solutions, what we call Branch A and Branch B.
We can then proceed to the next lowest power, $X^{N-1/2}$ and determine
$b_{N-2}$ from $a_{N-1}$. By continuing this process for $X^{n}$, for
$N/2<n\leq N$, we find

$b_{p}=\left\\{-3h\sqrt{2}\,a_{p+1},\,-3h\sqrt{2}\,a_{p+1}\,\frac{p+2}{p+1}\right\\}\qquad{\rm
for}\qquad 0\leq p\leq N-1\,.$ (25)

There are $N$ equations, for $X^{N}$, $X^{N-1/2}$,…$X^{(N+1)/2}$, and these
define the $N$ coefficients $b_{N-1}$, $b_{N-2}$,…$b_{0}$. Now let us jump to
the most negative power of $X$, i.e. $L\equiv n_{\rm min}<0$ and evaluate the
$b_{p}$ with $p<0$ working upward. One obtains basically the same expression
as Eq. (25), and again $p_{\rm min}=n_{\rm min}-1$:

$b_{p}=\left\\{-3h\sqrt{2}\,a_{p+1},\,-3h\sqrt{2}\,a_{p+1}\,\frac{p+2}{p+1}\right\\}\qquad{\rm
for}\qquad L-1\leq p<-1\,.$ (26)

Again there are two branches, and we must choose the same branch for the
negative powers as the positive powers. These determine the $L$ coefficients
$b_{L-1}$, $b_{L}$,…$b_{-2}$.

There are several equations for the powers $X^{n}$ with $-L/2\leq n\leq N/2$
and only a single parameter $b_{-1}$ left to determine. All these equations
must give a consistent solution for $b_{-1}$. Note that more and more terms
from Eq. (4) enter into these equations, but a consistent solution occurs for
each branch. For branch A, we need

$b_{-1}=h\sqrt{2}M-\frac{\lambda^{3}}{3h\sqrt{2}}-3h\sqrt{2}\,a_{0}\qquad[{\rm
Branch\ A}]\,,$ (27)

and for branch B, the requirements are

$b_{-1}=\\{{-h\sqrt{2}M,\ \rm arbitrary}\\},\quad
a_{0}=\frac{M}{2},\quad\lambda^{3}=3h^{2}M\qquad[{\rm Branch\ B}]\,.$ (28)

To summarize, all $b_{p}$ are determined by Eq. (25) (or equivalently Eq. 26),
except for $b_{-1}$ – special since for this order $K_{X}+2XK_{XX}=0$ – which
is given by Eq. (27) or (28) for the respective branch.

We can now try to sum up the series to obtain a functional relation. For
Branch A,

$\displaystyle K_{X}$ $\displaystyle=$ $\displaystyle\sum
b_{p}X^{p/2}=-3h\sqrt{2}\sum_{p\neq-1}a_{p+1}X^{p/2}+b_{-1}X^{-1/2}=-3h\sqrt{2}X^{-1/2}\sum_{p\neq-1}a_{p+1}X^{(p+1)/2}+b_{-1}X^{-1/2}$
(29) $\displaystyle=$ $\displaystyle-3h\sqrt{2}X^{-1/2}\sum_{n\neq
0}a_{n}X^{n/2}+b_{-1}X^{-1/2}$ $\displaystyle=$
$\displaystyle-3h\sqrt{2}X^{-1/2}g+h\sqrt{2}X^{-1/2}\left(M-\frac{\lambda^{3}}{6h^{2}}\right)\qquad{\rm[Branch\
A]}\,.$

For Branch B,

$\displaystyle K_{X}$ $\displaystyle=$
$\displaystyle-3h\sqrt{2}\sum_{p\neq-1}\frac{p+2}{p+1}\,a_{p+1}X^{p/2}+b_{-1}X^{-1/2}=-3h\sqrt{2}X^{-1/2}\sum_{p\neq-1}\frac{p+2}{p+1}\,a_{p+1}X^{(p+1)/2}+b_{-1}X^{-1/2}$
(30) $\displaystyle=$ $\displaystyle-3h\sqrt{2}X^{-1/2}\sum_{n\neq
0}\left(1+\frac{1}{n}\right)\,a_{n}X^{n/2}+b_{-1}X^{-1/2}$ $\displaystyle=$
$\displaystyle-3h\sqrt{2}X^{-1/2}\left(g-a_{0}\right)-\frac{3h\sqrt{2}}{2}\,X^{-1/2}\int
dX\,\frac{g-a_{0}}{X}+b_{-1}X^{-1/2}\qquad{\rm[Branch\ B]}\,.$

Recall that for Branch B, $\lambda^{3}=3h^{2}M$ and the constant part of $g$
is simply $a_{0}=M/2$. Thus we have derived our general solutions Eqs. (5) and
(6).

Substituting these back into the degeneracy Eq. (1), we find that Branch A
gives a zero for the first factor on the right hand side,
$K_{X}+2XK_{XX}+6h\dot{\phi}g_{X}$. This indicates the coefficient of
$\ddot{\phi}$ in the scalar field equation vanishes on shell. This is what we
called $(\checkmark)$ models in us2020a . Branch B does not zero out
coefficients of $\ddot{\phi}$ and gets a full $\checkmark$.

Note that if desired we can set $M=0$, removing the coupling to the Ricci
scalar and making the $G_{4}$ term standard. For Branch B this will also make
$\lambda^{3}=0$ and $a_{0}=0$.

## References

  * (1) S. Weinberg, The Cosmological Constant Problem, Rev. Mod. Phys. 61, 1 (1989)
  * (2) S. M. Carroll, The Cosmological constant, Living Rev. Rel. 4, 1 (2001) [arXiv:astro-ph/0004075]
  * (3) J. Martin, Everything You Always Wanted To Know About The Cosmological Constant Problem (But Were Afraid To Ask), Comptes Rendus Physique 13, 566 (2012) [arXiv:1205.3365]
  * (4) S. Nobbenhuis, Categorizing different approaches to the cosmological constant problem, Found. Phys. 36, 613 (2006) [arXiv:gr-qc/0411093]
  * (5) A. Padilla, Lectures on the Cosmological Constant Problem, arXiv:1502.05296
  * (6) S. Appleby, E.V. Linder, The Well-Tempered Cosmological Constant, JCAP 1807, 034 (2018) [arXiv:1805.00470]
  * (7) W.T. Emond, C. Li, P.M. Saffin, S-Y. Zhou, Well-Tempered Cosmology, JCAP 1905, 038 (2020) [arXiv:1812.05480]
  * (8) S. Appleby, E.V. Linder, The Well-Tempered Cosmological Constant: The Horndeski Variations, JCAP in press [arXiv:2009.01720]
  * (9) S. Appleby, E.V. Linder, The Well-Tempered Cosmological Constant: Fugue in B♭, JCAP in press [arXiv:2009.01723]
  * (10) G. Gubitosi, F. Piazza, and F. Vernizzi, The effective field theory of dark energy, JCAP 1302, 032 (2013) [arXiv:1210.0201]
  * (11) J.K. Bloomfield, E.E. Flanagan, M. Park, S. Watson, Dark energy or modified gravity? An effective field theory approach, JCAP 1308, 010 (2013) [arXiv:1211.7054]
  * (12) J. Gleyzes, D. Langlois, F. Piazza, F. Vernizzi, Essential building blocks of dark energy, JCAP 1308, 025 (2013) [arXiv:1304.4840]
  * (13) E. Bellini, I. Sawicki, Maximal freedom at minimum cost: linear large-scale structure in general modifications of gravity, JCAP 1407, 050 (2014) [arXiv:1404.3713]
  * (14) E.V. Linder, G. Sengör, S. Watson, Is the Effective Field Theory of Dark Energy Effective?, JCAP 1605, 053 (2016) [arXiv:1512.06180]
  * (15) B. Finelli, G. Goon, E. Pajer and L. Santoni, The Effective Theory of Shift-Symmetric Cosmologies, JCAP 1805, 060 (2018) [arXiv:1802.01580]
  * (16) L. Heisenberg, A systematic approach to generalisations of General Relativity and their cosmological implications, Phys. Rept. 796, 1 (2019) [arXiv:1807.01725]
  * (17) E. Silverstein, D. Tong, Scalar speed limits and cosmology: Acceleration from D-cceleration, Phys. Rev. D 70, 103505 (2004) [arXiv:hep-th/0310221]
  * (18) M. Alishahiha, E. Silverstein, D. Tong, DBI in the sky, Phys. Rev. D 70, 123505 (2004) [arXiv:hep-th/0404084]
  * (19) D. Mathis, A. Mousatov, G. Panagopoulos, E. Silverstein, A new branch of inflationary speed limits, arXiv:2010.00113


# Desingularization and p-Curvature of Recurrence Operators

Yi Zhou  Florida State UniversityTallahasseeUSA<EMAIL_ADDRESS>and  Mark
van Hoeij  Florida State UniversityTallahasseeUSA<EMAIL_ADDRESS>
(2022)

###### Abstract.

Linear recurrence operators in characteristic $p$ are classified by their
$p$-curvature. For a recurrence operator $L$, denote by $\chi(L)$ the
characteristic polynomial of its $p$-curvature. We can obtain information
about the factorization of $L$ by factoring $\chi(L)$. The main theorem of
this paper gives an unexpected relation between $\chi(L)$ and the true
singularities of $L$. An application is to speed up a fast algorithm for
computing $\chi(L)$ by desingularizing $L$ first. Another contribution of this
paper is faster desingularization.

linear recurrence equations, singularities, $p$-curvature, algorithms

††copyright: acmcopyright††journalyear: 2022††doi:
XXXXXXX.XXXXXXX††conference: International Symposium on Symbolic and Algebraic
Computation; July 04–07, 2022; Lille, France††price: 15.00††isbn:
978-1-4503-XXXX-X/18/06††ccs: Computing methodologies Algebraic algorithms

##  1\. Introduction

Singularities of linear difference operators can be divided into two groups,
true (i.e. non-removable) singularities, and apparant (i.e. removable)
singularities. Desingularization (detecting or removing apparant
singularities) can expedite various algorithms for difference or differential
equations. An early application (HomURL, ) appeared in
`DEtools[Homomorphisms]` in Maple 10. Other algorithms that benefit from
reducing the number of singularities include finding closed form solutions and
factoring, e.g. `LREtools[RightFactors]` in Maple 2021.

In characteristic $p$, difference operators can be classified by the so-called
$p$-curvature. Our main result gives a relation between $\chi(L)$, the
characteristic polynomial of the $p$-curvature of $L$, and the true
singularities of $L$. We prove that the denominator of $\chi(L)$ determines
the true singularities, including their multiplicities, up to shift
equivalence.

The algorithm from (Bostan:2014:FAC:2608628.2608650, ) computes $\chi(L)$,
multiplied by a denominator bound, by computing its $Z$-adic expansion. One
application of our theorem is that we can replace the denominator bound by the
exact denominator. This lowers the required $Z$-adic precision, which can
speed up the computation, see subsection 5.2.

We want desingularization to take less time than the time it saves in
applications. Then it is useful to compute a partial desingularization (where
the goal is to remove most apparant singularities, at a fraction of the cost
of a full desingularization). We give various algorithms for this in section
6.

Both authors were supported by NSF grants 1618657 and 2007959.

##  2\. Preliminaries

###  2.1. Desingularization

Let $F$ be a field. Let $P=F[x][y]$ and $D=F(x)[y]$. If $f\in P$, then $f$ is
called primitive if the gcd of its coefficients in $F[x]$ is 1. If $f\in
D-\\{0\\}$, then there is $c\in F(x)-\\{0\\}$, unique up to a factor in $F$,
for which $c^{-1}f\in P$ is primitive. The content of $f$, denoted ${\rm
Cont}(f)$, is this $c$, while the primitive part of $f$ is ${\rm
Prim}(f)=c^{-1}f\in P$. A version of Gauss’s lemma says ${\rm
Cont}(f_{1}f_{2})={\rm Cont}(f_{1}){\rm Cont}(f_{2})$.

Let $\tau$ be the shift-operator. If $r(x)$ is a rational function then the
product $\tau\cdot r(x)$ equals $r(x+1)\cdot\tau$. This product turns
$\mathcal{P}:=F[x][\tau]$ and $\mathcal{D}:=F(x)[\tau]$ into non-commutative
rings. The product corresponds to compositions of operators, where
$L=\sum_{i=0}^{n}a_{i}(x)\tau^{i}\in\mathcal{D}$ operates on $y(x)$ as
$L(y(x))=\sum_{i=0}^{n}a_{i}(x)y(x+i)$.

If $L\in\mathcal{D}-\\{0\\}$ we can define ${\rm Cont}(L)\in F(x)$ and ${\rm
Prim}(L)\in\mathcal{P}$ in the same way as before.

###### Definition 2.1.

An operator $L\in\mathcal{P}$ is called Gaussian if

$\forall_{A\in\mathcal{D}}\ AL\in\mathcal{P}\Longrightarrow A\in\mathcal{P}.$

If $L$ is Gaussian then $L$ is primitive. In the commutative case, the two
properties are equivalent by Gauss’s lemma. It is known that Gauss’s lemma
does not hold in the non-commutative case, which is illustrated in Example 2.3
below (see also (Jaroschek, , p. 27)).

An element $L=\sum_{i=0}^{n}a_{i}(x)\tau^{i}\in\mathcal{P}$ corresponds to a
recurrence relation $L(y(x))=0$, i.e.

(1) $a_{n}(x)y(x+n)+a_{n-1}(x)y(x+n-1)+\cdots+a_{0}(x)y(x)=0.$

So we can express $y(r)$ in terms of $y(r-1),\ldots,y(r-n)$ where $r=x+n$. The
expression is not defined when $r$ is a root of the denominator, which is
$a_{n}(r-n)$. Hence we define:

###### Definition 2.2.

Let $L=\sum_{i=0}^{n}a_{i}(x)\tau^{i}\in\mathcal{P}$ be primitive. If
$a_{n}\neq 0$ then define ${\rm ord}(L):=n$ and let $\mathfrak{lc}(L)$ be
$a_{n}(x-n)$ divided by its leading coefficient (to make it monic). The
_singularities_ of $L$ are the monic irreducible factors of $\mathfrak{lc}(L)$
in $F[x]$ (or equivalently, their roots in $\overline{F}$) with their
multiplicities.

###### Example 2.3.

Let

$L=x^{2}(x^{2}+1)\tau-(x+1)(x^{2}+2x+2)\in\mathbb{Q}[x][\tau].$

Substituting $x\mapsto x-{\rm ord}(L)$ in the leading coefficient we obtain
$\mathfrak{lc}(L)=(x-1)^{2}((x-1)^{2}+1)$. With repetition indicating
multiplicity, the singularities are $x-1$, $x-1$, $(x-1)^{2}+1$, or
equivalently $1$, $1$, $1\pm\sqrt{-1}$. Let

$A={\frac{1}{\left(x+1\right)\left({x}^{2}+2\,x+2\right)}}(10\,\tau+11\,{x}^{2}+15\,x+14).$

Then

$AL=10\left(x+1\right){\tau}^{2}+\left(11\,{x}^{3}-18\,{x}^{2}+35\,x-50\right)\tau-11\,{x}^{2}-15\,x-14.$

Now $\mathfrak{lc}(AL)=x+1-{\rm ord}(AL)=x-1$ (we divided by 10 to make it
monic). Compared to $L$, the singularity $(x-1)^{2}+1$ disappeared, as well as
one of the two copies of $x-1$.

###### Lemma 2.4 ((Jaroschek, , Theorem 4.1.7, Corollary 4.1.9) ).

Let $L\in\mathcal{P}$. For $k=0,1,2,\ldots$ let

(2) $\mathcal{I}_{k}=\\{0\\}\,\cup\,\\{\mathfrak{lc}(AL)\ :\
A\in\mathcal{D},AL\in\mathcal{P},{\rm ord}(A)=k\\}$

and let $\mathcal{I}_{\infty}$ be their union. Then
$\mathcal{I}_{0}\subseteq\mathcal{I}_{1}\subseteq\mathcal{I}_{2}\subseteq\cdots\subseteq\mathcal{I}_{\infty}$
are ideals in $F[x]$.

###### Proof.

To show $\mathcal{I}_{k}\subseteq\mathcal{I}_{k+1}$, take a non-zero
$a\in\mathcal{I}_{k}$. So $a=\mathfrak{lc}(AL)$ for some $A$ of order $k$.
Replacing $A$ by $\tau A$ shows that $a\in\mathcal{I}_{k+1}$. Clearly
$\mathcal{I}_{k}$ is closed under $F[x]$-multiplication, and it is not
difficult to show that it is closed under addition as well. ∎

###### Definition 2.5 (Essential parts, removable parts, (Jaroschek, ,
Definitions 4.1.8 and 4.1.10)).

With notations as in Lemma 2.4, let $\mathfrak{lc}_{k}(L)$ be the monic
generator of $\mathcal{I}_{k}$ for $k=0,1,2,\ldots,\infty$. Call
$\mathfrak{lc}_{k}(L)$ the _essential part of the leading coefficient at order
$k$_. Note that $\mathfrak{lc}_{0}(L)=\mathfrak{lc}(\mathrm{Prim}(L))$ and
$\mathfrak{lc}_{l}(L)$ divides $\mathfrak{lc}_{k}(L)$ if $l\geq k$. Let
$\mathfrak{rp}_{k}(L)=\frac{\mathfrak{lc}(L)}{\mathfrak{lc}_{k}(L)}\in F[x]$
and call it the _removable part of the leading coefficient at order $k$_.

Factors (or roots) of $\mathfrak{lc}(L)$ are divided into two (possibly
overlapping) sublists: Factors of $\mathfrak{lc}_{\infty}(L)$ are the _true
singularities_ of $L$. Factors of $\mathfrak{rp}_{\infty}(L)$ are the
_apparent singularities_. In 2.3 the true singularity is $x-1$ and the
apparant singularities are $x-1$, $(x-1)^{2}-1$.

###### Definition 2.6.

Let $L\in\mathcal{P}$ be primitive. We call $A\in\mathcal{D}$ a
_desingularizer_ if $AL\in\mathcal{P}$. Such $A$ is _trivial_ if
$A\in\mathcal{P}$ (that implies $\mathfrak{lc}(L)\ |\ \mathfrak{lc}(AL)$, so
no singularities were removed). A desingularizer $A$ is _optimal at order $k$_
if $\mathfrak{lc}(AL)=\mathfrak{lc}_{k}(L)$ and $\mathop{\rm ord}(A)\leq k$;
when $k=\infty$ say $A$ is _optimal_ since it removes all apparant
singularities while introducing no new singularities.

An optimal desingularizer at order $k$ exists because $\mathcal{I}_{k}$ is a
principal ideal; $F[x]$ is a PID.

###  2.2. LCLM method for desingularization

Desingularization of recurrence operators has been studied in (abramov1999, ),
(vanhoeij2006, ), (chen2013, ), (Jaroschek, ), (CHEN2016617, ) and (Zhang2016,
). Papers (abramov1999, ), (vanhoeij2006, ), (chen2013, ), (Jaroschek, ) aim
for full desingularization; (Zhang2016, ) focuses on desingularization over
$R[x]$ where $R$ is not a field. Here we describe the so-called _LCLM method_
, published in (CHEN2016617, ). LCLM stands for the least common left multiple
and GCRD the greatest common right divisor. The main result of (CHEN2016617, )
is restated below, where we focus on the recurrence case while the original
version also applies to other types of Ore operators.

###### Theorem 2.7 (Reformulation of Theorem 6 in (CHEN2016617, )).

Suppose $L\in\mathcal{P}$. Introduce new constants $c_{0},c_{1},\ldots,c_{k}$
that are algebraically independent over $F$. Denote
$A=c_{0}+c_{1}\tau+\cdots+c_{k}\tau^{k}$ and
$L^{\prime}=\mathrm{Prim}(\mathrm{LCLM}(L,A))\in
F(c_{0},c_{1},\ldots,c_{k})[x][\tau]$. Then
$\mathfrak{lc}(L^{\prime})=\mathfrak{lc}_{k}(L)f$ where $f\in
F(c_{0},c_{1},\ldots,c_{k})[x]$ has no non-trivial factor in $F[x]$.

The original form of Theorem 6 in (CHEN2016617, ):
_Let $q$ be an irreducible polynomial which appears with multiplicity $e$ in
$\mathfrak{lc}(L)$ and let $m\leq e$ be maximal such that
$q^{m}\mid\frac{\mathfrak{lc}(L)}{\mathfrak{lc}_{k}(L)}$ for $k\in\mathbb{N}$.
Let $A=c_{0}+c_{1}\tau+\cdots+c_{k}\tau^{k}$ in
$F(c_{0},\ldots,c_{k})[x][\tau]$, where $c_{0},\ldots,c_{k}$ are new constants
that are algebraically independent over $F$. Denote
$L^{\prime}=\mathrm{Prim}(\mathrm{LCLM}(L,A))$. Then the multiplicity of $q$
in $\mathfrak{lc}(L^{\prime})$ is $e-m$. _

The proof in (CHEN2016617, ) also holds when $e=m=0$, which results in Theorem
2.7. It was stated for the case $\mathrm{char}(F)=0$, but the proof is valid
for positive characteristic as well.

###### Remark 1.

Theorem 2.7 implies $\mathfrak{lc}_{k}(L)$ stays the same if $L$ is viewed as
an operator in $E(x)[\tau]$ where $E$ is a field extension of $F$, as the
field extension does not affect
$L^{\prime}=\mathrm{Prim}(\mathrm{LCLM}(L,A))$.

Theorem 2.7 implies the following desingularization algorithm.

Input : a primitive operator $L=\sum_{i=0}^{n}a_{i}\tau^{i}\in F[x][\tau]$ and
positive integer $k$

Output : $\mathfrak{lc}_{k}(L)$

1

2$A\leftarrow\sum_{i=0}^{k}c_{i}\tau^{i}$, where $c_{0},c_{1},\ldots,c_{k}$
are new constants that are algebraically independent over $F$;

3$L^{\prime}\leftarrow\mathrm{Prim}(\mathrm{LCLM}(A,L))$;

return _$\gcd(\mathfrak{lc}(L),\mathfrak{lc}(L^{\prime})$_ ;

Algorithm 1 LCLM_Method

The discussion following the main theorem in (CHEN2016617, ) states that in
characteristic 0, instead of new constants, we can let
$c_{0},c_{1},\ldots,c_{k}$ be random elements in $F$. In this case the
algorithm is Monte-Carlo, meaning it returns the desired result with a high
probability. The Monte-Carlo version is much faster since it avoids
computations in a transcendental extension of $F$. It was implemented (HomURL,
) with $k=1$ by the second author in 2004.

We refer to the algorithm where $k=1$ as _order-1 LCLM method_. The Monte-
Carlo version of order-1 LCLM method is useful in practice since it strikes a
good balance between benefit and cost. The LCLM computation is much faster for
$k=1$ than for larger $k$, and $\mathfrak{lc}_{1}(L)$ is often very close to
$\mathfrak{lc}_{\infty}(L)$. We will further speed up the algorithm in Section
6.

###  2.3. p-Characteristic polynomial

From here until Section 6, $F$ will be a field of characteristic $p$, where
$p$ is a prime number.

A general theory of linear difference equations in positive characteristic is
developed in (van2003galois, , Chapter 5). In
(Bostan:2014:FAC:2608628.2608650, ), $p$-characteristic polynomials of
recurrence operators (and differential operators) over $\mathbb{F}_{p}[x]$ are
studied and an algorithm for computing them is given. An algorithm for
computing $p$-characteristic polynomials of operators in $\mathbb{Z}[x][\tau]$
for a number of $p$ is presented in (pages, ), based on the algorithm from
(Bostan:2014:FAC:2608628.2608650, ). We will give more information about these
algorithms in subsection 2.4

Let $Z=x^{p}-x=x(x+1)\cdots(x+p-1)$. Clearly $Z$ is fixed by $\tau$ and hence
elements of $F(Z)$ are $\tau$-constant. In fact, $F(Z)$ _is_ the field of
$\tau$-constants; to see this, notice that $F(x)$ is a degree $p$ field
extension of $F(Z)$ and hence there is no proper intermediate field.

Let $\mathcal{N}:F(x)\rightarrow F(Z)$ denote the norm map of the field
extension $F(x)/F(Z)$. It is given by the formula

$\mathcal{N}:f(x)\mapsto f(x)f(x+1)\cdots f(x+p-1).$

Denote $T=\tau^{p}$. The center of $\mathcal{D}$ is $F(Z)[T]$. Since
$F(x)[T]\subseteq\mathcal{D}$, any $\mathcal{D}$-module is naturally an
$F(x)[T]$-module, that is, an $F(x)$-vector space equipped with an
$F(x)$-linear map.

###### Definition 2.8.

For a $\mathcal{D}$-module $M$, call the $F(x)$-linear map induced by $T$ the
_$p$ -curvature_ of $M$.

For an operator $L\in\mathcal{D}$, define its $p$-curvature to be that of
$\mathcal{D}/\mathcal{D}L$. Denote $\chi(L)\in F(x)[T]$ its characteristic
polynomial with $T$ as its variable and call it the _$p$ -characteristic
polynomial_ of $L$.

A characteristic polynomial is monic by definition so the leading coefficient
of $L$ is lost in $\chi(L)$. To reinsert it, denote
$\tilde{\chi}(L)=\mathcal{N}(\mathfrak{lc}(L))\chi(L)$. It is called the
_reduced norm_ of $L$ in (Bostan:2014:FAC:2608628.2608650, ).

###### Lemma 2.9.

Properties of $p$-characteristic polynomials.

  1. (i)

For $L\in\mathcal{D}$, $\chi(L)\in F(Z)[T]$.

  2. (ii)

For $L\in\mathcal{D}$, $\chi(L)\in\mathcal{D}L$.

  3. (iii)

For $L_{1},L_{2}\in\mathcal{D}$, $\chi(L_{1}L_{2})=\chi(L_{1})\chi(L_{2})$ and
$\tilde{\chi}(L_{1}L_{2})=\tilde{\chi}(L_{1})\tilde{\chi}(L_{2})$.

  4. (iv)

For $L_{1},L_{2}\in\mathcal{D}$, if $\,\mathrm{GCRD}(L_{1},L_{2})=1$, then

$\chi(\mathrm{LCLM}(L_{1},L_{2}))=\chi(L_{1})\chi(L_{2}).$

  5. (v)

For $L\in\mathcal{P}$, $\tilde{\chi}(L)\in F[Z][T]$ and
$\deg_{Z}(\tilde{\chi}(L))\leq\deg_{x}(L)$.

  6. (vi)

If $L\in F(Z)[T]$ then $\tilde{\chi}(L)=L^{p}$.

###### Proof.

All except item (iv) are proved in (Bostan:2014:FAC:2608628.2608650, , Section
3) for the case $F=\mathbb{F}_{p}$ and the proofs are valid for a general
field $F$ with positive characteristic. We now prove (iv). Denote
$L=\mathrm{LCLM}(L_{1},L_{2})$. If $\mathrm{GCRD}(L_{1},L_{2})=1$ then
$\mathcal{D}/\mathcal{D}L\cong\mathcal{D}/\mathcal{D}L_{1}\oplus\mathcal{D}/\mathcal{D}L_{2}$
as $\mathcal{D}$-modules (and hence as $F(x)[T]$-modules). Now (iv) follows
from the fact that characteristic polynomials are multiplicative on direct
sums. ∎

Lemma 2.9(iii) implies that an operator factors only when its
$p$-characteristic polynomial factors (as a polynomial in $F(Z)[T]$). In fact,
the $p$-characteristic polynomial tells us even more. See
(cluzeaufactorization, ) and (VANDERPUT1996367, ) for discussions on this
topic in the differential case. The $p$-characteristic polynomial is also
useful for testing or proving irreducibility of operators in
$\mathbb{Q}(x)[\tau]$ by reduction modulo $p$.

###  2.4. BCS algorithm and Pagès’ algorithm

Bostan, Caruso and Schost (2015) present an algorithm for computing the
$p$-characteristic polynomial of an operator in $\mathbb{F}_{p}[x][\tau]$,
called Xi_theta_d in (Bostan:2014:FAC:2608628.2608650, ). We refer to it as
the BCS algorithm. Their implementation in Magma is available at
https://github.com/schost.

The BCS algorithm takes a prime $p$ and a difference operator
$L\in\mathbb{F}_{p}[x][\tau]$ as its input and computes
$\tilde{\chi}(L)\in\mathbb{F}_{p}[Z][T]$ (making $\tilde{\chi}(L)$ monic gives
$\chi(L)$). The algorithm computes $\tilde{\chi}(L)$ in
$\mathbb{F}_{p}[[Z]][T]$ to precision $O(Z^{\deg_{x}(L)+1})$ which suffices by
Lemma 2.9(v).

For $L\in\mathcal{P}$, the first part of Lemma 2.9(v) implies that
$\mathcal{N}(\mathfrak{lc}(L))$ is a denominator bound for $\chi(L)$. The BCS
algorithm uses this bound to ensure that what it computes in
$\mathbb{F}_{p}[[Z]][T]$ is in $\mathbb{F}_{p}[Z][T]$, not just in
$\mathbb{F}_{p}(Z)[T]$. We will show that (partial) desingularization leads to
sharper denominator bounds. That reduces the required $Z$-adic precision,
speeding up the computation. In fact, our main result Theorem 3.1 says that
full desingularization gives the exact denominator.

For an operator $L\in\mathbb{Q}(x)[\tau]$, denote by $\chi_{p}(L)$ the
$p$-characteristic polynomial of its reduction modulo $p$. Pagès (2021) gives
an algorithm for computing $\chi_{p}(L)$ for a number of primes at the same
time, if $L\in\mathbb{Z}[x][\tau]$ has a leading coefficient in $\mathbb{Z}$
((pages, , Algorithm 3)). The algorithm is based on the BCS algorithm.

##  3\. Main Theorem and Corollaries

Let $\operatorname{denom}(\cdot)$ be the monic denominator of a rational
function, or of a polynomial with rational function coefficients. We will use
this notation in the cases $F[x]\subset F(x)$ and $F[Z][T]\subset F(Z)[T]$.

###### Theorem 3.1.

For $L\in F[x][\tau]$,
$\operatorname{denom}(\chi(L))=\mathcal{N}(\mathfrak{lc}_{\infty}(L))$.

The theorem quickly implies two corollaries, expressed in terms of the
following definition.

###### Definition 3.2.

Let $r_{1},r_{2}\in F(x)-\\{0\\}$. We say that $r_{1}$ and $r_{2}$ are _shift
equivalent_ , denoted $r_{1}\sim r_{2}$, if $\tau-\frac{r_{1}}{r_{2}}$ has a
non-zero solution in $F(x)$, in other words, if there exists $f\in
F(x)-\\{0\\}$ for which $\frac{r_{1}}{r_{2}}=\frac{\tau(f)}{f}$.

If $r(x)$ has a factor $q(x)$ in the numerator or denominator, and one
replaces $q(x)$ by its shift $q(x+1)$, then the result is shift-equivalent to
$r(x)$. Note that $r_{1}\sim r_{2}$ if and only if
$\mathcal{N}(r_{1})=\mathcal{N}(r_{2})$.

###### Corollary 3.3.

If $\mathcal{D}/\mathcal{D}L_{1}\cong\mathcal{D}/\mathcal{D}L_{2}$ for
$L_{1},L_{2}\in\mathcal{D}$, then $\mathfrak{lc}_{\infty}(L_{1})$ and
$\mathfrak{lc}_{\infty}(L_{2})$ are shift equivalent, so $L_{1}$ and $L_{2}$
have the same true singularities up to shifts.

###### Corollary 3.4.

For $L_{1},L_{2}\in\mathcal{D}$, if

  * •

$L=L_{1}L_{2}$, or

  * •

$L=\mathrm{LCLM}(L_{1},L_{2})$ and $\mathrm{GCRD}(L_{1},L_{2})=1$

then $\mathfrak{lc}_{\infty}(L)$ and
$\mathfrak{lc}_{\infty}(L_{1})\mathfrak{lc}_{\infty}(L_{2})$ are shift
equivalent.

These corollaries are not true in characteristic 0, so we did not expect
Theorem 3.1.

##  4\. Proof of the Main Theorem

This section is devoted to the proof of Theorem 3.1. We start with an easy
lemma.

###### Lemma 4.1.

For any $A,L\in\mathcal{D}$

$\operatorname{denom}(\chi(AL))=\operatorname{denom}(\chi(A))\operatorname{denom}(\chi(L)).$

###### Proof.

Lemma 2.9(iii) and Gauss’s lemma for $F(Z)[T]$ gives

$\mathrm{Cont}(\chi(AL))=\mathrm{Cont}(\chi(A))\mathrm{Cont}(\chi(L)).$

But $\chi(\cdot)$ is monic by definition, and the denominator of a monic
polynomial is the reciprocal of the content. ∎

###  4.1. Special case, Gaussian operators

###### Lemma 4.2.

Let $L\in\mathcal{P}$. The following are equivalent.

  1. 1.

$L$ is Gaussian, i.e. $\forall_{A\in\mathcal{D}}\
AL\in\mathcal{P}\Longrightarrow A\in\mathcal{P}$.

  2. 2.

Every desingularizer is trivial.

  3. 3.

${\rm Cl}(L)=\mathcal{P}L$, where ${\rm
Cl}(L):=\mathcal{D}L\bigcap\mathcal{P}$. (This is called the _Weyl closure_ in
(Tsai00weylclosure, ).)

  4. 4.

$\mathfrak{lc}(L)=\mathfrak{lc}_{\infty}(L)$, i.e. there are no apparant
singularities.

###### Proof.

Items 2 and 3 are reformulations of item 1, and immediately imply item 4. It
remains to show that item 4 implies item 1. Suppose that
$\mathfrak{lc}(L)=\mathfrak{lc}_{\infty}(L)$ and $AL\in\mathcal{P}$. To prove:
$A\in\mathcal{P}$.

By partial fraction decomposition, $A=A_{1}+A_{2}$ where $A_{1}\in\mathcal{P}$
and $A_{2}=\sum q_{i}\tau^{i}$ with the numerator of $q_{i}$ having lower
degree than its denominator. Since $AL$ and $A_{1}L$ are in $\mathcal{P}$,
their difference $A_{2}L$ is in $\mathcal{P}$ as well. If $A_{2}\neq 0$, then
the leading coefficient of $A_{2}L$ will have lower degree than
$\mathfrak{lc}(L)$, contradicting item 4. Thus $A_{2}=0$ and hence
$A\in\mathcal{P}$. ∎

Lemma 2.9(v) says that
$\tilde{\chi}(L)=\mathcal{N}(\mathfrak{lc}(L))\chi(L)\in F[Z][T]$ when
$L\in\mathcal{P}$, in other words

(3) $\operatorname{denom}(\chi(L))\mid\mathcal{N}(\mathfrak{lc}(L)).$

Here we give a sharper denominator bound.

###### Lemma 4.3.

For $L\in F[x][\tau]$,
$\operatorname{denom}(\mathrm{\chi}(L))\mid\mathcal{N}(\mathfrak{lc}_{\infty}(L)).$

###### Proof.

Let $A\in\mathcal{D}$ be an optimal desingularizer of $L$, then
$\mathfrak{lc}(AL)=\mathfrak{lc}_{\infty}(L)$. From Lemma 4.1 and Equation (3)
applied to $AL$,
$\operatorname{denom}(\chi(L))\mid\operatorname{denom}(\chi(AL))\mid\mathcal{N}(\mathfrak{lc}(AL))=\mathcal{N}(\mathfrak{lc}_{\infty}(L)).$
∎

Next we show that our denominator bound is exact for Gaussian operators. The
next section will prove the general case by exploiting the fact that any
operator has a Gaussian multiple.

###### Lemma 4.4.

If $L\in F[x][\tau]$ is Gaussian, then

$\operatorname{denom}(\mathrm{\chi}(L))=\mathcal{N}(\mathfrak{lc}_{\infty}(L)).$

###### Proof.

Denote $f=\mathrm{Prim}(\chi(L))\in F[Z][T]$. By Lemma 2.9(ii),
$f\in\mathcal{D}L$ so there exists $Q\in\mathcal{D}$ such that

(4) $QL=f.$

In fact $Q\in\mathcal{P}$ since $L$ is Gaussian. Lemma 2.9(v) says
$\tilde{\chi}(Q),\tilde{\chi}(L)\in F[Z][T]$. Applying $\tilde{\chi}$ to
Equation (4), and Lemma 2.9(vi), gives

(5) $\tilde{\chi}(Q)\tilde{\chi}(L)=\tilde{\chi}(f)=f^{p}.$

Now $f^{p}\in F[Z][T]$ is primitive since $f$ is primitive. Then Gauss’s lemma
implies $\tilde{\chi}(L)\in F[Z][T]$ is primitive. It follows that

$\operatorname{denom}(\chi(L))=\mathfrak{lc}(\tilde{\chi}(L))=\mathcal{N}(\mathfrak{lc}(L))=\mathcal{N}(\mathfrak{lc}_{\infty}(L))$

where the last equality comes from Lemma 4.2, part 4. ∎

###  4.2. Proof for the general case

###### Lemma 4.5.

Suppose $L\in\mathcal{P}$ and
$A=\sum_{i=0}^{k}\frac{n_{i}}{d_{i}}\in\mathcal{D}$ is a desingularizer of
$L$, where $\frac{n_{i}}{d_{i}}\in F(x)$ is in lowest terms for each $i$. Then
$\mathcal{N}(d_{k})\mid\mathcal{N}(\mathfrak{rp}_{\infty}(L))$.

###### Proof.

Let $n$ be the order of $L$. The definition of $\mathfrak{lc}_{\infty}$ and
the product of the leading terms of $A$ and $L$ gives

$\mathfrak{lc}_{\infty}(L)\mid\tau^{-k-n}(\frac{n_{k}}{d_{k}})\mathfrak{lc}(L)$

and hence $\tau^{-k-n}(\frac{n_{k}}{d_{k}})\mathfrak{rp}_{\infty}(L)\in F[x]$.
Since $\frac{n_{k}}{d_{k}}$ is a reduced fraction, we have
$\tau^{-k-n}(d_{k})\mid\mathfrak{rp}_{\infty}(L)$, which leads to

$\mathcal{N}(d_{k})=\mathcal{N}(\tau^{-k-n}(d_{k}))\mid\mathcal{N}(\mathfrak{rp}_{\infty}(L)).$

∎

###### Lemma 4.6.

Suppose $L\in\mathcal{P}$ and $A\in\mathcal{D}$ is an optimal desingularizer
of $L$. Then there exists a positive integer $N$ such that

$\operatorname{denom}(\chi(A))\mid(\mathfrak{rp}_{\infty}(L))^{N}.$

###### Proof.

Write $A=\sum_{i=0}^{k}\frac{n_{i}}{d_{i}}\tau^{i}$, where
$\frac{n_{i}}{d_{i}}\in F(x)$ is a reduced fraction for each $i$. We deduce
$n_{k}=1,d_{k}=\tau^{n+k}(\mathfrak{rp}_{\infty}(L))$ from the fact that
$\mathfrak{lc}(AL)=\mathfrak{lc}_{\infty}(L)$. Clearly $d_{0}d_{1}\cdots
d_{k}A\in\mathcal{P}$. Then by Equation (3)

$\operatorname{denom}(\chi(A))\mid\mathcal{N}(\mathfrak{lc}(d_{0}d_{1}\cdots
d_{k}A))=\mathcal{N}(d_{0}d_{1}\cdots d_{k-1}).$

Now we bound $\mathcal{N}(d_{i})$ in terms of $\mathfrak{rp}_{\infty}(L)$. Let

$A_{j}=d_{k}d_{k-1}\cdots d_{j+1}(\sum_{i=0}^{j}\frac{n_{i}}{d_{i}}\tau^{i})$

for $j=0,1,\ldots,k-1$. Notice that

$A_{j}-d_{k}d_{k-1}\cdots d_{j+1}A=-d_{k}d_{k-1}\cdots
d_{j+1}(\sum_{i=j+1}^{k}\frac{n_{i}}{d_{i}}\tau^{i})\in\mathcal{P}.$

This implies $A_{j}L\in\mathcal{P}$, or equivalently, $A_{j}$ is a
desingularizer of $L$. Apply Lemma 4.5 to $A_{j}$:

$\mathcal{N}(\operatorname{denom}(d_{k}d_{k-1}\cdots
d_{j+1}\frac{n_{j}}{d_{j}}))\mid\mathcal{N}(\mathfrak{rp}_{\infty}(L)).$

Notice that

$d_{j}\mid d_{k}d_{k-1}\cdots
d_{j+1}\cdot\operatorname{denom}(d_{k}d_{k-1}\cdots
d_{j+1}\frac{n_{j}}{d_{j}}).$

Therefore

$\mathcal{N}(d_{j})\mid\mathcal{N}(d_{k}d_{k-1}\cdots
d_{j+1})\cdot\mathcal{N}(\mathfrak{rp}_{\infty}(L)).$

Recall that $\mathcal{N}(d_{k})\mid\mathcal{N}(\mathfrak{rp}_{\infty}(L))$. By
downward induction on $j$, we conclude that

$\mathcal{N}(d_{j})\mid(\mathcal{N}(\mathfrak{rp}_{\infty}(L)))^{k-j+1}.$

∎

We are now ready to finish the proof for Theorem 3.1.

###### Proof of Theorem 3.1.

It remains to show that
$\mathcal{N}(\mathfrak{lc}_{\infty}(L))\mid\operatorname{denom}(\mathrm{\chi}(L))$
for any $L\in\mathcal{D}$. There exists a sufficiently large $k$ such that
$\mathfrak{lc}_{k}(L)=\mathfrak{lc}_{\infty}(L)$. Introduce new constants
$c_{0},\ldots,c_{k}$ that are algebraically independent over $F$ and denote
$E=F(c_{0},c_{1},\ldots,c_{k})$. Let

$L^{\prime}=\mathrm{Prim}(\mathrm{LCLM}(c_{k}\tau^{k}+\cdots+c_{0},L))\in
E[x][\tau].$

Theorem 2.7 says $\mathfrak{lc}(L^{\prime})=\mathfrak{lc}_{\infty}(L)f,$ where
$f\in E[x]$ has no non-trivial factor in $F[x]$. It follows from Definition
2.5 (see also Equation 2) that
$\mathfrak{lc}_{\infty}(L)\mid\mathfrak{lc}_{\infty}(L^{\prime})$ and hence
$\mathfrak{rp}_{\infty}(L^{\prime})\mid f$. Remark 1 guarantees
$\mathfrak{lc}_{\infty}(L)$ does not change as we shift from $F$ to $E$. Let
$A$ be an optimal desingularizer of $L^{\prime}$. By Lemma 2.9 ((iii) and
(iv)), we have

(6) $\chi(AL^{\prime})=\chi(A)\chi(c_{k}\tau^{n}+\cdots+c_{0})\chi(L).$

Lemma 2.9(vi) implies

$\operatorname{denom}(\chi(c_{k}\tau^{n}+\cdots+c_{0}))=1.$

Applying Lemma 4.1 to Equation (6) gives

$\operatorname{denom}(\chi(AL^{\prime}))=\operatorname{denom}(\chi(A))\operatorname{denom}(\chi(L)).$

Since $AL^{\prime}$ is Gaussian, we know from Lemma 4.4

$\operatorname{denom}(\chi(AL^{\prime}))=\mathcal{N}(\mathfrak{lc}_{\infty}(AL^{\prime})),$

which equals $\mathcal{N}(\mathfrak{lc}_{\infty}(L^{\prime}))$ since $A$ is an
optimal desingularizer of $L^{\prime}$. As a consequence,

(7)
$\mathcal{N}(\mathfrak{lc}_{\infty}(L))\mid\mathcal{N}(\mathfrak{lc}_{\infty}(L^{\prime}))=\operatorname{denom}(\chi(A))\operatorname{denom}(\chi(L)).$

Lemma 4.6 says $\operatorname{denom}(\chi(A))$ is a factor of $f^{N}$ for some
sufficiently large $N$, so $\operatorname{denom}(\chi(A))$ has no non-trivial
factor in $F[x]$. By taking only factors in $F[x]$ in Equation 7, we obtain
the desired result
$\mathcal{N}(\mathfrak{lc}_{\infty}(L))\mid\operatorname{denom}(\chi(L)).$

∎

##  5\. Application to Computations

In this section $F=\mathbb{F}_{p}$.

###  5.1. Algorithm

Let $L\in\mathcal{P}$ and $\alpha=\mathfrak{rp}_{k}(L)$. Theorem 3.1 implies
that $\mathcal{N}(\alpha)$, which is in $\mathbb{F}_{p}[Z]$, is a factor of
$\tilde{\chi}(L)$. Dividing this factor away reduces the degree bound from
Lemma 2.9(v) to

(8)
$\deg_{Z}(\,\mathcal{N}(\alpha)^{-1}\tilde{\chi}(L)\,)\leq\deg_{x}(L)-\deg_{x}(\alpha)$

which becomes an equality when $k$ is sufficiently large. However, we use
$k=1$ to minimize the time spent computing $\alpha$. The reduced degree bound
allows us to recover $\chi(L)$ from a lower precision $Z$-adic expansion. That
leads to the following algorithm.

Input : prime $p$ and $L\in\mathbb{F}_{p}[x][\tau]$

Output : $\mathrm{Prim}(\chi(L))\in\mathbb{F}_{p}[Z][T]$

1

2Pick $k\geq 1$ and compute $\mathfrak{lc}_{k}(L)$ and
$\alpha:=\mathfrak{rp}_{k}(L)\in\mathbb{F}_{p}[x]$. We use $k=1$ to minimize
the time spent in this step.

3Compute $\mathcal{N}(\alpha)\in\mathbb{F}_{p}[Z]$. Let $v$ be its $Z$-adic
valuation in $\mathbb{F}_{p}[[Z]]$ and let
$\beta=Z^{-v}\mathcal{N}(\alpha)\in\mathbb{F}_{p}[Z]$.
For computing $\mathcal{N}(\cdot)$ see Step 3 of Xi_theta_d in
(Bostan:2014:FAC:2608628.2608650, ).

4Let $d_{1}=\deg_{Z}(\beta)$ and $d:=\deg_{x}(L)$. Apply the BCS algorithm
with $d$ replaced by $d-d_{1}$ to $L$ to obtain $\tilde{\chi}(L)$ up to the
precision $O(Z^{d-d_{1}+1})$. Denote the result by $\chi_{1}$.

5Compute $\beta^{-1}$ in $\mathbb{F}_{p}[[Z]]$ up to $O(Z^{d-d_{1}-v+1})$.
This can be done by applying the extended Euclidean algorithm to $\beta$ and
$Z^{d-d_{1}-v+1}$.

6Compute $\beta^{-1}\cdot(Z^{-v}\chi_{1})$ in $\mathbb{F}_{p}[[Z]][T]$ up to
the precision $O(Z^{d-d_{1}-v+1})$. This gives
$\mathcal{N}(\alpha)^{-1}\chi_{1}\in\mathbb{F}_{p}[Z][T]$.
Return its primitive part (with respect to $T$).
Note: $\mathcal{N}(\alpha)^{-1}\chi_{1}$ and
$\mathcal{N}(\alpha)^{-1}\tilde{\chi}(L)$ agree to precision
$O(Z^{d-d_{1}-v+1})$ which suffices by (8).

Algorithm 2 Xi_p_desing

Step 3 is where we save CPU time over the original algorithm from
(Bostan:2014:FAC:2608628.2608650, ) if $d_{1}>0$. If $d_{1}=0$ then there is
no improvement in efficiency. However, as we will see in the following
section, the extra steps cost very little time.

###  5.2. Implementation and timings

Our Magma implementation of algorithm Xi_p_desing is available at
https://www.math.fsu.edu/~yzhou/magma/, together with experiments on a variety
of operators. One should load the implementation of
(Bostan:2014:FAC:2608628.2608650, ) at https://github.com/schost/pCurvature
(file pCurvature.mgm) prior to ours.

In the following table the data for two operators from OEIS ((A002777, ),
(A151329, )) is presented. Here $d_{1}$ is defined in the Step 3 of
Xi_p_desing; each running time is the average of ten runs.

OEIS index | order |  $x$-degree | $d_{1}$ | BCS | Xi_p_desing
---|---|---|---|---|---
A151329 | 9 | 18 | 10 | 17.2s | 9.6s
A002777 | 4 | 3 | 0 | 6.97s | 7.01s
Table 1. Timings for operators from OEIS. $p=27457$.

For the recurrence for OEIS A002777, we expect Xi_p_desing to be slower than
BCS since $d_{1}=0$. However, the running time difference between two
algorithms is nearly unnoticeable.

We also tested our algorithm on operators that are LCLMs of two operators
(implementations, ). Such operators tend to have many apparent singularities
and hence benefit more from our approach.

##  6\. Fast Algorithms for Desingularization at Order 1

###  6.1. First algorithm

In this section we present our first speedup of the order-1 LCLM method. We
used it for Step 1 of algorithm 2.

The order-1 LCLM method computes $L^{\prime}=\mathrm{LCLM}(L,\tau-c)$ where
$c$ is a new constant (or a random number in the Monte-Carlo version). To
speed this up, our idea is to obtain $\mathfrak{lc}_{1}(L)$ while only
computing a portion of $L^{\prime}$.

First we express of $L^{\prime}$ in terms of $c$ and coefficients of $L$.
Suppose $L=\sum_{i=0}^{n}a_{i}\tau^{i}\in F(x)[\tau]$. Then
$L^{\prime}=\sum_{i=0}^{n}c^{i}L_{i},$ where

(9) $L_{i}=a_{i}\tau L-\tau(a_{i-1})L=(a_{i}\tau-\tau(a_{i-1}))L,$

where $a_{i}=0$ for $i<0$ and $i>n$. Clearly this $L^{\prime}$ is a left
multiple of $L$. To verify it is also a left-multiple of $\tau-c$, use the
fact that the remainder of $\tau^{i}$ right-divided by $\tau-c$ is $c^{i}$. We
skip the tedious computation. As a result $L^{\prime}$ is an LCLM of $L$ and
$\tau-c$.

The order-1 LCLM method computes $L^{\prime}$ which amounts to compute all
$L_{i}$’s. The following proposition shows that one can provably obtain
$\mathfrak{lc}_{1}(L)$ from just a subset of the $L_{i}$’s.

###### Proposition 6.1.

Let $L=\sum_{i=0}^{n}a_{i}\tau^{i}\in F[x][\tau]$. Let $L_{i}$ be defined by
Equation 9, where $a_{i}=0$ for $i<0$ and $i>n$. If

(10) $\gcd(a_{i_{1}},a_{i_{2}},\ldots,a_{i_{k}})=1,$

then

$\gcd(\mathfrak{lc}_{0}(L_{i_{1}}),\mathfrak{lc}_{0}(L_{i_{2}}),\ldots,\mathfrak{lc}_{0}(L_{i_{k}}))=\mathfrak{lc}_{1}(L).$

The proof will be given in the next section. Note that there exist
$i_{1},i_{2},\ldots,i_{k}$ satisfying the gcd condition (Equation 10) if and
only if $L$ is primitive. The proposition immediately implies algorithm 3.

Input : a primitive operator $L=\sum_{i=0}^{n}a_{i}\tau^{i}\in F[x][\tau]$

Output : $\mathfrak{lc}_{1}(L)$

1

2Find $I\subset\\{0,1,\ldots,n\\}$ such that $a_{i}\neq 0$ for any $i\in I$
and $\gcd(a_{i}\mid i\in I)=1$. Note: the algorithm is still correct if we
allow $a_{i}=0$, but that $i$ is redundant since it does not affect the gcd at
all.

3Compute $L_{i}$ for $i\in I$ by Equation 9.

4Return $\gcd(\mathfrak{lc}_{0}(L_{i}):i\in I)$.

Algorithm 3 $\mathfrak{lc}_{1}$

###### Remark 2.

Computing $\mathfrak{lc}_{0}(L_{i})=\mathfrak{lc}(\mathrm{Prim}(L_{i}))$ is
the most time-consuming part in the algorithm, because $L_{i}$ has twice the
$x$-degree as $L$.

###  6.2. Proof

Always assume $L=\sum_{i=0}^{n}a_{i}\tau^{i}\in F[x][\tau]$ is primitive and
$L_{i}$ is defined by Equation 9.

###### Lemma 6.2.

There exists $b\in F[x]$ such that

$\mathrm{Cont}((\tau-b)L)=\tau^{n+1}(\mathfrak{rp}_{1}(L)).$

###### Proof.

Let $A\in\mathcal{D}$ be an optimal desingularizer of $L$ at order $1$. Then
$A=\frac{1}{d_{1}}\tau-\frac{n_{2}}{d_{2}}$, where
$d_{1}=\tau^{n+1}(\mathfrak{rp}_{1}(L))$ and $\frac{n_{2}}{d_{2}}\in F(x)$ is
a reduced fraction. Let $b=d_{1}\frac{n_{2}}{d_{2}}$. Observe that

(11) $bL=\tau\cdot L-d_{1}AL\in\mathcal{P}.$

Due to $L$ being primitive, $b$ has to be a polynomial. Since $A$ is an
optimal desingularizer of $L$ at order 1, $AL\in\mathcal{P}$ is primitive;
otherwise dividing out the content of $AL$ yields a more optimal
desingularizer. By rearranging Equation 11 we see that
$\frac{1}{d_{1}}(\tau-b)L=AL$ is primitive, which completes the proof. ∎

###### Theorem 6.3.

Let $C=(c_{0},c_{1},\ldots,c_{n+1})\in F^{n+2}$. Denote

$C_{1}=\sum_{i=0}^{n+1}c_{i}a_{i},\quad
C_{0}=\sum_{i=0}^{n+1}c_{i}\tau(a_{i-1}),\quad L^{\prime}=(C_{1}\tau-C_{0})L.$

Then

$\mathfrak{lc}_{1}(L)\mid\mathfrak{lc}_{0}(L^{\prime})\mid\tau^{-n-1}(C_{1})\mathfrak{lc}_{1}(L).$

###### Proof.

Assume $C_{1}\neq 0$ since otherwise it is trivial.

The relation $\mathfrak{lc}_{1}(L)\mid\mathfrak{lc}_{0}(L^{\prime})$
immediately follows from the definition of $\mathfrak{lc}_{1}$. By Lemma 6.2,
there exists $b\in F[x]$ such that
$\mathrm{Cont}((\tau-b)L)=\tau^{n+1}(\mathfrak{rp}_{1}(L))$. Since

$(\tau-b)L=\sum_{i=0}^{n+1}(\tau(a_{j-1})-ba_{j})\tau^{j},$

we have

(12) $\gcd(\tau(a_{j-1})-ba_{j}\mid
j=0,1,\ldots,n+1)=\tau^{n+1}(\mathfrak{rp}_{1}(L)).$

Notice that

$L^{\prime}=C_{1}(\tau-b)L+(bC_{1}-C_{0})L,$

and in particular

$bC_{1}-C_{0}=\sum_{i=0}^{n+1}bc_{i}a_{i}-\sum_{i=0}^{n+1}c_{i}\tau(a_{i-1})=\sum_{i=0}^{n+1}c_{i}(ba_{i}-\tau(a_{i-1}))$

is a multiple of $\tau^{n+1}(\mathfrak{rp}_{1})$ due to Equation 12. Hence
$\frac{1}{\tau^{n+1}(\mathfrak{rp}_{1})}L^{\prime}\in F[x][\tau].$ When
$C_{1}\neq 0$,

$\mathfrak{lc}(\frac{1}{\tau^{n+1}(\mathfrak{rp}_{1}(L))}L^{\prime})=\frac{1}{\mathfrak{rp}_{1}(L)}\tau^{-n-1}(C_{1})\mathfrak{lc}(L)=\tau^{-n-1}(C_{1})\mathfrak{lc}_{1}(L).$

Then we have

$\mathfrak{lc}_{0}(L^{\prime})=\mathfrak{lc}(\mathrm{Prim}(L^{\prime}))\mid\mathfrak{lc}(\frac{1}{\tau^{n+1}(\mathfrak{rp}_{1}(L))}L^{\prime})=\tau^{-n-1}(C_{1})\mathfrak{lc}_{1}(L).$

∎

###### Proof of Proposition 6.1.

In Theorem 6.3, setting $c_{i}=1$ for some $i$ and $c_{j}=0$ for any $j\neq i$
yields

$\mathfrak{lc}_{1}(L)\mid\mathfrak{lc}_{0}(L_{i})\mid\tau^{-n-1}(a_{i})\mathfrak{lc}_{0}(L).$

The desired result follows immediately. ∎

###  6.3. Desingularizing both leading and trailing coefficients

The variation in this section handles both leading and trailing singularities.
It uses only one $L_{i}$ (defined in Equation 9) without checking the $\gcd$
condition (Equation 10), since most apparant singularities are already
detected with one $L_{i}$.

In the algorithm, $\mathfrak{tc}_{1}$ denotes the _essential part of trailing
coefficient at order 1_ , which is the counterpart of $\mathfrak{lc}_{1}$ for
trailing coefficients.

Input : a primitive operator $L=\sum_{i=0}^{n}a_{i}\tau^{i}\in F[x][\tau]$
with $a_{0}a_{n}\neq 0$

Output : $l,t\in F[x]$ such that $\mathfrak{lc}_{1}(L)\mid
l\mid\mathfrak{lc}(L)$ and $\mathfrak{tc}_{1}(L)\mid l\mid\mathfrak{tc}(L)$

1

2$i\leftarrow\lfloor\frac{n}{2}\rfloor$

3$L_{i}\leftarrow(a_{i}\tau-\tau(a_{i-1}))L$

4$l,t\leftarrow\mathfrak{lc}_{0}(L_{i}),\mathfrak{tc}_{0}(L_{i})$

5$l,t\leftarrow\gcd(\tau^{-n}(a_{n}),l),\gcd(a_{0},t)$

return _$l,t$_

Algorithm 4 $\mathfrak{lc}$1_$\mathfrak{tc}$1

###  6.4. Examples and comparisons

We have implemented algorithm 3 and algorithm 4 in Maple and SageMath, and
done some experiments to compare the running time of our algorithm with the
order-1 LCLM method. All can be found at (implementations, ). Below we give an
experiment we did in Maple.

###### Example 6.4.

In this example the base field is $\mathbb{Q}$. We took random operators

$L_{1}=(26x^{4}+20)\tau^{11}-96x^{3}\tau^{9}+64x^{5}\tau^{8}+45x^{11}\tau^{4}-x^{2}\tau^{3},$
$L_{2}=-55x^{3}\tau^{7}+85x^{3}\tau^{4}+64x^{4}\tau^{3}+(-14x^{8}-20x^{4})\tau+79x,$

and then computed

$L=\mathrm{Prim}(\mathrm{LCLM}(L_{1},L_{2})).$

The $x$-degree of $L$ is $109$. We desingularize $L$ using three different
algorithms. For the LCLM method we used the Monte-Carlo version and randomly
choose $c=7$. The results are shown in the Table 2, where each time is the
average of ten runs.

algorithms | running time |  $x$-degree in output
---|---|---
Order-1 LCLM | 1.191s | 6
$\mathrm{LC}1$ | 0.055s | 6
$\mathfrak{lc}$1_$\mathfrak{tc}$1 | 0.092s | 6
Table 2. Comparison of different desingularization algorithms

One might expect $\mathfrak{lc}$1_$\mathfrak{tc}$1 to be slower than
$\mathrm{LC}1$ because it treats both the leading and trailing coefficient,
however, we expected it to be faster because it corresponds to taking just one
$L_{i}$ in $\mathrm{LC}1$.

##  7\. Future Work

###  7.1. Application to Pagès’ algorithm

Pagès’ Algorithm computes $\chi_{p}(L)$ for $L\in\mathbb{Z}[x][\tau]$ with
$\mathrm{lc}(L)\in\mathbb{Z}$, but with minor adjustments it applies to all
recurrence operators in $\mathbb{Z}[x][\tau]$. We expect desingularization to
be beneficial here as well.

###  7.2. Differential case

The desingularization improvement should also work for the differential case
or Ore operators. For a differential operator
$L=\sum_{i=0}^{n}a_{i}\partial^{i}\in F[x][\partial]$, we can write
$\mathrm{LCLM}(L,\tau-c)=\sum_{i=0}^{n}c^{i}L_{i}$, where

$L_{i}=(a_{i}\partial-(a_{i-1}+a_{i}^{\prime}))\,L.$

We expect that there should also be a differential analog of our main result,
Theorem 3.1.

## References

  * [1] S. A. Abramov, M. A. Barkatou, and M. van Hoeij.  Apparent singularities of linear difference equations with polynomial coefficients.  Applicable Algebra in Engineering, Communication and Computing, 17:117 – 133, 2006.
  * [2] Sergei Abramov.  Eg-eliminations.  Journal of Difference Equations and Applications, 5:393–433, 01 1999.
  * [3] Alin Bostan, Xavier Caruso, and Éric Schost.  A fast algorithm for computing the characteristic polynomial of the p-curvature.  In Proceedings of the 39th International Symposium on Symbolic and Algebraic Computation, ISSAC ’14, pages 59–66, New York, NY, USA, 2014. ACM.
  * [4] Shaoshi Chen, Maximilian Jaroschek, Manuel Kauers, and Michael F. Singer.  Desingularization explains order-degree curves for ore operators.  In Proceedings of the 38th International Symposium on Symbolic and Algebraic Computation, ISSAC ’13, pages 157–164, New York, NY, USA, 2013\. Association for Computing Machinery.
  * [5] Shaoshi Chen, Manuel Kauers, and Michael F. Singer.  Desingularization of ore operators.  Journal of Symbolic Computation, 74:617 – 626, 2016.
  * [6] Thomas Cluzeau.  Factorization of differential systems in characteristic p.  In Proceedings of the 2003 International Symposium on Symbolic and Algebraic Computation, ISSAC ’03, pages 58–65, New York, NY, USA, 2003. Association for Computing Machinery.
  * [7] OEIS Foundation Inc.  Entry A002777 in the on-line encyclopedia of integer sequences.  http://oeis.org/A002777.
  * [8] OEIS Foundation Inc.  Entry A151329 in the on-line encyclopedia of integer sequences.  http://oeis.org/A151329.
  * [9] Maximilian Jaroschek.  Removable Singularities of Ore Operators.  PhD thesis, RISC, Johannes Kepler University Linz, 2013.
  * [10] Raphaël Pagès.  Computing characteristic polynomials of p-curvatures in average polynomial time.  In Proceedings of the 2021 on International Symposium on Symbolic and Algebraic Computation, ISSAC ’21, pages 329–336, New York, NY, USA, 2021. Association for Computing Machinery.
  * [11] Harrison Tsai.  Weyl closure of a linear differential operator.  J. SYMBOLIC COMPUT, 29:4–5, 2000.
  * [12] M. van der Put and M.F. Singer.  Galois Theory of Linear Differential Equations.  Grundlehren der mathematischen Wissenschaften. Springer Berlin Heidelberg, 2003.
  * [13] Marius van der Put.  Reduction modulo p of differential equations.  Indagationes Mathematicae, 7(3):367–387, 1996.
  * [14] Mark van Hoeij.  Implementation of `DEtools[Homomorphisms]`, added to Maple in 2005. It uses `LCLM` to discard Apparant Singularities. http://www.math.fsu.edu/~hoeij/files/Hom.  2004\.
  * [15] Yi Zhang.  Contraction of ore ideals with applications.  In Proceedings of the ACM on International Symposium on Symbolic and Algebraic Computation, ISSAC ’16, pages 413–420, New York, NY, USA, 2016\. Association for Computing Machinery.
  * [16] Yi Zhou.  Implementations and examples.  http://www.math.fsu.edu/~yzhou/desing, 2021.


:
AABI 20235th Symposium on Advances in Approximate Bayesian Inference, 2023

# Amortised Inference in Neural Networks for Small-Scale Probabilistic Meta-
Learning

Matthew Ashman∗<EMAIL_ADDRESS>
Tommy Rochussen∗<EMAIL_ADDRESS>
Adrian Weller<EMAIL_ADDRESS>
University of Cambridge
* Equal contributions

##  1 Introduction

In many machine learning applications, well-calibrated posterior predictive
distributions are required for a number of closely-related datasets. Given
similarity between datasets, it is natural to wish to develop meta-learning
algorithms that utilise other datasets to reduce the computational complexity
and / or improve predictive performance when deploying models on newly-seen
datasets at test time. There have been a number of significant recent
developments in meta-learning for predictive distributions, most notably that
of the neural process (NP) family (Garnelo et al., 2018a, b; Foong et al.,
2020; Gordon et al., 2018, 2019). Despite the utility of these methods on
large-scale meta-datasets, they perform poorly in settings where the number of
datasets and the total number of datapoints is small. We argue that this is a
result of the large number of shared model parameters overfitting to the meta-
dataset. A natural solution is to remove these shared model parameters, and
instead train a meta-model to learn to approximate fully Bayesian inference
over task-specific model parameters.

Recently, Ober and Aitchison (2021) developed a variational approximation for
Bayesian neural networks (BNNs) based on using a set of inducing inputs to
construct a series of conditional distributions that accurately approximate
the conditionals of the true posterior distribution. Notably, the variational
distribution consists of the prior multiplied by a set of approximate
likelihoods for each inducing input. Our key insight is that these inducing
inputs can be replaced by the actual data, such that the variational
distribution consists of a set of approximate likelihoods for each datapoint.
This structure lends itself to amortised inference, in which the parameters of
each approximate likelihood are obtained by passing each datapoint through a
meta-model known as the inference network. By training this inference network
across related datasets, we can meta-learn Bayesian inference over task-
specific BNNs, addressing the challenge above.

##  2 Related Work

#### Neural processes

Our work is most similar to the NP family (Garnelo et al., 2018a, b), which
seeks to meta-learn predictive distributions either through maximisation of
the posterior predictive likelihood or variational inference (VI) (Foong et
al., 2020). Similar to our method, NPs utilise an encoder to create embeddings
for each datapoint. These embeddings are then aggregated to form a
distribution over a latent variable which is then sampled and passed, together
with a test datapoint, through a decoder. Volpp et al. (2020) propose the use
of Bayesian aggregation, in which embeddings of individual datapoints take the
form of approximate likelihoods which are multiplied together with the prior
to form an approximate posterior distribution over the latent variable. Whilst
these methods differ to ours in their use of shared model parameters, our
method can be reinterpreted as a member of the NP family in which the latent
variables are the parameters of the decoder. Through this perspective we can
train our model in an identical way to NPs.

#### Meta-learning neural networks

Meta-learning for neural networks had received a significant amount of
attention from the research community. Notable examples include MAML Finn et
al. (2017) and its extensions (Yoon et al., 2018; Antoniou et al., 2018),
which seek good parameter initialisation, and those which explicitly condition
on the dataset to obtain task specific parameters (Requeima et al., 2019;
Gordon et al., 2018). Whilst conceptually similar to our approach, these
methods differ in their use of shared model parameters—the task specific
parameters amount to a small subset of the overall model parameters. In
addition to requiring a large meta-dataset, this limits these methods to meta-
datasets in which the individual datasets are very similar. By contrast, our
approach does not use any shared model parameters but rather meta-learns
inference in a BNN. We discuss the relationship between our method and NPs in
Appendix A, demonstrating that they conceptually differ only in a change of
objective functions.

##  3 Background

In this section we review the GI-BNN of Ober and Aitchison (2021) and NP
family (Garnelo et al., 2018b). Throughout, let
$\boldsymbol{\mathbf{\Xi}}=\\{\mathcal{D}\\}$ denote a meta-dataset of $|\Xi|$
datasets, and $\mathcal{D}=\\{\mathbf{X},\mathbf{y}\\}$ denote a dataset
consisting of inputs $\mathbf{X}\in\mathbb{R}^{N\times D_{0}}$ and outputs
$\mathbf{y}\in\mathbb{R}^{N\times P}$.

###  3.1 Global Inducing Point Variational Posteriors for BNNs

Let $\mathbf{W}=\\{\mathbf{W}^{\ell}\\}_{l=1}^{L}$ denote the weights of an
$L$-layer neural network, such that
$\mathbf{W}^{\ell}\in\mathbb{R}^{D^{\ell-1}\times D^{\ell}}$ where $D^{\ell}$
denotes the dimensions of the $\ell$-th hidden layer, and $\psi(\cdot)$ denote
the element-wise activation function acting between layers. Ober and Aitchison
(2021) introduce the global inducing point variational approximation for BNNs,
in which the variational approximation to the posterior
$p(\mathbf{W}|\mathcal{D})$ is defined recursively as
$q_{\phi}(\mathbf{W})=\prod_{\ell=1}^{L}q_{\phi}(\mathbf{W}^{\ell}|\\{\mathbf{W}^{\ell^{\prime}}\\}_{\ell^{\prime}=1}^{\ell-1},\mathbf{U}^{0})$,
where

$q_{\phi}\left(\mathbf{W}^{\ell}|\\{\mathbf{W}^{\ell^{\prime}}\\}_{\ell^{\prime}=1}^{\ell-1},\mathbf{U}_{0}\right)\propto\prod_{d=1}^{D^{\ell}}p\left(\mathbf{w}^{\ell}_{d}\right)\underbrace{\mathcal{N}\left({\mathbf{v}^{\ell}_{d}};{\psi(\mathbf{U}^{\ell})\mathbf{w}^{\ell}_{d}},{\left[\boldsymbol{\mathbf{\Lambda}}^{\ell}_{d}\right]^{-1}}\right)}_{t^{\ell}_{d}(\mathbf{w}^{\ell}_{d})}.$
(1)

This mirrors the structure of the true posterior in the sense that it is
equivalent to the product of the prior and an approximate likelihood,
$t^{\ell}_{d}(\mathbf{w}^{\ell}_{d})$. The variational parameters $\phi$ of
this approximation are the parameters of each approximate likelihood,
$\mathbf{v}^{\ell}_{d}\in\mathbb{R}^{M}$ and
$\boldsymbol{\mathbf{\Lambda}}_{d}^{\ell}\in\mathbb{R}^{M\times M}$—which
themselves can be interpreted as pseudo observations—and the global inducing
locations, $\mathbf{U}_{0}\in\mathbb{R}^{M\times D_{0}}$, which are used to
define $\\{\mathbf{U}^{\ell}\\}_{\ell=1}^{L}$ according to

$\mathbf{U}^{1}=\mathbf{U}^{0}\mathbf{W}_{1},\quad\mathbf{U}^{\ell}=\psi(\mathbf{U}^{\ell-1})\mathbf{W}^{\ell}\quad\ell=2,\ldots,L.$
(2)

Optimisation of $\phi$ is achieved through maximisation of the evidence lower
bound (ELBO):

$\mathcal{L}_{\text{ELBO}}(\phi;\mathcal{D})=\mathbb{E}_{q_{\phi}(\mathbf{W})}\left[\log
p(\mathbf{y}|\mathbf{W},\mathbf{X})\right]-\mathrm{KL}\left[{q_{\phi}(\mathbf{W})}||{p(\mathbf{W})}\right].$
(3)

We refer to this variational approximation as pseudo-observation variational
inference for BNNs (POVI-BNN).

Ober and Aitchison (2021) demonstrate the efficacy of POVI-BNNs relative to
mean-field Gaussian variational approximations for BNNs, achieving state-of-
the-art performance on a number of regression and classification experiments.
Their effectiveness is further demonstrated by Bui (2022), who shows that the
estimate of the marginal likelihood provided by the POVI-BNN approximation is
close to the true value, indicating the approximation is close to the true
posterior.

##  4 Amortising Inference in Bayesian Neural Networks

In this section, we build upon the work of Ober and Aitchison (2021) to
develop an effective method of performing amortised inference in BNNs.

Consider the same variational approximation described in Section 3.1, except
with diagonal precision matrices $\boldsymbol{\mathbf{\Lambda}}^{\ell}_{d}$
and inducing locations $\mathbf{U}$ replaced by training locations
$\mathbf{X}$, such that

$q\left(\mathbf{W}^{\ell}|\\{\mathbf{W}^{\ell^{\prime}}\\}_{\ell^{\prime}=1}^{\ell-1},\mathcal{D}\right)\propto\prod_{d=1}^{D_{\ell}}p\left(\mathbf{w}^{\ell}_{d}\right)\prod_{n=1}^{N}\underbrace{\mathcal{N}\left({v^{\ell}_{d,n}};{x^{\ell}_{d,n}},{{\sigma_{d,n}^{\ell}}^{2}}\right)}_{t^{\ell}_{d,n}(\mathbf{w}^{\ell}_{d})}$
(4)

where

$\mathbf{x}^{1}_{n}=\mathbf{W}^{1}\mathbf{x}_{n},\quad\mathbf{x}^{\ell}_{n}=\mathbf{W}^{\ell}\psi(\mathbf{x}^{\ell-1}_{n})\quad\forall\ell=2,\ldots,L.$
(5)

This form of variational approximation enables the use of per-datapoint
amortised inference. Specifically, rather than treating the variational
parameters
$\\{\\{\\{v^{\ell}_{d,n},\sigma^{\ell}_{d,n}\\}_{n=1}^{N}\\}_{d=1}^{D_{\ell}}\\}_{\ell=1}^{L}$
as free-form parameters to optimise, we obtain them by passing each datapoint
$(\mathbf{x}_{n},\mathbf{y}_{n})$ through an inference network at each layer:

$\mathbf{v}^{\ell}_{n},\
\log\sigma^{\ell}_{n}=g^{\ell}_{\phi}(\mathbf{x}_{n},\mathbf{y}_{n})\quad\forall\ell=1,\ldots,L.$
(6)

The variational parameters now become those of the inference networks,
$\\{g^{\ell}_{\phi}\\}^{L}_{\ell=1}$. These variational parameters can be
shared across datasets, and can be learned using a meta-dataset $\Xi$ with the
training objective:

$\phi^{*}=\operatorname*{arg\,max}_{\phi}\mathcal{L}(\\{\phi^{\ell}\\}_{\ell=1}^{L};\Xi);\quad\mathcal{L}(\\{\phi^{\ell}\\}_{\ell=1}^{L};\Xi)=\frac{1}{|\Xi|}\sum_{i=1}^{|\Xi|}\mathcal{L}_{\text{ELBO}}(\\{\phi^{\ell}\\}_{\ell=1}^{L};\mathcal{D}_{i})$
(7)

where $\mathcal{L}_{\text{ELBO}}$ is defined in Equation 3. We refer to our
method as amortised pseudo-observation variational inference for BNNs (A-POVI-
BNN).

An important limitation of this approach is that performing stochastic
optimisation through mini-batching datapoints is not possible, as to compute
the $q(\mathbf{W}|\mathcal{D})$ we require passing the entire dataset
$\mathcal{D}$ through the inference network. Nonetheless, this limitation is
only significant for large datasets—provided the entire dataset can be passed
through the network at once (i.e. in the small to medium-sized dataset regime
which we consider here) this is not an issue. At test time, we can obtain an
approximate posterior $q(\mathbf{W}|\mathcal{D}_{*})$ with a single pass of
the dataset through the inference networks.

fig:regression_experiment [A-POVI-BNN ($|\Xi|=1$)]  [A-POVI-BNN ($|\Xi|=100$)]
[A-MFVI-BNN ($|\Xi|=1$)]  [A-MFVI-BNN ($|\Xi|=100$)]

Figure 1: Posterior predictive distributions for A-POVI-BNN (top), A-MFVI-BNN
(middle), and ConvCNP (bottom) after training on a meta-dataset of size
$|\Xi|=1$ (left) and $|\Xi|=100$ (right). Data points are shown as black dots,
the predictive mean is shown as a dark blue line and the 95% confidence
interval is shown as shaded blue.

[ConvCNP ($|\Xi|=1$)]  [ConvCNP ($|\Xi|=100$)]

##  5 Results and Discussion

In this section, we evaluate the performance of our model in the meta-learning
setting. We consider a synthetic meta-dataset consisting of samples from a GP
with a squared-exponential (SE) kernel. Each dataset consists of between 10
and 20 training datapoints. We compare A-POVI-BNNs to amortised mean-field
variational inference for BNNs (A-MFVI-BNN), which we detail in Appendix B,
and a ConvCNP (Gordon et al., 2019). All NNs used in the amortised BNN
architectures (both the model and inference network) consist of two layers of
50 hidden units and ReLU activation functions. The ConvCNP implementation and
architecture is identical to that provided in https://github.com/cambridge-
mlg/convcnp/blob/master/convcnp_regression.ipynb.

Figures 1 to 1 compare the posterior predictive distributions of each method
on an unseen dataset drawn from a GP with the same hyperparameters as those in
the meta-dataset. We evaluate the predictions of each method trained on a
meta-dataset of size $|\Xi|=1$ and $|\Xi|=100$. We see that only A-POVI-BNN
obtains a sensible predictive posterior on the unseen dataset in both cases.
The ConvCNP performs poorly for $|\Xi|=1$, which is unsurprising given the
large number of model parameters increasing it susceptibility to overfitting.
The A-MFVI-BNN performs significantly better for $|\Xi|=100$, yet, in both
cases the quality of its predictive posterior is poor relative to both the
A-POVI-BNN and ConvCNP.

Despite these results being very preliminary, they are encouraging and suggest
that A-POVI-BNNs may provide a more effective alternative to NPs when the size
of meta-datasets are small. We intend to explore the effectiveness of A-POVI-
BNNs in more diverse settings, such as image completion, in future work.

## References

  * Antoniou et al. (2018) Antreas Antoniou, Harrison Edwards, and Amos Storkey.  How to train your MAML.  _arXiv preprint arXiv:1810.09502_ , 2018.
  * Bui (2022) Thang D Bui.  Biases in variational Bayesian neural networks.  2022\.
  * Finn et al. (2017) Chelsea Finn, Pieter Abbeel, and Sergey Levine.  Model-agnostic meta-learning for fast adaptation of deep networks.  In _International conference on machine learning_ , pages 1126–1135. PMLR, 2017.
  * Foong et al. (2020) Andrew Foong, Wessel Bruinsma, Jonathan Gordon, Yann Dubois, James Requeima, and Richard Turner.  Meta-learning stationary stochastic process prediction with convolutional neural processes.  _Advances in Neural Information Processing Systems_ , 33:8284–8295, 2020.
  * Garnelo et al. (2018a) Marta Garnelo, Dan Rosenbaum, Christopher Maddison, Tiago Ramalho, David Saxton, Murray Shanahan, Yee Whye Teh, Danilo Rezende, and SM Ali Eslami.  Conditional neural processes.  In _International Conference on Machine Learning_ , pages 1704–1713. PMLR, 2018a.
  * Garnelo et al. (2018b) Marta Garnelo, Jonathan Schwarz, Dan Rosenbaum, Fabio Viola, Danilo J Rezende, SM Eslami, and Yee Whye Teh.  Neural processes.  _arXiv preprint arXiv:1807.01622_ , 2018b.
  * Gordon et al. (2018) Jonathan Gordon, John Bronskill, Matthias Bauer, Sebastian Nowozin, and Richard E Turner.  Meta-learning probabilistic inference for prediction.  _arXiv preprint arXiv:1805.09921_ , 2018.
  * Gordon et al. (2019) Jonathan Gordon, Wessel P Bruinsma, Andrew YK Foong, James Requeima, Yann Dubois, and Richard E Turner.  Convolutional conditional neural processes.  _arXiv preprint arXiv:1910.13556_ , 2019.
  * Ober and Aitchison (2021) Sebastian W Ober and Laurence Aitchison.  Global inducing point variational posteriors for Bayesian neural networks and deep Gaussian processes.  In _International Conference on Machine Learning_ , pages 8248–8259. PMLR, 2021.
  * Requeima et al. (2019) James Requeima, Jonathan Gordon, John Bronskill, Sebastian Nowozin, and Richard E Turner.  Fast and flexible multi-task classification using conditional neural adaptive processes.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Volpp et al. (2020) Michael Volpp, Fabian Flürenbrock, Lukas Grossberger, Christian Daniel, and Gerhard Neumann.  Bayesian context aggregation for neural processes.  In _International Conference on Learning Representations_ , 2020.
  * Yoon et al. (2018) Jaesik Yoon, Taesup Kim, Ousmane Dia, Sungwoong Kim, Yoshua Bengio, and Sungjin Ahn.  Bayesian model-agnostic meta-learning.  _Advances in neural information processing systems_ , 31, 2018.

##  Appendix A Relationship to Neural Processes

The NP family (Garnelo et al., 2018a, b) describes the class of meta-learning
models which use neural networks to map from datasets to posterior predictive
distributions in a single forward pass. Broadly speaking, NPs fall into two
categories—conditional NPs (CNPs) and latent NPs (LNPs). In both cases, each
observation in the context dataset $\mathcal{D}_{c}$ is processed separately
using an encoder,
$\mathbf{r}_{c}=e_{\theta_{e}}(\mathbf{x}_{c},\mathbf{y}_{c})\
\forall(\mathbf{x}_{c},\mathbf{y}_{c})\in\mathcal{D}_{c}$, where $\theta_{e}$
parameterises a NN. The local encodings $\\{\mathbf{r}_{c}\\}_{c=1}^{C}$ are
combined into a single representation $\mathbf{r}(\mathcal{D}_{c})$ using a
permutation invariant aggregation function, such as summation. This encoding
is then passed to the model decoder $d_{\theta_{d}}$, where $\theta_{d}$ also
parameterises a NN, such that
$d_{\theta_{d}}(\cdot;\mathbf{r}(\mathcal{D}_{c}))$ maps from target inputs
$\mathbf{X}_{T}$ to predictive distributions. CNPs assume the posterior
predictive distribution $p(\mathbf{y}_{T}|\mathcal{D}_{C},\mathbf{X}_{T})$
factorises as

$p(\mathbf{y}_{T}|\mathcal{D}_{C},\mathbf{X}_{T})=\prod_{t=1}^{T}p_{\theta}(\mathbf{y}_{t}|d_{\theta_{d}}(\mathbf{x}_{t};\mathbf{r}(\mathcal{D}_{c}))).$
(8)

LNPs differ to CNPs through their use of a stochastic latent variable
$\mathbf{z}\sim p_{\theta}(\mathbf{z}|\mathbf{r}(\mathcal{D}_{c}))$, rather
than a fixed embedding:

$p_{\theta}(\mathbf{y}_{T}|\mathbf{X}_{T},\mathcal{D}_{C})=\int
p_{\theta}(\mathbf{z}|\mathbf{r}_{C})\prod_{t\in\mathcal{T}}p_{\theta}(\mathbf{y}_{t}|d_{\theta_{d}}(\mathbf{x}_{t};\mathbf{z}))d\mathbf{z}.$
(9)

Importantly, doing so enables dependencies to be maintained across target
locations inducing by marginalisation over $\mathbf{z}$.

###  A.1 Training Neural Processes

Members of the NP family are trained through minimising the empirical loss of
the meta-dataset $\boldsymbol{\mathbf{\Xi}}$:

$\theta^{*}=\operatorname*{arg\,min}_{\theta}\mathcal{L}(\theta;\Xi);\quad\mathcal{L}(\theta;\Xi):=\frac{1}{M}\sum_{m=1}^{M}\mathcal{L}(\theta;\mathcal{D}_{m}).$
(10)

There are a variety of ways in which the loss function
$\mathcal{L}(\theta;\mathcal{D}_{m})$ can be constructed. Focusing on LNPs,
the two predominant methods are variational inference (VI) and maximum
likelihood (ML). VI training employs the loss function

$\mathcal{L}(\theta;\mathcal{D})=-\mathbb{E}_{p_{\theta}(\mathbf{z}|\mathcal{D}_{c}\cup\mathcal{D}_{T})}\left[\log
p_{\theta}(\mathbf{y}_{T}|\mathbf{X}_{T},\mathbf{z})\right]+\mathrm{KL}\left[{p_{\theta}(\mathbf{z}|\mathcal{D}_{c}\cup\mathcal{D}_{T})}||{p(\mathbf{z}|\mathcal{D}_{c})}\right]$
(11)

where $\mathcal{D}_{c}$ and $\mathcal{D}_{T}$ are two random non-overlapping
partitions of the dataset $\mathcal{D}$. This objective is motivated by the
interpretation of $p_{\theta}(\mathbf{z}|\mathcal{D})$ as an approximation to
the posterior distribution over $\mathbf{z}$ given $\mathcal{D}$. If
$p_{\theta}(\mathbf{z}|\mathcal{D})$ is a good approximation to the posterior,
then this objective becomes equivalent to the standard VI objective (ELBO):

$\mathcal{L}(\theta;\mathcal{D})=-\mathbb{E}_{p_{\theta}(\mathbf{z}|\mathcal{D})}\left[\log
p_{\theta}(\mathbf{y}_{T}|\mathbf{X}_{T},\mathbf{z})\right]+\mathrm{KL}\left[{p_{\theta}(\mathbf{z}|\mathcal{D})}||{p(\mathbf{z})}\right].$
(12)

Nonetheless, this objective potentially wastes modelling power on obtaining
consistent posterior approximations, rather than good posterior predictive
distributions. Instead, Foong et al. (2020) argue that the following ML
objective should be used instead:

$\mathcal{L}(\theta;\mathcal{D})=-\log
p_{\theta}(\mathbf{y}_{T}|\mathbf{X}_{T},\mathcal{D}_{c}).$ (13)

An important difference between the two objectives is that in the VI
objective, the parameters of the embedding function can be interpreted as
variational parameters, whereas in this ML objective all parameters are model
parameters.

###  A.2 Amortised Inference in BNNs as Neural Processes

There are distinct similarities between the model introduced in the previous
section and members of the latent neural process (LNP) family. In both cases
we

  1. 1.

compute per-datapoint embeddings of the dataset;

  2. 2.

perform a permutation invariant aggregation of embeddings to get a single
embedding for the entire dataset;

  3. 3.

use this embedding to define a distribution over a latent variable;

  4. 4.

sample this latent variable to define mappings from input location to
predictive distributions.

The only principal difference between LNPs and our model is the use of shared
model parameters in the decoder of LNPs. In our method, the latent variable is
the parameters of the decoder. Alternative models such as CNAPs (Requeima et
al., 2019) and ML-PIP (Gordon et al., 2018) are similar in this manner, except
that the per-dataset latent variables are only a small subset of the decoder
parameters.

Given these similarities, it is natural to consider training our model in an
identical manner to LNPs—in particular, using the ML objective. In doing so,
the parameters of the inference network no longer act as variational
parameters but rather become part of the model, and so increase the propensity
of the model to overfit. We leave exploration of this to future work.

##  Appendix B Amortised Mean-Field Variational Inference for BNNs

A popular variational approximation for BNNs is the mean-field Gaussian
approximation, given by

$q(\mathbf{W})=\prod_{\ell=1}^{L}\prod_{i=1}^{|\mathbf{W}^{\ell}|}\mathcal{N}\left({w_{i}};{\mu_{i}},{\sigma_{i}^{2}}\right)$
(14)

where $|\mathbf{W}^{\ell}|$ denotes the number of elements of weight matrix
$\mathbf{W}^{\ell}$. We cannot directly apply per-datapoint amortisation to
this form of the mean-field approximation as there are not per-datapoint
variational parameters. However, we can obtain per-datapoint variational
parameters by instead constructing the variational approximation as

$q(\mathbf{W})\propto p(\mathbf{W})\prod_{n=1}^{N}t_{n}(\mathbf{W})$ (15)

where each approximate likelihood $t_{n}(\mathbf{W})$ is an unnormalised mean-
field Gaussian distribution over $\mathbf{W}$. Provided $p(\mathbf{W})$ is
also mean-field Gaussian (which is typical), then $q(\mathbf{W})$ will be
mean-field Gaussian. The variational parameters of each $t_{n}(\mathbf{W})$
can then be obtained by passing each data pair
$(\mathbf{x}_{n},\mathbf{y}_{n})$ through an inference network, achieving per-
datapoint amortisation.