# Bounds for the estimation of matrix-valued parameters of a Gaussian random process

Lorena Leon Arencibia, Herwig Wendt, Jean-Yves Tourneret

<figure><img src="image_1.png" /></figure>

## To cite this version:

Lorena Leon Arencibia, Herwig Wendt, Jean-Yves Tourneret. Bounds for the estimation of matrix-valued parameters of a Gaussian random process. Signal Processing, 2023, 211, pp.109106. 10.1016/j.sigpro.2023.109106. hal-04254235

HAL Id: hal-04254235

https://hal.science/hal-04254235v1

Submitted on 23 Oct 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
**Goals, outline and contributions.** This paper, we consider $\Sigma_1, \Sigma_2$ to be unknown with inverse Wishart (IW) prior distributions. Our goal is to derive lower bounds for the MSE of estimators of $\Sigma_1, \Sigma_2$. Assuming the estimation to be conducted within a Bayesian formulation, Section 2 derives the BBs of $\Sigma_1$ and $\Sigma_2$ and analytically studies their properties. Section 3 studies the properties of the bounds in this framework using Monte Carlo simulations. The use of the proposed bounds for the parameters associated with the bivariate multifractal spectrum is finally illustrated in Section 3.3. The main contributions of this paper are i) the derivation of the BB for the above specified statistical model, which is a new theoretical result obtained from (4), (8), (9) and (11), ii) the derivation of a novel closed-form expression for computing non-trivial expectations involving Wishart random matrices, see Proposition 1, iii) the study of the analytic properties of the bounds (see Section 2) and iv) extensive numerical experiments and results that validate and illustrate the obtained theoretical expressions of the bounds (see Section 3).

## 2. Bayesian Cramér-Rao bound

This section derives the BB for the MSE of estimators of $\Sigma_1, \Sigma_2$, when $\Sigma_1$ and $\Sigma_2$ are assigned independent IW priors, i.e., for $i \in \{1, 2\}, \Sigma_i \sim \text{IW}(v_i, \Omega_i)$, with $v_i$ degrees of freedom ($v_i \in \mathbb{R}$ and $v_i > R+1$), and mean matrix $(v_i - R - 1)^{-1}\Omega_i$, where $\Omega_i$ is a real-valued p.d. scale matrix. To this end, we make use of the following novel results.

**Proposition 1.** Moments of the type $E[\text{WAWBW}]$.

If $\Sigma \sim \text{IW}(v, \Omega)$, then $W = \Sigma^{-1}$ has the Wishart distribution $W(v, \Delta = \Omega^{-1})$. Then, for any pair of real-valued symmetric matrices $(A, B)$:

$$ \mathbb{E}_{\Sigma} [ \begin{aligned} & \Sigma^{-1} A \Sigma^{-1} B \Sigma^{-1} \\ & = \Delta A \Delta B \Delta (v^3 + 2v^2 + v) + \Delta B \Delta A \Delta (v^2 + 3v) \\ & \quad + [\text{tr}(\Delta A)] \Delta B \Delta (v^2 + v) + [\text{tr}(\Delta B)] \Delta A \Delta (v^2 + v) \\ & \quad + \Delta [(v^2 + v)\text{tr}(\Delta A \Delta B) + v\text{tr}(\Delta A)\text{tr}(\Delta B)] \end{aligned} ] \quad (1) $$

The proof of Proposition 1 can be conducted using the approach detailed in [12]. Moreover, according to [13], for any real-valued symmetric matrix $A$, we have

$$ \mathbb{E}_{\Sigma} [\Sigma^{-1} A \Sigma^{-1}] = (v^2 + v) \Delta A \Delta + v \text{tr}(\Delta A) \Delta. \quad (2) $$

**Definitions.** Let $\theta \in \mathbb{R}^p$, with $p = R^2 + R$, the vector obtained by concatenating the vectors $\text{vec}_{\text{triu}}(\Sigma_1)$ and $\text{vec}_{\text{triu}}(\Sigma_2)$, where the matrix operator $\text{vec}_{\text{triu}}(A)$ returns the vector of the elements of the upper triangular part of $A$. Note that the first and the last $p/2$ elements of $\theta$, denoted as $\theta_1:_{\frac{p}{2}}$ and as $\theta_2:_{1+p}$, correspond to the main diagonal and all elements of $\Sigma_1$ and $\Sigma_2$ above the diagonal, respectively. In the following, the matrices $R_n$ will be denoted as $R_n(\theta)$ to emphasize the dependence of $R_n$ on $\theta$. The evaluation of the BBs of $\Sigma_1$ and $\Sigma_2$ requires to invert the posterior Fisher information matrix (PFIM) defined as [1]

