2

In MATLAB I have a vector x of length n, where n is usually O(10), and I would like to build a tall matrix A of size [n^m,m], where m is again 0(10). The matrix has a peculiar form: if n=4 and m=6, let

x=[x1; x2; x3; x4]

then A is

   x1 x1 x1 x1 x1 x1
   x2 x1 x1 x1 x1 x1
   x3 x1 x1 x1 x1 x1
   x4 x1 x1 x1 x1 x1
   x1 x2 x1 x1 x1 x1
   x2 x2 x1 x1 x1 x1
   x3 x2 x1 x1 x1 x1
   x4 x2 x1 x1 x1 x1
   x1 x3 x2 x1 x1 x1
   .               .         
   .               . 
   .               .
   x4 x4 x4 x4 x4 x4

In practice, each column is obtained by repeating the elements of x, with an increasing stride for each column. How can I do that? If possible, I'd prefer an efficient (vectorized) solution, because, as you can see, the number of rows of A increases exponentially with m.

EDIT: whoops, sorry! I forgot I also need to build a vector V of size [n^m,1], based on vector w having the same length of x

w=[w1; w2; w3; w4]

V is

   w1^6
   w2*w1^5
   w3*w1^5
     .
     .
     .
   w4^6

Hope the crappy graphics is clear enough. Anyway, V is a column vector of lenght n^m. Guess I could create a matrix B from w, in the same way one creates a matrix A from x, and then use prod(B,2)?

Nimantha
  • 6,405
  • 6
  • 28
  • 69
DeltaIV
  • 4,773
  • 12
  • 39
  • 86

5 Answers5

5

Use allcomb tool from MATLAB file-exchange to generate the possible combinations of the indices [1 2 3 4] and then use them to index into x -

v = repmat({1:numel(x)},1,m);
A = x(fliplr(allcomb(v{:})));

Also, it seems instead of using fliplr, you can use - allcomb(v{:},'matlab') instead.

For the edited part of the question, you can use a modified version of it -

V = prod(x(allcomb(v{:})),2)

Benchmarking

Please note that these are for runnable solutions posted here.

Benchmarking Code

%// Parameters and input x
n = 10; m = 6;num_runs = 20; x =  randi(9,n,1);

disp('-------- With allcomb')
tic
for runs = 1:num_runs
    v = repmat({1:numel(x)},1,m);
    A = x(fliplr(allcomb(v{:})));
end
toc,    clear v A

