\section{Introduction}

Non-Local Means (NLM) is a data-driven diffusion mechanism that was introduced by Buades et al. in \cite{BCM2005}. It has proved to be a simple yet powerful method for image denoising. In this method, a given pixel is denoised
using a weighted average of other pixels in the (noisy) image. In particular, given a noisy image $u = (u_i)$, the denoised image $\hat{u} = (\hat{u}_i)$ at pixel $i$ is computed using the formula
\begin{equation}
\label{NLM_formula}
 \hat{u}_i = \frac{\sum_{j} w_{ij} u_j}{\sum_{j} w_{ij} },
\end{equation}
where $w_{ij}$ is some weight (affinity) assigned to pixels $i$ and $j$. The sum in \eqref{NLM_formula} is ideally performed over the whole image. In practice, however, one restricts $j$ to a geometric neighborhood of $i$, e.g., to a
sufficiently large window of size $S \times S$ \cite{BCM2005}.

The central idea in \cite{BCM2005} was to set the weights using image patches centered around each pixel, namely as
\begin{equation}
\label{Weights}
 w_{ij} = \exp\Big(- \frac{1}{h^2}\lVert \mathbf{P}_i - \mathbf{P}_j \rVert^2  \Big),
\end{equation}
where $\mathbf{P}_i$ and $\mathbf{P}_j$ are the image patches of size $k \times k$ centered at pixels $i$ and $j$. Here, $\lVert \mathbf{P} \rVert$ is the Euclidean norm of patch $\mathbf{P}$ as a point in $\mathbf{R}^{k^2}$, and $h$ is a smoothing parameter.
Thus, pixels with similar neighborhoods are given larger weights compared to pixels whose neighborhoods look different. The algorithm makes explicit use of the fact that repetitive patterns appear in most natural images.
It is remarkable that the simple formula in \eqref{NLM_formula} often provides state-of-the-art results in image denoising \cite{BCM2010}. One outstanding feature of NLM is that, in comparison to other denoising techniques such as
Gaussian smoothing, anisotropic diffusion, total variation denoising, and wavelet regularization, the so-called  method noise (difference of the denoised and the noisy image) in NLM appears more
like white noise \cite{BCM2005,BCM2010}. We refer the reader to \cite{BCM2010} for a detailed review of the algorithm.

The rest of the letter is organized as follows. In Section \ref{S2}, we explain how the denoising performance of NLM can be improved in the vicinity of edges using the Euclidean median. Based on this observation, we propose a
new denoising algorithm in Section \ref{S3} and discuss its implementation. This is followed by some preliminary denoising results in Section \ref{S4}, where we compare our new algorithm with NLM. We conclude with some remarks
in Section \ref{S5}.

\section{Better robustness using Euclidean median}
\label{S2}

The denoising performance of NLM depends entirely on the reliability of the weights $w_{ij}$. The weights are, however, computed from the noisy image and not the clean image. Noise affects the distribution of weights,
particularly when the noise is large. By noise, we will always mean zero-centered white Gaussian noise with variance $\sigma^2$ in the rest of the discussion.

\begin{figure}[!htp]
  \centering
  \subfloat[Clean and noisy edge.]{\label{fig1}\includegraphics[width=0.7\linewidth]{fig1.eps}} \\
  \subfloat[Weights.]{\label{fig2}\includegraphics[width=0.7\linewidth]{fig2.eps}}
  \caption{Ideal edge in one dimension.}
    \label{figI}
\end{figure}

To understand the effect of noise and to motivate the main idea of the paper, we perform a simple experiment. We consider the particular case where the pixel of interest is close to an (ideal) edge. For this, we take a clean
edge profile of unit height and add noise ($\sigma = 0.2$) to it. This is shown in Fig. \ref{figI}. We now select a point of interest a few pixels to the right of the edge (marked with a star). The goal is to estimate its true value from its neighboring points using NLM.
To do so, we take 3-sample patches around each point ($k=3$), and a search window of $S=41$. The patches are shown as points in $3$-dimensions in Fig \ref{fig3}. The clean patch for the point of interest
is at $(1,1,1)$.

We now use \eqref{Weights} to compute the weights, where we set $h = 10\sigma$. The weights corresponding to the point of interest are shown in Fig. \ref{fig2}.
Using the noisy weights, we obtain an estimate of around $0.65$. This estimate has a geometric interpretation. It is the center coordinate of the Euclidean mean $\sum_{j} w_j \mathbf{P}_j /\sum_{j} w_j$,
where $w_j$ are the weights in Fig. \ref{fig2}, and $\mathbf{P}_j$ are the patches in Fig. \ref{fig3}. The Euclidean mean is marked with a purple diamond in Fig. \ref{figII}. Note that
the patches drawn from the search window are clustered around the centers $A=(0,0,0)$ and $B=(1,1,1)$. For the point of interest, the points around $A$ are the outliers, while the ones around $B$ are the inliers.
The noise reduces the relative gap in the corresponding weights, causing the mean to drift away from the inliers toward the outliers.

\begin{figure}[!htp]
  \centering
  \subfloat[3d patch space close to the edge.]{\label{fig3}\includegraphics[width=0.7\linewidth]{fig3.eps}}   \\
  \subfloat[2d projection (first $2$ coordinates).]{\label{fig4}\includegraphics[width=0.7\linewidth]{fig4.eps}}
  \caption{Outlier model of the patch space for the point of interest in Fig. \ref{fig1}. Due to its robustness to outliers, the Euclidean median behaves as a better estimator of the clean patch than the
  Euclidean mean (see text for explanation).}
  \label{figII}
\end{figure}

Note that the Euclidean mean is the minimizer of $\sum_{j} w_j \lVert \mathbf{P} - \mathbf{P}_j \rVert^2$ over all patches $\mathbf{P}$. Our main observation is that, if we instead compute the minimizer of $\sum_{j} w_j \lVert \mathbf{P} - \mathbf{P}_j \rVert$ over
all $\mathbf{P}$, and take the center coordinate of $\mathbf{P}$, then we get a much better estimate. Indeed, the denoised value turns out to be around $0.92$ in this case.
The above minimizer is called the Euclidean median (or the geometric median) in the literature \cite{HR2009}. We will often simply call it the median.
We repeated the above experiment using several noise realizations, and consistently got better results using the median. Averaged over $10$ trials, the denoised value using the mean and median were found to be $0.62$ and $0.93$, respectively.
Indeed, in Fig. \ref{figII}, note how close the median is to the true patch compared to the mean. This does not come as a surprise since it is well-known that the median is more robust to outliers than the mean. This fact has
a rather simple explanation in one dimension. In higher dimensions, note that the former is minimizer of (the square of) the weighted $\ell^2$ norm of the distances $\lVert \mathbf{P} - \mathbf{P}_j \rVert$, while the latter is the
minimizer of the weighted $\ell^1$ norm of these distances. It is this use of the $\ell^1$ norm over the $\ell^2$ norm that makes the Euclidean median more robust to outliers \cite{HR2009}.

\section{Non-Local Euclidean Medians}
\label{S3}

Following the above observation, we propose Algorithm \ref{algo1} below which we call the Non-Local Euclidean Medians (NLEM). We use $S(i)$ to denote the search window of size $S \times S$ centered at pixel $i$.

\begin{algorithm}
\caption{Non-Local Euclidean Medians}
\label{algo1}
\begin{algorithmic}
     \State \textbf{Input}: Noisy image $u = (u_i)$, and parameters $h, \lambda, S, k$.
      \State \textbf{Return}: Denoised image $\hat{u} = (\hat{u}_i)$.
     \State (1) Estimate noise variance $\sigma^2$, and set $h = \lambda \sigma$
     \State (2) Extract patch $\mathbf{P}_i \in \mathbf{R}^{k^2}$ at every pixel $i$.
     \State (3) For every pixel $i$, do
       \State \hspace{3mm}   (a)   Set $w_{ij} = \exp(-\lVert \mathbf{P}_i - \mathbf{P}_j \rVert^2/h^2)$ for every $j \in S(i)$.
        \State \hspace{3mm}   (b)  Find patch $\mathbf{P}$ that minimizes $\sum_{j \in S(i)} w_{ij} \lVert \mathbf{P} - \mathbf{P}_j \rVert$.
        \State \hspace{3mm}   (c)  Assign $\hat{u}_i$ the value of the center pixel in $\mathbf{P}$.
\end{algorithmic}
\end{algorithm}

The difference with NLM is in step $(3)$b which involves the computation of the Euclidean median. That is, given points $\boldsymbol{x}_1,\ldots,\boldsymbol{x}_n \in \mathbf{R}^d$ and weights $w_1,\ldots,w_n$, we are required to find $\boldsymbol{x} \in  \mathbf{R}^d$ that
minimizes the convex cost $\sum_{j=1}^n  w_j \lVert \boldsymbol{x} - \boldsymbol{x}_j \rVert$. There exists an extensive literature on the computation of the Euclidean median; see \cite{W1937,XY1997}, and the references therein. The simplest algorithm in this area is the
so-called Weiszfeld algorithm \cite{W1937,XY1997}. This is, in fact, based on the method of iteratively reweighted least squares (IRLS), which has received renewed interest in the compressed sensing
community in the context of $\ell^1$  minimization \cite{CW2008,DDFG2009}. Starting from an estimate $\boldsymbol{x}^{(k)}$, the idea is to set the next iterate as
\begin{equation*}
\boldsymbol{x}^{(k+1)} = \arg\ \min_{\boldsymbol{x} \in \mathbf{R}^d} \sum_{j=1}^n  w_j \frac{\lVert \boldsymbol{x} - \boldsymbol{x}_j \rVert^2}{\lVert \boldsymbol{x}^{(k)} - \boldsymbol{x}_j \rVert}.
\end{equation*}
This is a least-squares problem, and the minimizer is given by
\begin{equation}
\label{IRLS}
\boldsymbol{x}^{(k+1)} = \frac{\sum_j \mu^{(k)}_j \boldsymbol{x}_j}{\sum_j \mu^{(k)}_j },
\end{equation}
where $\mu^{(k)}_j = w_j / \lVert \boldsymbol{x}^{(k)} - \boldsymbol{x}_j \rVert$. Starting with an initial guess, one keeps repeating this process until convergence. In practice, one needs to address the situation when $\boldsymbol{x}^{(k)}$ gets close to some $\boldsymbol{x}_j$, which
causes $\mu^{(k)}_j$ to blow up. In the Weiszfeld algorithm, one keeps track of the proximity of  $\boldsymbol{x}^{(k)}$ to all the $\boldsymbol{x}_j$, and  $\boldsymbol{x}^{(k+1)}$ is set to be $\boldsymbol{x}_i$ if $\lVert \boldsymbol{x}^{(k)} - \boldsymbol{x}_i \rVert < \varepsilon$ for some $i$. It has been proved
by many authors that the iterates converge globally, and even linearly (exponentially fast) under certain conditions, e.g., see discussion in \cite{XY1997}.

Following the recent ideas in \cite{CW2008,DDFG2009}, we have also tried regularizing \eqref{IRLS} by adding a small bias, namely by setting $\mu^{(k)}_j = w_j / (\lVert \boldsymbol{x}^{(k)} - \boldsymbol{x}_j \rVert^2 + \varepsilon_k^2)^{1/2}$.
The bias $\varepsilon_k$ is updated at each iterate, e.g., one starts with $\varepsilon_0 = 1$ and gradually shrinks it to zero.
The convergence properties of a particular flavor of this algorithm are discussed in \cite{DDFG2009}. We have tried both the Weiszfeld algorithm and the one in \cite{CW2008}. Experiments show us that faster
convergence is obtained using the latter, typically within $3$ to $4$ steps. The overall computational complexity of NLEM is $k^2\cdot S^2 \cdot N_{\mathrm{iter}}$ per pixel, where $N_{\mathrm{iter}}$ is the average number of
IRLS iterations. The complexity of the standard NLM is $k^2\cdot S^2$ per pixel.

