How to generate all possible (10x10) matrices containing only 0s and 1s, along with other restrictions, for a Graph Theory application.

***EDIT TO REQUIREMENT 6 AND NEW REQUIREMENT ADDED***

6) Exactly 4 columns/rows must have degree 3.

7) No two vertices of degree 3 are adjacent to each other.

My goal:

To generate and save all matrices that meet specific requirements. Then compare each matrix to additional matrices that have been manually entered previously to check for specific similarities. I can add more detail if somebody thinks it would be helpful. I believe I have the comparison aspect of the code sorted out already, so I am waiting on the matrix generation portion. I need to do this for multiple sizes but I will focus this question on the 10x10 case.

Specific requirements:

1) Must be a 10x10 matrix (representing a graph on 10 vertices).

2) Must be symmetric (representing an adjacency matrix).

3) Have a diagonal of 0s (no loops).

4) Only 1s and 0s (simple graph).

5) The entire matrix must have exactly 48 1s (the graph has 24 edges).

6) Each column/row must have either 3 or 6 1s (each node as degree 3 or 6).

Application:

I am investigating a conjecture and believe I have come up with a possible solution which could break down the conjecture into smaller pieces and possibly prove some aspects. I want to use brute force to show if my idea works for a small specific case. Also having a base code in place could allow for future modifications to test other possible cases or ideas.

Ideas and thought process:

- I used the edges of a graph to manually input my comparison set. For example:

G9=graph([1 1 1 2 2 3 4 4 4 5 5 6 6 6 6 3 3 9 2 2 2 7 7 8],[2 3 4 3 4 4 5 6 7 6 7 7 3 9 10 9 10 10 7 8 9 8 9 9]);

- I think this is the only graph, up to isomorphism, which meets the previously listed requirements.

- My original thought was to create the possible matrices that satisfy the given conditions then compare them to my comparison set. I still think this is the best approach.

- I foolishly attempted to generate random matrices, completely overlooking the massive number of possibilities. Using a while loop, I first generated a random matrix that satisfied the first four requirements. Then in separate nested for statements I checked for requirement 5 using numedes() and requirement 6 using all(mod(degree())). That was a bad approach for several fairly obvious reasons, but I learned a lot through the process and it led me to the code that should do my final comparisons.

- This is the first time I have used Matlab so I am learning as I go. I have been working on this one code for nearly 2 weeks and do not know if what I have come up with is "good", but I am proud of what I have been able to do by myself. I have reached the point where I feel like I need some outside advice. I have looked through scores of previous questions but can not seem to find anything that I have been able to successfully work with. I am open to any suggestions and any level of help. A reference to a source, a function suggestion, another approach, or a complete solution with a "plug and play" code would be appreciated. I do not shy away from putting forth the effort to achieve my goals.

I appreciate any feedback.

[llueg]

My first instinct was to pose this as an optimization problem. To do this easily in Matlab, you should have the Optimization Toolbox installed. If you do not have that, modeling will be more complicated.

Below, I have basically formulated each of the requirements you mentioned as a constraint, and the matrix as a binary variable of size 10x10. I have added decision variables to switch between rowsums of 3 and 6.

I set the objective function of the optimization problem to 0, so that any feasible solution can be found as a result. The option @savemilpsolutions allows you to save all solutions that were found, however when I ran it, only one was returned. But, when I set different arbitrary Objective functions such as min(sum(z)) or max(sum(z)), different solutions were returned. So there seems to be more than one solution. I think the @savemilpsolutions option only keeps track of solutions within the Branch and Bound tree, not in the root preprocessing. But since the problem is solved very quickly, no nodes are explored. I am not sure how to find all feasible points, so maybe this answer is not very useful. But at least you can find a couple of feasible points by modifying the Objective function. I did not check for isomorphism in the solutions I found.

%Contraint Formulation
N = 10;
M = optimvar('M',N,N,'Type','integer','LowerBound',0,'UpperBound',1);

%1st Constraint: Symmetry
Constraints.Msym_constr = optimconstr((N-1)*N/2);
k = 1;
for i=1:N
    for j=i+1:N
            Constraints.Msym_constr(k) = M(i,j) == M(j,i);
            k = k+1;
    end
end

%2nd Constraint: Diagonal of 0s
Constraints.diag_constr = optimconstr(N);

for i=1:N
    Constraints.diag_constr(i) = M(i,i) == 0;
end

%3rd Contraint: 48 1s
Constraints.ones_contstr = sum(sum(M)) == 48;

