title: 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
body: <p>I am looking for a way to find same eigenvectors for 2 given matrices, this way I would make a joint diagonalisation. For this, I found out and tried to use qndiag (from <a href="https://github.com/pierreablin/qndiag.git" rel="nofollow noreferrer">https://github.com/pierreablin/qndiag.git</a> ) from the following function :</p>
  <pre><code>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| &lt; 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 &lt;pierre.ablin@inria.fr&gt;
  %          Alexandre Gramfort &lt;alexandre.gramfort@inria.fr&gt;
  %
  % 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 &lt; 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 &lt; 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 &lt; current_loss
              success = true;
              break
          end
          step = step / 2;
      end
      new_D = transform_set(M, D, false);
      delta = step * direction;
  end
  </code></pre>
  <p>I use it with the following script that you can test with downloading the 2 input matrices at the bottom of this post :</p>
  <pre><code>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'
  </code></pre>
  <p>But unfortunately, I get the following error :</p>
  <pre><code>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&gt;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);
  </code></pre>
  <p>I don't understand the utility of function <code>squeeze</code> the most important : why the function <code>eigs</code> returns <strong>only 6 values and not 7</strong> like in my data (the <strong>input matrices has 7x7 size</strong>).</p>
  <p>What might be wrong with this issue of array dimension and how can I fix it ?</p>
  <p>I put the 2 input files available here :</p>
  <p><a href="http://31.207.36.11/Fisher_GCsp_flat.txt" rel="nofollow noreferrer">Matrix Fisher_GCsp_flat.txt</a></p>
  <p><a href="http://31.207.36.11/Fisher_XC_GCph_WL_flat.txt" rel="nofollow noreferrer">Matrix Fisher_XC_GCph_WL_flat.txt</a></p>
  <p>You can test the above code that calls <code>qndiag</code> for these 2 matrices.</p>
  <h2>Update 1</h2>
  <p>To allow people interested to test quickly the code, I put a link of the archive:</p>
  <p><a href="http://31.207.36.11/Archive_Matlab_StackOver.tar" rel="nofollow noreferrer">Archive_Matlab_StackOver.tar</a></p>
  <p>You just have to untar and execute under Matlab the script <code>compute_joint_diagonalization.m</code> and you will see normally the above error (regarding the use of <code>eigs</code> and <code>squeeze</code> functions).</p>
  <p>It should help you understand the origin of this issue.</p>
  <h2>Update 2</h2>
  <p>If I replace <code>[p, d] = eigs(C_mean)</code> by  <code>[p, d] = eigs(C_mean,7)</code> , I get another error :</p>
  <pre><code>Index in position 1 exceeds array bounds (must not exceed 2).
  
  Error in qndiag&gt;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);
  </code></pre>
  <p>However, the dimensions of the 2 matrices used are 7x7 and should be correctly processed with <code>eigs(C_mean,7)</code>.</p>
  <h2>Update 3</h2>
  <p>The size of <code>op, D, M</code> and <code>k</code> are equal to (including after the error message) :</p>
  <pre><code>size(D) =
       2     7     7
  
  length(D) =
       7
     
  size(M) =
       7     7
  
  size(op) =
       2     7     7
  
  Index in position 1 exceeds array bounds (must not exceed 2).
  
  Error in qndiag&gt;transform_set (line 231)
              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);
  </code></pre>
  <p>Notice  that <code>k</code> varies from 1 to <code>length(D)=7</code>.</p>
  <p>Is there an issue which could appear with these dimensions ?</p>