octave-maintainers
[Top][All Lists]
Advanced

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

Re: Multiply matrix slices


From: Parsiad Azimzadeh
Subject: Re: Multiply matrix slices
Date: Sun, 22 Nov 2015 11:49:56 -0500

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.



reply via email to

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