# SClib, a hack for straightforward embedded C functions in Python

Esteban Fuentes∗†, Hector E. Martinez‡ * Corresponding author:
<EMAIL_ADDRESS>Electric Drives and Power Electronics group at the
Technical University of Munich‡ T30f group at the Physics Department of the
Technical University of MunichCopyright © 2014 Esteban Fuentes et al. This is
an open-access article distributed under the terms of the Creative Commons
Attribution License, which permits unrestricted use, distribution, and
reproduction in any medium, provided the original author and source are
credited. http://creativecommons.org/licenses/by/3.0/

###### Abstract

We present SClib, a simple hack that allows easy and straightforward
evaluation of C functions within Python code, boosting flexibility for better
trade-off between computation power and feature availability, such as
visualization and existing computation routines in SciPy.

We also present two cases were SClib has been used.

In the first set of applications we use SClib to write a port to Python of a
Schrödinger equation solver that has been extensively used the literature, the
resulting script presents a speed-up of about 150x with respect to the
original one. A review of the situations where the speeded-up script has been
used is presented. We also describe the solution to the related problem of
solving a set of coupled Schrödinger-like equations where SClib is used to
implement the speed-critical parts of the code. We argue that when using SClib
within IPython we can use NumPy and Matplotlib for the manipulation and
visualization of the solutions in an interactive environment with no
performance compromise.

The second case is an engineering application. We use SClib to evaluate the
control and system derivatives in a feedback control loop for electrical
motors. With this and the integration routines available in SciPy, we can run
simulations of the control loop a la Simulink. The use of C code not only
boosts the speed of the simulations, but also enables to test the exact same
code that we use in the test rig to get experimental results. Again,
integration with IPython gives us the flexibility to analyze and visualize the
data.

###### Index Terms:

embedded C code, particle physics, control engineering

##  1 Introduction

Embedding code written in oder languages is a common theme in the Python
context, the main motivation being speed boosting. Several alternatives exist
to achieve this, such as Cython [Cython], CFFI [CFFI], SWIG [SWIG], weave
[weave], among others. We present yet another alternative, which may be quite
close to CFFI than to the others. The motivation to write SClib grew out of
the urge to integrate C code, which was already written, into the Python
environment, minimizing the intervention of the code. Part of the resulting
work is briefly introduced later, in the engineering application section.

Nevertheless, embedding compiled code in Python will naturally have an impact
in performance, for instance, when the compiled code takes care of computer
intensive numerics. The first application we introduce (in particle physics),
leverages SClib in this sense, outsourcing the numerics to the compiled code
and using the Python environment for visualization.

##  2 SClib

At the core of SClib1 is ctypes [Hell], which actually does the whole work: it
maps Python data to C compatible data and provides a way to call functions in
DLLs or shared libraries. SClib acts as glue: it puts things together for the
user, to provide him with an easy to use
interface.\raisebox{10.00002pt}{\hypertarget{id7}{}}\hyperlink{id5}{1}\raisebox{10.00002pt}{\hypertarget{id7}{}}\hyperlink{id5}{1}footnotetext:
The code for SClib and example use are available at
<https://github.com/drestebon/SClib>

The requirements for SClib are very simple: call a function on an array of
numbers of arbitrary type and size and return the output of the function,
again of arbitrary type and size.

The resulting interface is also very simple: A library is initialized in the
Python side with the path to the DLL (or shared library) and a list with the
names of the functions to be called:

In [1]: import SClib as scIn [2]: lib = sc.Clib(’test.so’, [’fun’])The
functions are then available as members of the library and can be called with
the appropriate number of arguments, which are one dimensional arrays of
numbers. The function returns a list containing the output arrays of the
function:

In [3]: out, = lib.fun([0])In the C counterpart, the function declaration must
be accompanied with specifications of the inputs and outputs lengths and
types. This is accomplished with the helper macros defined in sclib.h:

#include <sclib.h>SCL_OL(fun, 1, 1); /* outputs lengths */SCL_OT(fun, 1, INT);
/* outputs types */SCL_IL(fun, 1, 1); /* inputs lengths */SCL_IT(fun, 1, INT);
/* inputs types */void fun(int * out, int * in) { *out = 42;}An arbitrary
number of inputs or outputs can be specified, for example:

#include <math.h>#include <sclib.h>SCL_OL(fun, 2, 1, 2);SCL_OT(fun, 2, INT,
FLOAT);SCL_IL(fun, 2, 1, 2);SCL_IT(fun, 2, INT, FLOAT);void fun(int * out0,
float * out1, int * in0, float * in1) { *out0 = 42*in0[0]; out1[0] =
in1[0]*in1[1]; out1[1] = powf(in1[0], in1[1]);}In the function declaration,
all the outputs must precede the inputs and must be placed in the same order
as in the SCL macros.

These specifications are processed during compilation time, but only the
number of inputs and outputs is static, the lengths of each component can be
overridden at run time:

In [4]: lib.INPUT_LEN[’fun’] = [10, 1]In [5]: lib.retype()In these use cases
the length of the arguments should be given to the function through an extra
integer argument.

In the function body, both inputs and outputs should be treated as one
dimensional arrays.

##  3 Application in Quarkonium Physics

###  3.1 Motivation

The Schrödinger equation is the fundamental equation for describing non-
relativistic quantum mechanical dynamics. For the applications we will present
in this section we will focus on the time-independent version which, in
natural units, is given by

$\left(-\frac{\nabla_{\mathbf{r}}^{2}}{2\mu}+V(\mathbf{r})\right)\psi(\mathbf{r})=E\psi(\mathbf{r}).$
(1)

It corresponds to an eigenvalue equation where the term inside the parenthesis
in l.h.s. is called the Hamiltonian operator, the value $E$, its eigenvalue,
is the measurable quantity (the energy) associated with it, $\mu$ is the
reduced mass of the system (it correspond the mass of the particle in one-
particle systems) and the wavefunction, $\psi(\mathbf{r})$ is the entity
containing all the information about the system, since its modulus squared
correspond to the probability density of a given measurement, it has to be
normalized to unity. The term $V({\mathbf{r}})$ in the Hamiltonian is called
the potential.

Since its discovery, the Schrödinger equation has played an important role in
our understanding of nature and it is present in almost every aspect of modern
physics. In this section we will review some cases where SClib has been used
to implement solutions of the computing problems associated with eq. (1) that
arise in the study of heavy
quarkonia2.\raisebox{10.00002pt}{\hypertarget{id9}{}}\hyperlink{id8}{2}\raisebox{10.00002pt}{\hypertarget{id9}{}}\hyperlink{id8}{2}footnotetext:
For a comprehensive review of the status and perspectives of the research in
heavy quarkonia we refer the reader to chapter four of [Bra14].

Quarkonium is a bound-state composed by a quark and its corresponding
antiquark. By heavy we mean states composed by charm and bottom quarks, called
charmonium and bottomonium respectively. Due to its large mass, the top quark
decays before forming a bound state. In heavy quarkonium the relative velocity
between the quark and antiquark inside of the bound-system is believed to be
small enough for the system to be considered, at least in a first
approximation, non-relativistic, making it suitable for being described by eq.
(1). Considering the equal mass case with a spherically symmetric potential
the angular part can be neglected (it correspond to the spherical harmonics)
and the relevant part of eq. (1) reduces to the one-dimensional equation given
by

$\left[-\frac{1}{m}\frac{d^{2}}{dr^{2}}+\frac{l(l+1)}{mr^{2}}+V(r)\right]y_{n,l}(r)=E_{n,l}y_{n,l}(r),$
(2)

where $r$ is the relative distance between the quark and the antiquark, $l$ is
the angular momentum quantum number, $m$ is the (anti)quark mass (for equal
mass $2\mu=m$), $y_{n,l}$ is called the reduced wavefunction and the
eigenvalue $E_{n,l}$ is interpreted as the binding energy of the bound-system,
where $n=0,1,2,\dots$ accounts for the number of nodes (radial excitations) of
the wavefunction. The mass of the quarkonium is then given by

$M=2m+E_{n,l}.$ (3)

The potential $V(r)$ describes the quark-antiquark interaction, it is a
function of $r$ and $\Lambda_{\rm QCD}$, the typical hadronic scale ($\sim
200\,{\rm MeV}$). For $r\Lambda_{\rm QCD}\ll 1$ (short-distance regime) the
potential may be evaluated perturbatively, but for $r\Lambda_{\rm QCD}\sim 1$
(long-distance regime) it cannot. To overcome this issue, models based on non-
relativistic reductions of phenomenological observations have been used to
describe heavy quarkonia, one these being the so-called Cornell potential
[Eich74], [Eich78], [Eich79]

