Thread Subject: Derangement: efficient full permutation

Dear readers,

An "derangement" is a permutation of the numbers from 1 to N such that no p(i)==i. If you use it to shuffle a vector, no element remains on the same position.
After some discussions about an efficient algorithm, e.g. in the comment section of the former version of RANDPERMFULL in the FEX (Jos, new version: #30189), I had the following idea:
Create the permutation by Knuth's well tested shuffle algorithm, and reject single elements, which would collide with the derangement condition.

function p = Derange(n)
p = 1:n;
for k = 1:n
   r = ceil(rand*n);
   while p(k) == r || p(r) == k % Reject non-derangement
      r = ceil(rand * n);
   end
   t = p(k);
   p(k) = p(r);
   p(r) = t;
end

For large N (> 10000) the most permutations are derangements and the above can be separated for efficiency:
1. Standard Knuth shuffle:
p = 1:n;
for k = 2:n
  r = ceil(rand * k);
  p(k) = p(r);
  p(r) = k;
end

2. Care for p(i)==i elements only:
for k = 1:n
  if p(k) == k
    ... and then the above Derangement method for this k here
  end
end

This is at least 5 times faster for K=1e5 than the SORT(RAND(1, K)) approach in RANDPERMFULL with a following complete rejection. If I use the C-Mex Shuffle(FEX: #27076) to get the Knuth-Shuffle, this is a further factor of 3.5 times faster. So I assume, the full C implementation would be efficient.

BUT: All who had participated in the former discussion (Jos, Derek, me, ...) had suggested algorithm with bugs. Therefore I hope, that somebody can confirm, that this method creates all derangements without a bias - or point me to the errors, before I publish the mexed version in the FEX.

Thanks, Jan

>
> BUT: All who had participated in the former discussion (Jos, Derek, me, ...) had suggested algorithm with bugs. Therefore I hope, that somebody can confirm, that this method creates all derangements without a bias - or point me to the errors, before I publish the mexed version in the FEX.

Your algorithm surely has BUG, e.g., n = 1 it generate a derangement 1, which is impossible.

Bruno

Dear Bruno,

> Your algorithm surely has BUG, e.g., n = 1 it generate a derangement 1, which is impossible.

Yup. I've posted a reduced version only. In the full version I reject N=1 and reply [2,1] for N=2.
Do you see limitations for N>=3 --- except for the necessary warnings for e.g. N > 2^48 or complex N?

Kind regards, Jan

"Jan Simon" wrote in message <ii91fn$5eo$1@fred.mathworks.com>...
> Dear Bruno,
>
> > Your algorithm surely has BUG, e.g., n = 1 it generate a derangement 1, which is impossible.
>
> Yup. I've posted a reduced version only. In the full version I reject N=1 and reply [2,1] for N=2.
> Do you see limitations for N>=3 --- except for the necessary warnings for e.g. N > 2^48 or complex N?

If you have to make a special case for small N, that's a strong evidence of biased in the algorithm. Otherwise you don't have to.

Bruno

"Jan Simon" wrote in message <ii91fn$5eo$1@fred.mathworks.com>...
> Dear Bruno,
>
> > Your algorithm surely has BUG, e.g., n = 1 it generate a derangement 1, which is impossible.
>
> Yup. I've posted a reduced version only. In the full version I reject N=1 and reply [2,1] for N=2.
> Do you see limitations for N>=3 --- except for the necessary warnings for e.g. N > 2^48 or complex N?
>
> Kind regards, Jan

Derange(3) sometime returns [1 2 3].

To me the algorithm is simple flawed.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ii920m$s53$1@fred.mathworks.com>...
> "Jan Simon" wrote in message <ii91fn$5eo$1@fred.mathworks.com>...
> > Dear Bruno,
> >
> > > Your algorithm surely has BUG, e.g., n = 1 it generate a derangement 1, which is impossible.
> >
> > Yup. I've posted a reduced version only. In the full version I reject N=1 and reply [2,1] for N=2.
> > Do you see limitations for N>=3 --- except for the necessary warnings for e.g. N > 2^48 or complex N?
> >
> > Kind regards, Jan
>
> Derange(3) sometime returns [1 2 3].
>
> To me the algorithm is simple flawed.
>
> Bruno

Hmm, I cannot replicate the bug observed by Bruno. Here "derange" produces derangements only. However it is biased:

% bootstrapping
clear R
N = 4 ;
for k=1000:-1:1,
   R(k,:) = derange(N) ;
end ;

%
OnlyDerangements = all(~any(bsxfun(@eq,R,1:N)))
[uR,~,ix] = unique(R,'rows') ;
n = histc(ix,1:size(uR,1)) ;
AllDerangements = (size(uR,1) == floor((prod(1:N)/exp(1)) + 0.5) )
[n(:) uR] % indeed derangements, but with different frequencies !!

... getting derangements is not so simple ...

Jos

"Jos (10584)" wrote in message <ii98b6$sb7$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message
>
> Hmm, I cannot replicate the bug observed by Bruno. Here "derange" produces derangements only. However it is biased:
>

Which Derange implementation? Jan spoke verbally his way of Derange [split in two parts etc] and I certainly missunderstood. Can someone (Jan?) post the exact code? Thanks.

Anyway I find strange the approach here. Knuth's method is rigorously backup by a proof of non-biased and why it's a certainly clever method. If Jan idea of alter Knuth's for derangement problem, then the method also needs to be backup up by a proof (by adapting Knuth's proof). May be the proof is hard, but at least Jan should have a strong feeling about the non-biased of his method. He should elaborate his though.

Bruno

Dear Bruno, dear Jos,

> If you have to make a special case for small N, that's a strong evidence of biased in the algorithm. Otherwise you don't have to.

I don't have to. I exclude the N=2 case only, because the rejection method would reject 50% of the solutions, which is a far to inefficient method to calculte [2, 1] dynamically.

> Derange(3) sometime returns [1 2 3].
> To me the algorithm is simple flawed.

Thanks! Exactly what I searched for.
But I'm surprised. I tried it 1e7 times and do not get results different from [2,3,1] and [3,1,2]. So either your RAND can reply more values than mine, or something other went wrong.

> Jos wrote:
> % indeed derangements, but with different frequencies !!
>... getting derangements is not so simple ...

Thanks! I was convinced it is a good idea to discuss this before publishing it in the FEX. I'm starting further investigations.

Jan

Dear Bruno,

> Which Derange implementation?

The one from the first post of this thread:
  function p = Derange(n)
  p = 1:n;
  for k = 1:n
     r = ceil(rand * n); % Knuth: ceil(rand * k)
     while p(k) == r || p(r) == k
        r = ceil(rand * n);
     end
     t = p(k);
     p(k) = p(r);
     p(r) = t;
  end

The idea was based on the "full rejected" approach:
  p = KnuthShuffle(n);
  while any(p == 1:n)
    p = KnuthShuffle(n);
  end
Here KnuthShuffle is either RANDPERM(n), Shuffle(n, 'index') or SORT(RAND(1, n)), see also Jos' RANDPERMFULL and Derek's DERrej in RPGLab.
My naive approach tried to reuse the elements with p~=1:n and swap just the rare elements colliding with the restriction. But Jos found out, that the index vectors, which are not created by Sattolo's modification of Knuth's shuffle, *are* created at least, but with the half probability than the others.

Afterwards I split the above method into two parts: 1. Standard Knuth shuffle, 2. process the colliding elements in an extra loop.

> Anyway I find strange the approach here. Knuth's method is rigorously backup by a proof of non-biased and why it's a certainly clever method. If Jan idea of alter Knuth's for derangement problem, then the method also needs to be backup up by a proof (by adapting Knuth's proof). May be the proof is hard, but at least Jan should have a strong feeling about the non-biased of his method. He should elaborate his though.

