[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: |
Mon, 23 Nov 2015 12:23:30 -0500 |
Thanks everyone for your input.
> Can you use `reshape` to transform the matrices from m x n x k to m x (n*k),
> multiply them, then reshape again?
I could not think of a way to implement your suggestion. I ended up
using the bsxfun approach detailed in a previous message.
Timing information below for anyone who's interested:
MATLAB financial toolbox results
================================
N Elapsed time (secs.)
----------------------------
16 0.543231
32 1.053423
64 2.167072
128 4.191894
256 8.361655
512 16.609718
1024 32.839757
octave-financial results
========================
N Elapsed time (secs.)
----------------------------
16 0.048691
32 0.064110
64 0.097092
128 0.162552
256 0.294098
512 0.568558
1024 1.136864
Code used for timing
====================
Asset1Price = 100.; Asset2Price = 90. ;
Volatility1 = 0.2 ; Volatility2 = 0.3 ;
Dividends1 = 0. ; Dividends2 = 0.005;
RiskFreeRate = 0.04;
Correlation = 0.5;
ExpiryTime = 1.;
Drift = drift ([0;0], [RiskFreeRate-Dividends1 0;0
RiskFreeRate-Dividends2]);
Diffusion = diffusion ([1;1], [Volatility1 0;0 Volatility2]);
SDE = sde (Drift, Diffusion, 'StartState', [Asset1Price;Asset2Price], ...
'Correlation', [1 Correlation;Correlation 1]);
M = 1000;
Ns = 2.^(4:10);
for N = Ns
t0 = clock ();
[Paths, ~, ~] = simulate (SDE, 1, 'NTRIALS', M, 'NSTEPS', N);
elapsed_time = etime (clock (), t0);
fprintf('%d\t%f\n', N, elapsed_time);
end
On Sun, Nov 22, 2015 at 3:17 PM Nicholas Jankowski <address@hidden> wrote:
>
> On Sun, Nov 22, 2015 at 12:32 PM, Alex Krolick <address@hidden> wrote:
>>
>> Can you use `reshape` to transform the matrices from m x n x k to m x (n*k),
>> multiply them, then reshape again?
>> ________________________________
>> From: Parsiad Azimzadeh
>> Sent: 11/22/2015 8:50 AM
>> To: Juan Pablo Carbajal
>> Cc: octave-maintainers
>> Subject: Re: Multiply matrix slices
>>
>> 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.
>>
>
>
> Just FYI here's the Stackoverflow thread where they discussed different
> approaches and timings. if you know the dimensions, unrolling the matrix
> multiply like I did for the two-by-two was an order of magnitude faster. It
> may be worth having explicit code paths for some lower order matrix
> dimensions, and in scales well to n-dimensions with the : indexing. They do
> discuss other approaches as well:
>
> http://stackoverflow.com/questions/6580656/matlab-how-to-vector-multiply-two-arrays-of-matrices
Re: Multiply matrix slices, Parsiad Azimzadeh, 2015/11/22