TITLE: Calculating a limit with constraints
QUESTION [1 upvotes]: Given the function $f(x)$,
$$ f(x,y,z,w) = \frac{x+iy}{\sqrt{|w+z|}} \text{.} $$
How do I calculate the limit
$$ \lim\limits_{w\rightarrow -z} f $$
under the constraint that the points $(x,y,z) \in \mathbb{R}^3$ lie on a sphere with radius $w$:
$$ x^2+y^2+z^2=w^2 \text{.} $$
Note that the constraint in the limit $w \rightarrow - z$ implies $x \rightarrow 0$ and $y \rightarrow 0$. Numerically, I obtained the following conjecture:
$$ \sqrt{z} (1+i) $$
How do I prove this analytically?

REPLY [2 votes]: I post this as an answer because it would be a bit long for a comment. But I am still not sure to really understand the question.
I don't think you can expect a limit in general.
If we choose $y=0$ and $x=\sqrt{w^2-z^2}$ then the conditions you set are satisfied and (if it makes sense)
$$\frac{x+iy}{\sqrt{|w+z|}}=\frac{\sqrt{w^2-z^2}\sqrt{|w-z|}}{\sqrt{|w+z|}\sqrt{|w-z|}}=\sqrt{|w-z|}\to\sqrt{2|z|}\,.$$
If we choose $x=0$ and $y=\sqrt{w^2-z^2}$ then the conditions you set are satisfied and (if it makes sense)
$$\frac{x+iy}{\sqrt{|w+z|}}=\frac{i\sqrt{w^2-z^2}\sqrt{|w-z|}}{\sqrt{|w+z|}\sqrt{|w-z|}}=i\sqrt{|w-z|}\to i\sqrt{2|z|}\,.$$


text:
\begin{document}
\title{Uniqueness of degree-one Ginzburg-Landau vortex in the unit ball in dimensions $N\geq 7$}

