Find the maximum determinant for each size Toeplitz matrix

For a fixed n, consider the n by n Toeplitz matrices with entries which are either 0 or 1. The aim is to find maximum determinant over all such Toeplitz matrices.

For each n from 1 upwards, output the maximum determinant over all n by n Toeplitz matrices with entries which are either 0 or 1. There should be one output per n which should have the maximum determinant and also an example matrix that reaches it.

Score

Your score is the largest n your code gets to in 2 minutes on my computer. To clarify a little, your code can run for 2 minutes in total, this is not 2 minutes per n.

Tie breaker

If two entries get the same n score then the winning entry will be the one that gets to the highest n in the shortest time on my machine. If the two best entries are equal on this criterion too then the winner will be the answer submitted first.

Languages and libraries

You can use any freely available language and libraries you like. I must be able to run your code so please include a full explanation for how to run/compile your code in linux if at all possible.

My Machine The timings will be run on my machine. This is a standard ubuntu install on an AMD FX-8350 Eight-Core Processor. This also means I need to be able to run your code.

For n = 1..10 the outputs should be 1,1,2,3,5,9,32,56,125,315

This sequence is not in OEIS and so the winning entry also gets to propose a new entry there.

Entries so far

• n=10 n=11 by Vioz in Python
• n=9 by Tyilo in C
• n=12 by Legendre in J
• n=10 by Tensibai in R
• n=14 by SteelRaven in C++
• n=14 by RetoKoradi in C++
• @AlexA. You are right and I have apologised. Luckily the two problems are very similar so he should easily be able to modify his code. – user9206 Jul 9 '15 at 16:47
• The solution by @Vioz comes up with a sequence that starts with 1, 1, 2, 3, 5, 9, 32. So the value for n=5 is different from what you list. Since all other values match, it looks like the solution is probably correct, and this is just a typo in the question? – Reto Koradi Jul 9 '15 at 18:56
• @RetoKoradi Thank you. Fixed. – user9206 Jul 9 '15 at 19:03
• Here are 10 possible binary Toeplitz matrices with maximum determinants for n = 1..10: ghostbin.com/paste/axkpa – Tyilo Jul 10 '15 at 3:55
• As an observation that may help others but I can't verify beyond 14. It appears respective means of the top row and the first column of the Toeplitz matrix are always 0.4 <= m <= 0.6 for the maximum determinant. – MickyT Jul 16 '15 at 1:52

This gets to n=14 in just under 1 minute on my machine. But since that's just a 2-core laptop, I hope that the 8-core test machine can finish n=15 in under 2 minutes. It takes about 4:20 minutes on my machine.

I was really hoping to come up with something more efficient. There has got to be a way to calculate the determinate of a binary matrix more efficiently. I wanted to come up with some kind of dynamic programming approach that counts the +1 and -1 terms in the determinant calculation. But it just hasn't quite come together so far.

Since the bounty is about to expire, I implemented the standard brute force approach:

• Loop over all possible Toeplitz matrices.
• Skip one of the two in each transposed matrix pair. Since the matrix is described by bitmask values, this is simple to do by skipping all values where the reverse of the bitmask is smaller than the bitmask itself.
• The determinate is calculated with a text book LR decomposition. Except for some minor performance tuning, the main improvement I made to the algorithm from my college numerical methods book is that I use a simpler pivot strategy.
• Parallelization is done with pthreads. Just using regular spacing for the values processed by each thread caused very bad load balancing, so I introduced some swizzling.

I tested this on Mac OS, but I used similar code on Ubuntu before, so I hope this will compile and run without a hitch:

1. Save the code in a file with a .cpp extension, e.g. optim.cpp.
2. Compile with gcc -Ofast optim.cpp -lpthread -lstdc++.
3. Run with time ./a.out 14 8. The first argument is the maximum n. 14 should finish in under 2 minutes for sure, but it would be great if you could try 15 as well. The second argument is the number of threads. Using the same value as the number of cores of the machine is normally a good start, but trying some variations could potentially improve the times.

