1
$\begingroup$

I have to solve the equality between 2 matrixes 12x12 containing a lot of symbolic variables and with which I perform inversion of matrix. There is only one unknown called "SIGMA_O".

My system is solved fastly when I take for example 2 matrices 2x2, the inversion is pretty direct.

Now, with the case of 2 matrixes 12x12, I need before actually to inverse a 31x31 matrix of symbolic variables (I marginalize after), since inversion takes a lot of time.

I would like to benefit from my GPU NVIDIA card to achieve this inversion faster but the GPU optimization is not supported currently for Symbolic arrays.

Below the script where you will find the line of inversion :

COV_ALL = inv(FISH_SYM)

and the entire code :

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);

I think that I have to get "known methods" to inverse a matrix, in order to avoid to apply a raw inversion with INV Matlab's function.

Maybe this technique could be exploited by GPU. If not this is no bad, the most important is to find an algorithm of inversion of symbolic matrix that allows to gain in runtime.

Have you got any suggestions about all the existing algorithms of inversion to use potentially ?

UPDATE: thanks to Federico, it might possible to reduce the runtime.

I did the follwing modifications :

% 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);

instead of "big" inversion with first inversion on 31x31 matrix:

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

Do you think my modifications are right ?

Improve this question
$\endgroup$
4
  • $\begingroup$ Before describing how you want to do things, can you explain what it is you want to do in mathematical terms? $\endgroup$
    – Wolfgang Bangerth
    Apr 23, 2021 at 16:29
  • $\begingroup$ @WolfgangBangerth. It is a little bit complicated to explain. First, the system of equations is just an equality with 2 matrixes 12x12. The first one is obtained quickly by removing columns/rows to get a first 12x12 matrix. The second matrix is more difficult to get : I have to inverse a 31x31 matrix, then on the inverse matrix, I marginalize by removing all nuisance terms, that is to say, by removing colums/rows to get a 12x12 matrix and I reinverse this latter to finaly have the second matrix equal to the first one described above. Hoping this is a little more clear for you. $\endgroup$
    – user14831
    Apr 23, 2021 at 18:09
  • 1
    $\begingroup$ Even if you can run it, chances are that the answer will be a long, unreadable mess. What do you really want to know about this inverse? Maybe there is a faster way to know it? $\endgroup$
    – Federico Poloni
    Apr 23, 2021 at 18:50
  • $\begingroup$ @FedericoPoloni the inversion is necessary since I am working with Fisher matrixes, and to know the associated covariance matrix, I need to invert it. What I am looking for is alternatives methods to Matlab INV functions like for example Gauss method or other methods which would be relatively easy to implement. My only unknown is SIGMA_0, Sorry, I have no choice, I am obliged to inverse. $\endgroup$
    – user14831
    Apr 23, 2021 at 21:18

1 Answer 1

Reset to default
1
$\begingroup$

The inverse of the (1,1) block of $$ \begin{bmatrix} A & B\\ C & D \end{bmatrix}^{-1} $$ is $A-BD^{-1}C$ (Schur complement). This is what you are trying to compute, if I understand correctly from your explanation ("marginalize" may be standard in your domain, but it is not standard linear algebra language). So at least you can reduce to inverting a $19\times 19$ matrix rather than a $31\times 31$ one followed by a $12\times 12$ one.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks and sorry for my specific language of "marginalize" : this simply means that I remove the useless rows/columns which doesn't modify the constraints given the fact that I am working with Fisher matrix. $\endgroup$
    – user14831
    Apr 23, 2021 at 22:13
  • $\begingroup$ And what do you think about if block $A$, $B$ are symmetric and $B$ is the symmetric bloc of $C$ ? I mean that my entire starting matrix is symmetric : is threre a specific method for this kinf of matrix to inverse ? And in your solution, you are talking about a 19x19 matrix : I don't have such a block of this size : how to deal with it ? $\endgroup$
    – user14831
    Apr 23, 2021 at 22:18
  • $\begingroup$ Could you write please a code snippet in order to better understand your method from my little script in the original post ? $\endgroup$
    – user14831
    Apr 23, 2021 at 22:26
  • $\begingroup$ I don't want to be boring but could you give me please help about your suggestion, I didn't grasp all the subtilities of the 2 inversions that you talk about and mostly how to adapt it in my code. $\endgroup$
    – user14831
    Apr 24, 2021 at 7:02
  • $\begingroup$ @youpilat13 your 13:31 block has size 19. On a symmetric matrix, the Schur complement is symmetric. No, sorry, I am not going to do the coding for you. $\endgroup$
    – Federico Poloni
    Apr 24, 2021 at 7:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.