5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Stupid question: your objective is to normalise your $x_i$ Or do you have specific constraint on the problem formulation? $\endgroup$ – gdupont Aug 10 '16 at 20:31
  • 2
    $\begingroup$ What do powers of unnormalized probabilities represent? (/how does this particular need arise? It seems rather an odd problem to be honest and I wonder if there might be a different approach that solves it better) Do you have positivity constraints on $\lambda$? $\endgroup$ – Glen_b -Reinstate Monica Aug 11 '16 at 0:53
  • $\begingroup$ The exponent of a probability can be interpreted as the number of independent trials. Without going into details, if $x$ are likelihoods of data, then $\lambda$ shapes the corresponding distribution. And yes, $\lambda$ should be greater than zero under normal circumstances. $\endgroup$ – fungs Aug 12 '16 at 12:20
3
$\begingroup$

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.

Histogram of y and plot of f with the search points marked

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`.
#------------------------------------------------------------------------------#
$\endgroup$
  • $\begingroup$ Thanks for your very comprehensive answer, I was thinking of the Newton-Raphson method but I wasn't sure I was missing something. Thanks for the starting point and log considerations! $\endgroup$ – fungs Aug 12 '16 at 12:23
0
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ Note that the $x_i$ do not some to one, it is a vector of probabilities but not a normalized probability vector. I corrected the posting so that there cannot be multiple ones in the vector (for there is no solution i this case). $\endgroup$ – fungs Aug 10 '16 at 17:55
  • $\begingroup$ @fungs edited, I proposed another solution $\endgroup$ – RUser4512 Aug 10 '16 at 22:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.