# Differential equations, splines and Gaussian processes

Lars Lau Raket

###### Abstract

We explore the connections between Green’s functions for certain differential
equations, covariance functions for Gaussian processes, and the smoothing
splines problem. Conventionally, the smoothing spline problem is considered in
a setting of reproducing kernel Hilbert spaces, but here we present a more
direct approach. With this approach, some choices that are implicit in the
reproducing kernel Hilbert space setting stand out, one example being choice
of boundary conditions and more elaborate shape restrictions.

The paper first explores the Laplace operator and the Poisson equation and
studies the corresponding Green’s functions under various boundary conditions
and constraints. Explicit functional forms are derived in a range of examples.
These examples include several novel forms of the Green’s function that, to
the author’s knowledge, have not previously been presented. Next we present a
smoothing spline problem where we penalize the integrated squared derivative
of the function to be estimated. We then show how the solution can be
explicitly computed using the Green’s function for the Laplace operator. In
the last part of the paper, we explore the connection between Gaussian
processes and differential equations, and show how the Laplace operator is
related to Brownian processes and how processes that arise due to boundary
conditions and shape constraints can be viewed as conditional Gaussian
processes. The presented connection between Green’s functions for the Laplace
operator and covariance functions for Brownian processes allows us to
introduce several new novel Brownian processes with specific behaviors.
Finally, we consider the connection between Gaussian process priors and
smoothing splines.

This paper was originally developed as part of the teaching material for a
graduate course in functional data analysis held in 2015 at Department of
Mathematical Sciences, University of Copenhagen.

> “Polynomials are wonderful even after they are cut into pieces, but the
> cutting must be done with care.”
>

##  1 Introduction

Splines are piecewise polynomials that are required to have certain smoothness
properties. The introductory quote, by one of the early pioneers in the
development of splines, encapsulates much of the feeling of working with
splines. Splines are truly wonderful and are applicable to a large body of
problems, but to understand the details of how they arise and how they are
connected to various branches of mathematics and statistics require care.

During the last half century, smoothing splines have become a fundamental tool
for estimating continuous functions from discretely observed data. In its
typical formulation, the smoothing spline problem involves finding the
function $\theta$ that minimizes a discrete Gaussian likelihood term that
measures the deviation from observed data plus a penalty term that measures
the roughness of the function. The conventional choice of roughness measure is
in terms of the integral of its squared second derivative. Let
$y_{1},\dots,y_{m}$ denote observed data with associated continuous covariate
values $t_{1},\dots,t_{m}$. In this note, we will think of the $t$-variables
as representing time and denote them as such, but in practice they can
represent any continuous covariate. The classical smoothing spline estimator
is defined as the minimizer

$\displaystyle\hat{\theta}=\operatorname*{arg\,min}_{\theta}\sum_{i=1}^{m}(y_{i}-\theta(t_{i}))^{2}+\lambda\int_{T}|\theta^{\prime\prime}(t)|^{2}\,\mathrm{d}t.$
(1)

It was shown by [1] that the function that minimizes the roughness term
$\int_{T}|\theta^{\prime\prime}(t)|^{2}\,\mathrm{d}t$ under the constraint of
interpolating $y_{1},\dots,y_{m}$ at time points $t_{1}<\dots<t_{m}$ is a
cubic spline with the so-called natural boundary contions
$\theta^{\prime\prime}(t_{1})=\theta^{\prime\prime}(t_{m})=0$. It easily
follows that the solution to (1) must be a natural cubic smoothing spline, in
fact, one can replace the quadratic data term with any reasonable pointwise
data term and the solution will still be a natural cubic smoothing spline. The
argument is simple: Given any candidate solution $\tilde{\theta}$, choose
$\hat{\theta}$ as the natural cubic spline that interpolates $\tilde{\theta}$
in the points $t_{1},\dots,t_{m}$. The pointwise data term will produce the
same measure of deviation for $\tilde{\theta}$ and $\hat{\theta}$, but the
roughness term for the interpolating cubic spline $\hat{\theta}$ will always
be less than or equal to the roughness term for $\tilde{\theta}$.

Since Holladay’s results, much work has been devoted to studying the smoothing
spline problem (1) and its generalizations. The theory and methodology for
smoothing splines matured greatly when Grace Wahba and collaborators
formulated and solved the problem in the setting of reproducing kernel Hilbert
spaces (RKHS) [2, 8]. While the RKHS approach is powerful, it is not easily
accessible in terms of the required level of mathematical skill and the
generality of the approach may lead to underappreciation of some of the finer
details of the considered problems. The avid reader may for example have
noticed how the smoothing spine problem (1) did not mention the space over
which the problem was considered. The classical RKHS approach will simply
choose the function space to minimize over as the least restrictive space, but
from the perspective of a statistician that choice may not align with the
knowledge of the problem at hand. For simplicity, the rest of this paper will
explore a slightly simpler smoothing spline problem, namely

$\displaystyle\hat{\theta}=\operatorname*{arg\,min}_{\theta}\sum_{i=1}^{m}(y_{i}-\theta(t_{i}))^{2}+\lambda\int_{T}|\theta^{\prime}(t)|^{2}\,\mathrm{d}t,$
(2)

where the roughness penalty is replaced by the integrated squared derivative
of the function. The solution to this problem will in general be a piecewise
linear function with knots in $t_{1},\dots,t_{m}$ as shown by [7], but we show
here how we may get complex behavior of the solution by imposing different
criteria on the space we are minimizing over. We even show that by imposing
sufficiently strict criteria, we may get pathological behavior of the spines,
for example such that they become infinitely smooth polynomials. Most of the
approaches and connections presented here will seamlessly transfer to the
classical smoothing spline problem (1) and the general L-spline problem. While
requiring some determination, explicit forms of various Green’s functions
defining the solutions to the classical spline problem can be derived (see
e.g. [6], rytgaard2016statistical.)

The approach taken here is related to and influenced by a number of previous
works. One of them is the work by [3] where the connection between Green’s
functions and smoothing splines is used to develop a framework for approximate
inference for mixed smoothing spline models. Along the way, [3] presents a
result [3, Theorem 1] that gives an explicit (although quite complex) form for
Green’s function associated with constant-coefficient linear differential
operators under various boundary conditions. This approach was extended to
data on multidimensional domains such as $[0,1]^{d}$ by [5]. In this paper,
the connection between Green’s function for various differential operators and
covariance functions for known Gaussian processes such as the Brownian sheet,
the tied-down Brownian sheet and the Matérn processes are explored in
examples. Finally, parts of the work presented here has been extended to the
original smoothing spline problem (1) in [6]. In this work, [6] presents
explicit forms for the Green’s functions under all combinations of Dirichlet
and Neumann boundary conditions along with a number of other extensions
involving shape restrictions and alternative likelihood terms.

##  2 Green’s functions and the Poisson equation

In this note we will consider the Poisson equation as our driving example. The
equation is given by

$\displaystyle-\partial_{t}^{2}f(t)=h(t)$ (3)

where $-\partial_{t}^{2}$ is the _Laplace operator_
$-\partial_{t}^{2}f=-f^{\prime\prime}$, $f\,:\,[0,1]\rightarrow\mathbb{R}$ is
twice differentiable in a suitable sense and $h\in L^{2}([0,1])$. This
differential equation is not too difficult to solve, an obvious idea is to
write the solutions $f$ as double integrals of $h$

$f(t)=-\int_{a_{1}}^{t}\int^{s_{1}}_{a_{2}}h(s_{2})\,\mathrm{d}s_{2}\,\mathrm{d}s_{1}.$

Here we are however not too interested in the specific solutions to the
differential equation, but rather in the general problem of inverting the
differential operator $-\partial_{t}^{2}$. Suppose an inverse operator
$\mathscr{G}$ exists, then it should hold that

$-\partial_{t}^{2}\mathscr{G}h(t)=h(t)$

and thus this operator provides a general method for solving the differential
equation (3).

###  2.1 Green’s functions

When a differential operator $\mathscr{L}$ is invertible its inverse
$\mathscr{G}$ will be an integral operator with kernel $\mathcal{G}$,

$\displaystyle\mathscr{G}f(t)=\int_{0}^{1}\mathcal{G}(s,t)f(s)\,\mathrm{d}s.$
(4)

The kernel $\mathcal{G}$ is referred to as the _Green’s function_ for
$\mathscr{L}$. Green’s functions have two important properties that we will
use extensively, but not prove here: firstly, Green’s functions are symmetric
$\mathcal{G}(s,t)=\mathcal{G}(t,s)$; and secondly, the Green’s function as a
function of the individual coordinates $\mathcal{G}(s,\,\cdot\,)$ typically
obey boundary conditions on the space of function under consideration. In the
cases we will consider here, the kernel will be a well-behaved function, but
it will sometimes have a more complex structure. In general the Green’s
function will be a _distribution_ or _generalized function_. The most basic
distribution is the Dirac delta at $t$, $\delta_{t}$, for which it holds that

$\displaystyle\int_{0}^{1}\delta_{t}(s)f(s)\,\mathrm{d}s=f(t).$ (5)

It is important to note here that the above definition of the Dirac delta
function only defines it as the point evaluation functional in combination
with an inner product. Thus, it does not always have the standard form of the
Dirac delta where one would think of it as a function with $\delta_{t}(s)=0$
when $s\neq t$ with an infinitely large value at $t$ corresponding to a point
mass.

######  Example 2.1.

Using distributions, it is possible to define derivatives of discontinuous
functions in a distributional sense. For example, suppose that the constant
function $1(t)=1$ is in the function space we are considering. Then for
$a\in[0,1]$

$\int_{0}^{t}\delta_{a}(s)\,\mathrm{d}s=\boldsymbol{1}_{[a,1]}(t)$

where $\boldsymbol{1}_{A}(t)$ is the indicator function that is 1 when $t\in
A$ and $0$ otherwise. Thus, in the distributional sense

$\partial_{t}\boldsymbol{1}_{[a,1]}(t)=\delta_{a}(t).$

$\circ$

By applying the differential operator $\mathscr{L}$ on both sides equation
(4), we get the important property of Green’s functions that

$f(t)=\int_{0}^{1}\mathscr{L}\mathcal{G}(s,t)f(s)\,\mathrm{d}s,$

which gives rise to an alternative definition of Green’s functions: a Green’s
function $\mathcal{G}$ for $\mathscr{L}$ is a distribution for which

$\displaystyle(\mathscr{L}\mathcal{G}(s,\,\cdot\,))(t)=\delta_{t}(s)$ (6)

or equivalently

$\mathscr{L}\mathcal{G}(s,\,\cdot\,)=\delta_{s}.$

Again, we note that this definition is in the distributional sense, thus,
depending on the function space we are considering, the representation of
$\delta_{t}$ may vary.

######  Example 2.2.

From Example 2.1 we see that $\mathcal{G}(s,t)=\boldsymbol{1}_{[s,1]}(t)$ is a
Green’s function for $\partial_{t}$. $\circ$

Since $\mathscr{L}$ is a linear operator on a function space $\mathscr{H}$, it
is invertible if and only if its kernel is trivial

$\ker(\mathscr{L})=\\{f\in\mathscr{H}\,:\,\mathscr{L}f=o\\}=\\{o\\}$

where $o(t)=0$ is the zero function.

######  Example 2.3.

The kernel for the Laplace operator on the space of twice continuously
differentiable functions $\mathscr{H}=\mathcal{C}^{2}([0,1])$ are the
functions $f$ for which

$-\partial_{t}^{2}f(t)=0\qquad\textrm{for all }t.$

Thus, the kernel consists of all affine functions

$\ker(-\partial_{t}^{2})=\\{f\in\mathscr{H}\,:\,f(t)=at+b\\}$

and the Laplace operator is not invertible on $\mathscr{H}$. $\circ$

As the previous example showed, we may need to impose extra conditions on the
function space $\mathscr{H}$ in order to identify a unique Green’s function
for $\mathscr{L}$. We may however be able to say something about its general
form in certain cases.

######  Example 2.4.

Any Green’s function for the Laplace operator on a subspace $\mathscr{H}$ of
$\mathcal{C}^{2}([0,1])$ must obey

$-\partial_{t}^{2}\mathcal{G}(s,t)=\delta_{s}(t).$

Assuming there exists $f\in\mathscr{H}$ such that $-\partial_{t}^{2}f=1$, we
can use Example 2.1 to write

$-\partial_{t}^{2}\mathcal{G}(s,t)=\partial_{t}\boldsymbol{1}_{[s,1]}(t).$

This means means that

$\partial_{t}\mathcal{G}(s,t)=-\boldsymbol{1}_{[s,1]}(t)+a_{1}(s),$

so that

$\displaystyle\mathcal{G}(s,t)=-(t-s)\boldsymbol{1}_{[s,1]}(t)+a_{1}(s)t+a_{0}(s).$
(7)

Because of the symmetry of $\mathcal{G}$, $a_{1}$ and $a_{0}$ will generally
be first-order affine functions. Again we note that this form is only in the
distributional sense, and the exact functional form may differ from this.
$\circ$

###  2.2 Boundary conditions

The typical solution to the problem of non-trivial kernels is to consider
linear subspaces of functions subject to boundary conditions. The two most
important types of boundary conditions are the _homogeneous Dirichlet boundary
conditions_

$f(0)=f(1)=0$

and the _homogeneous Neumann boundary conditions_

$f^{\prime}(0)=f^{\prime}(1).$

Given a function space, both of these types of boundary conditions will
produce linear subspaces of the original space. One can also consider higher-
order Neumann boundary conditions

$f^{(k)}(0)=f^{(k)}(1)=0,$

and combinations of boundary conditions, both in the form of requiring
multiple homogeneous boundary conditions to hold simultaneously, or by
specifying linear combinations of boundary conditions that should be zero.

######  Example 2.5.

Consider the Laplace operator $-\partial_{t}^{2}$ on the subspace
$\mathscr{H}$ of $\mathcal{C}^{2}([0,1])$ that arise from imposing homogeneous
Dirichlet boundary conditions. Since the kernel of the Laplace operator on
$C^{2}([0,1])$ consists of all affine functions, the kernel of
$-\partial_{t}^{2}$ on $\mathscr{H}$ is the affine functions for which the
boundary conditions hold. Thus $\ker(-\partial_{t}^{2})=\\{o\\}$ on
$\mathscr{H}$ and the Laplace operator is invertible.

The Green’s function is in a space of distributions associated to
$\mathscr{H}$, which means that $\mathcal{G}(s,\,\cdot\,)$ must obey the
boundary conditions on $\mathscr{H}$. Since the conditions of Example 2.4 hold
(see Exercise 1), the Green’s function is given by equation (7). By the
boundary contions we have that

$0=\mathcal{G}(s,0)=a_{0}(s)$

and

$0=\mathcal{G}(s,1)=-(1-s)+a_{1}(s)$

Thus,

$\displaystyle\mathcal{G}(s,t)=-(t-s)\boldsymbol{1}_{[s,1]}(t)+(1-s)t=s\wedge
t-st$ (8)

where $\wedge$ denotes the minimum. $\circ$

######  Example 2.6.

Consider again the Laplace operator, but this time on the subspace
$\mathscr{H}$ of $f\in\mathcal{C}^{2}([0,1])$ where $f$ obeys homogeneous
Neumann boundary conditions. We see that the kernel of $-\partial_{t}^{2}$
consists of all constant functions, and thus we need an extra condition to
ensure invertibility. A simple choice is to impose the restriction that
$f(0)=0$, which will make the kernel trivial. $\circ$

