[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Multiply matrix slices
From: |
Juan Pablo Carbajal |
Subject: |
Re: Multiply matrix slices |
Date: |
Sun, 22 Nov 2015 16:29:54 +0100 |
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
Re: Multiply matrix slices, Parsiad Azimzadeh, 2015/11/22