Sorry, Bruno. You spent a lot of time to offer improvements and corrections of my posts in the last days (and years). I hope you do not have the impression, that I want you to do my homework. Of course it is my turn to elaborate the method and I could have been more precise in the formulation, e.g. by asking: "Has anybody investigated this naive approach already and can give an advice if further work is a loss of time?"
However, in opposite to 99.9% of you other posts, the message that my algorithm is wrong because it fails for 1, was not helpful.

We had some failing Derangement algorithms here in CSSM, in the FEX and their comments and in the email contact between some CSSM-users. And now I've added my piece of junk. But I do not think, that it is a too strange approach to show an idea and ask for help. I have the impression, this is very common in CSSM...

Kind regards, Jan

> I have the impression, this is very common in CSSM...

Ah Jan, don't get me wrong, my intention is *not* to complain about someone does someone else homework. My complain is more about the scientific approach of designing an algorithm.

Let's assume Jos's bootstrapping does not show an obvious biased, then I believe one still have to show the algorithm works, otherwise to me it is remains a unreliable algorithm.

As we are all notices, generating permutation with constraints (such as derangement) is particularly tricky, since it related to one of the most difficult mathematics domain: combinatorics. So a proof of yes or no would be very hard. Still if it has to be done, then it must.

Not that because it's a common thing is CSSM that we need to follow it.

You are one of the most rigorous CSSM participant. Almost anything you wrote is sure and exact, and I appreciate a lot this quality. So keep it that way.

Enough about the derangement for today. Let's talk about something else more fun.

Bruno

Dear Bruno,

> Let's assume Jos's bootstrapping does not show an obvious biased, then I believe one still have to show the algorithm works, otherwise to me it is remains a unreliable algorithm.

I'm an experimental physicist. If the cheap bootstrapping reveals a bias, I do not start to think about an expensive proof. I know the fundamental role of prooves, but I keep in mind their limitations, when a prooved algorithm uses e.g. the cheap RNG of Matlab6, real hardware like an Intel processor with a long list of open bugs, and I've heared that even Matlab should have bugs.
Of course you cannot get a reliable program using unreliable algorithms. And you cannot get a reliable program just by implementing reliable algorithms.

And therefore programming combinatorics is fun for *me*.
 
> You are one of the most rigorous CSSM participant. Almost anything you wrote is sure and exact, and I appreciate a lot this quality. So keep it that way.

Am I rigorous?? After a long hesitating, I cannot decide if this is humorous or a kind of invitation. I'm not used that anybody suggests me how to behave. Anyway, I take it as friendly suggestion with a certain level of humor. I'm going to be very careful if I answer questions or publish functions in the FEX. And I feel free to ask what ever I'm interested in, as long as it concerns Matlab. ;-)

Kind regards, Jan

On Feb 1, 10:08 pm, "Jan Simon" <matlab.THIS_Y...@nMINUSsimon.de>
wrote:
> Dear Bruno,
>
> > Let's assume Jos's bootstrapping does not show an obvious biased, then I believe one still have to show the algorithm works, otherwise to me it is remains a unreliable algorithm.
>
> I'm an experimental physicist. If the cheap bootstrapping reveals a bias, I do not start to think about an expensive proof. I know the fundamental role of prooves, but I keep in mind their limitations, when a prooved algorithm uses e.g. the cheap RNG of Matlab6, real hardware like an Intel processor with a long list of open bugs, and I've heared that even Matlab should have bugs.
> Of course you cannot get a reliable program using unreliable algorithms. And you cannot get a reliable program just by implementing reliable algorithms.
>
> Kind regards,Jan

---------------------------

 Dear All,

 The Knuth Shuffle is not Knuth's but Richard Durstenfeld's
implementation of the Fisher-Yates method
 (see Algorithm 235, Communications of the ACM, Jan. 1964).
 It is provably correct and has complexity O(n). Hence it is optimal
up to a multiplicative constant.

 A few weeks ago I tried this rejection variation of Durstenfeld's
algorithm in the hope of generating
 derangements:

p = GRDrejr(n);
p = I(n); % Start with Identity permutation
for k = 2:n
    r = k;
    while r==k
        r = ceil(rand*k);
    end;
    t = p(k);
    p(k) = p(r);
    p(r) = t;
end; return; % GRDrejr

I tested this variation and found (experimentally) that it produces
derangements but with the same flaw that RANDPERMFULL had: all the
derangements are cyclic.

I should have been more careful and used part of Jan's test: while
r==k || p(r)==k

but this modification did not fix the flaw in GRDrejr

I then tested Jan's GRDsimon for various n=5,6,7,8,9,10 and found it
is flawed. For example:

>> target = GRDrej(9);[p,k,t] = FindGRD(target,2); % 2 limits the search to 2*n! generations.
Starting search for [2 3 7 6 9 5 4 1 8]
----------------------------------------------------------
 Target FOUND after k = 21492 iterations. Time = 2.964 secs.
 Fraction of Space searched (k/n!) = 0.059226 Rate = 7251 permutations
per second

 Target and last p generated:
     1 2 3 4 5 6 7 8 9
     2 3 7 6 9 5 4 1 8
     2 3 7 6 9 5 4 1 8

 GRDsimon generated number of permutations = 21492
 Generated number of derangements = 14415 Frac. of
total = 0.67071
 Expected number of derangements = 21492
 Generated number of cyclic permutations = 7398 Frac. of
total = 0.34422
 Expected number of cyclic permutations = 6492
 
-----------------------------------------------------------------------------------------------------------------

This shows that not all permutations generated by GRDsimon are
derangements. In this case GRDsimon was lucky and hit the target after
just 21492 generations. On this next test it was not so lucky.

>> target = GRDrej(10); [p,k,t] = FindGRD(target,2);

Starting search for [8 3 4 1 2 9 5 10 6 7]
 
----------------------------------------------------------------------------------------------------------
 Target NOT FOUND after k = 7257600 iterations. Time = 1051.1338 secs.
 Fraction of Space searched (k/n!) = 2 Rate = 6905 permutations per
second

 Target and last p generated:
     1 2 3 4 5 6 7 8 9 10
     8 3 4 1 2 9 5 10 6 7
     5 1 6 8 2 7 9 3 10 4

 GRDsimon generated number of permutations = 7257600
 Generated number of derangements = 4823421 Frac. of
total = 0.6646
 Expected number of derangements = 7257600
 Generated number of cyclic permutations = 2236580 Frac. of
total = 0.30817
 Expected number of cyclic permutations = 1972821
 
---------------------------------------------------------------------------------------------------------------------------

Repeated tests gave similar results.

Doing some "one-shot" tests for large n also showed up the same flaw:

