# Investigating techniques to optimise the layout of turbines in a windfarm
using a quantum computer

James Hancock  Centre for Mathematical Sciences, University of Plymouth
<EMAIL_ADDRESS>Matthew J. Craven  Centre for Mathematical
Sciences, University of Plymouth  Craig McNeile  Centre for Mathematical
Sciences, University of Plymouth  Davide Vadacchino  Centre for Mathematical
Sciences, University of Plymouth

###### Abstract

We study the optimal placement of wind turbines within a windfarm to maximize
the power produced by mapping the system to a Quadratic Unconstrained Binary
Optimisation (QUBO) problem. We investigate solving the resulting QUBO problem
using the Variational Quantum Eigensolver (VQE) on a quantum computer
simulator and compare the results to those from two classical optimisation
methods: simulated annealing and the Gurobi solver. The maximum grid size we
study is 4 $\times$ 4, which requires 16 qubits.

##  1 Introduction

With the climate crisis and goals to reduce our emissions, there is a clear
motivation for extracting the maximum amount of energy from renewable sources.
The use of renewable energy sources will minimize our negative effect on the
environment and lead to a more sustainable future.

One key renewable energy source is wind. There are many challenges to making
it an economically viable and reliable energy source [1]. In this paper, we
focus on using quantum computers to solve one specific problem towards
maximizing the energy extracted from the wind. There are currently limitations
on the available quantum computing hardware that we will discuss later. There
have been a number of reviews of the use of quantum computing to solve
problems in renewable energy [2, 3, 4]. There has been a recent review on
using quantum computers to solve optimisation problems [5].

Energy from the wind is extracted using wind turbines. Typically, wind
turbines are combined together to form windfarms. The spatial layout of the
turbines in the windfarm will change the maximum amount of energy that the
windfarm can produce. There is a long history of using optimization techniques
[6, 7] to find the best locations of the turbines in the windfarm to produce
the maximum amount of power.

The number of turbines in the windfarms increases to maximize the power
produced. For example, China is building a windfarm with 3500 turbines [8]. As
the number of turbines increases, the complexity of the optimization problem
grows exponentially. This motivates the study of improved optimization
algorithms for finding the best layout of turbines in a windfarm. Recently,
one form of the Windfarm Layout Optimization (WFLO) problem has been rewritten
as a Quadratic Unconstrained Binary Optimization (QUBO) problem [9] that can
be solved on a quantum computer.

There are currently two dominant types of quantum computers. One is based on
quantum circuits. For example, IBM and Google have constructed quantum
computers based on this paradigm. The current largest real-world quantum
computer is IBM’s Osprey [10], which has 433 qubits. One way to solve QUBO
problems using a quantum circuit computer is to use a Variational Quantum
Eigensolver (VQE).

Other types of quantum computers are quantum annealers (QA) [11]. The problem
that quantum annealers can solve are QUBO problems. D-WAVE’s Advantage is the
largest real-world QA, having 5000+ (although not all connected) qubits [12].
A new system in the coming years which promises to have 7000+ qubits [13] is
currently under development. Fujitsu’s Digital Annealer can also solve QUBO
problems using specialized hardware [14].

Another important practical issue is that quantum error correction is not
implemented on existing quantum computers. We are currently in the noisy
intermediate-scale quantum (NISQ) [6] era. Any algorithm should be resilient
to decoherence. Practical algorithms in the NISQ era typically run on hybrid
systems made from a combination of classical and quantum computers. We will
postpone the study of the effect of the additional problems, such as
decoherence and quantum errors from running on real quantum hardware to later
work.

In this paper, we investigate using quantum optimisation algorithms to solve
the WFLO problem mapped to a QUBO problem using the qiskit package [15] from
IBM, running on a classical simulator. The solution method uses the
Variational Quantum Eigensolver (VQE), so it is an appropriate method for the
noisy intermediate-scale quantum (NISQ) era of quantum computing [16]. This
work is a necessary first step to solving the WFLO problem using quantum
computers. The next steps are studying the introduction of quantum noise and
errors in the simulator, and finally using a real quantum computer.

Senderovich et al. [9] developed the formalism to convert a WFLO problem into
a QUBO problem. They compared solving the resulting QUBO problem using
classical methods and the Fujitsu Digital Annealer. Our contribution is to
investigate solving the QUBO problem mapped from a WFLO using a simulator of a
circuit-based quantum computer.

Apart from the importance of solving the WFLO problem in order to produce more
electricity from a windfarm, there is another motivation for this study. QUBO
problems are one of the main problems that can be solved by quantum annealers
such as those sold by D-Wave and Fujitsu [17]. Once our project is running on
quantum computing hardware, the solution of WFLO problems may be an
interesting test case for comparing the performance of adiabatic and circuit-
based quantum computers. Previous work has shown the difficulty in comparing
the performance of classical computer and adiabatic annealers in solving QUBO
problems [18, 19, 20, 21].

In this work we have used the circuit model of quantum computing as the basis
of the calculation. The next step in the project will be to run the code on
quantum computers and we will then be able to compare the performance on
circuit and quantum-based systems.

##  2 Windfarm Layout Optimisation

Choosing the placement of turbines within a windfarm, or Windfarm Layout
Optimisation (WFLO) [6] is selecting the locations of the turbines in order to
maximize the power output, with respect to certain constraints. These
constraints include the maximum number of turbines we can use, as well as the
minimum distance that we must have between them due to their rotors. We use
the simplifying assumption that the turbines can be placed only on a discrete
grid.

When a turbine is downstream (meaning that it is behind the other turbine with
respect to the wind direction) from another within its cone, it will have
reduced power output compared to if it is out in the open. This effect is
called the wake of a turbine. The key issue in the placement of the wind
turbines in the windfarm is to minimize the lost energy due to the wakes from
the turbines.

The Sum-of-Squares (SS) model for wake speeds [22] best captures the effect of
these wakes. However, the formulation leads to intractable optimisation
problems. To combat this, we examine using a simplified model, the Linear
Superposition of Wakes (LS) [23]. This allows us to map the problem to a
Quadratic Unconstrained Binary Optimisation (QUBO) problem, which can then in
turn be optimized using various quantum (and classical) methods.

##  3 Windfarm layout as a Quadratic Unconstrained Binary Optimisation problem

The Quadratic Unconstrained Binary Optimisation (QUBO) problem is defined as:

$\underset{x}{\text{argmin }}f_{Q}(x)\text{,}$ (1)

where,

$f_{Q}(x):=x^{T}Qx=\sum_{i,j=0}^{q}Q_{ij}x^{i}x^{j}$ (2)

$x\in\mathbb{B}^{q}=\\{0,1\\}^{q}$.

Many combinatorial optimization problems can be formulated as QUBO problems.
For example, the reconstruction of the tracks of charged particles can be
computed by solving a QUBO problem [24]. QUBO is a class of NP-hard problems
[25]. It is not expected that the solution of QUBO problems will be in the BQP
complexity class and hence that the problem can be solved in polynomial time
on quantum computers. QUBO problems are normally formulated in an
unconstrained way. However, constraints can be included by introducing
penalties in the form of large additive constants when the conditions are
violated. We follow the mapping of the WFLO problem to a QUBO problem that was
developed by A. Senderovich, et al in Ref. [9].

We will simplify the land that the farm is on to a square grid of length
$l_{grid}$. There are in total $l_{grid}^{2}$ total spots on which a turbine
could be placed. These will be labelled as in Fig. 1. The indicator of the
presence of a turbine on site $q$ is $x^{q}$ where
$q\in\\{1,...l_{grid}^{2}\\}$. Cartesian coordinates of the grid are shown in
Fig. 2.

There are two constraints that we impose on this simplified model of a
windfarm, which are:

  1. 1.

There can be no more than $m$ turbines on the grid.

  2. 2.

Each turbine must be at least distance $\xi$ apart from one another.

15913261014371115481216 Figure 1: Labelling of sites on a $l_{grid}=4$
windfarm grid.
(1,1)(2,1)(3,1)(4,1)(1,2)(2,2)(3,2)(4,2)(1,3)(2,3)(3,3)(4,3)(1,4)(2,4)(3,4)(4,4)
Figure 2: Coordinates of sites on a $l_{grid}=4$ windfarm grid.

A wake can be defined by three parameters, the angle $\alpha$ (degrees) that
the wind is coming from with respect to the west counted clockwise, $x$, the
maximum distance that is affected by the wake and $r$ is the radius with which
it spreads out per unit distance, away from the turbine causing the wake.
These physical factors are problem-specific. For our model, we say that if a
grid location is at all within the wake, it is entirely within its wake. We
denote a wake starting at position $i$ as $w(i,x,r;\alpha)$.

Figure 3: Turbine on location 15 or (1,4), $D=\\{0,12ms^{-1},1\\}$ Figure 4:
Turbine on location 23 or (2,2) and 34 or (3,3), $D=\\{0,12ms^{-1},1\\}$
Figure 5: Turbine on location 27 or (2,6), $D=$ Second Mosetti. Figure 6:
Turbine on location 35 or (3,4) and 76 or (7,5), $D=$ Second Mosetti.

To visualize the wakes from some turbines in the wind, they are plotted on a
$l_{grid}=10$ grid (that is bigger than used in the simulations) with the
parameters: $x=3$, $r=1$. Fig. 3 and Fig. 5 show a configuration of the wakes
from one turbine. Fig. 4 and Fig. 6 show a configuration of the wakes from two
turbines. If one of the turbines is in the wake of another turbine, then it
will see a reduced wind velocity and thus produce less power. The algorithm’s
goal is to place wind turbines in locations that have the windspeeds
maximized.

###  3.1 Model

To ensure that the model is more realistic, the system includes several wind
regimes. Each possible wind regime is defined as,
$d:=\\{\alpha_{d},v_{d},p_{d}\\}$. $\alpha_{d}$ is the angle at which the wind
comes in, $v_{d}$ is the free wind speed, i.e. the wind speed that powers a
turbine not in the wake of another. $p_{d}$ is the probability that we would
expect to encounter this wind regime. We denote the set of possible wind
regimes in our system as $D$. To clarify this, $D$ is a collection of $d$’s,
with $\sum_{d\in D}p_{d}=1$

The particular model that we are using is the Linear Superposition (LS) of
wakes. This model is useful as it can be transformed into a QUBO problem. The
power output for a system of $q$ sites within this model is computed as:

$E^{\text{LS}}=\sum_{d\in
D}\sum_{i=1}^{q}p_{d}\left[\frac{1}{3}v_{d}^{3}-\sum_{j\in
w_{i}}\frac{1}{3}\left(v_{d}^{3}-u_{ij}^{3}\right)\right],$ (3)

where $w_{i}$ is shorthand for the wake at $i$, and $u_{ij}$ is the reduced
windspeed at site $j$, caused by being in the wake of $i$. If $j$ is in the
wake of $i$, $u_{ij}$ can be calculated as in,

$u_{ij}=v_{d}\left(1-\frac{2a}{\left(1+\alpha_{T}(\delta/r)^{2}\right)^{2}}\right)$
(4)

this is the Jansen wake model [26], where $a$ in is the axial induction
factor, currently set to 0.1, and $\delta=||i-j||_{2}$. When we control the
distance which the wake can go back ($x$), as well as its radius per distance
($r$), we can calculate $\alpha_{T}$ using Eq. 5.

$\alpha_{T}=\frac{1}{x}\left(r-r_{t}\right),$ (5)

where $r_{t}$ is the turbine radius which must be set. For this work, we set
$r_{t}=0.33$ [27]. We note that $r_{t}$ and $a$ from Eq. 4 carry units, and
change in real systems dependent on the size and energy output of the
turbines. The values which have been chosen are arbitrary, they do not relate
to any real-life system - this is a simplified model.

The QUBO formulation for this problem is then:

$\underset{x}{\text{argmin }}(-f(x))$ (6) $f(x)=\sum_{d\in
D}\sum_{i=1}^{m}p_{d}\left[\frac{1}{3}v_{d}^{3}x^{i}x^{i}-\sum_{j\in
w_{i}}\frac{1}{3}\left(v_{d}^{3}-u_{ij}^{3}\right)x^{i}x^{j}\right]$ (7)

