What is the best way to determine an exponent $\lambda$ of a vector of probabilies $x$ so that the sum of the power-transformed vector equals one?
$\sum\limits_{i=1}^{|x|} x_i^\lambda = 1$
$x$ is a large sample and contains values between zero and one. Let's assume that there is at least one value different from zero and at most one value at one. I'm looking for a fast computational method to find $\lambda$. I'm also interested in approximate, numerical approaches.
There is no closed-form solution. Computation must be numerical.
Let $n$ be the length of the vector, set $y_i = \log(x_i)$, and define
$$f(\lambda) = \log\left(\sum_{i=1}^n e^{\lambda y_i}\right)$$
for $\lambda \ge 0$. The problem is to find the (unique positive) zero of $f$. Using the logarithm helps linearize the function near its zero.
The task, then, is to find a value of $\lambda$ close to that zero and then polish it, perhaps using a few Newton-Raphson iterations. These are easy to calculate because computing the derivative of $f$ requires almost no more effort than finding $f$ itself, where the work consists of computing the vector $z=(e^{\lambda y_i})$, for
$$\frac{d}{d\lambda}f(\lambda) = \frac{\sum_{i=1}^n y_i e^{\lambda y_i}}{\sum_{i=1}^n e^{\lambda y_i}} = \frac{y \cdot z}{1 \cdot z}$$
amounts to taking two quick dot products.
A general-purpose approach begins by bracketing the root between $0$ and an upper bound. One upper bound is found by replacing each $y_i$ by the largest of the $y_i$, giving
$$\lambda^{*} = \frac{-\log(n)}{\max(y)}.\tag{1}$$
A bisection search (halve the value of $\lambda^{*}$ until the value of $f$ exceeds $1$) does well even with challenging distributions of the $x_i$. The entire procedure typically requires five to 12 evaluations of $f$ for $n \le 10^6$, requiring much less than a second to complete. Here is an example.
A third of the values of $x$ are close to $1$ while two-thirds are tiny. This causes the bound $(1)$ to grossly overestimate the root. That initial estimate is shown on the right panel, where $\exp(f)$ is plotted on log-log axes: the dotted red line at the right marks its value. Consequently, six bisection steps are needed. They are shown as the regular progression of dotted lines with colors moving through the spectrum from red to yellow, green, and blue. The leftmost line marks the value of $\lambda$ where $f$ is found to exceed $0$. At this point, five Newton-Raphson iterations rapidly converge on the solution, which is computed to double-precision tolerances (in $\lambda$ and $f(\lambda)$). The final value is plotted in purple, exactly where $\exp(f)=1$.
The R code that produced these plots includes a general-purpose solution.
#
# Generate sample data.
#
n <- 1e6
set.seed(17)
a <- c(-10, -5, -0.5)
b <- c(-5, -3, 0)
x <- exp(runif(n, a, b))
#
# Find lambda.
#
system.time(
{
y <- log(x)
n <- length(x)
#
# Initial bisection search.
#
f <- function(lambda, y) {
z <- exp(lambda * y)
list(f=log(sum(z)), df=crossprod(y, z) / sum(z))
}
lambda <- -log(n) / max(y)
z <- c(lambda, rep(NA, iter.max))
j <- 1
while(f(lambda, y)$f < 0) { #$
lambda <- lambda/2
j <- j+1; z[j] <- lambda
}
#
# Final NR search.
#
tol <- 1e-12
tol2 <- 1e-14
iter.max <- 25
delta <- NA
converged <- -1
for (i in j:iter.max) {
#
# Evaluate the function and compare to zero.
#
converged <- 0
delta.old <- delta
lambda.old <- lambda
w <- f(lambda, y)
if (abs(w$f) < tol2) break
#
# Take a step.
#
converged <- 1
delta <- w$f / w$df # NR
if (delta > lambda) delta <- delta / 2 # Bisection
lambda <- lambda - delta
#
# Store the step and check for progress.
#
z[i+1] <- lambda
if (abs(delta) < tol * abs(lambda)) break
converged <- -1
}
z <- z[1:max(which(!is.na(z)))]
}
)
if (converged < 0) stop("Failed to converge.")
#
# Plot the search.
#
par(mfrow=c(1,2))
hist(y, breaks=100)
x.0 <- seq(min(z, na.rm=TRUE)/1.25, max(z, na.rm=TRUE)*1.25, length.out=101)
z.0 <- sapply(x.0, function(lambda) exp(f(lambda, y)$f))
plot(x.0, z.0, type="l", log="xy", xlab="lambda", ylab="",
main=paste("n =", n),
sub=paste(i, "evaluations"))
abline(h=1, lwd=2, col="Gray")
abline(v=z, lty=3, col=hsv(seq(0, 5/6, length.out=length(z)), 0.8, 0.8))
#
# Display some useful statistics, such as the last NR steps.
#
print(c(d.f=sum(exp(lambda * y)) - 1, d.loglambda=delta/lambda))
#
# The entire search is in `z`. The final value is stored in `lambda`.
#------------------------------------------------------------------------------#
First note that: $$\lambda \rightarrow f(\lambda)=\sum\limits_{i=1}^{|x|} x_i^\lambda = 1$$ is strictly decreasing with $\lambda$ (provided all the $x_i's$ are in $[0,1]$).
Then, $f(0)=|x|$ and $f\rightarrow |\{i|x_i=1\}|$ as $\lambda\rightarrow\infty$.
However, as the function is strictly decreasing, there is only one solution in $[0,\infty)$ and it exists if all the $x_i$'s are smaller than 1.
If your $x$ is a vector of probabiliies, $\lambda=1$ is a solution (and the only one).
Edit. As the vector of probabilities is not normalized, the above solution does not work.
using the inequality :
$$||x||_\infty^\lambda|x|\geq f(\lambda)$$
Looking for $\lambda$ such that :
$$1\geq||x||_\infty^\lambda|x| $$
you are sure that the solution $\lambda^*$ lies in $[0,-\frac{\log(|x|)}{\log(||x||_\infty)}]$ and this is a good starting point for a dichotmic search.
