问题描述
我正在使用 qndiag 库来尝试找到 2 个给定矩阵的对角线。
github 在这里：qndiag libray
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,)
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,n_features)
Set of matrices jointly diagonalized
B : array,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,_ = 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)
我正在使用这个 Python 脚本来计算这 2 个尽可能封闭的对角线：
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,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)
给定我的 2 个矩阵 7x7 FISH_sp 和 FISH_xc，我在最后打印 M0 和 M1 时出现错误：
{'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,B)
AttributeError: 'dict' object has no attribute 'T'
实际上这似乎是操作符转置“T”，它不能应用于字典。
如果是这样，如何将此列表转换为 numpy 数组以制作矩阵产品？
更新
为了检查函数qndiag的有效性，我再次尝试找到近似对角矩阵。
为此，我做到了：
# Print diagonal matrices
M0 = np.dot(np.dot(B.T,B)
也就是说，遵循维基百科的公式，这里 O=B 和 D 是对角矩阵。
M0 = O D O^T eq(1)
给出
D = O^T M0 O eq(2)
但最后一个公式 eq(2) 似乎不正确，我的意思是在我得到的约束级别（我正在使用 Fisher 形式主义）。
我是否使用了正确的公式？为什么只有 eq(1) 给出了相对较好的约束？最合乎逻辑的是公式 eq(2) 而不是 eq(1)。
解决方法
如 toy example 中所述，如果您更改此行，您应该能够运行您的代码
[D,B] = qndiag(C,None,1000,1e-3);
到
B,_ = qndiag(C,1e-3);
并删除
D0 = np.array(D[0])
D1 = np.array(D[1])
print(D0)
print(D1)