Yes, makes sense. ContinuousTransition node can indeed be helpful as long as you can frame your transition as $$x_t=Ax_{t-1}$$.

By the way, how does Linearization perform in your case?
Another question: doesn't the inference run without an initialization parameter in your case?

Thanks for your question!

It's not entirely fair to compare with the Kalman Filter since your current approach is smoothing on a full graph with 12,000 observations, as outlined by @albertpod. Each iteration requires a full forward and backward pass, involving many approximations. Specifically, it needs to call the multivariate Unscented transform 48,000 times per iteration.

If you want to adapt this to the filtering case (using only the forward pass), I recommend this material in addition to the one suggested by Albert. Filtering should indeed improve performance.

Additionally, this hidden example demonstrates the significant performance difference between filtering and smoothing. For a small dataset of univariate (!) 500 observations with no approximations, filtering shows a 30x performance improvement. Let us know if you need further assistance in adapting your example to the filtering case.

Another suggestion: you might use the dot function with a standard basis vector as the observation function. This approach can eliminate the need for the g function and the Unscented transform, improving accuracy and convergence speed, and removing the need for initializing messages and performing iterations. We have specialized implementations for dot that compute messages faster and more accurately, for example:

@model function ssm(y, f, g, tS, x0, u, t, Q, R, p)
    x_prior ~ x0
    c = ReactiveMP.StandardBasisVector(2, 1)
    x_prev = x_prior
    for i in eachindex(y)
        x[i] ~ MvNormal(μ=f(x_prev, u[:,i], p, t[i], tS), Σ=Q)
        y[i] ~ Normal(μ=dot(c, x[i]), var=R[1, 1])
        x_prev = x[i]
    end
end

prior_x = MvNormalMeanCovariance(x[1], diageye(2))

model = ssm(f=f, g=g, tS=tS, x0=prior_x, u=u_sim, t=t_sim, Q=Q, R=R, p = p)

ssm_meta = @meta begin 
    f() -> Unscented()
end

Additionally, avoid using the global variable p; instead, pass it as an argument to the model function (and any other variables I may have missed). Otherwise, Julia gets confused and checks its type every time it is accessed, which significantly slows down performance.

Also, consider checking Linearization instead of Unscented for potentially faster performance. (beware global variables, it makes it even slower when using Linearization because it utilizes ForwardDiff, which performs poorly with global variable access.)

By replacing g with dot and making p local to the model function, I was able to reduce inference times on my computer from 150 seconds per iteration to about 50 seconds. Not that bad if you consider full forward and full backward pass with 12,000 observations.

@time result = infer(
    model = model,
    data  = (y = first.(y), ), 
    options = (limit_stack_depth = 200, ), 
    # initialization = imessages,
    meta = ssm_meta,
    iterations = 1,
    showprogress = true,
)

The inference also required only one iteration and converged to the following:

result.posteriors[:x][1] # why is this a 1-element vector?

it depends on the number of iterations. Perhaps you ran it with a single iteration?

options = (limit_stack_depth = 10, ),

This setting is very low and can significantly slow down inference. It's recommended to set it as high as your computer allows; I typically use 500-600. However, in this example, I found 250 to be sufficient. This setting is usually not required for the filtering case.

How could I infer p.m?

To infer m, you can treat it as a standalone random variable, assign a prior to it, and infer it like any other variable. We do not support inference over components of a named tuple, so it should be a separate variable.

Please let us know if you have other question. I or my colleagues will be happy to help you.

I would also emphasize that while filtering is generally faster, it is often less accurate. The difference in accuracy should be evaluated on a case-by-case basis.

Additionally, as @albertpod mentioned, our current Unscented approximation technique requires some refactoring and is quite slow. We would be happy to assist if you choose RxInfer and decide to work on improving it. We are willing to guide you through the process.

Great, many thanks @Nimrais ! How long does it take for you if you don't compute the free energy? It's generally not necessary for the inference, only for checking if the algorithm has converged or not.

I would also assume that the R (which is the observational noise) has quite a big value. It was the same in the original formulation from @ValentinKaisermayer, but I can see that the scale of the data is somewhat (-20, 20) and the observational noise is set to 10. That might also explain the difference in the inferred mean. Does it infer better if you decrease R from 10 to let's say 1 or even 0.1?

@bvdmitri Yes, by reducing the R value, I can obtain improvement for the state estimation part, which makes sense as the noise is quite strong in the original formulation.

Here, I reduced the R to 1

However, the free energy does not look smooth. I think it's the approximation happening around the m edge that does not really perform well.

Regarding the free energy computation overhead, it does not seem to really contribute to the computation time: around 1.13 sec per iteration with the FE computation, and around 0.9 sec per iteration without.

@albertpod In the model with the inference for p using the Linearization() meta, the inference is diverging.

Hm, interesting. Do you measure the compilation time? @time should report these, and while compilation can take a significant portion of the runtime, it should typically occur only once. However, it might happen every time if Julia is confused for some reason. If @time doesn't show excessive compilation times, the issue might be with the model creation, although it shouldn't take that long. You could use callbacks to measure the exact time of each step in the procedure. For example you can use before_model_creation and after_model_creation callbacks to measure time between those. This would help us understand why it's taking so much time, even though the progress bar only reports 8 seconds.