Dla stałej n rozważmy macierze Toeplitza n przez n z wpisami, które są albo 0 albo 1. Celem jest znalezienie maksymalnego wyznacznika we wszystkich takich macierzach Toeplitz.

Zadanie

Dla każdego nz 1 w górę wyprowadzaj maksymalną determinantę dla wszystkich n na n macierzy Toeplitza z wpisami, które są albo 0 albo 1. Powinno być jedno wyjście, na nktóre powinien mieć maksymalną determinantę, a także przykładową macierz, która ją osiąga.

Wynik

Twój wynik jest najwyższy w ntwoim kodzie w ciągu 2 minut na moim komputerze. Aby to trochę wyjaśnić, Twój kod może działać łącznie przez 2 minuty, nie jest to 2 minuty na n.

Łamacz krawatów

Jeśli dwa zgłoszenia otrzymają ten sam nwynik, zwycięskim wejściem będzie ten, który osiągnie najwyższy nw najkrótszym czasie na mojej maszynie. Jeśli dwa najlepsze zgłoszenia są równe również w tym kryterium, zwycięzcą będzie odpowiedź przesłana jako pierwsza.

Języki i biblioteki

Możesz korzystać z dowolnego darmowego języka i bibliotek, które ci się podobają. Muszę być w stanie uruchomić Twój kod, więc proszę podać pełne wyjaśnienie, jak uruchomić / skompilować kod w systemie Linux, jeśli to w ogóle możliwe.

Moja maszyna Czasy zostaną uruchomione na moim komputerze. Jest to standardowa instalacja ubuntu na ośmiordzeniowym procesorze AMD FX-8350. Oznacza to również, że muszę być w stanie uruchomić Twój kod.

Małe odpowiedzi

Dla n = 1..10 wyniki powinny wynosić 1,1,2,3,5,9,32,56,125,315

Ta sekwencja nie znajduje się w OEIS, więc zwycięski wpis może również zaproponować tam nowy wpis.

Dotychczasowe wpisy

• n=10 n=11autor: Vioz w Python
• n=9autor: Tyilo w C
• n=12autor: Legendre in J.
• n=10autor: Tensibai w R.
• n=14autor: SteelRaven w C ++
• n=14autor: RetoKoradi w C ++

C ++ z pthreads

To osiąga n = 14 w niecałą minutę na moim komputerze. Ale ponieważ jest to tylko 2-rdzeniowy laptop, mam nadzieję, że 8-rdzeniowa maszyna testowa może ukończyć n = 15 w niecałe 2 minuty. Na moim komputerze zajmuje to około 4:20 minut.

Naprawdę miałem nadzieję wymyślić coś bardziej wydajnego. Tam dostał się sposób bardziej efektywnie obliczyć zdeterminowanej macierzy binarnej. Chciałem wymyślić jakieś dynamiczne podejście programistyczne, które uwzględnia warunki +1 i -1 w obliczaniu wyznacznika. Ale do tej pory po prostu nie doszło do siebie.

Ponieważ nagroda wkrótce wygaśnie, wdrożyłem standardowe podejście brutalnej siły:

• Zapętlić wszystkie możliwe macierze Toeplitz.
• Pomiń jeden z dwóch w każdej transponowanej parze macierzy. Ponieważ macierz jest opisywana przez wartości maski bitowej, jest to proste, pomijając wszystkie wartości, w których odwrotność maski bitowej jest mniejsza niż sama maska ​​bitowa.
• Wyznaczenie jest obliczane na podstawie rozkładu LR w podręczniku. Oprócz drobnych poprawek wydajności, głównym ulepszeniem algorytmu z podręcznika metod numerycznych w college'u jest to, że korzystam z prostszej strategii przestawnej.
• Równoległość odbywa się za pomocą pthreads. Samo stosowanie regularnych odstępów dla wartości przetwarzanych przez każdy wątek spowodowało bardzo złe równoważenie obciążenia, więc wprowadziłem trochę swizzlingu.

Testowałem to na Mac OS, ale wcześniej użyłem podobnego kodu na Ubuntu, więc mam nadzieję, że to się skompiluje i uruchomi bez żadnych problemów:

1. Zapisz kod w pliku z .cpprozszerzeniem, np optim.cpp.
2. Kompiluj z gcc -Ofast optim.cpp -lpthread -lstdc++.
3. Uruchom z time ./a.out 14 8. Pierwszy argument to maksimum n. 14 powinno zakończyć się w mniej niż 2 minuty, ale byłoby wspaniale, gdybyś mógł spróbować również 15. Drugi argument to liczba wątków. Korzystanie z tej samej wartości co liczba rdzeni maszyny jest zwykle dobrym początkiem, ale wypróbowanie niektórych odmian może potencjalnie poprawić czasy.

