Reconfigurable Intelligent Surface Enabled Spatial Multiplexing with Fully Convolutional Network

Bile Peng, Member, IEEE,
Ke Guan, Senior Member, IEEE,
Danping He, Member, IEEE,
Jan-Aike Termöhlen, Student Member, IEEE,
Cong Sun, Member, IEEE,
Tim Fingscheidt, Senior Member, IEEE,
and Eduard A. Jorswieck, Fellow, IEEE
B. Peng, J. Termöhlen T. Fingscheidt, and E. A. Jorswieck are with Institute for Communications Technology, TU Braunschweig, 38106 Braunschweig, Germany (e-mail: {peng,termoehlen,fingscheidt,jorswieck}@ifn.ing.tu-bs.de).
K. Guan and D. He are with State Key Laboratory of Rail Traffic Control and Safety, Beijing Jiaotong University, Beijing, China, 100044 (e-mail<EMAIL_ADDRESS>hedanping@bjtu.edu.cn).
C. Sun is with Beijing University of Posts and Telecommunications, Beijing, China, 100876 (e-mail: suncong86@bupt.edu.cn).

ris is an emerging technology for future wireless communication systems.
In this work, we consider downlink spatial multiplexing enabled by the ris for wsr maximization.
In the literature, most solutions use alternating gradient-based optimization,
which has moderate performance, high complexity, and limited scalability.
We propose to apply a fcn to solve this problem,
which was originally designed for semantic segmentation and depth estimation of images.
The rectangular shape of the ris and the spatial correlation of channels with adjacent ris antennas due to the short distance between them encourage us to apply it for the ris configuration.
We design a set of channel features that includes both cascaded channels via the ris and the direct channel.
In the bs, the differentiable mmse precoder is used for pretraining
and the wmmse precoder is then applied for fine-tuning, which is nondifferentiable, more complex, but achieves a better performance.
Furthermore, we propose a method to discretize the continuous phase shifts,
which significantly reduces the requirement for hardware complexity with acceptable performance loss.
Evaluation results show that the proposed solution has higher performance, good robustness against channel estimation error and allows for a faster evaluation than the baselines.
Hence it scales better to a large number of antennas,
advancing the ris one step closer to practical deployment.

Reconfigurable intelligent surface, weighted sum-rate maximization, fully convolutional network, WMMSE precoding.

§ INTRODUCTION

The ris is an emerging technology for next-generation wireless communication systems <cit.>.
It comprises many antennas on a surface.
Each antenna has the ability to receive a signal, process it without external power
(i.e., the signal cannot be amplified),
and reflect it.
With cooperation between the antennas,
the ris can realize complicated signal processing.
Due to its simple structure and good compatibility with the other components of wireless communication systems (e.g., precoding in the bs),
the ris is widely believed to be an essential part of the next-generation wireless communication systems and has been widely studied for optimizing
wsr <cit.>,
capacity <cit.>,
energy efficiency <cit.>, reliability <cit.>,
physical layer security <cit.>
and wireless power transfer <cit.>.

Among different applications of the ris, we focus on improving the wsr with spatial multiplexing in this work.
The high data rate is a major challenge for future wireless communication systems,
which can only be overcome by both acquiring broader bandwidth in higher frequencies
and improving the spectrum efficiency.
With a favorable propagation channel,
we can use the same resource block (a resource allocation unit in time and frequency domains) to serve multiple users and hence improve the spectrum efficiency.

Without ris, the propagation channel depends mainly on positions of transmitter and receivers and nearby scatterers,
which we cannot optimize.
If we deploy an ris in the environment,
we can optimize the environment by configuring the ris such that the propagation channel is more favorable for spatial multiplexing.

This is a challenging problem because of its high dimensionality due to the large number of ris antennas and the constraint that ris does not amplify the signal.
In the literature, the mmse precoder <cit.> and the wmmse precoder <cit.> are proposed for precoding in the bs without consideration of the ris,
stochastic successive convex approximation <cit.>,
majorization-maximization <cit.>,
ao <cit.>,
alternating direction method of multipliers (ADMM) <cit.>
are applied to jointly optimize the ris configuration and the bs precoding.
These proposed algorithms have achieved reasonable performance at the cost of high computational effort.
With similar problem formulations,
spatial multiplexing in ris-assisted uplink is considered and a gradient-based solution is proposed <cit.>.
The Riemannian manifold conjugate gradient (RMCG) and the Lagrangian method are applied to configure multiple ris and bs to serve users on the cell edge <cit.>.
A joint precoding scheme with an ao method is proposed for cell-free ris-aided communication <cit.>.
The sum-rate is maximized with the majorization-maximization (MM) algorithm for grouped users in a ris-aided communication system <cit.>.

The downlink snr is maximized with drl <cit.>.
The spectrum and energy efficiency are maximized with ris <cit.>.
The interference caused by secondary networks is mitigated with an ris <cit.>.
The bs configuration and the phase shifts of the ris are iteratively optimized with the ofdm transmission scheme <cit.>.
The robust transmission scheme design is addressed as well <cit.>.

Although the above-mentioned pioneering works have achieved reasonably good performances compared to random phase shifts and scenarios without ris,
there is still considerable space for performance improvement
and the computational complexities of these proposed algorithms are still too high for real-time application.
Besides, the high complexity constrains the number of ris antennas.
In most works mentioned above, the number of ris antennas is assumed to be less than 200,
which is far less than the vision of up to thousands of elements introduced in Ref. <cit.>.
In recent years, there have been attempts to apply machine learning to ris optimization for better scalability.
For example, meta-learning is used to optimize the spectrum efficiency but only one user per ris is considered in <cit.>,
a reflection beamforming codebook is learnt with drl for groups of users <cit.> (no spatial multiplexing is considered though),
the achievable rate is predicted with user location and beamforming in <cit.>.

In this work, we propose an approach based on fcn and wmmse precoding
that achieves better performance with less operating expenses.
The fcn is comprised of multiple convolutional layers and was first proposed in <cit.> and
is successfully applied for semantic segmentation and depth estimation <cit.>,
where an image is fed into the neural network as input and the category (for semantic segmentation, categories are, e.g., the car, the road, or the sky) or the depth (for depth estimation, which is the distance between the camera and the object) of each pixel of the image is estimated as output.
The first convolutional layer is able to extract the local features of the image
while the following convolutional layers can process information on higher levels.
Inspired by the rectangular shape of the ris
and the spatial correlation of the channels of adjacent ris antennas
due to the close distance between the antennas <cit.>,
we apply an fcn for the ris configuration.
In the literature, the fcn has been applied to object detection <cit.>,
visual tracking <cit.> and
receiver design <cit.> besides the above-mentioned applications.
To the authors' best knowledge, it has not yet been applied to ris configuration.

Our contributions in this paper are as follows:

* We apply an fcn for the ris configuration.
We show that the fcn is an efficient and scalable architecture for the considered wsr maximization problem.
By introducing dropout layers, we prevent overfitting and improve the generalization of the trained models.

* We design a trick to consider the direct channel without ris in the fcn,
which makes the proposed approach more applicable because a weak direct channel exists in most use cases.

* We design a training process with the wmmse precoder,
which is nondifferentiable and therefore incompatible with the gradient ascent optimization.
By alternatingly training the fcn and updating the precoding vectors, the fcn can work properly with the wmmse precoder.

* We propose a method to discretize the output phase shifts with low resolution (e.g., $\pi$ with binary, small codebook size) such that the solution can be applied to ris with simple and cheap discrete phase shifts. The proposed method uses the Frobenius norm as a penalty and aims to minimize the performance loss compared to the continuous phase shift.

In the following part of this work, Section <ref> formulates the problem,
Section <ref> describes the precoding techniques,
Section <ref> proposes our solution,
Section <ref> explains some issues of implementation and training and extends the method to discrete phase shifts,
Section <ref> presents the training and evaluation results and
Section <ref> concludes the work.

Notations:  Throughout this work, we use the following notations:
$\mb{I}_n$ is an identity matrix of size $n \times n$,
$\mb{0}_n$ is a zero matrix of size $n \times n$,
$\mb{A} \star \mb{B}$ denotes the 2-dimensional cross-correlation between $\mb{A}$ and $\mb{B}$,
$|z|$ and $\arg(z)$ are the amplitude and phase of complex number $z$, respectively,
$\mb{A}^+$ denotes the pseudoinverse of matrix $\mb{A}$,
$\bar{\theta}$ refers to the untrainable neural network parameter set $\theta$ (i.e., it is considered constant during training),
$\exp(\mb{A})$ is the elementwise exponential function.

§ PROBLEM FORMULATION

We consider an ris-aided communication with direct propagation path between bs and users.
Our objective is to serve multiple users with the same resource block and to maximize the wsr of the users.
The bs is assumed to have multiple antennas and can do precoding subject to the transmit power constraint $E_{Tr}$.
The ris is assumed to be rectangular with $H$ rows and $W$ columns of antennas.
Each ris antenna has the ability to receive signal, adjust its complex phase without changing its amplitude and transmit it.

We denote the precoding matrix as $\mb{V}$ of size $M\times U$, where $M$ is the number of bs antennas,
$U$ is the number of single-antenna users,
the channel from bs to ris as $\mb{H}$ of size $N\times M$,
where $N$ is the number of ris antennas ($N=WH$) and element $h_{nm}$ in row $n$ and column $m$ represents the channel gain from the $m$th antenna of the bs to the $n$th antenna of the ris.
The ris signal processing is denoted by the diagonal matrix $\boldsymbol{\Phi}$ of size $N\times N$,
where the element $\phi_{nn}$ in row $n$ and column $n$ is $e^{j\psi_n}$,
with $\psi_n$ being the phase shift of antenna $n$.
We consider both continuous and discrete phase shift in this work.
If the phase shift is continuous, $\psi_n \in [0, 2\pi)$.
If the phase shift is discrete, $\psi_n \in \boldsymbol{\Phi}_d$,
where $\boldsymbol{\Phi}_d$ is the set of predefined discrete phase shifts.
The channel matrix from ris to users is denoted as $\mb{G}$ of size $U\times N$,
where the element $g_{un}$ denotes the channel gain from ris antenna $n$ to user $u$.
The direct channel from bs to users is denoted as $\mb{D}$ of size $U\times M$,
where the element $d_{um}$ in row $u$ and column $m$ is the channel gain from the $m$th bs antenna to the $u$th user, the transmission is described as

\begin{equation}
\mb{y} = \left(\mb{G} \boldsymbol{\Phi} \mb{H} + \mb{D} \right) \mb{V} \mb{x} + \mb{n},
\label{eq:transmission_los}
\end{equation}

where $\mb{x}$ is the vector of transmitted symbols,
$\mb{y}$ is the vector of received symbols and $\mb{n}$ is the vector of thermal noise.
All three vectors have the same size of $U \times 1$.
The whole system model is illustrated in Fig. <ref>.

System model.

\begin{equation}
\mb{C} = (\mb{G} \boldsymbol{\Phi} \mb{H} + \mb{D}) \mb{V}
\label{eq:channel}
\end{equation}

and $c_{uv}$ being the element of $\mb{C}$ in row $u$ and column $v$, the objective is to maximize the wsr.
Therefore, the problem can be formulated as