$V(r)=\frac{a}{r}+kr,$ (4)

where $a$ and $k$ are parameters which need to be fixed by experimental (or
lattice) data of some observable. This potential incorporates two of the main
observed characteristics of the quark-antiquark interaction: at short
distances it exhibits a Coulombic behavior and in the long-distance regime the
interaction is dominated by a confinement phase.

Since the beginning of the last decade, non-relativistic effective field
theories (EFT), in particular non-relativistic QCD (NRQCD) [Cas85], [Bod94]
and potential NRQCD (pNRQCD) [Bra99], have become the state-of-the-art tools
for the study of heavy quarkonia (for review see [Bra04]). NRQCD is obtained
from QCD integrating out modes that scale like $m$, while pNRQCD is obtained
from NRQCD integrating out modes that scale like the quark
momentum3.\raisebox{10.00002pt}{\hypertarget{id19}{}}\hyperlink{id18}{3}\raisebox{10.00002pt}{\hypertarget{id19}{}}\hyperlink{id18}{3}footnotetext:
These EFT exploit the hierarchy of energy scales present in the bound-system.
If the relative velocity of the (anti)quark, $v$, is small, we have that
$mv^{2}(\sim E)\ll mv(\sim p)\ll m$, where $p$ is the momentum of the
particles and $E$ its kinetic energy. If one is interested in studying a
phenomena that happens at the scale $E$ (like the binding) it is more suitable
to integrate out degrees of freedoms with energies that scale like the other
two higher scales, this is the motivation behind pNRQCD. For a detailed
analysis of the scales present in heavy quarkonia we refer the reader to
[Bra04].

The physics of the modes that have been integrated out is encoded in Wilson
coefficients that must be calculated comparing at the same scale the results
(observables, Green functions) of the EFT, with the ones of QCD (for NRQCD) or
NRQCD (for pNRQCD). A key feature of pNRQCD is that it allows the relativistic
corrections to the quark-antiquark potential to be organized as an expansion
in powers of $1/m$. Up to second order $V(r)$ can be written as

$V(r)=V^{(0)}(r)+\frac{V^{(1/m)}(r)}{m}+\frac{V^{(1/m^{2})}(r)}{m^{2}},$ (5)

where $V^{(1/m)}$ and $V^{(1/m^{2})}$ are derived from QCD through the
matching procedure with NRQCD. The details about $V^{(1/m)}$ and
$V^{(1/m^{2})}$ and how they are obtained are beyond the scope of this
document, however, we can list some of their features:

  * •

They correspond to correlators that in the short-distance regime can be
computed in perturbation theory.

  * •

In the long-distance regime they can be computed in lattice QCD, however only
some of these correlators have been calculated.

  * •

Eq. (4) correspond, at least qualitatively, to the leading order $V^{(0)}$ in
eq. (5)

For the details about the derivation of the terms present in eq. (5) we refer
the reader to [Bra00] and [Pin00]. It is important to recall that, although it
can not be evaluated analytically in the whole range of $r$, eq. (5)
represents a model-independent expression for the quark-antiquark potential,
contrary to models like the one presented in eq. (4).

Using perturbation theory to include the relativistic corrections to the
potential, the expression for the bound-state mass reads

$\displaystyle M$ $\displaystyle=$ $\displaystyle
2m+E_{n,l}^{(0)}+\frac{\langle nl|V^{(1/m)}(r)|nl\rangle}{m}$ $\displaystyle+$
$\displaystyle\frac{\langle
nl|V^{(1/m^{2})}(r)|nl\rangle}{m^{2}}+\frac{1}{m^{2}}\sum_{m\neq
n}^{\infty}\frac{|\langle
nl|V^{(1/m)}|ml\rangle|^{2}}{E_{n,l}^{(0)}-E_{ml}^{(0)}},$

where $E_{il}^{(0)}$ comes from solving eq. (2) with $V(r)=V^{(0)}(r)$ and

$\langle
nl|f(r)|n^{\prime}l^{\prime}\rangle\propto\int_{0}^{\infty}dr\,y_{n,l}(r)f(r)y_{n^{\prime}l^{\prime}}(r),$

where the proportionality factor will depended on the corresponding quantum
numbers of the operators appearing in $V^{(1/m)}$ and $V^{(1/m^{2})}$.

Figure 1: Reduced wavefunctions $y_{n,l}(r)$ for two steps in the search of
the eigenvalue $E_{n=1\,l=1}$. For each step in the process to find the
eigenvalue the nodes of the $y_{n,l}$ are counted, only when the value of
$E_{n,l}$ corresponds to an eigenvalue $y_{n,l}$ is not divergent. In the plot
(a) $E_{n,l}=3.1\,m$ and $y_{n,l}$ fulfills the condition of having one node,
however, the accuracy in the value of $E_{n,l}$ is too low and the function
diverges. In the plot (b) $E_{n,l}=3.10952\,m$ so $y_{n,l}\rightarrow 0$ for
larger values of $r$. We have used the Cornell potential eq. (4) with
parameters $m=1$ $a=0.1$, $k=0.5m^{2}$, all dimensions defined in terms of the
mass.

###  3.2 Applications of SClib

The simplest computational problem related to eq. (2) is to find $E_{n,l}$ for
a given $n$ and $l$. Methods to solve this problem have been implemented since
long ago (see for instance [Fal85]), in a nutshell, the standard method
consist of applying two known constraints to the reduced wavefunction
$y_{n,l}$:

  * •

The number of nodes of $y_{n,l}(r)$ must be equal to $n$.

  * •

$y_{n,l}(r)$ has to be normalizable

$\int_{0}^{\infty}dr[y_{n,l}(r)]^{2}=1.$ (7)

In general $y_{n,l}(r)$ will diverge except when $E_{n,l}$ corresponds to an
eigenvalue. The procedure to find the eigenvalue consists in to perform a scan
of values of $E_{n,l}$ until $y_{n,l}(r)$ has $n$ nodes and converges for a
large enough value of $r$ (see Fig. 1). This implies that for each test value
of $E_{n,l}$ eq. (2) must be (numerically) solved. A popular4 Mathematica
[Mat9] implementation of this method has been available in [Luc98]. This
script has the advantage that the user can profit from the Mathematica built-
in functions to plot, integrate or store the resulting wavefunctions, however,
it has a very poor performance. With the goal of mimicking some of the
advantages of this script, but without compromising speed, we ported the
algorithm in [Luc98] to Python. The resulting script, SChroe.py5, uses SClib
to implement the speed-critical parts of the algorithm. In Schroe.py the
wavefunctions are stored as NumPy arrays [NumPy] so when the script is run
within IPython [IPy] together with SciPy [SciPy], NumPy and Matplotlib [Mplot]
the user can profit of the same or more flexibility as with the Mathematica
script plus a boosted speed. In table 1 we compare the performance of
SChroe.py against other implementations of the same
algorithm6.\raisebox{10.00002pt}{\hypertarget{id34}{}}\hyperlink{id24}{4}\raisebox{10.00002pt}{\hypertarget{id34}{}}\hyperlink{id24}{4}footnotetext:
The paper describing the script ranks fifth among the most cited papers (91
citations) of the International Journal of Modern Physics C with the last
citation from July 2014.
\raisebox{10.00002pt}{\hypertarget{id35}{}}\hyperlink{id28}{5}\raisebox{10.00002pt}{\hypertarget{id35}{}}\hyperlink{id28}{5}footnotetext:
Code available at <https://github.com/heedmane/schroepy/>
\raisebox{10.00002pt}{\hypertarget{id36}{}}\hyperlink{id33}{6}\raisebox{10.00002pt}{\hypertarget{id36}{}}\hyperlink{id33}{6}footnotetext:
Although the aim of this section is not to compare performance of Schrödinger
equation solvers, but to present an application in which SClib can improve the
speed of a known algorithm, we must mention that there are solvers that offer
better performance than the current version of SChroe.py. For instance, the
solver presented in [dftatom] implements a more sophisticated integration
method and allows refinements in the radial mesh. With these improvements the
dftatom solver can reach a speed-up of at least two orders of magnitude
compared to the current version of SChroe.py.

