Answered by:
Linear Dynamical System (Migrated from community.research.microsoft.com)
Question

tommy1802 posted on 06292010 9:10 PM
Hello,
I tried to build a Kalmanfilter for a linear dynamical system on Infer.NET. But I can not implement the MatrixTimesVector state transition A*x. (where x is the previous state).
My two questions are the following:
1.) Is it generally possible to make inference on a LDS with Infer.net?
2.) What is the mistake in the following code (this is just a first trial version, there are also no observation matrices implemented). I've marked the line that causes the problem with the state transition.

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using MicrosoftResearch.Infer.Models;
using MicrosoftResearch.Infer;
using MicrosoftResearch.Infer.Maths;
using MicrosoftResearch.Infer.Distributions;
namespace LDS
{
class Program
{
// True state transition matrix
static double[,] trueA =
{
{ 0.7, 0.0 },
{ 0.0, 0.7 }
};
// True state transition precision matrix
static double[,] trueB =
{
{ 1.0, 0.0 },
{ 0.0, 1.0 }
};
// True observation matrix
static double[,] trueC =
{
{ 1.0, 0.0 },
{ 5.8, 4.0 }
};
// True observation precision matrix
static double[,] trueD =
{
{ 10.0, 0.0 },
{ 0.0, 10.0 }
};
// True initial state mean
static double[] trueMu0 = { 0.0, 0.0 };
// True inital state precision
static double[,] trueB0 =
{
{ 1.0, 0.0 },
{ 0.0, 1.0 }
};
// 
// Main Program
// 
static void Main(string[] args)
{
// Defining dimensions
int nSamples = 10;
int dimObs = trueC.GetLength(0);
int dimState = trueC.GetLength(1);
// Generate Data
// 
// Set a stable random number seed for repeatbale runs
Rand.Restart(12347);
// Generate Samples from an LDS
double[,] data = generateData(nSamples);
//printArrayToConsole(data);
// Inference
// 
Vector vMu0 = new Vector(trueMu0);
PositiveDefiniteMatrix mB0 = new PositiveDefiniteMatrix(trueB0);
Matrix mA = new Matrix(trueA);
Variable<Matrix> vA = Variable.Constant<Matrix>(mA);
PositiveDefiniteMatrix mB = new PositiveDefiniteMatrix(trueB);
Matrix mC = new Matrix(trueC);
PositiveDefiniteMatrix mD = new PositiveDefiniteMatrix(trueD);
/*Variable<PositiveDefiniteMatrix> mD = Variable.WishartFromShapeAndScale(10.0,
PositiveDefiniteMatrix.IdentityScaledBy(2, 0.001));*/
Variable<Vector>[] vState = new Variable<Vector>[nSamples];
Variable<Vector>[] vObservations = new Variable<Vector>[nSamples];
vState[0] = Variable.VectorGaussianFromMeanAndPrecision(vMu0, mB0);
vObservations[0] = Variable.VectorGaussianFromMeanAndPrecision(vState[0], mD);
double[] currentData = new double[dimObs];
for (int j = 0; j < dimObs; j++)
currentData[j] = data[0,j];
vObservations[0].ObservedValue = new Vector(currentData);
for (int i = 1; i < nSamples; i++)
{
//vState = Variable.VectorGaussianFromMeanAndPrecision(vState[i  1], mB);
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!! Next Line causes problems !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
vState = Variable.VectorGaussianFromMeanAndPrecision(Variable.MatrixTimesVector(vA, vState[i1]), mB);
vObservations = Variable.VectorGaussianFromMeanAndPrecision(vState, mD);for (int j = 0; j < dimObs; j++)
currentData[j] = data[i, j];
vObservations.ObservedValue = new Vector(currentData);
}
InferenceEngine engine = new InferenceEngine(new VariationalMessagePassing());
//engine.ShowFactorGraph = true;
//engine.ShowTimings = true;
Console.WriteLine("State[" + Convert.ToString(nSamples  1) + "] = " + engine.Infer(vState[nSamples  1]));
Console.ReadLine();
}
/// Generate data from the true model
static double[,] generateData(int nSamples)
{
int dimObs = trueC.GetLength(0);
int dimState = trueC.GetLength(1);
double[,] data = new double[nSamples, dimObs];
Vector vMu0 = new Vector(trueMu0);
PositiveDefiniteMatrix mB0 = new PositiveDefiniteMatrix(trueB0);
Matrix mA = new Matrix(trueA);
PositiveDefiniteMatrix mB = new PositiveDefiniteMatrix(trueB);
Matrix mC = new Matrix(trueC);
PositiveDefiniteMatrix mD = new PositiveDefiniteMatrix(trueD);
Vector vStateTemp = new Vector(dimState);
Vector vObsTemp = new Vector(dimObs);
// Sample from initial state and corresponding observation
vStateTemp = VectorGaussian.Sample(vMu0, mB0);
vObsTemp = VectorGaussian.Sample(mC * vStateTemp, mD);
// Saving the first Observation
for (int j = 0; j < dimObs; j++)
data[0,j] = vObsTemp[j];
for (int i = 1; i < nSamples; i++)
{
// Sampling from new state and corresponding observation
vStateTemp = VectorGaussian.Sample(mA * vStateTemp, mB);
vObsTemp = VectorGaussian.Sample(mC * vStateTemp, mD);
// Saving the observations at the current time instant
for (int j = 0; j < dimObs; j++) data[i, j] = vObsTemp[j];
}
return data;
}
/// Print a 2D double array to the console
private static void printArrayToConsole(double[,] matrix)
{
for (int i = 0; i < matrix.GetLength(0); i++)
{
for (int j = 0; j < matrix.GetLength(1); j++)
Console.Write("{0,5:0.00}\t", matrix[i, j]);
Console.WriteLine();
}
}
}
}
Thanks
Thomas
PS: sorry for the uncolored code. I tried to usw Word import, but the color is also gone. And also some symbols are transfered into bulbs. The Editor does not offer code import. Or can someone tell me how to do this better.
Friday, June 3, 2011 5:56 PM
Answers

