How can I get randperm to return a permutation of a vector that has no entries at their original positions?

I want take a random permutation of a vector such that all entries of the vector move to a new location.

For example, if I have a vector [1,2,3,4,5], then the following permutations are acceptable:

[2,1,4,5,3], [3,1,5,2,4], [5,4,2,3,1], etc.

However, for me, the following vector is not acceptable:

[2,4,3,5,1]

because the "3" has remained in the same location.

The "randperm" function in MATLAB allows for some of the entries in the vector to stay in the same position. Is there some way to use randperm that stops it from doing this? Or is there some other function out there that I am missing? (I have also looked at the functions "datasample" and "randsample" but they also do not seem to allow for this).

[Paul]

Is there an upper bound on the length of the input vector?

[Stephen23]

This type of permutation is called a derangement:

https://en.wikipedia.org/wiki/Derangement

You can start by searching FEX:

https://www.mathworks.com/matlabcentral/fileexchange?q=derangement

and reading the descriptions of those submissions.

Answer 1

[Bruno Luong]

This will do what you ask for

https://www.mathworks.com/matlabcentral/fileexchange/30189-randpermfull?s_tid=srchtitle

Even better but need MEX

https://www.mathworks.com/matlabcentral/fileexchange/27076-shuffle?s_tid=srchtitle

Answer 2

[Bruno Luong]

Here is an implementation of a non-rejection method and unbiased random derangement:

function p = randder(n)
% p = randder(n)
% Generate a random derangement if length n
%
% INPUT:
%   n: scalar integer >= 2
% OUTPUT:
%   p: array of size (1 x n), such that 
%        unique(p) is equal to (1:n)
%        p(i) ~= i for all i = 1,2,....n.
%
% Base on: "Generating Random Derangement", Martinez Panholzer, Prodinger
% Remove the need inner loop by skrink index table J (still not ideal)
%
% See also: randperm

p = 1:n;
b = true(1,n);
m = n-1;
J = 1:m;
i = n;
u = n;
utab = 1:n;
qtab = (utab-1).*subfactorial(utab-2)./subfactorial(utab);
overflowed = ~isfinite(qtab);
qtab(overflowed) = 1./utab(overflowed);
x = rand(1,n);
r = rand(1,n);
while u>=2
    if b(i)
        k = ceil(x(i)*m);
        j = J(k);
        p([i j]) = p([j i]);
        if r(i) < qtab(u)
            b(j) = false;
            J(k:m-1) = J(k+1:m);
            m = m-1;
            u = u-1;
        end
        u = u-1;
    end
    i = i-1;
    if J(m)==i
        m = m-1;
    end
end
end % randder

%%
function D = subfactorial(n)
D = floor((gamma(n+1)+1)/exp(1));
end

It might be slower but it's a non rejection method, so at least the run-time is always predictable.

Test:

p=randder(10)

p =

     7     6     1     3     9     5    10     4     2     

Check validity and uniformity

m = 100000;
n = 6;
D=arrayfun(@(x) randder(n), 1:m, 'UniformOutput', false);

D=cat(1,D{:});
[U,~,J]=unique(D,'rows');

OK = all(U-(1:n),'all') && ...
     ~any(sort(U,2)-(1:n),'all'); % must be true
if OK
    fprintf('All derangements are valid\n');
end

% Check uniformity
close all
nbins = min(1000,size(U,1));
histogram(J,nbins);

[Bruno Luong]

Cleanup code posted in File Exchange

[Darcy Cordell]

@Bruno Luong Thanks for this answer. It's a great solution, but I do find that it is slower. For my application, I am running this random permutation millions of times, so the speed does matter. I did a comparison with N = 10^7 and found that randpermfull completes in about 1.1 s while your function randder takes about 20 s. For reference, randperm takes 0.60 s, so even randpermfull is going to almost double my runtime :(

[Bruno Luong]

You should download my FEX where I make it quite faster than the code I posted above.

Here is the result of three codes run on my Windows laptop; MATLAB R2020b Update 5:

N = 1e7;
ntest = 16;