Let me know if you have any problems building or running the code.

#include <stdint.h>
#include <cstdlib>
#include <iostream>

static int NMax = 14;

static int BarrierCount = 0;

static float* MaxDetA;
static uint32_t* MaxDescrA;

static inline float absVal(float val)
{
return val < 0.0f ? -val : val;
}

static uint32_t reverse(int n, uint32_t descr)
{
uint32_t descrRev = 0;
for (int iBit = 0; iBit < 2 * n - 1; ++iBit)
{
descrRev <<= 1;
descrRev |= descr & 1;
descr >>= 1;
}

return descrRev;
}

static void buildMat(int n, float mat[], uint32_t descr)
{
int iDiag;
for (iDiag = 1 - n; iDiag < 0; ++iDiag)
{
float val = static_cast<float>(descr & 1);
descr >>= 1;
for (int iRow = 0; iRow < n + iDiag; ++iRow)
{
mat[iRow * (n + 1) - iDiag] = val;
}
}

for ( ; iDiag < n; ++iDiag)
{
float val = static_cast<float>(descr & 1);
descr >>= 1;
for (int iCol = 0; iCol < n - iDiag; ++iCol)
{
mat[iCol * (n + 1) + iDiag * n] = val;
}
}
}

static float determinant(int n, float mat[])
{
float det = 1.0f;
for (int k = 0; k < n - 1; ++k)
{
float maxVal = 0.0f;
int pk = 0;
for (int i = k; i < n; ++i)
{
float q = absVal(mat[i * n + k]);
if (q > maxVal)
{
maxVal = q;
pk = i;
}
}

if (pk != k)
{
det = -det;
for (int j = 0; j < n; ++j)
{
float t = mat[k * n + j];
mat[k * n + j] = mat[pk * n + j];
mat[pk * n + j] = t;
}
}

float s = mat[k * n + k];
det *= s;

s = 1.0f / s;
for (int i = k + 1; i < n; ++i)
{
mat[i * n + k] *= s;
for (int j = k + 1; j < n; ++j)
{
mat[i * n + j] -= mat[i * n + k] * mat[k * n + j];
}
}
}

det *= mat[n * n - 1];

return det;
}

{

++BarrierCount;
{
}
else
{
BarrierCount = 0;
}

}

