octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Multiply matrix slices


From: Nicholas Jankowski
Subject: Re: Multiply matrix slices
Date: Sun, 22 Nov 2015 07:24:30 -0500



On Sun, Nov 22, 2015 at 7:21 AM, Nicholas Jankowski <address@hidden> wrote:
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

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



oh, I forgot it also works if one matrix is 2x2x1 and the other is 2x2xn, but only in Octave since that would use implicit broadcasting.  would have to rewrite with bxsfun to be matlab compatible.

reply via email to

[Prev in Thread] Current Thread [Next in Thread]