Python/Matlab - fsolve in Fisher context : Sensitivity to guess

Question

I have 2 Fisher matrices F1 and F2. I can diagonalise these 2 matrices which give D1 and D2 as diagonal matrices, P1 and P2 the respective passing matrices.

The eigen values are respectively eigen_1 and eigen_2.

For the moment, if I want to make the combination of these 2 Fisher informations (called also "cross-correlation), I can't build a "combinated" Fisher matrix directly by summing the 2 diagonal matrices since the linear combination of random variables is different between the 2 Fisher matrices after diagonalisation.

I have eigen vectors represented by P1 (with D1 diagonal) and P2 matrices (with D2 diagonal matrix).

That's why I think that I could perform a "global" combination of eigen vectors where I can respect the MLE (Maximum Likelihood Estimator on standard deviation) as each eigen value :

1/sigma_tau^2 = 1/sigma_1^2 + 1/sigma_2^2

because sigma_tau corresponds to the best estimator constraint from MLE method.

So, I thought a convenient linear combination of each eigen vectors P1 and P2 that could allow to achieve it would be under a new matrix P whose each column represents a new eigen global vector like this :

P = aP1 + bP2

I am working considering a and b as matrix 7x7, always with fsolve Matlab and fsolve Python. As you can guess, there are a lot of different solutions as function of the initial guess I provide for fsolve. I have not yet found a way to set the "correct" guess values for a and b:

1.Initially, I took : x0 = 0.5*ones(2*7*7,1) but I got constraints too really small (too optimistic).

Looking as a weighted sum with the formula : P = aP1 + bP2

So these initial guess values have to be chosen more accurately.

By considering a and b as matrices, my goal is to find a common basis P such that P = aP1 + bP2 and that the product P^-1 F P = D eq(1) with D a "new" diagonal matrices in which I could respect the MLE (Maximum Likelihood Estimator by taking the second derivative of $log(L) = 0$ see (eq 1))

There too, we could simply think that : Var(aP1 + bP2) = a^2 Var(P1) + b^2 Var(P2) + 2ab Covar(P1,P2) = a^2 Var(P1) + b^2 Var(P2) (eq 2) if P1 and P2 are not correlated. But you agree this relation is only valid for 2 random variables P1 and P2, and not for matrices as I want to consider it here.

So, I would like to find the equivalent of this equation (eq 2) that could define an equation potentially solved by Matlab fsolve.

Here the different attempts I did and tried at each time to use fsolve with them :

• D = (a*a)*D1 + (b*b)*D2; : that converges very quickly with a guess equal to x0 = ones(2*7*7,1); ; this formulation gives encouraging results (at least from a physical point of view) but is it correct from a mathematical point of view ? My initial motivation for this expression is like for random variables (and not necessarily right for matrices) : Var(aD1 + bD2) = a^2 Var(D1) + b^2 Var(D2) with D1and D2 uncorrelated

Ideally, maybe I am wrong, we should be able to write :

P^(-1) F1 P + P^(-1) F2 P = P^(-1) F P = D' (eq 3)

as a looking like with D1 + D2 = D (eq 4) even if (eq 3) is not equal to (eq 4). In this (eq 3), I am looking for an expression of D' as a function of unknown a and b and surely D1 and D2.

I am keeping to find this equation because if I can express a matricial equation involving P1, P2, D1, D2, F1 and F2 and from which I could build the wanted passing matrix P by finding the matrix a and b, then I could find the final Fisher matrix representing the combination of the 2 first initial Fisher matrices F1 and F2 with :

F = P D P^(-1) (eq 5)

I would qualify this Fisher matrix of "cross-correlated Fisher Matrix". After, a simple inversion of it gives the constraints on each random variables.

UPDATE 1:

1. For the definition of D matrix, I did :

   D = np.dot(np.identity(7)*np.diag(a)**2,D1) + np.dot(np.identity(7)*np.diag(b)**2,D2)

That is, I have only take the diagonal elements of a and b matrices : I recall, I have defined : P = a P1 + b P2.

2. In Python's fsolve I tried with the following guess for [a,b], i.e a flatten array [7*7*2]= 98 which is a vector containing 98 elements, first 49 elements of eigen vector P1 and secondly the 49 elements eigen vector from matrix P2.

   x0 = [eigenv_sp.ravel(), eigenv_xc.ravel()] x = fsolve(f, x0)

   I get more physical acceptable constraints (compared to x0 = 0.5*ones(2*7*7,1) case) but the contours of joint distribution are not good.

However, wit this guess, the Python's fsolve seems to converge quickly unlike to Matlab s fsolve which gives wrong results (like negative variance) with this 1x98 guess vector.

Anyone could give me advices about another guess vectors that I could take for a and b matrices (which are flattened for fsolve) ?