Distances and Riemannian metrics for multivariate spectral densities


Preliminaries on multivariate prediction

Consider a multivariate discrete-time, zero mean, weakly stationary stochastic process {{bf u}(k),~ kinmathbb{Z}} with {bf u}(k) taking values in mathbb{C}^{mtimes 1}. Let

 R_k=mathcal{E}left{{bf u}(ell) {bf u}^*(ell-k)  right} ~mbox{~for~} l,kinmathbb{Z}

denote the sequence of matrix covariances and dmu(theta) be the corresponding matricial power spectral measure for which

 R_k=int_{-pi}^{pi}  e^{-{rm j} ktheta}{frac{dmu(theta)}{2pi}}.

We will be concerned with the case of non-deterministic process of full rank with an absolutely continuous power spectrum. Hence, dmu(theta)=f(theta)dtheta with f(theta) being a matrix valued power spectral density (PSD) function. Further, log det f(theta) is assumed to be integrable throughout.

Geometry of multivariable process

Define L_{2}({bf u}) to be the closure of mtimes 1-vector-valued finite linear combinations of {{bf u}(k)} with respect to convergence in the mean:

 L_{2}({bf u}):=overline{left{sum_{rm finite} P_k{bf u}(-k) ;:; P_kinmathbb{C}^{mtimes m},; kinmathbb{Z}  right}}.

This space is endowed with both, a matricial inner product

 llbracket sum_k P_k{bf u}(-k), sum_k Q_k {bf u}(-k) rrbracket :=  mathcal{E} left{ left( sum_k P_k {bf u}(-k) right) left( sum_k Q_k{bf u}(-k)right)^* right},

as well as a scalar inner product

 langle sum_k P_k{bf u}(-k), sum_k Q_k {bf u}(-k) rangle= {rm trace} llbracket  sum_k P_k{bf u}(-k), sum_k Q_k {bf u}(-k) rrbracket .

It is standard to establish the Kolmogorov isomorphism between the “temporal” space L_2({bf u}) and the “spectral” space L_2(fdtheta),

 varphi ;:; L_2({bf u}) to L_{2}(fdtheta) ;:; sum_k P_k {bf u}(-k)mapsto sum_k P_kz^k,

with z=e^{{rm j}theta} for thetain[-pi, pi]. It is convenient to endow the latter space L_2(fdtheta) the matricidal inner product

 llbracket p, qrrbracket_{fdtheta}:=int_{-pi}^{pi} left(p(e^{{rm j} theta}){frac{fdtheta}{2pi}} q(e^{{rm j} theta})^*right).

We denote llbracket p, qrrbracket_{f}:=llbracket p, qrrbracket_{fdtheta} for simplicity. The corresponding scalar inner product

 langle p, qrangle_{f}:={rm trace} llbracket p, qrrbracket_{f}.

The least-variance linear prediction

     minleft{ {rm trace}mathcal{E}{{bf p} {bf p}^* } : {bf p}={bf u}(0)-sum_{k>0}P_k{bf u}({-k}),;P_kinmathbb{C}^{mtimes m} right}

can be expressed in the spectral domain

  min left{ langle p,prangle_f : p(z)=I-sum_{k>0}P_kz^k,;P_kinmathcal{C}^{mtimes m}  right}.

The minimizer coincides with that of

  min left{ llbracket p,prrbracket_f : p(z)=I-sum_{k>0}P_kz^k,;P_kinmathcal{C}^{mtimes m}  right},

where the minimum is sought in the positive-definite sense. Let p(z) denote the the minimizer of such a problem, then the minimal matrix of the above problem, denoted by Omega, is

  Omega:=llbracket p, prrbracket_f.

Spectral factorization and optimal prediction

For a non-deterministic process of full rank, the determinant of the error variance Omega is non-zero, and

  det Omega = exp{int_{-pi}^{pi} logdet f(theta) frac{dtheta}{2pi}}

this the is well-known Szegddot{o}-Kolmogorov formula. We consider only non-deterministic processes of full rank and hence we assume that

log det f(theta)in L_1[-pi, pi].
In this case, f(theta) admits a unique factorization
 f(theta)=f_+(e^{{rm j}theta})f_+(e^{{rm j}theta})^*,

