Untuk n tetap, pertimbangkan matriks n oleh n Toeplitz dengan entri yang 0 atau 1. Tujuannya adalah untuk menemukan penentu maksimum atas semua matriks Toeplitz tersebut.

Tugas

Untuk masing-masing ndari 1 ke atas, output penentu maksimum atas semua n oleh n Toeplitz matriks dengan entri yang 0 atau 1. Harus ada satu output per nyang harus memiliki penentu maksimum dan juga contoh matriks yang mencapainya.

Skor

Skor Anda adalah nkode Anda yang terbesar dalam 2 menit di komputer saya. Untuk sedikit memperjelas, kode Anda dapat berjalan selama 2 menit total, ini bukan 2 menit per n.

Tie breaker

Jika dua entri mendapatkan nskor yang sama maka entri yang menang akan menjadi yang tertinggi ndalam waktu singkat di komputer saya. Jika dua entri terbaik sama pada kriteria ini juga maka pemenangnya adalah jawaban yang dikirimkan terlebih dahulu.

Bahasa dan perpustakaan

Anda dapat menggunakan bahasa dan perpustakaan yang tersedia secara bebas yang Anda suka. Saya harus dapat menjalankan kode Anda jadi tolong sertakan penjelasan lengkap tentang cara menjalankan / kompilasi kode Anda di linux jika memungkinkan.

Mesin Saya Pengaturan waktu akan dijalankan pada mesin saya. Ini adalah instalasi ubuntu standar pada Prosesor Delapan Core AMD FX-8350. Ini juga berarti saya harus dapat menjalankan kode Anda.

Jawaban kecil

Untuk n = 1..10 output harus 1,1,2,3,5,9,32,56,125,315

Urutan ini tidak ada dalam OEIS dan entri yang menang juga akan mengusulkan entri baru di sana.

Entri sejauh ini

• n=10 n=11oleh Vioz in Python
• n=9oleh Tyilo di C
• n=12oleh Legendre in J
• n=10oleh Tensibai di R
• n=14oleh SteelRaven di C ++
• n=14oleh RetoKoradi di C ++

C ++ dengan pthreads

Ini mencapai n = 14 hanya dalam waktu kurang dari 1 menit pada mesin saya. Tetapi karena itu hanya laptop 2-inti, saya berharap bahwa mesin uji 8-inti dapat menyelesaikan n = 15 dalam waktu kurang dari 2 menit. Butuh sekitar 4:20 menit di mesin saya.

Saya benar-benar berharap untuk menghasilkan sesuatu yang lebih efisien. Sudah ada menjadi cara untuk menghitung penentu dari matriks biner lebih efisien. Saya ingin membuat semacam pendekatan pemrograman dinamis yang menghitung istilah +1 dan -1 dalam perhitungan determinan. Tapi sejauh ini belum cukup.

Karena hadiah akan segera kedaluwarsa, saya menerapkan pendekatan brute force standar:

• Simpulkan semua matriks Toeplitz yang mungkin.
• Lewati salah satu dari dua di setiap pasangan matriks yang dialihkan. Karena matriks dijelaskan oleh nilai bitmask, ini mudah dilakukan dengan melewatkan semua nilai di mana kebalikan dari bitmask lebih kecil dari bitmask itu sendiri.
• Yang ditentukan dihitung dengan dekomposisi LR buku teks. Kecuali untuk beberapa penyetelan kinerja kecil, peningkatan utama yang saya lakukan pada algoritma dari buku metode numerik kuliah saya adalah bahwa saya menggunakan strategi pivot yang lebih sederhana.
• Paralelisasi dilakukan dengan pthreads. Hanya menggunakan spasi reguler untuk nilai-nilai yang diproses oleh masing-masing thread menyebabkan penyeimbangan beban yang sangat buruk, jadi saya memperkenalkan beberapa swizzling.

Saya menguji ini di Mac OS, tapi saya menggunakan kode yang sama di Ubuntu sebelumnya, jadi saya harap ini akan dikompilasi dan dijalankan tanpa hambatan:

1. Simpan kode dalam file dengan .cppekstensi, mis optim.cpp.
2. Kompilasi dengan gcc -Ofast optim.cpp -lpthread -lstdc++.
3. Jalankan dengan time ./a.out 14 8. Argumen pertama adalah maksimum n. 14 harus selesai dalam waktu kurang dari 2 menit, tetapi akan lebih baik jika Anda bisa mencoba 15 juga. Argumen kedua adalah jumlah utas. Menggunakan nilai yang sama dengan jumlah inti mesin biasanya merupakan awal yang baik, tetapi mencoba beberapa variasi berpotensi meningkatkan waktu.

