Encuentre el determinante máximo para cada matriz de Toeplitz de tamaño

14

Para un n fijo, considere las matrices de Toeplitz n por n con entradas que son 0 o 1. El objetivo es encontrar el determinante máximo sobre todas esas matrices de Toeplitz.

Tarea

Para cada uno nde 1 en adelante, genere el determinante máximo sobre todas las matrices de n por n Toeplitz con entradas que sean 0 o 1. Debe haber una salida por la ncual debe tener el determinante máximo y también una matriz de ejemplo que lo alcance.

Puntuación

Su puntaje es el más grande nque alcanza su código en 2 minutos en mi computadora. Para aclarar un poco, su código puede ejecutarse durante 2 minutos en total, esto no es 2 minutos por n.

Desempate

Si dos entradas obtienen el mismo npuntaje, la entrada ganadora será la que obtenga la mayor puntuación nen el menor tiempo posible en mi máquina. Si las dos mejores entradas son iguales en este criterio también, entonces el ganador será la respuesta presentada primero.

Idiomas y bibliotecas

Puede utilizar cualquier idioma y bibliotecas disponibles gratuitamente que desee. Debo poder ejecutar su código, así que incluya una explicación completa sobre cómo ejecutar / compilar su código en Linux si es posible.

Mi máquina Los tiempos se ejecutarán en mi máquina. Esta es una instalación estándar de ubuntu en un procesador AMD FX-8350 de ocho núcleos. Esto también significa que necesito poder ejecutar su código.

Respuestas pequeñas

Para n = 1..10 las salidas deben ser 1,1,2,3,5,9,32,56,125,315

Esta secuencia no está en OEIS, por lo que la entrada ganadora también puede proponer una nueva entrada allí.

Entradas hasta ahora

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

fuente
@AlexA. Tienes razón y me he disculpado. Afortunadamente, los dos problemas son muy similares, por lo que debería poder modificar fácilmente su código.
La solución de @Vioz surge con una secuencia que comienza con 1, 1, 2, 3, 5, 9, 32. Entonces, el valor para n = 5 es diferente de lo que usted enumera. Como todos los demás valores coinciden, parece que la solución probablemente sea correcta, ¿y esto es solo un error tipográfico en la pregunta?
Reto Koradi
@RetoKoradi Gracias. Fijo.
Aquí hay 10 posibles matrices binarias de Toeplitz con determinantes máximos para n = 1..10: ghostbin.com/paste/axkpa
Tyilo
2
Como observación que puede ayudar a otros pero no puedo verificar más allá de 14. Parece que los medios respectivos de la fila superior y la primera columna de la matriz de Toeplitz son siempre 0.4 <= m <= 0.6 para el determinante máximo.
MickyT

Respuestas:

3

C ++ con pthreads

Esto llega a n = 14 en poco menos de 1 minuto en mi máquina. Pero dado que es solo una computadora portátil de 2 núcleos, espero que la máquina de prueba de 8 núcleos pueda terminar n = 15 en menos de 2 minutos. Tarda unos 4:20 minutos en mi máquina.

Realmente esperaba encontrar algo más eficiente. No ha llegado a ser una manera de calcular el determinante de una matriz binaria de manera más eficiente. Quería crear algún tipo de enfoque de programación dinámica que cuente los términos +1 y -1 en el cálculo determinante. Pero simplemente no se ha unido hasta ahora.

Como la recompensa está por expirar, implementé el enfoque estándar de fuerza bruta:

  • Pase sobre todas las posibles matrices de Toeplitz.
  • Omita uno de los dos en cada par de matriz transpuesta. Dado que la matriz se describe mediante valores de máscara de bits, esto es fácil de hacer omitiendo todos los valores donde el reverso de la máscara de bits es más pequeño que la propia máscara de bits.
  • El determinado se calcula con un libro de texto LR descomposición. Excepto por algunos ajustes menores de rendimiento, la mejora principal que hice al algoritmo de mi libro de métodos numéricos de la universidad es que uso una estrategia de pivote más simple.
  • La paralelización se realiza con pthreads. El solo uso de un espaciado regular para los valores procesados ​​por cada subproceso causó un equilibrio de carga muy malo, por lo que introduje algunos cambios.

Probé esto en Mac OS, pero usé un código similar en Ubuntu antes, así que espero que esto se compile y se ejecute sin problemas:

  1. Guarde el código en un archivo con una .cppextensión, por ejemplo optim.cpp.
  2. Compilar con gcc -Ofast optim.cpp -lpthread -lstdc++.
  3. Corre con time ./a.out 14 8. El primer argumento es el máximo n. 14 debería terminar en menos de 2 minutos, pero sería genial si pudieras probar 15 también. El segundo argumento es el número de hilos. Usar el mismo valor que el número de núcleos de la máquina es normalmente un buen comienzo, pero probar algunas variaciones podría mejorar los tiempos.

