text:
\begin{document}

\title{Derivative-free global minimization for a class of multiple minima problems}

\titlerunning{Derivative-free global minimization}

\author{Xiaopeng Luo \and Xin Xu \and Daoyi Dong}

\institute{\Letter~Xin Xu \\
\email{xu.permanent@gmail.com} \\
\at Department of Control and Systems Engineering, School of Management and Engineering, Nanjing University, Nanjing, 210008, China \\
Department of Chemistry, Princeton University, Princeton, NJ 08544, USA
\and
Xiaopeng Luo \\
\email{luo.permanent@gmail.com} \\
\at Department of Control and Systems Engineering, School of Management and Engineering, Nanjing University, Nanjing, 210008, China \\
Department of Chemistry, Princeton University, Princeton, NJ 08544, USA
\and
Daoyi Dong \\
\email{daoyidong@gmail.com} \\
\at School of Engineering and Information Technology, University of New South Wales, Canberra, 2600, Australia
}

\date{}

\maketitle

\begin{abstract}
  We prove that the finite-difference based derivative-free descent (FD-DFD) methods have a capability to find the global minima for a class of multiple minima problems. Our main result shows that, for a class of multiple minima objectives that is extended from strongly convex functions with Lipschitz-continuous gradients, the iterates of FD-DFD converge to the global minimizer $x_*$ with the linear convergence $\|x_{k+1}-x_*\|_2^2\leqslant\rho^k \|x_1-x_*\|_2^2$ for a fixed $0<\rho<1$ and any initial iteration $x_1\in\mathbb{R}^d$ when the parameters are properly selected. Since the per-iteration cost, i.e., the number of function evaluations, is fixed and almost independent of the dimension $d$, the FD-DFD algorithm has a complexity bound $\mathcal{O}(\log\frac{1}{\epsilon})$ for finding a point $x$ such that the optimality gap $\|x-x_*\|_2^2$ is less than $\epsilon>0$. Numerical experiments in various dimensions from $5$ to $500$ demonstrate the benefits of the FD-DFD method.

\keywords{multiple minima problem \and global minima \and nonconvex \and derivative-free \and convergence rate \and complexity}

\subclass{65K05 \and 68Q25 \and 90C26 \and 90C56}
\end{abstract}

\section{Introduction}
\label{DFD:s1}

Derivative-free descent (DFD) methods, also known as zero-order methods \cite{DuchiJ2015A_ZeroOrderCO,ShamirO2017A_ZeroOrderConvex} in the literature or bandit optimization in the machine learning literature \cite{HazanE2014A_BanditOptimization,ShamirO2017A_ZeroOrderConvex}, do not require the availability of derivatives. They are generally applied to instances where derivatives are unavailable or unreliable \cite{ConnA2009M_DerivativeFree,PoliakB1987M_Optimization}. DFD methods do not depend directly on gradient information, but some of them are related to gradient estimates, e.g., the finite-difference based derivative-free descent (FD-DFD) methods \cite{NesterovY2017A_GF,PoliakB1987M_Optimization}.

The FD-DFD method could be regarded as a smoothed extension of the gradient method because an FD-DFD descent direction is an unbiased estimate of smoothed gradient at the current iterate \cite{NemirovskiA1983M_optimization,NesterovY2017A_GF}. It also has some characteristics of the gradient method \cite{NesterovY2017A_GF} but its ``cognitive range'' is closely related to the smoothing parameter. One may expect that the FD-DFD method has certain global search capability when the smoothing and stepsize parameters are properly selected. Therefore, we attempt to analyze the convergence behavior of the FD-DFD method in global optimization under certain conditions, regardless of whether the derivative is available.

Recently, a regularized asymptotic descent (RAD) method \cite{LuoX2020A_RAD} was proposed to find the global minima with linear convergence and logarithmic work complexity for certain class of multiple minima problems. It is inspired by an asymptotic solution of the regularized minimization problems which is extended from the Pincus asymptotic solution formula \cite{PincusM1968A_AsymptoticSolution,PincusM1970A_AsymptoticSolution}. Under a mild assumption, the RAD iterates will converge to the global minimizer without being trapped in saddle points, local minima, or even discontinuities.  In this work, we will prove that the FD-DFD method also enjoys a similar convergence behavior under the same assumption.

We will see that this global linear convergence comes from the smoothing effect on objective gradients. Smoothing techniques have been extensively studied \cite{BurkeJ2020A_Max,DuchiJ2012A_SGsmoothing,NesterovY2017A_GF}, especially for nonsmooth objectives
\cite{ChenX2012A_Smoothing,NesterovY2005A_Smooth}. One may notice that a smoothing method is different from a random noise based approach, e.g., random search method \cite{BergstraJ2012M_RS_HyperParameter} or perturbed gradient method \cite{JinC2017A_EscapeSaddlePoints,PemantleR1990A_SaddlePoints}. The reason is that a smoothed gradient contains certain global information but a random noise is independent of the objective.

Specifically, we analyze the FD-DFD method for finding the global minima
\begin{equation}\label{RAD:eq:P}
  x_*=\arg\min_{x\in\mathbb{R}^d}f(x),
\end{equation}
where the objective function $f:\mathbb{R}^d\to\mathbb{R}$ satisfies the following assumption:
\begin{assum}\label{DFD:ass:A}
The objective function $f:\mathbb{R}^d\to\mathbb{R}$ satisfies there exist $x_*\in\mathbb{R}^d$ and $0<l\leqslant L<\infty$ such that for all $x\in\mathbb{R}^d$,
\begin{equation}\label{DFD:eq:A}
  f_*+\frac{l}{2}\|x-x_*\|_2^2\leqslant f(x)
    \leqslant f_*+\frac{L}{2}\|x-x_*\|_2^2.
\end{equation}
Hence, $f$ has a unique global minimizer $x_*$ with $f_*:=f(x_*)$.
\end{assum}

\begin{figure}[!h]
\centering
\includegraphics[width=0.4\textwidth]{fig1a}
\includegraphics[width=0.4\textwidth]{fig1b}
\caption{One-dimensional examples. Left: the objective $f(x)=x^2+x^2\cos(5\pi x)/2$ (solid) with lower bound $x^2/2$ (dashed) and upper bound $3x^2/2$ (dash-dotted). Right: the objective $f(x)=x^2-\cos(5\pi x)/2+1/2$ (solid) with lower bound $x^2$ (dashed) and upper bound $65x^2$ (dash-dotted).}
\label{DFD:fig:A}
\end{figure}

As shown in Fig. \ref{DFD:fig:A}, such a class of functions can be extended from strongly convex functions with Lipschitz-continuous gradients; however, it is not ruling out the possibility of multiple minima. The lower bound $f_*+\frac{l}{2}\|x-x_*\|_2^2$ guarantees the uniqueness of the global minima while the upper bound $f_*+\frac{L}{2}\|x-x_*\|_2^2$ controls the sharpness of the minima.

Under Assumption \ref{DFD:ass:A}, the Lipschitz-continuous objective $f$ has a unique global minima $x_*$ and possibly multiple local minima. Our goal here is to find this global minima $x_*$ without being trapped in saddle points or local minima. We prove that the FD-DFD method enjoys global linear convergence, i.e., $\|x_{k+1}-x_*\|_2^2=\rho^k\|x_1-x_*\|_2^2$, for finding the global minimizer $x_*$ when the parameters are properly selected (see Theorem \ref{DFD:thm:main2}).

The remainder of the paper is organized as follows. In Sect. \ref{DFD:s2}, we establish the convergence property and complexity bound
for the FD-DFD method. In Sect. \ref{DFD:s3}, we compare the characteristics of the FD-DFD method with the RAD method. In Sect. \ref{DFD:s4}, we demonstrate the benefits of the FD-DFD method by numerical experiments in various dimensions from $5$ to $500$. And finally, we draw conclusions in Sect. \ref{DFD:s5}.

\section{Derivative-free methods}
\label{DFD:s2}

\subsection{Gradient of Gaussian smoothing}

Let random vector $\xi$ have $d$-dimensional standard normal distribution, that is, $\xi\sim\mathcal{N}(0,I_d)$. Denote by $\mathbb{E}_\xi(h(\xi))$ the expectation of corresponding random variable. For any smoothing parameter $\sigma>0$, we consider
\begin{equation*}
  g_\sigma(x,\xi)=
  \frac{\big(f(x+\sigma\xi)-f_*\big)\xi}{\sigma}
\end{equation*}
as an unbiased estimate for the gradient of the smoothed objective function
\begin{equation}\label{DFD:eq:SO}
  f_\sigma(x)=\mathbb{E}_\xi\big[f(x+\sigma\xi)\big].
\end{equation}
In practice, $f_*$ and can be replaced with corresponding estimate.

Using the substitution $\xi=\frac{\theta-x}{\sigma}$, this gradient estimate can also be written as \begin{equation*}
  g_\sigma(x,\theta)=\frac{\big(f(\theta)-f_*\big)(\theta-x)}{\sigma^2},
  ~~~\theta\sim\mathcal{N}(x,\sigma^2I_d).
\end{equation*}
See \cite{NesterovY2017A_GF,PoliakB1987M_Optimization} for other different finite-difference based schemes.

Theorem \ref{SGFD:thm:stepsmoothing} establishes $\nabla f_\sigma(x)= \frac{1}{\sigma}\mathbb{E}_\xi[f(x+\sigma\xi)\xi]= \frac{1}{\sigma^2}\mathbb{E}_\theta[f(\theta)(\theta-x)]$, thus, we obtain the unbiasedness
\begin{equation*}
  \mathbb{E}_\xi[g_\sigma(x,\xi)]=
  \mathbb{E}_\theta[g_\sigma(x,\theta)]=\nabla f_\sigma(x).
\end{equation*}

\begin{theorem}\label{SGFD:thm:stepsmoothing}
Suppose that $f_\sigma:\mathbb{R}^d\to\mathbb{R}$ is the smoothed objective with a smoothing parameter $\sigma>0$ defined by \eqref{DFD:eq:SO}, then its gradient
\begin{equation}\label{DFD:eq:Gsmoothing}
  \nabla f_\sigma(x)=\frac{1}{\sigma}
  \mathbb{E}_\xi[f(x+\sigma\xi)\xi]
  =\frac{1}{\sigma^2}\mathbb{E}_\theta[f(\theta)(\theta-x)],
\end{equation}
where $\xi\sim\mathcal{N}(0,I_d)$ and $\theta\sim\mathcal{N}(x,\sigma^2I_d)$.
\end{theorem}
\begin{proof}
Using the substitution $\xi=\frac{\theta-x}{\sigma}$, we obtain
\begin{align*}
  f_\sigma(x)=\mathbb{E}_\xi\big[f(x+\sigma\xi)\big]
  =\frac{1}{(\sqrt{2\pi}\sigma)^d}\int_{\mathbb{R}^d}
  f(\theta)e^{-\frac{\|\theta-x\|_2^2}{2\sigma^2}}\ud\theta,
\end{align*}
and therefore,
\begin{align*}
  \nabla f_\sigma(x)\!=&\!\frac{1}{(\sqrt{2\pi}\sigma)^d}
  \!\!\int_{\mathbb{R}^d}\!f(\theta)\nabla_x
  e^{-\frac{\|\theta-x\|_2^2}{2\sigma^2}}\!\ud\theta \\
  \!=&\!\frac{1}{(\sqrt{2\pi}\sigma)^d}
  \!\!\int_{\mathbb{R}^d}\!f(\theta)\!
  \left(\!\frac{\theta\!-\!x}{\sigma^2}\!\right)
  e^{-\frac{\|\theta-x\|_2^2}{2\sigma^2}}\!\ud\theta.
\end{align*}
That is, $\nabla f_\sigma(x)=\frac{1}{\sigma^2} \mathbb{E}_\theta[f(\theta)(\theta-x)]$. Using the substitution $\theta=x+\sigma\xi$, we further obtain
\begin{equation*}
  \nabla f_\sigma(x)=\frac{1}{\sigma}\frac{1}{(\sqrt{2\pi})^d}
  \int_{\mathbb{R}^d}f(x+\sigma\xi)\xi
  e^{-\frac{\|\xi\|_2^2}{2}}\ud\xi
  =\frac{1}{\sigma}\mathbb{E}_\xi[f(x+\sigma\xi)\xi],
\end{equation*}
and the proof is complete. \qed
\end{proof}

We will see that the smoothed gradient $\nabla f_\sigma$ plays a key role in the subsequent analysis. At the same time, we also notice that the Gaussian smoothing does not make any significant changes to the bounds of the objective function $f$. Specifically, under Assumption \ref{DFD:ass:A}, we have,
\begin{equation*}
  f_*+\frac{l}{2}\|x-x_*\|_2^2\leqslant f(x)\leqslant
  f_*+\frac{L}{2}\|x-x_*\|_2^2,
\end{equation*}
together with
\begin{align*}
  &\frac{1}{(\sqrt{2\pi}\sigma)^d}\int_{\mathbb{R}^d}
  \|\theta-x_*\|_2^2e^{-\frac{\|\theta-x\|_2^2}{2\sigma^2}}\ud\theta \\
  =&\frac{1}{(\sqrt{2\pi}\sigma)^d}\sum_{i=1}^d\int_{\mathbb{R}^d}
  \big(\theta^{(i)}-x^{(i)}+x^{(i)}-x_*^{(i)}\big)^2
  e^{-\frac{\|\theta-x\|_2^2}{2\sigma^2}}\ud\theta \\
  =&\frac{1}{(\sqrt{2\pi}\sigma)^d}\sum_{i=1}^d\int_{\mathbb{R}^d}
  \big(\theta^{(i)}-x^{(i)}\big)^2
  e^{-\frac{\|\theta-x\|_2^2}{2\sigma^2}}\ud\theta
  +\sum_{i=1}^d\big(x^{(i)}-x_*^{(i)}\big)^2 \\
  =&d\sigma^2+\|x-x_*\|_2^2,
\end{align*}
we obtain the bounds
\begin{equation*}
  f_*+\frac{l}{2}\|x-x_*\|_2^2+\frac{d\sigma^2l}{2}
  \leqslant f_\sigma(x)\leqslant
  f_*+\frac{L}{2}\|x-x_*\|_2^2+\frac{d\sigma^2L}{2}.