\begin{equation}
\begin{array}[t]{ll}
\max_{\mb{V}, \boldsymbol{\Phi}} & %\begin{aligned}
f = \sum_{u=1}^U\alpha_u\log_2\left(1+\frac{c_{uu}^2}{\sum_{v\neq u}c_{uv}^2+\frac{1}{\rho}}\right)
\\
\text{s.t.} & \tr{\mb{V}\mb{V}^H} \leq E_{Tr}\\
& |\phi_{nn}|=1\\
& |\phi_{nn'}|=0 \text{ for } n \neq n',
\end{array}
\label{eq:problem}
\end{equation}

where $\alpha_u$ is weight of user $u$, $\alpha_u \in [0, 1]$ and $\sum_{u=1}^U\alpha_u = 1$, and $\rho$ is the tsnr,
which is the ratio between transmit power and noise power.

§ THE MMSE AND WMMSE PRECODER

In problem (<ref>), $\mb{H}$, $\mb{G}$ and $\mb{D}$ are given and we would like to optimize both $\mb{V}$ and $\boldsymbol{\Phi}$.
While the ris optimization problem (i.e., optimizing $\boldsymbol{\bs{\Phi}}$) is new, the precoding problem (i.e., optimizing $\mb{V}$) has been intensively studied in the literature.
Different precoding techniques, such as mrt, zf, mmse and wmmse are proposed.
Among them, the mmse precoder can minimize the mean squared error caused by both interference and noise with a closed form solution and is the optimal precoder for a single user <cit.>.
The iterative wmmse precoder introduces weights of the mse of different users in the multi-user mimo system.
It is proven that minimizing the sum of weighted mse for certain weights is equivalent to maximizing the wsr <cit.>.
However, the wmmse precoder does not have a closed-form solution and is only available as an iterative algorithm, which significantly increases its complexity compared to the mmse precoder.
Accordingly we cannot compute its derivative.
In the following, we briefly introduce the two precoders that will be applied jointly with the ris optimization.

§.§ The MMSE Precoder

In the mmse precoder, the error is defined as the difference between transmitted symbols $\mb{x}$ and received symbols $\mb{y}$ scaled by a proper factor $\beta^{-1}$ (because the received symbols are weaker due to the propagation loss).
Accordingly, the mse is defined as

\begin{equation}
\begin{aligned}
\mb{e}= & \mathbb{E}( || \beta^{-1}\mb{y}-\mb{x} ||^2_2 )\\
= & \mathbb{E}\left(||\beta^{-1} (\mb{C}\mb{V}\mb{x} + \mb{n}) - \mb{x}||_2^2\right).
\end{aligned}
\label{eq:error_mmse}
\end{equation}

The objective is to minimize the mse subject to the transmit power constraint:

\begin{equation}
\begin{array}{rl}
\min_{\beta, \mb{V}} & \mathbb{E}\left(||\beta^{-1} (\mb{C}\mb{V}\mb{x} + \mb{n}) - \mb{x}||_2^2\right) \\
\text{s.t.} & \tr{\mb{V}\mb{V}^H} \leq E_{Tr}. \\
\end{array}
\label{eq:mmse_formulation}
\end{equation}

When we assume $\mathbb{E}(\mb{x}\cdot \mb{x}^T)=\mb{I}_2$, $\mathbb{E}(\mb{n}\cdot \mb{n}^T)=1/\rho \cdot \mb{I}_2$ and $\mathbb{E}(\mb{n}\cdot \mb{x}^T)=\mb{0}_2$, we obtain a closed-form solution <cit.>

\begin{equation}
\mb{V} = \beta\left(\mb{C}^H\mb{C}+\frac{1}{\rho}\mb{I}\right)^{-1}\mb{C}^H
\label{eq:mmse}
\end{equation}

\begin{equation}
\beta=\sqrt{\frac{E_{Tr}}{\tr{\left(\mb{C}^H\mb{C}+\frac{1}{\rho}\mb{I}\right)^{-2}\mb{C}^H\mb{C}}}}.
\end{equation}

The mmse precoder is the optimal precoder to maximize the data rate for a single user.
However, it uses the same scalar $\beta^{-1}$ to scale all received symbols of multiple users.
If the received symbols at different users have very different signal strengths due to different channel gains,
the mmse precoder is forced to equalize the channel gains,
which results in a suboptimal performance.
Besides, it minimizes the sum of mean squared errors of all users.
If we would like to optimize the wsr,
the weights (i.e., $\alpha_u$ in (<ref>)) cannot be considered in the mmse precoder.
Therefore, the mmse precoder is not the optimal precoder to maximize the wsr with multiple users.
Nevertheless, it is still a good choice due to its satisfying performance and simple, closed form (and therefore differentiable) solution.

§.§ The WMMSE Precoder

Similar to the mmse precoder, we define the error of user $u$ as

\begin{equation}
e_u = \xi_u y_u - x_u,
\label{eq:wmmse_error}
\end{equation}

where $\xi_u$ is the scaling factor of user $u$ (corresponding to $\beta^{-1}$ in (<ref>) but can be different for different users),
$y_u$ is the received symbol of user $u$ (i.e., the $u$th element of $\mb{y}$ in (<ref>))
and $x_u$ is the transmitted symbol of user $u$ (i.e., the $u$th element of $\mb{x}$ in (<ref>)).

The wmmse precoder considers the following problem

\begin{equation}
\begin{array}{rl}
\min_{w_u, \mb{v}_u, \xi_u} & \sum_{u=1}^U\ \alpha_u (w_u \mathbb{E}(e_u e_u^H) - \log w_u)\\
\text{s.t.} & \sum_{u=1}^U \tr{\mb{v}_u \mb{v}_u^H} \leq E_{Tr},
\end{array}
\label{eq:wmmse_formulation}
\end{equation}

where $w_u$ is the weight of the mse of user $u$,
$U$ is the number of users,
$\mb{v}_u$ is the precoding vector for user $u$ of size $M \times 1$ (i.e., the $u$th column of $\mb{V}$ in (<ref>)).
It is proven that problem (<ref>) is equivalent to the problem of wsr maximization (<ref>) (when $\bs{\Phi}$ is fixed) in the sense that their optimal solutions of $w_u$ and $\mb{v}_u$ are the same when $\xi_u$ is carefully selected for all users $u$.

Unfortunately, there is no closed-form solution to (<ref>).
Instead, Algorithm <ref> is proposed to iteratively solve problem (<ref>),
where we define $\mb{c}_u$ as the channel to user $u$ of size $1 \times M$ (i.e., the $u$th row of $\mb{C}$),
$\mu^*$ is the Lagrange multiplier of the power constraint and
is chosen such that the transmit power constraint in (<ref>) is satisfied.

Iterative WMMSE Precoder

Initialize $\mb{v}_u$
$w_u' \leftarrow w_u$ for all users $u$
$\xi_u \leftarrow \left( \sum_{\nu=1}^U \mb{c}_u \mb{v}_\nu \mb{v}_\nu^H \mb{c}_u^H + 1/\rho \right)^{-1} \mb{c}_u \mb{v}_u$ for all users $u$
$w_u \leftarrow \left( 1 - \xi_u^H \mb{c}_u \mb{v}_u \right)^{-1}$ for all users $u$
$\mb{v}_u \leftarrow \alpha_u \left( \sum_{\nu=1}^U \alpha_\nu \mb{c}_\nu^H\xi_\nu w_\nu \xi_\nu^H \mb{c}_\nu + \mu^* \mb{I}_M \right)^{-1} \mb{c}_u^H \xi_u w_u$ for all users $u$
$\big| \sum_{u=1}^U |w_u| - \sum_{u=1}^U |w'_u| \big|\leq \epsilon$

The wmmse precoder is the optimal precoder to maximize the wsr (see Theorem 1 in <cit.> for details). However, its complexity is considerably higher than the mmse precoder. Besides, the iterative solution is not differentiable.

§ PROPOSED FCN BASED APPROACH

§.§ The FCN Approach

An fcn consists of several two-dimensional convolutional layers,
where each layer processes the feature maps of its preceding layer and extracts new feature maps from them.
For example, in semantic segmentation and depth estimation, usually an RGB image is used as input for the network.
The filters in the first layer operate on the three color channels of the image and extract feature maps for the following convolutional layer.
The filter kernels learned in the first layer are three-dimensional as they extract features from a three-dimensional input tensor.
Each filter kernel generates one output feature map.
The output of a semantic segmentation fcn in a certain task is a category for each pixel (e.g., road, car, pedestrian)
wheres the output of a depth segmentation fcn is the estimated depth for each pixel.
The channels $\mb{H}$ and $\mb{G}$ reveal spatial correlation due to the short distance between ris antennas.
This property shows similarity to images, where adjacent pixels form an object,
which can be captured by a two-dimensional convolutional filter,
thereby justifying the choice of an fcn architecture.
The optimal output of each ris antenna depends on both its own feature and features of other antennas,
which is the same to semantic segmentation and depth estimation[Note that the neural network is an optimized mapping from channel features to the ris phase shifts.
For different channels as inputs,
the phase shifts are different.
Therefore, the achieved data rate is not the ergodic rate.].
Encouraged by these analogies, we would like to compute $\bs{\Phi}$ with the fcn.
The elementary unit of an ris is an ris antenna.
Its features are the properties characterizing the wireless channels with respect to this ris antenna, which will be defined later in Section <ref>.
Other conventional neural network architectures, such as neural networks with dense layers would not work for such high-dimensional input and output.

As shown in Fig. <ref>(a), the input of the fcn is a three-dimensional array,
where the first dimension is the feature map index $k \in \{1, \dots, K \}$ with $K$ the number of features,
second and third dimensions are height $h \in \{ 1, \dots, H\}$ and width $w \in \{ 1, \dots, W\}$, respectively.
For example, the value of feature map index $k$ of the ris antenna at position $(h, w)$ is stored in coordinate $(k, h, w)$ in the three-dimensional array.
The output of the fcn is a two-dimensional array of the same height and width as the input,
but with only one feature map being the phase shift of the antenna,
as shown in Fig. <ref>(b).

[Input of the FCN]
[Output of the FCN]
Input and output format of the FCN. The input is a three-dimensional array with dimensions of channel features, height and width. Channel feature $k$ of RIS antenna at position $(h, w)$ is stored in coordinate $(k, h, w)$. The output is a two-dimensional array with dimensions of height and width. The phase shift of RIS antenna at position $(h, w)$ is stored in coordinate $(h, w)$.

Define $\mb{Q}_{\tau}$ as the $\tau$th feature map of the output of a convolutional layer.
It is computed as

\begin{equation}
\mb{Q}_{\tau} = b_{\tau} + \sum_{k=1}^K \mb{W}_{k\tau} \star \mb{P}_k,
\label{eq:conv2d}
\end{equation}

where $b_{\tau}$ is the bias of feature map with index $\tau$,
$\mb{W}_{k\tau}$ is the filter kernel which takes feature map with index $k$ as input and produces feature map with index $\tau$ as output,
$\mb{P}_k$ is the array of feature map with index $k$ of the input,
$K$ is the number of feature maps of the input.

When applied in a straightforward manner, the cross-correlation operation makes the size of $\mb{Q}_{\tau}$ smaller than the size of $\mb{P}_k$.
We use zero padding <cit.> such that the size of output equals the input size,
which is illustrated in Fig. <ref>.
Intuitively, the number of zeros added to each side of the input in the horizontal direction is

\begin{equation}
w_z = \frac{w_f - 1}{2}
\end{equation}

where $w_f$ is the filter size in the horizontal direction.
Note that $w_f$ has to be an odd number such that $w_z$ is an integer.
A similar relationship can be obtained for the vertical direction.
In this way, the sizes of input and output are kept identical.

[step=1cm] (0,0) grid (5, 5);
[step=1cm, dotted] (-1,-1) grid (6, 6);
[blue!20!white] (-1,3) rectangle (2,6);
[blue!50!white] (0,4) rectangle (1,5);
[step=1cm] (-1,3) grid (2, 6);

Illustration of filter and zero padding.
A filter of size $3\times 3$ (the light blue area) process information in the units it covers according to (<ref>) and stores the result in the dark blue unit in the output of the layer.
Since the filter size is $3\times 3$, the output size is reduced by 2 in both width and height.
In order to keep the output size identical to the input size,
we add zeros around the actual input (dotted units around the solid units),
which is referred to zero padding.
The zeros are also in accordance with reality: since the channel features are amplitudes and phases.
For units without an antenna (the dotted units),
the amplitude is 0 and the phase does not matter.

The size of the filters and the number of convolutional layers should be designed in such a way,
that each output elementary unit has its receptive field[The region of the input that produces the output] of the whole ris.
This idea is illustrated in Fig. <ref>.

Since every filter aggregates information from both sides in the horizontal direction by $(w_f - 1) / 2$,
the required number of layers is

\begin{equation}
N_{\text{layers}} = \frac{(w_f - 1) / 2}{W},
\label{eq:number_layers}
\end{equation}

such that the receptive field of a unit on one edge covers a unit on the opposite edge.
Equation (<ref>) can be extended to the vertical layer easily.

[blue!10!white] (1, 1) rectangle (8, 8);
[blue!30!white] (2, 2) rectangle (7, 7);
[blue!60!white] (3, 3) rectangle (6, 6);
[blue] (4,4) rectangle (5,5);
[step=1cm] (0,0) grid (9, 9);

Illustration of receptive field in different FCN layers. We consider the darkest unit in the middle and a filter size of $3\times 3$. After the first convolutional layer, the unit under consideration is in the receptive field of the units in the second darkest areas. After the second convolutional layer, the unit under consideration is in the receptive field of the units in the lighted areas. In this way, the information in the considered unit propagates to all units through 4 layers.

In order to prevent overfitting and to improve robustness, we put a dropout layer <cit.> between each convolutional layer,
which randomly turns off a certain portion of neurons.

After all layers, the output of the fcn $\bs{\Psi}$ is the array of phase shifts of all ris antennas of size $H\times W$.
The structure of the proposed fcn is shown in Fig. <ref>.
The computation of the fcn is denoted as $\bs{\Psi} = \p(\bs{\Gamma})$,
where $\theta$ is the parameters of the fcn (i.e., all biases and filters of all layers),
$\bs{\Gamma}$ is the input feature map,
$\bs{\Psi}=(\psi_n)$ is the two-dimensional array with linearly addressed elements $\psi_{n}, n=1, \dots, N$, $n$ is the index of the ris antenna.
For antenna at position $(h, w)$, we have $n = H \cdot (h - 1) + w$.

(gamma) [rectangle] $\bs{\Gamma}$;
[box, below of=gamma] (conv1) Convolutional 1;
[box, below of=conv1, yshift=-.2cm] (dropout1) Dropout 1;
[->] (gamma) – (conv1);
[->] (conv1) – (dropout1);
[rectangle, below of=dropout1, yshift=-.0cm] (dots) …;
[box, below of=dots, yshift=-.2cm] (dropout7) Dropout 7;
[box, below of=dropout7, yshift=-.2cm] (conv8) Convolutional 8;
[->] (dots) – (dropout7);
[->] (dropout7) – (conv8);
[rectangle, below of=conv8, yshift=-.2cm] (psi) $\mb{\Psi}$;
[->] (conv8) – (psi);

Structure of the proposed FCN for phase shift optimization.

§.§ Definition of Features

In problem (<ref>),
we assume that the csi at the ris controller is known, which is a common assumption in relevant works <cit.>.
Although the channel estimation involving an ris is challenging due to the large number of ris antennas,
there are multiple promising attempts to estimate the channel with sufficient accuracy in real-time, such as <cit.>.
Besides, the csi can also be used as intermediate information to derive an ris optimization method with an input easier to obtain than the csi itself,
such as user locations <cit.>.
Therefore, we believe that csi based ris optimization has practical importance in the near future.

Among the channel matrices,
$\mb{H}$ is assumed to be constant because both bs and ris are stationary,
$\mb{G}$ and $\mb{D}$ vary with the user positions and are required to compute $\bs{\Phi}$ and $\mb{V}$.
We would like to interpret $\mb{G}$ and $\mb{D}$ as two feature maps with element index $n$ referring to a certain ris antenna. Accordingly, we can use our fcn.
While each column of $\mb{G}$ can be injectively mapped to each ris antenna,
$\mb{D}$ is irrelevant to the ris. Therefore, we define $\mb{J}=\mb{D}\mb{H}^+$ and (<ref>) becomes

\begin{equation}
\mb{y} = \left(\mb{G} \boldsymbol{\Phi} + \mb{J}\right) \mb{H} \mb{V} \mb{x} + \mb{n}.
\label{eq:transmission_los_equivalent}
\end{equation}

Equation (<ref>) can be interpreted as follows:
the precoded signal $\mb{Vx}$ is transmitted through channel $\mb{H}$ to every ris antenna and then through channel $\mb{G}\bs{\Phi} + \mb{J}$ to the users.
Both $\mb{G}$ and $\mb{J}$ have $N$ columns and their columns can be injectively mapped to the ris antennas.
Therefore, the features of an ris element at position $(w, h)$ can be defined as

\begin{equation}
\begin{aligned}
\bs{\gamma}_{wh} = ( & |g_{1n}|, \arg(g_{1n}), \dots, |g_{Un}|, \arg(g_{Un}),\\
& |j_{1n}|, \arg(j_{1n}),\dots, |j_{Un}|, \arg(j_{Un})),
\end{aligned}
\label{eq:feature_definition}
\end{equation}

where $g_{k\tau}$ is the element in row $k$ and column $\tau$ of $\mb{G}$,
$j_{k\tau}$ is defined similarly to matrix $\mb{J}$.
In this way, we use 4$U$ features to characterize the wireless channels with $U$ users with respect to each ris antenna $n$, yielding the fcn input $\bs{\Gamma}=(\mb{G}, \mb{J})$.

§.§ Objective Functions

Our objective function $f$ in (<ref>) is a function of channels $\mb{G}$ and $\mb{D}$,
ris phase shifts $\bs{\Phi}$ and bs precoding $\mb{V}$.
When $\bs{\Phi}$ is determined,
$\mb{V}$ can be computed with either mmse precoder or wmmse precoder as described in Section <ref>,
which are denoted as $\mb{V}_{\text{MMSE}}(\mb{G}, \mb{D}, \bs{\Phi})$ and $\mb{V}_{\text{WMMSE}}(\mb{G}, \mb{D}, \bs{\Phi})$, respectively.
It is to note that $\mb{V}_{\text{MMSE}}$ is a differentiable function but $\mb{V}_{\text{WMMSE}}$ is not differentiable.

Since $f$ is a function of $\mb{G}$, $\mb{D}$, $\bs{\Phi}$ and $\mb{V}$,
we can express $f$ as $f(\mb{G}, \mb{D}, \bs{\Phi}, \mb{V})$. Therefore, for mmse precoding, the objective function is

\begin{equation}
\begin{aligned}
f\bigg(& \mb{G}, \mb{D}, \exp\big(j\p(\mb{G}, \mb{D})\big), \\
& \mb{V}_{\text{MMSE}}\Big(\mb{G}, \mb{D}, \exp\big(j\p(\mb{G}, \mb{D})\big)\Big)\bigg).
\end{aligned}
\label{eq:objective_mmse}
\end{equation}

For wmmse precoding, the objective function is

\begin{equation}
\begin{aligned}
f\bigg(& \mb{G}, \mb{D}, \exp\big(j\p(\mb{G}, \mb{D})\big), \\
& \mb{V}_{\text{WMMSE}}\Big(\mb{G}, \mb{D}, \exp\big(j\bs{\Psi}_{\bar{\theta}}(\mb{G}, \mb{D})\big)\Big)\bigg).
\end{aligned}
\label{eq:objective_wmmse}
\end{equation}

Please note that $\bs{\Psi}_{\bar{\theta}}$ in (<ref>) is untrainable, i.e., $\bar{\theta}$ is given and $\bs{\Psi}_{\bar{\theta}}(\mb{G}, \mb{D})$ is a constant.
Therefore, we do not need to compute the derivative of $\mb{V}_{\text{WMMSE}}$ since it does not contain trainable variable $\theta$.

We can then compute $\frac{\partial f}{\partial \theta}$ to do gradient ascent to optimize $\p$ such that $f$ is maximized for both (<ref>) and (<ref>).
The compositions of the objective functions is illustrated in Fig. <ref>.
The Adam optimizer <cit.> is chosen to perform the optimization.

[With MMSE precoder].45!

[With WMMSE precoder].45!

Composition of the objective function.
With the MMSE precoder, the input of the FCN $\bs{\Gamma}$ comprises $\mb{G}$ and $\mb{D}\mb{H}^+$.
The channel $\mb{C}$ is determined by $\mb{G}$, $\mb{D}$ and $\bs{\Phi}$, see (<ref>).
The MMSE precoding matrix $\mb{V}$ is computed with (<ref>).
The WSR is determined by $\mb{C}$ and $\mb{V}$ together.
With the WMMSE precoder, the structure of the objective function is similar,
but $\mb{V}$ is no more a function of $\mb{G}$, $\mb{D}$ and $\bs{\Phi}$,
but a constant which is updated periodically.
Therefore, we do not compute the derivative of $\mb{V}$.

Since the gradient ascent improves the objective monotonically and the objective wsr is bounded from above with a finite transmit power,
the iterative optimization must converge in the end of training.
In practice, we use stochastic gradient ascent instead of gradient ascent.
The gradient is computed for a randomly sampled batch of data rather than the whole data set.
After that, another batch with newly sampled data is used for the next gradient ascent step.
As a result,
we cannot guarantee monotonicity in the optimization.
However, the stochastic gradient ascent is a standard technique in machine learning and machine learning practice shows that if the data samples are identically distributed and the batch size is sufficiently large,
the stochastic gradient ascent shows solid convergence behavior,
which is also true in our case,
as will be shown in Section <ref>.

Note that the proposed approach is an unsupervised learning method.
In supervised learning, a target is provided for each input.
The optimizer tries to minimize the difference between the output of the neural network and the target (or label).
The loss function is usually defined as mse (for continuous output) or cross entropy (for categorical output).
By minimizing the loss function, the neural network learns to reproduce the given target.
In our problem, the target is unknown and therefore not given.
By maximizing the objective function (<ref>) or (<ref>),
the neural network learns an optimized mapping from channel features $\boldsymbol{\Gamma}$ to the ris phase shift $\boldsymbol{\Psi}$ that maximizes the wsr.
We believe that the unsupervised approach has a wider range of application since the target is unavailable for many difficult problems.
Besides, if the target is already available,
it is difficult to justify why we use machine learning to reproduce the target that is already known
(common reasons are, e.g., it is very time-consuming to obtain the target and the evaluation of the neural network is faster than the original algorithm <cit.>, or the purpose of machine learning is the predistortion of a distorted signal and the target (the original signal) is therefore known <cit.>).
Nevertheless, it is of great practical importance and also very exciting if the neural network finds a good solution by itself without supervision.

§ TRAINING PROCESS FORMULATION AND IMPLEMENTATION CONSIDERATIONS

§.§ The Two-Phase Training

As described in Section <ref>, the mmse precoder is simple and differentiable and the wmmse precoder has a higher performance at the cost of complexity and nondifferentiability.
During the training, we apply the mmse precoder and use (<ref>) as the objective function in the first phase of training for a fast training
and then switch to the wmmse precoder and use (<ref>) as the objective function to further improve the wsr in the second phase of training.

In the first phase, objective (<ref>) is a function composition and its derivative can be computed with the chain rule.
The mmse precoding matrix is computed with the current phase shifts $\bs{\Psi}$.
On the contrary, in the second phase, the precoding matrix $\mb{V}_{\text{WMMSE}}$ is nondifferentiable and is considered constant.
Since $\mb{V}_{\text{WMMSE}}$ changes with the channel, when $\bs{\Psi}$ changes because of the training, $\mb{V}_{\text{WMMSE}}$ is no more valid.
Therefore, we update $\mb{V}_{\text{WMMSE}}$ every 10 epochs during training.

In the iterative wmmse precoding algorithm (Algorithm <ref>), the initial precoding matrix is arbitrary subject to the transmit power constraint.
The iteration is supposed to continue until convergence.
In our training, however, it is desirable that
(a) we only run a few iterations such that the training time is reasonable,
(b) the updated precoding matrix is close to the precoding matrix from the previous update,
such that a drastic change is avoided and the training process is more fluent.
With these considerations, we use the previous precoding matrix as the initial precoding matrix in the next update and run 5 iterations per update.
In the first wmmse precoding matrix update, the mmse precoding matrix is applied.
Although the wmmse precoding matrix is not updated in real time as the mmse precoding matrix,
the change of the wmmse is smooth and the training process is stable,
as will be shown in Section <ref>.

Besides the above-mentioned complexity consideration,
another important reason for applying mmse precoder in the first phase of the training is that the wmmse precoder selects users by weighting them with different $w_u$ in Algorithm <ref>.
If the channel condition is poor and transmit power is insufficient following the water-filling principle <cit.>,
the wmmse precoder selects good users and allocate no power to poor users,
which is the situation for many data samples in the beginning of training since the ris is randomly initialized and cannot realize a good channel.
During the ris optimization, since the poor user does not have any transmit power,
the ris is optimized only for the good user and the poor user is ignored.
Therefore, the wmmse precoder does not allocate any transmit power to the poor user in the next iteration because the channel to the good user is improved by the ris but the channel to the poor user is not improved.
As a result, user selection under the initial poor channel quality may last until the end of training,
which results in a suboptimal performance.
On the contrary, the mmse precoder treats the two users equally
and the ris tries to improve the channel quality for both users.
In the end of the first phase,
channel quality is significantly improved compared to the initialization.
Therefore, user selection of the wmmse precoder abandons less users,
which results in a high wsr and better fairness.
The two-phase training is a reasonable choice especially when the initial channel quality is poor.

The proposed training is formulated in Algorithm <ref>.

FCN Training in two phases

Compute the gradient of the objective (<ref>) w.r.t. the neural network parameters
Perform a gradient ascent step with the Adam optimizer
WSR stops increasing
Compute WMMSE precoding vectors $\mathbf{v}_u$ for all $u$ according to Algorithm <ref>
Compute the gradient of the objective (<ref>) with current $\mathbf{v}_u$ w.r.t. the neural network parameters
Perform a gradient ascent step with the Adam optimizer
Preset number of iterations reached
WSR stops increasing

§.§ Solution for Discrete Phase Shifts

Although the above-described method optimizes the phase shifts of ris antennas such that the wsr is improved significantly,
a practical ris is expected to have discrete phase shifts with a high granularity (e.g., $\pi$) rather than continuous phase shifts in order to reduce the hardware complexity and cost.
The phase shift of each antenna is chosen from a finite set of discrete phases $\boldsymbol{\Phi}_d$,
rather than a continuous value.
In this section, we propose a penalty term to convert continuous phase shifts to discrete phase shifts while keeping the performance loss as low as possible.
The penalty term is defined as the smallest Frobenius norm between continuous phase shift and the closest feasible discrete phase defined by $\boldsymbol{\Phi}_d$:

\begin{equation}
p = \sqrt{\sum_n \left(\min_{\varphi_d \in \boldsymbol{\Phi}_d} \varphi_n - \varphi_d \right)}.
\label{eq:penalty}
\end{equation}

The objective is defined as

\begin{equation}
\sum_{u=1}^U \alpha_u \log_2\left( 1 + \frac{c_{uu}}{\sum_{v \neq u} c_{uv} + \frac{1}{\rho}} \right) - \pc p
\label{eq:discrete_objective}
\end{equation}

where $\pc$ is the factor for the penalty term. The training process with discrete phase shifts is formulated as Algorithm <ref>.

FCN Training with Discrete Phase Shifts

Compute the gradient of the objective (<ref>) w.r.t. the neural network parameters
Perform a gradient ascent step with the Adam optimizer
WSR stops increasing
Set $\pc = 0$
Compute WMMSE precoding vectors $\mathbf{v}_u$ for all $u$ according to Algorithm <ref>
Compute the gradient of the objective (<ref>) with current $\mathbf{v}_u$ w.r.t. the neural network parameters
Perform a gradient ascent step with the Adam optimizer
Preset number of iterations reached
WSR stops increasing
$\pc \leftarrow \pc + 0.05$
$p$ is smaller than a given threshold

§.§ Implementation

We employ PyTorch <cit.> as the machine learning framework,
which only supports real-valued inputs and parameters.
Therefore, the channel features (<ref>) use amplitude and phase to describe a complex channel gain.
On the contrary, since PyTorch 1.9, complex numbers are supported in the objective function,
which significantly reduces the complexity in implementation of (<ref>) and (<ref>).
After implementation of the objective function, the gradient is computed with PyTorch's differentiation engine Autograd <cit.> without human intervention.

§ TRAINING AND EVALUATION RESULTS

In this section, we present the training and evaluation results.
We apply a ray tracing channel simulator in an urban environment <cit.>,
where mobile users are supposed to have high demand for data rates.
We generate channel from the bs to the ris $\mb{H}$,
channels from the ris to the users $\mb{G}$ and the direct channels $\mb{D}$ of 5000 pairs of users,
whose positions are randomly generated with a minimum distance of 2 meters between them.
The scenario is shown in Fig. <ref>,
where a los propagation path between bs and users is assumed to be unavailable due to the blockage by a building in the middle.
A weak direct path is available through reflection on the building in the lower right corner.
An ris is equipped such that there are los propagation paths from bs to ris and from ris to users.
A building behind the bs enables a strong reflection path from bs to ris such that channel matrix $\mb{H}$ has a rank higher[The rank of a matrix product is less than or equal to the rank of every matrix in the product.
Since $\mb{G}$ is the channel matrix to two users at different positions,
it has a rank of 2.
$\bs{\Phi}$ is a diagonal matrix and therefore has full rank.
Therefore, $\mb{H}$ has to have a rank higher than 1,
otherwise we have to use the weak direct channel $\mb{D}$,
which leads to a poor performance.] than 1.

Note that we assume the channel from bs to ris $\mb{H}$ is a constant,
which holds because both positions of bs and ris are fixed.
However, it requires per-bs-ris-pair channel measurement and training before deployment and when surrounding environment changes (e.g., a new building is built nearby),
which should be plausible since per-cell optimization is a common practice <cit.>.

The considered scenario.
Positions of BS and RIS are given and positions of user pairs are generated randomly inside the yellow area.
The LoS propagation path between BS and users is blocked by the high building in the middle.
However, LoS propagation paths from BS to RIS and from RIS to users exist.
A building behind the BS enables a strong reflection path from BS to RIS.
A weak direct path is available through reflection on the building in the lower right corner.

The important parameters of scenario and model are presented in Table <ref>.
If multiple parameter values are considered, the default value is marked with a star (*),
which is always used without explicit specification.

Scenario and model parameters

Parameter     Value

Number of bs antennas     9

ris size     (16, 16)*, (32, 16), (32, 32)

Carrier frequency     5.8 GHz

Distance between adjacent antennas at bs     0.5 wave length

Distance between adjacent antennas at ris     0.25 wave length

Transmit snr     $10^{11}$*, $1.1\times 10^{11}$, …$10^{12}$

Weights of users     (0, 1), (0.25, 0.75), (0.5, 0.5)*, (0.75, 0.25), (1, 0)

Filter size     (5, 5) for ris size of (16, 16), (13, 13) otherwise

Zero padding size     (2, 2) for ris size of (16, 16), (6, 6) otherwise

Number of convolutional layers     8

Dropout rate     0.1 for ris size of (16, 16), 0.35 otherwise

Training rate with mmse precoder     $10^{-4}$

Epoches with mmse precoder     4000

Training rate with wmmse precoder     $2\times 10^{-6}$

Updating interval with wmmse precoder     10 Epoches

Epoches with wmmse precoder     4000

Batch size     256

Optimizer     ADAM

Number of data samples in training set     5000

Number of data samples in testing set     1024

Fig. <ref> shows the improvement of the wsr with two users, a tsnr of $10^{11}$ and user weights of (0.5, 0.5) in training (therefore, the sum rate is the wsr in the figure times 2).
As explained in Section <ref>,
we first use the mmse precoder to train the model for 4000 epochs,
then switch to the wmmse precoder for further training.
It can be observed that the training has improved the wsr significantly with the mmse precoder and the wmmse precoder outperforms the mmse precoder.

Improvement of WSR in our proposed two-phase training protocol with pretraining with the MMSE precoder and “booster” training with the WMMSE precoder.

Fig. <ref> shows the first 100 epochs of the training with the wmmse precoder (i.e., epochs 4000 - 4100 in Fig. <ref>).
Time points, where the precoding vectors are updated with the wmmse precoder are marked with the red dots.
At the beginning of the training with the wmmse precoder,
there is always a considerable performance improvement after the precoding vectors are updated.
The performance is then improved slightly due to the fcn improvement.
The improvement slows down as the training carries on until convergence.

Details to Fig. 6: Training with WMMSE precoder with update of precoding matrix every 10 epochs.

Fig. <ref> shows the wsr improvement with four users and two-phase training.
Since the weights of all four users are 0.25,
the sum rate is the wsr in the figure times 4.
Compare Fig. <ref> and Fig. <ref>,
we can observe that the scenario with four users achieves slightly better performance with longer training time.
The reason for the better performance is the higher potential to reuse the spectrum resource due to the larger number of users
whereas the reason for the longer training time is the more complicated objective function and neural architecture.
However, the advantage of four users compared to two users is only marginal.
As described in Section <ref>,
the wmmse precoder allocates no transmit power to poor users when the available transmit power is insufficient following the water-filling principle.
In order to verify whether the marginal advantage of four users is due to the transmit power constraint,
we train a model with a tsnr of $10^{13}$ while keeping other parameters unchanged[Note that this is not realistic assumption given a typical noise power.
The purpose of this assumption is only to validate whether the transmit power constraint prevents the four-user scenario from outperforming the two-user scenario significantly.].

Improvement of WSR in the two-phase training protocol for 4 users.

The comparison of achieved sum rate in the testing set (i.e., data that are not seen by the optimizer such that we can test the generality of the trained model) with equal user weights and different algorithms is shown in Fig. <ref>,
where we choose the bcd algorithm in <cit.> as a baseline.
It can be observed that all three setups achieve similar sum rates with the the tsnr of $10^{11}$
but the two setups with four users significantly outperform the setup with two users with the tsnr of $10^{13}$
and the proposed algorithm is better than the bcd algorithm.
This result validates that the transmit power is the constraint that prevent the four-user scenario from outperforming the two-user scenario significantly.

Sum rate with 2 and 4 users and different algorithms.

Fig. <ref> provides further insights into the transmit power constraint,
where the numbers of data samples with different numbers of users with a data rate of 0 are shown for both transmit powers.
As shown in Table <ref>, 1024 data samples are tested in total.
We can observe that two users have a data rate of 0 in more than half (577) of all data samples with the tsnr of $10^{11}$,
which makes the four-user scenario equivalent to the two-user scenario.
Therefore, the average performance difference between the two scenarios is very small.
On the contrary, no user has 0 data rate with the tsnr of $10^{13}$.
This is because the resource is fully reused and the sum rate outperforms the two-user scenario as a result.

Number of users with a data rate of 0 in testing of the four-user scenario with the proposed algorithm.

To compare the proposed solution with other algorithms,
we choose random ris phase shifts and phase shifts optimized with the bcd algorithm proposed as Algorithm 2 in <cit.> as two baselines.
The number of iterations is increased from 100 in the original paper to 5000 due to the larger number of ris antennas.
All three algorithms use the wmmse precoder.
As before, the fcn is evaluated with the testing data set.

Fig. <ref> shows the rate regions of the three algorithms.
Each curve is obtained with the user weights of (0, 1), (0.25, 0.75), (0.5, 0.5), (0.75, 0.25) and (1, 0) and
every point is the average of performances of 1024 channel samples.
We can see that both the bcd algorithm and the fcn improve the data rates of two users significantly from the random initialization.
The fcn outperforms the bcd algorithm clearly.
Besides, although training fcn and running bcd algorithm have time consumption in the same order of magnitude,
training the fcn is only done once before the application.
Testing fcn with 1024 channel samples takes merely a few seconds
while running the bcd algorithm takes more than 4 hours with a 12-core CPU
(channel samples are processed in parallel with the parfor loop in MATLAB).
This makes the proposed approach much more promising to be deployed in the near future
for real-time applications given hardware with comparable performance as of today.
Besides, Fig. <ref> also shows the achieved sum rate with random channel estimate error,
which is assumed to be an iid Gaussian random error.
The ratio between standard deviation of the Gaussian random error and the mean channel gain is denoted as $\gamma$.
We can observe that the fcn is robust against a channel estimate error
partially because of the dropout layers,
which add noise to inputs of each convolutional layer in the training and therefore enhance robustness of the neural work.
This result shows that the proposed solution does not require very precise channel estimation in order to work,
which would be very challenging in ris-aided systems.

Average rate regions of baselines and proposed algorithm.

Fig. <ref> shows the sum rate with the same user weights (0.5, 0.5) and different tsnr.
We can see that higher tsnr (i.e., transmit power) increases the data rate considerably
and the ris configured by the fcn outperforms random phase shifts and phase shifts optimized with the bcd algorithm for all tested tsnr.
It is also to note that the performance of the fcn is robust against random channel estimation errors.

In the evaluation, we have obtained a wsr of 6.87 bit/s/Hz for tsnr=$10^{11}$.
Compared to the wsr at the end of training (7.26 bit/s/Hz),
the wsr in the evaluation is roughly 5% less than the wsr in the training,
which suggests a low overfitting level and a good generalization of the trained model.

Sum rates of baselines and proposed algorithm with different TSNRs.

Fig. <ref> shows the sum rates of different approaches with different number of ris antennas.
It is to see that our proposed solution outperforms the two baselines with all tested numbers of antennas
and is robust against random channel estimation errors.
Besides, although the fcn requires more training time due to bigger filter sizes (see TABLE <ref>) and more ris antennas,
the evaluation of the fcn is within two minutes even for 1024 antennas on the author's two-year-old laptop with Intel's 8th generation i7 CPU (evaluation is run on the CPU) and 8 GB RAM,
as shown in TABLE <ref>.
In contrast, the bcd algorithm takes more than 3 days to optimize the ris configuration with 1024 antennas on the server with a 12-core CPU,
which makes it almost impractical for real-time deployment.

Sum rates of baselines and proposed algorithm with different numbers of RIS antennas.

Time consumption of different approaches

# of antennas     256     512     1024

Random     Seconds     Seconds     Seconds

BCD     5 hours     1 day     3 days

FCN training     5 hours     2 days     3 days

FCN evaluation     Seconds     Seconds     Minutes

Fig. <ref> shows the ecdf of the sum-rate with 256 ris antennas, a tsnr of $10^{11}$ and the user weights of (0.5, 0.5),
where we can confirm the advantage of the proposed algorithm again.
It is also to note that the ecdf of the proposed algorithm is flatter than the other two algorithms
and the outage probability that the wsr is less than a given threshold is higher with the proposed algorithm than with the bcd algorithm proposed in <cit.>,
suggesting that the performance depends more heavily on the channels.
The proposed algorithm only brings significantly advantage with favorable channels.
Therefore, user selection is crucially important for the best performance,
which remains an open problem.

ECDF of baselines and proposed algorithm.

Finally, we evaluate the performance with discrete phase shifts with a granularity of $\pi$.
Fig. <ref> shows the penalty (defined in (<ref>)) development during training with the wmmse precoder.
Algorithm <ref> effectively reduces the difference between the output phase shift and the predefined discrete phase shifts.
The wsr reduces as the fcn output is being discretized.
After training, we force the output phase to be one of the predefined phases and compare the result between the model without discretization (i.e., with Algorithm <ref>) and the model with discretization (i.e., with Algorithm <ref>),
which realize a sum rate of 5.72 bit/s/Hz and 6.05 bit/s/Hz, respectively.
This result suggests that the proposed discretization method in Section <ref> realizes only a marginal improvement compared to brutal discretization.
To the authors' best knowledge,
it remains an open problem to obtain the global optimum of a high dimensional combinatorial optimization problem.

Besides the simpler ris hardware,
another benefit of the discrete phase shift is the small amount of controlling signal.
With the default assumption of number of ris antennas
and two discrete phase shifts,
which can be represented by one bit,
the ris controlling signal has a data amount of 256 bit,
which can be easily transmitted within the next-generation wireless communication systems.

Development of penalty (blue) and WSR (orange) during training during the training phase with the WMMSE precoder. A lower penalty suggests that the output phase shifts are closer to the predefined discrete phase shifts. The WSR reduces as the FCN output is being discretized. The purpose of the training is to minimize the WSR loss during the discretization.

§ CONCLUSIONS

We proposed an fcn based solution for spatial multiplexing enabled by a ris in this paper.
The fcn is a widely applied neural network architecture for computer vision but has not yet been used for ris configuration.
The rectangular shape of the ris and the spatial correlation of channels on adjacent antennas of the ris encourage us to apply the fcn to the ris configuration.
We use the simple and differentiable mmse precoder for pretraining
and the iterative and complex wmmse precoder with periodically updated precoding vectors to fine-tune the model.
A set of channel features is designed to include both cascaded channels via the ris and the direct channel.
We use a ray-tracing simulator to generate channel models with spatial correlation.
Evaluation results show that the proposed solution achieves higher performance and also shows faster evaluation speed than the baselines.
Therefore, it can be better scaled to a large number of antennas.
The dropout layers between the convolutional layers of the fcn reduce overfitting.
A penalty term is introduced to discretize fcn output such that the proposed method can be applied to ris with discrete phase shifts.

It is shown that the wsr depends strongly on the wireless channels.
As future works, it is desirable to find a criterion to optimally decide which users are selected to share the same resource block,
which can significantly improve the average performance.
Furthermore, more advanced techniques of interference management,
such as nonorthogonal multiple access (NOMA) and rate splitting,
can also be combined with ris and may achieve a better performance than treating interference as noise.

The source code and channel data of this paper will be provided for public access if the article is accepted.

§ ACKNOWLEDGEMENT

The results with random ris phase shifts and with bcd algorithm are obtained with the open-source code generously shared by authors of <cit.> under <https://github.com/guohuayan/WSR_maximization_for_RIS_system>.

The authors would like to thank Mrs. X. Shan for providing the ray-tracing channel models,
Dr. K. L. Besser for the fruitful discussion on the combinarotial optimization
and Mr. R. Wang for assisting training and simulation.

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

journaltitleDigital Signal Processing
titleMassive MIMO is a reality—What is next?: Five promising research directions for antenna arrays

journaltitleIEEE Journal on Selected Areas in Communications
titleSmart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead

journaltitleIEEE Wireless Communications
titleHolographic MIMO surfaces for 6G wireless networks: Opportunities, challenges, and trends

journaltitleIEEE Transactions on Wireless Communications
titleMulticell MIMO communications relying on intelligent reflecting surfaces

journaltitleIEEE Transactions on Wireless Communications
titleWeighted sum-rate maximization for reconfigurable intelligent surface aided wireless networks

journaltitleIEEE Transactions on Signal Processing
titleA joint precoding framework for wideband reconfigurable intelligent surface-aided cell-free network

journaltitleIEEE Transactions on Wireless Communications
titleReconfigurable intelligent surfaces for energy efficiency in wireless communication

journaltitlearXiv preprint arXiv:2107.11852
titleReconfigurable Intelligent Surface Phase Hopping for Ultra-Reliable Communications

journaltitleIEEE Wireless Communications Letters
titleIntelligent reflecting surface aided multi-antenna secure transmission

booktitleICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)
titleChannel estimation and low-complexity beamforming design for passive intelligent surface assisted MISO wireless energy transfer

