#  A new iterative method for construction of the Kolmogorov-Arnold
representation

Michael Poluektov  Department of Mathematics, School of Science and
Engineering, University of Dundee, Dundee DD1 4HN, UK  Corresponding author,
email<EMAIL_ADDRESS>Andrew Polar  Independent software
consultant, Duluth, GA, USA

( DRAFT: August 28, 2024)

###### Abstract

The Kolmogorov-Arnold representation of a continuous multivariate function is
a decomposition of the function into a structure of inner and outer functions
of a single variable. It can be a convenient tool for tasks where it is
required to obtain a predictive model that maps some vector input of a black
box system into a scalar output. However, the construction of such
representation based on the recorded input-output data is a challenging task.
In the present paper, it is suggested to decompose the underlying functions of
the representation into continuous basis functions and parameters. A novel
lightweight algorithm for parameter identification is then proposed. The
algorithm is based on the Newton-Kaczmarz method for solving non-linear
systems of equations and is locally convergent. Numerical examples show that
it is more robust with respect to the section of the initial guess for the
parameters than the straightforward application of the Gauss-Newton method for
parameter identification.
Keywords: discrete Urysohn operator, generalised additive model,
Kolmogorov-Arnold representation, ridge function, model identification,
Kaczmarz method, Newton-Kaczmarz method.

## 1 Introduction

One of the typical data science tasks is constructing a predictive model that
maps some vector input into a scalar output. Such model can be viewed as a
multivariate function. In 1950s, Andrei Kolmogorov and Vladimir Arnold showed
that any continuous multivariate function can be represented by a composition
of functions of a single variable [1, 2]. Thus, the Kolmogorov-Arnold
representation is convenient tool for data modelling tasks, replacing the
problem of construction of a multivariate function in its entirety by building
a set of functions of a single variable. The representation has its use in
practical applications, for example, in the field of image processing [3, 4].

The original Kolmogorov’s proof was not constructive and did not provide
insights into choosing the underlying functions of the representation. There
is a series of works that gives constructive proofs of the Kolmogorov’s
theorem and proposes convergent algorithms for determining the functions [5,
6, 7, 8, 9, 10]. There are also papers, where various approximations of the
underlying functions and algorithms for their determination have been
suggested [11, 12, 13]. Recently, a lightweight iterative algorithm for
construction of the underlying functions has been proposed by the authors of
this paper [14]; however, the functions were limited to the class of
piecewise-linear functions and the algorithm was given without the convergence
proof.

The aim of the present paper is two-fold: (1) to propose a more general
decomposition of the underlying functions of the Kolmogorov-Arnold
representation via continuous basis functions and parameters; and (2) to
propose a novel lightweight algorithm, which is proven to be locally
convergent, for determining the parameters based on the input-output data.

As has been pointed out in [14], the Kolmogorov-Arnold representation can be
viewed as a tree of objects called the discrete Urysohn operators. This
property is useful, as the algorithms for construction of the discrete Urysohn
operator and the Kolmogorov-Arnold representation can have common elements.
Therefore, section 2 provides an overview of the Urysohn and the Kolmogorov-
Arnold models. Section 3 presents the novel material. First, the decomposition
of the underlying functions of the Urysohn model into continuous basis
functions and parameters is introduced. Then, the previously-proposed by the
authors of this paper identification algorithm for the Urysohn model [15] is
rewritten for such decomposition. Finally, the novel identification algorithm
for the Kolmogorov-Arnold model is proposed.

The algorithm is tested in section 4. It has already been shown in [14] that
the Kolmogorov-Arnold model in application to data modelling tasks performs
equivalently to neural networks in terms of accuracy, but can require less
computational time and resources for model identification. Therefore, the
numerical examples do not aim to compare it to other models, such as neural-
network-based ones, but rather aim to demonstrate the advantages of the
proposed identification technique.

## 2 Models

### 2.1 The Urysohn model

The Urysohn operator has originally been introduced within the theory of non-
linear integral equations [16]. It remained largely a theoretical concept
until 1970s, when it was picked up by researchers from the theory of automatic
control and was reintroduced as a continuous-time integral operator that
transforms function $x\left(t\right)$ to function $z\left(t\right)$ in the
following way [17]:

$z\left(t\right)=\int_{0}^{T}V\left(s,x\left(t-s\right)\right)\mathrm{d}s,$
(1)

where
$x:\left[-T,+\infty\right)\to\left[x_{\mathrm{min}},x_{\mathrm{max}}\right]$,
$z:\left[0,+\infty\right)\to\mathbb{R}$,
$V:\left[0,T\right]\times\left[x_{\mathrm{min}},x_{\mathrm{max}}\right]\to\mathbb{R}$
are continuous almost everywhere functions, $t\geq 0$, $T\geq 0$ and
$x_{\mathrm{min}},x_{\mathrm{max}}\in\mathbb{R}$. Function $V\left(s,x\right)$
is the kernel of the continuous-time Urysohn operator.

Since digital equipment collects continuous signals in a form of numeric
series, the discrete-time Urysohn operator has also been introduced [17]:

$z_{i}=\sum_{j=1}^{m}g^{j}\left(x_{i-j+1}\right),\quad i\in\mathbb{N},$ (2)

where $x_{i}\in\left[x_{\mathrm{min}},x_{\mathrm{max}}\right]$ is the series
of input values, $z_{i}\in\mathbb{R}$ is the series of output values, $m$ is
the memory depth of the operator and
$g^{j}:\left[x_{\mathrm{min}},x_{\mathrm{max}}\right]\to\mathbb{R}$ are
continuous almost everywhere functions. It can be seen that equation (2)
results from a numerical quadrature of equation (1). In the present paper,
superscripts are attached to quantities (e.g. functions or vectors) to denote
their numbers, while a quantity with subscript is a single scalar (e.g. an
element of a vector).

An interesting property of the discrete-time form of the Urysohn operator is
that it can be used not only for control systems (i.e. mapping series $x_{i}$
to series $z_{i}$), but also for a general task of data modelling. Given set
of data records $\left(X^{i},y_{i}\right)$, $i\in\left\\{1,\ldots,N\right\\}$,
where $X^{i}\in\mathbb{R}^{m}$ is the input of the $i$-th record,
$y_{i}\in\mathbb{R}$ is the output of the $i$-th record and $N$ is the number
of records, operator $\mathcal{U}:\mathbb{R}^{m}\to\mathbb{R}$ can be
introduced as [14]:

$\mathcal{U}\left(X^{i}\right)\equiv\sum_{j=1}^{m}g^{j}\left(X^{i}_{j}\right),$
(3)

where $X^{i}_{j}$ denote the $j$-th component of vector $X^{i}$. For
consistency of the notation, it is assumed that
$X^{i}_{j}\in\left[x_{\mathrm{min}},x_{\mathrm{max}}\right]$. The operator
notation can now be used to formulate the Urysohn model:

$\hat{y}_{i}=\mathcal{U}\left(X^{i}\right),$ (4)

where $\hat{y}_{i}$ is the calculated model output of the $i$-th record. It is
important to emphasise that the distinction is made between the actual
(recorded) output $y_{i}$ and the model output $\hat{y}_{i}$.

It is useful to highlight that the discrete Urysohn operator written for the
case of data modelling, equation (3), is mathematically equivalent to the so-
called generalised additive model (GAM). The latter has been introduced and
researched independently [18]. One possible reason for this is that the
discrete Urysohn operator has been introduced originally for dynamical
systems, equation (2), and most research has been done within the fields of
automatic control and signal processing. In the original definition of the
discrete Urysohn operator, the inputs are the discrete values of the same
signal; therefore, all functions are defined on the same domain and the
adjacent functions are expected to be similar. In the data modelling case,
equation (3), the functions can be defined on different domains and each
argument can be continuous or quantised.

### 2.2 The Kolmogorov-Arnold model

In [1, 2], it has been shown that any continuous multivariate function can be
represented by a composition of functions of a single variable. More
precisely, function $F:\mathbb{R}^{m}\to\mathbb{R}\in
C\left(\left[0,1\right]^{m}\right)$ can be represented as

$F\left(X_{1},X_{2},\ldots,X_{m}\right)=\sum_{k=1}^{2m+1}\varPhi^{k}\left(\sum_{j=1}^{m}f^{kj}\left(X_{j}\right)\right),$
(5)

where $f^{kj}:\mathbb{R}\to\mathbb{R}\in C\left[0,1\right]$ and
$\varPhi^{k}:\mathbb{R}\to\mathbb{R}\in C\left(\mathbb{R}\right)$. This
decomposition can be used for a general task of data modelling, where the
output is a continuous function of the inputs. Using the notation of the
discrete Urysohn operator and given set of data records
$\left(X^{i},y_{i}\right)$, intermediate variable
$\theta^{i}\in\mathbb{R}^{2m+1}$ is introduced for each data record and the
Kolmogorov-Arnold model for the input-output relationship [14] is then written
as

$\hat{y}_{i}=\mathcal{U}^{0}\left(\theta^{i}\right),\quad\theta^{i}_{k}=\mathcal{U}^{k}\left(X^{i}\right),\quad
k\in\left\\{1,\ldots,2m+1\right\\},$ (6)

where operator $\mathcal{U}^{0}$ contains functions $\varPhi^{k}$ and
operators $\mathcal{U}^{k}$ contain functions $f^{kj}$. Here, $\theta^{i}_{k}$
denote the $k$-th component of vector $\theta^{i}$. For consistency of the
notation, the domains of functions $f^{kj}$ have been extended to
$\left[x_{\mathrm{min}},x_{\mathrm{max}}\right]$ in equation (6), and the
functions are taken as continuous almost everywhere. This model can also be
interpreted as a tree of the discrete Urysohn operators with a particular
hierarchy, containing $1$ ‘root’ operator with $2m+1$ functions and $2m+1$
‘branch’ operators with $m$ functions each. Intermediate variable $\theta^{i}$
can then be interpreted as a hidden layer in the model.

It is useful to note that the Kolmogorov-Arnold representation can also be
viewed as a generalisation of the so-called ridge function [19]:

$F\left(X_{1},X_{2},\ldots,X_{m}\right)=\varPhi\left(\sum_{j=1}^{m}c_{j}X_{j}\right),$
(7)

where $c_{j}\in\mathbb{R}$ are some constants and
$\varPhi:\mathbb{R}\to\mathbb{R}$ is some function. The ridge function models
are used in some areas of data modelling and are even suggested as an
approximate replacement of the Kolmogorov-Arnold model and neural networks
[19].

## 3 Discretisation and identification

### 3.1 Building the Urysohn model

To handle the models computationally, the underlying functions must be
represented in some discretised form. For the Urysohn model, one option is to
use the following decomposition of functions $g^{j}$:

$g^{j}\left(x\right)=\sum_{l=1}^{n}U_{jl}\phi^{l}\left(x\right),$ (8)

where $U_{jl}$ are the model parameters and
$\phi^{l}:\left[x_{\mathrm{min}},x_{\mathrm{max}}\right]\to\mathbb{R}$ are the
basis functions with compact support. From equations (3) and (8), it can be
seen that the discrete Urysohn operator is defined by its parameters $U_{jl}$,
which can also be understood as elements of $m$-by-$n$ matrix $U$.