%4th Contraint: sum in rows = 6 or = 3 -> decision variable z
z = optimvar('z',N,'Type','integer','LowerBound',0,'UpperBound',1);

Constraints.rowsum_constr = optimconstr(N);

for i=1:N
   Constraints.rowsum_constr(i) = sum(M(i,:)) == z(i)*3 + (1-z(i))*6; 
end

%Create problem and solve
problem = optimproblem('ObjectiveSense','max');
problem.Objective = sum(z);
problem.Constraints = Constraints;

options = optimoptions('intlinprog','OutputFcn',@savemilpsolutions);
% The Matrix will be in x.M
[x,fval,exitflag,output]  = solve(problem,'Options',options);

[James Snead]

I truely appreciate the detailed response. I need to do some reading to better understand how and why each step works since I am unfamiliar with Matlab. Also, I need to create some more constraints. The reason I say that is because x.M is not a valid matrix for my problem even though it satisfies all the constraints I listed. The specific case my question refers to requires the graph to be decomposable into four 4-complete graphs. I plotted x.M and the edge set of a graph that does satisfy everything I need;

x.M=graph([0,0,0,0,0,0,0,1,1,1;0,0,0,1,1,1,0,1,1,1;0,0,0,0,1,0,0,1,1,0;0,1,0,0,1,0,1,1,1,1;0,1,1,1,0,1,1,0,0,1;0,1,0,0,1,0,0,0,1,0;0,0,0,1,1,0,0,1,0,0;1,1,1,1,0,0,1,0,0,1;1,1,1,1,0,1,0,0,0,1;1,1,0,1,1,0,0,1,1,0]);
G=graph([1 1 1 2 2 3 4 4 4 5 5 6 6 6 6 3 3 9 2 2 2 7 7 8],[2 3 4 3 4 4 5 6 7 6 7 7 3 9 10 9 10 10 7 8 9 8 9 9]);

figure(1);
subplot(1,2,1);
plot (x.M)
subplot(1,2,2);
plot (G)

If you look in the figure for x.M; 7 is connected to 4,5, and 8. According to my problem those connections imply: 4 must be connected to 5,7, and 8. 5 connected to 4,7, and 8. 8 connected to 4,5, and 7. In G, you can compare that to nodes 1, 2, 3, and 4.

Since G was not generated by your code, I will need to work on the details to see how to produce all feasible solutions. You have given me a lot to think about. I will work on this after work and give an update when I make some progress.

[James Snead]

I have made some progress but not much. This approach seems to give me few, if any, matrices, and trying to generate matrices in a loop gives me too many matrices; most of which are not valid. While using loops, I use H to denote the most recently generated matrix. I then use S=sum(H,1) to create an array of the column sums. Then I can use find() to see which elements of S have a sum of 3, which tells me which columns of H have a sum of 3. Since no two vertices of degree 3 are adjacent to each other, if S(i)=3 and S(j)=3 then H(i,j)=0 and H(j,i)=0. Do you know how to do something similiar using your method? Also, can we make it where exactly 4 columns have a sum of 3?

[Bruno Luong]

If you want to generate only degree 3 one idea is to organize the vertex by triangle latices on a sphere, then randomly permutted them.

Similarly If you want to generate only degree 6 then put the vertexes on a cube-like latices with in a 3-sphere, then randomly permutted them.

It should have a named finite group out there that reflects those ideas.

I have no idea how to mix both and not sure they generate all the possible graphes.

Answer 1

[Bruno Luong]

One simple thing you could do is move outt the while loop all the preparation step that compute variable that does change

graphs=0;
valid_length_Y=0;
valid_length_Y_y=0;
valid_graph=0;
valid_new_graph=0;
new_graph={}; % set of matrices that are valid and not isomorphic to G.
m = 10;
I = triu(ones(m),1);
n = sum(I(:));
I(I>0) = 1:n;
I = I+I';
ne = 48; % number of edges
lb = zeros(n,1);
ub = ones(n,1);
[~, r, y] = find(I);
Aeq = accumarray([r y], 1, [m+1 n]);
beq = zeros(m+1,1); % node degree
p6 = ne/3-m;
Aeq(m+1,:) = 1;
beq(m+1,:) = ne/2;
while graphs<2^25
    beq(1:m) = 3;
    beq(randperm(m,p6)) = 6;
    fdummy = randn(n,1);
    x = intlinprog(fdummy, 1:n, [], [], Aeq, beq, lb, ub);
    ...
end

The second thing is your way of growing cellarray is very bad. If you can't compute the bound of the size I would grow the array as this

