I am using the qndiag library to try to find a diagonalisation for 2 given matrices.
The github is here : qndiag libray
The function qndiag is defined like this (not entirely source) :
def qndiag(C, B0=None, weights=None, max_iter=1000, tol=1e-6,
lambda_min=1e-4, max_ls_tries=10, diag_only=False,
return_B_list=False, verbose=False):
"""Joint diagonalization of matrices using the quasi-Newton method
Parameters
----------
C : array-like, shape (n_samples, n_features, n_features)
Set of matrices to be jointly diagonalized. C[0] is the first matrix,
etc...
B0 : None | array-like, shape (n_features, n_features)
Initial point for the algorithm. If None, a whitener is used.
weights : None | array-like, shape (n_samples,)
Weights for each matrix in the loss:
L = sum(weights * KL(C, C')) / sum(weights).
No weighting (weights = 1) by default.
max_iter : int, optional
Maximum number of iterations to perform.
tol : float, optional
A positive scalar giving the tolerance at which the
algorithm is considered to have converged. The algorithm stops when
|gradient| < tol.
lambda_min : float, optional
A positive regularization scalar. Each eigenvalue of the Hessian
approximation below lambda_min is set to lambda_min.
max_ls_tries : int, optional
Maximum number of line-search tries to perform.
diag_only : bool, optional
If true, the line search is done by computing only the diagonals of the
dataset. The dataset is then computed after the line search.
Taking diag_only = True might be faster than diag_only=False
when the matrices are large (n_features > 200)
return_B_list : bool, optional
Chooses whether or not to return the list of iterates.
verbose : bool, optional
Prints informations about the state of the algorithm if True.
Returns
-------
D : array-like, shape (n_samples, n_features, n_features)
Set of matrices jointly diagonalized
B : array, shape (n_features, n_features)
Estimated joint diagonalizer matrix.
infos : dict
Dictionnary of monitoring informations, containing the times,
gradient norms and objective values.
References
----------
P. Ablin, J.F. Cardoso and A. Gramfort. Beyond Pham's algorithm
for joint diagonalization. Proc. ESANN 2019.
https://www.elen.ucl.ac.be/Proceedings/esann/esannpdf/es2019-119.pdf
https://hal.archives-ouvertes.fr/hal-01936887v1
https://arxiv.org/abs/1811.11433
"""
t0 = time()
n_samples, n_features, _ = C.shape
if B0 is None:
C_mean = np.mean(C, axis=0)
d, p = np.linalg.eigh(C_mean)
B = p.T / np.sqrt(d[:, None])
else:
B = B0
if weights is not None: # normalize
weights_ = weights / np.mean(weights)
else:
weights_ = None
D = transform_set(B, C)
I am using this Python script to compute these 2 diagonalisation as closed as possible :
import os, sys
import numpy as np
from qndiag import qndiag
# dimension
m=7
# number of matrices
n=2
# 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:m,0:m]
COV_XC = COV_XC_first[0:m,0:m]
# Invert to get Fisher matrix
FISH_sp = np.linalg.inv(COV_GCsp)
FISH_xc = np.linalg.inv(COV_XC)
# Drawing a random set of commuting matrices
C=np.zeros((n,m,m));
B=np.zeros((m,m));
C[0] = np.array(FISH_sp)
C[1] = np.array(FISH_xc)
# Perform operation of diagonalisation
[D, B] = qndiag(C, None, None, 1000, 1e-3);
D0 = np.array(D[0])
D1 = np.array(D[1])
print(D0)
print(D1)
print(B)
# Print diagonal matrices
M0 = np.dot(np.dot(B.T,C[0]),B)
M1 = np.dot(np.dot(B.T,C[1]),B)
print(M0)
print(M1)
Given my 2 matrices 7x7 FISH_sp and FISH_xc, I get an error with printing at the end M0 and M1 of kind :
{'t_list': [0.00012111663818359375, 0.00034308433532714844, 0.0004680156707763672, 0.0005850791931152344, 0.0007319450378417969, 0.0008790493011474609, 0.00098419189453125, 0.0010869503021240234, 0.0011870861053466797, 0.0012888908386230469, 0.0013890266418457031, 0.0014889240264892578, 0.0015878677368164062, 0.0016870498657226562, 0.0017862319946289062], 'gradient_list': [2.435835480314046, 8.032854098264083, 13.226556048661022, 9.681075100894695, 6.477682875227688, 3.1869761663221587, 2.0459590467438877, 7.102415981965997, 9.580245771870109, 3.4537238605601552, 3.9813687469559267, 2.1137748034714305, 1.3730779100371122, 0.04799779556789997], 'loss_list': [40.08624519813238, 39.65401446920329, 39.298969010821644, 38.83363862937428, 38.557138257558975, 38.3655948952275, 38.36418356814169, 38.165628179855645, 37.82921628860782, 37.80456354387957, 37.71472965598052, 37.641983813016495, 37.63102124815874, 37.630980901887284]}
Traceback (most recent call last):
File "compute_joint_diagonalization.py", line 38, in <module>
M0 = np.dot(np.dot(B.T,C[0]),B)
AttributeError: 'dict' object has no attribute 'T'
Actually this seems about the operator transpose "T" which can't be applied on a dict.
If it is that, how to convert this list into numpy arrays in order to make matricial product ?
Update
In order to check the validity of the function qndiag, I have tried to find again approximative diagonal matrices.
For this, I did :
# Print diagonal matrices
M0 = np.dot(np.dot(B.T,C[0]),B)
M1 = np.dot(np.dot(B.T,C[1]),B)
That is, following the formula from Wikipedia , with here O=B and D the diagonal matrices.
M0 = O D O^T eq(1)
which gives
D = O^T M0 O eq(2)
but this last formula eq(2) doesn't seem to be correct, I mean at the level of constraints I get with it (I am working with Fisher formalism).
Have I used the right formula? And why only eq(1) gives relative good constraints? The most logical would be the formula eq(2) and not eq(1).
question from:
https://stackoverflow.com/questions/65877295/attributeerror-dict-object-has-no-attribute-t-t-for-transposition 