\section{Experiments}
\label{S4}

We now present the result of some preliminary denoising experiments. For this, we have used the synthetic images shown in Fig. \ref{figIII} and some natural images. For NLEM, we computed the Euclidean median using the simple IRLS scheme described in \cite{CW2008}. For all experiments, we have set $S = 21, k = 7$, and $h = 10 \sigma$ for both algorithms.
These are the standard settings originally proposed in \cite{BCM2005}.

\begin{figure}[!htp]
  \centering
  \subfloat[\textit{Checker}.]{\includegraphics[width=0.5\linewidth]{ckb.eps}}
  \subfloat[\textit{Circles}.]{\includegraphics[width=0.5\linewidth]{circles.eps}}
  \caption{Synthetic grayscale test images, original size $256 \times 256$. Black is at intensity $0$, and white is at intensity $255$.}
  \label{figIII}
\end{figure}

First, we considered the \textit{Checker} image. We added noise with $\sigma =100$ resulting in a peak-signal-to-noise ratio (PSNR) of $8.18\text{dB}$. The PSNR improved to $17.94 \text{dB}$ after applying NLM, and with NLEM this
further went up to $19.45 \text{dB}$ (averaged over $10$ noise realizations).
This $1.5 \text{dB}$ improvement is perfectly explained by the arguments provided in Section \ref{S2}. Indeed, in Fig.\ref{fig6}, we have marked those pixels where the estimate from NLEM is significantly closer to the clean
image than that obtained using NLM. More precisely, denote the original image by $f_i$, the noisy image by $u_i$, and the denoised images from NLM and NLEM by $\hat{u}_i$ and $\hat{u}'_i$. Then the ``+'' in the figure denotes
pixels $i$ where $\lvert \hat{u}'_i - f_i \rvert < \lvert \hat{u}_i - f_i \rvert - 10$. Note that these points are concentrated around the edges where the median performs better than the mean.

\begin{figure}[!htp]
  \centering
   \subfloat[\textit{Checker} ($\sigma$=100).]{\label{fig6}\includegraphics[width=0.5\linewidth]{fig6.eps}}
  \subfloat[\textit{House} ($\sigma$=60).]{\label{fig7}\includegraphics[width=0.5\linewidth]{fig8.eps}}
  \caption{The ``+'' symbol marks pixels where the estimate returned by NLEM is significantly better than that returned by NLM.}
  \label{figIV}
\end{figure}

So what happens when we change the noise level? We went from $\sigma = 10$ to $\sigma = 100$ in steps of $10$. The plots of the corresponding PSNRs are shown in Fig. \ref{fig8}. At low noise levels ($\sigma < 30$), we see that
NLM performs as good or even better than NLEM. This is because at low noise levels the true neighbors in patch space are well identified, at least in the smooth regions. The difference between them is then mainly due to noise, and since the noise
is Gaussian, the least-squares approach in NLM gives statistically optimal results in the smooth regions. On the other hand, at low noise levels, the two clusters in Fig. \ref{figII} are well separated and hence the weights $w_{ij}$ for NLM are good enough to push
the mean towards the right cluster. The median and mean in this case are thus very close in the vicinity of edges. At higher noise levels, the situation completely reverses and NLEM performs consistently better than NLM. In Fig. \ref{fig8}, we see this phase
transition occurs at around $\sigma=30$. The improvement in PSNR is quite significant beyond this noise level, ranging between $0.5 \text{dB}$ and $2.1 \text{dB}$. The exact PSNRs are given in Table \ref{table1}.

\begin{figure}[!htp]
  \centering
   \subfloat[\label{fig8}]{\includegraphics[width=0.7\linewidth]{fig7.eps}}   \\
   \subfloat[\label{fig9}]{\includegraphics[width=0.7\linewidth]{ComparisonPSNR_FourMethods.eps}}
  \caption{Comparison of denoising performance for the \textit{Checker} image at noise $\sigma = 80$. NLM-kNN and NLEM-kNN refer to the
  respective modifications of NLM and NLEM where we only use the top $50\%$ of the weights.}
\label{figV}
\end{figure}

Next, we tried the above experiment on the \textit{Circles} image. The results are given in table \ref{table1}. The phase transition in thsi case occurs around $\sigma=25$.
NLM performs significantly better before this point, but beyond the phase transition, NLEM begins to perform better, and the gain in PSNR over NLM can be as large as $2.2\text{dB}$.
The method noise for NLM and NLEM obtained from a typical denoising experiment are shown in Fig. \ref{figVI}. Note that the method noise for NLEM appears more white (with less structural features) than that for NLM.

Finally, we considered some benchmark natural images, namely \textit{House}, \textit{Barbara}, and \textit{Lena}. The PSNRs obtained from NLM and NLEM for these images at different noise levels are shown in Table \ref{table1}.
The table also compares the Structural SIMilarity (SSIM) indices \cite{SSIM2004} obtained from the two methods. Note that a phase transition in SSIM is also observed for some of the images. In Fig. \ref{fig7}, we show the
pixels where NLEM does better (in the sense defined earlier) than NLM for the \textit{House} image at $\sigma = 60$. We again see that these pixels are concentrated around the edges.

\begin{figure}[!htp]
  \centering
\subfloat{\includegraphics[width=0.5\linewidth]{fig11.eps}}
{\includegraphics[width=0.5\linewidth]{fig13.eps}}
\caption{Comparison of the method noise for the \textit{Circles} image at $\sigma=80$.}
  \label{figVI}
\end{figure}

The improvements in PSNR and SSIM are quite dramatic for the synthetic images \textit{Checker} and \textit{Circles}. This is expected because they contain many edges and the edges have large transitions.
The improvement in PSNR and SSIM are less dramatic for the natural images. But, considering that NLM already provides top quality denoising results, this small improvement is already significant. We have also noticed that the performance of NLEM (and that of NLM) can be further improved by considering only the top $50\%$ of the $S^2$ neighbors with largest weights. This simple modification improves both NLM and
NLEM (and also their run time), while still maintaining their order in terms of PSNR performance. A comparison of the four different methods is given in Fig. \ref{fig9}.
The results provided here are rather incomplete, but they already provide a good indication of the denoising potential of NLEM at large noise levels.  While we have considered a fixed setting for the parameters, our general observation based on extensive experiments is that
NLEM consistently performs better than NLM beyond the phase transition, irrespective of the parameter settings. In future, we plan to investigate ways of further improving NLEM, and study the effect of the
parameters on its denoising performance.

\begin{sidewaystable}
\caption{Comparison of NLM and NLEM in terms of PSNR and SSIM, at noise levels $\sigma = 10, 20, \ldots, 100$ (Averaged over $10$ noise realizations). Same parameters used: $S=21, k = 7$, and $h = 10\sigma$.}
\centering
\begin{tabular}{l  c rrrrrrrrrr}

\hline

\bf{Image} & \bf{Method} &\multicolumn{10}{c}{\bf{PSNR (dB)}} \\

\hline
&NLM    &\bf{40.04}  &\bf{35.16}  &\bf{31.74}   &27.84     &25.21          &23.13       &21.39       &19.96       &18.84       &17.94 \\[-1ex]
\raisebox{1.5ex}{\textit{Checker}}
&NLEM   &39.73       &34.66       &31.66   &\bf{29.37} &\bf{26.71}     &\bf{24.76}  &\bf{23.22}  &\bf{21.82}  &\bf{20.50}  &\bf{19.45}  \\
\hline

&NLM    &\bf{37.31}  &\bf{34.67}       &31.79    &28.82  &26.46  & 24.69  &23.10  &21.58  &20.23  &19.03 \\[-1ex]
\raisebox{1.5ex}{\textit{Circles}}
&NLEM   &34.27      & 34.08  &\bf{32.33} &\bf{30.05}  &\bf{27.92}   &\bf{26.24}   &\bf{24.87}  &\bf{23.57} &\bf{22.36}  &\bf{21.16} \\

\hline

&NLM    &\bf{34.22}  &29.78        &26.88       &25.21        &24.07      &23.37      &22.78        &22.38      &22.06  &21.81 \\[-1ex]
\raisebox{1.5ex}{\textit{House}}
&NLEM   &33.96       &\bf{30.10}  &\bf{27.15}  &\bf{25.39}  &\bf{24.30}  &\bf{23.51}  &\bf{22.96}  &\bf{22.54}  &\bf{22.17}  &\bf{21.95} \\

\hline

&NLM    &\bf{32.37}   &27.39        &24.93   &23.52       &22.64        &22.04      &21.62       &21.29      &21.07     &20.88 \\[-1ex]
\raisebox{1.5ex}{\textit{Barbara}}
&NLEM   &32.11       &\bf{27.75} &\bf{25.26} &\bf{23.84}  &\bf{22.90}  &\bf{22.29}  &\bf{21.83}  &\bf{21.48} &\bf{21.20}  &\bf{21.01} \\

\hline

&NLM      &\bf{33.24}   &29.31   &27.40  &26.16  &25.24  &24.54  &24.04  &23.66  &23.34  &23.06 \\[-1ex]
\raisebox{1.5ex}{\textit{Lena}}
&NLEM     &33.15  &\bf{29.45}  &\bf{27.61}  &\bf{26.40}  &\bf{25.53}  &\bf{24.84}  & \bf{24.31}  &\bf{23.90} &\bf{23.53}  &\bf{23.24} \\
\end{tabular}

\begin{tabular}{l  c rrrrrrrrrr}
\hline
\hline

&  &\multicolumn{10}{c}{\bf{SSIM (\%)}} \\

\hline
&NLM    &\bf{99.41}      &\bf{98.66}  &97.59        &95.51         &92.32           &88.35       &83.37       &78.05       &72.31       &66.47 \\[-1ex]
\raisebox{1.5ex}{\textit{Checker}}
&NLEM   &99.36           &98.51       &\bf{97.60}   &\bf{96.30}    &\bf{ 94.03}     &\bf{91.18}  &\bf{87.81}  &\bf{84.04}  &\bf{79.59}  &\bf{74.51}  \\
\hline

&NLM    &\bf{96.31}  &\bf{94.02}       &\bf{91.43}    &88.66        &85.79        &82.49       &78.28       &73.28       & 68.13  &63.54 \\[-1ex]
\raisebox{1.5ex}{\textit{Circles}}
&NLEM   &95.28      &93.60             &91.20     &\bf{88.74}  &\bf{86.30}   &\bf{83.75}   &\bf{80.90}  &\bf{77.85}  &\bf{ 74.21}  &\bf{70.21} \\

\hline

&NLM    &\bf{86.90}        &\bf{81.80}       &76.93        &72.97       & 69.57      &66.86       &64.38        &62.24       &60.33       &58.55 \\[-1ex]
\raisebox{1.5ex}{\textit{House}}
&NLEM   &86.86             &81.25            &\bf{77.61}   &\bf{73.85}  &\bf{70.51}  &\bf{67.77}  &\bf{65.21}   &\bf{62.96}  &\bf{60.97}  &\bf{ 59.03} \\

\hline

&NLM    &89.95             &79.48        &71.32       &65.20       &60.79        &57.42      &54.64        &52.58      &50.78     &49.27 \\[-1ex]
\raisebox{1.5ex}{\textit{Barbara}}
&NLEM   &\bf{90.25}       &\bf{80.38}    &\bf{72.49} &\bf{66.48}  &\bf{62.04}    &\bf{58.57}  &\bf{55.64}  &\bf{53.45} &\bf{51.50}  &\bf{49.86} \\

\hline

&NLM      &86.95        &80.41      &76.45       &73.42        &70.80      &68.45        &66.46         &64.66    &62.94            &61.32 \\[-1ex]
\raisebox{1.5ex}{\textit{Lena}}
&NLEM     &\bf{87.06}  &\bf{80.57}  &\bf{76.77}  &\bf{73.97}  &\bf{71.50}  &\bf{69.21}  & \bf{67.20}  &\bf{65.32}  &\bf{ 63.52}  &\bf{61.81} \\