\end{equation*}

\subsection{Algorithms}

With an initial point $x_1$, a stepsize $\alpha>0$, an initial exploration radius $\lambda>0$, a fixed contraction factor $0<\rho<1$ and a number of function evaluations per-iteration $n\in\mathbb{N}$, the FD-DFD method is characterized by the iteration
\begin{equation}\label{DFD:eq:Ixk}
  x_{k+1}=x_k-\alpha g_k,
\end{equation}
where the gradient estimate
\begin{equation}\label{DFD:eq:gk}
  g_k=\frac{1}{n\sigma_k^2}\sum_{i=1}^n
  \bigg(f(\theta_{k,i})-\min_{1\leqslant j\leqslant n}
  f(\theta_{k,j})\bigg)(\theta_{k,i}-x_k),
\end{equation}
here, $\theta_{k,i}\sim\mathcal{N}(x_k,\sigma_k^2I_d)$ and $\sigma_k^2=\rho^{k}\lambda^{-1}$. To increase the stability of the iteration, we recommend to use the gradient estimate
\begin{equation}\label{DFD:eq:stablegk}
  \hat{g}_k=\frac{1}{n\hat{m}_k}\sum_{i=1}^n
  \bigg(f(\theta_{k,i})-\min_{1\leqslant j\leqslant n}
  f(\theta_{k,j})\bigg)(\theta_{k,i}-x_k),
\end{equation}
where
\begin{equation*}
  \hat{m}_k^2=\frac{1}{n}\sum_{i=1}^n \bigg(f(\theta_{k,i})-\min_{1\leqslant j\leqslant n}
  f(\theta_{k,j})\bigg)^2=\frac{1}{n}\sum_{i=1}^n
  \bigg(f(x_k+\sigma_k\xi_i)-\min_{1\leqslant j\leqslant n}
  f(\xi_j)\bigg)^2.
\end{equation*}
Here the random vector $\xi$ has $d$-dimensional standard normal distribution.

Note that, if $x_k\to x_*$ as $k\to\infty$, $\hat{m}_k^2$ can be viewed as an estimate for
\begin{equation*}
  m_k^2=\mathbb{E}_\xi\big(f(x_*+\sigma_k\xi)-f_*\big)^2.
\end{equation*}
Under Assumption \ref{DFD:ass:A}, i.e.,
\begin{equation*}
  \frac{l\sigma_k^2}{2}\|\xi\|_2^2\leqslant f(x_*+\sigma_k\xi)-f_*\leqslant
  \frac{L\sigma_k^2}{2}\|\xi\|_2^2,
\end{equation*}
we can obtain $m_k=\mathcal{O}(\sigma_k^2)$; and futher, $\hat{g}_k=\mathcal{O}(g_k)$. The difference between these two iterations based on $g_k$ and $\hat{g}_k$ is shown in Fig. \ref{DFD:fig:B}. We can see that when the stepsize parameter is properly selected, they have almost the same convergence behavior in the early stage while the estimate $\hat{g}_k$ is more stable than the estimate $g_k$ as the iteration progresses. Hence, we mainly consider the gradient estimate $\hat{g}_k$ in the following. Now we present the FD-DFD method as Algorithm \ref{DFD:alg:DFD}.

\begin{figure}[!h]
\centering
\includegraphics[width=0.325\textwidth]{fig2a}
\includegraphics[width=0.325\textwidth]{fig2b}
\includegraphics[width=0.325\textwidth]{fig2c}
\caption{Difference between two iterations. Left: the objective $f(x)=\|x\|_2^2-\frac{1}{2}\sum_{i=1}^2 \cos\big(5\pi x^{(i)}\big)+1,~ x\in\mathbb{R}^2$. Middle: the iteration $x_{k+1}=x_k-\alpha g_k$ with parameter settings $x_1=(1,-1)$, $\lambda=1/\sqrt{2}$, $\rho=0.9$, $n=5$ and $\alpha=0.5$. Right: the iteration $x_{k+1}=x_k-\alpha\hat{g}_k$ with parameter settings $x_1=(1,-1)$, $\lambda=1/\sqrt{2}$, $\rho=0.9$, $n=5$ and $\alpha=0.5$.}
\label{DFD:fig:B}
\end{figure}

\begin{algorithm}
\caption{FD-DFD Method}
\label{DFD:alg:DFD}
\begin{algorithmic}[1]
\STATE{Choose an initial iterate $x_1$ and preset parameters $\alpha>0$, $\lambda>0$, $\rho\in(0,1)$, $n\in\mathbb{N}$.}
\FOR{$k=1,2,\cdots$}
\STATE{Set the variance $\sigma_k=\rho^{\frac{k}{2}} \lambda^{-\frac{1}{2}}$.}
\STATE{Generate $n$ realizations $\{\theta_{k,i}\}_i^n$ of $\theta_k$ from $\mathcal{N}(x_k,\sigma_k^2I_d)$.}
\STATE{Compute a stochastic vector $\hat{g}_k$.}
\STATE{Set the new iterate as $x_{k+1}=x_k-\alpha\hat{g}_k$.}
\ENDFOR
\end{algorithmic}
\end{algorithm}

It follows from Chebyshev's inequality and $\mathbb{E}[g_k]=\nabla f_{\sigma_k}(x_k)$ that, for all $C>0$ and any $1\leqslant i\leqslant d$, with probability at least $1-\frac{1}{C^2}$, the $i$-th component of $g_k$ satisfies
\begin{equation}\label{DFD:eq:gki}
  \bigg|-g_k^{(i)}+\nabla^{(i)}f_{\sigma_k}(x_k)\bigg|\leqslant
  C\sqrt{\mathbb{V}\big[g_k^{(i)}\big]}.
\end{equation}
Furthermore, the $i$-th component of $x_{k+1}$ satisfies
\begin{equation}\label{DFD:eq:xki}
  x_{k+1}^{(i)}=x_k^{(i)}-\alpha g_k^{(i)}
  =x_k^{(i)}-\alpha\nabla^{(i)}f_{\sigma_k}(x_k).
\end{equation}

\subsection{Analyses}

To establish the expected linear convergence of the FD-DFD method, we first build upper bounds for $-\nabla^{(i)}f_{\sigma_k}(x_k)$ and $\mathbb{V}\big[g_k^{(i)}\big]$ by the following two lemmas.

\begin{lemma}\label{DFD:lem:gkbound}
Under Assumption \ref{DFD:ass:A}, suppose that $\sigma_k=\rho^{\frac{k}{2}}\lambda^{-\frac{1}{2}}$ with $0<\rho<1$ and there is an $M>0$ such that $\|x_k-x_*\|_2^2\leqslant\rho^kM$. Then for any $1\leqslant i\leqslant d$, the $i$-th component of $\nabla f_{\sigma_k}(x_k)$ satisfies the following inequality:
\begin{equation*}
  \bigg|-\nabla^{(i)}f_{\sigma_k}(x_k)
  +\frac{L+l}{2}\big(x_k^{(i)}-x_*^{(i)}\big)\bigg|
  \leqslant\frac{\rho^{\frac{k}{2}}(L-l)}{\sqrt{2\pi}}
  \left((d+2)\lambda^{-\frac{1}{2}}
  +\frac{2M}{\lambda^{-\frac{1}{2}}}\right);
\end{equation*}
further, when $\lambda^{-1}=\frac{2M}{d+2}$, this inequality can be improved as
\begin{equation*}
  \bigg|-\nabla^{(i)}f_{\sigma_k}(x_k)
  +\frac{L+l}{2}\big(x_k^{(i)}-x_*^{(i)}\big)\bigg|
  \leqslant\rho^{\frac{k}{2}}(L-l)\frac{\sqrt{(d+2)M}}{\sqrt{\pi}},
\end{equation*}
where $x^{(i)}$ be the $i$-th component of $x\in\mathbb{R}^d$.
\end{lemma}
\begin{proof}
For convenience we define $\varphi_k(x)=\prod_{i=1}^d\phi_k(x^{(i)})$, where
\begin{equation*}
  \phi_k(x^{(i)})=\frac{1}{\sqrt{2\pi}\sigma_k}
  \exp\bigg(\frac{-(x^{(i)}-x_k^{(i)})^2}{2\sigma_k^2}\bigg).
\end{equation*}
According to Theorem \ref{SGFD:thm:stepsmoothing}, for any $1\leqslant i\leqslant d$, we have
\begin{equation*}
  \nabla^{(i)}f_{\sigma_k}(x_k)
  =\frac{1}{\sigma_k^2}\int_{\mathbb{R}^d}f(x)
  \big(x^{(i)}-x_k^{(i)}\big)\varphi_k(x)\ud x.
\end{equation*}
Noting that
\begin{equation*}
  \int_{\mathbb{R}^d}\big(x^{(i)}-x_k^{(i)}\big)\varphi_k(x)\ud x=0,
\end{equation*}
it follows that
\begin{align*}
  \nabla^{(i)}f_{\sigma_k}(x_k)
  =&\frac{1}{\sigma_k^2}\int_{\mathbb{R}^d}\big(f(x)-f_*\big)
  \big(x^{(i)}-x_k^{(i)}\big)\varphi_k(x)\ud x \\
  =&\frac{1}{\sigma_k^2}\int_{x_k^{(i)}}^\infty\int_{\mathbb{R}^{d-1}}
  \big(f(x)-f_*\big)\big(x^{(i)}-x_k^{(i)}\big)
  \varphi_k(x)\ud x^{(-i)}\ud x^{(i)} \\
  &-\frac{1}{\sigma_k^2}\int_{-\infty}^{x_k^{(i)}}\int_{\mathbb{R}^{d-1}}
  \big(f(x)-f_*\big)\big(x_k^{(i)}-x^{(i)}\big)
  \varphi_k(x)\ud x^{(-i)}\ud x^{(i)},
\end{align*}
where $\ud x^{(-i)}=\ud x/\ud x^{(i)}$; together with Assumption \ref{DFD:ass:A}, i.e.,
\begin{equation*}
  \frac{l}{2}\|x-x_*\|_2^2\leqslant
  f(x)-f_*\leqslant\frac{L}{2}\|x-x_*\|_2^2,
\end{equation*}
this yields
\begin{align*}
  \nabla^{(i)}f_{\sigma_k}(x_k)
  \geqslant&\frac{l}{2\sigma_k^2}\int_{x_k^{(i)}}^\infty
  \int_{\mathbb{R}^{d-1}}\|x-x_*\|_2^2\big(x^{(i)}-x_k^{(i)}\big)
  \varphi_k(x)\ud x^{(-i)}\ud x^{(i)} \\
  &-\frac{L}{2\sigma_k^2}\int_{-\infty}^{x_k^{(i)}}
  \int_{\mathbb{R}^{d-1}}\|x-x_*\|_2^2
  \big(x_k^{(i)}-x^{(i)}\big)\varphi_k(x)\ud x^{(-i)}\ud x^{(i)} \\
  =&\frac{l}{2\sigma_k^2}\int_{\mathbb{R}^d}\|x-x_*\|_2^2
  \big(x^{(i)}-x_k^{(i)}\big)\varphi_k(x)\ud x \\
  &-\frac{L-l}{2\sigma_k^2}\int_{-\infty}^{x_k^{(i)}}
  \int_{\mathbb{R}^{d-1}}\|x-x_*\|_2^2
  \big(x_k^{(i)}-x^{(i)}\big)\varphi_k(x)\ud x^{(-i)}\ud x^{(i)}.
\end{align*}
Similarly, since
\begin{equation*}
  \int_{\mathbb{R}^d}\big(x^{(i)}-x_k^{(i)}\big)\varphi_k(x)\ud x=
  \int_{\mathbb{R}^d}\big(x^{(i)}-x_k^{(i)}\big)^3\varphi_k(x)\ud x=0,
\end{equation*}
it holds that
\begin{align*}
  &\frac{l}{2\sigma_k^2}\int_{\mathbb{R}^d}\|x-x_*\|_2^2 \big(x^{(i)}-x_k^{(i)}\big)\varphi_k(x)\ud x \\
  =&\frac{l}{2\sigma_k^2}\int_{\mathbb{R}^d}\big(x^{(i)}-x_*^{(i)}\big)^2
  \big(x^{(i)}-x_k^{(i)}\big)\varphi_k(x)\ud x \\
  =&\frac{l}{2\sigma_k^2}\!\int_\mathbb{R}\!\big[
  \big(x^{(i)}\!\!-\!x_k^{(i)}\big)^3\!\!+\!2
  \big(x_k^{(i)}\!\!-\!x_*^{(i)}\big)
  \big(x^{(i)}\!\!-\!x_k^{(i)}\big)^2\!\!+\!
  \big(x_k^{(i)}\!\!-\!x_*^{(i)}\big)^2\big(x^{(i)}\!\!-\!x_k^{(i)}\big)
  \big]\phi_k(x^{(i)})\ud x^{(i)} \\
  =&l\big(x_k^{(i)}-x_*^{(i)}\big).
\end{align*}
Hence, we obtain
\begin{equation*}
  -\nabla^{(i)}f_{\sigma_k}(x_k)\leqslant
  -l\big(x_k^{(i)}\!-\!x_*^{(i)}\big)+
  \frac{L-l}{2\sigma_k^2}\int_{-\infty}^{x_k^{(i)}}\!
  \int_{\mathbb{R}^{d-1}}\!\|x\!-\!x_*\|_2^2
  \big(x_k^{(i)}\!-\!x^{(i)}\big)\varphi_k(x)\ud x^{(-i)}\ud x^{(i)}.
\end{equation*}
Further, noting that
\begin{align*}
  &\int_\mathbb{R}\big(x^{(j)}-x_*^{(j)}\big)^2
  \phi_k(x^{(j)})\ud x^{(j)} \\
  =&\int_\mathbb{R}\big(x^{(j)}-x_k^{(j)}+x_k^{(j)}-x_*^{(j)}\big)^2
  \phi_k(x^{(j)})\ud x^{(j)} \\
  =&\int_\mathbb{R}\big(x^{(j)}-x_k^{(j)}\big)^2
  \phi_k(x^{(j)})\ud x^{(j)}
  +2\big(x_k^{(j)}-x_*^{(j)}\big)\int_\mathbb{R}
  \big(x^{(j)}-x_k^{(j)}\big)\phi_k(x^{(j)})\ud x^{(j)} \\
  &+\big(x_k^{(j)}-x_*^{(j)}\big)^2\int_\mathbb{R}
  \phi_k(x^{(j)})\ud x^{(j)}\\
  =&\sigma_k^2+\big(x_k^{(j)}-x_*^{(j)}\big)^2,
