Visualizing Probability Update Schemes

Unlike the last "visualizing" post, I don't have a real clear point to this one. This is more like "let's look at some pretty pictures".

I do think it helps to get intuition by actually looking at charts & graphs, rather than just always look at the end result score for something like compression.

We're going to look at binary probability estimation schemes.

Binary probability estimation is just filtering the history of bits seen in some way.

Each bit seen is a digital signal of value 0 or 1. You just want some kind of smoothing of that signal. Boom, that's your probability estimate, P(1). All this nonsense about Poisson and Bernoulli processes and blah blah, what a load of crap. It's just filtering.

For example, the "Laplace" estimator

P(1) = (n1 + c)/(n0 + n1 + 2c)

That's just the average of the entire signal seen so far. Add up all the bits, divide by number. (countless papers have been written about what c should be (c= 1? c = 1/2?), why not 1/4? or 3/4? Of course there's no a-priori reason for any particular value of c and in the real world it should just be tweaked to maximize results.)

That's the least history-forgetting estimator of all, it weights everything in the past equally (we never want an estimator where older stuff is more important). In the other direction you have lots of options - FIR and IIR filters. You could of course take the last N bits and average them (FIR filter), or take the last N and average with a weight that smoothly goes from 0 to 1 across the window (sigmoidal or just linear). You can of course take an IIR average, that's

average <- (average*N + last)/(N+1)

Which is just the simplest IIR filter, aka geometric decay, "exponential smoothing" blah blah.

Anyway, that's just to get us thinking in terms of filters. Let's look at some common filters and how they behave.

---

Green bars = probability of a 1 bit after the bit at bottom is processed.

Blue bars = log odds = codelen(0) - codelen(1)

---

Laplace : count n0 & n1 :

After 10 1's and 10 0's, P is back to 50/50 ; we don't like the way this estimator has no recency bias.

---

Geometric IIR with updshift 2 (very fast) :

When a fast geometric gets to the 010101 section is wiggles back and forth wildly. If you look at the codelen difference in that region you can see it's wasting something like 1.5 bits on each coding (because it's always wiggled in just the wrong way, eg. biased towards 0 right before a 1 occurs).

Note that geometric probabilities have an exponential decay shape. (specifically , 0.75^n , where 0.75 is from (1 - 1/4) and the 4 is because shift 2). HOWEVER the codelen difference steps are *linear*.

The geometric probability update (in the direction of MPS increasing probability) is very close to just adding a constant to codelen difference (logit). (this just because P *= lambda , so log(P) += log(lambda)). You can see after the 111 region, the first 0 causes a big negative step in codelen difference, and then once 0 becomes the MPS the further 0 steps are the same constant linear step.

---

Geometric IIR with updshift 5 (fast) :

Shift 5 is still fast by real world standards but looks slow compared to the crazy speed of geo 2. You can see that the lack of an initial adaptation range hurts the ramp up on the 1111 portion. That is, "geo 5" acts like it has 33 symbols in its history; at the beginning it actually has 0, so it has a bogus history of 33 times P = 0.5 which gives it a lot of intertia.

---

FIR 8-tap window flat weight :

Just take the last 8 and average. (I initialize the window to 01010 which is why it has the two-tap stair step in the beginning). In practice you can't like the probabilities get to 0 or 1 completely, you have to clamp at some epsilon, and in fact you need a very large epsilon here because over-estimating P(1) and then observing a 0 bit is very very bad (codelen of 0 goes to infinity fast as P->1). The "codelen difference" chart here uses a P epsilon of 0.01

---

bilateral filter :

It's just filtering, right? Why not? This bilateral filter takes a small local average (4 tap) and weights contribution of those local averages back through time. The weight of each sample is multiplied by e^-dist * e^-value_difference. That is, two terms (hence bilateral), weight goes down as you go back in time, but also based on how far away the sample is from the most recent one in value. ("sample" is the 4-tap local average)

What the bilateral does is that as you get into each region, it quickly forgets the previous region. So as you get into 000 it forgets the 111, and once it gets into 0101 it pretty solidly stabilizes at P = 0.5 ; that is, it's fast in forgetting the past when the past doesn't match the present (fast like geo 2), BUT it's not fast in over-responding to the 0101 wiggle like geo 2.

There are an infinity of different funny impulse responses you could do here. I have no idea if any of them would actually help compression (other than by giving you more parameters to be able to train to your data, which always helps, sort of).