\hline
\end{tabular}
\label{table1}
\end{sidewaystable}

\section{Discussion}
\label{S5}

The purpose of this note was to communicate the idea that one can improve (often substantially) the denoising results of NLM by replacing the $\ell^2$ regression on patch space by the more robust $\ell^1$ regression. This led us to propose the NLEM algorithm. The experimental results presented in this paper reveal two facts: (a) Phase transition phenomena -- NLEM starts to perform better (in terms of PSNR)
beyond a certain noise level, and (b) The bulk of the improvement comes from pixels close to sharp edges. The latter fact indicates that NLEM is better suited for denoising images
that have many sharp edges. This suggests that we could get similar PSNR improvements if we simply used NLEM in the vicinity of edges and the less expensive NLM elsewhere. Unfortunately, it is difficult to get a reliable
edge map from the noisy image. On the other hand, observation (b) suggests that, by comparing the denoising results of NLM and NLEM, one can devise a robust method of detecting edges at large noise levels. We note that Wang et al. have recently proposed a non-local median-based estimator in \cite{WSL}. This can be considered as a special case of NLEM corresponding to $k=1$, where single pixels (instead of patches) are used
for the median computation. On the other hand, some authors have proposed median-based estimators for NLM where the noisy image is median filtered before computing the weights \cite{CFFM2010,WSL}.
In fact, most of the recent innovations in NLM were concerned with better ways of computing the weights; e.g., see \cite{T2009,deVK2011}, and reference therein.
It would be interesting to see if the idea of robust Euclidean medians could be used to complement these innovations.


\section{Introduction}
The interference channel (IC) describes a network where multiple transmitters communicate with their intended receivers via a common medium. The characterization of the capacity region for a two-user IC is an open problem except for the strong and very strong interference cases \cite{Carleial:75IT,Sato:81IT,Han&Kobayashi:81IT}. To date, the largest achievable rate region is the celebrated Han-Kobayashi (HK) region that employs rate splitting at the transmitters and simultaneous decoding at the receivers \cite{Han&Kobayashi:81IT}. Not surprisingly, for those ICs whose capacity regions are completely characterized, it is without an exception that the capacity region coincides with the HK region.

However, the general HK region involves a time sharing variable that makes its evaluation intractable. For the Gaussian interference channel (GIC), another difficulty is the input distribution. A two-user GIC in its standard form can be represented as
\begin{equation}
\begin{array} {ccc}
Y_1&=&X_1+bX_2+Z_1, \\
Y_2&=&aX_1+X_2+Z_2,
\end{array}
\label{eq:GIC}
\end{equation}
where $X_1$ and $X_2$ are the input signals and are subject to respective power constraints $P_1$ and $P_2$; $Y_1$ and $Y_2$ are the received signals; $Z_1$ and $Z_2$ are Gaussian noises of unit variance and are independent of the inputs $X_1$ and $X_2$. This model is depicted in Fig.~\ref{Fig:f1}. While for all the cases where the capacity results are known for a GIC, the optimal input distribution is invariably Gaussian, it is not yet known (or proven) that such is the case for the general GIC.

There has been recent progress in obtaining computable subregion of the HK achievable region using Sato's non-na\"{i}ve frequency division \cite{Shang&Chen:07ISIT}. For the one-sided GIC (denoted as ZGIC) shown in Fig.~\ref{Fig:5}(a),  Motahari and Khandani established that such a non-na\"{i}ve frequency division scheme achieves the HK region with Gaussian input \cite{Motahari&Khandani:09IT}. Most recently, Costa introduced the so-called {\em noiseberg} scheme which uses water filling to achieve optimal power sharing between two orthogonal dimensions \cite{Costa:Noiseberg}. It turns out, as shown in the next section, that this simple noiseberg scheme achieves precisely the same HK region with Gaussian input.

\begin{figure}
\centering
\includegraphics{Fig1.eps}
\caption{ Two-user GIC}\label{Fig:f1}
\end{figure}

This paper focuses on GICs with mixed interference (MGIC) and with $ab\geq 1$, $a\leq 1$ and $b\geq 1$ (cf. Fig.~\ref{Fig:f1} and Eq.~(\ref{eq:GIC})). We describe a simple coding scheme that combines the noiseberg scheme with that of simultaneous decoding at the receiver experiencing strong interference. The obtained rate region is subsequently shown to coincide with the HK region with Gaussian input.

The rest of the paper is organized as follows. In Section \ref{sec:ZGIC}, we review the noiseberg scheme for the ZGIC and provide a proof of its equivalence to the HK region with Gaussian input. Section \ref{MGIC} describes the coding scheme for a class of MGIC and proves that the scheme achieves the HK region with Gaussian input. Section \ref{Conclusion} concludes this paper.

\section{Noiseberg Region for the ZGIC}\label{sec:ZGIC}

We consider the degraded GIC shown in Fig.~\ref{Fig:5}(b), which is equivalent to the ZGIC with $a<1$ in Fig.~\ref{Fig:5}(a)\cite{Costa:85IT}.

The noiseberg region, denoted by $\mathcal{R_{N}}$ and introduced by Costa in \cite{Costa:Noiseberg} for a ZGIC with weak interference ($a<1$ in Fig.~\ref{Fig:5}(a)) is the set of all nonnegative rate pairs $(R_1, R_2)$ satisfying
\@setsize\small{9pt}\viiipt\@viiipt\let\@listi\@listI
\begin{eqnarray}
R_{1} &\leq& \bar{\lambda}R_{1 \bar{\lambda}} + \lambda R_{1 \lambda}, \nonumber \\
R_{2}  &\leq& \bar{\lambda}R_{2 \bar{\lambda}}, \nonumber
\end{eqnarray}
\normalsize
where
\@setsize\small{9pt}\viiipt\@viiipt\let\@listi\@listI
\begin{eqnarray}
R_{1 \bar{\lambda}} &\leq& \frac{1}{2} \log \left(1+ \frac{P_{1A}}{\bar{\lambda}}\right) +\frac{1}{2} \log \left(1+ \frac{a^{2} \frac{P_{1C}}{\bar{\lambda}}}{1+a^{2}\frac{P_{1A}}{\bar{\lambda}}+ \frac{P_{2}}{\bar{\lambda}}}\right), \nonumber \\\label{eq:2}\\
R_{2 \bar{\lambda}}&\leq& \frac{1}{2} \log \left(1+ \frac{\frac{P_{2}}{\bar{\lambda}}}{1+a^{2}\frac{P_{1A}}{\bar{\lambda}}}\right), \label{eq:3}\\
R_{1 \lambda} &\leq& \frac{1}{2} \log \left(1+ \frac{P_{1B}}{\lambda}\right) \label{eq:1},
\end{eqnarray}
\normalsize
and the power limits $P_{1A}$, $P_{1B}$ and $P_{1C}$ are determined by two parameters $h$ and $\lambda$ such that
\@setsize\small{9pt}\viiipt\@viiipt\let\@listi\@listI
\begin{eqnarray}
\frac{P_{1A}}{\bar{\lambda}} &=& P_{1}-\frac{P_{2}\lambda}{a^{2}\bar{\lambda}}- \lambda \min \left\{h, \frac{1-a^2}{a^2}\right\} \nonumber \\ &&-\max\left \{0,h - \frac{1-a^2}{a^2}\right\}, \nonumber \\
\frac{P_{1B}}{\lambda} &=& P_{1}+\frac{P_{2}}{a^{2}}+ \bar{\lambda} \min\left\{h, \frac{1-a^{2}}{a^{2}}\right\},\nonumber \\
\frac{P_{1C}}{\bar{\lambda}} &=& \max\left\{0, h-\frac{1-a^{2}}{a^{2}}\right\}. \nonumber
\end{eqnarray}
\normalsize

\begin{figure}[!t]
\centering
\subfloat[ZGIC]{\includegraphics{FigZ.eps}}
\subfloat[Degraded GIC ]{\includegraphics{Fig5.eps}}
\caption{ZGIC}
\label{Fig:5}
\end{figure}

Costa \cite{Costa:Noiseberg} showed that $\mathcal{R_{N}}$ is achievable for the ZGIC with weak interference by a coding scheme that uses a two-band non-na\"{i}ve frequency division multiplexing (FDM) with water filling for optimal power allocation between the two subbands. The coding scheme, as well as its achievable region, involves two parameters $0\leq \lambda\leq 1$ and $h\geq 0$. They vary over the admissible region as shown in Fig. \ref{Fig:f2}, resulting in different transmission schemes depending on the values of the parameters. The parameter $\lambda$ determines how to allocate the frequency band.

\begin{itemize}
\item
The $\lambda$ subband is reserved for the communication between transmitter 1 and receiver 1.
\item
The $\bar{\lambda}$ subband is shared by both transceiver pairs and the corresponding coding scheme is determined by the other parameter $h$.
\end{itemize}

As the noise $Z_2$ does not affect the transmission of $X_1$, water filling allows the overall power level in the $\lambda$-subband to be raised above that of the $\bar{\lambda}$-subband, with part of the noise spectrum of $Z_2$ floating above the signal level. This phenomenon, i.e., the existence of difference in heights of power spectrum for the two subbands is referred to as the \emph{noiseberg}. The parameter h is defined as the height of total power density in the $\lambda$-subband above that of $X_2$'s power density in the $\bar{\lambda}$-subband. Different $h$ values divide the admissible region for the parameter pairs into two regions, each employing a different coding scheme in the $\bar{\lambda}$-subband:

\noindent {\bf Multiplex region} This corresponds to $h \leq \frac{1-a^2}{a^2}$. As shown in Fig.~\ref{Fig:f3}, $Z_2'$ prevents user 1's power from spilling over to the $\bar{\lambda}$-band thus no rate-splitting is involved. Receiver 2 decodes $W_1$ first, subtracts it and decodes $W_2$.

\noindent  {\bf Overflow region} This corresponds to $h > \frac{1-a^2}{a^2}$. As shown in Fig.~\ref{Fig:f4}, water-filling of user 1's power occurs as the power spills over from the $\lambda$-subband to the $\bar{\lambda}$-subband. The encoding scheme in the $\bar{\lambda}$ subband thus involves rate splitting for $W_1$: a common message $W_{1c}$ with power $P_{1c}$ decoded by both receivers and a private message $W_{1p}$ with power $P_{1A}$ decoded only by receiver 1. Receiver 2 decodes $W_{1c}$ first, subtracts it, and decodes $W_2$, all the while treating $W_{1p}$ as noise.

\begin{figure}
\centering
\includegraphics{Fig2.eps}
\caption{Admissible region for $(\lambda, h)$}\label{Fig:f2}
\end{figure}
\begin{figure}
\centering
\includegraphics{Fig3.eps}
\caption{Multiplex region}\label{Fig:f3}
\end{figure}
\begin{figure}
\centering
\includegraphics{Fig4.eps}
\caption{Overflow region}\label{Fig:f4}
\end{figure}

It is remarkable that such a simple transmission scheme turns out to achieve precisely the same HK region with time sharing and with Gaussian input.
\begin{theorem} \label{lemma:Noiseberg}
For the weak ZGIC, $\mathcal{R_{N}} = \mathcal{R_{HK}}$.
\end{theorem}
\begin{proof}
Motahari and Khandani showed that for the ZGIC, the non-na\"{i}ve FDM region, denoted by $\mathcal{R_{FDM}}$, is equivalent to $\mathcal{R_{HK}}$, whose boundary points can be characterized by the optimization problem \cite[Eq. (151)]{Motahari&Khandani:09IT}. It suffices to verify the equivalence between $\mathcal{R_{N}}$ and $\mathcal{R_{FDM}}$.