>> n = 10^6;t=tic; p=GRDsimon(n); t=toc(t); [IsDer(p) IsCyc(p)]

ans = 0 0

GRDsimon seems to generate derangements about 67% of the time.

Note that GRDrej, GRDmpp, and GRDsimon are all rejection algorithms
and have random running times.

Durstenfeld's Algorithm (my GRPfys) and Sattola's cyclic algorithm (my
GRCsat) have constant running times for a fixed n. I wonder why it is
so hard to find such an algorithm for derangements. Maybe some
computational complexity expert can answer this question.

Regards,

Derek O'Connor

Dear Derek,

Thanks for your answer.

> I then tested Jan's GRDsimon for various n=5,6,7,8,9,10 and found it
> is flawed. [...]
> This shows that not all permutations generated by GRDsimon are
> derangements.

No doubt, my approach is flawed. But I cannot reproduce a non-derangement permutation. Although this is surprising to me, it is not worth to discuss if my approach is "more correct" than you assume, as long as it is less than correct.

Kind regards, Jan

On Feb 5, 11:16 pm, "Jan Simon" <matlab.THIS_Y...@nMINUSsimon.de>
wrote:
> Dear Derek,
>
> Thanks for your answer.
>
> > I then testedJan'sGRDsimon for various n=5,6,7,8,9,10 and found it
> > is flawed. [...]
> > This shows that not all permutations generated by GRDsimon are
> > derangements.
>
> No doubt, my approach is flawed. But I cannot reproduce a non-derangement permutation. Although this is surprising to me, it is not worth to discuss if my approach is "more correct" than you assume, as long as it is less than correct.
>
> Kind regards,Jan


Jan, I'm re-posting this because after 6 hours it has tot appeared.

POSTED to CSSM Sunday, February 6, 2011 10:26

Dear Jan,

Mea culpa, mea culpa, mea maxima culpa!

Everything I said in my previous post is complete rubbish, but let
it stay there as a lesson.

In error I was using for k = 2:n. When I changed this to your for k =
1:n
then p = GRDsim(n) worked perfectly in all the tests I have done on it
so far.

Here are some typical results

 
-------------------------------------------------------------------------------
>> n=10; target = GRDrej(n); [p,k,t] = FindGRDrejr(target, 3);

Starting search for [6 3 1 9 4 8 5 10 2 7]

 Target FOUND after k = 1040194 iterations. Time = 23.0014 secs.
 Fraction of Space searched (k/n!) = 0.28665 Rate = 45224 permutations
per second

 Target and last p generated:
     1 2 3 4 5 6 7 8 9 10
     6 3 1 9 4 8 5 10 2 7
     6 3 1 9 4 8 5 10 2 7

 GRDsim generated number of permutations = 1040194
 Generated number of derangements = 1040194 Frac. of total
= 1
 Expected number of derangements = 1040194
 Generated number of cyclic permutations = 234379 Frac. of total
= 0.22532
 Expected number of cyclic permutations = 282755
 Number of random numbers generated = 11888239 Frac. of min.
= 1.1429
 
--------------------------------------------------------------------------------

As you can see, all your permutations are derangements and 22.5% of
them are cyclic which is about right. This suggests that your
method is sampling from the complete space of derangements and is
not biased.

A very significant feature of GRDsim is that it uses only 1.1429
times the random numbers that GRPfys uses (based on my limited
tests). This means that it is a very efficient Las Vegas rejection
algorithm. Remember that GRDrej uses 2.718 and GRDmpp uses 2 times
the random numbers that GRPfys uses, on average.

Here are timing tests

>>[times,ntimes] = TestDerangeTimes(0,[10^4 10^5 10^6],10);
 ------------------------------------------------------------------
  Four Derangement Generators. No. Samples = 10
 ------------------------------------------------------------------
            n GRDrej GRDsim GRDmpp GRDmex
  Times (secs) ----------------------------------------------------
        10000 0.0018628 0.0016665 0.0021822 0.0020541
       100000 0.016393 0.015044 0.020054 0.0072654
      1000000 0.48478 0.25515 0.29849 0.19957

  Normalized Times -----------------------------------------------
        10000 1.1177 1 1.3094 1.2326
       100000 2.2563 2.0707 2.7601 1
      1000000 2.4291 1.2785 1.4956 1
------------------------------------------------------------------

You can see that your GRDsim is fast and even manages to beat
GRDmex sometimes.

If GRDsim passes further tests then I believe you have invented a
very significant combinatorial algorithm. I wish I had the
mathematical ability to prove its correctness and analyze its
random running time. I suggest you contact Martinez et al.

Regards,

Derek O'Connor

Derek:
>
> As you can see, all your permutations are derangements and 22.5% of
> them are cyclic which is about right. This suggests that your
> method is sampling from the complete space of derangements and is
> not biased.
>

Humm I still don't quite manage to be convinced by the conclusion. Possibly because the code I use is wrong?

  function p = Derange(n) % Jan Simon's derangement method
  p = 1:n;
  for k = 1:n
     r = ceil(rand * n); % Knuth: ceil(rand * k)
     while p(k) == r || p(r) == k
        r = ceil(rand * n);
     end
     t = p(k);
     p(k) = p(r);
     p(r) = t;
  end

%% Script to test JS's derangement

m = 100000;
n = 3;
D = zeros(m,n);
for k=1:m
    D(k,:) = Derange(n);
end
[allD, ~, J] = unique(D,'rows');
P = accumarray(J(:), 1) / m;
fprintf('JS'' derangement with n = %d\n', n);
fprintf('number of experiences m = %d\n', m);
for k=1:size(allD,1)
    fprintf('P( %s )= %g\n', mat2str(allD(k,:)), P(k)) ;
end

% Here is the result I get:

JS' derangement with n = 3
number of experiences m = 100000
P( [2 3 1] )= 0.62399
P( [3 1 2] )= 0.37601

That number looks like a bias to me.

Bruno

Dear Derek, lovely newgroup!

This thread drives me crazy - and I am the originator.

1. I've posted a bold approach after a faulty test let me estimate the method to produce valid data.
2. Bruno finds a not existing bug in the code. (This is worth to mention, because Bruno fails very rarely, as we all know)
3. Jos found the obvious problems of my algorithm.
4. Bruno posts the not useful criticism, that my code fails for a derangement of length 1, which is an absolute irrelavent exception. Finally I have to face, that he calls me "rigorous". (This is worth to mention, because I have read just one other not useful post from him, and this concerned the trolling of Marco)
5. After some misunderstandings, Derek found out, that I'm a super hero.

What's the next step?
I should prove that there are no solutions for:
  a^n + b^n = c^n with n = 0.

Derek, let's dicuss this per email before we reveal too many awkward details.
Thanks for your interest, Jan

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <iin0v4$ic2$1@fred.mathworks.com>...
> Derek:
> >
> > As you can see, all your permutations are derangements and 22.5% of
> > them are cyclic which is about right. This suggests that your
> > method is sampling from the complete space of derangements and is
> > not biased.
> >
>
> Humm I still don't quite manage to be convinced by the conclusion. Possibly because the code I use is wrong?
>
-----snip ------------------
> % Here is the result I get:
>
> JS' derangement with n = 3
> number of experiences m = 100000
> P( [2 3 1] )= 0.62399
> P( [3 1 2] )= 0.37601
>
> That number looks like a bias to me.
>
> Bruno

