Générer de manière algorithmique tous les points de grille à l'intérieur d'un hypercube

8

Le problème vient directement des mathématiques computationnelles et peut être a déclaré ce qui suit:

Étant donné une matrice régulière MRd×d , trouver efficacement tous les vecteurs vZd tels que Mv1 , où Mv est la composante maximale de la vecteur en module.

Je donne un algorithme ci-dessous, mais il est très inefficace, car il a beaucoup plus d'étapes que le nombre de points de grille. C'est toujours supportable pour d=5 dans mon cas, mais ça échoue complètement pour d=6 , je n'ai pas beaucoup de temps; de plus, j'aimerais également travailler sur des dimensions plus élevées.


Mon algorithme actuel est basé sur ce qui suit (avertissement: contient plus de mathématiques): Tout d'abord, calculez M1 . Ensuite, on observe que vM1MvM1 . Cela signifie que nous pouvons calculer L=floorM1 puis essayer tous les vecteurs v{L,L+1,,L1,L}d ; il y en a exactement (2L+1)d . (Et là réside le problème: si d=6 et L=18 , j'obtiens 2.5E9 itérations, ce qui est des ordres de grandeur plus grand que le nombre de points que je pense qu'il y a.)

yo '
la source
1
Avez-vous envisagé de formuler votre problème sous forme de programme linéaire entier, puis de compter les solutions? Je pense qu'il existe des méthodes disponibles dans les solveurs IP (par exemple CPLEX) pour le faire, mais je ne sais pas à quelle vitesse elles sont en pratique. Une méthode pour ce faire est connue sous le nom d'algorithme de Barvinok.
Cornelius Brand
@ C.Brand Merci, je vais y jeter un œil (cependant, ça fait longtemps que je n'ai vu aucun LP / IP pour la dernière fois). Par hasard, avez-vous une idée si cela est présent dans SageMath?
yo
Quant à vos préoccupations, la formulation originale du problème est déjà (presque) un programme entier; la seule chose dont vous devrez vous occuper est l'utilisation de la norme (non linéaire) -norm. Mais je n'ai aucune idée si un algorithme pour compter les solutions IP est disponible dans SageMath.
Cornelius Brand

Réponses:

3

Voici une autre façon de voir le problème: Vous avez un réseau engendré par les colonnes de . Utilisez l'algorithme Lenstra – Lenstra – Lovász (LLL) pour obtenir une base réduite de ce réseau. Si vous remplacez par une nouvelle matrice formée par la sortie de LLL, alors les colonnes de généreront toujours le même réseau, mais les vecteurs de base seront plus proches d'être orthogonaux entre eux, et les entrées de devrait avoir une plus petite ampleur.MMMM1

À partir de là, il serait également utile de lier chaque composant de séparément: c'est-à-dire que vous pouvez lier le ème composantpar. (Soit dit en passant, le n'est pas correct; nous devons utiliser la somme des éléments sur chaque ligne, pas le maximum.)vi|vi|j=1d|(M1)ij|vM1

Pour des valeurs de jusqu'à environ 30, l'algorithme LLL se terminera pratiquement instantanément. De manière asymptotique, cela prend , donc cela ralentira pour les très grands , mais en même temps le nombre de points que nous devons vérifier augmente de façon exponentielle en , donc le temps d'exécution du LLL n'est pas vraiment le goulot d'étranglement. En revanche, les économies sur le nombre de points à contrôler peuvent être énormes. J'ai écrit du code GAP pour générer une matrice aléatoire régulière (stochastique) et comparer les limites sur les composants dedO(d6)ddMv que nous obtenons en utilisant la base d'origine, par rapport à la base réduite en LLL (Soit dit en passant, nous n'avons pas besoin de supposer que la matrice est régulière; j'ai fait cette restriction uniquement parce que c'était le cas dans votre demande):

d: = 8;
M: = IdentityMat (d);
pour i dans [1..d] faire
  pour j dans [1..d] do
    M [i] [j]: = aléatoire ([- 10 ^ 8..10 ^ 8]);
  od;
  M [i]: = M [i] / Sum (M [i]);
od;
L: = LLLReducedBasis (M) .basis;
MM: = M ^ -1 * 1,0;
LL: = L ^ -1 * 1.0;
pour i dans [1..d] faire
  pour j dans [1..d] do
    MM [i] [j]: = MM [i] [j] * SignFloat (MM [i] [j]);
    LL [i] [j]: = LL [i] [j] * SignFloat (LL [i] [j]);
  od;
od;

Impression ("Limites de la base d'origine:");
uns: = [1..d] * 0 + 1;
v: = MM * uns;
pour i dans [1..d] faire
  v [i]: = Int (Floor (v [i]));
  Imprimer (v [i]);
  Impression(" ");
od;
Imprimer ("\ n (");
Imprimer (produit (v * 2 + 1));
Imprimer ("points à vérifier) ​​\ n");

Print ("Limites pour la base LLL:");
v: = LL * ceux;
pour i dans [1..d] faire
  v [i]: = Int (Floor (v [i]));
  Imprimer (v [i]);
  Impression(" ");
od;
Imprimer ("\ n (");
Imprimer (produit (v * 2 + 1));
Imprimer ("points à vérifier) ​​\ n");

La sortie suivante (basée sur la graine aléatoire par défaut, avec ) n'est pas atypique:d=8

Limites de la base d'origine: 9 23 24 4 23 16 23 4 
(258370076349 points à vérifier)
Limites pour la base de LLL: 3 3 2 2 3 4 2 3 
(2701125 points à vérifier)

Edit : Ce problème est un cas particulier du problème général d'énumération des points de réseau dans les polytopes convexes, qui s'avère être un problème bien étudié, et il existe des algorithmes plus efficaces que celui décrit ci-dessus. Voir cet article pour une enquête.

Brent Kerby
la source
J'ai pensé à obtenir une bonne base de la grille d'abord, mais je ne sais pas à quel point elle est efficace dans les dimensions élevées? (Remarque: le problème n'est actuellement pas intéressant pour moi puisque nous avons réussi à le contourner complètement, mais je vais essayer de garder un œil sur cette question.)
yo '26
J'ai ajouté un peu plus à ma réponse, avec quelques informations sur l'efficacité de LLL, ainsi que des exemples de code et de sortie.
Brent Kerby