MATLAB: efficient construction of large matrices using elements of a small matrix

Question

I have a small matrix A of size (n+1)x(n+1) with n=O(10). I need to build N=(n+1)^m much larger matrices B_i, each one of size Nxm, where m=O(10) also. Luckily, I need to store just one B_i at the time, so memory is not an issue. To generate each B_i I have a matrix I of indices, of size Nxm. I contains all the possible m-ples of integers between 1 and n+1. For example, let n=4,m=4, then

I=[1 1 1 1; 
   2 1 1 1; 
   3 1 1 1; 
   4 1 1 1; 
   5 1 1 1; 
   1 2 1 1; 
       .
       .
   5 5 5 5]

The indices in row i of I are the column indices in A of the elements of the corresponding column of B_i. For example, consider i=3, i.e., B_3=C. Then column 1 of matrix C is assembled using elements from column 3 of A, while columns 2 to 4 use elements from column 1 of A. In the same way, indices from column j of I are the row indices in A of the elements of column j of B_i. Consider again B_3=C:

C(1,1) = A(1,3)
C(2,1) = A(2,3)
C(3,1) = A(3,3)
C(4,1) = A(4,3)
C(5,1) = A(5,3)
C(6,1) = A(1,3)
       .
       .
C(N,1) = A(5,3)

C(1,2) = A(1,1)
C(2,2) = A(1,1)
C(3,2) = A(1,1)
C(4,2) = A(1,1)
C(5,2) = A(1,1)
C(6,2) = A(2,1)
       .
       .
C(N,2) = A(5,1)

And so on. It's trivial to build each B_i with a double for loop. However, since I have to build N B_i, the resulting triple for loop takes forever. Can you show me some clever MATLAB trick to build these matrices super-quickly?

By the way, I don't really need the the B_i, but just the elementwise product of their columns, i.e.

PolyMD=prod(C,2)

which is an Nx1 column vector. As I said before, I store just one B_i at a time, so computing B_i first and then PolyMD is not a big issue. Still, if there's a fast and simple way to get PolyMD directly, I'd love to know.

EDIT: here is a full example to hopefully make the question more clear.

% input
n=2;
m=3;
A=[1.000000000000000e+00    -1.732050807568877e+00     1.414213562373095e+00; 
   1.000000000000000e+00     1.160285448625519e-16    -7.071067811865475e-01;
   1.000000000000000e+00     1.732050807568877e+00     1.414213562373095e+00]; 
I=[     1     1     1;
        2     1     1;
        3     1     1;
        1     2     1;
        2     2     1;
        3     2     1;
        1     3     1;
        2     3     1;
        3     3     1;
        1     1     2;
        2     1     2;
        3     1     2;
        1     2     2;
        2     2     2;
        3     2     2;
        1     3     2;
        2     3     2;
        3     3     2;
        1     1     3;
        2     1     3;
        3     1     3;
        1     2     3;
        2     2     3;
        3     2     3;
        1     3     3;
        2     3     3;
        3     3     3];

