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$.
Currently, 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).