This is the unconstrained version of the problem, i.e. the optimal solution
will always be a turbine on every site. We now must add the constraints in the
form of an energetic penalty.

$g(x;\lambda_{1},\lambda_{2})=\lambda_{1}\left(\sum_{i=1}^{l_{grid}^{2}}x^{i}-m\right)^{2}+\lambda_{2}\sum_{||i-j||_{2}<\xi}x^{i}x^{j}$
(8)

The term which is multiplied by $\lambda_{1}$ in Eq. 8 limits the number of
total turbines and the second term pertains to the minimum distance. The
values of $\lambda_{1},\lambda_{2}$ have to be large enough so that the
constraints are met. The full QUBO problem is:

$\underset{x}{\text{argmin }}(-f(x)+g(x;\lambda_{1},\lambda_{2}))$ (9)

The problem from Eq. 9 can then written in the form of a weight matrix Q from
Eq. 2:

$Q_{ij}=\begin{cases}-\frac{1}{3}\sum_{d}p_{d}v_{d}^{3}+\lambda_{1}(1-2nm)&\text{if
$i=j$}\\\ \\\
-\frac{1}{3}\sum_{d}p_{d}\left(v_{d}^{3}-u_{ij}^{3}\right)+2\lambda_{1}&\text{if
$i\neq j$}\end{cases}$ (10)

Here, we have neglected the proximity constraints term as we do not include
this in our simulations. This can often lead to dense matrices, with all
entries nonzero. We note here that $q=l_{grid}^{2}$ is the number of
variables/qubits.

##  4 Classical methods

As a baseline comparison to the use of quantum computers to solve QUBO
problems we also investigated two classical optimization techniques.

Gurobi is a set of solver packages that can be used in many different
programming languages to solve linear and quadratic optimization problems. The
package uses several different optimization techniques, choosing the best
depending on the problem at hand. Examples of the algorithms that are made use
of are: simplex and parallel barrier, further details can be found in Ref.
[28]. In this work, we used Gurobi to solve our QUBO problem, which is a
special case of quadratic optimization. Gurobi has been used as a comparison
to using quantum and digital annealers to solve different types of QUBO
problems [29, 30].

Simulated Annealing [31] (SA) is a standard approach to finding the optimum of
a function which has many local minima in a huge search space. It is
particularly useful for searching in discrete spaces [32]. The algorithm
searches the space via random fluctuations, controlled by a probability to
transfer between states which is a function of temperature. SA has been used
to find solutions to various NP-complete combinatorial problems such as the
Travelling Salesman Problem (TSP) [33], Minimum Linear Arrangement (MLA) [34]
and instances of the Packing Problem (PP) [35].

##  5 Variational Quantum Eigensolver

The Variational Quantum Eigensolver (VQE) is a hybrid quantum-classical
algorithm for finding the lowest eigenvalue of a Hamiltonian that is in the
form of a Pauli string [36]. For further details about the VQE see Ref. [37].
The VQE is an algorithm that has shown great promise in making use of current
NISQ computers. It uses a small parameterized quantum circuit to calculate the
expected value of the Hamiltonian with respect to the parameters ($\theta$). A
classical optimization routine is then used in order to best choose $\theta$
so that we minimize this value.

The VQE uses the inequality,

$E_{min}^{H}\geq\langle 0|U^{\dagger}(\theta)HU(\theta)|0\rangle,$ (11)

and the variational principle in order to find a tight upper bound on the
lowest eigenvalue, $E_{min}^{H}$, of a Pauli string matrix, $H$ which are
defined as follows,

$H=\sum_{\gamma}h_{\gamma}P_{\gamma}$ (12)
$P_{\gamma}=\bigotimes_{i=0}^{q}\sigma_{m_{i}}$ (13)

Where $m_{i}\in\\{0,3\\}$ tells us the Pauli matrix and $q$ is the number of
qubits. The two Pauli matrices that we make use of in this problem are,

$\sigma_{0}=\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}=I_{2},$ (14)

and

$\sigma_{3}=\begin{pmatrix}1&0\\\ 0&-1\end{pmatrix}$ (15)

The parameter values are then chosen by a classical optimization routine to
minimize this value.

$\underset{\theta}{\text{min }}\langle
0|U^{\dagger}(\theta)HU(\theta)|0\rangle$ (16)

There are two main features of the VQE that we must design and control in
order to make use of it efficiently. These are: the design of the
parameterized circuit (the ansatz) and the choice of optimization routine. Our
choice of ansatz must create states which have sufficient overlap with
solution space in order to be able to find the minimum. When using the VQE to
solve QUBO problems this is fairly trivial as one parameter on each qubit
would be sufficient (we are only looking for basis states). However, the
choice of optimization routine is certainly non-trivial and is vital in order
to best make use of the algorithm. In this work, we examine the effectiveness
of three different optimization routines, Powell’s gradient-free conjugate
search method (PO), Constrained Optimization by Linear Approximation (COByLA)
and Bayesian Optimization (BO).

####  5.0.1 Powell optimization

PO is a widely known standard optimization routine for finding the global
minimum of a function where we do not know/cannot calculate the derivatives.
It works well for the VQE (in small parameter spaces) but requires the system
to use a high shot 111Number of calls of the quantum circuit to gain one
expected value evaluation. count in order to not be too greatly affected by
the quantum measurement noise.

For this work, we make use of scipy.optimize.minimize’s option ‘powell’.

####  5.0.2 COByLA

COByLA is a trust region based, surrogate-assisted method for finding the
global minimum of a function. Similar to PO, COByLA requires sufficient shots
in order to not be too greatly affected by the noise; however, it is more
resilient than PO. COByLA has been used for the VQE several times within the
literature. Liu, et al. [38] used it as a test for their Layer-VQE (L-VQE)
approach to generating ansatz, and Tim Schwägerl, et al. [39] used it to
minimize a QUBO problem for reconstructing particle track results from a Large
Hadron Collider (LHC). COByLA has shown good promise for this black-box
optimization.

For this work, we make use of scipy.optimize.minimize’s option ‘COBYLA’.

####  5.0.3 Bayesian Optimization

BO is a surrogate-assisted method for finding the global minimum of a
function, which models the underlying function as a Gaussian Process
Regression (GPR) and then updates this model based on samples, according to
Bayes’ theorem [40, 41, 42, 43] We follow the work using BO for the VQE [43],
where the authors made use of the Noisy Expected Improvement acquisition
function. BO has also been used for other Variational Quantum Algorithms
(VQAs), such as solving a MAXCUT problem with the Quantum Approximate
Optimization Algorithm (QAOA) [44]. It has also been used for generative
modelling, where combinations of Pauli gates are used rather than
parameterized gates, to implement a training algorithm on the canonical Bars-
and-Stripes data set. Instead of iterations, unique circuits at every step are
used, see Ref. [45]. As well as these quantum examples, BO has been
extensively studied in the field of hyperparameter optimization for machine
learning (e.g. Ref. [46], [47], [48]).

For this work, we use a BO that we have coded. The definitions that we use are
the same as presented in Ref. [43]. It uses the periodic kernel, for $K$
parameters, defined as,

$k^{P}(\theta^{1},\theta^{2})=\sigma^{2}\prod_{i=1}^{K}\exp\left[\frac{-2}{l^{2}}\sin^{2}\left(\pi\left|\frac{\theta^{1}_{i}-\theta^{2}_{i}}{p}\right|^{2}\right)\right],$
(17)

where $\sigma$ is the expected variance of the GP, $l$ is the length scale
over which we would expect data points to be correlated and $p$ is the period
of the underlying function.

We make use of Expected Improvement (EI) as the acquisition function. If we
define our best guess so far to be $E_{min}=\min(E_{1},...E_{n})$, EI is
defined as,

$a_{EI}(\theta)=\mathbb{E}_{\theta}[\max(0,E_{min}-E(\theta))]$ (18)

Where $\mathbb{E}_{\theta}[\cdot]$ is evaluated over all the possible
$E(\theta)$ from the surrogate model.

Expected Improvement can be used in its explicit closed form which is,

$\displaystyle a_{EI}(\theta)=$
$\displaystyle\left(E_{min}-E(\theta)\right)\Phi\left(\frac{E_{min}-E(\theta)}{\Delta
E(\theta)}\right)$ (19) $\displaystyle+\Delta
E(\theta)\phi\left(\frac{E_{min}-E(\theta)}{\Delta E(\theta)}\right)$

Where $\Phi$ is the Normal distribution cumulative distribution function and
$\phi$ is the corresponding probability density function. We find the minimum
of $a_{EI}$ using Powell’s algorithm.

Bayesian optimization uses very few shots when compared to the other methods,
as it takes into account the noise on the measurement. It is, however,
extremely computationally expensive.

###  5.1 Mapping QUBO to the VQE

We can use the VQE to solve QUBO problems using the transformation

$x^{i}\mapsto\frac{1}{2}\left(\sigma^{i}_{0}+\sigma_{3}^{i}\right),$ (20)

where $\sigma_{3}^{i}$ denotes the Pauli-Z spin matrix acting on the $i^{th}$
qubit. We then find the groundstate eigenvalue of the Hamiltonian

$H=\sum_{i,j=1}^{q}Q_{ij}\frac{1}{4}\left((\sigma^{i}_{0}+\sigma_{3}^{i})(\sigma^{j}_{0}+\sigma_{3}^{j})\right).$
(21)

We then take the parameters which gave us this groundstate, measure the
produced state in the $\sigma_{3}^{\otimes q}$ basis, and the state with the
highest likelihood of being measured which meets our constraints is our
solution.

Due to the fact that Hamiltonians only consist of $\sigma_{0}$ and
$\sigma_{3}$ terms, they are diagonal, hence we only need one parameter per
qubit in order to map a sufficient amount of the space to receive the optimum.
However, to aid the classical optimizer in the noisy space of quantum
measurements, $q$ layers of rotation and entangling $CNOT$ gates will be used.

We note here that due to the density of the weight matrix $Q$ defined in Eq.
10, the Hamiltonian is very dense, often containing $q^{2}$ nonzero terms.
This means that this is a worst case in terms of stress on the VQE.

###  5.2 Dimensionallity Expressivity Analysis

When using variational quantum circuits, something that we often look at is
Dimensionallity Expressivity Analysis (DEA) [49]. This is a method by which we
can determine if the parameters of the gates are independent or redundant. The
methodology behind this can be found in Ref. [50].

DEA is more applicable when we are looking at more complex systems, i.e. when
we are not only searching for basis states and thus need to be able to span a
subset of the whole space.

If we do carry out this analysis on our circuit, we see that all the
parameters are non-redundant, as long as we do not start the optimization
process with any parameter values equal to $(2k+1)/4\cdot\pi$ and
$(2k)/4\cdot\pi$ for $k=\\{0,1\\}$. In order to ensure this, the parameter
values are uniformly randomly generated in the region $[0,2\pi)$, using
numpy.random.uniform.

Further details about DEA are included in Appen. A.

###  5.3 Conditional Value at Risk

A recent development in the usability of the VQE to solve combinatorial
optimization problems is to use the Conditional Value at Risk (CVaR) technique
from finance. Full details of CVaR-VQE and its advantages can be found in Ref.
[51]. In this reference, CVaR was used with regard to the sample mean, but in
this work, we apply it to the expected value.

For a Hamiltonian containing $K$ terms, we define the energy measurements to
be $H_{\gamma}$. In this notation, the coefficients and normalisation with
regards to the shot count have been included in $H_{\gamma}$. Our usual
objective function is,

$\text{VQE}\rightarrow\sum_{\gamma=1}^{K}H_{\gamma}$ (22)

If we now order the measurements so that $H_{k}\leq H_{k+1}$, the CVaR
adaptation is to then introduce a new parameter: $0<\alpha\leq 1$ such that,

$\text{CVaR-VQE}\rightarrow\sum_{\gamma=1}^{\lceil\alpha K\rceil}H_{\gamma}$
(23)

In this way, we truncate our measurements to only include the terms which
greatest contribute to our minimization. By changing this value we can more
quickly find the state which gives us our minimum.

Note that the minimum value of CVaR-VQE is not necessarily equal to the
minimum of the VQE, but the correct solution can be the same.

