[a]Aix Marseille Univ, CNRS/IN2P3, CPPM, Marseille, France [b]Aix Marseille
Univ, CNRS/IN2P3, CPPM, IPhU, Marseille, France<EMAIL_ADDRESS>
# Firmware implementation of a recurrent neural network for the computation of
the energy deposited in the liquid argon calorimeter of the ATLAS experiment

G. Aad,Corresponding author  T. Calvet,  N. Chiedde  R. Faure  E.M. Fortin  L.
Laatu  E. Monnier  N. Sur

###### Abstract

The ATLAS experiment measures the properties of particles that are products of
proton-proton collisions at the LHC. The ATLAS detector will undergo a major
upgrade before the high luminosity phase of the LHC. The ATLAS liquid argon
calorimeter measures the energy of particles interacting electromagnetically
in the detector. The readout electronics of this calorimeter will be replaced
during the aforementioned ATLAS upgrade. The new electronic boards will be
based on state-of-the-art field-programmable gate arrays (FPGA) from Intel
allowing the implementation of neural networks embedded in firmware. Neural
networks have been shown to outperform the current optimal filtering
algorithms used to compute the energy deposited in the calorimeter. This
article presents the implementation of a recurrent neural network (RNN)
allowing the reconstruction of the energy deposited in the calorimeter on
Stratix 10 FPGAs. The implementation in high level synthesis (HLS) language
allowed fast prototyping but fell short of meeting the stringent requirements
in terms of resource usage and latency. Further optimisations in Very High-
Speed Integrated Circuit Hardware Description Language (VHDL) allowed
fulfilment of the requirements of processing 384 channels per FPGA with a
latency smaller than 125 ns.

Calorimeter methods, Digital signal processing (DSP), Data processing methods.

##  1 Introduction

The ATLAS experiment [1] at the Large Hadron Collider (LHC) [2] measures the
properties of particles produced in proton-proton collisions at energies of
several teraelectronvolts [TeV] and a collision frequency of 40 MHz. In the
years 2026-2029 the LHC will undergo a major upgrade to increase its
instantaneous luminosity by a factor of 5-7 leading to the High luminosity LHC
(HL-LHC). During the same period the ATLAS detector will be upgraded to cope
with the increased luminosity of the HL-LHC. This upgrade is called the phase-
II upgrade. The readout electronics of the ATLAS liquid argon (LAr)
calorimeter will be replaced as part of the phase-II upgrade [3].

The LAr calorimeter measures the energy of particles produced in LHC
collisions. An excellent energy resolution and accurate detection of the
energy-deposit time are crucial to enhance the ATLAS physics discovery
potential at the HL-LHC. Currently, the transverse energy is computed using
optimal filtering algorithms [4] that assume a nominal pulse shape of the
electronic signal. Calorimeter electronic signals of up to 25 subsequent
collisions overlap and create distortions to the pulse shape. This increases
the difficulty of energy reconstruction and identification of the
corresponding proton-proton bunch crossing.

Up to 200 simultaneous proton-proton collisions are expected at the HL-LHC,
which will lead to a high rate of overlapping signals in a given calorimeter
channel. This will result in a significant energy degradation especially for
low time-gap between two consecutive pulses, as discussed in [5].

The use of modern energy reconstruction algorithms based on neural networks
(NNs) improves the resolution especially in the low time-gap region as shown
in [5]. Both convolutional neural networks (CNNs) and recurrent neural
networks (RNNs) are shown to outperform the optimal filtering algorithm. The
transverse energy reconstruction is performed in custom electronic boards
based on the latest state-of-the-art FPGAs [3]. Each FPGA should reconstruct
the energies of 384 independent channels. The computation is done on-the-fly
at the collision frequency of 40 MHz. Furthermore, the energy reconstruction
should be done within a latency of 125 ns. These requirements put very
stringent constraints on the firmware implementation of the energy
reconstruction algorithms. The implementation should balance serialisation of
the channels to save FPGA computation resources, and parallelisation to reduce
the implementation latency and keep up with the high bandwidth.

The electronic boards for the LAr calorimeter phase-II upgrade are currently
under development and will use the Agilex [6] family of Intel FPGAs. A
demonstrator board is already produced with Stratix 10 FPGAs [7]. In this
article we consider the implementation of a vanilla RNN [8] algorithm in a
Stratix 10 FPGA from Intel (part number 1SG280HU1F50E2VG). The implementation
is carried out initially in HLS for fast prototyping and optimisation of the
network parameters and the implementation architecture. After this initial
implementation, VHDL is used for the final optimisation, ensuring that the
requirements in terms of FPGA resource occupancy and latency are met.

##  2 Vanilla RNN

The neural network used in this article is the vanilla RNN with the sliding-
window approach described in [5]. At each bunch crossing, the network receives
as input five consecutive samples of the electronic pulse of one channel of
the LAr calorimeter and computes the corresponding deposited transverse
energy. The five samples correspond to a time window of five bunch crossings.
The computed transverse energy corresponds to a possible deposit at the second
bunch crossing. Four of the samples are around the pulse peak generated by the
deposited energy while the first sample is prior to the pulse allowing to
detect the presence of overlapping pulses from past energy deposits.

The network implementation is composed of five consecutive RNN cells followed
by a dense layer as shown in figure 1. Each cell process the input from one
sample corresponding to one bunch crossing. Each of the cells is composed of 6
blocks of computation as shown in figure 2: two addition blocks, one
multiplication of a scalar with a vector, one multiplication of a vector with
a matrix, and one activation function. The weights of the network obtained
from the training are stored in memory and are given as input to the RNN cells
along with the electronic pulse samples. The same weights are used for each of
the five network cells. The ReLU, described in equation 1, is used as the
activation function for its simple implementation in FPGAs.

$f(x)=\begin{cases}0&\text{if $x\leq 0$}\\\ x&\text{otherwise}\end{cases}$ (1)

Vanilla RNN 5 cellsEnergy reconstructedFirst cellSecond cellThird cellForth
cell Fifth cell Dense layerSignal from one LAr channel/8/8/8/8/8/1

Figure 1: Schematics of a Vanilla RNN with 5 cells and 8 internal dimensions
(size of the state vector) followed by a dense layer. Each cell calculates a
new state using the input data and the state of the previous cell. As the
first cell does not have a previous cell, the calculation is only based on the
input data. The dense layer computes the transverse energy using the state of
the last cell.

Figure 2: Schematics of the operations performed inside each of the RNN cells.
The recurrent weight multiplication multiplies the state vector from the
previous cell (St) with the recurrent kernel weight matrix (R).
Simultaneously, the LAr cell input (Xt+1) is multiplied by the kernel weight
vector (W) and added to the bias weight (B). The results from the two above
operations are added to create the internal vector Tt+1. The ReLU activation
function is applied on the elements of Tt+1 to create the state vector St+1.
The state vector size (n) is equal to 8.

##  3 High-Level Synthesis Implementation

The High Level Synthesis (HLS) is a design process that takes as input a
behavioural specification of a digital system and convert it to a register-
transfer level (RTL) realising the described functions. For the study
presented in this article, the Intel HLS [9] is used. The Intel HLS takes a
C++ code as input and produces RTL optimised for Intel FPGAs. It provides
several macros which are inlined in the C++ code that allow control over the
RTL implementation. The Intel HLS implements the standard C++ types but also
defines several other types, in particular arbitrary precision fixed-point
representations.

Neural network cell computations are made of additions and multiplications.
The aim of the firmware implementation is to reproduce with high precision the
software computation results while keeping low resource usage in the FPGA and
low latency.