Avíseme si tiene algún problema para crear o ejecutar el código.

#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
fuente
Hay una forma interesante de calcular el determinante de una matriz de enteros utilizando solo la aritmética de enteros: descomposición de LU en algún campo finito (básicamente mod un primo grande). No sé si esto sería más rápido.
lirtosiast
@ThomasKwa ¿Eso probablemente seguiría siendo O (n ^ 3)? Puede ser útil para matrices más grandes donde la precisión de coma flotante se convertiría en un problema. Realmente no busqué literatura. Bueno, hice una búsqueda rápida y encontré un artículo sobre el cálculo de los determinantes de las matrices de Toeplitz. Pero había demasiadas preguntas abiertas para dedicar tiempo a intentar implementarlo.
Reto Koradi
1
@Lembik Trataré de echarle un vistazo más tarde hoy. Lo cambié para manejar tamaños más grandes para su otro desafío relacionado ayer. No he podido superar los puntajes más altos para n = 30 hasta ahora, mis heurísticas están estancadas por debajo de 5 * 10 ^ 13.
Reto Koradi
1
@Lembik Consulte paste.ubuntu.com/11915546 para obtener el código y paste.ubuntu.com/11915532 para obtener resultados hasta n = 19.
Reto Koradi
1
@Lembik Los resultados hasta n = 20 están en paste.ubuntu.com/11949738 . Ahora enumeran todas las soluciones vinculadas, incluidos los atributos para ver rápidamente el valor de la diagonal y si son circulantes. Todos los máximos para m = 18,19,20 son matrices circulantes. Verifique los determinantes antes de publicarlos en cualquier lugar.
Reto Koradi
8

J

Actualización: código mejorado para buscar más de la mitad de los valores. Ahora calcula n=12cómodamente en 120 segundos (de 217 a 60 segundos).

Necesitará la última versión de J instalada.

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

