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

Initialize

Off[General::spell1];
Needs["ErrorBarPlots`"]

Last time

Neural population codes representing probability distributions

Generative modeling: Multivariate gaussian, mixtures

Drawing samples

Mixtures of gaussians

Will use mixture distributions in the next lecture on EM application to segmentation

Introduction to Optimal inference and task

Introduction to Bayesian learning

Interpolation using smoothness: 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 (SuperStar[y]), and the true ("hidden") depth estimates by y.

(Figure from Weiss (1999).

First-order smoothness

Earlier we saw that under specific assumptions, biologically plausible neural updating can be seen decrease the value of an energy (or "cost") function. One can also start off with an assumed cost function, determined say by a set of contraints, and use gradient descent to derive an update rule that minimizes the cost.  (See supplementary material in Lecture 20). However, such an update rule will not necessarily resemble a biologically plausible neural mechanism.

For an interpolation problem, we can write the energy or cost function by:

where (Subscript[w, k])= xs[[k]] is an "indicator function", and (Subscript[y, k]^( * ))= d, are the data values. The indicator function is 1 if there is data available, and zero otherwise. . 
Gradient descent gives the following local update rule:

λ is a free parameter that 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

Consider the problem of interpolating a set of points with missing data, marked by an indicator function with the following notation: 
(Subscript[w, k])= xs[[k]], y* = data, y=f.
We'll assume the true model is that f=y=j, where j=1 to size.  data is a function of the sampling process on  f=j

size=32;xs=Table[0,{i,1,size}];xs[[1]]=1;xs[[size]]=1;data=Table[N[j] xs[[j]],{j,1,size}];g3=ListPlot[Table[N[j],{j,1,size}],Joined->True,PlotStyle->{RGBColor[0,0.5,0]}];g2=ListPlot[data,Joined->False,PlotStyle->{Opacity[0.35],RGBColor[0.75,0.,0],PointSize[Large]}];

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

Show[{g2,g3},ImageSize->Medium]

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:

f0 = Table[0,{i,1,size}];
result=f0;
f=f0;
iter=25;

Now plot the interpolated function.

Animate[
result=Nest[Tf,f,iter];
f=result;
g1=ListPlot[result,Joined->False,AspectRatio->Automatic,PlotRange->{{0,size},{-1,size+1}},ImageSize->Medium];
Show[{g1,g2,g3,Graphics[{Text["Iteration="<>ToString[iter],{(size/2),(size/2)}]}]},PlotRange->{-1,size+1}],{j,1,10}]

Try starting with f = random values. Try various numbers of iterations.

Try different sampling functions xs[[i]].

Belief Propagation

Same interpolation problem, but now using belief propagation

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

Probabilistic generative model

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

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

The first term is the "data formation" model, i.e. how the data is directly influenced by the interaction of the underlying influences or causes: y with the sampling and noise. The second term reflects our prior assumptions about the smoothness of y, i.e. nearby y's are correlated, and in fact identical except for some added noise. So with no noise the prior reflects the assumption that lines are horizontal--all y's are the same.

Some theory

We'd like to know the distribution of the random variables at each node i, conditioned on all the data: I.e. we want the posterior
	p((Subscript[Y, i])=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((Subscript[Y, i])=u | all the data) be normally distributed: NormalDistribution[μi,σi].

Consider the ith unit. The posterior p((Subscript[Y, i])=u | all the data) =

p((Subscript[Y, i])=u| all the data) ∝ p((Subscript[Y, i])=u | data before i) p(data at i | (Subscript[Y, i])=u) p((Subscript[Y, i])=u | data after  i)

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

And the noise model for the data: 
						p(data at i | (Subscript[Y, i])=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((Subscript[Y, i])=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 (Subscript[Y, i])=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/(Subscript[σ, 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, at the i-1th location).

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-(Subscript[y, p])))^2)/(2 (Subscript[σ, D]^2)))-((((u-v))^2)/(2 (Subscript[σ, R]^2)))-((((v-(Subscript[μα, p])))^2)/(2 (Subscript[σα, 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]*Lp[v],{v,-Infinity,Infinity}]

Some uninspired Mathematica manipulations

To find an expression for the mode of the above calculated expression for α[u]

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:

μα←(((Subscript[μα, p]) (Subscript[σ, D]^2)+(Subscript[y, p]) (Subscript[σα, p]^2))/((Subscript[σ, D]^2)+(Subscript[σα, p]^2)))=(((((Subscript[μα, p]) (Subscript[σ, D]^2))/((Subscript[σα, p]^2) (Subscript[σ, D]^2)))+(((Subscript[y, p]) (Subscript[σα, p]^2))/((Subscript[σα, p]^2) (Subscript[σ, D]^2))))/(((Subscript[σ, D]^2)/((Subscript[σα, p]^2) (Subscript[σ, D]^2)))+((Subscript[σα, p]^2)/((Subscript[σα, p]^2) (Subscript[σ, D]^2)))))=((((Subscript[μα, p])/(Subscript[σα, p]^2))+((Subscript[y, p])/( (Subscript[σ, D]^2))))/((1/(Subscript[σα, p]^2))+(1/(Subscript[σα, p]^2))))

The update rule for the variance is:

((σα)^2)←(Subscript[σ, R]^2)+(1/((1/(Subscript[σ, D]^2))+(1/(Subscript[σα, p]^2))))

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

μβ←((((Subscript[μβ, a])/(Subscript[σβ, a]^2))+((Subscript[y, a])/( (Subscript[σ, D]^2))))/((1/(Subscript[σβ, a]^2))+(1/(Subscript[σα, a]^2))))

((σβ)^2)←(Subscript[σ, R]^2)+(1/((1/(Subscript[σ, D]^2))+(1/(Subscript[σβ, 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:

(Subscript[y, p])→xs[i-1]data[i-1]=(Subscript[w, i-1])(Subscript[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 (SuperStar[y])->data.(Subscript[w, k])->xs[[k]]). But now we have a principled way of assigning the relative amount of smoothing.

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

size=32;
μ0=1; 
μα=1;σα=100000; (*large uncertainty *)
μβ=1;σβ=100000;(*large*)
σ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

yfit=Table[{0,0},{i1,1,size}];
g1b=ErrorListPlot[{yfit}];
Dynamic[Show[{g1b,g2,g3,Graphics[{Text["Iteration="<>ToString[iter],{(size/2),size}]}]},PlotRange->{-50,50},Axes->{False,True}]]

Execute the next cell to run 31 iterations. The display is slowed down so that you can see the progression of the updates in the above graph.

Do[
Pause[.5];
µ[[i]]=((((xs[[i]] data[[i]])/(σD))+((µα[[i]])/(σα[[i]]))+((1. µβ[[i]])/(σβ[[i]])))/(((xs[[i]])/(σD))+(1/(σα[[i]]))+(1/(σβ[[i]]))));
σ[[i]]=(1./(((xs[[i]])/(σD))+(1/(σα[[i]]))+(1/(σβ[[i]]))));
µ[[j]]=((((xs[[j]] data[[j]])/(σD))+((µα[[j]])/(σα[[j]]))+((1. µβ[[j]])/(σβ[[j]])))/(((xs[[j]])/(σD))+(1/(σα[[j]]))+(1/(σβ[[j]]))));
σ[[j]]=(1./(((xs[[j]])/(σD))+(1/(σα[[j]]))+(1/(σβ[[j]]))));
nextj=j-1;
µα[[nextj]]=((((xs[[j]] data[[j]])/(σD))+((1. µα[[j]])/(σα[[j]])))/(((xs[[j]])/(σD))+(1/(σα[[j]]))));
σα[[nextj]]=σR+(1./(((xs[[j]])/(σD))+(1/(σα[[j]]))));
nexti=i+1;
µβ[[nexti]]=((((xs[[i]] data[[i]])/(σD))+((1. µβ[[i]])/(σβ[[i]])))/(((xs[[i]])/(σD))+(1/(σβ[[i]]))));
σβ[[nexti]]=σR+(1./(((xs[[i]])/(σD))+(1/(σβ[[i]]))));
j--;
i++;
iter++;
yfit=Table[{µ[[i1]],σ[[i1]]},{i1,1,size}];
g1b=ErrorListPlot[{yfit}];
,{size-1}];

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[RandomInteger[1],{i,1,size}];data=Table[N[Sin[((2 π j)/20)] xs[[j]]],{j,1,size}];g3b=ListPlot[Table[N[Sin[((2 π j)/20)]],{j,1,size}],Joined->True,PlotStyle->{RGBColor[0,0.5,0]}];g2b=ListPlot[data,Joined->False,PlotStyle->{RGBColor[0.75,0.,0]}];

Initialize

μ0=1; 
μα=1;σα=100000; (*large uncertainty *)
μβ=1;σβ=100000;(*large*)
σ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;

yfit=Table[{0,0},{i1,1,size}];
g1bb=ErrorListPlot[{yfit}];
Dynamic[Show[{g1bb,g2b,g3b},PlotRange->{-2,2},Axes->{False,True}]]

SINE WAVE DEMO: Belief Propagation Routine

Do[
Pause[0.2];
µ[[i]]=((((xs[[i]] data[[i]])/(σD))+((µα[[i]])/(σα[[i]]))+((1. µβ[[i]])/(σβ[[i]])))/(((xs[[i]])/(σD))+(1/(σα[[i]]))+(1/(σβ[[i]]))));
σ[[i]]=(1./(((xs[[i]])/(σD))+(1/(σα[[i]]))+(1/(σβ[[i]]))));
µ[[j]]=((((xs[[j]] data[[j]])/(σD))+((µα[[j]])/(σα[[j]]))+((1. µβ[[j]])/(σβ[[j]])))/(((xs[[j]])/(σD))+(1/(σα[[j]]))+(1/(σβ[[j]]))));
σ[[j]]=(1./(((xs[[j]])/(σD))+(1/(σα[[j]]))+(1/(σβ[[j]]))));
nextj=j-1;
µα[[nextj]]=((((xs[[j]] data[[j]])/(σD))+((1. µα[[j]])/(σα[[j]])))/(((xs[[j]])/(σD))+(1/(σα[[j]]))));
σα[[nextj]]=σR+(1./(((xs[[j]])/(σD))+(1/(σα[[j]]))));
nexti=i+1;
µβ[[nexti]]=((((xs[[i]] data[[i]])/(σD))+((1. µβ[[i]])/(σβ[[i]])))/(((xs[[i]])/(σD))+(1/(σβ[[i]]))));
σβ[[nexti]]=σR+(1./(((xs[[i]])/(σD))+(1/(σβ[[i]]))));
j--;
i++;
iter++;
yfit=Table[{µ[[i1]],σ[[i1]]},{i1,1,size}];
g1bb=ErrorListPlot[{yfit}];
,{size-1}]

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