Firstly, I want to give you a maximum of informations and precisions about my issue. If I can't manage to get the expected results, I will launch a bounty, maybe some experts or symply people who have been already faced to a similar problem will be able to help me

I have 2 covariance matrices known $\text{Cov}_X$ and $\text{Cov}_Y$ that I want to cross-correlate.

1) For this, I have performed a diagonalisation of each Fisher matrix $F_X$ (FISH_eigen_sp_flat in code) and $F_Y$ (FISH_eigen_xc_flat in code) associated of Covariance matrices $\text{Cov}_X$ and $\text{Cov}_Y$.

So, I have 2 different linear combinations that are uncorraleted, i.e just related by eigen values ($\dfrac{1}{\sigma_i^2}$) as respect of their combination.

Then, I get the diagonal matrices $D_X$ and $D_Y$. I can't build a "global" Fisher matrix directly by summing the 2 diagonal matrices since the linear combination of random variables is different between the 2 Fisher matrices.I have eigen vectors represented by $X$ (with $D_X$ diagonal) and $Y$ matrices (with $D_Y$ diagonal matrix)

That's why I think that I could perform a "global" combination of eigen vectors where I can respect the MLE estimator as eigen value :


because $\sigma_{\hat{\tau}}$ corresponds to the best estimator from MLE method (The only inconvenient in this approach for me is that I respect it only on a linear combination of random variables, and not just one, if someone could tell me if am wrong or not about this)

So, I thought a convenient linear combination that could allow to achieve it would be under the form :

$$M = X + Y + Z$$

with $$Z =aX + bY$$

where $a$ and $b$ have to be determined and such that :


Therefore, I am looking for a way to find $Z$ (not null) that could verify :


$$= \text{Var}(X)+\text{Var}(Y)+2\,\text{Cov}(X,Y)+\text{Var}(Z)+2\,\text{Cov}(X,Z)+2\,\text{Cov}(Y,Z)=\text{Var}(X)+\text{Var}(Y)$$

So, we can remove on both sides the terms $\text{Var}(X)+\text{Var}(Y)$ to get the equation :


It follows with solution under the form $Z=a\,X+b\,Y$ :




Finally, If I fix a value for $a$, I have to solve this equation like a second order equation with $b$ as unknown parameter (below b=x=unknown) :

$$\alpha\,x^2+\beta x + \delta$$

with :

  • $\alpha=\text{Var}(Y)$

  • $\beta=2\,\text{Cov}(X,Y)+2a\,\text{Cov}(X,Y)+2\,\text{Var}(Y)$

  • $\delta=2\,\text{Cov}(X,Y)+a^2\,\text{Var}(X)+2a\,\text{Var}(X)+2a\,\text{Cov}(X,Y)$

To avoid complex solutions of this second order equation, I have to respect :

$$\beta^2-4\,\alpha\delta > 0$$

2) We could for example take $a=-1$ and try to find the b unknown like this (more precision, I think that b unknow parameter is surely a matrix or a vector : if someone could confirm it since I take the scalar value $a=1$, this would fine to tell it). This choice of $a=-1$ may seem to be arbitral, it is not only the unique solution for this issue but I have not any error concerning my issue.

Therefore, by taking $a=-1$, I have implemented a function to compute varX, varY and CovXY in Python :

# Solution Z = aX + bY ; Taking a = -1 => solve unknown b represented by b_coef :
b_coef = np.zeros((7,7))
VarX = np.zeros((7,7))
VarY = np.zeros((7,7))
CovXY = np.zeros((7,7))
# Compute covariance between vectors of matrix
def compute_Cov(A,B):
  C = np.zeros((7,7))
  for i in range(7):
    C[0:7,i]= np.mean(A[0:7,i]*B[0:7,i]) - np.mean(A[0:7,i])*np.mean(B[0:7,i])
  return C

I can this way compute the variance of X (compute_Cov(X,X)), Y (compute_Cov(Y,Y) and covariance(X,Y) (compute_Cov(X,Y)).

Diagonalisation) By diagonalising the 2 Covariance matrices (inverse of Fisher matrices FISH_sp_flat and FISH_xc_flat) representing the variances on each parameters that I want to compute the standard deviations, I think I build a linear combination of these parameters that is uncorrelated :

This way, I can sum respectively for each combination the eigen values like this :

# Get eigen values (marginalized error since we handle Covariance error)
# and eigen vectors giving matrix "P" with convention : F = P D P^-1
eigen_sp_flat, eigenv_sp_flat = np.linalg.eig(np.linalg.inv(F_X))
eigen_xc_flat, eigenv_xc_flat = np.linalg.eig(np.linalg.inv(F_Y))

# Sum of FISH_eigen sp and xc
FISH_eigen_sp_flat = np.linalg.inv(np.diag(eigen_sp_flat))
FISH_eigen_xc_flat = np.linalg.inv(np.diag(eigen_xc_flat))

# MLE method (Maximum Likelihood Estimator) : Sum of FISH_eigen sp and xc : 1/sig^2 = 1/sig1^2 + 1/sig^2
FISH_eigen_sum = FISH_eigen_sp_flat + FISH_eigen_xc_flat

3) Once linear combination of parameters and Var(X), Var(Y), Cov(X,Y) computed, I build my final covariance matrix like this :

# Sum of passing matrix X+Y+Z with X = eigenv_sp_flat, Y = eigenv_xc_flat and Z to infer
eigenv_final = eigenv_sp_flat + eigenv_xc_flat + a*eigenv_sp_flat + b_coef*eigenv_xc_flat

# Covariance and Fisher matrix final
COV_final = np.dot(eigenv_final, np.linalg.inv(FISH_eigen_sum))
COV_final = np.dot(COV_final, np.linalg.inv(eigenv_final))
FISH_final = np.linalg.inv(COV_final)