We start by considering water filling in the two-band FDM applied to the degraded GIC shown in Fig.~\ref{Fig:5}(b). First, we split $W_1$ into private message $W_{1p}$ with power constraint $P_{1p}$ and common message $W_{1c}$ with power constraint $P_{1c}$ such that $P_{1p}+P_{1c}=P_1$. Power allocation into $\lambda$ and $\bar{\lambda}$ subbands is done in the following order. First, $P_{1p}$ is allocated to the two subbands in an arbitrary way. On top of that, $P_2$ is allocated to the two subbands via waterfilling. As $Y_2'$ sees additional noise $Z_2'$, $P_2$ is allocated on top of $Z_2'$ (see, e.g., Fig.~\ref{Fig:6}(d)). Finally, $P_{1c}$ is allocated to the two subbands, again, using waterfilling.
\begin{figure*}[!t]
\centering
\subfloat[]{\label{Fig:6a}\includegraphics{Fig6a.eps}}
\subfloat[]{\label{Fig:6b}\includegraphics{Fig6b.eps}}
\\\subfloat[]{\label{Fig:6c}\includegraphics{Fig6c.eps}}
\subfloat[]{\label{Fig:6d}\includegraphics{Fig6d.eps}}
\caption{Possible power allocation outcomes of the two-band non-na\"{i}ve FDM scheme with ZGIC}
\label{Fig:6}
\end{figure*}

Depending on $P_{1p}$ and its allocation between the two subbands, there are four possible power allocation outcomes of this scheme, as shown in Fig.~\ref{Fig:6}. Since the scenarios illustrated in Fig.~\ref{Fig:6}(a)(b)(c) are equivalent to noiseberg cases, it remains to argue that the power allocation scheme with flat spectrum top as in Fig.~\ref{Fig:6}(d) is not optimal. This is because the achievable rates under such a scheme are formed by convex combinations of points on the curve of associated broadcast channel capacity, as the flat top requires $\frac{P_{1c \bar{\lambda}}}{\bar{\lambda}} = \frac{P_{1c \lambda}}{\lambda}$. As the broadcast channel capacity curve is convex, we can only achieve the points on the chord, which are dominated by the points on the curve corresponding to the scheme with no frequency division. Thus they are not optimal.

Next we generalize this conclusion to three-band FDM. Similarly we argue that the power-bandwidth allocation schemes with two or more subbands sharing the same flat top are not optimal. If two subbands are occupied solely by the transmission of $X_{1}$, they can be merged into one and this reduces to one of the noiseberg scenarios.
\end{proof}

\section{Achievable Region for a Class of MGIC} \label{MGIC}
For the MGICs with $0<a \leq 1$, $b \geq 1$ and $ab \geq 1$, the HK region with Gaussian input can be simplified to be the set of all rate pairs $(R_1, R_2)$ satisfying
\@setsize\small{9pt}\viiipt\@viiipt\let\@listi\@listI
\begin{eqnarray}
R_1 &\leq& \sum\limits_{i=1}^{q} \lambda_i \left\{\frac{1}{2} \log (1 + P_{1i})\right\},\nonumber \\
R_2 &\leq&  \sum\limits_{i=1}^{q} \lambda_i \left\{\frac{1}{2} \log \left(1 + \frac{P_{2i}}{1+a^2 \alpha_i P_{1i}}\right)\right\}, \nonumber \\
R_1 + R_2 &\leq&  \sum\limits_{i=1}^{q} \lambda_i \left\{\frac{1}{2} \log \left(1+ \frac{P_{2i}+a^2\bar{\alpha}_i P_{1i}}{1+a^2\alpha_i P_{1i}}\right) \right. \nonumber \\
&&\left. + \frac{1}{2} \log (1+\alpha_i P_{1i})\right\},\nonumber
\end{eqnarray}
\normalsize
where $q \in \mathbb{N}$, $i \in \{1, \dots, q\}$, $0 \leq \alpha_i \leq 1$, , $\bar{\alpha}_i = 1-\alpha_i$, $\lambda_i \geq 0$, $\sum\limits_{i=1}^{q} \lambda_i = 1$, $\sum\limits_{i=1}^{q} \lambda_i P_{1i} \leq P_1$ and $\sum\limits_{i=1}^{q} \lambda_i P_{2i} \leq P_2$. We refer to $\mathcal{R_{HK}}$ as the above HK region with Gaussian input. The difficulty in using this region is largely due to the presence of the time sharing variable.

We now describe a simple transmission scheme for a MGIC with $ab\geq 1$. The scheme resembles the noiseberg scheme as it also utilizes the two-band non-na\"{i}ve FDM. Specifically,
in the $\lambda$-subband, only transmitter 1 transmits and receiver 1 decodes $W_1$ with a rate constraint defined in (\ref{eq:1}).

In the $\bar{\lambda}$-subband, transmitter 2 does not use rate splitting. Transmitter 1, on the other hand, employs two encoding schemes depending on the value of $h$. The corresponding decoding schemes are also different. We describe them in details below.

\noindent{\bf Multiplex region} This corresponds to $h \leq \frac{1-a^2}{a^2}$. Sequential decoding is used at both receivers. Receiver 1 first decodes $W_2$, which leads to the constraint
\@setsize\small{9pt}\viiipt\@viiipt\let\@listi\@listI
\begin{eqnarray}
R_{2\bar{\lambda}} & \le & \frac{1}{2} \log \left(1+ \frac{b^{2}P_{2}/\bar{\lambda}}{1+ P_{1A}/\bar{\lambda}}\right). \label{eq:4}
\end{eqnarray}
\normalsize
Subsequently, $X_2$ is subtracted from $Y_1$ and $W_1$ is decoded with constraint (\ref{eq:2}). Receiver 2 decodes $W_2$ with constraint (\ref{eq:3}). Comparing (\ref{eq:3}) and (\ref{eq:4}), we see that (\ref{eq:4}) is redundant. Therefore sequential decoding in the multiplex region achieves $\mathcal{R_{N}}$.

\noindent{\bf Overflow region} This corresponds to $h > \frac{1-a^2}{a^2}$. Receiver 1 employs simultaneous decoding while receiver 2 still uses sequential decoding. For receiver 1, let $S_{1}$, $T_{1}$ and $T_{2}$ be, respectively, the rates of transmitter 1's private message $W_{1p}$, common message $W_{1c}$ and $W_2$. Then $R_{1\bar{\lambda}} = S_{1}+T_{1}$ and $R_{2\bar{\lambda}} = T_{2}$. Evaluation of error probability will give us
\@setsize\small{9pt}\viiipt\@viiipt\let\@listi\@listI
\begin{eqnarray}
S_{1} &\leq& \frac{1}{2}\log \left(1+\frac{P_{1A}}{\bar{\lambda}}\right), \label{eq:5}\\
T_{1} &\leq& \frac{1}{2}\log \left(1+\frac{P_{1C}}{\bar{\lambda}}\right),\\
T_{2} &\leq& \frac{1}{2}\log \left(1+\frac{b^{2}P_2}{\bar{\lambda}}\right),\\
S_{1}+T_{1} &\leq& \frac{1}{2}\log \left(1+\frac{P_{1A}}{\bar{\lambda}}+\frac{P_{1C}}{\bar{\lambda}}\right),\\
S_{1}+T_{2} &\leq& \frac{1}{2}\log \left(1+\frac{P_{1A}}{\bar{\lambda}}+\frac{b^{2}P_2}{\bar{\lambda}}\right),\\
T_{1}+T_{2} &\leq& \frac{1}{2}\log \left(1+\frac{P_{1C}}{\bar{\lambda}}+\frac{b^{2}P_2}{\bar{\lambda}}\right),\\
S_{1}+T_{1}+T_{2} &\leq& \frac{1}{2}\log \left(1+\frac{P_{1A}}{\bar{\lambda}}+\frac{P_{1C}}{\bar{\lambda}}+\frac{b^{2}P_2}{\bar{\lambda}}\right).
\end{eqnarray}
\normalsize

As receiver 2 decodes $W_{1C}$ and $W_{2}$ sequentially, there will be two constraints
\@setsize\small{9pt}\viiipt\@viiipt\let\@listi\@listI
\begin{eqnarray}
T_{1} &\leq& \frac{1}{2} \log \left(1+ \frac{a^{2} \frac{P_{1C}}{\bar{\lambda}}}{1+a^{2}\frac{P_{1A}}{\bar{\lambda}}+ \frac{P_2}{\bar{\lambda}}}\right), \\
T_{2} &\leq& \frac{1}{2} \log \left(1+ \frac{ \frac{P_2}{\bar{\lambda}}}{1+a^{2}\frac{P_{1A}}{\bar{\lambda}}}\right). \label{eq:6}
\end{eqnarray}
\normalsize

Fourier-Motzkin elimination on (\ref{eq:5})-(\ref{eq:6}) gives us
\@setsize\small{9pt}\viiipt\@viiipt\let\@listi\@listI
\begin{eqnarray}
R_{1\bar{\lambda}} &\leq& \frac{1}{2}\log \left(1+\frac{P_{1A}}{\bar{\lambda}}+\frac{P_{1C}}{\bar{\lambda}}\right), \label{eq:7}\\
R_{1\bar{\lambda}} &\leq& \frac{1}{2}\log \left(1+\frac{P_{1A}}{\bar{\lambda}}\right) \nonumber \\ &&+  \frac{1}{2} \log \left(1+ \frac{a^{2} \frac{P_{1C}}{\bar{\lambda}}}{1+a^{2}\frac{P_{1A}}{\bar{\lambda}}+ \frac{P_2}{\bar{\lambda}}}\right), \label{eq:9}\\
R_{2\bar{\lambda}} &\leq& \frac{1}{2} \log \left(1+ \frac{ \frac{P_2}{\bar{\lambda}}}{1+a^{2}\frac{P_{1A}}{\bar{\lambda}}}\right), \label{eq:10}\\
R_{1\bar{\lambda}}+R_{2\bar{\lambda}} &\leq& \frac{1}{2}\log \left(1+\frac{P_{1A}}{\bar{\lambda}}+\frac{P_{1C}}{\bar{\lambda}}+\frac{b^{2}P_2}{\bar{\lambda}}\right), \label{eq:11}\\
R_{1\bar{\lambda}}+R_{2\bar{\lambda}} &\leq& \frac{1}{2}\log \left(1+\frac{P_{1A}}{\bar{\lambda}}+\frac{b^{2}P_2}{\bar{\lambda}}\right) \nonumber \\ &&+  \frac{1}{2} \log \left(1+ \frac{a^{2} \frac{P_{1C}}{\bar{\lambda}}}{1+a^{2}\frac{P_{1A}}{\bar{\lambda}}+ \frac{P_2}{\bar{\lambda}}}\right), \label{eq:12}\\
R_{1\bar{\lambda}} + 2R_{2\bar{\lambda}} &\leq& \frac{1}{2}\log \left(1+\frac{P_{1A}}{\bar{\lambda}}+\frac{b^{2}P_2}{\bar{\lambda}}\right) \nonumber \\&&+ \frac{1}{2}\log \left(1+\frac{P_{1C}}{\bar{\lambda}}+\frac{b^{2}P_2}{\bar{\lambda}}\right). \label{eq:8}
\end{eqnarray}
\normalsize

Then the achievable rate region in overflow region is
\@setsize\small{9pt}\viiipt\@viiipt\let\@listi\@listI
\begin{eqnarray}
R_{1} &\leq& \bar{\lambda}R_{1 \bar{\lambda}} + \lambda R_{1 \lambda}, \nonumber \\
R_{2}  &\leq& \bar{\lambda}R_{2 \bar{\lambda}}, \nonumber
\end{eqnarray}
\normalsize
where $R_{1 \bar{\lambda}}$, $R_{2 \bar{\lambda}}$ and $R_{1 \lambda}$ satisfy (\ref{eq:7})-(\ref{eq:8}) and (\ref{eq:1}). In the appendix we simplify this region and prove that it is equivalent to $\mathcal{R_{N}}$.

Combining the results in both the multiplex and overflow regions, we conclude that $\mathcal{R_{N}}$ is achievable for the MGICs with $ab>1$. Thus we have,

\begin{theorem} \label{theorem: Mixed IC}
For the MGICs with $ab>1$, $\mathcal{R_{N}}$ associated with the ZGIC obtained by removing the interfering link with gain $b$ is achievable.
\end{theorem}

\begin{corollary}
For the MGICs with $ab>1$, $\mathcal{R_{N}} = \mathcal{R_{HK}}$.
\end{corollary}
\begin{proof}
On the one hand, with Gaussian input, $\mathcal{R_{HK}}$ of the MGIC is a subset of that of the ZGIC. On the other hand, $\mathcal{R_{N}} \subseteq \mathcal{R_{HK}}$ for the MGIC in general. But Theorem 1 states that $\mathcal{R_{N}} = \mathcal{R_{HK}}$ for the ZGIC. This establishes the equivalence between $\mathcal{R_{N}}$ and $\mathcal{R_{HK}}$ for the MGICs with $ab>1$.
\end{proof}

This is consistent with \cite[Theorem 12]{Motahari&Khandani:09IT} that states for the MGICs with $ab>1$, $\mathcal{R_{HK}}$ is equivalent to that of the corresponding ZGIC.
\section{Conclusion} \label{Conclusion}
This paper established the equivalence between Costa's noiseberg region and the HK region with Gaussian input for the ZGIC. For the MGICs with $ab>1$, an achievable rate region was developed by rate splitting and a mixture of sequential and simultaneous decoding. By comparing the new region to the noiseberg region of the ZGIC, we proved that it is a simplification of the HK region for the MGICs with $ab>1$ and with Gaussian codebook.

\section*{Acknowledgment}
The authors would like to thank Prof. Max Costa for helpful discussions about the noiseberg region.


\section{Definitions, Notation, Statements}
We write $\C$ for the field of complex numbers. For every positive
integer $m$ let us  consider the affine space $\A^m=\C^m$ of all
monic complex polynomials of degree  $m$
$$u(x)=x^m+\sum_{i=0}^{m-1} a_i x^i$$ with coefficients
 $a=(a_0, \dots , a_{m-1})\in \C^m=\A^m$.
It is convenient to identify the tangent
  space
   $\C^m$  to $u(x)\in \A^m$ with the space of all polynomials
  $p(x)$ of degree $\le m-1$. Namely, to a polynomial $p(x)=\sum_{i=0}^{m-1} c_i x^i$
 one assigns the tangent vector $(c_0, \dots, c_{m-1}) \in \C^m$
 that corresponds to
``the tangency class  at $u(x)$ of the curve"
 $\epsilon \to u(x)+\epsilon \cdot p(x) \in \A^m$ \cite[Part II, Ch. III, Sect. 8, pp. 81--82]{Serre}.

Let $P_m\subset\A^m$ be the everywhere dense Zariski-open affine
subset
  that consists of all  polynomials  without multiple roots. Let $f(x)=x^m+\sum_{i=0}^{m-1} a_i x^i\in P_m$
and let us choose a root $\alpha$   of $f(x)$. Locally (with respect
to $a$), one may view
  $\alpha$  (using Implicit Function Theorem) as a holomorphic
(univalued) function in  $a=(a_0, \dots , a_{m-1})$.  We have
(\cite[Sect. 2]{ZarhinMatZametki2012})
$$d\alpha/da_i=-  [f^{\prime}(\alpha)]^{-1} \alpha^i.$$
(Since $\alpha$ is a simple root of $f(x)$, we have $f^{\prime}(\alpha)\ne 0$.)
We also have (ibid)
 $$d f^{\prime}(\alpha)/d a_i =i \alpha^{i-1} - [f^{\prime}(\alpha)]^{-1} \alpha^i  f^{\prime\prime}(\alpha)$$
 (of course, if $i=0$ then the first term disappears). Using these formulas, let us compute the differential
 $d\N: \C^m \to \C$ (at $f(x)$) of   locally defined holomorphic function
$$\N: P_{m} \to  \C, \ f(x)  \mapsto f^{\prime}(\alpha).$$

 It follows that $d\N$ sends the tangent vector $p(x)=\sum_{i=0}^{m-1} c_i x^i$ to the number
$$d\N(p(x))=\sum_{i=0}^{m-1} c_i \frac{d f^{\prime}}{d a_i}(\alpha)=p^{\prime}(\alpha)- [f^{\prime}(\alpha)]^{-1}p(\alpha)  f^{\prime\prime}(\alpha).$$

\begin{ex}
\label{exm}
Suppose that $m \ge 3$ and $f(x)=x^m-x$. Then $\alpha$ is either zero or $(m-1)$th  root of unity.
If $\alpha=0$ then $f^{\prime\prime}(0)=0$ and
$$d\N(p(x))=p^{\prime}(0)=c_1.$$
The {\sl gradient} of $\N$ at $f(x)=x^m-x$ (with respect to the root
$0$) is
$$Q_1(0)=(0,1, \dots ,0) \in \C^m.$$
If $\alpha^{m-1}=1$ then
$$f^{\prime}(\alpha)=m\alpha^{m-1}-1=m-1,$$
$$f^{\prime\prime}(\alpha)=m(m-1)\alpha^{m-2}=m(m-1)/\alpha,$$
and
$$d\N(p(x))=p^{\prime}(\alpha)-\frac{m p(\alpha)}{\alpha}.$$
The {\sl gradient}  of $\N$ at $f(x)=x^m-x$ (with respect to the
root $\alpha$) is
$$Q_1(\alpha)=\left(-\frac{m}{\alpha}, (1-m), (2-m)\alpha, \dots ,-\alpha^{m-2}\right) \in \C^n.$$
\end{ex}

Let $n \ge 2$ be an integer and $g(x) \in \C[x]$ a degree $n$ monic
polynomial with complex coefficients. For every positive integer $r$
we denote by $g^{\circ r}(x)$ the composition $g(\dots g(x))$ ($r$
times). Clearly, $g^{\circ r}(x)$ is a degree $n^r$ monic polynomial
with complex coefficients.
Let us consider the polynomial map
$$G: \C \to\C, \ z \mapsto g(z).$$
Clearly, the fixed points of $G$ are exactly the roots of $g(x)-x$
while the roots of  $g^{\circ r}(x)-x$ are exactly the points of
period (dividing) $r$.

\begin{ex}
If $g(x)=x^n$ then $g^{\circ r}(x)=x^{n^r}, \ g^{\circ r}(x)-x=
x^{n^r}-x$.
\end{ex}

 We write $Z_{n,r}\subset
\A^n$ for the everywhere dense Zariski-open  affine subset that
consists of all monic degree $n$ polynomials $g(x)$ such that
$g^{\circ r}(x)-x$ lies in $P_{n^r}$ (i.e., does not have multiple
roots). For example, $x^n \in Z_{n,r}$ for all $r$. Clearly, for
every positive integer $m$
$$Z_{m,1}= \{f(x)+x\mid f(x)\in P_m \}.$$
It is also clear that the holomorphic map
$$U_m: Z_{m,1} \to P_{m}, \ h(x)\mapsto h(x)-x$$
is a holomorphic isomorphism, whose tangent map
$$dU_n:\C^m \to \C^m$$
is the identity map at all points of $Z_{m,1}$.

 Let us consider a locally defined holomorphic function
$$\M^{r}: Z_{n,r}\to Z_{n^r,1} \to P_{n^r} \to  \C, \ g(x) \mapsto g^{\circ r}(x)
\overset{U_{n^r}}{\longmapsto}
  g^{\circ r}(x)-x \overset{\N}{\mapsto} [g^{\circ
r}(x)-x]^{\prime}(\alpha)$$ where $\alpha$ is a root of $g^{\circ
r}(x)-x$. We are going to discuss its differential $$d\M^{r}_{\mid
g(x)}:\C^n \to \C,$$ paying special attention to the computation of
the corresponding gradient
$$\grad(\M^{r})_{\mid g(x)}\in \C^n$$
 at the point $g(x)=x^n\in Z_{n,r}$. In what follows we denote this gradient by $Q_r(\alpha)$.
This notation is compatible with our previous notation for $Q_1(\alpha)$ in Example \ref{exm}.

\begin{rem}
Let us consider the locally defined multiplier function
$$\Mult^{r}: Z_{n,r}\to Z_{n^r,1} \to  \C, \ g(x) \longmapsto [g^{\circ r}]^{\prime}(\alpha).$$
Clearly, $\M^{r}(g)=\Mult^{r}(g)- 1$. It follows that the differentials
 $d\M^{r}$ and $d\Mult^{r}$ everywhere coincide. In other words
$$\grad(\M^{r} )_{\mid g(x)}=\grad(\Mult^{r})_{\mid g(x)} \ \forall g(x) \in Z_{n,r}.$$
\end{rem}

\begin{sect}
In order to state our main results, first notice that if $g(x) \in
\in Z_{n,r}$ and $\nu(n,r)$ is the number of of all orbits of length
$r$ for the map $z \mapsto g(z)$ then
$$\nu(n,r) \ge \frac{n}{r}$$
(see Subsection \ref{numberoforbits}). Second, let us consider a
positive integer $\ell$ and a sequence $\{r_1, \dots, r_{\ell}\}$ of
$\ell$ positive integers. Let $Z(n,\ell; r_1, \dots , r_{\ell})$ be
the intersection of all $Z_{n,r_i}$; it is a nonempty Zariski-open
affine subset in $\A^n$ that contains $g(x)=x^n$. Let $g(x)\in
Z(n,\ell; r_1, \dots , r_{\ell})$.
 For each $i$ pick a complex number $\beta_i$ that is a periodic point of $G: z\mapsto g(z)$
of {\bf exact period} $r_i$. Locally (with respect to  $g$), each
$\beta_i$ is a holomorphic function (in the coefficients of $g(x)$.)

Suppose that $\beta_1,  \dots \beta_{\ell}$ belong to {\bf distinct orbits} of $z\mapsto g(z)$. Let us consider the following $\ell$
locally defined holomorphic functions
$$\Mult_{\beta_i,r_i}: Z(n,\ell; r_1, \dots , r_{\ell})\to \C, \ g(x) \mapsto  [g^{\circ r_i}]^{\prime}(\beta_i).$$
Let  $Z^0(n,\ell; r_1, \dots , r_{\ell})$ be the set of all polynomials $g(x) \in Z(n,\ell; r_1, \dots , r_{\ell})$
such that the $\ell$-element set
$$\{\grad(\Mult_{\beta_i,r_i})_{\mid g(x)}\in \C^n \ \mid 1 \le i\le \ell\}$$
of gradients of $\Mult_{\beta_i,r_i}$'s at $g(x)$ is linearly
independent in $\C^n$  for every choice of $\{\beta_1, \dots,
\beta_{\ell}\}$. Clearly, $Z^0(n,\ell; r_1, \dots , r_{\ell})$  is
an open subset of $Z(n,\ell; r_1, \dots , r_{\ell})$ and therefore
of $\A^n$ in complex topology. However, this set may be empty; e.g.,
when $\ell\ge n$.
\end{sect}

The following statements are  main results of this paper.

\begin{thm}
 \label{open}
The set  $Z^0(n,\ell; r_1, \dots , r_{\ell})$  is  a Zariski-open subset of $Z(n,\ell; r_1, \dots , r_{\ell})$ and
therefore of $\A^n$.
\end{thm}

\begin{thm}
 \label{mainT}
Suppose that $n \ge 3$. Assume that $\sum_{i=1}^{\ell}r_i \le n$. If
$r_j=1$ for some $j$ then we assume additionally that
$\sum_{i=1}^{\ell}r_i < n$.

Then    $Z^0(n,\ell; r_1, \dots , r_{\ell})$ contains $g(x)=x^n$  and therefore is nonempty.
\end{thm}

Notice that under the notation and assumptions of Theorem
\ref{mainT}, if $r$ is a positive integer and $l(r)$ is the number
of $i$'s with $r_i=r$ then
$$l(r) \le \frac{n}{r} \le \nu(n,r).$$

Combining Theorems \ref{open} and \ref{mainT}, we obtain the
following statement.

\begin{cor}
\label{maincor} Suppose that $n \ge 3$. Assume that
$\sum_{i=1}^{\ell}r_i \le n$. If $r_j=1$ for some $j$ then we assume
additionally that $\sum_{i=1}^{\ell}r_i < n$.

Then $Z^0(n,\ell; r_1, \dots , r_{\ell})$ is a Zariski-open
everywhere dense  subset of $Z(n,\ell; r_1, \dots , r_{\ell})$ that
contains $g(x)=x^n$.
\end{cor}

\begin{ex}
 Suppose that $\ell = n-1$ and all $r_i=1$ (i.e., all the $\beta_i$ involved are  fixed points).
 It follows from results of \cite{ZarhinMatZametki2012,Rees}
that $$Z^0(n,n-1;1, \dots , 1)=Z(n,n-1;1, \dots ,1)=Z_{n,1}.$$
\end{ex}

\begin{rem}
\label{maprank}
 In the notation and assumptions of Corollary \ref{maincor}, let us consider the locally defined holomorphic map
$$Z(n,\ell; r_1, \dots , r_{\ell})\to \C^{\ell}$$
defined by the collection of functions
$\{\Mult_{\beta_i,r_i}\}_{i=1}^{\ell}$. Corollary \ref{maincor}
 asserts that this map has (maximal) rank $\ell$ on a nonempty Zariski-open
 subset
  $$Z^0(n,\ell; r_1, \dots , r_{\ell}) \subset Z(n,\ell; r_1, \dots , r_{\ell})$$
for every choice of periodic points $\{\beta_1, \dots ,
\beta_{\ell}\}$. It would be interesting to study its image. For the
case of fixed points (i.e., when all $r_i=1$), see \cite{Kulikov}.
\end{rem}

Notice that Remark \ref{maprank} gives a partial answer to a
question of Yu.S. Ilyashenko, who was interested in the case of two
orbits, in connection with \cite{Ilyash,Igor}.

\begin{rem}
In the case of two orbits,
 it turns out (see Examples \ref{n3r2}, \ref{rplus1} and Remark \ref{rplus2} below) that $g(x)=x^n$ does not belong
to $Z(n,2; r_1,r_2)$ if either $r_1=r_2=n-1$ or $r_1=1, r_2=n-1$. It would be interesting to find out whether in
these cases $Z(n,2; r_1,r_2)$ is empty.
\end{rem}

The paper is organized as follows.
In Section \ref{compute} we compute explicitly the differentials
$d\Mult_{\beta,r}$ at $g(x)=x^n$
where $\beta$ is a $(n^r-1)$th root of unity.
 This allows us to write down explicitly
the corresponding gradients $Q_r(\beta)\in\C^n$.
Now Theorem \ref{mainT} becomes equivalent to an assertion that
the corresponding set of vectors $\{Q_{r_i}(\beta_i)\}$ (and $Q_1(0)$ if one of $r_i$ is $1$)
 is linearly independent in $\C^n$.
We prove this assertion in Section \ref{ind}. Using standard
properties of finite maps (\cite[Ch. 1]{Sh},  \cite[Sect.
8]{Milne}), we prove Theorem \ref{open}  in Section \ref{global}.

\section{Computations of tangent maps}
\label{compute}

\begin{lem}
\label{tangentiter} Let us consider the holomorphic map $\Phi_{n,r}:
\A^n \to \A^{n^r}$ that sends a degree $n$ monic polynomial $g(x)$ to
the monic degree $n^r$ polynomial $g^{\circ r}(x)$. Then the tangent
map $d\Phi_{n,r}$ at $g(x)=x^n$ is as follows. It sends a tangent
vector $x^k$ (at the point $x^n$) to the tangent vector
$$p_{r,k}(x):=\sum_{i=1}^r n^{r-i} x^{n^r-n^i+n^{i-1}k}$$
(at the point $x^{n^r}$).
 In particular,
$$p_{r,0}(x):=\sum_{i=1}^r n^{r-i} x^{n^r-n^i}=n^{r-1}x^{n^r-n}+n^{r-2}x^{n^r-n^2}+\dots ,$$
$$p_{r,1}(x):=\sum_{i=1}^r n^{r-i} x^{n^r-n^i+n^{i-1}}=n^{r-1}x^{n^r-n+1}+n^{r-2}x^{n^r-n^2+n}+\dots $$
and $$\deg(p_{r,0})=n^r-n, \ \deg(p_{r,1})=n^r-n+1, \
\deg(p_{r,k})=n^r-n+k.$$
\end{lem}

\begin{proof}
Notice that $p_{1,k}(x)=x^k$ and for all positive integers $r$
$$p_{r+1,k}(x)=n x^{(n-1)n^r}p_{r,k}(x)+x^{k
n^r}.$$

{\sl Induction by} $r$.
 We need to prove that if $g^{[\epsilon]}(x)=x^n+\epsilon x^k$
then $[g^{[\epsilon]}]^{\circ r}(x)=x^{n^k}+\epsilon
p_{r,k}(x)+O(\epsilon^2)$. If $r=1$ then it is obvious. Assume that
this assertion is true for $r$ and let us check it for $r+1$. We
have
$$[g^{[\epsilon]}]^{\circ (r+1)}(x)=$$
$$\{[g^{[\epsilon]}]^{\circ r}(x)\}^n+\epsilon
\{[g^{[\epsilon]}]^{\circ r}(x)\}^k= (x^{n^r}+\epsilon
p_{r,k}(x)+O(\epsilon^2))^n+\epsilon (x^{n^r}+\epsilon
p_{r,k}(x)+O(\epsilon^2))^k=$$
$$x^{n^{r+1}}+\epsilon n
x^{(n-1)n^r}p_{r,k}(x)+\epsilon x^{k n^r}+O(\epsilon^2)=$$
$$x^{n^{r+1}}+\epsilon\{n x^{(n-1)n^r}p_{r,k}(x)+x^{k
n^r}\}+O(\epsilon^2)= x^{n^{r+1}}+\epsilon
p_{r,k+1}(x)+O(\epsilon^2).$$

\end{proof}

 Since all $p_{r,k}$ (for given $n$ and $r$) have distinct degrees, the
set $\{p_{r,0}, \dots , p_{r,n-1}\}$ is linearly independent. This
means that the rank of the tangent map to $\Phi_{n,r}$ at $g(x)=x^n$
is $n$, i.e. the tangent map at this point is injective and its
image coincides with
$$\oplus_{k=0}^{n-1}\C\cdot p_{r,k}.$$

\begin{sect}
\label{tangentMult} Suppose that $n^r \ge 3$.
Let us compute the differential
 $$d\Mult^{r}=d\M^{r}=d(\M \Phi_{n,r})=d\M \circ d \Phi_{n,r}$$
  at $g(x)=x^n \in Z_{n,r}$. Clearly,
  $$\Phi_{n,r}(x^n)=x^{n^r}\in P_m$$
  with $m=n^r$. Let $\alpha$ be a nonzero root of $x^m-x$, i.e., $\alpha^{n^r-1}=1$. Using Lemma \ref{tangentiter} and Example \ref{exm}, we obtain the following.
  The image
  $$q_{r,k}(\alpha):=d\Mult^{r}_{\mid g(x)=x^n}(x^k)$$
   of tangent vector $x^k$ to $g(x)=x^n \in Z_{n,r}$ is
  $$p_{r,k}^{\prime}(\alpha)-\frac{n^r p_{r,k}(\alpha)}{\alpha}=$$
 $$\alpha^{-1}\sum_{i=1}^r  (n^r-n^i+n^{i-1} k) n^{r-i}
\alpha^{n^r-n^i+n^{i-1}k}-\alpha^{-1}\sum_{i=1}^r n^r  n^{r-i}
\alpha^{n^r-n^i+n^{i-1}k}=$$
 $$\alpha^{n^r-1}\sum_{i=1}^r [n^{2r-i}-n^r+k
n^{r-1}-n^{2r-i}]\alpha^{-n^i+n^{i-1}k}=-(n^r-k n^{r-1})\sum_{i=1}^r
\alpha^{-n^i+n^{i-1}k}=$$
$$-(n-k)n^{r-1}\sum_{i=1}^r
\alpha^{-n^i+n^{i-1}k}=
-(n-k)n^{r-1}\sum_{i=1}^r \left(\frac{1}{\alpha^{n^{i-1}}}\right)^{n-k}.$$
In other words,
$$q_{r,k}(\alpha)=d\Mult^{r}_{\mid g(x)=x^n}(x^k)=-(n-k)n^{r-1}\sum_{i=1}^r \left(1/\alpha^{n^{i-1}}\right)^{n-k}.$$
Notice that
$$q_{r,k}(\alpha)=q_{r,k}(\alpha^n), \ q_{r,0}(\alpha)=n\cdot q_{r,n-1}(\alpha).$$
It follows that the {\sl gradient} of $\Mult^{r}$ at $g(x)=x^n$
(with respect  to $\alpha$) is
$$Q_{r}(\alpha)=(q_{r,0}(\alpha), q_{r,1}(\alpha), \dots, q_{r,n-1}(\alpha))=(n q_{r,n-1}(\alpha), q_{r,1}(\alpha), \dots, q_{r,n-1}(\alpha))\in \C^n.$$
Clearly,
$$Q_{r}(\alpha)=Q_{r}(\alpha^n)= \dots = Q_{r}(\alpha^{n^{r-1}}).$$

Let $\Oc(\alpha)=\{\alpha, \alpha^n, \dots , \alpha^{n^{r-1}}\}$ be the orbit of $\alpha$ with respect to $z\mapsto z^n$ and let $d(\alpha)$ be the cardinality of the set $\Oc(\alpha)$. Clearly, $d(\alpha)$ is a positive integer that divides $r$ and
$$\beta^{d(\alpha)}=\beta \ \forall \beta \in \Oc(\alpha).$$
It is also clear that
$$q_{r,k}(\alpha)=-\frac{r}{d(\alpha)}(n-k)n^{r-1}\sum_{\beta \in \Oc(\alpha)} \left(1/\beta\right)^{n-k}.$$
This implies that
$$Q_r(\alpha) =\frac{rn^{r}}{d(\alpha)n^{d(\alpha)}}\cdot Q_{d(\alpha)}(\alpha).$$
 \end{sect}

\begin{ex}
 \label{n3r2}
Suppose that $n=3$ and $r=2$. Then $n^r-1=8$. Let $\alpha$ be a $8$th root of unity that is not $\pm 1$.
Then $\alpha$ is a periodic point of exact period $2$ for $z\mapsto z^3$. We have
$$Q_2(\alpha)=
 -3^1\cdot \left(3\left[\frac{1}{\alpha}+\frac{1}{\alpha^3}\right],2\left[\frac{1}{\alpha^2}+\frac{1}{\alpha^6}\right],
      \left[\frac{1}{\alpha}+\frac{1}{\alpha^3}\right]\right).$$
So, if $\alpha$ is a primitive fourth root of unity, i.e.,
$$\alpha^2=-1, \ \alpha=\pm \mathbf{i},$$
then
$$\frac{1}{\alpha}+\frac{1}{\alpha^3}=0, \ \frac{1}{\alpha^2}+\frac{1}{\alpha^6}=-2$$ and
       $$Q_2(\alpha)=-3\cdot (0, -2,0)=(0,6,0).$$
 (Notice that $\mathbf{i}$ and $-\mathbf{i}$ lie in the same orbit.)

If $\alpha$ is a primitive $8$th root of unity then
$$\alpha^4=-1, \ 1+\frac{1}{\alpha^4}=0$$ and therefore
$$Q_2(\alpha)=
   -3\cdot \left(3\left[\frac{1}{\alpha}+\frac{1}{\alpha^3}\right],0,
      \left[\frac{1}{\alpha}+\frac{1}{\alpha^3}\right]\right)=
       -\frac{3}{\alpha^3}\cdot \left(3[\alpha^2+1],0, \alpha^2+1)\right)=