\end{align*}
it follows that
\begin{align*}
  &\frac{1}{\sigma_k^2}\sum_{j\neq i}\int_{-\infty}^{x_k^{(i)}}
  \int_{\mathbb{R}^{d-1}}\big(x^{(j)}-x_*^{(j)}\big)^2
  \big(x_k^{(i)}-x^{(i)}\big)\varphi_k(x)\ud x^{(-i)}\ud x^{(i)} \\
  =&\frac{1}{\sigma_k^2}\int_{-\infty}^{x_k^{(i)}}
  \big(x_k^{(i)}-x^{(i)}\big)\phi_k(x^{(i)})\ud x^{(i)}
  \cdot\sum_{j\neq i}
  \left[\sigma_k^2+\big(x_k^{(j)}-x_*^{(j)}\big)^2\right] \\
  \leqslant&\frac{1}{\sqrt{2\pi}\sigma_k}
  \left(\sigma_k^2d+\|x_k-x_*\|_2^2\right).
\end{align*}
Hence, it follows from $\sigma_k=\rho^{\frac{k}{2}}\lambda^{-\frac{1}{2}}$ and $\|x_k-x_*\|_2^2\leqslant\rho^kM$ that
\begin{equation*}
  \frac{1}{\sigma_k^2}\sum_{j\neq i}\int_{-\infty}^{x_k^{(i)}}
  \int_{\mathbb{R}^{d-1}}\big(x^{(j)}-x_*^{(j)}\big)^2
  \big(x_k^{(i)}-x^{(i)}\big)\varphi_k(x)\ud x^{(-i)}\ud x^{(i)}
  \leqslant\frac{\rho^{\frac{k}{2}}}{\sqrt{2\pi}}
  \left(\lambda^{-\frac{1}{2}}d+\frac{M}{\lambda^{-\frac{1}{2}}}\right).
\end{equation*}
Similarly, we have
\begin{align*}
  &\frac{1}{\sigma_k^2}\int_{-\infty}^{x_k^{(i)}}
  \int_{\mathbb{R}^{d-1}}\big(x^{(i)}-x_*^{(i)}\big)^2
  \big(x_k^{(i)}-x^{(i)}\big)\varphi_k(x)\ud x^{(-i)}\ud x^{(i)} \\
  =&\frac{1}{\sigma_k^2}\int_{-\infty}^{x_k^{(i)}}
  \big(x^{(i)}-x_*^{(i)}\big)^2\big(x_k^{(i)}-x^{(i)}\big)
  \phi_k(x^{(i)})\ud x^{(i)} \\
  =&\frac{1}{\sigma_k^2}\int_{-\infty}^{x_k^{(i)}}
  \big(x^{(i)}-x_k^{(i)}+x_k^{(i)}-x_*^{(i)}\big)^2
  \big(x_k^{(i)}-x^{(i)}\big)\phi_k(x^{(i)})\ud x^{(i)} \\
  =&\frac{1}{\sigma_k^2}\int_{-\infty}^{x_k^{(i)}}
  \big(x_k^{(i)}-x^{(i)}\big)^3\phi_k(x^{(i)})\ud x^{(i)}
  -2\frac{x_k^{(i)}-x_*^{(i)}}{\sigma_k^2}\int_{-\infty}^{x_k^{(i)}}
  \big(x_k^{(i)}-x^{(i)}\big)^2\phi_k(x^{(i)})\ud x^{(i)} \\
  &+\frac{\big(x_k^{(i)}-x_*^{(i)}\big)^2}{\sigma_k^2}
  \int_{-\infty}^{x_k^{(i)}}
  \big(x_k^{(i)}-x^{(i)}\big)\phi_k(x^{(i)})\ud x^{(i)} \\
  =&\frac{2\sigma_k}{\sqrt{2\pi}}-\big(x_k^{(i)}-x_*^{(i)}\big)
  +\frac{1}{\sqrt{2\pi}\sigma_k}\big(x_k^{(i)}-x_*^{(i)}\big)^2 \\
  \leqslant&-\big(x_k^{(i)}-x_*^{(i)}\big)
  +\frac{\rho^{\frac{k}{2}}}{\sqrt{2\pi}}
  \left(2\lambda^{-\frac{1}{2}}+\frac{M}{\lambda^{-\frac{1}{2}}}\right).
\end{align*}
Finally, we obtain
\begin{align*}
  -\nabla^{(i)}f_{\sigma_k}(x_k)
  \leqslant&-\frac{L+l}{2}\big(x_k^{(i)}-x_*^{(i)}\big)
  +\frac{\rho^{\frac{k}{2}}(L-l)}{\sqrt{2\pi}}
  \left((d+2)\lambda^{-\frac{1}{2}}
  +\frac{2M}{\lambda^{-\frac{1}{2}}}\right) \\
  \leqslant&-\frac{L+l}{2}\big(x_k^{(i)}-x_*^{(i)}\big)
  +\rho^{\frac{k}{2}}(L-l)\frac{\sqrt{(d+2)M}}{\sqrt{\pi}},
\end{align*}
where $\lambda^{-1}=\frac{2M}{d+2}$.

Similarly, we ca prove that
\begin{align*}
  -\nabla^{(i)}f_{\sigma_k}(x_k)
  \geqslant-\frac{L+l}{2}\big(x_k^{(i)}-x_*^{(i)}\big)
  -\rho^{\frac{k}{2}}(L-l)\frac{\sqrt{(d+2)M}}{\sqrt{\pi}},
\end{align*}
and the proof is complete.\qed
\end{proof}

\begin{lemma}\label{DFD:lem:Vkbound}
Under Assumption \ref{DFD:ass:A}, suppose that $\sigma_k=\rho^{\frac{k}{2}}\lambda^{-\frac{1}{2}}$ with $0<\rho<1$ and there is an $M>0$ such that $\|x_k-x_*\|_2^2\leqslant\rho^kM$. Then for any $1\leqslant i\leqslant d$, the variance of $g_k^{(i)}$ satisfies the following inequality:
\begin{equation*}
  \mathbb{V}\big[g_k^{(i)}\big]\leqslant
  \rho^k\frac{L^2}{n}\left(\frac{d+2}{\lambda}+M\right),
\end{equation*}
where $g_k^{(i)}$ is the $i$-th component of $g_k\in\mathbb{R}^d$.
\end{lemma}
\begin{proof}
According to Assumption \ref{DFD:ass:A}, for any $1\leqslant s\leqslant n$, we have
\begin{equation*}
  f(\theta_{k,s})-\min_{1\leqslant j\leqslant n}f(\theta_{k,j})
  \leqslant f(\theta_{k,s})-f_*\leqslant
  \frac{L}{2}\|\theta_{k,s}-x_*\|_2\leqslant
  \frac{L\sigma_k^2}{2}\|\xi_s\|_2+\frac{L}{2}\|x_k-x_*\|_2,
\end{equation*}
then for any $1\leqslant i\leqslant d$, together with the definition of $g_\sigma(x,\xi)$, we further obtain
\begin{align*}
  \mathbb{V}_\xi\big[g_\sigma^{(i)}(x,\xi)\big]
  \leqslant&\mathbb{E}_\xi\big[g_\sigma^{(i)}(x,\xi)\big]^2\\
  =&\mathbb{E}_\xi\bigg[\frac{\big(f(\theta_k)-\min f(\theta_k)\big)
  \xi^{(i)}}{\sigma_k}\bigg]^2 \\
  \leqslant&\mathbb{E}_\xi\bigg[\frac{L\sigma_k}{2}\|\xi_s\|_2\xi^{(i)}
  +\frac{L\|x_k-x_*\|_2}{2\sigma_k}\xi^{(i)}\bigg]^2 \\
  \leqslant&L^2\sigma_k^2\mathbb{E}_\xi
  \big[\|\xi\|_2\xi^{(i)}\big]^2
  +\frac{L^2\|x_k-x_*\|_2^2}{\sigma_k^2}
  \mathbb{E}_\xi\big[\xi^{(i)}\big]^2 \\
  \leqslant&L^2\rho^k\lambda^{-1}\mathbb{E}_\xi
  \big[\|\xi\|_2\xi^{(i)}\big]^2+L^2\rho^kM.
\end{align*}
Since
\begin{equation*}
  \mathbb{E}_\xi\big[\|\xi\|_2\xi^{(i)}\big]^2
  =\sum_{j=1}^d\frac{1}{(\sqrt{2\pi})^d}
  \int_{\mathbb{R}^d}\big(\xi^{(j)}\xi^{(i)}\big)^2
  e^{-\frac{\|\xi\|_2^2}{2}}\ud\xi=d+2,
\end{equation*}
it follows that
\begin{equation*}
  \mathbb{V}\big[g_k^{(i)}\big]\leqslant
  \rho^k\frac{L^2}{n}\left(\frac{d+2}{\lambda}+M\right),
\end{equation*}
and the proof is complete.\qed
\end{proof}

\begin{theorem}\label{DFD:thm:main1}
Under Assumption \ref{DFD:ass:A}, suppose that $\sigma_k=\rho^{\frac{k}{2}}\lambda^{-\frac{1}{2}}$ with $0<\rho<1$ and there is an $M>0$ such that $\|x_s-x_*\|_2^2\leqslant\rho^sM$ for all $1\leqslant s\leqslant k$. If the FD-DFD method (Algorithm \ref{DFD:alg:DFD}) is run with a stepsize parameter $\alpha>0$ such that
\begin{equation*}
  \rho_\alpha=2\left[1-\frac{\alpha(L+l)}{2}\right]^2<\rho<1,
\end{equation*}
then with probability at least $1-\frac{1}{C^2}$, the iterates of FD-DFD satisfy:
\begin{equation*}
  \big\|x_{k+1}-x_*\big\|_2^2\leqslant
  \rho_\alpha^k\big\|x_1\!-\!x_*\big\|_2^2
  +\rho^k\frac{2\rho_\alpha\alpha^2K_C^2}{\rho-\rho_\alpha},
\end{equation*}
where the constant
\begin{equation*}
  K_C=(L-l)\sqrt{\frac{(d+2)M}{\pi}}
  +L\frac{C}{\sqrt{n}}\sqrt{\frac{d+2}{\lambda}+M}
\end{equation*}
is independent of $k$ but depends on $C$.
\end{theorem}
\begin{proof}
According to \eqref{DFD:eq:gki}, Lemmas \ref{DFD:lem:gkbound} and \ref{DFD:lem:Vkbound}, it holds that for any $1\leqslant i\leqslant d$ and all $C>0$, with probability at least $1-\frac{1}{C^2}$,
\begin{equation*}
  \bigg|-g_k^{(i)}+\frac{L+l}{2}\big(x_k^{(i)}-x_*^{(i)}\big)\bigg|
  \leqslant\rho^{\frac{k}{2}}K_C,
\end{equation*}
Together with \eqref{DFD:eq:xki}, we obtain
\begin{align*}
  \Big|x_{k+1}^{(i)}-x_*^{(i)}\Big|
  =&~\Big|x_k^{(i)}-x_*^{(i)}-\alpha g_k^{(i)}\Big| \\
  =&~\Big|x_k^{(i)}-x_*^{(i)}-\frac{\alpha(L+l)}{2}
  \big(x_k^{(i)}-x_*^{(i)}\big)+\Big(-\alpha g_k^{(i)}
  +\frac{\alpha(L+l)}{2}\big(x_k^{(i)}-x_*^{(i)}\big)\Big)\Big| \\
  \leqslant&~\Big|x_k^{(i)}-x_*^{(i)}-\frac{\alpha(L+l)}{2}
  \big(x_k^{(i)}-x_*^{(i)}\big)\Big|+\rho^{\frac{k}{2}}\alpha K_C \\
  =&\left|1-\frac{\alpha(L+l)}{2}\right|\Big|x_k^{(i)}-x_*^{(i)}\Big|
  +\rho^{\frac{k}{2}}\alpha K_C,
\end{align*}
further, by the Arithmetic Mean Geometric Mean inequality, we have
\begin{align*}
  \big(x_{k+1}^{(i)}-x_*^{(i)}\big)^2
  \leqslant&2\left[1-\frac{\alpha(L+l)}{2}\right]^2
  \big(x_k^{(i)}-x_*^{(i)}\big)^2+\rho^k\cdot2\alpha^2K_C^2 \\
  =&\rho_\alpha\big(x_k^{(i)}-x_*^{(i)}\big)^2
  +\rho^k\cdot2\alpha^2K_C^2.
\end{align*}
Finally, by summing $i$ from $1$ to $d$, one obtains
\begin{equation*}
  \big\|x_{k+1}-x_*\big\|_2^2\leqslant\rho_\alpha
  \big\|x_k-x_*\big\|_2^2+\rho^k\cdot2\alpha^2K_C^2,
\end{equation*}
and doing it recursively, one further obtains,
\begin{align*}
  \big\|x_{k+1}\!-\!x_*\big\|_2^2\leqslant&
  \rho_\alpha^k\big\|x_1\!-\!x_*\big\|_2^2
  +\rho^k\left(1\!+\!\frac{\rho_\alpha}{\rho}\!+\!\cdots\!+\!
  \frac{\rho_\alpha^{k-1}}{\rho^{k-1}}\right)2\alpha^2K_C^2 \\
  \leqslant&\rho_\alpha^k\big\|x_1\!-\!x_*\big\|_2^2
  +\rho^k\frac{2\rho_\alpha\alpha^2K_C^2}{\rho-\rho_\alpha},
\end{align*}
and the proof is complete.\qed
\end{proof}