##  6 Design of the tests of the algorithms

We have tested the power produced by the different methods for grids of size
$l_{grid}=4$. This was the maximum system size that we could run on quantum
simulators using the resources available to us. By exhaustive search of all
the possibilities we found the maximum power for the grid with $l_{grid}=4$.
The WFLO problem has many local optima, so we run each algorithm 36 times from
different initial conditions. To check the scaling of the performance of the
different algorithms, we run the quantum simulators with $l_{grid}=2$ to
$l_{grid}=4$, and the classical optimizers with $l_{grid}=2$ to $l_{grid}=10$.

Justification for 36 being a sufficient sample size can be found in Appen. B.

###  6.1 Test model

In order to test these different algorithms, we will use the following
parameter values (the Mosetti benchmark cases are from Ref. [6]),

$\displaystyle D$ $\displaystyle=\text{Mosetti second benchmark case}$
$\displaystyle x$ $\displaystyle=1,\text{ }r=1.5,\text{ }m=4,\text{ }\xi=0$

where $x$ is the maximum distance that is affected by the wake and $r$ is the
radius with which it spreads out per unit distance away from the turbine
causing the wake. $m$ is the maximum number of turbines allowed, and $\xi$ is
the minimum distance allowed between turbines, in this work, we have set this
to zero so that turbines may be on adjacent sites. Due to our limited ability
to simulate large quantum systems, we will be looking at a very small grid
with $l_{grid}=4$. This means that there will be 16 possible positions for the
turbines, and so will need 16 qubits.

Mosetti’s second benchmark case is defined as:

$\displaystyle D$ $\displaystyle=\\{\\{10k,12ms^{-1},1/36\\}\\}_{k=0}^{36}$
(24)
$\displaystyle=\\{\\{0,12ms^{-1},0.028\\},...,\\{350,12ms^{-1},0.028\\}\\}$
(25)

This can be visualized as the distribution shown in Fig. 7. Here, each box
represents a windspeed of $12ms^{-1}$.

###  6.2 Computational details of simulations

The quantum simulations were run using IBM’s qiskit:0.39.2 package for python
[15].

Both PO and COByLA were utilized via scipy.optimize.minimize \- using the
methods ‘powell’ and ‘COBYLA’ respectively as part of the scipy library [27].
BO was coded up by for this work, using a heavily modified version of the code
used in [52], where we no longer make use of
sklearn.gaussian_process.GaussianProcess-
Regressor as this has problems when using the periodic kernel - and so utilize
the matrix form in order to carry out the calculations.

SA was implemented using the package pyqubo [53, 54]. This library has an
inbuilt annealer by using
neal.SimulatedAnnealingSampler() once the system has been turned into a QUBO
problem via model.to_qubo(). Pyqubo has a wrapper to the D-Wave Ocean SDK, so
we also ready to run on the D-wave systems.

Gurobi was implemented using the gurobipy package [28]. Gurobi was used with
an academic license through the University of Plymouth.

As well as the previously mentioned methods, we also investigated qiskit’s
built-in VQE QUBO solver. This uses their SamplingVQE function, and then
performs a MinimumEigenOptimizer over it in order to find the groundstate. The
ansatz used is the same as the one described in Sec. 5.1. This method acts as
a control, as it is the VQE method with no statistical noise present. We use
COByLA as the classical optimization routine for this method. Details about
this method can be found in Ref. [55].

All of the simulations for the results provided below were run using the
computational facilities of the High Performance Computing Centre, University
of Plymouth [56]. The CPUs were dual node Intel Xeon E5-2683v4 CPU each with
16 cores. Task farming was used to run the simulations using all the cores.

##  7 Results

###  7.1 How solutions are selected

For both the Gurobi solver and quantum annealing, the selection of the
solution is trivial as it will be the final binary string on which the
algorithm lands.

However, for the VQE this is less obvious. Once the optimization process has
finished, we are given a set of parameters which corresponds to the quantum
circuit gates. We then input these values back into the circuit, and measure
in the $\sigma_{3}$ basis on all qubits; this produces a distribution over all
possible basis states. From this distribution, we then take the state which
has the highest probability and only has four $|1\rangle$ terms. We only
accepted solutions which have exactly the correct number of turbines, which is
four for our test case.

###  7.2 Discussion about the degeneracy of solutions

As for many other QUBO problems, there are many degenerate optimal solutions.
For the case $l_{grid}=4$, there are 79 unique optimal solutions out of 1820
binary vectors which meet the constraints and 65536 total possible binary
strings of length 16. A list of all optimal solutions can be found in the
Appen. C. Each of these optimal solutions gives us a power output of 2304.0
kW.

Some of these degeneracies are due to the fact that there are optimal
solutions for the conditions in smaller system sizes. The ‘smaller’ solutions
can then be embedded in the larger systems. There exist no optimal solutions
for a system size $l_{grid}\leq 2$.

Figure 7: $2^{nd}$ Mosetti benchmark case wind regime distribution. Each bar
represents a windspeed of 12 $ms^{-1}$. Figure 8: Heatmap of $\langle
X\rangle$ values for the degenerate solutions with maximum power. All the
values sum to 4 (the number of turbines)

Another way in which we could study a structure in the optimal solutions is to
look at their average placement locations. If we label a solution matrix to be
$X_{*}^{i}$, this is a binary matrix which is representative of the solution
grid. We can map from solutions by using the labelling as in Fig. 2. We can
then calculate the average placement as:

