Out of memory error in combination combnk

Hi

I want to get all combinations of the n elements in v taken k at a time. I use combnk(v,k). For example combnk(1:121,10). But it is a huge matrix and MATLAB return an "out of memory" error.

I want to use each of the rows of the c in a loop. Is there any way to get each row of the c (one by one) in the loop to avoid out of memory error.

Thanks a lot

M. Rajaei

[Alejandro Casares]

Isn't there some way to modify nchoosek or another similar function to make it write the huge matrix to disk?

Then, the user could read it row by row and process each one independently, without overflowing memory.

A. Casares M.

[DGM]

If you can assume each combination of N numbers can be somehow packed into a single byte, then you're still going to need over 100 terabytes of disk space. That doesn't seem very practical either.

Answer 1

[Roger Stafford]

The total number of rows you would get from combnk(1:121,10) is 121!/10!/111! which is 1.2652e14 which is a very large number indeed. Even if you produced them in a loop one-at-a-time, the loop would require a very long time to finish up. Are you sure this is what you want to do?

A rather ugly way to generate them one-at-a-time would be in nested for-loops ten deep:

 for i1 = 1:121-9
  for i2 = i1+1:121-8
   for i3 = i2+1:121-7
    for i4 = i3+1:121-6
     ....
        for i10 = t9+1:121
         % Now process the set of values i1,i2,i3,i4,...,i10 as desired
        end
     ....
    end
   end
  end
 end

[Mohsen Rajaei]

Thank you Roger for your answer.

Answer 2

[Jos (10584)]

Questions

1. Do you need them all 1.2652e14 at the same time? (Impossible!)
2. Do you need only a random fraction of these? (Easy!)
3. Do you need specifically the K-th ones? (Do-able)

Answers

1. See Roger's answer
2. Here is some code:

V = 1:121 ;
K = 10 ;
N = 5 ; % I want five random combinations
[~, x] = sort(rand(N,numel(V)),2) ;
x = sort(x(:,1:K),2) ;
R = V(x) % N random combinations of K unique elements from V

3. upon request

[Mohsen Rajaei]

Thank you, Roger and Jos.

I should explain my problem more. I need to find roots of the equation det(Ax+B)=0. A and B are K-by-K matrices and B is diagonal. x is multiplied element by element with A.

For example for K=2:

f(x)=det(a11*x+b1,a12*x;a21*x,a22*x+b2)=0

The function f is a polynomial of the degree K. If I find the coefficients of this polynomial, I can find the roots of the equation by roots() in MATLAB. To find these coefficients, I wanted to know those combinations mentioned above.

[Roger Stafford]

Obviously the approach you proposed with 'combnk' is not going to work and you should therefore abandon it.

If you are seeking the real roots to your equation, I would suggest that you make use of matlab's 'fzero' function. You can first do a plot of it over a sufficiently broad span of x values so as to include all zero crossings and that will give you a set of very approximate roots which you can then use for the x0 estimates required by 'fzero' to make the roots accurate

[Roger Stafford]

How large a value of K did you have in mind for your problem? If you are determined to find the coefficients of your polynomial, it could be done as follows. For example, if K is equal to 4, the coefficient of x^2 could be found as the sum of all six possible products of two b's and the corresponding 2 x 2 determinant of the a's with the rows and columns that the two b's occupied removed:

 b1*b2*det([a33,a34;a43,a44]) +
 b1*b3*det([a22,a24;a42,a44]) +
 .... +
 b3*b4*det([a11,a12;a21,a22])   % Six terms altogether

Note however, that for large K this would need to be done using the appropriate for-loops. For example if K were equal to 11, the coefficient of x^5 would be formed from all possible 462 products of 6 of the 11 b's together with the associated reduced 5 x 5 determinant of a's, and this is too much to have to write out by hand. You could use 'combnk' for this purpose, but hopefully not for anything as large as K = 121.

[Mohsen Rajaei]

Thank you Roger,

I have tested 'fzero' function, but I had some problems with it. So I tried to use the approach I mentioned above (using combinations). But it seems that it is impossible to use this approach :)

K comes form (2M+1)(2N+1). for M=N=5, K=121. I need 4 or 5 first roots of the equation (and I know that they are less than 1).

My problem with 'fzero' is that, because of large dimensions of the matrix (K=121), the determinant for some x's becomes Inf and 'fzero' returns an error and exits running the code.

My code is something like this:

 f=@(x) det(A.*x+B)

 for i=linspace(0,1,100)
     if abs(f(i))<Inf
          r(n)=fzero(f,i)
          n=n+1;
     end
end

In my equation there is a parameter kx. For some kx's, for all 100 i's (first points to find root), f(i) becomes Inf and the code exits the loop without finding any root!

[Roger Stafford]

If you are having trouble with your determinant overflowing to infinity, you might try rescaling all the elements in A and B by the same scale factor, s, less than 1. If they all have the same rescaling, that won't affect the roots you obtain for x. Be careful not to make the scale factor too small or you will have trouble underflowing to zero. It could be a delicate adjustment.

Question: What are the significances of M and N in K = (2*M+1)*(2*N+1)? There may be some way of avoiding these huge 121 by 121 matrices which would involve only the sizes M and N. But we can't speculate on this unless you tell us that connection.

Answer 3

[Jozef Lukac]

You can generate a specific one combination corresponding to an index idx by this.

N = 6;
k = 4;

pasc = [ones(N,1) zeros(N,k-1)]; %pascal triangle cropped
for i=2:N
	for j=2:k
		pasc(i,j) = pasc(i-1,j-1) + pasc(i-1,j);
	end
end
combs_forCheck = nchoosek(1:N,k)
pasc
%% compute idx-th combination

idx = 7;


auxIdx = idx;
pascK = k; %aux. k index in pascal triang.
pascN = N; %aux N index in pascal triang.
currDig = 1; %current digit
possOut = zeros(1,k); %idx-th combination
while(pascK > 0)
	while(pascK>0 && pascN>0 && auxIdx <= pasc(pascN,pascK))
		%move up-left the pascal-trinagle matrix
		possOut(k-pascK+1) = currDig;
		currDig = currDig + 1;
		pascN = pascN - 1;
		pascK = pascK - 1;
	end
	if pascK<1
		break;
	end
	while(pascN>0 && auxIdx > pasc(pascN,pascK))
		%move up the pascal-trinagle matrix
		auxIdx = auxIdx - pasc(pascN,pascK);
		currDig = currDig + 1;
		pascN = pascN - 1;
	end
end
%possOut -- output combination
fprintf(1,'%d, ',possOut); fprintf(1,'\n');
fprintf(1,'%d, ',combs_forCheck(idx,:)); fprintf(1,'\n'); %just check