The following theorem states that when the parameters are properly selected,  the iterates of FD-DFD satisfy $\|x_{k+1}-x_*\|_2^2\leqslant\rho^k\big\|x_1-x_*\big\|_2^2$ for all $k\in\mathbb{N}$ in probability.
\begin{theorem}\label{DFD:thm:main2}
Under Assumption \ref{DFD:ass:A}, suppose that $\sigma_k=\rho^{\frac{k}{2}}\lambda^{-\frac{1}{2}}$ with $0<\rho<1$. If the FD-DFD method (Algorithm \ref{DFD:alg:DFD}) is run with a stepsize parameter $\alpha>0$ such that
\begin{equation*}
  \frac{\rho_\alpha}{\rho}\left[1+
  \frac{2\alpha^2K_C^2}{\rho-\rho_\alpha}\right]<1,
\end{equation*}
then with high probability, the iterates of FD-DFD satisfy for all $k\in\mathbb{N}$:
\begin{equation*}
  \|x_{k+1}-x_*\|_2^2\leqslant\rho^k\big\|x_1-x_*\big\|_2^2,
\end{equation*}
where $\rho_\alpha=2\left[1-\frac{\alpha(L-l)}{2}\right]^2$ and
\begin{equation*}
  K_C=(L-l)\sqrt{\frac{(d+2)\|x_1-x_*\|_2^2}{\pi\rho}}
  +L\frac{C}{\sqrt{n}}\sqrt{\frac{d+2}{\lambda}
  +\frac{\|x_1-x_*\|_2^2}{\rho}}.
\end{equation*}
\end{theorem}
\begin{proof}
Let $M=\rho^{-1}\|x_1-x_*\|_2^2$, that is, $\|x_1-x_*\|_2^2\leqslant\rho M$, which satisfies the condition of Theorem \ref{DFD:thm:main1}. Using Theorem \ref{DFD:thm:main1} and induction, one can deduce
\begin{equation*}
  \big\|x_{k+1}-x_*\big\|_2^2\leqslant
  \rho^{k+1}M\left(\frac{\rho_\alpha^k}{\rho^k}
  +\frac{2\rho_\alpha\alpha^2K_C^2}{\rho(\rho-\rho_\alpha)}\right).
\end{equation*}
Notice that for every $k\in\mathbb{N}$, it follows that
\begin{equation*}
  \frac{\rho_\alpha^k}{\rho^k}
  +\frac{2\rho_\alpha\alpha^2K_C^2}{\rho(\rho-\rho_\alpha)}
  <\frac{\rho_\alpha}{\rho}\left[1+
  \frac{2\alpha^2K_C^2}{\rho-\rho_\alpha}\right]<1,
\end{equation*}
thus, one can finally obtain
\begin{equation*}
  \big\|x_{k+1}-x_*\big\|_2^2\leqslant
  \rho^{k+1}M=\rho^k\big\|x_1-x_*\big\|_2^2,
\end{equation*}
and the proof is complete.\qed
\end{proof}

Since $n$ is independent of $k$, the following total work complexity bound for the FD-DFD method is immediate from Theorem \ref{DFD:thm:main2}.
\begin{corollary}[Complexity bound]
Suppose the conditions of Theorem \ref{DFD:thm:main2} hold. Then the
number of function evaluations of the FD-DFD (Algorithm \ref{DFD:alg:DFD}) required to achieve $\|x_k-x_*\|_2^2\leqslant\epsilon$ is $\mathcal{O}(\log(1/\epsilon))$.
\end{corollary}

\section{Comparison of RAD and FD-DFD}
\label{DFD:s3}

With an initial point $x_1$, three fixed parameters $\lambda>0$, $0<\rho<1$ and $n\in\mathbb{N}$, the RAD method \cite{LuoX2020A_RAD} is characterized by the iteration
\begin{equation}\label{DFD:eq:RAD}
  x_{k+1}=\frac{\sum_{i=1}^n\theta_{k,i} \exp[-m_k^{-1}(f(\theta_{k,i})-f_*)]}{\sum_{i=1}^n \exp[-m_k^{-1}(f(\theta_{k,i})-f_*)]}.
\end{equation}
With an additional stepsize parameter $\alpha>0$, the FD-DFD method is characterized by the iteration
\begin{equation}\label{DFD:eq:DFD}
  x_{k+1}=x_k-\frac{\alpha}{nm_k}\sum_{i=1}^n
  \bigg(f(\theta_{k,i})-\min_{1\leqslant j\leqslant n}
  f(\theta_{k,i})\bigg)(\theta_{k,i}-x_k),
\end{equation}
where $\theta_{k,i}\sim\mathcal{N}(x_k,\rho^{k}\lambda^{-1}I_d)$ and $m_k^2=\mathbb{E}[(f(\theta_{k,i})-f_*)^2]$. In practice, $f_*$ and $m_k$ should be replaced with corresponding estimates.

Both the RAD and FD-DFD method have their own characteristics. Since RAD is based on the asymptotic representation for the solution of regularized minimization, it does not require a stepsize parameter, or in other words, it can automatically obtain the optimal stepsize; however, the parameter $n$ in RAD algorithms increases as the dimension $d$ increases (See Fig. 3 in \cite{LuoX2020A_RAD}). In comparison, the parameter $n$ in FD-DFD algorithms is almost independent of $d$ because the variance term containing $n$ is insignificance for an appropriate $n$ (See Fig. \ref{DFD:fig:E2} and Theorem \ref{DFD:thm:main1}); but as a price, there is an additional stepsize parameter $\alpha$ and the choice of $\alpha$ directly affects whether the iterate sequence will converge to the global minimizer.

\section{Numerical experiments}
\label{DFD:s4}

We illustrate the performance of the FD-DFD algorithm by considering the revised Rastrigin function in $\mathbb{R}^d$ defined as
\begin{equation*}
  f(x)=\|x\|_2^2-\frac{1}{2}\sum_{i=1}^d\cos\big(5\pi x^{(i)}\big)
  +\frac{d}{2},~~\textrm{where $x^{(i)}$ be the $i$-th component of $x$}.
\end{equation*}
As shown in Fig. \ref{DFD:fig:A}, this revised Rastrigin function satisfies Assumption \ref{DFD:ass:A}. It has a unique global minima located at the origin and many local minima, e.g., the number of its local minima reaches $5^d$ in the hypercube $[-1,1]^d$.

Fig. \ref{DFD:fig:E2} shows the performance of FD-DFD algorithms for the revised Rastrigin function in various dimensions from $5$ to $500$. These experiments clearly demonstrate the global linear convergence of the FD-DFD method.

\begin{figure}[tbhp]
\centering
\subfigure{\includegraphics[width=0.325\textwidth]{fig3a}}
\subfigure{\includegraphics[width=0.325\textwidth]{fig3b}}
\subfigure{\includegraphics[width=0.325\textwidth]{fig3c}}
\subfigure{\includegraphics[width=0.325\textwidth]{fig3d}}
\subfigure{\includegraphics[width=0.325\textwidth]{fig3e}}
\subfigure{\includegraphics[width=0.325\textwidth]{fig3f}}
\subfigure{\includegraphics[width=0.325\textwidth]{fig3g}}
\subfigure{\includegraphics[width=0.325\textwidth]{fig3h}}
\subfigure{\includegraphics[width=0.325\textwidth]{fig3i}}
\caption{Performance of the FD-DFD method for the revised Rastrigin function in various dimensions, every initial iterate is randomly selected on a sphere of radius $\sqrt{d}$ centered at the origin, the parameter $\lambda=1/\sqrt{d}$, three different settings for the parameters $\rho$ and $\alpha$ are run independently for each plot.}
\label{DFD:fig:E2}
\end{figure}

In these experiments, every initial iterate is randomly selected on a $d$-dimensional sphere of radius $\sqrt{d}$ centered at the origin. Furthermore, the random vectors in each iteration are sequentially generated by a halton sequence with RR$2$ scramble type \cite{KocisL1997A_QMCscramble}. The algorithm is implemented in Matlab. The source code of the implementation is available at https://github.com/xiaopengluo/dfd.

\section{Conclusions}
\label{DFD:s5}

In this work we have analyzed that the finite-difference derivative-free descent (FD-DFD) method enjoys linear convergence for finding the global minima of a class of multiple minima functions. It also has a total work complexity bound $\mathcal{O}(\log\frac{1}{\epsilon})$ to find a point such that the gap between this point and the global minimizer is less than $\epsilon$. Numerical experiments in various dimensions demonstrate all the benefits.

\begin{acknowledgements}
We thank Prof. Herschel A. Rabitz for several discussions about global optimization for multiple minima problems.
\end{acknowledgements}

\bibliographystyle{spmpsci}
\bibliography{MReferences}

\end{document}


text:
\section{Implementations and proofs}\label{sec:algorithm}

In this section, we describe our implementations of the algorithms
in the previous sections and our proofs of the main results stated in the Introduction. As we have seen, there are two approaches to enumerating
superspecial Howe curves: (A) $(E_1,E_2)$-first
and (B) $C$-first.
The arguments in the previous sections show that (B) has an advantage in the complexity analysis.
Here we see that (B) is much superior to (A) also when we execute their implementations. Indeed, Theorems \ref{thm:main1} and \ref{thm:main2} in
the Introduction were obtained by Magma implementations based on (B) that were run on a PC with ubuntu 16.04 LTS OS at 3.40GHz CPU (Intel Core i7-6700) and 15.6 GB memory. The same result for $p\le 53$ was obtained by implementing the method (A) over Magma with an execution by the same PC. Although
it took 11871 seconds to obtain Theorem \ref{thm:main2} for $p\le 53$ by (A), the second strategy (B) finishes the enumeration for $p \le 199$
only in 924 seconds; see Table \ref{tableAB:2} for benchmark timing data for small $p$.
\begin{table}[ht]
\caption{Benchmark timing data for (A) Algorithm \ref{alg:main1}
and (B) Algorithm \ref{alg:main2}. All times shown are in seconds.}
\label{tableAB:2}
\vspace{-1ex}
\begin{center}
\renewcommand{\arraystretch}{0.9}
\begin{tabular}{rrrrrrrrrrr}
\toprule
  $p$ & \multicolumn{1}{c}{(A)} & \multicolumn{1}{c}{(B)} &\hbox to 1em{}& $p$ &  \multicolumn{1}{c}{(A)} &  \multicolumn{1}{c}{(B)} &\hbox to 1em{}& $p$ &  \multicolumn{1}{c}{(A)} &  \multicolumn{1}{c}{(B)}\\
\cmidrule{1-3}\cmidrule{5-7}\cmidrule{9-11}
  $5$ &  $0.02$ &  $0.08$ &&  $19$ &   $6.14$ &  $0.12$ && $41$ &  $1118.63$ &  $0.71$ \\
  $7$ &  $0.01$ &  $0.01$ &&  $23$ &  $27.59$ &  $0.21$ && $43$ &  $1423.26$ &  $0.80$ \\
 $11$ &  $0.17$ &  $0.04$ &&  $29$ & $114.70$ &  $0.31$ && $47$ &  $2686.17$ &  $1.03$ \\
 $13$ &  $0.76$ &  $0.05$ &&  $31$ & $193.82$ &  $0.34$ && $53$ &  $5678.32$ &  $1.46$ \\
 $17$ &  $3.92$ &  $0.09$ &&  $37$ & $617.23$ &  $0.54$ &&      &            &         \\
\bottomrule
\end{tabular}
\end{center}
\end{table}
The code for our implementations is available on the first author's web site
for (A) and on the third author's web site for (B).
In case (A), it is very costly to find Cartier--Manin matrices, and in addition to that there are many pairs $(E_1, E_2)$ of supersingular elliptic curves. This fact is consistent with the complexity analysis in Section \ref{subsec:alg1comp}.
On the other hand, the method (B) contains few intensive computations
and it enables us to find and enumerate superspecial Howe curves very efficiently.

The preceding remarks prove the computational results in Theorems~\ref{thm:main1} and~\ref{thm:main2}, and we are left to prove
the statement in Theorem~\ref{thm:main1} concerning primes $p\equiv 5\bmod 6.$
This fact is shown by using the Howe curve defined by $E_1\colon z^2y=x^3+y^3$ and $E_2\colon w^2y=x^3+ay^3$ with $a\in\{-1,1/4\}$.
Indeed, if $p\equiv 5 \bmod 6$, then
these two elliptic curves are supersingular and moreover
$y^2 = (x^3+1)(x^3+a)$ is superspecial.
This can be checked
by observing that the curve has two nonhyperelliptic involutions,
given by $(x,y)\mapsto (a^{1/3}/x, \pm a^{1/2}y/x^3)$, so that its Jacobian
is $(2,2)$-isogenous to a product of elliptic curves. For $a = -1$ we find that
these two curves are both isomorphic to the $j=0$ curve with CM by $-3$,
and for $a = 1/4$ we find that they are both isomorphic to the $j = -12288000$
curve with CM by $-27$. In both cases, these elliptic curves are supersingular
for primes $p\equiv 5\bmod 6$.

We remark that this curve for $a=1/4$ is isomorphic to the curve $X^3+Y^3+W^3=2YW+Z^2=0$ in $\bbP^3$ studied by the first author in \cite{K19}, by the correspondence $x=X$, $y=Y+W$, $z=\sqrt{-3/2}Z$ and $w=\sqrt{-3/4}(Y-W)$.

config: arxiv


TITLE: Distributions which are intractable to sample from?
QUESTION [5 upvotes]: I'm looking for an interesting family of probability distributions $P$ that is intractable to efficiently sample from.
I'm not sure what the right notion of intractable is, though I know the notion of efficiently samplable: the existence of a polynomial-time randomized Turing machine that can simulate samples from $P$.

Q1: Is a reasonable notion of intractable that if there exists such a randomized Turing machine, it can be derandomized to solve an NP-complete problem? If not, what is the right notion? Is there a good reference for this? (Or just some good search terms.)

There are a lot of tautological examples along the following lines: there is no algorithm to sample from [set of solutions to NP-complete problem union a symbol for no solution], because there is an obvious reduction to [that NP-complete problem] by taking the output of the sampler.

Q2: What are some more interesting examples of intractable sampling problems than these?

There are also examples perhaps like this: Given a set $X$, the goal is to sample uniformly from $P(P(X))$ (powerset). If you could do so with a polynomial time turing machine, you could produce exponentially many bits of entropy in only polynomially many coin flips, which is impossible (because $H(f(Y)) \leq H(Y)$ for any random variable $Y$ and deterministic function $f$).
In order to avoid this kind of example, I want the underlying set of the distribution to grow at most like $2^{p(|x|)}$, where $p$ is a polynomial in the size of the input.

