网格数据上的高效卡尔曼滤波器实现

Posted 2023-02-21

【中文标题】网格数据上的高效卡尔曼滤波器实现

【英文标题】：Efficient Kalman filter implementation on gridded data

【发布时间】：2013-09-17 16:45:28

【问题描述】：

我编写了一个非常简单的卡尔曼滤波器，它在时间序列上运行（包括数据间隙）。它工作得很好，但我碰巧有一个 数据立方体 数据（一个形状为 Nt, Ny, Nx 的数组），我想为数据立方体中的每个像素应用我的时间卡尔曼滤波器。我已经完成了明显的（在最后两个维度上循环），但这需要相当长的时间。

最终，我总是不得不从单个“像素”中提取数据，并构建相关的矩阵/向量，因此这个过程非常缓慢（请注意，每个单独的时间序列中的间隙是不同的，通常，所以是将状态与观察联系起来的 H 矩阵）。我不熟悉 cython，它可能会有所帮助（只是我不熟悉它）。

我只是想知道对问题的巧妙改写或巧妙的数据结构是否可以更有效地进行时间过滤。我宁愿这只使用 numpy/scipy，而不是 OpenCV，否则它会很麻烦，这取决于额外的包。

【问题讨论】：

“显而易见的”是不要循环。改用向量化的 NumPy 操作。

@larsmans 由于问题的性质，我看不出如何进行矢量化操作来解决这个问题（您需要构建矩阵等）

然后你需要先对程序进行概要分析以找出什么是昂贵的（Kern's line_profiler 非常适合），然后发布缓慢的代码部分并询问如何将它们矢量化。就目前而言，这个问题无法正确回答。

您能否发布您目前获得的代码，以便我们了解如何提供帮助？

【参考方案1】：

我制作了一个像这样的简单矢量化卡尔曼滤波器来处理电影帧。它非常快，但目前仅限于一维输入和输出，并且它不对任何滤波器参数进行 EM 优化。

import numpy as np

def runkalman(y, RQratio=10., meanwindow=10):
    """
    A simple vectorised 1D Kalman smoother

    y       Input array. Smoothing is always applied over the first 
            dimension

    RQratio     An estimate of the ratio of the variance of the output 
            to the variance of the state. A higher RQ ratio will 
            result in more smoothing.

    meanwindow  The initial mean and variance of the output are 
            estimated over this number of timepoints from the start 
            of the array

    References:
        Ghahramani, Z., & Hinton, G. E. (1996). Parameter Estimation for
        Linear Dynamical Systems.
        http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.55.5997

        Yu, B., Shenoy, K., & Sahani, M. (2004). Derivation of Extented
        Kalman Filtering and Smoothing Equations.
        http://www-npl.stanford.edu/~byronyu/papers/derive_eks.pdf
    """

    x,vpre,vpost = forwardfilter(y, RQratio=RQratio, meanwindow=meanwindow)
    x,v = backpass(x, vpre, vpost)
    x[np.isnan(x)] = 0.
    return x

def forwardfilter(y, RQratio=10., meanwindow=10):
    """
    the Kalman forward (filter) pass:

        xpost,Vpre,Vpost = forwardfilter(y)
    """

    y = np.array(y,copy=False,subok=True,dtype=np.float32)

    xpre = np.empty_like(y)
    xpost = np.empty_like(y)
    Vpre = np.empty_like(y)
    Vpost = np.empty_like(y)
    K = np.empty_like(y)

    # initial conditions
    pi0 = y[:meanwindow].mean(0)
    ystd = np.std(y,0)
    R = ystd * ystd
    Q = R / RQratio
    V0 = Q

    xpre[0] = xpost[0] = pi0
    Vpre[0] = Vpost[0] = V0
    K[0] = 0

    # loop forwards through time
    for tt in xrange(1, y.shape[0]):

        xpre[tt] = xpost[tt-1]
        Vpre[tt] = Vpost[tt-1] + Q

        K[tt] = Vpre[tt] / (Vpre[tt] + R)
        xpost[tt] = xpre[tt] + K[tt] * (y[tt] - xpre[tt])
        Vpost[tt] = Vpre[tt] - K[tt] * (Vpre[tt])

    return xpost,Vpre,Vpost

def backpass(x, Vpre, V):
    """ 
    the Kalman backward (smoothing) pass:

        xpost,Vpost = backpass(x,Vpre,V)
    """

    xpost = np.empty_like(x)
    Vpost = np.empty_like(x)
    J = np.empty_like(x)

    xpost[-1] = x[-1]
    Vpost[-1] = V[-1]

    # loop backwards through time
    for tt in xrange(x.shape[0]-1, 0, -1):
        J[tt-1] = V[tt-1] / Vpre[tt]
        xpost[tt-1] = x[tt-1] + J[tt-1] * (xpost[tt] - x[tt-1])
        Vpost[tt-1] = V[tt-1] + J[tt-1] * (Vpost[tt] - Vpre[tt]) * J[tt-1]

    return xpost,Vpost 

如果有人知道 Python 中支持多维输入/输出和 EM 参数优化的矢量化卡尔曼平滑器实现，我很想听听！我向 PyKalman 维护者提出了功能请求，但他们说清晰度比速度更重要。

【讨论】：

如果您有这样一个“清晰”算法循环数组索引并希望它快速运行，请查看 fortran 和 f2py（它是 scipy/numpy 宇宙的一部分）。该代码与您的 python 代码基本相同。我曾经为我的伙伴比较过不同的加速方法，其中也包括 fortran。 nbviewer.ipython.org/5774959

@P.R.我不确定您的意思-您链接到的代码与卡尔曼滤波器无关

否，但它展示了如何将循环的 numpy 代码转换为运行速度非常快的 Fortran 代码，同时不会产生奇怪的依赖关系。如果您不关心依赖关系，并且希望代码在几乎没有（实际上没有）重新编码的情况下更快，请查看该笔记本的 numba 部分。

@P.R.一维卡尔曼滤波器已经涉及大量的向量化计算（矩阵乘法等），因此并不特别适合 LLVM。在 Fortran 中重写部分或全部内容还需要大量工作。老实说，我认为在这种情况下最好的策略是仔细考虑每个中间变量的维度，并使用 Numpy 的广播规则将计算应用于输入数组的任何其他维度。

@ali_m 非常正确，事实证明我的问题实际上很容易使用广播来解决，因为一些矩阵有效地评估为标量（1x1 矩阵），大大简化了事情！