How to extract diagonal elements of multidimensional array ?

If I have a m-order n-dimensional tensor. How should I extract the diagonal elements ?

For example

% Generate a 3-order 4-dimensional tensor
rng('default')
A = rand(4,4,4)

A = 
A(:,:,1) =

    0.8147    0.6324    0.9575    0.9572
    0.9058    0.0975    0.9649    0.4854
    0.1270    0.2785    0.1576    0.8003
    0.9134    0.5469    0.9706    0.1419


A(:,:,2) =

    0.4218    0.6557    0.6787    0.6555
    0.9157    0.0357    0.7577    0.1712
    0.7922    0.8491    0.7431    0.7060
    0.9595    0.9340    0.3922    0.0318


A(:,:,3) =

    0.2769    0.6948    0.4387    0.1869
    0.0462    0.3171    0.3816    0.4898
    0.0971    0.9502    0.7655    0.4456
    0.8235    0.0344    0.7952    0.6463


A(:,:,4) =

    0.7094    0.6551    0.9597    0.7513
    0.7547    0.1626    0.3404    0.2551
    0.2760    0.1190    0.5853    0.5060
    0.6797    0.4984    0.2238    0.6991

The diagonal elements are A(1,1,1) = 0.8147, A(2,2,2) = 0.0357, A(3,3,3) = 0.7655 and A(4,4,4) = 0.6991.

I was hoping to have a tensor_diag function that takes a tensor A as an input parameter and returns a vector consisting of its diagonal elements.

[Mattia Marsetti]

Ciao shuang Yang

I don't know of any native Matlab function for this purpose. You can try to use the following custom function to solve your problem:

function out=get_tensor(v)
size_v=size(v);
if sum(size_v == size_v(1))<numel(size_v)
    error('the input vector is not a sqare matrix');
end

N = size_v(1);
out=zeros(N,1);
for i=1:N
    str='v(';
    for s=1:N
        str=[str 'i,'];
    end
    str(end:end+1)=');';
    out(i)=eval(str);
end
end

[Stephen23]

@Mattia Marsetti: your code throws an error on the example array:

rng('default')
A = rand(4,4,4)

A = 
A(:,:,1) =

    0.8147    0.6324    0.9575    0.9572
    0.9058    0.0975    0.9649    0.4854
    0.1270    0.2785    0.1576    0.8003
    0.9134    0.5469    0.9706    0.1419


A(:,:,2) =

    0.4218    0.6557    0.6787    0.6555
    0.9157    0.0357    0.7577    0.1712
    0.7922    0.8491    0.7431    0.7060
    0.9595    0.9340    0.3922    0.0318


A(:,:,3) =

    0.2769    0.6948    0.4387    0.1869
    0.0462    0.3171    0.3816    0.4898
    0.0971    0.9502    0.7655    0.4456
    0.8235    0.0344    0.7952    0.6463


A(:,:,4) =

    0.7094    0.6551    0.9597    0.7513
    0.7547    0.1626    0.3404    0.2551
    0.2760    0.1190    0.5853    0.5060
    0.6797    0.4984    0.2238    0.6991

get_tensor(A)

Index in position 4 exceeds array bounds. Index must not exceed 1.

Error in solution>get_tensor (line 19)
out(i)=eval(str);

function out=get_tensor(v)
size_v=size(v);
if sum(size_v == size_v(1))<numel(size_v)
    error('the input vector is not a sqare matrix');
end

N = size_v(1);
out=zeros(N,1);
for i=1:N
    str='v(';
    for s=1:N
        str=[str 'i,'];
    end
    str(end:end+1)=');';
    out(i)=eval(str);
end
end

Note that you could easily replace the evil EVAL with a cell array and a comma-separated list.

[Mattia Marsetti]

Thanks for the hint Stephen23, I'll try that

Regards

Answer 1

[the cyclist]

rng('default')

N = 4;
A = rand(N,N,N);

A(1:N^2+N+1:end)

ans = 1×4

0.8147 0.0357 0.7655 0.6991

[Bruno Luong]

Note that cyclist's indexing (1:N^2+N+1:end) is

sub2ind([N,N,N],1:N,1:N,1:N)

[the cyclist]

Nice -- and more intuitive about what is going on

[shuang Yang]

Thank you very much for your answer, although this is not the correct answer, but I should have guessed how to write it, the correct version of it should be like this

rng('default')
%% order = 2
N = 3;
A = rand(N,N)

A = 3×3

0.8147 0.9134 0.2785 0.9058 0.6324 0.5469 0.1270 0.0975 0.9575

A(1:N+1:end)

ans = 1×3

0.8147 0.6324 0.9575

%% order = 3
N = 4;
A = rand(N,N,N) % 

A = 
A(:,:,1) =

    0.9649    0.4854    0.9157    0.0357
    0.1576    0.8003    0.7922    0.8491
    0.9706    0.1419    0.9595    0.9340
    0.9572    0.4218    0.6557    0.6787


A(:,:,2) =

    0.7577    0.1712    0.0462    0.3171
    0.7431    0.7060    0.0971    0.9502
    0.3922    0.0318    0.8235    0.0344
    0.6555    0.2769    0.6948    0.4387


A(:,:,3) =

    0.3816    0.4898    0.7547    0.1626
    0.7655    0.4456    0.2760    0.1190
    0.7952    0.6463    0.6797    0.4984
    0.1869    0.7094    0.6551    0.9597


A(:,:,4) =

    0.3404    0.2551    0.9593    0.2575
    0.5853    0.5060    0.5472    0.8407
    0.2238    0.6991    0.1386    0.2543
    0.7513    0.8909    0.1493    0.8143

A(1:N^2+N+1:end)

ans = 1×4

0.9649 0.7060 0.6797 0.8143

%% order = 4
N = 4;
A = rand(N,N,N,N);
A(1:N^3+N^2+N+1:end)

ans = 1×4

0.2435 0.3692 0.0292 0.6569

A(1)

ans = 0.2435

A(2,2,2,2)

ans = 0.3692

A(3,3,3,3)

ans = 0.0292

%% order = 5
N = 7;
A = rand(N,N,N,N,N);
A(1:N^4+N^3+N^2+N+1:end)

ans = 1×7

0.6280 0.7118 0.6596 0.5435 0.6933 0.9827 0.9083

[the cyclist]

Ah, I read your question too quickly, and didn't make my solution general enough. Glad you found it.

[Stephen23]

@shuang Yang: you could generalize that:

A(1:sum(N.^(0:ndims(A)-1):end)

Answer 2

[Bruno Luong]

N = 7;
A = rand(N,N,N,N,N);

p=ndims(A);
N=length(A);

% Method 1: generalization of cyclist's answer
step = polyval(ones(1,p),N);
idx = 1:step:N^p;
A(idx)

ans = 1×7

0.4728 0.2804 0.8097 0.9442 0.5680 0.1666 0.9129

% Method 2
c = repmat({1:N}, [1,p]);
idx = sub2ind(size(A), c{:});
A(idx)

ans = 1×7

0.4728 0.2804 0.8097 0.9442 0.5680 0.1666 0.9129