######  Example 2.7.

We continue to consider the Laplace operator. This time on the subspace
$\mathscr{H}$ of functions $f$ subject to the mixed boundary conditions
$f(0)=f^{\prime}(1)=0$. Clearly the kernel of the Laplace operator is trivial,
so a unique Green’s function exists. Since there exists $f\in\mathscr{H}$ such
that $-\partial_{t}^{2}f(t)=1$ (see Exercise 1), we can use Example 2.4 to
give the form of the Green’s function. Imposing the boundary conditions on the
general form of the Green’s function (7) we get that

$0=\mathcal{G}(s,0)=a_{0}(s)$

and

$0=\partial_{t}\mathcal{G}(s,1)=-1+a_{1}(s).$

Thus, the Green’s function is

$\displaystyle\mathcal{G}(s,t)=-(t-s)\boldsymbol{1}_{[s,1]}(t)+t=s\wedge t.$
(9)

$\circ$

###  2.3 Finding Green’s functions using Fourier series

Recall that the _Fourier series_ of a function
$f\,:\,[0,1]\rightarrow\mathbb{R}$ is given by

$\displaystyle s_{f}(t)=\frac{1}{2}a_{0}+\sum_{i=1}^{\infty}a_{i}\cos(2i\pi
t)+\sum_{i=1}^{\infty}b_{i}\sin(2i\pi t)$ (10)

where

$a_{i}=2\int_{0}^{1}f(t)\cos(2i\pi t)\,\mathrm{d}t\qquad\text{and}\qquad
b_{i}=2\int_{0}^{1}f(t)\sin(2i\pi t)\,\mathrm{d}t.$

One idea for finding Green’s functions is to use a Fourier series
representation for them, and use the differential operator and boundary
conditions to identify the coefficients. We have that

$\displaystyle-\partial_{t}^{2}s_{f}(t)=4\pi^{2}\sum_{i=1}^{\infty}a_{i}i^{2}\cos(2i\pi
t)+4\pi^{2}\sum_{i=1}^{\infty}b_{i}i^{2}\sin(2i\pi t).$ (11)

The sine and cosine functions make up an orthogonal basis for $L^{2}([0,1])$
and we have that

$\int_{0}^{1}\cos(2i\pi t)\cos(2j\pi t)\,\mathrm{d}t=\begin{cases}1&\textrm{
if }i=j=0\\\ 1/2&\text{ if }i=j\neq 0\\\ 0&\text{ otherwise}\end{cases},$
$\int_{0}^{1}\sin(2i\pi t)\sin(2j\pi t)\,\mathrm{d}t=\begin{cases}1/2&\text{
if }i=j\neq 0\\\ 0&\text{ otherwise}\end{cases},$

and that the inner product between any sine and cosine function is zero. To
identify the coefficients for the solution $s_{f}$ of the differential
equation (3)

$-\partial_{t}^{2}s_{f}(t)=h(t)$

we can multiply both sides by respectively $\cos(j\pi t)$ and $\sin(j\pi t)$
and integrate. In the first case, we get

$4\pi^{2}a_{j}j^{2}\int_{0}^{1}\cos(2j\pi
t)^{2}\,\mathrm{d}t=\int_{0}^{1}\cos(2j\pi t)h(t)\,\mathrm{d}t.$

Solving this in $a_{j}$ gives

$a_{j}=\int_{0}^{1}\frac{1}{2j^{2}\pi^{2}}\cos(j\pi t)h(t)\,\mathrm{d}t,$

and a similar result holds for $b_{j}$. By rearranging, the solution $s_{f}$
thus has the form

$\displaystyle
s_{f}(t)=\frac{1}{2}a_{0}+\int_{0}^{1}\sum_{i=1}^{\infty}\frac{1}{2i^{2}\pi^{2}}(\cos(2i\pi
s)\cos(2i\pi t)+\sin(2i\pi s)\sin(2i\pi t))h(s)\,\mathrm{d}s.$ (12)

Comparing this to the definition of the Green’s function (4) it is not too
hard to see why this may be useful for finding Green’s functions.

######  Example 2.8 (Constant term).

From the definition of $a_{0}$, we see that $a_{0}=0$ when

$\int_{0}^{1}f(t)\,\mathrm{d}t=0.$

Thus, if we consider the linear subspace $\mathscr{H}$ of
$f\in\mathcal{C}^{2}([0,1])$ for which this holds, any Green’s function has
the form

$\displaystyle\mathcal{G}(s,t)=\sum_{i=1}^{\infty}\frac{1}{2i^{2}\pi^{2}}(\cos(2i\pi
s)\cos(2i\pi t)+\sin(2i\pi s)\sin(2i\pi t)).$ (13)

On the other hand, assume that $a_{0}\neq 0$. Using integration by parts
twice, the first term can be written on the form

$\displaystyle\frac{1}{2}a_{0}$ $\displaystyle=\int_{0}^{1}f(s)\,\mathrm{d}s$
$\displaystyle=[f(s)(2s+c_{1})]_{0}^{1}-[f^{\prime}(s)(t^{2}+c_{1}s+c_{0})]_{0}^{1}+\int_{0}^{1}(s^{2}+c_{1}s+c_{0})f^{\prime\prime}(t)\,\mathrm{d}s$
$\displaystyle=f(1)(2+c_{1})-f(0)c_{1}-f^{\prime}(1)(1+c_{1}+c_{0})+f^{\prime}(0)c_{0}$
$\displaystyle\quad-\int_{0}^{1}(s^{2}+c_{1}s+c_{0})h(s)\,\mathrm{d}s$

where it has been used that $f^{\prime\prime}(s)=-h(s)$. Thus, if we can make
the first four terms disappear, the Green’s function has the form

$\mathcal{G}(s,t)=\sum_{i=1}^{\infty}\frac{1}{2i^{2}\pi^{2}}(\cos(i\pi
s)\cos(i\pi t)+\sin(i\pi s)\sin(i\pi t))-(s^{2}+c_{1}(t)s+c_{0}(t)).$

$\circ$

We note once again that the found forms of the Green’s functions are only to
be understood in distributional sense, that is, they give the right results
under the integration. But as a function of $t$ the given form of a Green’s
function will typically not belong to $\mathscr{H}$. We end this section with
an example that demonstrates how to find the closed-form solution of the
cosine-series in the general form of the Green’s function (13).

######  Example 2.9 (Cosine series).

Consider the cosine series

$c(s,t)=\sum_{i=1}^{\infty}\frac{1}{2i^{2}\pi^{2}}\cos(2i\pi s)\cos(2i\pi t).$

From the product-to-sum formula, we get that

$\displaystyle\cos(2i\pi s)\cos(2i\pi
t)=\frac{1}{2}(\cos(2i\pi(s-t))+\cos(2i\pi(s+t))).$ (14)

We can thus rewrite the representation

$c(s,t)=\sum_{i=1}^{\infty}\frac{1}{4i^{2}\pi^{2}}\cos(2i\pi(s-t))+\sum_{i=1}^{\infty}\frac{1}{4i^{2}\pi^{2}}\cos(2i\pi(s+t)).$

Since the cosine function is even $\cos(2i\pi(s-t))=\cos(2i\pi|s-t|)$. On the
other hand, since the cosine function has period $2\pi$,
$\cos(2i\pi(s+t))=\cos(2i\pi(s+t-\boldsymbol{1}_{[1,2]}(s+t)))$.

Using some complex analysis and properties of dilogarithms, one can show that
for $s\in[0,1]$

$\displaystyle\sum_{i=1}^{\infty}\frac{1}{4i^{2}\pi^{2}}\cos(2i\pi
s)=\frac{1}{4}(s(s-1)+1/6).$ (15)

Putting it all together, we get that

$\displaystyle c(s,t)=\frac{1}{12}+\frac{1}{2}(s^{2}+t^{2}-s\wedge
t)-\frac{1}{2}\boldsymbol{1}_{[1,2]}(s+t)(s+t-1)$ (16)

which can also be written on the form

$c(s,t)=\begin{cases}\frac{1}{12}+\frac{1}{2}(s^{2}+t^{2}-s\wedge
t)&\textrm{for }0\leq s\leq 1-t\\\
\frac{1}{12}+\frac{1}{2}(s(s-1)+t(t-1)+1-s\wedge t)&\textrm{for }1-t<s\leq
1\end{cases}.$

$\circ$

###  2.4 Fourier bases and linear subspaces

The approach for finding Green’s functions for the Poisson equation presented
in the previous section is very general. To find the Green’s function
$\mathcal{G}$ for $-\partial_{t}^{2}$ on some subspace $\mathscr{H}$, we may
formulate the constraints of the subspace in terms of the Fourier basis
functions.

######  Example 2.10.

Consider the space $\mathscr{H}$ of $f\in\mathcal{C}^{2}([0,1])$ subject to
homogeneous Dirichlet boundary conditions. Imposing these restrictions on the
representation (10) gives that

$0=s_{f}(0)=\frac{1}{2}a_{0}+\sum_{i=1}^{\infty}a_{i}$

and

$0=s_{f}(1)=\frac{1}{2}a_{0}+\sum_{i=1}^{\infty}a_{i}.$

Thus, the representation of $f$ becomes

$s_{f}(t)=\sum_{i=1}^{\infty}a_{i}(\cos(2i\pi
t)-1)+\sum_{i=1}^{\infty}b_{i}\sin(2i\pi t),$

and the Green’s function can be written as

$\mathcal{G}(s,t)=\sum_{i=1}^{\infty}\frac{1}{2i^{2}\pi^{2}}((\cos(2i\pi
s)-1)(\cos(2i\pi t)-1)+\sin(2i\pi s)\sin(2i\pi t)).$

From Example 2.5 we know that this series must equal $s\wedge t-st$, but we
can also show that explicitly. The product-to-sum formulas for the sine
functions give that

$\displaystyle\sin(2i\pi s)\sin(i\pi
t)=\frac{1}{2}(\cos(2i\pi(s-t))-\cos(2i\pi(s+t))).$ (17)

