# \[ANN\] QuadraticKalman.jl - Quadratic State-Space Filtering

**URL:** <https://discourse.julialang.org/t/ann-quadratickalman-jl-quadratic-state-space-filtering/125538>\
**Category:** Package Announcements\
**Tags:** package\
**Created:** [February 4, 2025, 3:37pm UTC](https://discourse.julialang.org/t/ann-quadratickalman-jl-quadratic-state-space-filtering/125538 "2025-02-04T15:37:35Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![dleather](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dleather/32/43762_2.png) [@dleather](https://discourse.julialang.org/u/dleather)\
**Post date:** [February 4, 2025, 3:37pm UTC](https://discourse.julialang.org/t/ann-quadratickalman-jl-quadratic-state-space-filtering/125538/1 "2025-02-04T15:37:35Z")

</div>

# [ANN] QuadraticKalman.jl: Kalman Filtering for Quadratic State-Space Models

I’m excited to announce the first release (and my first package) of [QuadraticKalman.jl](https://github.com/dleather/QuadraticKalman.jl), a Julia package that implements Kalman filtering and smoothing for state-space models with quadratic measurement equations.

## Overview

QuadraticKalman.jl implements the _Quadratic Kalman Filter_ from [Monfort et. al. (2015, Journal of Econometrics)](https://www.sciencedirect.com/science/article/abs/pii/S0304407615000123) for the following state-space model with Gaussian noise:

State equation:

X\_{t} = \mu + \Phi X\_{t-1} + \Omega \epsilon\_t, \quad \epsilon\_t \sim \mathcal{N}(\mathbf{0},\mathbb{I}\_N)

Measurement Equation:

Y\_t = A + B X\_t + \alpha Y\_{t-1} + \sum\_{i=1}^M e\_i X\_t^\prime C\_i X\_t + D \eta\_t, \quad \eta\_t \sim \mathcal{N}(\mathbf{0}, \mathbb{I}\_M), e\_i := \mathbb{I}\_M[:,i]

The Quadratic Kalman Filter is an improvement over naively adding quadratic terms to the classic Kalman Filter by maintaining the structure of the covariance matrix in the linearized state-space. Although less flexible than the EKF, UKF, and Quadrature Kalman Filter, the QKF is computationally cheaper and more accurate when the measurement equation is known to be quadratic.

## Key Features

- Kalman filtering and smoothing for quadratic state-space models
- Extends original implementation with autoregressive measurement equations
- Automatic differentiation of negative log-likelihood using ForwardDiff.jl
- Visualization tools for filtered and smoothed states
- Efficient parameter-model conversion for optimization
- Automatic reparametrization for positive-definite covariance matrices
- 9x - 62x faster than equivalent R implementation
- Numerically stable implementation

## Future Planned Features

- A mutating version of the filter optimized for memory allocations
- Extend the filter for state-dependent noise
- Compatibility with ReverseDiff.jl

Edit: Changed R comparison after properly benchmarking

---

<div class="post-metadata">

**Author:** ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)\
**Post date:** [February 4, 2025, 7:26pm UTC](https://discourse.julialang.org/t/ann-quadratickalman-jl-quadratic-state-space-filtering/125538/2 "2025-02-04T19:26:52Z")

</div>

Congratulations, this looks very cool. The example docs look great! I like the nice plotting of the uncertainty.

> 8.8x faster than equivalent R implementation

Normally I’d expect R to be 100x slower, because this kind of iterative matrix mul-adds is usually pretty slow. Do you know if that version is calling a C library? The equations look well-suited to Julia, which should then be about same speed as C. Or is there something in the computation that’s bogging down Julia?

---

<div class="post-metadata">

**Author:** ![dleather](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dleather/32/43762_2.png) [@dleather](https://discourse.julialang.org/u/dleather)\
**Post date:** [February 4, 2025, 7:58pm UTC](https://discourse.julialang.org/t/ann-quadratickalman-jl-quadratic-state-space-filtering/125538/3 "2025-02-04T19:58:20Z")

</div>

Thanks @apo383! I appreciate the complements and feedback.

I agree 8.8x is a little slow. For the initial release I optimized the in-place version of the filter and smoother for readability, design, and AD-friendliness. The 8.8x result is for the in-place version. Once I get around to optimizing the mutating version that will likely be much faster. And then I will need to do another iteration to get that code to be AD-friendly which will likely mean trying to get it to play with Enzyme.jl. I tried it once through Enzyme and it caused a crash, but I’ve done some work previously on this and think I can get it to work.

The key bottleneck I faced is a robustness technique to correct for slightly non-PSD covariance update matrices where you replace the negative eigenvalues to 0 and reconstruct the matrix. Currently I maintain AD-friendliness through that routine by using ForwardDiff.jl which can handle `eigen(Symmetric(A))`.

I’d definitely welcome tips on how to write code that is _both_ AD-friendly and limits memory allocation.

---

<div class="post-metadata">

**Author:** ![apo383](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/apo383/32/11272_2.png) [@apo383](https://discourse.julialang.org/u/apo383)\
**Post date:** [February 4, 2025, 10:25pm UTC](https://discourse.julialang.org/t/ann-quadratickalman-jl-quadratic-state-space-filtering/125538/4 "2025-02-04T22:25:07Z")

</div>

That does sound like a likely bottleneck that probably swamps everything else. I briefly looked at the paper you referenced, and the algorithm equations look pretty straightforward, so the eigenvalue fixing seems like the only demanding issue.

I imagine the R code is using BLAS or something, but have no idea how they might be doing AD, which I thought was quite limited in R.

I briefly looked through `QuadraticKalman.jl` and saw a function `d_eigen`, but couldn’t find it used anywhere. Are `d_eigen` and AD are for future use? You may find `DifferentiationInterface.jl` helpful for trying out various approaches.

Probably a question for the original authors, but I wonder why they had to do the eigenvalue trick. A usual way to avoid that would be a square root filter, updating the covariance square root so that the covariance is always PDF. I’m imagining that they couldn’t formulate it in that way.

---

<div class="post-metadata">

**Author:** ![dleather](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dleather/32/43762_2.png) [@dleather](https://discourse.julialang.org/u/dleather)\
**Post date:** [February 5, 2025, 5:08am UTC](https://discourse.julialang.org/t/ann-quadratickalman-jl-quadratic-state-space-filtering/125538/5 "2025-02-05T05:08:13Z")

</div>

So I implemented some proper bench-marking and was underselling myself. For some reason in the 2-dimensional case the speedup is around 10x, but it is around 60x when looking at the scalar case and 25x in the 5-dimensional case.

[![](https://global.discourse-cdn.com/julialang/original/3X/e/f/efe98bedcb9892916c04e3ee46d8b0ccf0302af6.png) ](https://raw.githubusercontent.com/dleather/QuadraticKalman.jl/main/benchmarks/results/scaling_comparison.png)

I think a square-root version of the filter would have a much larger augmented state-space maybe posing a challenge for the R code. It would be interesting to see if Julia could handle it with reasonable speed.

And yes, `d_eig` is an artifact of an older version of the filter I tried to make work with Ezyme.jl.