I am currently preparing some ‘demo’ teaching sessions around a few subjects:

• Multithreaded & lock-less C++
• (Kalman) filtering
• Trading strategy backtesting, the technical part, i.e. /not/ the strategies generation process per se
• likely, a few other subjects… later.

While looking at Kalman filter litterature, I had an extra look at so-called tuning methods.

Let us put some background.

In its linear, discrete time form, Kalman filters can be described by the following equations

$$\begin{array}{lcl} x_{t+1} &=& A x_t + w \\ y &=& C x_t + v \\ \left( \begin{array}{c} w \\ v \end{array} \right) & \sim & N \left( 0, \left( \begin{array}{cc} Q_w & 0 \\ 0 & R_v \end{array} \right) \right) \end{array}$$

Now, the best case is when covariance matrices $ Q_w $ & $ R_v $ are known. In my last real-world usage of filters (spot-vol dynamics of EQ implied vol surfaces, if you ask), due to the way I used data, this was not a problem as everything is ‘observable’ eventually (as in: can be measured, not the matrix rank observability).

But what do you do if you cannot really observe empirical errors? One scenario I can imagine is when you are using some sensor and, well, cross checking would require… uhh, the same sensor and it all depends on some physical calibration of your embedded system and this is just a pain, at the very last…

One possible method consists in estimating a maximum likelihood. Without much ado, it can be shown we end up with a matrix optimization problem.

$$\begin{array}{l} \underset{R_w, Q_w}{min}{ ln \, det P_{R_w,Q_w} + Y' P^{-1}_{R_w, Q_w} Y } \\ s.t. Q_w, R_w \geq 0 \\ P_{R_w,Q_w} = \overset{N+K-1}{\underset{i=1}{\Sigma}} \mathbb{O}_i Q_w \mathbb{O}'_i + \overset{N}{\underset{j=1}{\Sigma}}\mathbb{I}_j R_v \mathbb{I}_j' \end{array}$$

One solution is to decompose $ P $, Cholesky or other. Note that $ P $ is quite sparse. To this end, I want to suggest, at least as a demo, some randomized matrix algorithm to find the Eigen decomposition. The little experience I have with generalized sparse matrix packages make it unsuitable for some embedded applications: too much code, and so, how do you embed that much code in your embedded solution? Who will audit the code?

So let us expose this algorithm.

We are dealing here with a ‘Follow the Perturbed Leader’, oracle-based eigenvector in an online setting.

We assume we have a distribution $ \mathcal{D} $ over $ \mathbb{S}_n $, the set of symmetrix matrices. More on a specific example /which can be evaluated ‘sparsely’ /.

We write $ A $, which is the matrix we want the eigenvectors of, as a sum of (possibly more sparse, possibly random) matrices $ A_\tau $. Kindly note that rescaling by a constant doesn’t change the result.

Lastly, we assume we have an ‘oracle’ $ EV(A) $ which will compute an ‘approximate’ leading eigenvector. Concretely, this is one variant of an iterative method such as $ Power $ or $ Lanczos $.

$$\mathbf{Algorithm}: Follow \, the \, Perturbed \, Leader \, (FPL)\\ \begin{array}{ll} 1: & Sample \, a \, matrix \, N \sim \mathcal{D}. \\ 2: & x_1 \leftarrow EV(n). \\ 3: & \mathbf{for} \, t= 1,2, \ldots \, \mathbf{do} \\ 4: & \quad Play \, x_t \, and \, observe \, A_t \\ 5: & \quad x_{t+1} \leftarrow EV \left( \Sigma_{\tau=1}^t A_\tau + N \right) \\ 6: & \mathbf{end \, for} \end{array}$$

One way of evaluating $ N \sim \mathcal{D} $ is to draw a vector $ v $ of i.i.d normals and compute $ N = c . vv’ $, with $c$ some constant. (this happen to be of the form $c = \sqrt{ \frac{T}{n} max \left\{ 1, ln \frac{T}{n} \right\} }$ ).

Assuming we wish to compute the approximate eigenvectors, values of $A$. On possible overall algorithm can be

$$\begin{array}{ll} 1: & \mathbf{for} \, i = 1, \ldots \, N \, \mathbf{do} \\ 2: & \quad Guess \, approximate \, 'leading' \, eigenvector \, l_i and \, eigenvalue \, \lambda_i \, of \, A \, via \, FPL \, algorithm \\ & \quad( Power \, method \, being \, completed \, with \, orthonormalization \, steps ) \\ 3: & \quad A \leftarrow A - \lambda_i l_i l_i' \\ 4: & \mathbf{end \, for} \end{array}$$

References

Identification of Disturbance Covariances Using Maximum Likelihood Estimation, Megan Zagrobelny &al, 2014

Online Learning of Eigenvectors, Dan Garber & al, JMLR 2015