The cosine term $((\cos(2i\pi s)-1)(\cos(2i\pi t)-1)$ equals

$\cos(2i\pi s)\cos(2i\pi t)+1-\cos(2i\pi s)-\cos(2i\pi t),$

so using the cosine product-to-sum formula (14), we get that the Green’s
function has the form

$\mathcal{G}(s,t)=\sum_{i=1}^{\infty}\frac{1}{2i^{2}\pi^{2}}(\cos(2i\pi(s-t))+1-\cos(2i\pi
s)-\cos(2i\pi t)).$

which, using the cosine-series result (15) gives the explicit form

$\mathcal{G}(s,t)=\frac{1}{2}\bigg{(}(|s-t|(|s-t|-1)+1/6)+\frac{1}{6}-(s(s-1)+1/6)-(t(t-1)+1/6)\bigg{)}.$

which can easily be manipulated to give the form

$\mathcal{G}(s,t)=s\wedge t-st.$

######  Example 2.11 (Zero coefficients).

Consider (the closure of) a linear subspace $\mathscr{H}$ of $C^{2k}$ spanned
by Fourier basis functions subject to the constraints

$a_{i_{1}}=0,\dots,a_{i_{\ell}}=0.$

Suppose that $-\partial_{t}^{2}$ is invertible on on $\mathscr{H}$. From (11),
we see that any admissible function $h$ on the left-hand side of the
differential equation (3) will also have zero coefficients for the
corresponding basis functions. Thus, due to the orthogonality all the
corresponding terms will disappear from the series inside the integral in
(12). $\circ$

######  Example 2.12 (Linear combinations of coefficients).

Consider (the closure of) a subspace $\mathscr{H}$ of $C^{2k}$ spanned by the
Fourier basis functions subject to a linear constraint on the coefficients

$\sum_{i=0}^{\infty}c_{i}a_{i}=0.$

where $c_{i}\in\mathbb{R}$. Suppose that $-\partial_{t}^{2}$ is invertible on
on $\mathscr{H}$, and for simplicity assume that $c_{0}=0$, so we have no
constraint on $a_{0}$. If we consider the general form of the Green’s function
(13), we see that it will generally not be in $\mathscr{H}$, since this would
require that

$\sum_{i=0}^{\infty}c_{i}\frac{1}{i^{2}}=0.$

To correct for this, consider the compensating term

$\mathcal{G}_{c}(s,t)=\sum_{i=0}^{\infty}c_{i}\frac{1}{2i^{2}\pi^{2}}\cos(2i\pi
s)\cos(2i\pi t).$
$-\partial_{t}^{2}\mathcal{G}_{c}(s,t)=\sum_{i=0}^{\infty}c_{i}\cos(2i\pi
s)\cos(2i\pi t).$

Using the Fourier series representation $s_{f}$ of a function
$f\in\mathscr{H}$, we see that

$\int_{0}^{1}-\partial_{t}^{2}\mathcal{G}_{c}(s,t)s_{f}(s)\,\mathrm{d}s=\frac{1}{2}\sum_{i=1}^{\infty}c_{i}a_{i}=0.$

Thus, we can subtract $\mathcal{G}_{c}$ from the general form of the Green’s
function without it having any effect in the distributional sense. $\circ$

###  2.5 Examples of Green’s functions

######  Example 2.13 (Balanced, periodic functions).

Consider the subspace $\mathscr{H}$ of functions $f\in\mathcal{C}^{2}([0,1])$
for which $\int_{0}^{1}f(t)\,\mathrm{d}t=0$. The kernel of the Laplace
operator on $\mathscr{H}$ is

$\ker(-\partial_{t}^{2})=\\{f\in\mathscr{H}:f(t)=a(t-1/2)\\},$

so we need one additional constraint. The least restricting constraint is that
$f(0)=f(1)$, so we restrict ourselves to those functions. The integral
constraint is equivalent to $a_{0}=0$ in the Fourier expansion (10), and all
the Fourier basis functions obey the constraint of periodicity (that the
endpoints are the same). Thus, the Green’s function is on the form

$\displaystyle\mathcal{G}(s,t)=\sum_{i=1}^{\infty}\frac{1}{2i^{2}\pi^{2}}(\cos(2i\pi
s)\cos(2i\pi t)+\sin(2i\pi s)\sin(2i\pi t)),$ (18)

which we can rewrite using the sine and cosine product-to-sum formulas as we
did in Example 2.10. This gives the the representation

$\mathcal{G}(s,t)=\sum_{i=1}^{\infty}\frac{1}{2i^{2}\pi^{2}}\cos(2i\pi(s-t)).$

Using the cosine series (15) we arrive at the closed form solution

$\mathcal{G}(s,t)=\frac{1}{2}|s-t|^{2}-\frac{1}{2}|s-t|+\frac{1}{12}.$

We note that this Green’s function is equivalent to the covariance function of
the periodic Brownian motion constructed by [4] (mumford2010pattern; Section
6.6.1). $\circ$

######  Example 2.14 (Odd functions).

Consider the space of functions $\mathscr{H}\subset\mathcal{C}^{2}([0,1])$ for
which $f(t)=-f(1-t)$ for $t\in[0,1/2]$. Again we see that the kernel for the
Laplace operator consists of the lines that are zero at $1/2$. Under the
oddness constraint, the additional constraints $f(0)=f(1)$, $f(0)=0$ and
$f(1)=0$ are all equivalent. Let $\mathscr{H}$ be the restriction to these
functions.

The oddness implies that $\int_{0}^{1}f(t)\,\mathrm{d}t=0$, so the Green’s
function must have a similar form to the previously found Green’s function
(18). The terms that are left out are exactly all the cosine functions, since
they are even. Using the product-to-sum formula for sine, the Green’s function
is found to be

$\mathcal{G}(s,t)=\sum_{i=1}^{\infty}\frac{1}{4i^{2}\pi^{2}}(\cos(2i\pi(s-t))-\cos(2i\pi(s+t))).$

We see that this must equal

$\displaystyle\mathcal{G}(s,t)$
$\displaystyle=\frac{1}{4}\bigg{(}|s-t|(|s-t|-1)-(s+t-\boldsymbol{1}_{[1,2]}(s+t))(s+t-\boldsymbol{1}_{[1,2]}(s+t)-1)\bigg{)}$
$\displaystyle=\begin{cases}\frac{1}{2}s\wedge t-st&\text{for }s+t\leq 1\\\
s\wedge t+\frac{1}{2}s\vee t-st-\frac{1}{2}&\text{for }s+t>1\end{cases}$

$\circ$

######  Example 2.15 (Mixed boundary conditions).

In this example we will consider the Laplace operator $-\partial_{t}^{2}$ on
function spaces $\mathscr{H}$ consisting of functions
$f\in\mathcal{C}^{2}([0,1])$ satisfying the mixed boundary conditions
$f(0)=f^{\prime}(1)=0$ and possibly other conditions. In the case that no
other conditions are available, we have already seen that the Green’s function
is

$\mathcal{G}(s,t)=s\wedge t.$

If we add the constraint that

$\int_{0}^{1}f(t)\,\mathrm{d}t=0,$

we see that the previous Green’s function no longer works, since

$\int_{0}^{1}\mathcal{G}(s,t)\,\mathrm{d}t=-\frac{s^{2}-2s}{2}.$

On the other hand, we see that any constant multiplied to $f$ will integrate
to zero, and thus, looking at the definition of the Green’s function it seems
natural to add a quadratic polynomial to the Green’s function (which will
become a constant after applying the Laplace operator). In particular, we
choose the polynomial $p_{s}(t)=a_{2}(s)t^{2}+a_{1}(s)t+a_{0}(s)$ that obey
the boundary conditions and integrate to $\frac{s^{2}-2s}{2}$. We easily see
that $a_{0}(s)=0$, and from the the boundary condition at 1,

$-2a_{2}(s)=a_{1}(s).$

We now simply need to determine $a_{2}(s)$ such that

$\frac{s^{2}-2s}{2}=\int_{0}^{1}p_{s}(t)\,\mathrm{d}t=\int_{0}^{1}a_{2}(s)t^{2}-2a_{2}(s)t\,\mathrm{d}t=-\frac{2}{3}a_{2}(s).$

We see that, $a_{2}(s)=-\frac{3}{4}(s^{2}-2s)$, and thus the Green’s function
becomes

$\mathcal{G}(s,t)=(s\wedge t)-\frac{3}{4}(s^{2}-2s)(t^{2}-2t).$

######  Example 2.16 (Dirichlet boundary).

In this example we will consider the Laplace operator $-\partial_{t}^{2}$ on
function spaces $\mathscr{H}$ consisting of functions
$f\in\mathcal{C}^{2}([0,1])$ satisfying homogeneous Dirichlet boundary
conditions $f(0)=f(1)=0$ and possibly other conditions. In the simplest case
we saw that the Green’s function was

$\displaystyle\mathcal{G}(s,t)=s\wedge t-st.$ (19)

Consider the additional constraint

$\int_{0}^{1}f(t)\,\mathrm{d}t=0.$

As in the previous example, it seems natural to substract a polynomial that is
second order in $s$ and $t$ from the Green’s function (19) to find a new
Green’s function. We see that

$\int_{0}^{1}(s\wedge t-st)\,\mathrm{d}t=\frac{1}{2}(1-s)s.$

From the boundary conditions, the polynomial must be of the form

$p_{s}(t)=a(s)t^{2}-a(s)t.$

We have that

$\int_{0}^{1}a(s)(t-1)t\,\mathrm{d}t=-\frac{1}{6}a(s)$

so $a(s)=3(1-s)s$. The Green’s function becomes

$\mathcal{G}(s,t)=s\wedge t-st-3(1-s)s(1-t)t.$

$\circ$

######  Example 2.17 (Second order polynomials).

Consider the constraint that $f^{\prime\prime\prime}(t)=0$. The subspace of
such functions are the second order polynomials. Now let $\mathscr{H}$ be the
space of second order polynomials subject to the mixed Boundary conditions
$f(0)=f^{\prime}(1)=0$. This means that any function $f\in\mathscr{H}$ has the
form

$f(t)=a(t-2)t,\qquad a\in\mathbb{R}.$

The Green’s function should obey that

$\int_{0}^{1}-\partial_{t}^{2}\mathcal{G}(s,t)f(t)\,\mathrm{d}t=f(s),$

but, using integration by parts twice, we get that the above equation is
equivalent to

$-2a\int_{0}^{1}\mathcal{G}(s,t)\,\mathrm{d}t=a(s-2)s.$

Thus,

$\int_{0}^{1}\mathcal{G}(s,t)\,\mathrm{d}t=-\frac{1}{2}(s-2)s,$

meaning that the Green’s function is

$\mathcal{G}(s,t)=\frac{3}{4}(s-2)s(t-2)t.$

$\circ$

######  Example 2.18 (Second order polynomial bridges).

Consider the same setup as in Example 2.17, but with the subspace constraint
that $\mathscr{H}$ is the space of second order polynomials subject to
homogeneous Dirichlet boundary conditions. This means that any function
$f\in\mathscr{H}$ has the form

$f(t)=a(t-1)t,\qquad a\in\mathbb{R}.$

The Green’s function should obey that

$\int_{0}^{1}-\partial_{t}^{2}\mathcal{G}(s,t)(a(t-1)t)\,\mathrm{d}t=a(s-1)s,$

but we have that

$\int_{0}^{1}\partial_{t}^{2}\mathcal{G}(s,t)t\,\mathrm{d}t=0$

and

$\int_{0}^{1}\partial_{t}^{2}\mathcal{G}(s,t)t^{2}\,\mathrm{d}t=2\int_{0}^{1}\mathcal{G}(s,t)\,\mathrm{d}t.$

Thus it must hold that

$\int_{0}^{1}\mathcal{G}(s,t)\,\mathrm{d}t=-\frac{1}{2}(s-1)s.$

Thus the Green’s function has the form

$\mathcal{G}(s,t)=3s(1-s)t(1-t).$

$\circ$

##  3 Smoothing splines

Let pairs of observation times and observations
$(t_{1},y_{1}),\dots,(t_{m},y_{m})$ be given. Assume that the observations are
generated according to the statistical model

$\displaystyle y_{i}=\theta(t_{i})+\varepsilon_{i}$ (20)

where the $\varepsilon_{i}$s are independent
$\mathcal{N}(0,\sigma^{2})$-variables, and $\theta\in\mathscr{H}$ where
$\mathscr{H}$ is some space of functions from the unit interval into the
reals. The negative log likelihood is

$\displaystyle\ell(\theta,\sigma^{2})=\frac{m}{2}\log\sigma^{2}+\frac{1}{2\sigma^{2}}\sum_{i=1}^{m}(y_{i}-\theta(t_{i}))^{2},$
(21)

but if $\mathscr{H}$ is a flexible space of functions, such as
$\mathcal{C}^{2}([0,1])$ there will be infinitely many functions that matches
the observations exactly, so the maximum likelihood estimator is not unique.
Furthermore, since we have noise on our observations, we do not want our
estimator of $\theta$ to be an interpolating function, but rather a function
that matches the trend of the observations while filtering out some of the
noise. One way of achieving this is to add a term to the likelihood (21) that
penalizes roughness. Here we will consider the integrated square derivative
magnitude as our choice of roughness measure

$P(\theta)=\int_{0}^{1}\|\theta^{\prime}(t)\|^{2}\,\mathrm{d}t.$

We can define the penalized log likelihood functional in $\theta$ as

$\displaystyle\ell_{\lambda}(\theta)=\sum_{i=1}^{m}(y_{i}-\theta(t_{i}))^{2}+\lambda
P(\theta)$ (22)

where $\lambda>0$ is a parameter that controls the weighting between data
fidelity and regularity. We will call the minimizer of this functional

$\hat{\theta}=\operatorname*{arg\,min}_{\theta\in\mathscr{H}}\ell_{\lambda}(\theta)$

a smoothing spline (although it will typically not be smooth). In the
following we will go through the necessary theory of identifying the smoothing
spline.

###  3.1 Calculus of variations

Minimization of functionals is the topic of calculus of variations. The
functional derivative $\delta F/\delta f$ of a functional $F\,:\,f\mapsto
F(f)$ in the present setting is defined by the equation

$\displaystyle\int_{0}^{1}\frac{\delta F}{\delta
f}(t)h(t)\,\mathrm{d}t=\frac{\mathrm{d}}{\mathrm{d}\epsilon}F(f+\epsilon
h)\Big{|}_{\epsilon=0}$ (23)

where $h$ is an arbitrary function living in the same function space as $f$.
As with conventional derivatives, a minimizer is found by solving for the
functional derivative equal to the zero function.

######  Example 3.1.

Let us find the functional derivative of the roughness measure $P$. Taking the
right-hand side of (23), we see that

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\epsilon}\int_{0}^{1}\|(\theta+\epsilon
h)^{\prime}(t)\|^{2}\,\mathrm{d}t\Big{|}_{\epsilon=0}$
$\displaystyle=\frac{\mathrm{d}}{\mathrm{d}\epsilon}\int_{0}^{1}(\theta^{\prime}(t)+\epsilon
h^{\prime}(t))^{2}\,\mathrm{d}t\Big{|}_{\epsilon=0}$
$\displaystyle=2\int_{0}^{1}\theta^{\prime}(t)h^{\prime}(t)\,\mathrm{d}t.$

Using integration by parts, this integral can be rewritten

$\displaystyle
2\int_{0}^{1}\theta^{\prime}(t)h^{\prime}(t)\,\mathrm{d}t=2(\theta^{\prime}(1)h(1)-\theta^{\prime}(0)h(0))-2\int_{0}^{1}\theta^{\prime\prime}(t)h(t)\,\mathrm{d}t.$
(24)

Thus, if we for example consider a space of functions satisfying homogeneous
Dirichlet boundary conditions $f(0)=f(1)=0$, the functional derivative is
$\delta P/\delta\theta=2\partial_{t}^{2}\theta$. To minimize the roughness
measure we thus need a function $\theta$ such that
$\partial_{t}^{2}\theta(t)=0$, which in the present case amounts to
$\theta(t)=0$. $\circ$

In the previous example we saw that minimization of functionals may be related
to differential equations. There are a number of interesting links between
splines, Gaussian processes and stochastic differential equations. In the
following we will find the penalized maximum likelihood estimator for
$\theta$. First, we need to identify the functional derivative of the data
term in the penalized log likelihood (22). Computing the right-hand side of
(23) we get

$-2\sum_{i=1}^{m}(y_{i}-\theta(t_{i}))h(t_{i}).$

As a consequence, the functional derivative must be of the form

$-2\sum_{i=1}^{m}(y_{i}-\theta(t))\delta_{t_{i}}(t).$

If we choose boundary conditions such that the boundary terms in (24)
disappear, we can write the functional derivative of $\ell_{\lambda}$ as

$-2\sum_{i=1}^{m}(y_{i}-\theta(t))\delta_{t_{i}}(t)-2\lambda\theta^{\prime\prime}(t).$

Suppose the Laplace operator $-\partial_{t}^{2}$ is invertible on
$\mathscr{H}$ with Green’s function $\mathcal{G}$. Then we can write
$\delta_{t_{i}}=-\partial_{t}^{2}\mathcal{G}(t_{i},t)$. Solving the for the
functional derivative equal to the zero function then leads to

$\frac{1}{\lambda}\sum_{i=1}^{m}(y_{i}-\theta(t))(-\partial_{t}^{2}\mathcal{G}(t_{i},t))=-\partial_{t}^{2}\theta(t).$

Consider a solution on the form

$\theta(t)=\sum_{i=1}^{m}c_{i}\mathcal{G}(t_{i},t),$

insertion in the differential equation gives

$\frac{1}{\lambda}\sum_{i=1}^{m}\bigg{(}y_{i}-\sum_{j=1}^{m}c_{j}\mathcal{G}(t_{j},t)\bigg{)}\delta_{t_{i}}(t)=\sum_{i=1}^{m}c_{i}\delta_{t_{i}}(t),$

but since $\delta_{t_{i}}$ is the point evaluation in $t_{i}$ we can
substitute $\mathcal{G}(t_{j},t)$ with $\mathcal{G}(t_{j},t_{i})$. The
coefficients $c_{i}$ can thus be found by solving a linear system of equations

$\displaystyle\sum_{j=1}^{m}c_{j}(\mathcal{G}(t_{j},t_{i})+\lambda\boldsymbol{1}_{\\{i\\}}(j))=y_{i},\qquad
i=1,\dots,m.$ (25)

Define the matrix of the Green’s function evaluated at all combinations of
observation points

$G=\begin{pmatrix}\mathcal{G}(t_{1},t_{1})&\cdots&\mathcal{G}(t_{1},t_{m})\\\
\vdots&\ddots&\vdots\\\
\mathcal{G}(t_{m},t_{1})&\cdots&\mathcal{G}(t_{m},t_{m})\end{pmatrix}.$

The linear system (25) can be written on the form

$(G-\lambda\mathbb{I}_{m})\boldsymbol{c}=\boldsymbol{y}$

where $\mathbb{I}_{m}$ is the $m\times m$ identity matrix and bold indicates
vectorization. Thus, the coefficients are given by

$\boldsymbol{c}=(G+\lambda\mathbb{I}_{m})^{-1}\boldsymbol{y}$

Finally, define the vector valued function

$G(t)=\begin{pmatrix}\mathcal{G}(t_{1},t)\\\ \vdots\\\
\mathcal{G}(t_{m},t)\end{pmatrix}.$

Putting everything together, we have shown that

$\displaystyle\hat{\theta}(t)=G(t)(G+\lambda\mathbb{I}_{m})^{-1}\boldsymbol{y}.$
(26)

In particular, if we want to evaluate $\hat{\theta}$ in all observation
points, we can compute

$\begin{pmatrix}\hat{\theta}(t_{1})\\\ \vdots\\\ \hat{\theta}(t_{m})\\\
\end{pmatrix}=G(G+\lambda\mathbb{I}_{m})^{-1}\boldsymbol{y}.$

##  4 Gaussian processes and differential equations

A Gaussian process $x$ is a continuous-time stochastic process for which all
finite-dimensional distributions follow multivariate normal distributions. In
this note, we will largely ignore the details of the underlying probability
space and think of $x$ as a random variable on $L^{2}([0,1])$ such that the
individual sample paths $x\,:\,[0,1]\rightarrow\mathbb{R}$ are considered
random functions.