The identification problem for the Urysohn model is finding parameters
$U_{jl}$ given set of records $\left(X^{i},y_{i}\right)$, such that the
discrepancy between $\hat{y}_{i}$ and $y_{i}$ is minimised. This problem and
its properties (e.g. the non-uniqueness of the solution) have been studied in
detail in [15] for piecewise-linear and piecewise-constant basis functions.

It can be seen that $\hat{y}_{i}$ is a linear function of the parameters (but
not a linear function of the inputs). Therefore, it is possible to assemble a
system of linear equations with respect to $U_{jl}$ and use it to reformulate
the identification problem. Thus, the model parameters are rearranged into
vector $Z$, the values of the basis functions are rearranged into matrix $M$,
rows of which are denoted as $M^{i}$, and the outputs are rearranged into
vector $Y$:

$\displaystyle
Z=\begin{bmatrix}U_{11}&\ldots&U_{1n}&U_{21}&\ldots&U_{2n}&\ldots&U_{m1}&\ldots&U_{mn}\end{bmatrix}^{\mathrm{T}},$
$\displaystyle
M^{i}=\begin{bmatrix}\phi^{1}\left(X^{i}_{1}\right)&\ldots&\phi^{n}\left(X^{i}_{1}\right)&\phi^{1}\left(X^{i}_{2}\right)&\ldots&\phi^{n}\left(X^{i}_{2}\right)&\ldots&\phi^{1}\left(X^{i}_{m}\right)&\ldots&\phi^{n}\left(X^{i}_{m}\right)\end{bmatrix},$
$\displaystyle Y=\begin{bmatrix}y_{1}&\ldots&y_{N}\end{bmatrix}^{\mathrm{T}}.$

For convenience, elements of vector $Y$ are denoted as $Y_{i}=y_{i}$. This
transforms the identification problem into the following minimisation problem:

$Z^{*}=\underset{Z}{\operatorname{argmin}}\left\|Y-MZ\right\|^{2},$ (9)

where $Z^{*}$ is the sought solution. In [15], it has been suggested to find
an approximate solution of this problem iteratively row-by-row (or record-by-
record) using the projection descend method [20, 21], also known as the
Kaczmarz method, which represents the following sequence:

$Z^{q+1}=Z^{q}+\mu\frac{Y_{i}-M^{i}Z^{q}}{\left\|M^{i}\right\|^{2}}{M^{i}}^{\mathrm{T}},$
(10)

where $\mu\in\left(0,2\right)$ is the regularisation parameter and $Z^{q}$ is
the approximation of the solution at iteration $q$. Index $i$ changes each
iteration. The method requires initial guess $Z^{0}$ and continues until some
convergence criteria are reached. It should be noted that in [15], the method
is presented in less general terms — it is written specifically for piecewise-
linear and piecewise-constant basis functions.

A number of important remarks should be made briefly regarding the method
above (for details see [15] and references therein). First, the method
considers system of linear equations $MZ=Y$. For an underdetermined system,
the method converges to a solution that is closest to the initial guess in
terms of the $L^{2}$-norm. For an overdetermined system, it produces a
solution that oscillates in a region around the true solution of minimisation
problem (9). The size of this region decreases with decreasing $\mu$. Second,
although one gets an approximate solution, the main advantage of the method is
that it does not require assembling matrix $M$ and keeping it in the memory,
as the method iterates row-by-row. Furthermore, in the case of the Uryshon
model, since functions $\phi^{l}$ have compact support, matrix $M$ is sparse;
therefore, only a small subset of elements of $Z^{q}$ is updated in equation
(10), which is computationally cheap. Third, a consequence of the sparsity of
$M$ is that the projection descend method converges relatively fast, as the
rows of $M$ are either orthogonal or close to orthogonal.

### 3.2 Building the Kolmogorov-Arnold model

As in the case of the Urysohn model, the underlying functions of the model are
decomposed as

$\displaystyle
f^{kj}\left(x\right)=\sum_{p=1}^{n}H_{kjp}\phi^{p}\left(x\right),$ (11)
$\displaystyle\varPhi^{k}\left(t\right)=\sum_{l=1}^{s}G_{kl}\psi^{l}\left(t\right),$
(12)

where $H_{kjp}$ and $G_{kl}$ are the model parameters,
$\phi^{p}:\left[x_{\mathrm{min}},x_{\mathrm{max}}\right]\to\mathbb{R}$ and
$\psi^{l}:\left[t_{\mathrm{min}},t_{\mathrm{max}}\right]\to\mathbb{R}$ are the
basis functions with compact support and
$t_{\mathrm{min}},t_{\mathrm{max}}\in\mathbb{R}$. Functions $\phi^{p}$ and
$\psi^{l}$ must be continuous, and functions $\psi^{l}$ must be differentiable
almost everywhere. Here, $n$ and $s$ are the numbers of the basis functions
used for the inner and the outer functions of the model, respectively.

The identification problem for the Kolmogorov-Arnold model is finding
parameters $H_{kjp}$ and $G_{kl}$ given set of records
$\left(X^{i},y_{i}\right)$, such that the discrepancy between $\hat{y}_{i}$
and $y_{i}$ is minimised. As above, this problem can be written as
minimisation problem

$\left\\{H_{kjp}^{*},G_{kl}^{*}\right\\}=\underset{H,G}{\operatorname{argmin}}\sum_{i=1}^{N}\left(y_{i}-\hat{y}_{i}\right)^{2},$
(13)

where $H_{kjp}^{*}$ and $G_{kl}^{*}$ are the sought parameters. In contrast to
the Urysohn model, output $\hat{y}_{i}$ is now a non-linear function of the
model parameters.

An approximate solution of problem (13) can be found using the Newton-Kaczmarz
method [22, 23, 24] for solving non-linear systems of equations. As outlined
in appendix A, this method can be understood as an iterative sequence of
alternating between equations’ linearisation and application of one step of
the Kaczmarz method. To make the paper self-contained, the steps of obtaining
the iterative scheme are outlined below.

First, the system of non-linear equations with respect to the model parameters
is written as

$\hat{y}_{i}-y_{i}=0,\quad\quad i\in\left\\{1,\ldots,N\right\\},$ (14)

which is expanded by substituting equations (11) and (12) into (5), leading to

$\sum_{k=1}^{2m+1}\sum_{l=1}^{s}G_{kl}\psi^{l}\left(\sum_{j=1}^{m}\sum_{p=1}^{n}H_{kjp}\phi^{p}\left(X^{i}_{j}\right)\right)-y_{i}=0,\quad\quad
i\in\left\\{1,\ldots,N\right\\}.$ (15)

To find the solution of system (15) iteratively, it is necessary to find the
expressions for the parameters at iteration $q+1$, knowing the parameters at
iteration $q$. The approximations of the model parameters at iteration $q$ are
denoted as $H_{kjp}^{q}$ and $G_{kl}^{q}$. It is assumed that initial guess
$H_{kjp}^{0}$ and $G_{kl}^{0}$ for the parameters is made. The increments of
the model parameters from iteration $q$ to iteration $q+1$ are denoted as

$\Delta H_{kjp}=H_{kjp}^{q+1}-H_{kjp}^{q},\quad\quad\quad\Delta
G_{kl}=G_{kl}^{q+1}-G_{kl}^{q}.$ (16)

The arguments of functions $\psi^{l}$ in system (15) are the components of
intermediate variable $\theta^{i}$, introduced in section 2.2. Hence, at
iteration $q$,

$\theta_{k}^{i}=\sum_{j=1}^{m}\sum_{p=1}^{n}H_{kjp}^{q}\phi^{p}\left(X^{i}_{j}\right).$
(17)

For brevity, it is useful to denote the entire first term of equation (15) at
iteration $q$ as

$E=\sum_{k=1}^{2m+1}\sum_{l=1}^{s}G_{kl}^{q}A_{kl},\quad\quad\quad
A_{kl}=\psi^{l}\left(\theta_{k}^{i}\right).$ (18)

Next, it is necessary to introduce the derivatives of basis functions
$\psi^{l}$,

$\eta^{l}\left(t\right)=\frac{\mathrm{d}\psi^{l}\left(t\right)}{\mathrm{d}t}.$
(19)

Now, system (15) is taken at iteration $q+1$ and is linearised with respect to
increments $\Delta H_{kjp}$ and $\Delta G_{kl}$ of the model parameters,
leading to

$E+\sum_{k=1}^{2m+1}\sum_{l=1}^{s}A_{kl}\Delta
G_{kl}+\sum_{k=1}^{2m+1}\sum_{j=1}^{m}\sum_{p=1}^{n}B_{kjp}\Delta
H_{kjp}+O\left(\left\|\Delta G\right\|^{2}+\left\|\Delta
H\right\|^{2}\right)-y_{i}=0,$ (20)

where additionally-introduced quantities are defined as

$B_{kjp}=\phi^{p}\left(X^{i}_{j}\right)\sum_{l=1}^{s}G_{kl}^{q}C_{kl},\quad\quad\quad
C_{kl}=\eta^{l}\left(\theta_{k}^{i}\right).$ (21)

For quantities $A_{kl}$, $B_{kjp}$, $C_{kl}$ and $E$, superscript/subscript
$i$ is omitted for simplicity, although the quantities depend on $X^{i}$,
hence, change with index $i$. Neglecting the second-order terms, gives a
linear system with respect to $\Delta H_{kjp}$ and $\Delta G_{kl}$. This
allows writing one step of the Kaczmarz method to find $\Delta H_{kjp}$ and
$\Delta G_{kl}$, taking the previous-step values as zeros, leading to

$\displaystyle
H_{kjp}^{q+1}=H_{kjp}^{q}+\mu\frac{y_{i}-E}{\zeta}B_{kjp},\quad\quad\quad
G_{kl}^{q+1}=G_{kl}^{q}+\mu\frac{y_{i}-E}{\zeta}A_{kl},$ (22)
$\displaystyle\zeta=\sum_{k=1}^{2m+1}\sum_{l=1}^{s}A_{kl}^{2}+\sum_{k=1}^{2m+1}\sum_{j=1}^{m}\sum_{p=1}^{n}B_{kjp}^{2},$
(23)

where $\mu\in\left(0,2\right)$ is the regularisation parameter. Sequence (22)
represents the iterative scheme for finding the model parameters, given
initial approximations $H_{kjp}^{0}$ and $G_{kl}^{0}$. As previously, index
$i$ changes each iteration. The iterations are performed until some
convergence criteria are reached, e.g. $\left|\hat{y}_{i}-y_{i}\right|$ is
sufficiently small for sufficiently large number of iterations (the exact
thresholds are decided by the user). As for the case of the Urysohn model,
system (14) can be underdetermined or overdetermined. An example of the latter
is a system with noise-affected recorded output $y_{i}$. As in section 3.1,
the method then gives an approximate solution, which oscillates in a region
around the true solution of minimisation problem (13) with the size of this
region depending on $\mu$.

It is important to note that the original Newton-Kaczmarz method is formulated
the systems of non-linear equations where functions are continuously
differentiable, as outlined in appendix A. Thus, if functions $\psi^{l}$ are
continuously differentiable, then the algorithm is locally convergent, as it
is the direct application of the Newton-Kaczmarz method to system (14), which
makes it distinct from the algorithm proposed previously in [14] that was
given without the convergence proof.

