MATLAB：找到一种算法来快速求逆一个大的符号变量矩阵

计算科学 matlab 矩阵 显卡 逆 矩阵方程
2021-12-19 01:53:49

我必须解决包含大量符号变量的 2 个 12x12 矩阵之间的相等性，并用它执行矩阵求逆。只有一个未知数叫做“ SIGMA_O”。

当我以 2 个 2x2 矩阵为例时，我的系统很快得到解决，反演非常直接。

现在，对于 2 个 12x12 矩阵的情况，我实际上需要先对一个 31x31 的符号变量矩阵求逆（之后我将其边缘化），因为求逆需要很多时间。

我想从我的 GPU NVIDIA 卡中受益，以更快地实现这种反转，但符号数组目前不支持 GPU 优化。

在脚本下方，您将找到 inversion 行：

COV_ALL = inv(FISH_SYM)

和整个代码：

clear;
clc;
%format long;

%parpool(64);

% 2 Fisher Matrixes symbolic : FISH_GCsp_SYM, : 1 cosmo params + 1 bias spectro put for common
%                              FISH_XC_SYM : 1 cosmo params + 2 bias photo correlated

% GCsp Fisher : 7 param cosmo and 5 bias spectro which will be summed
%FISH_GCsp_SYM = sym('sp_', [17,17], 'positive');
FISH_GCsp_SYM = load('Fisher_GCsp_FoM_52.76');
% Force symmetry for GCsp
%FISH_GCsp_SYM = tril(FISH_GCsp_SYM.') + triu(FISH_GCsp_SYM,1)

% GCph Fisher : 7 param cosmo + 3 I.A + 11 bias photo correlated 
%FISH_XC_SYM = sym('xc_', [21,21], 'positive');
FISH_XC_SYM = load('Fisher_OPT_GCph_WL_XC_F_N_XSAF_EXTENDED_FoM_1037.69');
% Force symmetry for GCph
%FISH_XC_SYM = tril(FISH_XC_SYM.') + triu(FISH_XC_SYM,1)

% Brutal Common Bias : sum of 7 cosmo param ans 5 bias spectro : FISH_ALL1 = first left matrix
FISH_ALL1 = sym('xc_', [12,12], 'positive');
% Sum cosmo
FISH_ALL1(1:7,1:7) = FISH_GCsp_SYM(1:7,1:7) + FISH_XC_SYM(1:7,1:7);
% Brutal sum of bias
FISH_ALL1(7:12,7:12) = FISH_GCsp_SYM(7:12,7:12) + FISH_XC_SYM(15:20,15:20);
% Adding cross-terms
FISH_ALL1(1:7,7:12) = FISH_GCsp_SYM(1:7,7:12) + FISH_XC_SYM(1:7,15:20);
FISH_ALL1(7:12,1:7) = FISH_GCsp_SYM(7:12,1:7) + FISH_XC_SYM(15:20,1:7);

% Adding new observable "O" terms
FISH_O_SYM = sym('o_', [5,2,2], 'positive');
%FISH_O_SYM = sym('o_', [2,2], 'positive');

% Bias fiducial from spectro
Bias_sp_fid = [1.43218954, 1.52439937, 1.63460801, 1.77880054, 1.92107801]

z_ph_fid = [ 0.9595, 1.087, 1.2395, 1.45, 1.688 ]
% General case
%z_ph_fid  = np.array([ 0.2095, 0.489, 0.619, 0.7335, 0.8445, 0.9595, 1.087, 1.2395, 1.45, 1.688, 2.15 ])
Bias_ph_fid = sqrt(z_ph_fid + 1)
% JUST TEST : CAUTION
%Bias_ph_fid = Bias_sp_fid


% Definition of sigma_o
SIGMA_O = sym('sigma_o_',[5], 'positive');

for i=1:5
  FISH_O_SYM(i,1,1) = 1/SIGMA_O(i)^2*4*Bias_sp_fid(i)^2/Bias_ph_fid(i)^4
  FISH_O_SYM(i,2,2) = 1/SIGMA_O(i)^2*4*Bias_sp_fid(i)^4/Bias_ph_fid(i)^6
  FISH_O_SYM(i,1,2) = 1/SIGMA_O(i)^2*(-4*Bias_sp_fid(i)^3/Bias_ph_fid(i)^5)
  FISH_O_SYM(i,2,1) = FISH_O_SYM(i,1,2)
end

% Force symmetry
%FISH_O_SYM = (tril(FISH_O_SYM.') + triu(FISH_O_SYM,1))
FISH_O_SYM

FISH_SYM = zeros(31,31,'sym');
FISH_BIG_GCsp = zeros(31,31,'sym');
FISH_BIG_XC = zeros(31,31,'sym');

% Block bias spectro + pshot and correlations;
FISH_BIG_GCsp(1:7,1:7) = FISH_GCsp_SYM(1:7,1:7);
FISH_BIG_GCsp(7:17,7:17) = FISH_GCsp_SYM(7:17,7:17);
FISH_BIG_GCsp(1:7,7:17) = FISH_GCsp_SYM(1:7,7:17);
FISH_BIG_GCsp(7:17,1:7) = FISH_GCsp_SYM(7:17,1:7);
% Block bias photo and correlations;
FISH_BIG_XC(1:7,1:7) = FISH_XC_SYM(1:7,1:7);
FISH_BIG_XC(21:31,21:31) = FISH_XC_SYM(11:21,11:21);
FISH_BIG_XC(1:7,21:31) = FISH_XC_SYM(1:7,11:21);
FISH_BIG_XC(21:31,1:7) = FISH_XC_SYM(11:21,1:7);
% Block I.A and correlations;
FISH_BIG_XC(18:20,18:20) = FISH_XC_SYM(8:10,8:10);
FISH_BIG_XC(1:7,18:20) = FISH_XC_SYM(1:7,8:10);
FISH_BIG_XC(18:20,1:7) = FISH_XC_SYM(8:10,1:7);

% Final summation
FISH_SYM = FISH_BIG_GCsp + FISH_BIG_XC;

% Add O observable

for i=1:5
  FISH_SYM(7+i,7+i) = FISH_SYM(7+i,7+i) + FISH_O_SYM(i,1,1);
  FISH_SYM(7+i,25+i) = FISH_SYM(7+i,25+i) + FISH_O_SYM(i,1,2);
  FISH_SYM(25+i,7+i) = FISH_SYM(25+i,7+i) + FISH_O_SYM(i,2,1);
  FISH_SYM(25+i,25+i) = FISH_SYM(25+i,25+i) + FISH_O_SYM(i,2,2);
end


%{
FISH_SYM(7+i,7+i) = FISH_SYM(7+i,7+i) + FISH_O_SYM(1,1);
FISH_SYM(7+i,25+i) = FISH_SYM(7+i,25+i) + FISH_O_SYM(1,2);
FISH_SYM(25+i,7+i) = FISH_SYM(25+i,7+i) + FISH_O_SYM(2,1);
FISH_SYM(25+i,25+i) = FISH_SYM(25+i,25+i) + FISH_O_SYM(2,2);
%}

% Force symmetry
FISH_SYM = (tril(FISH_SYM.') + triu(FISH_SYM,1))

% Marginalize FISH_SYM2 in order to get back a 2x2 matrix

% Using Schur complement formula
D = FISH_SYM(13:31,13:31)
inv_D = inv(D)
% Apply formula of Schur complement
COV_ALL = FISH_SYM(1:12,1:12) - FISH_SYM(1:12,13:31)*inv_D*FISH_SYM(13:31,1:12)
FISH_ALL2 = inv(COV_ALL);

%%%%%%%%%%%% BRUTAL METHOD %%%%%%%%%%%
% Invert to marginalyze
%COV_ALL = inv(FISH_SYM);
% Marginalize
%COV_ALL([13:31],:) = [];
%COV_ALL(:,[13:31]) = [];
%FISH_ALL2 = inv(COV_ALL);

%FISH_ALL1_final = FISH_ALL1(3:4,3:4)
%FISH_ALL2_final = FISH_ALL2(3:4,3:4)

FISH_ALL1_final = FISH_ALL1
FISH_ALL2_final = FISH_ALL2

% Matricial equation to solve
eqn = FISH_ALL1_final == FISH_ALL2_final;

% Solving : sigma_o unknown
[solx, parameters, conditions] = solve(eqn, [SIGMA_O(1), SIGMA_O(2), SIGMA_O(3), SIGMA_O(4),
                                             SIGMA_O(5)], 'ReturnConditions', true);

我认为我必须获得“已知方法”来反转矩阵，以避免使用INV Matlab 的函数应用原始反转。

也许这种技术可以被 GPU 所利用。如果不是这样也不错，最重要的是找到一种允许在运行时获得增益的符号矩阵求逆算法。

您对所有可能使用的现有反演算法有什么建议吗？

更新：感谢 Federico，它可能会减少运行时间。

我做了以下修改：

% Using Schur complement formula
D = FISH_SYM(13:31,13:31)
inv_D = inv(D)
% Apply formula of Schur complement
COV_ALL = FISH_SYM(1:12,1:12) - FISH_SYM(1:12,13:31)*inv_D*FISH_SYM(13:31,1:12)
FISH_ALL2 = inv(COV_ALL);

而不是在 31x31 矩阵上进行第一次反转的“大”反转：

COV_ALL = inv(FISH_SYM);
% Marginalize
COV_ALL([13:31],:) = [];
COV_ALL(:,[13:31]) = [];
FISH_ALL2 = inv(COV_ALL);

你觉得我的修改对吗？

1个回答

的 (1,1) 块的逆 是 ( Schur 补码)。如果我从您的解释中正确理解，这就是您要计算的内容（“边缘化”可能是您领域的标准，但它不是标准的线性代数语言）。所以至少你可以减少到反转一个的矩阵，而不是一个 ×31 的矩阵，然后是一个 ×12 的矩阵。

[ABCD]−1
A−BD−1C19×1931×3112×12

其它你可能感兴趣的问题
上一篇用神经网络改进函数逼近 下一篇LP优化和GLPK优化的区别