4 views (last 30 days)

Show older comments

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

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.

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.

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

Start Hunting!