Distances and Riemannian metrics for multivariate spectral densities

Consider a multivariate discrete-time, zero mean, weakly stationary stochastic process {{bf u}(k),~ kinmathbb{Z}}

with

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

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}}.$

$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},$

$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}.$

$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}$

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

$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

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 Szeg ddot{o}

-Kolmogorov formula. We consider only non-deterministic processes of full rank and hence we assume that

$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

with square-integrable radial limits. The spectral factorization presents an expression of the optimal prediction error in the form

$p(z)=f_+(0)f_+^{-1}(z).$

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$

$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

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

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

The mismatch between the two power spectra f_i

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

, 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

$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

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}. }$

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

${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. }$

$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}}.$

${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. }$

Examples

A scalar example

$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}$

$f_{tau}=f_0^{1-tau}f_1^{tau}.$

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

$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)