-----------------------------------------------

Dear Bruno,

Yes, for n = 3 it is biased. But...

I re-wrote your script as a function:

function P = TestBrunoGRDsim(n,ns);
D = zeros(ns,n);
for k=1:ns
    D(k,:) = Derange(n);
end;
[allD, ~, J] = unique(D,'rows');
P = accumarray(J(:), 1) / ns;
return;

Then I ran it, with ns = 10^4,

P = TestBrunoGRDsim(n,ns); semilogy(P,'.');

Go through n = 3,4,5,6,7,8,9,10,11,12,20 and watch what happens.

There is very odd behavior at n = 7, 8, 9, but for n = 10 you get a
horizontal line at 1/ns, with a handfull of dots above it.

For n = 11, 12, 20 ... you get a horizontal line at 1/ns = 10^(-4) only,
which means that GRDsim is generating derangements uniformly.

I assume (hope?) this uniformity continues for much higher values
of n, where such testing is impossible.


Regards,

Derek.

"Derek O'Connor" wrote in message <iin5l5$h5n$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message
>
> I assume (hope?) this uniformity continues for much higher values
> of n, where such testing is impossible.
>

I think that observation can be explained very well: for large n, derangement converges to permutation, and Jan Simon's converges to Knuth's. Though the small bias will be there.

Bruno

Dear Bruno,

> Humm I still don't quite manage to be convinced by the conclusion. Possibly because the code I use is wrong?

You are correct: The code you use is wrong. *And* it is the code I've posted.
The 2nd method I've described partially in my first post is not biased for N=2 and N=3, but for N>=4, and therefore it just of limited use...

Thanks for your comment, which is completely helpful again.

Jan


 

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <iin77n$rje$1@fred.mathworks.com>...
> "Derek O'Connor" wrote in message <iin5l5$h5n$1@fred.mathworks.com>...
> > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message
> >
> > I assume (hope?) this uniformity continues for much higher values
> > of n, where such testing is impossible.
> >
>
> I think that observation can be explained very well: for large n, derangement converges to permutation, and Jan Simon's converges to Knuth's. Though the small bias will be there.
>
> Bruno

Bruno,

I have tested GRDsim using your very useful TestBrunoGRDsim(n,ns) for n = 3.:20, and ns up to 10^7. In all tests GRDsim gave derangements, so it does not 'converge' to Durstenfeld's (Knuth) algorithm, which generates all possible permutations. Also, at n =
20, ns = 10^7, your test showed that the derangements were generated uniformly.

I was still intrigued by the strange behavior of GRDsim for n = 3:10 so I tested my GRDrej, which is 'provably correct', being based on GRPfys. I got exactly the same strange behavior!

Then I tested GRDmpp (Martinez et al.) and got the same behavior, converging to uniformity at about n = 12,13.

This implied that the strange frequency behavior was coming from GRPfys, so I tested it with TestBrunoGRDsim(n,ns), n = 3:12, ns up to 10^7. Of course, only 37% of the permutations generated are derangements. But again, I got exactly the same strange behavior.

So it seems that the source of the strange frequency behavior for small values of n, is due to the Durstenfeld-Fisher-Yates algorithm. I have no idea why this is so. These small, combinatorial algorithms are still a mystery to me.

Bruno, I think we may have been arguing at cross purposes: I don't care about the behavior of GRPfys(n), GRDrej(n), and GRDsim(n), for n < 20, but you do.

Regards,

Derek.

"Derek O'Connor" wrote in message <iio9b0$m0i$1@fred.mathworks.com>...
>
> So it seems that the source of the strange frequency behavior for small values of n, is due to the Durstenfeld-Fisher-Yates algorithm. I have no idea why this is so. These small, combinatorial algorithms are still a mystery to me.

That's odd. I think it merits to find out what is going on.

>
> Bruno, I think we may have been arguing at cross purposes: I don't care about the behavior of GRPfys(n), GRDrej(n), and GRDsim(n), for n < 20, but you do.
>

Yes I do but not for any specific application in mind. I just do because it shows the correctness of the algorithm, and the bias as we have seen will be more clearly seen for smaller n.

Bruno

Dear Derek,

> > Bruno wrote:
> > I think that observation can be explained very well: for large n, derangement converges to permutation, and Jan Simon's converges to Knuth's.

I agree: The difference between the standard shuffle and my derange algorithm is the on-the-fly rejection of a swap, if it produces a fix point. As Bruno says, for large N the chance for getting a fix point decreases and the rejection is unlikely to used ever. Or in other words: A long permutation is most likely a derangement and the special treatment of the rare exceptions will cause a decreasing bias only - but it will not vanish.

> I was still intrigued by the strange behavior of GRDsim for n = 3:10 so I tested my GRDrej, which is 'provably correct', being based on GRPfys. I got exactly the same strange behavior!

Two different methods show a strange behaviour? Then the behaviour could be hard coded in the test and not an effect of the algorithms.
 
I'd never enter an airplane, if it was designed by using a Monte-Carlo-method based on my Derangement algorithm.

Jan

"Jan Simon" wrote in message <iios8j$eno$1@fred.mathworks.com>...
> Dear Derek,
>
> > > Bruno wrote:
> > > I think that observation can be explained very well: for large n, derangement converges to permutation, and Jan Simon's converges to Knuth's.
>
> I agree: The difference between the standard shuffle and my derange algorithm is the on-the-fly rejection of a swap, if it produces a fix point. As Bruno says, for large N the chance for getting a fix point decreases and the rejection is unlikely to used ever. Or in other words: A long permutation is most likely a derangement and the special treatment of the rare exceptions will cause a decreasing bias only - but it will not vanish.
>
> > I was still intrigued by the strange behavior of GRDsim for n = 3:10 so I tested my GRDrej, which is 'provably correct', being based on GRPfys. I got exactly the same strange behavior!
>
> Two different methods show a strange behaviour? Then the behaviour could be hard coded in the test and not an effect of the algorithms.
>
> I'd never enter an airplane, if it was designed by using a Monte-Carlo-method based on my Derangement algorithm.
>
> Jan
---------------------------------------

Jan,

As the Americans say "Hang on one cotton-pickin' moment" . You say :

>"As Bruno says, for large N the chance for getting a fix point decreases and the rejection >is unlikely to used ever. Or in other words: A long permutation is most likely a >derangement and the special treatment of the rare exceptions will cause a decreasing >bias only - but it will not vanish."


Theorem: Let D(n) be the number of derangements in the set of n! permutations of length n. Then,

as n --> inf, D(n)/n! --> 1/e, independent of n.

In other words, on average, 1/e = 37% of all permutations generated by GRPfys will be derangements. What do you think the remaining 63% will be?

Yes, 63% of all permutations generated by GRPfys will have one or more fixed points, no matter how large n is.

>>RPG=@(n)GRPfys(n); n=10;nsamp=10^4; TestClassFreq(RPG,n,nsamp);

--------------- Class Frequency Test ----------------
n = 10 Perms Derange Cyclic HasFP
Actual 1.000 0.368 0.100 0.632000
Expect 1.000 0.368 0.100 0.632121
------------------------------------------------------

