Higher order Tensor computations

Question

I have a $D$-way tensor H of dimension $I \times I \times \dots \times I$ ($D$ times), that represent the coefficients of a polynomial. For better understanding, I provided an image of 3-way tensor and explained the image in the description below.

H is a 3-dimensional tensor of size 4x4x4 representing coefficients of a 3rd order polynomial. The element (111) gives the constant term in the polynomial, the sum of vectors (2:end,1,1), (1,2:end,1)', reshape((1,1,2:end),[],1) give the coefficients of first order monomials in the polynomial and the slices (2:end,2:end,1), (2:end,1,2:end), (1,2:end,2:end) give the coefficients of 2nd order polynomial and a sub-tensor of size 3x3x3 inside will give the coefficients of 3rd order monomials.

I am only interested in extracting the coefficients of nonlinear monomials of order (2,3,...,D) and arranging them in a vector of particular fashion. I am providing a following code snippet, where I do the computation for 3rd order tensor. But It has nested for loops and the number of for loops increases as $D$ increases, which is not desired. Can this computation be simplified using vectorization or any other tensor operation?

    pos = 1; start = 2;

if Memory > 0
    if Order > 1
        for i = start:I
            for j = i:I
                indices = num2cell(unique(perms([j,i,ones(1,Order-2)]),'rows'));
                [F,pos] = calcF(indices,F,pos,H);
            end
        end
    end
    if Order > 2
        for k = start:I
            for i = k:I
                for j = i:I
                    indices = num2cell(unique(perms([j,i,k,ones(1,Order-3)]),'rows'));
                    [F,pos] = calcF(indices,F,pos,H);
                end
            end
        end
    end
    if Order > 3
        for m = start:I
            for k = m:I
                for i = k:I
                    for j = i:I
                        indices = num2cell(unique(perms([j,i,k,m,ones(1,Order-4)]),'rows'));
                        [F,pos] = calcF(indices,F,pos,H);
                    end
                end
            end
        end
    end
    if Order > 4
        for n = start:I
            for m = n:I
                for k = m:I
                    for i = k:I
                        for j = i:I
                            indices = num2cell(unique(perms([j,i,k,m,n,ones(1,Order-5)]),'rows'));
                            [F,pos] = calcF(indices,F,pos,H);
                        end
                    end
                end
            end
        end
    end
end

function [F,pos] = calcF(indices,F,pos,H)
    for l = 1:size(indices,1)
        F(pos) = F(pos) + H(indices{l,:});
    end
    pos = pos + 1;
end

I am trying to sum the elements in a position in H and all possible unique permutations of this positional index and rearranging them in a vector form. The vector first contains the coefficients of 2nd order monomial, followed by 3rd order and so on. How can I extend this code for 4th and higher order polynomials without having to write many for loops?

Example: to find the coefficients of these monomials $\begin{bmatrix} x_1^2 & x_1x_2 & x_1u & x_2^2 & x_2u & u^2 & x_1^3 & x_1^2x_2 & x_1^2u & x_1x_2^2 & x_1x_2u & x_1u^2 & x_2^3 & x_2^2u & x_2u^2 & u^3\end{bmatrix}$, the following indices of H are summed up:

$$ \begin{bmatrix}H(221) + H(212) + H(122)\\ H(231) + H(321) + H(123) + H(132) + H(213) + H(312)\\ H(241) + H(124) + H(142) + H(214) + H(421) + H(412)\\ H(331) + H(313) + H(133)\\ H(341) + H(143) + H(134) + H(413) + H(431) + H(314)\\ H(441) + H(414) + H(144)\\ H(222)\\ H(223) + H(322) + H(232)\\ H(224) + H(422) + H(242)\\ H(233) + H(323) + H(332)\\ H(234) + H(342) + H(432) + H(324) + H(243) + H(423)\\ H(244) + H(442) + H(424)\\ H(333)\\ H(334) + H(343) + H(433)\\ H(344) + H(434) + H(443)\\ H(444)\end{bmatrix}$$

where the first 6 coefficients are of the 2nd order monomials and the next 10 are of 3rd order monomials.

For a simple 2 Dimensional case, it looks as below: $$y = \begin{pmatrix}\begin{bmatrix} 1 \\ x_1 \\ x_2 \\ u \end{bmatrix} \otimes \begin{bmatrix} 1 \\ x_1 \\ x_2 \\ u \end{bmatrix}\end{pmatrix}* h$$

$$y = rowVec\begin{pmatrix}\begin{bmatrix}1 & x_1 & x_2 & u \\ x_1 & x_1^2 & x_1x_2 & x_1u \\ x_2 & x_1x_2 & x_2^2 & x_2u \\ u & x_1u & x_2u & u^2 \end{bmatrix}\end{pmatrix} * colVec\begin{pmatrix}\begin{bmatrix}h_0 & h_1(0) & h_1(1) & h_1(2) \\ h_1(0) & h_2(00) & h_2(01) & h_2(02) \\ h_1(1) & h_2(10) & h_2(11) & h_2(12) \\ h_1(2) & h_2(20) & h_2(21) & h_2(22) \end{bmatrix}\end{pmatrix} $$ The $h$ here is a column vector, that can be reshaped into a $D$ dimensional tensor, where the sum of the elements h(2,3) and h(3,2) gives the coefficient of monomial $x_1x_2$.

I am looking for possible ways to make the code given work for any given Order. For now, the code works for Order = 5. For higher order, I have to add more nested for loops. But this is not an optimal solution.