Permuting elements in an symmetrical matrix with restricted outcomes

I would like to derive all possible permutations of a matrix A = [0 2 0 1; 2 0 1 0; 0 1 0 1; 1 0 1 0] by permuting the positions of elements with the following restrictions:

[1] Diagonal elements can only be "0" or "2" (i.e. the permuated matrix A = [0 2 0 1; 0 1 2 0; 0 1 0 1; 1 0 1 0] would NOT be valid as "1" is a diagonal element

[2] The positions and values (i.e. 0, 1 or 2) of elements (not including diagonal elements) must be "mirrored" across the diagonal (i.e. the permuted matrix A = [0 2 0 1; 2 0 1 0; 1 0 0 1; 0 1 1 0] is NOT valid

Can someone help me with a code for this ?

Thus far I have been generating all possible permutation matrices for this 4 x 4 matrix (4!) and iteratively multiplying the matrix A by each permutation matrix but I am not sure if I can derive all possible permutations using this method (not to mention it being extremely time consuming as eventually I will be permuting many more, larger matrices)

Answer 1

[Adam Danz]

Since the permuted matrices are all symmetric, you really only need to permute the lower (or upper) triangle of the matrix, excluding the diagonal, and then reflect the values. Since your matrix is 4x4, there are 6 values in the lower triangle excluding the diagonal. That results in 6! permutations (720). But there are also 16 variations of the diagonal (4 choices from 2 values [0,2] with replacement; if that's what you're looking for since it was unlcear).

so in total, there are 720 * 16 = 11520 permutations of matrix A.

Some of the lines below could be combined but I chose to keep them separated to increase readability. There may be additional shortcuts that I didn't think of but the code completes in less than 150 milliseconds.

The code requires Jos' permn() function from the file exchange in order to compute the permutations of diagonal elements with replacement.

See the inline comments for detials.

% Input matrix
A = [0 2 0 1; 2 0 1 0; 0 1 0 1; 1 0 1 0];

% Get index of lower triangle values (below diagonal)
trilMat = tril(reshape(1:numel(A),size(A)),-1); 
trilIdx = trilMat(trilMat>0);

% List all permutation of lower triangle indices
trilPerms = perms(trilIdx); 

% List accepted diagonal values 
diagPermitted = [0,2]; 

% List all permutation of diagonal values (with replacement)
% https://www.mathworks.com/matlabcentral/fileexchange/7147-permn
diagPerms = permn(diagPermitted,numel(diag(trilMat)));

% Loop through all lower-triangle-permutation and diagonal-permutations
diagMask = logical(eye(size(A)));
nPerms = size(trilPerms,1) * size(diagPerms,1); % total number of permutations
A_PERM = nan([size(A),nPerms]);                 % all permutations
for i = 1:size(trilPerms,1)
    for j = 1:size(diagPerms,1)
        c = (i-1)*size(diagPerms,1)+j;
        base = zeros(size(A)); 
        base(trilIdx) = A(trilPerms(i,:));      % permute lower triangle
        base = base + tril(base,-1).';          % reflect lower triangle 
        base(diagMask) = diagPerms(j,:);        % permute diagonal
        if ~issymmetric(base)
            % Sanity check: matrix is symmetric
            error('Matrix not symmetic. i = %d  j = %d',i,j)
        end
        A_PERM(:,:,c) = base;                   % store matrix
    end
end

A_PERM is a [4 x 4 x 11520] array containing all 11520 permutations of the 4x4 matrix A.

[addendum]

"I would like to restrict the resultant, permuted matrices to also have rows and columns that sum to 3,3,2,2, more specifically two of the rows (or columns) in the "accepted" permuted matrices must sum to 2 and the other two must sum to 3, the order of the row (column) sums does not matter."

It would be complicated to work that restriction into the permutation rules. Instead, keep all of the code above, produce all permutations, and then eliminate any that do no follow that rule. That just requires adding the the following two lines below at the end of the code above.

goodRows = ismember(squeeze(sort(sum(A_PERM,2))).',[2 2 3 3],'rows'); 
A_PERM(:,:,~goodRows) = []; 

To see the sum of each row,

sum(A_PERM,2)

[addendum II]

Extract the unique matrices (see the discussion below for details).

% Extract unique matrices.
% Each matrix is collapsed to a row vector and all
% row vector are combined into a mega-matrix that 
% contains all of the A_PERM data.  Then we'll  
% identify unique rows numbers.  
megaMat = squeeze(reshape(A_PERM,1,prod(size(A_PERM,[1,2])),size(A_PERM,3))).'; 
[~, unqRowIdx] = unique(megaMat,'rows');
A_PERM_UNQ = A_PERM(:,:,unqRowIdx); 

[Maxwell Day]

Thanks Adam, this is great!

How would I go about adding the following output restriction:

In the input matrix A = [0 2 0 1; 2 0 1 0; 0 1 0 1; 1 0 1 0] the sum of the first two rows (or coloumns) equals 3 and the sum of the 3rd and 4th rows (or columns) equals 2. We therfore have the row sums 3,3,2,2.

I would like to restrict the resultant, permuted matrices to also have rows and columns that sum to 3,3,2,2, more specifically two of the rows (or columns) in the "accepted" permuted matrices must sum to 2 and the other two must sum to 3, the order of the row (column) sums does not matter.

I.e. the permuted matrix = [0 2 1 1; 2 0 0 0; 1 0 0 1; 1 0 1 0] would NOT be valid as the first row (or column) sums to 4.

I would also like to display permuted matrices as individual 4 x 4 matrices rather than an array, is this possible?

[Adam Danz]

Thanks Adam, this is great!

Glad I could get things started!

How would I go about adding the following output restriction:

I was afraid this was going to happen.

Just to put this into perspective, I got interested in solving this last night and ended up spending more than a couple of hours constructing that solution (trial & error, tweeking, simplifying, etc).

This is why it's critical to make sure the initial description is as precise as possible because now we've got a good solution to a problem that doesn't exist.

I'll have to think about this more but after I'm done with some other work I'm doing now.

But before any more time is invested, there's another part of your question that worries me, "eventually I will be permuting many more, larger matrices". How large is the largex matrix? The perms function will definitely hit a limit if the input gets too large.

I would also like to display permuted matrices as individual 4 x 4 matrices rather than an array, is this possible?

I'm not sure what you mean by 'display'. Do you mean you want to print all 11520 matrices to the command window? I doubt that's what you want to do.

You've got to store them somewhere. The only other viable option is within a 11520 x 1 cell array where A_PERM{n} would contain the n_th matrix. But that would be much less efficient than storing them in a 4 x 4 x 11520 3D array where the n_th matrix is accessed by A_PERM(:,:,n).

[Adam Danz]

Actually, I just realized an easy solution; see updated answer "addendum".

[Maxwell Day]

Thanks you for your reply, I think we are very close to finalizing this code.

I have included the last two lines you wrote in the addendum to the intial code and it resulted in 336 permuted versions of matrix A. All of such matrices do in fact have the row sums [3,3,2,2], this is great!

However there are two issues that remain:

[1] Although none of the "accepted" matrices have diagonal values of "1" (which is good) they only have diagonal values of "0" instead of "0" or "2" despite specifying this, which I am assuming is what you did with the following lines;

diagPermitted = [0,2];

diagPerms = permn(diagPermitted,numel(diag(trilMat)));

Initially I though I improperly integrated the "permn" function, but the code seems to run fine so I dont think thats it.

[2] Some of the "accepted" matrices are not "symmetric" with respect to the off-diagonal elements, perhaps this is my fault in not being sufficiently clear while explaining what I meant.

the following matrix was one of the 336 that was generated;

P1 = [0 1 0 2; 1 0 1 0; 0 2 0 1; 1 0 1 0]

However this matrix should NOT be accepted as the element "2" occupies the row 1/column 4 position above the diagonal (or in the upper triangle) and the row 3/column 2 position below the diagonal (or in the lower triangle). If properly reflected (flipped) across the diagonal the element "2" in the lower triangle should occupy the row 4/column 1 position, which is simply the inverse of the position of the element "2" in the upper triangle.

The valid permutation of this matrix would read;

P1 = [0 1 0 2; 1 0 1 0; 0 1 0 1; 2 0 1 0] or....

P1' = [0 1 0 1; 1 0 2 0; 0 2 0 1; 1 0 1 0]

According to the calculations I have done by hand, after we finalize these restrictions [1] and [2], the number of "accepted" permutations should be drastically reduced.

I hope my explanations are clear and I apologize for my limited MatLab vocabulary.

To answer your earlier questions; the largest matrix I would be permuting would be 18 x 18.

[Adam Danz]

@Maxwell Day, I discovered an error in my answer. I was not reflecting the matrices properly. I've made the following 2 changes to my answer.

% 1) matrix is now reflected properly

% OLD
base(diagMask) = diagPerms(j,:);        % permute diagonal
base(upperMask) = base(lowerMask);      % copy to upper triangle

% NEW
base = base + tril(base,-1).';          % reflect lower triangle  
base(diagMask) = diagPerms(j,:);        % permute diagonal (not changed, but comes after reflection)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 2) The Sanity check was operating on the wrong variable
%    and the error message is now more informative.

% OLD
if ~issymmetric(A)
     % Sanity check: matrix is symmetric
     error('Matrix not symmetic.')
end

% NEW
if ~issymmetric(base)                                 % #UPDATED
     % Sanity check: matrix is symmetric
     error('Matrix not symmetic. i = %d  j = %d',i,j)  % #UPDATED
end

Now the result is 144 matrices and the new sanity check gives us confidence that none of them are non-symmetric. The matix you identified in the comment above, P1, no longer exists.

One thing to keep in mind is that there are plenty of duplicate matrices. This is to be expected since there are duplicate values in the input matrix A and we're performing permutations across all elements of that matrix (at least the lower triangle portion of it).

There are still no matrices that have diagonals other than [0 0 0 0].

Prior to removing the unaccepted matrices at the end, we've still got all 11520 permutations that are possible under the initial rules. When we eliminate the matrices whose row-sums do not contains two 2s and two 0, the matricies that contain a 2 in the diagonal all vanish. That tells me that it's not possible to combine the initial set of rules with the new row-sum rule and still have 2s in the diagonal. As an exercise, can you create such a matrix that follows the initial rules and the row-sum rule and still has a 2 in the diagonal? If you can, then our algorithm is flawed. But my guess is that it's not possible

"the largest matrix I would be permuting would be 18 x 18."

Not good news. An 18x18 matrix would result in permuting 153 values which would result in 2.0063e+269 permutation -- way beyond the capabilites of Matlab. Even permuting 15 values will freeze your computer (don't try it: perms(1:15)).

[Maxwell Day]

@Adam Danz, I just ran the updated code, very nice!

Yes you are correct, I should have realized that matrix A had two elements "2" and therfore 6 of the other elements would have to be "1" (=10 (3+3+2+2)), to maintain symmetry both "2"s would have to be diagonal elements and therefore force two rows to have a sum of 4 which is not allowed.

I tested the new code with a new adjacency matrix that contains only a single "2" and I do see permuted matrices with the diagonal element 2, therfore I am confident in saying this code works!

One last thing:

Is there any way to eliminate the duplicate permuted matrices after the code has been executed? I permuted the matrix Ar = [2 1 0 0; 1 0 1 1; 0 1 0 1; 0 1 1 0] and got 576 matrices when there are actually only 24 unique permuted matrices (according to what I've calculated by hand). Extracting these unique matrices would save me alot of time if possible.

Thank you very much for your time, I appreciate it!

[Adam Danz]

"Is there any way to eliminate the duplicate permuted matrices after the code has been executed?"

Yes, I've updated my answer - see addendum II. However, since the original problem changed again, I'm losing confidence that the my answer (which is a solution to the original problem that has been adapted a few times to fit changes to the problem) is the most efficient approach. Think about this: With the original input matrx A, 11520 permutations are created. Then we eliminate most of them due to the row-sum rule and we end up with 144 matrics. Then we eliminate duplicates and we end up with 12 final matrices. That's like cooking 11520 meals for 12 people.

Given all of the updates to the problem, there's likely a more efficient way to manage the permutations and it would probably be functional with your 18x18 matrix (or even larger ones). But that would require starting over at the conceptual level followed by several more hours of trial-and-error implementation.

I ran your new Ar example matrix and only found 12 unique matrices so please check your calculations that suggested you'd get 24 and check the list of unique matrices to see if there's a problem (I don't expect there to be a problem with that part of the code).

[Maxwell Day]

Everything works as expected, and yes there should be 12 unique matrices rather than 24 (12 of the 24 matrices are duplicates in my calculation). Although there may be a more efficient algorithm to produce the same result this code will work for now. Thank you.

When permuting matrices that contain different element combinations or matrices of a different size (i.e. 5 x 5 or 3 x 3), I am assuming i will need to make the following changes to the code:

[1] revise input matrix A

[2] Revise row sums that are allowed; I am assuming I can do this by simply changing the values [2 2 3 3] in the following lines;

goodRows = ismember(squeeze(sort(sum(A_PERM,2))).',[2 2 3 3],'rows');

A_PERM(:,:,~goodRows) = []

Is there anything else that will need to be changed in the code when permuting larger (or smaller) matrices or matrices with different elements? Anything in the addendum II code?

Out of curiosity, is it possible to draw (generate figures) of simple, undirected, non-weighted graphs from the permuted adjacency matrices? Here row 1 = node 1, row 2 = node 2....etc. and the elements would equal the number of edges between respective nodes. Matrix elements of "1" correspond to undirected edges between two different nodes, off-diagonal elements of "2" correspond to two undirected edges between respective nodes and diagonal elements of "2" correspond to self-loops to a single node.

Thanks again.

[Adam Danz]

1. Yes, the input matrix is A
2. Yes, that will have to change, too. Perhaps you could use the diagPermitted variable to programmatically create that vector in ismember().
3. Anything in the addendum II code?: nope, that's adaptive to the size of the input matrix.

I tested the code using a 3x3 input matrix and only 3 values in the vector mentioned in point #2 above and there were no errors.

[Maxwell Day]

Great! Thank you.

I will try usine the diagPermitted variable in the ismember() line

Let me know about drawing graphs as mentioned in my previous comment.

[Adam Danz]

Check out the examples of creating graphs here.

To access matrix number n:

A_PERM_UNQ(:,:,n)

If you'd like to create a graph for each matrix, you can show me the code needed to create a graph for 1 of the matrices and I can show you how to set that up in arrayfun() that will apply the code to each matrix.

[Maxwell Day]

Thus far I have simply used:

G = graph(A_PERM_UNQ(:,:,n))

To get the edges and nodes tables

then to draw the graph I use:

plot(G) 

Which is produces a simple, undirected graph. This is all I actually need I just have to draw one graph at a time.

[Adam Danz]

Try this (requires r2019b)

% Plot it
tiledlayout('flow')
arrayfun(@(n)plot(nexttile,graph(A_PERM_UNQ(:,:,n))),1:size(A_PERM_UNQ,3))

To set axis properties and such, see

• https://www.mathworks.com/help/matlab/ref/tiledlayout.html
• https://www.mathworks.com/help/matlab/ref/nexttile.html

[Maxwell Day]

Wow! very cool, thanks again.

[Maxwell Day]

@Adam Danz very sorry to keep pestering you but I may have found one small flaw in the code.

When permuting the matrix A = [0 2 0 1; 2 0 1 0; 0 1 0 1; 1 0 1 0] I recieve 12 valid permutations. However, we expect more as the elements "2" can occupy off-diagonal AND diagonal positions and still correspond to a symmterical matrix. (we discussed this issue before but I didnt realize it was possible at the time)

The code thus far generates all valid permutations where the element "2" occupies off-diagonal elements but not the symmetrical matrices that satisfy our symmetry rule and row sum rule where both elements "2" occupy diagonal positions such as;

A2 = [2 1 0 0; 1 0 1 0; 0 1 0 1; 0 0 1 2]

I believe this is the only possible permutation where "2" occupies diagonal elements, It is an easy enough problem to work out by hand after the code has been executed but an interesting problem nonetheless.

Let me know what you think. Thanks.

[Adam Danz]

Maxwell Day, it's really good you're testing this with a fine level of detail.

My first thought was to include the diagonal in the permutations. Currently we're only permuting the lower triangle exluding the diagonal. For a 4x4 matrix, that increase the number of values being permuted from 6 to 10 which doesn't sound like much but adds 3,628,080 additional permutations (10!-6!) and it froze my laptop completely. I'm afraid I can't put too much more time into this right now. I really hate letting problems that I've dipped into go unresolved. If you find that this problem is significant and you cannot make any progress on it over the next few days, please follow up with another comment so it comes back on my radar.

By the way, here is some code that might be useful to you. The first line converts the 4x4xn array into a 1xn cell array each containing 4x4 matrix. The following lines demonstrate how to locate a matrix within the array (if one exists).

% convert to cell array
C =  squeeze(mat2cell(A_PERM,4,4,ones(size(A_PERM,3),1)));

% Find matrix A2 within the array
A2 = [2 1 0 0; 1 0 1 0; 0 1 0 1; 0 0 1 2]; 
idx = cellfun(@(m)isequal(m,A2),C); % logical index
sub = find(idx)                     % subscript

[Maxwell Day]

Thanks @Adam Danz I will use these.

Permuting the diagonal elements in addition to the lower triangle isnt a big deal, I can solve for these "extra" permutations by hand.

One other issue when ploting the graphs has emerged:

I permuted the matrix Ax = [2 0 0 1; 0 0 2 0; 0 2 0 1; 1 0 1 0] and recieved 36 unique permuted matrices, good.

However, when I attempt to plot the graphs based on these matrices using the code:

tiledlayout('flow')
arrayfun(@(n)plot(nexttile,graph(A_PERM_UNQ(:,:,n))),1:size(A_PERM_UNQ,3))

It doesnt work and I recieve the error message:

Error using graph/subsref (line 8)
Indexing not supported.

Error in @(n)plot(nexttile,graph(A_PERM_UNQ(:,:,n)))

Do you know why this is happening? Is there an easy fix here?

Thanks

[Adam Danz]

No error when I try it. A_PERM_UNQ is a 4x4x36 array it produces the figure.

Answer 2

[Maxwell Day]

Hey @Adam Danz

The code is working great so far thanks again. I have another interesting problem if you are interested.

After executing the code and generating a set of unique (valid) permuted matrices (A1, A2...An) I would like to determine which of these matrices correspond to non-isomorphic (topologically distinct) graphs (G1,G2...Gn).

Two graphs G1 (from matrix A1) and G2 (from matrix A2) are isomorphic only if there exists a permutation matrix P such that:

(P)(A1)(P^-1) = A2

Here P^-1 is simple the inverted permutation matrix P.

Doing this for certain matrices is not possible as they will have a determinant of 0 and therfore not be invertible.

Is there a way to determine which matrices in a set generate isomorphic graphs by getting MatLab to determine if a permutation matrix P exists that satisfies the above equation? If this permutation matrix P exists we can then say G1 is isomorphic with G2 and eliminate either A1 or A2 (as they are effectively the same), if P did not exist we would keep both A1 and A2.

Such a code would have to apply the above equation to each unique combination of matrices in set. I.e. if we have 12 matrices there would be 66 unique combinations of two matrices as order doesnt matter and there is no reason to pair a matrix with itself.

A1-2,-3,-4,-5,-6,-7,-8,-9,-10,-11,-12

A2-3,-4,-5,-6,-7,-8,-9,-10,-11,-12

A3-4,-5,-6,-7,-8,-9,-10,-11,-12

etc.

Depending on the size of the matrix, this may be a very large and intensive calculation, let me know if this is possible, thanks.

-Max

[Adam Danz]

Something like this?

for i = 1:size(A_PERM_UNQ,3)
    A2 = A_PERM_UNQ(:,:,i) * A1 * A_PERM_UNQ(:,:,i)^(-1);
end

[Maxwell Day]

When I run this code I get:

Unrecognized function or variable 'A1'.

Do we have to specify what "A1" is?

The permutation matrix "P" or what seems to be A_PERM_UNQ in your code should refer to 1 of n permutation matrices possible for an n x n matrix where the total number of permutation matrices equals n!

i.e. for a 4 x 4 matrix there are 24 permutation matrices (P1-P24), if there are 12 unique matrices the equation should execute as follows:

A1-1 = (P1)(A_PERM_UNQ(:,:,1)(P1^-1)

A1-2 = (P2)(A_PERM_UNQ(:,:,1)(P2^-1)

....... A1-24 = (P24)(A_PERM_UNQ(:,:,1)(P24^-1)

Where A_PERM_UNQ(:,:,1) is the first of the 12 unique matrices

Here, if any of the calculated matrices A1-1 to A1-24 are identical to any of the 12 unique matrices (other than the A_PERM_UNQ(:,:,1) matrix obviously) those two matrices are topologically identical or isomorphic.

The calculation should then proceed with applying all 24 permutation matrices to the second unique matrix (A_PERM_UNQ(:,:,2) as follows:

A2-1 = (P1)(A_PERM_UNQ(:,:,2)(P1^-1)

A2-2 = (P2)(A_PERM_UNQ(:,:,2)(P2^-1)

.......A2-24 = (P24)(A_PERM_UNQ(:,:,2)(P24^-1)

Where A_PERM_UNQ(:,:,2) is the second of the 12 unique matrices

Here, if any of the calculated matrices A2-1 to A2-24 are identical to any of the 12 unique matrices (other than the A_PERM_UNQ(:,:,2) matrix obviously) those two matrices are topologically identical or isomorphic. This is what I am hoping MatLab can tell me.

The calculation would then continue with the third unique matrix (A_PERM_UNQ(:,:,3)) calculating matrices A3-1 to A3-24 by applying all 24 permutation matrices (P1-P24).

I hope this clarifies the problem.

[Adam Danz]

A_PERM_UNQ came from my answer - it contains the the permutations of A.

I don't know what A1 is supposed to be.

Answer 3

[Maxwell Day]

Hello @Adam Danz

I seem to be running into error related to exceeding the maximum array size with a 6 x 6 matrix, I was under the impression the code we had developed could handle matrices of this size.

I tried to permute the following matrix using the following code and recieved this error:

>>
% Input matrix
A = [0 1 1 1 1 0; 1 0 1 1 1 0; 1 1 0 0 0 0; 1 1 0 0 0 1; 1 1 0 0 0 0; 0 0 0 1 0 0];
% Get index of lower triangle values (below diagonal)
trilMat = tril(reshape(1:numel(A),size(A)),-1); 
trilIdx = trilMat(trilMat>0);
% List all permutation of lower triangle indices
trilPerms = perms(trilIdx); 
% List accepted diagonal values 
diagPermitted = [0,2]; 
% List all permutation of diagonal values (with replacement)
% https://www.mathworks.com/matlabcentral/fileexchange/7147-permn
diagPerms = permn(diagPermitted,numel(diag(trilMat)));
% Loop through all lower-triangle-permutation and diagonal-permutations
diagMask = logical(eye(size(A)));
upperMask = triu(true(size(A)),1); 
lowerMask = tril(true(size(A)),-1);
nPerms = size(trilPerms,1) * size(diagPerms,1); % total number of permutations
A_PERM = nan([size(A),nPerms]);                 % all permutations
for i = 1:size(trilPerms,1)
    for j = 1:size(diagPerms,1)
        c = (i-1)*size(diagPerms,1)+j;
        base = zeros(size(A)); 
        base(trilIdx) = A(trilPerms(i,:));      % permute lower triangle
        base = base + tril(base,-1).';          % reflect lower triangle  
base(diagMask) = diagPerms(j,:);        % permute diagonal (not changed, but comes after reflection)

if ~issymmetric(base)                                 % #UPDATED
     % Sanity check: matrix is symmetric
     error('Matrix not symmetic. i = %d  j = %d',i,j)  % #UPDATED
        end
        A_PERM(:,:,c) = base;                   % store matrix
    end
end
goodRows = ismember(squeeze(sort(sum(A_PERM,2))).',[1 3 3 3 3 3],'rows'); 
A_PERM(:,:,~goodRows) = []

Error using zeros
Requested 479001600x12 (42.8GB) array exceeds maximum array size preference. Creation of arrays
greater than this limit may take a long time and cause MATLAB to become unresponsive. See array size
limit or preference panel for more information.

Error in perms>permsr (line 49)
    P = zeros(nn*m,nn);

Error in perms (line 31)
    P = permsr(n);

Is there anyway around this? let me know, Thanks.

[Adam Danz]

I don't recall claiming that it would work on 6x6 arrays. With that matrix size, you've got 15 values in the lower triangle below the diagonal. That's 1.3077e+12 permutations which, as you can see, is beyond the maximum array size set in the preferences.

factorial(numel(trilIdx))
   1.3077e+12

The development of this algorithm changed a lot over time as the rules were changed and new rules were added. In the end, the algorithm is doing a heck of a lot of work and then toss most of it away. That's inefficient. There's likely a better solution. Could you summarize all of the rules as precisely as possible in a bulleted list? That may help in inspire an alternative.

[Maxwell Day]

Hi @Adam Danz, sorry for the delayed response.

The following is what I would like to accomplish:

For a given nxn square adjacency matrix (A) I would like to permute the positions of elements in A to find all non-identical (unique) permuted versions of A that are "valid". A "valid" matrix must satisfy the following requirements:

[1] A valid matrix must be symmetrical. The upper and lower triangles of the matrix must be symmetrical when flipped across the diagonal.

[2] Diagonal elements CAN ONLY BE 0 and 2. Because matrices will only ever have elements of 0, 1 and/or 2, this means that matrix elements of 1 will only ever occupy off-diagonal positions and matrix elements of 0 and 2 can occupy any matrix position.

[3] The sums of the rows or columns of the input matrix A must be maintained in no particular order. i.e. The matrix A = [0 1 0 0; 1 2 1 0; 0 1 0 1; 0 0 1 0] has the row (and column) sums 1, 4, 2, 1. A valid permuted version of this matrix must also have the row (column) sums 1, 4, 2, 1 in any order. Therfore a matrix with the sums 2, 1, 4, 1 would be valid.

[4] Permuted matrices with "floating vertices" are NOT valid. Because we are relating all permuted matrices to graphs (where each row (or column) represents a vertex (or node) and each element value represents the number of linkages (edges) between respective nodes it is necessary to draw each permuted matrix which we have been doing using the code;

>> % Plot it

tiledlayout('flow')

arrayfun(@(n)plot(nexttile,graph(A_PERM_UNQ(:,:,n))),1:size(A_PERM_UNQ,3))

In these graphical representations of each matrix it is required that NO vertices are "floating" or have no linkages connecting them to other vertices. i.e. in the matrix A = [0 0 1 1; 0 2 0 0; 1 0 0 2; 1 0 2 0] the element "2" in the second row has no connections to other vertices and is therfore floating. This matrix would have 4 vertices as there are 4 rows (or columns).

In some cases it is not a single isolated vertex that is floating but rather a group of vertices. i.e. In the matrix A = [2 1 0 0; 1 2 0 0; 0 0 2 1; 0 0 1 0] vertices 1 and 2 are connected to eachother but not to vertices 3 and 4. Vertices 3 and 4 are connected to eachother but not to vertices 1 and 2. Once drawn, you will see this matrix results in a graph that has two "floating" pairs of vertices, this is also NOT ALLOWED and is therfore not representative of a valid permuted matrix.

[5] After all valid permuted matrices that satisfy requirements [1] - [4] have been determined, I would like to determine which of those valid matrices correspond to topologically distinct (non-isomorphic) graphs.

After executing the code and generating a set of unique (valid) permuted matrices (A1, A2...An) I would like to determine which of these matrices correspond to non-isomorphic (topologically distinct) graphs (G1,G2...Gn).

Two graphs G1 (from matrix A1) and G2 (from matrix A2) are isomorphic only if there exists a permutation matrix P such that:

(P)(A1)(P^-1) = A2

Here P^-1 is simply the inverted permutation matrix P.

Doing this for certain matrices is not possible as they will have a determinant of 0 and therfore not be invertible.

Is there a way to determine which matrices in a set generate isomorphic graphs by getting MatLab to determine if a permutation matrix P exists that satisfies the above equation? If this permutation matrix P exists we can then say G1 is isomorphic with G2 and eliminate either A1 or A2 (as they are effectively the same), if P did not exist we would keep both A1 and A2 as they are non-isomorphic (topologically different)

Such a code would have to apply the above equation to each unique combination of matrices in set. I.e. if we have 12 matrices there would be 66 unique combinations of two matrices as order doesnt matter and there is no reason to pair a matrix with itself.

A1-2,-3,-4,-5,-6,-7,-8,-9,-10,-11,-12

A2-3,-4,-5,-6,-7,-8,-9,-10,-11,-12

A3-4,-5,-6,-7,-8,-9,-10,-11,-12

etc.

I am sure this piece of code will be much more difficult to get right, perhaps it is worth writing seperately from the code involving constraints [1] - [4] ?

Please let me know what you think and let me know if anything is unclear.

-Max

[Adam Danz]

The bottleneck is the perms() function that causes memory problems when the input vector exceeds ~10 elements (or even less).

1) Instead of computing all of the permutations at once (which will not be possible), you'll need to compute the permutations in chunks and fully analyze each chunk prior to generating the next chunk. Here are two FEX functions that you can check out (I haven't used either of them).

https://www.mathworks.com/matlabcentral/fileexchange/57429-nextperm-v-k

https://www.mathworks.com/matlabcentral/fileexchange/44887-nthperm

2) Most of the algorithm we already have will still be used. However, instead of building all possible matrix combinations and then eliminating the ones that don't fulfill all of the requirements, generate the permuted matrices one at a time within a loop and run each test on that matrix. If it doesn't pass all of the tests, toss it. Most will be tossed and this could still take a while to run. But it should avoid the memory restrictions imposed by perms().

[Maxwell Day]

Hey @Adam Danz

Makes sense, I am assuming I would only have to modify the order of operations in the following code:

>>

% Input matrix
A = [0 1 1 1 1 0; 1 0 1 1 1 0; 1 1 0 0 0 0; 1 1 0 0 0 1; 1 1 0 0 0 0; 0 0 0 1 0 0];
% Get index of lower triangle values (below diagonal)
trilMat = tril(reshape(1:numel(A),size(A)),-1); 
trilIdx = trilMat(trilMat>0);
% List all permutation of lower triangle indices
trilPerms = perms(trilIdx); 
% List accepted diagonal values 
diagPermitted = [0,2]; 
% List all permutation of diagonal values (with replacement)
% https://www.mathworks.com/matlabcentral/fileexchange/7147-permn
diagPerms = permn(diagPermitted,numel(diag(trilMat)));
% Loop through all lower-triangle-permutation and diagonal-permutations
diagMask = logical(eye(size(A)));
upperMask = triu(true(size(A)),1); 
lowerMask = tril(true(size(A)),-1);
nPerms = size(trilPerms,1) * size(diagPerms,1); % total number of permutations
A_PERM = nan([size(A),nPerms]);                 % all permutations
for i = 1:size(trilPerms,1)
    for j = 1:size(diagPerms,1)
        c = (i-1)*size(diagPerms,1)+j;
        base = zeros(size(A)); 
        base(trilIdx) = A(trilPerms(i,:));      % permute lower triangle
        base = base + tril(base,-1).';          % reflect lower triangle  
base(diagMask) = diagPerms(j,:);        % permute diagonal (not changed, but comes after reflection)
if ~issymmetric(base)                                 % #UPDATED
     % Sanity check: matrix is symmetric
     error('Matrix not symmetic. i = %d  j = %d',i,j)  % #UPDATED
        end
        A_PERM(:,:,c) = base;                   % store matrix
    end
end
goodRows = ismember(squeeze(sort(sum(A_PERM,2))).',[1 3 3 3 3 3],'rows'); 
A_PERM(:,:,~goodRows) = []

I am assuming all the perms() functions need to be replaced or changed. How do we do this to permute only a "chunk" of matrices ? Then after I am assuming the code would have to cycle through the "%Sanity check: matrix is symmetric" code to ensure symmetry and then the "goodrows" code to ensure the row sums are correct? This would have to happen for each chunk, how can we set this up?

-Max

[Adam Danz]

"I am assuming all the perms() functions need to be replaced or changed. How do we do this to permute only a "chunk" of matrices ?"

The perms() function will need replaced. In my previous comment I recommended two functions to check out on the file exchange (I've never used either of them - you can decided which one is best).

After you do that, "trilPerms" will not be a full list of permutations, I will only contain a subset of them, or maybe even just one permutation at a time, depending on which function you choose and how you use it. Perhaps it would be easiest just to produce one at a time. If you do it in small chunks it might (?) be faster but you'll still need to loop over single permutations.

You can probably still use the permn function since there are only 2 possible values.

The code is currently computing all permutation and looping through them starting at

for i = 1:size(trilPerms,1)
    for j = 1:size(diagPerms,1)

But now trilPerms will only be 1 permutation or a chunk of permutations. The permuted matrix is called 'base' and the current code is saving all of them. You'll change this. Instead, you'll perform all of the tests on "base" and only if it passes all of the tests will you store it in A_PERM. Otherwise, you won't store it and you'll move on to the next permutation.

It's a good idea to save a copy of your current code so you can go back and reference it. When you're done with the new version, run a small input matrix through both versions and make sure you're getting the same results.

Lastly, it would be a good investment of your time to go through the existing code you shared above, line by line, to make sure you understand exactly what's happening. It will be much easire to redesign it if you know what it's doing.

[Maxwell Day]

Hi @Adam Danz

Thanks for the recommendations, after reading through the NEXTPERM and NTHPERM codes I've decided to try both. Starting with the nthperm code as it allows me to loop the tests through a single permuted matrix whereas it seems like nextperm code can only do this in "blocks".

After attemping to revise the code above over and over it seems my lack notation knowledge is causing problems. The problems I think I am facing now is as follows;

How/what do I modify to stop the code from saving all permutations of the lower triangle after replacing perms() with nthperm code?

What ordering is required to loop the matrix tests through each matrix, one at a time?

I appriciate the time you have put into this thus far, if you can't afford to put in anymore time I understand.

Thanks again!

-Max

[Maxwell Day]

Hey @Adam Danz

Attempting to use the following code revised with the nthperm code you recommened as follows

>> % Input matrix
A = [0 1 1 1 1 0; 1 0 1 1 1 0; 1 1 0 0 0 0; 1 1 0 0 0 1; 1 1 0 0 0 0; 0 0 0 1 0 0];
% Get index of lower triangle values (below diagonal)
trilMat = tril(reshape(1:numel(A),size(A)),-1); 
trilIdx = trilMat(trilMat>0);
% List all permutation of lower triangle indices
trilPerms = nthperm(trilIdx 1); 
% List accepted diagonal values 
diagPermitted = [0,2]; 
% List all permutation of diagonal values (with replacement)
% https://www.mathworks.com/matlabcentral/fileexchange/7147-permn
diagPerms = permn(diagPermitted,numel(diag(trilMat)));
% Loop through all lower-triangle-permutation and diagonal-permutations
diagMask = logical(eye(size(A)));
upperMask = triu(true(size(A)),1); 
lowerMask = tril(true(size(A)),-1);
nPerms = size(trilPerms,1) * size(diagPerms,1); % total number of permutations
A_PERM = nan([size(A),nPerms]);                 % all permutations
for i = 1:size(trilPerms,1)
    for j = 1:size(diagPerms,1)
        c = (i-1)*size(diagPerms,1)+j;
        base = zeros(size(A)); 
        base(trilIdx) = A(trilPerms(i,:));      % permute lower triangle
        base = base + tril(base,-1).';          % reflect lower triangle  
base(diagMask) = diagPerms(j,:);        % permute diagonal (not changed, but comes after reflection)
if ~issymmetric(base)                                 % #UPDATED
     % Sanity check: matrix is symmetric
     error('Matrix not symmetic. i = %d  j = %d',i,j)  % #UPDATED
        end
        A_PERM(:,:,c) = base;                   % store matrix
    end
end
goodRows = ismember(squeeze(sort(sum(A_PERM,2))).',[1 3 3 3 3 3],'rows'); 
A_PERM(:,:,~goodRows) = []

getting this error code when I attempt to run it:

trilPerms = nthperm(trilIdx 1); 
                             ↑
Error: Invalid expression. Check for missing multiplication operator, missing or unbalanced
delimiters, or other syntax error. To construct matrices, use brackets instead of parentheses.

Let me know what you think. thanks!

[Adam Danz]

Looks like you're missing the comma

trilPerms = nthperm(trilIdx, 1)
                           ^

[Maxwell Day]

Ah, I see, still getting empty arrays for permutations I know should work, assuming this is due to incorrect revisions of the trilIdx code elsewhere

>> % Input matrix
A = [2 1 0 0; 1 0 1 1; 0 1 0 1; 0 1 1 0];
% Get index of lower triangle values (below diagonal)
trilMat = tril(reshape(1:numel(A),size(A)),-1); 
trilIdx = trilMat(trilMat>0);
% List all permutation of lower triangle indices
trilPerms = nthperm(trilIdx, 1); 
% List accepted diagonal values 
diagPermitted = [0,2]; 
% List all permutation of diagonal values (with replacement)
% https://www.mathworks.com/matlabcentral/fileexchange/7147-permn
diagPerms = permn(diagPermitted,numel(diag(trilMat)));
% Loop through all lower-triangle-permutation and diagonal-permutations
diagMask = logical(eye(size(A)));
upperMask = triu(true(size(A)),1); 
lowerMask = tril(true(size(A)),-1);
nPerms = size(trilPerms,1) * size(diagPerms,1); % total number of permutations
A_PERM = nan([size(A),nPerms]);                 % all permutations
for i = 1:size(trilPerms,1)
    for j = 1:size(diagPerms,1)
        c = (i-1)*size(diagPerms,1)+j;
        base = zeros(size(A)); 
        base(trilIdx) = A(trilPerms(i,:));      % permute lower triangle
        base = base + tril(base,-1).';          % reflect lower triangle  
base(diagMask) = diagPerms(j,:);        % permute diagonal (not changed, but comes after reflection)
if ~issymmetric(base)                                 % #UPDATED
     % Sanity check: matrix is symmetric
     error('Matrix not symmetic. i = %d  j = %d',i,j)  % #UPDATED
        end
        A_PERM(:,:,c) = base;                   % store matrix
    end
end
goodRows = ismember(squeeze(sort(sum(A_PERM,2))).',[2 2 3 3],'rows'); 
A_PERM(:,:,~goodRows) = []

A_PERM =

  4×4×0 empty double array

[Adam Danz]

goodRows = ismember(squeeze(sort(sum(A_PERM,2))).',[2 2 3 3],'rows');

goodRows (a column vector) is all False values. In other words, the [2 2 3 3] vector is not found anywhere in the permuration sums.

squeeze(sort(sum(A_PERM,2))).'