$$ F = \mathbb{E}_{z, \Sigma_1, \Sigma_2} \left[ - \frac{\partial^2 L(z, \Sigma_1, \Sigma_2)}{\partial \theta \partial \theta^T} \right], \quad (3) $$

where $L(z, \Sigma_1, \Sigma_2)$ is the joint log-likelihood of the model, which is twice differentiable with respect to (w.r.t.) $\theta$ and has a bounded support independent of $\theta$. These required regularity conditions ensure the existence of the BB. Equation (3) can be rewritten as

$$ \begin{aligned} F &= \mathbb{E}_{\Sigma_1, \Sigma_2} \left[ \mathbb{E}_{z | \Sigma_1, \Sigma_2} \left[ -\frac{\partial^2 L(z | \Sigma_1, \Sigma_2)}{\partial \theta \partial \theta^T} \right] - \frac{\partial^2 \pi_1(\Sigma_1)}{\partial \theta \partial \theta^T} - \frac{\partial^2 \pi_2(\Sigma_2)}{\partial \theta \partial \theta^T} \right] \\ &= \mathbb{E}_{\Sigma_1, \Sigma_2} \left[ F_{\theta} + F_{\Sigma_1} + F_{\Sigma_2} \right], \end{aligned} \quad (4) $$

where $L(z | \Sigma_1, \Sigma_2)$ is the log-likelihood of $z$ given $\Sigma_1$ and $\Sigma_2$, which can be expressed as

$$ L(z | \Sigma_1, \Sigma_2) = K - \ln \det \mathbf{R}(\theta) - z^H \mathbf{R}(\theta)^{-1} z, \quad (5) $$

with $K = -N \ln \pi$. Moreover, $\pi_i(\Sigma_i)$ is the log-prior of $\Sigma_i$ defined as

$$ \pi_i(\Sigma_i) = -((v_i + R + 1)/2) \ln \det \Sigma_i - \frac{1}{2} \text{tr}(\Omega_i \Sigma_i^{-1}) + \text{constant}, \quad (6) $$

where $\text{tr}(\cdot)$ denotes the trace operator. Note that functions (5) and (6) also satisfy the regularity conditions ensuring the existence of the BBs of $\Sigma_1$ and $\Sigma_2$.

**Computation of $F_{\theta}, F_{\Sigma_1}$ and $F_{\Sigma_2}$**. The matrix $F_{\theta}$ is known as the Fisher information matrix of $\theta$. Since $z$ is a zero mean complex circularly-symmetric Gaussian vector, the element of $F_{\theta}$ located at the $k$th row and $l$th column, denoted as $[F_{\theta}]_{kl}$ for $k, l \in \{1, \dots, p\}$, can be calculated as [14,15]

$$ [F_{\theta}]_{kl} = \text{tr} \left\{ \mathbf{R}^{-1}(\theta) \frac{\partial \mathbf{R}(\theta)}{\partial \theta_k} \mathbf{R}^{-1}(\theta) \frac{\partial \mathbf{R}(\theta)}{\partial \theta_l} \right\}. \quad (7) $$

Note that in general the computation of (7) needs $O(R^3N^3)$ operations because it requires to invert the matrix $\mathbf{R}(\theta)$. This computation can lead to a high computational cost for large values of $N$ and $R$. We propose to overcome this limitation by exploiting the block diagonal structure of $\mathbf{R}(\theta)$. This structure allows the inverse of the matrix $\mathbf{R}(\theta)$ to be computed using the inverses of the $R \times R$ diagonal blocks $\mathbf{R}_n(\theta)$ of $\mathbf{R}(\theta)$. Specifically, $\mathbf{R}_n^{-1}(\theta) = (\Sigma_1 g_1(n) + \Sigma_2 g_2(n))^{-1}$ and