new_graph = cell(1,1024);
valid_new_graph = 0;
while ...
    ...
    valid_new_graph=valid_new_graph+1;
    if  valid_new_graph > length(new_graph)
        new_graph(2*end) = {[]}; % grow it exponentially
    end
end
new_graph(cellfun('isempty', newgraph)) = [];

Edit: fix bug of new_graph(2*end) = {[]}

[James Snead]

I originally had something similar to what you have in your first suggestion. I was having several issues and I evenutally noticed if I moved everything inside the while loop the issues were resolved. When I was going through the details I noticed my mistake so I will fix that in a little bit.

I'm not completely sure what you mean by your second suggestion. I have so little experience that understanding the format and functions takes me longer than it probably should, but I will keep looking at. Is this what you where suggesting:

% G,H,I,and J are edge sets I used to test for specific errors
G=full(adjacency(graph([1 1 1 2 2 3 4 4 4 5 5 6 6 6 6 3 3 9 2 2 2 7 7 8],[2 3 4 3 4 4 5 6 7 6 7 7 3 9 10 9 10 10 7 8 9 8 9 9])));
% G error: none
H=full(adjacency(graph([1 1 1 1 1 1 2 2 3 3 4 4 4 5 5 5 5 6 6 6 6 7 7 8],[2 3 5 6 7 8 6 9 8 9 5 7 10 7 8 9 10 7 8 9 10 8 9 9])));
% H error: length(Y)
I=full(adjacency(graph([1 1 1 1 1 1 2 2 2 2 2 3 3 4 4 5 5 5 6 6 7 7 8 8],[2 5 6 7 8 10 3 4 5 7 8 8 10 5 10 7 8 9 7 10 9 10 9 10])));
% I error: length(Y{y})
J=full(adjacency(graph([1 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 5 6 7 7 8],[2 4 7 8 9 10 3 4 6 8 9 4 5 7 8 9 5 6 7 8 7 8 10 10])));
% J error: complete subgraph

k=4;
min_degree=k-1;
graphs=0;
valid_length_Y=0;
valid_length_Y_y=0;
valid_graph=0;
valid_new_graph=0;
new_graph = cell(1,1024);
new_graph={}; % set of matrices that are valid and not isomorphic to G.
nodes = k*(k+1)/2;
I = triu(ones(nodes),1);
n = sum(I(:));
I(I>0) = 1:n;
I = I+I';
degree_sum = k^2*(k-1); % twice the number of edges
lb = zeros(n,1);
ub = ones(n,1);
[~, r, y] = find(I);
Aeq = accumarray([r y], 1, [nodes+1 n]);
beq = zeros(nodes+1,1); % node degree
max_degree = degree_sum/min_degree-nodes;
Aeq(nodes+1,:) = 1;
beq(nodes+1,:) = degree_sum/2;
while graphs<2^25
    beq(1:nodes) = 3;
    beq(randperm(nodes,max_degree)) = 6;
    fdummy = randn(n,1);
    x = intlinprog(fdummy, ones(n,1), [], [], Aeq, beq, lb, ub);
    F={}; % set of complete subgraphs
    Y={}; % set of shared vertices
    f=1;
    y=1;
    pass=1;
    step=1;
    if ~isempty(x)
        x = round(x);
        A = zeros(size(I));
        A(I>0) = x(I(I>0));
        S=sum(A);
        for s=1:nodes       % s is the nodes of degree 3
           if S(s)==3
              a=A(:,s); % a is a column of degree 3
              [row,col,v] = find(a);
              F{f}=sort([row;s]); % F{f} is a complete subgraph of the graph
              f=f+1;
           end
        end
        for u=1:length(F)
            for v=u+1:length(F)
                if ~isempty(intersect(F{u},F{v}))
                    Y{y}=intersect(F{u},F{v});
                    y=y+1;
                end
            end
        end
        if length(Y)~=max_degree
            'Bad Matrix: length(Y)';
            graphs=graphs+1;
            pass=0;
        end
        if pass==1
            valid_length_Y=valid_length_Y+1;
            for y=1:length(Y)
                if length(Y{y})~=1
                    'Bad Matrix: length(Y{y})';
                    graphs=graphs+1;
                    pass=0;
                    continue
                end
            end
        end
        if pass==1
            valid_length_Y_y=valid_length_Y_y+1;
            for o=1:length(F)
                if pass==0
                    break
                end              
                for p=1:length(F)
                    if pass==0
                        break
                    end                    
                    for q=p+1:length(F)
                        if A(F{o}(p),F{o}(q))==1
                            continue
                        else
                            'Bad Matrix: Subgraphs are not complete';
                            graphs=graphs+1;
                            pass=0;
                            break
                        end
                    end
                end
            end
        end
        if pass==1
            valid_graph=valid_graph+1;
            'Valid Matrix';
            graphs=graphs+1;
            if isisomorphic(graph(G),graph(A))==false
                valid_new_graph=valid_new_graph+1;
                'New Graph';
                new_graph{valid_new_graph}=A;                
                if  valid_new_graph > length(new_graph)
                    new_graph(2*end) = {}; % grow it exponentially
                end
                new_graph(cellfun('isempty', new_graph)) = [];
            end
        end
    end