\author{Radu Ignat\thanks{Institut de Math\'ematiques de Toulouse \& Institut Universitaire de France, UMR 5219, Universit\'e de Toulouse, CNRS, UPS
IMT, F-31062 Toulouse Cedex 9, France. Email: Radu.Ignat@math.univ-toulouse.fr
}~, Luc Nguyen\thanks{Mathematical Institute and St Edmund Hall, University of Oxford, Andrew Wiles Building, Radcliffe Observatory Quarter, Woodstock Road, Oxford OX2 6GG, United Kingdom. Email: luc.nguyen@maths.ox.ac.uk}~, Valeriy Slastikov\thanks{School of Mathematics, University of Bristol, University Walk, Bristol, BS8 1TW, United Kingdom. Email: Valeriy.Slastikov@bristol.ac.uk}~ and Arghir Zarnescu\thanks{IKERBASQUE, Basque Foundation for Science, Maria Diaz de Haro 3,
48013, Bilbao, Bizkaia, Spain.}\,  \thanks{BCAM,  Basque  Center  for  Applied  Mathematics,  Mazarredo  14,  E48009  Bilbao,  Bizkaia,  Spain.
(azarnescu@bcamath.org)}\, \thanks{``Simion Stoilow" Institute of the Romanian Academy, 21 Calea Grivi\c{t}ei, 010702 Bucharest, Romania.}}

\date{}

\maketitle
\begin{abstract}
For $\eps>0$, we consider the Ginzburg-Landau functional for $\R^N$-valued maps defined in the unit ball $B^N\subset \R^N$ with the vortex boundary data $x$ on $\partial B^N$.
In dimensions $N\geq 7$, we prove that for every $\eps>0$, there exists a unique global minimizer $u_\eps$ of this problem; moreover, $u_\eps$ is symmetric {and of the form}
$u_\eps(x)=f_\eps(|x|)\frac{x}{|x|}$ for $x\in B^N$.  \\

\noindent {\it Keywords: uniqueness, symmetry, minimi{z}ers, Ginzburg-Landau.}

\noindent {\it MSC: 35A02, 35B06, 35J50.}
\end{abstract}

\section{Introduction and main results}

In this note, we consider the following Ginzburg-Landau type energy functional
$$
E_\eps(u)=\int_{B^N} \Big[\frac{1}{2}|\nabla u|^2+\frac{1}{2\eps^2}W(1 - |u|^2)\Big]\,dx,
$$
where $\eps>0$, $B^N$ is the unit ball in $\R^N$, {$N \geq 2$}, and the potential $W\in C^1((-\infty,1] {\blue ;} \R)$ satisfies

\begin{equation}
 W(0)=0,\, W(t)>0 \hbox{ for all } t\in (-\infty, 1]\setminus \{0\}, {\text{ and }} W \textrm{ is convex}.
 	\label{Eq:26VI18-E1}
\end{equation}
We investigate the global minimi{z}ers of the energy $E_\eps$ { in the set }
$$
\mcA:=\{ u\in H^1(B^N; \R^N):\,  u(x)=x \textrm{ on } \partial B^N = \bS^{N-1}\}.
$$
{The requirement that $u(x) = x$ on $\bS^{N-1}$ is sometimes referred {in the literature} as the vortex boundary condition.}

We note that in our analysis the convexity of $W$ needs not be strict; compare \cite{INSZ_AnnENS} where strict convexity is assumed.

The direct method in the calculus of variations yields the existence of a global minimizer $u_\eps$ of $E_\eps$ over $\mcA$ for all range of $\eps>0$. Moreover, any minimizer $u_\eps$ belongs to $C^1(\overline{B^N}; \R^N)$ and satisfies $|u_\eps| \leq 1$ and the system of PDEs (in the sense of distributions)
\be
\label{E-L}
-\Delta u_\eps=\frac1{\eps^2} u_\eps \, W'(1-|u_\eps|^2)\quad \textrm{ in } \, B^N.
\ee

The goal of this note is to give a short proof of the uniqueness and symmetry of the global minimizer of $E_\eps$ in $\mcA$ for \emph{all} $\eps>0$ in dimensions $N \geq 7$. We prove that,  in these dimensions, the global minimizer is unique and given by the unique radially symmetric critical point of $E_\eps$ defined by

\be
\label{def_sol_equa}
u_\eps(x)=f_\eps(|x|)\frac{x}{|x|} \quad \textrm{ for all } x\in B^N,
\ee

where the radial profile $f_\eps : [0,1]\to \R_+$ is the unique solution of

\be
\label{1}
\left\{\begin{array}{l}
-f''_\eps-\frac{N-1}{r}f'_\eps+\frac{N-1}{r^2}f_\eps=\frac{1}{\eps^2} f_\eps \,W'(1-f_\eps^2) \quad
\text{ for } r \in (0,1),\\
f_\eps(0)=0, f_\eps(1)=1.
\end{array}
\right.
\ee

Moreover, $f_\eps>0$ and $f'_\eps>0$ in $(0,1)$ (see e.g. \cite{ODE_INSZ}).

\medskip

\begin{theorem}
\label{GL_equator}
Assume that $W$ satisfies \eqref{Eq:26VI18-E1}. If $N\geq 7$, then for every $\eps>0$, $u_\eps$ given in \eqref{def_sol_equa} is the unique global minimizer of $E_\eps$  in $\mcA$.
\end{theorem}

To our knowledge, the question about the uniqueness of minimizers/critical points of $E_\eps$ in $\mcA$ for any $\eps > 0$ was raised in dimension $N = 2$ in the book of Bethuel, Brezis and H\'elein \cite[Problem 10, page 139]{BBH_book}, and in general dimensions $N \geq 2$ and also for the  blow-up limiting problem around the vortex (when the domain is the whole space $\R^N$ and by rescaling, $\eps$ can be assumed equal to $1$) in an article of Brezis \cite[Section 2]{Brezis99_PSPM}.

It is well known that uniqueness is present for large enough $\eps>0$ for any $N \geq 2$. Indeed,  for any $\eps > (W'(1)/\lambda_1)^{1/2}$ where $\lambda_1$ is the first eigenvalue of $-\Delta$ in $B^N$ with  zero Dirichlet boundary condition, $E_\eps$ is strictly convex in $\mcA$ and thus has a unique critical point in $\mcA$ {(that is the global minimizer of our problem)}.

For \emph{sufficiently small} $\eps> 0$ all results regarding uniqueness question available in the literature are in the affirmative.  In particular, we have:
\begin{enumerate}[(i)]
\item Pacard and Rivi\`ere \cite[Theorem 10.2]{Pacard_Riviere} showed in dimension $N = 2$ that, for small $\eps > 0$, $E_\eps$ has in fact a unique critical point in $\mcA$.

\item Mironescu \cite{Mironescu_symmetry} showed in dimension $N = 2$ that, when $B^2$ is replaced by $\RR^2$ and $\eps = 1$, a local minimizer of $E_\eps$ subjected to a degree-one boundary condition at infinity is unique (up to translation and suitable rotation). This was generalized to dimension $N = 3$ by Millot and Pisante \cite{mil-pis} and dimensions $N \geq 4$ by Pisante \cite{Pisante11-JFA}, also in the case of the blow-up limiting problem on $\RR^N$ and $\eps = 1$.
\end{enumerate}
These results should be compared to those for the limit problem on the unit ball obtained by sending $\eps \rightarrow 0$. In this limit, the Ginzburg-Landau problem `converges' to the harmonic map problem from $B^N$ to $\bS^{N-1}$. It is well known that, the vortex boundary condition gives rise to a unique minimizing harmonic map $x \mapsto \frac{x}{|x|}$ if $N\geq 3$; see Brezis, Coron and
Lieb~\cite{BrezisCoronLieb}  in dimension $N = 3$, J\"{a}ger and Kaul \cite{JagKaul} in dimensions $N \geq 7$, and Lin \cite{Lin-CR87} in dimensions $N \geq 3$.

We highlight that, in contrast to the above, our result holds for \emph{all} $\eps > 0$, provided that $N \geq 7$. The method of our proof deviates somewhat from that in the aforementioned works. In fact it is reminiscent of our recent work \cite{INSZ_AnnENS} on the (non-)uniqueness and symmetry of minimizers of the Ginzburg-Landau functionals for $\RR^M$-valued maps defined on $N$-dimensional domains, where $M$ is not necessarily the same as $N$. However we note that the results in \cite{INSZ_AnnENS} do not  directly apply to the present context, as in \cite{INSZ_AnnENS} it is required that $W$ be \emph{strictly convex}. Furthermore, a priori, it is not clear why non-strict convexity of the potential $W$ is sufficient to ensure uniqueness of global minimizers.

We exploit the convexity of $W$ to lower estimate the `excess' energy by a suitable quadratic energy which can be handled by the factorization trick \`a la Hardy. Indeed, the positivity of the excess energy is then related to the validity of a Hardy-type inequality, which explains our restriction of $N \geq 7$. This echoes our observation made in \cite{INSZ_AnnENS} that a result of J\"ager and Kaul \cite{JagKaul} on the minimality of the equator map in these dimensions is related to a certain inequality involving the sharp constant in the Hardy inequality.

We expect that our result remains valid in dimensions $2 \leq N \leq 6$, but this goes beyond the scope of this note and remains for further investigation.

\section{Proof of Theorem \ref{GL_equator}}

Theorem \ref{GL_equator} will be obtained as a consequence of a stronger result on the uniqueness of global minimizers of for the $\RR^M$-valued Ginzburg-Landau functional with $M \geq N$. By a slight abuse of notation, we consider the energy functional
$$
E_\eps(u)=\int_{B^N} \Big[\frac{1}{2}|\nabla u|^2+\frac{1}{2\eps^2}W(1 - |u|^2)\Big]\,dx,
$$
where $u$ belongs to
$$
\mcA:=\{ u\in H^1(B^N; \R^M):\,  u(x)=x \textrm{ on } \partial B^N = \bS^{N-1} \subset \RR^M\}.
$$

\begin{theorem}
\label{GL_equatorNM}
Assume that $W$ satisfies \eqref{Eq:26VI18-E1}. If $M \geq N\geq 7$, then for every $\eps>0$, $u_\eps$ given in \eqref{def_sol_equa} is the unique global minimizer of $E_\eps$  in $\mcA$.
\end{theorem}

When $W$ is strictly convex, the above theorem is proved in \cite{INSZ_AnnENS}; see Theorem 1.7. The argument therein uses the strict convexity in a crucial way.

\begin{proof}

The proof will be done in several steps. First, we consider the difference between the energies of the critical point $u_\eps$,  defined in \eqref{def_sol_equa},  and an arbitrary competitor $u_\eps+v$ and show that this difference is controlled from below by some quadratic energy functional $F_\eps(v)$. Second, we employ the positivity of the radial profile $f_\eps$ in \eqref{1} and apply the Hardy decomposition method in order to show that $F_\eps(v)\geq 0$, which proves in particular that $u_\eps$ is a global minimizer of $E_\eps$. Finally, we characterise the situation when this difference is zero and conclude to the uniqueness of the global minimizer $u_\eps$.

\bigskip
\par\noindent{\it Step 1: Lower bound for energy difference.}
For any  $v\in H^1_0(B^N; \R^M)$, we have
\begin{align*}
E_\eps(u_\eps+v)-E_\eps(u_\eps)
	&=\int_{B^N}\Big[ \nabla u_\eps\cdot\nabla v +\frac{1}{2}|\nabla v|^2\Big]\,dx\\
		&\qquad +\frac{1}{2\eps^2}\int_{B^N} \Big[W(1-|u_\eps+v|^2)-W(1-|u_\eps|^2)\Big]\,dx.
\end{align*}
Using the convexity of $W$, we have
\[
W(1-|u_\eps+v|^2)-W(1-|u_\eps|^2)\ge -W'(1-|u_\eps|^2)(|u_\eps+v|^2-|u_\eps|^2).
\]
The last two relations imply that
\begin{align*}
E_\eps(u_\eps+v)-E_\eps(u_\eps)
	&\ge \int_{B^N} \Big[\nabla u_\eps\cdot\nabla v-\frac{1}{\eps^2} W'(1-f_\eps^2)u_\eps\cdot v\Big]\,dx\nonumber\\
		&\qquad\qquad +\int_{B^N} \Big[ \frac{1}{2} |\nabla v|^2-\frac{1}{2\eps^2}W'(1-f_\eps^2)|v|^2\Big]dx.
\end{align*}
Moreover, by \eqref{E-L}, we obtain
\bea\label{rel:endif3}
E_\eps(u_\eps+v)-E_\eps(u_\eps) \ge  \int_{B^N} \Big[ \frac 12 |\nabla v|^2-\frac{1}{2\eps^2} W'(1-f_\eps^2) |v|^2\Big]\,dx =: \frac12 F_\eps(v)
\eea
for all  $v\in H^1_0({B^N};\R^M)$.

\bigskip
\par\noindent{\it Step 2: A rewriting of $F_\eps(v)$ using the decomposition $v=f_\eps w$ for every scalar test function $v\in C^\infty_c(B^N\setminus \{0\}; \R)$.}
We consider the operator
$$L_\eps:=\frac12 \nabla_{L^2}F_\eps=-\Delta-\frac1{\eps^2} W'(1-f_\eps^2).$$
Using the decomposition $$v=f_\eps w$$ for the scalar function $v\in C^\infty_c(B^N\setminus \{0\}; \R)$, we have (see e.g. \cite[Lemma~A.1]{INSZ3}):
\begin{align*}
F_\eps(v)=\int_{B^N} L_\eps v\cdot v\, dx&=\int_{B^N} w^2 L_\eps f_\eps\cdot f_\eps\, dx+\int_{B^N} f_\eps^2 |\nabla w|^2\, dx\\
& =\int_{B^N} f_\eps^2 \bigg( |\nabla w|^2-\frac{N-1}{r^2}w^2\bigg)\, dx,
\end{align*}
because \eqref{1} yields $L_\eps f_\eps\cdot f_\eps=-\frac{N-1}{r^2} f_\eps^2$ in ${B^N}$.

\bigskip
\par\noindent{\it Step 3: We prove that  $F_\eps(v)\geq 0$ for every scalar test function $v\in C^\infty_c(B^N\setminus \{0\}; \R)$.}
Within the notation $v=f_\eps w$ of Step 2 with $v,w \in  C^\infty_c(B^N\setminus \{0\}; \R)$, we use the decomposition $$w=\f g$$ with $\f=|x|^{-\frac{N-2}{2}}$ being the first eigenfunction of the Hardy's operator $-\Delta-\frac{(N-2)^2}{4|x|^2}$ in $\R^N\setminus\{0\}$ and $g\in C^\infty_c({B^N}\setminus \{0\}; \R)$. We compute
\begin{align*}
|\nabla w|^2=|\nabla \f|^2 g^2+|\nabla g|^2 \f^2+\frac12 \nabla (\f^2)\cdot \nabla (g^2).
\end{align*}
As $|\nabla \f|^2=\frac{(N-2)^2}{4|x|^2}\f^2$ and $\f^2$ is harmonic in ${B^N}\setminus\{0\}$, integration by parts yields
\begin{align}
\nonumber
F_\eps(v)&=\int_{B^N} f_\eps^2 \bigg( |\nabla g|^2 \f^2+\frac{(N-2)^2}{4r^2}\f^2g^2-\frac{N-1}{r^2}\f^2 g^2\bigg)\, dx-\frac12 \int_{{B^N}} \nabla (\f^2)\cdot \nabla (f_\eps^2)g^2\, dx\\
\nonumber &\geq \int_{B^N} f_\eps^2 |\nabla g|^2 \f^2\, dx+\bigg(\frac{(N-2)^2}{4}-(N-1)\bigg)\int_{{B^N}}\frac{f_\eps^2}{r^2}\f^2g^2\, dx\\
\label{23}
&\geq \bigg(\frac{(N-2)^2}{4}-(N-1)\bigg)\int_{{B^N}}\frac{v^2}{r^2}\, dx\geq 0,
\end{align}
where we have used $N\geq 7$ and $\frac12\nabla (\f^2)\cdot \nabla (f_\eps^2)=2\f \f' f_\eps f'_\eps\leq 0$ in ${B^N}\setminus \{0\}$.

\bigskip
\par\noindent{\it Step 4: We prove that $F_\eps(v)\geq 0$ for every $v\in H^1_0(B^N; \R^M)$ meaning that $u_\eps$ is a global minimizer of
$E_\eps$ over $\mcA$; moreover, $F_\eps(v)=0$ if and only if $v=0$.} Let $v\in H^1_0(B^N; \R^M)$. As a point has zero $H^1$ capacity in $\R^N$, a standard density argument implies the existence of a sequence $v_k\in C^\infty_c({B^N}\setminus \{0\}; \R^M)$ such that $v_k\to v$ in $H^1(B^N, \R^M)$ and a.e. in $B^N$. On the one hand, by
definition \eqref{rel:endif3} of $F_\eps$, since $W'(1-f_\eps^2)\in L^\infty$, we deduce that $F_\eps(v_k)\to F_\eps(v)$ as $k\to \infty$.
On the other hand, by \eqref{23} and Fatou's lemma, we deduce
\begin{align*}
\liminf_{k\to \infty} F_\eps(v_k)&\geq  \bigg(\frac{(N-2)^2}{4}-(N-1)\bigg)\liminf_{k\to \infty} \int_{{B^N}}\frac{v_k^2}{r^2}\, dx\\
&\geq \bigg(\frac{(N-2)^2}{4}-(N-1)\bigg) \int_{{B^N}}\frac{v^2}{r^2}\, dx.
\end{align*}
Therefore, we conclude that
$$F_\eps(v)\geq \bigg(\frac{(N-2)^2}{4}-(N-1)\bigg) \int_{{B^N}}\frac{v^2}{r^2}\, dx\geq 0,\quad \forall v\in H^1_0(B^N; \R^M),$$
implying by \eqref{rel:endif3} that $u_\eps$ is a minimizer of $E_\eps$ over $\mcA$.
Moreover, $F_\eps(v)=0$ if and only if $v=0$.

\bigskip
\par\noindent{\it Step 5: Conclusion.} We have shown that $u_\eps$ is a global minimizer. Assume that $\tilde u_\eps$ is  another global minimizer of $E_\eps$ over $\mcA$. If $v:=\tilde u_\eps-u_\eps$,
then $v\in H^1_0(B^N; \R^M)$ and by Steps 1 and 4, we have that $0=E_\eps(\tilde u_\eps)-E_\eps(u_\eps)\geq F_\eps(v)\geq 0$, which yields $F_\eps(v)=0$. Step 4 implies that $v=0$, i.e., $\tilde u_\eps=u_\eps$.
 \end{proof}

\begin{remark}
Recall that in the case $M \geq N \geq 7$, J\"{a}ger and Kaul \cite{JagKaul} proved the uniqueness of global minimizer for harmonic map problem
$$
\min_{u \in \mcA_*} \int_{B^N} |\nabla u|^2 \, dx,
$$
where $\mcA_* = \{  u\in H^1(B^N; \Sphere^{M-1}):\,  u(x)=x \textrm{ on } {\partial B^N} = \bS^{N-1} \subset \bS^{M-1} \}$. This can also be seen by the method above as observed in our earlier paper \cite{INSZ_AnnENS}. We give the argument here for readers' convenience: Take a perturbation $v \in H_0^1(B^N, \R^M)$ of the harmonic map $u_*(x) = \frac{x}{|x|}$ such that $|u_*(x) + v(x)|=1$ a.e. in $B^N$. Then, by \cite[Proof of Theorem~5.1]{INSZ_AnnENS},
\[
\int_{B^N} \Big[|\nabla (u_*+v)|^2 - |\nabla u_*|^2\Big]\,dx =\int_{B^N} \Big[|\nabla v|^2- |\nabla u_*|^2|v|^2\Big]\,dx= \int_{B^N}\Big[|\nabla v|^2- (N-1)\frac{|v|^2}{|x|^2}\Big]\,dx.
\]
Using Hardy's inequality in dimension $N$ we arrive at
$$
\int_{B^N} \Big[|\nabla (u_*+v)|^2 - |\nabla u_*|^2\Big]\,dx \geq \left(\frac{(N-2)^2}{4} - (N-1)\right) \int_{B^N} \frac{|v|^2}{|x|^2}\,dx.
$$
The result follows since $N \geq 7$.
\end{remark}

\section*{Acknowledgment.} R.I. acknowledges partial support by the ANR project ANR-14-CE25-0009-01. V.S. acknowledges  support by  the Leverhulme
grant RPG-2014-226. A.Z. was partially supported by a Grant of the Romanian National
Authority for Scientific Research and Innovation, CNCS-UEFISCDI, project number
PN-II-RU-TE-2014-4-0657; by the Basque Government through the BERC
2014-2017 program; and by the Spanish Ministry of Economy and Competitiveness
MINECO: BCAM Severo Ochoa accreditation SEV-2013-0323.

\bibliographystyle{acm}
\bibliography{paris.bib,LiquidCrystals.bib}

\end{document}

config: arxiv


TITLE: Non-Standard analysis and infinitesimal
QUESTION [3 upvotes]: Can someone please explain how Non Standard Analysis is used to justify infinitesimals?
I am not very clear about this but apparently it has something to do with hyperreals.

REPLY [4 votes]: The "infinitesimal", or a quantity which is smaller than any non-zero real number but still non-zero, is a non-rigorous idea that had been a weak link in early calculus, which extensively relied upon it. The Leibniz notation $dy/dx$, for example, is originally supposed to mean "an infinitesimal change in $y$ divided by an infinitesimal change in $x$". Well, it worked, but exactly why eluded the early mathematicians who developed calculus.
The infinitesimal problem was first resolved by looking at limits not as precise evaluations of a function at an infinitesimal value, but a somewhat more complicated definition (the $\epsilon$-$\delta$ definition) that made the concept rigorous and logically sound. However, people have always wondered if it is at all possible to actually extend the real numbers in a way (much like the extension to complex numbers) which makes it satisfy two properties:

It has to preserve certain properties of the real numbers familiar to us, and
It has to contain elements such that any "infinitesimal number" that arises in early calculus corresponds to a precisely-defined element.

Such an extension is called the hyperreals. It justifies the infinitesimal numbers because they are now well-defined: say one wishes to take the derivative of $x^2$ at $x=3$. The definition gives us
$$\lim_{\Delta x\to 0}\frac{(3+\Delta x)^2-3^2}{\Delta x}$$
Doing some quick algebra works this out to $6+\Delta x$ where $\Delta x$ is infinitesimal. In the standard real numbers this word has no precise meaning: it's simply a figure of saying that we happen to be able to manipulate as if it is zero, but without the restrictions (in the simplification above, for instance, we have to divide the numerator and denominator by $\Delta x$, which we cannot do if it is zero). In the hyperreals, however, $6+\Delta x$ is a specific element which is not equal to $6$. However, an equivalence relation called the standard part on a subset of the hyperreals (it is only defined on a subset because the hyperreals actually contain infinite elements as well) gives a homomorphism from the hyperreals to the reals, which can be considered as the hyperreals modulo this equivalence relation. So taking the standard part of $6+\Delta x$ easily shows that the real derivative is $6$. (In fact, any hyperreal construction necessarily admits infinitely many infinitesimals in the equivalence class of any real number.) This would give us property 2 as we discussed above.
What, then, about property 1? It turns out that the hyperreals can preserve many properties of the reals, but not all of them. The ones that they do preserve are the first-order properties. These are properties that can be expressed as a finite formula using quantifiers on the elements, but not on subsets. For instance, "for any two real numbers there is a real number strictly in between" is first-order, hence it's preserved in the hyperreals. On the other hand, the infinitary statement "for any real number $x$ there exists a natural number $n$ such that $x<\underbrace{1+1+\cdots+1}_n$" is not preserved in the hyperreals. (In particular, the field axioms are first-order, so the hyperreals form a field.) This preservation is called the "transfer theorem", a straightforward consequence of the standard construction of the hyperreals as an ultrapower of $\omega$ copies of $\mathbb{R}$ modulo some suitable ultrafilter; Łoś's theorem shows that the first-order properties are all preserved.


TITLE: Function with all but mixed second partial derivatives twice differentiable?
QUESTION [2 upvotes]: Let $f(x,y)$ be a a real valued function on an open subset of $\mathbf{R}^2$ with continuous partial derivatives $\frac{\partial^2 f}{\partial x^2}$ and $\frac{\partial^2}{\partial y^2}$. Is $f$ twice differentiable?

REPLY [3 votes]: Here is an explicit counterexample in $B_1$,
$$f(x,y) = x\, y \, \log(-\log(x^2+y^2))$$
Both $\partial_{xx} f$ and $\partial_{yy} f$ are continuous but $f$ is not twice differentiable at the origin.


TITLE: Choosing 4 integers in increasing order from a set of n integers
QUESTION [0 upvotes]: Suppose there is a set with n distinct integers, suppose I would like to choose a tuple of 4 integers with the following restriction

$(a,b,c,d)$
$d>c$, $c>b$, $c>a$

If in this case, could I say the number of tuples are $n\choose 4$ $\cdot 2$ because for every such tuple, there are 2 ways to arrange first 2 elements. Please correct me if I am wrong

REPLY [1 votes]: As has been stated in the comments, your approach and result are correct.


TITLE: Application of Rolle's theorem? Establish existence of $c\in(a,b)$ such that $f(c)+f'(c)=f(c)f'(c)$
QUESTION [36 upvotes]: Following my question on Meta, I post this as a new question; it was originally asked by Gmgfg, and closed due to the lack of context, background, and shown  effort.

Let $f\colon[a,b]\to\mathbb{R}$ be a function continuous on $[a,b]$ and differentiable on $(a,b)$ with $f(a)=f(b)=0$. Show that there exists $c\in(a,b)$ such that $f'(c)+f(c)=f(c)f'(c)$.

Now, given the assumptions this looks like it should be a straightforward application of Rolle's theorem. However:

as far as I can tell there is no simple auxiliary function $\Phi$ such that $\Phi'= f+f'-f\cdot f'$ to which one could apply Rolle's theorem. (If there is one, I failed to find it.)
indeed, while the RHS could come from $\left(\frac{f^2}{2}\right)'$; it's mostly the LHS which looks difficult to handle (and there is no clear advantage I can see in introducing an antiderivative $F$ of $f$ to have it be $(F+f)'$).

I verified the statement for some functions I could think of, such as $x\mapsto x(1-x)$ and $x\mapsto \sin \pi x$ on $[a,b]=[0,1]$. So, in that regard, it seems to hold at least against basic sanity checks.
However, what bothers me is the lack of "homogeneity." Usually, in things like that I'm used to saying that "without loss of generality, one can assume $[a,b]=[0,1]$." It does not appear to be the case here: if, given $f$, one defines $g\colon[0,1]\to\mathbb{R}$ by $g(x) = f((b-a)x+a)$, finding $c\in(0,1)$ such that $g'(c)+g(c)=g(c)g'(c)$ does not directly yield $c'$ such that $f'(c')+f(c')=f(c')f'(c')$ (but rather would give $c'$ such that $f'(c')+(b-a)f(c')=f(c')f'(c')$, if I'm not mistaken).

So, in short: is it true? And how to prove or disprove it — I'm at a loss, and embarrassed about it.

REPLY [37 votes]: Let
$$ g(x) = f(x)e^{x-f(x)}. $$
Then $g(a) = g(b) = 0$ since $f(a)=f(b)=0$, so $g'(c) = 0$ for some $c\in(a,b)$. However,
$$ g'(x) = f'(x)e^{x-f(x)} + f(x)e^{x-f(x)}(1-f'(x)) = e^{x-f(x)}(f'(x)+f(x)-f(x)f'(x)),$$
so $g'(c) = 0\implies f'(c)+f(c)-f(c)f'(c) = 0$, i.e. $f'(c)+f(c) = f(c)f'(c)$.


text:
Iterative methods are widely used in machine learning and stochastic optimization where the objective functions naturally induce \textit{stochastic processes}. For example, when an algorithm uses stochastic gradient descent (SGD) updates, the loss function forms a stochastic process. Therefore, to study the performance of a learning algorithm, it usually suffices to understand the \textit{behavior} of the corresponding stochastic process.

There have been many successes in providing theoretical guarantees for various learning algorithms. However, as the learning algorithms nowadays become increasingly complicated, it is more and more challenging to perform clean and tight theoretical analysis. Moreover, due to the lack of general principles for analysis, the existing theoretical studies are usually tailored to specific learning dynamics and hence are difficult to extend to other problems.

A main challenge of the analysis, sometimes known as a \textit{``chicken and egg'' phenomenon}~\cite{harvey2019tight}, lies in the circular entanglement between the improvement of the process and the stochasticity. ``Chicken and egg'' phenomenon is ubiquitous due to the non-linearity of iterative methods (see~\autoref{app:framework}). Consider a toy example: Let $\{X_t\}_{t\in\N}$ be a \textit{non-negative-valued} stochastic process with
\begin{equation}\label{eq:toy example}
X_t=X_{t-1}+\frac{N_t}{t+1}
\end{equation}
where $X_0=0.1$, $\Exp[N_t\, |\, X_{t-1}]=-X_{t-1}$, $\Var[N_t\, |\, X_{t-1}]\leq|X_{t-1}|/t$, and $|X_t|,|N_t|\leq1$ almost surely for all $t\geq1$. Note that the expectation of $N_t$ has a \textit{drifting term} $-X_{t-1}$ that drags the dynamic going to $0$ while the stochasticity of $N_t$ is the \textit{noise term}. To show that $X_t$ indeed converges to $0$, we need to show that the aggregated noise is small. However, the noise at time $t$ (i.e., $N_t$) depends on $X_{t-1}$ and hence this creates a circular entanglement between the improvement of the drifting term (toward $0$) and the bound for the noise term.

To overcome the ``chicken and egg'' problem, previous analysis usually either (i) only looks at the expectation and loses important local information (e.g.,~\cite{chen2019convergence,qiao2019rates,yu2019double,zou2019sufficient}), (ii) directly analyzes the process step-by-step which can cause overcomplication (e.g.,~\cite{AL17,allen2019convergence,richards2019optimal}), or (iii) focuses on a specific problem setting and hence is not easy to generalize.

\subsection{An overview on our framework}\label{sec:framework overview}
We propose an analysis framework with an attempt to address the lack of unifying theoretical analysis and the ``chicken and egg'' phenomenon. The framework is \textit{general} as it gives unifying proofs for various learning problems and different convergence guarantees. The framework is \textit{flexible} as it accommodates specialized techniques. To tackle the ``chicken and egg'' phenomenon, our framework provides a handy \textit{stopping time} toolbox to systematically reduce the effort of disentanglement. In the following, we outline the three-step recipe of the framework and apply on the previous toy example (\autoref{eq:toy example}) to show that $\Pr[\exists t\in\N, X_t>10\log(1/\delta)/t]<\delta$ for every $\delta\in(0,1)$.

\paragraph{Step 1: Recursion analysis.}
Find a recursion that dissects the drifting term and the noise term in the dynamic. Consider the toy example in~\autoref{eq:toy example}, we can choose the following recursion (note that there could be multiple choices of recursions that work). For every fixed $T\in\N$, we have
\begin{equation}\label{eq:toy example recursion}
X_T=\left(1-\frac{1}{T+1}\right)X_{T-1} + \frac{N_T+X_{T-1}}{T+1} \, .
\end{equation}
Furthermore, we have $X_{t} = \frac{1}{t+1}(X_0 + \sum_{t'=1}^{t} (N_{t'} + X_{t'-1}))$ for every $1\leq t\leq T$. Intuitively, $\frac{1}{t+1}X_0$ dominates the dyanmic of $\{X_t\}$ while $\sum_{t'=1}^{t}\left(N_{t'}+X_{t'-1}\right)$ is a minor/noise process.

\paragraph{Step 2: Moment and concentration analysis.}
Bound the moments of the minor term and apply a concentration inequality to show that the minor term is small. Since the minor term generally is still entangled with the process, we use \textit{stopping time} to get a tighter control on the moments. In the above toy example, we can use the stopping time $\tau$ for the event $\{X_t>\Lambda_t\}$, i.e., $\tau$ is the random variable defined as $\tau=\min_{X_t>\Lambda_t}t$, where $\{\Lambda_t\}_{t\in\N}$ is a sequence of non-negative threshold parameters that will be chosen later. By a standard application of concentration inequality, we have
\begin{equation}\label{eq:toy example concentration}
\Pr\left[\max_{1\leq t\leq T}\left|\sum_{t'=1}^{\min\{t,\tau\}}N_{t'}+X_{t'}\right|>4\sqrt{\sum_{t'=1}^T\frac{\Lambda_{t'}}{t'}\log\frac{1}{\delta}}+2\log\frac{1}{\delta}\right]<\delta,\ \ \forall \delta\in(0,1),\ \forall T\in\N \, .
\end{equation}

\paragraph{Step 3: Improvement analysis.}
The concentration bound in Step 2 only works for the \textit{stopped process} of the minor term (i.e., in~\autoref{eq:toy example concentration} the noise is summed up to $\min\{t,\tau\}$ instead of $t$). We use a \textit{pull-out lemma} (see~\autoref{lem:pull out}) which gives a sufficient condition to extend the concentration bound from the stopped process of the original minor term without introducing extra failure probability. Informally, the pull-out condition is of the form:
\begin{equation*}
\text{if the concentration bound holds for the noise process}\Rightarrow\text{the stopping time won't stop.}
\end{equation*}
In the toy example, the pull-out condition is
\begin{equation}\label{eq:toy example pull out condition}
\forall1\leq t\leq T,\ \left|\sum_{t'=1}^{t}N_{t'}+X_{t'}\right|\leq4\sqrt{\sum_{t'=1}^T\frac{\Lambda_{t'}}{t'}\log\frac{1}{\delta}}+2\log\frac{1}{\delta}\Rightarrow X_{t}\leq\Lambda_{t} \, .
\end{equation}
Finally, when we pick $\Lambda_t=\min\{1,10\log(1/\delta)/t\}$,~\autoref{eq:toy example pull out condition} holds. Thus, by the pull-out lemma,~\autoref{eq:toy example concentration} extends to the original minor process and hence combining with~\autoref{eq:toy example recursion} we have $\Pr[\exists t\in [T], X_t>10\log(1/\delta)/t]<\delta$. Since the choice of $T$ is arbitrary here, by the monotone convergence theorem, we have $\Pr[\exists t\in\N, X_t>10\log(1/\delta)/t]<\delta$ as desired.

We present the details and intuitions of our analysis framework in~\autoref{sec:framework}. We would like to stress in advance that the technical tools in our framework are basic and it is the \textit{composition} of these tools that makes the framework powerful and flexible. The goal of this paper is to provide a simple recipe for future analysis on a wide-ranged of learning algorithms.

\subsection{Examples of implementing the framework}\label{sec:examples}
To show the power and the flexibility of our framework, we analyze three different learning dynamics and either improve or match the state-of-the-art bounds. In the following, $T$ and $\delta$ stand for the number of iterations and the failure probability respectively.

\paragraph{A textbook example: SGD for strongly convex functions.}
We start with a well-known textbook learning algorithm: \textit{stochastic gradient descent (SGD)} algorithm for (smooth) strongly convex functions. We first give an expository proof that matches the state-of-the-art $O(T^{-1}(\log\delta^{-1}+\log\log T))$ convergence rate~\cite{hazan2014beyond,rakhlin2012making} for \textit{strong uniform convergence}. To show the flexibility of the framework, we further use the same analysis with a different stopping time to achieve the optimal $O(T^{-1}\log\delta^{-1})$ convergence rate~\cite{jain2019making,harvey2019tight} for the \textit{last iterate convergence}. The previous analysis for the two convergence guarantees were different but here we can analyze the two under a unifying analysis. See~\autoref{sec:sgd} and~\autoref{app:sgd} for details.

\paragraph{A classic non-convex example: Streaming $k$-PCA.}
Next, we move on to a classic dynamic with a non-convex structure: the \textit{streaming $k$-principle component analysis ($k$-PCA)} problem. Given i.i.d. samples from an unknown distribution, the goal is to use $O(d)$ space and output a vector that is close to the top-$k$ eigenspace of the distribution. For simplicity, we focus on the \textit{gap-dependent and local convergence setting}\footnote{There are other standard guarantees such as global convergence, gap-free convergence, exponential convergence etc. We leave it as an interesting future direction to apply our framework on these other guarantees and simplify/unify the previous proofs.}. The previous state-of-the-art analysis~\cite{AL17}\footnote{\cite{AL17} obtains $\log^5$ in their bound in the global gap-free setting. Using their analysis techniques on the local convergence would get $\log^3$ in the bound.} gets $O(T^{-1}(\log\delta^{-1}+\log^3(dT)))$ where $d$ is the dimension of the problem. We apply our framework and get $O(T^{-1}\log\delta^{-1})$ convergence rate. When $k=1$, this is the first high-probability result that matches the information-theoretic lower bound. See~\autoref{sec:pca} and~\autoref{app:pca} for details.

\paragraph{An active learning example: Solving stochastic linear bandit with SGD update.}
Finally, we consider a problem with active dynamic where the updates are adaptive and dependant on the whole history: solving stochastic linear bandit with SGD update \cite{jun2017scalable,korda2015fast}. This problem is not only useful for designing scalable bandit algorithms but also serves as an important intermediate step towards analyzing model-free algorithms in linear parameterized Markov decision processes (MDPs)~\cite{jin2018q,jin2019provably}. We analyze an algorithm where previous technique cannot analyze and improve the state-of-the-art regret from $O(d(T\log^2 (T) \log (T/\delta))^{1/2})$ in \cite{jun2017scalable} to $O(d(T\log^2 (T) \log (1/\delta))^{1/2})$. See~\autoref{sec:ql} and~\autoref{app:ql} for details.

\subsection{Related work}\label{sec:related work}
We focus on analyzing stochastic processes in learning algorithms which broadly appear in theoretical machine learning~\cite{moitra2018algorithmic,settles2009active,sutton2018reinforcement}, optimization theory~\cite{boyd2004convex,hazan2019introduction,shalev2014understanding}, statistical learning theory~\cite{hastie2009elements}, etc. There have been many beautiful results providing theoretical analysis on a wide range of important learning problems, e.g., principal component analysis~\cite{AL17}, non-negative matrix factorization~\cite{arora2016computing,lee1999learning}, topic models~\cite{arora2013practical,arora2012learning}, matrix completion~\cite{hardt2014understanding,jain2013low}, tensor decomposition~\cite{anandkumar2014tensor,ge2015escaping}, neural networks~\cite{allen2018convergence,du2018gradient,jacot2018neural}, continual learning~\cite{parisi2019continual,kirkpatrick2017overcoming}, etc. However, the lack of a unifying framework often makes the progress in analyzing frontier learning dynamics slow and sub-optimal. This paper attempts to propose a general framework for the future studies in new and complex learning dynamics.

The technical ingredients in the framework are standard, simple, and inspired by the stochastic approximation theory~\cite{Kushner_1997} and a recent analysis for streaming PCA~\cite{CW20}. For example, the stopping time technique or the martingale concentration have been widely applied in many other analysis~\cite{AL17,rakhlin2012making}. We emphasize that it is the \textit{composition} of these tools that makes our framework powerful and flexible, and the main contribution of this paper is to propose a general and streamlined recipe that future analysis of new and complex learning algorithms can easily adopt.


TITLE: Understanding the behaviour of an interacting Bose gas
QUESTION [3 upvotes]: From the Bose-Einstein distribution it follows that a non-interacting Bose gas condenses into the Bose-Einstein condensate below a certain critical temperature. What happens when interactions are introduced in a Bose gas is not dealt with in introductory Statistical mechanics courses. Hence my questions are pretty naive and basic.
How is(are) the interaction(s) quantitatively modelled in a Bose gas and how does it change the behaviour (compared to the noninteracting Bose gas) when the temperature is lowered? Is there a way to physical understand the change in the behaviour, if any?
As a minor comment, I have been informed that in case of Fermi gases, the role of interactions shifts the effective mass from $m\to m^*$ and the energy levels (effectively mapping the interacting system to a system of quasiparticles that still obey FD statistics). Does the similar thing happen here too?

REPLY [4 votes]: 1) First note that a non-interacting Bose gas is an idealization. If the gas was truly non-interacting, then it would be impossible to cool by evaporative cooling (or any other method that removes energy and requires the gas to re-equilibrate).
2) In a Bose gas the short range part of the interaction has to be repulsive (otherwise the gas will collape at low temperature). A typical model is a repulsive delta function
$$
V(x_1,x_2)=\frac{4\pi a}{m} \delta(x_1-x_2)
$$
controlled by the s-wave scattering length $a$. Indeed, for a dilute gas this is not a model but a systematic description of the low energy properties.
3) In a non-interacting gas Bose condensation takes place at the Einstein temperature
$$
 T_c = \frac{2\pi}{m}\left(\frac{n}{\zeta(3/2)}\right)^{2/3}.
$$
The leading shift due to (repulsive, $a>0$) interactions is
$$
\Delta T_c \simeq 1.3 an^{1/3} T_c
$$
which, even in a strongly interacting gas like helium, is not a large shift.
4) The systematic study of perturbation theory in $a$ goes back to Bogoliubov. He found, for example, that the dispersion relation of quasi-particle in a Bose condensed fluid is
$$
\epsilon_p = \frac{1}{2m}\sqrt{(p^2+8\pi an)^2-(8\pi an)^2}
$$
which smoothly interpolates between a Goldstone mode at low $p$, and non-interacting atoms at large $p$.
5) And indeed, as remarked below, you can take the interaction in (2) and treat in a mean field approximation. This leads to a non-linear Schroedinger equation (the Gross-Pitaevski equation) for the condensate wave function. This equation can be used to study cloud profiles, collective modes, etc.


TITLE: Double integral - finding the region
QUESTION [2 upvotes]: I need to find the volume of an object bounded by the following planes:
$$z=y+4$$
$$x^{2}+y^{2}=4$$
$$z=0$$
Is this true that the region of interest on the XY plane is the whole circle drew by the second formula above?

REPLY [3 votes]: If you do a little sketch, you get
$$-2\le x\le 2\;,\;\;-\sqrt{4-x^2}\le y\le\sqrt{4-x^2}\;,\;\;0\le z\le y+4\implies$$
$$\int\limits_{-2}^2dx\int\limits_{-\sqrt{4-x^2}}^{\sqrt{4-x^2}}dy\int\limits_0^{y+4}dz=\int\limits_{-2}^2dx\int\limits_{-\sqrt{4-x^2}}^{\sqrt{4-x^2}}(y+4)dy=$$
$$=\int\limits_{-2}^2 8\sqrt{4-x^2}=16\int\limits_{-2}^2\sqrt{1-\left(\frac x2\right)^2}\;dx=32\int\limits_{-\pi/2}^{\pi/2}\cos^2t\,dt$$
The last equality following from the substitution
$$\frac x2=\sin t\;,\;\;dx=2\cos t\,dt\;,\;\;\text{ so }$$
$$32\int\limits_{-\pi/2}^{\pi/2}\cos^2t\,dt=\left.32\frac{t+\sin t\cos t}2\right|_{-\pi/2}^{\pi/2}=16\pi$$


TITLE: On the nearest integer to $\zeta(1-1/B),B \ge 2$
QUESTION [1 upvotes]: Let $B \ge 2$ be integer and $[x]$ denote the nearest integer
to real $x$.
For $2 \le B \le 10^5$ computations with mpmath suggest:
$$ [\zeta(1-1/B)]=-B+1 \qquad (1)$$
Is (1) true for all $B \ge 2$?

REPLY [7 votes]: We can make the error mentioned by Wojowu in his comment explicit by using some results on the Laurent coefficients of the zeta function. There are a few results on this, but I'll just use Theorem 2 of this paper by Berndt (which gives a more general result for the Hurwitz zeta function). The relevant result is
Theorem. (Berndt) Write
$$
\zeta(s) = \frac{1}{s-1} + \sum_{n=0}^\infty \gamma_n (s-1)^n
$$
for the Laurent series expansion of $\zeta$ around $s=1$. Then for $n\geq 1$, we have
$$
\left|\gamma_n\right| \leq \frac{4}{n\pi^n}.
$$
For $B\geq 2$ an integer, write
$$
\zeta\left(1-\frac{1}{B}\right) = -B + \gamma + R(B),
$$
where $\gamma=\gamma_0\approx .577...$ is the usual Euler-Mascheroni constant and
$$
R(B) = \sum_{n=1}^\infty \gamma_n \left(-\frac{1}{B}\right)^n.
$$
Applying the triangle inequality and Berndt's result, we have
$$
\left|R(B)\right| \leq \sum_{n=1}^\infty \frac{4}{n(B\pi)^n} \leq \frac{4}{B}\sum_{n=1}^{\infty} \frac{1}{n\pi^n}.
$$
The sum on the right converges rapidly and is numerically $\approx 0.383...$. Thus
$$
\left| R(B) \right| \leq \frac{2}{B}.
$$
(By a more careful analysis, the constant 2 can be replaced with 1.) For $B \geq 26$, it follows that
$$
\frac{1}{2} < \gamma+ R(B) < 1.
$$
Therefore
$$
\left[ \zeta\left(1-\frac{1}{B} \right) \right] = -B+1
$$
for $B \geq 26$.  The remaining cases can be verified via computer.


TITLE: Double derangement
QUESTION [1 upvotes]: Let $[n]$ be a set with $n$ elements. Two derangements of this set, $\sigma$ and $\psi$ satisfy the next condition.

$\forall i\in [n], \sigma(i)\neq \psi(i)$

This condition makes $\sigma$ also a derangement of $\psi$. A double derangement of $[n]$ is an (ordered) pair {$\sigma;\psi$} such that $\forall i\in[n], \sigma(i)\neq\psi(i)$. How many double derangements can there be?
I have been thinking about this problem since I learned about derangements. I tried to solve this problem by finding a recurrence formula and using inclusion–exclusion principle, but figured out that it is too complicated to solve by this method. Are there any better ideas to solve this problem?

REPLY [0 votes]: I believe you are looking for
https://oeis.org/A000186.
The formula is given by
$a(n) = n!\sum_{k+j<=n} \dfrac{2^j}{j!}k!\binom{-3(k+1)}{n-k-j}$.
I checked with the following code:
from itertools import permutations
import math

def diff(x,y):
    return [a-b for a,b in zip(x,y)]
def derangements(s):
    'All deranged permutations of the integers 0..n-1 inclusive'
    return [perm for perm in permutations(s) if all(x!=0 for x in diff(perm,s)) ]

def t(n):
    c = 0
    ders = derangements(range(n))
    for i in ders:
        for j in ders:
            if all(x!=0 for x in diff(i,j)):
                c+=1
    return c
for i in range(7):
    print(t(i))


text:
\begin{document}

\title{A Dynamic Programming Approach To Length-Limited Huffman Coding}

\author{Mordecai~Golin,~\IEEEmembership{Member,~IEEE,}
        and~Yan~Zhang
\thanks{M. Golin and Y. Zhang are with the
Department of Computer Science \& Engineering,  Hong Kong UST,
ClearWater Bay, Kowloon, Hong Kong. Their research was partially supported by HK RGC CERG grants
HKUST6312/04E and 613105.}

}

\maketitle

\begin{abstract}

The ``state-of-the-art'' in Length Limited Huffman Coding algorithms is the
$\Theta(ND)$-time, $\Theta(N)$-space one of Hirschberg and Larmore, where $D\le N$ is the length restriction
on the code.  This is a very clever, very problem specific,
technique.  In this note we show that there is a simple Dynamic-Programming (DP) method
that solves the problem with the same time and space bounds.
The fact that there was an
$\Theta(ND)$ time DP algorithm was previously known; it is a
straightforward DP with the {\em Monge} property  (which permits an order
of magnitude speedup).  It was not
 interesting, though, because it also required
$\Theta(ND)$ space.

The main result of this paper is
the {\em  technique} developed for reducing the space. It is quite simple and applicable to many
other problems modeled by DPs with the Monge property.
We illustrate this with examples from web-proxy design and wireless mobile paging.
\end{abstract}

\begin{IEEEkeywords}
Prefix-Free Codes,  Huffman Coding,  Dynamic Programming,  Web-Proxies,  Wireless Paging, the Monge property.
\end{IEEEkeywords}

\IEEEpeerreviewmaketitle

\section{Introduction}
\label{sec:intro}
\IEEEPARstart{O}{}ptimal prefix-free coding, or \emph{Huffman coding}, is a standard compression technique.
Given an \emph{encoding alphabet}
$\Sigma = \{\sigma_1, \ldots, \sigma_r\}$, a {\em code} is just a set of words in $\Sigma^*.$
Given $n$ probabilities or nonnegative frequencies $\{p_i : 1 \le i \le n\}$, and associated code $\{w_1,w_2,\ldots, w_n\}$
the {\em cost} of the code is $\sum_{i=1}^n p_i |w_i|$ where $|w_i|$ denotes the length of
$w_i$.
A code is {\em prefix-free} if no codeword $w_i$ is a prefix of any other codeword $w_j$.
An {\em optimal} prefix-free code for $\{p_i : 1 \le i \le n\}$
is a prefix-free code that minimizes its cost among all prefix-free codes.

In \cite{Hu52}, Huffman gave the now classical $O(n \log n)$ time algorithm for solving this problem.
If the $p_i$'s are given in sorted order, Huffman's algorithm can be improved to $O(n)$ time
\cite{Va76}.
In this note we will always assume that the $p_i$'s are presorted and  that $p_1 \le p_2 \le \ldots \le p_n$.

In some applications, it is desirable that the length of all code words are bounded by a constant,
i.e., $|w_i| \le D$ where $D$ is given. The problem of finding the minimal cost prefix-free code among
all codes satisfying this length constraint is
the \emph{length-limited Huffman coding} (LLHC) problem, which we will consider here.
Fig.~\ref{fig:seq_ex} gives an example of inputs for which the Huffman code is {\em not} the same as the length-limited Huffman code.

The first algorithm for LLHC was due to Karp \cite{Ka61} in 1961; his
algorithm is based on integer linear programming (ILP), which, using standard ILP solving techniques,
leads to an exponential time algorithm.
Gilbert \cite{Gi71}
in 1971 was interested in this problem because of the issue of inaccurately known sources;
since the probabilities $p_i$'s are not known precisely, a set of codes with limited length will,
in some sense, be ``safe''.  The algorithm presented in \cite{Gi71} was an enumeration one
and therefore also runs in exponential time.
In 1972 Hu and Tan \cite{HuTa72} developed an $O(nD 2^D)$
time Dynamic Programming (DP) algorithm.
The first polynomial time algorithm, running in $O(n^2 D)$ time and using
$O(n^2 D)$ space, was presented by Garey in 1974 \cite{Ga74}.
Garey's algorithm was based on a DP formulation similar to that developed by
Knuth for deriving optimal binary search trees in \cite{Kn71} and hence only works for
binary encoding alphabets.
A decade later, Larmore \cite{La87} gave an algorithm running in $O(n^{3/2} D \log^{1/2} n)$ time
and using $O(n^{3/2} D \log^{-1/2} n)$ space.
This algorithm is a hybrid of \cite{HuTa72} and \cite{Ga74}, and therefore also only works for the
binary case.
This was finally improved by
Larmore and Hirschberg \cite{LaHi90} who gave a totally different algorithm
running in $O(nD)$ time and using $O(n)$ space.
In that paper, the authors first transform the length-limited Huffman coding problem to the
\emph{Coin Collector's} problem, a special type of Knapsack problem, and then, solve
the Coin Collector's problem by  what they name the \emph{Package-Merge} algorithm.
Their result is a very clever special case algorithm developed for this specific problem.

Theoretically, Larmore and Hirschberg's result was later superseded for the case\footnote{
$f(n) = \omega(g(n))$ if $\exists N,c > 0$ such that $\forall n >N$, $f(n) \ge g(n).$}
$D = \omega(\log n)$
by two algorithms based on
the \emph{parametric search} paradigm \cite{Me83}.
The algorithm by Aggarwal, Schieber and Tokuyama \cite{AgSc94+} runs in
$O(n \sqrt{D \log n} + n \log n)$ time and $O(n)$ space.  A later improvement by Schieber \cite{Sc98} runs
in $n 2^{O\left(\sqrt{\log D \log\log n}\right)}$ time and uses $O(n)$ space.

These algorithms are very complicated, though, and even for
$D = \omega(\log n)$, the Larmore-Hirschberg one is the one used in practice
\cite{TuMo96,TuMo95}.
For completeness, we point out that the algorithms of \cite{LaHi90,AgSc94+,Sc98} are all only
claimed for the binary ($r = 2$) case but they can be extended to work for the
non-binary ($r > 2)$ case using observations similar to those we provide in Appendix \ref{app:LLHCDP}
for the derivation of a DP for the generic $r$-ary LLHC problem.

Shortly after \cite{LaHi90} appeared, Larmore and Przytycka
\cite{LaPr91,larmore:1163}, in the context of parallel programming,
gave a simple dynamic programming
formulation for the binary Huffman coding problem.  Although their DP
was for regular Huffman coding and not the LLHC problem, we will see that it is quite
easy to modify their DP to model the LLHC problem.
It is then straightforward to show that their formulation also permits constructing
the optimal tree in $\Theta(nD)$ time by constructing a size $\Theta(nD)$ DP table.
This is done is Section \ref{Sec:DP}.
This straight DP
approach would not be as good as the Larmore-Hirschberg one, though, because, like
many DP algorithms, it requires maintaining the entire DP table to
permit backtracking to construct the solution, which would require $\Theta(nD)$ space. The main
result of this  note is the development of
 a simple technique (section \ref{Sec:space})
that permits reducing the DP space consumption down to $O(n)$,
thus matching the Larmore-Hirschberg performance with a straightforward DP model.
Our technique is not restricted to Length-Limited coding.  It can be used to reduce space
from $O(nD)$ to $O(n+D)$ in a variety of $O(nD)$ time DPs in the literature.  In Section
\ref{sec:further}
we illustrate with examples from the D-median on a line problem
(placing web proxies on a linear topology network)
\cite{Wo00} and
wireless paging \cite{973511}.

\section{The Dynamic Programming Formulation}
\label{Sec:DP}

Set $S_0 = 0$ and $S_m= \sum_{i=1}^m p_i$ for $1 \le m \le n$.
Larmore and Przytycka \cite{larmore:1163}
formulated the binary Huffman coding problem as a DP (\ref{eq:LM})
where
$H(0) = 0$ and for $0 < i < n$:
\begin{equation}
\label{eq:LM}
H(i) = \min_{\max\{0,2i-n\} \le j < i} \left(H(j) + S_{2i-j}\right).
\end{equation}
In this DP, $H(n-1)$ is the cost of the optimal
Huffman code.  Another version of this DP, generalized for unequal-cost
binary coding alphabets, also appeared in \cite{BrGo02+}.

It is straightforward to modify
(\ref{eq:LM}) to model the binary LLHC problem.  The resulting DP is

\begin{equation}
\label{eq:DP_Def}
H(d,i) =
\left\{
\begin{array}{@{\,}ll}
0 & d = 0, i = 0 \\
\infty & d = 0, 0 < i < n \\
\displaystyle{\min_{0 \le j \le i} \left(H(d-1,j) + c^{(d)}_{i,j}\right)} & d > 0, 0 \le i < n
\end{array}
\right.
\end{equation}

where $H(D,n-1)$ will denote the cost of the optimal length-limited Huffman code and
\begin{equation}
\label{eqn:c}
c_{i,j}^{(d)} =
\left\{
\begin{array}{@{\,}ll}
0 & i = j = 0 \\
S_{2i-j} & \max\{0,2i-n\} \le j < i \\
\infty & \textrm{otherwise.}
\end{array}
\right.
\end{equation}
In the next subsection we will see an interpretation of this DP
(which also provides an interpretation of (\ref{eq:LM})).
In order to make this note self-contained, a complete derivation of the
DP for the $r$-ary alphabet case is provided in Appendix \ref{app:LLHCDP}.

As far as running time is concerned,
(\ref{eq:LM}) appears to a-priori require
$O(n^2)$ time to fill in its corresponding DP table.
\cite{larmore:1163} used the inherent concavity of $S_m$ to reduce this time
down to $O(n)$ by transforming the problem to an instance of
the
Concave Least Weight Subsequence (CLWS) problem and using one of the
known $O(n)$ time algorithms, e.g., \cite{Wi88}, for solving that problem.

Similarly, (\ref{eq:DP_Def}) appears to a-priori require
$\Theta(n^2D)$ time to fill in its DP table.  We will see that we may again
use the concavity of $S_m$ to reduce this down by an order of
magnitude, to $O(nD)$ by using the SMAWK algorithm \cite{AgKl87+}
for finding row-minima of matrices as a subroutine.
Unlike the CLWS algorithms, the SMAWK one is very simple to code and
very efficient implementations are available in different packages, e.g.,
\cite{Ep02,NeoBio}.  In the conclusion to this note, after the application
of the technique becomes understandable, we will explain
why \cite{larmore:1163} needed to use the more complicated CLWS
routine to solve the basic DP while we can use the simpler SMAWK one.

The $O(nD)$ DP algorithm for solving the LLCH problem,
while seemingly never explicitly stated in the literature, was known as folklore.
Even though it is much simpler to implement than the $O(nD)$
Larmore and Hirschberg \cite{LaHi90}
Package-Merge algorithm it suffers from the drawback of requiring $\Theta(nD)$ space.
The main contribution of this note is the observation that its space can be reduced down to
$O(n+D)$ making it comparable with Package-Merge.
Note that since, for the LLHC problem we may trivially assume $D \le n$, this implies
a space requirement of $O(n)$.
Furthermore, our
space improvement will work not only for the LLHC problem but for all
DPs in form (\ref{eq:DP_Def}) where the $c^{(d)}_{i,j}$ satisfy a particular property.

\subsection{The meaning of The DP}

\begin{figure*}[!htbp]
\centering
\includegraphics{tree.eps}
\caption{Two trees and their corresponding sequences $\cali$ and codes.
The left tree has sequence $\cali_1 = (0,1,3,4,5,6).$ The right tree has
sequence $\cali_2 = (0,2,4,5,6).$ Note that, for both trees, $2 i_k - i_{k-1}$
is the number of leaves below level $k.$
For input frequencies $(p_1,\ldots,p_7) = (1,1,2,2,2,4,5,9)$.
The left tree is an  optimal Huffman code while
the right tree is an  optimal length-limited Huffman
code for $D = 4$.
Note that we allow padding sequences with initial {\bf 0}s, so the right tree could
also be represented by sequences $(0,0,2,4,5,6),$ $(0,0,0,2,4,5,6),$ etc..
}
\label{fig:seq_ex}
\end{figure*}

\begin{figure*}[!htbp]
\centering
\includegraphics{table.eps}
\caption{Solving the DP in equation  \ref{eq:DP_Def} for
$(p_1,\ldots,p_7) = (1,1,2,2,2,4,5,9)$ with $D=4.$
$H(d,i)$ is the value defined by (\ref{eq:DP_Def});  $J(d,i)$
is the index $j$ for which the value $H(d,i)$ in (\ref{eq:DP_Def})
is achieved. The circled entries yield the sequence $(0,2,4,5,6)$  (the $6$
comes from the fact that we are calculating $H(4,{\bf 6})$) which
is exactly the sequence $\cali_2$  from  Figure
\ref{fig:seq_ex}.  The righthand tree in Figure
\ref{fig:seq_ex} is therefore an  optimal length-limited Huffman
code for $D = 4$.
}
\label{fig:dp_table}
\end{figure*}

We quickly sketch the meaning of the DP (\ref{eq:DP_Def}) for the binary case.
Figures \ref{fig:seq_ex} and \ref{fig:dp_table} illustrate this sketch.
We note that in order to stress the parts important to our analysis,
our formalism is a bit
different than \cite{larmore:1163,BrGo02+}.  A complete derivation of the
DP for the $r$-ary case with the appropriate general versions
of the lemmas and observations stated below along with their proofs,
is provided in Appendix \ref{app:LLHCDP}.

It is standard that there is a $1-1$ correspondence between
binary prefix-free code with $n$ words and binary tree with $n$
leaves.
The set of edges from an internal node to its children are labeled by a {\bf 0} or {\bf 1}.
Each leaf corresponds to a code word, which is the concatenation of the characters on the
root-to-leaf path.
The cost of the code equals the \emph{weighted external path length} of the tree. So
we are really interested in finding a binary tree with minimum weighted external path length.

Denote the height of the tree by $h.$
The bottommost leaves are on level $0$; the root on level $h$.
Optimal assignments of the $p_i$'s to the leaves always assign
smaller valued $p_i$'s to leaves at lower levels.

A node in a binary tree is \emph{complete} if it has two
children and  a tree is \emph{complete} if all of its internal nodes are complete.
A min-cost tree must be complete, so we restrict ourselves to complete trees.
A complete tree $T$ of height $h$ can be completely
represented by a sequence $(i_0, i_1, \ldots, i_h)$,
where $i_k$ denotes the number of internal nodes at levels $\le k$.  Note that,
by definition, $i_0 = 0$,  $i_h = n-1$.  Also note that every level must contain at least
one internal node so $i_0 < i_1 < \cdots < i_h$.  Finally,
it is straightforward (see Appendix \ref{app:LLHCDP}) to show that the total number of leaves
on level $< k$ is
$2 i_k - i_{k-1}$, so $2 i_k - i_{k-1} \le n$ for all $k$.  For technical reasons,
because we will be dealing with trees having height \emph{at most} $h$
(but not necessarily {\em equal to} $h$),
we allow initial padding of the sequence by
{\bf 0}s so a sequence representing a tree will be of the form
$(i_0, i_1, \ldots, i_h)$ that has the following properties

\begin{definition}
\label{def:valid.list}
Sequence $(i_0, i_1, \ldots, i_h)$ is \emph{valid} if
\begin{itemize}
\item $\exists t> 0 $ such that
$i_0 = i_1 = \cdots = i_t = 0$,
\item $0 < i_{t+1} <  i_{t+2} < \cdots < i_h \le n-1$
\item $2 i_k - i_{k-1} \le n$ for all
$1 \le k \le h$.
\end{itemize}
A sequence is \emph{complete} if it is valid and $i_h = n-1$.
\end{definition}

We can rewrite the cost function for a tree in terms of its complete sequence.
\begin{lemma}
\label{lem:eqn:cost.tree}
If complete sequence $(i_0, i_1, \ldots, i_h)$ represents a tree, then the cost of the tree is
$\sum_{k=1}^h S_{2  i_k - i_{k-1}}$.
\end{lemma}
(Note that padding  complete sequences with initial {\bf 0}s does not change the cost of the sequence.)

We may mechanically extend this cost function to {\em all}  valid
 sequences as follows.
\begin{definition}
\label{def:cost}
For valid  $\cali = (i_0, i_1, \ldots, i_h)$, set
$$cost(\mathcal{I}) = \sum_{k=1}^h S_{2 i_k - i_{k-1}}.$$
$\cali$ is {\em optimal} if
$cost(\cali) = \min_{\cali'} cost(\cali')$ where the minimum is taken over all length $h$
sequences
$\cali' = (i'_0, i'_1, \ldots, i'_h)$ with $i'_h = i_h,$ i.e., all sequences of the same length that end with
the same value.
\end{definition}
Our goal is to find optimal trees by using the DP to optimize over valid sequences.  An
immediate issue is that not all complete sequences represent trees,  e.g.,
$\cali = (0,3,4,5)$ is complete for $n=6$ but,  by observation,  does not represent a tree.
The saving fact is  that even though not all complete sequences represent trees, all {\em optimal}
complete sequences represent trees.
\begin{lemma}
\label{lem:sequence.2.tree}
An optimal valid sequence ending in $i_h = n-1$ always represents a tree.
\end{lemma}

Thus,  to solve the LLHC problem of finding an
optimal tree of height $\le D$, we only  need to find an optimal valid sequence of length $h=D$
ending with $i_D=n-1$
(reconstructing the tree from the sequence can be done in $O(n)$ time).
In the  DP defined  by equations (\ref{eq:DP_Def}) and (\ref{eqn:c}), $H(d,j)$ clearly models
the recurrence for finding an optimal valid sequence $(i_0, i_1, \ldots, i_d)$ of length $d$ with $i_d=j$ so
this DP  solves the problem.

Note that, a-priori,   filling in the DP table $H(\cdot,\cdot)$ one entry at a time seems to  require
$O(n^2D)$ time.  We will now sketch the standard way of reducing this time down to $O(nD).$  Before doing
so we must  distinguish between the {\em value} problem and the {\em construction problem}.  The value
problem would be to calculate the value of $H(D,n-1).$  The construction problem would be to
construct an optimal valid {\em sequence}
$\cali = (I_1,I_2,\ldots,I_{D})$  with $I_{D}=n-1$
and  $cost(\cali) = H(D,n-1).$
This would require  backtracking through the DP table by
setting $I_0=0$, $I_{D} = n-1$ and finding  $I_1, I_2, \ldots I_{D-1}$ such that
\begin{equation}
\label{eq:const}
\forall 0 < d \le D,\quad H(d,I_d) = H(d-1,I_{d-1}) + c^{(d)}_{I_d,I_{d-1}}.
\end{equation}

\subsection{Solving the Value problem in $O(nD)$ time}

\begin{figure*}[!htbp]
\centering
\includegraphics{matrix.eps}
\caption{The matrices used for calculating the DP tables in Fig.~\ref{fig:dp_table}.
The shaded entries are the row minima. The row minima for $M^{(i)}$ are exactly the
row entries in the $H(d,i)$ table in Fig.~\ref{fig:dp_table}. The column indices of
the corresponding row minima are the $J(d,i)$ entries.}
\label{fig:row_ex}
\end{figure*}

\begin{definition}
An $n \times m$ matrix $M$ is
{\em Monge}\footnote{This property is sometimes alternatively defined by:
for  $0 \le i < i' < n$ and $0 \le j < j' <  m$
$M_{i,j} + M_{i',j'} \le M_{i',j} + M_{i,j'}$ but it is well known,
see, e.g., \cite{BuKl96+}, that this is equivalent to (\ref{eq:Monge def}).}

if for $0 \le i < n-1$ and $0 \le j <  m-1$
\begin{equation}
\label{eq:Monge def}
M_{i,j} + M_{i+1,j+1} \le M_{i+1,j} + M_{i,j+1}
\end{equation}
\end{definition}
The Monge property can be thought of as a discrete version of concavity.  It appears implicitly
in many optimization problems for  which it permits speeding up their solutions (\cite{BuKl96+})
provides a nice survey).
One of the classic techniques used is the SMAWK algorithm for finding row-minima.

Given an $n \times m$ matrix $M$, the {\em minimum} of row $i$, $i=1,\ldots,n$ is the entry of row $i$ that
has the smallest value; in case of ties, we take the rightmost entry.
Thus, a solution of the row-minima problem is a collection of indices $j(i)$, $i=1,\ldots,n$ such
that
{\small
\begin{equation}
\label{eq:res}
M_{i,j(i)} = \min_{0 \le j < m} M_{i,j}

\mbox{ and }
\hspace*{.02in}
j(i) = \max \{j \,:\, M_{i,j} = M_{i,j(i)}\}.
\end{equation}
}
Figure \ref{fig:row_ex} gives four examples of Monge matrices and their row minima.

At first glance it seems that we would have to examine all of the $mn$ entries in $M$ to find the
row minima but, \cite{AgKl87+} proved\footnote{Technically,
\cite{AgKl87+} proved their result for a larger class, the {\em totally-monotone matrices}.  But all
applications in the literature seem to be for Monge matrices.}
\begin{lemma} (The SMAWK algorithm \cite{AgKl87+})\\
\label{lem:SMAWK}
Let $M$ be a $n \times m$ Monge matrix such that entry $M_{i,j}$ can be calculated in
$O(1)$ time.
Then the row minima problem on $M$ can be solved in  $O(n + m)$ time.
\end{lemma}

The constant hidden by the $O(\,)$ is very small, around 2, and the algorithm is easy to code,  so it
is quite practical to use.

Note that the SMAWK algorithm doesn't have the time available to build the entire $n \times m$
 matrix.  Instead,  it searches through the matrix in a clever way,  constructing entries as needed.
One standard use of the SMAWK algorithm is  in the speedup of dynamic programs that have Monge properties.

\begin{definition}
A  DP in the form (\ref{eq:DP_Def}) is  {\em Monge} if,
for all $1 \le d \le D$ and $0 \le j \le i < n$,
\begin{equation}
\label{eq:def_Quad}
c^{(d)}_{i,j} + c^{(d)}_{i+1,j+1} \le c^{(d)}_{i+1,j} + c^{(d)}_{i,j+1}
\end{equation}
\end{definition}

{\small \em Note: In many DP applications, it is possible that for some $i,j$,
$c^{(d)}_{i,j} = \infty$.  The inequality in (\ref{eq:def_Quad}) treats
$\infty$ in the natural way, e.g., for any constant $c$;  $c < \infty$ and   $c + \infty = \infty$. Also,
$\infty + \infty = \infty$. The SMAWK algorithm permits the use of $\infty$ in this way.}

Now suppose that a DP  defined by (\ref{eq:DP_Def}) is Monge.
For $d=1,2,\ldots,D$ define matrix $M^{(d)}$ by
$$
M^{(d)}_{i,j} =
\left\{
\begin{array}{ll}
H(d-1,j) + c^{(d)}_{i,j} & \mbox{if $0 \le j \le i < n$} \\
\infty & \mbox{otherwise}
\end{array}
\right.
$$
Then, from (\ref{eq:def_Quad}), we have
{\footnotesize
\begin{eqnarray*}
\hspace*{-.15in} M^{(d)}_{i,j} + M^{(d)}_{i+1,j+1}
&=& H(d-1,j) + H(d-1,j+1) + c^{(d)}_{i,j} +c^{(d)}_{i+1,j+1}\\
&\le& H(d-1,j) + H(d-1,j+1) + c^{(d)}_{i+1,j} +c^{(d)}_{i,j+1}\\
&=& M^{(d)}_{i+1,j} + M^{(d)}_{i,j+1}
\end{eqnarray*}
}
and  $M^{(d)}$ is Monge.
Note that
\begin{eqnarray*}
H(d,i) &=&  \min_{0 \le j \le i} \left(H(d-1,j) + c^{(d)}_{i,j}\right)\\
&=& \min_{0 \le j \le i} M^{(d)}_{i,j}  =  \min_{0 \le j \le N} M^{(d)}_{i,j}.
\end{eqnarray*}
So, $H(d,i)$ are just the row-minima of $M^{(d)}$.
See Figure \ref{fig:row_ex}.
Since $M^{(d)}$ is Monge, we can use the SMAWK algorithm to, in $O(n)$ time, find {\em all} of its
row minima at one time.
More specifically,  let $J(d,i)$ and
$M^{(d)}_{i,J(d,i)}$ be the corresponding
values (\ref{eq:res}) returned when running SMAWK$(M^{(d)}).$
Then the algorithm for filling in the table is just to iteratively
run down the rows of the table, using SMAWK to fill in each  row by using knowledge of
the previous row:

\begin{figure}[h]

\par\noindent\underline{\tt Fill\_Table}\\

For $d=1$ to $D-1$\\
\hspace*{.1in} SMAWK $(M^{(d)})$\\
\hspace*{.1in} $\forall 0 \le i < n$ set   $H(d,i) = M^{(d)}_{i,J(d,i)}$
\caption{The $O(nD)$ algorithm for the value problem.}
\label{fig:alg_FT}
\end{figure}

Note that this algorithm uses $\Theta(nD)$ time, since, for each fixed $d$, the SMAWK algorithm
only uses $O(n)$ time.
Also note that if we're only interested in the final row, then the algorithm uses only
$O(n)$ space, since once row $d$ has been calculated,  the values from row $d-1$ can be thrown away.

We now return to the LLHC problem and show that it can be plugged into the above machinery.
\begin{lemma}
\label{lem:c.qi}
The $c_{i,j}^{(d)}$ defined in (\ref{eqn:c}) satisfy Monge property
(\ref{eq:def_Quad}).
\end{lemma}
\begin{IEEEproof}
If $i=j=0$ the righthand side of (\ref{eq:def_Quad}) is $\infty$, so
(\ref{eq:def_Quad}) is satisfied.

If $j+1 = i$ or $2(i+1)-n > j$, the righthand side of (\ref{eq:def_Quad}) is $\infty$, so
(\ref{eq:def_Quad}) is satisfied.

If $j+1 < i$ and $2(i+1)-n \le j$, (\ref{eq:def_Quad}) can be rewritten as
\begin{equation}
\label{eqn:c.qi.rewrite}
S_{2i-j} + S_{2(i+1)-(j+1)} \le  S_{2i-(j+1)} + S_{2(i+1)-j}
\end{equation}
It is easy to verify
\begin{eqnarray*}
S_{2i-j} + S_{2(i+1)-(j+1)} - S_{2i-(j+1)} - S_{2(i+1)-j}&&\\
\hspace*{.2in} = p_{2i-j} - p_{2i-j+2} \le 0&&
\end{eqnarray*}

Hence, (\ref{eqn:c.qi.rewrite}) holds.
\end{IEEEproof}

Thus,  from the discussion above,  we can find all of the $H(d,i)$ in $\Theta(nD)$ time.
In particular, $H(D,n-1)$ will be the cost of the optimal tree with height {\em at most} $D$
which is the required cost of the optimum $D$-limited code.

We have thus seen how to solve the value problem in $O(nD)$ time.
The difficulty is that {\em constructing}  the optimal tree associated with
$H(D,n-1)$  would require
finding the associated optimal valid sequence with $i_{D} = n-1.$  This would require solving
the construction problem by finding all indices $I_d$ in
(\ref{eq:const}).  The standard way of  solving this problem is to
 maintain  an   array storing the
$J(d,i)$ values returned by the algorithm.

Starting from   $H(D,n-1)$ and backtrack through the $j(\cdot,\cdot)$ array, constructing the
corresponding sequence by setting $I_{D} = n-1$ and $I_{d-1} = j(d,I_d).$
Unfortunately, this requires maintaining a size $\Theta(nD)$ auxiliary array,  which
requires too much space.

\section{ Solving the Construction Problem in $O(nD)$ time and $O(n+D)$ space}
\label{Sec:space}

\begin{figure}[!htbp]
\centering
\includegraphics{graph.eps}
\caption{The dropping-level graph associated with the example from Figures
\ref{fig:dp_table} and  \ref{fig:row_ex}.  The bold edges are the
minimum cost path from $(0,0)$ to $(6,4).$  Note that the $i$ coordinates
of the path are $(0,2,4,5,6)$ which is {\em exactly} the sequence of $J(d,i)$'s corresponding to optimal solution of the problem,  which is also the sequence
corresponding to the optimal tree.
}
\label{fig:drop_ex}
\end{figure}

Let $V$ be the grid nodes $(d,i)$ with $0 \le d \le D$ and $0 \le i <n$.
Consider the directed graph $G=(V,E)$ in which $(d,i)$ points  to all nodes immediately below it and to its  right, i.e.,
$$E = \{(\, (d,j),\, (d+1,i)\,) \mid  (d,j) \in V,\,  d < D,\,  j\le i \}$$
See Figure \ref{fig:drop_ex}.
Such graphs are sometimes called {\em dropping level-graphs} \cite{MuRa82}.
Now assign edge $(\, (d-1,j),\, (d,i)\,)$ the weight $c^{(d)}_{i,j}$.
The length of a path in $G$ will just be the sum of the weights of the edges in the path.
The important observation is that $H(d,i)$ in  DP (\ref{eq:DP_Def}) is simply

the length of the min-cost path from $(0,0)$ to $(d,i)$ in this weighted $G$. More
specifically,  the value
problem is to find the {\em length of a shortest path} and the construction problem is
to find an  {\em actual shortest path.}

A-priori,  finding such  a path seems to  require $O(nD)$ space.  There are two different
algorithms in the literature for reducing the space down to $O(n+D)$ in  related problems.

The first was for finding a
maximum common subsequence of two sequences. This reduced down to the problem of finding
a {\em max-length} path in something very similar to a dropping level-graph in which each vertex has
bounded indegree and bounded outdegree.  Hirschberg \cite{Hi75} developed an $\Theta(nD)$ time, $\Theta(n+D)$ space algorithm for this problem.  His algorithm was very influential in the bioinformatics community and its technique is incorporated into many later algorithms e.,g \cite{ChHa94+,MyMi88}. The techniques's performance  is very dependent upon the bounded degree of the vertices,  which is not true in our case.

The second, due to  Munro and Ramirez \cite{MuRa82},  was exactly for the problem of constructing
min-cost  paths in full dropping level-graphs.  Their algorithm ran in $\Theta(n^2D)$ time and $\Theta(n+D)$ space.
Their $\Theta(n^2D)$ time is too expensive for us.  We will now see how to reduce this down to
$\Theta(nD)$ using the Monge speedup while still maintaining the $\Theta(n+D)$ space.

The general problem  will be
to construct an optimal $u$-$w$ path in $G$ where $u=(d_u,i_u)$ is above and not to the left of  $w=(d_w,i_w),$
i.e., $d_u < d_w$ and $i_u \le i_w.$
Let $G(u,w)$ be the subgrid with upper-left corner  $u$ and lower-right corner
$w$ (with associated induced edges from $G$).
First note that, because $G$ is a dropping level-graph,
any optimal (min or max cost) $u$-$w$ path in $G$ must lie completely in $G(u,w)$.
Both algorithms  \cite{Hi75,MuRa82} start  from the same observation, which is  to  build the
path {\em recursively} i.e.,   by first (a)
 finding a point $v=(\bar d, \bar i)$ halfway (by link distance) on the optimal $u$-$w$ path
 in $G(u,w)$ and then (b) output the recursively constructed optimal $u$-$v$ path  in $G(u,v)$ and
 optimal $v$-$w$ path  in $G(v,w)$.

For dropping level-graphs, if $u=(d_1,i_1)$ and $w=(d_2,i_2)$ then the midlevel must be
$\bar d = \lfloor {(d_1 + d_2)}/2 \rfloor $.  Suppose that we had an algorithm $Mid(u,w)$ that returned
a  point $v=(\bar d, \bar i)$ on a shortest $u$-$w$ path in $G(u,w).$  Then, translated into our notation and
with appropriate termination conditions
 the construction algorithm can be written as:

\vspace*{.5in}

 \begin{figure}[h]

 \par\noindent\underline{{\tt Path}$(u,w)$}\\
 1. If $u=(d,j)$ and $w=(d+1,i)$ then\\
 2. \hspace*{.1in} output edge $(u,w)$\\
 3. Else if $u=(d,i)$ and $w=(d',i)$ then\\
 4.  \hspace*{.1in} Output  vertical path from $u$ to $w$\\
 5. Else\\
 6. \hspace*{.1in}  set $v=Mid(u,w)$\\
 7. \hspace*{.1in} {\tt Path}$(u,v)$; \ {\tt Path}$(v,w)$\\

  \label{fig:alg_Path}
  \caption{The algorithm for constructing a min-cost $u$-$w$ path.}
  \end{figure}

\begin{figure*}[!htbp]
\centering
\includegraphics{recursion.eps}
\caption{An illustration for finding the optimal path.
Here, $D = 8$ and there are $3$ levels of recursions.
The solid circles are the intermediate nodes found by the $Mid(u,v)$ procedures.
The first level of recursion finds the midpoint on level $4;$ the second level, the midpoints on
levels $2$ and $6$; the third the midpoints on levels $1,3,5,7.$  At that point all
subproblems are of height one and easily solvable. Note that each recursive call splits
a problem on a box of height $2^i$ into two problems on disjoint boxes of height $2^{i-1}.$}
\label{fig:path_illus}
\end{figure*}

(Figure \ref{fig:path_illus} illustrates this idea.)
 To solve the original problem we just call {\tt Path}$(u_0,w_0)$
 where $u_0=(0,0)$ and $w_0=(D,n-1).$
 Correctness follows from the fact that at each recursive call,
 the vertical distance $d_w-d_u$ decreases so the recursion must terminate.
 Furthermore, when the recursion terminates,  either (i)
 $u=(d,j)$ and $w=(d+1,j)$ so  the {\em only} $u$-$w$
 path in $G(u,w)$  is the edge $(u,w)$ or (ii) $u=(d,i)$ and $w=(d',i)$
 so  the {\em only} $u$-$w$
 path in $G(u,w)$  is the vertical path going down from $u$ to $w$.

 The efficiency of the resulting algorithm, both in time and space,  will depend upon how efficiently  $v=Mid(u,w)$ can be found.
 Note that with the exception of the calls of type
  $Mid(u,w)$, the  rest of the execution   of {\tt Path}$(u_0,w_0)$ (including all recursive calls)
 only requires  a total of $O(D)$ space,  since each recursive call  uses
 only $O(1)$ space and there are at most $O(D)$ such calls.  Thus,  if  $Mid(u,w)$ can be
  found  using $O(n+D)$ space,  then the entire procedure requires only $O(n+D)$ space.
  This is actually how both
 \cite{Hi75,MuRa82} achieve their  space bounds.  The two algorithms differ in how they calculate $v$.
   Although both their approaches  can be used  for our problem,  we will work with a modified version of that of \cite{MuRa82}, since it will be simpler to explain.

 We now describe how to use the SMAWK algorithm to find
 $Mid(u_0,w_0)$ in $O(nD)$ time and $O(n)$ space.  The  extension to general $Mid(u,w)$
 will follow later.  Recall that the procedure {\tt Fill\_Table} from Figure  \ref{fig:alg_FT}
  used the  fact that
 $H(\cdot,\cdot)$ was Monge and  the SMAWK algorithm
  to iteratively fill in the rows $H(d,\cdot)$,  for $d=1,2,\ldots,D.$  Given row $H(d-1,\cdot)$,
 the  procedure
 calculated $H(d,\cdot)$ in $O(n)$ time using SMAWK,   and then threw away $H(d-1,\cdot)$.

Consider an arbitrary node $(d,i)$ on level $d > \bar d.$
The shortest path from $u_0$ to $(d,i)$ must pass through {\em some} node on level $\bar d$.
We now modify {\tt Fill\_Table}
 to ``remember''  this  node.
More specifically, our algorithm will  calculate  auxiliary data  $pred(d,i)$.
\begin{itemize}
\item For $d < \bar d$, $pred(d,i)$ will be undefined.
\item For $d \ge \bar d$, $pred(d,i)$ will be an
index $j$ such that node $(\bar d,j)$  appears
on some  shortest path from $u_0$ to $(d,i)$.
\end{itemize}

So, when the procedure terminates,  $v=(\bar d,\,  pred(d, n-1))$ will be $Mid(u_0,w_0).$

By definition, on level $\bar d$,  we have  $pred(\bar d,i) =i.$

For $d > \bar d$ suppose $(d-1,j')$ is the immediate predecessor of $(d,i)$ on the shortest path
from $u_0$ to $(d,i)$.  Then (i) a  shortest path from $u_0$ to $(d-1,j')$ followed by
(ii) the edge from
$(d-1,j')$ to
$(d,i)$ is  (iii) a shortest
path from $u_0$ to  $(d,i)$; we may therefore set $pred(d,i) = pred(d-1,j').$

We can use this observation to  modify {\tt Fill\_Table}
to calculate  the $pred(d,\cdot)$ information.

\begin{figure}[h]
\par \noindent \underline{$Mid(u_0,w_0)$}\\
For $d=1$ to $\bar d$\\
\hspace*{.1in} SMAWK $(M^{(d)})$\\
\hspace*{.1in} $\forall 0 \le i < n$ set   $H(d,i) = M^{(d)}_{i,J(d,i)}$\\
$\forall 0 \le i < n$ set $pred(\bar d,i) = i;$\\
For $d=\bar d +1 $ to $D$\\
\hspace*{.1in} SMAWK $(M^{(d)})$\\
\hspace*{.1in} $\forall 0 \le i < n$,\ set   $H(d,i) = M^{(d)}_{i,J(d,i)}$\\
\hspace*{.1in} $\forall 0 \le i < n$,\  set $pred(d,i)  = pred\left(d-1,j(d,i)\right)$

\caption{Returns the midpoint, by link distance, on min-cost $u_0$-$w_0$ path.}
\label{fig:Midalg}
\end{figure}
  Note that $Mid(u_0,w_0)$ can throw away all
of the values $pred(d-1,\cdot)$ and $H(d-1,\cdot)$ after the values $pred(d,\cdot)$ and $H(d,\cdot)$ have been calculated, so it only uses  $O(n)$ space. Similarly to the analysis of
{\tt Fill\_Table}, it uses only $O(nD)$ time since each call to the SMAWK algorithm uses only
$O(n)$ time.

 So far,  we have only shown how to find $v=Mid(u_0,w_0)$.
 Note that the {\em only assumptions} we used were that $H(\cdot,\cdot)$ satisfies
 DP (\ref{eq:DP_Def}) and  is Monge, i.e.,
 the $c^{(d)}_{i,j}$ satisfy  (\ref{eq:def_Quad}).

Now suppose that we are given
$$u=(d_u,i_u),\   w=(d_w,i_w) \mbox{ with $d_u < d_w$ and $i_u \le i_w$}.$$
$G(u,w)$ is  a dropping level-graph on its own nodes so the cost of the shortest path from $u$ to any node $(d_u+d,i_u +i) \in G(u,w)$ is $\tilde H(d,i)$
defined by
{\small

\begin{equation}
\label{eq:DP_Def_2}
\tilde H(d,i) =
\left\{
\begin{array}{ll}
0 & \mbox{if $d = 0$, $i = 0$} \\
\infty & \mbox{if $d = 0$, $0 < i <N$} \\
\displaystyle{\min_{0 \le j \le i} \left(\tilde H(d-1,j) + \tilde c^{(d)}_{i,j}\right)} & \mbox{if $d>0$, $0 < N$}
\end{array}
\right.
\end{equation}
}
where $N=i_w-i_u +1$ and
$\tilde c^{(d)}_{i,j} = c^{(d)}_{i_u+ i,i_u + j}.$
Note that this new DP is exactly in the same form as (\ref{eq:DP_Def}),  just with a different $n$
and shifted $c^{(d)}_{i,j}.$
Since the original $c^{(d)}_{i,j}$ satisfy (\ref{eq:def_Quad}),
so do the $\tilde c^{(d)}_{i,j}$.  Thus (\ref{eq:DP_Def_2}) with the
$\tilde c^{(d)}_{i,j}$ is Monge as well.

Therefore,  we can
run exactly the same algorithm written in Figure \ref{fig:Midalg}  to find
the midpoint
 $v=(\bar d, \bar i)= Mid(u,w),$
 of the min-cost $u$-$w$ path in $O((d_w -d_u)N)$ time and
 $O((d_w -d_u)+ N)= O(D+n)$ space.

 As discussed  previously,  if $Mid(u,w)$ only requires $O(n+D)$ space, then  {\tt Path}$(u,w)$ only
 requires $O(n+D)$ space, so we have completed the space analysis.

 It remains to analyze running  time. Set
 $$Area(u,w) = (N-1) (d_w-d_u)
 $$
 to be the ``area'' of $G(u,w)$.
 Recall that line 3 of  {\tt Path}$(u,w)$implies  that $d_u \not = d_w$ when $Mid(u,w)$ is called.
 Therefore  $N \ge 1$ and the running time of
 $Mid(u,w)$ is
 $$O((d_w -d_u)N) = O(Area(u,w)).$$

We now analyze the {running time} of {\tt Path}$(u_0,w_0)$.
First consider the recursive calls
 when lines 1-4  occur, i.e., the recursion terminates.
 The {\em total} work performed by such calls is the total number of
 edges outputted.  Since an edge is outputted only once and
 the total path contains $D$ edges, the total work performed is $O(D).$

Next consider the calls when line 5-7  occur.  Since each such call returns a vertex $v$ on the path,
 there are only $D-1$ such calls so lines 6 and 7 are only called $O(D)$ times and their total work, with
 the exception of the call to $Mid(u,v)$, is $O(D).$

Finally consider the work performed by the $Mid(u,w)$ calls.  Partition the calls into levels.
\begin{itemize}
\item  Level 1 is the original call $Mid(u_0,w_0)$.
\item Level 2 contains the recursive calls directly
made by the level-1 call.
\item In general,  level $i$ contains the recursive calls  directly
made by the level-$(i-1)$ calls.
\end{itemize}

Note that if  $Mid(u,w)$ is a level $i$ call
with $u=(d_u,i_u)$ and $w=(d_w,i_w)$) then
$$\frac D {2^i} \le d_w - d_u <  \frac D {2^i} + 1.$$
Furthermore,  by induction,
if  $Mid(u,w)$ and $Mid(u',w')$ are two different  level $i$ calls, then horizontal ranges
$[d_u,d_w]$ and $[d_{u'},d_{w'}]$ are
{\em disjoint} except for possibly $d_w = d_{u'}$ or $d_u = d_{w'}$.

Fix $i.$  Let $(u_j,w_j)$ $j=1,\ldots t$ be the calls at level $i$.  The facts that each grid $G(u_j,w_j)$
has height $\le \frac D {2^i} + 1$ and that the horizontal ranges of the grids are disjoint implies
$$\sum_{j=1}^t Area(u_j,w_j) \le n \left(\frac D {2^i}+ 1\right).$$

Thus the total
of all level-$i$ calls is $O\left( n \left(\frac D {2^i} + 1\right)\right)$. Summing over the $\lceil \log D \rceil$
levels we get  that the total work performed by all of the $Mid(u,w)$ calls on line 6 is
$$O\left( \sum_i   n \left(\frac D {2^i} + 1\right)\right) = O(nD).$$
Thus,  the total work performed by {\tt Path}$(u_0,w_0)$ is $O(nD)$ and we are finished.

\section{Further Applications}
\label{sec:further}

We just saw how, in $\Theta(nD)$ time and $\Theta(n+D)$ space,
to solve the construction problem for any  DP in form  (\ref{eq:DP_Def})
that satisfies the Monge property (\ref{eq:def_Quad}).
$\Theta(nD)$ time was known previously;  the $\Theta(n+D)$ space bound, is
the new improvement.
There are many other DP problems
besides the binary LLHC that satisfy (\ref{eq:def_Quad}) and whose space
can thus be improved.
We illustrate with three examples.

\medskip

\par\noindent\underline{The $r$-ary LLHC problem:}\\
We have discussed the binary LLHC problem in which $|\Sigma|=2.$  The general $r$-ary alphabet case
with $N$ probabilities  is
still modeled by a DP in form  (\ref{eq:DP_Def}) but with $n = \frac {N-1} {r-1} + 1$.
 The only difference is that (\ref{eqn:c})
is replaced by
\begin{equation}
\label{eqn:c_r}
c_{i,j}^{(d)} =
\left\{
\begin{array}{ll}
S_{ri-j} & \textrm{if } \max\{0,ri-N\} \le j < i \\
\infty & \textrm{otherwise.}
\end{array}
\right.
\end{equation}
A full derivation of this DP is given in Appendix \ref{app:LLHCDP}.
The proof that the
$c_{i,j}^{(d)}$ satisfy the Monge property (\ref{eq:def_Quad}) is similar to  the proof of
Lemma \ref{lem:c.qi}.
Thus,  we can construct a solution to the $r$-ary LLHC problem in $\Theta(ND)$ time and $\Theta(N)$ space as well.

\medskip

\par\noindent\underline{$D$ medians on a line:}\\
We are given $n-1$ customers located on the positive real  line;  customer $i$ is at location $v_i.$
Without loss of generality, assume $v_1 < v_2 < \cdots < v_{n-1}.$  There are $D \le n$ service centers located
on the line and a customer is serviced by the closest service center to its left (thus we always assume
a service center at $v_0=0$).  Each customer has a service request $w_i >0.$  The cost of servicing
customer $i$ is $w_i$ times the distance to its service center.  In \cite{Wo00}, motivated by
the application of optimally placing web proxies on a linear topology network,
 Woeginger showed that this
problem could be modeled by a DP in form in form  (\ref{eq:DP_Def}) where
$$c^{(d)}_{i,j} = \sum_{l=j+1}^i w_l(v_l-v_{j+1})$$
and proved that these $c^{(d)}_{i,j}$ satisfy Monge property
(\ref{eq:def_Quad}).  He then used the SMAWK algorithm to construct a solution in $O(nD)$ time and $O(nD)$ space.
Using the technique we just described,  this can be reduced to $O(nD)$ time and $O(n)$ space.

We also mention that there is an undirected variant of this problem in which a node is serviced by its {\em closest}
service center looking both left and right.  There are many algorithms in the literature that
(explicitly or implicitly) use concavity to construct solutions for this problem
 in $O(nD)$ time using $O(n)$ space,
e.g., \cite{HaTa91,AuPa98+,FGZ06}. \cite{FGZ06} does this by using a DP formulation that is in the
DP form (\ref{eq:DP_Def})
and satisfies the Monge property (\ref{eq:def_Quad}) so
the technique in this paper can reduce the space  for this problem down to
$O(n)$ as well.

\medskip

\par\noindent\underline{Wireless Paging:}\\
The third application comes from wireless mobile paging. A user can be in one of $N$
different cells. We are given
a probability distribution in which $p_i$ denotes the probability that a
user will be in cell $i$ and want to minimize the bandwidth needed to send paging requests to
identify the cell in which the user resides.
This problem was originally conjectured to be NP-complete, but \cite{KGW04+}  developed a
DP algorithm for it.
The input of the problem is the $n$ probabilities $p_1 \ge p_2 \ge \cdots \ge p_n$ and an integer
$D\le n$ (corresponding to the number of paging rounds used).
The DP developed by \cite{KGW04+} is exactly in our DP  form (\ref{eq:DP_Def})
with
\begin{equation}
\label{eq:wmp.cd}
c^{(d)}_{i,j} = \left\{
\begin{array}{ll}
i \left(\sum_{\ell=j+1}^i p_{\ell}\right) & \textrm{if } d-1 \le j < i \\
\infty & \textrm{otherwise.}
\end{array}
\right.
\end{equation}
The goal is to compute $H(D,n)$, which will be the minimum expected bandwidth needed.
Solving the construction version of this DP  permits constructing the actual paging protocol that yields this minimum
bandwidth.

\cite{KGW04+} used the naive algorithm to solve the DP in $\Theta(n^2D)$ time and $\Theta(nD)$ space.
\cite{BGZ06} proved that the $c^{(d)}_{i,j}$ defined by (\ref{eq:wmp.cd}) satisfy the
the Monge property (\ref{eq:def_Quad}) and thus reduced the time to $\Theta(nD)$, but still
required  $\Theta(nD)$ space.
The algorithm in this paper permits improving the space complexity of constructing the protocol
down to $\Theta(n)$.

\section{Conclusion}
The standard approach to  solving the Length-Limited Huffman Coding (LLHC) problem is via the
special purpose Package-Merge algorithm of  Hirschberg and Larmore \cite{LaHi90} which runs in $O(nD)$
time and $O(n)$  space, where $n$ is the number of codewords and $D$ is the length-limit on the code.

In this note we point out that this problem can be solved in the same time and space using a straightforward
Dynamic Programming formulation.  We started by noting that it was known that
the LLHC problem could  be modeled using a DP in the form
{\small
\begin{equation}
\label{eq:H_final}
H(d,i) =
\left\{
\begin{array}{ll}
0 & \mbox{if $d = 0$, $i = 0$} \\
\infty & \mbox{if $d = 0$, $0 < i <n$} \\
\displaystyle{\min_{0 \le j \le i} \left(H(d-1,j) + c^{(d)}_{i,j}\right)} & \mbox{if $d>0$, $0 < i < n$}
\end{array}
\right.
\end{equation}
}

where $H(d,n)$ will denote the minimum cost of a code with longest word at most $d$ and the
$c^{(d)}_{i,j}$ are easily calculable constants. This implies an $O(n^2D)$ time $O(nD)$ space algorithm.
We then note that, using standard DP speedup techniques, e.g., the SMAWK algorithm,  the time could be
reduced down to $O(nD)$.  The main contribution of this paper is to note that, once the problem is
expressed in this formulation,   the space can be reduced down to $O(n)$ while maintaining the
time at $O(nD)$.  The space reduction developed for this problem was also shown to apply to other problems
in the literature that previously had been thought to require $\Theta(nD)$ space.

We conclude by noting that if we're only interested in solving the standard Huffman coding problem and not the LLHC one then DP (\ref{eq:H_final}) with
$c^{(d)}_{i,j}$ defined by (\ref{eqn:c_r}) collapses down to
\begin{equation}
H(i) = \min_{\max\{0,ri-N\} \le j < i} H(j) + S_{ri-j}.
\end{equation}
where $H(i)$ denotes the minimum cost of a ``valid sequence'' ending in $i$.
$H\left(\frac  {N-1} {r-1}\right)$ will be the cost
of an optimal complete sequence and solving the construction
problem for this DP will give this optimal sequence.  We can construct the code
from this optimal sequence in $O(N)$ time.

There is a subtle point here which should be mentioned.  The matrix $M$ defined by
$$
M_{i,j} = \left\{
\begin{array}{ll}
H(j) + S_{ri-j}  & \mbox{ if }  \max\{0,ri-N\} \le j < i\\
\infty  & \mbox{ otherwise}
\end{array}
\right.
$$
is Monge (the proof is similar to that of Lemma \ref{lem:c.qi}).  We can {\em not} use the SMAWK algorithm to
find its row minima and solve the problem,  though.  The reason is that,
as stated in Lemma \ref{lem:SMAWK},  the SMAWK algorithm
requires  being able to calculate any arbitrary requested entry $M_{i,j}$ in $O(1)$ time.  In
our current DP, though,  the $M_{i,j}$ are dependent upon the values $H(j)$
which are the  row-minima of other rows
in the same matrix!  Thus,  we have no way of calculating
$M_{i,j}$ in $O(1)$ time when required
and  the SMAWK algorithm can not be applied.  This is the reason why
Larmore and Przytycka  \cite{larmore:1163}   needed to use the more sophisticated
CLWS algorithm of  \cite{Wi88} to solve the binary ($r=2$) version of this
problem.  Other algorithms for more generalized versions of the
CLWS  have since appeared, e.g., \cite{LaSc91}, that could also be used to
solve this problem in $O(n)$ time, but they are also quite complicated.
To summarize,  by transforming $r$-ary
Huffman coding into a DP and using sophisticated tools such as
\cite{Wi88} or \cite{LaSc91} we can solve the problem in
$O(n)$ time.  This is not of practical interest, though,  since the simple, greedy, Huffman encoding algorithm is just as fast. Where the DP formulation helps is in the LLHC problem, exactly where the greedy procedure fails.  In that case we have the added practical benefit of being able to use the
simple SMAWK algorithm rather than the more complicated \cite{Wi88} or \cite{LaSc91}.

\bibliographystyle{IEEEtran}


TITLE: an example on Lie bialgebra
QUESTION [1 upvotes]: An example: $\mathfrak{sl}(2,\mathbb C)$
Let us consider the Lia algebra $\mathfrak g=\mathfrak{sl}(2,\mathbb C)$, with basis
  $$H=\begin{pmatrix}1&0\\0&-1\end{pmatrix},
\ X=\begin{pmatrix}0&1\\0&0\end{pmatrix},
\ Y=\begin{pmatrix}0&0\\1&0\end{pmatrix}$$
  and commutation relations,
  $$[H,X]=2X,\ [H,Y]=-2Y,\ [X,Y]=H.$$
  The dual vector space $\mathfrak g^*$ has the dual basis $H^*, X^*, Y^*$, where, by definition, $\langle H^*,H\rangle=1$, $\langle X^*,X\rangle=1$, $\langle Y^*,Y\rangle=1$, and all other duality brackets are $0$. We shall consider the following commutation relations on $\mathfrak g^*$,
  $$[H^*,X^*]_{\mathfrak g^*}=\frac14X^*,
\ [H^*,Y^*]_{\mathfrak g^*}=\frac14Y^*,
\ [X^*,Y^*]_{\mathfrak g^*}=0.$$

My question is: how to get the final equation?
Thank you very much!

REPLY [2 votes]: This seems to be $1.1$ here, in the book by Basil Grammaticos and K.M. Tamizhmani. The commutation relations are chosen such that $\mathfrak{g}\bowtie \mathfrak{g}^{*}$ is a Lie algebra having both $\mathfrak{g}$ and $\mathfrak{g}^*$ as subalgebras. The Lie bracket $[,]_{\mathfrak{g}^*}$ is recovered on page $117$ in an explicit way. In general, Lie bialgebra structures on $\mathfrak{sl}_2(\Bbb{C})$ are parametrised by the $3$-dimensional vector space $\Lambda^2(\mathfrak{sl}_2(\Bbb{C}))$.