journaltitleIEEE Wireless Communications Letters
titleWeighted sum power maximization for intelligent reflecting surface aided SWIPT

abstractWe examine and compare the different types of linear transmit processing for multiple input, multiple output systems, where we assume that the receive filter is independent of the transmit filter contrary to the joint optimization of transmit and receive filters. We can identify three filter types similar to receive processing: the transmit matched filter, the transmit zero-forcing filter, and the transmit Wiener filter. We show that the transmit filters are based on similar optimizations as the respective receive filters with an additional constraint for the transmit power. Moreover, the transmit Wiener filter has similar convergence properties as the receive Wiener filter, i.e., it converges to the matched filter and the zero-forcing filter for low and high signal-to-noise ratio, respectively. We give closed-form solutions for all transmit filters and present the fundamental result that their mean-square errors are equal to the errors of the respective receive filters, if the information symbols and the additive noise are uncorrelated. However, our simulations reveal that the bit-error ratio results of the transmit filters differ from the results for the respective receive filters. © 2005, IEEE. All rights reserved.
journaltitleIEEE Transactions on Signal Processing
titleLinear Transmit Processing in MIMO Communications Systems
/Users/peng/AppData/Local/Mendeley Ltd./Mendeley Desktop/Downloaded/Joham, Utschick, Nossek - 2005 - Linear Transmit Processing in MIMO Communications Systems.pdf:pdf
Linear transmit processing,Wiener filtering,linear precoding,multiple input,multiple output systems (MIMO),pre-equalization,prerake

