MATLAB Answers

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

Accepted Answer

Walter Roberson
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 Comments
Walter Roberson
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.

Sign in to comment.

More Answers (2)

petit
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
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

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!