text:
\begin{document}

\title{Some continued fractions for $\pi$ and $G$}
\date{\today}
\author {Amrik Singh Nimbran}
\address{B3-304, Palm Grove Heights, Ardee City, Gurgaon, Haryana, INDIA }
\email{amrikn622@gmail.com}

\author {Paul Levrie}
\address{Faculty of Applied Engineering, University of Antwerp, Groenenborgerlaan 171, B-2020 Antwerpen;
Department of Computer Science, KU Leuven, P.O. Box 2402, B-3001 Heverlee, BELGIUM}
\email{paul.levrie@cs.kuleuven.be}

\subjclass[2010]{11A55, 11J70, 30B70, 40A25}
\keywords{Infinite series; Continued fractions; Pi; Catalan's constant}

\begin{abstract}
We present here two classes of infinite series and the associated continued fractions involving $\pi$ and  Catalan's constant based on the work of Euler and Ramanujan. A few sundry continued fractions are also given.
\end{abstract}

\maketitle

\section{Introduction}

A continued fraction is an expression of the general form

\[
b_0+\cfrac{a_1}{b_1+ \cfrac{a_2}{b_2+ \cfrac{a_3}{b_3+ \dotsb}}}
\]~\\
with $a_1,a_2,\ldots$ and $b_0, b_1, b_2,\ldots$ real numbers. The following space-saving notation can be used:
\[
\newcommand{\cfplus}{\mathbin{\genfrac{}{}{0pt}{}{}{+}}}
b_0+\frac{a_1}{b_1}\cfplus\frac{a_2}{b_2}\cfplus\frac{a_3}{b_3}\cfplus\dotsb
\]
or the shorter notation
$ b_0 +\cfK_{n=1}^\infty \displaystyle\frac{a_n}{b_n},$
where the letter K comes from the German word \emph{Kettenbr\"{u}che} (Ketten - chain, Br\"{u}che - fraction) for a continued fraction.

It is not known when, where and by whom continued fractions were first used, but the notion seems to be quite old. For an overview, see \cite{brezinski}. Euclid's algorithm for finding the greatest common devisor of two integers in effect converts a fraction into a terminating continued fraction. Continued fractions were used in India in the 6th century by Aryabhata to solve linear Diophantine equations and in the 12th century by Bhaskaracharya for solving the `Pell' equation. Rafael Bombelli (1526--1572) and Pietro Antonio Cataldi (1548--1626), both of Bologna (Italy), gave continued fractions for $\sqrt{13}$ and $\sqrt{18}$ respectively. The first continued fraction for a number other than the quadratic irrationals is due to William Brouncker (1620-84) who converted the product $\displaystyle \frac{4}{\pi}=\prod_{n=1}^\infty \frac{(2n +1)^2}{2n(2n+2)}$ submitted to him by John Wallis (1616--1703) in 1655 (see \cite[p.182]{wallis}):
\[
\frac{4}{\pi}=1 +\cfK_{n=1}^\infty \frac{(2n-1)^2}{2}=1+\cfrac{1^2}{2+ \cfrac{3^2}{2+ \cfrac{5^2}{2+ \dotsb}}} .
\]

It was one of the sequence of continued fractions found by Brouncker.\cite[p.307]{stedall}

\section{Euler's continued fractions for $\pi$}

Leonhard Euler (1707--1783) laid the foundation of the theory of continued fractions in a dozen papers (six published in his life-time and six posthumously) and chapter 18 of his \textit{Introductio in Analysin Infinitorum} (1748). Of these, five papers written between 1737 and 1780, with Enestr\"{o}m index numbers: E071, E123, E522, E593 and E745, are of interest to us here. Brouncker's fraction occurs in E071 \cite[\S 4]{euler71}.
Euler explains the method of converting an infinite series into a continued fraction in E593 \cite{euler593}. He gives among others this formula in \cite[\S7]{euler123}:

\[
\int_{0}^1 \frac{x^{n-1}{dx}}{1+x^m}=\cfrac{1}{n+ \cfrac{n^2}{m+ \cfrac{(m+n)^2}{m +\cfrac{(2m+n)^2}{m+ \cfrac{(3m+n)^2}{m+ \dotsb}}}}} .
\]

Now $\displaystyle \int \frac{x {dx}}{1+x^4}=\frac{1}{2}\arctan x^2 +C.$ So taking $m=4, \, n=2,$ we get Brouncker's fraction in disguise:

\[
\frac{\pi}{8}=\cfrac{1}{2+ \cfrac{2^2}{4+ \cfrac{6^2}{4 +\cfrac{10^2}{4+ \cfrac{14^2}{4+ \dotsb}}}}} .
\]

We find in \cite[\S18]{euler593} Euler's Theorem I:

\textit{If such an infinite series will have been proposed
\[
s=\frac{1}{\alpha} -\frac{1}{\beta} +\frac{1}{\gamma} -\frac{1}{\delta} +\frac{1}{\epsilon} -\cdots,
\]
a continued fraction of the form
\[
\frac{1}{s}=\alpha +{\cfrac{\alpha \alpha}{\beta-\alpha+ \cfrac{\beta\beta}{\gamma-\beta+ \cfrac{\gamma\gamma}{\delta -\gamma+ \dotsb}}}}
\]
can always be formed from it.}

Taking the Leibnitz series $\displaystyle \sum_{n=1}^\infty \frac{(-1)^{n-1}}{2n-1}=\frac{\pi}{4},$ Euler deduced Brouncker's fraction.

We find this theorem as Theorem II in \cite[\S23]{euler593}:

\textit{If the proposed series is of this form
\[
s=\frac{1}{ab} -\frac{1}{bc} +\frac{1}{cd} -\frac{1}{de} +\frac{1}{ef} -\cdots,
\]
from this, the following continued fraction springs forth
\[
\frac{1}{a s}=b +{\cfrac{ab}{c-a+ \cfrac{bc}{d-b+ \cfrac{cd}{e -c+\cfrac{de}{f-d \dotsb}}}}} .
\]
}
Noting that $\displaystyle \sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)(2n+1)}=\frac{\pi}{4} -\frac{1}{2},$ he deduced
\[
\frac{4}{\pi -2}=3 +{\cfrac{1\cdot3}{4+ \cfrac{3\cdot5}{4+ \cfrac{5\cdot7}{4+ \cfrac{7\cdot9}{4+\dotsb}}}}}
\]
or
\[
\frac{\pi}{2}-1 =\cfrac{2}{3+ \cfrac{1\cdot3}{4+ \cfrac{3\cdot5}{4+ \cfrac{5\cdot7}{4+ \dotsb}}}} .
\]

Note that the convergents of this continued fraction $\displaystyle c_n=\frac{p_n}{q_n}$ satisfy the following relations:
\[ p_0=q_0=1, \quad
p_{n+1}=(2n+3) p_n +(-1)^n 2 (2n-1)!!, \quad q_n=(2n+1)!!,
\]
for $n=0, 1, 2, \ldots$,
where $n!!=n \cdot (n-2) \cdot (n-4) \cdot \ldots \cdot 1$ for $n$ odd, with $(-1)!!=1$.

This continued fraction is a special case of a more general one which can be deduced from a formula of Ramanujan giving a continued fraction for a product of quotients of gamma function values \cite[p.227]{perron}:

\noindent{\textbf{General Formula}}: \textit{Let
\[
P_{-1}=\frac{1}{2}, \, \ P_{2m}=\prod_{k=1}^{m} \frac{2k(2k+2)}{(2k+1)^2}, \, \  P_{2m+1}=P_{2m}\cdot \frac{2m+2}{2m+3} \quad (m=0,1,2,\ldots).
\]
Then we have
\[
\frac{\pi}{4 P_n} -1 =\cfrac{(-1)^{n+1} 2}{4(n+2)+(-1)^n+ \cfrac{1\cdot3}{4(n+2)+ \cfrac{3\cdot5}{4(n+2)+ \cfrac{5\cdot7}{4(n+2) +\dotsb}}}}
\]
}

It may be compared with \cite[p.227, eq(30a)]{perron}. Note that these continued fractions give $\pi$ as a combination of the partial products of Wallis's formula and a continued fraction. For instance, for $n=4$ we have:
\[
\frac{\pi}{4}=\frac{2}{3}\cdot\frac{4}{3}\cdot\frac{4}{5}\cdot\frac{6}{5} \cdot   [ 1- \cfrac{2}{25+ \cfrac{1\cdot3}{24+ \cfrac{3\cdot5}{24+ \cfrac{5\cdot7}{24+\ddots}}}} ] .
\]

Euler gave the following two continued fractions and a couple more in \cite[\S31, \S33]{euler123}:
\[
\frac{\pi}{2}=1 +\cfrac{1}{1+ \cfrac{1\cdot2}{1+ \cfrac{2\cdot3}{1+ \cfrac{3\cdot4}{1+ \dotsb}}}} \ ,
\]
\[
\frac{\pi}{2}=2- \cfrac{1}{2+ \cfrac{1^2}{2+ \cfrac{2^2}{2+ \cfrac{3^2}{2+ \dotsb}}}} .
\]

He obtained this fraction with partial denominators $4n$ from his general formula in \cite[\S 18]{euler522} and \cite[\S36]{euler593}:
\[
\frac{6\sqrt{3}} {\pi} =3 +\cfrac{3\cdot1^2}{8+\cfrac{3\cdot3^2}{12+ \cfrac{3\cdot5^2}{16+ \cfrac{3\cdot7^2} {20+ \dotsb}}}} .
\]

These continued fractions with partial denominators $3n-2$ and $5n-3$ are due to Glaisher \cite{glaisher73, glaisher76}.

\[
\frac{2}{\pi}=1-\cfrac{1\cdot1}{4- \cfrac{2\cdot3}{7- \cfrac{3\cdot5}{10- \cfrac{4\cdot7} {13- \dotsb}}}} \ ,
\]
\[
\frac{3\sqrt{3}}{\pi}=2-\cfrac{2(1\cdot1)}{7- \cfrac{2(2\cdot3)}{12- \cfrac{2(3\cdot5)}{17- \dotsb}}} .
\]
They were converted from the following series:
\[
\frac{\pi}{2}=1 +\frac{1}{3} +\frac{1\cdot2}{3\cdot5} +\frac{1\cdot2\cdot3}{3\cdot5\cdot7} +\frac{1\cdot2\cdot3\cdot4}{3\cdot5\cdot7\cdot9}+\cdots\ ,
\]
\[
\frac{2\pi}{3\sqrt{3}}=1 +\frac{1}{2\cdot3} +\frac{1\cdot2}{3\cdot4\cdot5} +\frac{1\cdot2\cdot3}{4\cdot5\cdot6\cdot7} +\frac{1\cdot2\cdot3\cdot4}{5\cdot6\cdot7\cdot8\cdot9}+\cdots.
\]

We will now obtain some new continued fractions for $\pi.$

\section{Series with linear factors and continued fractions}

If we define:
\[
y_k = \sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)(2n+1)(2n+3)\cdots(2n+2k-3)(2n+2k-1)}
\]
then the following recurrence relation with $k \in \mathbb{N}$ is easy to establish:
\[
y_k = \frac{1}{k} \cdot y_{k-1} -\frac{1}{2k \cdot 1\cdot3\cdot5\cdots(2k-3)\cdot (2k-1)}.
\]
Indeed, we have that:
\begin{align*}
2k y_ k & = 2k \sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)(2n+1)\cdots(2n+2k-3)(2n+2k-1)}\\
& = \sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)(2n+1)\cdots(2n+2k-5)(2n+2k-3)} \\
& \hspace{4cm} -\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n+1)\cdots(2n+2k-3)(2n+2k-1)} .
\end{align*}
By shifting the index of summation, the last sum on the right can be rewritten as:
\begin{align*}
& -\sum_{n=2}^\infty \frac{(-1)^{n-1}}{(2n-1)(2n+1)\cdots(2n+2k-5)(2n+2k-3)}\\
& = \frac{1}{1\cdot3\cdot5\cdots(2k-1)} -\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)(2n+1)\cdots(2n+2k-5)(2n+2k-3)}
\end{align*}
immediately leading to the expected result. See also \cite{nimbrie}.

Using this recurrence with $y_0 = \frac{\pi}{4}$ (Leibnitz's series), we get a class of series with an increasing number of factors in the denominator:
\begin{equation}
\frac{\pi}{8}-\frac{1}{3}=\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)(2n+1)(2n+3)} \ , \label{eq:a}
\end{equation}

\[
\frac{\pi}{24}-\frac{11}{90}=\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)(2n+1)(2n+3)(2n+5)} \ ,
\]

\[
\frac{\pi}{96}-\frac{2}{63}=\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)(2n+1)(2n+3)(2n+5)(2n+7)} .
\]