--------------- Class Frequency Test ----------------
n = 100 Perms Derange Cyclic HasFP
Actual 1.000 0.375 0.010 0.624700
Expect 1.000 0.368 0.010 0.632121
------------------------------------------------------


--------------- Class Frequency Test ----------------
n = 1000 Perms Derange Cyclic HasFP
Actual 1.000 0.371 0.001 0.629100
Expect 1.000 0.368 0.001 0.632121
------------------------------------------------------


I would certainly fly in a plane that uses GRDsim(n), n > 100.

Derek.

Jan,

I looked at (and traced through) your algorithm and found the reason why it is biased.

For N = 3, there are only 5 possible [k,r] "paths" through the code. I print them below with the resulting derangement beside. Each row in the brackets is one [k,r] selection through the FOR loop.

[1 2;2 2;3 2] 231
[1 2;2 3;3 3] 231
[1 3;2 1;3 3] 231
[1 2;2 2;3 1] 312
[1 3;2 3;3 3] 312

Thus you are eliminating future options in an biased way as you choose [k,r] early in the FOR loop. Note that if you draw [1 2] on the first iteration of the FOR loop, then 2 out of three times you will get one derangement instead of the other.

I also traced through for N = 4, and there are 56 [k r] paths possible, (I won't list them here!) but we see a similar pattern. In this case, there are 2 of the derangements which are produced by 9 paths, 5 which are produces by 6 paths and 2 that are produced by 4 paths.

Dear Matt,

> I looked at (and traced through) your algorithm and found the reason why it is biased.

Thanks! It is very helpful to find the cause of the bias, because this allows to estimate if the effects vanish fast enough for large N. But most of all I try to find a method to avoid the bias, of course.

Kind regards, Jan

"Jan Simon" wrote in message <ii8ut8$3fb$1@fred.mathworks.com>...
> Dear readers,
> An "derangement" is a permutation of the numbers from 1 to N such that no p(i)==i. If you use it to shuffle a vector, no element remains on the same position.
- - - - - - - -
  I just noticed this interesting thread the other day and couldn't resist throwing in my two cents. The following is a non-rejection scheme that is calculated to have a uniform distribution over all possible derangement permutations of 1:n, that being the objection raised in this thread about Jan's code.

  Obviously a penalty is being paid for achieving this uniformity. This penalty relative to other methods would be less if it were rewritten so that many random derangements could be requested in a single call, since the probability vector below would only have to be computed once.

  It should be noted that the algorithm is O(n) so that for large n it might be competitive with rejection methods using an O(n*log(n)) sort operation.

Roger Stafford
=============================================
function A = randderange(n)

% Column vector A is generated as a random derangement
% permutation of (1:n).' which is uniformly distributed
% among all possible such derangements. The algorithm
% does not involve rejection - all generated vectors are
% valid derangements. Roger Stafford - 2/14/11

if n < 2, error('n must be at least 2'), end

% Compute probability decision vector, P
P = zeros(n-1,1);
n0 = min(n-1,17);
s = 1;
d2 = 0;
for k = 1:n0
 d1 = d2;
 d2 = k*s;
 % At this point d1 = !k and d2 = !(k+1)
 s = d1 + d2;
P(k) = d1/s; % = !k/(!k+!(k+1)), approx. 1/(k+2) for large k
end
% Past k = 17, d1 and d2 suffer round-off errors
% so the following approximation is used instead.
P(n0+1:n-1) = 1./(n0+3:n+1);

% Generate random derangement in A
A = zeros(n,1);
B = (1:n-1).';
u = n; % Pointer to beginning of cycle
v = n; % Pointer to current end of cycle
b = true;
for k = n-1:-1:1
 if b % Force cycle to have at least two elements
  t = 0;
  b = false;
 else
  t = P(k); % Use P table for uniform distribution
 end
 if t <= rand % Decide whether to extent cycle further
  w = ceil(k*rand); % Random choice in B
  A(v) = B(w); % Insert it in A
  v = B(w); % Advance end pointer
  B(w) = B(k); % Fill in used space in B
 else % Close current cycle and begin another
  A(v) = u; % Complete cycle
  u = B(k); % Start a new cycle
  v = u;
  b = true; % Cycle must have at least one more element
 end
end
A(v) = u; % Complete the last cycle

return

Roger, just a quick question:

Do you decompose a permutation in random cycles and find a derangement as permutation such that no cycle is length-one?

Bruno

"Roger Stafford" wrote in message <ijcp3v$4i2$1@fred.mathworks.com>...
> "Jan Simon" wrote in message <ii8ut8$3fb$1@fred.mathworks.com>...
> > Dear readers,
> > An "derangement" is a permutation of the numbers from 1 to N such that no p(i)==i. If you use it to shuffle a vector, no element remains on the same position.
> - - - - - - - -
> I just noticed this interesting thread the other day and couldn't resist throwing in my two cents. The following is a non-rejection scheme that is calculated to have a uniform distribution over all possible derangement permutations of 1:n, that being the objection raised in this thread about Jan's code.
>
> Obviously a penalty is being paid for achieving this uniformity. This penalty relative to other methods would be less if it were rewritten so that many random derangements could be requested in a single call, since the probability vector below would only have to be computed once.
>
> It should be noted that the algorithm is O(n) so that for large n it might be competitive with rejection methods using an O(n*log(n)) sort operation.
>
> Roger Stafford
> ----- snip --------------------



I have done some quick tests on your method (which I call GRDsta) and so far so good. Here is a comparison of normalized times

 Normalized Times: Four Derangement Generators. No. Samples = 10 ------------------------------------------------------------------
            n GRDrej GRDsim GRDmpp GRDsta ----------------------------------------------------------------
       1e+005 1.1864 1 1.3028 1.3527
       1e+006 1.7055 1 1.1602 1.4618
       1e+007 1.8913 1 1.1876 1.3999

I have not analyzed it in detail, but my main criticism, so far, is that your method, GRDsta, uses 3 times as much memory as GRDrej and GRDsim: P = zeros(n-1,1); A = zeros(n,1); B = (1:n-1).';

Note that none of the other methods uses sorting, and that the provable expected running times for GRDrej and GRDmpp are 2.8*GRPfys and 2.0*GRPfys, respectively. I have shown experimentally that GRDsim is about 1.15*GRPfys. (See previous reply above).

Do you have a detailed expected running time analysis of GRDsta? I'm trying to save myself some work.

Regards,

Derek O'Connor.

"Roger Stafford" wrote in message <ijcp3v$4i2$1@fred.mathworks.com>...
> "Jan Simon" wrote in message <ii8ut8$3fb$1@fred.mathworks.com>...
> > Dear readers,
> > An "derangement" is a permutation of the numbers from 1 to N such that no p(i)==i. If you use it to shuffle a vector, no element remains on the same position.
> - - - - - - - -
> I just noticed this interesting thread the other day and couldn't resist throwing in my two cents. The following is a non-rejection scheme that is calculated to have a uniform distribution over all possible derangement permutations of 1:n, that being the objection raised in this thread about Jan's code.
>
> Obviously a penalty is being paid for achieving this uniformity. This penalty relative to other methods would be less if it were rewritten so that many random derangements could be requested in a single call, since the probability vector below would only have to be computed once.
>
> It should be noted that the algorithm is O(n) so that for large n it might be competitive with rejection methods using an O(n*log(n)) sort operation.
>
> Roger Stafford
> =============================================
> function A = randderange(n)

Hi Roger,
Isn't this just the Martinez algorithm? Below is version I deleted from the FEX since the rejection scheme proved much more efficient.
However, I think a few improvements can still be made on this taking advantage of the power of matlab
~ Jos

----

function R = derange_martinez(Nelements)

% The implementation is based on the algorithm described in the paper
% "Generating Random Derangements" by C. Martinez, A. Panholzer & H.
% Prodinger, published in the "Proceedings of the 10 ACM-SIAM Workshop on
% Algorithm Engineering and Experiments (ALENEX)" (Munro et al. (eds), pp
% 234-240; SIAM 2008).
% A copy of this paper can be found here:
% http://www.siam.org/proceedings/analco/2008/anl08_022martinezc.pdf