###  4.1 Brownian motion

The most important of all stochastic processes is the Brownian motion. We will
leave its construction to other courses, but give one definition of it. The
Brownian motion $x$ is a stochastic processes for which $x(0)=0$ almost
surely, and the increments

$x(t_{1})-x(0),\quad x(t_{2})-x(t_{1}),\quad\dots\quad x(t_{n})-x(t_{n-1})$

are independent, zero-mean normally distributed with variances $t_{i+1}-t_{i}$
for all choices of $n$ and $0<t_{1}<\dots<t_{n}$.

From this definition, we can derive some important properties of Brownian
motion. Clearly $x$ is a zero-mean process, that is

$\mathrm{E}[x(t)]=0\,\qquad\textrm{for all }t\in[0,1].$

We can also find the covariance function. Let $s\leq t$

$\mathrm{Cov}[x(s),x(t)]=\mathrm{E}[x(s)x(t)]-\mathrm{E}[x(s)]\mathrm{E}[x(t)]=\mathrm{E}[x(s)x(t)]$

Now, write $x(t)=x(t)-x(s)+x(s)$. Then

$\mathrm{Cov}[x(s),x(t)]=\mathrm{E}[x(s)(x(t)-x(s))]+\mathrm{E}[x(s)^{2}]=s.$

We thus conclude that for general $s,t\in[0,1]$

$\mathrm{Cov}[x(s),x(t)]=s\wedge t,$

and note that we found this expression as the Green’s function for the Laplace
operator in equation (9).

######  Example 4.1 (The Brownian bridge).

Let $x$ be a Brownian motion. Let $0\leq s\leq t\leq 1$. Since $x$ is a
Gaussian process, the distribution of $(x(s),x(t),x(1))$ is a multivariate
normal distribution with mean $\mu=(0,0,0)$ and covariance matrix $\Sigma$
given by evaluating the covariance function at all combinations og the points
$s$, $t$ and $1$. That is,

$\Sigma=\begin{pmatrix}s&s&s\\\ s&t&t\\\ s&t&1\end{pmatrix}.$

From the theory of the multivariate normal distribution, we know that if
$(\boldsymbol{x}_{1},\boldsymbol{x}_{2})$ is joint normal with mean
$(\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2})$ and covariance matrix

$\begin{pmatrix}\Sigma_{11}&\Sigma_{12}\\\
\Sigma_{21}&\Sigma_{22}\end{pmatrix},$

then the conditional distribution of
$\boldsymbol{x}_{1}|\boldsymbol{x}_{2}=\boldsymbol{a}$ is again multivariate
normal with mean

$\boldsymbol{\mu}_{1}+\Sigma_{12}\Sigma_{22}^{-1}(\boldsymbol{a}-\boldsymbol{\mu}_{2})$

and covariance matrix

$\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}.$

Now consider the condition $x(1)=0$. We can use the above formulas with
$\boldsymbol{x}_{1}=(x(s),x(t))$ and $\boldsymbol{x}_{2}=x(1)$ to get that the
conditional distribution of $(x(s),x(t))|x(1)=0$ is a zero-mean normal
distribution with covariance matrix

$\Sigma=\begin{pmatrix}s-s^{2}&s-st\\\ s-st&t-t^{2}\end{pmatrix}.$

This process is called the Brownian bridge and we see that its covariance
function must be $s\wedge t-st$. Again we note the connection to the Green’s
function for the Laplace operator found in Example 2.5. $\circ$

###  4.2 Gaussian priors for curves

In Bayesian statistics, one can choose prior probability distributions for
parameters in the model. The Bayesian equivalent to maximum likelihood
estimation is maximum a posteriori estimation, that is, we maximize the
posterior probability given the data and the prior. Assume that we are
interested in the parameter $\theta$ and have data $\boldsymbol{y}$. The
posterior probability (distribution) is given by Bayes’ theorem

$p(\theta\,|\,\boldsymbol{y})=\frac{p(\boldsymbol{y}\,|\,\theta)p(\theta)}{p(\boldsymbol{y})}$

where $p(\boldsymbol{y}\,|\,\theta)$ is the likelihood, $p(\theta)$ is the
prior and $p(\boldsymbol{y})$ is known as the evidence. The evidence may be
very hard to compute because we need to integrate out $\theta$ in the joint
distribution, but typically we can avoid computing this term altogether, since
it is just a normalization constant that is independent of $\theta$. The
maximum a posteriori estimate $\hat{\theta}_{\mathrm{MAP}}$ can be found by
minimizing

$-\log p(\theta\,|\,\boldsymbol{y})\propto-\log
p(\boldsymbol{y}\,|\,\theta)-\log p(\theta).$

Comparing this to the penalized log likelihood (22), we see that both methods
minimize the sum of the negative log likelihood and an extra term which in
both cases acts as a regularization term.

Consider again the statistical model (20) where we observe a function
$\theta\in\mathscr{H}$ with independent, identically distributed Gaussian
noise. We may choose a zero-mean Gaussian process with covariance function
$\sigma^{2}\tau^{2}\mathcal{G}(\,\cdot\,,\,\cdot\,)$ as our prior for
$\theta$. This means that any finite evaluation of $\theta$ in the vector of
points $\boldsymbol{s}$,
$\boldsymbol{\theta}=(\theta(s_{1}),\dots,\theta(s_{\ell}))$ follows a
multivariate normal distribution of the form

$p(\boldsymbol{\theta})=\frac{1}{\sqrt{2^{\ell}\pi^{\ell}\det
G(\boldsymbol{s},\boldsymbol{s})}}\exp\bigg{(}-\frac{1}{2\sigma^{2}\tau^{2}}\boldsymbol{\theta}^{\top}G(\boldsymbol{s},\boldsymbol{s})^{-1}\boldsymbol{\theta}\bigg{)}$

where

$G(\boldsymbol{s},\boldsymbol{s})=\begin{pmatrix}\mathcal{G}(s_{1},s_{1})&\dots&\mathcal{G}(s_{1},s_{\ell})\\\
\vdots&\ddots&\vdots\\\
\mathcal{G}(s_{\ell},s_{1})&\dots&\mathcal{G}(s_{\ell},s_{\ell})\end{pmatrix}.$

The joint distribution of $\boldsymbol{y}$ and $\boldsymbol{\theta}$ is zero-
mean multivariate normal with covariance matrix

$\begin{pmatrix}\sigma^{2}\tau^{2}G(\boldsymbol{t},\boldsymbol{t})+\sigma^{2}\mathbb{I}_{m\times
m}&\sigma^{2}\tau^{2}G(\boldsymbol{t},\boldsymbol{s})\\\
\sigma^{2}\tau^{2}G(\boldsymbol{s},\boldsymbol{t})&\sigma^{2}\tau^{2}G(\boldsymbol{s},\boldsymbol{s})\end{pmatrix}$

where $\boldsymbol{t}=(t_{1},\dots,t_{m})$ is the vector of observation
points. The maximum posterior estimate of $\theta$ in the points
$\boldsymbol{s}$ is the conditional expectation of $\boldsymbol{\theta}$ given
$\boldsymbol{y}$, which is given by

$\displaystyle\hat{\boldsymbol{\theta}}_{\mathrm{MAP}}=G(\boldsymbol{s},\boldsymbol{t})(G(\boldsymbol{t},\boldsymbol{t})+\tau^{-2}\mathbb{I}_{m\times
m})^{-1}\boldsymbol{y}.$ (27)

Comparing with the previously found smoothing spline (26), we see that the
maximum a posteriori estimate will be identical if we choose $\mathcal{G}$ as
the Green’s function for the Laplace operator on $\mathscr{H}$ and set
$\lambda=\tau^{-2}$.

###  4.3 Conditional Gaussian processes

In this section we will consider a zero-mean Gaussian process $x$ with
covariance function $\mathcal{G}$.

######  Example 4.2 (Neumann boundary conditions).

Consider the observations of $x$ at the points $s$, $t$, $1-\epsilon$ and $1$,
ordered increasingly. We see that $(x(s),x(t),(x(1)-x(1-\epsilon))/\epsilon)$
follow a zero-mean normal distribution with covariance

$\begin{pmatrix}\scriptstyle\mathcal{G}(s,s)&\scriptstyle\mathcal{G}(s,t)&\frac{\mathcal{G}(s,1)-\mathcal{G}(s,1-\epsilon)}{\epsilon}\\\\[1.99997pt]
\scriptstyle\mathcal{G}(t,s)&\scriptstyle\mathcal{G}(t,t)&\frac{\mathcal{G}(t,1)-\mathcal{G}(t,1-\epsilon)}{\epsilon}\\\\[1.99997pt]
\frac{\mathcal{G}(s,1)-\mathcal{G}(s,1-\epsilon)}{\epsilon}&\frac{\mathcal{G}(t,1)-\mathcal{G}(t,1-\epsilon)}{\epsilon}&\frac{\mathcal{G}(1,1)+\mathcal{G}(1-\epsilon,1-\epsilon)-2\mathcal{G}(1,1-\epsilon)}{\epsilon^{2}}\\\
\end{pmatrix}.$

Thus, the conditional distribution of $(x_{s},x_{t})$ given that
$(x(1)-x(1-\epsilon))/\epsilon=0$ is zero-mean normal with covariance

$\begin{pmatrix}\scriptstyle\mathcal{G}(s,s)&\scriptstyle\mathcal{G}(s,t)\\\\[1.99997pt]
\scriptstyle\mathcal{G}(t,s)&\scriptstyle\mathcal{G}(t,t)\end{pmatrix}-\frac{\scriptstyle\epsilon^{2}}{\scriptstyle\mathcal{G}(1,1)+\mathcal{G}(1-\epsilon,1-\epsilon)-2\mathcal{G}(1,1-\epsilon)}\begin{pmatrix}\scriptstyle\mathcal{G}_{\epsilon}^{\prime}(s,1)^{2}&\scriptstyle\mathcal{G}_{\epsilon}^{\prime}(s,1)\mathcal{G}_{\epsilon}^{\prime}(t,1)\\\\[1.99997pt]
\scriptstyle\mathcal{G}_{\epsilon}^{\prime}(s,1)\mathcal{G}_{\epsilon}^{\prime}(t,1)&\scriptstyle\mathcal{G}_{\epsilon}^{\prime}(t,1)^{2}\end{pmatrix}$

where
$\mathcal{G}_{\epsilon}^{\prime}(s,1)=\mathcal{G}_{\epsilon}^{\prime}(1,s)=\frac{\mathcal{G}(s,1)-\mathcal{G}(s,1-\epsilon)}{\epsilon}$.
$\circ$

######  Example 4.3.

Let $x$ be a Brownian motion. The covariance function is
$\mathcal{G}(s,t)=s\wedge t$. The second term in the formula for the
conditional variance from Example 4.2 is zero, since
$\mathcal{G}_{\epsilon}^{\prime}(s,1)=0$ for all $s\leq 1-\epsilon$. Thus the
covariance is unchanged, which is also what we would expect. This result
easily generalizes (see Exercise 12) to any future increment of any size,
which means that the Brownian motion is independent of future increments.
$\circ$

######  Example 4.4.

Let $x$ be a Brownian bridge. The covariance function is
$\mathcal{G}(s,t)=s\wedge t-st$. The formula for the conditional variance from
Example 4.2 are now

$\begin{pmatrix}s-s^{2}&t-st\\\
t-st&t-t^{2}\end{pmatrix}-\frac{\epsilon}{1-\epsilon}\begin{pmatrix}s^{2}&st\\\
st&t^{2}\end{pmatrix}.$

As $\epsilon\rightarrow 0$, the covariance just remains that of the Brownian
bridge. $\circ$

##  5 Exercises

  * 1.

Consider the Poisson equation (3) on the following three function spaces

$\displaystyle\mathscr{H}_{1}$
$\displaystyle=\\{f\in\mathcal{C}^{2}\,:\,f(0)=f(1)=0\\},$
$\displaystyle\mathscr{H}_{2}$
$\displaystyle=\\{f\in\mathcal{C}^{2}\,:\,f(0)=f^{\prime}(1)=0\\},$
$\displaystyle\mathscr{H}_{3}$
$\displaystyle=\\{f\in\mathcal{C}^{2}\,:\,f^{\prime}(0)=f^{\prime}(1)=0,{\textstyle\int_{0}^{1}f(t)\,\mathrm{d}t=0}\\}.$

Decide whether $-\partial^{2}_{t}f(t)$ can equal the constant function
$1(t)=1$ for $f$ belonging to each of these sets. Explain the relation to
examples 2.1 and 2.7.

  * 2.

Show that

$\mathscr{H}=\\{f\in\mathcal{C}^{2}([0,1])\,:\,{\textstyle\int_{0}^{1}f(t)\,\mathrm{d}t=0}\\}$

is a subspace of $\mathcal{C}^{2}([0,1])$.

  * 3.

Show that

$\mathscr{H}=\\{f\in\mathcal{C}^{2}([0,1])\,:\,f(0)=f(1)\\}$

is a subspace of $\mathcal{C}^{2}([0,1])$.

  * 4.

Show that

$\mathscr{H}_{1}=\\{f\in\mathcal{C}^{2}([0,1])\,:\,f(t)=f(1-t)\text{ for all
}t\in[0,1/2]\\}$

and

$\mathscr{H}_{2}=\\{f\in\mathcal{C}^{2}([0,1])\,:\,f(t)=-f(1-t)\text{ for all
}t\in[0,1/2]\\}$

are subspaces of $\mathcal{C}^{2}([0,1])$.

  * 5.

Let $\mathscr{H}$ be the linear subspace of functions in
$\mathcal{C}^{2}([0,1])$ for which $f(0)=f(1)=0$ and $f^{(4)}(t)=0$. Identify
the Green’s function for the Laplace operator on $\mathscr{H}$.

  * 6.

Let $\mathscr{H}$ be the space of functions in $\mathcal{C}^{2}([0,1])$ that
obey homogeneous Dirichlet boundary conditions. Express this constraint in
terms of Fourier coefficients.

  * 7.

Let $x$ be a Brownian motion. Define a new process $\tilde{x}(t)=x(t)-tx(1)$.
compute its mean and covariance function. What process is $\tilde{x}$?

  * 8.

Let $\tilde{x}$ be a backward Brownian motion, that is, $\tilde{x}(t)=x(1-t)$
for $t\in[0,1]$ where $x$ is a standard Brownian motion. Determine the
covariance function for $\tilde{x}$.

  * 9.

Let $x$ be the sum of independent forward and backward Brownian motions (see
previous exercise). Determine the covariance function for $x$. Determine
boundary conditions that this covariance function obeys.

  * 10.

Let $\tilde{x}$ be the sum of a forward and backward Brownian motion (see
previous exercises), but this time assume that the backward Brownian motion
arise from the forward Brownian motion, that is $\tilde{x}(t)=x(t)+x(1-t)$
where $x$ is the standard Brownian motion. Determine the covariance function
for $x$. Determine boundary conditions that this covariance function obeys.

  * 11.

Argue why it must hold that the conditional expectation maximizes
$\mathrm{E}[\theta\,|\,\boldsymbol{y}]$ the posterior
$p(\theta\,|\,\boldsymbol{y})$ when we assume that $\theta$ is a Gaussian
process.

  * 12.