-\frac{3(\alpha^2+1)}{\alpha^3}\cdot \big(3,0,1\big).$$
 Now if we put $\beta=\alpha^{-1}$ then  $\alpha$ and $\beta$ are primitive $8$th roots of unity that do not belong
to the same orbit while $Q_2(\alpha)$ and $Q_2(\beta)$ generate the same line $\C\cdot (3,0,1)$ in $\C^3$. This implies
that $Z(3,2;2,2)$ does {\sl not} contain $g(x)=x^3$.

\end{ex}

\begin{ex}
\label{rplus1} Suppose that $n=r+1$ and $r\ge 2$. Then for all
positive integers $i$
$$n^i=(1+r)^i=1+ i \cdot r^{1}+ \dots + {i\choose j}r^{j} + \dots +r^i.$$
It follows that $n^i$ is congruent to $1+ir$ modulo $r^2$. In
particular, $n^r-1$ is divisible by $r^2=(n-1)^2$. It also follows
that   $(n^i-n)/(n-1)=n^{i-1}$ is congruent to $i-1$ modulo $r$ and
therefore
$$n^i-n \equiv (i-1)r \bmod r^2.$$

Suppose that $\alpha$ is a {\sl primitive} $r^2$th root of unity. Then $\alpha^{n^r-1}=1$ and therefore
$$\alpha^{n^r}=\alpha,$$
i.e., $\alpha$ is a periodic point for the map $z\mapsto z^n$. Clearly, its period divides $r$. On the other hand,
 for all positive integers $i<r$ the power $n^i$ is {\sl not} congruent to $1$ modulo $r^2$ and
