Maximum of the function $x\mapsto -N\log(x)+\sum\log(x-i)$

Question

Fix $n\in\mathbb{N}$ and $N> n$. I am trying to find the maximum of the following function: $$ f:[n,\infty)\rightarrow\mathbb{R},\quad x\mapsto -N\log(x)+\sum_{i=0}^{n-1}\log\left(x-i\right) $$ This function arose in a maximum-likelihood problem I am studying, which is why I expect that there should be a maximum. Since $f$ is differentiable, we can find candidate solutions with the first-order-condition: $$ \frac{N}{x}=\sum_{i=0}^{n-1}\frac{1}{x-i} $$ The right-hand-side can be rewritten as: $$ \sum_{i=0}^{n-1}\frac{1}{x-i} =\sum_{j=x-n+1}^x\frac{1}{j} =\sum_{j=0}^x\frac{1}{j}-\sum_{j=0}^{x-n}\frac{1}{j} =H_{x}-H_{x-n} $$ Above, $H_x$ denothes the $x$-th harmonic number. Here I get stuck. Is there a way to find a closed-form solution for the maximum of $f$? Alternative approaches are also welcome!

Answer 1

$$\color{red}{\frac{N}{x}=\sum_{i=0}^{n-1}\frac{1}{x-i} }\qquad \implies N-n=\sum_{i=1}^{n-1}\frac{i}{x-i} $$ In this form, it looks very similar to the so-called Underwood equation $$\sum_{i=1}^n \frac{\alpha_i\, z_i}{ \theta-\alpha_i}=1-q$$

Our most recent work was published in $2014$ in this paper (you can also find it here) where we proposed rapid and robust solution methods using convex transformations.

Which is a shame is that we do not focus on the largest root (the one you want) but we wan adapt what was done earlier.

In order to remove the left asymptote, rewrite the original equation and consider that you look for the zero of function $$f(x)=(N-1)\Big(x-(n-1)\Big)-x \Bigg(1+\sum_{i=1}^{n-2}\frac{x-(n-1)}{x-i} \Bigg)$$ which is close to linearity when $x$ is large (this is good for any root-finding method).

The only (small) problem is that, before the root, this function goes to a minimum. So, by inspection, we need to find a value $x_*$ such that $f(x_*)>0$.

Expanding $f(x)$ as a series for infinitely large values of $x$,

$$f(x)=(N-n)\,x-\frac{(n-1) (2 N-n)}{2} +\frac{n(n-1) (n-2) }{6 x}+O\left(\frac{1}{x^2}\right)$$

$$\color{blue}{x_*=\frac{(n-1) (2 N-n)}{2 (N-n)}}$$ would be a good initial estimate. Moreover, since $f''(x) >0$, then, by Darboux theorem, we know that it is an overestimate of the solution which will be reached without any overshoot.

For a test, trying for $n=123$ and $M=234$, Newton iterates will be $$\left( \begin{array}{cc} k & x_k & \text{Objective function}\\ 0 & 189.59459 & -634.985\\ 1 & 162.06335 & -631.477\\ 2 & 159.31150 & -631.433\\ 3 & 159.25471 & -631.433\\ 4 & 159.25468 & -631.433\\ \end{array} \right)$$

What we can see, at least based on a small number of simulations, is that, when $M$ and $n$ are large, the area around the maximum of the objective function is flat. If not requiring for too much accuracy, the first iterate of Newton method (there is a simple closed form for it) could be sufficient.

For the worked example, for $150 \leq x \leq 170$, the objective function varies by less than $0.1$%.

Edit

We could do a bit better (still using the series expansion) and use for $x_*$ the largest solution of the quadratic equation $$6 (M-n)\, x^2-3(n-1) (2 M-n)\,x+n(n-1) (n-2)=0$$ Fot the worked case, it would give $x_*=173.919$ but this will only save one iteration.

Update

Considering the initial estimate proposed above $$x_*=\frac{(n-1) (2 N-n)}{2 (N-n)}$$ consider $N=kn$. For large values of $n$ $$\frac{x_*}n=\frac{2 k-1}{2 (k-1)}-\frac{2 k-1}{2 (k-1) n}+O\left(\frac{1}{n^3}\right) \sim \text{Constant}$$ This seems also by true for the solution.

Consider $k=2$ $$\left( \begin{array}{ccc} n & x_n & \frac{x_n} n \\ 10 & 11.6798 & 1.16798 \\ 20 & 24.2444 & 1.21222 \\ 30 & 36.7990 & 1.22663 \\ 40 & 49.3513 & 1.23378 \\ 50 & 61.9026 & 1.23805 \\ 60 & 74.4535 & 1.24089 \\ 70 & 87.0041 & 1.24292 \\ 80 & 99.5546 & 1.24443 \\ 90 & 112.105 & 1.24561 \\ 100 & 124.655 & 1.24655 \\ \end{array} \right)$$