###  3.1 Implementation of multiplications

Inside the Vanilla RNN cell there is one vector multiplication and one matrix
multiplication which reduces to several vector multiplications which in turn
reduce to several scalar multiplications and scalar additions. There is a
dedicated component inside the FPGA to perform scalar multiplications called
the Digital Signal Processing (DSP). The DSP can be used in three possible
modes. The first mode performs one 32$\times$32 bits multiplication in the
floating-point representation. The second mode performs one 27$\times$27 bits
multiplication in the fixed-point representation. The third mode can perform
two independent 19$\times$18 bits multiplications in the fixed-point
representation. As explained in 3.2, the number of multiplications of the FPGA
limits the number of calorimeter cells that can be handled by one FPGA. The
third mode is thus chosen since it allows doubling the available dedicated
multiplication resources on the FPGA.

###  3.2 Multiplexing

The number of multiplications inside the neural network depends on the size of
the state vector ($n$) and the number of network cells ($c$) following
equation 2 :

$n^{2}\times(c-1)+n\times(c+1)$ (2)

where $n=8$ and $c=5$. Therefore, 304 multiplications are needed for the
vanilla RNN used in this article. This leads to 116,736 multiplications that
are needed to handle the 384 channels processed in one FPGA. The Stratix 10
FPGA possesses only 11,520 multiplication blocks. Therefore, each computation
block must be utilised for several channels. The easiest way to do so is to
increase the DSP block computation frequency as a multiple of the input
frequency (40MHz). Several channels are then multiplexed using the same block.
This multiplexing is limited by the maximum frequency of the computation
block.

###  3.3 Implementation of arithmetic operations

Intel HLS implements two binary representations of numbers: a fixed point
representation and a floating point representation. The implementation of
floating point representation is relatively complex and uses more arithmetic
and logic blocks for the addition operations. Furthermore, the DSP blocks of
the Stratix 10 FPGA allows only one floating point multiplication instead of
two simultaneous fixed point multiplications as explained in section 3.1. The
fixed point representation is chosen in order to minimise the resource usage
in the FPGA.

The fixed point representation implemented in Intel HLS follows the
Algorithmic C (AC) Data-types defined by Mentor Graphics [10]. Four parameters
define this representation: the total number of bits, the number of bits of
the integer part, the quantisation type, and the treatment of the overflow.
Three different data categories are defined for the vanilla RNN
implementation: the input and output data, the neural network weights, and the
intermediate data which are the results of the internal computation inside the
neural network blocks. 19 bits are used for the width of the internal and
input/output categories, while 16 bits are used for the weights. This ensures
an efficient use of the DSP resources inside the FPGA while keeping a good
(¡0.1%) resolution of the firmware computation with respect to the software
computation that uses floating points arithmetic operations.

The overflow defines how the bits to the left of the most significant bit
(MSB) are lost due to saturation. The number of bits and the position of the
radix point are chosen to be able to represent the maximum value that could
occur in the network. Thus, no saturation detection is needed and a simple
drop of bits implementation is used.

The quantisation defines how the bits to the right of the least significant
bit (LSB) are lost. The loss of bits can occur during the conversion of the
floating points representing the input and weights. The inputs and weights are
given by a simulation111With real data the input will be a 16-bit integer
coming from an ADC and therefore this conversion will be implicitly done in
the detector electronics. of the detector and by the neural network training
software. These two use 32-bit signed floating point representations that need
to be converted to a fixed point representation in order to be used in the
firmware. The loss of bits can also occur inside the neural network internal
computations. The 37-bit output of the DSP is reduced to the number of bits
internally used. Two types of quantisation are implemented in the Intel HLS
compilation: truncation and rounding. Each of these types possess different
subtypes which are explained in [10].

Figure 3 shows a comparison of the resource usage, in terms of DSPs,
arithmetic lookup tables (ALUT), Flip Flops (FF), and random access memory
(RAM), and the latency of the firmware for different implementations of the
quantisation. The same quantisation is applied for all categories of data of
the network. Two modes are interesting, the default truncation (TRN) and the
default rounding (RND) which use lower resources and lower latency compared to
other quantisation modes. Figure 4 shows a comparison of the transverse energy
computed in firmware and the one computed in software, with a full floating
point implementation, for the different quantisation modes. All quantisation
modes show similar resolution with the exception of the TRN quantisation modes
which have large tails. The RND mode gives a good compromise among resource
usage, latency, and resolution.

Figure 3: Comparison of the resource usage (left) and the latency (right) of
the vanilla RNN firmware with different implementation of the quantisation in
Intel HLS. For each mode, the same quantisation is applied for all categories
of the data of the network. The different quantisation modes are described in
[10].

Figure 4: Resolution of the transverse energy (ET) computed in firmware with
respect to the one computed in software. The ET(firmware) is computed with
different implementations of the quantisation in Intel HLS. In each mode, the
same quantisation is applied for all categories of the data of the network.
The different quantisation modes are described in [10]. The different RND
modes give very similar results and their corresponding curves overlap. A
lower cut of 240 MeV is applied on $E_{\mathrm{T}}$(software) to remove low
energies below the 3$\sigma$ noise level as described in [5].

To further optimise the firmware implementation, a mix of quantisation
procedures are used for different data categories. Figure 5 shows the resource
usage and the latency while applying the TRN and RND quantisations on
different data categories. The rounding of the weights does not require any
additional resources in the FPGA since it can be done in software before
loading these weights into the FPGA. The same applies to the input data.
Rounding the output data induces a slight increase in the latency. Rounding
the internal data categories increases the needed resources and the latency
significantly. Figure 6 shows a comparison of the transverse energy computed
in firmware and the one computed in software depending on which of the data
categories is rounded. One can see that it is important to round the weights
and input/output categories, while rounding the internal data category does
not have any significant impact on the resolution. The root mean square (RMS)
of the TRN distribution is 0.2%. It decreases to 0.07% if all categories are
rounded (RND_IWD). The RMS becomes slightly worse (0.09%) if only the weights
and the input/output are rounded (RND_WD). Therefore, the TRN mode is used for
the internal data computation while the RND is used for the other data
categories. These optimisations allow significant improvement of the
resolution of the network at a low resource and latency cost.

Figure 5: Comparison of the resource usage (left) and the latency (right) of
the vanilla RNN firmware with different implementation of the quantisation for
the three data categories of the network computation. For each test, the
letters I, W, and D indicate that the RND is applied for the internal data
category, the weights, and the input/output data, respectively, while the TRN
mode is applied by default in all other categories.

Figure 6: Resolution of the transverse energy (ET) computed in firmware with
respect to the one computed in software. The ET(firmware) is computed with
different implementations of the quantisation for the three data categories of
the network. For each test, the letters I, W, and D indicate that the RND is
applied for the internal data category, the weights, and the input/output
data, respectively, while the TRN mode is applied by default in all other
categories. A lower cut of 240 MeV is applied on $E_{\mathrm{T}}$(software) to
remove low energies below the 3$\sigma$ noise level as described in [5].

###  3.4 Implementation of the neural network

Figure 7: Resource usage for ALUTs and FFs depending on the frequency for two
implementations of a matrix multiplication with and without chained DSPs. The
maximum frequency of an MLAB is 450 MHz. Above this frequency, there is a
significant increase of the logic resource usage for the chained DSP
implementation.

