On Sun, Nov 22, 2015 at 10:29 AM, Juan Pablo Carbajal <address@hidden> wrote: > On Sun, Nov 22, 2015 at 1:24 PM, Nicholas Jankowski <address@hidden> wrote: >> >> >> 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. > > check ndmultiply. does it do what you want? > http://sourceforge.net/p/octave/linear-algebra/ci/default/tree/inst/ndmult.m
It seems like ndmultiply might be robust enough to accomplish this. However, I think I might be able to rework the function using only multi-dimensional matrix-vector multiplies. In case this is useful to anyone:
# Multiply m-by-n matrix with nx1 vector over p slices A = rand(m, n, p); x = rand(n, p); b = sum(bsxfun(@times, A, permute(x, [3 1 2])),2);
This is presumably faster than the for-loop approach.