Find combination of 12 matrix rows, fulfilling certain criterion, from matrix of size 3432 *7

M=nchoosek(1:14,7) gives matrix of size 3432 * 7.

I want combinations of 12 rows from M, such that each of the 14 numbers in the 12-row combinations has frequency 6.

12 rows with 7 numbers are 84 numbers, and 6 *14 = 84. I think there are several 12-row combinations that fulfill the criterion.

How many such combinations are there and how can I get the row numbers of them?

Answer 1

[John D'Errico]

Finding all possible computations there are would be a computationally hard problem, certainly if you used brute force. It would be intractable at that point. Yes, there are probably many such combinations. However, you need to find only ONE combination, because all others could probably be generated from the first. (That is just conjecture on my part, because I'd need to think about it. And counting the number of possible solutions is an interesting question, probably solvable using group theory.)

Anyway, how would I solve it? First, DON'T use nchoosek. Instead, form it as a binary indicator matrix. We want to find the set of numbers with exactly 7 bits set to 1, from the numbers 0:(2^14-1)

B = dec2bin(0:2^14-1) - '0';
B(sum(B,2) ~= 7,:) = [];
size(B)

ans = 1×2

3432 14

As you can see, there are 3432 rows. Each row is one of the rows that nchoosek would generate.

B(1:5,:)

ans = 5×14

0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 1 0 1 1 1 1 1 1 0 0 0 0 0 0 1 1 0 1 1 1 1 1 0 0 0 0 0 0 1 1 1 0 1 1 1 1 0 0 0 0 0 0 1 1 1 1 0 1 1 1

Now, how can we find ONE solution? That part is also simple, or it should be. (Said, before I solve it.) I'll use a binary ILP solver, thus intlinprog. There will be 3432 binary variables. We want the solution where exactly 12 of them are set to 1, the rest zero. See how this will work.

help intlinprog

 INTLINPROG Mixed integer linear programming.
 
    X = INTLINPROG(f,intcon,A,b) attempts to solve problems of the form
 
             min f'*x    subject to:  A*x  <= b
              x                       Aeq*x = beq
                                      lb <= x <= ub
                                      x(i) integer, where i is in the index
                                      vector intcon (integer constraints)
 
    X = INTLINPROG(f,intcon,A,b) solves the problem with integer variables
    in the intcon vector and linear inequality constraints A*x <= b. intcon
    is a vector of positive integers indicating components of the solution
    X that must be integers. For example, if you want to constrain X(2) and
    X(10) to be integers, set intcon to [2,10].
 
    X = INTLINPROG(f,intcon,A,b,Aeq,beq) solves the problem above while
    additionally satisfying the equality constraints Aeq*x = beq. (Set A=[]
    and b=[] if no inequalities exist.)
 
    X = INTLINPROG(f,intcon,A,b,Aeq,beq,LB,UB) defines a set of lower and
    upper bounds on the design variables, X, so that the solution is in the
    range LB <= X <= UB. Use empty matrices for LB and UB if no bounds
    exist. Set LB(i) = -Inf if X(i) is unbounded below; set UB(i) = Inf if
    X(i) is unbounded above.
 
    X = INTLINPROG(f,intcon,A,b,Aeq,beq,LB,UB,X0) sets the initial point
    to X0.
 
    X = INTLINPROG(f,intcon,A,b,Aeq,beq,LB,UB,X0,OPTIONS) minimizes with
    the default optimization parameters replaced by values in OPTIONS, an
    argument created with the OPTIMOPTIONS function. See OPTIMOPTIONS for
    details.
 
    X = INTLINPROG(PROBLEM) finds the minimum for PROBLEM. PROBLEM is a
    structure with the vector 'f' in PROBLEM.f, the integer constraints in
    PROBLEM.intcon, the linear inequality constraints in PROBLEM.Aineq and
    PROBLEM.bineq, the linear equality constraints in PROBLEM.Aeq and
    PROBLEM.beq, the lower bounds in PROBLEM.lb, the upper bounds in
    PROBLEM.ub, the initial point in PROBLEM.x0, the options structure in
    PROBLEM.options, and solver name 'intlinprog' in PROBLEM.solver.
 
    [X,FVAL] = INTLINPROG(f,intcon,A,b,...) returns the value of the
    objective function at X: FVAL = f'*X.
 
    [X,FVAL,EXITFLAG] = INTLINPROG(f,intcon,A,b,...) returns an EXITFLAG
    that describes the exit condition. Possible values of EXITFLAG and the
    corresponding exit conditions are
 
      3  Optimal solution found with poor constraint feasibility.
      2  Solver stopped prematurely. Integer feasible point found.
      1  Optimal solution found.
      0  Solver stopped prematurely. No integer feasible point found.
     -1  Solver stopped by an output function or plot function.
     -2  No feasible point found.
     -3  Root LP problem is unbounded.
     -9  Solver lost feasibility probably due to ill-conditioned matrix.
 
    [X,FVAL,EXITFLAG,OUTPUT] = INTLINPROG(f,A,b,...) returns a structure
    OUTPUT containing information about the optimization process. OUTPUT
    includes the number of integer feasible points found and the final gap
    between internally calculated bounds on the solution. See the
    documentation for a complete description.
 
    See also LINPROG.

    Documentation for intlinprog
       doc intlinprog

N = 3432;
intcon = 1:N;
LB = zeros(N,1);
UB = ones(N,1);

f = ones(N,1); % this is pretty much irrelevant, since we just want to find a feasible solution

Aeq1 = ones(1,N);
beq1 = 12; % EXACTLY 12 of the rows will be chosen

Aeq2 = B';
beq2 = repmat(6,14,1); % this enforces that the set chosen will be perfectly balanced

[xsol,fval,exitflag] = intlinprog(f,intcon,[],[],[Aeq1;Aeq2],[beq1;beq2],LB,UB);

LP:                Optimal objective value is 12.000000.                                            


Optimal solution found.

Intlinprog stopped at the root node because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0. The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05.

The rows chosen were:

find(xsol)

ans = 12×1

876 925 1000 1321 1442 1606 1858 2069 2103 2431

B(logical(xsol),:)

ans = 12×14

0 1 0 0 0 1 1 1 1 1 1 0 0 0 0 1 0 0 1 0 1 1 1 0 0 1 1 0 0 1 0 0 1 1 1 1 1 0 0 1 0 0 0 1 1 0 0 1 0 0 1 1 0 1 0 1 0 1 1 0 1 0 1 1 0 1 0 0 0 1 0 1 1 1 0 1 0 1 0 1 1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 1 1 1 1 0 0 1 1 0 0 1 0 0 1 0 1 1 1 0 0 1 1 0 1 0 1 0 1 1 0 0 1 0 1 1 0 0 0 0 0 0 1 1 1 1

You can now easily enough convert that into an output of the form nchoosek would generate.

Can you generate other random solutions? This should work, to find a different randomly selected set:

f2 = rand(N,1);
[xsol,fval,exitflag] = intlinprog(f2,intcon,[],[],[Aeq1;Aeq2],[beq1;beq2],LB,UB);

LP:                Optimal objective value is 0.054465.                                             

Cut Generation:    Applied 3 strong CG cuts,                                                        
                   and 1 zero-half cut.                                                             
                   Lower bound is 0.057734.                                                         

Heuristics:        Found 2 solutions using total rounding.                                          
                   Upper bound is 0.215474.                                                         
                   Relative gap is 12.98%.                                                         

Heuristics:        Found 2 solutions using diving.                                                  
                   Upper bound is 0.058546.                                                         
                   Relative gap is 0.08%.                                                          

Branch and Bound:

   nodes     total   num int        integer       relative                                          
explored  time (s)  solution           fval        gap (%)                                         
       3      0.45         4   5.854572e-02   0.000000e+00                                          

Optimal solution found.

Intlinprog stopped because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 0. The intcon variables are integer within tolerance,
options.IntegerTolerance = 1e-05.

find(xsol)

ans = 12×1

34 113 405 1198 1258 1617 1985 2148 2352 3024

Counting the number of possible solutions, or finding all such possible solutions is another question, as I suggest, probavbly in the realm of group theory, but I have work at home for now. ;-)

[Ernesto]

John D'Errico thanks, but I thinkk your answer is wrong.

My row numbering is:

1: 1 2 3 4 5 6 7

2: 1 2 3 4 5 6 8

.

.

3432: 8 9 10 11 12 13 14

If you check your two solutions I think you won't get a frequency of 6 for the 14 numbers, assuming my number ordering.

Two examples I have found:

176

326

514

687

1146

1437

1996

2287

2746

2919

3107

3257

****

1153

1355

1458

1553

1581

1620

1813

1852

1880

1975

2078

2280

****

Answer 2

[David Goodmanson]

Hi Ernesto,

Once you have a single 12x7 solution you can make many more by random element swapping of elements in that matrix. This will eventually create every possible 12x7, although not in a deterministic fashion. The idea is that for two rows of the 12x7, if element b in row p is not contained in row q, and if element c in row q is not contained in row p, then you can swap elements and get a new solution. The new rows will still have distinct elements and the rule that the number of ones = 6, twos = 6 etc. is not affected.

For each 1x7 row in the matrix, the code below can calculate the row index in the 3432x7 matrix created by nchoosek. Since the total number of rows is only 3432, I chose to use a simple method to determine the row index rather than use some iterative process involving euclidean remainders or the like.

% construct a solution
a = [1:7]+[0:5]';
b = mod([5:-1:1]-[1:6]',12)+1;
c = repmat(13:14,6,1);
n1 = sort([a;[b c]],2)

% 1000 swaps (way more than you really need) 
n = n1;
for k = 1:1000
    r1 = randi(12);
    r2 = randi(12);
    c1 = randi(7);
    c2 = randi(7);
    e1 = n(r1,c1);
    e2 = n(r2,c2);
    if ~ismember(e1,n(r2,:)) & ~ismember(e2,n(r1,:))
        n(r1,c1) = e2;
        n(r2,c2) = e1;
    end
end
n = sort(n,2)    % the result

% find the row index for 3432x7 nchoosek, using the 6th row of n as an example
n(6,:)
ind = fun(n(6,:))
% check
nch = nchoosek(1:14,7);  
nch(ind,:)            % agrees with n(6,:)
return

function ind = fun(rowvec)
% row index of a given nchoosek(14,7) row vector
nch = nchoosek(1:14,7);
nchx =  nch*(17.^(7:-1:1))';
nx = rowvec*(17.^(7:-1:1))';
ind = find(nx==nchx);
end

You can also create another 12x7 matrix nnew from n with with

p = randperm(14);
nnew = p(n)  

This changes the elements more quickly than the method above, but it merely maps old numbers into new ones and does not change the basic stucture of n.