journaltitleIEEE Transactions on Communications
titleGeneralized linear precoder and decoder design for MIMO channels using the weighted MMSE criterion

journaltitleIEEE Transactions on Wireless Communications
titleWeighted sum-rate maximization using weighted MMSE for MIMO-BC beamforming design

journaltitleIEEE Transactions on Signal Processing
titleAn iteratively weighted MMSE approach to distributed sum-utility maximization for a MIMO interfering broadcast channel

booktitle2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)
titleAchievable rate maximization by passive intelligent mirrors

journaltitlearXiv preprint arXiv:2110.01700
titleOn the Maximum Achievable Sum-rate of the RIS-aided MIMO Broadcast Channel

booktitle2021 IEEE International Conference on Communications Workshops (ICC Workshops)
titleTwo-user SINR Region for Reconfigurable Intelligent Surface Aided Downlink Channel

journaltitleIEEE Wireless Communications Letters
titleOn Spatial Multiplexing Using Reconfigurable Intelligent Surfaces

journaltitleIEEE Wireless Communications Letters
titleWeighted sum-rate maximization for multi-IRS aided cooperative transmission

journaltitleIEEE Transactions on Signal Processing
titleIntelligent reflecting surface aided multigroup multicast MISO communication systems

journaltitleIEEE Wireless Communications Letters
titleDeep reinforcement learning based intelligent reflecting surface optimization for MISO communication systems