The computation inside the neural network requires 304 multiplications and 231
additions. The multiplications are implemented inside the DSPs, and the
additions are initially implemented using ALUTs and FFs. The DSP allows
summing the output of the two multiplications internally. The usage of such
functionality reduces the number of additions implemented in the ALUTs and FFs
from 231 to 131. Furthermore, the DSP component contains an additional adder
that can take an external input. It is possible to chain two DSPs to sum their
outputs by using the output of the first DSP in the additional adder of the
second DSP. By doing this, it is necessary to synchronise the DSPs so that
their results arrive at the same time to the additional adder of the second
DSP. Therefore, the second DSP in the chain should be shifted by 1 clock
cycle. The same procedure is applied to more DSPs to build a full chain. In
such case, all additions can be implemented in DSPs.

To perform the timing shift, each input of the DSPs needs an additional level
of registers to delay the results. This is done by combining several ALUTs to
create a memory logic array block (MLAB) to implement a first in first out
(FIFO) memory. However, the MLAB frequency is limited to 450 MHz in read-
during-write mode needed for FIFO implementation. At higher frequencies the
synchronisation cannot be implemented in MLAB, in such case the delays are
implemented in basic ALUTs and FFs which increase the number of needed logic
elements significantly. At high frequency, more ALUTs and FFs are needed to
synchronise the DSPs in chained mode than to implement the additions in a non-
chained mode. Figure 7 shows the FPGA resource usage as function of the
frequency for DSPs used in chained or non-chained mode for one matrix
multiplication. The chained mode is advantageous below 450 MHz. Above this
frequency, it is more advantageous to use logic elements to do the additions
than to chain DSPs. This will be the option retained for the implemented RNN
since we seek higher frequencies to increase the multiplexing and thus reduce
the overall resource usage.

###  3.5 Results of the HLS Implementation

The designed RNN HLS code is compiled in Intel HLS and Quartus [11]. The
results of these two compilations are given in table 1. For one implemented
network the design can run at 455 MHz which a priory allows multiplexing up to
11 channels per network for a data input rate of 40 MHz. However, with a
multiplexing of 10 the maximum frequency reached is 393 MHz while the needed
frequency is 400 MHz for the firmware to run without timing violations. The
maximum frequency is reduced when applying multiplexing due to the additional
weights that are needed by the network to perform the computation for several
channels. Up to 37 networks can fit in the FPGA which leads to the usage of
100% of the available DSP resources and most of the logic resources. The logic
resources are given in terms of ALUTs and FFs in the HLS report and adaptive
logic modules (ALM) for the Quartus report. The ALMs are the actual physical
resources in the FPGA and each ALM can be configured to be used as two ALUTs
or 4 FFs. The HLS design does not allow reaching the required 384 channels
processed in one FPGA, even with full utilisation of resources in the FPGA. In
practice only part of these resources will be available for the transverse
energy computation as discussed in section 4.4. Additional optimisations are
needed to fulfil the requirements of the LAr calorimeter. Furthermore, the
latency of this HLS firmware is 277.5 ns which is significantly larger than
the required 125 ns. The HLS design was optimised for the maximum possible
frequency which leads to additional registers added by the compiler to meet
the required timing constraints. This in turn increases the latency of the
design.

Table 1: Resource usage and the maximum frequency (Fmax) given by the Intel
HLS and Quartus compiler reports for one and 37 implemented networks in the
FPGA. The multiplexing is set to 10 so that each network instance can process
10 independent channels.

| Networks | Multiplexing | Channels | ALUTs | FFs | ALMs | DSPs | Memory | Fmax
---|---|---|---|---|---|---|---|---|---
| 1 | 10 | 10 | 0,6% | 0.7% | - | 3% | 0,3% | -
HLS | 37 | 10 | 370 | 22.9% | 25.9% | - | 100% | 11% | -
| 1 | 10 | 10 | - | - | 2,4% | 3% | 0,3% | 455 MHz
Quartus | 37 | 10 | 370 | - | - | 90% | 100% | 11% | 393 MHz

##  4 VHDL Implementation

The HLS implementation adds an additional level of abstraction that allows
fast and efficient optimisation of the network parameters and firmware
implementation. However, this additional level of abstraction prevents some
finer optimisations that are possible in VHDL. These finer optimisations allow
to meet the specifications, which is not possible for the HLS implementation.
That is why VHDL is used for the final optimisation and placement of the RNN
firmware implementation.

###  4.1 Reuse of common computations between RNN cells

As shown in figure 2, the Kernel Weight Multiplication and the Bias Addition
depend only on the input of the NN and do not depend on the past network
state. However, these two computations done in the second cell at a given time
$t$ are identical to the computations in the first cell at time $t-1$. More
generally, these computations in cell $n$ at time $t$ are identical to the
computations in the first cell at time $t-n-1$. One can do these two
computations one time at the first cell and propagate the output to the other
four cells at the proper time. Moreover, the first cell does not need the
Recurrent Weight Matrix Multiplication since there is no previous RNN state at
this level. Therefore, this computation is removed from the first cell. These
two optimisations allow a reduction in the number of used DSPs by 10%. The ALM
usage is also reduced by 21% due to the removal of the duplicated Bias
Additions.

###  4.2 Placement constraints

Several instances of the neural network are needed to process all 384
channels. In a given compilation each instance have a different placement
shape. Moreover, these shapes change between compilations due to the
randomisation in Quartus. This complicates the optimisation of the timing
critical paths that is needed to reach higher frequencies and thus higher
multiplexing.

Placement constraints are used to force the same placement shape of the
implemented neural network. Thus, all instances of the neural network have the
same shape which simplifies the optimisation of the critical paths. Moreover,
the placement of the 5 network cells is optimised to minimise the distance
between connected cells. The shape of the neural network is shown in figure 8.
The first cell is placed in the middle since it is connected to all other four
cells as described in section 4.1. The other cells are placed around the first
cell and ordered to reduce the distance between consecutive cells. The dense
layer is placed next to the fifth cell.

Second cellThird cellForth cellFifth cellDense layerFirst cell

DSP Line

DSP Line

Memory Line

Figure 8: Schematic view of the optimised placement of the Vanilla RNN showing
the 5 cells and the dense layer with respect to the memory and DSP lines
inside the FPGA.

All cells use the same set of recurrent weights except for the first cell that
does not have a recurrent block as explained in section 4.1. These weights are
stored in memory (M20K blocks) and are directly connected to the DSPs to make
the matrix multiplications. To reduce the mean distance between the M20k block
and each DSP the recurrent weights are duplicated. Each set of weights is used
for two cells instead of four.

These placement constraints allow increasing the maximum frequency from 434
MHz to 492 MHz with 28 instances of the RNN implemented in the FPGA.

###  4.3 Incremental compilation

Quartus compilation is composed of three parts: Analysis and Synthesis,
Fitting, and Timing Analysis. The Analysis and Synthesis will translate the
VHDL code into RTL. The RTL code creates a high-level representation of a
circuit. The Fitter will place and route the design into the FPGA. It will
determine the required resources and the wiring between the different
components. The Timing Analysis will determine the maximum frequency that the
FPGA can reach with a given design.

Quartus can divide a firmware design into multiple partitions. It also
provides the possibility to preserve the results of a given compilation at
different steps of the process for each partition. To further increase the
maximum frequency the firmware design is partitioned in a way that each
partition corresponds to one neural network instance. A sequence of
compilations are performed and each partition that reaches the target
frequency is preserved while the others are recompiled. Four compilations are
needed for a target frequency of 560 MHz and 28 instances of the RNN in the
firmware. 18 out of the 28 instances reach the target frequency at the first
compilation. This allows to run the RNN with a multiplexing of 14 and reach
the required number of channels to be processed in one FPGA.