Following Euler's method, we obtain
\[
\frac{8}{3\pi -8}=5 +{\cfrac{1\cdot5}{6+ \cfrac{3\cdot7}{6+ \cfrac{5\cdot9}{6+ \cfrac{7\cdot11}{6+ \dotsb}}}}} \ ,
\]

\[
\frac{24}{15\pi -44}=7 +{\cfrac{1\cdot7}{8+ \cfrac{3\cdot9}{8+ \cfrac{5\cdot11}{8+ \cfrac{7\cdot13}{8+ \dotsb}}}}} \ ,
\]

\[
\frac{96}{105\pi -320}=9 +{\cfrac{1\cdot9}{10+ \cfrac{3\cdot11}{10+ \cfrac{5\cdot13}{10+ \cfrac{7\cdot15}{10+ \dotsb}}}}} .
\]

It may be noted that the numerator of the left hand side is $4(k!)$ while the multiple of $\pi$ is $(2k-1)!!,$ where $(2k-1)$ is the constant in the largest factor $2n+2k-1$ in the denominator of the series. These continued fractions are all special cases of:

\begin{theorem*}
Let $s=\displaystyle \sum_{n=1}^\infty \frac{(-1)^{n-1}}{\prod_{j=1}^k{(2n+2j-3)}},$ then we have
\[
\frac{1}{(2k-1)!! \, s}=2k+1 +{\cfrac{1\cdot(2k+1)}{2k+2+ \cfrac{3\cdot(2k+3)}{2k+2+ \cfrac{5\cdot(2k+5)}{2k+2+ \cfrac{7\cdot(2k+7)}{2k+2+ \dotsb}}}}} \ .
\]
\end{theorem*}
To illustrate how the continued fractions are derived from the corresponding series (Euler's theorem I), we give an example (see also \cite{levrie}).
We will convert the following series given in \cite[p.269, Ex.109(c)]{knopp}:
\begin{equation}
\sum_{n=1}^\infty \frac{(-1)^{n-1}}{2n(2n+1)(2n+2)}=\frac{\pi-3}{4} \label{eq:b}
\end{equation}
into
\[
\frac{1}{\pi -3}=6 +{\cfrac{1^2}{6+ \cfrac{3^2}{6+ \cfrac{5^2}{6+ \dotsb}}}} ,
\]
the ``new" continued fraction derived in 1999 by Lange \cite{lange} though it already occurs in a 1988 paper of Castellanos \cite{castellanos} who rightly ascribed it to Euler. This continued fraction is related to Brouncker's.

Let us define $z_n = (-1)^n (\frac{\pi-3}{4}-\sum_{k=1}^n \frac{(-1)^{k-1}}{2k(2k+1)(2k+2)})$. From this it follows that:
\[
z_{n}+z_{n-1} =\frac{1}{2n(2n+1)(2n+2)} , \ \mbox{and} \
z_{n+1}+z_{n} =\frac{1}{(2n+2)(2n+3)(2n+4)} .
\]
Dividing both equations, we find:
\[
\frac{z_{n+1}+z_n}{z_n+z_{n-1}}=\frac{n(n+1)(2n+1)}{(n+1)(n+2)(2n+3)}=\frac{n(2n+1)}{(n+2)(2n+3)} .
\]

We now get rid of the denominators:
\[
(2n+3)(n+2) z_{n+1}+[(2n+3)(n+2)-(2n+1) n]z_n = (2n+1) n z_{n-1}
\]
or
\[
(2n+3)(n+2) z_{n+1}+6 (n+1) z_n = (2n+1) n z_{n-1} .
\]

If we use the transformation $w_n = (n+1) z_n$, we have:
\[
(2n+3) w_{n+1}+6w_n = (2n+1) w_{n-1} \Rightarrow (2n+3)\frac{w_{n+1}}{w_n} +6=\frac{2n+1}{\displaystyle \frac{{w_n}}{w_{n-1}}}
\]
giving the relation that generates the continued fraction on taking $n=1, 2, 3, \cdots$
\[
\frac{{w_n}}{w_{n-1}}=\frac{(2n+1)}{6 +(2n+3)\displaystyle \frac{w_{n+1}}{w_n}}.
\]
The value of the continued fraction obtained in this way is given by:
\[
\frac{w_1}{w_0} = \frac{2z_1}{z_0}=2 \cdot \frac{- (\frac{\pi-3}{4}-\frac{1}{2\cdot 3\cdot 4})}{\frac{\pi-3}{4}} .
\]
Some manipulations lead to the desired form.
Other continued fractions for $\pi$ can be obtained in a similar way. We derived this interesting series by combining (\ref{eq:a}) and (\ref{eq:b}):

\begin{equation}
\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)2n(2n+1)(2n+2)(2n+3)}=\frac{10-3\pi}{72}
\end{equation}
which gives the second convergent for $\pi:$
\[
\pi=\frac{22}{7} -24\sum_{n=2}^\infty \frac{(-1)^{n}}{(2n+1)(2n+2)(2n+3)(2n+4)(2n+5)}.
\]

We converted the preceding series into the following continued fraction:

\[
\frac{6}{10 -3\pi}=10 +{\cfrac{1\cdot5}{10+ \cfrac{3\cdot7}{10+ \cfrac{5\cdot9}{10+ \cfrac{7\cdot11}{10+ \cfrac{9\cdot13}{10+ \dotsb}}}}}}
\]
which may be compared with the continued fraction derived earlier (with $a_k=10, \; \forall \; k\in \mathbb{N}$) and one given by Euler in \cite[\S12]{euler745} and also by Osler \cite{osler}:

\[
\frac{16}{\pi}=5 +{\cfrac{1^2}{10+ \cfrac{3^2}{10+ \cfrac{5^2}{10+ \cfrac{7^2}{10+ \cfrac{9^2}{10+ \dotsb}}}}}}
\]

Osler \cite{osler} (also in \cite[p.226, eq(28a)]{perron}) gives the following two classes of continued fractions related to Brouncker's. They were already known to John Wallis \cite{wallis}. Let
\[
P_0=1, \ \mbox{and} \ P_n=\prod_{k=1}^n \frac{(2k-1)(2k+1)}{(2k)^2}.
\]
Then for $n=0, 1, 2, 3, \cdots$ we have:
\begin{gather*}
(4n+1) +\cfrac{1^2}{2(4n+1)+ \cfrac{3^2}{2(4n+1)+ \cfrac{5^2}{2(4n+1)+ \dotsb}}}=\frac{(2n+1)}{P_n}\frac{4}{\pi},\\
(4n+3) +\cfrac{1^2}{2(4n+3)+ \cfrac{3^2}{2(4n+3)+ \cfrac{5^2}{2(4n+3)+ \dotsb}}}=(2n+1)P_{n} \pi.
\end{gather*}

These formulas are special cases of the known formula
\[
\frac{4 \, \displaystyle \Gamma\left(\frac{x+3 +y}{4}\right) \Gamma\left(\frac{x+3 -y}{4}\right)}{\displaystyle \Gamma\left(\frac{x+1 +y}{4}\right) \Gamma\left(\frac{x+1- y}{4}\right)}=x +\cfrac{1^2 -y^2}{2x+ \cfrac{3^2 -y^2}{2x+ \cfrac{5^2 -y^2}{2x+ \dotsb}}}
\]
valid for either $y$ an odd integer and $x$ any complex number or $y$ any complex number and $\Re(x)>0.$ Originally due to Euler \cite[\S67]{euler123}, it occurs in an inverted form in \cite[p.140]{bernt2}.

\section{Series with quadratic factors and continued fractions}

Ramanujan's Notebook II \cite[pp.151, 153]{bernt2} contains these two continued fractions for  Catalan's constant $\displaystyle G=\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)^2}=0.915965594177\dots$:

\[
2G=2 -{\cfrac{1}{3+ \cfrac{2^2}{1+ \cfrac{2^2}{3+ \cfrac{4^2}{1+ \cfrac{4^2}{3+ \dotsb}}}}}}
\]
and
\[
2G=1 +{\cfrac{1^2}{\frac{1}{2}+ \cfrac{1\cdot2}{\frac{1}{2}+ \cfrac{2^2}{\frac{1}{2}+ \cfrac{2\cdot3}{\frac{1}{2}+ \cfrac{3^2}{\frac{1}{2}+ \dotsb}}}}}}
\]

We find in Ramanujan's Manuscript Book 1 \cite[Ch.XIV, p.107, entry 14]{Ramanujan1} (and in \cite[p.123, entry 16]{bernt2}) this formula with $m, n \notin \mathbb{Z^{-}}$:

\[
\sum_{k=1}^\infty \frac{(-1)^{k+1}}{(m+k)(n+k)}=\cfrac{1}{m+n+1+mn+ \cfrac{(m+1)^2 (n+1)^2}{m+n+3+ \cfrac{(m+2)^2 (n+2)^2}{m+n+5+ \cfrac{(m+3)^2 (n+3)^2}{m+n+7+ \dotsb}}}}
\]

Setting $m=n=-\frac{1}{2}$ and doing a little manipulation, the last formula yields a fraction given in \cite{bowman}:
\begin{equation}
G=\cfrac{1}{1+ \cfrac{1^4}{8+ \cfrac{3^4}{16+ \cfrac{5^4} {24+ \dotsb}}}} .\label{eq:bowman}
\end{equation}
Using a generalisation of the series defining $G$, we will obtain new continued fractions for this constant and for $\pi$. In the \textit{Appendix} we prove that $y_k$ defined by
\[
y_k = \sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)^2(2n+1)^2(2n+3)^2 \cdots (2n+2k-3)^2}
\]
satisfies the recurrence
\[
y_k=\frac{10k^2+8k+1}{2(2k+1)!!^2} -4k(k+1)^3 y_{k+2} .
\]
Since $y_1 = G$ and $y_2 = \frac{1}{2} - \frac{\pi}{8}$ (which can be proved using a partial fraction expansion and telescoping), the previous recurrence gives for $y_4$, $y_6$ and $y_8$:

\[
\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)^2(2n+1)^2(2n+3)^2(2n+5)^2}=\frac{2!^3}{4!^3 1!}\pi -\frac{7}{4050}
\]

\[
\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)^2(2n+1)^2\cdots(2n+9)^2}=-\frac{3!^3}{6!^3 2!}\pi +\frac{41}{44651250}
\]

\[
\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)^2(2n+1)^2\cdots(2n+13)^2}=\frac{4!^3}{8!^3 3!}\pi -\frac{14789}{134221791453750}
\]
and for $y_3$, $y_5$ and $y_7$:
\[
\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)^2(2n+1)^2(2n+3)^3}=-\frac{1!!^3 1!}{2!^4 2^1} G + \frac {19}{576}
\]
\[
\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)^2(2n+1)^2(2n+3)^2(2n+5)^2(2n+7)^3}=\frac{3!!^3 2!}{4!^4 2^2} G - \frac {3919}{108380160}
\]
\[
\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)^2(2n+1)^2\cdots(2n+11)^2}=-\frac{5!!^3 3!}{6!^4 2^3} G + \frac {22133579}{2549361475584000} .
\]
Note that an \emph{odd number} of consecutive odd squares in the denominator gives a series for $G$ while an \emph{even number} of these squares gives a series for $\pi.$

In general we found the following sums for the series related to $\pi$: $(k=0, 1, 2, \cdots)$
\[
\sum_{n=1-k}^\infty \frac{(-1)^{n-1}}{\prod_{j=0}^{2k+1} (2n+2j-1)^2}=(-1)^{k+1} \left(\frac{(k+1)!^3}{(2k+2)!^3 k!} \pi - \frac{1}{2(2k+1)!!^4}\right).
\]

For the series above the corresponding continued fractions are given in the following theorem:
\begin{theorem*}
Let
\[S
=\displaystyle \sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)^2(2n+1)^2(2n+3)^2 \cdots (2n+(2k-1))^2}.\]
Then,
\[
 \frac{1}{(2k-1)!!^2 \, S}=(2k+1)^2 +\cfrac{1^2 (2k+1)^2}{(2k+3)^2-1^2+ \cfrac{3^2 (2k+3)^2}{(2k+5)^2-3^2+ \cfrac{5^2 (2k+5)^2}{(2k+7)^2-5^2+ \dotsb}}} .
\]
\end{theorem*}

Using the above formulas and the theorem we get:

\subsection{Continued fractions for $\pi$}

The values $k=1, \, 3, \, 5$ give

