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. Thanks a lot and regards,

Sergio

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
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
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
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
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

Luis Mendo Luis Mendo · Accepted Answer · 2014-08-31T15:21:17

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

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

1 Answers