###  4.4 Results of the VHDL implementation

The final firmware implemented in VHDL contains 28 instances of the vanilla
RNN each with a multiplexing of 14 which allows covering 392 channels. The
results from this firmware are summarised in table 2. This firmware can run at
561 MHz while the required frequency for a multiplexing of 14 is 560 MHz. The
final requirements on the resource usage are not available since the full
final phase-II firmware is not yet available. However, we require a resource
usage of the neural network block similar to the transverse energy computation
block used for the phase-I upgrade [12] of the LAr calorimeter which is
currently in operation. This block uses about 70% of the DSPs and 30% of the
logic. The optimised RNN resource usage is within these requirements as shown
in table 2. The latency of the firmware is 65 clock cycles corresponding to
116 ns which is within the 125 ns required latency.

Table 2: Resource usage and the maximum frequency (Fmax) given by the Quartus compiler report for 28 implemented networks in the FPGA. The multiplexing is set to 14 so that each network can process 14 independent channels. Type | Networks | Multiplexing | Channels | ALM | DSP | Memory | Fmax
---|---|---|---|---|---|---|---
Specification | - | - | 384 | 30% | 70% | 30% | 560 MHz
Quartus | 28 | 14 | 392 | 18% | 66% | 16% | 561 MHz

##  5 Conclusion

The phase-II upgrade of the LAr calorimeter allows the unique opportunity to
implement neural networks in FPGAs in order to improve the computation of the
transverse energy deposited in the calorimeter. This report presents the
implementation of a vanilla RNN in Stratix 10 FPGAs from Intel. The firmware
is first developed in HLS for fast prototyping and optimisation of the network
architecture. The implemented network in firmware matches the computation of
the deposited transverse energy in software with a resolution better than 1‰
after the optimisation of the bit width of the fixed point representation and
the quantisation of the arithmetic operations. However, the HLS implementation
could neither reach the required frequency to allow the implementation of 384
channels per FPGA nor the required latency. VHDL is used to further optimise
the HLS implementation and to add constraints on the placement allowing to
reach the required specifications of the firmware. The final result is a
firmware containing 28 instances of the RNN capable of running at 560 MHz with
a multiplexing of 14 and a latency of 65 clock cycles (116 ns). This firmware
requires less than 70% of the available DSPs and less than 20% of the
available logic elements in the Stratix 10 FPGA.

##  6 Acknowledgements

The project leading to this publication has received funding from Excellence
Initiative of Aix-Marseille Université - A*MIDEX, a French “Investissements
d’Avenir” programme, AMX-18-INT-006 and from the French “Agence National de la
Recherche”, ANR-20-CE31-0013. This work received support from the French
government under the France 2030 investment plan, as part of the Excellence
Initiative of Aix-Marseille University - A*MIDEX (AMX-19-IET-008 - IPhU). This
work would not be possible without the support of the Institut für Kern- und
Teilchenphysik of the Technische Universität Dresden.

## References

  * [1] ATLAS Collaboration, The ATLAS Experiment at the CERN Large Hadron Collider, https://doi.org/10.1088/1748-0221/3/08/S08003JINST 3 (2008) S08003.
  * [2] L. Evans and P. Bryant, LHC machine, https://doi.org/10.1088/1748-0221/3/08/S08001JINST 3 (2008) S08001.
  * [3] ATLAS Collaboration, Technical Design Report for the Phase-II Upgrade of the ATLAS LAr Calorimeter, https://cds.cern.ch/record/2285582CERN-LHCC-2017-018.
  * [4] W.E. Cleland and E.G. Stern, Signal processing considerations for liquid ionization calorimeters in a high rate environment, https://doi.org/10.1016/0168-9002(94)91332-3NIM A 338 (1994) 467-497.
  * [5] G. Aad et al, Artificial Neural Networks on FPGAs for Real-Time Energy Reconstruction of the ATLAS LAr Calorimeters, https://doi.org/10.1007/s41781-021-00066-yComput. Softw. Big. Sci. 5, 19 (2021).
  * [6] Intel Corporation, Intel Agilex FPGA, available on the https://www.intel.com/content/www/us/en/products/details/fpga/agilex.htmlIntel Agilex product website.
  * [7] Intel Corporation, Intel Stratix 10 FPGA, available on the https://www.intel.com/content/www/us/en/products/details/fpga/stratix/10/sx.htmlIntel Stratix 10 product website.
  * [8] J.L. Elman, Finding Structure in Time, https://doi.org/10.1207/s15516709cog1402_1Cognitive Science 14, 2 (1990) 179-211.
  * [9] Intel Corporation, Intel High Level Synthesis Compiler, available on the https://www.intel.com/content/www/us/en/software/programmable/quartus-prime/hls-compiler.htmlIntel HLS compiler website.
  * [10] Mentor Graphics Algorithmic C (AC) data types, provided under the Apache license as part of the https://www.intel.com/content/www/us/en/docs/programmable/683310/19-1/arbitrary-precision-math-support.htmlIntel HLS Compiler.
  * [11] Intel Corporation, Intel Quartus Prime Software, available on the https://www.intel.com/content/www/us/en/products/details/fpga/development-tools/quartus-prime.htmlIntel FPGA Development Tools website.
  * [12] ATLAS Collaboration, ATLAS Liquid Argon Calorimeter Phase-I Upgrade : Technical Design Report, https://cds.cern.ch/record/1602230CERN-LHCC-2013-017.


# Sum of squares decomposition of positive polynomials with rational
coefficients

Santiago Laplagne<EMAIL_ADDRESS>Instituto de Cálculo, FCEN, Universidad
de Buenos Aires - Ciudad Universitaria, Pabellón I - (C1428EGA) - Buenos
Aires, Argentina

###### Abstract.

We present an example of a strictly positive polynomial with rational
coefficients that can be decomposed as a sum of squares of polynomials over
${\mathbb{R}}$ but not over ${\mathbb{Q}}$. This answers an open question by
C. Scheiderer posed as the second question in [14, Section 5.1]. We verify
that the example we construct defines a nonsingular projective hypersurface,
giving also a positive answer to the third question posed in that section.

##  1\. Introduction

Given a polynomial $f\in{\mathbb{R}}[x_{1},\dots,x_{n}]$, deciding if it is
always positive or non-negative is a basic problem in real algebra, with
important applications in optimization and other areas of mathematics. This is
in general a difficult problem, and a usual relaxation is the problem of
deciding if a polynomial can be decomposed as a sum of squares of polynomials,
which provides a non-negativity certificate. This problem can be tackled using
semidefinite programming, for which many numerical algorithms and fast
implementations exist. B. Sturmfels asked the question if every polynomial
with rational coefficients that is a sum of squares of polynomials allows
always a decomposition with rational coefficients. This question is relevant
in the problem of computing exact non-negativity certificates.

In [14], C. Scheiderer provides the first examples of polynomials with
rational coefficients that can be decomposed as a sum of squares with real
coefficients (${\mathbb{R}}$-sos) but cannot be decomposed as a sum of squares
of polynomials with rational coefficients (${\mathbb{Q}}$-sos). The
polynomials obtained in his construction always have common real roots, and
the existence of these common roots is central to prove that there is no
rational decomposition. A different class of examples have been constructed in
[9] and [4]. The strategy used there also produces polynomials with real
roots. A natural question, posed by C. Scheiderer as part of [14, Open
Question 5.1], is whether there exist a polynomial over ${\mathbb{Q}}$ that
can be decomposed as a sum of squares of polynomials over ${\mathbb{R}}$ but
not over ${\mathbb{Q}}$ that is strictly positive, that is, it has no
nontrivial real zeros. This problem is also discussed in [6, Open Question
1.10], where the authors comment on the difficulties related to this problem.
In Section 3 we present an example that gives a positive answer to this
question.