\[
\frac{8}{4 -\pi}=3^2+ \cfrac{1^2 3^2}{(5^2-1^2)+ \cfrac{3^2 5^2}{(7^2-3^2)+ \cfrac{5^2 7^2}{(9^2-5^2)+ \cfrac{7^2 9^2}{(11^2-7^2)+ \dotsb}}}}\ ,
\]
\[
\frac{576}{75\pi -224}=7^2+ \cfrac{1^2 7^2}{(9^2-1^2)+ \cfrac{3^2 9^2}{(11^2-3^2)+ \cfrac{5^2 11^2}{(13^2-5^2)+ \cfrac{7^2 13^2}{(11^2-7^2)+ \dotsb}}}}
\]
and
\[
\frac{25600}{20992 -6615\pi}=11^2+ \cfrac{1^2 11^2}{(13^2-1^2)+ \cfrac{3^2 13^2}{(15^2-3^2)+ \cfrac{5^2 15^2}{(17^2-5^2)+ \cfrac{7^2 17^2}{(19^2-7^2)+ \dotsb}}}} .
\]

\subsection{Continued fractions for G}

The values $k=0, \, 2, \, 4$ give

\[
\frac{1}{G}=1^2+ \cfrac{1^4}{(3^2-1^2)+ \cfrac{3^4}{(5^2-3^2)+ \cfrac{5^4}{(7^2-5^2)+ \cfrac{7^4}{(9^2-7^2)+ \dotsb}}}}\ ,
\]
(this is the same one as (\ref{eq:bowman}))

\[
\frac{64}{19 -18G}=5^2+ \cfrac{1^2 5^2}{(7^2-1^2)+ \cfrac{3^2 7^2}{(9^2-3^2)+ \cfrac{5^2 9^2}{(11^2-5^2)+ \cfrac{7^2 11^2}{(13^2-7^2)+ \dotsb}}}}
\]
and
\[
\frac{49152}{22050G -19595}=9^2+ \cfrac{1^2 9^2}{(11^2-1^2)+ \cfrac{3^2 11^2}{(13^2-3^2)+ \cfrac{5^2 13^2}{(15^2-5^2)+ \cfrac{7^2 15^2}{(17^2-7^2)+ \dotsb}}}}.
\]

\section{Other continued fractions for $G$}
While looking for continued fractions for $G$ related to (\ref{eq:bowman}) in the same way that Lange's continued fraction for $\pi$ is related to Brouncker's, we found these:
\begin{equation}
\frac{2^5}{6G -1}=7+ 3\cdot \cfrac{1^4}{3(3^2-1^2)+ \cfrac{3^4}{3(5^2-3^2)+ \cfrac{5^4}{3(7^2-5^2)+ \cfrac{7^4}{3(9^2-7^2)+ \dotsb}}}} \label{eq:G}
\end{equation}
\[
\frac{2^{13}}{82G -19}=145+ 41\cdot \cfrac{1^4}{5(3^2-1^2)+ \cfrac{3^4}{5(5^2-3^2)+ \cfrac{5^4}{5(7^2-5^2)+ \cfrac{7^4}{5(9^2-7^2)+ \dotsb}}}}
\]
\[
\frac{2^{17}}{882G -\frac{713}{3}}=229+ 49\cdot \cfrac{1^4}{7(3^2-1^2)+ \cfrac{3^4}{7(5^2-3^2)+ \cfrac{5^4}{7(7^2-5^2)+ \cfrac{7^4}{7(9^2-7^2)+ \dotsb}}}}
\]
and it goes on like this. We'll prove (\ref{eq:G}). This continued fraction is related to the following series:
\[
\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(4(n-1)^2+3)(4n^2+3)(2n-1)^2}
\]
as can be checked using Euler's method to transform a series into a continued fraction. Now the sum of this series can be calculated like this:
\begin{align*}
& \sum_{n=1}^\infty \frac{(-1)^{n-1}}{(4(n-1)^2+3)(4n^2+3)(2n-1)^2} \\
& = \sum_{n=1}^\infty {(-1)^{n-1}}\left(-\frac{1}{32 (4 (n-1)^2+3)}-\frac{1}{32 (4n^2+3)}+\frac{1}{16(2n-1)^2} \right)\\
& = -\frac{1}{32} \sum_{n=1}^\infty \frac{(-1)^{n-1}}{4 (n-1)^2+3} -\frac{1}{32} \sum_{n=1}^\infty \frac{(-1)^{n-1}}{4n^2+3} + \frac{1}{16}\sum_{n=1}^\infty \frac{(-1)^{n-1}}{(2n-1)^2}\\
& = -\frac{1}{3\cdot 32}+\frac{1}{16}G
\end{align*}
as a consequence of telescoping and the definition of $G$.\\
The series needed to prove the value of the other 2 continued fractions for $G$ are:
\[
\sum_{n=1}^\infty \frac{(-1)^{n-1}}{r(n-1) r(n) (2n-1)^2}
\]
with $r(n) = 16n^4+88n^2+41$ and $r(n)=64n^6+1168n^4+3628n^2+1323$ respectively.


text:
In general, the expectations involved in \eqref{eq:SoftConvex} cannot be computed analytically. Thus, we need to resort to sampling based approaches in order to solve these problems. This has been studied extensively in recent years in the context of machine learning, where stochastic gradient methods and variants have been shown to be efficient, particularly in the context training machine learning algorithms with huge amounts of data. We now describe stochastic gradient methods for solving \eqref{eq:SoftConvex} and adapt the convergence guarantees of stochastic gradient methods to our setting.

\subsection{Stochastic Gradient Methods with Convergence Guarantees}\label{sec:GradUpdate}
In this section, we will derive gradients of the convex objective function  \eqref{eq:SoftConvex}.  We will assume that the function $f$ is differentiable at all $\theta \in \R^k$. In order to get unbiased gradient estimates, we exponentiate the objective \eqref{eq:SoftConvex} to get:
\[\Gb\br{\theta}=\ExP{\noise \sim \Gauss{0}{\Sigma}}{\expb{\risk f\br{\theta+\noise}+\frac{\risk}{2}\tran{\theta}R\theta}}.\]
Since $f\br{\theta}$ is differentiable for all $\theta$, so is $\expb{\risk f\br{\theta+\noise}+\frac{1}{2}\tran{\theta}\br{\risk R}\theta}$. Further, suppose that
\[\ExP{}{\expb{2\risk f\br{\theta+\noise}}\norm{ \nabla f\br{\theta+\noise} + R \theta}^2}\]
exists and is finite for each $\theta$. Then, if we differentiate $\Gb\br{\theta}$ with respect to $\theta$, we can interchange the expectation and differentiation to get
\begin{align*}
\ExP{}{\expb{\risk f\br{\theta+\noise}+\frac{\risk}{2}\tran{\theta}R\theta}\br{\risk \nabla f\br{\theta+\noise}+\risk R \theta}}
\end{align*}
Thus, we can sample $\noise \sim \Gauss{0}{\Sigma}$ and get an unbiased estimate of the gradient
\begin{align}
\risk \expb{\risk f\br{\theta+\noise}+\frac{\risk}{2}\tran{\theta}R\theta}\br{\nabla f\br{\theta+\noise}+ R \theta}\label{eq:UnbiasedGrad1}
\end{align}
which we denote by $\hat{\nabla}\Gb\br{\theta,\noise}$. As in standard stochastic gradient methods, one saves on the complexity of a single iteration by using a single (or a small number of) samples to get a gradient estimate while still converging to the global optimum with high probability and in expectation, because over multiple iterations one moves along the negative gradient ``on average''.
\begin{algorithm}[Stochastic Gradient Method]
 \begin{algorithmic}
 \State $\theta \gets 0$
\For {$i = 1,\ldots,T$}
    \State $\noise \sim \Gauss{0}{\Sigma},\theta \gets  \Proj{\C}{\theta-\eta_i\hat{\nabla}\Gb\br{\theta,\noise}}$
\EndFor
\end{algorithmic}
 \caption{Stochastic Gradient Method for \eqref{eq:SoftConvex}}\label{alg:StochGrad}
\end{algorithm}

\begin{comment}
\begin{algorithm}[Derivative-Free Stochastic Gradient Method]
 \begin{algorithmic}
 \State $\theta \gets 0$
\For {$i = 1,\ldots,T$}
    \State $\ygu \sim \Gauss{\theta}{\Sigma}$
    \State $\delta \theta \gets \Gcal\br{\ygu}\expb{\frac{1}{2}\tran{\theta}\br{\risk R}\theta}\br{\inve{\Sigma}\br{\ygu-\theta}+\risk R\theta}$
    \State $\theta \gets  \Proj{\C}{\theta-\eta_i\delta \theta}$
\EndFor
\end{algorithmic}
 \caption{Derivative-Free Stochastic Gradient Method for \eqref{eq:SoftConvex}}\label{alg:DFStochGrad}
\end{algorithm}
\end{comment}

From standard convergence theory for stochastic gradient methods \citep{BubeckNotes}, we have:
\begin{corollary}
Suppose that $\ExP{}{\norm{\hat{\nabla}\Gb\br{\theta,\noise}}^2}\leq \zeta^2$, \\$\C$ is contained in a ball of radius $\Rad{\C}$ and \\$ \risk R \succeq \inve{\Sigma}$. Run algorithm \ref{alg:StochGrad} for $T$ iterations with \\$\eta_i=\frac{\Rad{\C}}{\zeta}\sqrt{\frac{1}{2i}}$ and define $\hat{\theta}=\frac{1}{T}\sum_{i=1}^T \theta_i$, $\opt{\Gb}=\argmin_{\theta\in \C}\Gb\br{\theta}$. Then, we have
\[\ExP{}{\Gb\br{\hat{\theta}}}-\opt{\Gb}\leq \Rad{\C}\zeta \sqrt{\frac{1}{2T}}\]
\end{corollary}

In the following theorem, we prove a convergence rate for algorithm \ref{alg:StochGrad}.
\begin{theorem}\label{thm:AlgConv}
Suppose that $R=\kappa I, \Sigma =\sigma^2 I,\risk \kappa \succeq \frac{1}{\sigma^2}$, $\C$ is contained in a sphere of radius $\Rad{\C}$ and that for all $\theta \in \C$ $\fpert{g}{\Sigma}{\theta} \leq \mup$. Also, suppose that for some $\beta<\powb{\risk}{-1}$:
\[\frac{1}{2\risk}\logb{\ExP{\noise\sim\Gauss{0}{\Sigma}}{\expb{2\frac{\risk\sigma^2}{\beta}\norm{\nabla f\br{\theta+\noise}}^2}}}\leq \gamma^2\]
Define $\delta=\sqrt{\frac{\beta \gamma^2 }{\sigma^2\br{1-\risk\beta}}}+\kappa\Rad{\C}$. Then, we can choose
$\zeta \leq \risk^2\delta^2\expb{2\risk(\mup+\gamma^2)+\frac{\risk\beta}{1-\risk\beta}-\sigma^2\kappa}$.
\end{theorem}

\begin{remark}
The convergence guarantees are in terms of the exponentiated objective $\Gb\br{\theta}$. We can convert these into bounds on $\logb{\Gb\br{\theta}}$ as follows:
\begin{align*}
& \ExP{}{\logb{\Gb\br{\hat{\theta}}}} \leq \logb{\ExP{}{\Gb\br{\hat{\theta}}}}  \leq \logb{\opt{\Gb}+\frac{\zeta\Rad{\C}}{\sqrt{2T}}}
\end{align*}
where the first inequality follows from concavity of the $\log$ function. Subtracting $\logb{\opt{\Gb}}$, we get
\begin{align*}
&\ExP{}{\logb{\Gb\br{\hat{\theta}}}}-\logb{\opt{\Gb}} \leq \logb{1+\frac{\zeta\Rad{\C}}{\sqrt{2T}\opt{\Gb}}}.
\end{align*}
\end{remark}

\begin{comment}
\subsection{Heuristic Approach in Log-Space}\label{sec:AlgLog}

A similar approach can be taken in the log-space where one optimizes $\frac{1}{\risk}\logb{\Gb\br{\theta}}$. We have proved that this is a convex optimization problem (theorem \ref{thm:ConvexityMain}) and further, this leads to numerically better behaved algorithms since the exponentiation can blow up or diminish values numerically. However, the problem is that it is not simple to get an unbiased estimate of the gradient of $\frac{1}{\risk}\logb{\Gb\br{\theta}}$ objective:
\[\nabla \frac{1}{\risk}\logb{\Gb\br{\theta}}=\frac{\ExP{}{\expb{\risk f\br{\theta+\noise}}}\br{ \nabla f\br{\theta+\noise}}}{\ExP{}{\expb{\risk f\br{\theta+\noise}}}}+ R\theta\]
Since the gradient is a ratio of expectations, taking a single instantiation of the noise $\noise \sim \Gauss{0}{\Sigma}$ does not give an unbiased estimate of the gradient. Instead, one must sample from a reweighted distribution given by
\[Q\br{\noise} = \frac{P\br{\noise}\expb{\risk f\br{\theta+\noise}}}{\ExP{\noise \sim P}{\expb{\risk f\br{\theta+\noise}}}}\]
where $P$ is the Gaussian density of $\Gauss{0}{\Sigma}$. This distribution in general is difficult to compute or sample from.

