# How to rearrange a Matrix elements such that all diagonal elements are non zero?

I'm trying to rearrange a matrix such that all the diagonal elements are non zero.I'm not sure how to proceed. Any suggestion are highly appreciated.

A={
{a, 0, 0, 0, 0},
{0, 0, 0, d1, 0},
{0, 0, 0, 0, e1},
{0, b, 0, d2, 0},
{0, 0, c, 0, e2}
};
Diagonal[A]


Results in

{a, 0, 0, d2, e2}


I used a solution given by "corey979" for the question "Writing a program to swap rows and columns [duplicate]". To swap any two rows i and j of a general m×n matrix:

  swapRij[mat_, {i_, j_}] := Block[{mat1 = mat},
mat1[[{i, j}]] = mat[[{j, i}]];
mat1 ]


If i'm using this function then I need to manually check the position and then send those row numbers as an input to that function.

It would be very difficult to arrange a big matrix.I have a matrix of size 30 x 30. I wanted to automate it.

Output i'm expecting is

   A={
{a, 0, 0, 0, 0},
{0, b, 0, d2, 0},
{0, 0, c, 0, e2},
{0, 0, 0, d1, 0},
{0, 0, 0, 0, e1}
}


Where the diagonal elements are

{a, b, c, d1, e1}


I'm adding my actual matrix as well.

    {{m1m1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, {0, m1m1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0},
{0, 0, m1m1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0}, {0, 0, 0, m4m4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, m4m20, m4m21, 1, 0, 0, 0, m4m26, m4m27, 1, 0, 0},
{0, 0, 0, 0, m5m5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, m5m19, 0, m5m21, 0, 1, 0, m5m25, 0, m5m27, 0, 1, 0}, {0, 0, 0, 0, 0, m6m6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, m6m19, m6m20, 0, 0, 0, 1, m6m25, m6m26, 0, 0, 0, 1},
{0, 0, 0, 0, 0, 0, m7m7, m7m8, m7m9, m7m10, m7m11, m7m12, 0, 0, 0, 0, 0, 0, m7m19, m7m20, m7m21, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, m8m7, m8m8, m8m9, m8m10, m8m11, m8m12, 0, 0, 0, 0, 0, 0, m8m19, m8m20, m8m21, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, m9m7, m9m8, m9m9, m9m10, m9m11, m9m12, 0, 0, 0, 0, 0, 0, m9m19, m9m20, m9m21, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, m10m13, m10m14, m10m15, m10m16, m10m17, m10m18, 0, 0, 0, 0, 0, 0, m10m25, m10m26, m10m27,
0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, m11m13, m11m14, m11m15, m11m16, m11m17, m11m18, 0, 0, 0, 0, 0, 0, m11m25, m11m26, m11m27, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, m12m13, m12m14, m12m15, m12m16, m12m17, m12m18, 0, 0, 0, 0, 0, 0,
m12m25, m12m26, m12m27, 0, 0, 0}, {0, 0, 0, 0, 0, 0, m13m7, m13m8, m13m9, m13m10, m13m11, m13m12, 0, 0, 0, 0, 0, 0, m13m19, m13m20, m13m21, m7m19, m7m20, m7m21, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, m14m7, m14m8, m14m9, m14m10, m14m11, m14m12, 0, 0, 0, 0, 0, 0, m14m19, m14m20, m14m21, m8m19, m8m20, m8m21, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, m15m7, m15m8, m15m9, m15m10, m15m11, m15m12, 0, 0, 0, 0, 0, 0, m15m19, m15m20, m15m21,
m9m19, m9m20, m9m21, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, m16m13, m16m14, m16m15, m16m16, m16m17, m16m18, 0, 0, 0, 0, 0, 0, m16m25, m16m26, m16m27, m10m25, m10m26, m10m27},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, m17m13, m17m14, m17m15, m17m16, m17m17, m17m18, 0, 0, 0, 0, 0, 0, m17m25, m17m26, m17m27, m11m25, m11m26, m11m27}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, m18m13, m18m14, m18m15, m18m16, m18m17, m18m18, 0, 0, 0, 0, 0,
0, m18m25, m18m26, m18m27, m12m25, m12m26, m12m27}, {0, 0, 0, m7m19, m7m20, m7m21, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, m8m19, m8m20, m8m21, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, m9m19, m9m20, m9m21, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, m10m25, m10m26, m10m27, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, m11m25, m11m26, m11m27, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, m12m25, m12m26, m12m27, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 0, 0, 0, m5m19, m6m19, m7m19, m8m19, m9m19, m13m19, m14m19, m15m19, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 1, 0, m4m20, 0, m6m20, m7m20, m8m20, m9m20, m13m20, m14m20, m26m12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, m4m21, m5m21, 0, m7m21, m8m21, m9m21, m13m21, m14m21, m15m21, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{1, 0, 0, 0, m5m25, m6m25, 0, 0, 0, 0, 0, 0, m10m25, m11m25, m12m25, m16m25, m17m25, m18m25, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 0, m4m26, 0, m6m26, 0, 0, 0, 0, 0, 0, m10m26, m11m26, m12m26, m16m26, m17m26, m29m18, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0}, {0, 0, 1, m4m27, m5m27, 0, 0, 0, 0, 0, 0, 0, m10m27, m11m27, m12m27, m16m27, m17m27, m18m27, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}}