Show the general version of the result about conditioning on future increments
for the Brownian motion mentioned in Example 4.3. You can either use the
conditional distribution, or use the definiton of Brownian motion.

  * 13.

Give an interpretation of the result in Example 2.17 in terms of the Brownian
motion.

## References

  * [1] John C Holladay  “A smoothest curve approximation”  In _Mathematical tables and other aids to computation_ 11.60 JSTOR, 1957, pp. 233–243
  * [2] George Kimeldorf and Grace Wahba  “Some results on Tchebycheffian spline functions”  In _Journal of Mathematical Analysis and Applications_ 33.1 Elsevier, 1971, pp. 82–95
  * [3] Bo Markussen  “Functional data analysis in an operator-based mixed-model framework”  In _Bernoulli_ 19, 2013, pp. 1–17
  * [4] David Mumford and Agnès Desolneux  “Pattern theory: the stochastic analysis of real-world signals”  CRC Press, 2010
  * [5] Lars Lau Raket and Bo Markussen  “Approximate inference for spatial functional data on massively parallel processors”  In _Computational Statistics & Data Analysis_ 72, 2014, pp. 227 –240  DOI: 10.1016/j.csda.2013.10.016
  * [6] Helene Charlotte Rytgaard  “Statistical models for robust spline smoothing”, 2016
  * [7] Isaac Jacob Schoenberg  “Spline interpolation and best quadrature formulae”  In _Bulletin of the American Mathematical Society_ 70.1, 1964, pp. 143–148
  * [8] Grace Wahba  “Spline Models for Observational Data”  Society for IndustrialApplied Mathematics, 1990


#  2D discrete Yang-Mills equations on the torus

Volodymyr Sushch
Koszalin University of Technology
Sniadeckich 2, 75-453 Koszalin, Poland
<EMAIL_ADDRESS>
###### Abstract

In this paper, we introduce a discretization scheme for the Yang-Mills
equations in the two-dimensional case using a framework based on discrete
exterior calculus. Within this framework, we define discrete versions of the
exterior covariant derivative operator and its adjoint, which capture
essential geometric features similar to their continuous counterparts. Our
focus is on discrete models defined on a combinatorial torus, where the
discrete Yang-Mills equations are presented in the form of both a system of
difference equations and a matrix form.

Key words and phrases: Yang-Mills equations, discrete exterior calculus,
discrete operators, combinatorial torus, difference equations

Math. Subj. Clas.: 39A12, 39A70, 81T13

##  1 Introduction

Discrete models of mathematical physics problems that preserve a geometric
structure are crucial for the successful numerical simulation of differential
equations. In [22], a discrete version of the de Rham-Hodge theory in the
specific case of a two-dimensional torus has been discussed in detail. In this
paper, motivated by this and considering the Yang-Mills theory as nonlinear
generalization of Hodge theory, we construct a discrete counterpart of Yang-
Mills equations in the two-dimensional case. Our discretization scheme is
based on the geometric approach proposed by Dezin in [11]. This approach is
characterized by the construction of a discrete exterior calculus on a complex
of real-valued cochains. In author’s previous papers [17, 18, 19, 20, 21],
this technique has been developed on matrix-valued complexes of cochains, and
discrete analogues of the Yang-Mills equations in $\mathbb{R}^{n}$ and in
Minkowski space have been studied. In [17, 19], we focused on the construction
of gauge invariant discrete models. Depending on the definition of a discrete
version of the Hodge star, we discuss two types of discrete Yang-Mills
equations. In [20, 21], we study discrete analogues of instanton and anti-
instanton solutions.

There are various approaches to discretizing Yang-Mills theories, with
numerous papers dedicated to this subject. For instance, see [5, 6, 14, 15,
16] and the references therein. Many of these approaches rely on lattice
discretization schemes. However, in the lattice formulation, maintaining the
geometrical properties of the original gauge theory poses challenges. Two
dimensional Yang-Mills theories have been investigated from several points of
view in [2, 3, 9, 13, 23].

Recently, from a numerical simulation perspective, interest in the
discretization of the Yang-Mills equations has been fueled by the development
of methods that preserve the geometric and topological structures of their
continuum counterparts. Constraint-preserving numerical methods for the Yang-
Mills equations, based on the finite element exterior calculus [1], have been
developed and investigated in [7, 8, 4]. In [10], within framework of a
discrete de Rham method, another structure-preserving numerical scheme to
approximate the Yang-Mills equations has been constructed.

Let’s provide brief reminders of the main definitions of smooth Yang-Mills
theory. Although our focus in this paper is on the two dimensional torus, we
will describe the Yang–Mills action in the more general context of compact
Riemannian manifolds of arbitrary dimension. Let $M$ be a compact connected
Riemannian manifold. For simplicity, we assume that a principal $G$-bundle $P$
over $M$ is trivial and that the structure group $G$ is the Lie group $SU(2)$.
A connection on $P$ is defined as the $su(2)$-valued differential 1-form $A$,
where $su(2)$ is the Lie algebra of the Lie group $SU(2)$. The curvature
2-form $F$ of the connection $A$ is defined as follows

$F=dA+A\wedge A,$ (1.1)

where $d$ is the exterior derivative and $\wedge$ is the exterior product. The
exterior covariant derivative $d_{A}$ of the connection $A$ acts on an
$su(2)$-valued $r$-form $\Omega$ as

$d_{A}\Omega=d\Omega+A\wedge\Omega+(-1)^{r+1}\Omega\wedge A.$ (1.2)

Denote by $\delta_{A}$ the formal adjoint operator of $d_{A}$ with respect to
the $L^{2}$-norm on the sections of $P$. It is given by $\delta_{A}=\pm\ast
d_{A}\ast$, where $\ast$ is the Hodge star. See [2] for details. The Yang-
Mills equations in terms of differential forms can be represented as

$\delta_{A}F=0,$ (1.3)

or, equivalently, as

$d_{A}\ast F=0.$ (1.4)

Equations (1.3) or (1.4) form a system of second order nonlinear differential
equations. A solution $A$ of (1.3) or (1.4) is referred to as the Yang-Mills
connection. Furthermore, by virtue of the Bianchi identity $d_{A}F=0$, the
Yang-Mills connection satisfies the following Laplace-type equation

$(d_{A}\delta_{A}+\delta_{A}d_{A})F=0.$ (1.5)

This explains why the Yang-Mills theory can be interpreted as nonlinear
generalization of Hodge theory. See, for example, [12] or [2] for details.

Our purpose is to propose a geometric structure-preserving discrete model of
the Yang-Mills equations on a two-dimensional torus. The technique developed
in [22] to describe a discrete analogue of the de Rham-Hodge theory also
allows us to construct and investigate discrete analogues of Equations (1.3)
and (1.4) on a combinatorial torus. We show that discrete analogues of the
operators $d_{A}$ and $\delta_{A}$ exhibit properties similar to those in the
continuous case. Furthermore, we present a system of difference equations
explicitly, which are equivalent to discrete Yang-Mills equations.

The paper is organized as follows. In Section 2, we briefly review the concept
of a combinatorial model of $\mathbb{R}^{2}$, revisiting relevant material
from [22]. We introduce a discrete version of matrix-valued differential forms
(a complex of $su(2)$-valued cochains) and define analogues of the main
exterior calculus operations on them. Section 3 establishes a discrete
analogue of the Yang-Mills equations in the two-dimensional case.
Additionally, we discuss the properties of discrete counterparts of the
covariant exterior derivative and its adjoint operator. In Section 4, we apply
the proposed discretization scheme to a combinatorial torus. We present the
discrete Yang-Mills equations in terms of a system of difference equations and
also provide the matrix form of this system.

##  2 Background on a discrete model

The construction of a combinatorial model of the two-dimensional Euclidean
space $\mathbb{R}^{2}$ is presented in detail in [22]. We will only recall the
main definitions and key operations that are most relevant to this paper. A
combinatorial model of $\mathbb{R}^{2}$ consists of a two-dimensional chain
complex $C(2)=C_{0}(2)\oplus C_{1}(2)\oplus C_{2}(2)$ generated by the 0-, 1-,
and 2-dimensional basis elements $\\{x_{k,s}\\}$,
$\\{e_{k,s}^{1},e_{k,s}^{2}\\}$, and $\\{V_{k,s}\\}$, respectively, where
$k,s\in{\mathbb{Z}}$. The boundary operator $\partial:C_{r}(2)\rightarrow
C_{r-1}(2)$ is defined by

$\displaystyle\partial x_{k,s}=0,\qquad\partial e_{k,s}^{1}=x_{\tau
k,s}-x_{k,s}\qquad\partial e_{k,s}^{2}=x_{k,\tau s}-x_{k,s},$
$\displaystyle\partial V_{k,s}=e_{k,s}^{1}+e_{\tau k,s}^{2}-e_{k,\tau
s}^{1}-e_{k,s}^{2},$ (2.1)

where $\tau$ denotes the shift by one to the right, i.e., $\tau k=k+1$. The
definition (2) is extended to arbitrary chains by linearity.

Let us consider an object dual to the chain complex $C(2)$. The dual complex
$K(2)=K^{0}(2)\oplus K^{1}(2)\oplus K^{2}(2)$ is defined as in [22]. It has a
structure similar to $C(2)$ and consists of cochains with coefficients
belonging to the algebra $su(2)$. For the construction of $su(2)$-valued
cochains see also [17]. Denote by $\\{x^{k,s}\\}$, $\\{e^{k,s}_{1},\
e^{k,s}_{2}\\}$, and $\\{V^{k,s}\\}$ the basis elements of $K^{0}(2)$,
$K^{1}(2)$, and $K^{2}(2)$, respectively. Then the cochains $\Phi\in
K^{0}(2)$, $\Omega\in K^{1}(2)$, and $\Psi\in K^{2}(2)$ can be written as

$\Phi=\sum_{k,s}\Phi_{k,s}x^{k,s},\quad\Omega=\sum_{k,s}(\Omega^{1}_{k,s}e_{1}^{k,s}+\Omega^{2}_{k,s}e_{2}^{k,s}),\quad\Psi=\sum_{k,s}\Psi_{k,s}V^{k,s},$
(2.2)

where $\Phi_{k,s},\Omega^{1}_{k,s},\Omega^{2}_{k,s},\Psi_{k,s}\in su(2)$ for
any $k,s\in{\mathbb{Z}}$. As in [22], we will call cochains forms (discrete
forms). Note the difference between $K(2)$ here and the cochain complex in
[22]. The complex $K(2)$ consists of $su(2)$-valued cochains while in [22]
cochains are real-valued.

For discrete forms (2.2), the pairing is defined with the basis elements of
$C(2)$ according to the rule

$\langle x_{k,s},\ \Phi\rangle=\Phi_{k,s},\ \langle e_{k,s}^{1},\
\Omega\rangle=\Omega^{1}_{k,s},\ \langle e_{k,s}^{2},\
\Omega\rangle=\Omega^{2}_{k,s},\ \langle V_{k,s},\Psi\rangle=\Psi_{k,s}.$
(2.3)

As in [22], we consider the coboundary operator $d^{c}:K^{r}(2)\rightarrow
K^{r+1}(2)$. The coboundary operator can be interpreted as a discrete analogue
of the exterior differential. For the forms (2.2), we have

$d^{c}\Phi=\sum_{k,s}(\Phi_{\tau k,s}-\Phi_{k,s})e_{1}^{k,s}+(\Phi_{k,\tau
s}-\Phi_{k,s})e_{2}^{k,s},$ (2.4) $d^{c}\Omega=\sum_{k,s}(\Omega_{\tau
k,s}^{2}-\Omega_{k,s}^{2}-\Omega_{k,\tau s}^{1}+\Omega_{k,s}^{1})V^{k,s},$
(2.5)

and $d^{c}\Psi=0$.

Finally, let’s recall the definitions of the $\cup$ product and the star
operator on $K(2)$, following [22]. For the basis elements of $K(2)$, the
$\cup$ product is defined as follows

$x^{k,s}\cup x^{k,s}=x^{k,s},\qquad x^{k,s}\cup e^{k,s}_{1}=e^{k,s}_{1},\qquad
x^{k,s}\cup e^{k,s}_{2}=e^{k,s}_{2},$ $x^{k,s}\cup V^{k,s}=V^{k,s},\qquad
e^{k,s}_{1}\cup x^{\tau k,s}=e^{k,s}_{1},\qquad e^{k,s}_{2}\cup x^{k,\tau
s}=e^{k,s}_{2},$ $V^{k,s}\cup x^{\tau k,\tau s}=V^{k,s},\qquad e^{k,s}_{1}\cup
e^{\tau k,s}_{2}=V^{k,s},\qquad e^{k,s}_{2}\cup e^{k,\tau s}_{1}=-V^{k,s},$

with the product being zero in all other cases. This operation is extended to
arbitrary forms by linearity, where the form coefficients are multiplied as
matrices. In [11, Ch.3, Proposition 2], it is shown that for real-valued
discrete forms, the discrete counterpart of the Leibniz rule is valid. For any
matrix-valued forms, this rule holds as well, and we have

$d^{c}(\Omega\cup\Phi)=d^{c}\Omega\cup\Phi+(-1)^{r}\Omega\cup d^{c}\Phi,$
(2.6)

where $r$ is the degree of $\Omega$. The star operator
$\ast:K^{r}(2)\rightarrow K^{2-r}(2)$ is defined by the rule

$\ast x^{k,s}=V^{k,s},\quad\ast e^{k,s}_{1}=e^{\tau k,s}_{2},\quad\ast
e^{k,s}_{2}=-e^{k,\tau s}_{1},\quad\ast V^{k,s}=x^{\tau k,\tau s}.$ (2.7)

Again, it’s extended to arbitrary forms by linearity. The operation $\ast$ has
properties similar to the Hodge star operator. Thus, it can be regarded as a
discrete analogue of the Hodge star operator.

Let’s now shift our focus to the inner product on $K(2)$. We’ll start by
recalling that the real Lie algebra $su(2)$ has a basis given by

$E_{1}=\begin{bmatrix}i&0\\\ 0&-i\end{bmatrix},\quad
E_{2}=\begin{bmatrix}0&1\\\ -1&0\end{bmatrix},\quad
E_{3}=\begin{bmatrix}0&i\\\ i&0\end{bmatrix},$

where $i$ is the imaginary unit. Any matrix $A\in su(2)$ can be expressed as

$A=a_{1}E_{1}+a_{2}E_{2}+a_{3}E_{3}=\begin{bmatrix}a_{1}i&a_{2}+a_{3}i\\\
-a_{2}+a_{3}i&-a_{1}i\end{bmatrix},$ (2.8)

where $a_{1},a_{2},a_{3}\in\mathbb{R}$. The algebra $su(2)$ is endowed with an
inner product defined as

$(A,B)=-\operatorname{tr}(AB).$

Here, $\operatorname{tr}$ denotes the trace of the matrix product. It’s easy
to show that this forms a real symmetric bilinear form, and in fact, it’s
positive definite.

Let’s denote by $V$ the two-dimensional finite chain with unit coefficients,
defined by

$V=\sum_{k=1}^{N}\sum_{s=1}^{M}V_{k,s}.$ (2.9)

Now, we define the inner product of discrete forms on $V$ as follows

$(\Phi,\ \Omega)_{V}=-\frac{1}{2}\operatorname{tr}\langle V,\
\Phi\cup\ast\Omega\rangle,$ (2.10)

