143 views (last 30 days)

Show older comments

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).

Stephen Cobeldick
on 27 Feb 2021

This type of permutation is called a derangement:

You can start by searching FEX:

and reading the descriptions of those submissions.

Bruno Luong
on 27 Feb 2021

Edited: Bruno Luong
on 27 Feb 2021

This will do what you ask for

Even better but need MEX

Bruno Luong
on 28 Feb 2021

Edited: Bruno Luong
on 28 Feb 2021

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
on 3 Mar 2021

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

Jeff Miller
on 26 Feb 2021

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.

Image Analyst
on 27 Feb 2021

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

David Goodmanson
on 28 Feb 2021

Edited: David Goodmanson
on 28 Feb 2021

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?

David Goodmanson
on 28 Feb 2021

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.

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!