A sharper question posed also by C. Scheiderer is whether there exist such
polynomials defining a nonsingular hypersurface. In Section 4, we show that
our example gives also a positive answer to this question.

##  2\. The cone of sums of squares

We set some notation and recall basic results that will be used in our
construction (see [2, Chapter 3 and 4] for details and proofs). Let $H_{n,d}$
be the vector space of homogeneous polynomials in $n$ variables of degree $d$
with coefficients in ${\mathbb{R}}$. Let $\Sigma_{n,2d}\subset H_{n,2d}$ be
the cone of polynomials that can be written as a sum of squares of polynomials
in $H_{n,d}$ and let $\partial\Sigma_{n,2d}$ be its boundary.

The dual cone $\Sigma_{n,2d}^{*}$ of $\Sigma_{n,2d}$ is the set of all linear
functionals in the dual space $H_{n,2d}^{*}$ that are nonnegative on
$\Sigma_{n,2d}$. A polynomial $f\in\Sigma_{n,2d}$ is in
$\partial\Sigma_{n,2d}$ if and only if there exists $\ell\in\Sigma_{n,2d}^{*}$
such that $\ell(f)=0$. For an element $\ell\in H_{n,2d}^{*}$, we denote
$B_{\ell}$ the associated bilinear form $B_{\ell}(p,q)=\ell(pq)$ for $p,q\in
H_{n,d}$. $B_{\ell}$ is positive semidefinite if and only if
$\ell\in\Sigma_{n,2d}^{*}$. We can also associate to $\ell$ a quadratic form
$Q_{\ell}(p)=\ell(p^{2})$. If $f=p_{1}^{2}+\dots+p_{s}^{2}$ is a sum of
squares and $\ell\in\Sigma_{n,2d}^{*}$, then $\ell(f)=0$ implies
$Q_{\ell}(p_{i})=0$ for all $1\leq i\leq s$. For computations, it is
convenient to fix a monomial basis ${\mathcal{M}}$ of $H_{n,d}$ and represent
$Q_{\ell}$ in the coordinates of ${\mathcal{M}}$ as a matrix
$\mathsfit{Q}_{\ell}\in{\mathbb{R}}^{N\times N}$, where $N=\dim H_{n,d}$.
Since $B_{\ell}$ is positive semidefinite, $Q_{\ell}(p)=0$ implies
$(p)_{{\mathcal{M}}}\in\ker(\mathsfit{Q}_{\ell})$. Moreover
$p\in\ker(Q_{\ell})$ implies $B_{\ell}(p,q)=0$ for all $q\in H_{n,d}$ (see [1,
Lemma 2.6]).

Let $m$ be the vector of all monomials in ${\mathcal{M}}$. A polynomial $f\in
H_{n,2d}$ is a sum of squares if and only if there is a positive semidefinite
matrix $A\in{\mathbb{R}}^{N\times N}$ such that $f=m^{T}Am$. We call any such
matrix a Gram matrix of $f$. If $f$ is in the boundary of the SOS cone then
any Gram matrix of $f$ will be singular, and conversely, if $f$ is in the
interior of the SOS cone then there exists an invertible Gram matrix for $f$.

By the following result, we know that any strictly positive polynomial over
${\mathbb{Q}}$ which is ${\mathbb{R}}$-sos but not ${\mathbb{Q}}$-sos must lie
in the boundary of the SOS cone.

######  Theorem 2.1 (Theorem 1.2, [8]).

Let $f\in{\mathbb{Q}}[x_{1},\dots,x_{n}]$. If there is an invertible Gram
matrix for $f$ with real entries, then there is a Gram matrix for $f$ with
rational entries.

Constructions for strictly positive polynomials in the boundary of the SOS
cone have been found by several authors. In [1], G. Blekherman gives a general
construction for forms of degree 6 in 3 variables and degree 4 in 4 variables.
In [3], the authors extend the latter construction to more general varieties
and in [12], the authors use those ideas to build concrete examples for
polynomials of any degree and number of variables.

##  3\. Main construction

In our construction we glue together two different polynomials. The first
polynomial is a polynomial with rational coefficients that can be decomposed
as a sum of squares of polynomial over ${\mathbb{R}}$ but not over
${\mathbb{Q}}$, which has real zeroes (that is, it is not strictly positive).
The second polynomial is a strictly positive polynomial in the boundary of the
SOS cone.

###  3.1. A polynomial with no rational decomposition

The example we use as one of the blocks was obtained in [4] by J. Capco, S.
Laplagne and C. Scheiderer, simplifying the construction in [9]. We recall the
example here. It is a sum of squares of 3 polynomials in 4 variables with
coefficients in ${\mathbb{Q}}(\sqrt[3]{2})$. Noting $\alpha$ the real root of
$X^{3}-2$, the 3 polynomials are

$\displaystyle p_{1}$
$\displaystyle=(-4\alpha^{2}+4\alpha-2)x_{0}^{2}+x_{1}^{2}-2x_{1}x_{2}+2x_{2}x_{3}-2x_{3}^{2},$
$\displaystyle p_{2}$
$\displaystyle=(-4\alpha^{2}-4\alpha+6)x_{0}^{2}+x_{1}^{2}+2x_{1}x_{2}+2x_{2}x_{3}+2x_{3}^{2},$
$\displaystyle p_{3}$ $\displaystyle=4\alpha
x_{0}x_{1}+4x_{0}x_{2}+4\alpha^{2}x_{0}x_{3},$

and the sum of squares is

$f=p_{1}^{2}+p_{2}^{2}+p_{3}^{2}=40x_{0}^{4}+8x_{0}^{2}x_{1}^{2}+32x_{0}^{2}x_{1}x_{2}+64x_{0}^{2}x_{1}x_{3}+16x_{0}^{2}x_{2}^{2}+16x_{0}^{2}x_{2}x_{3}+32x_{0}^{2}x_{3}^{2}+2x_{1}^{4}+8x_{1}^{2}x_{2}^{2}+8x_{1}^{2}x_{2}x_{3}+16x_{1}x_{2}x_{3}^{2}+8x_{2}^{2}x_{3}^{2}+8x_{3}^{4},$

a polynomial of degree 4 with rational coefficients.

This polynomial is not strictly positive, the real roots are the common real
roots of $p_{1},p_{2},p_{3}$, with approximate affine representatives

$r_{1}=(1,-1.817,0.452,1.158)\text{ and }r_{2}=(0,0,1,0)$

(the coordinates of $r_{1}$ are in an algebraic extension of ${\mathbb{Q}}$ of
degree 12). For our construction, it is only important that there are no real
roots with $x_{2}=0$, which can be verified by hand.

We will now prove that the polynomial $f$ does not allow a rational
decomposition. The strategy used here is different from the original proof by
the authors of [4] and we will use this new proof to extend the example to our
main problem.

######  Proposition 3.1.

Let $f\in{\mathbb{Q}}[x_{0},x_{1},x_{2},x_{3}]$ be as above, then $f$ cannot
be decomposed as a sum of squares with rational coefficients.

See [10, Worksheet A] for the computations in Maple [13]. To prove the claim,
we study first which polynomials can appear in a sum of squares decomposition
of $f$. We seek to compute a linear form $\ell\in\Sigma_{4,4}^{*}$ such that
$\ell(f)=0$. Since $\ell$ is required to be non-negative on squares, the
associated quadratic form