end
graphs
valid_length_Y
valid_length_Y_y
valid_graph
valid_new_graph

There is a general equation for the upper bound of new_graph{}, but generating random matrices gives me duplicates. If I can figure out a way to systematically generate all matrices, then I would have an unpper bound. In that case, would I just replace 1024 with the bound:

new_graph = cell(1,1024); % current
new_graph = cell(1,upper bound) % updated

Also, occationally your matrix code generates a matrix with nodes of degree 7 such that sum(A)=54. My code saves it in new_graph{} so I can easily notice the error and disregard it. I'm not concerned with it since it has only happened about 10 times out of over 10 billion iterations, but I am curious if you know why it might be doing that.

[Bruno Luong]

"Also, occationally your matrix code generates a matrix with nodes of degree 7 such that sum(A)=54."

I put it as a bug of intlinprog, since I impose the sum is 48 in the last row of Aeq. If it gives different result then the fault is from intlinprog. If you have such case please store the input of intlinprog and post it here for us to look at, or report the bug to TMW.

Sorry I believe I might make a bug you should grow the array. The correct statement is:

new_graph(2*end) = {[]};

The statement

new_graph(cellfun('isempty', new_graph)) = [];

should be OUTSIDE the while-loop.

You might also prefer to store the random A matrices in 3D array 10 x 10 x n rather than cell array. It needs a big avalilability of a contiguous memory block, but overall the memory requirement is less than CELL. You can also use fancy stuff such as TALLARRAY, MATFILE object to store the results on Disk.

[James Snead]

If it happens again I'll post it here and report it to TMW. I made the two adjustments and the code is working almost twice as fast as my code. What I have now is enough for this case so I am going to modify it to fit a few more cases. Once I do that I'll come back to this case and work on better storage. I am hoping to have everything completed within a week.

I do have one other question. Right now I need 4 nodes of degree 3 and 6 nodes of degree 6. If I need x nodes of degree 3, y nodes of degree 6, and z nodes of degree 9 could I simply expand the same method that I am currently using or would I need to figure out another approach? At this point, I am reasonably confident I can work through everything after the matrix is generated. The matrix generating portion is what I am still struggling with.

[Bruno Luong]

Yes I'm confident the method can be generalized with three degrees 3/6/9.

Answer 2

[Bruno Luong]

Here is working code that generate a random adjadcent matrix when allowing

• 4 nodes of degree 3
• 6 nodes of degree 6

to be compatible with # of edges of 48

m = 10;
I = triu(ones(m),1);
n = sum(I(:));
I(I>0) = 1:n;
I = I+I';
ne = 48; % number of edge
fdummy = randn(n,1);
lb = zeros(n,1);
ub = ones(n,1);
[~, r, c] = find(I);
Aeq = accumarray([r c], 1, [m+1 n]);
beq = 3 + zeros(m+1,1); % node degree
p6 = ne/3-m;
beq(randperm(m,p6)) = 6;
Aeq(m+1,:) = 1;
beq(m+1,:) = ne/2;
x = intlinprog(fdummy, 1:n, [], [], Aeq, beq, lb, ub);
if ~isempty(x)
    x = round(x);
    A = zeros(size(I));
    A(I>0) = x(I(I>0));
    % Check result
    A
    sum(A(:))
    sum(A)
end

Edit: shorter code (remove for-loop)

[James Snead]

I am sure there are far more efficient ways to do this, but I'll show you what I've done so far. I used your code to generate the random matrix. Then, I checked that matrix for a few requirements to see if it was a valid matrix. If it was a valid matrix, I checked to see if it was isomorphic to a matrix which I know to be valid. If it is not isomorphic then it is stored into a set of other valid matrices. I placed everything in a while loop and added counters to track the number of matrices that fail at specific checkpoints. Your edited code makes a significant difference in run time once I increase the number of matrices I check. My current issues are how long it takes to check a large number of matrices and the fact I am relying on randomly generated matrices. The extended time is not a big deal since I can let it run while I do other things. Using randomly generated matrices will never guarantee that I have checked every possible matrix. Do you have any suggestions for improving this code in any way? Here is my current code:

% G,H,I,and J are edge sets I used to test for specific errors
G=full(adjacency(graph([1 1 1 2 2 3 4 4 4 5 5 6 6 6 6 3 3 9 2 2 2 7 7 8],[2 3 4 3 4 4 5 6 7 6 7 7 3 9 10 9 10 10 7 8 9 8 9 9])));
% G error: none
H=full(adjacency(graph([1 1 1 1 1 1 2 2 3 3 4 4 4 5 5 5 5 6 6 6 6 7 7 8],[2 3 5 6 7 8 6 9 8 9 5 7 10 7 8 9 10 7 8 9 10 8 9 9])));
% H error: length(Y)
I=full(adjacency(graph([1 1 1 1 1 1 2 2 2 2 2 3 3 4 4 5 5 5 6 6 7 7 8 8],[2 5 6 7 8 10 3 4 5 7 8 8 10 5 10 7 8 9 7 10 9 10 9 10])));
% I error: length(Y{y})
J=full(adjacency(graph([1 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 5 6 7 7 8],[2 4 7 8 9 10 3 4 6 8 9 4 5 7 8 9 5 6 7 8 7 8 10 10])));
% J error: complete subgraph

graphs=0;
valid_length_Y=0;
valid_length_Y_y=0;
valid_graph=0;
valid_new_graph=0;
new_graph={}; % set of matrices that are valid and not isomorphic to G.

while graphs<2^25
    m = 10;
    I = triu(ones(m),1);
    n = sum(I(:));
    I(I>0) = 1:n;
    I = I+I';
    ne = 48; % number of edges
    fdummy = randn(n,1);
    lb = zeros(n,1);
    ub = ones(n,1);
    [~, r, y] = find(I);
    Aeq = accumarray([r y], 1, [m+1 n]);
    beq = 3 + zeros(m+1,1); % node degree
    p6 = ne/3-m;
    beq(randperm(m,p6)) = 6;
    Aeq(m+1,:) = 1;
    beq(m+1,:) = ne/2;
    x = intlinprog(fdummy, ones(n,1), [], [], Aeq, beq, lb, ub);
    F={}; % set of complete subgraphs
    Y={}; % set of shared vertices
    f=1;
    y=1;
    pass=1;
    step=1;
    if ~isempty(x)
        x = round(x);
        A = zeros(size(I));
        A(I>0) = x(I(I>0));
        S=sum(A);
        for s=1:m       % s is the nodes of degree 3
           if S(s)==3
              a=A(:,s); % a is a column of degree 3
              [row,col,v] = find(a);
              F{f}=sort([row;s]); % F{f} is a complete subgraph of the graph
              f=f+1;
           end
        end
        for u=1:length(F)
            for v=u+1:length(F)
                if ~isempty(intersect(F{u},F{v}))
                    Y{y}=intersect(F{u},F{v});
                    y=y+1;
                end
            end
        end
        if length(Y)~=p6
            'Bad Matrix: length(Y)';
            graphs=graphs+1;
            pass=0;
        end
        if pass==1
            valid_length_Y=valid_length_Y+1;
            for y=1:length(Y)
                if length(Y{y})~=1
                    'Bad Matrix: length(Y{y})';
                    graphs=graphs+1;
                    pass=0;
                    continue
                end
            end
        end
        if pass==1
            valid_length_Y_y=valid_length_Y_y+1;
            for o=1:length(F)
                if pass==0
                    break
                end              
                for p=1:length(F)
                    if pass==0
                        break
                    end                    
                    for q=p+1:length(F)
                        if A(F{o}(p),F{o}(q))==1
                            continue
                        else
                            'Bad Matrix: Subgraphs are not complete';
                            graphs=graphs+1;
                            pass=0;
                            break
                        end
                    end
                end
            end
        end
        if pass==1
            valid_graph=valid_graph+1;
            'Valid Matrix';
            graphs=graphs+1;
            if isisomorphic(graph(G),graph(A))==false
                valid_new_graph=valid_new_graph+1;
                'New Graph';
                new_graph{valid_new_graph}=A;
            end
        end
    end
end
graphs
valid_length_Y
valid_length_Y_y
valid_graph
valid_new_graph

[Bruno Luong]

Stephan finds the bug I made,

The correct call must be:

x = intlinprog(fdummy, 1:n, [], [], Aeq, beq, lb, ub);