Pour un n fixe, considérons les matrices de Toeplitz n par n avec des entrées qui sont soit 0 soit 1. Le but est de trouver le déterminant maximum sur toutes ces matrices de Toeplitz.

Tâche

Pour chacun à npartir de 1, émettez le déterminant maximum sur toutes les matrices Toeplitz n par n avec des entrées qui sont soit 0 ou 1. Il devrait y avoir une sortie pour nlaquelle doit avoir le déterminant maximum et également un exemple de matrice qui l'atteint.

But

Votre score est le plus élevé atteint par nvotre code en 2 minutes sur mon ordinateur. Pour clarifier un peu, votre code peut s'exécuter pendant 2 minutes au total, ce n'est pas 2 minutes par n.

Tie break

Si deux entrées obtiennent le même nrésultat alors l'entrée gagnante sera celle qui obtient le plus haut ndans le temps le plus court sur ma machine. Si les deux meilleures entrées sont également égales sur ce critère, le gagnant sera la réponse soumise en premier.

Langues et bibliothèques

Vous pouvez utiliser n'importe quelle langue et bibliothèque librement disponibles. Je dois être en mesure d'exécuter votre code, veuillez donc inclure une explication complète sur la façon d'exécuter / compiler votre code sous Linux si possible.

Ma machine Les horaires seront exécutés sur ma machine. Il s'agit d'une installation ubuntu standard sur un processeur AMD FX-8350 à huit cœurs. Cela signifie également que je dois pouvoir exécuter votre code.

Petites réponses

Pour n = 1..10, les sorties devraient être 1,1,2,3,5,9,32,56,125,315

Cette séquence n'est pas dans OEIS et donc l'entrée gagnante peut également y proposer une nouvelle entrée.

Entrées jusqu'à présent

• n=10 n=11par Vioz en Python
• n=9par Tyilo en C
• n=12par Legendre en J
• n=10par Tensibai dans R
• n=14par SteelRaven en C ++
• n=14par RetoKoradi en C ++

C ++ avec pthreads

Cela arrive à n = 14 en un peu moins d'une minute sur ma machine. Mais comme il ne s'agit que d'un ordinateur portable à 2 cœurs, j'espère que la machine de test à 8 cœurs pourra terminer n = 15 en moins de 2 minutes. Cela prend environ 4:20 minutes sur ma machine.

J'espérais vraiment trouver quelque chose de plus efficace. Il a obtenu d'être un moyen de calculer la déterminée d'une matrice binaire de manière plus efficace. Je voulais trouver une sorte d'approche de programmation dynamique qui compte les termes +1 et -1 dans le calcul des déterminants. Mais cela ne s'est tout simplement pas encore réalisé jusqu'à présent.

Puisque la prime est sur le point d'expirer, j'ai implémenté l'approche standard de la force brute:

• Boucle sur toutes les matrices Toeplitz possibles.
• Sautez l'un des deux dans chaque paire de matrices transposées. Étant donné que la matrice est décrite par des valeurs de masque binaire, cela est simple à faire en ignorant toutes les valeurs où l'inverse du masque binaire est plus petit que le masque binaire lui-même.
• Le déterminé est calculé avec une décomposition LR du manuel. À l'exception de quelques ajustements mineurs des performances, la principale amélioration que j'ai apportée à l'algorithme de mon livre de méthodes numériques est que j'utilise une stratégie de pivot plus simple.
• La parallélisation se fait avec pthreads. Le simple fait d'utiliser un espacement régulier pour les valeurs traitées par chaque thread a causé un très mauvais équilibrage de charge, j'ai donc introduit un peu de swizzling.

J'ai testé cela sur Mac OS, mais j'ai déjà utilisé du code similaire sur Ubuntu, donc j'espère que cela se compilera et fonctionnera sans accroc:

