Following the advise of @Federico Poloni
on a previous post, one suggested, to find a basis of common eigen vectors between 2 matrices, to simply do :
Generate 2 random scalars
x
andy
Compute the eigen vectors and values of the quantity
x*A + b*Y
by doing in Matlab :eig(x * A + y * B)
But I would have thought that eigen vectors obtained would be surely dependent of x
and y
scalar random values, wouldn't they ?
Any clarification is welcome, I have serious doubts about this method, what do you think about ?
UPDATE 1:
@Federico Polini
. Thanks for your suggestion. Here is below how I understood the way to implement it. Surely some improvements have to be done, especially the random process that produces differents constraints from the 2 combined Fisher matrix :
import os,sys
import numpy as np
# seed the pseudorandom number generator
from numpy.random import seed
from numpy.random import rand
N = 7
Nsq2 = 2*N*N
# Load spectro and WL+GCph+XC
FISH_GCsp = np.loadtxt('Fisher_GCsp_flat.txt')
FISH_XC = np.loadtxt('Fisher_XC_GCph_WL_flat.txt')
# Marginalizing over uncommon parameters between the two matrices
COV_GCsp_first = np.linalg.inv(FISH_GCsp)
COV_XC_first = np.linalg.inv(FISH_XC)
COV_GCsp = COV_GCsp_first[0:N,0:N]
COV_XC = COV_XC_first[0:N,0:N]
# Invert to get Fisher matrix
FISH_sp = np.linalg.inv(COV_GCsp)
FISH_xc = np.linalg.inv(COV_XC)
# Search for common build eigen vectors between FISH_sp and FISH_xc
D1,V1 = np.linalg.eig(FISH_sp)
D2,V2 = np.linalg.eig(FISH_xc)
# Loop to generate a basis of common eigen vectors
V = np.zeros((7,7))
D = np.zeros((7,7))
# seed random number generator : based on current system time
seed(None)
# generate some random numbers
coef = rand(2)
C = coef[0]*FISH_sp + coef[1]*FISH_xc
D, V = np.linalg.eig(C)
# Check if we get diagonal matrix
eigen_final = np.diag(coef[0]*D1 + coef[1]*D2)
#eigen_final = D
# Final sum
FISH_final = np.dot(np.dot(V,eigen_final),np.linalg.inv(V))
# Save Fisher_final
np.savetxt('Fisher_final.txt',FISH_final)
# Print constraints and FoM
status = os.system('./printFoM_and_Constraints_from_FisherTotallySAF Fisher_final.txt OPT GCph 0 0 F N')
The following Archive contains all the file necessary to execute :
$ python compute_common_eigen_vectors_basis.py
But I don't get the same constraints for multiple executions, the FoM (figure of merit) is also always changing.
The constraints can vary from a factor 0.5 up to 2-3 between the multiple executions (roughly for FoM between 700 and 2000 whereas ~ 1400-1700 is expected)
What can I do to have fixed values, even if I do multiple executions ?
I used seed(None)
to have a really different random but I wonder if it is a good idea, what do you think about this ?
Maybe I have made a mistake in my Python script above.
OR maybe I have to do a lot of executions and take the average of all constraints and FoM ? This would be interesting maybe ...
UPDATE 2: I tried to check if the eigen vectors come from V = eig(x * A + y * B)
could check the relation :
Diagonal_Matrix_A = V^-1 A V
Diagonal_Matrix_B = V^-1 B V
I talk here about only one random generation of x
and y
with numpy.random
library functions, that is, 2 randoms numbers generated comprised between 0 and 1 and following a flat PDF (if I don't say nonsense things, correct me if necessary).
So I decided to generate a lot of random pair (x,y)
, for example for a large sample of 1000 pairs (x,y)
and make a synthesis by computing the average elements of matrix V
eigen vectors used in the relation (see also code snippet above) :
D, V = np.linalg.eig(C)
# Check if we get diagonal matrix
eigen_final = np.diag(coef[0]*D1 + coef[1]*D2)
#eigen_final = D
# Final sum
FISH_final = np.dot(np.dot(V,eigen_final),np.linalg.inv(V))
But the results are not as pertinent as expected. Indeed, like I said, I get the following average constraints :
Omega_m_mean = 0.008009373151819661
Omega_b_mean = 0.0021941755866995267
w0_mean = 0.05820248972906542
wa_mean = 0.0316912395233744
h_mean = 0.004471179162394616
ns_mean = 0.012590924583566296
sig8_mean = 0.005018408809039552
with a mean V
passing matrix called V_mean
that is expected to produce, with A
and B
matrices, diagonal matrices by the relation :
Diagonal_Matrix_A = V^-1 A V
Diagonal_Matrix_B = V^-1 B V
But unfortunately, this is not the case, the 2 matrices Diagonal_Matrix_A
and Diagonal_Matrix_B
are not diagonal :
For A
matrix (FISH_sp
) :
Checking diagonal : inv(V_mean)*FISH_sp*V_mean
[[ 2.12855597e+06 1.95764387e+06 1.84943650e+05 -1.66914168e+05
-2.23929461e+05 -1.51995849e+04 -6.38824270e+01]
[ 1.95764387e+06 2.03080752e+06 2.61872127e+05 -5.37480105e+03
-2.51793418e+05 -4.88211121e+03 6.60750409e+01]
[ 1.84943650e+05 2.61872127e+05 1.92233641e+06 7.96864527e+04
-8.73830084e+03 -1.37696109e+03 -1.58946314e+02]
[-1.66914168e+05 -5.37480105e+03 7.96864527e+04 1.60081106e+05
-1.52294236e+04 9.51239018e+03 4.79120267e+02]
[-2.23929461e+05 -2.51793418e+05 -8.73830084e+03 -1.52294236e+04
4.33096664e+04 -1.61417441e+03 2.41535313e+02]
[-1.51995849e+04 -4.88211121e+03 -1.37696109e+03 9.51239018e+03
-1.61417441e+03 1.39146896e+03 -4.63011124e+01]
[-6.38824270e+01 6.60750409e+01 -1.58946314e+02 4.79120267e+02
2.41535313e+02 -4.63011124e+01 1.87869422e+01]]
and for B
(FISH_xc
) :
Checking diagonal : inv(V)*FISH_xc*V
[[ 2.92703219e+07 -2.88904692e+06 -2.72935692e+05 2.46328187e+05
3.30470077e+05 2.24312066e+04 9.42762532e+01]
[-2.88904692e+06 4.09024426e+06 -3.86465013e+05 7.93201086e+03
3.71591080e+05 7.20491025e+03 -9.75120638e+01]
[-2.72935692e+05 -3.86465013e+05 8.85586342e+04 -1.17599480e+05
1.28957884e+04 2.03208830e+03 2.34569406e+02]
[ 2.46328187e+05 7.93201086e+03 -1.17599480e+05 5.40042526e+05
2.24752419e+04 -1.40381721e+04 -7.07074946e+02]
[ 3.30470077e+05 3.71591080e+05 1.28957884e+04 2.24752419e+04
4.18669716e+05 2.38216240e+03 -3.56452399e+02]
[ 2.24312066e+04 7.20491025e+03 2.03208830e+03 -1.40381721e+04
2.38216240e+03 1.23971873e+04 6.83301435e+01]
[ 9.42762532e+01 -9.75120638e+01 2.34569406e+02 -7.07074946e+02
-3.56452399e+02 6.83301435e+01 9.04483695e+01]]
I would have thought that the fact of generating a large sample of coefficient pairs and take the average for all elements of passing matrix V
would allow to gain in precision the determination of the common eigen vectors basis between A
and B
but the passing matrix represented by the matrix V
is not a valid common eigen vectors basis, by checking if V
checks the relation of diagonal matrices production (what I expected for Diagonal_Matrix_A
and Diagonal_Matrix_B
but these 2 matrices are not diagonal).
Moreover, the mean Figure of Merit (FoM) has a too low value on this large sample : FoM = 1011.59
whereas we expect rather a FoM around 1400-1800
(coming from other approximative studies of probes combination). Indeed, the single matrix B (FISH_xc) has a FoM = 1030
; so by combining FISH_xc
with FISH_sp (the spectro Fisher)
must make increase more the final FoM up to 1400-1800
as I said before. The FoM of FISH_xc
is about 55
.
Now, which strategy could anyone suggest me to find this such wanted common eigen vectors basis ?
Remark: If I multiply, for example by 10, the random of pairs (x,y)
, like this :
# generate some random numbers
coef = 10*rand(2)
# print('coef = ', coef)
C = coef[0]*FISH_sp + coef[1]*FISH_xc
D, V = np.linalg.eig(C)
Then, the mean FoM is too big (order of FoM = 10000
). Finally, how to adjust or find the appropriate coefficient in front of rand(2)
in order to get consistent constraints and FoM ?
I am opened to any kind of improvements, provided the method has a way to be justified rationnaly.
PS: Here the version of the last code used in UPDATE 1 to compute means (FoM and constraints) from a large sample (1000) of generated pairs (x,y)
:
import os,sys
import numpy as np
# seed the pseudorandom number generator
from numpy.random import seed
from numpy.random import rand
# Capture output of printFom
import subprocess
N = 7
Nsq2 = 2*N*N
# Load spectro and WL+GCph+XC
FISH_GCsp = np.loadtxt('Fisher_GCsp_flat.txt')
FISH_XC = np.loadtxt('Fisher_XC_GCph_WL_flat.txt')
# Marginalizing over uncommon parameters between the two matrices
COV_GCsp_first = np.linalg.inv(FISH_GCsp)
COV_XC_first = np.linalg.inv(FISH_XC)
COV_GCsp = COV_GCsp_first[0:N,0:N]
COV_XC = COV_XC_first[0:N,0:N]
# Invert to get Fisher matrix
FISH_sp = np.linalg.inv(COV_GCsp)
FISH_xc = np.linalg.inv(COV_XC)
# Search for common build eigen vectors between FISH_sp and FISH_xc
D1,V1 = np.linalg.eig(FISH_sp)
D2,V2 = np.linalg.eig(FISH_xc)
# Loop to generate a basis of common eigen vectors
sizeSample = 1000
V = np.zeros((7,7,sizeSample))
V_sum = np.zeros((7,7))
V_mean = np.zeros((7,7))
D = np.zeros((7,7,sizeSample))
# FoM array
FoM = np.zeros(sizeSample)
# Constraints array
Omega_m = np.zeros(sizeSample)
Omega_b = np.zeros(sizeSample)
w0 = np.zeros(sizeSample)
wa = np.zeros(sizeSample)
h = np.zeros(sizeSample)
ns = np.zeros(sizeSample)
sig8 = np.zeros(sizeSample)
for i in range(sizeSample):
print(i)
# seed random number generator : based on current system time
seed(None)
# generate some random numbers
coef = rand(2)
# print('coef = ', coef)
C = coef[0]*FISH_sp + coef[1]*FISH_xc
D, V = np.linalg.eig(C)
#print('eigen values of first vector = ', D)
#print('Checking diagonal : inv(V)*FISH_sp*V')
#print(np.dot(np.dot(np.linalg.inv(V),FISH_sp),V))
#print('Checking diagonal : inv(V)*FISH_xc*V')
#print(np.dot(np.dot(np.linalg.inv(V),FISH_xc),V))
# Check if we get diagonal matrix
eigen_final = np.diag(coef[0]*D1 + coef[1]*D2)
# Final sum
FISH_final = np.dot(np.dot(V,eigen_final),np.linalg.inv(V))
# Save Fisher_final
np.savetxt('Fisher_final.txt',FISH_final)
# Print constraints and FoM
#status = os.system('/Users/henry/bin/printFoM_and_Constraints_from_FisherTotallySAF Fisher_final.txt OPT GCph 0 0 F N')
FoM[i] = subprocess.Popen("/Users/henry/bin/printFoM_and_Constraints_from_FisherTotallySAF Fisher_final.txt OPT GCph 0 0 F N | grep FoM | awk '{print $3}'",
shell=True,
stdout=subprocess.PIPE,
universal_newlines=True).communicate()[0]
Omega_m[i] = subprocess.Popen("/Users/henry/bin/printFoM_and_Constraints_from_FisherTotallySAF Fisher_final.txt OPT GCph 0 0 F N | grep -v FoM | grep wm | awk '{print $3}'",
shell=True,
stdout=subprocess.PIPE,
universal_newlines=True).communicate()[0]
Omega_b[i] = subprocess.Popen("/Users/henry/bin/printFoM_and_Constraints_from_FisherTotallySAF Fisher_final.txt OPT GCph 0 0 F N | grep -v FoM | grep wb |awk '{print $3}'",
shell=True,
stdout=subprocess.PIPE,
universal_newlines=True).communicate()[0]
w0[i] = subprocess.Popen("/Users/henry/bin/printFoM_and_Constraints_from_FisherTotallySAF Fisher_final.txt OPT GCph 0 0 F N | grep -v FoM | grep w0 | awk '{print $3}'",
shell=True,
stdout=subprocess.PIPE,
universal_newlines=True).communicate()[0]
wa[i] = subprocess.Popen("/Users/henry/bin/printFoM_and_Constraints_from_FisherTotallySAF Fisher_final.txt OPT GCph 0 0 F N | grep -v FoM | grep wa | awk '{print $3}'",
shell=True,
stdout=subprocess.PIPE,
universal_newlines=True).communicate()[0]
h[i] = subprocess.Popen("/Users/henry/bin/printFoM_and_Constraints_from_FisherTotallySAF Fisher_final.txt OPT GCph 0 0 F N | grep -v FoM | grep -v ph | grep h | awk '{print $3}'",
shell=True,
stdout=subprocess.PIPE,
universal_newlines=True).communicate()[0]
ns[i] = subprocess.Popen("/Users/henry/bin/printFoM_and_Constraints_from_FisherTotallySAF Fisher_final.txt OPT GCph 0 0 F N | grep -v FoM | grep ns | awk '{print $3}'",
shell=True,
stdout=subprocess.PIPE,
universal_newlines=True).communicate()[0]
sig8[i] = subprocess.Popen("/Users/henry/bin/printFoM_and_Constraints_from_FisherTotallySAF Fisher_final.txt OPT GCph 0 0 F N | grep -v FoM | grep s8 | awk '{print $3}'",
shell=True,
stdout=subprocess.PIPE,
universal_newlines=True).communicate()[0]
FoM_sum = 0
Omega_m_sum = 0
Omega_b_sum = 0
w0_sum = 0
wa_sum = 0
h_sum = 0
ns_sum = 0
sig8_sum = 0
## Mean of all constraints and FoM
for i in range(sizeSample):
FoM_sum += FoM[i]
Omega_m_sum += Omega_m[i]
Omega_b_sum += Omega_b[i]
w0_sum += w0[i]
wa_sum += wa[i]
h_sum += h[i]
ns_sum += ns[i]
sig8_sum += sig8[i]
for j in range(7):
for k in range(7):
V_sum[j,k] += V[j,k]
FoM_mean = FoM_sum/sizeSample
Omega_m_mean = Omega_m_sum/sizeSample
Omega_b_mean = Omega_b_sum/sizeSample
w0_mean = w0_sum/sizeSample
wa_mean = wa_sum/sizeSample
h_mean = h_sum/sizeSample
ns_mean = ns_sum/sizeSample
sig8_mean = sig8_sum/sizeSample
V_mean = V_sum/sizeSample
print('eigen values of first vector = ', D)
print('Checking diagonal : inv(V_mean)*FISH_sp*V_mean')
print(np.dot(np.dot(np.linalg.inv(V_mean),FISH_sp),V_mean))
print('Checking diagonal : inv(V)*FISH_xc*V')
print(np.dot(np.dot(np.linalg.inv(V_mean),FISH_xc),V_mean))
print('')
print('FoM Mean = ', FoM_mean)
print('')
print('Omega_m_mean = ', Omega_m_sum/sizeSample)
print('Omega_b_mean = ', Omega_b_sum/sizeSample)
print('w0_mean = ', w0_sum/sizeSample)
print('wa_mean = ', wa_sum/sizeSample)
print('h_mean = ', h_sum/sizeSample)
print('ns_mean = ', ns_sum/sizeSample)
print('sig8_mean = ', sig8_sum/sizeSample)