It should also be noted that in this section, the continuity requirements for
$\psi^{l}$ have been relaxed. As has been found out by the extensive testing
of the algorithm, in practice, this does not noticeably impact the convergence
of the method, although care must be taken when $\theta_{k}^{i}$ falls exactly
at a point of discontinuity of $\eta^{l}$. The relaxation of the continuity
requirements allows using the piecewise-linear basis functions, which is
convenient from a practical point of view.

#### 3.2.1 Piecewise-linear basis functions

In practice, it is convenient to use the piecewise-linear basis functions. To
define the functions, equally-spaced nodes are introduced as

$x_{p}=x_{\mathrm{min}}+\frac{p-1}{n-1}\left(x_{\mathrm{max}}-x_{\mathrm{min}}\right),\quad
p\in\left\\{1,\ldots,n\right\\}.$ (24)

The same is done for the intermediate variable:

$t_{l}=t_{\mathrm{min}}+\frac{l-1}{s-1}\left(t_{\mathrm{max}}-t_{\mathrm{min}}\right),\quad
l\in\left\\{1,\ldots,s\right\\}.$ (25)

As introduced previously, $n$ and $s$ are the numbers of the inner and the
outer basis functions, respectively. For convenience, the nodes are combined
in sets $P=\left\\{x_{p}\right\\}$ and $Q=\left\\{t_{l}\right\\}$. The basis
functions are continuous, are linear between the nodes, are exactly $1$ at a
given node and are exactly $0$ at all other nodes:

$\phi^{p}\left(x\right)=\begin{cases}0,&\text{if }x\in P\text{ and }x\neq
x_{p},\\\ 1,&\text{if }x=x_{p},\\\
\text{linear},&\text{otherwise},\end{cases}\quad\quad\quad\psi^{l}\left(t\right)=\begin{cases}0,&\text{if
}t\in Q\text{ and }t\neq t_{l},\\\ 1,&\text{if }t=t_{l},\\\
\text{linear},&\text{otherwise}.\end{cases}$ (26)

#### 3.2.2 Initial guess for the parameters

In general, initial approximations $H_{kjp}^{0}$ and $G_{kl}^{0}$ can be
arbitrary. However, the outputs of the inner functions constitute the
intermediate variable, the components of which are discretised according to
equation (25). With ‘bad’ initial approximations of the parameters, the
intermediate variable might significantly drift during the identification
process. In the previous method [14], it was suggested to update limits
$t_{\mathrm{min}}$ and $t_{\mathrm{max}}$ of the intermediate variable and its
nodal positions $t_{l}$. A better approach is to initialise the parameters in
such a way that this updating procedure is not necessary.

The extensive testing of the method with the piecewise-linear basis functions
showed that in practice, $t_{\mathrm{min}}=y_{\mathrm{min}}$ and
$t_{\mathrm{max}}=y_{\mathrm{max}}$ can be taken, where $y_{\mathrm{min}}$ and
$y_{\mathrm{max}}$ are the minimum and the maximum values of output $y_{i}$,
respectively. The model parameters are then initialised as uniformly-
distributed random numbers:

$H_{kjp}^{0}\sim\operatorname{unif}\left(\frac{y_{\mathrm{min}}}{m},\frac{y_{\mathrm{max}}}{m}\right),\quad\quad
G_{kl}^{0}\sim\operatorname{unif}\left(\frac{y_{\mathrm{min}}}{2m+1},\frac{y_{\mathrm{max}}}{2m+1}\right),$
(27)

where $2m+1$ in the denominator corresponds to the number of addends in the
model111If the number of addends is decreased, then the denominator in
equation (27) should be modified accordingly..

#### 3.2.3 Prevention of overfitting

The common way of the overfitting detection is the usage of the training and
the validation datasets. The validation data is not used for the training and
is typically referred to as the unseen data. The complexity of the Kolmogorov-
Arnold model is determined by the numbers of the basis functions, $n$ and $s$.
By modifying these parameters, fitting the model to the training dataset and
verifying the model’s performance using the validation dataset, the optimal
complexity can be established.

Furthermore, as was shown in [14], it is not necessary to use the ‘full’ model
all the time. The ‘full’ Kolmogorov-Arnold model consists of $2m+1$ addends
(corresponding to functions $\varPhi^{k}$), as seen from equation (5).
Depending on the complexity of the modelled system and the desired accuracy,
the number of addends can be decreased, which gives the third parameter
responsible for the complexity of the model that should be fine-tuned using
the training and the validation datasets.

## 4 Numerical examples

In this section, two numerical examples are presented to show the efficiency
of the proposed solution technique. After having represented unknown functions
$f^{kj}$ and $\varPhi^{k}$ via basis functions $\phi^{p}$ and $\psi^{l}$ and
parameters $H_{kjp}$ and $G_{kl}$, various methods can be used to determine
the parameters. The straightforward choice is the Gauss-Newton (GN) method.
Therefore, it is necessary to investigate how it performs compared to the
Newton-Kaczmarz (NK) method for the considered identification problem, which
is done in the first example.

### 4.1 Example 1: Identification of a ridge function

To compare the GN and the NK methods, the ridge function model, equation (7),
is taken. It has similar structure to the Kolmogorov-Arnold model, but with
one addend and with inner functions being linear. Reasons for this choice will
be clarified later.

Since Newton-type methods require an initial guess, it is important that the
identification technique is robust with respect to the choice of the guess and
does not get stuck at local minima. Therefore, this example aims to compare
the dependence of the methods’ convergence on the initial guess.

The example consists in assuming a certain model, generating the data (i.e.
generating the inputs and calculating the corresponding outputs using the
assumed model) and then trying to identify the model parameters using the
data. The advantage of such approach is the availability of the exact solution
of the identification problem. The code of this example is available at the
GitHub repository222https://github.com/andrewpolar/RidgeIdentM.

As before, outer function $\varPhi$ is represented via the basis functions and
the parameters:

$\varPhi\left(t\right)=\sum_{l=1}^{s}G_{l}\psi^{l}\left(t\right).$ (28)

To avoid any possible convergence issues, smooth basis functions are taken:

$\psi^{l}\left(t\right)=\exp\left(-2\left(t-t_{l}\right)^{2}\right),$ (29)

where $t_{l}$ are the nodal positions.

To generate the data, a model with $m=5$ inputs is taken. The following
parameters are assumed:

$c=\begin{bmatrix}-0.7&2.5&-1.2&0.8&1.6\end{bmatrix},\quad\quad
G=\begin{bmatrix}2.1&-0.9&0.7\end{bmatrix},\quad\quad t_{l}=l-0.5,$ (30)

where $l\in\left\\{1,2,3\right\\}$. To create the datasets of $N=400$ records,
uniformly-distributed random inputs from $x_{\mathrm{min}}=0$ to
$x_{\mathrm{max}}=1$ are generated,
$X^{i}_{j}\sim\operatorname{unif}\left(0,1\right)$, and outputs $y_{i}$ are
calculated.

Having the datasets, the identification problem is solved using the NK and the
GN methods. The latter method is summarised in appendix B. In the NK method,
regularisation parameter $\mu=0.1$ is used and $10^{4}$ iterations are
performed. In the GN method, tolerance $\delta=10^{-12}$ is set and the
maximum number of iterations is taken to be $100$.

To solve the identification problem, the initial guess for the parameters is
taken as the exact solution plus a random perturbation, i.e.

$c_{j}^{0}=c_{j}+\alpha e_{j}^{1},\quad\quad G_{l}^{0}=G_{l}+\alpha
e_{l}^{2},$ (31)

where $e_{j}^{1},e_{l}^{2}\sim\operatorname{unif}\left(-0.5,0.5\right)$ are
uniformly-distributed random perturbations. For small $\alpha$, the initial
guess is close to the true solution and the methods are expected to converge,
while for large $\alpha$, the methods might not converge.

The normalised root-mean-square error (RMSE),

$E_{\mathrm{RMSE}}=\frac{1}{y_{\mathrm{max}}-y_{\mathrm{min}}}\sqrt{\frac{1}{N}\sum_{i=1}^{N}\left(y_{i}-\hat{y}_{i}\right)^{2}},$
(32)

is used as a measure to compare the methods. After the methods solve the
identification problem and produce the model parameters for given $\alpha$,
error $E_{\mathrm{RMSE}}$ is calculated. Since the generation of the datasets
and the initial guess for the parameters involve randomisation, error
$E_{\mathrm{RMSE}}$ can significantly differ from one run333One execution of
the code is implied by the ‘run’. to another; therefore, its average behaviour
should be considered. In particular, an illustrative characteristic of the
performance of the methods is the average proportion of runs, in which
$E_{\mathrm{RMSE}}$ is below some threshold. Such characteristic is more
useful than a direct comparison of average $E_{\mathrm{RMSE}}$ produced by the
GN and the NK methods because it is not sensitive to the outliers444As
observed during the testing, in some rare runs, parameter values diverge
significantly, leading to a very large error..

The average percentages of runs where the RMSE was below $5\%$, $10\%$ and
$20\%$ for different amplitude of the perturbations of the initial guess,
parameter $\alpha$, are shown in table 1 for both methods. As expected, as
$\alpha$ increases, the percentage of accurate runs drops. More importantly,
it can be seen that the GN method consistently results in fewer runs with the
error below a given threshold than the NK method. Furthermore, the accuracy of
the GN method drops more rapidly than of the NK method. For example, for
$\alpha=1.2$, more than $90\%$ of the runs of the NK method produced errors
less than $20\%$, while only about a half of the runs of the GN method had
such accuracy; for $\alpha=2.8$, which is close to the completely random
initial guess, since such value of $\alpha$ is greater than any parameter
value in the assumed model, about $20\%$ of the runs of the NK method produced
errors less than $10\%$, while only about $1\%$ of the runs of the GN method
had such accuracy.