clear p;
tic;
for i=1:ntest
    p=randperm(N);
end
t_randperm = toc/ntest % 0.4993 seconds.

clear p;
tic;
for i=1:ntest
    p=randpermfull(N);
end
t_randpermfull = toc/ntest % 1.6284 seconds. The theory said it takes 2.71828*t_randperm in average

clear p;
tic;
for i=1:ntest
    p=randder(N);
end
t_randder = toc/ntest % 2.6069, slowest but very predictable

% Shuffle with MEX implementation is the fatest, but the speed varies as with randpermfull
clear p;
tic;
for i=1:ntest
    p = Shuffle(N, 'derange');
end
t_rShuffle = toc/ntest % Elapsed time is 0.7444

[Bruno Luong]

To illustrate the stability/instability of non-rejection vs rejection methods, here is the script and the results of run on 32 random generation requested:

clear
N = 1e7;
ntest = 32;

dfuns = {@(N) randpermfull(N), ...
         @(N) randder(N), ...
         @(N) Shuffle(N, 'derange')};

     
for k=1:length(dfuns)
    t = zeros(1,ntest);
    df = dfuns{k};
    for i=1:ntest
        tic
        p = df(N);
        t(i) = toc;
    end
    fprintf('%s \n', func2str(df));
    fprintf('\tmeant(t)         = %f\n', mean(t));
    fprintf('\tmax(t)           = %f\n', max(t));
    fprintf('\tmin(t)           = %f\n', min(t));
    fprintf('\t(max-min)/mean   = %f\n', (max(t)-min(t))/mean(t));
end

The results is

>> benchderangement
@(N)randpermfull(N) 
	meant(t)         = 1.666878
	max(t)           = 5.088765
	min(t)           = 0.519352
	(max-min)/mean   = 2.741300
@(N)randder(N) 
	meant(t)         = 2.642569
	max(t)           = 2.747349
	min(t)           = 2.598667
	(max-min)/mean   = 0.056264
@(N)Shuffle(N,'derange') 
	meant(t)         = 0.719942
	max(t)           = 1.726542
	min(t)           = 0.280354
	(max-min)/mean   = 2.008756

Another run

    >> benchderangement
@(N)randpermfull(N) 
	meant(t)         = 1.518073
	max(t)           = 5.090753
	min(t)           = 0.507698
	(max-min)/mean   = 3.018995
@(N)randder(N) 
	meant(t)         = 2.613948
	max(t)           = 2.688238
	min(t)           = 2.404195
	(max-min)/mean   = 0.108665
@(N)Shuffle(N,'derange') 
	meant(t)         = 1.055848
	max(t)           = 4.574268
	min(t)           = 0.278108
	(max-min)/mean   = 4.068919

RANDER is only method that the run time is consistent and predictible.

With rejection methods there is a chance - perhaps tiny - where the run time just blew up.

[Paul]

What is the Shuffle() function? Is it a FEX?

[Stephen23]

@Paul:

https://www.mathworks.com/matlabcentral/fileexchange/27076-shuffle

[Bruno Luong]

I run de tic/toc up to 256 runs for each method and plot the histogram. Here is the results (one can see a hint of the power law (1-1/e)^p for runtime for rejection methods):

>> benchderangement
@(N)randpermfull(N) 
	mean(t)          = 1.432840
	max(t)           = 7.470663
	min(t)           = 0.463656
	(max-min)/mean   = 4.890291
@(N)randder(N) 
	mean(t)          = 1.744887
	max(t)           = 2.063773
	min(t)           = 1.632292
	(max-min)/mean   = 0.247283
@(N)Shuffle(N,'derange') 
	mean(t)          = 0.592406
	max(t)           = 2.253336
	min(t)           = 0.229715
	(max-min)/mean   = 3.415937

Answer 3

[Jeff Miller]

I don't think randperm can do that by itself, but I think this would work for an even number of items in the original vector:

orig = 1:10;        % an example for the original vector of items
nvec = numel(orig);
halfn = nvec/2;
perm1 = randperm(nvec);
final = zeros(size(orig));
for ivec=1:halfn
    swap1pos = perm1(ivec);
    swap2pos = perm1(ivec+halfn);
    final(swap1pos) = orig(swap2pos);
    final(swap2pos) = orig(swap1pos);
end
sum(final==orig)

If you have an odd number of items in the original vector, handle one item specially and use this method with the remaining ones.

Would not be surprised if someone has a better solution, though.

Answer 4

[Image Analyst]

Just keep looping until there are no matches, like this:

n = 5;
originalVector = 1 : n;
maxIterations = 10000;
loopCounter = 1;
while loopCounter < maxIterations
	newVector = randperm(n);
	matches = newVector == originalVector;
	if ~any(matches)
		break; % Break out of loop if there are no matches.
	end
	loopCounter = loopCounter + 1;
end
% Print out newVector to command window.
fprintf('Found answer after %d iterations.\n', loopCounter);
newVector

[Steven Lord]

Setting maxIterations to 10000 is probably a bit overkill given the probability of a random permutation of a set being a derangement.

[Image Analyst]

Yeah, usually it exited the loop with less than 4 iterations, so it doesn't really matter how big it is - it should never get that high. But I always use a failsafe with while loops because often users have errors where they'll get an infinite loop and this will make sure that doesn't happen. If it never did break out of the loop, then the failsafe should kick them out of the loop in a fraction of a second even with ten thousand max iterations. Then, to be super robust, I usually tell the user after the loop if the loopCounter == maxIterations so they know that they did not exit the loop normally.

Yeah, it's a brute force method but it's very fast. Not sure if the derangement code in the File Exchange is more elgant - I didn't look at it.

[Bruno Luong]

"Yeah, it's a brute force method but it's very fast. Not sure if the derangement code in the File Exchange is more elgant - I didn't look at it."

So you claim a brute force method is the most elegant?

Jan's code is amazing piece of code, as usual with him.

[Image Analyst]

I did not claim my code was elegant at all, much less the "most elegant".

And yes, Jan does write very amazing code. And you write very clever code also!

[Darcy Cordell]

I think this is often the fastest solution for large N. But for small N, I think randpermfull is faster.

Answer 5

[David Goodmanson]

Hi Darcy,

the methods I have seen here seem to involve trying randperm and rejecting the result if an element remains in the same location. Here is a method that uses the cycle structure of the permutation and does not allow any 1-cycles (element stays where it is). Randperm is called once. The code uses the fact that if you have n elements, and do a chain of randomly chosen elements starting with a given element, the odds that you obtain a k-cycle is 1/n for every k.

I don't know how the speed is compared to the rejection method, but this code is not slow, taking half a second for 10 million elements.

n = 100
tic
q = randperm(n);
p = zeros(1,n);
nrem = n;           % number of remaining elements
cycstruct = [];     % cycle structure (just the lengths)
while nrem > 0
  if nrem == 2;
    cyc = 2;
  else
    cyc = randi(nrem-2)+1;  % cycle length
    if cyc == nrem-1;       % unallowed cycle length
       cyc = nrem;
    end
  end
  cycstruct = [cycstruct cyc];
  nrem = nrem-cyc;
  ind = q(1:cyc);
  q(1:cyc) = [];
  p(ind) = circshift(ind,-1);
end  
toc
cycstruct  
[1:n; p];                % p is the result
any(p==1:n)              % check if any element stays at home
any(diff(sort(p))~=1)    % any duplicated elements?

[Paul]

@David Goodmanson

Can you comment on the distribution of the results of this procedure (even though the OP didn't speciify a requirement on such)? I ran the code many times with n = 4, for which there are 9 possible outcomes, but the results were not uniformly distributed among those 9 possibilities.

[David Goodmanson]

Hi Paul,

I agree with what you are saying. I had thought that the probability of getting a 4-cycle was the same as getiing a 2-cycle (and hence a pair of 2-cycles) but actually there are six cases of one 4-cycle and three cases of a pair of 2-cycles. So I will go look at that, but meanwhile Bruno has been on the job.