booktitleICC 2020-2020 IEEE International Conference on Communications (ICC)
titleReconfigurable intelligent surface for MISO systems with proportional rate constraints

booktitle2020 IEEE 21st International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)
titleResource allocation for intelligent reflecting surface-assisted cognitive radio networks

journaltitleIEEE Transactions on Communications
titleIntelligent reflecting surface meets OFDM: Protocol design and rate maximization

booktitle2020 IEEE Wireless Communications and Networking Conference (WCNC)
titleIRS-enhanced wideband MU-MISO-OFDM communication systems

journaltitleIEEE Transactions on Signal Processing
titleA framework of robust transmission design for IRS-aided MISO communications with imperfect cascaded channels

booktitleICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)
titleMeta-Learning for 6G Communication Networks with Reconfigurable Intelligent Surfaces

journaltitlearXiv preprint arXiv:2109.14909
titleLearning reflection beamforming codebooks for arbitrary RIS and non-stationary channels

journaltitleIEEE Open Journal of the Communications Society
titleA deep learning based modeling of reconfigurable intelligent surface assisted wireless communications for phase shift configuration

booktitleProceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)
titleFully convolutional networks for semantic segmentation

journaltitlestanford. edu
titleFully convolutional network for depth estimation and semantic segmentation

journaltitleIEEE Journal on Selected Areas in Communications
titleSpatially-stationary model for holographic MIMO small-scale fading

booktitleAdvances in neural information processing systems
titleR-FCN: Object detection via region-based fully convolutional networks

booktitleProceedings of the IEEE International Conference on Computer Vision (ICCV)
titleVisual tracking with fully convolutional networks

journaltitleIEEE Transactions on Wireless Communications
titleDeepRx: Fully convolutional deep learning receiver

booktitle2017 International Conference on Engineering and Technology (ICET)
titleUnderstanding of a convolutional neural network

journaltitlearXiv preprint arXiv:1207.0580
titleImproving neural networks by preventing co-adaptation of feature detectors

journaltitleIEEE Journal on Selected Areas in Communications
titleRobust and secure wireless communications via intelligent reflecting surfaces

booktitle2019 IEEE Globecom Workshops (GC Wkshps)
titleIntelligent reflecting surface for downlink non-orthogonal multiple access networks

journaltitleIEEE Journal on Selected Areas in Communications
titleReconfigurable intelligent surface aided NOMA networks

journaltitleIEEE Transactions on Communications
titlePower efficient IRS-assisted NOMA

journaltitlearXiv preprint arXiv:1906.02360
titleIntelligent reflecting surface assisted wireless communication: Modeling and channel estimation

journaltitleIEEE Wireless Communications Letters
titleCascaded channel estimation for large intelligent metasurface assisted massive MIMO

journaltitleIEEE access
titleEnabling large intelligent surfaces with compressive sensing and deep learning

booktitle2021 Joint European Conference on Networks and Communications & 6G Summit (EuCNC/6G Summit)
titleChannel estimation and hybrid architectures for RIS-assisted communications

journaltitleIEEE Journal of Selected Topics in Signal Processing
titleRIS-aided joint localization and synchronization with a single-antenna receiver: Beamforming design and low-complexity estimation

journaltitlearXiv preprint arXiv:1412.6980
titleAdam: A method for stochastic optimization

journaltitleIEEE Transactions on Signal Processing
titleA globally optimal energy-efficient power control framework and its efficient implementation in wireless interference networks

booktitleGLOBECOM 2020-2020 IEEE Global Communications Conference
titleResidual neural networks for digital predistortion

journaltitleIEEE Transactions on Information Theory
titleOptimum power allocation for parallel Gaussian channels with arbitrary input distributions

Curran Associates, Inc.

booktitleAdvances in Neural Information Processing Systems 32
titlePyTorch: An Imperative Style, High-Performance Deep Learning Library

journaltitleIEEE Communications Surveys & Tutorials
titleThe design and applications of high-performance ray-tracing simulation platform for 5G and beyond wireless communications: A tutorial

booktitle2008 IEEE International Conference on Communications
titleMMSE optimization with per-base-station power constraints for network MIMO systems

booktitle2012 IEEE 23rd International Symposium on Personal, Indoor and Mobile Radio Communications-(PIMRC)
titleBase station selection and per-cell codebook optimization for CoMP with joint processing

booktitle2011 IEEE 73rd Vehicular Technology Conference (VTC Spring)
titleVertical antenna tilt optimization for LTE base stations


# Roots of Polynomials and Umbilics of Surfaces

Brendan Guilfoyle  Brendan Guilfoyle
School of STEM
Munster Technological University
Tralee
Co. Kerry
Ireland<EMAIL_ADDRESS>and  Wilhelm Klingenberg  Wilhelm
Klingenberg
Department of Mathematical Sciences
University of Durham
Durham DH1 3LE
United Kingdom<EMAIL_ADDRESS>
###### Abstract.

For certain polynomials we relate the number of roots inside the unit circle
with the index of a non-degenerate isolated umbilic point on a real analytic
surface in Euclidean 3-space.

In particular, for $N>0$ we prove that for a certain ($N+2$)-real dimensional
family of complex polynomials of degree $N$, the number of roots inside the
unit circle is less than or equal to $1+N/2$. This bound is established as
follows.

From the polynomial we construct a convex real analytic surface containing an
isolated umbilic point, such that the index of the umbilic point is determined
by the number of roots of the polynomial that lie inside the unit circle. The
bound on the number of roots then follows from Hamburger’s bound on the index
of an isolated umbilic point on a convex real analytic surface.

###### Key words and phrases:

Zeros, complex polynomial, umbilic point, convex surface

###### 2020 Mathematics Subject Classification:

Primary 30C15; Secondary 53A05

##  1\. Introduction

Almost two centuries ago Cauchy established that all zeros of a complex
polynomial must lie inside a circle of a radius that can be determined by the
moduli of the polynomial’s coefficients [1]. Over the intervening years a vast
number of relations have been found between polynomial roots lying inside a
circle or annulus, and the coefficients of the polynomial [2] [9] [11].

The determination of the location of the zeros of a given polynomial relative
to the unit circle in the complex plane arises in questions of stability in
many areas of applied mathematics, including numerical integration, difference
equations, control theory and economic growth models [3] [10] [12] [13].

In this paper we present a relationship between the number of zeros contained
inside the unit circle of certain polynomials and the index of isolated
umbilic points on certain real analytic convex surfaces in ${\mathbb{R}}^{3}$.
A deep result of Hamburger [8] bounds the index of an isolated umbilic point
on any real analytic convex surface, which yields a bound on the number of
zeros inside the unit circle.

Umbilic points have many subtle properties - they can control the evolution of
surfaces under curvature flows [7], while Hamburger’s index bound does not
hold for even small perturbations of the ambient Euclidean metric [4].

The specific polynomials considered in this paper have relations between their
coefficients that arise from a Lagrangian condition in the space of oriented
affine lines in ${\mathbb{R}}^{3}$.

In more detail, let $P_{N}$ be a polynomial over ${\mathbb{C}}$ of degree $N$:

$\displaystyle P_{N}(\zeta)$ $\displaystyle=\prod_{n=1}^{N}(\zeta-\zeta_{n})$
$\displaystyle=\zeta^{N}+N\Delta_{1}^{N}\zeta^{N-1}+{\textstyle{\frac{N(N-1)}{2}}}\Delta_{2}^{N}\zeta^{N-2}+...+\;\binom{N}{n}\Delta_{n}^{N}\zeta^{N-n}+...+\Delta_{N}^{N},$

where the weighted symmetric functions of the roots
$\\{\zeta_{n}\\}_{n=1}^{N}$ are defined by

${\textstyle{\binom{N}{n}}}\Delta_{n}^{N}=\sum_{1\leq i_{1}<...<i_{n}\leq
N}\zeta_{i_{1}}\zeta_{i_{2}}...\zeta_{i_{n}}.$

Consider the following conditions on the coefficients of $P_{N}$:

(1)
$|\Delta_{N-2}^{N}|=1\qquad\qquad\Delta_{n}^{N}=\overline{\Delta_{N-2-n}^{N}}\Delta_{N-2}^{N}\qquad{\mbox{
for}}\;1\leq n\leq l-1,$

where $N=2l+1$ or $N=2l+2$. Note that in either case, these are $N-2$ real
conditions on the coefficients of $P_{N}$, which give rise to a real $N+2$
dimensional family of degree $N$ polynomials.

The second of conditions (1) is complex and does not depend solely on the
moduli of the polynomial coefficients. Moreover, our approach is complimentary
to Cauchy’s in that the arguments of the coefficients play a dominant role.

The main result of this paper is:

Main Theorem:

Let $P_{N}$ be a polynomial satisfying the relations (1) and suppose that none
of the roots $\zeta_{n}$ lies on the unit circle. Then the number of roots
(counted with multiplicity) of $P_{N}$ inside the unit circle is less than or
equal to $1+N/2$.

We prove this result by showing in Theorem 3.3 that, given a polynomial
$P_{N}$ satisfying conditions (1), there exists a convex real analytic surface
in ${\mathbb{R}}^{3}$ with an isolated umbilic point which has index
$I=K-N/2$, where $K$ is the number of zeros of $P_{N}$ lying inside the unit
circle. Thus, the local umbilic index bound of [8], which states that $I\leq
1$, yields the stated result.

In the next section we give the background on the construction of convex
surfaces, using Lagrangian surfaces in the space of oriented affine lines in
${\mathbb{R}}^{3}$ \- further details can be found in [5]. Section 3 contains
the proof of the Main Theorem and in the final section we illustrate our
techniques by constructing explicit examples.

##  2\. Real Analytic Convex Surfaces in ${\mathbb{R}}^{3}$

######  Definition 2.1.

Let ${\mathbb{L}}$ be the set of oriented (affine) lines in Euclidean 3-space
${\mathbb{R}}^{3}$.

######  Definition 2.2.

Let $\Phi:TS^{2}\rightarrow{\mathbb{L}}$ be the map that identifies
${\mathbb{L}}$ with the tangent bundle to the unit 2-sphere in Euclidean
${\mathbb{R}}^{3}$, by parallel translation. This bijection gives
${\mathbb{L}}$ the structure of a differentiable 4-manifold [6].

Let ($\xi,\eta$) be holomorphic coordinates on $TS^{2}$, where $\xi$ is
obtained by stereographic projection from the south pole onto the plane
through the equator, and we identify ($\xi,\eta$) with the vector

$\eta\frac{\partial}{\partial\xi}+\overline{\eta}\frac{\partial}{\partial\overline{\xi}}\in
T_{\xi}S^{2}.$

######  Proposition 2.3.

[5] The map $\Phi$ takes ($\xi,\eta$)$\in TS^{2}$ to the oriented line given
by

(2)
$z=\frac{2(\eta-\overline{\eta}\xi^{2})+2\xi(1+\xi\overline{\xi})r}{(1+\xi\overline{\xi})^{2}}$
(3)
$t=\frac{-2(\eta\overline{\xi}+\overline{\eta}\xi)+(1-\xi^{2}\overline{\xi}^{2})r}{(1+\xi\overline{\xi})^{2}},$

where $z=x^{1}+ix^{2}$, $t=x^{3}$, ($x^{1},x^{2},x^{3}$) are Euclidean
coordinates on ${\mathbb{R}}^{3}={\mathbb{C}}\oplus{\mathbb{R}}$ and $r$ is an
affine parameter along the line.

An oriented line congruence (or line congruence for short) is a surface
$\Sigma\subset{\mathbb{L}}$. An oriented surface $S$ in ${\mathbb{R}}^{3}$
gives rise to a line congruence by its oriented normal lines. Conversely a
line congruence is called Lagrangian if it is orthogonal to a (possibly
singular) surface in ${\mathbb{R}}^{3}$. The term refers to the canonical
symplectic structure on ${\mathbb{L}}$ which vanishes on such a surface - see
[6] for further details of the geometric background.

######  Proposition 2.4.

The oriented normals of a convex surface in ${\mathbb{R}}^{3}$ form a surface
$\Sigma\subset{\mathbb{L}}$ which is the graph of a local section of the
bundle $\pi:{\mathbb{L}}\rightarrow S^{2}$.

Such congruences are given by $\eta=F(\xi,\bar{\xi})$, where $F$ is a complex-
valued function on an open neighbourhood of $S^{2}$. For such local sections,
a necessary and sufficient condition for the line congruence to be Lagrangian
is:

######  Proposition 2.5.

The graph of a local section of $\pi:{\mathbb{L}}\rightarrow S^{2}$ is
Lagrangian iff the defining function $F$ satisfies

(4)
$\partial\left(\frac{F}{(1+\xi\bar{\xi})^{2}}\right)=\bar{\partial}\left(\frac{\bar{F}}{(1+\xi\bar{\xi})^{2}}\right),$

where, here and throughout, we denote
${\textstyle{\frac{\partial}{\partial\xi}}}$ by $\partial$. In particular,
there exists a real function $r:{\mathbb{C}}\rightarrow{\mathbb{R}}$, called
the support function, satisfying

(5) $\bar{\partial}r=\frac{2F}{(1+\xi\bar{\xi})^{2}}.$

Given a support function $r$ and the associated $F$ defined by (5), equations
(2) and (3) yield a parameterised surface in ${\mathbb{R}}^{3}$. Thus we
obtain a surface in 3-space parameterised by the inverse of its Gauss map.
Note that $r$ is defined up to an additive constant - replacing $r$ by $r+C$
moves the surface along its normal lines to the parallel surface.

For a convex real analytic surface the graph function $F$ is real analytic and
we can expand it in a power series about $\xi=0$:

(6) $F=\sum_{n,m=0}^{\infty}A_{nm}\xi^{n}\bar{\xi}^{m},$

for $A_{nm}\in{\mathbb{C}}$.

######  Proposition 2.6.

The Lagrangian conditions (4) restrict the coefficients of the power series by

(7) $A_{10}=\overline{A}_{10}$ (8)
$(n+1)A_{n+1\>0}=\overline{A}_{1n}-2\overline{A}_{0\;n-1}$ (9)
$(n+1)A_{n+1\>m}+(n-2)A_{n\>m-1}=(m+1)\overline{A}_{m+1\;n}+(m-2)\overline{A}_{m\;n-1},$

where $n,m\geq 1$.

###### Proof.

Suppose that

$F=\sum_{n,m=0}^{\infty}A_{nm}\xi^{n}\bar{\xi}^{m},$

and differentiating

$\partial F=\sum_{n,m=0}^{\infty}nA_{nm}\xi^{n-1}\bar{\xi}^{m}.$

Then

$\displaystyle(1+\xi\bar{\xi})\partial F-2\bar{\xi}F$
$\displaystyle=\sum_{n,m=0}^{\infty}nA_{nm}\xi^{n-1}\bar{\xi}^{m}+(n-2)A_{nm}\xi^{n}\bar{\xi}^{m+1}$
$\displaystyle=\sum_{n=1}^{\infty}nA_{n0}\xi^{n-1}+\sum_{m=1}^{\infty}A_{1m}\bar{\xi}^{m}+\sum_{n=2}^{\infty}\sum_{m=1}^{\infty}nA_{nm}\xi^{n-1}\bar{\xi}^{m}$
$\displaystyle\qquad-\sum_{m=1}^{\infty}2A_{0m-1}\bar{\xi}^{m}+\sum_{n,m=1}^{\infty}(n-2)A_{nm-1}\xi^{n}\bar{\xi}^{m}$
$\displaystyle=A_{10}+\sum_{m=1}^{\infty}(m+1)A_{m+1\>0}\xi^{m}+\sum_{m=1}^{\infty}(A_{1m}-2A_{0m-1})\bar{\xi}^{m}$
$\displaystyle\qquad+\sum_{n,m=1}^{\infty}[(n+1)A_{n+1m}+(n-2)A_{nm-1}]\xi^{n}\bar{\xi}^{m}.$

The result then follows from noting that the condition (4) is equivalent to

$(1+\xi\bar{\xi})\partial
F-2\bar{\xi}F=(1+\xi\bar{\xi})\bar{\partial}\bar{F}-2\xi\bar{F},$

and comparing terms in the above expression. ∎

We now rearrange these conditions into a more convenient form.

######  Proposition 2.7.

The Lagrangian conditions can be rewritten

$\begin{matrix}A_{nn-1}=\bar{A}_{nn-1}&{\mbox{I}}\\\
2(n-2)A_{0n-2}-2(n-1)A_{1n-1}+2nA_{2n}=n(n+1)\bar{A}_{n+1\;1}&II\\\
(n+2)A_{n+2\;2}+(n-1)A_{n+1\;1}=3\bar{A}_{3\;n+1}&III\\\
(n+2)A_{n+2\;m+1}+(n-1)A_{n+1\;m}=(m+2)\bar{A}_{m+2\;n+1}+(m-1)\bar{A}_{m+1\;n}&IV\\\
nA_{n0}=\bar{A}_{1n-1}-2\bar{A}_{0n-2}.&V\\\ \end{matrix},$

for $n>1,m<n$.

From here on, we assume that $S$ is a convex real analytic surface in
${\mathbb{R}}^{3}$ and the oriented normals are given by
$\xi\rightarrow(\xi,\eta=F(\xi,\bar{\xi}))$, with power series expansion (6)
satisfying I to V of Proposition 2.7.

We are interested in the umbilic points on $S$, that is, the points where the
second fundamental form of $S$ has a double eigenvalue. In our formalism,
these points can be characterised as follows:

######  Proposition 2.8.

[5] A point $p$ on $S$ is umbilic iff $\bar{\partial}F(\gamma)=0$, where
$\gamma$ is the oriented normal to $S$ at $p$.

Suppose that $S$ has an isolated umbilic point situated at the point with
$\xi=0$, and define the polar coordinates $(R,\theta)$ by $\xi=Re^{i\theta}$.

######  Definition 2.9.

The index $I$ of the isolated umbilic point at 0 is

$I=\lim_{R\rightarrow 0}{\textstyle{\frac{1}{8\pi
i}}}\int_{0}^{2\pi}\frac{\partial}{\partial\theta}\ln\left(\frac{\bar{\partial}F}{\partial\bar{F}}\right)d\theta.$

######  Example 2.10.

For the triaxial ellipsoid with semi-axes $a_{1}$, $a_{2}$ and $a_{3}$ the
support function is

$r=\left[a_{1}\left(\frac{\xi+\overline{\xi}}{1+\xi\overline{\xi}}\right)^{2}-a_{2}\left(\frac{\xi-\overline{\xi}}{1+\xi\overline{\xi}}\right)^{2}+a_{3}\left(\frac{1-\xi\overline{\xi}}{1+\xi\overline{\xi}}\right)^{2}\right]^{\scriptstyle{\frac{1}{2}}},$

and the Lagrangian section is

$F=\frac{a_{1}(\xi+\overline{\xi})(1-\xi^{2})+a_{2}(\xi-\overline{\xi})(1+\xi^{2})-2a_{3}\xi(1-\xi\overline{\xi})}{2\left[a_{1}(\xi+\overline{\xi})^{2}-a_{2}(\xi-\overline{\xi})^{2}+a_{3}(1-\xi\overline{\xi})^{2}\right]^{\scriptstyle{\frac{1}{2}}}}.$

For $a_{1}=a_{2}\neq a_{3}$ the surface is the rotationally symmetric
ellipsoid (the $x^{3}-$axis being the axis of symmetry). For
$a_{1}=a_{2}=a_{3}$ the surface is the round 2-sphere about the origin of
radius $\sqrt{a_{1}}$ which is totally umbilic.

For the rotationally symmetric ellipsoid $a_{1}=a_{2}=a\neq a_{3}$

$r=\frac{[4aR^{2}+a_{3}(1-R^{2})^{2}]^{\scriptstyle{\frac{1}{2}}}}{1+R^{2}}\qquad\qquad
F=\frac{(a-a_{3})R(1-R^{2})e^{i\theta}}{[4aR^{2}+a_{3}(1-R^{2})^{2}]^{\scriptstyle{\frac{1}{2}}}},$

in polar coordinates $\xi=Re^{i\theta}$.

Computing the derivative of $F$

