Quick way to find a common basis of eigenvectors between 2 matrices : valid or not?

Question

Following the advise of @Federico Polonion a previous post, one suggested, to find a basis of common eigen vectors between 2 matrices, to simply do :

Generate 2 random scalars x and y
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_Bbut 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)

Why downvote ? Initially, this question has been suggested by someone other than me. At first sight, I would have tought it was deserved to be asked, I mean, I am just looking for a solution, so this downvote is not justified. — youpilat13, Jan 18 at 19:14
If the two matrices are dense, it suffices to transform the second matrix into the eigensystem of the first and check whether it becomes diagonal. Then there is no need to disentagle the x-y multiplied eigenvalues, which might be tricky. — davidhigh, Jan 19 at 6:56
@davidhigh And what if it doesn't become diagonal? If A and B have multiple eigenvalues, then this problem is non-trivial. For instance, take $A = V diag(1,1,2,2) V^{-1}$ and $B = V diag(1,2,1,2) V^{-1}$ with a random-generated $V$, and try to recover the common eigenvectors. — Federico Poloni, Jan 19 at 9:12
@FredericoPoloni: Ok, I misunderstood the question, but then the basic idea is to expand the second matrix into eigenspaces of the first, and diagonalize again. Which, however, might not always be convenient, particularly if the matrices are sparse. I've used the proposed method to find eigenspaces corresponding to L-S symmetry in quantum mechanics, and it indeed simplified things, but I also disliked the work-around character, where it's not easy to obtain the eigenvalue-pairs. — davidhigh, Jan 19 at 15:01
@youpilat13 Sorry, I don't have the time nor the expertise to debug your code. — Federico Poloni, Jan 19 at 21:38

Federico Poloni · Accepted Answer · 2021-01-18 20:50:58Z

Yes, the eigenvectors found with this method may depend on $x$ and $y$, but no, it doesn't matter in practice.

If $A$ and $B$ share a basis of common eigenvectors, then $$ A = V\operatorname{diag}(\alpha_1,\dots,\alpha_n)V^{-1}, \quad B = V\operatorname{diag}(\beta_1,\dots,\beta_n)V^{-1}. $$ If $A$ has all distinct eigenvalues, then $V$ is unique, up to taking multiples of each column, so eig(A) returns that basis. Problems ensue if $A$ has repeated eigenvalues, though, for instance, $\alpha_i=\alpha_j$ for two indices $i,j$, because then eig can return an arbitrary basis of $span(v_i,v_j)$ and nothing guarantees that it will be also a basis of eigenvectors for $B$: indeed, if $\beta_i=\beta_j$ then any nontrivial linear combination of $v_i,v_j$ is not an eigenvalue of $B$. As an extreme case, think about $A=I$, but any case with repeated eigenvalues will cause non-uniqueness. The same can be said, of course, about eig(B).

So, a workaround is computing an eigenbasis for $$ C = xA + yB = V\operatorname{diag}(x\alpha_1 + y\beta_1, \dots, x\alpha_n + y\beta_n)V^{-1}. $$ Even if both $A$ and $B$ have repeated eigenvalues, the eigenvalues of $C$ are distinct for almost any choice of $x$ and $y$, unless $\alpha_i=\beta_i$ and $\alpha_j=\beta_j$ for a pair of indices $i,j$ (or more than two of them, or more than one pair). Hence ties are broken, and $V$ is the unique solution to eig(C), up to column multiples.

If $\alpha_i=\beta_i$ and $\alpha_j=\beta_j$ for a pair of indices, then any linear combination of $v_i$ and $v_j$ is an eigenvector of $C$, and eig(C) may return an arbitrary basis of $span(v_i, v_j)$. But then those are also eigenvalues of $A$ and $B$, so it doesn't matter!

This analysis is just sketched, and one should speak of eigenspaces instead of pairs of indices to be more general, but I hope this makes things clearer.

Could you take a look please at my UPDATE to see what I have scripted with your method, rather from what I have understood about your method. — youpilat13, Jan 19 at 12:40
@youpilat13 That seems about right, at a first glance. Of course, the best way to confirm it is checking in a few examples if V actually diagonalizes both A and B. :) — Federico Poloni, Jan 19 at 12:43
Unfortunately, as I said in my UPDATE 2, V vectors don't diagonalize A and B matrices, at least from my last tests. — youpilat13, Jan 19 at 21:38
@youpilat13 Do your A and B have a basis of common eigenvectors? Have you tried some artificially generated examples that do? — Federico Poloni, Jan 19 at 21:39
If I understand correctly what you are asking, there is a ton of literature on this problem, known as joint diagonalization. There are methods for finding an orthogonal eigenvector basis (which is the easier case) as well as for the non-orthogonal case. Here is one paper: pdfs.semanticscholar.org/56e0/…. — Amit Hochman, Jan 23 at 19:46

Stack Exchange Network

current community

your communities

more stack exchange communities

Quick way to find a common basis of eigenvectors between 2 matrices : valid or not?

1 Answer 1

Your Answer

Not the answer you're looking for? Browse other questions tagged matrix eigenvalues matrix-factorization basis-set randomized-algorithms or ask your own question.

Linked

Hot Network Questions

Quick way to find a common basis of eigenvectors between 2 matrices : valid or not?

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged matrix eigenvalues matrix-factorization basis-set randomized-algorithms or ask your own question.

Linked

Related

Hot Network Questions