Log-Sum-Exp Trick to Prevent Numerical Underflow

标签：

引自： http://wittawat.com/log-sum_exp_underflow.html

Multiplying a series of terms p₁p₂?pn where 0≤p_i≤1 can easily result in a numerical underflow. This is especially the case when dealing with probabilistic models which involve such multiplication all the time. One solution commonly used is to work with the log of probabilities instead:

logp₁?p_n=∑_ilogp_i,

which turns multiplication into a summation. Since a summation does not decrease the magnitude of the result, the underflow problem can be avoided.

Unfortunately, this little trick alone does not work when we have a sum of products of probabilities as in the case of, for example, the backward algorithm used in HMM:

β_k(z_k)=∑_zk+1 β_k+1(z_k+1)p(x_k+1 | z_k+1)p(z_k+1 | z_k).

The exact meaning of each term is irrelevant. What is important here is that as we successively apply this recurrence relation (i.e., as k gets lower), β_k(z_k) will underflow. The solution here is the same as before. That is, we work with logβ_k(z_k) instead. So, by taking the log on both sides as well as putting explog in the sum, we have:

logβ_k(z_k)=log∑_zk+1 exp logβ_k+1(z_k+1)p(x_k+1 | z_k+1)p(z_k+1 | z_k)

               →log∑_iexpa_i.

The reason to apply explog on the RHS is to make the log version of β_k+1 so the recurrence relation can be used. We have defined a_i which absorbs the log term on the RHS. Note that typically ai?0. So we still have the underflow when calculating expa_i. The remedy this time is to define

b=max_ia_i,  a_i = logβ_k+1(z_k+1)p(x_k+1 | z_k+1)p(z_k+1 | z_k)

and calculate logexp(b)∑_iexp(a_i−b) instead.

log∑_iexpa_i=logexp(b)∑_iexp(a_i−b)=b+log∑_iexp(a_i−b)

Since b<0 (因为 a_i = logβ_k+1(z_k+1)p(x_k+1 | z_k+1)p(z_k+1 | z_k)), a_i−b is actually positive. So, underflow does not occur when computing exp(a_i−b). We have avoided the underflow problem by using the RHS of the last line.

对比下面这段代码：

def logsumexp(A):
    k = max(A)
    return log(sum( exp(i-k) for i in A ))+k

def calc_feat(x, i, l, r):
    return { ("T", l, r): 1, ("E", r, x[i]): 1 }
def calc_e(x, i, l, r, w, e_prob):
    if (i, l, r) not in e_prob:
        e_prob[i,l,r] = dot(calc_feat(x, i, l, r), w)
    return e_prob[i,l,r]

def calc_f(x, i, l, w, e, f):   #### forward
    if (i, l) not in f:
        if i == 0:
            f[i,0] = 0
        else:
            prev_states = (range(1, len(tagids)) if i != 1 else [0])
            f[i,l] = logsumexp([
                calc_f(x, i-1, k, w, e, f) + calc_e(x, i, k, l, w, e)
                    for k in prev_states])
    return f[i,l]   # 重复利用log值

Soft-Max Function

Just a side note, log∑_iexpa_i （分母）is called a soft max function because

log∑_iexpa_i≈max_ia_i.

The approximation error tends to get larger as the magnitude of a_i gets smaller, however.

Source: (ML 14.10) Underflow and the log-sum-exp trick by mathematicalmonk

Log-Sum-Exp Trick to Prevent Numerical Underflow

标签：

原文地址：http://www.cnblogs.com/Dminer/p/4340716.html