$\alpha$ | $0.4$ | $0.8$ | $1.2$ | $1.6$ | $2.0$ | $2.4$ | $2.8$
---|---|---|---|---|---|---|---
Gauss-Newton (GN) method
err. $<5\%$ | $93.6\pm 1.3$ | $59.0\pm 5.2$ | $31.4\pm 1.7$ | $12.4\pm 2.1$ | $5.0\pm 2.1$ | $2.0\pm 1.6$ | $0.4\pm 0.5$
err. $<10\%$ | $93.8\pm 1.1$ | $59.0\pm 5.2$ | $32.6\pm 2.3$ | $13.6\pm 2.3$ | $5.6\pm 1.9$ | $2.2\pm 1.9$ | $1.2\pm 0.8$
err. $<20\%$ | $99.2\pm 0.4$ | $80.0\pm 4.8$ | $55.0\pm 6.4$ | $33.8\pm 2.9$ | $21.8\pm 4.8$ | $11.0\pm 3.1$ | $7.4\pm 2.9$
% of conv. | $99.4\pm 0.5$ | $83.2\pm 5.5$ | $62.4\pm 7.2$ | $40.8\pm 3.6$ | $29.8\pm 7.2$ | $22.2\pm 2.9$ | $19.6\pm 5.6$
RMSE | $0.8\pm 0.2$ | $4.3\pm 0.1$ | $8.3\pm 1.1$ | $11.9\pm 1.0$ | $15.1\pm 2.8$ | $22.1\pm 2.5$ | $24.7\pm 3.1$
Newton-Kaczmarz (NK) method
err. $<5\%$ | $98.4\pm 0.9$ | $78.2\pm 3.0$ | $49.4\pm 4.7$ | $35.8\pm 4.1$ | $24.8\pm 2.6$ | $16.0\pm 5.1$ | $13.0\pm 4.3$
err. $<10\%$ | $99.8\pm 0.4$ | $95.0\pm 2.0$ | $78.0\pm 7.2$ | $56.2\pm 3.6$ | $42.2\pm 4.3$ | $30.2\pm 5.8$ | $20.4\pm 4.5$
err. $<20\%$ | $100\pm 0.0$ | $98.8\pm 1.1$ | $90.8\pm 3.8$ | $75.4\pm 1.7$ | $64.4\pm 3.2$ | $48.0\pm 5.3$ | $34.8\pm 7.2$
Table 1: The average number of runs (executions of the code) out of $100$,
where the normalised root-mean-square error was below $5\%$, $10\%$ and $20\%$
for different amplitude of the perturbations of the initial guess, parameter
$\alpha$, for the GN method (lines 3-5) and for the NK method (lines 9-11).
The average number of runs out of $100$, where the convergence criteria of the
GN method was fulfilled for different $\alpha$ (line 6). The average
normalised root-mean-square error, calculated only for the converged runs of
the GN method, for different $\alpha$ (line 7). In each cell, the mean and the
standard deviation across $5$ ensembles, each consisting of $100$ runs, are
given.

The percentages of accurate runs of the methods might seem to be relatively
small, but this example is specifically chosen to be challenging for the
methods. In particular, the exact values of $G_{l}$ given in equation (30) are
such that outer function $\varPhi$ is highly non-linear and non-monotonous.
Furthermore, since the model is the ridge function, there is no redundancy —
it consists of only one addend that must be identified precisely, otherwise
the error is large. This means that the convergence of the methods to local
minima (with the identified parameters that are different from the exact model
parameters) is insufficient in terms of accuracy.

In the comparison above, all runs of the GN method were used, even if the
convergence criteria was not fulfilled and the iterations were stopped after
the maximum number of iterations has been reached. An alternative approach is
that the runs where the convergence criteria is not fulfilled are disregarded.
However, this implies that if the convergence criteria is fulfilled, then the
result should be considered to be accurate. In this case, first, the
percentage of runs with the ‘converged’ results drops with $\alpha$, as shown
in line 6 of table 1. Second, the average RMSE, calculated only for the
‘converged’ runs, still increases with $\alpha$ and reaches almost $25\%$ for
$\alpha=2.8$, as shown in line 7 of table 1. This means that the runs converge
to local minima, providing insufficiently accurate result.

The percentages in table 1 were calculated based on $5$ ensembles, each
consisting of $100$ runs. Within each ensemble, the number of runs with the
error below a given threshold was calculated; then the mean and the standard
deviation across the ensembles were calculated and are given in the tables.
The same was done for the percentage of the converged runs of the GN method
and the corresponding RMSE.

### 4.2 Example 2: Approximation of a non-linear function

The purpose of this example is to demonstrate that the Kolmogorov-Arnold model
can indeed approximate a non-linear function and that the proposed method
solves the identification problem. The code of this example is available at
the GitHub repository555https://github.com/andrewpolar/UStandardM. The ‘full’
model is taken with piecewise-linear basis functions.

The following non-linear function of $m=5$ inputs is taken:

$\begin{split}&y_{i}=\frac{2+2X_{3}^{i}}{3\pi}\left(\operatorname{arctan}\left(20\left(X_{1}^{i}-\frac{1}{2}+\frac{X_{2}^{i}}{6}\right)\exp\left(X_{5}^{i}\right)\right)+\frac{\pi}{2}\right)+\\\
&\vphantom{y}+\frac{2+2X_{4}^{i}}{3\pi}\left(\operatorname{arctan}\left(20\left(X_{1}^{i}-\frac{1}{2}-\frac{X_{2}^{i}}{6}\right)\exp\left(X_{5}^{i}\right)\right)+\frac{\pi}{2}\right).\end{split}$
(33)

The datasets consisting of $N=10^{4}$ records for the training and
$N_{\mathrm{val}}=2000$ records for the validation have been generated. This
involved taking uniformly-distributed random inputs
$X^{i}_{j}\sim\operatorname{unif}\left(0,1\right)$ and calculating outputs
$y_{i}$ using equation (33). Using the data, the Kolmogorov-Arnold models with
$n=5$ inner and $s=7$ outer piecewise-linear basis functions were identified.
The models were identified using the training data only, while the RMSE was
calculated using the validation data only.

Since the identification method is iterative, it is interesting to study the
convergence rate with respect to the number of iterations. The number of
iterations should be large, but each iteration is computationally cheap, as
each iteration of the Newton-Kaczmarz method operates only with one record.
Therefore, the algorithm goes through all the training records consecutively
multiple times. After each pass666In machine learning literature, one ‘pass’
through the training data is also called ‘epoch’. through the training data,
the RMSE is calculated using the validation data. Thus, the RMSE can be
plotted as a function of the number of passes through the data, which is
illustrated in figure 1. Three different values of the regularisation
parameter were used; for each case, $10$ runs with different data and
different initial guesses for the parameters were performed.

Figure 1: The dependence of the RMSE on the number of passes through the data
for different values of $\mu$.

It can be seen that for any $\mu$, the method converges approximately
logarithmically with respect to the number of passes, i.e. there is
approximately linear dependence between $\log E_{\mathrm{RMSE}}$ and $\log
N_{\mathrm{pass}}$. However, smaller $\mu$ leads to higher absolute value of
$E_{\mathrm{RMSE}}$ for the same number of passes. The practical purpose of
decreasing $\mu$ is to filter out the noise, which is absent in this example.
A more detailed study of the noise filtering abilities of the Kaczmarz method
was performed for the case of the Urysohn model in [15]. Finally, it should be
mentioned that the model obtained using the proposed identification technique
is accurate, with the RMSE of around $0.5\%$ for $\mu=1$ after $500$ passes,
as seen from figure 1.

## 5 Conclusions

In this paper, the Kolmogorov-Arnold model has been considered. It is a
convenient tool for data modelling tasks, which consist of mapping the input
vectors into the output scalars. It has been shown that the underlying
functions constituting the model can be represented via the basis functions
and the model parameters. To determine these parameters, an iterative
identification algorithm that is based on the Newton-Kaczmarz method has been
proposed. Compared to the previously-proposed algorithm [14] that was not
backed up by the proof of the convergence, the present algorithm rests upon
the local convergence of the Newton-Kaczmarz method.

The proposed algorithm has been tested using the numerical examples. First, it
was compared to the straightforward application of the Gauss-Newton method to
identify the parameters of the ridge function model. A challenging example has
been selected, where the non-linearity of the outer function creates
convergence difficulties when the initial guess for the parameters is far from
their exact values. It has been shown that the algorithm based on the Newton-
Kaczmarz method is much less sensitive to the selection of the initial guess;
therefore, it is more robust from a practical point of view. The second
example has demonstrated that a challenging non-linear function can be
approximated accurately by the Kolmogorov-Arnold model with the piecewise-
linear basis functions, and the model parameters can be obtained using the
proposed algorithm. Finally, it should be mentioned that many more examples of
the application of the Urysohn and the Kolmogorov-Arnold models to various
synthetic and real datasets, as well as all codes of the examples, are
available at the website of the authors777http://www.ezcodesample.com/ and the
GitHub repository888https://github.com/andrewpolar.

##  Appendix A The Newton-Kaczmarz method

For non-linear system of equations $L\left(Z\right)=0$, where
$L:D\to\mathbb{R}^{N}$ is continuously differentiable in open set
$D\subset\mathbb{R}^{r}$, the Newton-Kaczmarz method [22, 23, 24] consists in
iterations

$Z^{q+1}=Z^{q}-\mu\frac{L^{i}\left(Z^{q}\right)}{\left\|\nabla
L^{i}\left(Z^{q}\right)\right\|^{2}}{\nabla
L^{i}\left(Z^{q}\right)},\quad\quad\nabla
L^{i}\left(Z^{q}\right)=\left.\begin{bmatrix}\frac{\partial L^{i}}{\partial
Z_{1}}&\ldots&\frac{\partial L^{i}}{\partial
Z_{r}}\end{bmatrix}^{\mathrm{T}}\right|_{Z=Z^{q}},$ (34)

where function $L^{i}$ denotes the $i$-th component of vector function $L$,
scalar $Z_{j}$ denotes the $j$-th component of vector $Z$, parameter
$\mu\in\left(0,2\right)$ is introduced for regularisation and vector $Z^{q}$
is the approximation of the solution at iteration $q$. Index $i$ changes each
iteration. The method requires an initial guess and continues until some
convergence criteria are reached.

This method can also be interpreted as an iterative scheme, where at each
step, a single equation is linearised and one step of the Kaczmarz method [20,
21] is performed. Indeed, to write an iterative scheme, one must find an
expression for $Z^{q+1}$, assuming that $Z^{q}$ is known. Considering single
equation $L^{i}\left(Z^{q+1}\right)=0$, the linearisation can be performed:

$L^{i}\left(Z^{q+1}\right)=L^{i}\left(Z^{q}\right)+\nabla
L^{i}\left(Z^{q}\right)^{\mathrm{T}}\Delta Z+O\left(\left\|\Delta
Z\right\|^{2}\right)=0,\quad\quad\Delta Z=Z^{q+1}-Z^{q}.$ (35)

When the second-order term is neglected, linear equation

$M^{i}\Delta Z=Y_{i}$ (36)

is obtained, where $M^{i}=\nabla L^{i}\left(Z^{q}\right)^{\mathrm{T}}$ and
$Y_{i}=-L^{i}\left(Z^{q}\right)$. Equation (36) with respect to $\Delta Z$
describes a hyperplane in $\mathbb{R}^{r}$. Using the key idea of the Kaczmarz
method (the projection descend), the origin is projected onto this hyperplane
to obtain the solution:

$\Delta Z=\frac{Y_{i}}{\left\|M^{i}\right\|^{2}}{M^{i}}^{\mathrm{T}}.$ (37)

Adding regularisation parameter $\mu$ to limit the magnitude of change of $Z$
over one step, leads to equation (34).

##  Appendix B Identification of the ridge function using the Gauss-Newton
method

The ridge function model with the outer function represented via the
parameters and the basis functions is given by

$\hat{y}_{i}=\sum_{l=1}^{s}G_{l}\psi^{l}\left(\sum_{j=1}^{m}c_{j}X_{j}^{i}\right),$
(38)

which is obtained by substituting equation (28) into (7). The sum of squares
of the residuals is introduced as

$S=\frac{1}{2}\sum_{i=1}^{N}\left(\hat{y}_{i}-y_{i}\right)^{2}.$ (39)