In other words, knowing a solution for a "small" $n$ allows to generate a much better estimate for a much larger $n$.

For example, using the value for $n=100$ and using the approximation for $n=1234$ gives an estimate of $124.655\times \frac {1234}{100}=1538.24$ (instead of $1849.5$) while the exact solution is $1547.83$.

In terms of the polygamma function, using $x_0=x_*$, the first iterate of Newton method is

$$x_1=x_*-\frac{f(x_*)}{f'(x_*)}$$ with $$f'(x)=(N-2)+(2x-n+1) \,\,(\psi ^{(0)}(n-x-1)-\psi ^{(0)}(1-x))-$$ $$x (x-n+1)\,\, (\psi ^{(1)}(n-x-1)-\psi ^{(1)}(1-x))$$

but $x_1$ is not as good as @Jason Zimba's simple and elegant approximation.

Answer 2

I'm not sure if such a heuristic result will be helpful, but for large $n$, I believe the maximum of $f$ occurs fairly close to $\hat{x} = n\left(1+\frac{n}{N}\,W(-\frac{N}{n}\,e^{-\frac{N}{n}})\right)^{-1}$, where $W$ is the Lambert W-function. Here are a few values for comparison:

$$\left.\matrix{n & N & \hat{x} & {\rm Numerical}\cr 100 & 200 & 125.5 & 124.655 \cr 300 & 305 & 9200.55 & 9170.19 \cr 300 & 700 & 345.593 & 344.87 \cr 1000 & 3000 & 1063.29 & 1062.68 }\right.$$

The "Numerical" values are the outputs of Mathematica's NMaximize command applied to $f(x) = -N\log{x}+\sum_{i=0}^{n-1}{\log{(x-i)}}$. That is, ${\tt x /. NMaximize[\{f[x, n, bigN], x > n\}, x][[2, 1]]}$. Note that the NMaximize command takes a long time to execute on my laptop: about 2 minutes for the case $n=1000$, $N=3000$.

The expression for $\hat{x}$ results from crudely replacing the sum $\sum_{i=0}^{n-1}{\frac{1}{x-i}}$ with the integral $\int_0^{n}{\frac{di}{x-i}}$ in the condition for vanishing derivative. This results in the equation $\frac{N}{x} = \log{\frac{x}{x-n}}$. Defining $z = \frac{1}{x}$, we have $e^{Nz}(1-nz) = 1$ which shows that $z$ is close to $\frac{1}{n}$ (so that the product of a large and small factor is unity). In general, unraveling such transcendental equations leads to the Lambert W-function. (Mathematica calls it the ${\tt ProductLog}$ function.)

The expression for $\hat{x}$ is a little messy. When $\frac{N}{n}\gg 1$, an approximation is $\hat{x} \approx n(1+e^{-\frac{N}{n}})$. For example,

$$\left.\matrix{n & N & n(1+e^{-\frac{N}{n}}) & \hat{x} & {\rm Numerical}\cr 40 & 80 & 45.4 & 50.2 & 49.3512 \cr 100 & 500 & 100.674 & 100.703 & 100.13}\right.$$

Update: For completeness, I'll show how to express the solution of $e^{Nz}(1-nz)=1$ in terms of the Lambert W-function. Define $u=1-nz$, then the equation says $u\,e^{\frac{N}{n}(1-u)}=1$. Multiplying both sides by $e^{-\frac{N}{n}}$, we have $u\,e^{-\frac{N}{n}u}=e^{-\frac{N}{n}}$. Multiplying both sides by $-\frac{N}{n}$, we have $(-\frac{N}{n}u)\,e^{(-\frac{N}{n}u)}=-\frac{N}{n}\,e^{-\frac{N}{n}}$. Now by definition, the Lambert W-function $W(y)$ has defining property $W\,e^W = y$. So from $(-\frac{N}{n}u)\,e^{(-\frac{N}{n}u)}=-\frac{N}{n}\,e^{-\frac{N}{n}}$, we have $-\frac{N}{n}u = W(-\frac{N}{n}\,e^{-\frac{N}{n}})$. Re-expressing this in terms of $z$, $-\frac{N}{n}(1-nz) = W(-\frac{N}{n}\,e^{-\frac{N}{n}})$. This equation is easily solved for $z$, and hence for $x$, leading to the expression for $\hat{x}$ given above.