REPLY [3 votes]: I'll expand my comment to an answer. Many combinatorial structures in graphs are actually NP-hard to sample from.
The earliest example I can think of is JVV86 (Thm 5.1), which shows that there is no polynomial-time algorithm to randomly generate a cycle in a directed graph, unless NP=RP.
Similar results along this line are NP-hardness for sampling independent sets, colourings etc. See e.g. DFJ02, Sly10, GSV15. Note that the corresponding search problems are all easy.
The problem can also be defined by "soft" constraints. E.g., sampling from an anti-ferromagnetic Ising model is NP-hard. See JS93 (Thm 14).
Results of this sort are usually stated as NP-hard approximate counting problems, due to the standard self-reductions JVV86.


TITLE: Does the concept of a positive null net make sense for all directed sets?
QUESTION [1 upvotes]: A common strategy to prove things is to use null sequences and I was wondering if it is possible to generalise the concept to aritrary nets.
Let $I$ be a directed set. Is there a net $\left( x_\alpha \right)_{\alpha \in I}$ of positive real numbers such that $lim_\alpha x_\alpha = 0$?
This is certainly equivalent to the existence of a monotonic and final map from $I$ to $\mathbb{R}$ and perhaps even to $\mathbb{N}$.
My gut feeling tells me that such maps do not exist for "too large" $I$, i.e that there are totally ordered sets $I$ for which such monotonic maps must eventually be constant.
But my intuition might be tainted because I would attempt to construct such a map by taking a strictly monotonic and cofinal map from $\mathbb{N}$ to $I$ (which most certainly does not always exist) and then to "fill the gaps" (which also does not work unless $I$ is totally ordered).
Thus, I could imagine that there are non-constructive ways to tackle this problem by starting with some arbitrary surjective map from $I$ to $\mathbb{N}$.

REPLY [2 votes]: Quite simply, let $I=(\omega_1,\le)$ be the first uncountable ordinal, let $(x_i)_{i\in I}$ be a null net and call $\alpha_n=\min\{i\in I\,:\, \forall j\ge i, x_j\le1/n\}$. Since $\operatorname{cof}\omega_1=\omega_1>\omega$, we have that $\beta=\sup_{n\in\Bbb N}\alpha_n<\omega_1$. By definition, $x_\gamma$ cannot be positive for any $\gamma\ge\beta$.


TITLE: Data transmission via laser illuminating the ceiling
QUESTION [0 upvotes]: I've heard a data transmission technique proposed more than once wherein a base station would fire a laser at the ceiling of an indoor space. The pulsing of this laser would be observed by one or more omnidirectional receivers (mounted on wireless devices such as phones) which would decode the signal, resulting in short-range high-bandwidth data transmission system, at least in one direction. Basically an optical fiber without the fiber. The intention would be to achieve a multi-gigabit connection that could stream video.
Does this work and if not, why? I can't find any discussion of it online and the few similar systems I've found invariably involve rotating the emitter and receiver to face each other, suggesting there is some problem with this approach.

REPLY [2 votes]: This is called optical wireless communication. It has been demonstrated many times in the lab, going back to at least the 1990's. There have been attempts at commercialization, for example under the name Li-Fi but none have been especially successful.


text:
\subsection{Volume-preserving involutive functions lead to Metropolis--Hastings proposals satisfying detailed balance}
\label{sec:validity-proof}

The proof of detailed balance for Hybrid Monte Carlo relies on the fact that Hamiltonian dynamics composed with negating momentum is involutive and volume-preserving. It is also known that volume-preserving involutive functions lead to Metropolis--Hastings proposals satisfying detailed balance (see, e.g., \citep{marco-proof}), but we were unable to find a published argument. For completeness we describe here a procedure for using any volume-preserving involutive function as a proposal, and prove its correctness.