---

mixed :

Linear mix of geo 2 (super fast) and geo 5 (still quite fast). mixing weight starts at 0.5 and is updated with online gradient descent. The mixer here has an unrealistically fast learning rate to exaggerate the effect.

The weight shown is the weight of geo 2, the faster model. You can see that in the 111 and 000 regions the weight of geo 2 shoots way up (because geo 5 is predicting P near 0.5 still), and then in the 0101 region the weight of geo 2 gradually falls off because the slow geo 5 does better in that region.

The mixer immediately does something that none of the other estimators did - when the first 0 bit occurs, it takes a *big* step down, almost down to 50/50. It's an even faster step than geo 2 did on its own. (same thing with the first 1 after the 000 region).

Something quite surprising popped out to me. The probability steps in the 111 and 000 regions wind up linear. Note that both the underlying estimators (geo 2 and geo 5) has exponential decay curving probabilities, but the interaction with the mixing weight cancels that out and we get linear. I'm not sure if that's a coincidence of the particular learning rate, but it definitely illustrates to us that a mixed probability can behave unlike any of its underlying experts!

Visualizing Binary Probability Updates

Some time ago I noted that the standard way we use binary probabilities to model statistics has some odd quirks : 06-12-14 - Some LZMA Notes

I thought I would make some pictures to make that more intuitive.

What I have here is a 4 bit symbol (alphabet of 16). It is coded with 4 binary probabilities in a tree. That is :

first code a bit for sym >= 8 or not (is it in [0-7] or [8-15])
then go to [0-3] vs [4-7] (for example)
then [4-5] vs [6-7]
lastly [6] vs [7]

One way you might implement this is something like :

U32 p0[16];

// sym is given
sym |= 16; // set a place-value marker bit

for 4 times
{
  int ctx = sym >> 4; int bit = (sym>>3)&1;
  arithmetic_code( p0[ctx] , bit );
  sym <<= 1;
}

and note that only 15 p0's are actually used, p0[0] is not accessed; p0[1] is the root probability for [0-7] vs [8-15] , p0[2] is [0-3] vs [4-7] , etc.

The standard binary probability update is exponential decay (the simplest geometric IIR filter) :

if ( bit )
{
    P0 -= P0 * lambda;
}
else
{
    P0 += (1.0 - P0) * lambda;
}

So what actually happens to the symbol probabilities when you do this?

Something a bit strange.

When you update symbol [5] (for example), his probability goes up. But so does the probability of his neighbor, [4]. And also his next neighbors [6 and 7]. And so on up the binary tree.

Now the problem is not the binary decomposition of the coding. Fundamentally binary decomposition and full alphabet coding are equivalent and should be able to get the same results if you form the correct correspondence between them.

(aside : For example, if you use a normal n0/n1 count model for probabilities, AND you initialize the counts such that the parents = the sum of the children, then what you get is : visualizing_binary_probability_updates_n0_n1.html - only the symbols you observe increase in probability. Note this is a binary probability tree, not full alphabet, but it acts exactly the same way.)

The problem (as noted in the previous LZMA post) is that the update rate at the top levels is too large.

Intuitively, when a 5 occurs you should update the binary choice [45] vs [67] , because a 5 is more likely, that pair is more likely. The problem is that [4] did not get more likely, and the probability of the group [45] represents the chance of either occuring. One of those things did get more likely, but one did not. The 4 in the group should act as a drag that keeps the group [45] probability from going up so much. Approximately, it should go up by half the speed of the leaf update.

The larger the group, the greater the drag of symbols that did not get more likely. eg. in [0-7] vs [8-15] when a 5 occurs, all of 0123467 did not get more likely, so 7 symbols act as drag and the rate should be around 1/8 the speed.

(see appendix at end for the exact speeds you should use if you wanted to adjust only one probability)

Perhaps its best to just look at the charts. These are histograms of probabilities (in percent). It starts all even, then symbol 5 occurs 20 X, then symbol 12 occurs 20 X. The probabilities are updated with the binary update scheme. lambda here is 1.0/8 , which is rather fast, used to make the visualization more dramatic.

What you should see : when 5 first starts getting incremented, the probabilities of its neighbors go way up, 4, and 6-7, and the whole 0-7 side. After 5 occurs many times, the probabilities of its neighbors goes down. Then when 12 starts beeing seen, the whole 8-15 side shoots up.

Go scroll down through the charts then we'll chat more.