Beri tahu saya jika Anda memiliki masalah dalam membangun atau menjalankan kode.

#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

Pembaruan: Peningkatan kode untuk mencari lebih dari setengah nilai. Sekarang menghitung dengan n=12nyaman dalam 120 detik (turun dari 217s ke 60s).

Anda akan membutuhkan J versi terbaru yang diinstal.

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

Jalankan ini dan bunuh ketika dua menit sudah habis. Hasil saya (MBP 2014 - 16GB 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│
└────┴──┴─────────────┴────┘

Total waktu berjalan = 61.83d.

Hanya untuk bersenang-senang

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

Ini membutuhkan waktu sekitar 210 detik sendiri.

Python 2

Ini adalah solusi yang sangat mudah, dan mungkin tidak akan memenangkan kontes. Tapi hei, itu berhasil!

Saya akan memberikan gambaran singkat tentang apa yang sebenarnya terjadi.

1. Saya pertama kali menghasilkan setiap baris awal yang mungkin untuk n. Misalnya, ketika n=2, ini akan menghasilkan panjang array 2 ^{n + 1} , di mana setiap baris panjang 2n-1. Ini akan terlihat seperti ini: [[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. Kemudian, untuk masing-masing baris awal yang mungkin, saya memutarnya nkali dan memotong nitem pertama untuk menghasilkan matriks yang sesuai, dan gunakan scipyuntuk menghitung determinan, semua sambil melacak nilai maksimum. Pada akhir ini, saya cukup mencetak maksimum, selisih n1, dan terus berjalan hingga 10 menit berlalu.

Untuk menjalankan ini, Anda akan perlu berhemat menginstal .

Sunting 1: Mengubah bagaimana baris awal dibuat dengan menggunakan itertools.produk sebagai gantinya, terima kasih Sp3000!

Sunting 2: Penyimpanan yang dihapus dari baris-baris awal yang mungkin untuk peningkatan kecepatan minimal.

Sunting 3: Diubah scipyagar memiliki kontrol lebih besar atas cara detkerjanya.

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

Berikut ini beberapa contoh output pada mesin rumah saya (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

C ++

Bruteforce dengan penggunaan OpenMP untuk paralelisasi dan optimalisasi sederhana untuk menghindari evaluasi determinan untuk matriks yang ditransformasikan.

$ 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

Kompilasi dengan:

$ clang -Ofast 52851.c -o 52851

Jalankan dengan:

$ ./52851

Dapat menampilkan penentu maksimum n = 1..10dalam ~ 115 detik di komputer saya.

Program ini hanya mendapatkan determinan setiap kemungkinan ukuran matriks Toeplitz biner n , namun setiap determinan matriks ukuran 5x5atau lebih kecil akan di-cache menggunakan memoisasi.

Pada awalnya saya salah berasumsi bahwa setiap submatrix dari sebuah matriks Toeplitz juga akan menjadi sebuah matriks Toeplitz, jadi saya hanya perlu memotret 2^(2n-1)nilai daripada 2^(n^2)untuk masing-masing n. Saya membuat program sebelum menyadari kesalahan saya, jadi pengajuan ini hanyalah perbaikan dari program itu.

#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

Anda harus menginstal R dan paket yang terdaftar install.packages("package_name")

Tidak mendapatkan kurang dari 2 menit pada komputer saya dengan versi ini (Saya sudah mencoba dengan modifikasi paralel)

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

Panggilan dan keluaran:

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

Tolok ukur pada mesin saya:

> 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

Untuk informasi, untuk rentang 1:11, dibutuhkan 285 detik.

PARI / GP, n = 11

Ini adalah kekuatan kasar tetapi memanfaatkan det(A^T) = det(A). Saya hanya mempostingnya untuk menunjukkan betapa mudahnya melewati transpos. Bit terendah b1memegang sel kiri atas, dan bit lainnya memegang sisa baris atas. b2memegang sisa kolom kiri. Kami hanya menegakkan 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);
  )
}

Mengenai menghitung faktor-faktor penentu Toeplitz dalam O(n^2)waktu: Dalam penelitian terbatas saya, saya terus berlari ke persyaratan bahwa semua anak di bawah umur utama harus nol untuk agar algoritma bekerja, yang merupakan hambatan utama untuk tugas ini. Jangan ragu untuk memberi saya petunjuk jika Anda tahu lebih banyak tentang ini daripada saya.