What is the complexity of the permute() fonction for multiway arrays?

Some functions to manipulate arrays, such as the "reshape" function, seem to run in O(1). On the other hand, the "permute" function runs quite slowly for larger arrays, especially multiway arrays.

Is it possible to bound the numerical complexity of the "permute" function? I guess my question relates to how "permute" is internally implemented.

[James Tursa]

reshape for full variables is probably more like O(0), since it doesn't involve any data movement.

[Walter Roberson]

Constant time is referred to as O(1) by convention. Work is being done, just independent of the size of the input.

... Though didn't you recently point out that some 2D arrays have their size stored differently?

[James Tursa]

I never was any good at telling jokes ... (my point was simply that the time for the reshape is small).

The size difference stuff is in the noise for the purposes of this discussion. The caveat about being "full" is significant, though, since in general reshaping a sparse matrix will almost always involve data movement.

[Bruno Luong]

Waoh, reshape on SPARSE never cross my mind who needs to do this.

But good point James about data movement, but I just think more about it: data under Pr/Pi can stay as it is, just Jc and Ir must be rebuilt, no? So strictly speaking there is no data movement me think.

[James Tursa]

The history of sparse storage makes this reply twofold.

R2008a and earlier:

Shared data copies of sparse matrices only shared the pr/pi pointers ... they did not share the ir/jc indexing. In theory, then, one could reshape with completely new indexing but still share the pr/pi pointers. I don't have an R2008a or earlier version at the moment to check if this does actually happen ... I just point out that it is physically possible. I will check this later.

R2008b and later:

Shared data copies of sparse matrices share pr/pi and ir/jc pointers. So any reshape that causes indexing to change (i.e. new ir/jc data) will require that the pr/pi cannot be shared either ... a deep copy will be required.

In fact, this is the reason my SHAREDCHILD code does not handle certain sparse operations for R2008a and earlier. It would have been extra code to write for 10+ year old versions of MATLAB so I just disallowed it instead of writing the code. I'm too busy writing new code for R2018a complex interleaved memory ...

[Jeremy Cohen]

Actually I can mention that, when computing tensor decompositions such as PARAFAC (a data-mining model), it is standard to reshape and permute arrays several times. Also, Sparse tensors are really a common thing ^^

[Bruno Luong]

I just check with debug format, RESHAPE sparse will result a new Pr, even if nothing is explicit shared.

   >> S = sparse(rand(1,2))

S =

Structure address = d92c4250
m = 2
n = 1
pr = 21159d80

     (1,1)       0.3437
     (1,2)       0.8925

>> S = reshape(S,[2 1])

S =

Structure address = d92c7f90
m = 2
n = 1
pr = 20a03700

     (1,1)       0.3437
     (2,1)       0.8925

Whereas for FULL matrix, Pr address will be unchanged if I do the same reshaping.

I believe for SPARSE, JIT accelerator could optimized further so that

• Pr data is fixed as well as data behind it,
• Ir pointer could be unchanged but new data is inplace updated, but
• Jc pointer must be reallocated and updated.

Obviously TMW fill it's not worth effort to optimize that way for sparse matrix.

Close this out of topic. Sorry.

Answer 1

[Bruno Luong]

The complexity of PERMUTE(A,P) is O(numel(A)) regardless p, but with the constant O depends on p, for example when p=1,2...,d, it will be significantly faster due to linear memory access, but still O(n).

I have no idea how it's implemented internally. There can be so many variants for such problem.

That's why if you can trade permute by reshape, don't hesitate.

Timing and test code

function t = timepermute

p = perms(1:3);
n3 = 10:10:100;
t = zeros(size(p,1),length(n3));
for j = 1:length(n3)
    A=rand(1000,1000,n3(j));
    for i=1:size(p,1)
        pi = p(i,:);
        tic
        B = permute(A,pi);
        t(i,j) = toc;
    end
end
s = arrayfun(@(i) sprintf('p = %s',mat2str(p(i,:))), 1:size(p,1),'unif',0);
close all
plot(n3*size(A,1)*size(A,2),t);
xlabel('n');
ylabel('permute time [s]');
legend(s,'location','northwest')
grid on

Answer 2

[Steven Lord]

Calling reshape on an array doesn't change the order of the elements, just the size/shape of the array, so there's no need to move the elements around in memory.

Calling permute needs to change the order of the elements (unless the permutation vector starts with 1:ndims(arrayBeingPermuted)) so there are memory manipulation operations involved.

>> r = reshape(1:10000, [100 100]);
>> r2 = reshape(r, [50 200]);
>> r(4836), r2(4836)
ans =
        4836

ans =
        4836

>> p2 = permute(r, [2 1]);
>> r(4836), p2(4836)
ans =
        4836

ans =
        3549

If you're trying to compare the performance of those two operations, you're comparing apples and oranges. Actually, a better analogy would be asking your friends to call you by a different name (low effort, reshape) versus filing the necessary paperwork with the government to change your name legally on your driver's license, passport, and other official documents (higher effort, permute.)

[Jeremy Cohen]

I completely agree with what you just explained, which is why I mentionned reshape is O(1). I'm just curious how the "permute" operator scales with the dimensions of the arrays.

[Walter Roberson]

permute() touches each input entry exactly once. It is not implemented as a multiplication by a permutation matrix.

Answer 3

[Walter Roberson]

permute is probably O(n) but is probably badly affected by cache line effects.

[Jeremy Cohen]

I'm sorry but this is not quite precise. Given a multidimensional array of dimensions n_1, n_2, ..., n_d, do you mean the complexity is probably around O(\prod_i {n_i}) ?

[Walter Roberson]

O(number of elements of input matrix)