Tips on generating a subset of linspace(0,1,k)^n for "large" n and k.

Hi all,

I am interested in producing all the elements of a particular set, which I will call B. It is defined as follows:

First, let A:=linspace(0,1,k) where .

Then, define , where .

For "small" n and k (i.e. less than 10 and 11, respectively), I am able to produce every vector in B by first using ndgrid to produce and then using find. However, when k=10, this method stops working (computer runs out of memory). An alternative method that avoids generating first would be appreciated.

[Jan]

Please post your code This is better than a description as text. If the result exhausts your RAM, other methods to find the results will not help in any way.

[Frederick Papazyan]

Yes, I should and I will in the future. My apologies!

And just to clarify: the result was not exhausting my ram, a (likely uneeded) intermediate step was (producing A^n)

Answer 1

[Jan]

See: FEX: VChooseKR (C-Mex) :

k = 14;
n = 12;  % <= k
A = linspace(0, 1, k);
P = VChooseKR(uint8(k:-1:1), n);
B = A(P);  % ??? is this really needed?

Now the output P is an UINT8 array. Applying it as index by B=A(P) creates B as double array, which needs 8 times more memory. If RAM is the limiting factor, creating B explicitly is the wrong approach. Stay at P for the further processing.

n = 11 and k = 10 produces a small output only, therefore I cannot understand, where in your code the problem appears. The shown code takes 0.004 sec on my computer and the output P has 1.85MB only, while A(P) would have 14.8 MB.

For n=14, k=12 the code needs 0.44 seconds (i7, 16 GB RAM, Matlab 2018b).

[Walter Roberson]

Would this be in decreasing order across the rows as requested, or would it be in increasing order?

If it is in increasing order the obvious fix would be A = linspace(1,0,k)

[Jan]

@Walter Roberson: The indices are produced in decreasing order by the posted code:

P = VChooseKR(uint8(k:-1:1), n)

This would create the increasing order:

P = VChooseKR(uint8(1:k), n)

Changing the order of A would be efficient also.

[Frederick Papazyan]

Thank you so very much Jan and Walter! I really do appreciate the (large amount) of work you two put in to answering my question.

Answer 2

[Walter Roberson]

This can certainly still run out of memory.

The number of outputs is k! / (k-n)! and they take n columns each, and need to be double precision for the final output. The intermediate calculation of indices requires k! / (k-n)! * n but they can be uint8 unless you are going for k >= 256 and n fairly small.

I ran k=18 n=14 on my system; B came out 124062000 x 14 for 13 gigabytes, and took about 5 1/2 minutes. I had not yet made the change to use uint8 indices so it would have required I as large as B was, and would have required building the parts, T, with a total memory as large -- so the peak memory demand would have been on the order of 40-ish gigabytes. The change here to use uint8 indices should reduce that significantly.

tic
k = 14;
n = 12;
A = linspace(0,1,k);

I = uint8(1:k-n+1).';
for N = 2 : n
    maxidx = k-n+N;
    nI = size(I,1);
    T = cell(nI,1);
    for J = 1 : nI
        H = uint8(I(J,1):maxidx).';
        T{J} = [H, repmat(I(J,:), length(H),1)];
    end
    I = vertcat(T{:});
end
B = A(I);
toc

Elapsed time is 3.834880 seconds.

size(B)

ans = 1×2

1931540 12

B(1:10,:)

ans = 10×12

0 0 0 0 0 0 0 0 0 0 0 0 0.0769 0 0 0 0 0 0 0 0 0 0 0 0.1538 0 0 0 0 0 0 0 0 0 0 0 0.2308 0 0 0 0 0 0 0 0 0 0 0 0.3077 0 0 0 0 0 0 0 0 0 0 0 0.3846 0 0 0 0 0 0 0 0 0 0 0 0.4615 0 0 0 0 0 0 0 0 0 0 0 0.5385 0 0 0 0 0 0 0 0 0 0 0 0.6154 0 0 0 0 0 0 0 0 0 0 0 0.6923 0 0 0 0 0 0 0 0 0 0 0

[Walter Roberson]

Using uint8 indices sped up the calculation a fair bit too; 18/14 only took about 4:40 instead of around 5:30

[Jan]

I'm confused. x1>=x2>=x3...>=xn means, that all x can have the same value. For k05 and n=3 your indices I contain [1,1,1], but not [5,5,5]. Is this correct?

H = uint8(I(J,1):maxidx) creates a double vector and casts it afterwards. To get the UINT8 data directly:

H = (uint8(I(J,1)):uint8(maxidx)).'

[Walter Roberson]

Ah, I think I programmed for strict > not for >=

tic
k = 14;
n = 12;
A = linspace(0,1,k);
maxidx = uint8(k);
I = (uint8(1):maxidx).';
for N = 2 : n
    nI = size(I,1);
    T = cell(nI,1);
    for J = 1 : nI
        H = (I(J,1):maxidx).';
        T{J} = [H, repmat(I(J,:), length(H),1)];
    end
    I = vertcat(T{:});
end
B = A(I);
toc

Elapsed time is 21.067155 seconds.

size(B)

ans = 1×2

5200300 12

B(1:10,:)

ans = 10×12

0 0 0 0 0 0 0 0 0 0 0 0 0.0769 0 0 0 0 0 0 0 0 0 0 0 0.1538 0 0 0 0 0 0 0 0 0 0 0 0.2308 0 0 0 0 0 0 0 0 0 0 0 0.3077 0 0 0 0 0 0 0 0 0 0 0 0.3846 0 0 0 0 0 0 0 0 0 0 0 0.4615 0 0 0 0 0 0 0 0 0 0 0 0.5385 0 0 0 0 0 0 0 0 0 0 0 0.6154 0 0 0 0 0 0 0 0 0 0 0 0.6923 0 0 0 0 0 0 0 0 0 0 0

[Jan]

Yes, now the indices contain [1,1,1] and [5,5,5] for k=5 and n=3. The number of outputs is (K+N-1 over N).