{

float* mat = new float[NMax * NMax];

for (int n = 1; n <= NMax; ++n)
{
uint32_t descrRange(1u << (2 * n - 1));
float maxDet = 0.0f;
uint32_t maxDescr = 0;

for (uint32_t descrBase = 0;
descrBase + descrInc < descrRange;
{
uint32_t descr = descrBase + descrInc;
descrInc = (descrInc + 1) % ThreadCount;

if (reverse(n, descr) > descr)
{
continue;
}

buildMat(n, mat, descr);
float det = determinant(n, mat);
if (det > maxDet)
{
maxDet = det;
maxDescr = descr;
}
}

// Let main thread output results.
}

delete[] mat;

return 0;
}

static void printMat(int n, float mat[])
{
for (int iRow = 0; iRow < n; ++iRow)
{
for (int iCol = 0; iCol < n; ++iCol)
{
std::cout << " " << mat[iRow * n + iCol];
}
std::cout << std::endl;
}

std::cout << std::endl;
}

int main(int argc, char* argv[])
{
if (argc > 1)
{
NMax = atoi(argv[1]);
if (NMax > 16)
{
NMax = 16;
}
}

if (argc > 2)
{
}

{
}

float* mat = new float[NMax * NMax];

for (int n = 1; n <= NMax; ++n)
{

float maxDet = 0.0f;
uint32_t maxDescr = 0;

{
{
}
}

std::cout << "n = " << n << " det = " << maxDet << std::endl;
buildMat(n, mat, maxDescr);
printMat(n, mat);

}

delete[] mat;

delete[] MaxDetA;
delete[] MaxDescrA;

return 0;
}

• There's an interesting way of calculating the determinant of an integer matrix using only integer arithmetic: LU decomposition in some finite field (basically mod a large prime). I don't know if this would be faster. – lirtosiast Jul 19 '15 at 3:56
• @ThomasKwa That would probably still be O(n^3)? It might be helpful for larger matrices where floating point precision would otherwise become an issue. I didn't really look for literature. Well, I did a quick search, and found a paper about calculating determinants of Toeplitz matrices. But there were too many open questions for me to commit the time to try and implement it. – Reto Koradi Jul 19 '15 at 4:08
• @Lembik I'll try and take a look at it later today. I changed it to handle larger sizes for your other related challenge yesterday. Couldn't beat the highest scores for n=30 so far, my heuristics are stuck below 5*10^13. – Reto Koradi Jul 20 '15 at 17:11
• @Lembik See paste.ubuntu.com/11915546 for code and paste.ubuntu.com/11915532 for results up to n=19. – Reto Koradi Jul 21 '15 at 17:23
• @Lembik The results up to n=20 are at paste.ubuntu.com/11949738. They now list all the tied solutions, including attributes to quickly see the value of the diagonal and whether they are circulant. All of the maxima for m=18,19,20 are circulant matrices. Please double check the determinants before publishing them anywhere. – Reto Koradi Jul 27 '15 at 23:31

J

Update: Improved code to search over half the values. Now calculates n=12 comfortably within 120 seconds (down from 217s to 60s).

#!/usr/bin/jconsole

dim =: -:@>:@#
take =: i.@dim
rotstack =: |."0 1~ take
toep =: (dim (|."1 @: {."1) rotstack)"1
det =: -/ . * @: toep
ps =: 3 : ',/(0 1 ,"0 1/ ,.y)'
canonical =: #. >: [: #. |. " 1

lss =: 3 : 0
shape =. (2^y), y
shape $,>{;/(y,2)$0 1
)

ls =: (canonical@:lss) # lss
ans =: >./ @: det @: ls @: <: @: +:

display =: 3 : 0
echo 'n = ';y;'the answer is';ans y
)
display"0 (1 + i.13)
exit''


Run this and kill when two minutes are up. My results (MBP 2014 - 16GB of RAM):

┌────┬─┬─────────────┬─┐
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬─┐
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬─┐
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬─┐
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬─┐
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬─┐
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬──┐
└────┴─┴─────────────┴──┘
┌────┬─┬─────────────┬──┐
└────┴─┴─────────────┴──┘
┌────┬─┬─────────────┬───┐
└────┴─┴─────────────┴───┘
┌────┬──┬─────────────┬───┐
└────┴──┴─────────────┴───┘
┌────┬──┬─────────────┬────┐
└────┴──┴─────────────┴────┘
┌────┬──┬─────────────┬────┐
└────┴──┴─────────────┴────┘


Total run time = 61.83s.

Just for fun

┌────┬──┬─────────────┬────┐
└────┴──┴─────────────┴────┘


This took approximately 210 seconds on its own.

• Note to testers: n = 12 requires roughly 18 GiB of memory. – Dennis Jul 11 '15 at 16:49
• This is a very nice improvement. The output is slightly buggy however. For me using j64-804 it outputs n =1 twice so it is out by one for ever more. – user9206 Jul 12 '15 at 19:53
• @Lembik Ah that's right. I just updated the code; could you try running again? Thanks! (I have set it to calculate up to n=13. You can change the 13 in the second-to-last line to have it calculate up whatever you wish.) – Legendre Jul 12 '15 at 22:07
• I ran it again and it still gets to 12. – user9206 Jul 14 '15 at 6:42
• @Lembik Hmm .. are you saying it gets to 12 within the time limit and gets to 13 some time after that (which is what I expect), or that it never gets to 13 (i.e. program halts after 12)? – Legendre Jul 14 '15 at 6:46

Python 2

This is a very straightforward solution, and probably won't win the contest. But hey, it works!

I'll give a quick overview of what exactly is happening.