The Gauss-Newton method consists in minimisation of $S$ with respect to the
model parameters using the Newton-Raphson method and neglecting the second-
order derivative terms. However, to obtain the most accurate Hessian matrix of
$S$ and to maximise the performance of the method, the second-order
derivatives are not neglected here. Thus, although the method deviates from
the classical Gauss-Newton algorithm, for the purposes of this paper, it is
still referred to as such.

The minimisation of $S$ gives

$\frac{\partial S}{\partial G_{l}}=\sum_{i=1}^{N}r_{i}\frac{\partial
r_{i}}{\partial G_{l}}=0,\quad\frac{\partial S}{\partial
c_{j}}=\sum_{i=1}^{N}r_{i}\frac{\partial r_{i}}{\partial c_{j}}=0,\quad
l\in\left\\{1,\ldots,s\right\\},\quad j\in\left\\{1,\ldots,m\right\\},$ (40)

with the following additionally-introduced quantities:

$\displaystyle
r_{i}=\hat{y}_{i}-y_{i},\quad\quad\quad\theta_{i}=\sum_{j=1}^{m}c_{j}X_{j}^{i},\quad\quad\quad\eta^{l}=\frac{\mathrm{d}\psi^{l}}{\mathrm{d}t},$
(41) $\displaystyle\frac{\partial r_{i}}{\partial
G_{l}}=\psi^{l}\left(\theta_{i}\right),\quad\quad\quad\frac{\partial
r_{i}}{\partial
c_{j}}=X_{j}^{i}\sum_{l=1}^{s}G_{l}\eta^{l}\left(\theta_{i}\right).$ (42)

The evaluation of the second derivatives of $S$ leads to

$\displaystyle\frac{\partial^{2}S}{\partial G_{l}\partial
G_{p}}=\sum_{i=1}^{N}\frac{\partial r_{i}}{\partial G_{l}}\frac{\partial
r_{i}}{\partial G_{p}},\quad\quad\quad\frac{\partial^{2}S}{\partial
G_{l}\partial c_{j}}=\sum_{i=1}^{N}\left(\frac{\partial r_{i}}{\partial
G_{l}}\frac{\partial r_{i}}{\partial
c_{j}}+r_{i}\frac{\partial^{2}r_{i}}{\partial G_{l}\partial c_{j}}\right),$
(43) $\displaystyle\frac{\partial^{2}S}{\partial c_{j}\partial
c_{k}}=\sum_{i=1}^{N}\left(\frac{\partial r_{i}}{\partial c_{j}}\frac{\partial
r_{i}}{\partial c_{k}}+r_{i}\frac{\partial^{2}r_{i}}{\partial c_{j}\partial
c_{k}}\right),\quad\quad\quad p\in\left\\{1,\ldots,s\right\\},\quad
k\in\left\\{1,\ldots,m\right\\},$ (44)

with the following additionally-introduced quantities:

$\frac{\partial^{2}r_{i}}{\partial G_{l}\partial
c_{j}}=X_{j}^{i}\eta^{l}\left(\theta_{i}\right),\quad\quad\quad\frac{\partial^{2}r_{i}}{\partial
c_{j}\partial
c_{k}}=X_{j}^{i}X_{k}^{i}\sum_{l=1}^{s}G_{l}\xi^{l}\left(\theta_{i}\right),\quad\quad\quad\xi^{l}=\frac{\mathrm{d}^{2}\psi^{l}}{\mathrm{d}t^{2}}.$
(45)

Using the obtained relations, equations (40) are solved using the Newton-
Raphson method, consisting in the following iterative scheme:

$Z^{q+1}=Z^{q}-\mu\left.M^{-1}L\right|_{Z=Z^{q}},$ (46)

where $Z$ is the vector containing the model parameters, with superscript
indicating the iteration number, $M$ and $L$ are the matrix and the vector
containing the first and the second derivatives of $S$ with respect to the
model parameters, respectively, and $\mu$ is the regularisation parameter
added for the robustness of the scheme999In the example of section 4.1, the
initial value of $\mu$ is $0.1$. As soon as $\left\|Z^{q+1}-Z^{q}\right\|<0.1$
is fulfilled, the value of $\mu$ is changed $1$.. These quantities are defined
as

$\displaystyle
Z=\begin{bmatrix}G_{1}&\ldots&G_{s}&c_{1}&\ldots&c_{m}\end{bmatrix}^{\mathrm{T}},$
(47) $\displaystyle L=\begin{bmatrix}\frac{\partial S}{\partial
G_{1}}&\ldots&\frac{\partial S}{\partial G_{s}}&\frac{\partial S}{\partial
c_{1}}&\ldots&\frac{\partial S}{\partial c_{m}}\end{bmatrix}^{\mathrm{T}},$
(48) $\displaystyle M_{l,p}=\frac{\partial^{2}S}{\partial G_{l}\partial
G_{p}},\quad\quad\quad M_{l,s+j}=M_{s+j,l}=\frac{\partial^{2}S}{\partial
G_{l}\partial c_{j}},\quad\quad\quad M_{s+j,s+k}=\frac{\partial^{2}S}{\partial
c_{j}\partial c_{k}},$ (49)

where comma is added between the subscripts indicating the components of
matrix $M$ for clarity. The iterations are performed until either the maximum
number of iterations is reached or $\left\|Z^{q+1}-Z^{q}\right\|<\delta$ is
fulfilled101010In the example of section 4.1, a more strict convergence
criteria, also requiring $\left\|L^{q}\right\|<\delta$, were investigated
additionally (results are not shown). However, no significant change of the
results was observed..

## References

  * [1] V. I. Arnold.  On functions of three variables.  Doklady Akademii Nauk SSSR, 114(4):679–681, 1957.
  * [2] A. N. Kolmogorov.  On the representation of continuous functions of many variables by superposition of continuous functions of one variable and addition.  Doklady Akademii Nauk SSSR, 114(5):953–956, 1957.
  * [3] D. Bryant.  Analysis of Kolmogorov’s superposition theorem and its implementation in applications with low and high dimensional data.  PhD thesis, University of Central Florida, 2008.
  * [4] X. Liu.  Kolmogorov superposition theorem and its applications.  PhD thesis, Imperial College London, 2015.
  * [5] D. A. Sprecher.  A numerical implementation of Kolmogorov’s superpositions.  Neural Networks, 9(5):765–772, 1996.
  * [6] D. A. Sprecher.  A numerical implementation of Kolmogorov’s superpositions II.  Neural Networks, 10(3):447–457, 1997.
  * [7] M. Köppen.  On the training of a Kolmogorov network.  In J. R. Dorronsoro, editor, Artificial Neural Networks — ICANN 2002, pages 474–479. Springer Berlin Heidelberg, 2002.
  * [8] J. Braun and M. Griebel.  On a constructive proof of Kolmogorov’s superposition theorem.  Constructive Approximation, 30(3):653–675, 2009.
  * [9] J. Actor and M. J. Knepley.  An algorithm for computing Lipschitz inner functions in Kolmogorov’s superposition theorem.  arXiv:1712.08286, 2017.
  * [10] J. Actor.  Computation for the Kolmogorov superposition theorem.  Master’s thesis, Rice University, 2018.
  * [11] M. Nees.  Approximative versions of Kolmogorov’s superposition theorem, proved constructively.  Journal of Computational and Applied Mathematics, 54(2):239–250, 1994.
  * [12] B. Igelnik and N. Parikh.  Kolmogorov’s spline network.  IEEE Transactions on Neural Networks, 14(4):725–733, 2003.
  * [13] M. Coppejans.  On Kolmogorov’s representation of functions of several variables by functions of one variable.  Journal of Econometrics, 123(1):1–31, 2004.
  * [14] A. Polar and M. Poluektov.  A deep machine learning algorithm for construction of the Kolmogorov-Arnold representation.  Engineering Applications of Artificial Intelligence, 99:104137, 2021\.
  * [15] M. Poluektov and A. Polar.  Modelling non-linear control systems using the discrete Urysohn operator.  Journal of the Franklin Institute, 357:3865–3892, 2020.
  * [16] M. A. Krasnoselskii.  Topological methods in the theory of nonlinear integral equations.  Pergamon, 1964.
  * [17] V. V. Krylov.  Modeling the internal structure of dynamical systems from input-output relationships (abstract realization theory) 2.  Automation and Remote Control, 45(3):277–288, 1984.
  * [18] T. J. Hastie and R. J. Tibshirani.  Generalized additive models.  Chapman and Hall/CRC, 1990.
  * [19] V. E. Ismailov.  A review of some results on ridge function approximation.  Azerbaijan Journal of Mathematics, 3(1):3–51, 2013.
  * [20] S. Kaczmarz.  Angenäherte auflösung von systemen linearer gleichungen.  Bulletin international de l’Acadámie polonaise des sciences et des lettres. Classe des sciences mathématiques et naturelles. Série A, Sciences mathématiques, 35:355–357, 1937.
  * [21] R. P. Tewarson.  Projection methods for solving sparse linear systems.  The Computer Journal, 12(1):77–80, 1969.
  * [22] K.-H. Meyn.  Solution of underdetermined nonlinear equations by stationary iteration methods.  Numerische Mathematik, 42(2):161–172, 1983.
  * [23] J. M. Martinez.  Solving systems of nonlinear equations by means of an accelerated successive orthogonal projections method.  Journal of Computational and Applied Mathematics, 16(2):169–179, 1986.
  * [24] J. M. Martinez and R. J. B. de Sampaio.  Parallel and sequential Kaczmarz methods for solving underdetermined nonlinear equations.  Journal of Computational and Applied Mathematics, 15(3):311–321, 1986.


# Commissioning of a clinical pencil beam scanning proton therapy unit for
ultrahigh dose rates (FLASH)

Konrad P. Nesteruk<EMAIL_ADDRESS>Michele Togno  Martin Grossmann
Paul Scherrer Institute, Center for Proton Therapy, Villigen, Switzerland
Anthony J. Lomax  Paul Scherrer Institute, Center for Proton Therapy,
Villigen, Switzerland
ETH Zurich, Department of Physics, Villigen, Switzerland  Damien C. Weber
Paul Scherrer Institute, Center for Proton Therapy, Villigen, Switzerland
University Hospital Zurich, Department of Radiation Oncology, Switzerland
University Hospital Bern, Department of Radiation Oncology, Switzerland
Jacobus M. Schippers  Paul Scherrer Institute, Division Large Research
Facilities, Villigen, Switzerland  Sairos Safai  David Meer  Serena Psoroulas
Paul Scherrer Institute, Center for Proton Therapy, Villigen, Switzerland

###### Abstract