$n$ |  $E_{n,l=1}\,\,[m]$ |  schroe.nb [Luc98] |  Python |  SChroe.py
---|---|---|---|---
0 |  2.15789 |  98.88 |  25.46 |  0.66
1 |  3.10952 |  124.14 |  30.95 |  0.75
2 |  3.93850 |  135.68 |  35.32 |  0.84
20 |  13.5995 |  370.0 |  88.04 |  1.99
TABLE II: Time in seconds taken to compute the eigenvalues and reduced
wavefunctions for the Cornell potential eq. (4). The column Python correspond
to the implementation of the algorithm in Python without SClib. The parameters
of the potential are the same as in Fig. 1. All the scripts were tested in the
same machine, a notebook with a 2.4 Ghz core i5 processor (dual core) and 8 GB
of RAM.

In [Bra14] SChroe.py has been used to evaluate the relativistic corrections to
the mass spectrum of quarkonium in the long-distance regime. In that paper the
relativistic corrections $V^{(1/m)}$ and $V^{(1/m^{2})}$ appearing in (3.1)
were evaluated assuming the hypothesis that in the long-distance regime the
interaction between the quark and the antiquark can be described by a string.
In Fig. 2 we show some of the energy levels (masses) corresponding to the
string spectrum. It is noteworthy to mention that all the numerical
calculations and plots of that paper were done with IPython using the SciPy
library.

Figure 2: Long-range energy levels of the first triplet quarkonium state. The
lines are calculated from eq. (3.1) using the relativistic corrections derived
from the string hypothesis [Bra14a]. The leading order (LO) correspond to eq.
(4) setting $a=0$ and $k=1$ (in the plot labeled $\sigma$) and $m=3\sqrt{k}$.
This plot shows the relative size of the next-to-leading-order (NLO)
correction (the term proportional to $1/m$ in the r.h.s. of eq. (3.1)) and the
newly computed next-to-next-to-leading-order (NNLO) corrections (the terms
proportional to $1/m^{2}$). For more details see [Bra14a].

An application in which the speed of SChroe.py plays an important role is
fixing the parameters of the potential given some experimental input. For
instance, consider the problem of finding the parameters $a$ and $k$ of eq.
(4) together with $m$, given the experimental values of the masses of three
different quarkonium states. If relativistic corrections are included, in
order to find the parameters we must solve a system of three equations like
eq. (3.1). For each probe value of $(a,k,m)$ we have to find the eigenvalues
and reduced wavefunctions of eq. (2) and then with these values evaluate the
sums and integrals in (3.1). A parameter fixing of this type was necessary to
implement in [Bra14b]. The implementation has been carried out using SChroe.py
together with a mixture of C and SciPy functions using SClib to link both
environments7.\raisebox{10.00002pt}{\hypertarget{id44}{}}\hyperlink{id43}{7}\raisebox{10.00002pt}{\hypertarget{id44}{}}\hyperlink{id43}{7}footnotetext:
Some of the code will be available once the paper appears online.

Another related computational problem that arises from the study of heavy
quarkonium hybrids, bound-states composed by a quark-antiquark pair plus an
exited gluon, is to solve a system of $N$ Schrödinger-like coupled equations.
Explicitly the system to solve reads

$\left(-\frac{\delta_{ij}}{m}\frac{d^{2}}{dr^{2}}+V_{ij}(r,l)\right)u_{j,(n,l)}(r)=E_{n,l}\,u_{i,(n,l)}(r),$
(8)

where $i=1,2,..N$ and the angular momentum dependence has been included in the
potential matrix. A method to solve this equation for the case $N=2$ has been
implemented in [Ber14]. The method relies on an extension of the nodal theorem
[Ama95] and convergence conditions for the components of the vector
wavefunction $u_{j,(n,l)}(r)$. The extension of the nodal theorem states that
the number of nodes of the determinant of the matrix $U_{n,l}(r)$, whose
columns are $N$ lineal-independent solutions of eq. (8), is equal to $n$. The
procedure then consist in a scan of values $E_{n,l}$; in each step the set of
equations (8) is solved and the nodes of $|U_{n,l}(r)|$ are counted for a
large enough interval of $r$. As in the one-dimensional case, if $E_{n,l}$
approached to an eigenvalue the components of $u_{j,(n,l)}$ converge for large
$r$. In the solution presented in [Ber14] the performance-intensive parts of
the implementation rely on C functions linked to the IPython interface trough
SClib.

As an example of the application of the method implemented in [Ber14], in Fig.
3 we show the results for the search of the first two eigenvalues and
wavefunctions with the matrix potential given by

$V_{ij}(r,l)=\begin{pmatrix}\frac{l(l+1)+2}{mr^{2}}+F_{0}(r)&-\frac{2\sqrt{l(l+1)}}{mr^{2}}\\\
-\frac{2\sqrt{l(l+1)}}{mr^{2}}&\frac{l(l+1)}{mr^{2}}+F_{1}(r)\end{pmatrix}$
(9)

where

$F_{i}(r)=\ln(a_{i}+b_{i}r).$

Figure 3: Solutions for the components of the vector wavefunction $u_{n,l}(r)$
for the first two eigenvalues ($l=1$) of eq. (8) with the matrix potential
given in (9). We have used $m=1$, $a_{0}=1$, $b_{0}=0.5$, $a_{1}=2$ and
$b_{1}=0.1$. The eigenvalues are $E_{n=0,l=1}=1.01727\,m$ for Fig. (a) and
$E_{n=1,l=1}=1.18789\,m$ for Fig. (b).

In all the applications described in this section the combination of SClib and
the SciPy library within an interactive environment like IPython provided a
powerful framework based entirely on open source software for solving problems
that require a high performance and visualization tools.

##  4 Application in Control Engineering

Figure 4: General scheme of a control system.

Most control systems have the structure depicted in Fig 4. $G$ is the plant,
it represent the natural phenomena we wish to control. We usually describe it
using ordinary differential equations:

$G:\;\left\\{\begin{array}[]{rl}\frac{dx}{dt}&=f(x,u,d)\\\
y&=c(x,u,d).\end{array}\right.$ (10)

$x$ represents the internal state of the plant and $y$ its output (the
measurements). $d$ is an independent variable, usually not measurable, named
the perturbation and $u$ is the actuation: the degree of freedom used by the
controller $C$ to achieve the control goal $r$. In general the controller is a
function of the measurements and the reference $r$:

$C:\;u=\pi(y,r),$

but it also may comprise internal states. They are commonly used to
reconstruct the state $x$ out of the history of $y$ and $u$. The latter
systems are called state observers and the whole is called feedback control.

We use SClib to put together a simulator for these kind of systems. Both the
system derivatives $f(\cdot)$ and the control $\pi(\cdot)$ are written in C
and are evaluated using SClib. As stated before, the system state represents a
natural phenomena, therefore it is natural to describe it as a continuous time
variable, as eq. (10) suggests. To calculate the system state we have to solve
this equation. In our simulator this is achieved using numerical methods,
namely the integration routines available in _scipy.integrate_. On the other
hand, the controller is usually implemented in a real-time computer, which can
only sample $y$ at a fixed interval (called $h$): it is a discrete-time
system. This means, that the simulator only needs to evaluate $\pi(\cdot)$ at
given times.

Traditional controllers took the form of linear filters, which could even be
implemented using analog circuitry. As control techniques and requirements
advance, more complex controllers are devised. Many modern control techniques
are based on optimization methods. Time-optimal controllers, for example,
require the solution of an usually very complex optimization problem, to find
a control $u$ that leads the system state $x$ towards its target $r$ in
minimum time [Gru11]:

$u^{*}=\pi^{*}(x)=\underset{\pi\in U,\,x\in
X}{\operatorname{argmin}}\left\\{T_{x}(u)\right\\}.$ (11)

Here $T_{x}(u)$ is the time required to lead $x$ towards its target and $X$
and $U$ are the regions where we want $x$ and $u$ to be confined, they
constitute the constraints for the control problem. These kind of controllers
require exhaustive computation and it is natural to implement them in C.

For motivation, we present the results for a minimum-time control strategy for
a relatively simple and well known problem, the double integrator [Fu13]:

$\frac{d}{dt}\left[\begin{array}[]{c}x_{0}\\\
x_{1}\end{array}\right]=\left(\begin{array}[]{c}\nicefrac{{u}}{{\tau_{0}}}\\\
\nicefrac{{x_{0}}}{{\tau_{1}}}\end{array}\right).$ (12)

The relevance of this system lays in that it models many mechanical systems:
$u$, $x_{0}$ and $x_{1}$ may represent acceleration, speed and position, for
example.

Fig. 5 presents a minimum time control strategy for this system.

Figure 5: Time optimal control for the double integrator considering
$\tau_{0}=\tau_{1}=5$, $u\in[-1,1]$, $h=1$ and $x\in[-1,1]\times\mathbb{R}$.

The form of $\pi(x)$ for this case reveals its non-linear nature.

Fig. 6 presents the trajectory developed by the state using this control
strategy and random initial conditions.

Figure 6: Time optimal trajectories for the double integrator, with random
initial conditions.

These results were obtained using SClib and the devised simulator. The example
code provided with SClib is ready to reproduce them.

The main advantage we obtained from this work was that, since we were using a
Linux based real time system in our test rig, we could use exactly the same
code for the simulations and the experimental tests. Another feature of this
work is that it effectively replaces Simulink in all of our use cases using
only free software.

##  5 Final Remark

We hope the applications of SClib scope beyond the ones listed in this paper
since we believe it provides a simple but powerful way to boost Python
performance.

##  6 Acknowledgments

H.M. acknowledges financial support from DAAD and the TUM Graduate School
during the realization of this work.

## References

  * [Cython] Stefan Behnel, Robert Bradshaw, Lisandro Dalcín, Mark Florisson, Vitja Makarov, Dag Sverre Seljebotn. _Cython is an optimising static compiler for both the Python programming language and the extended Cython programming language (based on Pyrex)._ http://www.cython.org/
  * [CFFI] Armin Rigo and Maciej Fijalkowski _C Foreign Function Interface for Python._ http://cffi.readthedocs.org/en/release-0.8/
  * [SWIG] SWIG developers. _SWIG: Simplified Wrapper and Interface Generator_ http://swig.org/
  * [weave] Ralf Gommers. _Weave provides tools for including C/C++ code within Python code._ https://github.com/scipy/weave
  * [Hell] Heller. _The ctypes module._ , https://docs.python.org/3.4/library/ctypes.html#module-ctypes
  * [Bra14] N. Brambilla, S. Eidelman, P. Foka, S.Gardner, A. S. Kronfeld, M. G. Alford, R. Alkofer and M. Butenschoen et al., _QCD and Strongly Coupled Gauge Theories: Challenges and Perspectives_ arXiv:1404.3723
  * [Eich74] E. Eichten, K. Gottfried, T. Kinoshita, J. B. Kogut, K. D. Lane and T. M. Yan, _The Spectrum of Charmonium,_ Phys. Rev. Lett. 34 , 369 (1975) Erratum-ibid. 36 , 1276 (1976)
  * [Eich78] E. Eichten, K. Gottfried, T. Kinoshita, K. D. Lane and T. M. Yan, _Charmonium: The Model,_ Phys. Rev. D 17 , 3090 (1978) Erratum-ibid. D 21 , 313 (1980)
  * [Eich79] E. Eichten, K. Gottfried, T. Kinoshita, K. D. Lane and T. M. Yan, _Charmonium: Comparison with Experiment,_ Phys. Rev. D 21 , 203 (1980).
  * [Cas85] W. E. Caswell and G. P. Lepage, _Effective Lagrangians for Bound State Problems in QED, QCD, and Other Field Theories,_ Phys. Lett. B 167 , 437 (1986).
  * [Bod94] G. T. Bodwin, E. Braaten and G. P. Lepage, _Rigorous QCD analysis of inclusive annihilation and production of heavy quarkonium,_ Phys. Rev. D 51 , 1125 (1995) Erratum-ibid. D 55 , 5853 (1997)
  * [Bra04] N. Brambilla, A. Pineda, J. Soto and A. Vairo, _Effective field theories for heavy quarkonium,_ Rev. Mod. Phys. 77 , 1423 (2005)
  * [Bra99] N. Brambilla, A. Pineda, J. Soto and A. Vairo, _Potential NRQCD: An Effective theory for heavy quarkonium,_ Nucl. Phys. B 566 , 275 (2000)
  * [Bra00] N. Brambilla, A. Pineda, J. Soto and A. Vairo, _The QCD potential at O(1/m),_ Phys. Rev. D 63 , 014023 (2001)
  * [Pin00] A. Pineda and A. Vairo, _The QCD potential at O(1/m^2): Complete spin dependent and spin independent result,_ Phys. Rev. D 63 , 054007 (2001) Erratum-ibid. D 64 , 039902 (2001)
  * [Fal85] P. Falkensteiner and H. Grosse and F. Schoeberl and P. Hertel Comput. Phys. Comm. 34 , 287 (1985)
  * [Luc98] W. Lucha and F. F. Schoberl, _Solving the Schrödinger equation for bound states with Mathematica 3.0,_ Int. J. Mod. Phys. C 10 , 607 (1999)
  * [dftatom] Čertík, O., Pask, J. E., Vackář, J. (2013). dftatom: A robust and general Schrödinger and Dirac solver for atomic structure calculations. Computer Physics Communications, 184(7), 1777–1791.
  * [Mat9] Wolfram Research, Inc. Mathematica Version 9.0 (2012)
  * [Bra14a] N. Brambilla, M. Groher, H. E. Martinez and A. Vairo, _Effective string theory and the long-range relativistic corrections to the quark-antiquark potential,_ arXiv:1407.7761
  * [Bra14b] N. Brambilla, H. E. Martinez and A. Vairo, TUM-EFT 40/13, In preparation.
  * [SciPy] Eric Jones and Travis Oliphant and Pearu Peterson and others http://www.scipy.org/ (2001–)
  * [NumPy] Stefan van der Walt, S. Chris Colbert and Gaël Varoquaux. The NumPy Array: A Structure for Efficient Numerical Computation, Computing in Science & Engineering, 13 , 22-30 (2011)
  * [Ber14] M. Berwein and H. E. Martinez, TUM-EFT 48/14, In preparation.
  * [Ama95] H. Amann and P. Quittner, _A nodal theorem for coupled systems of Schrödinger equations and the number of bound states,_ Journal of Mathematical Physics 36 , 4553 (1995), doi:10.1063/1.530907.
  * [Mplot] John D. Hunter. Matplotlib: A 2D Graphics Environment, Computing in Science & Engineering, 9 , 90-95 (2007)
  * [IPy] Fernando Perez and Brian E. Granger. IPython: A System for Interactive Scientific Computing, Computing in Science & Engineering, 9 , 21-29 (2007)
  * [Gru11] L. Gruene and J. Pannek, _Nonlinear Model Predictive Control: Theory and Algorithms,_ Springer-Verlag, 2011.
  * [Fu13] E. Fuentes, D. Kalise, J. Rodriguez, and R. Kennel _Cascade-free predictive speed control for electrical drives,_ Industrial Electronics, IEEE Transactions on , vol. PP, no. 99, pp. 1–1, 2013.


# A visualisation for conveying the dynamics of iterative eigenvalue
algorithms over PSD matrices

Ran Gutin
Department of Computer Science, Imperial College London
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We propose a new way of visualising the dynamics of iterative eigenvalue
algorithms such as the QR algorithm, over the important special case of PSD
(positive semi-definite) matrices. Many subtle and important properties of
such algorithms are easily found this way. We believe that this may have
pedagogical value to both students and researchers of numerical linear
algebra. The fixed points of iterative algorithms are obtained visually, and
their stability is analysed intuitively. It becomes clear that what it means
for an iterative eigenvalue algorithm to “converge quickly” is an ambiguous
question, depending on whether eigenvalues or eigenvectors are being sought.
The presentation is likely a novel one, and using it, a theorem about the
dynamics of general iterative eigenvalue algorithms is proved. There is an
accompanying video series, currently hosted on Youtube, that has certain
advantages in terms of fully exploiting the interactivity of the
visualisation.

##  1 Simple iterative eigenvalue algorithms

The (naive) QR algorithm [2, 8, 3] evaluated on matrix $M$ begins with setting
$M_{0}=M$ and then repeating the following two steps until convergence:

  1. 1.

Find the QR decomposition $M_{n}=Q_{n}R_{n}$.

  2. 2.

Set $\text{}M_{n+1}=R_{n}Q_{n}$.

The thing to note is that $M_{n+1}=Q_{n}^{T}M_{n}Q_{n}$. We only concern
ourselves with PSD matrices111 A PSD (positive semi-definite) matrix is a
symmetric $\mathbb{R}$-matrix whose eigenvalues are all non-negative.
Equivalently, it is a matrix $M$ for which $v^{T}Mv\geq 0$ for all vectors
$v$., so the fixed points of the iteration above are all diagonal matrices
with non-negative real entries.

We’ve also investigated a variant of the LR algorithm evaluated on a PSD
matrix $M$ that begins with setting $M_{0}=M$ and then repeating:

  1. 1.

Find the Cholesky decomposition $M_{n}=L_{n}L_{n}^{T}$.

  2. 2.

Set $M_{n+1}=L_{n}^{T}L_{n}$.

Observe that all $M_{n}$ are similar to each other because
$M_{n+1}=L_{n}^{-1}M_{n}L_{n}$. Each $M_{n}$ is also PSD. Therefore by the
spectral theorem, the stronger fact follows that all $M_{n}$ are orthogonally
similar to each other.

Note that while there are more modern and sophisticated versions of the QR
algorithm [3], the original one is still widely taught, and remains simpler
than its more recent improvements. We think there are conceptual advantages to
better understanding the original one.

##  2 Account of visualisation

In this paper, we present a visualisation of the QR algorithm. Our
visualisation makes many properties of the naive QR algorithm, and its various
improvements, intuitive. Our visualisation may be beneficial to students and
researchers. Our visualisation is also qualitative, and using it we can show
that many features of the QR algorithm are in fact universal to iterative
eigenvalue algorithms, and not just specific to the QR algorithm.

The QR algorithm is an iterative eigenvalue algorithm. By iterative, we mean
that it employs a function $f$ from a set to itself, and evaluates it
repeatedly on some starting value $x$ to produce the sequence:
$x,f(x),f(f(x)),f(f(f(x)))\ldots$ until it converges close enough to a fixed
point of $f$. It’s easy to verify that for the QR algorithm the fixed points
are matrices whose eigenvalues are easily found. The issue is that in general,
function iterations have very complicated dynamics. For some functions $f$ and
starting values $x$, the iteration can diverge, or even be chaotic. When
working over $\mathbb{R}$, the usual visualisation that’s employed is
sometimes called the cobweb plot, and using it one can develop a good intution
for the dynamics of some instances of function iteration. We cannot use cobweb
plots to understand the dynamics of iterative eigenvalue algorithms however,
because even in the smallest non-trivial case, the $2\times 2$ case, the
function $f$ maps $\mathbb{R}^{4}\rightarrow\mathbb{R}^{4}$.

We focus only on PSD matrices (positive semi-definite). It can be argued that
this special case is sufficient for computing SVDs (Singular Value
Decompositions), eigendecompositions of symmetric matrices, and
eigendecompositions of orthogonal matrices.222Briefly: For symmetric matrices,
one can perform a shift $M^{\prime}=M+\mu I$ for some easily chosen
$\mu\in\mathbb{R}$ which makes $M^{\prime}$ PSD. For an orthogonal matrix $M$,
we can use the matrix logarithm (or Cayley Transform) as a first step towards
making a matrix $M^{\prime}$ which is PSD. SVD is reducible to
eigendecomposition of symmetric matrices in various ways. The PSD case is
especially amenable to visualisation. We won’t try to justify further our
focus on the PSD case.

The visualisation is based on the one-to-one correspondence between $n\times
n$ PSD matrices, and ellipsoids centred at the origin of $\mathbb{R^{n}}$. The
correspondence is given by $M\mapsto\\{Mv\mid v\in S^{n}\\}$ where $S^{n}$ is
the origin-centred hypersphere in $\mathbb{R^{n}}$. Every PSD matrix
corresponds to a unique ellipsoid, and every ellipsoid (that is origin-
centred) corresponds to a unique PSD matrices. We therefore sometimes talk
about ellipsoids and their semi-axes, instead of about matrices and their
eigen-values/vectors. We go as far as to say “the ellipsoid” instead of “the
PSD matrix” sometimes.333 This all follows from the spectral theorem.

|  |
---|---|---
|  |

Input PSD matrix
---
1 iteration QR
1 iteration LR

Table 1:

Our input ellipse is in blue. Our output ellipse for 1 iteration is in red. We
also visualise the LR algorithm in green for comparison. We visualise the
naive QR algorithm.

We make some immediate observations:

  1. 1.

The algorithm always converges for PSD matrices.

  2. 2.

In every iteration, the angle between the large semi-axes and the $x$-axis
diminishes.

  3. 3.

The large semi-axis is in the same pair of quadrants for the output ellipse as
for the input ellipse.

  4. 4.

The fixed points only occur when the semi-axes are aligned with the coordinate
axes. For a given ellipse, this can happen two ways.

  5. 5.

If the large semi-axis is aligned with the $x$-axis, then this fixed point is
stable.

  6. 6.

If the large semi-axis is aligned with the $y$-axis, then this fixed point is
unstable. Being close to this fixed point is undesirable, because 1 iteration
moves you away from it, into the direction of the other fixed point. The
closer you are to the unstable fixed point, the more slowly you move away from
it.

  7. 7.

The rotation of the ellipse gets slower the closer it is to being a circle.
This presents unsolvable difficulties for finding eigenvectors in general as
opposed to eigenvalues.

  8. 8.

Opposite to point 7, the rotation of the ellipse gets faster the closer it is
to being a degenerate line segment. This can be harnessed to speed up the
algorithm.

  9. 9.

If the ellipse $M$ gets scaled to $M^{\prime}=\lambda M$, then the angle of
rotation remains the same.

  10. 10.

The input ellipse is congruent (in the sense of Euclidean geometry and not of
linear algebra) to the output ellipse.

Observation 7 has profound origins. It is a consequence of eigenvector
instability. When the geometric multiplicity of some eigenvalue $\lambda$ of a
matrix $M$ is 2 or more, then an infinitesimal perturbation of $M$ can
violently change the eigenspaces. Therefore, the slowness of rotation in the
near-circular case is a manifestation of an unsolvable difficulty which is
intrinsic to the eigendecomposition problem. Changing to a different algorithm
won’t make it go away.

Further discussion of observation 7 requires us to know what a nearly diagonal
matrix looks like in the ellipse model:

  1. i.

A near-circle.

  2. ii.

An ellipse whose semi-axes are nearly coordinate-axis aligned.

These two scenarios (i and ii) have little overlap. We can find eigenvalues in
both scenarios (because of the Gershgorin circle theorem), but not
eigenvectors, for which we require scenario ii. The eigenvector problem is
uncomputable in general.

This uncomputability is perhaps not a “big deal”. It often suffices to find
the eigendecomposition of a matrix close by to the input matrix, even though
this may violently change the eigenspaces. A discussion of why this can ever
be acceptable is left to a footnote.444 To understand why, we recall the
distinction between forwards numerical stability and backwards numerical
stability, and the difference in applications of each. An algorithm
$\widehat{f}$ for computing a function $f$ is forwards-stable if
$\widehat{f}(x)\approx f(x)$ for all $x$. An algorithm $\widehat{f}$ for
computing a function $f$ is backwards stable if given $x$ there is an
$\widehat{x}\approx x$ such that $\widehat{f}(x)\approx f(\widehat{x})$.
Clearly, forwards stability is desirable, but we’ve shown that for the
eigendecomposition problem it is unattainable. Backwards stability remains
attainable, but what’s the point of it? We explain. Let $f:R{\rightarrow}R$ be
some continuous function. Eigendecomposition enables us to extend $f$ to the
set of PSD matrices by assuming that $f(PDP^{-1})=Pf(D)P^{-1}$ is true. We can
only compute eigendecompositions in a backwards stable way, which means that
we can make $M\approx PDP^{-1}$ true but not $M=PDP^{-1}$. This is fine
though, because by the continuity of $f$ we have that $f(M)\approx
Pf(D)P^{-1}$, so the approximation error from computing $f$ via the
eigendecomposition of $M$ can be made negligible. Be careful though, because
while $D$ may be close to the eigenvalue matrix of $M$ (because eigenspectra
vary continuously) the same is not true for $P$, which may be far from any
true eigenvector matrix of $M$.

The fix for observation 6 is to make the function being iterated
discontinuous. This is obvious in 2 dimensions thanks to the visualisation,
but it is also true in all dimensions. A more general version of this is the
statement of theorem 1. A continuous choice of function $f$ might be nicer,
but continuity can’t hold everywhere.555For the sake of conceptual
understanding, the function $f$ may be engineered to be continuous over a
subspace of its domain over which it has a globally attractive fixed point,
while the discontinuities outside of this subspace may serve to “kick” the
input into the well-behaved subspace. The Wilkinson shift happens to be
discontinuous, so it provides the needed discontinuities.

Observation 8 appears to be a (kind of) converse to observation 7, but it’s
difficult to verify under what conditions an iterative eigenvalue algorithm
behaves this way. In the case of the QR and LR algorithms, this behaviour may
be directly verified, but it would be interesting to demonstrate that it is
true for a large class of other algorithms.

|  |
---|---|---

Input PSD matrix
---
1 iteration QR

Table 2: $3\times 3$ PSD matrices. Observe that when pancaked, the ellipsoid
aligns 1 semi-axis in 1 iteration.

Observation 8 can be combined with observation 9 to motivate shifting, and to
demonstrate the limits of shifting and how they may be bypassed through the
use of deflation. Shifting is the map $M^{\prime}=M+\mu I$ for some arbitrary
$\mu\in\mathbb{R}$. The effect this has is to increase the lengths of all the
semi-axes by $\mu$. To better understand this effect, we exploit observation 9
to scale down the ellipse so that its largest semi-axis has length exactly
equal to $1$. We then observe that the effect of shifting is that it either
fattens or thins the ellipse. The fattest case is a perfect circle, which by
observation 7 is the worst-case scenario. The thinnest case is a line segment,
for which convergence is instant. In higher dimensions, the extreme cases are
flat pancakes or perfect hyperspheres. Therefore, shifting is pancaking!
Pancakes are, of course, the best; and blobs are, of course, the worst. But
once you’ve got a perfect pancake, shifting runs out of steam as a speeding-up
technique, so we may employ deflation to destroy the dimension of the
ellipsoid that has been reduced to $0$, so that we may pancake the remaining
dimensions, and rekindle the speed-up.

The above analysis of observation 8 shows that iterative eigenvalue algorithms
in which observation 8 holds can be used as improvers of cruder but faster
eigenvalue estimation technique. If an eigenvalue estimation technique can
crudely estimate the smallest eigenvalue, then this can be used as the value
of the shift $\mu$, which can be used then to pancake the ellipsoid
effectively, and which in turn will accelerate the iterative eigenvalue
algorithm, and which may in turn make the crude estimation technique more
accurate in subsequent iterations, and so on.

##  3 Account of generalisation to $n$ dimensions

It may not be obvious that the conclusions generalise straightforwardly to PSD
matrices in $n$ dimensions, but they do:

  * $\bullet$

The algorithm converges for all PSD matrices.

  * $\bullet$

The stable fixed points are precisely the diagonal matrices of the form
$\operatorname{diag}(\lambda_{n},\lambda_{n-1},\ldots,\lambda_{1})$ where all
the $\lambda_{i}$’s are non-negative reals, and
$\lambda_{n}\geq\lambda_{n-1}\geq\ldots\geq\lambda_{1}$. The other diagonal
matrices are unstable fixed points, with convergence slowing down the smaller
the angle an ellipsoid makes with them.

  * $\bullet$

The output of each iteration is orthogonally similar (as a matrix) to its
input. In terms of Euclidean geometry, it is a congruent ellipsoid. (Be aware
that “congruence” in Euclidean geometry means something different from in
linear algebra).

  * $\bullet$

At eigenvalue clashes, the eigenvectors are unstable, but the algorithm may
still converge quickly to a nearly diagonal matrix. In which case, only the
eigenvalues and not the eigenvectors may be obtained under such circumstances.

  * $\bullet$

Shifting can be understood as “pancaking” the ellipsoid (once the semi-axes
are scaled so that the largest semi-axis has length $1$). Like in the 2D case,
this speeds up convergence.

  * $\bullet$

Deflation can be thought of as continuing the “pancaking” (shifting) once one
of the semi-axes has been fully pancaked.

These conclusions can be checked by doing the ellipsoid visualisation in 3
dimensions. But this isn’t necessary, as the 2D case is suggestive enough.

##  4 Remark about the dual purpose of Wilkinson shifts, and the search for
better shifting strategies

Wilkinson shifts are curiously enough not continuous. Wilkinson shifts appear
to both “kick” the ellipsoid away from the unstable fixed point (if it’s too
close to it) and to “pancake” the ellipsoid. Both of these features speed up
convergence. One might wonder though if these two properties: The kicking away
from the unstable fixed point, and the pancaking of the ellipsoid, can be
achieved separately. In particular, the shift $\mu$ can perhaps be made into a
continuous function, with it serving only to “pancake” the ellipsoid. This
might simplify the search for superior shifting strategies if we can narrow
the search down to only continuous functions, with the discontinuity needed to
do the “kicking” from the fixed point coming from elsewhere.

##  5 Axiomatic development of the theory of iterative eigenvalue algorithms

One advantage of the visualisation is that it is mostly qualitative. In
particular, we see that the LR algorithm exhibits the same qualitative
behaviour as the QR algorithm. We exploit this to show that much of the
behaviour and theory of the QR algorithm is not caused by the QR decomposition
as such, but is intrinsic to solving the eigenvalue estimation problem by the
use of fixed-point iteration.

Consider the following axioms for a fixed-point iteration algorithm:

  1. 1.

It consists of iterating some function $f$ over PSD matrices.

  2. 2.

$f(M)$ is orthogonally similar to $M$.

  3. 3.

The sequence $(f^{n}(M))_{n\in\mathbb{N}}$ converges for every $M$. It
furthermore converges to a fixed point of $f$.

  4. 4.

A fixed point can only be a diagonal matrix.

  5. 5.

All fixed points are attractive. (This makes convergence fast).

We consider these 5 axioms to be highly desirable. Intriguingly, these rule
out the possibility that $f$ can be continuous everywhere. This is the
statement of theorem 1. Morally, the theorem holds because given a non-scalar
matrix $M$, the topological space of matrices orthogonally similar to $M$ is
connected, but has some “holes” in it. The influence these “holes” have is
that they prevent axioms 1 to 5 being realised when $f$ is continuous, even
though these 5 conditions are highly desirable. A work-around might be to make
$f$ continuous over some neighbourhood of each fixed point, and use the
discontinuity of $f$ outside these neighbours to make $f(x)$ only map into
these neighbourhoods.

###### Theorem 1

Assuming axioms 1 to 5, the function $f$ being iterated must have
discontinuities in the set of PSD matrices.

Proof  We restrict the domain and codomain of $f$ to the set of matrices
orthogonally similar to some PSD matrix $M$, where $M$ is not a multiple of
the identity matrix. The fact that we can do this follows from condition 2. We
show that $f$ has at least two fixed points under continuity: Assume it only
has one fixed point $x$. It must be attractive by condition 5. By condition 3,
all points attract to $x$. This results in the space of matrices orthogonally
similar to $M$ being a contractible space, which it isn’t. We get a
contradiction. Therefore $f$ has at least two fixed points orthogonally
similar to $M$.

Now we show that under continuity, one of these two fixed points is not
attractive, which contradicts conditions 3 and 5. The space of matrices
orthogonally similar to some arbitrary matrix $M$ is connected. The set of
points which attract to some attractive fixed point is always open. The basins
of attraction of the attractive fixed points are disjoint, so by topological
connectivity there are points not in their union, which are the points which
don’t attract to an attractive fixed point. This contradicts conditions 3 and
5. $\Box$

##  6 Future work regarding visualisations

We ask the following questions. Forgive us if some happen to be easy:

  * •

Can a visualisation technique for iterative eigenvalue algorithms be developed
for non-symmetric matrices? The dynamics of the QR algorithm in the case of
non-symmetric matrices (especially featuring shift policies) is not as well
understood as in the symmetric case [1].

  * •

Can a visualisation technique help to better understand Hessenberg matrices
and their uses for finding eigenvalues and eigenvectors? We haven’t succeeded
yet in producing anything useful for this. Even symmetric Hessenberg matrices
would be interesting.

It’s important that visualisations be legible. Whether a visualisation is
legible or not is perhaps subjective.

##  7 Novelty

We first put a visualisation of the QR algorithm onto the Wikipedia page on
the QR algorithm on August 2021. We did so under a pseudonym. Afterwards, we
left it for some time. We then later presented some of the material in this
paper in a Youtube video series [4], also under a pseudonym. Based on a
private correspondence with Professor David Watkins [5] (who has written
numerous pedagogical papers on iterative eigenvalue algorithms like the QR
algorithm [6, 7]) the presentation is likely a novel one. The emphasis is
likely novel. Videos and lectures are an important accompaniment when trying
to present material like the one here.

## References

  * [1] Jess Banks, Jorge Garza-Vargas, and Nikhil Srivastava.  Global convergence of hessenberg shifted qr i: dynamics.  ArXiv preprint arXiv:2111.07976, 2021\.
  * [2] John GF Francis.  The qr transformation a unitary analogue to the lr transformation—part 1.  The Computer Journal, 4(3):265–271, 1961.
  * [3] Gene H. Golub and Charles F. van Loan.  Matrix Computations.  JHU Press, Fourth edition, 2013.
  * [4] Ran Gutin.  Visual qr algorithm tutorial.  Youtube video series.
  * [5] Ran Gutin and David Watkins.  Private correspondence.
  * [6] David S. Watkins.  Some perspectives on the eigenvalue problem.  SIAM Review, 35(3):430–471, 1993.
  * [7] David S. Watkins.  The qr algorithm revisited.  SIAM Review, 50(1):133–145, 2008\.
  * [8] James H Wilkinson.  Global convergene of tridiagonal qr algorithm with origin shifts.  Linear Algebra and its Applications, 1(3):409–420, 1968\.


0 Research  Left: view from inside the colon in VR. Right: a colormap view
indicating the areas already examined in blue and not yet seen in red.

# Development of an Immersive Virtual Colonoscopy Viewer
for Colon Growths Diagnosis

João Serras
IST - ULisboa  e-mail<EMAIL_ADDRESS>Anderson Maciel
IST - ULisboa e-mail<EMAIL_ADDRESS>INESC-ID  Soraia
Paulo
IST - ULisboa e-mail<EMAIL_ADDRESS>INESC-ID  Andrew
Duchowski
Clemson University  e-mail<EMAIL_ADDRESS>Regis Kopper
UNC Greensboro  e-mail<EMAIL_ADDRESS>Catarina Moreira
Queensland University of Technology  e-mail<EMAIL_ADDRESS>Joaquim Jorge
IST - ULisboa e-mail<EMAIL_ADDRESS>INESC-ID

###### Abstract

Desktop-based virtual colonoscopy has been proven to be an asset in the
identification of colon anomalies. The process is accurate, although time-
consuming. The use of immersive interfaces for virtual colonoscopy is
incipient and not yet understood. In this work, we present a new design
exploring elements of the VR paradigm to make the immersive analysis more
efficient while still effective. We also plan the conduction of experiments
with experts to assess the multi-factor influences of coverage, duration, and
diagnostic accuracy.

Immersive Virtual ColonoscopyVisualizationVisualization techniques; Immersive
Virtual ColonoscopyVisualizationVisualization design and evaluation methods
Human-centered computingHuman computer interaction (HCI)nteraction
paradigmsVirtual reality; Applied computingLife and medical
sciencesComputational biologyImaging

Introduction

In 2020, Colorectal Cancer (CRC) was ranked third in terms of incidence and
second as the cause of cancer-related deaths [8]. The field of medical imaging
has made significant advancements in recent years, offering new and improved
ways to diagnose and treat various conditions. One such innovation is the
development of CT colonography, also known as virtual colonoscopy. Unlike
other CT-based exams, the diagnosis of colon polyps and diverticles can hardly
be done solely on sliced image volumes. Instead, the internal surface of the
colon is reconstructed in 3D and displayed on the radiologist monitor from the
perspective of an observer at the colon mid-line. This, combined with a very
wide field of view projection, allows for tunnel-like navigation (see fig. 1)
that makes colon growths apparent, permitting experts to measure them.

Virtual reality has been applied to scientific and data visualization over two
decades with exciting results [11, 3]. Several features of VR, such as
presence, natural look around, navigation, and manipulation, among others, are
responsible for a more focused and efficient visualization task when the data
has multiple dimensions.

In this context, we propose the development of an immersive virtual
colonoscopy viewer for the diagnosis of colon growths. Such viewer utilizes
state-of-the-art technology to create a highly realistic and interactive 3D
representation of the colon, allowing doctors to identify and examine growths
in greater detail than ever before. This approach may improve the accuracy
and/or efficiency of colon growth diagnoses while offering a more comfortable
and less invasive alternative to traditional colonoscopy procedures and
current desktop implementations.

The translation of a desktop viewer for an immersive one sounds quite
straightforward. Nevertheless, very few works in the literature approach
diagnostic colon visualization with an immersive interface, and none provide
consistent guidelines for designing this visualization and the related
interactive actions.

In this paper, we present preliminary results from the development of a VR app
that allows an immersed user to enter the virtual colon and systematically
explore its whole length for clinically relevant findings. Our VR app also
benefits from eye tracking as well as a set of VR metaphors with the goal of
making the radiologist’s job at least as accurate and less time-consuming to
help make the procedure more widely available for patients.

Figure 1: Desktop-based virtual colonoscopy in a 2x2 view configuration.

##  1 Background and Related Work

Nowadays, Total Colonoscopy (TC) is the most accurate method to detect
colorectal cancer. Despite being so, it entails many disadvantages and
complications, such as discomfort, need for preparations, risk of rupture, and
many others due to its invasiveness, since it requires a scope to be inserted
through the rectum.

Virtual Colonoscopy (VC) was created to solve some of TC disadvantages.
Through the use of CT scans and desktop programs, radiologists navigate
through a reconstructed colon in a 2-way path: from the rectum to the cecum
and from the cecum to the rectum to avoid missing anomalies in folds. Despite
solving invasiveness, it brings its own disadvantages. Navigation is
restricted to a single view, causing the level of uncertainty regarding the
types of polyps to be higher, which in turn leads to a higher referral rate to
colonoscopy and superior cost. While VC is able to detect most colon
anomalies, it is highly time-consuming, taking much longer than a conventional
total colonoscopy. Moreover, the exam only takes into account the inside of
the colon; extracolonic lesions may still exist. Thus the need for further
exams increases the cost, and on top of that, false positives can occur.

Finally, Immersive Virtual colonoscopy(IVC) was proposed. A 3D model of the
colon is reconstructed from CT scans, similar to VC. However, the resulting 3D
colon surface is viewed within a virtual reality environment. IVC allows
radiologists to freely explore the inside and outside of the colon using
natural head and body motion, better analyze volumetric data with stereopsis
and parallax cues, and still have the same functionalities as VC. Since the
radiologist is not bound to a single path and view, many possibilities are
created. Measures can be taken directly in 3D. Despite all these
possibilities, radiologists are not interested in fancy VR interaction. The
major problem we face in designing a VR interface for IVC is its efficiency.
Thus, the major challenge we face is to reduce the duration of the exam
keeping the fidelity high. Later in the paper, we describe how we designed the
interactions so that the exam can be made more efficiently.

###  1.1 Anomalies

The colon exams described have the same objective: finding and categorizing
anomalies in the colon. Anyone can develop these anomalies that are
distinguished in 5 types:

Adenomatous (tubular adenoma) \- Most common type of colon polyp and mostly
benign. They can be transformed into other cancerous types although this
process can take several years and can be prevented with regular screening;

Serrated \- This type of polyp can become cancerous or not depending on its
size and location. If large, typically flat and found in the upper colon, they
are precancerous. This kind, due to its characteristics, is hard to locate.

Hyperplastic \- If a serrated polyp is small and found in the lower colon,
they are known as Hyperplastic; they are rarely malignant.

Inflammatory \- Occurs mostly in people who have inflammatory bowel. Also
known as pseudopolyps. They are a product of chronic inflammation in the
colon. Inflammatory polyps are benign and often do not turn into cancerous
types.

Villous Adenoma (Tubulovillous Adenoma) \- Commonly sessile (flat), this type
of polyp has a high risk of becoming cancerous. Small polyps can be removed
during a colonoscopy while larger ones may require additional surgery to
complete their removal.

All types of Polyps should be identified but adenomas, in particular, are
found to turn into cancerous polyps more often. When they are flat, they can
be hard to identify. Furthermore, serrated polyps are a strong factor to CRC
as the polyp may be too small and pass undetected, especially when
preparations were not done properly[8].

With the necessities of the radiologists in mind, previous works studied how
VR impacts the diagnosis [12]. Then, a few papers endeavored to improve their
experience using spatially aware displays [4] and immersive VR. Firstly, to
help navigation in VR, previous papers have shown that World in Miniature and
teleport are great ways to reduce time consumption and improve orientation
while inside cave-like environments such as a colon [2, 10]. Secondly, a
comparison of navigation methods through a path inside the colon has shown
that, despite ”Fly-over” having the best accuracy, ”Fly-through” was the most
liked [9].

Regarding the interface and interactions, there are many interaction and
selection techniques available and studied. In this project, we decided to use
ray-based techniques [1, 5]. A similar system was previously developed [7],
but at the time the HMD used did not have enough resolution and received
significant cybersickness complaints.

Another technology we rely on is eye-tracking. There are previous works
regarding eye-tracking tools for medical purposes. Studies have shown that
they bring significant benefits as visualization aids. They have also
suggested and explained metrics in order to evaluate such tools [6].

##  2 IVC Design

In VC radiologists are used to be able to navigate the colon through its
center line at different speeds, bookmark suspected regions, measure regions,
list areas measured, and inspect the CT scans of a suspected region to define
the diagnosis. Sometimes a droplet of mucus is shaped like a polyp. Density
values available in CT images are used to disambiguate.

In our IVC design, we tried to support the same tasks but with VR elements and
metaphors.

The VE contains an infinite horizon landscape where the virtual colon is
floating. Upon entrance into the VE, the user sees the colon from outside in
an anatomical position. They also see the avatars of the two controllers they
are handling, with labels indicating how to operate. With the left hand, the
user can rotate the colon model, and with the right hand, they can point to a
location on the colon where they want to start the analysis. Fig. 2
illustrates this first view.

Figure 2: Virtual environment outside the colon showing the controls.

The radiologist generally navigates from the rectum to the colon due to
training with conventional techniques. So, let us suppose that they pointed
toward the rectum and teleported there. From there, the radiologist navigates
the colon using the fly-through technique through the mid-line. Five different
velocities are chosen naturally with the joystick of the right hand. The
motion direction is the tangent of the mid-line, but the radiologist can look
around and, when looking back, the navigation direction is inverted. Small
natural head and body movements are still allowed, which helps quickly
approach the colon wall to see details or tilt the head sideways to better
grasp depth features.

Several cues help with self-orientation. First, anatomical elements allow the
specialist to recognize the portion of the colon they are seeing. Second, we
provide an arrow object that always points vertically. Finally, a world in
miniature (WIM) depicting an external view of the colon can be toggled and
appears attached to the left hand (fig. 3). It is possible to teleport to
other locations either pointing on the WIM or pointing at the colon itself.
The teleport destination is the closest position in the mid-line path to the
first ray intersection. The line shown at the mid-line path is originally
white but will be painted green at the segments already visited.

Figure 3: World in miniature attached to the hand.

An additional feature, without parallel in the desktop app, is the ability to
compute the time the user spent looking at each point on the colon wall. Using
the center of the viewport as a reference, we project a set of vectors toward
the model and add time to the points intercepted. When eye tracking is
available, which is the case with our prototype, the viewport center approach
is replaced by the respective eye point in the viewport, which increases the
reliability of the times computed. Time is then mapped to colors, and a
heatmap can be consulted by the radiologist if they want to make sure nothing
was left out. The heatmap can be applied either for the WIM (fig. 4) or the
colon model itself (fig. Development of an Immersive Virtual Colonoscopy
Viewer for Colon Growths Diagnosis). Moreover, the view time for each point is
valuable information for future validation of the VR interface and for meta-
analysis of reading patterns.

Figure 4: Coverage heatmap indicating on the world in miniature for how long
each area of the colon was observed.

The interface allows one to point and bookmark regions for later analysis and
categorization. A list of bookmarks is available where a selection teleports
the user to the respective location for quick access. Besides bookmarking,
measurement is possible in 3D. While in desktop apps the measures are taken on
the three planes, we offer a tool to select any two points on the colon and
measure the distance between them (see fig. 5).

Figure 5: Distance measurement tool to size the anomalies.

A final and crucial feature is to allow the analysis of the abnormality
location on the CT scan. Instead of ”splitting the screen”, which is usual in
desktop apps, to show three canonical planes, we propose to use the 3D space
for that. In our current implementation, CT images can be seen on a floating
panel upon pointing to a location (see fig. 6). Moving the pointing ray causes
the images to update for the new location. However, this is a critical feature
that we wish to investigate experimentally. We will consider other approaches,
such as registering the images with the model and constructing arbitrary
planes interactively.

Figure 6: Slice View

##  3 Final Comments

The current development was exposed to experts for informal feedback. Several
redesigns occurred in this iterative process. The next step is to detail a
protocol for formal experimentation. The main goal is to study coverage,
accuracy, and efficiency among conditions. We hypothesize that the immersive
interface allows for a reduction in the time consumption for the exam while
the false-positives and false-negatives rates remain similar to the desktop
approach.

A group of radiologists from more than one hospital will be recruited to
perform an IVC on colons generated from public data with ground truth. They
will perform the colonoscopy with only one run (rectum to cecum), two runs
(rectum to cecum and cecum to rectum), and also on a conventional desktop VC
with two runs (rectum to cecum and cecum to rectum) as it is the usual setup
in clinical settings. In all tests, the area covered, the time consumed, and
the characteristics of the anomalies discovered will be compared.

We recognize that a between-groups protocol is more suitable for this analysis
to avoid memorization among conditions. Still, some of the variables will be
tested within groups due to the restricted number of professional radiologists
with experience in VC. Initially, we plan to expose the participants to all
conditions but with different colons, randomizing the order of colons and
interfaces.

The system will be installed on Pico Neo 3 Pro Eye HMD. Every action will be
recorded for pattern analysis, and their reactions will be taken into account
for discussion over the quantitative results. Action patterns and exam
outcomes will also be crossed with subjective comfort and cybersickness.

###### Acknowledgements.

This work is partly supported by CNPq-Brazil project 311251/2020-0, by
Fundação para a Ciência e Tecnologia (Portuguese Foundation for Science and
Technology) through grants 2022.09212.PTDC (XAVIER), UIDB/50021/2020, and
Carnegie Mellon Portugal grant SFRH/BD/151465/2021 under the auspices of the
UNESCO Chair on AI&XR.

## References

  * [1] F. Argelaguet and C. Andujar.  A survey of 3d object selection techniques for virtual environments.  Computers & Graphics, 37(3):121–136, 2013. doi: 10 . 1016/j . cag . 2012 . 12 . 003
  * [2] D. Coffey, N. Malbraaten, T. B. Le, I. Borazjani, F. Sotiropoulos, A. G. Erdman, and D. F. Keefe.  Interactive slice wim: Navigating and interrogating volume data sets using a multisurface, multitouch vr interface.  IEEE Transactions on Visualization and Computer Graphics, 18(10):1614–1626, 2012. doi: 10 . 1109/TVCG . 2011 . 283
  * [3] A. Fonnet and Y. Prie.  Survey of immersive analytics.  IEEE Transactions on Visualization & Computer Graphics, 27(03):2101–2122, mar 2021. doi: 10 . 1109/TVCG . 2019 . 2929033
  * [4] J. G. Grandi and A. Maciel.  Interactive spatially aware visualization of medical images.  International Journal of Computer Assisted Radiology and Surgery, 5:347–348, 2010.
  * [5] T. Grossman and R. Balakrishnan.  The design and evaluation of selection techniques for 3d volumetric displays.  In Proceedings of the 19th Annual ACM Symposium on User Interface Software and Technology, UIST ’06, p. 3–12. Association for Computing Machinery, New York, NY, USA, 2006. doi: 10 . 1145/1166253 . 1166257
  * [6] E. Helbren.  Eye-tracking the moving medical image: Development and investigation of a novel investigational tool for CT Colonography.  PhD thesis, 04 2019.
  * [7] S. Mirhosseini, I. Gutenko, S. Ojal, J. Marino, and A. Kaufman.  Immersive virtual colonoscopy.  IEEE Transactions on Visualization and Computer Graphics, 25(5):2011–2021, 2019. doi: 10 . 1109/TVCG . 2019 . 2898763
  * [8] U. of Michigan Health.  Colon and rectal polyps.  Last accessed 10 January 2022.
  * [9] S. F. Paulo, D. Medeiros, P. B. Borges, J. Jorge, and D. S. Lopes.  Improving camera travel for immersive colonography.  In 2020 IEEE Conference on Virtual Reality and 3D User Interfaces Abstracts and Workshops (VRW), pp. 748–749, 2020. doi: 10 . 1109/VRW50115 . 2020 . 00225
  * [10] R. Stoakley, M. J. Conway, and R. Pausch.  Virtual Reality on a WIM: Interactive Worlds in Miniature.  In Proceedings of the SIGCHI Conference on Human Factors in Computing Systems, CHI ’95, p. 265–272. ACM Press/Addison-Wesley Publishing Co., USA, 1995. doi: 10 . 1145/223904 . 223938
  * [11] A. van Dam, A. Forsberg, D. Laidlaw, J. LaViola, and R. Simpson.  Immersive vr for scientific visualization: a progress report.  IEEE Computer Graphics and Applications, 20(6):26–52, 2000. doi: 10 . 1109/38 . 888006
  * [12] J. Venson, J. A. Berni, C. E. da Silva Maia, A. M. da Silva, M. C. d’Ornellas, and A. Maciel.  A case-based study with radiologists performing diagnosis tasks in virtual reality.  In InMEDINFO 2017: Precision Healthcare through Informatics: Proceedings of the 16th World Congress on Medical and Health Informatics, vol. 245, p. 244. IOS Press, 2018.