where $\Phi$ and $\Omega$ are discrete forms of the same degree. For forms of
different degrees, the product (2.10) is set equal to zero. From (2.3) and
(2.7), using the definition of the $\cup$ product, we obtain

$(\Phi,\
\Omega)_{V}=-\frac{1}{2}\operatorname{tr}\sum_{k=1}^{N}\sum_{s=1}^{M}\Phi_{k,s}\Omega_{k,s}=\sum_{k=1}^{N}\sum_{s=1}^{M}\sum_{\alpha=1}^{3}\varphi_{k,s,\alpha}\omega_{k,s,\alpha},$
(2.11)

where $\Phi_{k,s}=[\varphi_{k,s,\alpha}]$ and
$\Omega_{k,s}=[\omega_{k,s,\alpha}]$ are the matrices of the form (2.8). Thus,
the norm is given by

$\|\Phi\|^{2}=(\Phi,\
\Phi)_{V}=\sum_{k=1}^{N}\sum_{s=1}^{M}\sum_{\alpha=1}^{3}(\varphi_{k,s,\alpha})^{2}.$
(2.12)

Now, let’s define the adjoint operator of $d^{c}$ with respect to the inner
product (2.10). As in [22, Proposition 2], we denote the adjoint operator of
$d^{c}$ as $\delta^{c}$ which satisfies the following relation.

######  Proposition 2.1.

Let $\Phi\in K^{r}(2)$ and $\Omega\in K^{r+1}(2)$, $r=0,1$. Then we have

$(d^{c}\Phi,\ \Omega)_{V}=-\frac{1}{2}\operatorname{tr}\langle\partial V,\
\Phi\cup\ast\Omega\rangle+(\Phi,\ \delta^{c}\Omega)_{V},$ (2.13)

where

$\delta^{c}\Omega=(-1)^{r+1}\ast^{-1}d^{c}\ast\Omega$ (2.14)

and $\ast^{-1}$ is the inverse of $\ast$.

###### Proof.

Taking into account (2.6) and (2.10), the proof coincides with that [22,
Proposition 2]. ∎

It is clear that the operator $\delta^{c}:K^{r+1}(2)\rightarrow K^{r}(2)$
defined by (2.14) serves as a discrete counterpart to the codifferential
$\delta$ (the operator adjoint of $d$). Note that by (2.7), we have

$\ast^{-1}x^{k,s}=V^{\sigma k,\sigma s},\quad\ast^{-1}e^{k,s}_{1}=-e^{k,\sigma
s}_{2},\quad\ast^{-1}e^{k,s}_{2}=e^{\sigma
k,s}_{1},\quad\ast^{-1}V^{k,s}=x^{k,s}.$

Here, $\sigma$ denotes the shift by one to the left, i.e., $\sigma k=k-1$.
Combining these relations with (2.4) and (2.5) for the forms (2.2), we obtain
$\delta^{c}\Phi=0$, and

$\delta^{c}\Omega=\sum_{k=1}^{N}\sum_{s=1}^{M}(\Omega_{\sigma
k,s}^{1}-\Omega_{k,s}^{1}+\Omega_{k,\sigma s}^{2}-\Omega_{k,s}^{2})x^{k,s},$
(2.15) $\delta^{c}\Psi=\sum_{k=1}^{N}\sum_{s=1}^{M}(\Psi_{k,s}-\Psi_{k,\sigma
s})e_{1}^{k,s}-(\Psi_{k,s}-\Psi_{\sigma k,s})e_{2}^{k,s}.$ (2.16)

######  Remark 2.2.

The inner product is defined over $V$ in the form (2.9), but the relation
(2.13) includes terms with the components $\Phi_{0,s}$, $\Phi_{\tau N,s}$,
$\Omega_{0,s}$, $\Omega_{\tau N,s}$, $\Phi_{k,0}$, $\Phi_{k,\tau M}$,
$\Omega_{k,0}$ and $\Omega_{k,\tau M}$, which are not specified by the
definition of the form $\Phi$, and $\Omega$ on $V$. Hence, these components
must be defined additionally.

Assume the components of the form (2.2) satisfy the following conditions

$\displaystyle\Phi_{0,s}=\Phi_{N,s},\quad\Phi_{\tau
N,s}=\Phi_{1,s},\quad\Phi_{k,0}=\Phi_{k,M},\quad\Phi_{k,\tau M}=\Phi_{k,1},$
$\displaystyle\Psi_{0,s}=\Psi_{N,s},\quad\Psi_{\tau
N,s}=\Psi_{1,s},\quad\Psi_{k,0}=\Psi_{k,M},\quad\Psi_{k,\tau M}=\Psi_{k,1},$
$\displaystyle\Omega^{1}_{0,s}=\Omega^{1}_{N,s},\quad\Omega^{1}_{\tau
N,s}=\Omega^{1}_{1,s},\quad\Omega^{1}_{k,0}=\Omega^{1}_{k,M},\quad\Omega^{1}_{k,\tau
M}=\Omega^{1}_{k,1},$
$\displaystyle\Omega^{2}_{0,s}=\Omega^{2}_{N,s},\quad\Omega^{2}_{\tau
N,s}=\Omega^{2}_{1,s},\quad\Omega^{2}_{k,0}=\Omega^{2}_{k,M},\quad\Omega^{2}_{k,\tau
M}=\Omega^{2}_{k,1}$ (2.17)

for any $k=1,2,...,N$ and $s=1,2,..,M$. We denote by $H^{r}(V)$ the Hilbert
space generated by the inner product (2.10) of $r$-forms satisfying conditions
(2). For the operators

$d^{c}:H^{r}(V)\rightarrow H^{r+1}(V),\qquad\delta^{c}:H^{r+1}(V)\rightarrow
H^{r}(V),$

the relation (2.13) takes the form

$(d^{c}\Phi,\ \Omega)_{V}=(\Phi,\ \delta^{c}\Omega)_{V}.$ (2.18)

We’ll omit the proof here. See [22, Proposition 3] for details. However, it’s
worth noting that

$\partial V=\sum_{k=1}^{N}e^{1}_{k,1}+\sum_{s=1}^{M}e^{2}_{\tau
N,s}-\sum_{k=1}^{N}e^{1}_{k,\tau M}-\sum_{s=1}^{M}e^{2}_{1,s},$

and for any $\Phi\in H^{r}(V)$, $\Omega\in H^{r+1}(V)$, we have

$\langle\partial V,\ \Phi\cup\ast\Omega\rangle=0.$

The operator

$\Delta^{c}=d^{c}\delta^{c}+\delta^{c}d^{c}:H^{r}(V)\rightarrow H^{r}(V)$

serves as a discrete analogue of the Laplacian. Furthermore, $\Delta^{c}$ is
self-adjoint, i.e.,

$(\Delta^{c}\Phi,\ \Omega)_{V}=(\Phi,\ \Delta^{c}\Omega)_{V},$

where $\Phi$, $\Omega\in H^{r}(V)$. This follows immediately from (2.18).

##  3 Discrete Yang-Mills equations

In this section, we establish a discrete counterpart of the Yang-Mills
equations for the model considered in the previous section. We also discuss
some properties of operators generated by the discrete Yang-Mills equations.
Let the 1-form

$A=\sum_{k=1}^{N}\sum_{s=1}^{M}(A^{1}_{k,s}e_{1}^{k,s}+A^{2}_{k,s}e_{2}^{k,s}),$
(3.1)

where $A^{1}_{k,s},A^{2}_{k,s}\in su(2)$, be a discrete analogue of the
$su(2)$-valued connection 1-form. The discrete analogue of the curvature form
(1.1) is defined by

$F=d^{c}A+A\cup A.$ (3.2)

Using (2.5) and the definition of $\cup$, the components $F_{k,s}\in su(2)$ of
the 2-form

$F=\sum_{k=1}^{N}\sum_{s=1}^{M}F_{k,s}V^{k,s}$

can be represented as

$F_{k,s}=A^{2}_{\tau k,s}-A^{2}_{k,s}-A^{1}_{k,\tau
s}+A^{1}_{k,s}+A^{1}_{k,s}A^{2}_{\tau k,s}-A^{2}_{k,s}A^{1}_{k,\tau s}$ (3.3)

for any $k,s$. Now, let’s consider the operator $d^{c}_{A}:H^{r}(V)\rightarrow
H^{r+1}(V)$ defined by

$d^{c}_{A}\Omega=d^{c}\Omega+A\cup\Omega+(-1)^{r+1}\Omega\cup A,$ (3.4)

where $A\in H^{1}(V)$ and $\Omega\in H^{r}(V)$. Since $A$ represents the
discrete connection 1-form, we refer to $d^{c}_{A}$ as the discrete analogue
of the covariant exterior derivative (1.2). It is obvious that the equation
$d^{c}_{A}F=0$ is an identity (a discrete analogue of the Bianchi identity).
The discrete counterpart of the Yang-Mills equations (1.4) can be given as

$d^{c}_{A}\ast F=0,$ (3.5)

where $F$ is the discrete curvature form (3.2). The form $\ast F$ is a 0-form,
and applying (2.7) yields

$\ast F=\sum_{k=1}^{N}\sum_{s=1}^{M}F_{k,s}x^{\tau k,\tau
s}=\sum_{k=1}^{N}\sum_{s=1}^{M}F_{\sigma k,\sigma s}x^{k,s},$

where $F_{k,s}$ is given by (3.3). Then Equation (3.5) can be represented as

$\langle e^{1}_{k,s},\ d_{A}^{c}\ast F\rangle=F_{k,\sigma s}-F_{\sigma
k,\sigma s}+A^{1}_{k,s}F_{k,\sigma s}-F_{\sigma k,\sigma s}A^{1}_{k,s}=0,$
$\langle e^{2}_{k,s},\ d_{A}^{c}\ast F\rangle=F_{\sigma k,s}-F_{\sigma
k,\sigma s}+A^{2}_{k,s}F_{\sigma k,s}-F_{\sigma k,\sigma s}A^{2}_{k,s}=0$

for any $k=1,2,...,N$ and $s=1,2,...,M$.

The operator $d_{A}^{c}$ behaves like the coboundary operator $d^{c}$ under
the $\cup$ operation. It follows in particular that we have the following:

######  Proposition 3.1.

Let $\Omega\in H^{r}(V)$. Then we have

$d^{c}_{A}(\Omega\cup\Phi)=d^{c}_{A}\Omega\cup\Phi+(-1)^{r}\Omega\cup
d^{c}_{A}\Phi.$ (3.6)

###### Proof.

Let $\Omega\in H^{r}(V)$ and $\Phi\in H^{p}(V)$, where $r=0,p=1$ or $r=1,p=0$.
Then by (2.6) and (3.4), we obtain

$\displaystyle
d_{A}^{c}(\Omega\cup\Phi)=d^{c}(\Omega\cup\Phi)+A\cup\Omega\cup\Phi+(-1)^{r+p+1}\Omega\cup\Phi\cup
A$ $\displaystyle=d^{c}\Omega\cup\Phi+(-1)^{r}\Omega\cup
d^{c}\Phi+A\cup\Omega\cup\Phi+(-1)^{r+p+1}\Omega\cup\Phi\cup A$
$\displaystyle=d^{c}\Omega\cup\Phi+A\cup\Omega\cup\Phi+(-1)^{r+1}\Omega\cup
A\cup\Phi$ $\displaystyle+(-1)^{r}(\Omega\cup d^{c}\Phi+\Omega\cup
A\cup\Phi+(-1)^{p+1}\Omega\cup\Phi\cup A)$
$\displaystyle=d^{c}_{A}\Omega\cup\Phi+(-1)^{r}\Omega\cup d^{c}_{A}\Phi.$

In cases where $\Omega\cup\Phi\in H^{2}(V)$, both sides of the relation (3.6)
are zero. Therefore, Equation (3.6) holds for all cases, as desired. ∎

It should be noted that under conditions (2), the following 2-forms

$\Phi=\sum_{k=1}^{N}\sum_{s=1}^{M}\Phi_{k,s}V^{k,s},\
\tilde{\Phi}=\sum_{k=1}^{N}\sum_{s=1}^{M}\Phi_{\tau k,s}V^{k,s},\
\hat{\Phi}=\sum_{k=1}^{N}\sum_{s=1}^{M}\Phi_{k,\tau s}V^{k,s}$

take the same value over $V$, i.e.,

$\langle V,\ \Phi\rangle=\langle V,\ \tilde{\Phi}\rangle=\langle V,\
\hat{\Phi}\rangle=\sum_{k=1}^{N}\sum_{s=1}^{M}\Phi_{k,s}.$ (3.7)

######  Proposition 3.2.

Let $\Phi\in H^{r}(V)$ and $\Omega\in H^{2-r}(V)$. Then we have

$\operatorname{tr}\langle V,\ \Phi\cup\Omega\rangle=\operatorname{tr}\langle
V,\ \Omega\cup\ast\ast\Phi\rangle.$ (3.8)

###### Proof.

Let’s start with forms $\Phi\in H^{2}(V)$ and $\Omega\in H^{0}(V)$. These
forms are

$\Phi=\sum_{k=1}^{N}\sum_{s=1}^{M}\Phi_{k,s}V^{k,s},\
\Omega=\sum_{k=1}^{N}\sum_{s=1}^{M}\Omega_{k,s}x^{k,s}.$

By the definition of the $\cup$ product, we have

$\Phi\cup\Omega=\sum_{k=1}^{N}\sum_{s=1}^{M}\Phi_{k,s}\Omega_{\tau k,\tau
s}V^{k,s}$

and thus,

$\operatorname{tr}\langle V,\
\Phi\cup\Omega\rangle=\sum_{k=1}^{N}\sum_{s=1}^{M}\operatorname{tr}(\Phi_{k,s}\Omega_{\tau
k,\tau s}).$

Since

$\ast\ast\Phi=\sum_{k=1}^{N}\sum_{s=1}^{M}\Phi_{k,s}V^{\tau k,\tau
s}=\sum_{k=1}^{N}\sum_{s=1}^{M}\Phi_{\sigma k,\sigma s}V^{k,s},$

we calculate

$\Omega\cup\ast\ast\Phi=\sum_{k=1}^{N}\sum_{s=1}^{M}\Omega_{k,s}\Phi_{\sigma
k,\sigma s}V^{k,s}.$

By (2), it is true that for any 2-form we have

$\langle V,\
\Phi\rangle=\sum_{k=1}^{N}\sum_{s=1}^{M}\Phi_{k,s}=\sum_{k=1}^{N}\sum_{s=1}^{M}\Phi_{\sigma
k,\sigma s}=\langle V,\ \ast\ast\Phi\rangle.$

Hence,

$\displaystyle\operatorname{tr}\langle V,\
\Omega\cup\ast\ast\Phi\rangle=\sum_{k=1}^{N}\sum_{s=1}^{M}\operatorname{tr}(\Omega_{k,s}\Phi_{\sigma
k,\sigma s})=\sum_{k=1}^{N}\sum_{s=1}^{M}\operatorname{tr}(\Omega_{\tau k,\tau
s}\Phi_{k,s})$
$\displaystyle=\sum_{k=1}^{N}\sum_{s=1}^{M}\operatorname{tr}(\Phi_{k,s}\Omega_{\tau
k,\tau s})=\operatorname{tr}\langle V,\ \Phi\cup\Omega\rangle.$

