The most computationally intensive part of my application is the triangularisation of a matrix, $\mathbf{S} = \operatorname{triag}\left(\mathbf{A}\right)$, such that $\mathbf{S}\mathbf{S}^T = \mathbf{A} \mathbf{A}^T$.

For this application $\mathbf{A} \in \Re^{n \times 2(n+v)}$, where $n$ is generally close to the same value as $v$.

Currentl, I solve via a QR decomposition, $\mathbf{A}^T = \mathbf{Q} \mathbf{R}$:

$\mathbf{A} \mathbf{A}^T = (\mathbf{R}^T \mathbf{Q}^T) (\mathbf{Q} \mathbf{R})$

$ = \mathbf{R}^T \mathbf{Q}^T \mathbf{Q} \mathbf{R} $

$ = \mathbf{R}^T \mathbf{R}$

$ = \mathbf{S} \mathbf{S}^T$

Which gives the form that I want.

Is there a more efficient way to perform this operation? Even a modest saving over QR will provide a performance boost as this one operation is profiled as 30% of my total computing time.

The $\operatorname{triag}()$ function is the heart of the Divided Difference Filter, which is similar to the more commonly known Unscented Kalman Filter. The problem is mentioned explicitly in the paper that introduces the Cubature Kalman Filter.

In short, it calculates four "divided difference" matrices, based on the non-linear process model and state covariance Cholesky factors. These four matrices are horizontally concatenated, and then triangularised - this forms the predicted state covariance:

$ \mathbf{S}(k+1|k) = \operatorname{triag}\left(\begin{bmatrix} \mathbf{S}_{xx}^{(2)} & \mathbf{S}_{xv}^{(1)} & \mathbf{S}_{xx}^{(2)} & \mathbf{S}_{xv}^{(2)} \end{bmatrix} \right)$

where $\mathbf{S}_{pq}^n$ above are the "divided difference" matrices calculated above. The predicted state covariance Cholesy factor $\mathbf{S}(k+1|k)$ is then directly used to generate "divided difference" matrices of the predicted measurements.

A similar process follows of triangularisation of the predicted measurements, but as my measurement vector is small compared to my state vector, this calculation (and the subsequent Kalman Gain/Update) equations are not computationally intensive (about 4% of the process time).

  • $\begingroup$ What is your ultimate speed goal? Suppose you could cut the QR time in half, then you'd only be saving 15% CPU time. Is that really enough to achieve what you hope to achieve? $\endgroup$ May 7 '13 at 15:07
  • $\begingroup$ @MichaelC.Grant, the algorithm runs on a real-time embedded system where the filter must compete with plenty of other processes. There is no option to "throw more hardware at the problem" as there is a substantial install base. Any real CPU saving (even if only a couple of percent) gives us options, whether it be for more advanced algorithms or more features elsewhere in the system. Or, we'll simply claim the CPU for a latency reduction, which is a worthy goal in and of itself. $\endgroup$
    – Damien
    May 7 '13 at 23:07
  • $\begingroup$ Fair enough. Then I'm going to cast my lot with Brian's answer below. If you think $AA^T$ will be reasonably well conditioned, then it might save you some cycles to compute its Cholesky factor. $\endgroup$ May 10 '13 at 18:49

If $AA^{T}$ is well conditioned, then you might be able to get away with using the Cholesky factorization of $AA^{T}$ rather than using the QR factorization. Depending on the size of $A$ this could be vastly faster than computing the QR factorization of $A$. If you told us something about the size of $A$ we could address that question.

However, if $AA^{T}$ is badly conditioned, then using the Cholesky factorization could be a bad idea because the Cholesky factorization will be less accurate. It would help to know something about the condition number of $AA^{T}$.

You haven't told us why you want $S$ and what you do with it after computing this factorization. It's quite possible that there's a much more efficient way to solve the higher level problem without ever computing this factorization of $AA^{T}$. However, we can't make any suggestions unless you provide more background on what you're trying to accomplish. In other words, a good answer to this question is going to depend on a lot of context that you haven't provided yet.

  • $\begingroup$ I'll update my question shortly. $\endgroup$
    – Damien
    May 5 '13 at 4:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.