% C = LOGSUMEXP(A, B) % % Computes LOG(EXP(A)+EXP(B)), with attention to numerical precision. % % A naive computation would overflow if A or B were very positive, or if % both A and B were very negative, yielding +Inf or -Inf respectively. % More subtly, if A or B were moderately negative, the naive computation % would underflow (lose relative precision). Moderately negative means % that exp(A) or exp(B) is small but does not itself cause underflow. % % Example of positive overflow: % >> log(exp(-1e5)+exp(1e5)) % ans = Inf % >> logsumexp(-1e5,1e5) % ans = 100000 % % Example of negative overflow: % >> log(exp(-1e5)+exp(-1e5)) % ans = -Inf % >> logsumexp(-1e5,1e5) % ans = -9.9999e+04 % % Example of underflow: % >> log(exp(-1e2)+exp(0)) % ans = 0 % >> logsumexp(-1e2,0) % ans = 3.7201e-44 function c = logsumexp(a, b) % ensure that A is the larger argument mask = (a < b); s = a; a(mask) = b(mask); b(mask) = s(mask); % factor out exp(A), computing A + log(1 + exp(B-A)) c = a + log1p(exp(b-a));