1. Enregistrez le code dans un fichier avec une .cppextension, par exemple optim.cpp.
2. Compilez avec gcc -Ofast optim.cpp -lpthread -lstdc++.
3. Courez avec time ./a.out 14 8. Le premier argument est le maximum n. 14 devrait se terminer en moins de 2 minutes, bien sûr, mais ce serait bien si vous pouviez aussi essayer 15. Le deuxième argument est le nombre de threads. Utiliser la même valeur que le nombre de cœurs de la machine est normalement un bon début, mais essayer certaines variations pourrait potentiellement améliorer les délais.

Faites-moi savoir si vous rencontrez des problèmes pour créer ou exécuter le code.

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

static int NMax = 14;
static int ThreadCount = 4;

static pthread_mutex_t ThreadMutex;
static pthread_cond_t ThreadCond;
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;
}

static void threadBarrier()
{
    pthread_mutex_lock(&ThreadMutex);

    ++BarrierCount;
    if (BarrierCount <= ThreadCount)
    {
        pthread_cond_wait(&ThreadCond, &ThreadMutex);
    }
    else
    {
        pthread_cond_broadcast(&ThreadCond);
        BarrierCount = 0;
    }

    pthread_mutex_unlock(&ThreadMutex);
}

static void* threadFunc(void* pData)
{
    int* pThreadIdx = static_cast<int*>(pData);
    int threadIdx = *pThreadIdx;

    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;

        uint32_t descrInc = threadIdx;
        for (uint32_t descrBase = 0;
             descrBase + descrInc < descrRange;
             descrBase += ThreadCount)
        {
            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;
            }
        }

        MaxDetA[threadIdx] = maxDet;
        MaxDescrA[threadIdx] = maxDescr;

        threadBarrier();
        // Let main thread output results.
        threadBarrier();
    }

    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)
    {
        ThreadCount = atoi(argv[2]);
    }

    MaxDetA = new float[ThreadCount];
    MaxDescrA = new uint32_t[ThreadCount];

    pthread_mutex_init(&ThreadMutex, 0);
    pthread_cond_init(&ThreadCond, 0);

    int* threadIdxA = new int[ThreadCount];
    pthread_t* threadA = new pthread_t[ThreadCount];

    for (int iThread = 0; iThread < ThreadCount; ++iThread)
    {
        threadIdxA[iThread] = iThread;
        pthread_create(threadA + iThread, 0, threadFunc, threadIdxA + iThread);
    }

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

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

        float maxDet = 0.0f;
        uint32_t maxDescr = 0;

        for (int iThread = 0; iThread < ThreadCount; ++iThread)
        {
            if (MaxDetA[iThread] > maxDet)
            {
                maxDet = MaxDetA[iThread];
                maxDescr = MaxDescrA[iThread];
            }
        }

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

        threadBarrier();
    }

    delete[] mat;

    delete[] MaxDetA;
    delete[] MaxDescrA;

    delete[] threadIdxA;
    delete[] threadA;

    return 0;
}

J

Mise à jour: code amélioré pour rechercher plus de la moitié des valeurs. n=12Calcule maintenant confortablement en 120 secondes (de 217 à 60).

Vous aurez besoin de la dernière version de J installée.

#!/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''

Exécutez ceci et tuez quand deux minutes sont écoulées. Mes résultats (MBP 2014 - 16 Go de RAM):

┌────┬─┬─────────────┬─┐
│n = │1│the answer is│1│
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬─┐
│n = │2│the answer is│1│
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬─┐
│n = │3│the answer is│2│
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬─┐
│n = │4│the answer is│3│
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬─┐
│n = │5│the answer is│5│
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬─┐
│n = │6│the answer is│9│
└────┴─┴─────────────┴─┘
┌────┬─┬─────────────┬──┐
│n = │7│the answer is│32│
└────┴─┴─────────────┴──┘
┌────┬─┬─────────────┬──┐
│n = │8│the answer is│56│
└────┴─┴─────────────┴──┘
┌────┬─┬─────────────┬───┐
│n = │9│the answer is│125│
└────┴─┴─────────────┴───┘
┌────┬──┬─────────────┬───┐
│n = │10│the answer is│315│
└────┴──┴─────────────┴───┘
┌────┬──┬─────────────┬────┐
│n = │11│the answer is│1458│
└────┴──┴─────────────┴────┘
┌────┬──┬─────────────┬────┐
│n = │12│the answer is│2673│
└────┴──┴─────────────┴────┘

