Introduction to Neural Networks
U. Minn. Psy 5036
Daniel Kersten
Belief Propagation

Initialize

Read in Statistical Add-in packages:

Off[General::spell1];
<< Statistics`NormalDistribution`

<< Statistics`ContinuousDistributions`
<< Statistics`MultinormalDistribution`

<< Graphics`MultipleListPlot`

Interpolation using smoothness revisited: Gradient descent

For simplicity, we'll assume 1-D as in the lecture on sculpting the energy function. In anticipation of formulating the problem in terms of a graph that represents conditional probability dependence, we represent observable depth cues by \(TraditionalForm\`\(y\^*\)\), and the true ("hidden") depth estimates by y.

(Figure from Weiss (1999).

First-order smoothness

Recall that the energy or cost function is given by:

where (Subscript[w, k])= xs[[k]] is the indicator function, and \(TraditionalForm\`\(y\_k\%*\)\)= d, are the data values.
Gradient descent gives the following local update rule:

As before, λ controls the degree of smoothness, i.e. smoothness at the expense of fidelity to the data.
Gauss-Seidel: η[k_]:=1/(λ+xs[[k]]);
Successive over-relaxation (SOR): η2[k_]:=1.9/(λ+xs[[k]]);

A simulation: Straight line with random missing data points

Make the data

We return to the problem of interpolating a set of points with missing data, marked by an indicator function xs.

size = 32;
xs = Table[0, {i,1,size}]; xs[[1]]=1;xs[[size]]=1;(*xs[[size/2]]=1;*)
data = Table[N[j] xs[[j]],
{j, 1, size}];
g3 = ListPlot[Table[N[j],{j,1,size}],PlotJoined->True,
    DisplayFunction->Identity,PlotStyle->{RGBColor[0,.5,0]}];
g2 = ListPlot[data,PlotJoined->False,
PlotStyle->{RGBColor[.75,.0,0]}, Prolog-> AbsolutePointSize[5],
DisplayFunction->Identity];

The green line shows the a straight line connecting the three data points. The red dots on the abscissa mark the points where data is missing.

Show[g2,g3,DisplayFunction->$DisplayFunction];

Let's set up two matrices, Tm and Sm such that the gradient of the energy is equal to:
 Tm . f - Sm . f.  
 Sm will be our filter to exclude non-data points. Tm will express the "smoothness" constraint.

Sm = DiagonalMatrix[xs];
Tm = Table[0,{i,1,size},{j,1,size}];
For[i=1,i<=size,i++,Tm[[i,i]] = 2];
Tm[[1,1]]=1;Tm[[size,size]]=1; (*Adjust for the boundaries*)
For[i=1,i<size,i++, Tm[[i+1,i]] = -1];
For[i=1,i<size,i++, Tm[[i,i+1]] = -1];

Check the update rule code for small size=10:

Clear[f, d, λ]
(λ*Tm . Array[f, size] - 
          Sm . ((Array[d, size]) - Array[f, size])) // 
      MatrixForm

Run gradient descent

Clear[Tf,f1];
dt = 1; λ=2;
Tf[f1_] := f1 - dt*(1/(λ+xs))*(Tm.f1 - λ*Sm.(data-f1));

We will initialize the state vector to zero, and then run the network for iter iterations:

iter=256;
f = Table[0,{i,1,size}];
result = Nest[Tf,f,iter];

Now plot the interpolated function.

g1 = ListPlot[result,PlotJoined->True, AspectRatio->Automatic,PlotRange->{{1,size},{1,size}}, DisplayFunction->Identity];
Show[{g1,g2,g3,Graphics[{Text["Iteration="<>ToString[iter],{size/2,size/2}]}]}, DisplayFunction->$DisplayFunction,PlotRange->{0,40}];

Belief Propagation

Same interpolation problem, but now using belief propagation

Example is taken from Yair Weiss.(Weiss, 1999)

Probabilistic generative model

y[[i]] = xs[[i]] data[[i]]  + dnoise,  
    dnoise~N[0, (Subscript[σ, D])]

y[[i + 1]] = y[[i]] + znoise,  
    znoise~N[0, (Subscript[σ, R])]

Some theory

We'd like to know the distribution of each node i, conditioned on all the data: I.e. we want the posterior
	p(yi=u |all the data) 
If we could find this, we'd be able to: 1) say what the most probable value of the y value is, and 2) give a measure of confidence

Let p(yi=u |all the data) be normally distributed: NormalDistribution[μi,σi].

Consider the ith unit. The posterior p(yi=u|all the data) =

p(yi=u|all the data) ∝ p(yi=u|data before i) p(data at i| yi=u) p(yi=u|data after  i)

Suppose that p(yi=u|data before i) is also gaussian:  
						p(yi=u|data before i) = α[u] ∼ NormalDistribution[μα,σα]
and so is probability conditioned on the data after i:
						p(yi=u|data after  i)= β[u] ∼ NormalDistribution[μβ,σβ]

And the noise model for the data: 
						p(data at i| yi=u) = L[u]∼ NormalDistribution[yp, (Subscript[σ, D])]
						yp=data[[i]]

So in terms of these functions, the posterior probability of the ith unit taking on the value u can be expressed as proportional to a product of the three factors:

p(yi=u|all the data) ∝ α[u]*L[u]*β[u]

αudist = 
        NormalDistribution[μα, σα];
α[u] = 
        PDF[αudist, 
          u];

Ddist = 
        NormalDistribution[(Subscript[y, p]), (Subscript[σ, D])];
L[u] = PDF[Ddist, u];

βudist = 
        NormalDistribution[μβ, σβ];
β[u] = 
        PDF[βudist, u];

α[u]*L[u]*β[u]

This just another gaussian distribution on yi=u. What is its mean and variance? Finding the root enables us to complete the square to see what the numerator looks like. In particular, what the mode (=mean for gaussian) is.

Solve[-D[-((((u - μα))^2)/(2 σ(α^2))) - ((((u - μβ))^2)/(2 σ(β^2))) - ((((u - (Subscript[y, p])))^2)/(2 (Subscript[σ, D]^2))), 
            u] == 0, u]

The update rule for the variance is:

\(σ\^2 -> 
      1\/σα\^2 + 1\/σβ\^2 + 
        1\/σ\_D\%2\)

How do we get μα, μβ , σα, μβ?
We express the probability of the ith unit taking on the value u in terms of the values of the neighbor  before, conditioning on what is known (the observed measurements), and marginalizing over what isn't (the previous "hidden" node value, v).

We have three terms to worry about that depend on nodes in the neighborhood preceding i:

α[u]=(∫_(-∞)^∞)(Subscript[α, p])[v]*S[u]*
            L[v] dv∝(∫_(-∞)^∞)\(\(e\^\(\(-\(\((v - y\_p)\)\^2\/\(2\ σ\_D\%2\)\)\) - \((u - v)\)\^2\/\(2\ σ\_R\%2\) - \((v - μα\_p)\)\^2\/\(2\ σα\_p\%2\)\)\) dv\)

(Subscript[α, p]) = (Subscript[α, i - 1]) .  S[u] is our smoothing term,  
    or transition probability :  S[u] = 
      p (u | v) .  
          L[] is the likelihood of the previous data node,  
    given its hidden node value,  v.

Rdist = 
        NormalDistribution[v, (Subscript[σ, R])];
S[u] = PDF[Rdist, u];

αvdist = 
        NormalDistribution[μ(Subscript[α, p]), σ(Subscript[α, p])];
(Subscript[α, p])[v] = 
        PDF[αvdist, 
          v];

Lp[v] = PDF[Ddist, v];

Integrate[(Subscript[α, p])[v]*S[u]*L[v], {v, -Infinity, 
        Infinity}]

Some uninspired Mathematica manipulations

D[-(((((u - μ(Subscript[α, p])))^2) (Subscript[σ, D]^2) + μ(Subscript[α, p]^2) (Subscript[σ, R]^2) + (u^2) σ(Subscript[α, p]^2) + 
              (Subscript[y, p]^2) ((Subscript[σ, R]^2) + σ(Subscript[α, p]^2)) - 
              2 (Subscript[y, p]) (μ(Subscript[α, p]) (Subscript[σ, R]^2) + 
                    u σ(Subscript[α, p]^2)))/(2 ((Subscript[σ, R]^2) σ(Subscript[α, p]^2) + (Subscript[σ, D]^2) ((Subscript[σ, R]^2) + σ(Subscript[α, p]^2))))), u]

Solve[-% == 0, u]

Simplify[(((μ(Subscript[α, p]) (Subscript[σ, D]^2))/((Subscript[σ, R]^2) σ(Subscript[α, p]^2) + (Subscript[σ, D]^2) ((Subscript[σ, R]^2) + σ(Subscript[α, p]^2)))) + (((Subscript[y, p]) σ(Subscript[α, p]^2))/((Subscript[σ, R]^2) σ(Subscript[α, p]^2) + (Subscript[σ, D]^2) ((Subscript[σ, R]^2) + σ(Subscript[α, p]^2)))))/((Subscript[σ, D]^2)*σ(Subscript[α, p]^2))]

Simplify[(((Subscript[σ, D]^2)/((Subscript[σ, R]^2) σ(Subscript[α, p]^2) + (Subscript[σ, D]^2) ((Subscript[σ, R]^2) + σ(Subscript[α, p]^2)))) + σ((Subscript[α, p]^2)/((Subscript[σ, R]^2) σ(Subscript[α, p]^2) + (Subscript[σ, D]^2) ((Subscript[σ, R]^2) + σ(Subscript[α, p]^2)))))/((Subscript[σ, D]^2)*σ(Subscript[α, p]^2))]

(((μ(Subscript[α, p]) (Subscript[σ, D]^2) + (Subscript[y, p]) σ(Subscript[α, p]^2))/((Subscript[σ, D]^2) (Subscript[σ, R]^2) σ(Subscript[α, p]^4) + (Subscript[σ, D]^4) σ(Subscript[α, p]^2) ((Subscript[σ, R]^2) + σ(Subscript[α, p]^2)))))/((((Subscript[σ, D]^2) + σ(Subscript[α, p]^2))/((Subscript[σ, D]^2) (Subscript[σ, R]^2) σ(Subscript[α, p]^4) + (Subscript[σ, D]^4) σ(Subscript[α, p]^2) ((Subscript[σ, R]^2) + σ(Subscript[α, p]^2)))))

So we now have rule that tells us how to update the α(u)=p(yi=u|data before i), in terms of the mean and variance parameters of the previous node:

\(μα ← \(μα\_p\ σ\_D\%2 + y\_p\ σα\_p\%2\)\/\(σ\_D\%2 + σα\_p\%2\) = \(\(\(μα\_p\ σ\_D\%2\)\/\(σα\_p\%2\ σ\_D\%2\) + \(y\_p\ σα\_p\%2\)\/\(σα\_p\%2\ σ\_D\%2\)\)\/\(σ\_D\%2\/\(σα\_p\%2\ σ\_D\%2\) + σα\_p\%2\/\(σα\_p\%2\ σ\_D\%2\)\) = \(μα\_p\/σα\_p\%2 + y\_p\/\(\(\ \)\(σ\_D\%2\)\)\)\/\(1\/σα\_p\%2 + 1\/σα\_p\%2\)\)\)

The update rule for the variance is:

\(σα\^2 ← σ\_R\%2 + 
        1\/\(1\/σ\_D\%2 + 1\/σα\_p\%2\)\)

A similar derivation gives us the rules for μα, σ(α^2)

μβ←\(\(μβ\_a\/σβ\_a\%2 + y\_a\/\(\(\ \)\(σ\_D\%2\)\)\)\/\(1\/σβ\_a\%2 + 1\/σα\_a\%2\)\)

\(σβ\^2 ← σ\_R\%2 + 
        1\/\(1\/σ\_D\%2 + 1\/σβ\_a\%2\)\)

Where the subscript index  p (for "previous", i.e. unit i-1) is replaced by a (for "after", i.e. unit i+1).

Recall that sometimes we have data and sometimes we don't. So replace:

\(y\_p → 
        xs[i - 1] 
          data[i - 1] = \(w\_\(i - 1\)\) \(y\_\(i - 
              1\)\%*\)\)

And similarly for (Subscript[y, a]).

Summary of update rules

The ratio, (((((Subscript[σ, D])/(Subscript[σ, R]))))^2)plays the role of λ above. If (Subscript[σ, D]^2)>>(Subscript[σ, R]^2), there is greater smoothing. If (Subscript[σ, D]^2)<<(Subscript[σ, R]^2), there is more fidelity to the data. (Recall \(TraditionalForm\`\(y\^*\) -> \(\(data\)\(.\)\)\)(Subscript[w, k]) -> xs[[k]])

We'll follow Weiss, and  also make a (hopefully not too confusing) notation change to avoid the square superscripts for (Subscript[σ, D]^2)->(Subscript[σ, D]), (Subscript[σ, R]^2)->(Subscript[σ, R]).

A simulation: Belief propagation for interpolation with missing data

Initialization

μ0 = 1; 
μα = 1; σα = 
      100000;    
μβ = 
      1; σβ = 100000;   
σR = 
      4.0; σD = 1.0;
µ = Table[μ0, {i, 1, size}];
σ = 
        Table[σα, {i, 1, size}];
µα = 
        Table[μ0, {i, 1, size}];
σα = 
        Table[σα, {i, 1, size}];
µβ = 
        Table[μ0, {i, 1, size}];
σβ = 
        Table[σβ, {i, 1, size}];
iter = 0;
i = 1; j = size;

The code below implements the above iterative equations, taking care near the boundaries. The plot shows the estimates of (Subscript[y, i])=  µ, and the error bars show ±(Subscript[σ, i]).

Belief Propagation Routine: Execute this cell "manually" for each iteration

µ[[i]] = ((xs[[i]]/σD*data[[i]] + (1/σ)α[[i]]*µα[[i]] + (1.0/σ)β[[i]]*µβ[[i]])/(xs[[i]]/σD + (1/σ)α[[i]] + (1/σ)β[[i]]));
σ[[i]] = 
        (1.0/(xs[[i]]/σD + (1/σ)α[[i]] + (1/σ)β[[i]]));
µ[[j]] = ((xs[[j]]/σD*data[[j]] + (1/σ)α[[j]]*µα[[j]] + (1.0/σ)β[[j]]*µβ[[j]])/(xs[[j]]/σD + (1/σ)α[[j]] + (1/σ)β[[j]]));
σ[[j]] = 
        (1.0/(xs[[j]]/σD + (1/σ)α[[j]] + (1/σ)β[[j]]));


nextj = j - 1;
µα[[nextj]] = ((xs[[j]]/σD*data[[j]] + (1.0/σ)α[[j]]*µα[[j]])/(xs[[j]]/σD + (1/σ)α[[j]]));
σα[[nextj]] = σR + 
          (1.0/(xs[[j]]/σD + (1/σ)α[[j]]));

nexti  =  i + 1;
µβ[[nexti]] = ((xs[[i]]/σD*data[[i]] + (1.0/σ)β[[i]]*µβ[[i]])/(xs[[i]]/σD + (1/σ)β[[i]]));
σβ[[nexti]] = σR + 
          (1.0/(xs[[i]]/σD + (1/σ)β[[i]]));

j--;
i++;

iter++;

yfit = 
        Table[{µ[[i1]], ErrorBar[σ[[i1]]]}, {i1, 1, 
            size}];
g1b = 
        MultipleListPlot[yfit, 
          DisplayFunction -> Identity];
Show[{g1b, g2, g3, 
          Graphics[{Text["Iteration=" <> ToString[iter], {size/2, 
                  size}]}]}, DisplayFunction -> $DisplayFunction, 
        PlotRange -> {-40, 40}];

Expectation Maximization (EM) -- A preview & simulation

We've studied the problem of interpolation given missing data. We motivated the problem by the visual phenomenon of surface completion. Another aspect of surface perception is our ability to take noisy data (e.g. depth cues), and not only interpolate the data, but also decide which of several surfaces the data belong using stereo data (Madarasmi et al., 1993; Kersten & Madarasmi, 1995) or optic flow. A number of ways have been proposed to deal with this problem, but the technique of most general application is expectation-maximization--a general statistical technique developed in the 1970's that appears in various forms in many algorithms, including belief propagation. The algorithm has been applied to the surface estimation problem, e.g. from optic flow (Jepson, 1993; Weiss, 1997). We go to Weiss again for a nice simple demo.
Consider Generative Model 1.

Generative models

Two lines with slopes and intercepts (a1,b1) and (a2,b2).

ndist = NormalDistribution[0, 1];

Generative model 1

y1[x_] :=  -2*x  + 40 + 
          Random[ndist];
y2[x_] := 3*x + 1 + Random[ndist];
data = 
        Table[{x, If[Abs[x - 3] < 1.5, y1[x], y2[x]]}, {x, 0, 
            6,  .1}];
{x, y} = Transpose[data];

Generative model 2

y1[x_] :=  -2*x  + 15;
y2[x_] := 3*x + 1;
data = 
        Table[{x, If[Random[] <  .5, y1[x], y2[x]]}, {x, 0, 
            6,  .1}];
{x, y} = Transpose[data];

gdata = ListPlot[data];

Which data belong together? We'll assume that we know there are two linear models, but we don't know their parameters--i.e. we don't know their slopes and intercepts.

EM algorithm

Initialize parameters to random values

σ =  .1;
{a1, b1, a2, b2} = Table[10*Random[], {4}];

E-step

Compute residuals  r1, r2, the error in the predicted and actual y values under each of the two models.

r1 = a1*x  +  b1 -  y;
r2 = a2*x  +  b2 -  y;

Using the residuals, compute weights. We'll assing these weights to data in the M-step later.

w1 = 
        Exp[-r1^2/σ]/(Exp[-r1^2/σ] + 
              Exp[-r2^2/σ]);
w2 = 
        Exp[-r2^2/σ]/(Exp[-r1^2/σ] + 
              Exp[-r2^2/σ]);

M-step

With above, weights compute weighted linear regression:

{a1, b1} = 
        Inverse[{{w1 . (x x), w1 . x}, {w1 . x, 
                Apply[Plus, w1]}}] . {w1 . (x y), 
            w1 . y};
{a2, b2} = 
        Inverse[{{w2 . (x x), w2 . x}, {w2 . x, 
                Apply[Plus, w2]}}] . {w2 . (x y), w2 . y};

gfit = 
        Plot[{a1*x + b1, a2*x + b2}, {x, 0, 6}, 
          DisplayFunction -> Identity];
Show[{gfit, 
          Graphics[{PointSize[0.03], 
              Transpose[{Hue[#1] & /@ (w1/2), Point /@ data}]}]}, 
        PlotRange -> All, 
        DisplayFunction -> $DisplayFunction];

Run the descent algorithm using successive over-relaxation (SOR): η2[k_]:=1.9/(λ+xs[[k]]). 
How does convergence compare with Gauss-Seidel?

Run Belief Propation  using: σR=1.0; σD=4.0; How does fidelity to the data compare with the original case 
(σR=4.0; σD=1.0).

BP with missing sine wave data

Generate sine wave with missing data

size = 64;
xs = Table[Random[Integer,1], {i,1,size}];
data = Table[N[Sin[2 Pi (1/20) j] xs[[j]]],
{j, 1, size}];
g3b = ListPlot[Table[N[Sin[2 Pi (1/20) j]], {j, 1, size}],
		PlotJoined->True,DisplayFunction->Identity,
	PlotStyle->{RGBColor[0,.5,0]}];
g2b = ListPlot[data,PlotJoined->False,
PlotStyle->{RGBColor[.75,.0,0]}, Prolog-> AbsolutePointSize[5],
DisplayFunction->Identity];
Show[g2b,g3b,DisplayFunction->$DisplayFunction];

Initialize

μ0 = 1; 
μα = 1; σα = 
      100000;    
μβ = 
      1; σβ = 100000;  
σR =  .5; σD =  .1;
µ = Table[μ0, {i, 1, size}];
σ = 
        Table[σα, {i, 1, size}];
µα = 
        Table[μ0, {i, 1, size}];
σα = 
        Table[σα, {i, 1, size}];
µβ = 
        Table[μ0, {i, 1, size}];
σβ = 
        Table[σβ, {i, 1, size}];
iter = 0;
i = 1; j = size;

SINE WAVE DEMO: Belief Propagation Routine: Execute this cell "manually" for each iteration

µ[[i]] = (xs[[i]]/σD*data[[i]] + 
              (1/σ)α[[i]]*µα[[i]] + 
              (1.0/σ)β[[i]]*µβ[[i]])/(xs[[i]]/σD + (1/σ)α[[i]] + 
              (1/σ)β[[i]]);
σ[[i]] = 
        (1.0/(xs[[i]]/σD + (1/σ)α[[i]] + (1/σ)β[[i]]));
µ[[j]] = (xs[[j]]/σD*data[[j]] + 
              (1/σ)α[[j]]*µα[[j]] + 
              (1.0/σ)β[[j]]*µβ[[j]])/(xs[[j]]/σD + (1/σ)α[[j]] + 
              (1/σ)β[[j]]);
σ[[j]] = 
        (1.0/(xs[[j]]/σD + (1/σ)α[[j]] + (1/σ)β[[j]]));


nextj = j - 1;
µα[[nextj]] = ((xs[[j]]/σD*data[[j]] + (1.0/σ)α[[j]]*µα[[j]])/(xs[[j]]/σD + (1/σ)α[[j]]));
σα[[nextj]] = σR + 
          (1.0/(xs[[j]]/σD + (1/σ)α[[j]]));

    
nexti  =  i + 1;
µβ[[nexti]] = ((xs[[i]]/σD*data[[i]] + (1.0/σ)β[[i]]*µβ[[i]])/(xs[[i]]/σD + (1/σ)β[[i]]));
σβ[[nexti]] = σR + 
          (1.0/(xs[[i]]/σD + (1/σ)β[[i]]));

    
j--;
i++;

    
iter++;

    
yfit = 
        Table[{µ[[i1]], ErrorBar[σ[[i1]]]}, {i1, 1, 
            size}];
g1bb = 
        MultipleListPlot[yfit, 
          DisplayFunction -> Identity];
Show[{g1bb, g2b, g3b}, DisplayFunction -> $DisplayFunction, 
        PlotRange -> {-2, 2}];

Run EM with Generative Model 2. Increase the additive noise. How does attribution accuracy change? ("attribution" means assigning a point to its correct line)

References

Applebaum, D. (1996). Probability and Information . Cambridge, UK: Cambridge University Press.
Frey, B. J. (1998). Graphical Models for Machine Learning and Digital Communication. Cambridge, Massachusetts: MIT Press.
Jepson, A., & Black, M. J. (1993). Mixture models for optical flow computation. Paper presented at the Proc. IEEE Conf. Comput. Vsion Pattern Recog., New York.
Kersten, D. and P.W. Schrater (2000), Pattern Inference Theory: A Probabilistic Approach to Vision, in Perception and the Physical World, R. Mausfeld and D. Heyer, Editors. , John Wiley & Sons, Ltd.: Chichester. (pdf)
Kersten, D., & Madarasmi, S. (1995). The Visual Perception of Surfaces, their Properties, and Relationships. DIMACS Series in Discrete Mathematics and Theoretical Computer Science, 19, 373-389.
Madarasmi, S., Kersten, D., & Pong, T.-C. (1993). The computation of stereo disparity for transparent and for opaque surfaces. In C. L. Giles & S. J. Hanson & J. D. Cowan (Eds.), Advances in Neural Information Processing Systems 5. San Mateo, CA: Morgan Kaufmann Publishers.
Pearl, Judea. (1997) Probabilistic Reasoning in Intelligent Systems : Networks of Plausible Inference. (amazon.com link)
Ripley, B. D. (1996). Pattern Reco gnition and Neural Networks. Cambridge, UK: Cambridge University Press.
Weiss Y. (1999) Bayesian Belief Propagation for Image Understanding submitted to SCTV 1999. (gzipped postscript 297K) 
Weiss, Y. (1997). Smoothness in Layers: Motion segmentation using nonparametric mixture estimation. Paper presented at the Proceedings of IEEE conference on Computer Vision and Pattern Recognition.
Yuille, A., Coughlan J., Kersten D.(1998) (pdf)

For notes on Graphical Models, see:http://www.cs.berkeley.edu/~murphyk/Bayes/bayes.html