locked
Linear Dynamical System (Migrated from community.research.microsoft.com) RRS feed

  • Question

  • tommy1802 posted on 06-29-2010 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[i-1]), 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 2-D 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 07-07-2010 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

    Friday, June 3, 2011 5:56 PM

All replies

  • John Guiver replied on 07-01-2010 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 2-D 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[i-1]), 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 time-slice 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 07-01-2010 8:28 AM

    Hi John,

    Thank you for your reply.

    John Guiver:
    The answer to (1) depends on what inferences you are trying to make

    For 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[i-1]), 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 PS-section 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 07-07-2010 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

    Friday, June 3, 2011 5:56 PM