$\bar{\partial}F=\frac{2a(a-a_{3})R^{2}(1+R^{2})e^{2i\theta}}{[4aR^{2}+a_{3}(1-R^{2})^{2}]^{\scriptstyle{\frac{3}{2}}}},$

and so the umbilics lie at $R=0$ and $R\rightarrow\infty$.

To compute the index of the umbilic at $R=0$ we find

$I=\lim_{R\rightarrow 0}{\textstyle{\frac{1}{8\pi
i}}}\int_{0}^{2\pi}\frac{\partial}{\partial\theta}\ln\left(\frac{\bar{\partial}F}{\partial\bar{F}}\right)d\theta={\textstyle{\frac{1}{2\pi}}}\int_{0}^{2\pi}d\theta=1.$

On the other hand, for the triaxial ellipsoid $a_{1}\neq a_{2}\neq a_{3}$, a
similar calculation shows that there are 4 umbilic points which lie at
$\theta=0$ and

$R=\pm\left[-{{\frac{(a_{1}+a_{2})a_{3}-2a_{1}a_{2}}{(a_{1}-a_{2})a_{3}}}}\pm\left[\left({{\frac{(a_{1}+a_{2})a_{3}-2a_{1}a_{2}}{(a_{1}-a_{2})a_{3}}}}\right)^{2}-1\right]^{\scriptstyle{\frac{1}{2}}}\right]^{\scriptstyle{\frac{1}{2}}}.$

Denote either of these values by $R=R_{0}$. A rotation

$\xi\rightarrow\frac{\xi-R_{0}}{1+R_{0}\xi}\qquad\qquad
F\rightarrow\frac{(1+R_{0}^{2})F}{(1+R_{0}\xi)^{2}},$

moves this umbilic to the origin $R=0$ and a lengthy computation yields

$I=\lim_{R\rightarrow 0}{\textstyle{\frac{1}{8\pi
i}}}\int_{0}^{2\pi}\frac{\partial}{\partial\theta}\ln\left(\frac{\bar{\partial}F}{\partial\bar{F}}\right)d\theta={\textstyle{\frac{1}{4\pi}}}\int_{0}^{2\pi}d\theta={\textstyle{\frac{1}{2}}}.$

For any section $F$, define the complex-valued functions $G_{k}(\theta)$ by
the polar decomposition of $\bar{\partial}F$:

$\bar{\partial}F=\sum_{k=0}^{\infty}G_{k}(\theta)e^{-ik\theta}R^{k}.$

Then

$G_{k}(\theta)=\sum_{n=0}^{k}(k-n+1)A_{n\;k-n+1}e^{2in\theta},$

and, in particular, $G_{0}=A_{01}=\bar{\partial}F(0)$, which vanishes iff
$\xi=0$ is an umbilic point of $S$.

######  Definition 2.11.

An isolated umbilic point at $\xi=0$ is said to be non-degenerate of order $N$
if

$G_{k}(\theta)=0\qquad{\mbox{for}}\;k<N\qquad\qquad{\mbox{and}}\qquad\qquad
G_{N}(\theta)\neq 0\qquad{\mbox{for all}}\;\theta.$

For such an umbilic point, define the degree $N$ polynomial

$P_{N}=\sum_{n=0}^{N}(N-n+1)\frac{A_{n\;N-n+1}}{A_{N1}}\zeta^{n}.$

Here and throughout we assume that $A_{N1}\neq 0$.

######  Proposition 2.12.

For a non-degenerate isolated umbilic point of order N, the polynomial $P_{N}$
does not have any zeros on the unit circle $|\zeta|=1$ and the index of the
umbilic point is

$I=K-N/2,$

where $K$ is the number of zeros of $P_{N}$ that lie inside the unit circle.

###### Proof.

On the unit circle $\zeta=e^{2i\theta}$ we have

$P_{N}=\frac{G_{N}(\theta)}{A_{N1}}\neq 0,$

which proves the first statement. For the second statement

$\displaystyle I$ $\displaystyle=\lim_{R\rightarrow
0}{\textstyle{\frac{1}{8\pi
i}}}\int_{0}^{2\pi}\frac{\partial}{\partial\theta}\ln\left(\frac{\bar{\partial}F}{\partial\bar{F}}\right)d\theta$
$\displaystyle={\textstyle{\frac{1}{8\pi
i}}}\int_{0}^{2\pi}\frac{\partial}{\partial\theta}\ln\left(\frac{G_{N}(\theta)e^{-iN\theta}}{\bar{G}_{N}(\theta)e^{iN\theta}}\right)d\theta$
$\displaystyle={\textstyle{\frac{1}{8\pi
i}}}\int_{0}^{2\pi}\frac{\partial}{\partial\theta}\ln\left(\frac{P_{N}}{\bar{P}_{N}}\right)d\theta-N/2$
$\displaystyle=K-N/2,$

as claimed. ∎

However, we also have

######  Theorem 2.13.

[8] The index of an isolated umbilic point on a convex surface in
${\mathbb{R}}^{3}$ is less than or equal to 1.

As a corollary

######  Corollary 2.14.

For a non-degenerate umbilic point of order $N$ on a real analytic surface,
the number of zeros of the polynomial $P_{N}$ that lie inside the unit circle
is less than or equal to $1+N/2$.

Of course the polynomial $P_{N}$ is not arbitrary: its coefficients satisfy
certain relations which are contained in the Lagrangian conditions I to V of
Proposition 2.7. To prove our Main Theorem we reverse the preceding argument:
given a polynomial $P_{N}$ satisfying the order $N$ Lagrangian conditions, we
construct a convex surface in ${\mathbb{R}}^{3}$ with isolated umbilic point
for which $P_{N}$ is the lowest order term in $R$ of $\bar{\partial}F$. This
we do in the next section.

##  3\. Proof of the Main Theorem

Our first task is to determine the restrictions placed on the coefficients of
the polynomial $P_{N}$ by the Lagrangian conditions I to V. To understand
these conditions, define the Fourier decomposition of $\bar{\partial}F$ about
the origin by

$\bar{\partial}F=\sum_{k=-\infty}^{\infty}\gamma_{k}(R)e^{ik\theta},$

where

$\gamma_{k}(R)=\frac{1}{2\pi}\int_{0}^{2\pi}\bar{\partial}F\;e^{-ik\theta}d\theta.$

In Figure 1 we decompose each of these Fourier coefficients in powers of $R$
and indicate the Lagrangian conditions placed on the coefficients $A_{nm}$ by
enclosing related terms by a box. The four Lagrangian conditions I to IV each
has a distinctive shape which enables us to trace this net of relations and
assist in calculations.

Condition V does not introduce any relationship on the coefficients of
$\bar{\partial}F$ as it involves $A_{n0}$, which is annihilated by
$\bar{\partial}$. As such, we can view condition V as defining the
coefficients $A_{n0}$ in the full power series after we have determined the
other coefficients.

######  Proposition 3.1.

For the polynomial $P_{N}$ associated with a non-degenerate umbilic point of
order $N$ the Lagrangian conditions I to IV mean that the coefficients must
satisfy

(10) $(n+2)A_{n+2\;N-n-1}=(N-n)\bar{A}_{N-n\;n+1}\qquad{\mbox{for}}\;1\leq
n\leq l-1$ (11) $2A_{2\;N-1}=N\bar{A}_{N1},$

where $N=2l+1$ or $N=2l+2$.

###### Proof.

Given $N$, the $R^{N}$ column of Figure 1 lists the lowerest order non-zero
coefficients of $\bar{\partial}F$, namely those that enter into the definition
of $P_{N}$. The Lagrangian condition that these coefficients must satisfy are
indicated by the boxes in which they lie. Equations (10) and (11) are obtained
by inspection of the boxes that contain only these non-zero terms and lower
order terms, which by assumption vanish. ∎

We now show that these relations are precisely the equations (1) given in the
Introduction.

######  Corollary 3.2.

For the polynomial $P_{N}$ associated with a non-degenerate umbilic point of
order $N$ the Lagrangian conditions I to IV mean that the symmetric functions
$\Delta^{N}_{n}$ must satisfy

$|\Delta_{N-2}^{N}|=1\qquad\qquad\Delta_{n}^{N}=\Delta_{N-2}^{N}\overline{\Delta_{N-2-n}^{N}}\qquad{\mbox{
for}}\;1\leq n\leq l-1,$

where $N=2l+1$ or $N=2l+2$.

###### Proof.

We have the definitions

$P_{N}=\sum_{k=0}^{N}(k+1)\frac{A_{N-k\;k+1}}{A_{N1}}\zeta^{N-k}=\sum_{k=0}^{N}\binom{N}{k}\;\Delta_{k}^{N}\zeta^{N-k},$

and so we have the relations

(12) $\binom{N}{k}\Delta_{k}^{N}=(k+1)\frac{A_{N-k\;k+1}}{A_{N1}}.$

For $k=N-2$ this reads

${\textstyle{\frac{N(N-1)}{2}}}\Delta^{N}_{N-2}=(N-1)\frac{A_{2\;N-1}}{A_{N1}}.$

Thus, with the aid of equation (11) we have

$\Delta^{N}_{N-2}=\frac{2A_{2\;N-1}}{NA_{N1}}=\frac{\bar{A}_{N1}}{A_{N1}}.$

Therefore $|\Delta^{N}_{N-2}|=1$, as claimed.

On the other hand, utilising equation (10) we deduce that

$\binom{N}{k}\Delta_{k}^{N}=(k+1)\frac{A_{N-k\;k+1}}{A_{N1}}={\textstyle{\frac{(k+1)(k+2)}{N-k}}}\frac{A_{k+2\;N-k-1}}{A_{N1}}.$

Now applying equation (12) again with $k$ replaced by $N-k-2$, we get

$\binom{N}{k}\Delta_{k}^{N}={\textstyle{\frac{(k+1)(k+2)}{(N-k)(N-k-1)}}}\binom{N}{N-k-2}\frac{\bar{A}_{N1}}{A_{N1}}\bar{\Delta}^{N}_{N-k-2},$

which, after some cancellation, reads

$\Delta_{k}^{N}=\Delta^{N}_{N-2}\overline{\Delta_{N-k-2}^{N}},$

as claimed. ∎

Now, given a polynomial $P_{N}$ satisfying relations (1) we seek to construct
a convex real analytic surface in ${\mathbb{R}}^{3}$ with umbilic point at
$\xi=0$. We do this by first using the coefficients defined by $P_{N}$ as
terms of order $N$ in the power series of $\bar{\partial}F$, and then adding
further terms of higher order in such a way that the full Lagrangian condition
holds. By moving to higher order terms we ensure that the index of the umbilic
point resides in the lowest order term, namely $P_{N}$.

Such an extension will be far from unique, but we seek to extend in such a way
that minimises the number of non-zero terms introduced. Finally, we must
ensure that the resulting umbilic point at $\xi=0$ is indeed isolated, and
hence we can apply Hamburger’s result to prove the stated Theorem.

Once again Figure 1 is useful as a guide to the equations we are dealing with.
Suppose then that we have non-zero terms in the $R^{N}$ column of
$\bar{\partial}F$ (namely, those given by $P_{N}$) and all terms of lower
degree are zero. We will see that it is only necessary to go to terms of order
$R^{N+2}$ to complete the Lagrangian conditions. In this we identify 6
situations of overlapping equations that must dealt with separately.

Case 1: IV IV IV

Thus we have the Lagrangian condition

$(n+2)A_{n+2\;m+1}=(m+2)\bar{A}_{m+2\;n+1},$

and the two higher order equations

$(n+3)A_{n+3\;m+2}+nA_{n+2\;m+1}=(m+3)\bar{A}_{m+3\;n+2}+m\bar{A}_{m+2\;n+1}$
$(n+1)A_{n+3\;m+2}=(m+1)\bar{A}_{m+3\;n+2}.$

The first and third of these, rearranged to

$A_{m+2\;n+1}=\frac{n+2}{m+2}\bar{A}_{n+2\;m+1}\qquad\qquad
A_{m+3\;n+2}=\frac{n+1}{m+1}\bar{A}_{n+3\;m+2},$

and substituted in the second, yield

$A_{n+3\;m+2}=\frac{m+1}{m+2}A_{n+2\;m+1},$

and so

$A_{m+3\;n+2}=\frac{n+1}{m+2}\bar{A}_{n+2\;m+1}.$

With these choices the power series is truncated.

Case 2: III IV IV

First we have the Lagrangian condition at order $N$:

$(N-1)A_{N-1\;2}=3\bar{A}_{3\;N-2},$

and two order $N+2$ conditions:

$NA_{N\;3}+(N-3)A_{N-1\;2}=4\bar{A}_{4\;N-1}+\bar{A}_{3\;N-2}$
$(N-2)A_{N\;3}=2\bar{A}_{4\;N-1}.$

Combining these as before we get

$A_{N-1\;2}=\frac{3}{N-1}\bar{A}_{3\;N-2}\qquad
A_{4\;N-1}=\frac{N-2}{N-1}A_{3\;N-2}\qquad
A_{N\;3}=\frac{2}{N-1}\bar{A}_{3\;N-2}.$

Case 3: II III IV