Purpose: The purpose of this work was to provide a flexible platform for FLASH
research with protons by adapting a former clinical pencil beam scanning
gantry to irradiations with ultrahigh dose rates.
Methods: PSI Gantry 1 treated patients until December 2018. We optimized the
beamline parameters to transport the 250 MeV beam extracted from the PSI COMET
accelerator to the treatment room, maximizing the transmission of beam
intensity to the sample. We characterized a dose monitor on the gantry to
ensure good control of the dose, delivered in spot-scanning mode. We
characterized the beam for different dose rates and field sizes for
transmission irradiations. We explored scanning possibilities in order to
enable conformal irradiations or transmission irradiations of large targets
(with transverse scanning).
Results: We achieved a transmission of 86 % from the cyclotron to the
treatment room. We reached a peak dose rate of 9000 Gy/s at 3 mm water
equivalent depth, along the central axis of a single pencil beam. Field sizes
of up to 5x5 mm2 were achieved for single spot FLASH irradiations. Fast
transverse scanning allowed to cover a field of 16x1.2 cm2. With the use of a
nozzle-mounted range shifter we are able to span depths in water ranging from
19.6 to 37.9 cm. Various dose levels were delivered with a precision within
less than 1 %. Conclusions: We have realized a proton FLASH irradiation setup
able to investigate continuously a wide dose rate spectrum, from 1 to 9000
Gy/s in a single spot irradiation as well as in the pencil beam scanning mode.
As such, we have developed a versatile test bench for FLASH research.

Ultrahigh dose rates, FLASH, proton therapy, gantry, pencil beam scanning

##  I Introduction

One of the limiting factors in delivering radiation therapy to particularly
aggressive or radiation resistant tumors is the limitation on the amount of
dose healthy tissues can withstand. In recent years, several studies have
indicated that ultrahigh dose rates might result in reduced toxicities to
healthy tissues while keeping the same tumor control as for treatments with
standard dose rate levels. In 2014, the so called FLASH effect has been
demonstrated for the first time by irradiating lung tumors in mice with 4.5
MeV electrons Favaudon2014 . A reduced pulmonary fibrosis was observed with
electrons for a dose rate of 60 Gy/s when compared to a conventional dose rate
of 0.03 Gy/s. Whole brain irradiation of mice with 4.5 MeV electrons showed a
unique memory sparing for dose rates above 100 Gy/s Montay2017 . Almost no
skin toxicity has been observed in mini-pigs irradiated with electrons, as
well as in cat patients with locally advanced squamous cell carcinoma of the
nasal planum Vozenin2019 . Eventually, in Switzerland, the first ever FLASH
treatment of a human patient has been performed in 2018 at the Lausanne
University Hospital (CHUV) on a prototype linac for FLASH radiotherapy with
electrons Bourhis2019 .

Although the majority of the promising studies of the FLASH effects was
performed with electrons, it is conceptually possible that this FLASH effect
could be observed with proton irradiation. Several centers in the world, as
well as main vendors in the field, started to investigate the possibility of
reaching FLASH dose rates with existing treatment units Patriarca2018 ;
Beyreuther2019 ; Buonanno2019 ; Darafsheh2020 ; Diffenderfer2020 . Based on
the studies with photons and electrons, it has been concluded that the
threshold for the FLASH effect is at least 40 Gy/s, at least a factor 10
higher than conventional dose rates. However, the definition of the dose rate
is prone to ambiguities, as the beam has its microstructure and the average
dose rate may be substantially different from the maximum dose rate in a short
pulse. The definition of the dose rate is also difficult for pencil beam
scanning, due to the sequential character of the delivery technique. Since the
mechanism of FLASH has not been fully explained, a potential influence of
pauses in dose delivery and fractionation is not understood yet. Therefore,
exact conditions for the FLASH effect to occur remain unknown, which reflects
in non-consistent results of studies performed with protons so far. Thus,
flexible FLASH test benches, able to provide different dose rates and
irradiation conditions, are required.

In this paper, we report on the commissioning of the former clinical unit
Gantry 1 at the Paul Scherrer Institute for FLASH research with protons. PSI
Gantry 1 was the first gantry in the world to deliver routinely pencil beam
scanning (PBS) to cancer patients. The objective of this work was to transform
the in-house developed gantry into a versatile platform for FLASH effect
studies. As such, the modified gantry will be able to provide various dose
rates and allow for single spot irradiations as well as pencil beam scanning
with ultrahigh instantaneous dose rates for biological samples and potentially
patients. Dose delivery in both transmission and conformal (Bragg peak) modes
will be possible.

##  II Materials and Methods

###  II.1 Pencil beam scanning Gantry 1

Figure 1: Gantry 1 at an angle of -90 degrees seen from the treatment room.
Highlighted are spot-scanning axes and motions, the ”isocenter” defined for
this paper, and a Faraday cup used for experiments.

Gantry 1 treated patients from 1996 until the end of 2018. It was the first
facility worldwide that used the spot-scanning technique Pedroni1995 ; Lin2009
. With its radius of 2 m, it is still the most compact gantry in the world
competing advantageously with low footprint gantries sold by vendors. The
compactness was achieved by mounting the patient table off-axis (excentric
gantry) as well as by placing a scanning magnet upstream the last dipole.
Since the gantry is excentric and the table moves vertically together with the
gantry rotation, there is no true isocenter. In this paper, whenever we refer
to the isocenter, we have in mind a moving isocenter defined as a rotation-
invariant distance between the nozzle and the axis of the table at its
reference position (Fig. 1).

Magnetic scanning is realized in one transverse direction (U) by means of a
fast scanning magnet, so-called sweeper, with a dead time of a few ms between
spots. Scanning in depth (S) is achieved with a range shifter in the gantry
nozzle, which consists of 39 polystyrene (PS) plates of 4.53 mm thickness and
one PS plate of half thickness. They are moved pneumatically in and out of the
beam with a dead time of 60 ms per plate. In view of the FLASH applications,
this feature gives a unique opportunity to transport an un-degraded proton
beam to the isocenter in order to maximize the beam intensity, hence the dose
rate, while keeping the possibility of conformal (Bragg peak) irradiations.
The slowest scanning direction (T) is performed by moving a treatment table.
An overview of the gantry, as seen from the treatment room, is presented in
Fig. 1.

Figure 2: Schematic diagram of the Gantry 1 beamline. Selected components
relevant to this paper are highlighted.

###  II.2 Beam transport for ultrahigh dose rates

Proton beam used for treatments is produced by the COMET cyclotron
Schippers2007 . The energy and intensity of the extracted beam are 250 MeV and
up to 1 $\mu$A, respectively. Energy modulation is performed with a degrader
system, installed right after the cyclotron. In the case of Gantry 1 the
previously mentioned range shifter was used for an additional fine modulation
of the energy pre-selected by the degrader. However, the degrader introduces
massive intensity losses due to the scattering occuring in the degrader
material. At low energies, the fraction of the beam transported to the
treatment room is lower than 1 %. In order to reach ultrahigh dose rates, we
had to optimize the beamline parameters to transport the 250 MeV beam produced
by the PSI COMET cyclotron to the treatment room with minimum losses. This
energy has never been used for treatments, although the gantry was design to
accept energies up to 300 MeV. Therefore, a completely new configuration of
beamline magnet settings, so-called beamline tune, had to be found.

The beamline (Fig. 2), including the gantry, is 44 m long and consists of
about 40 magnets - dipoles and quadrupoles. Beam transport can be simulated in
numerous codes for beam optics calculations. We used a simulation by means of
TRANSPORT transport1 ; transport2  to find initial magnet settings and as a
guide for further experimental fine tuning. We considered the following
criteria for the desired beam optics:

  * •

beam size (2-sigma) always much smaller than apertures of magnets and
collimators

  * •

point-to-point imaging between the location of the degrader (upstream) and the
gantry coupling point

  * •

limited beam tilt at the entrance of the gantry

  * •

beam waist in the gantry at the same location in both bending (dispersive) and
non-bending planes

  * •

limited beam divergence of the beam extracted from the gantry

  * •

round beam spot at a certain distance from the gantry

  * •

minimized dispersion and its derivative at the isocenter.

The experimental tuning was based on numerous beam profilers and intensity
monitors located along the beamline. It was realized in two stages. The goal
of the first stage was to transport the beam to the gantry coupling point with
minimum loses. At the second stage, we optimized the transport through the
gantry to achieve the maximum transmission, a round beam spot at a defined
distance from the gantry nozzle, and a limited beam divergence in air. For
these optimizations, we used a scintillating foil and a in-room camera for
online tuning, and a scintillator-CCD detector for precise measurements of the
beam size and its envelope in air. The beam transmission to the gantry
coupling point was determined as a ratio between coupling-point and cyclotron-
extracted beam intensities, measured simultaneously by means of intensity
monitors installed along the beamline. The beam current at the isocenter was
calculated based on the total proton charge measured by a Farady cup (see Fig.
1) and measured delivery time. The transmission from the coupling point to the
isocenter was defined as a ratio between the calculated Faraday cup current
and the current measured at the coupling point.

###  II.3 Dose monitoring in the FLASH regime

Gantry 1 is equipped with two parallel-plane ionization chambers placed in the
nozzle to control the dose to be delivered Lin2009 . The first chamber
(monitor 1) is the main monitor which controls dose delivery with an ion
collection time of 90 $\mu$s, while the second chamber (monitor 2) with a
collection time of 350 $\mu$s is a backup element and it was part of the
therapy verification system. Assuming a FLASH irradiation time of the order of
10 ms, monitor 1 guarantees a precision of at least 1 %.

Compared to the clinical setting, the instantaneous dose rate we expect in
FLASH experiments is 100 to 1000 times higher. Ionization chambers at such
high dose rates are affected by significant ion recombination losses and
therefore, a thorough characterization of these kind of monitors is needed
dosimetry_preprint . In order to provide an accurate dose monitoring, we have
performed such a calibration of the nozzle monitor against a Faraday cup (see
Fig. 1). We used the same Faraday cup to validate the precision of the dose
delivery in terms of reproducibility. The prediction of the absolute dose with
different dose rates was based on a model described in dosimetry_preprint .

###  II.4 Delivery of different dose rates

For the definition of dose rate we assume the beam produced by our cyclotron
is continuous, as the RF frequency of the COMET cyclotron is 72.85 MHz. As
such, every 14 ns a 0.8 ns long pulse is delivered. Therefore, in a single
pulse the dose rate is 17.5 times higher than under the assumption of
continuous beam. However, a single pulse is so short that the dose delivered
within it is extremely low. Even for a FLASH dose rate of 1000 Gy/s, the dose
per pulse is of the order of 10 $\mu$Gy. Also in terms of ion recombination,
the dose delivery can be considered continuous Goma2014 .

####  II.4.1 Single spot irradiations

Single spot transmission irradiations are used to deliver homogeneously a dose
to small targets, such as cell lines. In order to obtain homogeneity in depth
we use the flat part of the Bragg peak curve on the central axis where the
dose, the beam size, and hence the dose rate, is typically constant for a few
centimeters, as depicted in Fig. 3 and described in dosimetry_preprint . For
the transverse homogeneity, we define the field size to be the 95 % isodose
surface for a given Gaussian pencil beam.

Figure 3: Comparison of Bragg peak curves - integral and on the central axis
for a beam energy of 230 MeV. The flat part of the curve with the maximum dose
rate is highlighted.

The instantaneous dose rate can be varied by adjusting the beam current and
beam size. The former is realized by a vertical deflector in the COMET
cyclotron. As such, we can vary the cyclotron current between 0.1 and 800 nA.
Up to 20 % higher beam currents are technically achievable. However, they are
not clinically commissioned and thus were not used in this study. The
variation of beam size can be achieved by choosing a different number of range
shifter plates in the gantry nozzle. Due to multiple Coulomb scattering the
beam size increases and hence the dose rate decreases in a quadratic manner.
This can also be used to increase the field size with 95 % homogeneity at the
cost of decreased dose rate.

