title: MATLAB: efficient construction of large matrices using elements of a small matrix
body: <p>I have a small matrix <code>A</code> of size <code>(n+1)</code>x<code>(n+1)</code> with <code>n</code>=<code>O(10)</code>. I need to build <code>N=(n+1)^m</code> much larger matrices <code>B_i</code>, each one of size <code>N</code>x<code>m</code>, where <code>m</code>=<code>O(10)</code> also.  Luckily, I need to store just one <code>B_i</code> at the time, so memory is not an issue. To generate each <code>B_i</code> I have a matrix <code>I</code> of indices, of size <code>N</code>x<code>m</code>. <code>I</code> contains all the possible <code>m</code>-ples of integers between <code>1</code> and <code>n+1</code>. For example, let <code>n=4</code>,<code>m=4</code>, then</p>
  <pre><code>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] 
  </code></pre>
  <p>The indices in row <code>i</code> of <code>I</code> are the column indices in <code>A</code> of the elements of the corresponding column of <code>B_i</code>. For example, consider <code>i=3</code>, i.e., <code>B_3=C</code>. Then column <code>1</code> of matrix <code>C</code> is assembled using elements from column <code>3</code> of <code>A</code>, while columns <code>2</code> to <code>4</code> use elements from column <code>1</code> of <code>A</code>. In the same way, indices from column <code>j</code> of <code>I</code> are the row indices in <code>A</code> of the elements of column <code>j</code> of <code>B_i</code>. Consider again <code>B_3=C</code>:</p>
  <pre><code>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)
  </code></pre>
  <p>And so on. It's trivial to build each <code>B_i</code> with a double for loop. However, since I have to build <code>N</code> <code>B_i</code>, the resulting triple <code>for</code> loop takes forever. Can you show me some clever MATLAB trick to build these matrices super-quickly?</p>
  <p>By the way, I don't really need the the <code>B_i</code>, but just the elementwise product of their columns, i.e.</p>
  <pre><code>PolyMD=prod(C,2)
  </code></pre>
  <p>which is an <code>N</code>x<code>1</code> column vector. As I said before, I store just one <code>B_i</code> at a time, so computing <code>B_i</code> first and then <code>PolyMD</code> is not a big issue. Still, if there's a fast and simple way to get <code>PolyMD</code> directly, I'd love to know.</p>
  <p>EDIT: here is a full example to hopefully make the question more clear.</p>
  <pre><code>% 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];
  </code></pre>