$Q_{\ell}:{\mathbb{R}}[x_{0},x_{1},x_{2},x_{3}]_{2}\rightarrow{\mathbb{R}}$

must be positive semidefinite and, es explained in Section 2, any polynomial
in a sum of squares decomposition of $f$ must be in the kernel of $Q_{\ell}$.

If we can construct a form as above with no non-trivial rational polynomials
in the kernel, then we can assure that $f$ has no rational decomposition. To
achieve this, we will construct a form with kernel as small as possible.

We follow a similar strategy as in [12, Propositions 5.2 and 5.8]. The form
$\ell$ must satisfy $\ell(p_{i}q)=0$ for any $q\in H_{4,2}$. It is defined
over $H_{4,4}$, which is a 35-dimensional space. The dimension of $H_{4,2}$ is
10, so for each $p_{i}$ we get 10 restrictions for $\ell$ . This would be a
total number of $30$ restrictions, but noting that $p_{i}p_{j}=p_{j}p_{i}$, we
are imposing at most $10+9+8=27$ restrictions. Doing the computations in Maple
we obtain that the space $\mathcal{E}$ of linear functionals satisfying all
the above conditions is indeed a $35-27=8$-dimensional space.

We fix a monomial basis

${\mathcal{M}}_{{\bm{x}}}=\\{x_{0}^{2},x_{0}x_{1},x_{0}x_{2},x_{0}x_{3},x_{1}^{2},x_{1}x_{2},x_{1}x_{3},x_{2}^{2},x_{2}x_{3},x_{3}^{2}\\}$

for the 10-dimensional space $H_{4,2}$ and a basis
$\\{\ell_{1},\dots,\ell_{8}\\}$ for $\mathcal{E}$, and we write a general
element of $\mathcal{E}$ as $\ell=t_{1}\ell_{1}+\dots+t_{8}\ell_{8}$,
$t_{i}\in{\mathbb{R}}$, $1\leq i\leq 8$. Let
$\mathsfit{Q}({\bm{t}})=\mathsfit{Q}(t_{1},\dots,t_{8})$, be the matrix
associated to the quadratic form $Q_{\ell}$ for the monomial basis
${\mathcal{M}}_{{\bm{x}}}$ (the entries of $\mathsfit{Q}({\bm{t}})$ are linear
forms in $\\{t_{1},\dots,t_{8}\\}$). We compute values of $t_{1},\dots,t_{8}$
so that the resulting quadratic form is positive semidefinite and has kernel
of low dimension. To achieve this, we compute a numerical approximation using
SEDUMI [15] and deduce from it exact values. SEDUMI uses interior point
methods, so in general the solution will be a solution of maximum rank. The
solution found has rank 4 and kernel of dimension 6. Inspecting the numerical
solution, we fix the following values:

$t_{1}:=0,t_{4}:=0,t_{7}:=0,t_{5}:=1,t_{8}:=1,t_{3}:=1/2.$

After fixing these values we can again check with SEDUMI that there exist
values for the remaining two variables such that the resulting quadratic form
if positive semidefinite. Let $\tilde{\mathsfit}{Q}$ be the new approximation
obtained with SEDUMI. In this case, it is not easy to deduce exact values for
the two remaining unknowns from the approximations given by SEDUMI. Computing
the rank of the principal submatrices of $\tilde{\mathsfit}{Q}$, we observe
that $\tilde{\mathsfit}{Q}_{[2,3]\times[2,3]}$ has rank 1 and
$\tilde{\mathsfit}{Q}_{[5,6,7]\times[5,6,7]}$ has rank 2. Hence we force the
determinants of the submatrices $\mathsfit{Q}({\bm{t}})_{[2,3]\times[2,3]}$
and $\mathsfit{Q}({\bm{t}})_{[5,6,7]\times[5,6,7]}$ to be 0. Solving in Maple
the equations we obtain

$t_{6}:=\beta,t_{2}:=-\frac{1}{2}\alpha-\frac{1}{2}\beta,$

where $\alpha$ is a root of $X^{3}-2$ as before and $\beta$ is a root of
$X^{2}+\alpha^{2}X+(1-\alpha^{2})$ (and we can verify
$\beta\not\in{\mathbb{Q}}(\alpha)$).

Applying these substitutions we obtain a matrix
$\mathsfit{Q}_{{\bm{x}}}\in{\mathbb{Q}}(\alpha,\beta)^{10\times 10}$ and we
can verify by symbolic computations that this matrix is positive semidefinite.
Its kernel has dimension 6, so every polynomial in a sum of squares
decomposition of $f$ is a linear combination of 6 fixed polynomials.

We can group the generators of the kernel of $\mathsfit{Q}_{{\bm{x}}}$ in two
sets supported in complementary sets of monomials. The first set is

$\begin{array}[]{r*{10}{l}}u_{1}&=(0,&0,&1,&-\alpha\beta,&0,&0,&0,&0,&0,&0)_{{\mathcal{M}}_{{\bm{x}}}},\\\
u_{2}&=(0,&1,&0,&\alpha+\beta,&0,&0,&0,&0,&0,&0)_{{\mathcal{M}}_{{\bm{x}}}}\end{array}$

(note that
$p_{3}=(0,4\alpha,4,4\alpha^{2},0,0,0,0,0,0,0)_{{\mathcal{M}}_{{\bm{x}}}}=4o_{1}+4\alpha
o_{2}$).

The second set is

$\begin{array}[]{r*{10}{l}}u_{3}&=(-\alpha+2,&0,&0,&0,&\dfrac{\alpha-4}{4\alpha^{2}-2},&0,&0,&0,&0,&1)_{{\mathcal{M}}_{{\bm{x}}}},\\\
u_{4}&=(-2\alpha^{2}+1,&0,&0,&0,&\dfrac{1}{2},&0,&0,&0,&1,&0)_{{\mathcal{M}}_{{\bm{x}}}},\\\
u_{5}&=(\alpha^{2}+2\alpha+2\beta+2,&0,&0,&0,&\dfrac{\alpha^{2}-4\alpha}{4\alpha^{2}-2},&0,&1,&0,&0,&0)_{{\mathcal{M}}_{{\bm{x}}}},\\\
u_{6}&=(-\alpha+2,&0,&0,&0,&\dfrac{-\alpha+4}{4\alpha^{2}-2},&1,&0,&0,&0,&0)_{{\mathcal{M}}_{{\bm{x}}}}.\end{array}$

We now prove that there exists no non-trivial combination of these vectors
such that the resulting vector has rational coordinates.

######  Lemma 3.2.

With notations as above, let $W=\langle
u_{1},u_{2},u_{3},u_{4},u_{5},u_{6}\rangle_{{\mathbb{R}}}\subset{\mathbb{R}}^{10}$.
Then $W\cap{\mathbb{Q}}^{10}=\\{0\\}$.

###### Proof.

Since $\\{u_{1},u_{2}\\}$ and $\\{u_{3},\dots,u_{6}\\}$ have complementary
non-zero coordinates, we analyze separately the spaces generated by each of
these sets.

We show first that there is no non-trivial linear combination of $u_{1}$ and
$u_{2}$ with rational entries. If $c_{1}u_{1}+c_{2}u_{2}$ has rational
coordinates, by looking at the second and third coordinates, we immediately
see that $c_{1},c_{2}\in{\mathbb{Q}}$. Suppose now
$c_{1}(-\alpha\beta)+c_{2}(\alpha+\beta)=q\in{\mathbb{Q}}$. Then