Here, we use the fact that $\operatorname{tr}(AB)=\operatorname{tr}(BA)$ for
any two matrices $A$ and $B$ of appropriate sizes.

Similarly, we derive (3.8) for $\Phi\in H^{0}(V)$ and $\Omega\in H^{2}(V)$.

Finally, let’s consider the forms $\Phi\in H^{1}(V)$ and $\Omega\in H^{1}(V)$
given by

$\Phi=\sum_{k=1}^{N}\sum_{s=1}^{M}(\Phi^{1}_{k,s}e_{1}^{k,s}+\Phi^{2}_{k,s}e_{2}^{k,s}),\
\Omega=\sum_{k=1}^{N}\sum_{s=1}^{M}(\Omega^{1}_{k,s}e_{1}^{k,s}+\Omega^{2}_{k,s}e_{2}^{k,s}).$

In this case, the $\cup$ product is

$\Phi\cup\Omega=\sum_{k=1}^{N}\sum_{s=1}^{M}(\Phi^{1}_{k,s}\Omega^{2}_{\tau
k,s}-\Phi^{2}_{k,s}\Omega^{1}_{k,\tau s})V^{k,s}.$

This gives

$\operatorname{tr}\langle
V,\Phi\cup\Omega\rangle=\sum_{k=1}^{N}\sum_{s=1}^{M}\operatorname{tr}(\Phi^{1}_{k,s}\Omega^{2}_{\tau
k,s}-\Phi^{2}_{k,s}\Omega^{1}_{k,\tau s}).$

Using (2.7), we get

$\ast\ast\Phi=-\sum_{k=1}^{N}\sum_{s=1}^{M}(\Phi^{1}_{\sigma k,\sigma
s}e_{1}^{k,s}+\Phi^{2}_{\sigma k,\sigma s}e_{2}^{k,s}).$ (3.9)

Then we have

$\Omega\cup\ast\ast\Phi=\sum_{k=1}^{N}\sum_{s=1}^{M}(\Omega^{2}_{k,s}\Phi^{1}_{\sigma
k,s}-\Omega^{1}_{k,s}\Phi^{2}_{k,\sigma s})V^{k,s}.$

By (3.7), we obtain

$\displaystyle\operatorname{tr}\langle V,\
\Omega\cup\ast\ast\Phi\rangle=\sum_{k=1}^{N}\sum_{s=1}^{M}\operatorname{tr}(\Omega^{2}_{k,s}\Phi^{1}_{\sigma
k,s}-\Omega^{1}_{k,s}\Phi^{2}_{k,\sigma s})$
$\displaystyle=\sum_{k=1}^{N}\sum_{s=1}^{M}\operatorname{tr}(\Omega^{2}_{\tau
k,s}\Phi^{1}_{k,s}-\Omega^{1}_{k,\tau
s}\Phi^{2}_{k,s})=\sum_{k=1}^{N}\sum_{s=1}^{M}\operatorname{tr}(\Phi^{1}_{k,s}\Omega^{2}_{\tau
k,s}-\Phi^{2}_{k,s}\Omega^{1}_{k,\tau s})$
$\displaystyle=\operatorname{tr}\langle V,\ \Phi\cup\Omega\rangle.$

∎

######  Proposition 3.3.

Let $\Phi\in H^{r}(V)$ and $\Omega\in H^{2-r}(V)$. Then

$(\Phi,\ \ast^{-1}\Omega)_{V}=(\ast\Phi,\ \Omega)_{V}.$

###### Proof.

From the definition of the inner product (2.10) and by Proposition 3.2, it
follows that

$\displaystyle(\Phi,\
\ast^{-1}\Omega)_{V}=-\frac{1}{2}\operatorname{tr}\langle V,\
\Phi\cup\ast\ast^{-1}\Omega\rangle=-\frac{1}{2}\operatorname{tr}\langle V,\
\Phi\cup\Omega\rangle$ $\displaystyle=-\frac{1}{2}\operatorname{tr}\langle V,\
\Omega\cup\ast\ast\Phi\rangle=(\Omega,\ \ast\Phi)_{V}=(\ast\Phi,\
\Omega)_{V}.$

∎

######  Corollary 3.4.

$\operatorname{tr}\langle V,\ \Phi\cup\Omega\rangle=\operatorname{tr}\langle
V,\ \ast\Phi\cup\ast\Omega\rangle.$

######  Proposition 3.5.

Let $\Phi\in H^{r}(V)$ and $\Omega\in H^{r+1}(V)$, $r=0,1$. Then we have

$(d_{A}^{c}\Phi,\ \Omega)_{V}=(\Phi,\ \delta_{A}^{c}\Omega)_{V},$ (3.10)

where

$\delta_{A}^{c}\Omega=\delta^{c}\Omega+\ast^{-1}(\ast\Omega\cup\ast\ast
A+(-1)^{r+1}A\cup\ast\Omega).$ (3.11)

###### Proof.

This is a fairly straight forward proof. Using (2.18), (3.4), and (3.8), we
have

$\displaystyle(d_{A}^{c}\Phi,\
\Omega)_{V}=-\frac{1}{2}\operatorname{tr}\langle V,\
d_{A}^{c}\Phi\cup\ast\Omega\rangle$
$\displaystyle=-\frac{1}{2}\operatorname{tr}\langle V,\
(d^{c}\Phi+A\cup\Phi+(-1)^{r+1}\Phi\cup A)\cup\ast\Omega\rangle$
$\displaystyle=(d^{c}\Phi,\ \Omega)_{V}-\frac{1}{2}\operatorname{tr}\langle
V,\
A\cup\Phi\cup\ast\Omega\rangle-(-1)^{r+1}\frac{1}{2}\operatorname{tr}\langle
V,\ \Phi\cup A\cup\ast\Omega\rangle$ $\displaystyle=(\Phi,\
\delta^{c}\Omega)_{V}-\frac{1}{2}\operatorname{tr}\langle V,\
\Phi\cup\ast\Omega\cup\ast\ast
A\rangle-(-1)^{r+1}\frac{1}{2}\operatorname{tr}\langle V,\ \Phi\cup
A\cup\ast\Omega\rangle$ $\displaystyle=(\Phi,\
\delta^{c}\Omega)_{V}-\frac{1}{2}\operatorname{tr}\langle V,\
\Phi\cup\ast\ast^{-1}(\ast\Omega\cup\ast\ast
A+(-1)^{r+1}A\cup\ast\Omega)\rangle$ $\displaystyle=(\Phi,\
\delta^{c}\Omega)_{V}+(\Phi,\ \ast^{-1}(\ast\Omega\cup\ast\ast
A))_{V}+(-1)^{r+1}(\Phi,\ \ast^{-1}(A\cup\ast\Omega))_{V}$
$\displaystyle=(\Phi,\ \delta^{c}\Omega+\ast^{-1}(\ast\Omega\cup\ast\ast
A+(-1)^{r+1}A\cup\ast\Omega))_{V}=(\Phi,\ \delta_{A}^{c}\Omega)_{V}.$

∎

Let us consider the following operator

$\Delta_{A}^{c}=d_{A}^{c}\delta_{A}^{c}+\delta_{A}^{c}d_{A}^{c}:H^{r}(V)\rightarrow
H^{r}(V).$

This operator serves as a discrete analogue of the Laplace-type operator
(1.5).

######  Proposition 3.6.

For any $r$-form $\Phi\in H^{r}(V)$, we have

$(\Delta_{A}^{c}\Phi,\ \Phi)_{V}\geq 0$

and $(\Delta_{A}^{c}\Phi,\ \Phi)_{V}=0$ if and only if $\Delta_{A}^{c}\Phi=0$.

###### Proof.

By (3.10), one has

$\displaystyle(\Delta_{A}^{c}\Phi,\ \Phi)_{V}=(d_{A}^{c}\delta_{A}^{c}\Phi,\
\Phi)_{V}+(\delta_{A}^{c}d_{A}^{c}\Phi,\ \Phi)_{V}$
$\displaystyle=(\delta_{A}^{c}\Phi,\ \delta_{A}^{c}\Phi)_{V}+(d_{A}^{c}\Phi,\
d_{A}^{c}\Phi)_{V}=\|\delta_{A}^{c}\Phi\|^{2}+\|d_{A}^{c}\Phi\|^{2},$

where $\|\cdot\|$ is the norm given by (2.12). From this, if
$(\Delta_{A}^{c}\Phi,\ \Phi)_{V}=0$, then $\|\delta_{A}^{c}\Phi\|^{2}=0$ and
$\|d_{A}^{c}\Phi\|^{2}=0$. It gives $\delta_{A}^{c}\Phi=0$ and
$d_{A}^{c}\Phi=0$. Hence,

$\Delta_{A}^{c}\Phi=d^{c}_{A}\delta_{A}^{c}\Phi+\delta_{A}^{c}d_{A}^{c}\Phi=0.$

∎

######  Corollary 3.7.

$\Delta_{A}^{c}\Phi=0$ if and only if $\delta_{A}^{c}\Phi=0$ and
$d_{A}^{c}\Phi=0$.

######  Proposition 3.8.

The operator $\Delta_{A}^{c}:H^{r}(V)\rightarrow H^{r}(V)$ is self-adjoint,
i.e.,

$(\Delta_{A}^{c}\Phi,\ \Omega)_{V}=(\Phi,\ \Delta_{A}^{c}\Omega)_{V}.$

###### Proof.

By (3.10), it is obvious. ∎

######  Remark 3.9.

Taking into account (2.14), the operator $\delta_{A}^{c}\Phi$ can be expressed
as

$\delta_{A}^{c}\Phi=(-1)^{r}\ast^{-1}(d^{c}\ast\Phi+(-1)^{r}\ast\Phi\cup\ast\ast
A+A\cup\ast\Phi),$ (3.12)

where $\Phi\in H^{r}(V)$. On the other hand, by (3.4), it is straightforward
to calculate that

$\delta_{A}^{c}\Phi=(-1)^{r}\ast^{-1}d_{A}^{c}\ast\Phi+\ast^{-1}(\ast\Phi\cup(\ast\ast
A+A)).$ (3.13)

It’s worth noting that in the continuous case, we have
$\delta_{A}=(-1)^{r}\ast^{-1}d_{A}\ast$. However, this doesn’t hold true in
our discrete model. More precisely, it follows from (2.7) that the operation
$\ast\ast$ corresponds to a shift, meaning $\ast\ast A$ takes the form (3.9),
unlike in the continuous case where $\ast\ast A=-A$. Consequently, the
discrete Yang-Mills equations (3.5) does not follow immediately from the
equation $\Delta_{A}^{c}F=0$.

The discrete counterpart of Equation (1.3) has the form

$\delta_{A}^{c}F=0,$ (3.14)

where $F$ is given by (3.2). Again, using (2.16) and (3.9), from (3.11),
Equation (3.14) can be represented in terms of difference equations as follows

$F_{k,\sigma s}-F_{\sigma k,\sigma s}+A^{1}_{k,s}F_{k,\sigma s}-F_{\sigma
k,\sigma s}A^{1}_{\sigma k,\sigma s}=0,$ $F_{\sigma k,s}-F_{\sigma k,\sigma
s}+A^{2}_{k,s}F_{\sigma k,s}-F_{\sigma k,\sigma s}A^{2}_{\sigma k,\sigma s}=0$

for any $k=1,2,...,N$ and $s=1,2,...,M$. Comparing the difference
representations of Equations (3.5) and (3.14), it can be noted that only the
last terms of the difference equations slightly differ in both cases.

##  4 Discrete Yang-Mills equations on a combinatorial torus

In this section, following [22], we delve into an example of our discrete
model on a combinatorial torus in detail. Recall that the torus can be
regarded as the topological space obtained by taking a rectangle and
identifying each pair of opposite sides with the same orientation. As was
already discussed in [22], let’s consider the partitioning of the plane
$\mathbb{R}^{2}$ by the straight lines $x=k$ and $y=s$, where
$k,s\in\mathbb{Z}$. We denote an open square bounded by these lines as
$V_{k,s}$, with vertices labeled as $x_{k,s},\ x_{\tau k,s},\ x_{k,\tau s}$,
$x_{\tau k,\tau s}$, where $\tau k=k+1$. Further, we define $e_{k,s}^{1}$ and
$e_{k,s}^{2}$ as the open intervals $(x_{k,s},\ x_{\tau k,s})$ and $(x_{k,s},\
x_{k,\tau s})$, respectively. These geometric objects can be associated with
the combinatorial objects we have previously discussed. We identify the
collection $V_{k,s}$ with $V$ given by (2.9) and let $N=M=2$. If we identify
the points and the intervals on the boundary of $V$ as follows

$\displaystyle x_{1,1}=x_{3,1}=x_{1,3}=x_{3,3},\qquad x_{1,2}=x_{3,2},\qquad
x_{2,1}=x_{2,3},$ $\displaystyle e_{1,1}^{1}=e_{1,3}^{1},\qquad
e_{2,1}^{1}=e_{2,3}^{1},\qquad e_{1,1}^{2}=e_{3,1}^{2},\qquad
e_{1,2}^{2}=e_{3,2}^{2},$ (4.1)

we obtain the geometric object which is homeomorphic to the torus. For a
visual representation, see [22, Figure 1]. Let $C(T)$ denote the complex
corresponding to this geometric object, serving as a combinatorial model of
the torus. Similarly, we denote by $K(T)$ the complex of cochains over $C(T)$.
From (4), it follows that forms on the combinatorial torus satisfy the
conditions (2).

Let us now proceed to the construction of the discrete Yang-Mills equations on
the combinatorial torus. The discrete connection 1-form $A\in K^{1}(T)$ has
the form

$A=\sum_{k=1}^{2}\sum_{s=1}^{2}(A^{1}_{k,s}e_{1}^{k,s}+A^{2}_{k,s}e_{2}^{k,s}).$

Using the definition of the $\cup$ product and (2.5), we have

$\displaystyle A\cup
A=(A^{1}_{1,1}A^{2}_{2,1}-A^{2}_{1,1}A^{1}_{1,2})V^{1,1}+(A^{1}_{2,1}A^{2}_{1,1}-A^{2}_{2,1}A^{1}_{2,2})V^{2,1}$
$\displaystyle+(A^{1}_{1,2}A^{2}_{2,2}-A^{2}_{1,2}A^{1}_{1,1})V^{1,2}+(A^{1}_{2,2}A^{2}_{1,2}-A^{2}_{2,2}A^{1}_{2,1})V^{2,2}$

and

$\displaystyle
d^{c}A=(A^{1}_{1,1}-A^{1}_{1,2}+A^{2}_{2,1}-A^{2}_{1,1})V^{1,1}+(A^{1}_{2,1}-A^{1}_{2,2}-A^{2}_{2,1}+A^{2}_{1,1})V^{2,1}$
$\displaystyle+(A^{1}_{1,2}-A^{1}_{1,1}+A^{2}_{2,2}-A^{2}_{1,2})V^{1,2}+(A^{1}_{2,2}-A^{1}_{2,1}+A^{2}_{1,2}-A^{2}_{2,2})V^{2,2}.$

Then the discrete curvature form (3.2) on $C(T)$ can be represented as

$F=F_{1,1}V^{1,1}+F_{2,1}V^{2,1}+F_{1,2}V^{1,2}+F_{2,2}V^{2,2},$