$$ \frac{\partial \mathbf{R}_n(\theta)}{\partial \theta_l} = \mathbf{B}_{n,l} = \begin{cases} J_{1,l} g_1(n) & \text{if } l \in \{1, \dots, \frac{p}{2}\}, \\ J_{2,l} g_2(n) & \text{if } l \in \{\frac{p}{2} + 1, \dots, p\}, \end{cases} $$

where $J_{i,l} = \frac{\partial \Sigma_i}{\partial \theta_l}$ does not depend on $\theta$, hence $\frac{\partial^2 \Sigma_i}{\partial \theta_k \partial \theta_l} = 0$. Thus, (7) can be rewritten as

$$ [F_{\theta}]_{kl} = \sum_{n=1}^{N} \text{tr} \left\{ \mathbf{R}_n^{-1}(\theta) \mathbf{B}_{n,k} \mathbf{R}_n^{-1}(\theta) \mathbf{B}_{n,l} \right\}, \quad (8) $$

which can be computed with $O(R^3N)$ operations. Overall, the values of $R$ considered in this paper are relatively small $R \ll N$ and lead to feasible computational times.

On the other hand, the following result is obtained

$$ \begin{aligned} [F_{\Sigma_i}]_{kl} &= - \frac{\partial^2 \pi_i(\Sigma_i)}{\partial \theta_k \partial \theta_l} \\ &= ((v_i + R + 1)/2) \text{tr}(\Sigma_i^{-1} J_{i,k} \Sigma_i^{-1} J_{i,l}) \\ &\quad -(1/2) \text{tr}(\Omega_i [\Sigma_i^{-1} J_{i,k} \Sigma_i^{-1} J_{i,l} \Sigma_i^{-1} + \Sigma_i^{-1} J_{i,l} \Sigma_i^{-1} J_{i,k} \Sigma_i^{-1}]). \end{aligned} \quad (9) $$

**Expectations.** This section explains how to compute the expectations in (4). Calculating the expectation $\mathbb{E}_{\Sigma_1, \Sigma_2} [F_{\Sigma_1} + F_{\Sigma_2}]$ reduces to determining the expectation of $-\frac{\partial^2 \pi_1(\Sigma_1)}{\partial \theta_1:_{\frac{p}{2}} \partial \theta_1:_{\frac{p}{2}}}$ w.r.t.

$$ \Sigma_1 \text{ and the expectation of } -\frac{\partial^2 \pi_2(\Sigma_2)}{\partial \theta_1:_{\frac{p}{2}} \partial \theta_2:_{1+p} \partial \theta_2:_{1+p}} \text{ w.r.t. } \Sigma_2, \text{ since the} $$

other terms equal 0. Both expectations have a closed-form expression that can be determined using the matrix expectations (1) and (2). The challenge here is to compute the expectation $\mathbb{E}_{\Sigma_1, \Sigma_2} [F_{\theta}]$, which involves calculating the expectation of the expression $\mathbf{R}_n^{-1} \mathbf{B}_{n,k} \mathbf{R}_n^{-1} \mathbf{B}_{n,l}$ with respect to $\Sigma_1, \Sigma_2$ for all $n = 1, \dots, N$.
This computation is possible provided we can compute the expectation

$$ \mathbb{E}_{\Sigma_1, \Sigma_2} \left[ \left( (a\Sigma_1 + b\Sigma_2)^{-1} A (a\Sigma_1 + b\Sigma_2)^{-1} B \right) \right] \quad (10) $$

for $a, b \in \mathbb{R}_+$, $\Sigma_1 \sim \text{IW}(\nu_1, \Omega_1)$, $\Sigma_2 \sim \text{IW}(\nu_2, \Omega_2)$ and any pair of symmetric matrices $A$ and $B$. If $a$ or $b$ equal zero, (10) can be calculated using (2). Otherwise, we propose to approximate (10) numerically via a Monte Carlo algorithm. Given $a$, $b$, $A$ and $B$, we can generate a large number $S$ of samples $\{\Sigma_1^{(s)}, \Sigma_2^{(s)}\}_{s=1}^S$ according to IW distributions, compute $E^{(s)} = \left[ (a\Sigma_1^{(s)} + b\Sigma_2^{(s)})^{-1} A (a\Sigma_1^{(s)} + b\Sigma_2^{(s)})^{-1} B \right]$ and approximate (10) by the average of the matrices $\{E^{(s)}\}_{s=1}^S$.