Thank you,

• Please remove MatrixForm[] or write (A= ....)//MatrixForm to work with your matrix. MatrixForm[] is only for display. Commented May 2, 2020 at 6:44
• And please link to the question that you reference so that we can easily look it up (but in this case, Michael has already said what the problem is.) Commented May 2, 2020 at 6:46
• Please. Could you tell us the reason why the diagonalization is necessary? Commented May 2, 2020 at 15:51
• @Cesareo 1. This will help me to check my main equations from which this matrix is formed. 2.I was going to perform Gauss elimination method to this symbolic matrix. Commented May 3, 2020 at 1:20

The question can be rephrased as:

Can we assign an index to each row such that 1) the row contains a nonzero entry at that index and 2) no other row is assigned that index?

Each row has a set of what we might call relevant indices, that is, the set of indices at which it is nonzero. (E.g. the list {1,0,0,1,0} has relevant indices {1,4}.) Armed with our sets of relevant indices, we can again rephrase the question as:

Can we choose an element (index) from each set such that we do not choose the same element from any other set?

Turns out: this is a bad problem to brute-force. We find that your matrix leads to 31 990 087 039 589 244 179 841 024 possibilities to check if we try to use, say, Outer! When I tried to run code based on Outer for this matrix, it crashed my kernel.

So let's not use Outer or otherwise check each tuple. Instead, let's rephrase the problem as a graph problem, which Mathematica has quite the toolkit for. We want to draw an edge from each given set of relevant indices to each of its members. (This is a bipartite graph, with sets of relevant indices (arising from rows) on one side and indices (corresponding to columns) on another.)

Our problem is now:

Can we choose a set of edges in this graph such that every vertex belongs to some edge, and no two edges share a vertex?

This is known as a matching, or independent edge set, which Mathematica has a built-in function for: FindIndependentEdgeSet.

We can clean up our representation of things a bit. Instead of using sets of relevant indices and their elements, since sets of relevant indices each correspond to a given row, we can simply ask for a vertex representing row i to be connected to a vertex representing column j (one of the indices in the set) whenever M[[i,j]] is nonzero. We can do this by

A = Position[M, Except[0], {2}, Heads->False]


We can then turn an index pairs into an edge by

edge[{a_, b_}] := DirectedEdge[{"row", a}, {"column", b}]


(It's important we distinguish the indices representing rows from those representing columns somehow, or Mathematica will think they're the same vertex.)

And we can get the whole graph by

G = edge /@ A


Then matching = FindIndependentEdgeSet[G] finds a maximum edge set.

But is it enough? Can we, in fact, reorder the rows or columns to get a diagonal with no zeros for this matrix?

Yes! Length[matching] is 30, and M is a 30 x 30 matrix. So we have a corresponding row for each column, and vice versa.

Now to apply it: we can switch rows or columns, whichever we like. The point is that now for each row, we have in matching an edge that lands on a unique column such that that row and column in M has a nonzero entry, and we want to simply move that column to the diagonal...

or we can read it vice versa, which is easier: for each column j, matching gives a unique row i such that the entry of M at that row and column is nonzero, and we want to move that row so that it appears at index j. In other words, we want to re-order our rows so that row i of M is row j of our new matrix (call it M1). So we want a list L that has i in the jth place; then we can say M1 = M[[L]] (which works just like the swapping function you have above, but all at once).

To make this list we can call upon SortBy to reorder matching and put all the edges in the right place, and then extract their row values.