therefore $\alpha^{n^i}\ne \alpha$. It follows that $\alpha$ has
exact period $r$.

The number $\gamma:=\alpha^{1-n}=\alpha^{-r}$ is a primitive $r$th
root of unity. For each integer $k$ with $0\le k\le n-1$ the number
$\delta=\gamma^{n-k}$ is an $r$th root of unity. Clearly, $\delta
\ne 1$ if and only if $n-k \ne n-1$, i.e. $k \ne 1$.  In particular,
if $k \ne 1$ then $\sum_{i=1}^r \delta^i=0$.

We have
$$q_{r,k}(\alpha)=-(n-k)n^{r-1} \sum_{i=1}^r
\left(\alpha^{-n^i}\right)^{n-k}=$$
$$-(n-k)n^{r-1} \cdot\alpha^{n(k-n)}\sum_{i=1}^r
\left(\alpha^{n-n^i}\right)^{n-k}=-(n-k)n^{r-1}\alpha^{n(k-n)}\sum_{i=1}^r
\delta^{i-1}=$$
$$-(n-k)n^{r-1}\delta^{-1}\alpha^{n(k-n)}\sum_{i=1}^r \delta^{i}=
0$$ if $k \ne 1$. On the other hand, if $k=1$ then $\delta=1$ and
$q_{r,1}(\alpha)=-r(n-1)n^{r-1}\gamma$. It follows that
$$Q_r(\alpha)=-(0,r(n-1)n^{r-1}\gamma, 0, \dots ,0)=-r(n-1)n^{r-1}\gamma \cdot Q_1(0)\in \C^n.$$
This implies that $Z^0(r+1,2;r,1)$ does {\sl not} contain
$g(x)=x^n$.
\end{ex}