The MSE of any estimator $(\hat{\Sigma}_1, \hat{\Sigma}_2)$ of $(\Sigma_1, \Sigma_2)$ is defined as the trace of the error covariance matrix MSE $= \text{tr}(\mathbb{E}_{\hat{\theta}}[(\hat{\theta} - \theta)(\hat{\theta} - \theta)^T])$, where $\hat{\theta}$ is the vector obtained by concatenating the vectors $\text{vec}_{\text{triu}}(\hat{\Sigma}_1)$ and $\text{vec}_{\text{triu}}(\hat{\Sigma}_2)$. Finally, the inverse of the PFM (3) denoted as $[\mathbf{F}]^{-1}$ yields the desired lower bound for the MSE of any estimator $(\hat{\Sigma}_1, \hat{\Sigma}_2)$ of $(\Sigma_1, \Sigma_2)$, when these matrices are IW distributed:

$$ \text{MSE} \ge \text{BB} = \text{tr}([\mathbf{F}]^{-1}). \quad (11) $$

**Analytic properties.** Assuming that $a$ or $b$ are zero and $\Delta$ is diagonal, (11) can be computed in closed-form. In particular, the entries of $\mathbf{F}$ are given by

$$ [\mathbf{F}]_{kl} = \gamma_1 \text{tr}(\mathbf{\Delta}_J_k) \text{tr}(\mathbf{\Delta}_J_l) + \gamma_2 \text{tr}(\mathbf{\Delta}_J_k \mathbf{\Delta}_J_l) \quad (12) $$

where $\gamma_1 = \frac{1}{2}(\nu^2 + \nu(2N - R + 3))$ and $\gamma_2 = \frac{1}{2}(\nu^3 + \nu^2(2N - R + 6) + \nu(2N + 9 - R))$. When $\Delta$ is the identity matrix, the first term of the PFM is a matrix with zero entries except for an $R \times R$ diagonal block with non-zero entries whereas the second term is a diagonal matrix. In this case, it can be shown that the PFM has $R - 1$ eigenvalues (ev.) equal to $\gamma_2$, $R(R+1)/2 - R = (R^2 - R)/2$ ev. equal to $2\gamma_2$ and one ev. equal to $R\gamma_1 + \gamma_2$. The trace of $[\mathbf{F}]^{-1}$ with the above assumptions, denoted as aBB (for approximate BB), is

$$ \text{aBB} = \text{tr}([\mathbf{F}]^{-1}) = \frac{R(R-1)}{4\gamma_2} + \frac{R((R-1)\gamma_1 + \gamma_2)}{\gamma_2^2 + R\gamma_1\gamma_2}. \quad (13) $$

This shows that the bound behaves asymptotically as $\nu^{-3}$ as $\nu \to +\infty$, and as a second order polynomial in $\nu$ as $\nu \to R + 1$. Moreover, the behavior is asymptotically linear in $R$, and a second order polynomial in $R$ as $R \to 1$. As expected, the asymptotic decay with sample size is $N^{-1}$.

### 3. Numerical illustrations

In this section, extensive numerical simulations are used to study the properties of the BB for the MSE of any estimator of the pair of matrices $\Lambda = (\Sigma_1, \Sigma_2)$ in the above probabilistic framework, and compare the bounds against the MSEs of Bayesian estimators.

#### 3.1. Monte Carlo simulations

**Estimation algorithms.** We consider the maximum a posteriori (MAP) and minimum mean square (MMSE) estimators, defined by:

$$ \hat{\Lambda}^{\text{MAP}} = \underset{\Lambda}{\operatorname{argmax}} p(\Lambda | \mathbf{z}), \quad (14) $$

$$ \hat{\Lambda}^{\text{MMSE}} = \mathbb{E}[\Lambda | \mathbf{z}]. \quad (15) $$

In order to compute (14) and (15), we use a Gibbs sampler [16] to generate a large collection $N_{\text{mc}}$ of samples distributed according the posterior distribution equation to be recalled here. After

