How to eliminate equal matrices that are created randomly inside loop?

The code segment I'm working on is given below:

        NphaseSteps = 6;
        phases = exp( 2*pi*1i * (0:(NphaseSteps-1))/NphaseSteps );
        i = 1;
        while i <= 10  %number of iterations
            ind = randi([1 NphaseSteps],10,10);
            inField{i} = phases(ind);
            save('inField.mat', 'inField')
            i = i + 1;
        end

Now, what I want is to keep track of these randomly created matrices "inField{i}" and eliminate the ones that are equal to each other. I know that I can use "if" condition but since I'm new to programming I don't know how to use it more efficiently so that it doesn't take too much time. So, I need your help for a fast working program that does the job. Thanks in advance.

Answer 1

[Walter Roberson]

Adapting from my https://www.mathworks.com/matlabcentral/answers/339203-permutation-of-an-array#comment_452147

A = phases;
n = 10; m = 10;  %each output will be 10 x 10

rows_needed = 10;  %you need 10 of them
cols_needed = n*m;

permidx = [];
more_needed = rows_needed;
while more_needed > 0
  trial_idx = randi([1 NphaseSteps], more_needed, cols_needed);
  permidx = unique([permidx; trial_idx], 'rows', 'stable');
  more_needed = rows_needed - size(permidx,1);
end

randperms = A(permidx);

inField = cellfun(@(M) reshape(M, m, n), num2cell(randperms, 2), 'uniform', 0);

[Bruno Luong]

Sorry but Walter's code is unreliable: if cols_needed > NphaseSteps^(n*m), e.g. cols_needed = 20, NphaseSteps=2, n=2,m=2, the while loop never end. It's unable to capture it's impossible realisation.

If cols_needed <= NphaseSteps^(n*m) but by small amount, it's still possible to end, but it might take very long time to end. Try with n = 4; m = 4; NphaseSteps = 2; rows_needed = 60000;

[Walter Roberson]

For NPhaseSteps = 6, n = m = 10, rows_needed could approach 6E77 .

If you ask for more unique matrices than can be made from NPhaseSteps^(n*m) then you are going to have a problem no matter what approach you use.

There are other approaches, including creating all NPhaseSteps^(n*m) possible arrangements of indices, and then picking rows_needed of them ... and getting stuck if there are not as many as rows_needed unique variations.

There are techniques to generate random permutations without repetition and without generating them all, but they are a bit tricky to implement efficiently.

[Bruno Luong]

Two parameters randperm(N,K) exists for few versions back calling it is not tricky at all.

Now in the case of m=n=10 and NPhaseSteps = 6, N is larger than 64-bit integer, and one cannot even call randperm; but essentially there is virtually 0 chance of collision, the original code works just fine without even checking it.

If you want reliable generic case, RANDPERM has to be used.

[Walter Roberson]

The code I linked to, that I adapted this from, implemented a multiple row version of randperm for efficiency. randperm as-is is only for one row at a time. The task in this question does not call for randperm() because it permits duplicates within any one row -- it has to, because it wants 10*10 = 100 values drawn from the range 1 to 6.

When I say it is difficult to generate random permutations without repetition, I am not talking about any one permutation: I am talking about multiple distinct permutations over sets of objects, and I am talking about methods that deliberately avoid or deliberately reduce the technique of generate-and-reject-duplicates -- techniques that have a high probability of making progress even when most of the possible candidates have already been generated, that have sub-exponential waiting time as the occupancy nears completion.

Somewhat recently there was a discussion of the implementation of randperm() in which one of the Mathworks staff (Steven maybe?) indicated that an algorithm was used for rand(N,K) that was not the previous algorithm of sort(rand) and pick the first K entries. Unfortunately I do not recall the name of the mathematical technique involved.. "dekekind shuffle" or something like that.

I guess you could do,

idx = randperm(NPhaseSteps^(m*n), rows_needed);
idxnum = dec2base(idx, NPhaseSteps, m*n) - '0' + 1;    %more work is needed if nPhaseSteps exceeds 10
randperms = A(idxnum);
inField = cellfun(@(M) reshape(M, m, n), num2cell(randperms, 2), 'uniform', 0);

and count upon randperm() being efficiently implemented for this case. But this will not actually work as randperm() can only handle 2^53-1 for the first parameter, not something in the range between 2^258 and 2^259.

[Bruno Luong]

They probably use Fisher–Yates or Knuth shuffle algorithm for randperm. TMW might realizes the drawback of the old one-parameter randperm (based on SORT) when there is a lot of discussion about lacking of such efficient algorithm in the case k<<n. Jan have a nice FEX implementation of such shuffle.

Your algo is certainly not sub-exponential in all case since it stucks forever for some circumstances (already given above).

[Walter Roberson]

The code I gave in this original Answer was never intended to have to deal with the case where the number of rows being requested exceeded the number of unique rows possible. It would not be difficult to put in

assert(rows_needed <= NPhaseSteps^(m*n), 'more rows requested than uniquely exist')

It would probably start to get seriously inefficient at roughly 90% of NPhaseSteps^(m*n)

Of course, 90% of 6^(10*10) would far exceed the amount of memory anyone is likely to have.

For NPhaseSteps up to 255 and modifying to

    trial_idx = randi([1 NphaseSteps], more_needed, cols_needed, 'uint8');

then the permidx matrix would not be the major memory roadblock: the major memory roadblock would be A(permidx) as that would be 8 times as large since A is double precision. The num2cell of that would take even more memory as cell is not especially efficient.

If we postulate that the user has 64 gigabytes of memory, and making the simplification of assuming that the num2cell version only takes the same amount of memory as the A(permidx) array, then as there would temporarily be two arrays of maximum size, each array could be no more than 32 gigabytes. The permidx matrix would be 1/8th as large, so 4 gigabytes. Each row would be 100 bytes... so you could request on the order of 2^25 rows before running out of memory in that circumstance.

There are algorithms that would generate directly into cells rather than generating into memory and needing to duplicate that into a cell, but I worry about the cost of the uniquifying step. Hmmmm, maybe if you used a hash table of NPhaseSteps^(m*n) bits... That would only require about 1/8th of the total baryons in the universe for the hash table.