where

$\displaystyle
F_{1,1}=A^{1}_{1,1}-A^{1}_{1,2}+A^{2}_{2,1}-A^{2}_{1,1}+A^{1}_{1,1}A^{2}_{2,1}-A^{2}_{1,1}A^{1}_{1,2},$
$\displaystyle
F_{2,1}=A^{1}_{2,1}-A^{1}_{2,2}-A^{2}_{2,1}+A^{2}_{1,1}+A^{1}_{2,1}A^{2}_{1,1}-A^{2}_{2,1}A^{1}_{2,2},$
$\displaystyle
F_{1,2}=A^{1}_{1,2}-A^{1}_{1,1}+A^{2}_{2,2}-A^{2}_{1,2}+A^{1}_{1,2}A^{2}_{2,2}-A^{2}_{1,2}A^{1}_{1,1},$
$\displaystyle
F_{2,2}=A^{1}_{2,2}-A^{1}_{2,1}+A^{2}_{1,2}-A^{2}_{2,2}+A^{1}_{2,2}A^{2}_{1,2}-A^{2}_{2,2}A^{1}_{2,1}.$

By (2.7) and (4), we obtain

$\ast F=F_{2,2}x^{1,1}+F_{1,2}x^{2,1}+F_{2,1}x^{1,2}+F_{1,1}x^{2,2}.$

Consequently, we can calculate the following:

$\displaystyle d^{c}\ast
F=(F_{1,2}-F_{2,2})e_{1}^{1,1}+(F_{2,2}-F_{1,2})e_{1}^{2,1}+(F_{2,2}-F_{2,1})e_{2}^{1,2}$
$\displaystyle+(F_{2,1}-F_{2,2})e_{2}^{1,1}+(F_{1,1}-F_{2,1})e_{1}^{1,2}+(F_{2,1}-F_{1,1})e_{1}^{2,2}+$
$\displaystyle+(F_{1,2}-F_{1,1})e_{2}^{2,2}+(F_{1,1}-F_{1,2})e_{2}^{2,1},$
(4.2) $\displaystyle\ast F\cup
A=F_{2,2}A^{1}_{1,1}e_{1}^{1,1}+F_{1,2}A^{1}_{2,1}e_{1}^{2,1}+F_{2,1}A^{2}_{1,2}e_{2}^{1,2}+F_{2,2}A^{2}_{1,1}e_{2}^{1,1}$
$\displaystyle+F_{2,1}A^{1}_{1,2}e_{1}^{1,2}+F_{1,1}A^{1}_{2,2}e_{1}^{2,2}+F_{1,1}A^{2}_{2,2}e_{2}^{2,2}+F_{1,2}A^{2}_{2,1}e_{2}^{2,1},$
(4.3)

and

$\displaystyle A\cup\ast
F=A^{1}_{1,1}F_{1,2}e_{1}^{1,1}+A^{1}_{2,1}F_{2,2}e_{1}^{2,1}+A^{2}_{1,2}F_{2,2}e_{2}^{1,2}+A^{2}_{1,1}F_{2,1}e_{2}^{1,1}$
$\displaystyle+A^{1}_{1,2}F_{1,1}e_{1}^{1,2}+A^{1}_{2,2}F_{2,1}e_{1}^{2,2}+A^{2}_{2,2}F_{1,2}e_{2}^{2,2}+A^{2}_{2,1}F_{1,1}e_{2}^{2,1}.$
(4.4)

Combining (4), (4), and (4), we obtain the representation of $d_{A}^{c}\ast F$
on the combinatorial torus. Then Equation (3.5) can be expressed in terms of
difference equations as

$\displaystyle\langle e^{1}_{1,1},\ d_{A}^{c}\ast
F\rangle=F_{1,2}-F_{2,2}+A^{1}_{1,1}F_{1,2}-F_{2,2}A^{1}_{1,1}=0,$
$\displaystyle\langle e^{1}_{2,1},\ d_{A}^{c}\ast
F\rangle=F_{2,2}-F_{1,2}+A^{1}_{2,1}F_{2,2}-F_{1,2}A^{1}_{2,1}=0,$
$\displaystyle\langle e^{2}_{1,2},\ d_{A}^{c}\ast
F\rangle=F_{2,2}-F_{2,1}+A^{2}_{1,2}F_{2,2}-F_{2,1}A^{2}_{1,2}=0,$
$\displaystyle\langle e^{2}_{1,1},\ d_{A}^{c}\ast
F\rangle=F_{2,1}-F_{2,2}+A^{2}_{1,1}F_{2,1}-F_{2,2}A^{2}_{1,1}=0,$
$\displaystyle\langle e^{1}_{1,2},\ d_{A}^{c}\ast
F\rangle=F_{1,1}-F_{2,1}+A^{1}_{1,2}F_{1,1}-F_{2,1}A^{1}_{1,2}=0,$
$\displaystyle\langle e^{1}_{2,2},\ d_{A}^{c}\ast
F\rangle=F_{2,1}-F_{1,1}+A^{1}_{2,2}F_{2,1}-F_{1,1}A^{1}_{2,2}=0,$
$\displaystyle\langle e^{2}_{2,2},\ d_{A}^{c}\ast
F\rangle=F_{1,2}-F_{1,1}+A^{2}_{2,2}F_{1,2}-F_{1,1}A^{2}_{2,2}=0,$
$\displaystyle\langle e^{2}_{2,1},\ d_{A}^{c}\ast
F\rangle=F_{1,1}-F_{1,2}+A^{2}_{2,1}F_{1,1}-F_{1,2}A^{2}_{2,1}=0.$

Similarly, taking into account (3.13), Equation (3.14) is equivalent to the
following system of difference equations

$\displaystyle F_{1,2}-F_{2,2}+A^{1}_{1,1}F_{1,2}-F_{2,2}A^{1}_{2,2}=0,$
$\displaystyle F_{2,2}-F_{1,2}+A^{1}_{2,1}F_{2,2}-F_{1,2}A^{1}_{1,2}=0,$
$\displaystyle F_{2,2}-F_{2,1}+A^{2}_{1,2}F_{2,2}-F_{2,1}A^{2}_{2,1}=0,$
$\displaystyle F_{2,1}-F_{2,2}+A^{2}_{1,1}F_{2,1}-F_{2,2}A^{2}_{2,2}=0,$
$\displaystyle F_{1,1}-F_{2,1}+A^{1}_{1,2}F_{1,1}-F_{2,1}A^{1}_{2,1}=0,$
$\displaystyle F_{2,1}-F_{1,1}+A^{1}_{2,2}F_{2,1}-F_{1,1}A^{1}_{1,1}=0,$
$\displaystyle F_{1,2}-F_{1,1}+A^{2}_{2,2}F_{1,2}-F_{1,1}A^{2}_{1,1}=0,$
$\displaystyle F_{1,1}-F_{1,2}+A^{2}_{2,1}F_{1,1}-F_{1,2}A^{2}_{1,2}=0.$ (4.5)

Finally, we give a matrix form of Equations (3.5) and (3.14) on the
combinatorial torus. Let’s introduce the following row vectors

$[F]=[F_{1,1}\ F_{2,1}\ F_{1,2}\ F_{2,2}],\quad[A]=[A^{1}_{1,1}\ A^{1}_{2,1}\
A^{2}_{1,2}\ A^{2}_{1,1}\ A^{1}_{1,2}\ A^{1}_{2,2}\ A^{2}_{2,2}\
A^{2}_{2,1}],$ $[x]=[x^{1,1}\ x^{2,1}\ x^{1,2}\
x^{2,2}],\quad[e]=[e_{1}^{1,1}\ e_{1}^{2,1}\ e_{2}^{1,2}\ e_{2}^{1,1}\
e_{1}^{1,2}\ e_{1}^{2,2}\ e_{2}^{2,2}\ e_{2}^{2,1}],$ $[V]=[V^{1,1}\ V^{2,1}\
V^{1,2}\ V^{2,2}].$

Denote by $[\cdot]^{T}$ the corresponding column vector. In this notation, the
1-form $A$ can be rewritten as

$A=[e][A]^{T}=[A][e]^{T}.$

By trivial computation, one finds that

$[\ast F]=[F]S,\quad\ast F=[F]S[x]^{T}$

and

$d^{c}\ast F=[e]DS[F]^{T},\quad\quad\ast F\cup A=[F]SD_{2}I_{A}[e]^{T},\quad
A\cup\ast F=[e]I_{A}D_{1}S[F]^{T},$

where

$D=\begin{bmatrix}-1&1&0&0\\\ 1&-1&0&0\\\ 1&0&-1&0\\\ -1&0&1&0\\\ 0&0&-1&1\\\
0&0&1&-1\\\ 0&1&0&-1\\\ 0&-1&0&1\end{bmatrix},\quad
S=\begin{bmatrix}0&0&0&1\\\ 0&0&1&0\\\ 0&1&0&0\\\ 1&0&0&0\end{bmatrix},$
$D_{1}=\begin{bmatrix}0&1&0&0\\\ 1&0&0&0\\\ 1&0&0&0\\\ 0&0&1&0\\\ 0&0&0&1\\\
0&0&1&0\\\ 0&1&0&0\\\ 0&0&0&1\end{bmatrix},\quad
D_{2}=\begin{bmatrix}1&0&0&1&0&0&0&0\\\ 0&1&0&0&0&0&0&1\\\ 0&0&1&0&1&0&0&0\\\
0&0&0&0&0&1&1&0\end{bmatrix},$
$I_{A}=\begin{bmatrix}A^{1}_{1,1}&0&0&0&0&0&0&0\\\
0&A^{1}_{2,1}&0&0&0&0&0&0\\\ 0&0&A^{2}_{1,2}&0&0&0&0&0\\\
0&0&0&A^{2}_{1,1}&0&0&0&0\\\ 0&0&0&0&A^{1}_{1,2}&0&0&0\\\
0&0&0&0&0&A^{1}_{2,2}&0&0\\\ 0&0&0&0&0&0&A^{2}_{2,2}&0\\\
0&0&0&0&0&0&0&A^{2}_{2,1}\end{bmatrix}.$

Thus, the discrete Yang-Mills equation (3.5) can be written in matrix form as

$[e]DS[F]^{T}+[e]I_{A}D_{1}S[F]^{T}-[F]SD_{2}I_{A}[e]^{T}=0$

or, equivalently,

$DS[F]^{T}+I_{A}D_{1}S[F]^{T}-([F]SD_{2}I_{A})^{T}=[0]^{T}.$

Similarly, Equation (3.14) or the system (4) can be represented as follows

$DS[F]^{T}+I_{A}D_{1}S[F]^{T}+([F]SD_{2}I_{\ast\ast A})^{T}=[0]^{T},$

where

$[\ast\ast A]=[-A^{1}_{2,2}\ -A^{1}_{1,2}\ -A^{2}_{2,1}\ -A^{2}_{2,2}\
-A^{1}_{2,1}\ -A^{1}_{1,1}\ -A^{2}_{1,1}\ -A^{2}_{1,2}].$

##  5 Conclusions

This study introduces a discrete model of the two-dimensional Yang-Mills
equations that preserves geometric structure. The properties of operators
generated by the discrete Yang-Mills equations have been thoroughly
investigated. Special attention is given to discrete models on the
combinatorial two-dimensional torus. In this context, a system of difference
equations equivalent to the discrete Yang-Mills equations was constructed, and
the matrix form of these equations was also presented. The author encourages
readers to apply the proposed discretization scheme to develop numerical
implementation techniques.

## References

  * [1] Arnold, D. Finite Element Exterior Calculus; SIAM, 2018.
  * [2] Atiyah, M.F.; Bott, R. The Yang-Mills equations over Riemann surfaces. Philos. Trans. R. Soc. Lond., Ser. A 1983, 308, 523–615.
  * [3] Aharony, O.; Kundu, S.; Sheaffer, T. A string theory for two dimensional Yang-Mills theory I. arXiv 2024, arXiv:2312.12266v2.
  * [4] Berchenko-Kogan, Y.; Stern, A. Charge-conserving hybrid methods for the Yang-Mills equations. SMAI J. Comput. Math. 2021, 7, 97–119.
  * [5] Catterall, S. A geometrical approach to $N=2$ super Yang-Mills theory on the two dimensional lattice. J. High Energy Phys. 2004, 11, 006.
  * [6] Catterall, S.; Jha, R. G.; Schaich, D.; Wiseman, T. Testing holography using lattice super-Yang–Mills theory on a 2-torus. Phys. Rev. D 2018, 97, 086020.
  * [7] Christiansen, S. H.; Winther, R. On constraint preservation in numerical simulations of Yang-Mills equations. SIAM J. Sci. Comput. 2006, 28(1), 75–101.
  * [8] Christiansen, S. H.; Halvorsen, T. G. A simplicial gauge theory. J. Math. Phys. 2012, 53(3).
  * [9] Driver, B. K.; Gabriel, F.; Hall, B. C.; Kemp, T. The Makeenko-Migdal equation for Yang-Mills theory on compact surfaces. Comm. Math. Phys. 2017, 352(3), 967–978.
  * [10] Droniou, J.; Oliynyk, T. A.; Qian, J. J. A polyhedral discrete de Rham numerical scheme for the Yang–Mills equations. Journal of Computational Physics. 2023, 478, 111955.
  * [11] Dezin, A.A. Multidimensional analysis and discrete models; CRC Press, Boca Raton, 1995.
  * [12] Freed, D.; Uhlenbeck, K. Instantons and Four-Manifolds; Springer–Verlag, 1984.
  * [13] Gross, A.; King, C.; Sengupta, A. Two dimensional Yang-Mills theory via stochastic differential equations. Annals of Physics 1989, 194, 65-112.
  * [14] Joseph, A. Supersymmetric Yang–Mills theories with exact supersymmetry on the lattice. Int. J. Mod. Phys. A 2011, 26, 5057.
  * [15] Kondo, K.-I.; Shibata, A.; Shinohara, T.; Murakami, T.; Kato, S.; Ito, S. New descriptions of lattice SU(N) Yang–Mills theory towards quark confinement. Phys. Lett. B 2008, 669, 107–118.
  * [16] Hanada, M.; Kanamori, I. Lattice study of two-dimensional $N(2,2)$ super Yang-Mills at large-$N$. Phys. Rev. D 2009, 80, 065014.
  * [17] Sushch, V. Gauge-invariant discrete models of Yang-Mills equations. Math. Notes 1997, 61(5), 621–631.
  * [18] Sushch, V. Discrete model of Yang-Mills equations in Minkowski space. Cubo 2004, 6(2), 35–50.
  * [19] Sushch, V. A gauge-invariant discrete analog of the Yang-Mills equations on a double complex. Cubo 2006, 8(3), 61–78.
  * [20] Sushch, V. Self-dual and anti-self-dual solutions of discrete Yang-Mills equations on a double complex. Cubo 2010, 12(3), 99–120.
  * [21] Sushch, V. Instanton-anti-instanton solutions of discrete Yang-Mills equations. Math. Bohemica 2012, 137(2), 219–228.
  * [22] Sushch, V. 2D discrete Hodge–Dirac operator on the torus. Symmetry 2022, 14, 1556.
  * [23] Witten, E. On quantum gauge theories in two dimensions.Comm. Math. Phys. 1991, 141(1), 153–209.