oseiskar / simdkalman

Python Kalman filters vectorized as Single Instruction, Multiple Data
https://simdkalman.readthedocs.io/
MIT License
176 stars 36 forks source link

Vectorized online update #13

Closed simejanko closed 4 years ago

simejanko commented 4 years ago

I didn't see any indication in the documentation and examples of whether vectorized online update of multiple data streams is possible.

Example:

import numpy as np
import simdkalman

N_STREAMS = 3

kf = simdkalman.KalmanFilter(
    state_transition = np.array([[1,1],[0,1]]),
    process_noise = np.diag([0.1, 0.01]),
    observation_model = np.array([[1,0]]),
    observation_noise = 1.0)

initial_means = np.zeros((N_STREAMS, 2))
initial_covariances = np.tile(np.eye(2), (N_STREAMS, 1, 1)) # shape = (N_STREAMS, 2, 2)
observations = np.random.normal(size=(N_STREAMS, 1))

kf.update(initial_means, initial_covariances, observations)

This throws an assertion error but i'm not sure whether i'm doing something wrong or this kind of use is not supported.

oseiskar commented 4 years ago

Yes, this is possible, but you have to use the primitives module, which works a bit differently: you do not create a KalmanFilter object, but use the primitive functions directly with the previous state and the relevant matrices (A, H, Q, or R). Then they return you the new state. The primitive functions are also a bit more picky about the matrix shapes.

Your example would be modified as follows

import numpy as np
from simdkalman.primitives import update

N_STREAMS = 3

state_transition = np.array([[1,1],[0,1]])
process_noise = np.diag([0.1, 0.01])
observation_model = np.array([[1,0]])
observation_noise = np.array([[1.0]]) # NOTE: wrap as matrix

initial_means = np.zeros((N_STREAMS, 2, 1)) # NOTE: (N, 1) -> (N, 2, 1)
initial_covariances = np.tile(np.eye(2), (N_STREAMS, 1, 1)) # shape = (N_STREAMS, 2, 2)
observations = np.random.normal(size=(N_STREAMS, 1, 1)) # NOTE: (N, 1) -> (N, 1, 1)

means, covariances = update(
    initial_means,
    initial_covariances,
    observation_model,
    observation_noise,
    observations)
simejanko commented 4 years ago

This is exactly what I needed. Thank you.