# Maximum of determinant of $n\times n$-matrix with elements 0 and 1

Generalizing https://math.stackexchange.com/questions/3265627/largest-value-of-a-third-order-determinant-whose-elements-are-0-or-1 I'd like to propose two related problems

(a) find the maximum value of the determinant of an $$n\times n$$-matrix whose elements are $$0$$ or $$1$$ (christianed a 0-1-matrix here).

(b) find all possible values of the determinants of a 0-1-matrix of order $$n$$.

I found the solutions for small $$n= 2,...,5$$ by "brute force" in Mathematica by just listing the values of the determinants all possible $$2^{n^2}$$ matrices.

On my PC I couldn't go beyond $$n=5$$ due to memory shortage.

I have two questions

1) Is somebody aware of an analytic solution?

2) Can you improve the Mathematica-code which I shall show soon in a self answer? (This part is intentionally delayed for letting the reader find own solutions).

• – AccidentalFourierTransform Jun 22 '19 at 12:29
• @ AccidentalFourierTransform Thank you for the hint to a standard problem. I shall look into it, particularly with respect to problem (b). – Dr. Wolfgang Hintze Jun 22 '19 at 21:30

Introduction

A comment to this OP pointed out (https://oeis.org/A003432) that this problem is a classic. In fact, the French mathematician Jacques Hadamard proposed an equivalent problem in 1893.

Many results have been obtained in the meantime, but even to this date the problem is not yet solved completely. (see results)

It might still be interesting to attempt a naive approach using Mathematica which can raise interesting questions and can hopefully be optimized using feedback from the readers (see Mathematica code).

Mathematica code

My code for an exhaustive enumeration of all determinants of a (0,1)-matrix of dimension $$n$$ is (example $$n=5$$, the case $$n=6$$ leads to memory overflow on my PC)

AbsoluteTiming[
With[{n = 5},
Clear[a];
m = Table[a[i, j], {i, 1, n}, {j, 1, n}]; (* the matrix *)
d = Det[m];(* the determinant *)
s = Flatten[Table[{a[i, j], 0, 1}, {i, 1, n}, {j, 1, n}], 1];(*
the index operations *)
t = Flatten[Table[d, Evaluate[Sequence @@ s]], n^2 - 1];(*
the list of all determinants *)
]
]

(* Out {122.248, Null} *)


The detailed spectrum (i.e. not just the values but also their frequency) of values of the determinants (problem (b)) is given by

Tally[t]

*( Out[100]= {{0, 21040112}, {1, 4851360}, {-1, 4851360}, {-2,
1213920}, {2, 1213920}, {3, 144720}, {-3, 144720}, {-4, 43560}, {4,
43560}, {5, 3600}, {-5, 3600}} *)


Graph of the detailed spectrum

Hence the detailed spectrum looks close to exponential.

Maximum value

mx = Max[t]

(* Out[102]= 5*)


Find matrices correponding to a given determinant (example mx)

Position of mx

p = Position[t, mx] // Flatten;
Length[p]

(* Out[104]= 3600 *)


Obtain matrix as a flat list from position

mp = PadLeft[#, 5^2] & /@ (IntegerDigits[# - 1, 2] & /@ p);

Length[mp]

(* Out[107]= 3600 *)


Partition list as matrix

mm = Partition[#, 5] & /@ mp;


and show them

mf = MatrixForm /@ mm;


Examples

Take[mf[[1 ;; 5]]]


$$\left( \begin{array}{ccccc} 0 & 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & 0 & 1 \\ 0 & 1 & 1 & 1 & 0 \\ 1 & 0 & 1 & 1 & 0 \\ 1 & 1 & 0 & 0 & 1 \\ \end{array} \right),\left( \begin{array}{ccccc} 0 & 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & 0 & 1 \\ 1 & 0 & 1 & 1 & 0 \\ 1 & 1 & 0 & 0 & 1 \\ 0 & 1 & 1 & 1 & 0 \\ \end{array} \right),\left( \begin{array}{ccccc} 0 & 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & 0 & 1 \\ 1 & 1 & 0 & 0 & 1 \\ 0 & 1 & 1 & 1 & 0 \\ 1 & 0 & 1 & 1 & 0 \\ \end{array} \right),\left( \begin{array}{ccccc} 0 & 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & 1 & 0 \\ 0 & 1 & 1 & 0 & 1 \\ 1 & 1 & 0 & 1 & 0 \\ 1 & 0 & 1 & 0 & 1 \\ \end{array} \right),\left( \begin{array}{ccccc} 0 & 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & 1 & 0 \\ 1 & 0 & 1 & 0 & 1 \\ 0 & 1 & 1 & 0 & 1 \\ 1 & 1 & 0 & 1 & 0 \\ \end{array} \right)$$

There have beeen proposals for reducing memory space by using symmetry porperties (see other solutions).

Monte-Carlo calculations

We can replace the systematic exhaustive procedure by using random values for the $$a(i,j)$$ and set up the following simple code

detRnd[n_, nn_] := Module[{mr, r := RandomInteger[{0, 1}]},
r := RandomInteger[{0, 1}]; (* define random number *)
mr := Table[r, {i, 1, n}, {j, 1, n}]; (*
form a matrix of random elements *)
Table[Det[mr], {i, 1, nn}] (* collect nn determinants *)]


Runs spitting out the spectrum for various small n are

Tally[detRnd[3, 10^3]] // Sort
Max[#[[1]] & /@ %]

(* Out[28]= {{-2, 6}, {-1, 163}, {0, 666}, {1, 161}, {2, 4}}

Out[29]= 2 *)

Tally[detRnd[4, 10^4]] // Sort
Max[#[[1]] & /@ %]

(* Out[41]= {{-3, 9}, {-2, 188}, {-1, 1534}, {0, 6594}, {1, 1466}, {2,
201}, {3, 8}}

Out[42]= 3 *)

Tally[detRnd[5, 10^4]] // Sort
Max[#[[1]] & /@ %]

(* Out[22]= {{-5, 1}, {-4, 12}, {-3, 52}, {-2, 354}, {-1, 1466}, {0,
6310}, {1, 1399}, {2, 343}, {3, 43}, {4, 18}, {5, 2}}

Out[23]= 5 *)


Now n=6, first try

Tally[detRnd[6, 10^4]] // Sort
Max[#[[1]] & /@ %]

(* Out[43]= {{-6, 7}, {-5, 5}, {-4, 63}, {-3, 119}, {-2, 527}, {-1,
1327}, {0, 5818}, {1, 1357}, {2, 572}, {3, 125}, {4, 63}, {5,
10}, {6, 6}, {8, 1}}

Out[44]= 8 *)


We find 8, but Max = 9 is correct. Hence we increase nn:

Tally[detRnd[6, 5*10^4]] // Sort
Max[#[[1]] & /@ %]

(* Out[49]= {{-9, 1}, {-8, 2}, {-7, 3}, {-6, 27}, {-5, 48}, {-4,
334}, {-3, 538}, {-2, 2786}, {-1, 6614}, {0, 29138}, {1, 6753}, {2,
2811}, {3, 549}, {4, 305}, {5, 55}, {6, 26}, {7, 5}, {8, 4}, {9, 1}}

Out[50]= 9 *)


Very nice, but for n=7 we must give up:

Tally[detRnd[7, 5*10^4]] // Sort
Max[#[[1]] & /@ %]

(* Out[56]= {{-16, 1}, {-14, 1}, {-12, 1}, {-11, 1}, {-10, 8}, {-9,
12}, {-8, 45}, {-7, 23}, {-6, 182}, {-5, 170}, {-4, 834}, {-3,
1084}, {-2, 3675}, {-1, 5982}, {0, 25921}, {1, 5976}, {2, 3795}, {3,
1030}, {4, 786}, {5, 161}, {6, 186}, {7, 46}, {8, 46}, {9,
15}, {10, 8}, {11, 3}, {12, 7}, {13, 1}}

Out[57]= 13 *)


13 is fine but it is very far from the correct Max = 32.

We have 2^(n^2) matrices, and it is a nice exercise to find out the number of runs to get the maximum determinant.

Maximize, space saving but very slow

maxDetMaximize[n_] := Module[{},
m = Table[a[i, j], {i, 1, n}, {j, 1, n}];
max = Maximize[{Det[m],
Table[a[i, j] == 0 || a[i, j] == 1, {i, 1, n}, {j, 1, n}]},
m // Flatten, Integers];
Print[max[[1]]];
mm = Table[a[i, j], {i, 1, n}, {j, 1, n}] /. max[[2]]]


Runs for n=2 and n=3 are quick, but n=4 runs very long, and I had to quit it.

Results

§1 We start by quoting from /4/

The Hadamard maximal determinant problem asks when a matrix of a given order with entries -1 and +1 has the largest possible determinant.

Despite well over a century of work by mathematicians, beginning with Sylvester's investigations of 1867, the question remains unanswered in general.

A discouraging statement of the experts. But let us carry on a bit.

A related (and even more difficult) problem is the determinant spectrum problem which asks, not just for the maximal determinant, but for the complete set of values taken by the determinant function.

This corresponds to my problem (b) (but is misses the question of multiplicity of determinants).

§ 2 Here's a list of largest determinants taken from OEIS /2/

Hadamard maximal determinant problem: largest determinant of a (real) {0,1}-matrix of order n.

a = (1, 1, 1, 2, 3, 5, 9, 32, 56, 144, 320, 1458, 3645, 9477, 25515, 131072, 327680, 1114112, 3411968, 19531250, 56640625)


Listplotting this (notice that the by convention the first element in the list is $$a(0)$$).

Looking carefully we can see, starting at $$n=3$$, consecutive groups of $$4$$ values each which climb up a staircase. This is probably related to the discovery proved by Hadamard "that that a(n) <= 2^(-n)*(n+1)^((n+1)/2), with equality if and only if a Hadamard matrix of order n+1 exists. It is believed that a Hadamard matrix of order n exists if and only if n = 1, 2 or a multiple of 4" (/1/).

The last statement indicates the special role of the number $$4$$

A Hadamard-matrix $$H$$ is defined as consisting of elements $$-1$$ and $$+1$$ which has the property that $$H.H^{T}= n E$$.

There's a Mathematica command HadamardMatrix[2^k] which works, however, only for powers of $$2$$ as indicated.

References

Here are some basic references, with my comments

/1/ https://oeis.org/A003432, Hadamard maximal determinant problem: largest determinant of a (real) {0,1}-matrix of order n.

/2/ http://mathworld.wolfram.com/HadamardsMaximumDeterminantProblem.html - A (-1,1)-matrix having a maximal determinant is known as a Hadamard matrix (Brenner and Cummings 1972). - The Mathematica-Code attached to this article does not run in my environment.

/4/ http://www.indiana.edu/~maxdet/: The Hadamard maximal determinant problem (Will Orrick's page)

/5/ http://www.indiana.edu/~maxdet/spectrum.html: Spectrum of the determinant function: Solution of part of my problem (b), lists all possible values of the determinante of (0,1)-matices for given $$n$$.

/6/ https://tomas.rokicki.com/newrec.html: "These are my new records for maximal determinants"

/7/ https://maths-people.anu.edu.au/~brent/maxdet/ : Richard Brent's web site, The Hadamard maximal determinant problem, contains a wealth of information about the problem, including data collected in numerous exhaustive and non-exhaustive searches.

/9/ https://en.wikipedia.org/wiki/Hadamard%27s_maximal_determinant_problem among other things equivalence of {-1,1}- und {0,1}-matrices

/10/ https://mathoverflow.net/questions/39786/maximum-determinant-of-0-1-valued-n-times-n-matrices: expert discussion 2010

/11/ https://mathoverflow.net/questions/18547/number-of-unique-determinants-for-an-nxn-0-1-matrix: expert discussion 2010

• You might use symmetry properties of the determinant (e.g. interchange of columns) to generate matrices in canonical order, vastly reducing search space – mikado Jun 23 '19 at 7:46
• Every row and every column needs at least one 1 or the Det is zero, so you can simplify generating the space by assuming the diagonal has all one's (you can always swap columns to force this). – MikeY Jun 24 '19 at 16:15
• Thanks for your proposals. They might improve the code for finding the maximum by reducing memory requirements, but don't forget part (b) of my question which asks for the spectrum of possible determinant values). – Dr. Wolfgang Hintze Jun 26 '19 at 9:28

Even though the sequence's first few terms are known and won't be extended by this solution, here's a way to throw Mathematica's integer math magic at the problem. This seems to use much less memory:

f[n_Integer?Positive] := Module[{v},
(* make a list of n^2 unique variables *)
v = Table[Unique[], n^2];
(* maximize the determinant under the given constraints *)
First@Maximize[{Det[Partition[v, n]], And @@ Thread[0 <= v <= 1]}, v, Integers]]


test:

f[1] // AbsoluteTiming
(*    {0.064168, 1}    *)

f[2] // AbsoluteTiming
(*    {0.045458, 1}    *)

f[3] // AbsoluteTiming
(*    {0.03024, 2}    *)

f[4] // AbsoluteTiming
(*    {2.25245, 3}    *)


### Better solution using the comment by @MikeY, and without much memory use:

Every row and every column needs at least one 1 or the Det is zero, so you can simplify generating the space by assuming the diagonal has all one's (you can always swap columns to force this).

First, define a function that, given a list of matrix elements, computes the determinant of the corresponding matrix with all ones on the diagonal:

q[n_Integer?Positive] := Module[{v, k, M},
(* make a list of n(n-1) unique variables *)
v = Table[Unique[], n (n - 1)];
(* make an nxn matrix with 1 on the diagonal *)
k = 0;
M = Table[If[i == j, 1, v[[++k]]], {i, n}, {j, n}];
(* compile the calculation of the determinant *)
Compile[Evaluate[{#, _Integer} & /@ v],
Evaluate[FullSimplify[Det[M]]],
CompilationTarget -> "C", RuntimeOptions -> "Speed"]]


Loop over all possible combinations of $$\{0,1\}$$ for the matrix elements, without ever constructing a list of these combinations:

h[n_Integer?Positive] := Module[{f, m = 0, x},
(* compile the determinant function *)
f = q[n];
(* loop over binary combinations *)
Do[If[(x = f @@ IntegerDigits[i, 2, n (n - 1)]) > m,
m = x; Echo[IntegerDigits[i, 2, n (n - 1)], m]],
{i, 0, 2^(n (n - 1)) - 1}];
(* return the highest determinant found *)
m]


Test:

h[4] // AbsoluteTiming
(*    {0.618653, 3}    *)

h[5] // AbsoluteTiming
(*    {20.4476, 5}    *)

• @ Roman Your code surely is compact and elegant. But (1) it didn't run unless I modified it toTable[Unnique[],{i,1,n^2}], and (2) when I tried to run f[5] my PC became trapped in a loop, and I had to switch the power off .Summary, nice but dangerous (at least on my PC). What is your result for f[5]? Can you go beyond n=5? – Dr. Wolfgang Hintze Jun 25 '19 at 22:07
• What is your \$Version? Is your Unnique a typo? – Roman Jun 25 '19 at 22:11
• @ Roman: "Even though the analytic solution is known". Can you be more precise, please? I am not aware of an analytic solution, and I doubt that thre is one. Compare my self solution, section results. – Dr. Wolfgang Hintze Jun 25 '19 at 22:11
• Sorry, yes it is a typo here in the comment. In the code it was correct, and the code ran without problems for n=2,3,4. – Dr. Wolfgang Hintze Jun 25 '19 at 22:13
• What I meant is: "even though the sequence has an OEIS number", sorry for the lack of precision. – Roman Jun 25 '19 at 22:13

Running with @Roman's answer, and noting that all rows and columns need at least two 1's in them, I assumed that the fixed part of the matrix is banded with 1's on the diagonal and the first band above (no guarantee this will give a valid answer!),

f2[n_Integer?Positive] := Module[{v, M, k},
(*make a list of n(n-1) unique variables*)
v = Table[Unique[], (n - 1)^2];
(*make an nxn matrix with 1 on the diagonal and first upper band*)
k = 0;
M = SparseArray[{
{i_, i_} -> 1,
{i_, j_} :> 1 /; i - j == -1,
{i_, j_} :> v[[++k]] /; i - j > 0,
{i_, j_} :> v[[++k]] /; i - j < -1
}, {n, n}];
(*maximize the determinant under the given constraints*)
First@Maximize[{Det[M], And @@ Thread[0 <= v <= 1]}, v, Integers]]


Running @Roman's code

f[5] // AbsoluteTiming
(*   {67.344, 5}   *)

f2[5] // AbsoluteTiming
(*   {3.61407, 5}   *)


It still bogs down on f2[6], but maybe a path forward? Tried a tridiagonal matrix, but that gave wrong answers.

EDIT

On the off chance that an optimal matrix from the $$n$$ dimension would form a kernel of the $$n+1$$ problem, I wrote this. It takes in a previous optimum and augments it with $$2n+1$$ unknowns and solves that, using @Roman's process. Works super fast and give the right answer for $$n=6$$ when given a good answer for $$n=5$$, but in general does not work. Shown for FYI.

f3[mat_] := Module[{n = Length@mat, v, vv, M, k, rules},
(*make a list of 2n+1 unique variables*)
v = Table[Unique[], 2 n + 1];

k = 0;
M = SparseArray[{
{n + 1, n + 1} :> v[[++k]],
{n + 1, j_} :> v[[++k]] /; j <= n,
{i_, n + 1} :> v[[++k]] /; i <= n,
{i_, j_} :> mat[[i, j]] /; i <= n && j <= n
}, {n + 1, n + 1}];

(*maximize the determinant under the given constraints*)
rules = Last@Maximize[{Det[M], And @@ Thread[0 <= v <= 1]}, v, Integers];

(* this is kludgy right here. I wanted to do M/.rules, but that doesn't work! *)

vv = v /. rules;
k = 0;

(* rebuild the matrix, bleh *)
SparseArray[{
{n + 1, n + 1} :> vv[[++k]],
{n + 1, j_} :> vv[[++k]] /; j <= n,
{i_, n + 1} :> vv[[++k]] /; i <= n,
{i_, j_} :> mat[[i, j]] /; i <= n && j <= n
}, {n + 1, n + 1}]

]


Starting from an optimal for $$n=5$$

$$\left( \begin{array}{ccccc} 1 & 1 & 0 & 1 & 0 \\ 0 & 1 & 1 & 0 & 1 \\ 0 & 0 & 1 & 1 & 0 \\ 0 & 0 & 0 & 1 & 1 \\ 1 & 0 & 1 & 0 & 1 \\ \end{array} \right)$$

It is fast and gives a right answer

(opt6 = f3[opt5]) // Det // AbsoluteTiming
(*   {0.0363656, 9}   *)


But fails in the next iteration (Det should be 32)

(opt7 = f3[opt6]) // Det // AbsoluteTiming
(*   {0.13524, 18}   *)


A one-line and very inefficient solution is of course

MaxDet[n_Integer] := Det /@ (Partition[#, n] & /@ Tuples[{0, 1}, n^2]) // Max
Table[MaxDet@n,{n,4}]//AbsoluteTiming
(*{1.6056, {1, 1, 2, 3}}*)


As @Roman points out, we can assume the matrix has all ones in the diagonal to reduce the search space

MaxDet2[n_Integer] :=  Det /@ Select[Partition[#, n] & /@ Tuples[{0, 1}, n^2], 1 == Times @@ Diagonal@# &] // Max
Table[MaxDet@n,{n,4}]//AbsoluteTiming
(*{0.331908, {1, 1, 2, 3}}*)


Computing MaxDet2[5] will still take a long time though:

MaxDet2@5//AbsoluteTiming
(*{252.608, 5}*)


We can also exploit other symmetries, e.g. the fact that the determinant is invariant under matrix transposition:

MinTranspose[m_] /; ! OrderedQ@{m, Transpose@m} := Transpose@m;
MinTranspose[m_] /; OrderedQ@{m, Transpose@m} := m;
MaxDet3[n_Integer] :=  Det /@ DeleteDuplicatesBy[Select[Partition[#, n] & /@ Tuples[{0, 1}, n^2], 1 == Times @@ Diagonal@# &], MinTranspose] // Max


But the benefit seems to be somewhat limited

MaxDet3@5//AbsoluteTiming
(*{248.654, 5}*)


I was hoping that one could reduce even further the search space by assuming that a maximal-determinant (n+1) x (n+1) matrix can be constructed using an n x n maximal one and adding 2n new elements in the first row and first column:

MatDesc[m_] :=
Flatten[
Table[Prepend[Transpose@Prepend[m, row], Prepend[col, 1]]
, {row, Tuples[{0, 1}, Length@m]}
, {col, Tuples[{0, 1}, Length@m]}]
, 1]

MaxMat[1] = {{{1}}};
MaxMat[n_Integer] := MaxMat[n] =
DeleteDuplicatesBy[MaximalBy[Flatten[MatDesc /@ MaxMat[n - 1], 1], Det], MinTranspose]


Unfortunately, this construction breaks down already for 7 x 7 matrices:

Det@First@MaxMat@7 // AbsoluteTiming
(*{107.684, 18}*)


Edit: Regardless of the approach, as we are dealing with small matrices it is always convenient to pre-compile the determinant expression (most of the time is spent in Det and not generating matrices):

CDet[n_Integer] := CDet[n] =
Block[{mat, det, var},
mat = Table[x[i, j], {i, n}, {j, n}];
det = Expand@Det@mat;
var = {#, _Integer} & /@ Flatten@mat;
Compile[Evaluate@var, Evaluate@det, CompilationTarget -> "C", RuntimeOptions -> "Speed"]
]


Note that Expand above makes a big difference. This allows even a very naive approach to reach n = 5

MaxCDet[n_Integer] := CDet[n] @@@ Tuples[{0, 1}, n^2] // Max
MaxCDet@5 // AbsoluteTiming
(*{144.214, 5}*)


Following @Roman again, we can then use base-2 representations, binary operations and an iterative approach (to avoid using too much memory):

Block[{n = 5, max = 0, id, i},
id = FromDigits[Flatten@IdentityMatrix@n, 2];
i = id;
While[i < 2^n^2,
max = Max[max, CDet[n] @@ IntegerDigits[i, 2, n^2]];
i = BitOr[i + 1, id];
];
max
] // AbsoluteTiming
(*{5.88838, 5}*)


Note: Of course, finding examples of small matrices with maximal determinants is a lot easier, since they are not too uncommon once you discard most of the null-determinant matrices. You can then just expect to find a maximal-determinant matrix if you try enough random instances

RandMat[n_Integer] := IdentityMatrix@n + RandomInteger[{0, 1}, {n, n}] /. 2 -> 1
MaxDet4[k_Integer][n_Integer] := Max@Table[Det@RandMat@n, {k}]

MaxDet4[10000]@5 // AbsoluteTiming
(*{0.941558, 5}*)
MaxDet4[100000]@6 // AbsoluteTiming
(*{13.7177, 9}*)


But of course this fails once they become too rare:

MaxDet4[1000000]@7 // AbsoluteTiming
(*{185.917, 18}*)