1. I first generate every possible starting row for n. For example, when n=2, this will generate an array length 2n+1, where each row is length 2n-1. It would look like this: [[0,0,0],[0,0,1],[0,1,0],[0,1,1],[1,0,0],[1,0,1],[1,1,0],[1,1,1]].
2. Then, for each of those possible starting rows, I rotate it n times and slice off the first n items to generate the appropriate matrix, and use scipy to calculate the determinant, all while keeping track of the maximum value. At the end of this, I simply print out the maximum, increment n by 1, and keep going until 10 minutes has passed.

To run this, you will need scipy installed.

Edit 1: Changed how initial rows were built by using itertools.product instead, thanks Sp3000!

Edit 2: Removed storage of possible starting rows for a minimal improvement in speed.

Edit 3: Changed to scipy to have more control over how det worked.

from scipy import linalg
from collections import deque
from time import time
from itertools import product

c=1
t=time()
while 1:
m=0
for d in product(range(2),repeat=2*c-1):
a=deque(d)
l=[d[0:c]]
for x in xrange(c-1):
a.rotate(1)
l+=[list(a)[0:c]]
m=max(m,linalg.det(l,overwrite_a=True,check_finite=False))
print m,'in',time()-t,'s'
c+=1


Here's some sample output on my home machine (i7-4510U, 8GB RAM):

1.0 in 0.0460000038147 s
1.0 in 0.0520000457764 s
2.0 in 0.0579998493195 s
3.0 in 0.0659999847412 s
5.0 in 0.0829999446869 s
9.0 in 0.134999990463 s
32.0 in 0.362999916077 s
56.0 in 1.28399991989 s
125.0 in 5.34999990463 s
315.0 in 27.6089999676 s
1458.0 in 117.513000011 s

• Thank you but I think you have answered an old version of the question I am afraid. It is now about Toeplitz matrices and the time limit is 2 minutes. – user9206 Jul 9 '15 at 16:19
• I see so much golfed Python on this site that I often forget that for general purposes it's actually a readable language. – Alex A. Jul 9 '15 at 18:25
• This could probably be sped up significantly, because it doesn't take advantage of the fact that it's a binary matrix. – lirtosiast Jul 9 '15 at 19:11
• @ThomasKwa If I'm honest, I have no idea how to take advantage of that :P – Kade Jul 9 '15 at 19:32
• Quote from numpy documentation: "The determinant is computed via LU factorization using the LAPACK routine z/dgetrf." I looked at dgetrf, and it says it uses double precision; depending on OP's GPU single precision could be faster. – lirtosiast Jul 10 '15 at 19:29

C++

Bruteforce with use of OpenMP for parallelization and simple optimization to avoid evaluation of determinant for transposed matrices.

$lscpu ... Model name: Intel(R) Core(TM) i5-2410M CPU @ 2.30GHz ...$ g++ -O2 toepl.cpp -fopenmp
$timeout 2m ./a.out 1 1 2 1 3 2 4 3 5 5 6 9 7 32 8 56 9 125 10 315 11 1458 12 2673 13 8118 14 22386  #include <cmath> #include <algorithm> #include <iostream> #include <vector> using namespace std; void updateReverses(vector < int > & reverses) { int reversesCnt = reverses.size(); for(int i = 0; i < reversesCnt; ++i){ reverses[i] <<= 1; reverses.push_back(reverses[i] | 1); } } const double eps = 1e-9; double determinant(vector < vector < double > > & matrix) { int n = matrix.size(); double det = 1; if(n == 1) return matrix[0][0]; for(int i = 0; i < n; ++i){ int p = i; for(int j = i + 1; j < n; ++j) if(fabs(matrix[j][i]) > fabs(matrix[p][i])) p = j; if(fabs(matrix[p][i]) < eps) return 0; matrix[i].swap(matrix[p]); if(i != p) det *= -1; det *= matrix[i][i]; matrix[i][i] = 1. / matrix[i][i]; for(int j = i + 1; j < n; ++j) matrix[i][j] *= matrix[i][i]; for(int j = i + 1; j < n; ++j){ if(fabs(matrix[j][i]) < eps) continue; for(int k = i + 1; k < n; ++k) matrix[j][k] -= matrix[i][k] * matrix[j][i]; } } return det; } int main() { vector < int > reverses(1, 0); reverses.reserve(1 << 30); updateReverses(reverses); for(int n = 1;; ++n){ double res = 0; int topMask = 1 << (2 * n - 1); vector < vector < double > > matrix(n, vector < double > (n)); #pragma omp parallel for reduction(max:res) firstprivate(matrix) schedule(dynamic,1<<10) for(int mask = 0; mask < topMask; ++mask){ if(mask < reverses[mask]) continue; for(int i = 0; i < n; ++i) for(int j = 0; j < n; ++j) matrix[i][j] = (mask >> (i - j + n - 1)) & 1; res = max(res, determinant(matrix)); } cout << n << ' ' << res << endl; updateReverses(reverses); updateReverses(reverses); } }  • Looks like you may soon be making your first OEIS entry unless someone comes up with a clever idea :) – user9206 Jul 16 '15 at 0:23 C Compile with: $ clang -Ofast 52851.c -o 52851