% desired output for B_2 (I skip B_1 since it's trivial, and B_3,....B_27 for brevity!)    
B_2=[ -1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
     1.160285448625519e-16     1.000000000000000e+00     1.000000000000000e+00;
     1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
    -1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
     1.160285448625519e-16     1.000000000000000e+00     1.000000000000000e+00;
     1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
    -1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
     1.160285448625519e-16     1.000000000000000e+00     1.000000000000000e+00;
     1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
    -1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
     1.160285448625519e-16     1.000000000000000e+00     1.000000000000000e+00;
     1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
    -1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
     1.160285448625519e-16     1.000000000000000e+00     1.000000000000000e+00;
     1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
    -1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
     1.160285448625519e-16     1.000000000000000e+00     1.000000000000000e+00;
     1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
    -1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
     1.160285448625519e-16     1.000000000000000e+00     1.000000000000000e+00;
     1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
    -1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
     1.160285448625519e-16     1.000000000000000e+00     1.000000000000000e+00;
     1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
    -1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00;
     1.160285448625519e-16     1.000000000000000e+00     1.000000000000000e+00;
     1.732050807568877e+00     1.000000000000000e+00     1.000000000000000e+00];

It's not clear in your description how `B_i` differs from `B_j`. Your construction of `C` does not depend on `i` — Luis Mendo, Aug 31 '14 at 15:18
C depends on i because the i-th row of I gives you the column indices in A. Example with B_2 and B_3: `B_2(1,1)=A(1,2)` while `B_3(1,1)=A(1,3)`, and that's because the second row of I is [2 1 1 1], while the third row is [3 1 1 1]. — DeltaIV, Aug 31 '14 at 15:21
I think I begin to understand what you want. But if `A` is m x (n+1), so 4x5 in your example, `A((5,1)` is not defined. What's the size of `A`? — Luis Mendo, Aug 31 '14 at 15:40
It seems like `A` has to be (at least) `(n+1)` x `(n+1)`, because for example `B_5(5,1)` is defined to be `A(5,5)`. — Luis Mendo, Aug 31 '14 at 15:49
Sorry.@LuisMendo ! I made a typing mistake. As you correctly noticed, `A` is `(n+1)` x `(n+1)`. Need some sleep now, but tomorrow I'll check your solution and if doesn't work, I'll post a full example. Thanks a lot — DeltaIV, Aug 31 '14 at 21:44

Luis Mendo · Accepted Answer · 2014-09-01T18:56:36.123

2

Maybe this is what you want:

n = 4;
m = 4;
N = (n+1)^m;
A = rand(n+1,n+1); %// example data
I = fliplr(dec2base(0:(n+1)^m-1,n+1)-'0'+1); %// combinations of column indices
for ii = 1:N %// it's better not to use "i" as a variable name
    B_i = A(bsxfun(@plus, I, (I(ii,:)-1)*(n+1))); %// linear indexing into A
    result_i = prod(B_i,2); %// "Poly" is a reserved word. Better not use it.
end

edited Sep 01 '14 at 18:56

answered Aug 31 '14 at 15:21

Luis Mendo

110,752
13
76
147

Hi @LuisMendo. I have already the matrix I of indices. BTW, you and others helped me generate it in [link] (http://stackoverflow.com/questions/25245105/matlab-create-a-large-matrix-by-repeating-elements-of-a-vector-with-increasing). Can I skip the fifth line in your code, then? – DeltaIV Aug 31 '14 at 21:50
Sure. I used that line only to generate an example for me to work with – Luis Mendo Aug 31 '14 at 22:11
Nope, sorry, it doesn't work. As promised, I edited my question with a worked out example. Thanks again for the help and the patience! – DeltaIV Sep 01 '14 at 12:56
@DeltaIV Try replacing `J` by `I` (see edited answer) – Luis Mendo Sep 01 '14 at 13:30
excellent! Now it works perfectly. I gained a 3x factor in speed (which may be even more for larger problem). Thanks a lot! Just a last question: I understand that bsxfun is used to add the row vector `(I(ii,:)-1)*(n+1)` to the matrix I. Why does this return the right elements from A? Is there a simple, intuitive explanation? – DeltaIV Sep 01 '14 at 15:12
1

@DeltaIV `I` is the row index you want to use, whereas the column index is `I(ii,:)`. Those two indices are combined with `bsxfun` into a [linear index](http://www.mathworks.es/es/help/matlab/math/matrix-indexing.html) that allows accessing the desired elements of `A`. So `bsxfun` here does the same as `sub2ind`, but it's faster for two reasons: 1) It avoids some overhead `sub2ind` has. 2) `sub2ind` would have required `I(ii,:)` to be `repmat`ted to match the size of `I` – Luis Mendo Sep 01 '14 at 15:40
1

@DeltaIV As to why I use `(I(ii,:)-1)*(n+1))`: matrix elements in Matlab are stored in [column-major order](http://en.wikipedia.org/wiki/Row-major_order). So the `(...-1)*(n+1)` when creating the linear index means "jump to another column", where `n+1` is the numer of rows – Luis Mendo Sep 01 '14 at 15:44
@DeltaIV Glad I explained myself! Linear indexing in Matlab is actually a very powerful tool – Luis Mendo Sep 01 '14 at 18:55

MATLAB: efficient construction of large matrices using elements of a small matrix

1 Answers1