disp('-------- With bsxfun')
tic
for runs = 1:num_runs
    A = x(floor(mod(bsxfun(@rdivide, (0:n^m-1).', n.^[0:m-1] ), n)+1)); %//'
end
toc,    clear A

disp('-------- With ttable')
tic
for runs = 1:num_runs
    I = ttable(n*ones(1,m));
    A = x(I);
end
toc,    clear I A

disp('-------- With arrayfun')
tic
for runs = 1:num_runs
    A = cell2mat(arrayfun(@(i)...
        (repmat(reshape(repmat(x',n^(i-1),1),[],1),n^(m-i),1)),1:m,'uni',0));
end
toc

Results

-------- With allcomb
Elapsed time is 6.544981 seconds.
-------- With bsxfun
Elapsed time is 11.547062 seconds.
-------- With ttable
Elapsed time is 15.729932 seconds.
-------- With arrayfun
Elapsed time is 4.319048 seconds.
Divakar
  • 218,885
  • 19
  • 262
  • 358
  • @Divakar, allcomb seems fast! Also, it builds an index matrix, which may be useful for vector V (see my edit, and apologies if I forgot to write it before!). I will try it out. Thanks a lot!!! – DeltaIV Aug 11 '14 at 15:12
  • 1
    @DeltaIV Exactly! Just use that index array and then use `prod` for the edit part - `prod(fliplr(allcomb(v,v,v,v,v,v)),2)` I think (Not tested). In fact you don't need to use `fliplr` then! Good luck! – Divakar Aug 11 '14 at 15:13
  • @Divakar, Ouch! I have an issue, though...m (the number of dimensions) is a parameter in my code. So I cannot hardcode `allcomb(x,x,x,x,x,x)` in my code, because I have no idea how many copies of x's I will need to multiply, until runtime! Is there a solution? – DeltaIV Aug 11 '14 at 15:25
  • @DeltaIV Check out the edited code! Good call on that, now it's generic! Also incorporated the solution to the EDIT into the solution. – Divakar Aug 11 '14 at 15:40
  • @Divakar Great! BTW, `allcomb` has an argument `matlab` which makes the first column change fastest,so using `allcomb(v{:},'matlab')` I got rid of the call to `fliplr`. – DeltaIV Aug 11 '14 at 15:44
  • @DeltaIV WOW!! Awesomeness in that tool! So you are all set it seems! :) – Divakar Aug 11 '14 at 15:45
4

One-liner based only on built-in functions (namely, mod and the very powerful bsxfun):

result = x(floor(mod(bsxfun(@rdivide, (0:n^m-1).', n.^[0:m-1] ), n)+1));
Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
  • This should be the fastest one. (I have been looking for your answer and here it is!) – Yvon Aug 11 '14 at 14:35
  • @Yvon :-) Yes, it's vectorized and looks efficient, so it should be very fast. In fact, I was thinking about doing some benchmarking... but most other answers involve downloading functions and I can't be bothered – Luis Mendo Aug 11 '14 at 14:38
  • 2
    @LuisMendo Posted some benchmarking results in my solution. – Divakar Aug 11 '14 at 15:00
  • 1
    Don't understand what it does, but respect is due for one-liners using bsxfun! +1 – DeltaIV Aug 11 '14 at 15:16
  • 1
    @LuisMendo Something with `bsxfun and @rdivide..`. You have my +1 too! – Divakar Aug 11 '14 at 15:21
  • @DeltaIV `bsxfun` applies an operation to different combinations of values spread along different dimensions. Try `bsxfun(@plus, [1 2 3], [10 20].')` to grasp what it does (note the transpose in the second vector) – Luis Mendo Aug 11 '14 at 18:53
  • @LuisMendo, thanks a lot for the answer! I know that `bsxfun` is used to apply an operation between two different-sized arrays, expanding them if needed. I still don't get the meaning of your one-liner. For starters, why the dot in `(0:n^m-1).'`? `(0:n^m-1)'` seems to generate the same vector. Why the `mod` function? If you used `ceil` instead of `floor`, would that get rid of the two -1? The main problem for me is that I don't understand what the code is doing. – DeltaIV Aug 12 '14 at 10:10
  • @LuisMendo `.'` is transpose, whereas `'` is conjugate transpose. See a discussion [here](http://stackoverflow.com/q/25150027/2586922). `mod` is applied to get only numbers from `1` to `n`. As for `ceil`, I think it wouldn't work at the "edge" cases. It's true that the code in one line is cryptic. I suggest you see intermediate results one by one (such as `mod(bsxfun(@rdivide, (0:n^m-1).', n.^[0:m-1] ), n)` etc) to understand what the code does – Luis Mendo Aug 12 '14 at 11:06
2

Try this: (all are bulit-in functions)

A = cell2mat(arrayfun(@(i)(repmat(reshape(repmat(x',n^(i-1),1),[],1),n^(m-i),1)),1:m,'UniformOutput',0))

Explanation:

n = 2;
m = 4;
x = (1:n)';
A = [];
for i = 1:m
%// temp1 is (n^(i-1)) x n matrix with each row equal to x' 
    temp1 = repmat(x',n^(i-1),1); 
%// temp2 is (n^(i-1))*n x 1 column vector with corresponding elements of temp1
    temp2 = reshape(temp1,[],1);
%// temp3 is a (n^(m-i))*(n^(i-1))*n x 1, i.e n^m x 1 column vector with elements of temp2 repeated n^(m-i) times
    temp3 = repmat(temp2,n^(m-i),1);
%// A is appending temp3 into its ith column
    A = cat(2,A,temp3);
end

For The EDIT part:

You can do what you said i.e do prod(B,2) where B is a matrix calculated using above code

Nishant
  • 2,571
  • 1
  • 17
  • 29
  • Not getting the desired results, could you check back again? – Divakar Aug 11 '14 at 14:51
  • @Divakar I cross checked for `n=5; m=6; x= randi(10,5,1)` and my results are same as that of Luis. In you benchmarking answer you have taken `x` as a row vector while OP has taken it as a column vector and so have I. That may be the reason – Nishant Aug 11 '14 at 16:04
  • Exactly that was the issue! Incorporating your results in a bit! – Divakar Aug 11 '14 at 16:04
  • 1
    Congratulations! Yours is the fastest for the main question! You had my +1 already. – Divakar Aug 11 '14 at 16:42
  • @Nishant, fast code, and thanks for the reply to the edit. I don't understand what the code does. Can you explain? – DeltaIV Aug 12 '14 at 10:11
  • @DeltaIV see the explanation, feel free to ask if you still have doubts. – Nishant Aug 12 '14 at 14:27
1

I think the generalized truth table function from the file exchange will help you

try (not tested):

I = ttable(n*ones(1,m));
x(I);
Divakar
  • 218,885
  • 19
  • 262
  • 358
Dan
  • 45,079
  • 17
  • 88
  • 157
0

Okay, there are tricky one-liner solutions for this. There is, however, a file on matlab file exchanging providing the solution that you are looking at if the order doesn't matter. See here combn. It basically makes use of the same tricks as the one-liner presented by other answers. If the order matters, you might have to do some sorting aftwerwards or adjust the source code directly.

Nras
  • 4,251
  • 3
  • 25
  • 37