4

I would like to compute the product of the next n adjacent elements of a matrix. The number n of elements to be multiplied should be given in function's input. For example for this input I should compute the product of every 3 consecutive elements, starting from the first.

[p, ind] = max_product([1 2 2 1 3 1],3);

This gives [1*2*2, 2*2*1, 2*1*3, 1*3*1] = [4,4,6,3].

Is there any practical way to do it? Now I do this using:

for ii = 1:(length(v)-2)
    p = prod(v(ii:ii+n-1));
end

where v is the input vector and n is the number of elements to be multiplied.

in this example n=3 but can take any positive integer value.

Depending whether n is odd or even or length(v) is odd or even, I get sometimes right answers but sometimes an error.
For example for arguments:

v = [1.35912281237829 -0.958120385352704 -0.553335935098461 1.44601450110386 1.43760259196739 0.0266423803393867 0.417039432979809 1.14033971399183 -0.418125096873537 -1.99362640306847 -0.589833539347417 -0.218969651537063 1.49863539349242 0.338844452879616 1.34169199365703 0.181185490389383 0.102817336496793 0.104835620599133 -2.70026800170358 1.46129128974515 0.64413523430416 0.921962619821458 0.568712984110933] 
n = 7

I get the error:

Index exceeds matrix dimensions.
Error in max_product (line 6)  
p = prod(v(ii:ii+n-1));

Is there any correct general way to do it?

Dev-iL
  • 23,742
  • 7
  • 57
  • 99
esem
  • 153
  • 1
  • 12
  • I suspect you don't get a wrong answer, you get an index out of bounds error, correct? – beaker Sep 04 '16 at 14:52
  • Right I get Index exceeds matrix dimensions. Error in max_product (line 6) p=prod(v(ii:ii+n-1)); – esem Sep 04 '16 at 14:53
  • That's because the limits on the range for `ii` are incorrect. Hint: When `ii` is at its maximum value, `ii+n-1` should give you the last element of `v`. – beaker Sep 04 '16 at 14:55
  • @beaker that does not really help. I get the error long before the ii reaches its maximum value. If ii=1:15, I get the error already from the 10th iteration. thats right my question is how to organize the range of ii. making it dependent on length(v) is wrong, trying to find the cases when n is odd or even is also wrong. – esem Sep 04 '16 at 15:36
  • What is `n` in this case? Why is this not in the original question? Specifically, what is `10+n-1`? Is it greater than the length of `v`? – beaker Sep 04 '16 at 15:48
  • Incidentally, I suspect you also want something like `p(ii) = prod...` rather than overwriting the value each iteration. – beaker Sep 04 '16 at 15:56
  • n= is the number of elements to be multiplied. (in my example 3). it means how many consecutive elements should be taken for multiplication. It can't be longer than v. 10+n-1 is my way of looping over elements and multiplying. But for some arguments I get error – esem Sep 04 '16 at 16:00
  • 2
    There is only one mistake in your code. You should change the beginning of your for loop to `for ii = 1:(length(v)-n+1)` and then it should work. – Erfan Sep 04 '16 at 16:53
  • @erfan that is wonderful. Thank you. yes that was the mistake. please post it as an answer I will accept it. – esem Sep 04 '16 at 17:04

6 Answers6

7

Based on the solution in Fast numpy rolling_product, I'd like to suggest a MATLAB version of it, which leverages the movsum function introduced in R2016a.

The mathematical reasoning is that a product of numbers is equal to the exponent of the sum of their logarithms:

enter image description here

A possible MATLAB implementation of the above may look like this:

function P = movprod(vec,window_sz)
  P = exp(movsum(log(vec),[0 window_sz-1],'Endpoints','discard'));
  if isreal(vec)   % Ensures correct outputs when the input contains negative and/or
    P = real(P);   %   complex entries.
  end
end

Several notes:

  1. I haven't benchmarked this solution, and do not know how it compares in terms of performance to the other suggestions.
  2. It should work correctly with vectors containing zero and/or negative and/or complex elements.
  3. It can be easily expanded to accept a dimension to operate along (for array inputs), and any other customization afforded by movsum.
  4. The 1st input is assumed to be either a double or a complex double row vector.
  5. Outputs may require rounding.
Community
  • 1
  • 1
Dev-iL
  • 23,742
  • 7
  • 57
  • 99
5

Update

Inspired by the nicely thought answer of Dev-iL comes this handy solution, which does not require Matlab R2016a or above:

out = real( exp(conv(log(a),ones(1,n),'valid')) )

The basic idea is to transform the multiplication to a sum and a moving average can be used, which in turn can be realised by convolution.