a burn-in period, where the first $N_{\text{bi}}$ samples are discarded, the Bayesian estimators $\hat{\Lambda}^{\text{MAP}}$ and $\hat{\Lambda}^{\text{MMSE}}$ are approximated using the last generated samples.

**Simulation setup.** Unless otherwise stated $R = 2$, $N = 2^8$, $\Omega_1 = \Omega_2 = \mathbb{I}_R$ ($R \times R$ identity matrix) and $\nu_1 = \nu_2 = 80$. Without loss of generality, we use the functions $g_1(n) = 2\pi \cos^2(x[n]) + 0.1$ and $g_2(n) = 2\pi \sin^2(x[n]) + 0.1$, where $x$ is the vector of $N$ components whose values have been generated in the interval $[0,2]$, equidistanted with a distance of $2/(N-1)$. In all cases, the sample MSE of the estimators is computed as the average of the trace of the error covariance matrix over 1000 independent realizations. Gibbs samplers are run with $N_{\text{mc}} = 1000$ and $N_{\text{bi}} = 500$, and (10) is approximated (when needed) as described before using $S = 200$.

#### 3.2. Performance analysis

**Performance vs. sample sizes.** Fig. 1 displays the BB, its approximation aBB, and the MSEs of the MMSE and MAP estimators, for various sample sizes, where $\Sigma_1, \Sigma_2$ are random matrices with IW prior distributions. The following comments are appropriate:

* The BB decreases as $N^{-1}$ when $N \to \infty$ and to a constant when $N \to 0$.

* The BB and its approximation aBB are asymptotically close but tend to different constants for small sample size.

* **BB vs. MSE:** The MSEs of both estimators are approaching the BBs when $N$ increases - the more data, the tighter the bounds.

* **MMSE vs. MAP:** Overall, the MMSE estimator has better performance than the MAP estimator, in particular for small sample size. This result was expected since the MMSE estimator minimizes the MSE.

**Performance vs. degrees of freedom.** Fig. 2 compares the MSEs of the MMSE and MAP estimators and the BBs for various degrees of freedom, $\nu_1 = \nu_2 \in \{10, 15, 20, 25, \dots, 120\}$. We can observe that the BB decreases when $\nu_1$ and $\nu_2$ increase; indeed, in that case, the priors are more informative. The approximation aBB is very similar to BB and predicts that this decay is of order $\nu^{-3}$. Moreover, the values of the MSEs for both the MAP and MMSE estimators are observed to be significantly larger than the lower BB for small values for $\nu_1$ and $\nu_2$ (uninformative priors), but very close to the bound for large values of $\nu_1, \nu_2$ (informative priors).

**Performance vs. number of components.** Fig. 3 displays the MSEs of the MMSE and MAP estimators and the BBs when the number of components $R$ - thus, the number of parameter $p$ - is varied, specifically $R \in \{1, 2, 3, \dots, 10\}$. We can observe that: 1) **BB vs. R:** The BB increases with increasing values for $R$ and is very tightly approximated by aBB, thus suggesting an asymptotically linear behavior in $R$. 2) **BB vs. MSE:** The values taken by the MSE and the BB are very similar for a small number of components/parameters. For large values of $R$, the bound is slightly less tight. Since the sample size is fixed here, this behavior is coherent. Indeed, we would expect that larger sample sizes are required to converge to the BB when more parameters are estimated.

#### 3.3. Application to a multivariate multifractal analysis

Finally, the theoretical results of this paper are applied to a practical example related to multivariate multifractal analysis [4,5,17]. In particular, the model and estimation framework of [7] are considered for bivariate time series ($R = 2$). In this context, the Fourier transform of the logarithm of wavelet leaders, defined as nonlinear and nonlocal transformations of wavelet coefficients, approximately obeys the data model considered in this paper, where the elements of $\Sigma_1$ are directly related to the multifractality of the data, i.e., $\Sigma_1 = -[c_{20}, c_{11}; c_{11}, c_{02}]$ and thus $\theta_{1:3} = -(c_{20}, c_{11}, c_{02})$, where $c_{20}, c_{02} < 0$ are related to the widths
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Comparison between the sample MSEs of the MAP and MMSE estimators averaged over 1000 independent realizations versus the BB for sample sizes $N \in \{2^5, 2^6, \dots, 2^{12}\}$, $v_1 = v_2 = 80$ and $R = 2$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Comparison between the sample MSEs of the MAP and MMSE estimators averaged over 1000 independent realizations versus the BB, varying the degrees of freedom $v_1 = v_2$, for $R = 2$ and $N = 2^8$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Comparison between the sample MSEs of the MAP and MMSE estimators averaged over 1000 independent realizations versus the BB, varying the number of components $R$, for $N=2^8$, $\Omega_1=\Omega_2=I_R$ and $v_1=v_2=80$.</figcaption></figure>