$\beta=\frac{q-c_{2}\alpha}{c_{2}-c_{1}\alpha},$

but $\beta$ has degree 6 over ${\mathbb{Q}}$ and the right hand size is in
${\mathbb{Q}}(\alpha)$, which is an extension of ${\mathbb{Q}}$ of degree 3
(note that the denominator does not vanish because $c_{1},c_{2}$ cannot be
both zero). Hence, $\\{c_{1}u_{1}+c_{2}u_{2}\\}\cap{\mathbb{Q}}^{10}=\\{0\\}$.

We show now that there is no non-trivial linear combination

$c_{3}u_{3}+c_{4}u_{4}+c_{5}u_{5}+c_{6}u_{6}$

such that the resulting vector has rational coordinates. Otherwise, looking at
the 6th, 7th, 9th and 10th coordinates, we see that $c_{i}\in{\mathbb{Q}}$,
$3\leq i\leq 6$. If we look at the first coordinate, $\beta$ only appears in
$u_{5}$ and $\beta\not\in{\mathbb{Q}}(\alpha)$, hence $c_{5}=0$. Similarly
$c_{4}=0$ because after setting $c_{5}=0$, $\alpha^{2}$ only appears in
$u_{4}$ and $\\{\alpha,\alpha^{2}\\}$ are linearly independent over
${\mathbb{Q}}$. We now look at $c_{3}u_{3}+c_{6}u_{6}$. To cancel $\alpha$ in
the first coordinate, it must be $c_{6}=-c_{3}$, but this gives a non-rational
5th coordinate in $c_{3}u_{3}-c_{3}u_{6}=c_{3}\frac{\alpha-4}{2\alpha^{2}-1}$.
Hence, $\\{c_{3}u_{3}+\dots+c_{6}u_{6}\\}\cap{\mathbb{Q}}^{10}=\\{0\\}$. ∎

This finishes the proof of Proposition 3.1.

######  Remark 3.3.

We observe that the $8$th coordinate (corresponding to $x_{2}^{2}$) is zero in
all the 6 vectors in the kernel. That is, $x_{2}^{2}$ cannot appear in any
polynomial in an SOS decomposition of $f$, which is clear since $x_{2}^{4}$ is
not a monomial in $f$.

###  3.2. A strictly positive polynomial in the boundary of the SOS cone

The second polynomial that we will use as a building block is an example of a
strictly positive polynomial in the boundary of $\Sigma_{4,4}$. The concrete
example we use was obtained in [5], following the construction in [1].

We take

$\displaystyle q_{1}$ $\displaystyle=y_{0}^{2}-y_{3}^{2},$ $\displaystyle
q_{2}$ $\displaystyle=y_{1}^{2}-y_{3}^{2},$ $\displaystyle q_{3}$
$\displaystyle=y_{2}^{2}-y_{3}^{2},$ $\displaystyle q_{4}$
$\displaystyle=-y_{0}^{2}-y_{0}y_{1}-y_{0}y_{2}+y_{0}y_{3}-y_{1}y_{2}+y_{1}y_{3}+y_{2}y_{3},$

and $g=q_{1}^{2}+q_{2}^{2}+q_{3}^{2}+q_{4}^{2}$. Following [5], we know that
the SOS decomposition of $g$ is unique up to orthogonal transformations. In
particular, there is a form $\ell\in\Sigma_{4,4}^{*}$ such that the kernel of
the associated quadratic form $Q_{\ell}$ is $\langle
q_{1},q_{2},q_{3},q_{4}\rangle$. For the monomial basis

${\mathcal{M}}_{{\bm{y}}}=\\{y_{0}^{2},y_{0}y_{1},y_{0}y_{2},y_{0}y_{3},y_{1}^{2},y_{1}y_{2},y_{1}y_{3},y_{2}^{2},y_{2}y_{3},y_{3}^{2}\\},$

the matrix associated to $Q_{\ell}$ (see also [11]) is

$\mathsfit{Q}_{{\bm{y}}}=\begin{pmatrix}6&-1&-1&1&6&-1&1&6&1&6\\\
-1&6&-1&1&-1&-1&1&-1&1&-1\\\ -1&-1&6&1&-1&-1&1&-1&1&-1\\\
1&1&1&6&1&1&-1&1&-1&1\\\ 6&-1&-1&1&6&-1&1&6&1&6\\\
-1&-1&-1&1&-1&6&1&-1&1&-1\\\ 1&1&1&-1&1&1&6&1&-1&1\\\
6&-1&-1&1&6&-1&1&6&1&6\\\ 1&1&1&-1&1&1&-1&1&6&1\\\
6&-1&-1&1&6&-1&1&6&1&6\end{pmatrix}.$

###  3.3. Combination of the two examples.

We define

$r=x_{2}^{2}-y_{2}^{2}$

and $h({\bm{x}},{\bm{y}})=f+g+r^{2}$.

######  Theorem 3.4.

The polynomial

$h(x_{0},x_{1},x_{2},x_{3},y_{0},y_{1},y_{2},y_{3})=f(x_{0},x_{1},x_{2},x_{3})+g(y_{0},y_{1},y_{2},y_{3})+r(x_{2},y_{2})^{2}$

is a strictly positive polynomial with rational coefficients that can be
decomposed as the sum of squares of polynomials in
${\mathbb{Q}}(\sqrt[3]{2})[{\bm{x}},{\bm{y}}]$ but cannot be decomposed as the
sum of squares of polynomials in ${\mathbb{Q}}[{\bm{x}},{\bm{y}}]$.

###### Proof.

We first verify that $h$ is strictly positive. If $g({\bm{x}},{\bm{y}})=0$, by
Section 3.2, it must be $y_{i}=0$ for $0\leq i\leq 3$. Now $r(x_{2},y_{2})=0$
and hence $x_{2}=0$. But $f$ has no non trivial root with $x_{2}=0$.

We now prove that $h$ does not allow a rational decomposition. Our strategy
will be to construct a form $\ell\in\Sigma_{8,4}^{*}$ that vanishes in $h$ and
such that the kernel of the associated quadratic form contains no non-trivial
polynomial with rational coefficients in the ${\bm{x}}$-variables.

We proceed in a similar way as in Section 3.1. The linear form $\ell$ must
satisfy $\ell(pw)=0$ for all
$p\in\\{p_{1},p_{2},p_{3},q_{1},q_{2},q_{3},q_{4},r\\}$ and all $w\in
H_{8,2}$. In this case, the space of forms satisfying these restrictions has
dimension 70. We add as additional restrictions that the associated quadratic
form $Q_{\ell}$ restricted to the ${\bm{x}}$-variables is $Q_{{\bm{x}}}$ and
restricted to the ${\bm{y}}$-variables is $Q_{{\bm{y}}}$. This reduces the
dimension of the affine space of forms to 62.