$\langle X\rangle=\frac{1}{\text{\\#solutions}}\sum_{i}X_{*}^{i}$ (26)

The results from this can be seen in Fig. 8, note here that the sum of these
should equal four as this is the maximum number of turbines, but the numbers
have been limited to three d.p. for ease of viewing. From this, we see that
when we sum over the optimal solutions which satisfy the constraints, that on
average they want to be near the outsides, with even more emphasis on being in
the corners. This makes sense physically, with the wind regime used because
this maximizes the distance between any two turbines and thus minimizes any
possible wake effects.

###  7.3 Discussion of the results

We found that the Powell optimization method in the simulations produced
results in a reasonable time for $l_{grid}=2$ and $l_{grid}=3$ cases. However,
the Powell optimizer was unable to produce results in a reasonable amount of
time for $l_{grid}=4$. This was due to that it requires a large shot count,
and also requires many function evaluations in order to be effective. We only
present results using the two surrogate model based methods, COByLA and BO,
and the two classical algorithms.

Fig. 9 reports the results for the power from the different algorithms
investigated for the $l_{grid}\leq 4$ case. Each of the methods was run 36
times and the boxplots show the spread of these results. The $x$-axis shows
the names of the different methods, the number under the name of the VQE
method is the CVaR $\alpha$ value. ‘Gurobi SS’ is the Gurobi method where each
run was started with the same default guess that the algorithm chooses,
‘Gurobi RS’ is the method with random starting guesses.

Figure 9: Each value along the $x$-axis represents a the different methods.
The box plots show the power obtained from running the problems 36 times. The
red dashed line is the optimal solution power output, found through an
exhaustive list search. The Gurobi method may look empty, but there is a line
at the optimal - it was able to achieve this every run. The circle points
which are not within the boxes are outlier points.

From Fig. 9 we can draw a few conclusions. The Gurobi method with the default
starting guess outperforms all the other methods, consistently reaching the
optimal solution every time. This is also true for the noiseless COBYLA-VQE
method from qiskit. We also see that the different VQE methods provide similar
quality of results to one another, which are also comparable to SA. The Gurobi
method using random restarts performs the worst out of all the methods.

We see that each of the COBYLA-CVaR methods is able to reach the optimal
solution; we also see that using $\alpha<1$ has little effect on the quality
of results. We see that the Bayesian-CVaR has very comparable results.

Tab. 1 shows the average power output of each of the methods as a percentage
of the optimal power output. We see that each of these methods does perform
very well on average, with the lowest average output being $93.3\%$ of the
maximum.

Method |  | Average output result
---
(% of optimal)
COByLA-1.00 | 95.5
COByLA-0.75 | 94.7
COByLA-0.50 | 94.5
COByLA-0.25 | 95.2
Bayes-1.00 | 93.3
Bayes-0.75 | 94.2
Bayes-0.50 | 94.6
Bayes-0.25 | 93.4
Table 1: Average output power from solutions of different VQE-based methods as
a percentage of the optimal power output.

###  7.4 Timing the results

As well as looking at the quality of results, we are also interested in how
long the algorithms take to get their results.

In order to test this, we have used Python’s time package, to take a total
time measurement from start to end of the whole process. The results from this
can be shown in Tab. 2

$l_{grid}$ | Gurobi [s] | Simulated Annealing [s] | COByLA-CVaR [s] | BO-CVaR [s]
---|---|---|---|---
2 | - | - | 19.74 | 250.11
3 | 0.09 | 0.04 | 561.78 | 5741.43
4 | 1.11 | 0.85 | 11546.10 | 55454.77
5 | 3.91 | 3.31 | - | -
6 | 14.51 | 11.49 | - | -
7 | 107.11 | 106.62 | - | -
8 | 374.72 | 377.49 | - | -
9 | 983.27 | 985.57 | - | -
10 | 2240.78 | 2220.44 | - | -
Table 2: Average time in seconds taken for the different algorithms

From these results, we can easily see that the simulated quantum system takes
far more time than the other two algorithms. However, what is of more interest
to us is the scalability of these results, i.e. how much slower they get as we
increase the system size. We can study this by looking at the coefficients of
LogLog plots of the timings.

Figure 10: LogLog plots of the time taken for the different methods over different system sizes. Method | Gradient | Intercept
---|---|---
Gurobi | 4.24 | -5.18
Simulated Annealing | 4.51 | -5.67
CVaR-VQE using COByLA | 4.57 | -1.50
CVaR-VQE using BO | 3.89 | 0.04
Table 3: Coefficients of linear fits for LogLog timings

The time to solve the QUBO problem is expected to depend on the system size
exponentially. However, as we only have three data points, we parameterize the
time to solution as a polynomial with coefficient $\alpha$.

$\mbox{Time to solve}\propto(V)^{\alpha}$

Where $V=l_{grid}^{2}$ is the system volume. Fig. 10 shows the log of time
taken to solution versus the log of the system volume. The timings are in
Table 3. The results from all four algorithms in Fig. 10 have similar slopes.

We find that: Gurobi $\sim O(V^{4.24})$, SA $\sim O(V^{4.51})$, COByLA $\sim
O(V^{4.57})$ and Bayesian $\sim O(V^{3.89})$. These values are calculated by
linearly fitting the LogLog data from timing each algorithm twelve times and
taking the mean. This is interesting, as while Gurobi and SA are the true
algorithms, the quantum method is a simulation, and thus we would expect that
the real quantum circuit would reduce this time.

Very interestingly, using Bayesian optimization even on the simulator scales
better than all the other methods, having the lowest polynomial scaling; this
does come with the caveat that the pre-factor is much larger than the other
systems. We can look at this more closely, and by searching for the
intersection point of the Bayesian and Gurobi timing lines, this is shown in
Fig. 11. We can estimate that the simulated Bayesian-CVaR method would be
faster than the Gurobi method for a system volume of $\approx 10^{15}$, which
is a square of sides length $\approx 31622776$.

Figure 11: Extrapolation of time taken for Gurobi and Bayesian Optimization in
order to find the point at which Bayesian would be faster.

###  7.5 Conclusions

The VQE-based solver of the WFLO problem finds good solutions when a
sufficient number of measurements have been done. The Gurobi optimizer always
finds the optimal solution and it outperforms the VQE-based method. This is
not surprising as Gurobi has been developed to specifically solve linear and
quadratic programming problems, of which QUBO is a special case. We have only
investigated small system sizes and it will be interesting to compare the
results from the VQE algorithm running on a quantum computer. Many
benchmarking studies that compare different algorithms to solve the WFLO
problem use a grid with 100 points, which would require 100 qubits in our
approach. Solving the WFLO on a quantum computer would require many more than
100 qubits when error correction is used.

It is interesting to see that the scaling of the time taken with system size
had similar exponents between the Gurobi and the VQE-based method running on
the quantum simulator. The VQE method uses a mixture of classical and quantum
computing. From our results, the performance of the VQE method critically
depends on the performance of the optimizer that runs on the classical
computer and this will have to be tuned for good performance.

For WFLO problems, it has been found to be useful to use multi-objective
optimisation [7] to, for example, maximize the power and minimize the cost.
The QUBO formalism has recently been extended to include multiple objectives
[57]. This would be an interesting area to explore in the future.

We are encouraged by the results we have obtained on the quantum computing
simulator to study the algorithms on real quantum computers. We will next
study the resilience of the implemented methods against quantum errors on the
simulators. It will be interesting to compare the performance on circuit-based
quantum computers to those from adiabatic quantum systems.

###  7.6 Acknowledgements

We thank Julian Stander for discussions about Bayesian optimisation. This work
was carried out using the computational facilities of the High Performance
Computing Centre, University of Plymouth – https://www.plymouth.ac.uk/about-
us/university-structure/faculties/science-engineering/hpc.

## References

  * [1] P. Veers et al.,  Science 366, eaau2027 (2019).
  * [2] A. Giani and Z. Eldredge,  SN Computer Science 2, 393 (2021).
  * [3] C. Berger et al.,  arXiv preprint arXiv:2107.05362 (2021).
  * [4] A. Ajagekar and F. You,  Renewable and Sustainable Energy Reviews 165, 112493 (2022).
  * [5] A. Abbas et al.,  arXiv preprint arXiv:2312.02279 (2023).
  * [6] G. Mosetti, C. Poloni, and B. Diviacco,  Journal of Wind Engineering and Industrial Aerodynamics 51, 105 (1994).
  * [7] P. L. Manikowski, D. J. Walker, and M. J. Craven,  Journal of Marine Science and Engineering 9, 1376 (2021).
  * [8] Vu Phong Energy Group JSC,  The world’s biggest wind farms,  https://vuphong.com/the-worlds-biggest-wind-farms/, 2022 (accessed 01/11/2023).
  * [9] A. Senderovich, J. Zhang, E. Cohen, and J. C. Beck,  IEEE Access 10, 78044 (2022).
  * [10] J. Gambetta,  Quantum-centric supercomputing: The next wave of computing,  https://research.ibm.com/blog/next-wave-quantum-centric-supercomputing,  Accessed: 27-06-2023.
  * [11] A. Rajak, S. Suzuki, A. Dutta, and B. K. Chakrabarti,  Philosophical Transactions of the Royal Society A 381, 20210417 (2023).
  * [12] C. McGeoch and P. Farré,  D-Wave systems (2022).
  * [13] C. McGeoch, P. Farré, and K. Boothby,  D-Wave systems (2022).
  * [14] S. Matsubara et al.,  Digital annealer for high-speed solving of combinatorial optimization problems and its applications,  in 2020 25th Asia and South Pacific Design Automation Conference (ASP-DAC), pp. 667–672, IEEE, 2020.
  * [15] Qiskit contributors,  Qiskit: An open-source framework for quantum computing, 2023.
  * [16] J. Preskill,  Quantum 2, 79 (2018).
  * [17] E. Zahedinejad and A. Zaribafiyan,  arXiv preprint arXiv:1708.05294 (2017).
  * [18] T. Albash and D. A. Lidar,  Physical Review X 8, 031016 (2018).
  * [19] A. Pearson, A. Mishra, I. Hen, and D. A. Lidar,  npj Quantum Information 5, 107 (2019).
  * [20] R. Yaacoby et al.,  Physical Review E 105, 035305 (2022).
  * [21] V. Martin-Mayor and I. Hen,  Scientific reports 5, 15324 (2015).
  * [22] P. Y. Zhang, D. A. Romero, J. C. Beck, and C. H. Amon,  Solving wind farm layout optimization with mixed integer programming and constraint programming,  in Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems, pp. 284–299, Berlin, Heidelberg, 2013, Springer.
  * [23] S. Donovan,  Proceedings of the Annual ORSNZ Conference 40 (2005).
  * [24] A. Crippa et al.,  (2023), arXiv:2304.01690.
  * [25] G. Kochenberger et al.,  Journal of combinatorial optimization 28, 58 (2014).
  * [26] M. Khanali, S. Ahmadzadegan, M. Omid, F. Keyhani Nasab, and K. W. Chau,  International Journal of Energy and Environmental Engineering 9, 399 (2018).
  * [27] T. do Cuoto, B. Farias, A. Carlos G.C. Diniz, and M. Vinicius G. de Morais,  10th World Congress on Structural and Multidisciplinary Optimization (2013).
  * [28] Gurobi Optimization LLC,  Gurobi optimizer reference manual,  https://www.gurobi.com, 2023,  Accessed: 27-06-2023.
  * [29] O. Şeker, N. Tanoumand, and M. Bodur,  Applied Soft Computing 127, 109367 (2022).
  * [30] K. Nakano et al.,  arXiv preprint arXiv:2207.03069 (2022).
  * [31] S. Kirkpatrick, C. D. Gelatt Jr, and M. P. Vecchi,  science 220, 671 (1983).
  * [32] M. Pincus,  Operations Research 18, 1225 (1970).
  * [33] Z. Wang, X. Geng, and Z. Shao,  Journal of Computational and Theoretical Nanoscience 6, 1680 (2009).
  * [34] E. Rodriguez-Tello, J.-K. Hao, and J. Torres-Jimenez,  Computers & Operations Research 35, 3331 (2008).
  * [35] K. Tole, R. Moqa, J. Zheng, and K. He,  Computers & Industrial Engineering 176, 109004 (2023).
  * [36] A. Peruzzo et al.,  Nature Communications 5 (2014).
  * [37] J. Tilly et al.,  Physics Reports 986, 1 (2022).
  * [38] X. Liu et al.,  IEEE Transactions on Quantum Engineering 3, 1 (2022).
  * [39] T. Schwägerl et al.,  Frontiers of Physics 18, 21308 (2023).
  * [40] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik,  New J. Phys. 18, 023023 (2016), arXiv:1509.04279.
  * [41] N. Klco et al.,  Phys. Rev. A 98, 032331 (2018), arXiv:1803.03326.
  * [42] R. Santagati et al.,  Science advances 4, eaap9646 (2018).
  * [43] G. Iannelli and K. Jansen,  arXiv preprint arXiv:2112.00426 (2021).
  * [44] J. S. Otterbach et al.,  arXiv preprint arXiv:1712.05771 (2017).
  * [45] D. Zhu et al.,  Science advances 5, eaaw9918 (2019).
  * [46] J. Snoek, H. Larochelle, and R. P. Adams,  Advances in neural information processing systems 25 (2012).
  * [47] V. Nguyen,  Bayesian optimization for accelerating hyper-parameter tuning,  in 2019 IEEE second international conference on artificial intelligence and knowledge engineering (AIKE), pp. 302–305, IEEE, 2019.
  * [48] J. Bergstra, R. Bardenet, Y. Bengio, and B. Kégl,  Advances in neural information processing systems 24 (2011).
  * [49] L. Funcke, T. Hartung, K. Jansen, S. Kühn, and P. Stornati,  Quantum 5, 422 (2021).
  * [50] L. Funcke et al.,  Dimensional expressivity analysis, best-approximation errors, and automated design of parametric quantum circuits, 2021, arXiv:2111.11489.
  * [51] P. K. Barkoutsos, G. Nannicini, A. Robert, I. Tavernelli, and S. Woerner,  Quantum 4, 256 (2020).
  * [52] J. Brownlee,  How to implement bayesian optimization from scratch in python,  https://machinelearningmastery.com/what-is-bayesian-optimization/,  Accessed: 27-06-2023.
  * [53] M. Zaman, K. Tanahashi, and S. Tanaka,  IEEE Transactions on Computers (2021).
  * [54] K. Tanahashi, S. Takayanagi, T. Motohashi, and S. Tanaka,  Journal of the Physical Society of Japan 88, 061010 (2019).
  * [55] Q. O. D. Team,  Improving variational quantum optimization using cvar,  https://qiskit.org/ecosystem/optimization/tutorials/08_cvar_optimization.html,  Accessed: 27-11-2023.
  * [56] HPC team UoP,  High performance computing (hpc),  https://www.plymouth.ac.uk/about-us/university-structure/faculties/science-engineering/hpc,  Accessed: 06-11-2023.
  * [57] M. Ayodele, R. Allmendinger, M. López-Ibáñez, and M. Parizy,  Multi-objective qubo solver: Bi-objective quadratic assignment problem,  in Proceedings of the Genetic and Evolutionary Computation Conference, pp. 467–475, 2022.

## Appendix

##  Appendix A Dimensionallity Expressivity Analysis

A way in which we can analyse the system and try to improve our results is to
carry out expressivity analysis on the circuit matrix. A parameter
$\theta_{k}$ is considered ‘redundant’ if a small change of $\theta_{k}$,
keeping all others constant produces a state which can also be obtained by
keeping $\theta_{k}$ constant and adjusting all other $\theta_{j}$
accordingly. In order to do this analysis, we are using the method that is
described by Lena Funcke, et al in Ref. [50]. The method is as follows:

  1. 1.

This can be checked by considering the real partial Jacobians $J_{k}$ of $C$
($C$ is the matrix of the circuit applied to the state $|0\rangle^{\otimes
q})$:

$J_{k}(\theta)=\begin{pmatrix}Re(\partial_{1}C)...Re(\partial_{k}C)\\\
Im(\partial_{1}C)...Im(\partial_{k}C)\end{pmatrix}$

Here it is key to understand that $\partial_{k}C$ is itself a $2^{q}\times 1$
vector, and thus $J_{k}$ is a $2^{q+1}\times k$ matrix, where we start with
$k=1$ for the first parameter.

  2. 2.

We then check the rank of the matrix $J_{k}$ for each $k$, and if adding a new
parameter does not increase the rank, the parameter must be redundant.

In order to efficiently check the rank of $J_{k}$, we consider the matrix
$S_{k}=J_{k}^{*}J_{k}$. Thus if we check that $S_{k}$ is invertible we know
that all the parameters are independent (i.e., not redundant).

##  Appendix B Justification for thirty-six samples

When gathering results for this work, we chose 36 samples for each of the
methods. In this section we will look at one particular method, being COBYLA
with the CVaR $\alpha=0.25$. We take a much larger sample (284) and compare
the averages. By doing this we can see that 36 samples shows similar features
and is thus sufficient to draw conclusions from.

The initial 36 samples taken were:

        2136.54, 2118.96, 2202.69, 2220.27, 2220.27, 2220.27, 2268.84,
        2220.27, 2286.42, 2220.27, 2136.54, 2118.96, 2118.96, 2220.27,
        2101.38, 2202.69, 2220.27, 2286.42, 2202.69, 2286.42, 2304.0,
        2118.96, 2286.42, 2118.96, 2118.96, 2136.54, 2136.54, 2268.84,
        2118.96, 2118.96, 2101.38, 2286.42, 2220.27, 2286.42, 2202.69,
        2118.96

The extra 284 samples taken were:

        2251.25, 2202.69, 2185.11, 2185.11, 2136.54, 2136.54, 2286.42,
        2118.96, 2268.84, 2304.0, 2202.69, 2136.54, 2220.27, 2268.84,
        2101.38, 2118.96, 2268.84, 2220.27, 2220.27, 2202.69, 2202.69,
        1933.92, 2118.96, 2202.69, 2268.84, 2220.27, 2220.27, 2220.27,
        2268.84, 2220.27, 2185.11, 2136.54, 2118.96, 2202.69, 2251.25,
        2304.0, 2052.82, 2118.96, 2286.42, 2136.54, 2118.96, 2220.27,
        2118.96, 2118.96, 2268.84, 2304.0, 2101.38, 2268.84, 2286.42,
        2035.23, 2304.0, 2202.69, 2118.96, 2304.0, 2202.69, 2017.65,
        2268.84, 2220.27, 2136.54, 2304.0, 2202.69, 2118.96, 2268.84,
        2136.54, 2136.54, 2202.69, 2118.96, 2035.23, 2202.69, 2304.0,
        2118.96, 2220.27, 2136.54, 2185.11, 2118.96, 2220.27, 2251.25,
        2220.27, 2304.0, 2202.69, 2202.69, 2220.27, 2268.84, 2304.0,
        2185.11, 2136.54, 2202.69, 2118.96, 2185.11, 2220.27, 2118.96,
        2118.96, 2268.84, 2220.27, 2185.11, 2136.54, 2251.25, 2304.0,
        2268.84, 2202.69, 2220.27, 2185.11, 2202.69, 2118.96, 2220.27,
        2220.27, 2220.27, 2286.42, 2268.84, 2220.27, 2185.11, 2202.69,
        2101.38, 2202.69, 2220.27, 2202.69, 2304.0, 2304.0, 2220.27,
        2202.69, 2220.27, 2286.42, 2220.27, 2118.96, 2118.96, 2220.27,
        2202.69, 2035.23, 2118.96, 2202.69, 2286.42, 2118.96, 2101.38,
        2118.96, 2202.69, 2304.0, 2136.54, 2304.0, 2185.11, 2202.69,
        2286.42, 2101.38, 2202.69, 2220.27, 2220.27, 2220.27, 2220.27,
        2220.27, 2136.54, 2220.27, 2101.38, 2220.27, 2286.42, 2118.96,
        2220.27, 2220.27, 2286.42, 2220.27, 2136.54, 2136.54, 2268.84,
        2101.38, 2202.69, 2101.38, 2035.23, 2220.27, 2136.54, 2202.69,
        2136.54, 2286.42, 2185.11, 2286.42, 2220.27, 2118.96, 2118.96,
        2220.27, 2202.69, 2202.69, 2118.96, 2118.96, 2118.96, 2202.69,
        2220.27, 2202.69, 2220.27, 2202.69, 2286.42, 2118.96, 2202.69,
        2136.54, 2220.27, 2268.84, 2118.96, 2251.25, 2286.42, 2136.54,
        2136.54, 2286.42, 2286.42, 2220.27, 2202.69, 2286.42, 2136.54,
        2220.27, 2202.69, 2220.27, 2118.96, 2202.69, 2304.0, 2220.27,
        2202.69, 2304.0, 2286.42, 2304.0, 2286.42, 2101.38, 2220.27,
        2101.38, 2220.27, 2202.69, 2202.69, 2220.27, 2035.23, 2304.0,
        2304.0, 2202.69, 2118.96, 2185.11, 2202.69, 2220.27, 2118.96,
        2101.38, 2202.69, 2268.84, 2304.0, 2202.69, 2017.65, 2118.96,
        2286.42, 1933.92, 2304.0, 2118.96, 2202.69, 2286.42, 2220.27,
        2136.54, 2136.54, 2185.11, 2017.65, 2185.11, 2220.27, 2202.69,
        2118.96, 2136.54, 2118.96, 2220.27, 2017.65, 2251.25, 2136.54,
        2202.69, 2035.23, 1933.92, 2136.54, 2017.65, 2220.27, 2118.96,
        2101.38, 2220.27, 2286.42, 2118.96, 2286.42, 2220.27, 2286.42,
        2220.27, 2286.42, 2202.69, 2268.84, 2101.38, 2251.25, 2220.27,
        2035.23, 2286.42, 2202.69, 2202.69

The sample mean of the 36 initial samples is: 2193.1, and the sample mean of
the extra is: 2192.4. We can see from these that these are similar values
($0.03\%$ difference), and so we are able to draw conclusions about the
method’s effectiveness on average by only taking 36 samples. This is very
useful as minimizing time and resource waste is an important factor when
carrying out simulations.

##  Appendix C Degeneracy of solutions

Below is a list of all possible optimal solutions for the problem defined in
Sec. 6.1.

     (0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1),
     (0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1),
     (0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0),
     (0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1),
     (0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1),
     (0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0),
     (0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1),
     (0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1),
     (0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0),
     (0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1),
     (0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1),
     (0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0),
     (0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0),
     (0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0),
     (0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1),
     (0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1),
     (0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0),
     (0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0),
     (0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0),
     (0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0),
     (0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1),
     (0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1),
     (0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0),
     (0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0),
     (0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0),
     (0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0),
     (0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1),
     (0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1),
     (0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0),
     (0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1),
     (0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1),
     (0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0),
     (0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1),
     (0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1),
     (0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0),
     (0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0),
     (0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0),
     (0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0),
     (0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1),
     (0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0),
     (0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1),
     (0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0),
     (0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0),
     (0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0),
     (1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1),
     (1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1),
     (1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0),
     (1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1),
     (1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1),
     (1, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 0),
     (1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1),
     (1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1),
     (1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0),
     (1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1),
     (1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0),
     (1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1),
     (1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1),
     (1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0),
     (1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0),
     (1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0),
     (1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0),
     (1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1),
     (1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0),
     (1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1),
     (1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0),
     (1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0),
     (1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0),
     (1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1),
     (1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1),
     (1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0),
     (1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0),
     (1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0),
     (1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0),
     (1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1),
     (1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0),
     (1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1),
     (1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0),
     (1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0),
     (1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0)


# Boosting neural video codecs by exploiting hierarchical redundancy

Reza Pourreza, Hoang Le, Amir Said, Guillaume Sautière, Auke Wiggers
Qualcomm AI Research333 Qualcomm AI Research is an initiative of Qualcomm
Technologies, Inc.
{pourreza, hoanle, asaid, gsautie<EMAIL_ADDRESS>
###### Abstract

In video compression, coding efficiency is improved by reusing pixels from
previously decoded frames via motion and residual compensation. We define two
levels of hierarchical redundancy in video frames: 1) first-order: redundancy
in pixel space, _i.e_., similarities in pixel values across neighboring
frames, which is effectively captured using motion and residual compensation,
2) second-order: redundancy in motion and residual maps due to smooth motion
in natural videos. While most of the existing neural video coding literature
addresses first-order redundancy, we tackle the problem of capturing second-
order redundancy in neural video codecs via predictors. We introduce generic
motion and residual predictors that learn to extrapolate from previously
decoded data. These predictors are lightweight, and can be employed with most
neural video codecs in order to improve their rate-distortion performance.
Moreover, while RGB is the dominant colorspace in neural video coding
literature, we introduce general modifications for neural video codecs to
embrace the YUV420 colorspace and report YUV420 results. Our experiments show
that using our predictors with a well-known neural video codec leads to 38%
and 34% bitrate savings in RGB and YUV420 colorspaces measured on the UVG
dataset.

##  1 Introduction

Video traffic is predicted to reach more than 82% of all internet traffic in
2022 [6]. As video resolution and framerate demands have steadily increased
over the years, there is a strong need to develop algorithms that enable
transmission of video data at low bit cost. Standard video codecs such as AVC
[34] and HEVC [27] were created exactly for this purpose.

The goal of a codec is to trade off the number of bits spent to transmit a
message and the distortion between the original data and its reconstruction.
In the video setting, standard codecs keep the bitrate low by using previously
transmitted information as context. For example, to make a prediction for a
current frame, the previously transmitted frame can be used. Instead of
transmitting the current frame, the motion between the previous reconstruction
and the current frame may be transmitted. This allows making an initial
prediction for the current frame using motion compensation. Transmitting the
motion typically has a lower bit cost than transmitting the frame, as motion
of pixels belonging to the same object are strongly correlated. Similarly, an
additive correction or _residual_ can be transmitted to refine a receiver-side
frame. If the residual is sparse, the cost of transmitting a residual will be
much lower than transmitting the frame in isolation.

Standard codecs take this idea one step further and make _motion vector
predictions_. As motion vectors are typically similar across timesteps due to
smooth object motion in video, an initial prediction can be made based on
previous motion. This means the codec only needs to transmit corrections to
this prediction, usually at low bit cost.

Figure 1:  First- and second-order redundancy in video, here shown for
(backward) optical flow. First-order redundancy exists between frames in
pixel-space, second-order redundancy exists between optical flow maps.

To distinguish these two approaches, we distinguish between two types of
redundancy in videos. First-order redundancy is redundancy in the pixel
domain. Inter-frame codecs exploit this using motion compensation and residual
coding. Second-order redundancy is the redundancy in flow and residual maps
due to smooth motions in videos, which can be exploited by making motion
vector predictions. These are shown schematically for motion in Fig. 1.

Most modern neural video codecs exploit only first-order redundancy using a
two-step process [17, 2, 13, 11]: a _flow model_ transmits motion information
(optical flow) and _residual model_ transmits residual maps. There are works
that exploit second-order redundancy as well. For example, ELF-VC [23]
predicts the flow, and Liu _et al._[16] use flow and residual prediction in
latent space.

In this work, we exploit second-order redundancy via pixel-space flow and
residual predictors. We design generic models that predict flow and residual
using the history available in a frame buffer. This leads to easy-to-interpret
predictions compared to predictions in latent space, makes it possible to
combine the predictors with many existing neural codecs, and delivers
significant performance improvement.

Furthermore, we provide evaluations for multiple input colorspaces. Most
literature on neural video codecs only provides solutions and evaluations for
the RGB colorspace, with some exceptions [7, 24]. However, standard codecs
were designed to operate with YUV420 input. YUV420 is a color format that
subsamples the chroma channels of the YUV colorspace which is closer to how
humans perceive images and leads to bandwidth saving in compression
applications. To gauge progress in learned video compression compared to
standard codecs, we train our codecs on both RGB and YUV420 input formats. We
demonstrate that neural codecs can excel in both domains, and that only minor
architectural modifications are needed to handle the difference in input
resolution between colorspaces.

In summary, the contributions of this paper are:

  1. 1.

We define the framework of hierarchical redundancy in video in order to
categorize neural codecs, and observe that most neural codecs only address the
first type of redundancy.

  2. 2.

To exploit second-order redundancy, we introduce generic flow and residual
predictors that can be paired with existing inter-frame codecs and deliver
strong performance improvements.

  3. 3.

We show that simple modifications enable codecs to run on YUV420 inputs, and
report both RGB and YUV420 performance.

##  2 Related work

###  2.1 Learned compression

The focus of this work is a so-called _learned codec_ , meaning each component
is a learned model. Neural network-based codecs have been successfully applied
to data compression in many domains, including the image [30, 28, 3] and video
[17, 9, 2, 13, 23, 11] domain. In these systems, a neural encoder takes the
data $\mathbf{x}$ and produces a quantized latent variable $\mathbf{z}$, and a
neural decoder produces a reconstruction $\widehat{\mathbf{x}}$ given this
latent. A neural prior or context model is used to learn the distribution of
latent variables $p(z)$. Using this prior and an entropy coding algorithm,
Shannons source-coding theory tells us that the latent $\mathbf{z}$ can be
losslessly compressed using $-\log p(\mathbf{z})$ bits.

Neural codecs are typically trained to minimize a rate-distortion loss
consisting of two terms:

$\mathcal{L}_{RD}=\mathbb{E}_{\mathbf{x}\sim
p(x)}\left[\beta\mathcal{L}_{rate}(\mathbf{z})+\mathcal{L}_{distortion}(\mathbf{x},\widehat{\mathbf{x}})\right].$
(1)

The rate term corresponds to the number of bits needed to transmit the
quantized latent variable $\mathbf{z}$, and a distortion term corresponds to
the distance between the reconstruction $\widehat{\mathbf{x}}$ and the ground
truth $\mathbf{x}$.

Neural codecs promise several advantages over handcrafted codecs. First, as
they learn to identify redundancies from example data, they have shown a
strong ability to specialize to a domain [9] or even a single datapoint [32,
25, 31]. For example, if the codec will only be used to code animated content,
it can easily be finetuned to this domain. Second, neural codecs benefit from
advances in general-purpose neural hardware, where most video codecs require
dedicated hardware to run in real time on device. Third, they have shown to be
able to hallucinate desirable textures, leading to improved perceptual quality
compared to traditional codecs in user studies [19, 18, 38].

###  2.2 Neural video compression

In the video setting, neural codec design has been inspired by techniques from
handcrafted codecs. Where early works used frame interpolation [35] or
predicted entire blocks of frames jointly [9], more recent neural video codecs
exploit similarity between frames using motion compensation and residual
coding. Having neural networks learn to perform these steps end-to-end, as
introduced by [17], has led to major bitrate savings in both the low latency
[2, 24, 13, 23, 11], and streaming setting [22, 14].

As consecutive frames often contain continuous motion, both motion information
and residuals are predictable given previous frames. Recent neural video
codecs exploit this fact, for example, by adding components that predict the
flow [23] or predict the motion and residual latents [16].

In this work, we use simple predictors for both flow and residual. By only
incorporating context from the past two timesteps, and avoiding statefulness
and recurrent components, we avoid instabilities due to aggregating error.
Unlike ELF-VC [23], our codec requires no normalization layers anywhere in the
network, which makes the codec easier to deploy to hardware [21, 1]. Unlike
Liu [16], we make predictions in pixel-space. This makes predictions easy to
interpret, and enables combining our predictors with different base models.

###  2.3 YUV420 colorspace

Most neural video codecs are trained to maximize PSNR or MS-SSIM in the RGB
colorspace. However, when visual quality according to human observers is
important, other colorspaces may be better suited. Distance in the YUV
colorspace is known to match human perception more closely than distance in
the euclidean RGB space [29].

Most standard codecs were designed to maximize PSNR in the YUV420 input
domain. This is a subsampling of YUV, where chroma components are subsampled
by $4\times$ to reduce bitrate and to meet bandwidth constraints. For a
$H\times W$ input, this subsampling results in one $H\times W$ luma channel,
and two $\frac{H}{2}\times\frac{W}{2}$ chroma channels.

There exist a few learned video codecs that support the YUV420 color format
[7, 24], or have been trained to optimize metrics in the YUV domain [23].
However, as of yet, no learned video codec has shown competitive PSNR or MS-
SSIM in the YUV420 domain.

##  3 Method

###  3.1 Base codec

In this work, we introduce flow and residual predictors to capture second-
order redundancy in the low-delay neural video coding setting. Although the
introduced predictors can potentially be added to any neural video codec,
here, we build upon the well-known video codec scale-space flow (SSF) [2] due
to the simple and efficient network design and the straightforward training
schedule. The block-diagram of SSF in shown is Figure 2, where
$\widehat{\mathbf{f}}_{t}$ indicates scale-space flow and
$\operatorname{Warp}$ indicates scale-space warping.

Figure 1 (a) shows two consecutive video frames, and the associated optical
flow and residual maps are shown in Figure 1 (b) and (c), respectively. Of
course the flow and residual for the first timestep depend on a previous
frame, which is omitted here for brevity. As can be seen from these figures,
optical flow maps and residuals are much less detailed than video frames,
making them easier to compress. Additionally, due to the smooth motion in this
video, optical flow and residual maps are highly correlated across the two
time steps.

Figure 2: Base inter-frame codec architecture performing motion compensation
and residual coding. W indicates a scale-space warp.  Figure 3: Architecture
for our inter-frame model. W+ indicates a scale-space warp followed by a
summation.

###  3.2 Predictors

Here, we build a neural video codec on top of SSF where we keep the intra-
frame codec as well the _flow-block_ and _residual-block_ in the inter-frame
codec unchanged. We add flow and residual predictors to the base inter-frame
codec as shown in Figure 3.

Inter-frame coding is done via a four-step process as follows: 1) _flow
prediction_ , 2) _flow compression_ , 3) _residual prediction_ , 4) _residual
compression_. These four steps are explained below:

  1. 1.