However, one can take a heuristic approach of approximate the numerator and denominator separately (akin to importance sampling) to compute
\[\frac{\sum_{k=1}^{M} \expb{\risk f\br{\theta+\noise^{(k)}}}\br{ \nabla f\br{\theta+\noise^{(k)}}}}{\sum_{k=1}^{M} \expb{\risk f\br{\theta+\noise^{(k)}}}}+ R\theta\]
where $\noise^{(k)} \sim \Gauss{0}{\Sigma}, k=1,\ldots,M$. As $M \to \infty$, this converges to the true gradient (it is consistent estimator) although for finite $M$ it is not an unbiased gradient estimate. However, one can use this in conjunction with a standard stochastic gradient method - it would enjoy better numerical stability and smaller variance since random values are not exponentiated but normalized while forming gradient estimates. In recent work \citep{StochasticProximal}, the authors study stochastic versions of proximal gradient algorithms and arrive at stronger convergence guarantees under appropriate assumptions, even with biased gradient estimates. Investigating how these results apply in our setting is a direction for future work.

\begin{remark}[Low-Rank Functions]

There are some special cases where one does not need to resort to sampling at all in order to solve optimization problems arising from our convexification. For example, in the supervised learning application described in section \ref{sec:Supervised}, the noise added is a scalar, and even if the expectation is computed by numerical integration, it can be evaluated efficiently. This can be generalized to the general case of \emph{low-rank functions} studied in \citep{mittal2013fptas}. A function $g$ is a said to have rank $k$ if it can be written as:
\[g\br{\theta}=h\br{\tran{a_1}\theta,\tran{a_2}\theta,\ldots,\tran{a_k}\theta}\]
In our setting, this means that we can convexify the problem by adding only $k$-dimensional noise. If $h\br{\ygu}=f\br{\ygu}+\frac{1}{2}\tran{\ygu}R\ygu$:
\begin{align*}
&\ExP{\noise_i \sim \Gauss{0}{\sigma_i}}{\expb{\risk f\br{\tran{a_1}\theta+\noise_1,\ldots,\tran{a_k}\theta+\noise_k}}}\\
& \quad +\sum_{i=1}^k r_i\powb{\tran{a_i}\theta}{2}.
\end{align*}
As long as $\risk r_i \geq \frac{1}{\sigma_i^2}$ for $1\leq i \leq k$, the resulting problem is convex and further, the optimal solution can be found using a standard deterministic optimization algorithm (like LBFGS) as long as $k$ is small (since the expectation can be computed explicitly by numerical integration).
\end{remark}
\end{comment}

config: arxiv


text:
\begin{document}

\twocolumn[
\icmltitle{Non-Convex Exact Community Recovery \\ in Stochastic Block Model}

\icmlsetsymbol{equal}{*}

\begin{icmlauthorlist}
	\icmlauthor{Peng Wang}{equal,to}
	\icmlauthor{Zirui Zhou}{equal,goo}
	\icmlauthor{Anthony Man-Cho So}{to}
\end{icmlauthorlist}

\icmlaffiliation{to}{Department of Systems Engineering and Engineering Management, The Chinese University of Hong Kong, Shatin, NT, Hong Kong}
\icmlaffiliation{goo}{Huawei Technologies Canada Co., Ltd., Burnaby, Canada}

\icmlcorrespondingauthor{Zirui Zhou}{zirui.zhou@huawei.com}

\icmlkeywords{Machine Learning, ICML}

\vskip 0.3in
]

\printAffiliationsAndNotice{\icmlEqualContribution}

\begin{abstract}
Learning community structures in graphs that are randomly generated by stochastic block models (SBMs) has received much attention lately. In this paper, we focus on the problem of exactly recovering the communities in a binary symmetric SBM, where a graph of $n$ vertices is partitioned into two equal-sized communities and the vertices are connected with probability $p = \alpha\log(n)/n$ within communities and $q = \beta\log(n)/n$ across communities for some $\alpha>\beta>0$. We propose a two-stage iterative algorithm for solving this problem, which employs the power method with a random starting point in the first stage and turns to a generalized power method that can identify the communities in a finite number of iterations in the second stage. It is shown that for any fixed $\alpha$ and $\beta$ such that $\sqrt{\alpha} - \sqrt{\beta} > \sqrt{2}$, which is known to be the information-theoretic limit for exact recovery, the proposed algorithm exactly identifies the underlying communities in $\tilde{O}(n)$ running time with probability tending to one as $n\rightarrow\infty$. We also present numerical results of the proposed algorithm to support and complement our theoretical development.
\end{abstract}

\section{Introduction}\label{sec:intro}
Learning community structures in graphs is a fundamental task in machine learning and computer science. It has found broad applications in physics \cite{fortunato2010community, newman2004finding}, biology \cite{cline2007integration}, and social science \cite{girvan2002community}, to name a few. In research on community estimation, the stochastic block model (SBM), which is a generative model for random graphs that admits community structures, has become increasingly popular as a platform for validating theoretical ideas and comparing numerical algorithms. In particular, substantial advances have been made in the past decade on understanding the statistical limits for detection and recovery of communities in graphs that are generated by SBMs, and on developing computationally tractable methods that can detect and recover the underlying communities at their corresponding statistical limits; see, for example, \citet{mossel2014consistency,abbe2015community,abbe2016exact,abbe2017entrywise,chen2016statistical,gao2017achieving}.

In this paper, we focus on the problem of exactly recovering the communities in the binary symmetric SBM (also known as the planted bisection model). Specifically, given $n$ nodes that are partitioned into two equal-sized clusters, a graph on these $n$ nodes is randomly generated such that each pair of nodes are connected with probability $p$ if they are in the same cluster and with probability $q$ if not, where $p>q>0$. Then, we aim to recover the underlying clusters using only the adjacency matrix of the generated graph. For this problem, the regime where $p = \alpha\log(n)/n$ and $q = \beta\log(n)/n$ for some $\alpha>\beta>0$ is of particular interest as it possesses a sharp threshold for exact recovery. Indeed, from an information-theoretic point of view, it has been proved that recovering the clusters with high probability is impossible if $\sqrt{\alpha} - \sqrt{\beta}<\sqrt{2}$ and is possible if $\sqrt{\alpha} - \sqrt{\beta} > \sqrt{2}$~\cite{abbe2016exact,mossel2014consistency}. Moreover, a number of computationally tractable algorithms have been shown to achieve this threshold for exact recovery~\cite{abbe2016exact,abbe2017entrywise,hajek2016achieving,gao2017achieving}. More precisely, these algorithms have running time that is polynomial in $n$ and for any fixed $\alpha$ and $\beta$ such that $\sqrt{\alpha} - \sqrt{\beta} > \sqrt{2}$, they successfully identify the underlying clusters with probability tending to one as $n\rightarrow\infty$. However, these algorithms rely on solving a semidefinite programming (SDP) or eigenvector computation, whose time complexities are polynomial but generally not linear in $n$. A major goal of this paper is to propose an algorithm that runs in time nearly linear in $n$ and achieves exact recovery at the information-theoretic limit.

\subsection{Related Works}
For the binary symmetric SBM, whether it is possible to recover the underlying communities with high probability depends on the scaling of $p$, $q$, and $p-q$. When $p = a/n$ and $q = b/n$ for some $a>b>0$ (also known as the sparse regime), it is impossible to recover the communities because the graph is unconnected with high probability \cite{decelle2011asymptotic}. In the regime where $p = \alpha\log(n)/n$ and $q = \beta\log(n)/n$ for some $\alpha>\beta>0$, \citet{abbe2016exact} and \citet{mossel2014consistency} independently showed that there exists a sharp threshold for exact recovery. Specifically, they proved that from an information-theoretic point of view, exact recovery is impossible if $\sqrt{\alpha} - \sqrt{\beta}<\sqrt{2}$, while it is possible and can be achieved by the maximum likelihood (ML) estimator if $\sqrt{\alpha} - \sqrt{\beta} > \sqrt{2}$. Computing the ML estimator is, however, computationally intractable in general. Indeed, by encoding the underlying clusters into a vector $\bxs\in\{\pm1\}^n$ such that $x_i^* = 1$ if node $i$ is in the first cluster and $x_i^* = -1$ otherwise, the ML estimator of $\bxs$ is the solution of the following problem:
\begin{equation}
\label{eq:MLE}
\max \left\{\bx^T\bA\bx: \, \mathbf{1}_n^T\bm{x}=0, \, x_i = \pm 1, \, i = 1,\dots,n \right\},
\end{equation}
where $\bA$ is the adjacency matrix of the graph and $\bo_n$ is the all-one vector of dimension $n$. Problem \eqref{eq:MLE} is equivalent to the problem of finding minimum bisection of a graph, which is known to be NP-hard in the worst case \cite{garey1974some}. Nevertheless, the above threshold provides the information-theoretic limit for exact recovery in the binary symmetric SBM, which is useful for benchmarking recovery algorithms.

Over the past decades, many computationally tractable algorithms have been proposed that can be applied to the problem of exact community recovery in the binary symmetric SBM. Some of them are proposed as heuristic combinatorial algorithms for the minimum bisection problem with an average-case performance guarantee and some others tackle the community recovery problem based on convex or non-convex relaxations of the ML estimation problem \eqref{eq:MLE}. In the comparison of these algorithms, of particular interest are the following two factors: (i) computational efficiency and (ii) conditions on $p$ and $q$ that guarantee exact recovery. Moreover, in light of the above discussion, a computationally efficient algorithm that can achieve exact recovery at the information-theoretic limit would be highly desirable. Prior to the work \citet{abbe2016exact}, there were two polynomial-time algorithms that allow for exact recovery in the regime $p = \alpha\log(n)/n$ and $q = \beta\log(n)/n$, which are respectively proposed in \citet{boppana1987eigenvalues} and \citet{mcsherry2001spectral}. However, \citet{boppana1987eigenvalues} needs $(\alpha - \beta)^2/(\alpha+\beta) > 72$ and \citet{mcsherry2001spectral} needs $(\alpha - \beta)^2/(\alpha+\beta) > 64$ for exact recovery, both of which do not meet the information-theoretic limit $\sqrt{\alpha} - \sqrt{\beta} > \sqrt{2}$. The first computationally tractable algorithm that can achieve exact recovery at this limit was given in \citet{abbe2016exact}, which employs the polynomial-time algorithm in \citet{massoulie2014community} for obtaining a partial recovery, followed by a local improvement procedure. Later, \citet{hajek2016achieving} and \citet{bandeira2018random} independently proved that the solution of a semidefinite relaxation of Problem \eqref{eq:MLE} can identify the underlying communities exactly with high probability if $\sqrt{\alpha} - \sqrt{\beta} > \sqrt{2}$, confirming a conjecture raised in \citet{abbe2016exact}. As SDPs can be solved in polynomial time, this also provides a computationally tractable algorithm that achieves exact recovery at the information-theoretic limit. More recently, by conducting an entrywise eigenvector analysis of random matrices with low expected rank, \citet{abbe2017entrywise} proved that the vanilla eigenvector-based algorithm achieves exact recovery at the information-theoretic limit. The dominating computational cost in this algorithm is to compute the eigenvector associated with the second largest eigenvalue of the adjacency matrix $\bA$, which is known to be polynomial time computable.

We would also like to mention some algorithms for the considered problem with extremely low computational complexity. \citet{condon2001algorithms} proposed a simple combinatorial algorithm for minimum bisection problem, whose per iteration cost is linear in $n$. To achieve exact community recovery with high probability, it needs $p - q\geq n^{-1/2 + \epsilon}$ with some $\epsilon>0$. Recently, \citet{bandeira2016low} proposed to apply the Burer-Monteiro decomposition~\cite{burer2003nonlinear} to the semidefinite relaxation of Problem \eqref{eq:MLE}, which results in solving an optimization problem defined on a smooth manifold. They showed that all second-order stationary points of the problem, which can be computed efficiently by the Riemannian trust-region method~\citep{boumal2018global}, correspond to the underlying communities with high probability as long as $(p-q)/\sqrt{p+q}\geq cn^{-1/6}$ for some constant $c>0$. Despite their low computational complexity, both approaches require much stronger conditions on $p$ and $q$ for ensuring exact recovery. In particular, they cannot even allow for exact recovery in the regime $p = \alpha\log(n)/n$ and $q = \beta\log(n)/n$ for some $\alpha>\beta>0$.

