Re: Multiply matrix slices

[Nicholas Jankowski]

On Sun, Nov 22, 2015 at 2:05 AM, Parsiad Azimzadeh <address@hidden> wrote:

A triply-nested loop in an octave-financial routine I wrote renders it fairly slow: http://sourceforge.net/p/octave/financial/ci/default/tree/inst/@sde/simByEuler.m.

This could be fixed through a _native_ routine to multiply slices. That is, if A and B are 3-dimensional real arrays, C = multiply_slices(A,B) would give:

C(:, :, 1) == A(:, :, 1)*B(:, :, 1),

...,

C(:, :, k) == A(:, :, k)*B(:, :, k).

I am guessing no such routine exists in the core. I may also be overlooking other possible solutions. Has anyone run into a similar problem/have advice?

Thanks,

Parsiad

Are your matrices always a predictable size?  I was doing an eigenvector based program over the summer that worked with n-D arrays of 2x2 matrices. Had to do the 'multiply by slice' problem you mentioned quite a bit. The fastest way to do it, since my matrix size was fixed, was to 'unroll the multiply and do vectorized multiplication of each element. I don't think there is a good vectorized way to do it without knowing the size apriori, though. Found a speed comparison somewhere on stackexchange I think.

Code below, and m-fie attached. It works for any set of 2 x 2 x n x m x p x  ...  as long as they are the same size.  It's far from 'Octave core' worthy, though.


------------------------------
function C = twobytwoarraymult(A,B)
   %either A and B should be 2x2 matrices. either or both can be a 2x2xn matrix array,
   %but if both are arrays, they must have the same n.  Output
   %will be matrix multiplication of 2x2 matrices along 3rd dimentsion in A and B.
   %either A or B will be broadcast in 3rd dimension if it is only a single matrix,
   %otherwise it needs to be the same size as A.
   %will work for higher dimensions than three. the : will recast them as 3D arrays for
   %the multiply, and C will be returned with the intended multidimensional size.
   
  
   if (ndims(A)>=ndims(B))
    C = zeros(size(A));
  else
    C = zeros(size(B));
  end
 
  C(1,1,:) = A(1,1,:).*B(1,1,:) + A(1,2,:).*B(2,1,:);
  C(1,2,:) = A(1,1,:).*B(1,2,:) + A(1,2,:).*B(2,2,:);
  C(2,1,:) = A(2,1,:).*B(1,1,:) + A(2,2,:).*B(2,1,:);
  C(2,2,:) = A(2,1,:).*B(1,2,:) + A(2,2,:).*B(2,2,:);
end

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

twobytwoarraymult.m
Description: Text Data