Run with:

\$ ./52851


Can output the maximum determinant for n = 1..10 in ~115 seconds on my computer.

The program is just getting the determinant every possible binary Toeplitz matrix of size n, however every determinant of matrices of size 5x5 or smaller will be cached using memoization.

At first I wrongly assumed that every submatrix of a Toeplitz matrix will also be a Toeplitz matrix, so I only needed to memoize 2^(2n-1) values instead of 2^(n^2) for each n. I made the program before realizing my mistake, so this submission is just a fix of that program.

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <limits.h>
#include <string.h>

#define ELEMENTS(x) (sizeof(x) / sizeof(*x))

int *dets[6];

void print_matrix(int n, int c) {
for(int row = 0; row < n; row++) {
for(int col = 0; col < n; col++) {
int j = n - 1 - row + col;
int val = !!(c & (1 << j));
printf("%d ", val);
}
puts("");
}
}

int det(int n, uint8_t *m) {
if(n == 1) {
return m[0];
}

int i = 0;

if(n < ELEMENTS(dets)) {
for(int j = 0; j < n * n; j++) {
i *= 2;
i += m[j];
}

int v = dets[n][i];
if(v != INT_MIN) {
return v;
}
}

int v = 0;

uint8_t *sub = malloc((n - 1) * (n - 1));

for(int removed = 0; removed < n; removed++) {
if(m[removed]) {
uint8_t *p = sub;
for(int row = 1; row < n; row++) {
for(int col = 0; col < n; col++) {
if(col == removed) {
continue;
}

*p = m[col + row * n];

p++;
}
}

v += (removed % 2 == 0? 1: -1) * det(n - 1, sub);
}
}

free(sub);

if(n < ELEMENTS(dets)) {
dets[n][i] = v;
}
return v;
}

int main(void) {
for(int i = 2; i < ELEMENTS(dets); i++) {
int combinations = 1 << (i * i);
dets[i] = malloc(combinations * sizeof(**dets));
for(int j = 0; j < combinations; j++) {
dets[i][j] = INT_MIN;
}
}

puts("1: 1");

for(int n = 2; n < 65; n++) {
int vars = 2 * n - 1;
size_t combinations = 1 << vars;

int best = -1;
int max = -1;

uint8_t *sub = malloc((n - 1) * (n - 1));

for(int c = 0; c < combinations; c++) {
int d = 0;
for(int i = 0; i < n; i++) {
if(c & (1 << (n - 1 + i))) {
uint8_t *p = sub;
for(int row = 1; row < n; row++) {
for(int col = 0; col < n; col++) {
if(col == i) {
continue;
}

int j = n - 1 - row + col;
*p = !!(c & (1 << j));

p++;
}
}
d += (i % 2 == 0? 1: -1) * det(n - 1, sub);
}
}

if(d > max) {
max = d;
best = c;
}
}

free(sub);

printf("%d: %d\n", n, max);
//print_matrix(n, best);
}

return 0;
}

