`smooth` throws a convert `MethodError` with PDMats ≥ 0.11.40; the version cap hides a type invariant that isn't enforced

[yebai]

GeneralisedFilters.smooth is unusable for linear Gaussian models with recent PDMats. filter on the
same model is fine. Released v0.4.2 has no PDMats upper bound, so this is what a user gets from the
registry today; main avoids it only through the PDMats = "0.11.35 - 0.11.39" cap added in #182.

The underlying cause is that PDMat(Symmetric(Σ)) and PDMat(Σ) are no longer the same type, so the
Kalman stages stopped agreeing on their state type. I think the cap is worth replacing with a fix,
both to unblock newer PDMats and because the invariant is currently only maintained by accident.

Reproducer

With GeneralisedFilters v0.4.2 and PDMats v0.11.41:

using GeneralisedFilters, PDMats, Random

model = create_homogeneous_linear_gaussian_model(
    [0.0],                          # μ0
    PDMat(reshape([1.7], 1, 1)),    # Σ0
    reshape([0.8], 1, 1),           # A
    [0.0],                          # b
    PDMat(reshape([0.5], 1, 1)),    # Q
    reshape([1.0], 1, 1),           # H
    [0.0],                          # c
    PDMat(reshape([0.3], 1, 1)),    # R
)
ys = [[0.4], [-0.7], [1.1], [0.2]]

_, ll = GeneralisedFilters.filter(model, KalmanFilter(), ys)   # ok, log p(y) = -5.565380648337504
GeneralisedFilters.smooth(Random.default_rng(), model, KalmanSmoother(), ys)

ERROR: MethodError: Cannot `convert` an object of type
  MvNormal{Float64,PDMat{Float64,Symmetric{Float64, Matrix{Float64}},Cholesky{…}},Vector{Float64}} to an object of type
  MvNormal{Float64,PDMat{Float64,Matrix{Float64},Cholesky{…}},Vector{Float64}}
Stacktrace:
 [1] setindex!(A::Vector{FullNormal{…PDMat{Float64, Matrix{Float64}…}}}, x::FullNormal{…PDMat{Float64, Symmetric{…}…}}, i::Int64)
 [2] (::StateCallback)(…, ::PostUpdateCallback; …)      # kalman.jl, filtered_states[iter] = deepcopy(state)
 [8] filter(…; callback::StateCallback, …)
[10] smooth(…)

It fails on the first step, and for any prior covariance type — swapping PDMat(Σ0) for
PDMat(Symmetric(Σ0)) just moves the failure from filtered_states to proposed_states.

Cause

Between PDMats 0.11.39 and 0.11.41, X_A_Xt began returning a Symmetric, and PDMat began keeping
a Symmetric argument in its type parameters instead of collapsing it to the parent array type. (The
cap's upper bound suggests 0.11.40 is the change; I tested 0.11.39 and 0.11.41 only.)

PDMats	PDMat(M)	PDMat(Symmetric(M))
0.11.39	PDMat{Float64,Matrix}	PDMat{Float64,Matrix} — same type
0.11.41	PDMat{Float64,Matrix,Cholesky}	PDMat{Float64,Symmetric{…},Cholesky}

The two Kalman stages then diverge, for a reason that is easy to miss:

• update forms X_A_Xt(Σ, I - K*H) + X_A_Xt(R, K), a sum of two Symmetrics, which stays
  Symmetric; PDMat(Symmetric(·)) now preserves that, giving PDMat{…,Symmetric{…}}.
• predict forms X_A_Xt(Σ, A) + Q, and because Q is a PDMat the sum falls back to a plain
  matrix, giving PDMat{…,Matrix}.

So a predicted state and a filtered state have different types. StateCallback sizes both of its
caches from the initial state:

callback.proposed_states = Vector{T}(undef, N)   # T = typeof(initial state)
callback.filtered_states = Vector{T}(undef, N)

and no single concrete T can hold both. filter never trips this because it stores nothing — the
callback is the only thing that keeps states, which is why only smooth breaks.

The callback isn't really at fault: it relies on the stages agreeing on a state type, which they did
until PDMats changed under it.

The cap works, but it has costs

PDMats = "0.11.35 - 0.11.39" on main does prevent the crash, and CI presumably drove it: the
existing "Kalman filter StaticArrays" test fails under 0.11.41 on main with

Expression: state.Σ isa PDMat{Float64, SMatrix{D, D, Float64, D * D}}
 Evaluated: … isa PDMat{Float64, SMatrix{2, 2, Float64, 4}, C} where C<:Cholesky{Float64}

Two costs worth weighing. The cap blocks co-installation with anything that wants newer PDMats, which
is a real constraint for a package meant to be used alongside others. And the invariant that made the
code correct is still unstated and unenforced, so the next equivalent upstream change lands the same
way.

Worth noting separately that smooth itself has no direct test on the linear Gaussian path — the
"Kalman smoother" tests hand-roll the backward pass with backward_smooth in a fold, so the only
thing that catches this is the StaticArrays type assertion.

Suggested fix

Make the stages agree at the source, so Vector{T} becomes correct as written and stays concretely
typed. The canonical form is a plain array, which is what initialise and predict already produce,
so predict needs no change:

"""
    pd_symmetric(Σ) -> PDMat

Wrap a covariance as a `PDMat` whose stored matrix is a plain array of the same kind the arithmetic
produced, averaging the two triangles and unwrapping a `Symmetric` rather than keeping it. …
"""
pd_symmetric(Σ::AbstractMatrix) = PDMats.PDMat((Σ + Σ') / 2)
pd_symmetric(Σ::Symmetric) = PDMats.PDMat(identity.(Σ))

then replace the six PDMat(Symmetric(·)) sites with it — four in algorithms/kalman.jl
(_apply_jitter_and_wrap, backward_smooth, backward_initialise, backward_predict) and two in
GFTest (make_pd, and its copy in GFTest/models/linear_gaussian.jl, which otherwise makes the
test models' own priors Symmetric-backed). That is +30/−7 including the docstring, and afterwards no
Symmetric( remains in src outside the helper, so the invariant lives in one place.

Two details behind the choices: unwrapping via broadcasting keeps a statically sized input statically
sized, where Matrix(Σ) would densify the SMatrix models; and averaging the triangles matches the
idiom already used in GFTest.make_pd.

I have a patch ready if it's useful.

Verification

Running the existing Kalman testitems in both PDMats versions, with the cap relaxed to "0.11.35":

	PDMats 0.11.39	PDMats 0.11.41
unfixed	8/8 pass	7/8 — StaticArrays fails
fixed	8/8 pass	8/8 pass

Also under 0.11.41 with the fix: the four algorithms/discrete.jl testitems pass (these do call
smooth), and smooth end to end matches joint-Gaussian conditioning on the model above to 1e-12
(p(x₁|y₁:₄) mean 0.20523634768, variance 0.207044369872, log p(y) −5.565380648338).

What I did not run: algorithms/csmc.jl, the JET type-stability items, Aqua, and the
Turing/Mooncake/CUDA integrations — I could not resolve those test deps locally, so
integrations/kalman_rrule.jl in particular is unverified, and a custom AD rule is where a covariance
type change seems most likely to bite.

Environment

Julia 1.12.6, macOS (aarch64). GeneralisedFilters v0.4.2 from the registry with PDMats v0.11.41 for
the reproducer; main (v0.5.0, at #182) with the cap relaxed for the fix and the test runs.