To extract column values, we can apply the rule DirectedEdge[_,{"column",b_}] :> b, and to extract row values, DirectedEdge[{"row",a_},_] :> a.

To sort by the column values and then extract row values:

reindexlist[matching_] :=
( SortBy[matching, Replace[#,DirectedEdge[_,{"column",b_}] :> b]& ]
/. DirectedEdge[{"row",a_},_] :> a )


Then

M1 = M[[reindexlist[matching]]]


And indeed, Diagonal[M1] yields

{1, m1m1, 1, m8m19, m5m5, m10m27, m7m20, m8m8, m7m9, 1, m14m21, 1,
m10m26, m11m25, m11m15, m12m16, 1, 1, m9m19, m6m20, m13m21, 1, m9m20,
m8m21, 1, m10m26, 1, m10m25, m12m26, m11m27}


So we're done! :)

For convenience, here's the code in a block, except for your matrix:

M = (*your matrix here*);

A = Position[M, Except[0], {2}, Heads -> False];

edge[{a_, b_}] := DirectedEdge[{"row", a}, {"column", b}];

G = edge /@ A;

matching = FindIndependentEdgeSet[G];

reindexlist[matching_] :=
(SortBy[matching, Replace[#, DirectedEdge[_, {"column", b_}] :> b] &]
/. DirectedEdge[{"row", a_}, _] :> a);

M1 = M[[reindexlist[matching]]]


Update: Inspired by @thorimur's excellent answer, an alternative way to get a maximal matching:

ClearAll[swapRows]
swapRows = #[[SparseArrayMaximalBipartiteMatching[# /.
Except[List, _Symbol] -> 1][[All, 1]]]] &;


Examples:

MatrixForm @ swapRows[A]


With mat as the 30X30 matrix in OP:

Diagonal @ swapRows[mat]

 {1, 1, 1, m4m21, m5m19, m6m26, m7m7, m8m8, m9m9, 1, 1, 1, m10m13,
m11m14, m12m15, 1, 1, 1, m5m19, m4m20, m13m21, m8m19, m9m20, 1, 1, 1,
1, m10m25, m12m26, m11m27}

MatrixForm[MapIndexed[If[Equal @@ #2, Style[#, Red], #] &, swapRows[mat], {2}]]


SeedRandom[1]
rm = RandomChoice[{5, 1, 1, 1} -> {0, a, b, c}, {10, 10}];

Row[MatrixForm /@ {rm,
ReplacePart[#, {i_, i_} :> Style[#[[i, i]], Red, Bold]] & @ swapRows[rm]}]


SeedRandom[1]
im = RandomSample[IdentityMatrix[15]];

Row[MatrixForm /@ {im,
ReplacePart[#, {i_, i_} :> Style[#[[i, i]], Red, Bold]] & @ swapRows[im]}]


The following methods work for the input matrix A in OP:

A = {{a, 0, 0, 0, 0}, {0, 0, 0, d1, 0}, {0, 0, 0, 0, e1}, {0, b, 0, d2, 0},
{0, 0, c, 0, e2}};

B = SortBy[LengthWhile[#, # == 0 &] &] @ A;

MatrixForm @ B


Also

B2 = SortBy[-FromDigits[# /. Except[List, _Symbol] -> 1] &] @ A;


and

B3 = SortBy[FirstPosition[#, Except[0], Heads -> False] &] @ A;

B == B2 == B3

True

• hi, thanks for the response. Would you kindly check the matrix which i have provided. i'm Sorry that I edited the information after your response.It would be very kind of you ,to provide a solution to above edited code. Commented May 2, 2020 at 7:23

EDIT: The method below only works for upper triangular matrices, I need to modify it to work for all matrices!

One quick way is as follows:

ClearAll[getDiagonal];
getDiagonal[list_List] := getDiagonal[list, {}];
getDiagonal[{}, a_] := a;
getDiagonal[list_List, results_List] := With[{
chosen = DeleteCases[list, {0, ___} | {}]},
Which[Flatten[list] == {},
results,
chosen == {},
getDiagonal[(#[[2 ;;]] & /@ list), Join[results, {0}]],
True,
getDiagonal[(#[[2 ;;]] & /@ DeleteCases[list, chosen[[1]]]),
Join[results, {First[chosen[[1]]]}]]
]
];


It is a recursive approach:

1. Given a matrix M, start with an empty list (call L)
2. Look for a row that starts with a nonzero element in M (call that row R)
3. Add the first element of R to L (add 0 to L if no R is found), call new L L'
4. Remove R from M and call resultant matrix M' (M'=M if no R is found)
5. Remove first column of M' (call resultant matrix M'')
6. If M'' does not have any column left, return L'. Otherwise go to step 1 with M'' as M and L' as L

For OP's matrices:

A={{a,0,0,0,0},{0,0,0,d1,0},{0,0,0,0,e1},{0,b,0,d2,0},{0,0,c,0,e2}};
B={{m1m1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0},{0,m1m1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0},{0,0,m1m1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0},{0,0,0,m4m4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,m4m20,m4m21,1,0,0,0,m4m26,m4m27,1,0,0},{0,0,0,0,m5m5,0,0,0,0,0,0,0,0,0,0,0,0,0,m5m19,0,m5m21,0,1,0,m5m25,0,m5m27,0,1,0},{0,0,0,0,0,m6m6,0,0,0,0,0,0,0,0,0,0,0,0,m6m19,m6m20,0,0,0,1,m6m25,m6m26,0,0,0,1},{0,0,0,0,0,0,m7m7,m7m8,m7m9,m7m10,m7m11,m7m12,0,0,0,0,0,0,m7m19,m7m20,m7m21,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,m8m7,m8m8,m8m9,m8m10,m8m11,m8m12,0,0,0,0,0,0,m8m19,m8m20,m8m21,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,m9m7,m9m8,m9m9,m9m10,m9m11,m9m12,0,0,0,0,0,0,m9m19,m9m20,m9m21,0,0,0,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,m10m13,m10m14,m10m15,m10m16,m10m17,m10m18,0,0,0,0,0,0,m10m25,m10m26,m10m27,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,m11m13,m11m14,m11m15,m11m16,m11m17,m11m18,0,0,0,0,0,0,m11m25,m11m26,m11m27,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,m12m13,m12m14,m12m15,m12m16,m12m17,m12m18,0,0,0,0,0,0,m12m25,m12m26,m12m27,0,0,0},{0,0,0,0,0,0,m13m7,m13m8,m13m9,m13m10,m13m11,m13m12,0,0,0,0,0,0,m13m19,m13m20,m13m21,m7m19,m7m20,m7m21,0,0,0,0,0,0},{0,0,0,0,0,0,m14m7,m14m8,m14m9,m14m10,m14m11,m14m12,0,0,0,0,0,0,m14m19,m14m20,m14m21,m8m19,m8m20,m8m21,0,0,0,0,0,0},{0,0,0,0,0,0,m15m7,m15m8,m15m9,m15m10,m15m11,m15m12,0,0,0,0,0,0,m15m19,m15m20,m15m21,m9m19,m9m20,m9m21,0,0,0,0,0,0},{0,0,0,0,0,0,0,0,0,0,0,0,m16m13,m16m14,m16m15,m16m16,m16m17,m16m18,0,0,0,0,0,0,m16m25,m16m26,m16m27,m10m25,m10m26,m10m27},{0,0,0,0,0,0,0,0,0,0,0,0,m17m13,m17m14,m17m15,m17m16,m17m17,m17m18,0,0,0,0,0,0,m17m25,m17m26,m17m27,m11m25,m11m26,m11m27},{0,0,0,0,0,0,0,0,0,0,0,0,m18m13,m18m14,m18m15,m18m16,m18m17,m18m18,0,0,0,0,0,0,m18m25,m18m26,m18m27,m12m25,m12m26,m12m27},{0,0,0,m7m19,m7m20,m7m21,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,m8m19,m8m20,m8m21,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,m9m19,m9m20,m9m21,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,m10m25,m10m26,m10m27,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,m11m25,m11m26,m11m27,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,0,m12m25,m12m26,m12m27,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0},{1,0,0,0,m5m19,m6m19,m7m19,m8m19,m9m19,m13m19,m14m19,m15m19,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,1,0,m4m20,0,m6m20,m7m20,m8m20,m9m20,m13m20,m14m20,m26m12,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,1,m4m21,m5m21,0,m7m21,m8m21,m9m21,m13m21,m14m21,m15m21,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{1,0,0,0,m5m25,m6m25,0,0,0,0,0,0,m10m25,m11m25,m12m25,m16m25,m17m25,m18m25,0,0,0,0,0,0,0,0,0,0,0,0},{0,1,0,m4m26,0,m6m26,0,0,0,0,0,0,m10m26,m11m26,m12m26,m16m26,m17m26,m29m18,0,0,0,0,0,0,0,0,0,0,0,0},{0,0,1,m4m27,m5m27,0,0,0,0,0,0,0,m10m27,m11m27,m12m27,m16m27,m17m27,m18m27,0,0,0,0,0,0,0,0,0,0,0,0}};


It yields following results:

getDiagonal[A]


{a, b, c, d1, e1}

getDiagonal[B]


{m1m1, m1m1, m1m1, m4m4, m5m5, m6m6, m7m7, m8m8, m9m9, m13m10, m14m11, m15m12, m10m13, m11m14, m12m15, m16m16, m17m17, m18m18, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

Performance of the code seems to be fine:

{RepeatedTiming[getDiagonal[A]][[1]], RepeatedTiming[getDiagonal[B]][[1]]}


{0.000065, 0.0016}

A "brute force" solution with Genetic Algorithms

Given a symbolic matrix, first we convert to a zero-one's matrix in which the ones represent non null elements. This is done as follows. Given M we obtain M0

{n, n} = Dimensions[M]
M0 = Table[If[NumericQ[M[[i, j]]] && M[[i, j]] == 0, 0, 1], {i, 1, n}, {j, 1, n}]


After that the fitness is calculated as the diagonal sum for the resulting transformed matrix, after a change in rows followed for a change in columns. This can be observed in the module fitnessFunction. The crossover operation is implemented as a single point crossover as can be observed in the module doSingleCrossover. The script can be optimized but it was left as it is as a means to show easily the GA procedures.

Clear[recover]
recover[M0_, bestIndividual_] := Module[{Mopt = {}, Mopt0, i},
For[i = 1, i <= length, i++, AppendTo[Mopt, M0[[bestIndividual[[1, i]]]]]];
Mopt0 = Transpose[Mopt];
Mopt = {};
For[i = 1, i <= length, i++, AppendTo[Mopt, M0[[bestIndividual[[2, i]]]]]];
Return[Mopt]
]

Clear[doMutation];
doMutation[{stringh_, stringv_}] := Module[{tempstring, i, ind1, ind2, atom, choice},
choice = RandomInteger[1];
If[choice == 1, tempstring = stringh, tempstring = stringv];
If[Random[] < mutationRate, ind1 = RandomInteger[{1, length}];
ind2 = RandomInteger[{1, length}];
atom = tempstring[[ind1]];
tempstring[[ind1]] = tempstring[[ind2]];
tempstring[[ind2]] = atom];
If[choice == 1, Return[{tempstring, stringv}], Return[{stringh, tempstring}]]
]

Clear[fitnessFunction];
fitnessFunction[{listh_, listv_}] := Module[{n = Length[M0], Mdum = {}, i, j, sum = 0, Mdum0, rowi},
For[i = 1, i <= n, i++, rowi = M0[[listh[[i]]]];
AppendTo[Mdum, rowi]];
Mdum0 = Transpose[Mdum];
Mdum = {};
For[i = 1, i <= n, i++, rowi = M0[[listv[[i]]]];
AppendTo[Mdum, rowi]];
Return[Total[Diagonal[Mdum]]]
]

Clear[doSingleCrossover];
doSingleCrossover[{stringh1_, stringv1_}, {stringh2_, stringv2_}] :=
Module[{cuth, cutv, temph1, temph2, tempv1, tempv2},
cuth = RandomInteger[{1, length}]; cutv = RandomInteger[{1, length}];
temph1 = Join[Take[stringh1, cuth], Drop[stringh2, cuth]];
temph2 = Join[Take[stringh2, cuth], Drop[stringh1, cuth]];
tempv1 = Join[Take[stringv1, cutv], Drop[stringv2, cutv]];
tempv2 = Join[Take[stringv2, cutv], Drop[stringv1, cutv]];
Return[{{temph1, tempv1}, {temph2, tempv2}}]
]

Clear[doCumSumOfFitness];
doCumSumOfFitness := Module[{temp}, temp = 0.0;Table[temp += popFitness[[i]], {i, popSize}]]

Clear[doSingleSelection];
doSingleSelection := Module[{rfitness, ind},
rfitness = RandomReal[{0, cumFitness[[popSize]]}];
ind = 1;
While[rfitness > cumFitness[[ind]], ind++];
Return[ind]
]

Clear[selectPair];
selectPair := Module[{ind1, ind2}, ind1 = doSingleSelection;
While[(ind2 = doSingleSelection) == ind1];
{ind1, ind2}
]

Clear[pickRandomPair];
pickRandomPair := Module[{ind1, ind2}, ind1 = RandomInteger[{1, popSize}];
While[(ind2 = RandomInteger[{1, popSize}]) == ind1];
{ind1, ind2}
]

Clear[exchangeString];
exchangeString[ind_, newstring_, newF_] := Module[{}, popStrings[[ind]] = newstring;
popFitness[[ind]] = newF
]

Clear[renormalizeFitness];
renormalizeFitness[fitness0_List] :=
Module[{minF, maxF, a, b, fitness = fitness0, i}, minF = Min[fitness];
maxF = Max[fitness];
a = 0.5*maxF/(maxF + minF);
b = (1 - a)*maxF;
Map[a # + b &, fitness]
]

Clear[bestDet]
bestDet := Module[{bestFitness = -1, i, ibest = 1},
For[i = 1, i <= popSize, i++,
If[popFitness[[i]] > bestFitness, bestFitness = popFitness[[i]];
ibest = i]];
If[bestFitness > bestOfAll, bestOfAll = bestFitness;
bestIndividual = popStrings[[ibest]]];
Return[popStrings[[ibest]]]
]

Clear[doInitialize];
doInitialize := Module[{i},
popFitness = Table[fitnessFunction[popStrings[[i]]], {i, popSize}];
popFitness = renormalizeFitness[popFitness];
cumFitness = doCumSumOfFitness;
listOfCumFitness = {cumFitness[[popSize]]};
historyOfPop = {bestDet}
]

Clear[updateGenerationSync];
updateGenerationSync := Module[{parentsid, children, ip}, parentsid = {};
Do[AppendTo[parentsid, selectPair], {popSize/2}];
children = {};
Do[AppendTo[children,
doSingleCrossover[popStrings[[parentsid[[ip, 1]]]],
popStrings[[parentsid[[ip, 2]]]]]], {ip, popSize/2}];
popStrings = Flatten[children, 1];
popStrings = Map[doMutation, popStrings];
popFitness = Map[fitnessFunction, popStrings];
popFitness = renormalizeFitness[popFitness];
cumFitness = doCumSumOfFitness
]


and now the main program

SeedRandom[4];
bestOfAll = -1;
popSize = 600;(*should be even*)
numberOfEpochs = 300;
mutationRate = 0.007;
n = Length[M0];
length = n;
popStrings = Table[{RandomSample[Table[i, {i, 1, n}]],
RandomSample[Table[i, {i, 1, n}]]}, {popSize}];
doInitialize;

Do[updateGenerationSync;
AppendTo[historyOfPop, bestDet];
AppendTo[listOfCumFitness,
cumFitness[[popSize]]], {numberOfEpochs}
];

ListLinePlot[Map[fitnessFunction, historyOfPop], PlotRange -> All]
bestIndividual
fitnessFunction[bestIndividual]

recover[M, bestIndividual] // Diagonal

(* {1, m1m1, 1, m8m19, m5m21, m6m20, m15m7, m9m8, m9m19, m14m10, m14m21, m15m21, m18m13, m17m14, m12m26, m10m16, m10m17, 1, m8m19, m4m20, m14m21, m8m19, m9m20, m7m21, m16m25, m16m26, m4m27, 1, m10m26, m10m27}*)


NOTE

This matrix has null determinant. Follows the fitness evolution plot, and the best individual.

(* {{25, 22, 27, 15, 12, 12, 16, 18, 5, 8, 27, 7, 19, 22, 19, 27, 22, 14, 30, 26, 18, 19, 12, 10, 25, 13, 29, 28, 7, 6}, {28, 2, 27, 20, 27, 26, 15, 9, 25, 14, 27, 27, 18, 17, 29, 10, 10, 24, 8, 4, 14, 14, 15, 13, 16, 16, 4, 4, 16, 16}}} *)

• Note that null determinant does not imply that one of the diagonal elements must be 0 under row/column reordering! The determinant is only the product of the diagonal elements when the matrix is diagonal; there's plenty of opportunity for cancellation from the other elements otherwise. (Consider the matrix {{1,1},{1,1}}`, which clearly has determinant 0, but has all nonzero diagonal elements.) Commented May 4, 2020 at 6:55