\begin{rem}
\label{rplus2}
Suppose that $r >2$ and $n=r+1$.  Example \ref{rplus1} tells us
 that if $\alpha$ and $\beta$ are
 primitive $r^2$th roots of unity then $Q_r(\alpha)$ and $Q_r(\beta)$
generate the same line $\C \cdot (0,1, \dots , 0)$ in $\C^n$.
 Since $r>2$, the number $\varphi(r^2)$ of primitive $r^2$th
roots of unity is strictly greater than $r$. (Here $\varphi$ is the
Euler function.) In particular, we may choose such $\alpha$ and
$\beta$ from different orbits (of length $r$) of the map $z \mapsto
z^n$. It follows that $Z^0(r+1,2;r,r)$ does {\sl not} contain
$g(x)=x^n$.
\end{rem}

\section{Linear independence}
\label{ind}
As was already pointed out,
Theorem \ref{mainT} is an immediate corollary of the following statement.

\begin{thm}
\label{bign}
 Let   $\ell$ be a positive integer. Let $\{r_1, \dots r_{\ell}\}$ be a sequence of $\ell$ positive integers.
  Let $\{\alpha_1, \dots, \alpha_{\ell}\}$ be a sequence of distinct complex numbers such that
  $$\alpha_i^{n^{r_i}-1}=1 \ \forall i=1, \dots, \ell.$$
  Assume that $\{\alpha_1, \dots, \alpha_{\ell}\}$
 belong to different orbits of the map $z \mapsto z^n$. Then:

  \begin{itemize}
  \item[(i)]
  the set of $\ell$ vectors
 $\{Q_{r_1}(\alpha_1), \dots , Q_{r_{\ell}}(\alpha_{\ell})\}$ in $\C^n$ is linearly independent if
 $n \ge \sum_{i=1}^{\ell}  d(\alpha_i)$. In particular, if $\sum_{i=1}^{\ell} r_i \le n$ then the $\ell$-tuple
 $\{Q_{r_1}(\alpha_1), \dots , Q_{r_{\ell}}(\alpha_{\ell})\}$  is linearly independent in $\C^n$.
 \item[(ii)]
 If $n \ge 2+ \sum_{i=1}^{\ell}  d(\alpha_i)$ then the $(\ell+1)$-tuple
 $\{Q_1(0); Q_{r_1}(\alpha_1), \dots , Q_{r_{\ell}}(\alpha_{\ell})\}$  is a linearly independent set in $\C^n$.
In particular, if $\sum_{i=1}^{\ell} r_i < n-1$ then the
$(\ell+1)$-tuple $\{Q_1(0); Q_{r_1}(\alpha_1), \dots ,
Q_{r_{\ell}}(\alpha_{\ell})\}$ is a linearly independent set in
$\C^n$.
 \end{itemize}
 \end{thm}

\begin{proof}[Proof of Theorem \ref{bign}]

In the course of the proof we will use the following elementary statement that will be proven at the end of this section.

\begin{lem}
\label{newton}
Let $d$ be a positive integer,  $S$ a set of $d$ nonzero complex numbers.
Let $c:S \to \C$ be a function such that for all positive integers $u=1, \dots ,d$
 $$\sum_{\beta \in S}\frac{c(\beta)}{\beta^u}=0.$$
Then $c(\beta)=0$ for all $\beta\in S$.

\end{lem}

Let us continue to prove Theorem \ref{bign}. Replacing each $r_i$ by
$d(\alpha_i)$ we may and will assume that $r_i=d(\alpha_i)$, i.e.,
the orbit $\Oc(\alpha_i)$ of $\alpha_i$ consists of $r_i$ distinct
elements (for all $i$). We also assume that $n \ge
\sum_{i=1}^{\ell}r_i$.

Let $\{c_1,  \dots c_{\ell}\}$ be a sequence of $\ell$ complex numbers  such that
 $$\sum_{i=1}^{\ell} c_i Q_{r_i}(\alpha_i)=0.$$
 Let $S\subset \C$ be the (disjoint) union of all $\Oc(\alpha_i)$, which consists of $\left(\sum_{i=1}^{\ell} r_i\right)$ elements. Let us define a complex valued function
 $c$ on $S$ that assigns to $\alpha \in \Oc(\alpha_i)$ the complex number
 $$c(\alpha):=n^{r_i-1}c_i.$$
 Then we obtain for all $k=0,1, \dots , n-1$
 $$0=\sum_{i=1}^{\ell} c_i q_{r_i,k}(\alpha_i)=-(n-k)\sum_{\alpha\in S}c(\alpha) (1/\alpha)^{n-k}.$$
 This implies that
 $$\sum_{\alpha\in S}c(\alpha) (1/\alpha)^{u}=0$$
 for all positive integers  $u=1, \dots, n$.

 It follows from Lemma \ref{newton} applied to $d=\sum_{i=1}^{\ell}r_i$  that
 all $c(\alpha)=0$. Since all $n^{r_i-1}\ne 0$, we conclude that all $c_i=0$. This proves (i).

 Now assume that $n \ge 2+\sum_{i=1}^{\ell}r_i$. We are going to prove (ii).
 Let $\{c_0, c_1,  \dots c_{\ell}\}$ be a sequence of $(\ell+1)$ complex numbers  such that
 $$c_0 Q_1(0)+\sum_{i=1}^{\ell} c_i Q_{r_i}(\alpha_i)=0.$$
 We have
 $$ -c_0 Q_1(0)=\sum_{i=1}^{\ell} c_i Q_{r_i}(\alpha_i).$$
 Recall that all the coordinates of $Q_1(0)$ except the second one do vanish. This implies that
 $$0=\sum_{i=1}^{\ell} c_i q_{r_i,k}(\alpha_i)=-(n-k)\sum_{\alpha\in S}c(\alpha) (1/\alpha)^{n-k}$$ for all $k=0, \dots n-1$ except $k=1$. It
 follows
 that
 $$\sum_{\alpha\in S} c(\alpha) (1/\alpha)^{u}=0$$
 for all positive integers integers $u=1, \dots, n-2$. Since $n-2\ge d$, the same arguments with Lemma \ref{newton} as above prove that $c_i=0$ for all positive integers $i$ and therefore
 $-c_0 Q_1(0)=0$, i.e., $c_0=0$.
\end{proof}

