How to build a Cross-correlated Covariance matrix by solving an equation with Covariance and Variance expression of an unknown random variable?

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 :

$$\dfrac{1}{\sigma_{\hat{\tau}}^{2}}=\dfrac{1}{\sigma_1^2}+\dfrac{1}{\sigma_2^2}\quad(1)$$

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 :

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

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

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

$$= \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 :

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

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

$$2\,\text{Cov}(X,Y)+a^2\,\text{Var}(X)+b^2\,\text{Var}(Y)+2ab\,\text{Cov}(X,Y)+2\,\text{Cov}(X,aX+bY)+2\,\text{Cov}(Y,aX+bY)$$

$$=2\,\text{Cov}(X,Y)+a^2\,\text{Var}(X)+b^2\,\text{Var}(Y)+2ab\,\text{Cov}(X,Y)+2\,\text{Cov}(X,aX)+2\,\text{Cov}(X,bY)$$

$$=2\,\text{Cov}(X,Y)+a^2\,\text{Var}(X)+b^2\,\text{Var}(Y)+2ab\,\text{Cov}(X,Y)+2a\,\text{Var}(X)+2b\,\text{Cov}(X,Y)+2a\,\text{Cov}(X,Y)+2b\,\text{Var}(Y)=0\quad(2)$$

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.

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

• @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 Nov 25 '20 at 3:01
• $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$. Nov 25 '20 at 7:50
• 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 Nov 25 '20 at 9:37
• 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 ? Nov 25 '20 at 9:55
• Here is the Julia code for the example in the answer. gofile.io/d/AfK04v Nov 25 '20 at 12:17