speeding up function

[Carlo de Falco]

Hi all,

I often store sets of matrices and vectors as part of NDArrays, e.g.

for ii=1:Ni
  for jj=1:Nj
    v(ii, :, jj) = randn (2, 1);
    m(ii, :, :, jj) = randn (2, 2);
  endfor
endfor

so, to be able to perform matrix-vector operations like

for ii=1:Ni
  for jj=1:Nj
    detm(ii, jj) = det (m(ii,:,:,jj));
  endfor
endfor

or

for ii=1:Ni
  for jj=1:Nj
    mv2(ii, :, jj) = squeeze (m(ii,:,:,jj)) * (v(ii,:,jj).^2)(:);
  endfor
endfor

without for cycles, I produced the funcion below.
Although it works it seems to be a bit slow.
Does anyone on the list see any obvious bottleneck I coud remove to  
speed it up?

Thanks in advance,
c.

-----------------------------------------------------------------------------------------------


% MATARRAYFUN: Apply a given functions to a set of matrices
%              concatenated into an array.
%
%      res = matarrayfun (fun, mat_array, dims, [mat_arrya2,  
dims2, ...])
%
% INPUTS:
%
%    fun:               function handle to apply to each matrix
%    mat_array:         n-dimensional array
%    dims:              dimensions into mat_arry that represent the
%                       rows and columns of the matrices. mat_array
%                       will be permuted so that these dimensions come  
first
%                       in the result
%    mat_array2, dims2: if fun takes more inputs, they can be passed  
as extra
%                       arguments to matarrayfun
%
% OUTPUTS:
%
%    res:         result of the application of fun to mat_array
%                 if the option 'scalar' is set to false res has
%                 the same size as mat_array,
%                 otherwise ndims (res) =  ndims (mat_array) - 2
%                 and size (res, i) = size (mat_array, i+2) for any i
%
% Copyright (C) 2010 Carlo de Falco
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.

% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.
% Author: Carlo de Falco <cdf AT users.sourceforge.net>


function  res = matarrayfun (fun, ma, dims, varargin)

dim_rep = setdiff (1:ndims(ma), dims);
s       = size (ma);
num_rep = s(dim_rep);
mat_size= s(dims);
if (numel(num_rep) == 1)
  num_rep = [1 num_rep];
end
in = {};
if ~isempty (varargin)
  in = cellfun (@(x, y) permute(x, [y, setdiff(1:ndims(x), y)]),  
varargin(1:end-1), varargin(2:end), 'uniformoutput', false);
  in = cellfun (@(x, y) squeeze (num2cell (x, 1:numel(y))), in,  
varargin(2:end), 'uniformoutput', false);
end
res = cellfun (fun, reshape (num2cell (permute (ma, [dims, dim_rep]),  
1:numel(dims)), num_rep), in{:}, 'uniformoutput', false);
res = squeeze (cell2mat (reshape (res, [ones(1, ndims (res{1})),  
num_rep])));
end

%!shared a, b
%!test
%! a = ones (7, 5, 4, 4, 3, 2);
%! res = matarrayfun (@det, a, [3, 4]);
%! assert (size (res), [7 5 3 2]);
%! assert (any (res(:)), false);
%!test
%! res = matarrayfun (@transpose, a, [3, 4]);
%! assert (res, ones (4, 4, 7, 5, 3, 2));
%!test
%! b = ones (7, 5, 3, 2, 4);
%! res = matarrayfun (@mtimes, a, [3, 4], b, 5);
%! assert (res, 4*ones (4, 7, 5, 3, 2));
%!test
%! aa = randn (4, 4); bb =  randn (4, 1);
%! a(1, 1, :, :, 1, 1) = aa;
%! b(1, 1, 1, 1, :) = bb;
%! res = matarrayfun (@mtimes, a, [3, 4], b, 5);
%! assert (res (:, 1, 1, 1, 1), (aa*bb)(:))
%!test
%! aa = randn (4, 4, 100); bb =  randn (4, 100);
%! res = matarrayfun (@mtimes, aa, [1, 2], bb, 1);
%! assert (size (res), [4 100])
%!test
%! aa = randn (4, 4, 100);
%! res = matarrayfun (@(x) inv(x'), aa, [1, 2]);
%! for ii=1:100
%!   assert (squeeze(res (:,:,ii)), inv(squeeze(aa(:,:,ii))'))
%! endfor
%!test
%! aa = randn (4, 4, 3, 2, 3);
%! res = matarrayfun (@(x) inv(x)', aa, [1, 2]);
%! for ii=1:3, for jj=1:2, for kk=1:3
%!   assert ((res (:,:,ii, jj, kk)), inv(squeeze(aa(:,:,ii, jj, kk)))')
%! endfor, endfor, endfor