Flow prediction: as the name suggests, we predict the flow for the current
time step $\mathbf{f}^{p}_{t}$ using a history of decoded frames,
$\widehat{\mathbf{x}}_{t-1}$ and $\widehat{\mathbf{x}}_{t-2}$, and the
previous flow $\widehat{\mathbf{f}}_{t-1}$, as follows:

$\mathbf{f}^{p}_{t}=\operatorname{Flow-
Pred}(\widehat{\mathbf{x}}_{t-2},\widehat{\mathbf{x}}_{t-1},\widehat{\mathbf{f}}_{t-1})$
(2)

Here, $\operatorname{Flow-Pred}$ is a lightweight network that directly
predicts $\mathbf{f}^{p}_{t}$ using the provided inputs. The predicted flow
$\mathbf{f}^{p}_{t}$ is then used to generate a prediction of the current
frame . This operation is done only on the transmitter side.

$\mathbf{x}_{t}^{p}=\operatorname{Warp}(\widehat{\mathbf{x}}_{t-1},\mathbf{f}^{p}_{t})$
(3)

  2. 2.

Flow compression: in this step, the predicted flow $\mathbf{f}^{p}_{t}$ is
corrected using delta flow $\Delta\widehat{\mathbf{f}}_{t}$ that is measured
and transmitted via $\operatorname{Flow-AE}$:

$\displaystyle\Delta\widehat{\mathbf{f}}_{t}$
$\displaystyle=\operatorname{Flow-AE}(\mathbf{x}_{t}^{p},\mathbf{x}_{t})$ (4)
$\displaystyle\widehat{\mathbf{f}}_{t}$
$\displaystyle=\Delta\widehat{\mathbf{f}}_{t}+\operatorname{Warp}(\mathbf{f}^{p}_{t},\Delta\widehat{\mathbf{f}}_{t})$
(5)

Here, $\operatorname{Flow-AE}$ is a hyperprior network identical to the one in
SSF. The motion compensated frame $\mathbf{x}^{w}_{t}$ is generated using the
corrected flow $\widehat{\mathbf{f}}_{t}$:

$\mathbf{x}_{t}^{w}=\operatorname{Warp}(\widehat{\mathbf{x}}_{t-1},\widehat{\mathbf{f}}_{t})$
(6)

  3. 3.

Residual prediction: we predict the residual for the current time step
$\mathbf{r}^{p}_{t}$ based on the motion compensated frames
$\mathbf{x}^{w}_{t}$ and $\mathbf{x}^{w}_{t-1}$ as well as the previous
decoded frame $\widehat{\mathbf{x}}_{t-1}$:

$\mathbf{r}^{p}_{t}=\operatorname{Res-
Pred}(\mathbf{x}_{t}^{w},\mathbf{x}_{t-1}^{w},\widehat{\mathbf{x}}_{t-1})$ (7)

Here, $\operatorname{Res-Pred}$ is a lightweight network that directly
predicts $\mathbf{r}^{p}_{t}$ using the provided inputs. Given that
$\mathbf{r}_{t}=\mathbf{x}_{t}-\mathbf{x}^{w}_{t}$, the rationale is that
$\operatorname{Res-Pred}$ should be able to predict $\mathbf{r}_{t}$ by
looking at $\widehat{\mathbf{x}}_{t-1}$ and $\mathbf{x}^{w}_{t-1}$. Next,
$\mathbf{r}^{p}_{t}$ is applied to $\mathbf{x}^{w}_{t}$:

$\mathbf{x}_{t}^{r}=\mathbf{x}_{t}^{w}+\mathbf{r}^{p}_{t}$ (8)

  4. 4.

Residual compression: in this step, the predicted residual
$\mathbf{r}^{p}_{t}$ is corrected using delta residual
$\Delta\widehat{\mathbf{r}}_{t}$ that is transmitted via $\operatorname{Res-
AE}$::

$\Delta\widehat{\mathbf{r}}_{t}=\operatorname{Res-
AE}(\mathbf{x}_{t}-\mathbf{x}_{t}^{r})$ (9)

Here, $\operatorname{Res-AE}$ is a hyperprior network identical to the one in
SSF. Finally, the decoded frame $\widehat{\mathbf{x}}_{t}$ is generated as
follows:

$\widehat{\mathbf{x}}_{t}=\mathbf{x}_{t}^{r}+\Delta\widehat{\mathbf{r}}_{t}$
(10)

The following points are worth emphasizing:

  * •

All the used flows are scale-space flow and all the warp operations are scale-
space warp [2].

  * •

Both $\operatorname{Flow-Pred}$ and $\operatorname{Res-Pred}$ have hourglass
architectures. The decoders of $\operatorname{Flow-AE}$ and
$\operatorname{Res-AE}$ are conditioned on the bottleneck features of
$\operatorname{Flow-Pred}$ and $\operatorname{Res-Pred}$, respectively.

  * •

We always encode the first frame in a sequence as an intra-frame and the
subsequent frames as inter-frames. For the first inter-frame $\mathbf{x}_{1}$,
note that $\widehat{\mathbf{x}}_{t-2}$, $\widehat{\mathbf{f}}_{t-1}$, and
$\mathbf{x}_{t-1}^{w}$ are not available. They are set as follows:

$\displaystyle\widehat{\mathbf{x}}_{t-2}$
$\displaystyle=\mathbf{x}_{t-1}^{w}=\widehat{\mathbf{x}}_{t-1}$ (11)
$\displaystyle\widehat{\mathbf{f}}_{t-1}$ $\displaystyle=\textbf{0}$ (12)

  * •

To construct the total flow $\Delta\widehat{\mathbf{f}}_{t}$ in step 2, we
choose to warp the predicted flow, which we observed to lead to better
performance than additive operations. Other schemes have been used as well,
such as letting by the $\operatorname{Flow-AE}$ output a scale-and-shift [23].

Figure 4:  Illustration of video (a), flow (b) and residual (c), flow
prediction and delta flow (d, e), and residual prediction and delta residual
(f, g). Screenshot from Tango video from Netflix Tango in Netflix El
Fuente11footnotemark: 1 [36]

Figure 1 (d) and (f) show the predicted flow and residual, respectively. The
associated delta flow and delta residual are shown in Figure 1 (e) and (g). As
can be seen here, the predicted flow and residual maps are very close to the
actual maps, leaving much less correlated delta flow and delta residual.

Figure 5: a) Network architecture for the RGB colorspace, b) YUV420 input
head, c) YUV420 output head. c* and t* indicate the number of filters in the
convolutional and transposed-convolutional layers, respectively.

###  3.3 RGB and YUV420 architecture

Most literature on neural video codecs provides solutions for the RGB
colorspace. Since the R, G, and B channels share the same resolution, working
with RGB videos is trivial as the networks can conveniently concatenate them
in the network input and output. However, standard codecs were designed to
operate with YUV420 input mainly due to the resemblance with human vision.
Since the chroma channels in the YUV420 colorspace are in a subsampled space,
some modifications to the network architecture are needed to accommodate
YUV420 videos. Figure 5 (a) shows the general architecture we used in the
predictors and the autoencoders of our network for the RGB colorspace. The
input is parsed by an encoder input head E1 and passed through an encoder
backbone E2. The decoder consists of a backbone D2 and a decoder output head
D1. At the bottleneck, the predictors simply pass the encoder outputs directly
to the decoder. Components that transmit information will quantize the latents
before passing these to the decoder.

To enable the codec with accommodating YUV420 inputs/outputs, we update the
input/output heads E1/D1 in all the components of the network as shown in
Figure 5 (b)/(c). E1 and D1 handle the luma (Y) and chroma (U and V) channels
separately, and only the luma channels are down- or upsampled. Lastly, for the
YUV420 codec, warping is applied for the luma and chroma channels separately.
This ensures warping is always performed on the highest available spatial
resolution.

##  4 Experiments

###  4.1 Network architecture

All of the components in our codec follow the architecture shown in Figure 5.
The the number of filters in convolutional and transposed-convolutional layers
is specified by c1, …, c4 and t1, …, t4, respectively. More details are
provided in Table 1. In the main autoencoders – including $\operatorname{Flow-
AE}$, $\operatorname{Res-AE}$, and $\operatorname{Image-AE}$ (used for intra-
frame coding) – the prior model follows the SSF design.

