0
$\begingroup$

I am trying to create a compiled version of this function:

cKalmanTimeLoop = Compile[{{Y, _Real, 2}, {X, _Real, 2}, {a, _Real, 2}, {b1, _Real, 
2}, {Ve, _Real}, {Vw, _Real, 2}, {Covb1, _Real, 2}},
Module[{Yhat, b, , k, n, e, t, Q, K, Covb, Zscore
(************************************************************)
(************************************************************)
(**UPDATES THE KALMAN FILTER PARAMETERS FOR EACH TIMESTEP**)
(************************************************************)
(*X=k x n matrix of k prices data series for t=1,...,n;
note that X(t)[[1,All]]=1 for all t;
a=k x k State Transition matrix estimated in initialization;
b=k x n matrix of state vectors b(t)={b1,...,bk},Ve=
scalar price error variance,estimated in initialization;
Vw=k x k state noise covariance matrix,estimated in initialization;
Covb[[t]]=k x k state covariance matrix at time t=1,...,n;*)
(***********************************************************)
{k, n} = Dimensions@X;
K = 0.;
e = ConstantArray[0., n];
Q = ConstantArray[0., n];
b = ConstantArray[0., {k, n}];
Covb = ConstantArray[0., {n, k, k}];
b[[All, 1]] = b1;
Covb[[1]] = Covb1;
Do[
(*************************************************)
(*update state prediction and state covariance matrix*)
If[t > 1, b[[All, t]] = a.b[[All, t - 1]];
 Covb[[t]] = a.Covb[[t - 1]].Transpose[a] + Vw;
 Covb[[t]] = DiagonalMatrix@Diagonal@Covb[[t]];];
(*************************************************)
(*Measurement Prediction Equation*)
Yhat = X[[All, t]].b[[All, t]];
(*************************************************)
(*Measurement Prediction Error*)
e[[t]] = Y[[t]] - Yhat;
(*************************************************)
(*Measurement Covariance Prediction*)
Q[[t]] = X[[All, t]].Covb[[t]].X[[All, t]] + Ve;
(*************************************************)
(*Kalman Gain*)
K = Covb[[t]].X[[All, t]]/Q[[t]];
(*************************************************)
(*State Update*)
b[[All, t]] = b[[All, t]] + K*e[[t]];
(*************************************************)
(*State Covariance Update*)
Covb[[t]] = Covb[[t]] - DiagonalMatrix[K*X[[All, t]]].Covb[[t]];
, {t, 1, n}];

 Zscore = e/Sqrt[Q];
{b, Covb, K, e, Q, Zscore}], {{Covb, _Real, 3}}, 
Parallelization -> True, CompilationTarget -> "C"]

But I cant seem to get past the following error message:

Compile::cplist: Compile`$28.Covb[[t]] should be a tensor of type Integer, Real, or Complex; evaluation will use the uncompiled function.

I assume that I am making some elementary syntax error, but I can't spot it.

Also, any suggestions for other possible speed-ups to the code would be gratefully received.

$\endgroup$
1
  • 1
    $\begingroup$ I debugged this code. How to check (test)? $\endgroup$ Dec 14, 2018 at 18:08

1 Answer 1

4
$\begingroup$

Tha Kalman gain is a list, not a scalar here. So, try this (note that I put all variables in lower case, so minor changes):

cKalmanTimeLoop = Compile[
    {
       {y, _Real, 2}, 
       {x, _Real, 2}, 
       {a, _Real, 2}, 
       {b1, _Real,2}, 
       {ve, _Real}, 
       {vw, _Real, 2}, 
       {covb1, _Real, 2}
    },
    Module[
        { yhat, b, k, n, e, t, q, capk, covb, zscore},
        (************************************************************)
        (************************************************************)
        (**UPDATES THE KALMAN FILTER PARAMETERS FOR EACH TIMESTEP**)
        (************************************************************)
        (*X=k x n matrix of k prices data series for t=1,...,n;
          note that X(t)[[1,All]]=1 for all t;
          a=k x k State Transition matrix estimated in initialization;
          b=k x n matrix of state vectors b(t)={b1,...,bk},
          Ve= scalar price error variance,estimated in initialization;
          Vw=k x k state noise covariance matrix,estimated in initialization;
          Covb[[t]]=k x k state covariance matrix at time t=1,...,n;*)
        (***********************************************************) 
        {k, n} = Dimensions@x;
        capk = Table[0., {n}]; (* needs to be a vector of length n *)
        e = Table[0., {n}];
        q = Table[0., {n}];
        b = Table[0., {k}, {n}];
        covb = Table[0., {n}, {k}, {k}];
        b[[All, 1]] = b1;
        covb[[1]] = covb1;
        Do[
            (*************************************************)
            (*update state prediction and state covariance matrix*)
            If[ t > 1,
                b[[All, t]] = a.b[[All, t - 1]];
                covb[[t]] = a.covb[[t - 1]].Transpose[a] + vw;
                covb[[t]] = DiagonalMatrix@Diagonal@covb[[t]];
            ];
            (*************************************************)
            (*Measurement Prediction Equation*)
            yhat = x[[All, t]].b[[All, t]];
            (*************************************************)
            (*Measurement Prediction Error*)
            e[[t]] = y[[t]] - yhat;
            (*************************************************)
            (*Measurement Covariance Prediction*)
            q[[t]] = x[[All, t]].covb[[t]].x[[All, t]] + ve;
            (*************************************************)
            (*Kalman Gain*)
            (* NOTE: THIS GIVES A LIST - NOT A SCALAR!!! *)
            capk = covb[[t]].x[[All, t]]/q[[t]];
            (*************************************************)
            (*State Update*)
            b[[All, t]] = b[[All, t]] + capk*e[[t]];
            (*************************************************)
            (*State Covariance Update*)
            covb[[t]] = covb[[t]] - DiagonalMatrix[capk*x[[All, t]]].covb[[t]];
            ,
            {t, 1, n}
       ];
       zscore = e/Sqrt[q];
       {b, covb, capk, e, q, zscore}
   ],
   {{covb,_Real,3}},
   Parallelization -> True, 
   CompilationTarget -> "C"
]

You can check with CompilePrint which needs the Package CompiledFunctionTools.

$\endgroup$
3
  • 2
    $\begingroup$ Note that there is still code that calls the main-evaluator. I'd have a look at Diagonal and DiagonalMatrix. These two seem to be one of the reasons and they can be replaced by injected functions. $\endgroup$
    – halirutan
    Dec 14, 2018 at 20:17
  • $\begingroup$ Thanks gwr. This works well. Having said that, it isn't any faster than the WL code, which comes as a surprise (to me). $\endgroup$ Dec 15, 2018 at 8:43
  • $\begingroup$ @JonathanKinlay There is still MainEvaluator being called as halirutan pointed out - that is the reason for example for me using Table instead of ConstantArray. So there is still some speed gains to be made by optimizing the code. $\endgroup$
    – gwr
    Dec 15, 2018 at 8:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.