John Guiver replied on 07072010 12:27 PM
Yes sorry  I see that now. Somehow when I copied and pasted into VisualStudio, it neglected to paste the lightbulbs.
As regards VMP versus EP. You could probably do an EP version by explicitly coding the matrix multiply as EP does support product of two Gaussian random variables. But it needs to do an expensive quadrature to calculate the messages. Let me know how you get along if you decide to do this.
John
 Marked as answer by Microsoft Research Friday, June 3, 2011 5:56 PM
Friday, June 3, 2011 5:56 PM
All replies

John Guiver replied on 07012010 5:18 AM
Hi Thomas
The answer to (1) depends on what inferences you are trying to make. With fixed A,B,C,D, you can make inferences about hidden states. Also you can make predictions about future states and outputs. You should also be able to make inferences about system and observation noise. If you structure it correctly, you also may be able to make inferences about individual matrix entries, though you would need to use the MatrixMultiply factor (or explicitly implement the matrix multiplication) which uses 2D variable arrays over type double, rather than a random variable over a matrix type; you will need to use variational message passing for those inferences.
The problem lines in your code are due to the fact that you are missing the indices on vState and vObservations. You should have:
vState[i] = Variable.VectorGaussianFromMeanAndPrecision(
Variable.MatrixTimesVector(vA, vState[i1]), mB);
vObservations[i] = Variable.VectorGaussianFromMeanAndPrecision(vState[i], mD);
You will also need to use expectation propagation for your current implementation which uses the MatrixTimesVector factor, as VMP does not support this factor. If you want to use VMP, you should use the MatrixMultiply factor instead.In general, you will should think carefully about how you structure your model. For example, you should think about an 'online' single timeslice version where the priors for state vector come from posteriors from the previous state vector. On the other hand, you may also have scenarios where you want to update the inferred state vector for a whole bunch of previous states  for example when observations are missing or delayed.
John
Friday, June 3, 2011 5:56 PM 
tommy1802 replied on 07012010 8:28 AM
Hi John,
Thank you for your reply.
John Guiver:The answer to (1) depends on what inferences you are trying to makeFor now, I would like to implement a Kalmanfilter. That means I would like to make inference about the last state in the chain given all the observations. But my future plans are to learn some of the Parameters A,B,C,D by defining prios over them. Like in "Variational Bayes Linear Dynamical Systems" (Beal 2003) .
John Guiver:The problem lines in your code are due to the fact that you are missing the indices on vState and vObservations. You should have:
vState[i] = Variable.VectorGaussianFromMeanAndPrecision(
Variable.MatrixTimesVector(vA, vState[i1]), mB);
vObservations[i] = Variable.VectorGaussianFromMeanAndPrecision(vState[i], mD);I have implemented the code exactly like you suggestet. In my code
= [i] : it automatically made a bulb out of it > compare my PSsection about beautiprint of code.
But you are right the problem was the inference engine. With "Expectation Propagation" it works.
So you mean, if I would like to infer some of the parameters, I must change the code to use variational message passing?
Thank you,
Thomas
Friday, June 3, 2011 5:56 PM 
John Guiver replied on 07072010 12:27 PM
Yes sorry  I see that now. Somehow when I copied and pasted into VisualStudio, it neglected to paste the lightbulbs.
As regards VMP versus EP. You could probably do an EP version by explicitly coding the matrix multiply as EP does support product of two Gaussian random variables. But it needs to do an expensive quadrature to calculate the messages. Let me know how you get along if you decide to do this.
John
 Marked as answer by Microsoft Research Friday, June 3, 2011 5:56 PM
Friday, June 3, 2011 5:56 PM