Module |  c1 & t2 |  c2 & t3 |  c3 & t4 | c4
---|---|---|---|---
$\operatorname{Image-AE}$ | 128 | 128 | 128 | 192
$\operatorname{Flow-AE}$ | 128 | 128 | 128 | 192
$\operatorname{Res-AE}$ | 128 | 128 | 128 | 192
$\operatorname{Flow-Pred}$ | 16 | 32 | 64 | 128
$\operatorname{Res-Pred}$ | 16 | 32 | 64 | 128
Table 1: Network architecture details. For all components, t1 depends on the number of output channels. | Module | Parameters | KMACs / pixel
---|---|---|---
RGB | $\operatorname{Image-AE}$ | 9.414M | 30.4% | 198.2 | 30.6%
$\operatorname{Flow-AE}$ | 10.038M | 32.4% | 210.2 | 32.4%
$\operatorname{Res-AE}$ | 10.028M | 32.5% | 207.8 | 32.0%
$\operatorname{Flow-Pred}$ | 0.748M | 2.4% | 16.2 | 2.5%
$\operatorname{Res-Pred}$ | 0.748M | 2.4% | 16.2 | 2.5%
YUV420 | $\operatorname{Image-AE}$ | 9.496M | 30.4% | 208.0 | 30.5%
$\operatorname{Flow-AE}$ | 10.124M | 32.4% | 225.7 | 33.0%
$\operatorname{Res-AE}$ | 10.111M | 32.3% | 217.6 | 31.9%
$\operatorname{Flow-Pred}$ | 0.751M | 2.4% | 16.3 | 2.4%
$\operatorname{Res-Pred}$ | 0.750M | 2.4% | 15.3 | 2.2%
Table 2: Parameters and MACs per pixel of our RGB (top) and YUV420 (bottom)
codecs.

We report the computational complexity of our codec in Table 2, where we show
number of parameters and the number of multiply-and-accumulate operations
(MACs) per module. The MAC count is normalized by number of input pixels to
enable comparison across input modalities and resolutions. We also show the
corresponding percentage to put these numbers in context.

We emphasize that for each I-frame, only the $\operatorname{Image-AE}$ is
used, and this model is not used for P-frames. For example, for a sequence of
length 10, the total number of MACs for the RGB codec would be $1\cdot
198.2+9\cdot(210.2+207.8+16.2+16.1)$ KMACs per pixel.

Using MACs as a proxy for computational complexity, we observe that receiver-
side predictors incur a small computational cost of only 2.5% of the entire
model. Additionally, for each P-frame, both predictors incur $12\times$ fewer
MACs than the P-frame autoencoders. Although the YUV420 input has fewer
dimensions than the RGB input due to the subsampled chroma components, the
fact that the YUV420 codec processes luma and chroma channels using separate
layers results in a higher parameter count and MAC count than that of the RGB
codec.

###  4.2 Dataset

Training is performed on the Vimeo90K dataset [37], a diverse set of 89,800
video sequences of resolution $256\times 448$ in RGB. We evaluate on three
common video compression benchmark datasets: UVG [20], MCL-JCV [33], and class
B of HEVC (Common Test Conditions) [5], all available in raw YUV420 format.
UVG contains 7 full-HD (_i.e_., $1920\times 1080$) videos at framerates of up
to 120 fps. MCL-JCV contains 30 full-HD videos at 30 fps. HEVC class B
contains 5 full-HD videos at various framerates.

###  4.3 Training

We train codecs with batch size of 8, where each sequences has $4$ frames of
size $256\times 256$, resulting in one I-frame and three P-frames. We optimize
our codec on $\operatorname{MSE}$ and $\operatorname{MS-SSIM}$. Specifically,
we initially train all models for 1,000,000 gradient updates on
$\operatorname{MSE}$, then trained the $\operatorname{MS-SSIM}$ models for
extra 200,000 gradient updates on $\operatorname{MS-SSIM}$, and finally
finetuned all the models for 100,000 gradient updates. Training is done on
$256\times 256$ patches with learning rate $10^{-4}$, and finetuning is done
on $256\times 384$ patches with learning rate $10^{-5}$, using the Adam
optimizer.

For the RGB codec training, Video90k is used in the native RGB colorspace, and
both $\operatorname{MSE}$ and $\operatorname{MS-SSIM}$ are averaged over the
R, G, and B channels.

For the YUV420 codec training, Video90k is converted to YUV420 using ffmpeg,
and we weigh the $\operatorname{MSE}$ and $\operatorname{MS-SSIM}$ for the
luma component by a factor $\frac{6}{8}$ and the chroma components by a factor
$\frac{1}{8}$, consistent with the evaluation procedure used for standard
codecs [26].

###  4.4 Evaluation

We evaluate all models in terms of rate-distortion performance. The
$\operatorname{PSNR}$ metric in RGB weighs all color-channels equally. For
YUV420 input, the luma channel is weighed by a factor $\frac{6}{8}$, and both
chroma channels by a factor $\frac{1}{8}$. This is consistent with the way
this metric is used in standard codecs [26]. As there is no commonly agreed-
upon way to measure $\operatorname{MS-SSIM}$ for YUV420 inputs, we simply
follow the same procedure as for $\operatorname{PSNR}$, re-weighting the
$\operatorname{MS-SSIM}$ for the luma and chroma components.

As the evaluation datasets are available in raw YUV420 format, they are used
in the native colorspace for the YUV420 codec evaluation, see Figure 6 (a).
However, evaluating in RGB requires a conversion from YUV420 to RGB using
ffmpeg as shown in Figure 6 (b).

Most standard codecs were designed to work in the YUV420 input domain. To
enable a fair comparison to standard codecs, we always feed them YUV420 inputs
and always operate them in the YUV420 colorspace. The YUV420 output is
converted to the target colorspace (if needed) as shown in Figure 6 (a) and
(b) for comparison with our codecs in YUV420 and RGB, respectively.

Figure 6: Evaluation procedure for YUV420 (a) and RGB (b). Color-space
conversions are performed using ffmpeg.

###  4.5 Results

|  UVG |  HEVC-B |  MCL-JCV
---|---|---|---
Average |  -38.01% |  -26.51% |  -16.48%
Table 3: BD rate savings for RGB PSNR with with respect to the SSF base model. Lower is better. |  UVG |  HEVC-B |  MCL-JCV
---|---|---|---
Y |  -28.34% |  -28.50% |  -21.85%
U |  -47.32% |  -44.38% |  -34.63%
V |  -55.68% |  -43.72% |  -36.51%
Average |  -34.13% |  -32.39% |  -25.28%
Table 4: BD-rate savings for YUV420 PSNR with respect to the SSF-YUV model.
For the average results, the luma and chroma BD-rates are weighed by
$\frac{6}{8}$ and $\frac{1}{8}$ respectively. Lower is better.

Figure 7: RGB rate-distortion curves for the UVG, HEVC Class-B, and MCL-JCV
datasets.

We measure rate-distortion performance on the UVG, MCL-JCV and HEVC-B
datasets. We use $\operatorname{PSNR}$ and $\operatorname{MS-SSIM}$ as
distortion measures, and measure bitrate in bits per pixel. We compare our RGB
and YUV420 codecs to the reference implementation of HEVC called HM-16.25 [10]
as well the the ffmpeg implementation of HEVC and AVC. We always use Group of
Pictures size “infinite“, even if some neural baselines use different GOP
size. HM is used in the low-delay-P configuration, and ffmpeg is used with all
the default parameters but with B-frames disabled. Please see appendix 6.2 for
the used HM and ffmpeg commands. We compare our RGB codec with representative
neural video coding methods including DVC [17], SSF [2], M-LVC [15], FVC [13],
ELF-VC [23], RDAE [9], FRAE [8], and C2F [12]. The RGB and YUV420 rate-
distortions results are shown in Figures 7 and 8, respectively. Here, SSF-YUV
is a re-implementation of SSF with the required YUV420 modifications.

In RGB colorspace, we outperform almost all of the compared methods in terms
of both $\operatorname{PSNR}$ and $\operatorname{MS-SSIM}$, except the very
recent C2F [12]. In YUV420 colorspace, classical codecs dominate on YUV
$\operatorname{PSNR}$, but our codec obtains much better $\operatorname{MS-
SSIM}$ on all the datasets.

Many modern neural codecs use motion compensation and residual coding, and we
emphasize that our method is often applicable to those codecs as well. Here,
we applied our method to the SSF architecture and therefore measured the
performance improvement with respect to a (reimplemented) SSF and SSF-YUV
baseline in Tables 3 and 4, respectively. We see that performance improvements
are substantial, _e.g_., more than 30% on the UVG dataset. Hence, although we
are not generating state-of-the-art results in the RGB colorspace, our method
has the potential to be added to a stronger base method and improve the state-
of-the-art.

Figure 8: YUV420 rate-distortion curves for the UVG, HEVC Class-B, and MCL-JCV
datasets.

###  4.6 Ablations

To separate the effect of each component of our method, we perform ablation
studies. We incrementally add the flow predictor, residual predictor, and
connections between the bottleneck of the predictors and autoencoders
($+$conditions). The corresponding Bjontegaard-Delta [4] (BD)-rate gains are
shown in Table 5. Each addition leads to strong BD-rate gains, and the flow
predictor has the biggest effect for the high-framerate UVG dataset.

Component | UVG | HEVC-B | MCL-JCV
---|---|---|---
base model | 0.0% | 0.0% | 0.0%
+flow predictor | -29.60% | -9.24% | -9.88%
+res predictor | -33.73% | -20.37% | -10.34%
+conditions | -38.01% | -26.51% | -16.48%
Table 5: BD rate savings for RGB PSNR with respect to the SSF base model.
Lower is better.

Additionally, we measure bitrate savings for the flow and residual bitstreams
individually to demonstrate that both $\operatorname{Flow-AE}$ and
$\operatorname{Res-AE}$ benefit from our method. To measure BD-rate per
component, we use the output PSNR as distortion and use the individual
bitstream sizes per component as rate. The BD-rate gains versus PSNR for the
UVG, HEVC class B, and MCL-JCV datasets are shown in Figure 9. Here, it can be
observed that the predictors deliver significant rate savings on both flow and
residual bitstreams across all PSNR values.

For completeness, we list two training techniques that did not work for us in
early experiments and were therefore abandoned. Training with auxiliary losses
for reconstructions resulting from flow and residual prediction,
$\mathbf{x}_{t}^{p}$ and $\mathbf{x}_{t}^{r}$, did not work better than using
only a distortion loss defined on the reconstructions
$\widehat{\mathbf{x}}_{t}$. Temporal loss modulation as in [23] resulted in
$<1\%$ rate savings, so we omitted it for simplicity.

Figure 9: BD-rate savings versus RGB PSNR for flow and residual bitstreams and
in-total for a SSF baseline. Lower is better.

##  5 Conclusion

In this work, we demonstrate that improved use of receiver-side information
substantially improves the rate-distortion performance of neural video codecs.
By having the decoder predict the flow and residual at a current timestep
based on previously transmitted information, our codec only needs to transmit
a correction, saving bits in the process. This results in Bjontegaard-Delta
rate savings of up to 30% compared to a Scale-Space Flow baseline on common
video benchmark datasets. The required architectural modification is
straightforward, requires only lightweight components, and is applicable to
other neural video codecs as well.

Furthermore, we train our model on both the RGB and YUV420 input domains with
only minor architectural modifications, and show strong performance in both
settings compared to neural and standard baselines. As standard codecs were
designed to optimize YUV420 PSNR, our codec does not outperform HM on this
metric in the low bitrate regime, but substantially outperform the HM codec on
YUV420 MS-SSIM. Ablations show that both the flow and residual predictors
contribute to the final performance.

#### Impact statement

It is possible that learned codecs exacerbate biases present in the training
data. Nevertheless, we believe that improved video coding efficiency has a net
positive impact on the world by reducing bandwidth and storage space needs,
and by improving visual quality for video applications such as video
conferencing.