• It looks like you're calculating the determinant using expansion by minors; this has O(n!) complexity, so you may be better off using a different algorithm. – lirtosiast Jul 10 '15 at 19:17
• @ThomasKwa I didn't know that there were faster algorithms, so yeah this solution is pretty bad. – Tyilo Jul 11 '15 at 10:35
• You may want to look into using LU Decomposition to find the determinant of a matrix. It's O(n^3), I believe, though can be made faster with some interesting algorithms. I believe most of the builtins used here generally use a variant of decomposition to perform determinants. – BrainSteel Jul 11 '15 at 16:13
• @BrainSteel, yeah I looked at it, but I might as well go for a O(n^2) algorithm if I'm updating my answer. – Tyilo Jul 11 '15 at 16:22
• According to a casual Wikipedia search, the determinant of a Toeplitz matrix can be determined in O(n^2). But I think the bottleneck of the problem is searching among the O(4^n) many 0-1 n by n matrices. – Legendre Jul 13 '15 at 8:34

R

You'll have to install R and the packages listed with install.packages("package_name")

Didn't get under 2 mins on my machine with this version (I've to try with a parallel modification)

library(pracma)
library(stringr)
library(R.utils)
library(microbenchmark)

f <- function(n) {
#If n is 1, return 1 to avoid code complexity on this special case
if(n==1) { return(1) }
# Generate matrices and get their determinants
dets <- sapply(strsplit(intToBin( 0:(2^n - 1)), ""), function(x) {
sapply( strsplit( intToBin( 0:(2^(n-1) - 1) ), ""),
function(y) {
det(Toeplitz(x,c(x[1],y)))
})

})
#Get the maximum determinant and return it
res <- max(abs(dets))
return(res)
}


Call and output:

> sapply(1:10,f)
[1]   1   1   2   3   5   9  32  56 125 315


Benchmark on my machine:

> microbenchmark(sapply(1:10,f),times=1L)
Unit: seconds
expr      min       lq     mean   median       uq      max neval
sapply(1:10, f) 66.35315 66.35315 66.35315 66.35315 66.35315 66.35315     1


For information, for a 1:11 range, it takes 285 seconds.

PARI/GP, n=11

This is brute force but taking advantage of det(A^T) = det(A). I'm only posting it to demonstrate how easy it is to skip transposes. The lowest bit of b1 holds the top left cell, and the other bits hold the rest of the top row. b2 holds the rest of the left column. We simply enforce b2 <= (b1>>1).

{ for(n=1,11,
res=0;
for(b1=0,2^n-1,
for(b2=0,b1>>1,
res=max(res,matdet(matrix(n,n,i,j,bittest(if(i>j,b2>>(i-j-1),b1>>(j-i)),0))));
)
);
print(n" "res);
)
}


Regarding computing Toeplitz determinants in O(n^2) time: In my limited research, I've kept running into a requirement that all leading principal minors must be nonzero in order for the algorithms to work, which is a major obstacle for this task. Feel free to give me pointers if you know more about this than I do.

• Did you see this paper? scienpress.com/upload/JAMB/Vol%201_1_4.pdf . It wasn't clear to me what the complexity is. There seemed to be quite a lot of terms for the n=5 example. – Reto Koradi Jul 19 '15 at 14:37
• @RetoKoradi Yes I saw that. It seems the complexity is not polynomial, considering that e.g. e_{k+1} has 4 times the number of components as e_k. There are many omissions in the paper. An invertible matrix has an LU decomposition iff all leading principal minors are nonzero. (Notice the denominators, e.g. a_0 -- implicitly they are guaranteed to be nonzero.) Uniqueness comes from L being unit triangular. The author also didn't mention numerical stability. In case the link becomes unavailable, the paper is "On Calculating the Determinants of Toeplitz Matrices" by Hsuan-Chu Li (2011). – Mitch Schwartz Jul 19 '15 at 21:05