Znajdź maksymalną determinantę dla każdego rozmiaru macierzy Toeplitz

14

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 ++

źródło
@AlexA. Masz rację, a ja przeprosiłam. Na szczęście oba problemy są bardzo podobne, więc powinien łatwo móc modyfikować swój kod.
Rozwiązanie @Vioz zawiera sekwencję rozpoczynającą się od 1, 1, 2, 3, 5, 9, 32. Zatem wartość dla n = 5 jest inna niż na liście. Ponieważ wszystkie pozostałe wartości są zgodne, wygląda na to, że rozwiązanie jest prawdopodobnie poprawne, a to tylko literówka w pytaniu?
Reto Koradi,
@RetoKoradi Dziękujemy. Naprawiony.
Oto 10 możliwych binarnych macierzy Toeplitz z maksymalnymi determinantami dla n = 1..10: ghostbin.com/paste/axkpa
Tyilo,
2
Jako obserwacja, która może pomóc innym, ale nie mogę zweryfikować powyżej 14. Wydaje się, że odpowiednie środki górnego rzędu i pierwszej kolumny macierzy Toeplitza wynoszą zawsze 0,4 <= m <= 0,6 dla wyznacznika maksymalnego.
MickyT,

Odpowiedzi:

3

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;
}
Reto Koradi
źródło
Istnieje ciekawy sposób obliczenia wyznacznika macierzy liczb całkowitych przy użyciu tylko arytmetyki liczb całkowitych: dekompozycja LU w pewnym polu skończonym (w zasadzie mod dużej liczby pierwszej). Nie wiem czy to by było szybsze.
lirtosiast
@ThomasKwa To prawdopodobnie nadal byłby O (n ^ 3)? Może to być pomocne w przypadku większych matryc, w których precyzja zmiennoprzecinkowa stałaby się problemem. Tak naprawdę nie szukałem literatury. Cóż, przeprowadziłem szybkie wyszukiwanie i znalazłem artykuł na temat obliczania wyznaczników macierzy Toeplitza. Ale było zbyt wiele otwartych pytań, aby poświęcić czas na ich wdrożenie.
Reto Koradi,
1
@Lembik Spróbuję to dzisiaj obejrzeć. Zmieniłem go, aby obsługiwać większe rozmiary dla twojego innego pokrewnego wyzwania wczoraj. Jak dotąd nie udało mi się pokonać najwyższych wyników dla n = 30, moja heurystyka utknęła poniżej 5 * 10 ^ 13.
Reto Koradi,
1
@Lembik Patrz paste.ubuntu.com/11915546 kodu i paste.ubuntu.com/11915532 na wyniki do N = 19.
Reto Koradi,
1
@Lembik Wyniki do n = 20 znajdują się na stronie paste.ubuntu.com/11949738 . Wyświetlają teraz listę wszystkich powiązanych rozwiązań, w tym atrybuty, aby szybko zobaczyć wartość przekątnej i czy są w obiegu. Wszystkie maksima dla m = 18,19,20 są matrycami krążącymi. Proszę dokładnie sprawdzić determinanty przed opublikowaniem ich w dowolnym miejscu.
Reto Koradi
8

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.

Legendre
źródło
1
Uwaga dla testerów: n = 12wymaga około 18 GiB pamięci.
Dennis
To bardzo fajna poprawa. Wyjście jest jednak nieco wadliwe. Dla mnie, używając j64-804, wypisuje n = 1 dwa razy, więc jest jeszcze jeden.
@Lembik Ah to prawda. Właśnie zaktualizowałem kod; mógłbyś spróbować ponownie? Dzięki! (Ustawiłem go na obliczanie do n=13. Możesz zmienić 13w przedostatnim wierszu, aby obliczyć, co chcesz.)
Legendre
Uruchomiłem go ponownie i nadal osiąga 12.
@Lembik Hmm .. czy mówisz, że dotrze do 12 w wyznaczonym terminie i do 13 po jakimś czasie (czego się spodziewam), czy też nigdy nie osiąga 13 (tzn. Program zatrzymuje się po 12)?
Legendre
4

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
Kade
źródło
Dziękuję, ale myślę, że odpowiedziałeś na starą wersję pytania, obawiam się. Teraz chodzi o matryce Toeplitz, a limit czasu wynosi 2 minuty.
4
Na tej stronie widzę tak bardzo golfowego Pythona, że ​​często zapominam, że dla celów ogólnych jest to właściwie czytelny język.
Alex A.,
Prawdopodobnie można to znacznie przyspieszyć, ponieważ nie korzysta z faktu, że jest to macierz binarna.
lirtosiast
@ThomasKwa Jeśli mam być szczery, nie mam pojęcia, jak to wykorzystać: P
Kade
Cytat z dokumentacji Numpy: „Wyznacznik jest obliczany przez faktoryzację LU przy użyciu procedury LAPACK z / dgetrf.” Spojrzałem na dgetrf i mówi, że używa podwójnej precyzji; w zależności od GPU OP pojedyncza precyzja może być szybsza.
lirtosiast
4

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);
  }
}
SteelRaven
źródło
Wygląda na to, że wkrótce będziesz robił swój pierwszy wpis OEIS, chyba że ktoś
2

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;
}
Tyilo
źródło
Wygląda na to, że obliczasz wyznacznik za pomocą ekspansji przez nieletnich; ma to O(n!)złożoność, więc lepiej użyć innego algorytmu.
lirtosiast
@ThomasKwa Nie wiedziałem, że istnieją szybsze algorytmy, więc tak, to rozwiązanie jest dość złe.
Tyilo,
Możesz zajrzeć za pomocą dekompozycji LU, aby znaleźć wyznacznik macierzy. To O(n^3), jak sądzę, choć mogą być wykonane szybciej ciekawych algorytmów. Wierzę, że większość wbudowanych tu wbudowanych generalnie używa wariantu dekompozycji do wykonania wyznaczników.
BrainSteel
@BrainSteel, tak, spojrzałem na to, ale równie dobrze mogę skorzystać z O(n^2)algorytmu, jeśli aktualizuję swoją odpowiedź.
Tyilo,
Według przypadkowego wyszukiwania w Wikipedii wyznacznik macierzy Toeplitza można określić w O(n^2). Ale myślę, że wąskim gardłem problemu jest wyszukiwanie wśród O(4^n)wielu 0-1 nwedług nmacierzy.
Legendre,
2

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.

Tensibai
źródło
1

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.

Mitch Schwartz
źródło
Widziałeś ten papier? scienpress.com/upload/JAMB/Vol%201_1_4.pdf . Nie było dla mnie jasne, jaka jest złożoność. Wydawało się, że istnieje wiele terminów dla przykładu n = 5.
Reto Koradi,
@RetoKoradi Tak, widziałem to. Wydaje się, że złożoność nie jest wielomianowa, biorąc pod uwagę, że np. e_{k+1}Ma 4-krotność liczby składników jak e_k. W dokumencie jest wiele pominięć. Odwracalna matryca ma rozkład LU, jeśli wszystkie wiodące nieletnie są niezerowe. (Zwróć uwagę na mianowniki, np. a_0- domyślnie gwarantuje się, że są niezerowe). Wyjątkowość wynika z tego, że L jest jednostką trójkątną. Autor nie wspomniał także o stabilności numerycznej. W przypadku, gdy łącze stanie się niedostępne, artykuł brzmi „O obliczaniu determinantów matematycznych Toeplitza” Hsuan-Chu Li (2011).
Mitch Schwartz,