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 ?


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$
    – youpilat13
    Nov 25 '20 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
    Nov 25 '20 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$
    – youpilat13
    Nov 25 '20 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$
    – youpilat13
    Nov 25 '20 at 9:55
  • 1
    $\begingroup$ Here is the Julia code for the example in the answer. gofile.io/d/AfK04v $\endgroup$
    – balaji
    Nov 25 '20 at 12:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.