Ejecute esto y mate cuando hayan transcurrido dos minutos. Mis resultados (MBP 2014 - 16 GB 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│
└────┴──┴─────────────┴────┘

Tiempo total de ejecución = 61.83s.


Solo por diversión

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

Esto tomó aproximadamente 210 segundos por sí solo.

Legendre
fuente
1
Nota para los probadores: n = 12requiere aproximadamente 18 GiB de memoria.
Dennis
Esta es una muy buena mejora. Sin embargo, la salida es ligeramente defectuosa. Para mí, usando j64-804, genera n = 1 dos veces, por lo que está fuera por uno para siempre.
@Lembik Ah, eso es correcto. Acabo de actualizar el código; ¿podrías intentar correr de nuevo? ¡Gracias! (Lo configuré para que calcule hasta n=13. Puede cambiar el 13en la penúltima línea para que calcule lo que desee.)
Legendre
Lo volví a ejecutar y todavía llega a 12.
@Lembik Hmm ... ¿estás diciendo que llega a 12 dentro del límite de tiempo y llega a 13 en algún momento después de eso (que es lo que espero), o que nunca llega a 13 (es decir, el programa se detiene después de 12)?
Legendre
4

Python 2

Esta es una solución muy sencilla, y probablemente no ganará el concurso. Pero oye, ¡funciona!

Daré un resumen rápido de lo que está sucediendo exactamente.

  1. Primero genero todas las filas iniciales posibles para n. Por ejemplo, cuándo n=2, esto generará una longitud de matriz 2 n + 1 , donde cada fila es longitud 2n-1. Se vería así: [[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. Luego, para cada una de esas posibles filas iniciales, lo giro nveces y corto los primeros nelementos para generar la matriz adecuada, y lo uso scipypara calcular el determinante, todo mientras hago un seguimiento del valor máximo. Al final de esto, simplemente imprimo el máximo, incremente nen 1, y continúo hasta que hayan pasado 10 minutos.

Para ejecutar esto, necesitará scipy instalado.

Edición 1: se modificó la forma en que se crearon las filas iniciales mediante el uso de itertools.product, ¡gracias Sp3000!

Edición 2: Se eliminó el almacenamiento de posibles filas iniciales para una mejora mínima en la velocidad.

Edición 3: cambiado para scipytener más control sobre cómo detfuncionó.

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

Aquí hay algunos resultados de muestra en mi máquina doméstica (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
Kade
fuente
Gracias pero creo que ha respondido una versión anterior de la pregunta, me temo. Ahora se trata de matrices Toeplitz y el límite de tiempo es de 2 minutos.
44
Veo tanto Python golfizado en este sitio que a menudo olvido que, para fines generales, en realidad es un lenguaje legible.
Alex A.
Esto probablemente podría acelerarse significativamente, porque no aprovecha el hecho de que es una matriz binaria.
lirtosiast
@ThomasKwa Si soy honesto, no tengo idea de cómo aprovechar eso: P
Kade
Cita de la documentación numpy: "El determinante se calcula mediante la factorización LU utilizando la rutina LAPACK z / dgetrf". Miré a dgetrf, y dice que usa doble precisión; dependiendo de la GPU de OP, la precisión individual podría ser más rápida.
lirtosiast
4

C ++

Fuerza bruta con el uso de OpenMP para paralelización y optimización simple para evitar la evaluación de determinantes para matrices transpuestas.

$ 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
fuente
Parece que pronto estarás haciendo tu primera entrada de OEIS a menos que a alguien se le
2

C

Compilar con:

$ clang -Ofast 52851.c -o 52851

Corre con:

$ ./52851

Puede generar el determinante máximo n = 1..10en ~ 115 segundos en mi computadora.

El programa solo está obteniendo el determinante de todas las posibles matrices binarias de tamaño de Toeplitz n, sin embargo, cada determinante de matrices de tamaño 5x5o menor se almacenará en caché mediante la memorización.

Al principio, supuse erróneamente que cada submatriz de una matriz de Toeplitz también sería una matriz de Toeplitz, por lo que solo necesitaba memorizar 2^(2n-1)valores en lugar de 2^(n^2)cada uno n. Hice el programa antes de darme cuenta de mi error, por lo que esta presentación es solo una solución de ese programa.


#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
fuente
Parece que estás calculando el determinante usando la expansión de menores; esto tiene O(n!)complejidad, por lo que es mejor que uses un algoritmo diferente.
lirtosiast
@ThomasKwa No sabía que había algoritmos más rápidos, así que sí, esta solución es bastante mala.
Tyilo
Es posible que desee considerar el uso de la descomposición LU para encontrar el determinante de una matriz. Es O(n^3), creo, aunque se puede hacer más rápido con algunos algoritmos interesantes. Creo que la mayoría de los componentes incorporados aquí generalmente usan una variante de descomposición para realizar determinantes.
BrainSteel
@BrainSteel, sí, lo miré, pero también podría buscar un O(n^2)algoritmo si estoy actualizando mi respuesta.
Tyilo
Según una búsqueda casual de Wikipedia, el determinante de una matriz de Toeplitz se puede determinar en O(n^2). Pero creo que el cuello de botella del problema está buscando entre los O(4^n)muchos 0-1 npor nmatrices.
Legendre
2

R

Tendrás que instalar R y los paquetes listados con install.packages("package_name")

No obtuve menos de 2 minutos en mi máquina con esta versión (tengo que probar con una modificación paralela)

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

Llamada y salida:

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

Punto de referencia en mi máquina:

> 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

Para información, para un rango de 1:11, toma 285 segundos.

Tensibai
fuente
1

PARI / GP, n = 11

Esta es la fuerza bruta pero que se aprovecha det(A^T) = det(A). Solo lo estoy publicando para demostrar lo fácil que es omitir las transposiciones. El bit más bajo b1contiene la celda superior izquierda, y los otros bits sostienen el resto de la fila superior. b2sostiene el resto de la columna izquierda. Simplemente hacemos cumplir 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);
  )
}

Con respecto al cálculo de los determinantes de Toeplitz a O(n^2)tiempo: en mi investigación limitada, me he topado con el requisito de que todos los principales principales menores de edad no sean cero para que los algoritmos funcionen, lo cual es un obstáculo importante para esta tarea. No dudes en darme consejos si sabes más sobre esto que yo.

Mitch Schwartz
fuente
¿Viste este artículo? scienpress.com/upload/JAMB/Vol%201_1_4.pdf . No estaba claro para mí cuál es la complejidad. Parecía haber bastantes términos para el ejemplo n = 5.
Reto Koradi
@RetoKoradi Sí, lo vi. Parece que la complejidad no es polinómica, considerando que, por ejemplo, e_{k+1}tiene 4 veces el número de componentes que e_k. Hay muchas omisiones en el documento. Una matriz invertible tiene una descomposición LU si todos los menores principales principales son distintos de cero. (Observe los denominadores, por ejemplo a_0, implícitamente se garantiza que no son cero.) La unicidad proviene de que L es unidad triangular. El autor tampoco mencionó la estabilidad numérica. En caso de que el enlace no esté disponible, el documento es "Sobre el cálculo de los determinantes de las matrices de Toeplitz" por Hsuan-Chu Li (2011).
Mitch Schwartz