In order to characterize the beam and estimate the corresponding dose rate we
measured beam sizes by means of a CCD-scintillator detector and depth-dose
curves in water by means of a large sensitive area, plane-parallel Bragg peak
ionization chamber. The dose rates were determined from the measured delivered
protons by means of a Faraday cup and recorded delivery time. The time is
recorded by the control system and corrected offline for the reaction time of
an upstream kicker magnet which switches the beam off when the requested dose
is reached. As such, an accuracy better than 50 $\mu$s is achieved.

####  II.4.2 Pencil beam scanning

We explored scanning possibilities in order to enable conformal irradiations
or transmission irradiations of large targets (transverse scanning).

Since the upstream degrader is not used anymore, the only possibility to scan
in depth is to use the range shifter in the gantry nozzle. We verified the
minimum reachable energy with all the range shifter plates inserted by
performing a range measurement in water.

For scanning in the transverse U direction (the bending direction), we
evaluated the performance of the scanning magnet (sweeper) at an energy of 250
MeV, as this energy has never been used for treatments. For the main FLASH-
optimized tune, a new scanning magnet calibration, so-called sweeper map, was
defined to obtain equal spacing between spots. Additionally, we explored the
possibility of fast scanning in the T direction (the non-dispersive direction)
with the use of a steering magnet usually utilized for the first spot
correction (beam centering). This magnet is not as fast as the sweeper. With a
dead time between consecutive spots of the order of 100 ms it would, however,
guarantee much faster scanning in T than the table moving. Eventually, a 2D
scanning map has been determined. For all the transverse scanning
measurements, a CCD-scintillator detector was used.

##  III Results

###  III.1 FLASH-optimized beam transport

Figure 4: Beam envelope corresponding to the main FLASH beamline tune.
Negative values of beam size correspond to the dispersive plane (U) and
positive values to the non-dispersive plane (T). Continuous line - result of
simulations with TRANSPORT; red dots - beam profile measurements.

We characterized two different beamline tunes to be used for FLASH
experiments. The main FLASH-optimized tune is shown in Fig. 4. The measured
beam envelope is compared with the one simulated in TRANSPORT with the final
magnet settings. Due to the design of the compact gantry, there is no beam
profile monitor in the gantry. The observed discrepancies are mostly due to
the phase space definition. The TRANSPORT simulation is based on the phase
space measurements performed in 2006. We updated the phase space (highlighted
in the figure) on the basis of upstream beam profile measurements. However,
the beam emittance was assumed to be the same as measured in 2006. The
agreement we reached is satisfactory for our purpose and the simulation in
TRANSPORT served as a useful guide for the beamline tuning, providing an
important qualitative description of our beam transport. The tune presented in
Fig. 4 provides the highest transmission to the isocenter and the minimum
eccentricy of the beam spot, even when the range shifter is not used. The beam
divergence in air was measured to be ($0.5\pm 0.3$) mrad and ($2.0\pm 0.3$)
mrad for the dispersive (U) and non-dispersive (T) plane, respectively. At a
distance of 41 cm from the nozzle, 5 cm behind the isocenter, the beam spot
was found to be perfectly round with a 1$\sigma$ size of ($3.5\pm 0.1$) mm.
The tune reproducibility has been assessed based on several measurements of
beam profiles in different days. The measured beam sizes along the beamline
were found to be in agreement within 5 % and those in air within 3 %.

The transmission to the gantry coupling point was measured to be ($86\pm 1$) %
for cyclotron currents exceeding 500 nA. A weak dependence on the beam current
was observed, as different settings of the vertical deflector slightly affect
the initial phase space and hence the beam transport. The minimum transmission
was measured to be ($82\pm 1$) % for the lowest cyclotron currents, 1 nA and
below. The transmission from the coupling point to the isocenter was
determined to be 100 % and independent of the beam current.

###  III.2 Ultrahigh dose rates

####  III.2.1 Single spot transmission irradiations

With the main FLASH tune we have achieved field sizes ranging from 2.3x2.3 mm2
for a beam energy of 250 MeV to 5x5 mm2 for a beam energy of 170 MeV (all
range shifter plates inserted). The corresponding maximum peak dose rates are
$(3600\pm 200)$ Gy/s and $(700\pm 30)$ Gy/s, respectively. The mentioned peak
dose rates correspond to 3.5 cm water equivalent depth and a cyclotron beam
current of 800 nA. For a given field size, we are able to continuously vary
the dose rate from the level of the conventional PBS dose rate (order of 1
Gy/s) to the maximum dose rate corresponding to the given field size.

The absolute maximum we have achieved was $(9300\pm 500)$ Gy/s at 3 mm water
equivalent depth. This dose rate was reached with another FLASH tune, which
provided a narrower pencil beam. However, for this tune the beam spot is not
round, when no range shifter plate is used. At a distance of 56 cm from the
gantry nozzle, 10 cm behind the isocenter, 1$\sigma$ sizes of the beam spot
were measured to be ($2.27\pm 0.05$) mm and ($1.80\pm 0.05$) mm in the U and T
directions, respectively. As such, the maximum field size with this beam spot
for transmission irradiations is 1.1x1.1 mm2. The characterization of the beam
spot in air for this maximum dose rate is depicted in Fig. 5.

Figure 5: Characterization of a beam spot in air. Left: Beam spot of 2.3/1.8
mm (U/T \- 1$\sigma$) in air recorded with a scintillator-CCD detector. Right:
Corresponding dose rate distribution in water at 3 mm water equivalent depth.

The uncertainty of the dose rate consists of the model uncertainty, dose
delivery uncertainty and the beam size uncertainty. The beam size uncertainty
contributes most to the total uncertainty. In our characterization studies,
the uncertainty due to the delivery time measurement was negligible due to 3
to 4 orders of magnitude longer delivery times than the precision of the time
measurement.

Figure 6: Efficiency of monitor 1 in the gantry nozzle as a function of
cyclotron beam current.

In order to deliver different dose levels with a wide range of dose rates we
used the previously described monitor 1 calibrated against dose rate
independent Faraday cup. For all cyclotron beam currents up to 5 nA we did not
observe any efficiency drop of monitor 1 due to ion recombination and the
ratio between monitor units (MU) of monitor 1 and the total charge $q$
measured with the Faraday cup remained constant. Also between 5 and 10 nA the
efficiency drop was very subtle, always below 1 %. For cyclotron currents
exceeding 10 nA, the monitor efficiency, normalized to the MU$/Q$ ratio
corresponding to a cyclotron current of 1 nA, was observed to drop rapidly
with the increasing current, as shown in Fig. 6. The corresponding calibration
curve had to be divided in two regions - moderate current range (10-300 nA)
and high current range (300-800 nA), as the monitor response presents a kink
at about 300 nA and another polynomial had to be fitted to the latter
cyclotron current range. With this calibration we were able to deliver
different doses with different dose rates with a precision always better than
1 %, as verified with the Faraday cup.

####  III.2.2 Pencil beam scanning mode

The range shifter enables scanning in depth between 19.6 g$\cdot$cm-2 and 37.9
g$\cdot$cm-2 (R80 ranges of the proton beam at 80 % distal fall-off). The
measured ranges (Fig. 7) correspond to the energies 250 MeV and 170 MeV,
respectively. The former is achieved when no range shifter (RS) plate is used,
the latter corresponds to the use of all the 40 RS plates.

Figure 7: Bragg peak curves measured with a range scanner for two minimum and
maximum beam energies. Figure 8: Transverse pencil beam scanning pattern
acquired with a CCD-scintillator detector.

The maximum range of the transverse scanning with the minimum beam spot
distortion was found to be 16 cm ($\pm 8$ cm) and 1.2 cm ($+0.5$/$-0.7$ cm) in
the U and T directions, respectively. In Fig. 8 a scanning pattern acquired
with a CCD-scintillator detector is shown. The pattern is based on a
previously measured 2D scanning map. The spots are separated by 2 cm in the U
direction and by 1.2 cm (maximum range) in the T direction.

##  IV Discussion

In this paper, we report the commissioning of the former treatment unit Gantry
1 at PSI for FLASH experiments. In order to maximize the dose rate in the
treatment room, we optimized the beam transport for 250 MeV beam, and decided
to use only the in-room range shifter for the energy degradation. We achieved
a transmission of $86\%$, i.e. we could reach more than 680 nA beam
irradiation at isocenter. Higher beam intensity could be potentially be
achieved, as already shown in a proof-of-principle study by Busold et al
Busold2017 . Moreover, losses have been observed only upstream of the
treatment room, and no loss has been observed on the gantry, therefore
satisfying the radiation protection requirements of our facility (where the
shielding of the treatment room is thinner than the one of the cyclotron
bunker, which can sustain higher beam losses). In order to minimize the losses
we refrained from using collimators which help define better the phase space
and make the beam transport less sensitive to cyclotron phase space
variations. This contributes to our beam size uncertainty and thus the total
uncertainty of the dose rate.

Such high beam currents can then be shaped to a dose and dose-rate
distribution using either multiple- or single-spot deliveries. The latter
provides the highest instantaneous dose rates available at our facility: we
could achieve more than 9000 Gy/s, which matches dose-rate levels achieved at
synchrotron facilities such as the European Synchrotron Research Facility
(ESRF) - 8000 Gy/s and 16000 Gy/s Serduc2006 ; Braeuer-Krisch2015  \- and
electron machines like the eRT6 electron linac at Lausanne University Hospital
(CHUV). The gantry control system does not limit the length of these single-
spot irradiations, therefore we can deliver all possible doses in both single-
and multiple-spots deliveries. The minimum pause between two spots corresponds
to 2 ms, but larger pauses can be introduced by the control system, to help
investigating the impact of these variables on irradiation modalities such as
clinical spot scanning. This minimum time can be further decreased down to 100
$\mu$s by using the vertical deflector in the cyclotron and following some
modifications of the control system.

The dose can be precisely delivered based on the nozzle monitor, which we
calibrated against the Faraday cup. However, we observed some day-to-day
instabilities in the calibration. Therefore, we apply a day-specific
correction for our calibration by measuring the monitor-1-to-Faraday-cup ratio
(MU$/Q$) for low and high dose rates, which takes a few minutes. We observed
an increase ratio over time since the original calibration in the MU$/Q$
ratio. We have not understood this effect yet. One of the explanations could
be temperature and pressure changes, for which we do not apply any correction.
The maximum deviations we have observed so far were up to 4 % and can be
completely neutralized by the mentioned day-specific correction.

Even though the energy can be modified easily inside the treatment room, the
presence of the range shifter at the end of the beamline causes a clear
correlation between final energy and dose rate. The beam size increases as
function of number of range shifter plates and distance from the gantry exit,
and consequently the dose rate drops. At the minimum energy achievable with
our system, 700 Gy/s could be reached at the isocenter. Such performance could
be improved in case a different material for the range shifters is used, or in
case the samples are placed closer to the gantry exit.