of the marginal multifractal spectra, and $c_{11}$ quantifies the joint multifractality. The matrix $\Sigma_2$ is an adjustment parameter that essentially subsumes the short-lag autocorrelation of log-wavelet leaders. Moreover, one can define a multifractal correlation parameter as $\rho_{mf} = -\frac{c_{11}}{\sqrt{c_{20}c_{02}}}$ (with $-1 \le \rho_{mf} \le 1$), that quantifies a dependence beyond linear correlation between the time series. Note that the BB of $\rho_{mf}$ is obtained from the BB of $\Sigma_1$ using the functional invariance principle [1].

**Simulation study.** This section considers 2000 independent copies of $2^{10} \times 2$ time series of a canonical multifractal model process, i.e., a bivariate multifractal random walk (bMRW) [17–19], to compute the sample MSE of the MMSE estimator, and the BB, for different multifractal parameter settings, controlled by $\Sigma_1$. In a first experiment, $\Sigma_1$ is generated using $\Omega_1 = [0.5, 0; 0, \omega]$, with $0.37 \le \omega \le 1.2$ and $v_1 = 10$, leading to realistic expected values

for the multifractal parameters, i.e., $-c_{20} = 0.05$ and $-c_{02} = \{0.037, \dots, 0.12\}$. In a second experiment $\Omega_1 = [0.4, \gamma; \gamma, 0.4]$ with $\gamma$ tuned such that $0 \le \rho_{mf} \le 0.8$ in average. The parameters of $\Sigma_2$ cannot be controlled by the bMRW synthesis and are thus unknown, and we set $\Omega_2 = \Omega_1$ and $v_2 = v_1$. Results for the two experiments are presented in Fig. 4 (top and bottom row, respectively). They indicate that the derived BBs provide good indications for the variations of the observed MSE of the multifractal parameter estimates. In particular, they show that: 1) the MSE of the estimator of $c_{20}$ does not depend on $c_{02}$, which is expected since $c_{20}$ is a marginal parameter of the first data component that is independent of $c_{02}$; 2) the MSE of the estimator of $c_{02}$ increases with $c_{02}$. Indeed, $c_{02}$ controls the variance of the marginal likelihood of the second data component; 3) the MSE of the estimator of $c_{11}$ also increases with $c_{02}$ because $\rho_{mf}$ is held fixed so that $c_{11}$ and
<figure><img src="image_6.png" /><figcaption>Fig. 4. Sample MSE for multifractal parameters $c_{20}$, $c_{02}$, $c_{11}$ as a function of $c_{02}$ (top) and sample MSE for multifractal correlation $\rho_{mf}$ (bottom).</figcaption></figure>

thus its covariance also increase; 4) the MSE of $\rho_{mf}$ decreases in a non-trivial way when this parameter increases.

## 4. Conclusion

This paper derived and studied Bayesian Cramér-Rao lower bounds for the MSE of estimators of two symmetric positive definite matrices whose sum is the covariance matrix of a zero mean complex circularly-symmetric Gaussian model. To calculate these Bayesian bounds, a novel closed-form expression for a non-trivial expectation involving Wishart random matrices was provided. The properties of these Bayesian bounds were studied analytically. Various numerical simulations were used to validate the theoretical results and study the properties of the proposed Bayesian bounds. The practical interest of the bounds derived in this paper was finally illustrated for the estimation of the parameters of the bivariate multifractal spectrum.

## Declaration of Competing Interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

### CRediT authorship contribution statement