if Nelements < 2
    % No derangements exist, so return []
    R = [] ;
elseif Nelements == 2,
    % Trivial derangement, keep shape
    R = [2 1] ;
else
    % The main algorithm, in short:
    % An array R is filled with the values 1:N. Working from the
    % right, we change the element R(k) with a random element R(j), where j
    % is in the interval [1 k-1]. So, an element can move to the left
    % multiple times, but it can move to the right of the array only once!
    % This algorithm ensures that we get a proper derangement, and also
    % that each possible derangement can be obtained.
    % However, we need to tweak the chance that an item can move to the
    % left more than once to get the same probability for each possible
    % derangement. This is done by flagging the element at position j so
    % that it will not move again, with a certain probability. This
    % probability (DNratio in the code below) is derived from the number of
    % possible derangements.
    % For details and the formal proofs, see Martinez (2008)
    
    % initialize array of indices
    R = 1:Nelements ; % A in the description above
    % All elements are allowed to swap positions
    ElementCanSwap = true(1,Nelements) ; % the "flag"
    
    % An auxiliary array to retrieve index numbers fast (similar to FIND)
    IDX0 = 1:Nelements ;
    
    % Make a list of random values. This value
    % is compared to probability that be calculated before the loop.
    randomValue = rand(1,Nelements) ;
    
    % The number of derangements D for a set of n elements is given by the
    % formulae: D(n) = floor(n!/e + 0.5)
    % See http://en.wikipedia.org/wiki/Derangement
    % The values increases rapidly, so we only calculate it up to a limit
    DNlimit = min(25,Nelements) ;
    DN = floor((cumprod(1:DNlimit)/exp(1)) + 0.5) ;
    
    % In the present algorithm we only use the ratio between D(n) and
    % D(n-2). For larger N the ratio "D(n) / D(n-2)" is closely approximated
    % by "1 / n(n-1)". Then this is multiplied by (n-1). For large N, this
    % becomes "(n-1) / n(n-1)" which equals "1/n".
    % Therefore we can initialize D(n) directly as ...
    DNratio = 1 ./ IDX0 ;
    % ... and we corerct for smaller N by using the formula:
    x = 3:DNlimit;
    % DNratio(x) = (x-1) .* DN(x-2) ./ DN(x) ;
    DNratio(x) = (x-1) .* DN(x-2) ./ DN(x) ;
    % the first two elements are not used
    
    % The value NtoSwap holds the number of positions that can
    % be modified.
    NtoSwap = Nelements ;
    
    % We start at the end of the array.
    for currentIDX = Nelements:-1:2,
        if NtoSwap > 1
            % Check if the current element may change
            if ElementCanSwap (currentIDX)
                
                % The element at position k will not move again in future
                % iterations of the loop.
                ElementCanSwap(currentIDX) = false ;
                
                % The list of possible positions that can swap elements with
                % the current position (k). Since these positions are all within the
                % range [1, K-1], we always get a derangement!
                PossibleIDX = IDX0(ElementCanSwap) ;
                % select a random position from this set
                j = ceil(numel(PossibleIDX) * rand(1)) ;
                randomIDX = PossibleIDX(j) ;
                
                % swap the elements at position k and the random position
                tmp = R(currentIDX) ;
                R(currentIDX) = R(randomIDX) ;
                R(randomIDX) = tmp ;
                
                % Randomly, the element at position ri will also not move
                % again. The random process ensures that elements will
                % move a random number of times (but at least once, namely
                % just now), before settling on their final position. With
                % this, we change the probability that a particular
                % derangement is obtained. Since DNratio is carefully
                % designed, each derangement becomes equally likely.
                if randomValue(currentIDX) < DNratio(NtoSwap)
                    ElementCanSwap(randomIDX) = false ;
                    NtoSwap = NtoSwap - 1 ;
                end
                
                % The element at position k will not move again
                NtoSwap = NtoSwap - 1 ;
            end
        else
            break ; % Nothing more to change, so break the loop
        end % end of if
    end % end of backward for-loop
    
    % The indices 1:N have been deranged. Now re-position the elements of V
    % accordingly.
end

"Derek O'Connor" wrote in message <ijdr94$dog$1@fred.mathworks.com>...
> I have done some quick tests on your method (which I call GRDsta) and so far so good. Here is a comparison of normalized times
>
> Normalized Times: Four Derangement Generators. No. Samples = 10 ------------------------------------------------------------------
> n GRDrej GRDsim GRDmpp GRDsta ----------------------------------------------------------------
> 1e+005 1.1864 1 1.3028 1.3527
> 1e+006 1.7055 1 1.1602 1.4618
> 1e+007 1.8913 1 1.1876 1.3999
>
> I have not analyzed it in detail, but my main criticism, so far, is that your method, GRDsta, uses 3 times as much memory as GRDrej and GRDsim: P = zeros(n-1,1); A = zeros(n,1); B = (1:n-1).';
>
> Note that none of the other methods uses sorting, and that the provable expected running times for GRDrej and GRDmpp are 2.8*GRPfys and 2.0*GRPfys, respectively. I have shown experimentally that GRDsim is about 1.15*GRPfys. (See previous reply above).
>
> Do you have a detailed expected running time analysis of GRDsta? I'm trying to save myself some work.
>
> Regards,
>
> Derek O'Connor.
- - - - - - - -
  Thanks for running those comparative tests, Derek. It would appear that the version I presented didn't compare very favorably with the other versions you tested. That has inspired me to revise the algorithm so that all but a much smaller part could be vectorized. It is still designed to yield a uniform distribution and on my machine it does run considerably faster but that is with an older matlab version. I would be pleased if you would try it out and see how it compares.

  It does call on 'randperm' which I asserted earlier was O(n*log(n)) because of sorting. However I realize now there are O(n) algorithms for accomplishing the same thing without doing a sort. For example, an efficiently coded mex equivalent to the following would be O(n) and presumably able to outdo a sort(rand(1,n)) at some point as n increases:

A = zeros(1,n); B = 1:n;
for k = n:-1:1, r = ceil(k*rand); A(k) = B(r); B(r) = B(k); end

I am hoping Mathworks has written their current version of 'randperm' so as to be O(n) in some such manner.

  Here is my revised derange code. I will put off explaining its logic until we see if it is of some use.