Our goal is to use a volume-preserving involutive function $f_\theta$ to construct a Markov process with $P_S(\phi)$ as a stationary distribution. To do this, we will find a transition such that the transition probabilities $P_M(\phi \mapsto \phi')$ satisfy the detailed balance condition:

\begin{align*}
P_S(\phi) P_M(\phi \mapsto \phi') = P_S(\phi') P_M(\phi' \mapsto \phi)
\end{align*}

We follow the original derivation of Hybrid Monte Carlo \citep{duane1987hybrid}, since the structure of the proof is similar.

In order to make a transition, we do the following.
\begin{enumerate}
\item Introduce an auxiliary random variable $\pi$ with probability density $P_G$.
\item Propose a transition drawn from $P_H((\phi,\pi) \mapsto (\phi',\pi'))$ according to the volume-preserving involutive function $f_\theta$:

\begin{align*}
P_H((\phi,\pi) \mapsto (\phi',\pi')) = \delta[(\phi',\pi') - f_\theta((\phi,\pi))].
\end{align*}

\item Accept or reject that transition according to the Metropolis--Hastings acceptance criterion
\begin{align*}
P_A&((\phi,\pi) \mapsto (\phi',\pi')) \\
&= \min\Bigl(1,\frac{P_S(\phi')P_G(\pi')P_H((\phi',\pi') \mapsto (\phi,\pi))}{P_S(\phi)P_G(\pi)P_H((\phi,\pi) \mapsto (\phi',\pi'))}\Bigr).
\end{align*}
\item Marginalize over the auxiliary variable $\pi$.
\end{enumerate}

Formally, we define our transition probability by

\begin{align*}
P_M&(\phi \mapsto \phi') \\
&= \int \Bigl(
    P_G(\pi)\, P_H((\phi,\pi) \mapsto (\phi',\pi')) \\
& \qquad \qquad P_A((\phi,\pi) \mapsto (\phi',\pi')) \Bigr)
    d\pi d\pi'.
\end{align*}

We now show that this transition satisfies the detailed balance condition.

\begin{lemma}
\label{lem:involutive-unchanged}
Applying $f_\theta$ within a Dirac $\delta$ distribution leaves the $\delta$ unchanged:

\begin{align*}
\delta[x-y] \,=\, \delta[f_\theta(x)-f_\theta(y)]
\end{align*}
\end{lemma}

\begin{proof}
For arbitrary $F$ we have
    \[
        F(y) = \int_{\Omega} \: \delta[x-y]F(x) \, dx
    \]
and
\begin{align*}
    \int_{\Omega}\: \delta[f_\theta(x)-&f_\theta(y)]F(x) \,dx \\
=\, &\int_{f_\theta(\Omega)} \: \delta[u-f_\theta(y)]F(f_\theta^{-1}(u)) \, \frac{du}{|\det f_\theta'(x)|}\\
=\, &\int_{f_\theta(\Omega)}\: \delta[u-f_\theta(y)]F(f_\theta^{-1}(u)) \,du \\ & \hspace{2cm} \text{(since $f_\theta$ preserves volume)} \\
=\, &F(f_\theta^{-1}(f_\theta(y))) \\
=\, &F(y),
\end{align*}

and so

\[
        \int_{\Omega} \: \delta[x-y]F(x)\,dx \,=\, \int_{\Omega} \: \delta[f(x)-f(y)]F(x)\, dx.
\]

Hence $\delta[x-y]	= \delta[f(x)-f(y)]$.
\end{proof}

\begin{lemma}
\label{lem:PH-symmetric}
$P_H$ is symmetric:

\begin{align*}
P_H((\phi,\pi) \mapsto (\phi',\pi')) = P_H((\phi',\pi') \mapsto (\phi,\pi)).
\end{align*}
\end{lemma}

\begin{proof}
We have
\begin{align*}
    P_H((\phi,\pi&) \mapsto (\phi',\pi')) \\
&= \delta[(\phi',\pi') - f_\theta((\phi,\pi))] \\
&= \delta[f_\theta((\phi',\pi')) - f_\theta \circ f_\theta((\phi,\pi))] \\ & \hspace{2cm} (\text{by
Lemma~\ref{lem:involutive-unchanged}}) \\
&= \delta[f_\theta((\phi',\pi')) - (\phi,\pi)] \\ & \hspace{2cm} \text{(since $f_\theta$ is involutive)} \\
&= \delta[(\phi,\pi) - f_\theta((\phi',\pi'))] \\ & \hspace{2cm} (\text{since} \: \delta[x] = \delta[-x]) \\
&=P_H((\phi',\pi') \mapsto (\phi,\pi)),
\end{align*}

establishing the lemma.
\end{proof}

\begin{lemma}
$P_A$ satisfies the following simpler form:
\begin{align*}
P_A((\phi&,\pi) \mapsto (\phi',\pi')) = \min\Bigl(1,\frac{P_S(\phi')P_G(\pi')}{P_S(\phi)P_G(\pi)}\Bigr).
\end{align*}
\end{lemma}

\begin{proof}
Apply Lemma~\ref{lem:PH-symmetric} to the definition of $P_A$.
\end{proof}

\begin{lemma}
\label{lem:PSPGPA-symmetric}
$P_SP_GP_A$ is symmetric:
\begin{align*}
P_S(\phi)P_G(\pi)P_A&((\phi,\pi) \mapsto (\phi',\pi')) \\
= \: &P_S(\phi')P_G(\pi')P_A((\phi',\pi') \mapsto (\phi,\pi)).
\end{align*}
\end{lemma}

\begin{proof}
We have
\begin{align*}
P_S(\phi)P_G(\pi)&P_A((\phi,\pi) \mapsto (\phi',\pi')) \\
&= P_S(\phi)P_G(\pi) \: \min\Bigl(1,\frac{P_S(\phi')P_G(\pi')}{P_S(\phi)P_G(\pi)}\Bigr) \\
&=  \: \min\bigl(P_S(\phi)P_G(\pi),P_S(\phi')P_G(\pi')\bigr) \\
&= P_S(\phi')P_G(\pi') \: \min\Bigl(\frac{P_S(\phi)P_G(\pi)}{P_S(\phi')P_G(\pi')},1\Bigr) \\
&= P_S(\phi')P_G(\pi')P_A((\phi',\pi') \mapsto (\phi,\pi)),
\end{align*}
as desired.
\end{proof}

\begin{theorem}
The Markov chain defined by transitions $P_M$ has $P_S$ as a stationary distribution.
\end{theorem}

\begin{proof}
It suffices to show that $P_M$ satisfies the detailed balance condition

\begin{align}
P_S(\phi) P_M(\phi \mapsto \phi') = P_S(\phi') P_M(\phi' \mapsto \phi).
\label{eq:det-bal}
\tag{$\star$}
\end{align}

We have
\begin{align*}
P_S&(\phi) P_M(\phi \mapsto \phi') \\
    &=\int \Bigl(
    P_S(\phi)P_G(\pi) P_H((\phi,\pi) \mapsto (\phi',\pi')) \\
    & \qquad \qquad P_A((\phi,\pi) \mapsto (\phi',\pi')) \Bigr)
    d\pi d\pi' \\
    &=\int \Bigl(
    P_S(\phi')P_G(\pi') P_H((\phi,\pi) \mapsto (\phi',\pi')) \\
    & \qquad \qquad P_A((\phi',\pi') \mapsto (\phi,\pi)) \Bigr)
    d\pi d\pi' \\
& \hspace{3cm} \text{(by Lemma~\ref{lem:PSPGPA-symmetric})} \\
    &=\int \Bigl(
    P_S(\phi')P_G(\pi') P_H((\phi',\pi') \mapsto (\phi,\pi)) \\
    & \qquad \qquad P_A((\phi',\pi') \mapsto (\phi,\pi)) \Bigr)
    d\pi d\pi' \\
& \hspace{3cm} \text{(by Lemma~\ref{lem:PH-symmetric})} \\
&=P_S(\phi') P_M(\phi' \mapsto \phi),
\end{align*}
and so \eqref{eq:det-bal} holds.
\end{proof}

We have established that for Markov processes generated by our transition, the desired distribution is stationary, by showing that it satisfies detailed balance.
This implies, when the chain is ergodic, that the MCMC procedure eventually generates samples which are arbitrarily close in total variation distance to the posterior distribution.


TITLE: Row and Column Sum-Preserving Matrix Tranformations
QUESTION [1 upvotes]: I have an $\infty \times n$ matrix with entries $a_{i,j}\ge 0$, with $1\le i<\infty, 1\le j\le n$ and the sum of the entries in the matrix is 1: $\sum_{i,j}a_{i,j}=1$.
What's the set of transformations (not necessarily linear) that I can apply that will provide another matrix of nonnegative entries that will preserve the row and column sums of the original matrix?
For a particular problem, I'm looking for a transform (if one exists) that maps $a_{i,j}\mapsto 0$ when $\frac{j}{\gcd(i,j)}$ is compositie while not altering the row and column sums: For each $j$, $\sum_{i\ge 1}a_{i,j} \mapsto
\sum_{i\ge 1}a_{i,j}$ and for each $i$, $\sum_{1\le j\le n}a_{i,j} \mapsto \sum_{1\le j\le n}a_{i,j}$.

REPLY [1 votes]: Here's a version of an algorithm I've used for a similar problem finding an $m \times n$ matrix with prescribed row and column sums and $0$'s in some specified locations. It's an iteration that converges to an answer. Since there's no closed form you'd have to be satisfied with a program (easy to write in any language). That could cause problems since you have infinitely many rows. But the tail must be  small since the matrix elements sum to $1$. Maybe it can be managed.
Given your starting $a_{i,j}$, record and remember the row and column sums. Then set all the entries that must be $0$ to $0$. Now loop on the rows, rescaling the nonzero values left so that the row sums are right. Then do the same for the columns.
Of course that will destroy the row sums. So rescale the rows again. Continue alternately rescaling rows and columns.
Most of the time this process will converge to a matrix with the right row and column sums.
I say "most of the time" because it may fail when the specified $0$ entries disconnect the matrix, in the following sense. Construct the bipartite graph with one vertex set for the rows and one for the columns. Join two vertices when their intersection is not one of the entries that must be $0$. Hope that graph is connected. If it's not you may run into problems. I haven't checked connectedness for the particular zero set in your problem.

Addition, in answer to comment.
I'm not surprised that any such transformation must be nonlinear. The only algebraic invariant that makes sense is that it's homogeneous: you can scale all the matrix entries and the corresponding sums by the same positive amount and the answer scales the same way.
The algorithm is meant to provide a reasonable answer to a restricted transportation problem. Suppose you have quantities $r_i$ at $m$ sources and want to move a total of $c_j$ to each of $j$ sinks. Of course the $r_i$ and $c_j$ must have the same sum, which you can assume is $1$. Then each possible routing scheme $A$ is an $m \times n$ matrix with the prescribed row an column sums: ship $A_{ij}$ units from source $i$ to sink $j$. The set of possible schemes is a convex polytope in matrix space: a transportation polytope. In operations research there's usually a cost function involved in shipping from source $i$ to sink $j$ and a question of optimization; I was't interested in that aspect so made shipping free. Then in some sense the assignment $A_{ij} = r_i c_j$ is the most natural solution assuming some kind of independence in routing choices. The most extreme solutions correspond to vertices of the polytope and in turn to spanning trees in the bipartite graph defined above.
The algorithm finds what I think is the most natural solution when some routes are forbidden: some matrix entries that must be $0$.
There are probably theorems that codify the notion of "most natural"; I didn't need them so never looked for them.
Google finds many links to papers on transportation polytopes.
I hope this background helps.


text:
\begin{document}

\title{Spherical parking functions, uprooted trees, \\
and yet another way to count $n^n$}

\author{Anton Dochtermann}
\address{Department of Mathematics, Texas State University, San Marcos}
\email{dochtermann@txstate.edu}

\date{\today}

\begin{abstract}
Parking functions are a widely studied class of combinatorial objects, with connections to several branches of mathematics.  The number of parking functions of length $n$ is given by $(n+1)^{n-1}$, which by Cayley's formula is equal to the number of spanning trees of the complete graph $K_n$.   Parking functions can be identified with the standard monomials of $M_n$, a certain monomial ideal in the polynomial ring $S = {\mathbb K}[x_1, \dots, x_n]$ with a set of generators indexed by the nonempty subsets of $[n] = \{1,2,\dots,n\}$.

In this note we study $M^{(n-2)}_n$, a subideal of $M_n$ obtained by allowing only generators corresponding to \emph{proper} subsets of $[n]$. The set of standard monomials of $M^{(n-2)}_n$ consists of the usual parking functions as well as a collection of sequences we term `spherical parking functions' of length $n$, which we show has cardinality $(n-1)^{n-1}$. The proof is recursive and relies on known enumerative results regarding a certain class of `uprooted trees'.   As a corollary we obtain a seemingly new combinatorial identity for $n^n$, reminiscent of a result of Benjamin and Juhnke.  We discuss combinatorial refinements of our result and a conjectural relationship between the degree of spherical parking functions and inversions in uprooted trees.

\end{abstract}
\maketitle

\section{Introduction}\label{sec:intro}

A \defi{parking function of size $n$} is a sequence $(a_1, a_2, \dots, a_n)$ of nonnegative integers such that its rearrangement $c_1 \leq c_2 \leq \cdots \leq c_n$ satisfies $c_i < i$.  This seemingly innocent construction turns out to have connections and applications to many areas of mathematics.  In \cite{KonWei} it is shown that there are $(n+1)^{n-1}$ parking functions of size $n$, which by Cayley's formula is equal to the number of spanning trees of the complete graph $K_{n+1}$. This correspondence can be generalized to the case of arbitrary graphs $G$ in the context of sequences of integers known as \emph{$G$-parking functions} \cite{Gab}.

On the algebraic side parking functions of length $n$ can be identified with the standard monomials of the \emph{parking function ideal} $M_n$, a monomial ideal living in the polynomial ring $S = {\mathbb K}[x_1, x_2, \dots, x_n]$ for some fixed filed ${\mathbb K}$.  In \cite{PosSha} a parking function ideal $M_G$ is defined for an arbitrary graph $G$ on $n+1$ vertices, but we will only need the case $G = K_{n+1}$ is the complete graph. A generating set of monomials for $M_n$ (in fact for any $M_G$) is indexed by all nonempty subsets of $[n] = \{1,2,\dots, n\}$ according to the following construction.  For any $\emptyset \neq \sigma \subseteq [n]$  define the monomial $m_\sigma$ by
\begin{equation}
m_\sigma =  \prod_{i \in \sigma} x_i^{n-|\sigma|+1}.
\end{equation}
The ideal $M_n$ is then by definition the ideal (minimally) generated by these monomials:
\[M_n = \langle m_\sigma: \emptyset \neq \sigma \subseteq [n] \rangle.\]
In this paper we study a subideal of $M_n$ generated by monomials corresponding to the \emph{proper} subsets of $[n]$.

\begin{defn}
The $M_n^{(n-2)}$ is the ideal in $S = {\mathbb K}[x_1, \dots, x_n]$ given by
\[M_n^{(n-2)} = \langle m_\sigma: \emptyset \neq \sigma \subsetneq [n] \rangle.\]
\end{defn}

The generators of $M_n^{(n-2)}$ then correspond to all subsets of size at most $n-1$, which can be thought of as the $(n-2)$-skeleton of a simplex (hence the notation).  In \cite{DocOne} the $k$-skeleton ideals of $G$-parking function ideals are introduced, where standard monomials and minimal free resolutions are studied for the case $k=1$.

\begin{example}
For $n = 4$ we have
\[M_4^{(2)} = \langle x_1^4, x_2^4, x_3^4, x_4^4, x_1^3x_2^3, x_1^3x_3^3, x_1^3x_4^3, x_2^3x_3^3, x_2^3x_4^3, x_3^3x_4^3, x_1^2x_2^2x_3^2, x_1^2x_2^2x_4^2, x_1^2x_3^2x_4^2, x_2^2x_3^2x_4^2 \rangle\]
\end{example}

We briefly recall some basic notions in commutative algebra that may be unfamiliar to the reader (see \cite{MilStu} for more details). Suppose $I \in S$ is an ideal such that $S/I$ is finite dimensional as a ${\mathbb K}$-vector space.  Then we say that $I$ is an \emph{Artinian} ideal and a basis for $S/I$ is called the set of \defi{standard monomials} of $I$ (determined by a chosen term order).  In the case that $I$ is monomial the standard monomials consist of those monomials $m$ such that $m$ is not divisible by any generators of $I$.

One can check that the standard monomials of $M_n$ are naturally in correspondence with the parking functions of length $n$, so that $\dim_{\mathbb K} S/M_n = (n+1)^{n-1}$.  In \cite{PosSha} this approach is used to define $G$-parking function ideals for an arbitrary graph $G$, where the ${\mathbb K}$-dimension of $S/M_G$ is equal to the number of spanning trees of $G$.  Our main result concerns the standard monomials of $M_n^{(n-2)}$.

\begin{thm}\label{thm:main}
The number of standard monomials of $M_{n}^{(n-2)}$ is given by
\begin{equation} \label{eq:codimoneskelstd}
(n+1)^{n-1} + (n-1)^{n-1}.
\end{equation}
\end{thm}

Since the ideal $M_n^{(n-2)}$ is contained in $M_n$ we have that the set of standard monomials of $M_n$ are contained in those of $M_n^{(n-2)}$.  Hence the $(n+1)^{n-1}$ term appearing in Equation \ref{eq:codimoneskelstd} can be seen to correspond to the set of parking functions of length $n$.  Define a \defi{spherical parking function} of length $n$ to be (the exponent vector of) a standard monomial of $M_n^{(n-2)}$ that is \emph{not} a parking function (recall that the $(n-2)$- skeleton of an $(n-1)$-dimensional simplex is a sphere).

\begin{cor}\label{cor:spherical}
The number of spherical parking functions of length $n$ is given by $(n-1)^{n-1}$.
\end{cor}

The set of spherical parking functions also has a straightforward combinatorial description as a sequence $(a_1, a_2, \dots, a_n)$ of nonnegative integers whose rearrangement $c_1 \leq c_2 \leq \cdots \leq c_n$ satisfies $c_1 = 1$  and $c_i < i$ for all $i = 2, 3,\dots, n$.

\begin{example}
Again for $n=4$ we have that $M_4^{(2)}$ contains $5^3 + 3^3 = 152$ standard monomials.  Among these, $125$ correspond to parking functions of length $n$ (for instance the parking function $(1,0,1,2)$ corresponds to the monomial $m = x_1x_3x_4^2$).  The other 27 standard monomials (the spherical parking functions) can be described as follows
\[x_1x_2x_3x_4, \hspace{.15  in} x_1^2x_2x_3x_4 \times 4, \hspace{.15  in} x_1^3x_2x_3x_4 \times 4, \hspace{.15  in}x_1^2x_2^2x_3x_4 \times 6, \hspace{.15  in} x_1^3x_2^2x_3x_4 \times 12.\]
\end{example}

The value $(n-1)^{(n-1)}$ appearing in (\ref{eq:codimoneskelstd}) also has combinatorial interpretations. By an \defi{uprooted tree} on $[n] = \{1, 2, \dots, n\}$ we mean a rooted tree with vertex set $[n]$ in which the root is larger than all of its immediate descendants.  In \cite{ChaDulGui} it is shown that the number of uprooted trees on $[n]$ is given by $(n-1)^{(n-1)}$, although do not know of a direct bijection between these object and the spherical parking functions (see Section \ref{sec:inv} for further discussion).  However, an explicit formula for the standard monomials of $M_{n}^{(n-2)}$ coming directly from the definition provides `another way of counting' $(n-1)^{n-1}$ in the spirit of \cite{BenJuh}.

\begin{cor} \label{cor:main}
For any integer $n \geq 1$ we have
\begin{equation} \label{eq:choose}
(n-1)^{n-1} \; = \sum\limits_{\substack{0 \leq k_1 \leq 1 \\ 0 \leq k_1 + k_2 \leq  2 \\ \cdots \\ 0 \leq k_1 + \cdots + k_{n-2} \leq n-2}} {n \choose k_1}{n-k_1 \choose k_2} \cdots {n-(k_1 + \dots + k_{n-3}) \choose k_{n-2}},
\end{equation}
where $n > 1$ is an integer and $k_1, k_2, \dots, k_{n-2}$ are nonnegative integers.
\end{cor}

To illustrate the identity, for $n=4$ we have the possible values of $k_1$, $k_2$ given by
\begin{center}
\begin{tabular}{|c|c|}
\hline
$k_1$ & $k_2$ \\
\hline
0 & 0 \\
0&1 \\
0&2 \\
1&0 \\
1&1 \\
\hline
\end{tabular}
\end{center}

and the formula becomes
\smallskip
\begin{align*}
3^3 &= {4 \choose 0}{4-0 \choose 0} + {4 \choose 0}{4-0 \choose 1} + {4 \choose 0}{4-0 \choose 2} + {4 \choose 1}{4-1 \choose 0} + {4 \choose 1}{4-1 \choose 1} \\[8 pt]
&= 1 + 4 + 6 + 4 + 12 = 27.
\end{align*}

We note that the formula for $n^n$ involves $C_n$ many terms in the summation, where
\[C_n = \frac{1}{n+1} {2n \choose n}\]
\noindent
is the $nth$ \emph{Catalan number}. In \cite{BenJuh} Benjamin and Juhnke established a similar looking identity:

\begin{equation} \label{eq:another}
(n-1)^{n-1} \; = \sum\limits_{\substack{0 \leq k_1 \leq 1 \\ 0 \leq k_1 + k_2 \leq  2 \\ \cdots \\ 0 \leq k_1 + \cdots + k_{n-2} \leq n-2}} \frac{(n-1)!}{k_1!k_2! \cdots k_{n-2}!},
\end{equation}

\noindent
where $n > 1$ is an integer and $k_1, k_2, \dots, k_{n-2}$ are nonnegative integers.  In \cite{DWC} this formula was generalized to an identity involving $n^m$, where $m<n$.  Note that Equation \ref{eq:another} involves a summation over the same indexing set as Equation \ref{eq:choose}, but one can check that the terms in the summation are not the same (even up to reordering).  For instance if $n=4$ Equation \ref{eq:another} becomes

\begin{align*}
3^3 &= \frac{3!}{0!0!}+\frac{3!}{0!1!}+ \frac{3!}{0!2!}+ \frac{3!}{1!0!} + \frac{3!}{1!1!} \\[8 pt]
&= 6 + 6 + 3 + 6 + 6 = 27.
\end{align*}

In the next section we provide proofs of our results. We end with a discussion regarding other combinatorial refinements of our identities, as well as conjectural connection between the degree of spherical parking functions and certain inversions on uprooted trees.

\section{Proofs and discussion}

In this section we prove Theorem \ref{thm:main} and Corollary \ref{cor:main}.  Recall that the ideal $M_n^{(n-2)}$ is generated by all monomials $m_\sigma$ where $\sigma \subseteq [n]$ and $1 \leq |\sigma | \leq n-1$.  In particular $M_n^{(n-2)}$ has the same generators as $M_n$ \emph{except} for the single monomial $m_{\{1,2,\dots,n\}} = x_1x_2 \cdots x_n$.

\begin{proof}[Proof of Theorem \ref{thm:main}]

Since $M_n^{(n-2)} \subset M_n$ we have that the set of standard monomials of $M_n$ are contained in those of $M_n^{(n-2)}$.  We know that there are $(n+1)^{n-1}$ standard monomials of $M_n$ so we need to show that there are $(n-1)^{(n-1)}$ monomials $m$ such that  $m$ is divisible by $x_1x_2 \cdots x_n$  and yet not divisible by any generators of $M_n^{(n-2)}$.  Such monomials have the form
\[x_1^{a_1} x_2^{a_2} \cdots x_n^{a_n} \text{, where $a_i \geq 1$ for all $i$},\]
\noindent
with certain conditions on the values of $a_i$. Let  ${\mathcal M}_n^+$ denote the set of such monomials. Hence ${\mathcal M}_n^+$ naturally corresponds to the set of spherical parking functions of length $n$, but here we prefer the monomial perspective.   As we have seen it suffices to show that $|{\mathcal M}_n^+| = (n-1)^{n-1}$.

Note that an element of ${\mathcal M}_n^+$ is obtained by taking the monomial $x_1x_2 \cdots x_n$, choosing a set $S \subseteq [n]$, and then raising each variable corresponding to a element of $S$ by a nonzero amount.  Let $F(n,s)$ denote the number of standard monomials in ${\mathcal M}_{n+1}^+$ that correspond to raising the exponent of all variables in a chosen (and fixed) set of $n-s$ variables, fixing the exponents of all other variables to be one.  We claim that $F(n,s)$ satisfies the recursion
\begin{equation}\label{eq:recursion}
F(n,s) = \sum_{j=0}^{n-s} {n-s \choose j} F(n-1, s+j-1),
\end{equation}
with $F(1,1) = 1$ and $F(n,0) = 0$ for $n \geq 1$.

We have $F(1,1) = 1$ since if we have only two variables, we get a single contribution to the set ${\mathcal M}_2^+$ (standard monomials of $M_2^{(1)}$ that are elements of $M_2$) by `changing zero variables', namely from the monomial $x_1x_2$.  In fact $F(n,n) = 1$ for all $n \geq 1$.

Similarly we have $F(n,0) = 0$ for $n\geq 1$ since if we have $n+1$ variables and raise the exponent of any $n$ variables we end up with a monomial in $M_{n+1}^{(n-2)}$.  This follows from the fact that $M_{n+1}^{(n-2)}$ contains the monomial $m_S$ for all $S \subseteq [n]$ with $|S| \leq n$.

Next, to establish the recursion, suppose $m$ is a standard monomial in ${\mathcal M}_{n+1}^+$ resulting from raising the exponents in some set of $t = n-s$ variables.  Without loss of generality suppose that the set is $\{1,2, \dots, t\}$ so that
\[m = x_1^{a_1}x_2^{a_2} \cdots x_t^{a_t}x_{t+1} \cdots x_{n+1},\]
\noindent
where each $a_i > 1$.  Suppose $j$ of these variables have exponent equal to 2 (so that $t-j$ variables have exponent 3 or larger).  Take the subset of variables consisting of those with exponent 3 or larger and subtract one from each of those exponents to get a new monomial $m^\prime$ (on that subset of variables, shifting indices if necessary) that is now a standard monomial of $M_n^{(n-2)}$ with the property that $t-j= n-s-j$ variables have exponent (strictly) larger than 1.  The number of such standard monomials is $F(n-1, s+j-1)$.  If we add up ${n-s \choose j} F(n-1, s+j-1)$ for all possible choices of $j$ variables with exponent equal to 2 we get $F(n,s)$ (see Example \ref{ex:recursion} for an illustration).  The recursion formula (\ref{eq:recursion}) follows.

Now, in \cite{Tak} it is shown that $F(n,s)$ counts the number of forests with $n$ labeled vertices consisting of $s$ distinct trees such that $s$ specified vertices belong to distinct trees.  There it is shown that
\[F(n,s) = s n^{n-s-1}.\]
In \cite{ChaDulGui} the authors show that the number of uprooted trees on vertex set $[n+1]$ is given by $n^n$, and they use this description to prove a number of combinatorial interpretations .  By considering the parameter `degree $s$ of the root', they obtain the identity:
\begin{equation} \label{eq:treesbydeg}
n^n = \sum_{s=1}^n {n+1 \choose s+1} s n^{n-s-1}.
\end{equation}

We have seen that $s n^{n-s-1}$ counts the number of monomials in ${\mathcal M}_{n+1}^+$ that correspond to raising the exponent of all variables in a chosen (and fixed) set of $n-s$ variables (while fixing all other variables to have exponent one).  There are ${{n+1} \choose {n-s}} = {n+1 \choose s+1}$ ways to choose which variables to raise, and hence our result follows.
\end{proof}

\begin{example}\label{ex:recursion}
To aid the reader we provide an illustration of the identity in Equation \ref{eq:recursion} for the case
\[F(4,2) = F(3,1) + 2 F(3,2) + F(3,3).\]
\noindent
On the left hand side we are considering monomials with $4+1 = 5$ variables. Assuming as above that the set of variables with raised coefficients is $\{x_1, x_2\}$ we have
\begin{equation}
\begin{split}
|\{33111,34111,43111,32111,23111,42111,24111,22111\}| \\
= |\{2211,2311,3211\}| + 2|\{2111, 3111\}| + |\{1111\}|,
\end{split}
\end{equation}

where monomials are represented by their exponent vectors.
\end{example}

We next turn to the proof of Corollary \ref{cor:main}.

\begin{proof}[Proof of Corollary \ref{cor:main}]
Our strategy will be to show that the right hand side of the identity naturally counts the elements of ${\mathcal M}_n^+$. Theorem \ref{thm:main} will then give the result. Note that if $x_1^{a_1} x_2^{a_2} \cdots x_n^{a_n} \in {\mathcal M}_n^+$, then no exponent can satisfy $a_i \geq n$ since $x_i^n$ is a generating monomial of $M_G^{(n-1)}$.  We can have at most one exponent satisfying $a_i = n-1$, at most two exponents satisfying $a_j = n-2$, and so on.  Let $k_1$ denote the number of exponents $a_i$ such that $a_i = n-1$, $k_2$ the number of exponents $a_j$ such that $a_j = n-2$, etc.  Hence the number of standard monomials of the form $x_1^{a_1} x_2^{a_2} \cdots x_n^{a_n}$ is given by the expression on the right hand side of (\ref{eq:choose}), which by Theorem \ref{thm:main} is equal to $(n-1)^{n-1}$.
\end{proof}

\section{Inversions and further questions.}\label{sec:inv}

We end with a brief discussion regarding some refinements of $n^n$ via combinatorial parameters. As we have seen Equation (\ref{eq:treesbydeg}) taken from \cite{ChaDulGui} enumerates uprooted trees via the parameter `degree $s$ of the root' (see Figure \ref{fig:trees_degree}).  In our context this corresponds to enumerating elements of ${\mathcal M}_{n+1}^+$ in terms of `raising the exponent of $n-s$ variables' .

\begin{figure}[h]
\includegraphics[scale = .475]{trees_degree.pdf}
\caption{Uprooted trees on $[4]$ by degree of the root, $27 = 1 + 8 + 18$.}
\label{fig:trees_degree}
\end{figure}

In the context of parking functions, a natural parameter to consider is the \emph{degree} (sum of the entries). Kreveras \cite{Kre} showed that the number of length $n$ parking functions of degree ${n+1 \choose 2} - k$ is equal to the number of spanning trees of $K_{n+1}$ with $k$ \emph{inversions}.  To describe this notion suppose $T$ is a tree on vertex set $\{1, 2,\dots,n\}$ and fix a root at 1.  Then an \defi{inversion} in $T$ is a pair $(\alpha, \beta)$ of vertices in $T$ such that $\beta$ is a descendant of $\alpha$ and $\alpha > \beta$. Kreweras established his result by relating parking functions to the \emph{Tutte polynomial} of $K_n$, and hence to the \emph{external activity} of trees.  A bijection between the set of trees with $k$ inversions and the set of trees with $k$ externally active edges was given by Beissinger in \cite{Bei}.  A geometric interpretation of parking functions counted by degree is also provided by the so-called `Pak-Stanley' labeling of regions in the Shi arrangement \cite{Sta}.

Many of these results can be generalized to arbitrary graphs.  If we let $T_G(x,y)$ denote the Tutte polynomial of a (connected) graph $G$, Merino \cite{Mer} showed that $T_G(1,y) = \sum_{i=0}^r a_i y^i$, where $r-a_i$ is the number of $G$-parking functions of degree $i$, and $r$ is the circuit rank of $G$.  As a corollary we see that the number of $G$-parking functions of degree $i$ is given by the number of spanning trees of $G$ with external activity $r - i$.  A  bijective proof of this result was provided by Cori and Le Borgne in \cite{CorLeB}.  Gessel \cite{Ges} has generalized the notion of inversion in the context of \emph{$\kappa$-inversions} for spanning trees of an arbitrary graph $G$.  He shows that the coefficient $a_i$ of $T_G(x,y)$  is also given by the number of spanning trees of $G$ with $\kappa$-number $i$.  In \cite{PerYanYu} Perkinson, Yang, and Yu give a bijection between $G$-parking functions and spanning trees that preserves degree and the number of $\kappa$-inversions.

In our context we can enumerate elements of ${\mathcal M}_{n}^+$ according to degree.  As we have seen elements in ${\mathcal M}_{n}^+$ are in bijection with \emph{uprooted trees}, and so degree gives a a new parameter on such objects.   Small examples indicate a possible connection to inversions that is analogous to the classical case.  If $T$ is a birooted tree on $[n]$ with root $r$, define a \defi{surface inversion} to be a pair of \emph{non-root} vertices $(\alpha, \beta)$ such that $\beta$ is a descendant of $\alpha$ and $\alpha > \beta$.  We then have the following conjectural correspondence.

\begin{conj}\label{conj:trees}
The number of spherical parking functions of length $n$ of degree ${n \choose 2} - k +1$ is equal to the number of birooted trees on vertex set $[n]$ with $k$ surface inversions.
\end{conj}

The correct proof of this conjecture would presumably involve a bijective proof of Theorem \ref{thm:main} that preserves degree and surface inversions.

\begin{figure}[h]
\includegraphics[scale = .475]{trees_oneinversion.pdf}
\caption{Uprooted trees on $[4]$ with one surface inversion.}
\label{fig:trees_inversions}
\end{figure}

\begin{example}
As an example of the above discussion, the relevant monomials in ${\mathcal M}_4^+$ (described as spherical parking functions)  are $1123$ (all 12 permutations), $1122$ (all six permutations), $1113$ (all four permutations), $1112$ (all four permutations), and $1111$.  Theorem \ref{eq:treesbydeg} gives the identity
\[27 = 18 + 8 + 1,\]
since 18 monomials have two variables raised, 8 have a single variable raised, and one has no variable raised.  In terms of labeled rooted trees discussed in \cite{ChaDulGui} this corresponds to the number of trees with particular root degree, see Figure \ref{fig:trees_degree}.

If we instead count using the \emph{degree} statistic then we get
\[27 = 12 + 10 + 4 + 1.\]
\noindent
Conjecture \ref{conj:trees} states that this the same as counting uprooted trees in terms of the number of surface inversions.  See Figure \ref{fig:trees_inversions} for the list of uprooted trees on vertex set $[4]$ with one surface inversion.
\end{example}

\begin{rem}
After this paper was prepared we found a connection between our objects of study and the \emph{${\bf u}$-vector parking functions} introduced by Yan in \cite{YanGen}.  For this suppose ${\bf u} = (u_1, u_2 , \dots u_n) \in {\mathbb N}^n$ is a vector of nonnegative integers.  Then, a sequence $(a_1, a_2, \dots, a_n)$ of nonnegative integers is a \defi{${\bf u}$-parking function} if its rearrangement $c_1 \leq c_2 \leq \cdots \leq c_n$ satisfies
\[c_j < \sum_{k=1}^j u_k\]
for all $1 \leq j \leq n$.

Let $PF({\bf u})$ denote the set of all ${\bf u}$-parking functions. Observe that the usual parking functions are recovered for the case ${\bf u} = (1,1, \dots, 1)$.  One can check that the standard monomials of $M_n^{(k)}$ are naturally in one-to-one correspondence with the set of ${\bf u}_{n,k}$ vector parking functions, where
\[{\bf u}_{n,k} = (n-k, \underbrace{0, 0 \dots, 0}_\text{$n-k-1$ times}, \underbrace{1,1, \dots, 1}_\text{$k$ times}).\]
\noindent
By results of Yan \cite{Yan} the number of such monomials  is given by
\[\sum_{j=0}^k {n \choose j} (k + 1 - j)(k+1)^{j-1}(n-k)^{n-j}\]
\noindent
By Theorem \ref{thm:main}, for $k=n-2$ this formula simplifies to  $(n+1)^{n-1} + (n-1)^{n-1}$ although we do not see a clear reduction.
\end{rem}

We remark that for an arbitrary graph $G$ (with specified root vertex), the standard monomials of the $k$-skeleton ideal $M_G^{(k)}$ provide a natural blending of the ${\bf u}$-parking functions with the $G$-parking functions from \cite{PosSha}.  One could hope for a generalization of Conjecture \ref{conj:trees} that extends to such objects, involving a generalization of the $\kappa$-inversions of \cite{Ges} and perhaps the burning algorithm described in \cite{PerYanYu} and \cite{GayHop}.  As far as we know these notions have not been explored.


TITLE: Finding joint distribuion of $X,Y$ given joint moment generating function
QUESTION [0 upvotes]: Let $M_{X,Y}(t_1,t_2) = {t_1^2+t_2^2\over 2}$ be the joint moment generating function of $X,Y$. Is it possible to find the joint density function of $X,Y$?
I´m having a really difficult time in this problem. I don´t know hot to compute the joint density. I would really appreciate any ideas or suggestions.

REPLY [0 votes]: The MGF is usually defined by $M_(X_1,X_2)(t_1,t_2)=Ee^{t_1X_1+t_2X_2}$. With this definition we have  $M_(X_1,X_2)(0,0)=1$ so your function is not a  MGF.
I suppose you want to consider $M_(X_1,X_2)(t_1,t_2)=Et_1^{X_1}t_2^{X_2}$. I will answer the question taking this as the definition.
Solution based on guess work: $(X,Y)$ takes just two values $(0,2)$  and$(2,0)$ with probability $\frac 1 2 $ each.


text:
\begin{document}
\title{A one-relator group with long lower central series}
\author{Roman Mikhailov}
\address{Chebyshev Laboratory, St. Petersburg State University, 14th Line, 29b,
Saint Petersburg, 199178 Russia and St. Petersburg Department of
Steklov Mathematical Institute} \email{rmikhailov@mail.ru}
\thanks{This research is supported by the Chebyshev
Laboratory  (Department of Mathematics and Mechanics, St.
Petersburg State University)  under RF Government grant
11.G34.31.0026, and by JSC Gazprom Neft, as well as by the RF
Presidential grant MD-381.2014.1.}
\begin{abstract}
A one-relator group with lower central series of length $\omega^2$
is constructed. This answers a problem of G. Baumslag.
\end{abstract}
\maketitle

\section{Introduction} For a group $G$, the lower central series are defined
inductively as follows: $$\gamma_1(G)=G,\
\gamma_{\alpha+1}(G)=[\gamma_\alpha(G),G]$$ and
$\gamma_\tau(G)=\cap_{\alpha<\tau}\gamma_\alpha(G)$ for a limit
ordinal $\tau$. The smallest ordinal $\alpha$, such that
$\gamma_{\alpha}(G)=\gamma_{\alpha+1}(G)$ is called the {\it
length of the lower central series of $G$}. First examples of
finitely presented groups with the lower central series of length
greater than $\omega$ were constructed by J. Levine \cite{Levine}.
T. Cochran and K. Orr constructed examples of 3-manifold groups
with lower central series of length greater than $\omega$ in
\cite{CochranOrr}.

G. Baumslag asked the following question (\cite{Baumslag}, Problem
10): {\it Is the lower central series of a 1-relator group of
length at most $\omega n$ for some finite ordinal $n$?}

Consider the following one-relator group
$$
G=\langle a,b\ |\ a^{b^2}=aa^{3b}\rangle
$$
This group is a semidirect product $F_2\rtimes \mathbb Z$, where a
generator of the cyclic group acts on $F_2=F(x,y)$ as the
automorphism

\begin{align*}
& x\mapsto y\\
& y\mapsto xy^3
\end{align*}
The main result is the following
\begin{theorem}
The group $G$ has the lower central series of length $\omega^2$.
\end{theorem}

The group $G$ provides an example of a one-relator group with
length greater than $\omega n$ for any finite $n$, thus answering
the above problem of G. Baumslag.

\section{Proof of theorem 1}
Let $\Gamma$ be a group and $M$ a $\mathbb Z[\Gamma]$-module. $M$
is called {\it residually nilpotent} if
$$
\cap_kM\Delta^k(\Gamma)=0,
$$
where $\Delta^k(\Gamma)$ is the $k$th power of the augmentation
ideal $\Delta(\Gamma):=\ker\{\mathbb Z[\Gamma]\to \mathbb Z\}.$
\begin{lemma}
Let $M$ be a finitely generated free abelian group of rank $m\geq
1$ and an infinite cyclic group with generator $t$ acts on $M$ as
a matrix $A\in GL_m(\mathbb Z)$. Suppose that the product of any
collection of eigenvalues of $A-Id$ is not $\pm 1$. Then the
module $M$ is residually nilpotent $\mathbb Z[\langle
t\rangle]$-module.
\end{lemma}
\begin{proof}
Let $\Delta$ be the augmentation ideal of $\mathbb Z[\langle
t\rangle]$. Observe that, for any $n\geq 1$, the power $\Delta^n$
is the ideal generated by $(1-t)^n.$ Consider the submodule
$N:=M\Delta^\omega=\cap_{k\geq 1}M\Delta^k$ and assume that $N\neq
0.$ The Krull intersection theorem implies that $NI=N$, hence
$N(1-t)=N$ and the restriction of $(1-t)$ on $N$ defines a
bijection on $N$. Hence $det(A-Id)=\pm 1=\beta_1\dots \beta_k$,
where $(\beta_1,\dots, \beta_k)$ is some collection of eigenvalues
of $A-Id$. This gives a needed contradiction. Hence $N=0$.
\end{proof}

Now we prove theorem 1. As usual, if $x,y,a_1,\dots,a_{k+1}$ are
elements of a group $G$ we set $[x,y] =x^{-1}y^{-1}xy$,
$x^y=y^{-1}xy$ and define
$$[a_1,a_2,\dots, a_{k+1}]=[[a_1,\dots,a_{k}], a_{k+1}]\ (\, k>1).$$ For $n\geq 1$, define $T_n$ as a normal
closure in $G$ of all elements of the form $[x_1,\dots,x_n],$
where all $x_i$-s are either $a$ or $a^b$. We have
$$
H_1=\langle a\rangle^G=F(a,a^b),\ H_2=\langle [a,a^b]\rangle^G,\
H_3=\langle [a,a^b,a], [a,a^b,a^b]\rangle^G
$$
etc. Observe that
$$
H_{n}=\gamma_n(H_1).
$$
We claim that, for all $n\geq 1,$
\begin{equation}\label{omegan}
\gamma_{n\omega}(G)\subseteq H_{n+1}
\end{equation}

First we prove (\ref{omegan}) for $n=1$. Consider the quotient of
$G$ by the normal closure of the element $[a,a^b]:$
$$
H=\langle a,b\ |\ a^{b^2}=aa^{3b}, [a,a^b]=1\rangle
$$
The group $H$ is the metabelian polycyclic group $(\mathbb Z\oplus
\mathbb Z)\rtimes \mathbb Z$, where the generator of the cyclic
group $\mathbb Z$ acts on $\mathbb Z\oplus \mathbb Z$ as the
matrix
\begin{equation}\label{matr}
U:=\begin{pmatrix}
0 & 1 \\
1 & 3
\end{pmatrix}
\end{equation}
The eigenvalues of the matrix $U-Id$ are $
\frac{3\pm\sqrt{13}}{2}-1.$ By lemma 2, the group $H$ is
residually nilpotent. Lets show now that $[a,a^b]\in
\gamma_\omega(G)$. Observe that the element $a$ is a generalized
3-torsion element, i.e. for every $n$ there exists $t(n)$, such
that $a^{3^{t(n)}}\in\gamma_n(G)$. Indeed, since
\begin{equation}\label{3tor} a^3=[a,b^2]^{b^{-1}},\end{equation}
$$
a^{3^k}\in\gamma_{k+1}(G),\ k\geq 1.
$$
Therefore, the element $[a,a^b]$ also is a generalized 3-torsion
element. On the other hand, taking the commutator with $a^b$ of
the both sides of the relation, we get
$$
[a^{b^2},a^b]=[a,a^b]^{a^{3b}},
$$
hence $[a^b,a]=[a,a^b]^{a^{3b}b^{-1}}$ and
\begin{equation}\label{relation}
[a^b,a]^2=[a,a^b,a^{3b}b^{-1}]\end{equation}  Therefore $[a,a^b]$
is a generalized 2-torsion element. This implies that $[a,a^b]\in
\gamma_\omega(G)$. Since $H$ is residually nilpotent, we conclude
that $\gamma_\omega(G)=\langle [a,a^b]\rangle^G,$
$G/\gamma_\omega(G)=H$ and the statement (\ref{omegan}) is proved
for $n=1$. It is easy to see at this point that the group $G$ has
length greater than $\omega$. Indeed, since $H_2(H)=\mathbb Z/2$
and $H_2(G)=0$, we get
$$
H_2(H)=H_2(G/\gamma_\omega(G))\simeq
\gamma_{\omega}(G)/\gamma_{\omega+1}(G)\simeq \mathbb Z/2
$$
and
$$
[a,a^b]\in\gamma_{\omega}(G)\setminus\gamma_{\omega+1}(G),\
[a,a^b]^2\in\gamma_{\omega+1}(G).
$$

Now suppose that the statement (\ref{omegan})  holds for a given
$n$. Here we prove this statement for $n+1$.

We have an inclusion
$$
\gamma_{n\omega+k}(G)\subseteq [H_{n+1},\underbrace {G,\dots,
G}_k],\ k\geq 1
$$
Define the group $T:=G/H_{n+2}$. By the assumption of induction,
the term $\gamma_{n\omega}(T)$ lies in the quotient
$H_{n+1}/H_{n+2}$, the term $\gamma_{n\omega+k}(T)$ is a subgroup
of the quotient $([H_{n+1},\underbrace {G,\dots,
G}_k]H_{n+2})/H_{n+2}$ and can be presented as a product of
conjugates of the elements
$$
[\alpha, \underbrace {b,\dots, b}_k].H_{n+2},\ \alpha\in H_{n+1}.
$$
Consider the quotient $H_{n+1}/H_{n+2}$ as a $\mathbb Z[\langle
b\rangle]$-module. If the module $H_{n+1}/H_{n+2}$ is residually
nilpotent $\mathbb Z[\langle b\rangle]$-module, i.e.
$$
0=\cap_k (H_{n+1}/H_{n+2})\Delta^k,
$$
where $\Delta$ is the augmentation ideal of $\mathbb Z[\langle
b\rangle]$, then the intersection
$$
\cap_k([H_{n+1},\underbrace {G,\dots, G}_k]H_{n+2})/H_{n+2}
$$
is trivial in $T$ and the statement (\ref{omegan}) will follow for
$n+1$. Denote $M:=H_1/H_2.$ Since $H_1$ is a free group and
$H_i=\gamma_i(H_1),\ i\geq 1,$ the quotient $H_{n+1}/H_{n+2}$ is
the $n+1$-st Lie power $L^{n+1}(M)$ by Magnus-Witt theorem (see
\cite{Magnus}, \cite{Witt}). Lie powers are subgroups of tensor
powers $L^{n+1}(M)\subseteq \otimes^{n+1}(M)$ and therefore,
$$
\cap_k L^{n+1}(M)\Delta^k\subseteq \cap_k \otimes^{n+1}(M)\Delta^k
$$
Here we consider the tensor powers of $M$ as $\mathbb Z[\langle
b\rangle]$-modules with the diagonal action of $\langle b\rangle$.

For given matrices $A$ and $B$, the eigenvalues of their tensor
product $A\otimes B$ consist of all possible products of the
eigenvalues of $A$ and $B$. The eigenvalues of the matrix $U$ are
$\alpha_1=\frac{3+\sqrt{13}}{2},\ \alpha_2=\frac{3-\sqrt{13}}{2}.$
Hence, the eigenvalues of the the tensor power $U^{\otimes n+1}$
are $\pm \alpha_1^l$ and $\pm \alpha_2^s$ for some $l,s\geq 1$.
Suppose that one can present $\pm 1$ as a product of elements of
the form $\pm\alpha_1^l\pm 1$ and $\pm\alpha_2^s\pm 1$. Then the
product of their conjugates in the sense of quadratic
irrationality also is $\pm 1$. Since $\alpha_1^l\pm 1$ is the
conjugate quadratic irrationality to $\alpha_2^l\pm 1$, one can
present $\pm 1$ as a product of integers of the form
$$
(\alpha_1^l-1)(\alpha_2^l-1)\ \text{and}\
(\alpha_1^s+1)(\alpha_2^s+1)
$$
for different $l,s\geq 1$. Since $(\alpha_1-1)(\alpha_2-1)=-3$
divides $(\alpha_1^l-1)(\alpha_2^l-1)$ for any $l\geq 1$, the
absolute values of the terms $(\alpha_1^l-1)(\alpha_2^l-1)$ are
greater than 1. For an odd $s$, the product
$(\alpha_1^s+1)(\alpha_2^s+1)$ divides
$(\alpha_1+1)(\alpha_2+1)=3.$ For an even $s$, the terms
$(\alpha_1^s+1)$ and $(\alpha_2^s+1)$ are greater than 1. We get a
contradiction with assumption that one can present $\pm 1$ as a
product of elements of the form $\pm\alpha_1^l\pm 1$ and
$\pm\alpha_2^s\pm 1$.

We conclude that, for any $n$, the matrix $U^{\otimes n+1}-Id$
satisfies the conditions of lemma 2 and therefore the $(n+1)$st
tensor power $M^{\otimes n+1}$ is a residually nilpotent $\mathbb
Z[\langle b\rangle]$-module. This implies that the Lie power
$L^{n+1}(M)$ also is a residually nilpotent $\mathbb Z[\langle
b\rangle]$-module and the inductive step is done,
$$
\gamma_{(n+1)\omega}(G)\subseteq H_{n+2},
$$
proving (\ref{omegan}) for all $n$.

Since $H_1$ is a free group, it is residually nilpotent, and we
get
$$
\gamma_{\omega^2}(G)=\cap_n\gamma_{n\omega}(G)\subseteq
\cap_nH_{n+1}=\cap_n\gamma_{n+1}(H_1)=1.
$$
Since $[a^b,a]\in \gamma_{\omega}(G),$ we have
$$
1\neq [[a^b,a],a,\underbrace {[a^b,a],\dots, [a^b,a]}_k]\in
\gamma_{(k+1)\omega}(G),\ k\geq 1.
$$
This element clearly is non-trivial, since it is a basic
commutator in the basic elements $a,a^b$ of the free group
$H_1=F(a,a^b)$. Hence, the lower central series length of $G$ is
$\omega^2$ and theorem 1 is proved.

\begin{remark}
Observe that one can easily use the proof of theorem 1 to prove
the following more general result. Let $F_n$ be a free group and
an element $t$ acts on $F_n$ as an automorphism $\phi:F_n\to F_n$.
Define $\Phi=F_n\rtimes \langle t\rangle$. Suppose that any tensor
power $(F_{ab})^{\otimes k}\ (k\geq 1)$ is a residually nilpotent
$\mathbb Z[\langle t\rangle]$-module. Then
$\gamma_{\omega^2}(\Phi)=1$. In particular, if $\Phi$ is not
residually nilpotent, then the lower central series length of
$\Phi$ is exactly $\omega^2$.
\end{remark}

We add the next remark here due to its applications in the theory
of localizations \cite{Mikha}.
\begin{remark} The $2\omega$th lower central series
term of $G$ can be described as follows
\begin{equation}\label{2o}\gamma_{2\omega}(G)=\langle
[a,a^b,a],[a,a^b,a^b]\rangle^G\end{equation}
\end{remark}

We proved before that there is an inclusion
$\gamma_{2\omega}(G)\subseteq H_3.$ To see the converse inclusion,
observe that $a$ and $a^b$ are generalized 3-torsion elements,
hence $$[a,a^b,a]^{3^l}, [a,a^b,a^b]^{3^l}\subseteq
\gamma_{\omega+l}(G),\ l\geq 1.$$ However, the same is true for
powers of 2, due to (\ref{relation}). Therefore, $H_3\subseteq
\gamma_{2\omega}(G)$ and description (\ref{2o}) follows.

\vspace{.5cm}\noindent {\it Acknowledgements.} The author thanks
G. Baumslag, S.O. Ivanov, F. Petrov and K. Orr for discussions
related to the subject of the paper.