**Lorena Leon:** Conceptualization, Writing – review & editing.
**Herwig Wendt:** Conceptualization, Writing – review & editing.
**Jean-Yves Tourneret:** Conceptualization, Writing – review & editing.

### Data availability

No data was used for the research described in the article.

### References

[1] H.L.V. Trees, Detection, Estimation, and Modulation Theory, Part I, John Wiley & Sons, Inc., 2001.

[2] O. Besson, S. Bidon, J.-Y. Tourneret, Bounds for estimation of covariance matrices from heterogeneous samples, IEEE Trans. Signal Process. 56 (7) (2008) 3357–3362.

[3] S. Jaffard, P. Abry, H. Wendt, Irregularities and scaling in signal and image processing: multifractal analysis, World scientific publishing, Singapore, 2015, pp. 31–116.

[4] S. Jaffard, S. Seuret, H. Wendt, R. Leonarduzzi, S. Roux, P. Abry, Multivariate multifractal analysis, Appl. Comput. Harmon. Anal. 46 (3) (2019) 653–663.

[5] S. Jaffard, S. Seuret, H. Wendt, R. Leonarduzzi, P. Abry, Multifractal formalisms for multivariate analysis, Proc. R. Soc. A 475 (2229) (2019).

[6] H. Wendt, S. Combrexelle, Y. Altmann, J.-Y. Tourneret, S. McLaughlin, P. Abry, Multifractal analysis of multivariate images using gamma Markov random field priors, SIAM J. Imag. Sci. (SIIMS) 11 (2) (2018) 1294–1316.

[7] L. Len, H. Wendt, J.-Y. Tourneret, P. Abry, Bayesian estimation for the parameters of the bivariate multifractal spectrum, in: Proc. European Signal Process. Conf. (EUSIPCO), 2021, pp. 1930–1934. Dublin, Ireland

[8] L. Leon, H. Wendt, J.-Y. Tourneret, P. Abry, A Bayesian framework for multivariate multifractal analysis, IEEE Trans. Signal Process. 70 (2022) 3663–3675.

[9] D.K. Duvenaud, H. Nickisch, C. Rasmussen, Additive Gaussian processes, Adv. Neural Inf. Process. Syst. 24 (2011).

[10] J. Hensman, N. Durrande, A. Solin, Variational fourier features for gaussian processes, J. Mach. Learn. Res. 18 (1) (2017) 5537–5588.

[11] L. Cheng, S. Ramchandran, T. Vatanen, N. Lietzén, R. Lahesmaa, A. Vehtari, H. Lähdesmäki, An additive Gaussian process regression model for interpretable non-parametric analysis of longitudinal data, Nat. Commun. 10 (1) (2019) 1–11.

[12] P. Graczyk, G. Letac, H. Massam, The hyperoctahedral group, symmetric group representations and the moments of the real wishart distribution, J. Theor. Probab. 18 (1) (2005) 1–42.

[13] T. Holgersson, J. Pielaszkiewicz, A collection of moments of the Wishart distribution, in: Recent Developments in Multivariate and Random Matrix Analysis, Springer International Publishing, 2020, pp. 147–162.

[14] R. Frehlich, Cramér-rao bound for Gaussian random processes and applications to radar processing of atmospheric signals, IEEE Trans. Geosci. Remote Sens. 31 (6) (1993) 1123–1131.

[15] B. Porat, B. Friedlander, Computation of the exact information matrix of Gaussian time series with stationary random components, IEEE Trans. Acoust., Speech, Signal Process. 34 (1) (1986) 118–130.

[16] C.P. Robert, G. Casella, Monte Carlo statistical methods, Springer, New York, USA, 2005.

[17] H. Wendt, R. Leonarduzzi, P. Abry, S. Roux, S. Jaffard, S. Seuret, Assessing cross-dependencies using bivariate multifractal analysis, in: Proc. IEEE Int. Conf. Acoust., Speech, and Signal Proces. (ICASSP), 2018. Calgary, Canada

[18] E. Bacry, J. Delour, J. Muzy, Multifractal random walk, Phys. Rev. E 64 (2001).

[19] J.F. Muzy, J. Delour, E. Bacry, Modelling fluctuations of financial time series: from cascade process to stochastic volatility model, Eur. Phys. J. B 17 (2000).