On the other hand, the relationship between range shifters and beam size
allows an easy way of adapting the beam conditions to the size of the sample
to irradiate, without the production of complex single or double-scattering
systems. We have tested different combinations of number of range shifter
plates and distances from the gantry nozzle, simulating different sample
sizes. This makes our experimental setup particularly interesting for
experiments with cells, as we could adapt flexibly to the requirements of
external users.

Thanks to its flexibility, PSI Gantry 1 will give a chance to address a few
important questions concerning the FLASH effect mechanism and to define
appropriate dose delivery conditions. Different studies of the FLASH effect
showed different results, from no protective effect of FLASH irradiations to a
significant reduction of toxicity to healthy tissues. However, there are
several aspects which have to be taken into account while comparing various
results obtained with different beam parameters. One of the most important is
the difference between the average and instantaneous dose rate. Depending on
the machine type and dose delivery technique, the two dose rates may be very
different. For instance, the Kinetron and Oriatron eRT6 electron linacs, used
at Orsay and Lausanne University Hospital (CHUV), respectively, allow an
instantaneous dose rate in a 1-2 $\mu$s pulse to exceed 107 Gy/s Favaudon2014
; Lansonneur2019 ; Jaccard2018 ; Petersson2017 . However, the average dose
rate under typical treatment conditions would be ”only” up to a few kGy/s and
would involve pauses between pulses of at least 5 ms. With these machines the
most promising results have been achieved so far. However, if the
instantaneous dose rate plays the key role in triggering of the FLASH effect,
it could well be that reaching an average dose rate of the order of 100 Gy/s
or even 1000 Gy/s will not be sufficient unless large enough instantaneous
dose rate is provided. On the other hand, studies in different animal and
cells models have shown FLASH effects at such average dose rates. Dose-rate
thresholds observed are not fully consistent among different studies, though,
and this prompted some authors to suggest that other factors might play a role
- not only the instantaneous dose rate, but also the total irradiation time as
well as the minimum dose within a (macro-)pulse Bourhis2019b . Different
machines are currently able to test at best only a part of the parameter
space, as highlighted in a recent overview of the so-called FLASH-compatible
machines Esplen2020 ; therefore it is of fundamental importance to develop
flexible irradiation facilities, allowing to test in the same conditions and
different animal models different combinations of dose, dose rate and pulse
structure.

##  V Conclusions and outlook

We have successfully converted the PSI pencil beam scanning Gantry 1 into a
test bench for FLASH experiments. We are able to deliver a wide range of dose
rates from 1 to 9000 Gy/s, which enables detailed studies of FLASH with
protons. Moreover, we can conduct transmission as well as conformal
irradiations using only local energy modulation. Field sizes larger than 5x5
mm2 can be irradiated by means of fast transverse pencil beam scanning.
Although the gantry has not been designed to provide scanning in both
transverse directions, we have proven that it is possible to include the
second direction for fast scanning of small fields. This feature can be used
for pre-clinical studies of FLASH with small animals.

Teaming with colleagues from CHUV and Varian Medical Systems, we are currently
carrying out our first biological experiments at Gantry 1.

###### Acknowledgements.

We thank our collaborators at PSI (M. Schippers, M. Eichin, U. Rechsteiner, B.
Rohrer, M. Egloff) and CHUV (M.-C. Vozenin, V. Grilj, C. Bailat). This work is
partially funded by the Swiss National Science Foundation (grant No. 190663).

## Conflict of interest statement

The authors have no conflict to disclose.

## References

  * (1) V. Favaudon, L. Caplier, V. Monceau, F. Pouzoulet, M. Sayarath, C. Fouillade, M.-F. Poupon, I. Brito, P. Hupé, J. Bourhis, J. Hall, J.-J. Fontaine, and M.-C. Vozenin,  Ultrahigh dose-rate FLASH irradiation increases the differential response between normal and tumor tissue in mice,  Science Translational Medicine 6, 245ra93 (2014).
  * (2) P. Montay-Gruel, K. Petersson, M. Jaccard, G. Boivin, J.-F. Germond, B. Petit, R. Doenlen, V. Favaudon, F. Bochud, C. Bailat, J. Bourhis, and M.-C. Vozenin,  Irradiation in a flash: Unique sparing of memory in mice after whole brain irradiation with dose rates above 100Gy/s,  Radiotherapy and Oncology 124, 365 – 369 (2017),  15th International Wolfsberg Meeting 2017.
  * (3) M.-C. Vozenin, P. De Fornel, K. Petersson, V. Favaudon, M. Jaccard, J.-F. Germond, B. Petit, M. Burki, G. Ferrand, D. Patin, H. Bouchaab, M. Ozsahin, F. Bochud, C. Bailat, P. Devauchelle, and J. Bourhis,  The Advantage of FLASH Radiotherapy Confirmed in Mini-pig and Cat-cancer Patients,  25, 35–42 (2019).
  * (4) J. Bourhis, W. J. Sozzi, P. G. Jorge, O. Gaide, C. Bailat, F. Duclos, D. Patin, M. Ozsahin, F. Bochud, J.-F. Germond, R. Moeckli, and M.-C. Vozenin,  Treatment of a first patient with FLASH-radiotherapy,  Radiotherapy and Oncology 139, 18 – 22 (2019),  FLASH radiotherapy International Workshop.
  * (5) A. Patriarca, C. Fouillade, M. Auger, F. Martin, F. Pouzoulet, C. Nauraye, S. Heinrich, V. Favaudon, S. Meyroneinc, R. Dendale, A. Mazal, P. Poortmans, P. Verrelle, and L. De Marzi,  Experimental Set-up for FLASH Proton Irradiation of Small Animals Using a Clinical System,  International Journal of Radiation Oncology*Biology*Physics 102, 619 – 626 (2018).
  * (6) E. Beyreuther, M. Brand, S. Hans, K. Hideghéty, L. Karsch, E. Leßmann, M. Schürer, E. R. Szabó, and J. Pawelke,  Feasibility of proton FLASH effect tested by zebrafish embryo irradiation,  Radiotherapy and Oncology 139, 46 – 50 (2019),  FLASH radiotherapy International Workshop.
  * (7) M. Buonanno, V. Grilj, and D. J. Brenner,  Biological effects in normal cells exposed to FLASH dose rate protons,  Radiotherapy and Oncology 139, 51 – 55 (2019),  FLASH radiotherapy International Workshop.
  * (8) A. Darafsheh, Y. Hao, T. Zwart, M. Wagner, D. Catanzano, J. F. Williamson, N. Knutson, B. Sun, S. Mutic, and T. Zhao,  Feasibility of proton FLASH irradiation using a synchrocyclotron for preclinical studies,  Medical Physics 47, 4348–4355 (2020).
  * (9) E. S. Diffenderfer et al.,  Design, Implementation, and in Vivo Validation of a Novel Proton FLASH Radiation Therapy System,  International Journal of Radiation Oncology*Biology*Physics 106, 440 – 448 (2020).
  * (10) E. Pedroni, R. Bacher, H. Blattmann, T. Böhringer, A. Coray, A. Lomax, S. Lin, G. Munkel, S. Scheib, U. Schneider, and A. Tourovsky,  The 200-MeV proton therapy project at the Paul Scherrer Institute: Conceptual design and practical realization,  Medical Physics 22, 37–53 (1995).
  * (11) S. Lin, T. Boehringer, A. Coray, M. Grossmann, and E. Pedroni,  More than 10 years experience of beam monitoring with the Gantry 1 spot scanning proton therapy facility at PSI,  Medical Physics 36, 5331–5340 (2009).
  * (12) J. Schippers, R. Dölling, J. Duppich, G. Goitein, M. Jermann, A. Mezger, E. Pedroni, H. Reist, and V. Vrankovic,  The SC cyclotron and beam lines of PSI’s new protontherapy facility PROSCAN,  Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 261, 773 – 776 (2007),  The Application of Accelerators in Research and Industry.
  * (13) K. L. Brown, D. C. Carey, F. C. Iselin, and F. Rothacker,  TRANSPORT: a computer program for designing charged-particle beam-transport systems,  CERN Yellow Reports: Monographs, CERN, Geneva, 1980,  Also publ. as SLAC and FERMILAB.
  * (14) U. Rohrer,  Graphic Transport Framework, 2007.
  * (15) S. Safai and et al,  Dosimetry and biological experiments with proton beams at ultra-high dose rates (FLASH): experimental setup and first dosimetry results,  Physics in Medicine and Biology (2020),  in preparation.
  * (16) C. Gomà, S. Lorentini, D. Meer, and S. Safai,  Proton beam monitor chamber calibration,  Physics in Medicine and Biology 59, 4961–4971 (2014).
  * (17) S. Busold and H. Röcken,  Beam Intensity Modulation Capabilities for Varian’s ProBeam® Isochronous Cyclotron,  Proceedings of the 21st Int. Conf. on Cyclotrons and their Applications Cyclotrons2016 , THP05 (2017).
  * (18) R. Serduc, P. Vérant, J.-C. Vial, R. Farion, L. Rocas, C. Rémy, T. Fadlallah, E. Brauer, A. Bravin, J. Laissue, H. Blattmann, and B. van der Sanden,  In vivo two-photon microscopy study of short-term effects of microbeam irradiation on normal mouse brain microvasculature,  International journal of radiation oncology, biology, physics 64, 1519—1527 (2006).
  * (19) E. Bräuer-Krisch et al.,  Medical physics aspects of the synchrotron radiation therapies: Microbeam radiation therapy (MRT) and synchrotron stereotactic radiotherapy (SSRT),  Physica medica : PM : an international journal devoted to the applications of physics to medicine and biology : official journal of the Italian Association of Biomedical Physics (AIFB) 31, 568—583 (2015).
  * (20) P. Lansonneur, V. Favaudon, S. Heinrich, C. Fouillade, P. Verrelle, and L. De Marzi,  Simulation and experimental validation of a prototype electron beam linear accelerator for preclinical studies,  Physica Medica 60, 50 – 57 (2019).
  * (21) M. Jaccard, M. Durán, K. Petersson, J. Germond, P. Liger, M.-C. Vozenin, J. Bourhis, F. Bochud, and C. Bailat,  High dose-per-pulse electron beam dosimetry: Commissioning of the Oriatron eRT6 prototype linear accelerator for preclinical use,  Medical Physics 45, 863–874 (2018).
  * (22) K. Petersson, M. Jaccard, J.-F. Germond, T. Buchillier, F. Bochud, J. Bourhis, M.-C. Vozenin, and C. Bailat,  High dose-per-pulse electron beam dosimetry — A model to correct for the ion recombination in the Advanced Markus ionization chamber,  Medical Physics 44, 1157–1167 (2017).
  * (23) J. Bourhis, P. Montay-Gruel, P. Goncalves Jorge, C. Bailat, B. Petit, J. Ollivier, W. Jeanneret-Sozzi, M. Ozsahin, F. Bochud, R. Moeckli, J.-F. Germond, and M.-C. Vozenin,  Clinical translation of FLASH radiotherapy: Why and how?,  Radiotherapy and Oncology 139, 11 – 17 (2019).
  * (24) N. M. Esplen, M. S. Mendonca, and M. Bazalova-Carter,  Physics and biology of ultrahigh dose-rate (FLASH) radiotherapy: a topical review,  Physics in Medicine and Biology (2020).