This is a funny thing. Most people in the data compression community think of binary coding symbols this way as just a way to encode an alphabet using binary arithmetic coding. They aren't thinking about the fact that what they're actually doing is a strange sort of group-probability update.

In fact, in many applications if you take a binary arithmetic coder like this and replace it with an N-ary full alphabet coder, results get *worse*. Worse !? WTF !? It's because this weird group update thing that the binary coder is doing is often actually a good thing.

You can imagine scenarios where that could be the case. In some types of data, when a new symbol X suddenly starts occuring (when it had been rare before), then that means (X-1) and (X+2) may start to be seen as well. We're getting some kind of complicated modeling that novel symbols imply their neighbors novel probability should go up. In some type of data (such as BWT post-MTF) the probabilities act very much in binary tree groups. (see Fenwick for example). In other types of data that is very bit structured (such as ascii text and executable code), when a symbol with some top 3 bits occurs, then other symbols with those top bits are also more likely. That is, many alphabets actually have a natural binary decomposition where symbol groups in the binary tree do have joint probability.

This is one of those weird things that happens all the time in data compression where you think you're just doing an implementation choice ("I'll use binary arithmetic coding instead of full alphabet") but that actually winds up also doing modeling in ways you don't understand.

The visuals :

5 x 1

5 x 2

5 x 3

5 x 4

5 x 5

5 x 6

5 x 7

5 x 8

5 x 9

5 x 10

5 x 11

5 x 12

5 x 13

5 x 14

5 x 15

5 x 16

5 x 17

5 x 18

5 x 19

5 x 20

12 x 1

12 x 2

12 x 3

12 x 4

12 x 5

12 x 6

12 x 7

12 x 8

12 x 9

12 x 10

12 x 11

12 x 12

12 x 13

12 x 14

12 x 15

12 x 16

12 x 17

12 x 18

12 x 19

12 x 20

---

Appendix :

How to do a binary probability update without changing your neighbor :

Consider the simplest case : 4 symbol alphabet :

[0 1 2 3]

A = P[0 1] vs [2 3]

B = P[1] vs [1]

P(0) = A * B
P(1) = A * (1 - B)

now a 0 occurs

alpha = update rate for A
beta = update rate for B

A' = A + (1-A) * alpha
B' = B + (1-B) * beta

we want P(1)/P(23) to be unchanged by the update

(if that is unchanged, then P(1)/P(2) and P(1)/P(3) is unchanged)

that is, the increase to P(0) should scale down all other P's by the same factor

P(1)/P(23) = A * (1-B) / (1-A)

so require :

