Low-rank multi-parametric covariance identification

Abstract

We propose a differential geometric approach for building families of low-rank covariance matrices, via interpolation on low-rank matrix manifolds. In contrast with standard parametric covariance classes, these families offer significant flexibility for problem-specific tailoring via the choice of “anchor” matrices for interpolation, for instance over a grid of relevant conditions describing the underlying stochastic process. The interpolation is computationally tractable in high dimensions, as it only involves manipulations of low-rank matrix factors. We also consider the problem of covariance identification, i.e., selecting the most representative member of the covariance family given a data set. In this setting, standard procedures such as maximum likelihood estimation are nontrivial because the covariance family is rank-deficient; we resolve this issue by casting the identification problem as distance minimization. We demonstrate the utility of these differential geometric families for interpolation and identification in a practical application: wind field covariance approximation for unmanned aerial vehicle navigation.

A Computations required for the variable projection method

We detail here the computations for the main steps of the variable projection method proposed in Sect. 4.2.2. Remember that the corresponding surface, $\varphi ^{\mathrm {LG}}(t_1,t_2) = Y_{\varphi }(t_1,t_2) Y_{\varphi }(t_1,t_2)^\top \!$, where $Y_{\varphi }$ is obtained by composition of two geodesics, respectively along the $t_1$ and $t_2$ variables:

$$\begin{aligned} Y_{\varphi }(t_1, t_2)&{:}{=}(1-t_2) Y_{1-2}(t_1) + t_2 Y_{3-4}(t_1)Q(t_1)^\top \!\,, \\ Y_{1-2}(t_1)&{:}{=}(1-t_1) Y_1 + t_1 Y_2 Q_{1-2}^\top \!\,, \\ Y_{3-4}(t_1)&{:}{=}(1-t_1) Y_3 + t_1 Y_4 Q_{3-4}^\top \!\,. \end{aligned}$$

1.1 A.1 Computation of the partial derivative with respect to $t_1$

Computing the partial derivative of the function $(t_1,t_2) \mapsto \varphi ^{\mathrm {LG}} (t_1,t_2)$ with respect to $t_1$ yields:

$$\begin{aligned} {\frac{\partial \varphi ^{\mathrm {LG}}}{\partial t_1} (t_1,t_2)} = {\frac{\partial Y_{\varphi }}{\partial t_1} (t_1,t_2)} \ Y_{\varphi } (t_1,t_2)^\top \!+ Y_{\varphi } (t_1,t_2) \left( {\frac{\partial Y_{\varphi }}{\partial t_1} (t_1,t_2)} \right) ^\top \!, \end{aligned}$$

with

$$\begin{aligned} {\frac{\partial Y_{\varphi }}{\partial t_1} (t_1,t_2)} = (1-t_2) \dot{Y}_{1-2}(t_1) + t_2 \dot{Y}_{3-4}(t_1) Q(t_1)^\top \!+ t_2 Y_{3-4}(t_1) \dot{Q}(t_1)^\top \!. \end{aligned}$$

The values of $\dot{Y}_{1-2}(t_1)$ and $\dot{Y}_{3-4}(t_1)$ are independent of $t_1$:

$$\begin{aligned} \dot{Y}_{1-2}(t_1) = - Y_1 + Y_2 Q_{1-2}^\top \!\, , \qquad \dot{Y}_{3-4}(t_1) = - Y_3 + Y_4 Q_{3-4}^\top \!\,, \qquad \forall t_1. \end{aligned}$$

The value of $\dot{Q}(t_1)$ can be obtained as follows. Recall from the geodesic definition that $Q(t_1)$ is the orthogonal factor of the polar decomposition of the matrix $M(t_1) = Y_{1-2}(t_1)^\top \!Y_{3-4}(t_1),$ which means that there exists a symmetric positive definite matrix $H(t_1)$ such that $M(t_1) = H(t_1)Q(t_1)$. Then,

$$\begin{aligned} \dot{M}(t_1) = \dot{H}(t_1) Q(t_1) + H(t_1) \dot{Q}(t_1), \end{aligned}$$

(A.1)

where $\dot{H}(t_1)$ is a symmetric matrix, and $\dot{Q}(t_1)$ is of the form $\dot{Q}(t_1) = \varOmega (t_1) Q(t_1)$, with $\varOmega (t_1)$ a skew-symmetric matrix. Right-multiplying this expression by $Q(t_1)^\top \!$ yields:

$$\begin{aligned} \dot{M}(t_1) Q(t_1)^\top \!= \dot{H}(t_1) + H(t_1) \varOmega (t_1), \end{aligned}$$

(A.2)

while left-multiplying the transpose of Eq. (A.1) by $-Q(t_1)$ yields:

$$\begin{aligned} - Q(t_1) \dot{M}(t_1)^\top \!= - \dot{H}(t_1) + \varOmega (t_1) H(t_1). \end{aligned}$$

(A.3)

Now, summing equations (A.2) and (A.3) yields:

$$\begin{aligned} \dot{M}(t_1) Q(t_1)^\top \!- Q(t_1) \dot{M}(t_1)^\top \!= H(t_1) \varOmega (t_1) +\varOmega (t_1) H(t_1). \end{aligned}$$