The Lagrangian condition at order $N$ is

$2A_{2\;N-1}=N\bar{A}_{N\;1},$

and the other conditions are

$(N+1)A_{N+1\;2}+(N-2)A_{N\;1}=3\bar{A}_{3\;N}$
$(N-1)A_{N+1\;2}=\bar{A}_{3\;N}.$

Solving these we have

$A_{N\;1}=\frac{2}{N}\bar{A}_{2\;N-1}\qquad
A_{3\;N}=\frac{N-1}{N}A_{2\;N-1}\qquad
A_{N+1\;2}=\frac{1}{N}\bar{A}_{2\;N-1}.$

Case 4: II

Here we have exhausted the Lagrangian condition, but must accommodate the tail
of the polynomial $P_{N}$. That is

$-NA_{1\;N}+(N+1)A_{2\;N+1}=0,$

or

$A_{2\;N+1}=\frac{N}{N+1}A_{1\;N}.$

The other term is

$(N+1)A_{0\;N+1}-(N+2)A_{1\;N+2}=0,$

or

$A_{1\;N+2}=\frac{N+1}{N+2}A_{0\;N+1}.$

Case 5: V

As mentioned earlier, this Lagrangian condition does not enter into the
relations enjoyed by the coefficients of $\bar{\partial}F$, but is taken as a
definition of the non-vanishing holomorphic part of $F$. In particular, we
have

$A_{N+1\;0}=\frac{1}{N+1}\bar{A}_{1\;N},$

and

$(N+3)A_{N+3\;0}=\bar{A}_{1\;N+2}-2\bar{A}_{0\;N+1}=\frac{N+1}{N+2}\bar{A}_{0\;N+1}-2\bar{A}_{0\;N+1},$

or

$A_{N+3\;0}=-\frac{1}{N+2}\bar{A}_{0\;N+1}.$

Case 6: I

When $N=2l+2$ is even, we have the Lagrangian condition

$A_{l+2\;l+1}=\bar{A}_{l+2\;l+1}.$

As this involves only one coefficient of the power series it is not strictly
necessary to relate it to any other term to maintain the Lagrangian condition.
However, in order to ensure that the umbilic on the surface is isolated, we
impose the final following condition

$A_{l+3\;l+2}=\frac{l}{l+2}A_{l+2\;l+1}.$

Combining all of the preceding results

######  Theorem 3.3.

Given a polynomial $P_{N}$ satisfying the relations (1), choose
$A_{N1}\in{\mathbb{C}}^{*}$ and define the complex numbers
$\\{A_{n\;N-n+1}\\}_{n=0}^{N-1}$ by

$A_{N-n\;n+1}={\textstyle{\frac{1}{n+1}}}\binom{N}{n}A_{N1}\Delta_{n}^{N},$

for $n=1,2,...N$.

Then, for $N=2l+1$ the following function defines a convex real analytic
surface $S$ in ${\mathbb{R}}^{3}$:

$F=\sum_{n=0}^{l+1}A_{n\;N-n+1}(1+{\textstyle{\frac{N-n+1}{N-n+2}}}\xi\bar{\xi})\xi^{n}\bar{\xi}^{N-n+1}+\bar{A}_{n\;N-n+1}({\textstyle{\frac{n}{N-n+2}}}+{\textstyle{\frac{n-1}{N-n+2}}}\xi\bar{\xi})\xi^{N-n+2}\bar{\xi}^{n-1},$

with support function

$r=\frac{2}{1+\xi\bar{\xi}}\sum_{n=0}^{l+1}\left({\textstyle{\frac{A_{n\;N-n+1}}{N-n+2}}}\xi^{n}\bar{\xi}^{N-n+2}+{\textstyle{\frac{\bar{A}_{n\;N-n+1}}{N-n+2}}}\xi^{N-n+2}\bar{\xi}^{n}\right)+C.$

For $N=2l+2$ the function is

$\displaystyle F$
$\displaystyle=\sum_{n=0}^{l+1}A_{n\;N-n+1}(1+{\textstyle{\frac{N-n+1}{N-n+2}}}\xi\bar{\xi})\xi^{n}\bar{\xi}^{N-n+1}+\bar{A}_{n\;N-n+1}({\textstyle{\frac{n}{N-n+2}}}+{\textstyle{\frac{n-1}{N-n+2}}}\xi\bar{\xi})\xi^{N-n+2}\bar{\xi}^{n-1}$
$\displaystyle\qquad\qquad\qquad+A_{l+2\;l+1}(1+{\textstyle{\frac{l+1}{l+2}}}\xi\bar{\xi})\xi^{l+2}\bar{\xi}^{l+1},$

and the support function is

$\displaystyle r$
$\displaystyle=\frac{2}{1+\xi\bar{\xi}}\sum_{n=0}^{l+1}\left({\textstyle{\frac{A_{n\;N-n+1}}{N-n+2}}}\xi^{n}\bar{\xi}^{N-n+2}+{\textstyle{\frac{\bar{A}_{n\;N-n+1}}{N-n+2}}}\xi^{N-n+2}\bar{\xi}^{n}\right)$
$\displaystyle\qquad\qquad\qquad+\frac{2}{(1+\xi\bar{\xi})(l+2)}A_{l+2\;l+1}\xi^{l+2}\bar{\xi}^{l+2}+C.$

In both cases, $S$ has an isolated umbilic point at $\xi=0$ which is non-
degenerate of order $N$.

###### Proof.

We will show that $F$ and $r$ satisfy (5), which implies the Lagrangian
condition (4). For $N=2l+1$ we compute:

$\displaystyle\bar{\partial}r$
$\displaystyle=\frac{2}{1+\xi\bar{\xi}}\sum_{n=0}^{l+1}\left(A_{n\;N-n+1}\xi^{n}\bar{\xi}^{N-n+1}+{\textstyle{\frac{n}{N-n+1}}}\bar{A}_{n\;N-n+1}\xi^{N-n+2}\bar{\xi}^{n-1}\right)$
$\displaystyle\qquad\qquad\qquad\qquad-\frac{2\xi}{(1+\xi\bar{\xi})^{2}}\sum_{n=0}^{l+1}\left(\frac{A_{n\;N-n+1}}{N-n+2}\xi^{n}\bar{\xi}^{N-n+2}+\frac{\bar{A}_{n\;N-n+1}}{N-n+2}\xi^{N-n+2}\bar{\xi}^{n}\right)$
$\displaystyle=\frac{2}{(1+\xi\bar{\xi})^{2}}\sum_{n=0}^{l+1}\left[A_{n\;N-n+1}\left(1+{\textstyle{\frac{N-n+1}{N-n+2}}}\xi\bar{\xi}\right)\xi^{n}\bar{\xi}^{N-n+1}\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\left.+\bar{A}_{n\;N-n+1}\left({\textstyle{\frac{n}{N-n+2}}}+{\textstyle{\frac{n-1}{N-n+2}}}\xi\bar{\xi}\right)\xi^{N-n+2}\bar{\xi}^{n-1}\right]$
$\displaystyle=\frac{2F}{(1+\xi\bar{\xi})^{2}},$

as claimed.

The proof for $N=2l+2$ is the same with the added computation

$\bar{\partial}\left(\frac{2A_{l+2\;l+1}\xi^{l+2}}{(1+\xi\bar{\xi})(l+2)}\bar{\xi}^{l+2}\right)=\frac{2{A}_{l+2\;l+1}}{(1+\xi\bar{\xi})^{2}}(1+{\textstyle{\frac{l+1}{l+2}}}\xi\bar{\xi})\xi^{l+2}\bar{\xi}^{l+1}.$

For $N>0$ clearly $\xi=0$ is an umbilic point. To show that it is isolated we
compute

$\displaystyle\bar{\partial}F$
$\displaystyle=(1+\xi\bar{\xi})\sum_{n=-1}^{l}(N-n)A_{n+1\;N-n}\xi^{n+1}\bar{\xi}^{N-n-1}+{\textstyle{\frac{n(n+1)}{N-n+1}}}\bar{A}_{n+1\;N-n}\xi^{N-n+1}\bar{\xi}^{n-1}$
$\displaystyle=(1+\xi\bar{\xi})\sum_{n=0}^{N}(N-n+1)A_{n\;N-n+1}\xi^{n}\bar{\xi}^{N-n}$
$\displaystyle=(1+R^{2})G_{N}(\theta)e^{-iN\theta}R^{N}.$

By assumption, the polynomial $P_{N}$ has no zeros on the unit circle, and so
$G_{N}(\theta)\neq 0$. Thus $R=0$ is an isolated umbilic point and is non-
degenerate of order $N$. ∎

##  4\. Examples

We now demonstrate direct proofs of the Main Theorem for orders $N=1,2,3,4$.
For order $N\geq 5$ there is no direct method.

N=1,2:

In these cases the statement is easily seen to be trivially true.

N=3:

Let $\zeta_{1},\zeta_{2}$ and $\zeta_{3}$ be the roots of the cubic polynomial
satisfying the single Lagrangian condition

$|\Delta_{1}^{3}|=1.$

Thus

${\textstyle{\frac{1}{3}}}|\zeta_{1}+\zeta_{2}+\zeta_{3}|=1$

and we conclude that $|\zeta_{i}|\geq 1$ for some $i$. Thus the number of
roots inside the unit circle must be less than or equal to 2, as claimed.

N=4:

Let $\zeta_{1},\zeta_{2}$,$\zeta_{3}$ and $\zeta_{4}$ be the roots of the
quartic polynomial satisfying the Lagrangian conditions

$|\Delta_{2}^{4}|=1\qquad\qquad\Delta_{1}^{4}=\Delta_{2}^{4}\bar{\Delta}_{1}^{4}.$

The first of these conditions implies that

${\textstyle{\frac{1}{6}}}|\zeta_{1}\zeta_{2}+\zeta_{1}\zeta_{3}+\zeta_{1}\zeta_{4}+\zeta_{2}\zeta_{3}+\zeta_{2}\zeta_{4}+\zeta_{3}\zeta_{4}|=1,$

and we conclude that $|\zeta_{i}\zeta_{j}|\geq 1$ for some $i$ and $j$. This
means that $|\zeta_{i}|\geq 1$ for some $i$, and hence the number of roots
inside the unit circle must be less than or equal to 3, as claimed.

N=5:

In the case of the quintic we are considering (after substitution of the
Lagrangian conditions) the family of polynomials:

$P_{5}=\zeta^{5}+5\alpha
e^{iA}\zeta^{4}+10\bar{\alpha}\zeta^{3}+10e^{iA}\zeta^{2}+\beta\zeta+\gamma,$

for $\alpha,\beta,\gamma\in{\mathbb{C}}$ and $A\in{\mathbb{R}}$.

The Main Theorem states that such a polynomial has at most 3 roots inside the
unit circle. By elementary arguments, as in the previous cases, it is not hard
to show that there are at most 4 roots inside the unit circle, but reducing
this to 3 roots is extremely difficult by elementary methods.

Acknowledgements: The authors would like to thank Gerhard Schmeisser for
helpful discussions during the original evolution of this work.

## References

  * [1] A.L. Cauchy, Exercices de mathématique, Œuvres 2.9 (1829) 122.
  * [2] B. Datt and N.K. Govil, On the location of the zeros of a polynomial, J. Approx. Theory 24.1 (1978) 78–82.
  * [3] R.J. Duffin, Algorithms for classical stability problems, Siam Rev. 11.2 (1969) 196–213.
  * [4] B. Guilfoyle, On isolated umbilic points, Comm. Anal. Geom. 28.8 (2020) 2005–2018.
  * [5] B. Guilfoyle and W. Klingenberg, Generalised surfaces in ${\mathbb{R}}^{3}$, Math. Proc. R. Ir. Acad. 104A (2004) 199–209.
  * [6] B. Guilfoyle and W. Klingenberg, An indefinite Kähler metric on the space of oriented lines, J. London Math. Soc. 72.2 (2005) 497–509.
  * [7] B. Guilfoyle and W. Klingenberg, Evolving to non-round Weingarten spheres: integer linear Hopf flows, Partial Differ. Equ. Appl. 2 (2021) 72.
  * [8] H. Hamburger, Beweis einer Carathéodory Vermutung I II III, Ann. Math. 41 (1940) 63–86, Acta. Math. 73 (1941) 175–228 and 229–332.
  * [9] P. Kumar and R. Dhankhar, On the location of zeros of polynomials, Complex Anal. Oper. Theory 16.1 (2022) 1–13.
  * [10] J.J.H. Miller, On the location of zeros of certain classes of polynomials with applications to numerical analysis, IMA J. Appl. Math. 8.3 (1971) 397–406.
  * [11] Q.I. Rahman and G. Schmeisser, Analytic theory of polynomials, London Mathematical Society Monographs, New Series 26, The Clarendon Press, Oxford University Press, Oxford, (2002).
  * [12] R.D. Richtmyer and K.W. Morton, Difference methods for initial-value problems, Malabar (1994).
  * [13] P.A. Samuelson, Conditions that the roots of a polynomial be less than unity in absolute value, Annals of Mathematical Statistics 12.3 (1941) 360–364.