Daj mi znać, jeśli masz problemy z budowaniem lub uruchomieniem kodu.

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

jot

Aktualizacja: Ulepszony kod do wyszukiwania ponad połowy wartości. Teraz oblicza n=12wygodnie w ciągu 120 sekund (z 217 do 60 sekund).

Trzeba będzie najnowsza wersja J zainstalowany.

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

Uruchom to i zabij, gdy miną dwie minuty. Moje wyniki (MBP 2014 - 16 GB pamięci 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│
└────┴──┴─────────────┴────┘

Całkowity czas pracy = 61,83s.

Dla żartu

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

Samo to zajęło około 210 sekund.

Python 2

To bardzo proste rozwiązanie i prawdopodobnie nie wygra konkursu. Ale hej, to działa!

Dam szybki przegląd tego, co się właściwie dzieje.

1. Najpierw generuję każdy możliwy wiersz początkowy dla n. Na przykład, kiedy n=2wygeneruje to tablicę o długości 2 ^{n + 1} , gdzie każdy wiersz ma długość 2n-1. To będzie wyglądać następująco: [[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. Następnie, dla każdego z tych możliwych początkowych wierszy, obracam go nrazy i odcinam pierwsze nelementy, aby wygenerować odpowiednią macierz, i używam scipydo obliczenia wyznacznika, wszystko z zachowaniem wartości maksymalnej. Na koniec po prostu drukuję maksimum, przyrost no 1 i kontynuuję, aż upłynie 10 minut.

Aby to uruchomić, musisz zainstalować scipy .

Edycja 1: Zmieniono sposób budowania początkowych wierszy za pomocą itertools.product, dzięki Sp3000!

Edycja 2: Usunięto przechowywanie możliwych rzędów początkowych dla minimalnej poprawy prędkości.

Edycja 3: Zmieniono na scipy mieć większą kontrolę nad sposobem detdziałania.

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

Oto kilka przykładowych danych wyjściowych na moim komputerze domowym (i7-4510U, 8 GB pamięci 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 z wykorzystaniem OpenMP do równoległości i prostej optymalizacji w celu uniknięcia oceny wyznacznika dla transponowanych matryc.

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

do

Połącz z:

$ clang -Ofast 52851.c -o 52851

Biegnij z:

$ ./52851

Może wyprowadzić maksymalną determinantę dla n = 1..10 w ciągu ~ 115 sekund na moim komputerze.

Program pobiera wyznacznik każdej możliwej binarnej macierzy Toeplitza wielkości n, jednak każdy wyznacznik macierzy wielkości5x5 lub mniejszej będzie buforowany przy użyciu zapamiętywania.

Na początku błędnie założyłem, że każda submatrix macierzy Toeplitz będzie również macierzą Toeplitz, więc musiałem tylko zapamiętać 2^(2n-1)wartości zamiast 2^(n^2)dla każdej n. Zrobiłem program, zanim zdałem sobie sprawę z mojego błędu, więc to zgłoszenie jest tylko naprawą tego programu.

#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

Musisz zainstalować R i pakiety wymienione w install.packages("package_name")

Ta wersja nie miała mniej niż 2 minuty na moim komputerze (muszę spróbować z równoległą modyfikacją)

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

Wywołanie i wyjście:

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

Benchmark na mojej maszynie:

> 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

Aby uzyskać informacje, dla zakresu 1:11 potrzeba 285 sekund.

PARI / GP, n = 11

To brutalna siła, ale wykorzystująca det(A^T) = det(A). Zamieszczam go tylko po to, aby pokazać, jak łatwo można pominąć transpozycję. Najniższy bit b1zawiera lewą górną komórkę, a pozostałe bity zawierają resztę górnego rzędu. b2zawiera resztę lewej kolumny. Po prostu egzekwujemy 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);
  )
}

Jeśli chodzi o obliczanie wyznaczników Toeplitza w O(n^2)czasie: w moich ograniczonych badaniach ciągle napotykałem wymóg, aby wszyscy główni główni nieletni byli niezerowi, aby algorytmy działały, co jest główną przeszkodą w tym zadaniu. Daj mi wskazówki, jeśli wiesz o tym więcej niż ja.