Using cumulative sum to create a symmetric n-dimensional array

I am trying to create a symmetric tensor of order 2, 3 and more using the below code.

n = 5;
mask = tril(ones(n));
posVecF = 0;
Order = 3; %order of tensor 
for j = 2:Order
    A0 = bsxfun(@times,ones(repelem(n,j)),mask);
    posMatrix = cumsum(A0(:));
    posMatrix = reshape(posMatrix,repelem(n,j));
    posMatrix(A0==0) = 0; % Make lower triangular tensor
    posMatrix = posMatrix+posMatrix'-diag(diag(posMatrix)); % Mirror the lower half of the triangular tensor
end

The expected and calculated 'posMatrix' for order 2 (j=2) using the below code is

1	2	3	4	5
2	6	7	8	9
3	7	10	11	12
4	8	11	13	14
5	9	12	14	15

For order 3 (j = 3), I expect my 'posMatrix' to be a 5x5x5 tensor

val(:,:,1) =
1     2     3     4     5
2     6     7     8     9
3     7    10    11    12
4     8    11    13    14
5     9    12    14    15
val(:,:,2) =
2     6     7     8     9
6    16    17    18    19
7    17    20    21    22
8    18    21    23    24
9    19    22    24    25
val(:,:,3) =
3     7    10    11    12
7    17    20    21    22
10    20    26    27    28
11    21    27    29    30
12    22    28    30    31
val(:,:,4) =
4     8    11    13    14
8    18    21    23    24
11    21    27    29    30
13    23    29    32    33
14    24    30    33    34
val(:,:,5) =
5     9    12    14    15
9    19    22    24    25
12    22    28    30    31
14    24    30    33    34
15    25    31    34    35

But this is not what I get. There is an issue in the line posMatrix = cumsum(A0(:)); for higher order and posMatrix = posMatrix+posMatrix'-diag(diag(posMatrix)); should be adjusted for any orders. I am struggling to find a solution for this problem here.

Answer 1

[Shubham]

Hi Neuling,

To create a symmetric tensor of orders 2, 3, and higher, adjustments are needed in your code to correctly handle the symmetry across all dimensions and orders. The primary challenges include:

1. Generalizing Symmetry for Higher Orders: The operation posMatrix = posMatrix+posMatrix'-diag(diag(posMatrix)); is specific to matrices (2D arrays) and does not generalize to tensors of higher dimensions.
2. Correct Cumulative Sum for Positioning: The line posMatrix = cumsum(A0(:)); computes the cumulative sum linearly, which doesn't account for the multidimensional structure properly when reshaped back for higher orders.

To address these issues, let's revise your approach. We need a method that:

• Constructs the lower triangular part correctly for any order of tensor.
• Mirrors this lower triangular part across all dimensions to achieve symmetry.

Here's an approach to achieve this:

function posTensor = createSymmetricTensor(n, Order)
    % Initialize the tensor with zeros
    posTensor = zeros(repelem(n, Order));
    
    % Generate indices for the lower triangular part (including diagonal)
    ind = tril(true(n));
    [row, col] = find(ind);
    
    % Initialize position counter
    posCounter = 1;
    
    % Iterate through each element in the lower triangular part to assign positions
    for i = 1:length(row)
        % Generate a template tensor for current position
        tempTensor = zeros(n);
        tempTensor(row(i), col(i)) = 1;
        
        % Expand the template tensor to higher dimensions if needed
        if Order > 2
            tempTensor = repmat(tempTensor, repelem(1, Order-2));
        end
        
        % Find all permutations of the indices to mirror the positions across dimensions
        permsRow = perms([row(i)*ones(1,Order-2), row(i), col(i)]);
        uniquePermsRow = unique(permsRow, 'rows');
        
        % Assign positions to all symmetric locations
        for j = 1:size(uniquePermsRow, 1)
            indices = num2cell(uniquePermsRow(j, :));
            posTensor(indices{:}) = posCounter;
        end
        
        posCounter = posCounter + 1;
    end
end

This function attempts to create a symmetric tensor by explicitly handling the symmetry for tensors of any order. However, note that the approach to generating all permutations and then finding unique rows can be computationally intensive for large tensors or higher orders. The approach needs to be optimized for specific needs or constraints.

To test this function for orders 2 and 3:

n = 5;
Order = 2; % Change this to 3, 4, etc., to generate tensors of different orders
posTensor = createSymmetricTensor(n, Order);
disp(posTensor);

Please note: The provided solution is a conceptual approach and might require adjustments based on the exact requirements for symmetry in higher dimensions. The complexity of generating symmetric tensors increases significantly with the order, and special attention is needed to ensure all symmetric positions are correctly identified and assigned.