## References

  * [1] Low-Power Computer Vision: Improve the Efficiency of Artificial Intelligence (1st ed.).  Chapman and Hall/CRC, 2022.
  * [2] Eirikur Agustsson, David Minnen, Nick Johnston, Johannes Balle, Sung Jin Hwang, and George Toderici.  Scale-space flow for end-to-end optimized video compression.  In Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, 2020.
  * [3] Eirikur Agustsson, Michael Tschannen, Fabian Mentzer, Radu Timofte, and Luc Van Gool.  Generative adversarial networks for extreme learned image compression.  In Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, 2019.
  * [4] Gisle Bjøntegaard.  Calculation of average psnr differences between rd-curves.  2001\.
  * [5] Frank Bossen et al.  Common test conditions and software reference configurations.  JCTVC-L1100, 12(7), 2013.
  * [6] Cisco.  Report: Where does the majority of internet traffic come from?  2017\.
  * [7] Hilmi E. Egilmez, Ankitesh K. Singh, Muhammed Coban, Marta Karczewicz, Yinhao Zhu, Yang Yang, Amir Said, and Taco S. Cohen.  Transform network architectures for deep learning based end-to-end image/video coding in subsampled color spaces.  IEEE Open Journal of Signal Processing, 2:441–452, 2021.
  * [8] Adam Goliński, Reza Pourreza, Yang Yang, Guillaume Sautière, and Taco S. Cohen.  Feedback recurrent autoencoder for video compression.  ACCV, 2020.
  * [9] Amirhossein Habibian, Ties van Rozendaal, Jakub M Tomczak, and Taco S Cohen.  Video compression with rate-distortion autoencoders.  In IEEE International Conference on Computer Vision, 2019.
  * [10] HM.  Hevc offical test model.  https://hevc.hhi.fraunhofer.de,2021, 2021.
  * [11] Zhihao Hu, Guo Lu, Jinyang Guo, Shan Liu, Wei Jiang, and Dong Xu.  Coarse-to-fine deep video coding with hyperprior-guided mode prediction.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 5921–5930, 2022.
  * [12] Zhihao Hu, Guo Lu, Jinyang Guo, Shan Liu, Wei Jiang, and Dong Xu.  Coarse-to-fine deep video coding with hyperprior-guided mode prediction.  In Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, 2022.
  * [13] Zhihao Hu, Guo Lu, and Dong Xu.  FVC: A new framework towards deep video compression in feature space.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 1502–1511, 2021.
  * [14] Théo Ladune, Pierrick Philippe, Wassim Hamidouche, Lu Zhang, and Olivier Déforges.  Conditional coding for flexible learned video compression.  ICLR neural compression workshop, 2021.
  * [15] Jianping Lin, Dong Liu, Houqiang Li, and Feng Wu.  M-lvc: Multiple frames prediction for learned video compression.  Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, 2020.
  * [16] Chao Liu, Heming Sun, Jiro Katto, Xiaoyang Zeng, and Yibo Fan.  Learned video compression with residual prediction and loop filter.  Aug. 2021.
  * [17] Guo Lu, Wanli Ouyang, Dong Xu, Xiaoyun Zhang, Chunlei Cai, and Zhiyong Gao.  Dvc: An end-to-end deep video compression framework.  In Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, 2019.
  * [18] Fabian Mentzer, Eirikur Agustsson, Johannes Ballé, David Minnen, Nick Johnston, and George Toderici.  Neural video compression using gans for detail synthesis and propagation.  arXiv preprint arXiv:2107.12038, 2021.
  * [19] Fabian Mentzer, George Toderici, Michael Tschannen, and Eirikur Agustsson.  High-fidelity generative image compression.  Neural Information Processing Systems, 2020.
  * [20] Alexandre Mercat, Marko Viitanen, and Jarno Vanne.  Uvg dataset: 50/120fps 4k sequences for video codec analysis and development.  In Proceedings of the 11th ACM Multimedia Systems Conference, MMSys ’20, page 297–302, New York, NY, USA, 2020. Association for Computing Machinery.
  * [21] Markus Nagel, Marios Fournarakis, Rana Ali Amjad, Yelysei Bondarenko, Mart van Baalen, and Tijmen Blankevoort.  A white paper on neural network quantization, 2021.
  * [22] Reza Pourreza and Taco S Cohen.  Extending neural p-frame codecs for b-frame coding.  IEEE International Conference on Computer Vision, 2021.
  * [23] Oren Rippel, Alexander G Anderson, Kedar Tatwawadi, Sanjay Nair, Craig Lytle, and Lubomir Bourdev.  ELF-VC: Efficient Learned Flexible-Rate Video Coding.  Neural Information Processing Systems, 2021.
  * [24] Ankitesh K Singh, Hilmi E Egilmez, Reza Pourreza, Muhammed Coban, Marta Karczewicz, and Taco S Cohen.  A combined deep learning based End-to-End video coding architecture for YUV color space.  Apr. 2021.
  * [25] Yannick Strümpler, Janis Postels, Ren Yang, Luc Van Gool, and Federico Tombari.  Implicit neural representations for image compression.  arXiv preprint arXiv:2112.04267, 2021.
  * [26] Jacob Ström, Kenneth Andersson, Rickard Sjöberg, Andrew Segall, Frank Bossen, Gary Sullivan, Jens-Rainer Ohm, and Alexis Tourapis.  Working practices using objective metrics for evaluation of video coding efficiency experiments (draft 4).  Input document JVET-T2016-v2, Joint Video Exploration Team (JVET) of ITU-T SG 16 WP 3 and ISO/IEC JTC 1/SC 29, Oct. 2020.
  * [27] Gary J. Sullivan, Jens-Rainer Ohm, Woo-Jin Han, and Thomas Wiegand.  Overview of the high efficiency video coding (hevc) standard.  IEEE Transactions on Circuits and Systems for Video Technology, 22(12):1649–1668, 2012.
  * [28] Lucas Theis, Wenzhe Shi, Andrew Cunningham, and Ferenc Huszár.  Lossy image compression with compressive autoencoders.  International Conference on Learning Representations, 2017.
  * [29] M Tkalcic and J F Tasic.  Colour spaces: perceptual, historical and applicational background.  In The IEEE Region 8 EUROCON 2003. Computer as a Tool., volume 1, pages 304–308 vol.1, Sept. 2003.
  * [30] George Toderici, Damien Vincent, Nick Johnston, Sung Jin Hwang, David Minnen, Joel Shor, and Michele Covell.  Full resolution image compression with recurrent neural networks.  In Proceedings of the IEEE conference on Computer Vision and Pattern Recognition, pages 5306–5314, 2017.
  * [31] Ties van Rozendaal, Johann Brehmer, Yunfan Zhang, Reza Pourreza, and Taco S Cohen.  Instance-adaptive video compression: Improving neural codecs by training on the test set.  arXiv preprint arXiv:2111.10302, 2021.
  * [32] Ties van Rozendaal, Iris A. M. Huijben, and Taco S. Cohen.  Overfitting for fun and profit: Instance-adaptive data compression, 2021\.
  * [33] Haiqiang Wang, Weihao Gan, Sudeng Hu, Joe Yuchieh Lin, Lina Jin, Longguang Song, Ping Wang, Ioannis Katsavounidis, Anne Aaron, and C.-C. Jay Kuo.  Mcl-jcv: A jnd-based h.264/avc video quality assessment dataset.  In IEEE International Conference on Image Processing, pages 1509–1513, 2016.
  * [34] Thomas Wiegand, Gary J Sullivan, Gisle Bjontegaard, and Ajay Luthra.  Overview of the h. 264/avc video coding standard.  IEEE Transactions on circuits and systems for video technology, 13(7):560–576, 2003.
  * [35] Chao-Yuan Wu, Nayan Singhal, and Philipp Krähenbühl.  Video compression through image interpolation.  Proceedings of the European Conference on Computer Vision, 2018\.
  * [36] Xiph.org.  Xiph.org video test media [derf’s collection].  https://media.xiph.org/video/derf/.  Accessed: 2020-02-21.
  * [37] Tianfan Xue, Baian Chen, Jiajun Wu, Donglai Wei, and William T Freeman.  Video enhancement with task-oriented flow.  International Journal of Computer Vision, 127(8):1106–1125, 2019\.
  * [38] Ren Yang, Luc Van Gool, and Radu Timofte.  Perceptual learned video compression with recurrent conditional GAN.  arXiv preprint arXiv:2109.03082, Sept. 2021.

##  6 Supplementary material

###  6.1 Commands for standard codecs

We use the following commands to obtain compressed videos for standard codecs.
For ffmpeg, we use the medium preset, disable B-frames and otherwise use
default settings. For HM, we use HM-16.25 with the LowDelay-P config, which
can be obtained from https://vcgit.hhi.fraunhofer.de/jvet/HM/-/tags/HM-16.25.

⬇

# For ffmpeg x264

ffmpeg -y -f rawvideo \

-pix_fmt yuv420p \

-s:v <width>x<height> \

-r <framerate> \

-i <input.yuv> \

-c:v libx264 \

-preset medium \

-crf <crf> \

-x264-params bframes=0 \

<output>

# For ffmpeg x265

ffmpeg -y -f rawvideo \

-pix_fmt yuv420p \

-s:v <width>x<height> \

-r <framerate> \

-i <input.yuv> \

-c:v libx265 \

-preset medium \

-crf <crf> \

-x265-params bframes=0 \

<output>

# For HM-16.25 LowDelayP

./bin/TAppEncoderStatic -c \

./cfg/encoder_lowdelay_P_main.cfg \

-i <input.yuv> \

–InputBitDepth=8 \

-wdt <width> \

-hgt <height> \

-fr <framerate> \

-f <numframes> \

-q <qp> \

-o <output>

###  6.2 Detailed computational complexity

We report the number of model parameters and multiply-and-accumulate (MAC)
operations for the RGB and YUV codec in Table 6, separating the numbers for
the sender (“send”) and receiver (“recv”) side. MACs are normalized with
respect to the number of input pixels to enable comparisons across input
resolutions and input modalities. For each group of five rows, the shown
percentage is computed with respect to that group. All percentages are
rounded, and some columns do not sum to 100% as a result.

Note that in order to transmit data, the modules are run in their entirety due
to the design of the SSF codec. Compute for the sender side is therefore
identical to the total compute as reported in Table 2 in the main text, and we
repeat those numbers here for completeness.

On the receiver side, the predictors are run in their entirety, but only the
decoder and hyper-decoder of each of the hyperprior modules has to be run. We
see that the computational complexity of our predictors is low on the receiver
side still: 7.4% of total parameters, 6.2% of total MACs for the RGB receiver,
and only 5.6% for the YUV420 codec.

As a minor detail, for the SSF model, the residual decoder does not have to be
run on the sender side for the very last frame in a Group of Pictures. As this
has a negligible effect on compute when coding with a large GoP size we do not
take this into account here.

| Module | Parameters | KMACs / pixel
---|---|---|---
RGB send | $\operatorname{Image-AE}$ | 9.414M | 30.4% | 198.2 | 30.6%
$\operatorname{Flow-AE}$ | 10.038M | 32.4% | 210.2 | 32.4%
$\operatorname{Res-AE}$ | 10.028M | 32.5% | 207.8 | 32.0%
$\operatorname{Flow-Pred}$ | 0.748M | 2.4% | 16.2 | 2.5%
$\operatorname{Res-Pred}$ | 0.748M | 2.4% | 16.2 | 2.5%
RGB recv | $\operatorname{Image-AE}$ | 5.790M | 28.8% | 158.9 | 30.1%
$\operatorname{Flow-AE}$ | 6.410M | 31.9% | 168.6 | 31.9%
$\operatorname{Res-AE}$ | 6.410M | 31.9% | 168.5 | 31.9%
$\operatorname{Flow-Pred}$ | 0.748M | 3.7% | 16.2 | 3.1%
$\operatorname{Res-Pred}$ | 0.748M | 3.7% | 16.2 | 3.1%
YUV420 send | $\operatorname{Image-AE}$ | 9.496M | 30.4% | 208.0 | 30.5%
$\operatorname{Flow-AE}$ | 10.124M | 32.4% | 225.7 | 33.0%
$\operatorname{Res-AE}$ | 10.111M | 32.3% | 217.6 | 31.9%
$\operatorname{Flow-Pred}$ | 0.751M | 2.4% | 16.3 | 2.4%
$\operatorname{Res-Pred}$ | 0.750M | 2.4% | 15.3 | 2.2%
YUV420 recv | $\operatorname{Image-AE}$ | 5.880M | 28.9% | 171.1 | 29.9%
$\operatorname{Flow-AE}$ | 6.500M | 31.9% | 188.8 | 33.0%
$\operatorname{Res-AE}$ | 6.490M | 31.9% | 180.7 | 31.6%
$\operatorname{Flow-Pred}$ | 0.751M | 3.7% | 16.3 | 2.9%
$\operatorname{Res-Pred}$ | 0.750M | 3.7% | 15.3 | 2.7%
Table 6: Parameters and MACs per pixel of our RGB and YUV420 codecs.