------
function p = derange2(n) % n >= 2
m = min(n,18);
p = fliplr(floor(cumprod(1:m)/exp(1)+.5));
b = rand(1,n+1) < [1,1./(n+1:-1:m+2),p(2:m)./(p(2:m)+p(1:m-1)),1];
for k = 2:n+1
 b(k) = b(k)&(~b(k-1));
end
f = find(b);
q = 2:n+1;
q(f(2:end)-1) = f(1:end-1);
p = randperm(n);
p(p) = p(q);
-------

Roger Stafford

"Roger Stafford" wrote in message <ijifmf$6uv$1@fred.mathworks.com>...

>
> I am hoping Mathworks has written their current version of 'randperm' so as to be O(n) in some such manner.

Hi Roger,

In the mean time you might use Jan's Shuffle on FEX, which I believe O(n).

I'm wait for Jos's Mex implementation of Martinez method (A real deal IMO) to make a benchmark.

Bruno

"Roger Stafford" wrote in message <ijifmf$6uv$1@fred.mathworks.com>...
> "Derek O'Connor" wrote in message <ijdr94$dog$1@fred.mathworks.com>...

> Thanks for running those comparative tests, Derek. It would appear that the version I presented didn't compare very favorably with the other versions you tested. That has inspired me to revise the algorithm so that all but a much smaller part could be vectorized. It is still designed to yield a uniform distribution and on my machine it does run considerably faster but that is with an older matlab version. I would be pleased if you would try it out and see how it compares.
>
> It does call on 'randperm' which I asserted earlier was O(n*log(n)) because of sorting. However I realize now there are O(n) algorithms for accomplishing the same thing without doing a sort. For example, an efficiently coded mex equivalent to the following would be O(n) and presumably able to outdo a sort(rand(1,n)) at some point as n increases:
>-- snip ----
> I am hoping Mathworks has written their current version of 'randperm' so as to be O(n) in some such manner.
>
> Here is my revised derange code. I will put off explaining its logic until we see if it is of some use.
> ---- snip -----------------------
>
> Roger Stafford

I have run your old and new versions of randderange alongside others which are in RPGLab (along with notes on their implementation and testing).

http://www.mathworks.com/matlabcentral/fileexchange/30101-rpg-lab

For consistency I have renamed your randderange and randderange2 as GRDsta1 and GRDsta2. Here is a brief description of each generator:

GRPfys: Pure Matlab implementation of Durstenfeld's permutation generator which is an O(n) implementation of the Fisher-Yates Shuffle.

GRDrej: Pure Matlab derangement generator and uses GRPfys in a rejection loop.

GRDmex: GRDrej using Jan Simon's mexed version of GRPfys, called GRPmex.

GRDsim: Jan Simon's new rejection derangement generator.

GRDmpp: Martinez, et al. rejection derangement generator.

GRDsta1: R. Stafford's randderange.

GRDsta2: R. Stafford's randderange2, using Matlab's randperm.

GRDstamex: GRDsta2, using GRPmex.

R2008a 64-bit, Windows7 64-bit, 2xQuad Xeon 2.33GHz, 16GB ram.

Seven Derangement Generators. No. Samples = 10. Normalized Times

n GRDrej GRDmex GRDsim GRDmpp GRDsta1 GRDsta2 GRDstamex

10^5 7.39 1 4.12 5.12 5.91 6.95 3.78
10^6 21.44 1 3.33 4.08 5.37 5.94 3.37
10^7 9.58 1 3.36 3.99 5.01 5.21 3.22

Avg 12.81 1 3.60 4.40 5.43 6.04 3.46

Comparing mexed with pure Matlab functions is not really 'fair'. However the same applies to the use of sort in Matlab's randperm. Using sort is a slick but lazy way of generating random permutations, because it is O(nlogn) and uses extra memory. I think Matlab should implement GRPfys, which is O(n) and uses just 1 array, p(1,n);

Also, comparing times is difficult because most(all?) the functions are Las Vegas algorithms, i.e., their running times are random, and so the table above changes each time the test is run. I need to run it with larger number of samples, but that takes time.

Regards,

Derek.

"Derek O'Connor" wrote in message <iio9b0$m0i$1@fred.mathworks.com>...

> I was still intrigued by the strange behavior of GRDsim for n = 3:10 so I tested my GRDrej, which is 'provably correct', being based on GRPfys. I got exactly the same strange behavior!

Sorry for bringing back an old thread, but I check GRDrej (rejection Durstenfeld based) (with small n) and did not notice any strange behavior.

Roger's code also behave correctly, and it's a non-rejection code!

Martinez is a no-no for small n. It proves that Roger's method is not at all equivalent (so that answers Jos's earlier question).

I won't be bother anymore with Martinez's method, which plays in another category.

I have my preferable derangement code now. Thanks Roger.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ijk9se$gcb$1@fred.mathworks.com>...
> "Derek O'Connor" wrote in message <iio9b0$m0i$1@fred.mathworks.com>...
>
> > I was still intrigued by the strange behavior of GRDsim for n = 3:10 so I tested my GRDrej, which is 'provably correct', being based on GRPfys. I got exactly the same strange behavior!
>
> Sorry for bringing back an old thread, but I check GRDrej (rejection Durstenfeld based) (with small n) and did not notice any strange behavior.
>
> Roger's code also behave correctly, and it's a non-rejection code!
>
> Martinez is a no-no for small n. It proves that Roger's method is not at all equivalent (so that answers Jos's earlier question).
>
> I won't be bother anymore with Martinez's method, which plays in another category.
>
> I have my preferable derangement code now. Thanks Roger.
>
> Bruno
---------------------------------

Bruno,

I would not get too excited, at least not yet.

Here is a more accurate set of average times(secs) for n = 10^6 and nsamp = 100;

GRDsim : 0.23744
GRDmpp : 0.28213
GRDsta1: 0.37914
GRDsta2: 0.41287
GRDrej : 0.54378

GRDmex : 0.10005

Also, remember that Roger's method uses at least twice as much memory as GRDrej and GRDsim, which becomes a problem for n > 10^8, with 16GB ram.