The issue is that I get a better FoM (1389) (Figure of Merit = 1/area(w0,wa) which is the inverse area of 1 C.L contour for the joint parameters (w0,wa)) compared to a simple sum betwen the 2 starting Fisher matrix (1235): I would expect a more significant improvement and some anomalies appears also in the fnal constraints (for example, for a random variable, it doesn't respect the equation(1) ( $\dfrac{1}{\sigma_{\hat{\tau}}^{2}}=\dfrac{1}{\sigma_1^2}+\dfrac{1}{\sigma_2^2}\quad(1)$), i.e I have no gain for $\sigma_{\hat{\tau}}^{2}$ compared to the smallest individual $\sigma$ ($\sigma_1$ or $\sigma_2$).

Anyone could confirm me if my approach is correct, especially the calculation of a new basis that cheks V(X+Y+Z) = Var(X)+Var(Y) which seems for me essential ?

I hope that I have been enough clear in my explanations : the issue is simple and a little tricky in the same time.

EDIT 1: I realized that the condition of orthogonality ($F.T^T\neq I_d$) on the final building of eigen vectors is not respected in my code.

So, this characteristic of orthogonality is essential if I want to check the condition $\text{Var}(X+Y+Z)=\text{Var}(X)+\text{Var}(Y)$ where $Z =ax + bY$.

Therefore, I think this condition of orthogonality adds a complementary condition on the choice (and so the computation) of coefficient $a$ and $b$ into $Z =ax + bY$.

That's why I have asked another question about this criterion of orthogonality on : Building an orthogonal basis from the sum of different orthogonal basis.

Hope this will help you

Anyone could see how to translate this condition to determine a single value for $a$ and for $b$ ?

Don't hesitate to ask further information.

Any help/fix/suggestion is welcome. Sorry if this is a little bit long to read.

EDIT 2 : I have slightly modified the function compute_Cov : Does it make sense to compute variance $\text{Var}(X)$, $\text{Var}(Y)$, and $\text{CoVar}(X,Y)$ with $X$, $Y$ covariance matrices like this ? :

# Compute covariance between vectors of matrix
def compute_Cov(A, B):
  C = np.zeros(7)
  for i in range(7):
    C[i]= np.cov(A[0:7,i]*B[0:7,i])
  return C

VarX = compute_Cov(X,X)
VarX = compute_Cov(Y,Y)
CoVarX = compute_Cov(X,Y)

I have serious doubts about what I do since the beginning, anyone could see clearer ?

  • $\begingroup$ You can see a reformulation of my problem on scicomp.stackexchange.com/questions/36373/… . I restarted from zero by trying to explain better my goal. $\endgroup$
    – user359328
    Commented Nov 25, 2020 at 16:15

1 Answer 1


You are trying to find $u, v, \text{and }w$ such that:

$$ u^2 + \sigma_X u v + \sigma_Y u w + \rho_{XY} \sigma_X \sigma_Y = 0 \\ u \geq 0 \\ -1 \leq v, w \leq +1, $$

where $u = \sigma_Z$, $v = \rho_{X,Z}$ and $w = \rho_{Y,Z}$. Once you know $u, v, \text{and }w$, it is straightforward to find $a \text{ and } b$ in $Z = aX + bY$.

This needs to be solved numerically. You might have to preclude edge cases ($\rho_{XY} = \pm 1$).

An example: $\sigma_1 = \sigma_X = 0.1, \sigma_1 = \sigma_Y = 0.2$ and $\rho_{XY}$ varies from -0.9 to 0.9. $a$ and $b$ from minimising $(u^2 + \sigma_X u v + \sigma_Y u w + \rho_{XY} \sigma_X \sigma_Y)^2$ subject to the constraints look like this. The optimal objective values are of the order of $10^{-15}$, so practically $0$.

enter image description here

  • $\begingroup$ @bataji . Thanks for your answer : how do you find the relation above : $$u^2 + \sigma_X u v + \sigma_Y u w + \rho_{XY} \sigma_X \sigma_Y = 0 \\ u \geq 0 \\ -1 \leq v, w \leq +1,$$ ? A detailled answer is welcome, Regards $\endgroup$
    – user359328
    Commented Nov 25, 2020 at 3:01
  • 1
    $\begingroup$ $u^2 + \sigma_X u v + \sigma_Y u w + \rho_{XY} \sigma_X \sigma_Y = 0$ comes from condition $\sigma(X+Y+Z)^2 = \sigma(X)^2 + \sigma(Y)^2$. $\endgroup$
    – balaji
    Commented Nov 25, 2020 at 7:50
  • $\begingroup$ Thanks for your plot. Could you give me the code that allowed you to produce this plot, whatever it is on python, Matlab or others, I will adapt myself but I really need it to have a starting point given that it's a numerical stuff. If you could, I would be very grateful. Regards $\endgroup$
    – user359328
    Commented Nov 25, 2020 at 9:37
  • $\begingroup$ Which values did you take for $u, v , w$ in your plot ? and how can I infer $a$ and $b$ from these values ? your goal is to minimize $$u^2 + \sigma_X u v + \sigma_Y u w + \rho_{XY} \sigma_X \sigma_Y = 0 \\ u \geq 0 \\ -1 \leq v, w \leq +1,$$ , isn't it ? $\endgroup$
    – user359328
    Commented Nov 25, 2020 at 9:55
  • 1
    $\begingroup$ Here is the Julia code for the example in the answer. gofile.io/d/AfK04v $\endgroup$
    – balaji
    Commented Nov 25, 2020 at 12:17

You must log in to answer this question.