with f_+(e^{{rm j}theta})in mathcal{H}_2^{mtimes m}(mathbb{D}) and f_+(0)=Omega^{frac12}, where {mathbb D}:={z: |z|<1}, and mathcal{H}_2(mathbb{D}) denotes the Hardy space of functions which are analytic in the unit disk mathbb{D} with square-integrable radial limits. The spectral factorization presents an expression of the optimal prediction error in the form


Comparison of PSD's

Prediction errors and innovations processes

Consider two processes { {bf u}_i(k), ~iin{1,2 } } with f_i and p_i=I-sum_{k>0} P_{i,k}z^k the corresponding PSD's and the optimal prediction filters, respectively. Let

 {bf p}_{i,j}(k):= {bf u}_{i}(k)-P_{j,1}{bf u}_i(k-1)-P_{j,2}{bf u}_i(k-2)-ldots

for i,jin{1,2 } be an innovations process. Then, from stationarity,

 llbracket {bf p}_{i,i}(k), {bf p}_{i,i}(k) rrbracket =Omega_i,


 llbracket {bf p}_{i,j}(k), {bf p}_{i,j}(k) rrbracket =:Omega_{i,j}geq Omega_i. hspace*{3cm}(1)

The color of innovations and PSD mismatch

We normalize the innovation processes as follows:

 {bf h}_{i,j}(k) =  Omega_j^{-frac12} {bf p}_{i,j}(k),  ~ kinmathbb{Z},

then the power spectral density of the process {{bf h}_{i,j}(k) } is

 f_{{bf h}_{i,j}} =f_{j+}^{-1}f_i f_{j+}^{-*}.

The mismatch between the two power spectra f_i, f_j can be quantified by the distance of f_{{bf h}_{i,j}} to the identity. To this end, we define

 begin{array}{rl} {rm D}_1(f_1,f_2):=&{rm trace} left( f_{j+}^{-1}f_i f_{j+}^{-*}-Iright)+ {rm trace}left( f_{i+}^{-1}f_j f_{i+}^{-*}-I right) frac{dtheta}{2pi}   =& int_{-pi}^{pi} | f_1^{-1/2}f_2^{1/2}- {f_1^{1/2}}{f_2^{-1/2}}|_{rm Fr}^2 frac{dtheta}{2pi}. end{array}

The following choices of divergence measures lead to the same Riemannian structure as the above one:

 begin{array}{l} sum_{i,j} int_{-pi}^{pi} |f_{j+}^{-1}f_i f_{j+}^{-*}- I|_{rm Fr}^2  frac{dtheta}{2pi}  sum_{i,j} int_{-pi}^{pi} |(f_{j+}^{-1}f_i f_{j+}^{-*})^{frac12}- I|_{rm Fr}^2  frac{dtheta}{2pi}   int_{-pi}^{pi} left({rm trace} (f_{2+}^{-1}f_1 f_{2+}^{-*})-logdet(f_{2+}^{-1}f_1 f_{2+}^{-*})-mright) frac{dtheta}{2pi} int_{-pi}^{pi} |log(f_{1+}^{-1}f_2 f_{1+}^{-*})|_{rm Fr}^2  frac{dtheta}{2pi}. end{array}

These are the Frobenius distance, the generalized Hellinger distance, the multivariate Itakuta-Saito distance, the log-spectral deviation between (f_{j+}^{-1}f_i f_{j+}^{-*})^{frac12} and I, respectively.

Suboptimal prediction and PSD mismatch

As was given in (1), Omega_{i,j}geq Omega_i. The above equality holds if and only if p_i=p_j. Thus a mismatch between the two spectral densities can be quantified by the strength of the above inequality. Since Omega_i^{-frac12}Omega_{i,j}Omega_i^{-frac12}geq I, we consider

 {rm D}_2(f_i,f_j):= logdetleft(Omega_i^{-frac12}Omega_{i,j}Omega_i^{-frac12}right)

as a “divergence measure” between the two PSD's. The following options lead to the same Riemannian structure:

 begin{array}{l} frac{1}{m}{rm trace}(Omega_i^{-frac12}Omega_{i,j}Omega_i^{-frac12})-1 det(Omega_i^{-frac12}Omega_{i,j}Omega_i^{-frac12})-1. end{array}

Riemannian structure on multivariate spectra

Consider the following class of mtimes m PSD's

 mathcal{F}:= {f | f(theta)>0, forall thetain[-pi, pi], df(theta)/dtheta mbox{~continuous} }.