Do you still believe that GRDsim (J. Simon's new method) is wrong? I still cannot find a counter-example.


Regards,

Derek.

"Derek O'Connor" wrote in message <ijlldg$181$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ijk9se$gcb$1@fred.mathworks.com>...

>
> Do you still believe that GRDsim (J. Simon's new method) is wrong? I still cannot find a counter-example.

Hi Derek,

I have no reason to change my opinion since the last time. I think Jan's method has flaws.

Roger's method is the only one that gives unbiased results (even for small value of n), has deterministic running time and relatively fast. It needs to be plugged with an O(n) RANDPERM engine, but that's easy.

Bruno

"Derek O'Connor" wrote in message <ijlldg$181$1@fred.mathworks.com>...

>
> Also, remember that Roger's method uses at least twice as much memory as GRDrej and GRDsim, which becomes a problem for n > 10^8, with 16GB ram.
>

I just take a closer look of Roger's algorithm, and it looks like it is doable to work with a single array (in a C Mex coding). It is also possible to plug a O(N) shuffle engine for it.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ijln5n$o67$1@fred.mathworks.com>...
> "Derek O'Connor" wrote in message <ijlldg$181$1@fred.mathworks.com>...
>
> >
> > Also, remember that Roger's method uses at least twice as much memory as GRDrej and GRDsim, which becomes a problem for n > 10^8, with 16GB ram.
> >
>
> I just take a closer look of Roger's algorithm, and it looks like it is doable to work with a single array (in a C Mex coding). It is also possible to plug a O(N) shuffle engine for it.
>
> Bruno

----------------------

Bruno,

I have just noticed that I omitted the average time for GRDsta2mex, which is Roger's second method with randperm replaced by Jan's GRPmex.

GRDsta2mex: 0.21258 secs, which is more that twice GRDmex, i.e., GRDrej with GRPmex.

Can you tell me exactly what is wrong with Jan's new derangement algorithm and give me counter-examples that I can test or examine.

Please remember, I am not interested in n < 20. I am interested in generating random permutations of various types with n >> 10^6.

Derek

"Derek O'Connor" wrote in message <ijm2l7$pco$1@fred.mathworks.com>...

>
> I have just noticed that I omitted the average time for GRDsta2mex, which is Roger's second method with randperm replaced by Jan's GRPmex.
>
> GRDsta2mex: 0.21258 secs, which is more that twice GRDmex, i.e., GRDrej with GRPmex.
>
> Can you tell me exactly what is wrong with Jan's new derangement algorithm and give me counter-examples that I can test or examine.
>
> Please remember, I am not interested in n < 20. I am interested in generating random permutations of various types with n >> 10^6.

Derek, why arguing back and forth I think I have clearly elaborated my position.

You are free to trust Jan's algorithm at larger n with experimental testing. That's after all your choice.

Here are the facts that I stick with:
1) It's unproven algorithm.
2) It's biased at small n.

Bruno

Dear Bruno,

> > Can you tell me exactly what is wrong with Jan's new derangement algorithm ...

> 1) It's unproven algorithm.
> 2) It's biased at small n.

3) Matt Fig's argument shows that there is a bias for large n also.
My algorithm starts from an equally distributed permutation and swaps the conflicting elements. But to get to the equally distributed set of derangements, a different number of swaps is needed. But my algorithm applies the minimal number of swaps ever. Therefore the result cannot be equally distributed.

Although you, Bruno, might have a higher requirements to the formal composition of a proof, I estimate that this algorithm *is* proven - to create biased derangements for small and large n.

Kind regards, Jan

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ijd8j0$8uq$1@fred.mathworks.com>...
> Roger, just a quick question:
> Do you decompose a permutation in random cycles and find a derangement as permutation such that no cycle is length-one?
> Bruno

"Jos (10584)" wrote in message <ijdtco$32b$1@fred.mathworks.com>...
> Hi Roger,
> Isn't this just the Martinez algorithm? ......
> ~ Jos
- - - - - - - - - -
  This is in answer to Bruno's and Jos' questions Tuesday concerning the first algorithm I presented in this thread. Both this algorithm and the Martinez algorithm are similar in that they establish partial cycles, one link at a time. However they differ in that mine only has one partial cycle being developed at any given time whereas the Martinez algorithm may have a number of partially completed cycles being developed concurrently. As it turns out, this makes no difference whatever to the probability situation involved. The decision to complete a cycle depends only on the number of remaining links yet to be established, and this is the same for the two algorithms. The probability sequences used are therefore exactly the same in both cases. For example they would both terminate with .... 44/309, 9/53, 2/11, 1/3, and 0 probabilities for completing cycles as the number of elements
to be linked-in decreases to zero.

  Another difference is that I use a separate vector, B, as a source for supplying additional elements for linking into the A vector, whereas the Martinez algorithm has swapping in place using only the one vector. This latter saves the memory required for the second vector but requires a "repeat ... until" kind of processing to avoid "marked" elements at each step. Fortunately the number of marked elements is typically in a minority, so the "repeat" action does not usually have to take many steps.

  After reading Jan's original article I became sufficiently intrigued by the problem to independently work out the basics of the above algorithm without reading any further. After later learning of the Martinez algorithm I finally decided there were enough differences between the two to warrant putting it into this thread for others to consider in spite of apparent similarities.

  The theory for each of these algorithms can be considered as firmly established as to generating uniformity. The second "vectorized" method I presented, though seemingly similar, actually is a somewhat different procedure, and I have not been able to work out a rigorous proof of uniformity - only an intuitive indication. However all the tests I have made for small n show that it produces the same kind of results that occur in the first algorithm. I did a million runs for the n = 6 case with both algorithms, and both distributions looked much as one would expect for a statistically uniform 265-valued distribution.

Roger Stafford

"Roger Stafford" wrote in message <ijniue$9so$1@fred.mathworks.com>...
> .......
> The theory for each of these algorithms can be considered as firmly established as to generating uniformity. The second "vectorized" method I presented, though seemingly similar, actually is a somewhat different procedure, and I have not been able to work out a rigorous proof of uniformity - only an intuitive indication. However all the tests I have made for small n show that it produces the same kind of results that occur in the first algorithm. I did a million runs for the n = 6 case with both algorithms, and both distributions looked much as one would expect for a statistically uniform 265-valued distribution.
>
> Roger Stafford
- - - - - - - - -
  Yesterday I stated that I had "only an intuitive indication" that the vectorized algorithm which I called 'derange2' yields uniformity over all possible derangements. After thinking about the matter overnight I now realize that this algorithm actually does possess a firm theoretical justification for its uniformity which I will now attempt to explain.

  The vector p obtained with "p = randperm(n);" in the code is to be regarded as being in cyclic notation with the divisions between cycles indicated by the b vector. For example suppose p and b at this point in the code happen to be:

 p = 4 7 2 8 5 9 6 1 3
 b = 1 0 0 1 0 0 0 1 0 1

Then this p would represent the cyclic notation (4 7 2)(8 5 9 6)(1 3) for a permutation. The vector q generated from this b will be

 q = 2 3 1 5 6 7 4 9 8

and this in turn serves in "p(p) = p(q);" to translate the cyclic notation interpretation of p into the corresponding permutation vector:

 p = 3 4 1 7 9 8 2 5 6.

  If we suppose that the initial values 4, 8, and 1 in each of the above cycles had been obtained, not at random, but somehow by the same deterministic rules as that used in what I called the 'randderange' algorithm (namely the rightmost remaining element at B(k) at the k-th step) we would have exactly the same procedure with the same numerical results (assuming corresponding 'rand' values to be the same) even though carried out by rather different coding. The crucial point to realize here (and one that didn't occur to me yesterday) is that it makes absolutely no difference in the statistical behavior of this type of processing how the initial element in each successive cycle is to be chosen - be it the largest element of those remaining; the smallest; chosen by randderange; chosen at random by 'randperm'; or chosen by any other wild scheme, biased or unbiased - as long as it is some
selection from among the remaining unused elements. The reason for this is that after closing one cycle, the remaining cycles in a cyclic notation can always be reordered and their contents rotated so that an arbitrary element can be chosen to be next in the list without altering the resulting permutation. That of course is not true for subsequent elements in a given cycle where different choices would yield different permutations.

  From this argument we can conclude that 'derange2' does indeed yield a uniform distribution of derangements, just as 'randderange' does.

Roger Stafford

Got you Roger, awesome.

Bruno

I have started a new thread for another derangement generator here:

http://www.mathworks.com/matlabcentral/newsreader/view_thread/304729

Derek O'Connor