\subsection{Our Contribution}\label{sec:oc}
In this work, we propose a two-stage iterative algorithm for exact community recovery in the binary symmetric SBM. Our algorithm is based on the following regularized version of the ML estimation problem \eqref{eq:MLE}:
\begin{equation}
\label{eq:rMLE}
\max \left\{\bx^T\bB\bx: \, x_i = \pm 1, \, i = 1,\dots,n \right\},
\end{equation}
where $\bB = \bA - \rho\bm{E}_n$ and $\rho = \bo_n^T\bA\bo_n/n^2$. Specifically, the proposed iterative algorithm starts with a vector that is chosen randomly with uniform distribution over the unit sphere of $\mathbb{R}^n$, employs the standard power method (PM) in the first stage for approximating the dominant eigenvector of $\bB$, and then applies a generalized power method (GPM) to Problem \eqref{eq:rMLE} in the second stage that can exactly identify the underlying communities in a finite number of iterations. We show that in the regime $p = \alpha\log(n)/n$ and $q = \beta\log(n)/n$ for some $\alpha>\beta>0$, the proposed algorithm achieves exact recovery at the information-theoretic limit and enjoys a time complexity that scales nearly linear in $n$. More precisely, for any fixed $\alpha$ and $\beta$ such that $\sqrt{\alpha} - \sqrt{\beta} > \sqrt{2}$, we prove that the following event happens with probability tending to one as $n\rightarrow\infty$: the proposed algorithm terminates at the ground truth $\bxs$ or $-\bxs$\footnote{Note that $\bx^*$ and $-\bx^*$ represent the same community structure up to a global flip of the labels.} in $O(\log n/\log\log n)$ iterations of the PM and $O(\log n/\log\log n)$ iterations of the GPM, with $O(n\log n)$ computational complexity in each iteration of the PM and the GPM. At the heart of our result is a careful analysis on the performance of the PM and the GPM when applied to our problem. In particular, we show the following results hold with high probability:
\begin{itemize}
\item[(i)] The sequence of iterates generated by the first stage of our algorithm converges linearly to a dominant eigenvector of $\bB$ and the ratio in defining the linear rate of convergence tends to $0$ as $n\rightarrow\infty$.
\item[(ii)] The iterates generated by the second stage of our algorithm possess a contraction property once their Euclidean distance to the ground truth is no greater than a threshold that scales linearly with $\sqrt{n/\log n}$.
\item[(iii)] The second stage of our algorithm exhibits finite termination as it admits a one-step convergence property once the generated iterate has Euclidean distance no greater than a constant threshold to the ground truth.
\end{itemize}
We also conduct experiments on synthetic data sets and compare the performance of the proposed algorithm to some existing ones. Numerical results support our theoretical development and also show the advantage of the proposed algorithm in terms of computational efficiency.

The rest of this paper is organized as follows. In section \ref{sec:main-result}, we introduce the proposed two-stage algorithm for exact community recovery and present the main results of this paper regarding its recovery performance and computational efficiency. In sections \ref{sec:pm} and \ref{sec:gpm}, we respectively analyze the properties of the power method and the generalized power method that are respectively used in the two stages of the proposed algorithm. In section \ref{sec:pf-main}, we provide the proofs of the main results. Numerical results of the proposed algorithm are reported in section \ref{sec:num}.

