Squeeze : Unable to perform assignment because the size of the left side is 1-by-7-by-7 and the size of the right side is 6-by-6

Hello,

I am trying to use qndiag (from https://github.com/pierreablin/qndiag.git ) from the following function :

function [D, B, infos] = qndiag(C, varargin)
% Joint diagonalization of matrices using the quasi-Newton method
%
% The algorithm is detailed in:
%
%    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
%
% The function takes as input a set of matrices of size `(p, p)`, stored as
% a `(n, p, p)` array, `C`. It outputs a `(p, p)` array, `B`, such that the
% matrices `B * C(i,:,:) * B'` are as diagonal as possible.
%
% There are several optional parameters which can be provided in the
% varargin variable.
%
% Optional parameters:
% --------------------
% 'B0'                        Initial point for the algorithm.
%                             If absent, a whitener is used.
% 'weights'                   Weights for each matrix in the loss:
%                             L = sum(weights * KL(C, C')).
%                             No weighting (weights = 1) by default.
% 'maxiter'                   (int) Maximum number of iterations to perform.
%                             Default : 1000
%
% 'tol'                       (float) A positive scalar giving the tolerance at
%                             which the algorithm is considered to have converged.
%                             The algorithm stops when  |gradient| < tol.
%                             Default : 1e-6
%
% lambda_min                  (float) A positive regularization scalar. Each
%                             eigenvalue of the Hessian approximation below
%                             lambda_min is set to lambda_min.
%
% max_ls_tries                (int), Maximum number of line-search tries to
%                             perform.
%
% return_B_list               (bool) Chooses whether or not to return the list
%                              of iterates.
%
% verbose                     (bool) Prints informations about the state of the
%                             algorithm if True.
%
% Returns
% -------
% D : Set of matrices jointly diagonalized
% B : Estimated joint diagonalizer matrix.
% infos : structure containing monitoring informations, containing the times,
%     gradient norms and objective values.
%
% Example:
% --------
%
%  [D, B] = qndiag(C, 'maxiter', 100, 'tol', 1e-5)
%
% Authors: Pierre Ablin <[e-mail]>
%          Alexandre Gramfort <[e-mail]>
%
% License: MIT

% First tests

if nargin == 0,
    error('No signal provided');
end

if length(size(C)) ~= 3,
    error('Input C should be 3 dimensional');
end

if ~isa (C, 'double'),
  fprintf ('Converting input data to double...');
  X = double(X);
end

% Default parameters

C_mean = squeeze(mean(C, 1));
[p, d] = eigs(C_mean);
p = fliplr(p);
d = flip(diag(d));
B = p' ./ repmat(sqrt(d), 1, size(p, 1));

max_iter = 1000;
tol = 1e-6;
lambda_min = 1e-4;
max_ls_tries = 10;
return_B_list = false;
verbose = false;
weights = [];

% Read varargin

if mod(length(varargin), 2) == 1,
    error('There should be an even number of optional parameters');
end

for i = 1:2:length(varargin)
    param = lower(varargin{i});
    value = varargin{i + 1};
    switch param
        case 'B0'
            B0 = value;
        case 'max_iter'
            max_iter = value;
        case 'tol'
            tol = value;
        case 'weights'
            weights = value / mean(value(:));
        case 'lambda_min'
            lambda_min = value;
        case 'max_ls_tries'
            max_ls_tries = value;
        case 'return_B_list'
            return_B_list = value;
        case 'verbose'
            verbose = value;
        otherwise
            error(['Parameter ''' param ''' unknown'])
    end
end

[n_samples, n_features, ~] = size(C);

D = transform_set(B, C, false);
current_loss = NaN;

% Monitoring
if return_B_list
    B_list = []
end

t_list = [];
gradient_list = [];
loss_list = [];

if verbose
    print('Running quasi-Newton for joint diagonalization');
    print('iter | obj | gradient');
end

for t=1:max_iter
    if return_B_list
        B_list(k) = B;
    end

    diagonals = zeros(n_samples, n_features);
    for k=1:n_samples
        diagonals(k, :) = diag(squeeze(D(k, :, :)));
    end

    % Gradient
    if isempty(weights)
        G = squeeze(mean(bsxfun(@rdivide, D, ...
                                reshape(diagonals, n_samples, n_features, 1)), ...
                         1)) - eye(n_features);
    else
        G = squeeze(mean(...
                bsxfun(@times, ...
                        reshape(weights, n_samples, 1, 1), ...
                        bsxfun(@rdivide, D, ...
                               reshape(diagonals, n_samples, n_features, 1))), ...
                    1)) - eye(n_features);
    end
    g_norm = norm(G);
    if g_norm < tol
        break
    end

    % Hessian coefficients
    if isempty(weights)
        h = mean(bsxfun(@rdivide, ...
                        reshape(diagonals, n_samples, 1, n_features), ...
                        reshape(diagonals, n_samples, n_features, 1)), 1);
    else
        h = mean(bsxfun(@times, ...
                        reshape(weights, n_samples, 1, 1), ...
                        bsxfun(@rdivide, ...
                                reshape(diagonals, n_samples, 1, n_features), ...
                                reshape(diagonals, n_samples, n_features, 1))), ...
                 1);
    end
    h = squeeze(h);

    % Quasi-Newton's direction
    dt = h .* h' - 1.;
    dt(dt < lambda_min) = lambda_min;  % Regularize
    direction = -(G .* h' - G') ./ dt;

    % Line search
    [success, new_D, new_B, new_loss, direction] = ...
        linesearch(D, B, direction, current_loss, max_ls_tries, weights);
    D = new_D;
    B = new_B;
    current_loss = new_loss;

    % Monitoring
    gradient_list(t) = g_norm;
    loss_list(t) = current_loss;
    if verbose
        print(sprintf('%d  - %.2e - %.2e', t, current_loss, g_norm))
    end
end

infos = struct();
infos.t_list = t_list;
infos.gradient_list = gradient_list;
infos.loss_list = loss_list;

if return_B_list
    infos.B_list = B_list
end

end

function [op] = transform_set(M, D, diag_only)
    [n, p, ~] = size(D);
    if ~diag_only
        op = zeros(n, p, p);
        for k=1:length(D)
            op(k, :, :) = M * squeeze(D(k, :, :)) * M';
        end
    else
        op = zeros(n, p);
        for k=1:length(D)
            op(k, :) = sum(M .* (squeeze(D(k, :, :)) * M'), 1);
        end
    end
end

function [v] = slogdet(A)
    v = log(abs(det(A)));
end

function [out] = loss(B, D, is_diag, weights)
    [n, p, ~] = size(D);
    if ~is_diag
        diagonals = zeros(n, p);
        for k=1:n
            diagonals(k, :) = diag(squeeze(D(k, :, :)));
        end
    else
        diagonals = D;
    end
    logdet = -slogdet(B);
    if ~isempty(weights)
        diagonals = bsxfun(@times, diagonals, reshape(weights, n, 1));
    end
    out = logdet + 0.5 * sum(log(diagonals(:))) / n;
end

function [success, new_D, new_B, new_loss, delta] = linesearch(D, B, direction, current_loss, n_ls_tries, weights)
    [n, p, ~] = size(D);
    step = 1.;
    if current_loss == NaN
        current_loss = loss(B, D, false);
    end
    success = false;
    for n=1:n_ls_tries
        M = eye(p) + step * direction;
        new_D = transform_set(M, D, true);
        new_B = M * B;
        new_loss = loss(new_B, new_D, true, weights);
        
        if new_loss < current_loss
            success = true;
            break
        end
        step = step / 2;
    end
    new_D = transform_set(M, D, false);
    delta = step * direction;
end

I use it with the following script :

clc; clear

m=7  % dimension
n=2   % number of matrices

% Load spectro and WL+GCph+XC
FISH_GCsp = load('Fisher_GCsp_flat.txt');
FISH_XC = load('Fisher_XC_GCph_WL_flat.txt');

% Marginalizing over uncommon parameters between the two matrices
COV_GCsp_first = inv(FISH_GCsp);
COV_XC_first = inv(FISH_XC);
COV_GCsp = COV_GCsp_first(1:m,1:m);
COV_XC = COV_XC_first(1:m,1:m);
% Invert to get Fisher matrix
FISH_sp = inv(COV_GCsp);
FISH_xc = inv(COV_XC);
% Drawing a random set of commuting matrices
C=zeros(n,m,m);
B0=zeros(n,m,m);
C(1,:,:) = FISH_sp
C(2,:,:) = FISH_xc

%[D, B] = qndiag(C, 'max_iter', 1e6, 'tol', 1e-6);
[D, B] = qndiag(C);

% Print diagonal matrices
B*C(1,:,:)*B'
B*C(2,:,:)*B'

But unfortunately, I get the following error :

Unable to perform assignment because the size of the left side is 1-by-7-by-7 and the size of the
right side is 6-by-6.

Error in qndiag>transform_set (line 224)
            op(k, :, :) = M * squeeze(D(k, :, :)) * M';

Error in qndiag (line 128)
D = transform_set(B, C, false);

Error in compute_joint_diagonalization (line 32)
[D, B] = qndiag(C);

I don't understand the utility of function squeeze since it seems that it removes a dimension to the op array (6x6 instead of 7x7).

How to fix it ?

Any help is welcome, I put the 2 input files in attachment.

Best regards

Answer 1

[Walter Roberson]

squeeze is not removing the original size.

The code uses eigs(). By default eigs returns information about 6 eigenvalues. your original data is something by 7 x 7 but the helper matrix becomes 6x6 because of eigs and it is the 6 that gets used to initialize the output size but the array being squeezed is 7.

[petit]

@Walter Roberson .ok, thanks for your quick answer. So, what's the suggestion you do to fix this error ?

Regards

[Walter Roberson]

Try changing

[p, d] = eigs(C_mean);

to

[p, d] = eig(C_mean);

unless you will be dealing with large sparse matrices.

If you will be dealing with large sparse matrices then more careful examination would be needed.

Answer 2

[petit]

@Walter Roberson

I tried your solution by using

[p, d] = eig(C_mean) 

or even :

[p, d] = eigs(C_mean,7)

But I still get the following error :

......
......
C(:,:,7) =

   1.0e+07 *

  Columns 1 through 4

  -0.007447517743280   0.005110669785473   0.000281814650722  -0.000229202051452
   1.134581616163971  -0.481981952270819  -0.032067199913283  -0.007354731002997

  Columns 5 through 7

  -0.019695822064900  -0.001780244439019   0.004619441234695
   0.136899307353765   0.238102166567297   0.873581186080297

Index in position 1 exceeds array bounds (must not exceed 2).

Error in qndiag>transform_set (line 224)
            op(k, :, :) = M * squeeze(D(k, :, :)) * M';

Error in qndiag (line 128)
D = transform_set(B, C, false);

Error in compute_joint_diagonalization (line 27)
[D, B] = qndiag(C);

Does it suggest to you at first sight the origin of this error ?

Answer 3

[petit]

@Walter Roberson

Hello,

the previous error into qndiag has been fixed :

Unable to perform assignment because the size of the left side is 1-by-7-by-7 and the size of the
right side is 6-by-6.
Error in qndiag>transform_set (line 224)
            op(k, :, :) = M * squeeze(D(k, :, :)) * M';
Error in qndiag (line 128)
D = transform_set(B, C, false);
Error in compute_joint_diagonalization (line 32)
[D, B] = qndiag(C);

I should replace 2 times :

for k=1:length(D)

by

for k=1:n

into qndiag.m.

By the way, I get constraints (by inversing the Fisher matrix) too bad (I mean, standard deviations are too high for each parameters) :

% 2 Fisher matrices to jointly diagonalize
C(1,:,:) = FISH_sp;
C(2,:,:) = FISH_xc;

% Calling qndiag
[D, B] = qndiag(C, 'max_iter', 1e9, 'tol', 1e-6);

% Check diagonal property for M1
M1 = B*FISH_sp*B'
% Check diagonal property for M2
M2 = B*FISH_xc*B'

% Summing the 2 diagonal matrices
FISH_final = M1 + M2;

% Back to starting space 
FISH_final = B'*FISH_final*B; 

So finally, If I inerse the FISH_final matrix, I get :

wm +/- 357.220530291
wb +/- 771.539306179
w0 +/- 20.0317070287
wa +/- 10.7752879009
h +/- 708.447530242
ns +/- 128.089711461
s8 +/- 169.678103149

FoM = 0.00779787914018

The FoM (Figure of Merit) expected is about ~ 1200 - 1600, so I am far away from these results.

The goal is to combine (do cross-correlations) between the first Fisher matrix (FISH_sp) and the second one (FISH_xc).

So, I wanted to jointly diagonalize the 2 matrices, i.e to find a similar eigenvectors basis, then I sum these 2 diagonal matrices and finally, I get back to parameters starting space by doing :

% Back to starting space 
FISH_final = B'*FISH_final*B; 

If someone could tell see at first sight where the bad results could come from or simply if the method is valid, this would fine to tell it.

Best regards