Symbolic multidimensional arrays with changable dimension

Consider an integer n, p a 1x4 symbolic array and c a 4x4 symbolic matrix. Here Latin letters are denoting integer indexes that run from 1 to 4. I would like to compute the following n-multidimensional array (let us call it M)

.

How can I compute this quantity?

Here is an example. If we have . We can solve this with a for loop such that

M = sym(zeros(4, 4, 4));
syms x y z w
p = [x y z w];
c = [x y z w; x y z w; x y z w; x y z w];
for i = 1:4
    for j = 1:4
        for k = 1:4
            M(i,j,k) = p(i)*c(i,j)*c(j,k);
        end
    end
end

How can I do this symbolically when n is not known a priori, but it is choosen by the user?

Thank you.

[Bruno Luong]

Your decription is not clear.

What are indexes from 1 to 4? a, b, c, d, ...? What are the size of the arrays p and c?

Then you did not tell how the product is related to M. Then there are new indexes a1, a2, ..., an. What are the relation the new indexes with a, b, c, d, ...

A small example would not be useless.

[Simone]

Thank you. I have modified my question to improve clarity.

[Bruno Luong]

Still a bad description, Should read:

"Consider an integer n, p a 1x4 array and c a 4x4 matrix. Here Latin letters are denoting integer indexes that run from 1 to 4. I would like to compute the following n-multidimensional array (let us call it M)"

Answer 1

[Bruno Luong]

p = [1 2 3 4];
c = [1 2 3 4; 1 2 3 4; 1 2 3 4; 1 2 3 4];

n=3;
M = p(:); % EDIT, bug fix
for k=1:n-1
    cc = reshape(c, [ones(1,k-1) size(c)]);
    M = M.*cc;
end
M

M = 
M(:,:,1) =

     1     2     3     4
     2     4     6     8
     3     6     9    12
     4     8    12    16


M(:,:,2) =

     2     4     6     8
     4     8    12    16
     6    12    18    24
     8    16    24    32


M(:,:,3) =

     3     6     9    12
     6    12    18    24
     9    18    27    36
    12    24    36    48


M(:,:,4) =

     4     8    12    16
     8    16    24    32
    12    24    36    48
    16    32    48    64

[Simone]

Is there a reason for this not working when c is initialized in the same way but with symbolic expressions?

[Bruno Luong]

Yes it uses auto-expansion not available in symbolic I guess (I don't have symbolic toolbox).

[Simone]

Thank you for your answer, it provides a useful insight. However, I have to do this symbolically. I updated the question to specify this point.

[Bruno Luong]

To avoid auto-expansion you could use repmat(), not sure if symbolic support thos command

p = [1 2 3 4];
c = [1 2 3 4; 1 2 3 4; 1 2 3 4; 1 2 3 4];

n = 3;
M = p;
q = length(M);
m = 1;
cc = reshape(c,[1 q q]);
for k=1:n-1
    M = reshape(M,[m,q]);
    M = repmat(M,[1 1 q]) .* repmat(cc,[m 1 1]);
    m = m*q;
end
if n == 1 % EDIT to make it works for n=1
    sz = [q 1];
else
    sz = q+zeros(1,n);
end
M = reshape(M, sz);
M

M = 
M(:,:,1) =

     1     2     3     4
     2     4     6     8
     3     6     9    12
     4     8    12    16


M(:,:,2) =

     2     4     6     8
     4     8    12    16
     6    12    18    24
     8    16    24    32


M(:,:,3) =

     3     6     9    12
     6    12    18    24
     9    18    27    36
    12    24    36    48


M(:,:,4) =

     4     8    12    16
     8    16    24    32
    12    24    36    48
    16    32    48    64

[Simone]

Thank you. I checked your solution and it works with symbolic arrays and matrices, except for for which it gives the error "New array must have the same number of elements as the original array.".

[Bruno Luong]

It works for me see the code I give and RUN.

Something you change that is not right.

[Simone]

I copied and pasted your code, clearing my workspace. If I replace with it gives me "Number of elements must not change. Use [] as one of the size inputs to automatically calculate the appropriate size for that dimension." as an error using reshape.

Answer 2

[Walter Roberson]

You have the implicit question of how you can do n nested loops where n is entered by the user.

There are two major ways:

First way:

You can use ndgrid to build cell arrays of elements that will be used in the calculations, or cell arrays of indices. The number of entries build can be dynamic:

n = 3;
T1 = cell(1,n);
[T1{:}] = ndgrid(sym(1:4));  %these will be individually multidimensional arrays
T2 = cellfun(@(M) M(:), T1, 'uniform', 0);  %these will be individually column vectors
T2 = horzcat(T2{:});

whos T2

  Name       Size            Bytes  Class    Attributes

  T2        64x3                 8  sym                

T2 will now be an array in which each row contains a combination of elements, and all possible combinations are listed. You could, for example, take

prod([p(T2(:,1)), c(T2(:,2:end))], 2)

... except if you are indeed going to use the T2 values as indices, you would not use the sym() in the ndgrid call.

This approach requires storing those intermediate multidimensional arrays, and that can add up to a lot of temporary memory.

Second way:

You can use what has sometimes been called an "odometer" pattern of programming. I gave general code for that in https://www.mathworks.com/matlabcentral/answers/623358-get-a-combination-of-unique-paths-for-given-pair-of-numbers#comment_1082638 -- and in there I do not assume that the values being selected from have the same datatype (you can be a bit more efficient if the values are all numeric.)

The odometer programming pattern at each point changes the last element of a vector to the next value. If it was already the last value possible, it rolls that position over to the first value for it, and then moves to the left and increments that value; if that was not already the last then it is done. If that one was the last possible for its position, it rolls it over to the first valid value for that position, and moves one left... and so on.

This only requires keeping storage for the current configuration, together with a vector that defines the valid values at each position. It does not require keeping multidimensional arrays.

The odometer pattern is not as fast as vectorization... when you have enough memory to do the vectorization. When you do not have enough memory to do the vectorization, then the odometer pattern is faster, in the sense that clearly eventually it will get through all of the patterns, and getting through all the patterns at some point is "faster" than not being able to get through all of the possibilities because you ran out of memory.

[Simone]

Thank you for your answer. I understand that your first method gives all the possible combinations of indexes that would come from an n-nested loop. However, not only the number of for loop is parametrized by n, but also the number of multiplied c matrices does (we have actually matrices multiplied with that particular choice of indexing).

How do I compute M with this method? I tried to apply your code but

prod([p(T2(:,1)), c(T2(:,2:end))], 2)

gives me "CAT arguments dimensions not consistent.".