As a result, the term $\dot{Q}_1(t_1)$ can be obtained by solving a Sylvester equation. Moreover, since $H(t_1)$ is always positive definite (except in the set of zero measure corresponding to low-rank matrices $Y_{1-2}(t_1)^\top \!Y_{3-4}(t_1)$), the solution to the Sylvester equation is unique ($H(t_1)$ and -$H(t_1)$ have no common eigenvalues).

1.2 A.2 Computation of $t_2^*(t_1)$

The first-order optimality condition

$$\begin{aligned} \left. \frac{\partial f}{\partial t_2}\right| _{(t_1 t_2^*(t_1))} = 0 \end{aligned}$$

implies that the optimal value $t_2^*(t_1)$ corresponding to an arbitrary value $t_1$ is the solution to a cubic equation:

$$\begin{aligned} s_1(t_1) t_2^3(t_1) + s_2(t_1) t_2^2(t_1) + s_3(t_1) t_2(t_1) + s_4(t_1) = 0, \end{aligned}$$

(A.4)

with

$$\begin{aligned}&s_1(t_1) = 2 \mathrm {tr}\left( R^2 \right) = 2 \sum _i \sum _j R_{ij}^2, \\&s_2(t_1) = 3 \mathrm {tr}\left( R S \right) = 3 \sum _i \sum _j R_{ij}S_{ij}, \\&s_3(t_1) = 2 \mathrm {tr}\left( R T \right) + \mathrm {tr}\left( S \right) = 2 \sum _i \sum _j R_{ij}T_{ij} + S_{ij}^2 ,\\&s_4(t_1) = 2 \mathrm {tr}\left( S T \right) = 2 \sum _i \sum _j S_{ij}T_{ij}. \end{aligned}$$

The matrices R, S, and T arising in those expressions are defined as:

$$\begin{aligned} R&= Y_{1-2}Y_{1-2}^\top \!+ Y_{3-4}Y_{3-4}^\top \!- \left( Y_{1-2} Q Y_{3-4}^\top \!+ Y_{3-4} Q^\top \!Y_{1-2}^\top \!\right) , \\ S&= \left( Y_{1-2} Q Y_{3-4}^\top \!+ Y_{3-4} Q^\top \!Y_{1-2}^\top \!\right) - 2 Y_{1-2}Y_{1-2}^\top \!\;, \\ T&= Y_{1-2}Y_{1-2}^\top \!- \widehat{C}, \end{aligned}$$

with all these matrices depending on $t_1$. Observe, however, that for a fixed value of $t_1$, the function $t_2 \rightarrow f(t_1,t_2)$ might not be convex; hence, the condition (A.4) might have several (up to three) real solutions. In that case, we compute the value of the cost function at those solutions, and we choose the one corresponding to the smallest value of f.

1.3 A.3 Gradient descent for the univariate cost function

We are now looking for the derivative of the cost function $\tilde{f}(t_1) = f(t_1, t_2^*(t_1))$, with respect to the variable $t_1$, in order to be able to apply a steepest descent method to that problem. Using the notation $\tilde{f} = F \circ \tilde{\varphi }^{\mathrm {LG}}$, with $\tilde{ \varphi }^{\mathrm {LG}}(t_1) = \varphi ^{\mathrm {LG}}(t_1,t_2^*(t_1))$, we have:

$$\begin{aligned} \dot{\tilde{f}}(t_1) = DF[ \dot{\tilde{\varphi }}^{\mathrm {LG}} (t_1)] = 2 \mathrm {tr}\left( \dot{\tilde{\varphi }}^{\mathrm {LG}}(t_1) (\tilde{\varphi }^{\mathrm {LG}}(t_1) - \widehat{C})^\top \! \right) . \end{aligned}$$

The derivative $\dot{\tilde{\varphi }}^{\mathrm {LG}}(t_1)$ is given by:

$$\begin{aligned} \dot{\tilde{\varphi }}^{\mathrm {LG}}(t_1) = \dot{Y}_{\tilde{\varphi }}(t_1) Y_{\tilde{\varphi }}(t_1)^\top \!+ Y_{\tilde{\varphi }} \dot{Y}_{\tilde{\varphi }}(t_1)^\top \!. \end{aligned}$$

Using the chain rule,

$$\begin{aligned} \dot{Y}_{\tilde{\varphi }}(t_1) = {\frac{\partial Y_{\varphi }}{\partial t_1}(t_1,t_2^*(t_1))} + {\frac{\partial Y_{ \varphi }}{\partial t_2}(t_1,t_2^*(t_1))} \dot{t}_2^*(t_1). \end{aligned}$$

By the definition of $t_2^*(t_1)$, the term $\frac{\partial Y_{\varphi }}{\partial t_2}(t_1,t_2^*(t_1))$ is equal to zero. As a result, $\dot{Y}_{\tilde{\varphi }}(t_1) = \frac{\partial Y_{\varphi }}{\partial t_1}(t_1,t_2^*(t_1))$, which has been computed earlier.