Let mathcal{D} be a class of admissible perturbations of Delta for fin mathcal{F} with

 mathcal{D}:= {Delta mid Delta(theta)=Delta(theta)^*, d Delta(theta)/dtheta mbox{~continuous} }.

It was shown in [1] that for PSD's in mathcal{F} and Delta in mathcal{D} the distance {rm D}_1 induces the following Riemannian metric

 {color{blue} {rm g}_{1,f}(Delta):= int_{-pi}^{pi} |f^{-1/2} Delta f^{-1/2} |^2_{rm Fr}frac{dtheta}{2pi}. }

In particular, for (f,Delta)inmathcal{F}times mathcal{D} and epsilon >0. Then for epsilon sufficiently small

 {rm D}_1(f,f+epsilonDelta) ={rm g}_{1,f}(epsilonDelta)+O(epsilon^3).

The geodesic path connecting two spectra f_0, f_1in mathcal{F} is

 {color{blue} f_{tau}=f_0^{1/2}(f_0^{1/2} f_1f_0^{1/2})^tau f_0^{1/2}, mbox{~for~} 0leqtauleq1. }

The geodesic distance is

 d_{{rm g}_1}(f_0, f_1)=sqrt{int_{-pi}^{pi} |log f_0^{-1/2}f_1f_0^{-1/2} |_{rm Fr}^2 frac{dtheta}{2pi}}.

For the distance {rm D}_2, the corresponding Riemannian metic is

 {color{blue} {rm g}_{2,f}(Delta):= {rm trace}int_{-pi}^{pi} hspace*{-2pt}(f^{-1}_+Delta f^{-*}_+)^2frac{dtheta}{2pi}-{rm trace}big(int_{-pi}^{pi} hspace*{-2pt}f^{-1}_+Delta f^{-*}_+frac{dtheta}{2pi}big)^2. }

The associtated geodesic path for {rm g}_{2,f}(Delta) is still unknown.


A scalar example

Consider the two power spectral densities

 f_i(theta)=frac{1}{|a_i(e^{{rm j}theta})|^2}, ;iin {0,1 },


 begin{array}{rl} a_0(z)=&(z^2-1.96cos(frac{pi}{5})+0.98^2)(z^2-1.7cos(frac{pi}{3})+0.85^2)(z^2-1.8cos(frac{2pi}{3})+0.9^2), a_1(z)=&(z^2-1.96cos(frac{2pi}{15})+0.98^2)(z^2-1.5cos(frac{7pi}{30})+0.75^2)(z^2-1.8cos(frac{5pi}{8})+0.9^2). end{array}

The two power spectra, log f_0(theta) and log f_1(theta) , are shown in the following figure.


We evaluate tauin left{frac{1}{3}, frac{2}{3}, frac{4}{3} right} , the corresponding spectra are shown in the following figure.


A multivariable example

Consider the two matrix-valued power spectral densities

 f_0=left[       begin{array}{cc}         1 & 0          0.1e^{{rm j}theta} & 1        end{array}     right]left[       begin{array}{cc}         frac{1}{|a_0(e^{{rm j}theta})|^2} & 0          0 & 1        end{array}     right]left[       begin{array}{cc}         1 & 0.1e^{-{rm j}theta}         0 & 1        end{array}     right]
 f_1=left[       begin{array}{cc}         1 & 0.1e^{{rm j}theta}          0 & 1        end{array}     right]left[       begin{array}{cc}         1 & 0          0 & frac{1}{|a_1(e^{{rm j}theta})|^2}        end{array}     right]left[       begin{array}{cc}         1 & 0          0.1e^{-{rm j}theta} & 1        end{array}     right].

The following two figures show the two spectra where the (1,2) block is the log|f_i(1,2)| and the (2,1) block is arg(f_i(2,1)), for iin {0,; 1}.


We compute the geodesic connecting f_0 and f_1 as

 f_tau=f_0^{frac{1}{2}}(f_0^{-frac{1}{2}}f_1f_0^{-frac{1}{2}})^tau, mbox{~for~} 0leqtauleq1.

The geodesic is shown in the following figure.




[1] X. Jiang, L. Ning and T.T. Georgiou, “Distances and Riemannian metrics for multivariate spectral densities,” IEEE Trans. on Signal Processing, to appear 2012. http://arxiv.org/abs/1107.1345