Setting the monomial basis
${\mathcal{M}}_{{\bm{xy}}}=\\{x_{0}^{2},x_{0}x_{1},x_{0}x_{2},x_{0}x_{3},x_{1}^{2},x_{1}x_{2},x_{1}x_{3},$
$x_{2}^{2}$,
$x_{2}x_{3},x_{3}^{2},y_{0}^{2},y_{0}y_{1},y_{0}y_{2},y_{0}y_{3},y_{1}^{2},y_{1}y_{2},y_{1}y_{3},y_{2}^{2}$,
$y_{2}y_{3}$, $y_{3}^{2}$, $x_{0}y_{0}$, $x_{0}y_{1}$, $x_{0}y_{2}$,
$x_{0}y_{3}$, $x_{1}y_{0}$,
$x_{1}y_{1},x_{1}y_{2},x_{1}y_{3},x_{2}y_{0},x_{2}y_{1},x_{2}y_{2},x_{2}y_{3},x_{3}y_{0},x_{3}y_{1},x_{3}y_{2},x_{3}y_{3}\\}$,
we get a matrix $\mathsfit{Q}$ associated to the quadratic form whose entries
are linear affine in 62 variables. Setting the 62 variables to 0, we get an
exact matrix $\mathsfit{Q}_{{\bm{xy}}}$ and we can verify in Maple that it is
positive semidefinite. The kernel of $\mathsfit{Q}_{{\bm{xy}}}$ has dimension
14. Computing generators in Maple, we see that the kernel can be generated by
the 6 polynomials $\\{u_{1},\dots,u_{6}\\}$ generating the kernel of
$\mathsfit{Q}_{{\bm{x}}}$, plus the polynomials $q_{1},q_{2},q_{3},q_{4},r$
and 3 new polynomials:

$\displaystyle s_{1}$ $\displaystyle=x_{0}y_{0}+x_{0}y_{1},$ $\displaystyle
s_{2}$ $\displaystyle=x_{1}y_{0}+x_{1}y_{1},$ $\displaystyle s_{3}$
$\displaystyle=x_{3}y_{0}+x_{3}y_{1}.$

We conclude that every polynomial in an SOS decomposition of $h$ is a linear
combination of these 14 polynomials.

Assume now that there is a rational SOS decomposition of $h$ and let $p$ be a
polynomial in that decomposition, then $p$ is a linear combination of the 14
generators of the kernel of $\mathsfit{Q}_{{\bm{xy}}}$. We have already seen
that there is no non-trivial polynomial with rational coefficients in the
kernel of $\mathsfit{Q}_{{\bm{x}}}$, and we verify that the monomials
appearing in $q_{1},q_{2},q_{3},q_{4},r,s_{1},s_{2},s_{3}$ are all different
from the monomials in the polynomials $u_{1},\dots,u_{6}$, hence the
coefficients of $u_{1},\dots,u_{6}$ in $p$ must all be zero. This implies that
$p$ has no term $x_{0}^{2}$. Since this is true for all polynomials in a
rational decomposition of $h$, the monomial $x_{0}^{4}$ cannot appear in $h$,
a contradiction since $x_{0}^{4}$ does occur in $h$.

∎

######  Remark 3.5.

If we attempt to construct other examples, numerical experiments using SEDUMI
suggest that if we replace $r=x_{2}^{2}-y_{2}^{2}$ by other polynomials such
as $x_{2}^{2}-y_{1}^{2}$ or $x_{2}^{2}-2y_{2}^{2}$, the resulting polynomials
are also examples of strictly positive polynomials that are ${\mathbb{R}}$-SOS
but not ${\mathbb{Q}}$-SOS. However, if we use instead
$r=x_{2}^{2}+y_{2}^{2}$, the resulting polynomial is in the interior of the
SOS cone, and hence it is a ${\mathbb{Q}}$-SOS.

##  4\. A nonsingular projective variety

In [14, Section 5.1], C. Scheiderer asks if there exists a ${\mathbb{Q}}$-form
which is an ${\mathbb{R}}$-SOS but not a ${\mathbb{Q}}$-SOS such that the
hypersurface it defines is nonsingular. This implies that the polynomial has
no real zeros, because any real zero of a non-negative polynomial is a
singular point of the hypersurface, and that the polynomial is irreducible,
because the intersection of the hypersurfaces defined by the factors is in the
singular locus. Hence, as the author mentions, being nonsingular is a common
sharpening of both properties. That is, the motivation for this question is to
find an example with as little structure as possible, so that it cannot be
detected by a priori methods (e.g. computing the singularities) if the
polynomial can be decomposed over ${\mathbb{Q}}$.

Using Singular [7] (see [10, Worksheet B]), we can verify by a Groebner basis
computation that the hypersurface $\\{h=0\\}$, where $h$ is the polynomial we
defined in the previous section, contains no singular points, that is:

$\sqrt{\left\langle f,\frac{\partial f}{\partial{x_{0}}},\dots,\frac{\partial
f}{\partial{y_{3}}}\right\rangle}=\langle
x_{0},x_{1},x_{2},x_{3},y_{0},y_{1},y_{2},y_{3}\rangle,$

hence giving a positive answer to that question.

## References

  * [1] Grigoriy Blekherman.  Nonnegative polynomials and sums of squares.  In Semidefinite optimization and convex algebraic geometry, volume 13 of MOS-SIAM Ser. Optim., pages 159–202. SIAM, Philadelphia, PA, 2013.
  * [2] Grigoriy Blekherman, Pablo A. Parrilo, and Rekha R. Thomas, editors.  Semidefinite Optimization and Convex Algebraic Geometry.  Society for Industrial and Applied Mathematics, Philadelphia, PA, 2012\.
  * [3] Grigoriy Blekherman, Rainer Sinn, and Mauricio Velasco.  Do sums of squares dream of free resolutions?  SIAM J. Appl. Algebra Geom., 1(1):175–199, 2017.
  * [4] Jose Capco, Santiago Laplagne, and Claus Scheiderer.  Exact polynomial sum of squares decompositions.  Preprint, 2023.
  * [5] Jose Capco and Claus Scheiderer.  Two remarks on sums of squares with rational coefficients.  In Algebra, logic and number theory, volume 121 of Banach Center Publ., pages 25–36. Polish Acad. Sci. Inst. Math., Warsaw, 2020.
  * [6] Lynn Chua, Daniel Plaumann, Rainer Sinn, and Cynthia Vinzant.  Gram spectrahedra.  In Ordered algebraic structures and related topics. International conference at CIRM, Luminy, France, October 12–16, 2015. Proceedings, pages 81–105. Providence, RI: American Mathematical Society (AMS), 2017.
  * [7] Wolfram Decker, Gert-Martin Greuel, G. Pfister, and H. Schönemann.  Singular 4.3.0 — A computer algebra system for polynomial computations.  http://www.singular.uni-kl.de, 2022.
  * [8] Christopher J. Hillar.  Sums of squares over totally real fields are rational sums of squares.  Proc. Amer. Math. Soc., 137(3):921–930, 2009.
  * [9] Santiago Laplagne.  Facial reduction for exact polynomial sum of squares decomposition.  Math. Comp., 89(322):859–877, 2020.
  * [10] Santiago Laplagne.  Maple and singular scripts code for this paper.  https://github.com/slap/sums-of-squares, 2023.
  * [11] Santiago Laplagne.  On the pythagoras number for polynomials of degree 4 in 5 variables.  Revista de la Unión Matemática Argentina, 2023.  In press.
  * [12] Santiago Laplagne and Marcelo Valdettaro.  Strictly positive polynomials in the boundary of the sos cone, 2023.  arXiv:2012.05951v2.
  * [13] Maple.  Version 2015.1.  Maple is a trademark of Waterloo Maple Inc., 2015.
  * [14] Claus Scheiderer.  Sums of squares of polynomials with rational coefficients.  J. Eur. Math. Soc. (JEMS), 18(7):1495–1513, 2016.
  * [15] Jos F. Sturm.  Using sedumi 1.02, a matlab toolbox for optimization over symmetric cones.  Optimization Methods and Software, 11(1-4):625–653, 1999.