Durée totale d'exécution = 61,83 s.

Juste pour le fun

┌────┬──┬─────────────┬────┐
│n = │13│the answer is│8118│
└────┴──┴─────────────┴────┘

Cela a pris environ 210 secondes à lui tout seul.

Python 2

Il s'agit d'une solution très simple et ne gagnera probablement pas le concours. Mais bon, ça marche!

Je vais donner un bref aperçu de ce qui se passe exactement.

1. Je génère d'abord toutes les lignes de départ possibles pour n. Par exemple, lorsque n=2cela générera une longueur de tableau 2 ^{n + 1} , où chaque ligne est de longueur 2n-1. Il ressemblerait à ceci: [[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. Ensuite, pour chacune de ces lignes de départ possibles, je la fais pivoter nfois et je coupe les premiers néléments pour générer la matrice appropriée, et utiliser scipypour calculer le déterminant, tout en gardant une trace de la valeur maximale. À la fin de cela, j'imprime simplement le maximum, incrémente nde 1 et continue jusqu'à ce que 10 minutes se soient écoulées.

Pour l'exécuter, vous aurez besoin de scipy installé.

Edit 1: Modification de la façon dont les lignes initiales ont été construites en utilisant itertools.product à la place, merci Sp3000!

Edit 2: Suppression du stockage des lignes de départ possibles pour une amélioration minimale de la vitesse.

Edit 3: Changé pour scipyavoir plus de contrôle sur le detfonctionnement.

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

Voici quelques exemples de sortie sur ma machine domestique (i7-4510U, 8 Go de 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

C ++

Bruteforce avec l'utilisation d'OpenMP pour la parallélisation et l'optimisation simple pour éviter l'évaluation du déterminant pour les matrices transposées.

$ 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);
  }
}

C

Compiler avec:

$ clang -Ofast 52851.c -o 52851

Courir avec:

$ ./52851

Peut sortir le déterminant maximum pendant n = 1..10~ 115 secondes sur mon ordinateur.

Le programme obtient simplement le déterminant de chaque matrice binaire Toeplitz de taille possible n, mais chaque déterminant de matrices de taille 5x5ou plus petit sera mis en cache à l'aide de la mémorisation.

Au début, j'ai supposé à tort que chaque sous-matrice d'une matrice Toeplitz serait également une matrice Toeplitz, donc je n'avais besoin que de mémoriser des 2^(2n-1)valeurs au lieu de 2^(n^2)pour chacune n. J'ai créé le programme avant de réaliser mon erreur, donc cette soumission n'est qu'un correctif de ce programme.

#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;
}

R

Vous devrez installer R et les packages répertoriés avec install.packages("package_name")

Je n'ai pas eu moins de 2 minutes sur ma machine avec cette version (je dois essayer avec une modification parallèle)

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)
}

Appel et sortie:

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

Benchmark sur ma 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

Pour information, pour une plage de 1:11, cela prend 285 secondes.

PARI / GP, n = 11

C'est de la force brute mais en profitant det(A^T) = det(A). Je ne le poste que pour montrer à quel point il est facile de sauter les transpositions. Le bit le plus bas de b1contient la cellule supérieure gauche et les autres bits contiennent le reste de la ligne supérieure. b2contient le reste de la colonne de gauche. Nous appliquons simplement 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);
  )
}

En ce qui concerne le calcul des déterminants de Toeplitz dans le O(n^2)temps: dans mes recherches limitées, j'ai continué à me heurter à une exigence selon laquelle tous les principaux mineurs principaux doivent être différents de zéro pour que les algorithmes fonctionnent, ce qui est un obstacle majeur à cette tâche. N'hésitez pas à me donner des conseils si vous en savez plus que moi.