\begin{proof}[Proof of Lemma \ref{newton}]
This Lemma is a variant of well-known classical results (e.g., see \cite{Yomdin}).
Let us consider the rational function
$$X(t)=\sum_{\beta \in S} \frac{c(\beta)}{\beta - t}=
\sum_{\beta \in S} \frac{c(\beta)/\beta}{1-\frac{t}{\beta}}.$$
Clearly,
$$X(t)=\frac{Q(t)}{\prod_{\beta\in S}(\beta- t)}$$
where $Q(t)$ is a polynomial, whose degree does not exceed $d-1$. (Recall that $d=\#(S)$.)
 For each positive integer $u$, the number $\sum_{\beta\in S} c(\beta)/ \beta^u$ is
 the $(u-1)$th coefficient of the Taylor power series of $X(t)$ at the origin
  (see \cite[Ch. 1, Sect. 2]{Mc}). It follows that $X(t)$ has a zero of order $\ge d$ at the origin.
This implies that $Q(t)$ is divisible by $t^{d}$ and therefore $Q(t)=0$, i.e. $X(t)=0$.
Since $-c(\beta)$ is the residue of $X(t)$ at $t=\beta$ for all $\beta \in S$, we conclude that $c(\beta)=0$.
\end{proof}

\begin{rem}
 One may give even more elementary proof of Lemma \ref{newton},
using the nondegeneracy of the Vandermonde matrix of size $d \times d$ for {\sl distinct} numbers
$\{1/\beta \mid \beta \in S\}$.
\end{rem}

\section{Openness in Zariski topology}
\label{global}

The aim of this Section is to prove Theorem \ref{open}.
We will need the following well known  easy statement.

\begin{lem}
\label{div}
 Let $n \ge 2$ be an integer and $g(x)\in \C[x]$ is a monic degree $n$ polynomial.
Suppose that $g(x)-x$ has a multiple root say, $\alpha$. Then for all positive integers $r$ the complex number
$\alpha$ is a multiple root of $g^{\circ r}(x)-x$.
\end{lem}

\begin{proof}
 We have
$$g(\alpha)=\alpha, \ g^{\prime}(\alpha)=1.$$
It follows easily that
$$g^{\circ r}(\alpha)=\alpha, \ [g^{\circ r}]^{\prime}(\alpha)=1.$$
This means that
$$[g^{\circ r}(x)-x](\alpha)=0, \ [g^{\circ r}(x)-x]^{\prime}(\alpha)=0.$$
In other words, $\alpha$ is a multiple root of $g^{\circ r}(x)-x$.
\end{proof}

\begin{cor}
\label{divcor} Let $m$ be a positive integer that divides $r$.
Suppose that $g^{\circ m}(x)-x$ has a multiple root say, $\alpha$.
Then $\alpha$ is a multiple root of $g^{\circ r}(x)-x$.
\end{cor}

\begin{sect}
\label{numberoforbits}
 Let $g(x)\in Z_{n,r}$. If $m$ is a positive
divisor of $r$ then Corollary \ref{divcor} implies that $g(x) \in
Z_{n,m}$. The number of periodic points of exact period $m$ for
$z\mapsto g(z)$ is
$$\nu_n(m)=\sum_{m^{\prime}\mid m} \mu\left(\frac{m}{m^{\prime}}\right)n^{m^{\prime}}$$
where $\mu$ is the M\"obius function \cite[pp. 74--75]{Milnor}. In particular, the number of orbits of length $m$ equals
$$d(n, m)=\frac{\nu_n(m)}{m}$$
and therefore $\nu_n(m)$ is divisible by $m$. The explicit formula
for $\nu_n(m)$ implies that $\nu_n(m)$ is also divisible by $n$
(ibid). It follows that $\nu_n(m)$ is divisible by $nm/(n,m)$ where
$(n,m)$ is the greatest common divisor  of $n$ and $m$. On the other
hand, the number $\nu_n(m)$ is always positive. Indeed, $\nu_n(1)=n$
and for $1<m\le 5$ we have $\nu_n(m)=n^m-n^{m/p}>0$ where $p$ is the
only prime divisor of $m$. Now assume that $m>5$, i.e, $m \ge 6$.
Notice that the points of exact period $m$ are exactly the roots of
$g^{\circ m}(x)-x$ that are not roots of $g^{\circ (m/p)}(x)-x$ for
any prime divisor $p$ of $m$. Since the number of prime divisors of
$m$ does not exceed $\log_2(m)$,
$$\nu_n(m) \ge n^{m} - n^{m/2}\log_2(m).$$
It is easy to check that (under our assumptions on $n$ and $m$) we
have $n^{m/2}\ge 2^{m/2}>\log_2(m)$ and therefore $\nu_n(m)$ is
positive.

Since  $\nu_n(m)$ is divisible by $\frac{nm}{(n,m)}$,
$$\nu_n(m) \ge \frac{nm}{(n,m)}, \ d(n,m) \ge \frac{n}{(n,m)} \ge
\frac{n}{m}.$$
\end{sect}

\begin{sect}
For each positive divisor $m$ of $r$ we pick a $d(n, m)$-element set $S_m$ and consider the corresponding
 $d(n, m)$-dimensional coordinate space $\C^{S_m}$ of all $\C$-valued functions on $S_m$.

 Let us consider the Zariski-closed subset
$$\hat{Z}_{n,r} \subset Z_{n,r} \times \prod_{m\mid r}\C^{S_m}$$
that is cut out by the following equations imposed on
$$\{g; \phi_m:S_m \to \C, \ m\mid r \}\in Z_{n,r} \times \prod_{m\mid r}\C^{S_m}.$$

For each $s \in S_m$ the complex number $\phi_m(s)$ is a periodic
point, whose period divides  $m$, with respect to  $z \mapsto g(z)$,
i.e. $g^{\circ m}(\phi_m(s))=\phi_m(s)$. In addition, we require
that
$$g^{\circ r}(x)-x=
\prod_{m\mid r} \left(\prod_{s \in S_m} (x-\phi_m(s))
\prod_{i=1}^{m-1} \left(x-g^{\circ i}(\phi_m(s))\right)\right).$$ In
other words, the coefficients of both polynomials in $x$ do
coincide. Notice that
$$g^{\circ r}(\phi_m(s))-\phi_m(s)=0$$
on $\hat{Z}_{n,r}$. In particular, all the coordinate functions
$\phi_m(s)$ on $\hat{Z}_{n,r}$ are integral over the polynomial ring
$\C[\A^n]=\C[a_0, \dots , a_{n-1}]$, which is generated by the
coefficients of $g(x)=x^n + \sum_{i=0}^{n-1} a_i x^i$.

Recall that $g^{\circ r}(x)-x$ has no multiple roots. It follows
that each map $\phi_m:S_m \to \C$ is injective, its image consists
of elements of exact period $m$ while distinct elements of $S_m$ go
under $\phi_m$ to distinct orbits of length $m$; in addition, every
orbit of length $m$ contains exactly one element of $\phi_m(S_m)$.
On the other hand, for any choice of an element $\zeta$ in every
orbit of length $m$ (for each divisor $m$ of $r$) there is (exactly
one) point of $\hat{Z}_{n,r}$ that lies ``above" $g(x)$ and such
that the corresponding $\phi_m(S_m)$ consists of these $\zeta$.

By construction, the projection map of affine varieties
$\hat{Z}_{n,r} \to Z_{n,r}$ is surjective. In addition, this map
 is finite, because the ring of regular functions
 $\C[\hat{Z}_{n,r}]$ on $\hat{Z}_{n,r}$ is generated over
 $\C[Z_{n,r}]\supset \C[\A^n]=\C[a_0, \dots , a_{n-1}]$ by the coordinate functions
$\phi_m(s)$ that are integral over $\C[\A^n]$.

 Now one may ``lift'' $\Mult^r=\Mult_{r,\beta}$ to globally defined
functions on $\hat{Z}_{n,r}$. Namely, for each $s\in S_r$ the
function
$$\overline{\Mult}_{r,s}:\hat{Z}_{n,r} \to \C, \{g; \phi_m:S_m \to \C, \ m\mid r \}\mapsto \Mult_{r,\phi_r(s)}(g)$$
is a globally defined regular function.  If $\phi_m(s)$ and $\beta$
lie in the same orbit of length $m$ then this function coincides
with the composition of projection map $\hat{Z}_{n,r} \to Z_{n,r}$
and $\Mult_{r,\beta}$.
 It is also clear that the vector function
$$\overline{\grad}(\Mult_{r,s}) :\hat{Z}_{n,r} \to \C^n, \ \{g; \phi_m:S_m \to \C^n, \ m\mid r \}\mapsto
\grad(\Mult_{r,\phi_r(s)})_{\mid g(x)}$$ is a regular map that
coincides with the composition of projection map $\hat{Z}_{n,r} \to
Z_{n,r}$  and $\grad(\Mult_{r,\beta}): Z_{n,r}\to \C^n$  with
$\beta=\phi_r(s)$.

Let $\ell$ be a positive integer that does not exceed $d(n,r)$. If
$D$ is an $\ell$-element subset of $S_r$ let us consider the the
subset $X_D$ of points $v \in \hat{Z}_{n,r}$ such that the
collection of $\ell$ vectors $\{\overline{\grad}(\Mult_{r,s})(v)\mid
s\in D\}$ is linearly dependent in $\C^n$. Clearly, $X_D$ is a
Zariski-closed subset in $\hat{Z}_{n,r}$. It follows that the union
$X$ of all $X_D$ (where $D$ runs through all $\ell$-element subsets
of $S_r$) is also closed in  $\hat{Z}_{n,r}$. The finiteness of the
projection map implies that the image $\bar{X}$ of $X$ in $Z_{n,r}$
is also Zariski-closed (\cite[Ch. 1, Sect. 5.3]{Sh}). On the other
hand, one may easily check that $\bar{X}$ is the complement of
$Z^0(n,\ell; r, r, \dots ,r)$ in $Z(n,\ell; r,r \dots, r)=Z_{n,r}$.
It follows that $Z^0(n,\ell; r,r \dots, r)$ is Zariski-open in
$Z_{n,r}$. This proves Theorem \ref{open} in the case when $r_1=r_2=
\dots = r_{\ell}$.
\end{sect}

\begin{sect}
Now let us consider the general case. Let $d$ be the number of distinct elements in
 the sequence $\{r_1, \dots ,r_{\ell}\}$ and
$R$ the corresponding $d$-element set of positive integers. For each
$r\in R$ we denote by $l(r)$ the number of $i$ with $r_i=r$. If
there is $r$ with $l(r)> d(n,r)$ then $Z(n,\ell; r_1, r_2, \dots
,r_{\ell})$ is empty. So further we assume that $l(r)\le d(n,r)$ for
all $r\in R$.

We write $\hat{Z}_{n,r}^{\prime}$ for the preimage of $Z(n,\ell;
r_1, r_2, \dots ,r_{\ell})$ in $\hat{Z}_{n,r}$. The natural
 regular map
$$\hat{Z}_{n,r}^{\prime}\to Z(n,\ell; r_1, r_2, \dots ,r_{\ell})$$
  is finite, because the base change of a finite map is finite \cite[Sect. 8, Prop. 8.22]{Milne}. Let us consider the fiber product
$\hat{Z}_{n,R}$ of all $\hat{Z}_{n,r}^{\prime}$ (where $r$ runs
through $R$) over $Z(n,\ell; r_1, r_2, \dots ,r_{\ell})$. If we
write $\pi$ for the natural regular map $$\hat{Z}_{n,R}\to Z(n,\ell;
r_1, r_2, \dots ,r_{\ell})$$ then $\pi$ is finite, because the
composition  of finite maps is also finite \cite[Sect. 8, Prop. 8.4
and 8.22]{Milne}. We denote by $\pi_r$
 the natural finite regular map
$$\pi_r:\hat{Z}_{n,R} \to \hat{Z}_{n,r}^{\prime}.$$
Now for each $r \in R$ pick a $l(r)$-element subset $D_r$ in $S_r$
and consider the subset $X_{\{D_r\mid r\in R\}}$ of points $v \in
\hat{Z}_{n,R}$ such that the collection of $\left(\sum_{r\in
R}l(r)\right)$ vectors
$$\{\overline{\grad}(\Mult_{r,s_r})(v_r)\mid s_r\in D_r, \ r \in R\}$$
 is linearly dependent in $\C^n$.
Here $v_r=\pi_r(v)\in \hat{Z}_{n,r}^{\prime}$. Clearly,
$X_{\{D_r\mid r\in R\}}$ is a closed algebraic subvariety in
$\hat{Z}_{n,R}$. Now let us consider the union $Y$ of all such
$X_{\{D_r\mid r\in R\}}$ for all choices of $\{D_r\mid r\in R\}$. It
is also clear that $Y$ is also a closed algebraic subvariety in
$\hat{Z}_{n,R}$. Since  $\pi$ is finite,
  the image $\pi(Y)$ is a closed algebraic subvariety in
$Z(n,\ell; r_1, r_2, \dots ,r_{\ell})$. On the other hand, one may easily check
that $Z^{0}(n,\ell; r_1, r_2, \dots ,r_{\ell})$ is the complement of $\pi(Y)$ in $Z(n,\ell; r_1, r_2, \dots ,r_{\ell})$.
It follows that $Z^{0}(n,\ell; r_1, r_2, \dots ,r_{\ell})$ is Zariski-open in $Z(n,\ell; r_1, r_2, \dots ,r_{\ell})$.
This ends the proof of Theorem \ref{open}.
\end{sect}

{\bf Acknowledgements.} I am grateful to Yulij S. Ilyashenko,
Vladimir L. Popov and Victor S. Kulikov  for  stimulating
discussions. My special thanks go to Tatiana Bandman, whose comments
helped to improve the exposition.

 Part of this work was done during my stay at
Weizmann Institute of Science Department of Mathematics in May--June
of 2012 and at Centre Interfacultaire Bernoulli (\'Ecole
Polytechnique F\'ed\'erale de Lausanne) in July--August of 2012:
 I am grateful to both of them  for the hospitality.