A' * (1-B') / (1-A') = A * (1-B) / (1-A)

solve for alpha [...algebra...]

alpha = A * beta / ( 1 - (1 - A) * beta )

that's as opposed to just using the same speed at all levels, alpha = beta.

In the limit of small beta, (slow update), this is just alpha = A * beta.

The update rate at the higher level is decreased by the probability of the updated subsection.

In the special case where all the probabilities start equal, A = 0.5, so this is just alpha = beta / 2 - the update rates should be halved at each level, which was the intuitive rule of thumb that I hand-waved about before.

Whirlwind Tour of Mixing Part 3

On to online gradient descent.

First we need to talk about how our experts are mixed. There are two primary choices :

Linear mixing :

M(x) = Sum{ w_i P_i(x) } / Sum{ w_i }

or = Sum{ w_i P_i(x) } if the weights are normalized

Geometric mixing :

M(x) = Prod{ P_i(x) ^ c_i } / Sum(all symbols y){ Prod{ P_j(y) ^ c_j } }

(the coefficients c_i do not need to be normalized)

P = probability from each expert
M = mixed probability

require

Sum(all symbols y){ M(y) } = 1

Linear mixing is just a weighted linear combination. Geometric mixing is also known as "product of expert" and "logarithmic opinion pooling".

Geometric mixing can be expressed as a linear combination of the logs, if normalization is ignored :

log( M_unnormalized ) = Sum{ c_i * log( P_i ) }

In the large alphabet case, there is no simple logarithmic form of geometric M (with normalization included).

In the special case of a binary alphabet, Geometric mixing has a special simple form (even with normalization) in log space :

M = Prod{ P_i ^ c_i } / ( Prod{ P_i ^ c_i } + Prod{ (1 - P_i) ^ c_i } )

M = 1 / ( 1 + Prod{ ((1 - P_i)/P_i) ^ c_i } )

M = 1 / ( 1 + Prod{ e^c_i * log((1 - P_i)/P_i) } )

M = 1 / ( 1 + e^ Sum { c_i * log((1 - P_i)/P_i) } )

M = sigmoid( Sum { c_i * logit(P_i) } )

which we looked at in the previous post, getting some background on logit and sigmoid.

We looked before about how logit space strongly preserves skewed probabilities. We can see the same thing directly in the geometric mixer.

Obviously if an expert has a P near 0, that multiplies into the product and produces an output near 0. Zero times anything = zero. Less obviously if an expert has a P near 1, the denomenator normalization factor becomes zero in the other terms (that aren't the numerator), so the mix becomes 1, regardless of what the other P's.

So assume we are using these mixtures. We want to optimize (minimize) code length. The general procedure for gradient descent is :

the code length cost of a previous coding event is -log2( M(x) )
evaluated at the actually coded symbol x

(in compression the total output file length is just the sum of all these event costs)

To improve the w_i coefficients over time :

after coding symbol x, see how w_i could have been improved to reduce -log2( M(x) )

take a step of w_i in that direction

the direction is the negative gradient :

delta(w_i) ~= - d/dw_i { -log2( M(x) ) }

delta(w_i) ~= (1/M) * d/dw_i{ M(x) }

(~ means proportional; the actual delta we want is scaled by some time step rate)

so to find the steps we should take we just have to do some derivatives.

I won't work through them in ascii, but they're relatively simple.

(the geometric mixer update for the non-binary alphabet case is complex and I won't cover it any more here. For practical reasons we will assume linear mixing for non-binary alphabets, and either choice is available for binary alphabets. For the non-binary geometric step see Mattern "Mixing Strategies in Data Compression")

When you work through the derivatives what you get is :

Linear mixer :

delta(w_i) ~= ( P_i / M ) - 1


Geometric mixer (binary) :

delta(w_i) ~= (1 - M) * logit(P_i)

Let's now spend a little time getting familiar with both of them.

First the linear mixer :

put r_i = P_i(x) / M(x)

ratio of expert i's probability to the mixed probability, on the symbol we just coded (x)

Note : Sum { w_i * r_i } = 1  (if the weights are normalized)

Your goal is to have M = 1

That is, you gave 100% probability to the symbol that actually occured

The good experts are the ones with high P(x)  (x actually occured)

The experts that beat the current mix have P > M

Note that because of Sum { w_i * r_i } = 1

Whenever there is an expert with r_i > 1 , there must be another with r_i < 1

The linear mixer update is just :

delta(w_i) ~= r_i - 1

experts that were better than the current mix get more weight
experts worse than the current mix get less weight

I think it should be intuitive this is a reasonable update.

Note that delta(w) does not preserve normalization. That is :

Sum{ delta(w_i) } != 0 , and != constant

however

Sum{ w_i * delta(w_i) } = Sum{ w_i * P_i / M } - Sum { w_i } = Sum {w_i} - Sum{w_i} = 0

the deltas *do* sum to zero when weighted.

We'll come back to that later.

If we write the linear update as :

delta(w_i) ~= ( P_i - M ) / M

it's obvious that delta is proportional to (1/M). That is, when we got the overall mixed prediction M right, M is near 1, the step is small. When we got the step wrong, M is near 0, the step we take can be very large. (note that P_i does not have to also be near zero when M is near zero, the w_i for that expert might have been small). In particular, say most of the experts were confused, they predict the wrong thing, so for the actual coded event x, their P(x) is near zero. One guy said no, what about x! His P(x) was higher. But everyone ignores the one guy piping up (his weight is low). The mixed M(x) will be very low, we make a huge codelen (very bad log cost), the step (1/M) will be huge, and it bumps up the one guy who got the P(x) right.

Now playing with the geometric update :

Geometric mixer (binary) :

delta(w_i) ~= (1 - M) * logit(P_i)

this is evaluated for the actually coded symbol x (x = 0 or 1)

let's use the notation that ! means "for the opposite symbol"
eg. if you coded a 0 , then ! means "for a 1"

eg. !P = 1-P and !M = 1-M

(in the previous post I put up a graph of the relationship between L and !L, codelens of opposing bits)

We can rewrite delta(w_i) as :

delta(w_i) ~= !M * ( !L_i - L_i )

recalling logit is the codelen difference

that's

delta(w_i) ~= (mixed probability of the opposite symbol) * (ith model's codelen of opposite - codelen of coded symbol)

I find this form to be more intuitive; let's talk through it with words.

The step size is proportional to !M , the opposite symbol mixed probability. This is a bit like the (1/M) scaling in the linear mixer. When we got the estimate right, the M(x) of the coded symbol will go to 1, so !M goes to zero, our step size goes to zero. That is, if our mixed probability M was right, we take no step. If it ain't broke don't fix it. If we got the M very wrong, we thought a 0 bit would not occur but then it did, M is near zero, !M goes to one.

While the general direction of scaling here is the same as the linear mixer, the scaling is very different. For example in the case where we were grossly wrong, M -> 0 , the linear mixer steps like (1/M) - it can be a huge step. The geometric mixer steps like (1-M) , it hits a finite maximum.

The next term (!L - L) means weights go up for experts with a large opposite codelen. eg. if the codelen of the actually coded symbol was low, that's a good expert, then the codelen of the opposite must be high, so his weight goes up.

Again this is similar to the linear mixer, which took a step ~ P(x) , so experts who were right go up. Again the scaling is different, and in the opposite way. In the linear mixer, if experts get the P right or wrong, the step just scales in [0,1] ; it varies, but not drastically. Conversely in the geometric case because step is proportional to codelen, it goes to infinity as P gets extremely wrong.

Say you're an expert who gets it all wrong and guesses P(x) = 0 (near zero) then x occurs. In the linear case the step is just (-1). In the geometric case, L(x) -> inf , !L -> 0 , your weight takes a negative infinite step. Experts who get it grossly wrong are not gradually diminished, they are immediately crushed. (in practice we clamp weights to some chosen finite range, or clamp step to a maximim)

Repeating again what the difference is :

In linear mixing, if the ensemble makes a big joint mistake, the one guy out of the set who was right can get a big weight boost. If the ensemble does okay, individuals experts that were way off do not get big adjustments.

In geometric mixing, individual experts that make a big mistake take a huge weight penalty. This is only mildly scaled by whether the ensemble was good or not.

Another funny difference :

With geometric mixing, experts that predict a 50/50 probability (equal chance of 0 or 1) do not change w at all. If you had a static expert that just always said "I dunno, 50/50?" , his weight would never ever change. But he would also never contribute to the mixture either, since as noted in the previous post equivocal experts simply have no effect on the sum of logits. (you might think that if you kept seeing 1 bits over and over, P(1) should be near 100% and experts that keep saying "50%" should have their weight go to crap; that does not happen in geometric mixing).

In contrast, with linear mixing, 50/50 experts would either increase or decrease in weight depending on if they were better or worse than the previous net mixture. eg. if M was 0.4 on the symbol that occurred, the weight of an 0.5 P expert would go up, since he's better than the mix and somebody must be worse than him.

Note again that in geometric mixing (for binary alphabets), there is no explicit normalization needed. The w_i's are not really weights, they don't sum to 1. They are not in [0,1] but go from [-inf,inf]. It's not so much "weighting of experts" as "adding of experts" but in delta-codelen space.

Another funny thing about geometric mixing is what the "learning rate" really does.

With the rate explicitly in the steps :

Linear mixer :

w_i += alpha * ( ( P_i / M ) - 1 )

(then w_i normalized)


Geometric mixer (binary) :

c_i += alpha * ( (1 - M) * logit(P_i) )

(c_i not normalized)

In the linear case, alpha really is a learning rate. It controls the speed of how fast a better expert takes weight from a worse expert. They redistribute the finite weight resource because of normalization.

That's not the case in the geometric mixer. In fact, alpha just acts as an overall scale to all the weights. Because all the increments to w over time have been scaled by alpha, we can just factor it out :

no alpha scaling :

c_i += ( (1 - M) * logit(P_i) )

factor out alpha :

M = sigmoid( alpha * Sum { c_i * logit(P_i) } )

The "learning rate" in the geometric mixer really just scales how the sum of logits is stretched back to probability. In fact you can think of it as a normalization factor more than a learning rate.

The geometric mixer has yet another funny property that it doesn't obviously pass through the best predictor. It doesn't even obviously pass through the case of a single expert. (with linear mixing it's obvious that crap predictors weight goes to zero, and a single expert would just pass through at weight=1)

aside you may ignore; quick proof that it does :

say you have only one expert P
say that P is in fact perfectly correct
0 bits occur at rate P and P is constant

initially you have some mixing coefficient 'c' (combining alpha in with c)

The mixed probability is :

M = sigmoid( c * logit(P) )

which is crap.  c is applying a weird scaling we don't want.

c should go to 1, because P is the true probability.

Does it?

delta(c) after one step of symbol 0 is (1-M) * logit(P)
for symbol 1 it's M * logit(1-P)

since P is the true probability, the probabilistic average step is :

average_delta(c) = P * ( (1-M) * logit(P) ) + (1-P) * ( (M) * logit(1-P) )

logit(1-P) = - logit(P)

average_delta(c) = P * logit(P) - M * logit(P) = (P - M) * logit(P)

in terms of l = logit(P)
using P = sigmoid(l)

average_delta(c) = ( sigmoid(l) - sigmoid( c * l ) ) * l

this is a differential equation for c(t)

you could solve it (anyone? anyone?)

but even without solving it's clear that it does go to c = 1 as t -> inf

if c < 1 , average_delta(c) is > 0 , so c goes up
if c > 1 , average_delta(c) is < 0 , so c goes down

in the special case of l small (P near 50/50), sigmoid is near linear,
then the diffeq is easy and c ~= (1 + e^-t)

Note that this was done with alpha set to 1 (included in c). If we leave alpha separate, then c should converge to 1/alpha in order to pass through a single predictor. That's what I mean by alpha not really being a "learning rate" but more of a normalization factor.

So anyway, despite geometric mixing being a bit odd, it works in practice; it also works in theory (the gradient descent can be proved to converge with reasonable bounds and so on). In practice it has the nice property of not needing any divides for normalization (in the binary case). Of course it doesn't require evaluation of logit & sigmoid, which are typically done by table lookup.

---

Bonus section : Soft Bayes

The "Soft Bayes" method is not a gradient descent, but it is very similar, so I will mention it here. Soft Bayes uses a weighted linear combination of the experts, just like linear gradient descent mixing.

Linear mixer :

d_i = ( P_i / M ) - 1

w_i += rate * d_i

Soft Bayes :

w_i *= ( 1 + rate * d_i )

Instead of adding on the increment, Soft Bayes multiples 1 + the increment.

Multiplying by 1 + delta is obviously similar to adding delta and the same intuitive arguments we made before about why this works still apply. Maybe I'll repeat them here :

General principles we want in a mixer update step :

The step should be larger when the mixer produced a bad result

  the worse M(x) predicted the observed symbol x, the larger the step should be
  this comes from scaling like (1/M) or (1-M)

The step should scale up weights of experts that predicted the symbol well

  eg. larger P_i(x) should go up , smaller P_i(x) should go down

  a step proportional to (P-M) does this in a balanced way (some go up, some go down)
  but a step of (P - 0.5) works too
  as does logit(P)

We can see a few properties of Soft Bayes :

Soft Bayes :

w_i *= ( 1 + rate * d_i )

w_i' = w_i + w_i * rate * d_i

w_i += rate * w_i * d_i

This is the same as linear gradient descent step, but each step d_i is scaled by w_i

Soft Bayes is inherently self-normalizing :

Sum{w_i'} = Sum{w_i} + Sum{w_i * rate * d_i}

Sum{ w_i * d_i } = Sum{ ( w_i * P_i / M ) - w_i } = ( M/M - 1 ) * Sum{ w_i } = 0

therefore :

Sum{w_i'} = Sum{w_i} 

the sum of weights is not changed

if they started normalized, they stay normalized

The reason that "soft bayes" is so named can be seen thusly :

Soft Bayes :

w_i *= ( 1 + rate * d_i )

d_i = ( P_i / M ) - 1

w_i *= (1-rate) + rate * ( P_i / M )

at rate = 1 :

w_i *= P_i / M

ignoring normalization this is 

w_i *= P_i

which is "Beta Weighting" , which is aka "naive Bayes".

So Soft Bayes uses "rate" as a lerp factor to blend in a naive Bayes update. When rate is 1, it's "fully Bayes", when rate is lower it keeps some amount of the previous weight lerped in.

Soft Bayes seems to work well in practice. Soft Bayes could have the bad property that if weights get to zero they can never rise, but in practice we always clamp weights to a non-zero range so this doesn't occur.

A quick aside on why Beta weighting is the naive Bayesian solution :

ignoring normalization through this aside
(also ignoring priors and all the proper stuff you should do if you're Bayesian)
(see papers for thorough details)

Mix experts :

M = Sum{ w_i * P_i }

P_i(x) is the probability of x given expert i ; P(x|i)

The w_i should be the probability of expert i, given past data , that is :

P(x|past) = Sum{ P(i|past) * P(x|i) }

P(i|past) is just ~= P(past|i)
  (ignoring prior P(i) and normalizing terms P(past) that factor out in normalization)

Now *if* (big if) the past events are all independent random events, then :

P(past|i) = P(x0|i)*P(x1|i)*P(x2|i)... for past symbols x0,x1,...

since P("ab") = P(a)*P(b) etc.

which you can compute incrementally as :

w_i *= P_i

the Beta Weight update.

In words, the weight for the expert should be the probability of that expert having made the data we have seen so far. That probability is just the P_i that expert has predicted for all past symbols, multiplied together.

Now this is obviously very wrong. Our symbols are very much not independent, they have strong conditional probability, that's how we get compression out of modeling the stream. The failure of that assumption may be why true beta weighting is poor in practice.

---

This Whirlwind Tour has been brought to you by the papers of Mattern & Knoll, PAQ of Mahoney, and the radness of RAD.

Small note about practical application :

You of course want to implement this with integers and store things like weights and probabilities in fixed point. I have found that the exact way you handle truncation of that fixed point is very important. In particular you need to be careful about rounding & biasing and try not to just lose small values. In general mixing appears to be very "fiddly"; it has a very small sweet spot of tuned parameters and small changes to them can make results go very bad.

Whirlwind Tour of Mixing Part 2

Before we dive into online gradient descent, let's look at a couple of building blocks.

In many cases we've been talking about mixers that work with arbitrary alphabets, but let's specialize now to binary. In that case, each expert just makes one P, say P(0), and implicitly P(1) = 1 - P(0).

We will make use of the logit ("stretch") and logistic sigmoid ("squash"), so let's look at them a bit first.

logit(x) = ln( x / (1-x) )

sigmoid(x) = 1 / (1 + e^-x )  (aka "logistic")

logit is "log odds" of a probability

logit takes [0,1] -> [-inf,inf]

sigmoid takes [-inf,inf] -> [0,1]

logit( sigmoid(x) ) = 1

I'm a little bit fast and loose with the base of the logs and exponents sometimes, because the difference comes down to scaling factors which wash out in constants in the end. I'll try to be not too sloppy about it.

You can see plots of them here for example . The "logistic sigmoid" is an annoyingly badly named function; "logistic" is way too similar to "logit" which is very confusion, and "sigmoid" is just a general shape description, not a specific function.

What's more intuitive in data compression typically is the base-2 variants :

log2it(P) = log2( P / (1-P) )

log2istic(t) = 1 / (1 + 2^-t )

which is just different by a scale of the [-inf,inf] range. (it's a difference of where you measure your codelen in units of "bits" or "e's").

In data compression there's an intuitive way to look at the logit. It's the difference of codelens of symbol 0 and 1 arising from the probability P.

log2it(P) = log2( P / (1-P) )

= log2( P ) - log2( 1 - P )

if P = P0

then L0 = codelen of 0 = - log2(P0)

L1 = - log2(1-P0)

log2it(P0) = L1 - L0

which I like to write as :

log2it = !L - L

where !L means the len of the opposite symbol

Working with the logit in data compression is thus very natural and has nice properties :

1. It is (a)symmetric under exchange of the symbol alphabet 0->1. In that case P -> 1-P , so logit -> - logit. The magnitude stays the same, so mixers that act on logit behave the same way. eg. you aren't favoring 0 or 1, which of course you shouldn't be. This is why using the *difference* of codelen between L0 and L1 is the right thing to do, not just using L0, since the 0->1 exchange acts on L0 in a strange way.

Any mixer that we construct should behave the same way if we swap 0 and 1 bits. If you tried to construct a valid mixer like that using only functions of L0, it would be a mess. Doing it in the logit = (!L-L) it works automatically.

2. It's linear in codelen. Our goal in compression is to minimize total codelen, so having a parameter that is in the space we want makes that natural.

3. It's unbounded. Working with P's is awkward because of their finite range, which requires clamping and theoretically projection back into the valid cube.

4. It's centered on zero for P = 0.5 and large for skewed P. We noted before that what we really want is to gather models that are skewed, and logit does this naturally.

Let's look at the last property.

Doing mixing in logit domain implicitly weights skewed experts, as we looked at before in "functional weighting". With logit, skewed P's near 0 or 1 can stretched out to -inf or +inf. This gives them very high implicit "weight", that is even after mixing with a P near 0.5, the result is still skewed.

Mix two probabilities with equal weight :

0.5 & 0.9

linear : 0.7
logit  : 0.75

0.5 & 0.99

linear : 0.745
logit  : 0.90867

0.5 & 0.999

linear : 0.74950
logit  : 0.96933

Logit blending is sticky at skewed P.

A little plot of what the L -> !L function looks like :

Obviously it is symmetric around (1,1) where both codelens are 1 bit, P is 0.5 ; as P gets skewed either way one len goes to zero and the other goes to infinity.

logistic sigmoid is the map back from codelen delta to probability. It has the nice property for us that no matter what you give it, the output is a valid normalized probability in [0,1] ; you can scale or screw up your logistics in whatever way you want and you still get a normalized probability (which makes it codeable).

Linear combination of logits (codelen deltas) is equal to geometric (multiplicative) mixing of probabilities :

M = mixed probability

linear combination in logit space, then transform back to probability :

M = sigmoid( Sum{ c_i * logit( P_i ) } )

M = 1 / ( 1 + e ^ - Sum{ c_i * logit( P_i ) } )

e ^ logit(P) = P/(1-P)
e ^ -logit(P) = (1-P)/P
e ^ - c * logit(P) = ((1-P)/P) ^ c

M = 1 / ( 1 + Prod{ (( 1 - P_i ) / P_i)^c_i } )

M = Prod{ P_i ^ c_i } / ( Prod{ P_i ^ c_i } + Prod{ ( 1 - P_i )^c_i } )

What we have is a geometric blend (like a geometric mean, but with arbitrary exponents, which are the mixing factors), and the denomenator is just the normalizing factor so that M0 and M1 some to one.

Note that the blend factors (c_i) do not need to be normalized here. We get a normalized probability even if the c_i's are not normalized. That is an important property for us in practice.

This is a classic machine learning method called "product of experts". In the most basic form, all the c_i = 1, the powers are all 1, the probabilities from each expert are just multiplied together.

I'm calling the c_i "coefficients" not "weights" to emphasize the fact that they are not normalized.

Note that while we don't need to normalize the c_i , and in practice we don't (and in the "product of experts" usage with all c_i = 1 they are of course not normalized), that has weird consequences and let us not brush it under the table.

For one thing, the overall scale of c_i changes the blend. In many cases we hand-waved about scaling because normalization will wipe out any overall scaling factors of weights. In this case the magnitude scale of the c_i absolutely matters. Furthermore, the scale of the c_i has nothing to do with a "learning rate" or blend strength. It's just an extra exponent applied to your probabilities that does some kind of nonlinear mapping to how they are mixed. e.g. doing c_i *= 2 is like mixing squared probabilities instead of linear probabilities.

For another, with unnormalized c_i - only skewed probabilities contribute. Mixing in more 50/50 experts does *nothing*. It does not bias your result back towards 50/50 at all!

We can see this both in the logit formulation and in the geometric mean formulation :

M = sigmoid( Sum{ c_i * logit( P_i ) } )

Now add on a new expert P_k with some coefficient c_k

M' = sigmoid( Sum{ c_i * logit( P_i ) + c_k * logit( P_k ) } )

if P_k = 50% , then logit(P_k) = 0

M' = M

even if c_k is large.

(if you normalized the c's, the addition of c_k would bring down the coffecient from the other contributions)

The Sum of logits only adds up codelen deltas

Experts with large codelen deltas are adding their scores to the vote
small codelen deltas simply don't contribute at all

In the geometric mean :

Q = Prod{ P_i ^ c_i }

M = Q / ( Q + !Q )

(where ! means "of the opposing symbol; eg. P -> (1-P) )

Add on another expert with P = 0.5 :

Q' = Q * 0.5 ^ c_k

!Q' = !Q * 0.5 ^ c_k

M' = Q * 0.5 ^ c_k / ( Q * 0.5 ^ c_k + !Q * 0.5 ^ c_k)

the same term multiplied around just factors out

M' = M

This style of mixing is like adding experts, not blending them. Loud shouty experts with skewed codelen deltas add on.

Next up, how to use this for gradient descent mixing in compression.