\emph{Notation.}  Let $\R^n$ be the $n$-dimensional Euclidean space and $\|\cdot\|_2$ be the Euclidean norm. We write matrices in capital bold letters like $\bA$, vectors in bold lower case like $\bm{a}$, and scalars as plain letters. We denote by $a_{ij}$ the $(i,j)$-th element of $\bm{A}$ and $x_i$ the $i$-th element of $\bx$. We use $\bo_n$ and $\bm{E}_n$ to denote the all-one vector and all-one matrix of dimension $n$, respectively. Given a vector $\bd$, we denote by $\Diag(\bd)$ the diagonal matrix with $\bd$ on its diagonal. For any positive integer $n$, let $[n]$ denote the set $\{1,\ldots,n\}$; for any discrete set $T$, let $|T|$ denote the number of elements in $T$. We use $\mathbf{Bern}(p)$ to denote the Bernoulli random variable with mean $p$. For any $\bv\in\mathbb{R}^n$, $\bv/|\bv|$ denotes the vector of $\mathbb{R}^n$ defined as
\begin{equation}
\label{eq:frac}
\left(\frac{\bv}{|\bv|}\right)_i = \left\{ \begin{array}{ll}
1, \quad \ \ \;\mbox{if} \ v_i \geq 0, \\
-1, \quad \mbox{otherwise},
\end{array}\right. \quad i = 1,\dots,n.
\end{equation}

\section{Main Results}\label{sec:main-result}
Given $n$ nodes that are partitioned into two equal-sized clusters, we denote by $\bx^*$ the vector that encodes its true community structure, {\it e.g.}, for every $i\in[n]$, $x_i^* = 1$ if the node $i$ belongs to the first cluster and $x_i^*=-1$ if it belongs to the second one. Given $\bx^*$, the binary symmetric SBM generates an undirected graph of $n$ vertices, whose adjacency matrix $\bA$ is randomly generated as follows.
\begin{model}\label{model:SBM}
	The elements $\{a_{ij}: 1\leq i \leq j\leq n\}$ of $\bA$ are generated independently by
	\begin{align*}\label{SBM}
	a_{ij}\sim \left\{
	\begin{aligned}
	\mathbf{Bern}(p),\quad & \text{if}\ \ x_i^*x_j^*=1,\\
	\mathbf{Bern}(q),\quad & \text{if}\ \ x_i^*x_j^*=-1,
	\end{aligned}
	\right.
	\end{align*}
	where
	$$ p = \frac{\alpha\log n}{n} \quad \mbox{and} \quad q = \frac{\beta\log n}{n} $$
	for some constants $\alpha>\beta>0$. Besides, we have $a_{ij} = a_{ji}$ for all $1\leq j<i\leq n$.
\end{model}
Then, the problem of community recovery is to develop efficient methods that can find $\bx^*$ or $-\bx^*$ with high probability, given the adjacency matrix $\bA$ generated from Model~\ref{model:SBM}. Our proposal for tackling this problem is a two-stage algorithm summarized in Algorithm \ref{alg:TSA}. It starts with a vector $\by^0$ that is chosen randomly with uniform distribution over the unit sphere of $\mathbb{R}^n$, employs the standard power iteration $N$ times in the first stage (lines 4-6 of Algorithm \ref{alg:TSA}), and then applies a generalized power iteration to Problem \eqref{eq:rMLE} in the second stage until a fixed point of the generalized power iteration is reached (lines 8--13 of Algorithm \ref{alg:TSA}). We use the name generalized power iteration due to the fact that the operation $\bB\bx^{k-1}/|\bB\bx^{k-1}|$ (see \eqref{eq:frac} for its definition) is essentially computing a projection of $\bB\bx^{k-1}$ to the feasible set of Problem \eqref{eq:rMLE}.
\begin{algorithm}[t]
	\caption{A Two-Stage Algorithm for Exact Recovery}
	\begin{algorithmic}[1]
		\STATE \textbf{Input:} adjacency matrix $\bA$, positive integer $N$
		\STATE set $\rho \leftarrow \bo_n^T\bA\bo_n/n^2$ and $\bB \leftarrow \bA - \rho\bm{E}_n$
		\STATE choose $\by^0$ randomly with uniform distribution over the unit sphere of $\mathbb{R}^n$
		\FOR{$k=1,2,\ldots, N$}
		\STATE set $\by^k\leftarrow \bB\by^{k-1}/\|\bB\by^{k-1}\|_2$
		\ENDFOR
		\STATE set $\bx^0\leftarrow \sqrt{n}\by^{N}$
		\FOR{$k=1,2,\dots$}
		\STATE set $\bx^k\leftarrow \bB\bx^{k-1}/|\bB\bx^{k-1}|$
		\IF{$\bx^k = \bx^{k-1}$}
		\STATE terminate and return $\bx^k$
		\ENDIF
		\ENDFOR
	\end{algorithmic}
	\label{alg:TSA}
\end{algorithm}

We next present the main result of this paper, which shows that Algorithm \ref{alg:TSA} achieves exact recovery at the information-theoretic limit and also provides explicit iteration complexity bounds for Algorithm \ref{alg:TSA} to exactly recover the underlying communities.
\begin{thm}
\label{thm:main-thm}
Let $\bA$ be randomly generated by Model \ref{model:SBM}. If $\sqrt{\alpha} - \sqrt{\beta} > \sqrt{2}$, then the following statement holds with probability at least $1-n^{-\Omega(1)}$: Algorithm \ref{alg:TSA} finds $\bx^*$ or $-\bx^*$ in $O(\log n/\log\log n)$ power iterations and $O(\log n/\log\log n)$ generalized power iterations.
\end{thm}

Note that the computational cost of both the power iteration and the generalized power iteration is dominated by the one for computing the matrix-vector product $\bB\bv$ for some $\bv\in\mathbb{R}^n$. Since $\bB = \bA - \rho\bm{E}_n$ and $\bm{E}_n = \bo_n\bo_n^T$ is a rank-one matrix, the cost of computing $\bB\bv$ is mainly determined by the number of non-zero entries in $\bA$. Also, as $\bA$ is randomly generated by Model \ref{model:SBM}, one can show by a simple concentration argument that the number of non-zero entries in $\bA$ is, with high probability, in the order of $n\log n$. These, together with the iteration bounds provided in Theorem \ref{thm:main-thm}, imply that the time complexity of Algorithm \ref{alg:TSA} is with high probability nearly linear in $n$.
\begin{coro}
\label{coro:comp-time}
Let $\bA$ be randomly generated by Model \ref{model:SBM}. If $\sqrt{\alpha} - \sqrt{\beta} > \sqrt{2}$, then with probability at least $1-n^{-\Omega(1)}$, Algorithm \ref{alg:TSA} finds $\bx^*$ or $-\bx^*$ in $O(n\log^2 n)$ time complexity.
\end{coro}

\section{Analysis of the Power Method}\label{sec:pm}
In this section, we analyze the performance of the PM with a random initial point that is employed in the first-stage of Algorithm \ref{alg:TSA}. In particular, we shall characterize the convergence rate of the iterates generated by the PM to a dominant eigenvector of $\bB$. To begin, we present two lemmas that shall be used in studying the spectral property of $\bB$. The proof of Lemma \ref{lem:rho-bound} can be found in the appendix, while the proof of Lemma \ref{lem:A-2norm} is omitted because it follows directly from \citet[Theorem 5.2]{lei2015consistency}.
\begin{lemma}\label{lem:rho-bound}
	Suppose that $\bA$ is generated by Model \ref{model:SBM} and $\rho = \bo_n^T\bA\bo_n/n^2$. Then, it holds with probability at least $1-2n^{-\frac{1}{2(\alpha+\beta+1)}}$ that
	\be\label{result:lem-rho-bound}
	\left |\rho - \frac{p+q}{2} \right| \le \frac{\log n}{n^{3/2}}.
	\ee
\end{lemma}

\begin{lemma}\label{lem:A-2norm}
	Let $\bA$ be generated by Model	\ref{model:SBM}. Then, there exist constants $c_1\ge 1$ and $c_2>0$, whose values depend on $\alpha$ and $\beta$, such that
	\begin{equation}
	\label{eq:A-2norm}
	\|\bA - \bE[\bA]\|_2 \leq c_1\sqrt{\log n}
	\end{equation}
	holds with probability at least $1 - c_2n^{-3}$.
\end{lemma}

Equipped with Lemmas \ref{lem:rho-bound} and \ref{lem:A-2norm}, and by applying Weyl's inequality, it is not hard to obtain the following result regarding the magnitudes of eigenvalues of $\bB$. The proof of this result can be found in the appendix.

\begin{lemma}\label{lem:eig-bound}
Suppose that $\bA$ is randomly generated by Model \ref{model:SBM} and $\bB=\bA-\rho\bm{E}_n$ with $\rho=\bo_n^T\bA\bo_n/n^2$. Let $\lambda_1\ge \lambda_2 \ge\dots \ge \lambda_n$ be the eigenvalues of $\bB$. Then, for all sufficiently large $n$, it holds with probability at least $1-2n^{-\frac{1}{2(\alpha+\beta+1)}} - c_2n^{-3}$ that
\begin{equation}\label{eq:lambda1}
\lambda_1\ge \frac{\alpha-\beta}{3}\log n
\end{equation}
and
\begin{equation}\label{eq:lambdai}
 |\lambda_i| \le 2c_1\sqrt{\log n}, \quad i=2,\dots,n,
\end{equation}
where $c_1$ and $c_2$ are the constants in Lemma \ref{lem:A-2norm}.
\end{lemma}
Recall that the PM in Algorithm \ref{alg:TSA} starts with a random initial point. To characterize its convergence rate, we need the following result on the uniform distribution over the unit sphere of $\mathbb{R}^n$, whose proof follows similar arguments as those in \citet{KW92} and is provided in the appendix.
\begin{lemma}\label{lem:unit-sphere}
	Suppose that $\bb\in\mathbb{R}^n$ is randomly generated from the uniform distribution over the unit sphere of $\mathbb{R}^n$. Then, for all $n\geq 8$, it holds that
	\begin{equation}
	\label{eq:unif-sphere}
	\P \left( \sum_{i=2}^n \left(\frac{b_i}{b_1} \right)^2 \leq \frac{n^2}{2} \right) \geq 1 - 2n^{-1/2}.
	\end{equation}
\end{lemma}
Now we are ready to characterize the convergence rate of the PM used in the first-stage of Algorithm \ref{alg:TSA}.
\begin{prop}\label{prop:rate-power}
Suppose that $\bA$ is randomly generated by Model \ref{model:SBM}. Let $\{\by^k\}_{k\geq 0}$ be the sequence generated in the first-stage of Algorithm \ref{alg:TSA}. Then, for all sufficiently large $n$, it holds with probability at least $1- n^{-\Omega(1)}$ that
\be\label{result:lem-rate-power}
\min_{ s\in \{\pm 1\} } \| \by^{k} - s\bu_1 \|_2 \le  n\cdot \left( \frac{6c_1}{(\alpha-\beta)\sqrt{\log n}} \right)^k
\ee
for all $k\geq 0$, where $\bu_1$ is an eigenvector of $\bB$ associated with the largest eigenvalue and $c_1$ is the constant in Lemma \ref{lem:A-2norm}.
\end{prop}
\begin{proof}
Let $\bB = \bU\bLam\bU^T$ be the eigenvalue decomposition of $\bB$, where $\bLam = \Diag(\lambda_1,\dots,\lambda_n)$ with $\lambda_1\ge\dots\ge\lambda_n$ being the eigenvalues of $\bB$ and $\bU = [\bu_1,\dots,\bu_n]$ are the associated eigenvectors of $\bB$. Suppose that \eqref{eq:lambda1} and \eqref{eq:lambdai} hold, which, according to Lemma \ref{lem:eig-bound}, happens with probability at least $1-2n^{-\frac{1}{2(\alpha+\beta+1)}} - c_2n^{-3}$ for all sufficiently large $n$. Let us define $\bb = \bU^T\by^0$. We claim that if $b_1 >0$, then
\begin{equation}
\label{eq:eig-vec-conv}
 \|\by^k - \bu_1\|_2 \leq \left(\frac{6c_1}{(\alpha-\beta)\sqrt{\log n}}\right)^k\cdot \sqrt{2\sum_{i=2}^n\left(\frac{b_i}{b_1}\right)^2}.
\end{equation}
Indeed, since $\{\by^k\}_{k\geq 0}$ is generated by the power method and $\bB = \bU\bLam\bU^T$, we have
\begin{equation}\label{eq:uy}
\langle\by^k, \bu_1\rangle = \frac{\bu_1^T\bB^k\by^0}{\|\bB^k\by^0\|_2} = \frac{\bu_1^T\bU\bLam^k\bb}{\|\bLam^k\bb\|_2} = \frac{\lambda_1^kb_1}{\sqrt{\sum_{i=1}^n \lambda_i^{2k}b_i^2}}.
\end{equation}
Notice that \eqref{eq:lambda1} and $\alpha>\beta$ imply that $\lambda_1>0$. Hence, if $b_1>0$, then $\lambda_1^kb_1 >0$ for all $k\geq 0$, which yields that $\langle\by^k, \bu_1\rangle >0$ for all $k\geq 0$. This, together with $\|\by^k\|_2 = \|\bu_1\|_2 = 1$, leads to
\begin{align*}
1 & = \|\by^k\|_2 \leq \|\by^k - \langle \by^k,\bu_1\rangle\bu_1\|_2 + \|\langle \by^k,\bu_1\rangle\bu_1\|_2 \\
& = \|\by^k - \langle \by^k,\bu_1\rangle\bu_1\|_2 + \langle \by^k,\bu_1\rangle,
\end{align*}
where the last equality follows from $\|\bu_1\|_2 = 1$ and $\langle\by^k, \bu_1\rangle >0$. In addition, by Cauchy-Schwarz inequality, we have $\langle\by^k,\bu_1\rangle\leq \|\by^k\|_2\|\bu_1\|_2 = 1$. Thus, we obtain
\begin{equation}\label{eq:1-uy}
	0 \leq 1 - \langle \by^k, \bu_1 \rangle = \|\by^k - \langle \by^k,\bu_1\rangle\bu_1\|_2.
\end{equation}
It then follows that
\begin{align*}
& \|\by^k - \bu_1\|_2^2\\
=\ & \|\by^k - \langle\by^k,\bu_1\rangle\bu_1 + \langle\by^k,\bu_1\rangle\bu_1 -\bu_1\|_2^2 \\
=\ & \|\by^k - \langle\by^k,\bu_1\rangle\bu_1\|_2^2 + \|\langle\by^k,\bu_1\rangle\bu_1 -\bu_1\|_2^2 \\
=\ & \|\by^k - \langle\by^k,\bu_1\rangle\bu_1\|_2^2 + (1 - \langle\by^k,\bu_1\rangle)^2 \\
\leq\ & 2\|\by^k - \langle\by^k,\bu_1\rangle\bu_1\|_2^2,
\end{align*}
where the inequality is due to \eqref{eq:1-uy}. In addition, by $\|\by^k\|_2 = \|\bu_1\|_2 = 1$ and \eqref{eq:uy}, we have
\begin{align*}
& \|\by^k - \langle\by^k,\bu_1\rangle\bu_1\|_2^2 = 1 - \langle\by^k,\bu_1\rangle^2 = \frac{\sum_{i=2}^n\lambda_i^{2k}b_i^2}{\sum_{i=1}^n\lambda_i^{2k}b_i^2}\\
& \ \leq \frac{\sum_{i=2}^n\lambda_i^{2k}b_i^2}{\lambda_1^{2k}b_1^2} \leq \left(\frac{\bar{\lambda}}{\lambda_1}\right)^{2k}\sum_{i=2}^n\left(\frac{b_i}{b_1}\right)^2,
\end{align*}
where $\bar{\lambda}:=\max\{|\lambda_i|: i=2,\dots,n\}$. These, together with \eqref{eq:lambda1} and \eqref{eq:lambdai}, establish the claim in \eqref{eq:eig-vec-conv}. By the same arguments, one can show that if $b_1<0$, then
\begin{equation}\label{eq:copy}
 \|\by^k + \bu_1\|_2 \leq \left(\frac{6c_1}{(\alpha-\beta)\sqrt{\log n}}\right)^k\cdot \sqrt{2\sum_{i=2}^n\left(\frac{b_i}{b_1}\right)^2}.
\end{equation}
Since $\by^0$ is chosen according to the uniform distribution over the unit sphere of $\mathbb{R}^n$, which is known to be orthogonally invariant, we have that $\bb = \bU^T\by^0$ also follows the uniform distribution over the unit sphere of $\mathbb{R}^n$. It then follows from Lemma \ref{lem:unit-sphere} that
\begin{equation}\label{eq:bb}
	\sum_{i=2}^n \left(\frac{b_i}{b_1} \right)^2 \leq \frac{n^2}{2}
\end{equation}
holds with probability at least $1 - 2n^{-1/2}$. The desired result then follows from \eqref{eq:eig-vec-conv}, \eqref{eq:copy}, \eqref{eq:bb}, and the union bound.
\end{proof}
One can observe from Proposition \ref{prop:rate-power} that with high probability, the sequence generated in the first-stage of Algorithm \ref{alg:TSA} converges at least linearly to the eigenvector of $\bB$ associated with its largest eigenvalue. Moreover, for any fixed $\alpha$ and $\beta$, \eqref{result:lem-rate-power} shows that the ratio in the linear rate of convergence tends to $0$ as $n\rightarrow \infty$.

\section{Analysis of the Generalized Power Method}\label{sec:gpm}
In this section, we analyze the performance of the GPM that is used in the second-stage of Algorithm \ref{alg:TSA}. In particular, we shall show that with high probability, the sequence of iterates generated by the GPM converges to the ground truth $\bxs$ or $-\bxs$ in a finite number of iterations, provided that its initial iterate is in a suitable neighborhood of $\bxs$ or $-\bxs$. To begin, we present the following two lemmas that will be used for establishing the contraction property of the GPM. Their proofs can be found in the appendix.
\begin{lemma}\label{lem:contrac-GPM}
Suppose that $\bA$ is randomly generated by Model \ref{model:SBM} and $\bB=\bA-\rho\bm{E}_n$ with $\rho=\bo_n^T\bA\bo_n/n^2$. For any $\bx\in\mathbb{R}^n$ with $\|\bx\|_2 = \sqrt{n}$, it holds with probability at least $1 - n^{-\Omega(1)}$ that
\begin{align}
\|\bB\bx - \bB\bx^*\|_2 & \leq \left( \frac{\log n}{\sqrt{n}} + c_1\sqrt{\log n}\right)\|\bx - \bxs\|_2 \nonumber \\
& \quad + \frac{(\alpha-\beta)\log n}{4\sqrt{n}}\|\bx - \bxs\|_2^2. \label{result:lem-contrac-GPM}
\end{align}
\end{lemma}
\begin{lemma}\label{lemma-2}
Let $\delta>0$ be arbitrarily chosen. Then, for any $\bu,\ \bv \in\R^n$ with $|v_i| \ge \delta$ for all $i=1,2,\dots,n$, the following inequality holds:
\be\label{result-2:prop-2}
\left\| \frac{\bu}{|\bu|} - \frac{\bv}{|\bv|} \right \|_2 \le \frac{2 \|\bu-\bv\|_2}{\delta}.
\ee
\end{lemma}
We are now ready to establish the following result, which implies that the GPM used in the second-stage of Algorithm~\ref{alg:TSA} possesses a contraction property.
\begin{lemma}
\label{lem:contrac-prop}
Suppose that $\bA$ is randomly generated by Model \ref{model:SBM} and $\bB=\bA-\rho\bm{E}_n$ with $\rho=\bo_n^T\bA\bo_n/n^2$. For any fixed $\alpha>\beta>0$ such that $\sqrt{\alpha} - \sqrt{\beta} > \sqrt{2}$, there exists a constant $\gamma>0$ such that the following event happens with probability at least $1 - n^{-\Omega(1)}$: for all $\bx\in\mathbb{R}^n$ such that $\|\bx\|_2 = \sqrt{n}$ and
\begin{equation}
\label{eq:neighbor}
\|\bx - \bxs\|_2 \leq \frac{6\sqrt{2}c_1}{\alpha - \beta}\cdot\sqrt{\frac{n}{\log n}},
\end{equation}
it holds that
\begin{equation}
\label{eq:GPM-conv}
\left\|\frac{\bB\bx}{|\bB\bx|} - \bxs\right\|_2 \leq \frac{8c_1}{\gamma\sqrt{\log n}}\cdot\|\bx - \bxs\|_2.
\end{equation}
\end{lemma}
\begin{proof}
Since $\bB = \bA - \rho\bm{E}_n$ and $\bm{E}_n\bxs = \bzz$, we obtain $\bB\bxs = \bA\bxs$. By $\sqrt{\alpha} - \sqrt{\beta} > \sqrt{2}$, there exists a postive constant $\gamma$ such that $(\sqrt{\alpha} - \sqrt{\beta})^2/2 - \gamma\log(\alpha/\beta)/2>1$. These, together with \citep[Lemma 8]{abbe2017community}, yield that
\begin{equation}
\label{eq:opt-tail-bd-B}
\min\left\{ x_i^*(\bB\bxs)_i : i=1,2,\dots,n\right\} \geq \gamma\log n
\end{equation}
holds with probability at least $1 - n^{-\Omega(1)}$. Notice that \eqref{eq:opt-tail-bd-B} implies $\bB\bx^*/|\bB\bx^*| = \bx^*$ and
\begin{equation}
\label{eq:opt-tail-bd-B1}
\min\left\{ \left|(\bB\bxs)_i\right| : i=1,2,\dots,n\right\} \geq \gamma\log n.
\end{equation}
It then follows from \eqref{eq:opt-tail-bd-B1}, Lemma \ref{lem:contrac-GPM}, and Lemma \ref{lemma-2} that for all $\bx\in\mathbb{R}^n$ satisfying $\|\bx\|_2 = \sqrt{n}$ and \eqref{eq:neighbor},
\begin{align*}
& \left\|\frac{\bB\bx}{|\bB\bx|} - \bxs\right\|_2 = \left\|\frac{\bB\bx}{|\bB\bx|} - \frac{\bB\bxs}{|\bB\bxs|}\right\|_2 \leq \frac{2\|\bB\bx - \bB\bxs\|_2}{\gamma\log n} \\

& \leq \left(\frac{(\alpha-\beta)}{2\gamma\sqrt{n}}\|\bx - \bxs\|_2 + \frac{2}{\gamma\sqrt{n}} + \frac{2c_1}{\gamma\sqrt{\log n}}\right)\|\bx - \bxs\|_2 \\
& \leq \left(\frac{3\sqrt{2} c_1}{\gamma\sqrt{\log n}} + \frac{2}{\gamma\sqrt{n}} + \frac{2c_1}{\gamma\sqrt{\log n}}\right)\cdot\|\bx - \bxs\|_2 \\
& \leq \frac{8c_1}{\gamma\sqrt{\log n}}\|\bx - \bx^*\|_2
\end{align*}
for all $n$ sufficiently large.
\end{proof}

From Algorithm \ref{alg:TSA}, we observe that all the iterates $\{\bx^k\}$ generated by the GPM satisfy $\|\bx^k\|_2 = \sqrt{n}$. Indeed, we have from line 7 of Algorithm \ref{alg:TSA} that $\|\bx^0\|_2 = \sqrt{n}$ and $\|\bx^k\|_2 = \sqrt{n}$ for all $k\geq 1$ due to the fact that $\bx^k\in\{\pm 1\}^n$ for all $k\geq 1$.
Thus, Lemma \ref{lem:contrac-prop} implies that for all sufficiently large $n$, the iterates generated by the GPM possess a contraction property once their Euclidean distance to the ground truth is no greater than a threshold that scales linearly with $\sqrt{n/\log n}$. This is summarized in the following proposition.
\begin{prop}\label{prop:rate-GPM}
Suppose that $\bA$ is randomly generated by Model \ref{model:SBM}. Let $\alpha>\beta>0$ be fixed such that $\sqrt{\alpha} - \sqrt{\beta} > \sqrt{2}$. Suppose that the $\bx^0$ in Algorithm \ref{alg:TSA} satisfies $\|\bx^0\|_2 = \sqrt{n}$ and
$$\|\bx^0 - \bx^*\|_2 \leq \frac{6\sqrt{2} c_1}{\alpha - \beta}\sqrt{\frac{n}{\log n}},$$
where $c_1$ is the constant in Lemma \ref{lem:A-2norm}. Then, for all sufficiently large $n$, it holds with probability at least $1-n^{-\Omega(1)}$ that
\begin{align}\label{eq:contract-GPM}
\|\bx^k-\bx^*\|_2 \le \|\bx^0 - \bx^*\|_2\left( \frac{8c_1}{\gamma\sqrt{\log n}} \right)^k,
\end{align}
where $\gamma>0$ is any constant such that $(\sqrt{\alpha} - \sqrt{\beta})^2/2 - \gamma\log(\alpha/\beta)/2>1$.
\end{prop}

Furthermore, the following result indicates that the GPM exhibits finite termination. More precisely, it admits a one-step convergence property once the generated iterate has Euclidean distance no greater than a constant threshold to the ground truth.

\begin{lemma}
\label{lem:local-one-step}
Suppose that $\bA$ is randomly generated by Model \ref{model:SBM} and $\bB=\bA-\rho\bm{E}_n$ with $\rho=\bo_n^T\bA\bo_n/n^2$. For any fixed $\alpha>\beta>0$ such that $\sqrt{\alpha} - \sqrt{\beta} > \sqrt{2}$, the following event happens with probability at least $1 - n^{-\Omega(1)}$: for all $\bx\in\{\pm1\}^n$ such that $\|\bx - \bxs\|_2 \leq 2$, it holds that
\begin{equation}
\label{eq:one-step}
\frac{\bB\bx}{|\bB\bx|} = \bxs.
\end{equation}
\end{lemma}
\begin{proof}
Notice that to prove \eqref{eq:one-step}, it suffices to show that $x_i^*(\bB\bx)_i >0$ for all $i = 1, 2, \dots, n$. Let $\bx\in\{\pm1\}^n$ with $\|\bx - \bxs\|_2 \leq 2$ be arbitrarily chosen. Then, one can observe that $\bx = \bxs$ or $\bx = \bxs \pm 2\bev_l$ for some $l = 1,2,\dots,n$, where $\bev_l$ is the vector of $\mathbb{R}^n$ with all entries being $0$ except the $l$-th entry being $1$. If $\bx = \bxs$, then it follows from \eqref{eq:opt-tail-bd-B} that $x_i^*(\bB\bx)_i >0$ for all $i = 1, 2, \dots, n$. If $\bx = \bxs \pm 2\bev_l$ for some $l = 1,2,\dots,n$, then for any $i=1,2,\dots,n$,
\begin{equation}
\label{eq:case-ii}
\begin{aligned}
x_i^*(\bB\bx)_i & = x_i^*(\bB\bxs \pm 2\bB\bev_l)_i \\
& = x_i^*(\bB\bxs)_i \pm 2x_i^*(\bB\bev_l)_i \\
& \geq \gamma\log n \pm 2x_i^*(\bB\bev_l)_i,
\end{aligned}
\end{equation}
where the inequality follows from \eqref{eq:opt-tail-bd-B}. Notice that
$ \bB\bev_l = \bA\bev_l - \rho\bm{E}_n\bev_l$. Hence, one has $(\bB\bev_l)_i = A_{il} - \rho$. It then follows from Model \ref{model:SBM} and \eqref{result:lem-rho-bound} that
$$ (\bB\bev_l)_i \leq 1 - \frac{(\alpha+ \beta)\log n}{2n} + \frac{\log n}{n^{3/2}}$$
and
$$ (\bB\bev_l)_i \geq - \frac{(\alpha+ \beta)\log n}{2n} - \frac{\log n}{n^{3/2}}. $$
This, together with $|x_i^*| = 1$ and \eqref{eq:case-ii}, implies that for all sufficiently large $n$, $x_i^*(\bB\bx)_i >0$ holds for all $i = 1, 2, \dots, n$.
\end{proof}

\section{Proofs of Main Results}\label{sec:pf-main}

In this section, we provide the proofs of our main results Theorem \ref{thm:main-thm} and Corollary \ref{coro:comp-time}. Recall that the iterates of the PM converge to the eigenvector $\bu_1$ of $\bB$ associated with its largest eigenvalue, not to the ground truth $\bxs$ or $-\bxs$. As such, we need the following lemma that bounds the Euclidean distance between these two vectors. The proof of this lemma can be found in the appendix.
\begin{lemma}\label{lem:dist-eig-truth}
	Suppose that $\bA$ is randomly generated by Model \ref{model:SBM} and $\bB=\bA-\rho\bm{E}_n$ with $\rho=\bo_n^T\bA\bo_n/n^2$. Let $\bu_1$ be the eigenvector of $\bB$ associated with its largest eigenvalue. Then, it holds with probability at least $1 - c_2n^{-3}$ that
	\be\label{dist:eig-truth}
	\min_{ s\in \{\pm 1\} } \left \|\bu_1 - \frac{s\bx^*}{\sqrt{n}} \right \|_2 \le \frac{3\sqrt{2}c_1}{(\alpha-\beta)\sqrt{\log n}},
	\ee
	where $c_1,c_2$ are the constants in Lemma \ref{lem:A-2norm}.
\end{lemma}
Now we are ready to provide the proof of Theorem \ref{thm:main-thm}.
\begin{proof}[Proof of Theorem \ref{thm:main-thm}]
	Suppose that \eqref{result:lem-rate-power} and \eqref{eq:contract-GPM} hold simultaneously, which, according to Propositions \ref{prop:rate-power} and \ref{prop:rate-GPM}, happens with probability at least $1-n^{-\Omega(1)}$. In the first stage, we require the power iterations to output a $\by\in\mathbb{R}^n$ such that $\|\by\|_2 = 1$ and $\|\by-\bu_1\|_2\le 3\sqrt{2}c_1/({(\alpha-\beta)\sqrt{\log n}})$. In view of Proposition \ref{prop:rate-power}, this can be achieved in $N_p$ power iterations, where $N_p$ is the smallest integer satisfying
	\begin{align}\label{eq-1:main-thm}
	n \left( \frac{6c_1}{(\alpha-\beta)\sqrt{\log n}} \right)^{N_p} \le \frac{3\sqrt{2}c_1}{(\alpha-\beta)\sqrt{\log n}}.
	\end{align}
	Thus, we have
	\[
	N_p = \left\lceil\frac{\log\left((\alpha-\beta)n\sqrt{\log n}\right) - \log(3\sqrt{2}c_1)}{\log\left((\alpha-\beta)\sqrt{\log n}\right) - \log(6c_1)} \right\rceil,
	\]
	which is roughly $O(\log n / \log\log n)$.
	Then, the second-stage of Algorithm \ref{alg:TSA} starts with $\bx^0 = \sqrt{n}\by$, which satisfies
	\begin{align*}
	\|\bx^0 - \bx^*\|_2 & \le \|\bx_0-\sqrt{n}\bu_1\|_2 + \|\sqrt{n}\bu_1 - \bx^*\|_2 \\
	& \le  \frac{6\sqrt{2}c_1\sqrt{n}}{(\alpha-\beta)\sqrt{\log n}},
	\end{align*}
	where the second inequality follows from \eqref{dist:eig-truth} and \eqref{eq-1:main-thm}. Thus, in the second stage, according to Proposition \ref{prop:rate-GPM} and Lemma \ref{lem:local-one-step}, we have that the GPM terminates at $\bx$ or $-\bxs$ in at most $N_g$ iterations, where $N_g$ is the smallest integer satisfying
	\begin{align*}
	\frac{6\sqrt{2}c_1\sqrt{n}}{(\alpha-\beta)\sqrt{\log n}}\left( \frac{8c_1}{\gamma\sqrt{\log n}} \right)^{N_g} \le 2.
	\end{align*}
	Thus, we have
	\[
	N_g = \left\lceil\frac{\log(3\sqrt{2}c_1\sqrt{n}) - \log((\alpha-\beta)\sqrt{\log n})}{\log(\gamma\sqrt{\log n}) - \log(8c_1)}\right\rceil,
	\]
	which is also roughly $O(\log n/\log\log n)$.
\end{proof}
As for the per-iteration cost of Algorithm \ref{alg:TSA}, it mainly depends on the cost of computing the matrix-vector product $\bm{B}v$ for some $\bv\in \R^n$ with $\bm{B}=\bm{A}-\rho \bo_n\bo_n^T$. Note that the cost of computing $\bB\bv$ is determined by the sparsity of $\bm{A}$. As such, we need the following lemma that bounds the number of non-zero entries in $\bA$. The proof of this lemma can be found in the appendix.
\begin{lemma}\label{lem:A-sparsity}
Let $\bA$ be randomly generated by Model \ref{model:SBM} with $\alpha,\ \beta>0$. It holds with probability at least $1 - n^{-\Omega(1)}$ that the number of non-zero entries in $\bA$ is less than $2(\alpha+\beta)n\log n$.
\end{lemma}
Armed with the above results, we can prove Corollary \ref{coro:comp-time}.
\begin{proof}[Proof of Corollary \ref{coro:comp-time}]
According to Algorithm \ref{alg:TSA} and Lemma \ref{lem:A-sparsity}, the time complexity of each iteration is $2(\alpha+\beta)n\log n$. By Theorem \ref{thm:main-thm}, we have the total iteration number of Algorithm \ref{alg:TSA} is $O(\log n)$. Thus, the total time complexity of Algorithm \ref{alg:TSA} is $O(n\log^2 n)$ with probability at least $1-n^{-\Omega(1)}$.
\end{proof}

\begin{figure*}[t]
	\begin{minipage}[b]{0.245\linewidth}
		\centering
		\centerline{\includegraphics[width=\linewidth]{fig-GPM-1.eps}}

		\centerline{(a)}\medskip
	\end{minipage}

	\begin{minipage}[b]{0.245\linewidth}
		\centering
		\centerline{\includegraphics[width=\linewidth]{fig-SDP-1.eps}}

		\centerline{(b)}\medskip
	\end{minipage}
	\begin{minipage}[b]{0.245\linewidth}
		\centering
		\centerline{\includegraphics[width=\linewidth]{fig-MGD-1.eps}}

		\centerline{(c)}\medskip
	\end{minipage}
	\begin{minipage}[b]{0.245\linewidth}
		\centering
		\centerline{\includegraphics[width=\linewidth]{fig-SC-1.eps}}

		\centerline{(d)}\medskip
	\end{minipage}
	\vskip -0.2in
	\caption{Phase transition: the $x$-axis is $\beta$, the $y$-axis is $\alpha$, and darker pixels represent lower empirical probability of success. The red curve is the information-theoretic threshold $\sqrt{\alpha}-\sqrt{\beta}=\sqrt{2}$.}
	\label{fig-1}
\end{figure*}

\begin{figure*}[t]
\begin{center}
	\begin{minipage}[b]{0.28\linewidth}
		\centering
		\centerline{\includegraphics[width=\linewidth]{f-1000.eps}}

		\centerline{(a)}\medskip
	\end{minipage}

	\begin{minipage}[b]{0.28\linewidth}
		\centering
		\centerline{\includegraphics[width=\linewidth]{f-5000.eps}}

		\centerline{(b)}\medskip
	\end{minipage}
	\begin{minipage}[b]{0.28\linewidth}
		\centering
		\centerline{\includegraphics[width=\linewidth]{f-10000.eps}}

		\centerline{(c)}\medskip
	\end{minipage}
\end{center}
	\vskip -0.2in
	\caption{Convergence performance: the $x$-axis is number of iterations, the $y$-axis for GPM is $\|\bx^k{\bx^k}^T - \bxs{\bxs}^T\|_F$, and the $y$-axis for MGD is $\|\bm{Q}^k{\bm{Q}^k}^T - \bxs{\bxs}^T\|_F$, where $\bx^k$ and $\bm{Q}^k$ are the iterates generated in the $k$-th iteration of GPM and MGD, respectively.}
	\label{fig-2}
\end{figure*}

\section{Experimental Results}\label{sec:num}

In this section, we conduct numerical experiments on synthetic data sets to test the performance of the proposed two-stage approach for community recovery in the binary symmetric SBM and compare it with several existing methods. For ease of reference, we denote our two-stage approach simply by GPM in this section. Our codes are implemented in MATLAB R2019a. The experiments are conducted on a PC with 16GB memory and Intel(R) Core(TM) i5-8600 3.10GHz CPU.

\subsection{Phase Transition and Computational Efficiency}

We first conduct the experiment of phase transition to test the recovery performance of our approach GPM and compare it with the SDP-based approach in \citet{amini2018semidefinite}, the manifold optimization (MFO) based approach in \citet{bandeira2016low}, and the spectral clustering (SC) approach in \citet{abbe2017entrywise}. In the implementation, we use alternating direction method of multipliers (ADMM) to solve the SDP as suggested in \citet{amini2018semidefinite}, manifold gradient descent (MGD) method to solve the MFO, and the Matlab function \textsf{eigs} for computing the eigenvector that is needed in the SC approach.  We choose $n = 300$ in the experiment and let $\alpha$ and $\beta$ in Model \ref{model:SBM} vary from $0$ to $30$ and $0$ to $10$, with increments $0.5$ and $0.4$, respectively. For every pair of $\alpha$ and $\beta$, we generate $40$ instances and calculate, for all the tested methods, the ratio of exact recovery. The simulation results are presented in Figure \ref{fig-1}. We can observe that all the methods admit a phase transition that is close to the information-theoretic limit, with the recovery performance of the GPM and the SC slightly better than the other two. This supports our theoretical development that the proposed algorithm can achieve exact recovery at the information-theoretic limit. Besides, as an indicator of computational efficiency, we record in Table \ref{table-1} the total CPU time consumed by every approach for completing the phase transition experiments. It can be observed that the GPM slightly outperforms the SC in terms of CPU time, and both of them are substantially computationally more efficient than the SDP and the MGD.
\vskip -0.4cm
\begin{table}[!htbp]
\caption{Total CPU times (in seconds) consumed by the test approaches in the experiment of phase transition.}
\label{table-1}
\vskip 0.15in
\begin{center}
\begin{small}
\begin{tabular}{lcccr}
\toprule
Methods & GPM & SDP &  MGD & SC\\
\midrule
Time (s) & 28 & 9193 &  1063 & 124 \\
\bottomrule
\end{tabular}
\end{small}
\end{center}
\vskip -0.2in
\end{table}

\subsection{Convergence Performance}

Next we conduct the experiments of convergence performance to test the number of iterations needed by our approach GPM to exactly identify the underlying communities. For comparison, we also test the convergence performance of MGD, which is an iterative algorithm that has similar per-iteration cost to GPM. In the experiments, we use $\alpha=10$ and $\beta=2$ and generate graphs of dimensions $n=1000$, $5000$, and $10000$ by Model \ref{model:SBM}. We use $\|\bx^k{\bx^k}^T-\bx^*{\bx^*}^T\|_F$ (resp. $\|\bm{Q}^k{\bm{Q}^k}^T-\bx^*{\bx^*}^T\|_F)$ as the measure of distance from $\bx^k$ (resp. $\bm{Q}^k$) to the ground truth $\bxs$ or $-\bxs$, where $\bx^k$ (resp. $\bm{Q}^k$) is the iterate generated in the $k$-th iteration of GPM (resp. MGD). In Figure \ref{fig-2}, we plot this measure of distance to the ground truth against the iteration number for both GPM and MGD. It is observed that GPM exhibits a finite termination phenomenon and converges to the ground truth much faster than MGD. This also corroborates our theoretical development that our approach identifies the underlying communities in at most $O(\log n/\log\log n)$ number of PM and GPM iterations.

\section{Concluding remarks}

In this work, we consider the problem of exact community recovery in the binary symmetric SBM. We propose a two-stage iterative algorithm to solve it, which employs the power method firstly and then turns to the generalized power method. We show that the proposed method can achieve exact recovery at the information-theoretic limit within $\tilde{O}(n)$ time complexity as $n$ goes to infinity. Numerical experiments on synthetic data sets demonstrate that the proposed approach has strong recovery performance and is highly efficient.

As suggested in our numerical simulations, applying Algorithm \ref{alg:TSA} to Problem \eqref{eq:rMLE} seems to have a significantly larger one-step convergence region than the constant region in Lemma \ref{lem:local-one-step}, which can be further improved in future research. Also, it is worthy to study whether the proposed
approach can be extended to more general SBMs.

\section*{Acknowledgements}

The first and third authors are supported in part by the Hong Kong Research Grants Council (RGC) General Research Fund (GRF) project CUHK 14208117 and in part by the CUHK Research Sustainability of Major RGC Funding Schemes project 3133236. The second author is supported in part by the National Natural Science Foundation of China (NSFC) project 11901490 and in part by an HKBU Start-up Grant. Most of the work of the second author was done when he was affiliated with the Department of Mathematics of Hong Kong Baptist University.

\bibliographystyle{icml2020}


TITLE: Numerical solution of the Laplace equation in semi-circle and quarter-circle
QUESTION [1 upvotes]: Problem: Find the numerical solution of $\Delta u = 0$ in semi-circle $0 \leq r \leq R$, $0 \leq \phi \leq \pi$ with boundary conditions $$u_{\big|r=R} = T$$
$$u_{\big|\phi=0}=u_{\big|\phi=\pi}=0$$
where $T>0.$
After that, notice that the solution is symmetric with respect to the vertical axis and try to solve the problem only on the quarter-circle $0 \leq \phi \leq \frac{\pi}{2}$.
I managed to do the first part. I used the grid given with
$$(r_i,\phi_j) = (i \cdot \Delta r, j \cdot \Delta \phi),$$
where $\Delta r = \frac{R}{M+1}, \Delta \phi = \frac{\pi}{N+1}$ for $M,N \in \mathbb{N}.$ Also, I noticed the relation between the solution in $5$ nodes and then I've found the equation:
$$\big(1-\frac{\Delta r}{2r_i}\big)u_{i-1,j} - 2 \big[ 1 + \big( \frac{\Delta r}{r_i \Delta \phi}\big)^2 \big] u_{i,j} + \big( 1 + \frac{\Delta r}{2r_i} \big) u_{i+1,j} + \big( \frac{\Delta r}{r_i \Delta \phi} \big) ^2 u_{i,j-1} + \big( \frac{\Delta r}{r_i \Delta \phi} \big) ^2 u_{i,j+1} = 0$$
and implemented everything in Matlab successfully.
I have a problem with the part concerning the quarter-circle. I understand that $\phi \mapsto \pi - \phi$ doesn't change the solution but I'm not sure why.
Most importantly, I can't reach the boundary condition for $\phi = \frac{\pi}{2}.$ I know that for fixed $r$ the solution has the biggest value here and it gets bigger for bigger $r$, but I'm not able to find the equation analogue to the one I computed.
Any hint would be appreciated.

REPLY [1 votes]: I have a problem with the part concerning the quarter-circle. I understand that $\phi\mapsto\pi-\phi$ doesn't change the solution but I'm not sure why.

The whole problem has a $\phi$-symmetry wrt $\phi=\frac{\pi}{2}$:
the domain and boundary condition is obvious but the same holds true for the Laplacian.
This follows by comparing the Laplacian wrt $(r,\phi)$ and $(r,\phi':=\pi-\phi)$
\begin{align}
  \Delta_{(r,\phi)} &= r^{-1} \partial_r (r \partial_r ) + r^{-2} \partial_\phi^2
  \\
  \Delta_{(r,\phi')} &= r^{-1} \partial_r (r \partial_r) + r^{-2} \partial_{\phi'}^2
  \\
  &= r^{-1} \partial_r (r \partial_r) + r^{-2} (-\partial_\phi)^2
  \\
  &= r^{-1} \partial_r (r \partial_r) + r^{-2} \partial_\phi^2
  \\
  &= \Delta_{(r,\phi)}
\end{align}
where we used $\partial_{\phi'} = - \partial_\phi \,.$
Thus, there is no possibility for asymmetry in $\phi$.

Most importantly, I can't reach the boundary condition for $\phi=\frac{\pi}{2}$.

I am not quite sure if you understood the idea correct as there is no boundary at $\phi=\frac{\pi}{2}$.
You still solve for the whole (!) semi-circle but you do that by only using variables (grid nodes) for $ 0 \leq \phi \leq \frac{\pi}{2} $.
Your linear equation system will contain equations for each grid node of the quarter circle and whenever there appears a grid in your discretized equation from the other quarter circle just use the reflected one.
You will also need to assign Dirichlet boundary conditions at $(\phi=0, 0\leq r \leq R)$ and at $(r=R, 0\leq \phi \leq \frac{\pi}{2})$.


TITLE: Simplify $\sum \limits_{k=0}^{n} \binom{n}{k} 2^{\sqrt{k}}$
QUESTION [10 upvotes]: Can this sum be simplified: $\sum \limits_{k=0}^{n} \binom{n}{k} 2^{\sqrt{k}}$
Or at least is there a simple fairly tight upperbound?
EDIT
So I think this sum is more easily bounded than I previously thought:
Clearly, $$\sum_{k=0}^{n} \binom{n}{k} 2^{\sqrt{k}} \cdot \leq 2^{\sqrt{n}} \sum_{k=0}^{n} \binom{n}{k} = 2^{n+\sqrt{n}} .$$
Also, along the same lines as Shai Covo's answer, $\sum \limits_{k=0}^{n} \binom{n}{k} 2^{\sqrt{k}} \geq \binom{n}{n/2} 2^{\sqrt{n/2}}$. The central binomial coefficient: $\binom{n}{n/2}$ is at least $\frac{2^n}{\sqrt{2n}}$, hence $$\sum_{k=0}^{n} \binom{n}{k} 2^{\sqrt{k}} \geq 2^{\sqrt{n/2}} \cdot \frac{2^n}{\sqrt{2n}} = \frac{2^{n + \sqrt{n/2}}}{2n}$$
So $\sum \limits_{k=0}^{n} \binom{n}{k} 2^{\sqrt{k}}$ is $O(2^{n + \sqrt{n}})$ and $\Omega(\frac{2^{n+\sqrt{n/2}}}{n})$
What about expressions of the form:
$\sum \limits_{k=0}^{n} \binom{n}{k} a^{\sqrt{k}} b^{\sqrt{n-k}}$?

REPLY [2 votes]: The mass in the sum is approximately proportional to a (shifted) Gaussian centered at $k = n/2 + C\sqrt{n}$ and with standard deviation $A \sqrt{n}$ for explicitly calculable constants $C$ and $A$.  This implies asymptotics of the form $(M + o(1))2^{n + \sqrt{n/2}}$ where $M$ is another computable constant.
The same is true for the $a,b$ version of the problem.
[edit: my calculations give $M = 2^{(\ln 2)/8} = 1.0618966...$.
[edit-2: similar calculations give $M=(a/b)^{\ln (a/b)}/8$ and asymptotics $(M+o(1))2^n (ab)^\sqrt{n/2}$ for the $a,b$ sum.]