Old answers

This is one way using gallery to get a circulant matrix and indexing the relevant part of the resulting matrix before multiplying the elements:

a = [1 2 2 1 3 1]
n = 3

%// circulant matrix
tmp = gallery('circul', a(:))
%// product of relevant parts of matrix
out = prod(tmp(end-n+1:-1:1, end-n+1:end), 2)

out =

     4
     4
     6
     3

More memory efficient alternative in case there are no zeros in the input:

a = [10 9 8 7 6 5 4 3 2 1]
n = 2

%// cumulative product
x = [1 cumprod(a)] 
%// shifted by n and divided by itself
y = circshift( x,[0 -n] )./x 
%// remove last elements 
out = y(1:end-n) 

out =

    90    72    56    42    30    20    12     6     2
Community
  • 1
  • 1
Robert Seifert
  • 25,078
  • 11
  • 68
  • 113
  • the vector a and n can be changed. they are input arguments. If i give a=[10 9 8 7 6 5 4 3 2 1] and n=2 i should get j=[90 72 56 42 30 20 12 6 2]. In your function I get something very different – esem Sep 04 '16 at 15:18
  • I have a look at it later – Robert Seifert Sep 04 '16 at 16:29
  • @AstrikSahakyan I made a correction. Now it should work. – Erfan Sep 04 '16 at 16:49
  • @erfan thanks for the correction, I made another change to avoid the flipud – Robert Seifert Sep 04 '16 at 20:28
  • @AstrikSahakyan I added a more memory efficient solution you should consider. – Robert Seifert Sep 04 '16 at 20:47
  • @thewaywewalk Your 2nd solution only works if the input contains no zeros (i.e. replace `10` with `0` and see what happens). – Dev-iL Sep 05 '16 at 12:46
  • @Dev-iL thanks for the comment, that lies in the nature of `cumprod`. I will add that. – Robert Seifert Sep 05 '16 at 12:48
  • @thewaywewalk nice solution it helped me reduce my code by 70%. The only problem I encounter is for arguments `a=[-0.289733526202033 -0.424117761081419 -0.109225062928504 0.437558086660525 -1.07901784704407 -1.40123341459366 0.99784653114119 0.0410647599704143 -0.136283604829852 -1.07514168008808 2.10545308817616 1.4640416843157 -0.223038336113035 0.580459996777139 2.45201246105665 -0.0361913086341127 1.41307815119679] and n=3` I receive a vector of complex numbers. Any idea why? – esem Sep 06 '16 at 21:04
  • @AstrikSahakyan Why? Because the logarithm of a negative number is a complex. But as `exp` inverts the calculation, the imaginary part is actually zero. If it bothers you, you can just cast the result to be `real` -> have a look at my edit. – Robert Seifert Sep 06 '16 at 21:47
2

Your approach is correct. You should just change the for loop to for ii = 1:(length(v)-n+1) and then it will work fine.

If you are not going to deal with large inputs, another approach is using gallery as explained in @thewaywewalk's answer.

Community
  • 1
  • 1
Erfan
  • 1,897
  • 13
  • 23
2

I think the problem may be based on your indexing. The line that states for ii = 1:(length(v)-2) does not provide the correct range of ii.

Try this:

function out = max_product(in,size)
    size = size-1;                 % this is because we add size to i later
    out = zeros(length(in),1)      % assuming that this is a column vector
    for i = 1:length(in)-size
        out(i) = prod(in(i:i+size));
    end

Your code works when restated like so:

for ii = 1:(length(v)-(n-1))
    p = prod(v(ii:ii+(n-1)));
end

That should take care of the indexing problem.

Adriaan
  • 17,741
  • 7
  • 42
  • 75
Alea Kootz
  • 913
  • 4
  • 11
  • Please do not beg for acceptance of your answer. The OP will choose whichever answer suits him/her best. – Adriaan Sep 15 '16 at 12:29
1

using bsxfun you create a matrix each row of it contains consecutive 3 elements then take prod of 2nd dimension of the matrix. I think this is most efficient way:

max_product = @(v, n) prod(v(bsxfun(@plus, (1 : n), (0 : numel(v)-n)')), 2);
p = max_product([1 2 2 1 3 1],3)

Update: some other solutions updated, and some such as @Dev-iL 's answer outperform others, I can suggest fftconv that in Octave outperforms conv

rahnema1
  • 15,264
  • 3
  • 15
  • 27
1

If you can upgrade to R2017a, you can use the new movprod function to compute a windowed product.

CKT
  • 781
  • 5
  • 6