# 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

4 views (last 30 days)
petit on 26 Jan 2021
Edited: petit on 28 Jan 2021
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 <pierre.ablin@inria.fr>
% Alexandre Gramfort <alexandre.gramfort@inria.fr>
%
% 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 = [];
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 = [];
loss_list = [];
if verbose
print('Running quasi-Newton for joint diagonalization');
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
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
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.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
% 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

Walter Roberson on 26 Jan 2021
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.
##### 2 CommentsShowHide 1 older comment
Walter Roberson on 26 Jan 2021
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.

petit on 26 Jan 2021
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 ?

petit on 28 Jan 2021
Edited: petit on 28 Jan 2021
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