mathematical - scipy optimize minimize in python

Una forma rápida de encontrar una respuesta de todo cero

Question

Para cada matriz de longitud n + h-1 con valores de 0 y 1, me gustaría comprobar si existe otra matriz de longitud n distinta de cero con valores de -1,0,1 para que todos los productos internos h cero. Mi manera ingenua de hacer esto es

import numpy as np
import itertools
(n,h)= 4,3
for longtuple in itertools.product([0,1], repeat = n+h-1):
bad = 0
for v in itertools.product([-1,0,1], repeat = n):
if not any(v):
continue
if (not np.correlate(v, longtuple, ''valid'').any()):
bad = 1
break
if (bad == 0):
print "Good"
print longtuple

Esto es muy lento si configuramos n = 19 y h = 10 que es lo que me gustaría probar.

Mi objetivo es encontrar una sola matriz "Buena" de longitud n+h-1 . ¿Hay alguna manera de acelerar esto para que n = 19 y h = 10 sea ​​factible?

El enfoque ingenuo actual toma 2 ^ (n + h-1) 3 ^ (n) iteraciones, cada una de las cuales toma aproximadamente n tiempo. Eso es 311,992,186,885,373,952 iteraciones para n = 19 y h = 10 cual es imposible.

Nota 1 Se modificó la convolve para correlate modo que el código considere la forma correcta.

10 de julio de 2015

El problema aún está abierto sin ninguna solución lo suficientemente rápida para n=19 y h=10 aún.

Answer 1

A continuación se muestra un algoritmo que reduce la complejidad de 3 ^ n a 3 ^ {nh}.

Deje que v_1, v_2, .., v_h sean los vectores a los que necesita ser ortogonal.

Considera el espacio vectorial (Z / 3Z) ^ n. Sea v''_1, .., v''_h las inclusiones naturales de v_1, .., v_h en este espacio.

Ahora sea w un vector con coeficientes en {-1,0,1}, y sea w ''el vector de (Z / 3Z) ^ n obtenido al ver naturalmente w como un vector de (Z / 3Z) ^ n. Entonces, una condición necesaria para que w tenga un producto escalar cero con v_1, .., v_h (en R) es que w ''tiene un producto escalar cero (en (Z / 3Z) ^ n) con v''_1, .., v ''_h.

Ahora puede determinar fácilmente el w ''que tiene cero producto escalar con v''_1, .., v''_h. Formarán un espacio de tamaño 3 ^ {nh}. Luego, debe verificar, para cada uno de ellos, si el w asociado fue en realidad ortogonal a todos los v_i.

Answer 2

Aquí hay un enfoque que lo reduce a O(n*h*3^(n/2 + 1)) . Eso escala mal, pero es lo suficientemente bueno para su caso de uso.

Iterar a través de todas las posibilidades para la primera mitad del vector. Cree un diccionario de diccionarios de ... de diccionarios de matrices cuyas claves son el valor de cada producto interno cambiado, y cuyo valor final es la matriz de las primeras mitades del vector que dio lugar a esa secuencia.

Ahora itere a través de todas las posibilidades para la segunda mitad del vector. A medida que calcula cada uno de sus productos internos, recorra el diccionario anidado y vea si hay las primeras mitades correspondientes cuya contribución al producto interno aún se cancela. Si recorres todo el camino hasta el final, puedes juntar la primera mitad que encontraste con la segunda mitad que también encontraste y tienes una respuesta.

¡No olvides ignorar la respuesta que es todo 0s!

Answer 3

Considere el siguiente enfoque de "reunirse en el medio".

Primero, replantea la situación en la formulación de matriz proporcionada por leekaiinthesky.

A continuación, tenga en cuenta que solo tenemos que considerar los vectores "cortos" s de la forma {0,1}^n (es decir, los vectores cortos que solo contienen 0 y 1) si cambiamos el problema para encontrar una matriz Hankel H de hxn de 0 y 1 es tal que Hs1 nunca es igual a Hs2 para dos vectores cortos diferentes de 0 y 1. Esto se debe a que Hs1 = Hs2 implica H(s1-s2)=0 que implica que hay un vector v de 1''s, 0''s y -1''s, es decir s1-s2 , tal que Hv = 0 ; a la inversa, si Hv = 0 para v en {-1,0,1}^n , entonces podemos encontrar s1 y s2 en {0,1}^n tal que v = s1 - s2 por Hs1 = Hs2 tanto, Hs1 = Hs2 .

Cuando n=19 hay solo 524,288 vectores s en {0,1}^n para probar; hash los resultados Hs y si el mismo resultado ocurre dos veces, H no es bueno y prueba con otra H En términos de memoria este enfoque es bastante factible. Hay 2^(n+h-1) matrices de Hankel H para probar; cuando n=19 y h=10 eso es 268,435,456 matrices. Eso es 2^38 pruebas, o 274,877,906,944, cada una con aproximadamente nh operaciones para multiplicar la matriz H y el vector s , aproximadamente 52 billones de operaciones. Eso parece factible, ¿no?

Dado que ahora solo está tratando con 0 y 1, no con -1, también podría acelerar el proceso mediante el uso de operaciones de bits (desplazamiento y contar 1).

Actualizar

Implementé mi idea en C ++. Estoy usando operaciones de bits para calcular productos de puntos, codificando el vector resultante como un entero largo, y utilizando unordered_set para detectar duplicados, obteniendo una salida temprana de un vector largo dado cuando se encuentra un vector duplicado de productos de puntos.

Obtuve 00000000010010111000100100 para n = 17 y h = 10 después de unos minutos, y 000000111011110001001101011 para n = 18 y h = 10 en un tiempo más. Estoy a punto de ejecutarlo para n = 19 y h = 10.

#include <iostream>
#include <bitset>
#include <unordered_set>
/* Count the number of 1 bits in 32 bit int x in 21 instructions.
* From /Hackers Delight/ by Henry S. Warren, Jr., 5-2
*/
int count1Bits(int x) {
x = x - ((x >> 1) & 0x55555555);
x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
x = (x + (x >> 4)) & 0x0F0F0F0F;
x = x + (x >> 8);
x = x + (x >> 16);
return x & 0x0000003F;
}
int main () {
const int n = 19;
const int h = 10;
std::unordered_set<long> dotProductSet;
// look at all 2^(n+h-1) possibilities for longVec
// upper n bits cannot be all 0 so we can start with 1 in pos h
for (int longVec = (1 << (h-1)); longVec < (1 << (n+h-1)); ++longVec) {
dotProductSet.clear();
bool good = true;
// now look at all n digit non-zero shortVecs
for (int shortVec = 1; shortVec < (1 << n); ++shortVec) {
// longVec dot products with shifted shortVecs generates h results
// each between 0 and n inclusive, can encode as h digit number in
// base n+1, up to (n+1)^h = 20^10 approx 13 digits, need long
long dotProduct = 0;
// encode h dot products of shifted shortVec with longVec
// as base n+1 integer
for(int startShort = 0; startShort < h; ++startShort) {
int shortVecShifted = shortVec << startShort;
dotProduct *= n+1;
dotProduct += count1Bits(longVec & shortVecShifted);
}
auto ret = dotProductSet.insert(dotProduct);
if (!ret.second) {
good = false;
break;
}
}
if (good) {
std::cout << std::bitset<(n+h-1)>(longVec) << std::endl;
break;
}
}
return 0;
}

Segunda actualización

El programa para n = 19 y h = 10 se ejecutó durante dos semanas en segundo plano en mi computadora portátil. Al final, acaba de salir sin imprimir ningún resultado. Salvo algún tipo de error en el programa, parece que no hay vectores largos con la propiedad que desea. Sugiero buscar razones teóricas por las que no existen vectores tan largos. Quizás algún tipo de argumento de conteo funcionará.

Answer 4

Esta es solo una respuesta parcial, ya que todavía parece demasiado lento para verificar el caso n=19, h=10 (y en este caso no es posible que existan vectores "buenos").

Aquí hay una implementación del algoritmo de verificación que describe @vib, y el uso de muestreo aleatorio en todos los vectores 2^(n+h-1) , en Mathematica.

TestFullNullSpace[A_, ns_] := Module[{dim, n},
{dim, n} = Dimensions[ns];
For[i = 1, i < 3^dim, i++,
testvec = Mod[IntegerDigits[i, 3, dim].ns, 3] /. {2 -> -1};
If[Norm[A.testvec] == 0,
Return[False]];
];
Return[True];
]
n = 17;
h = 10;
Do[
v = Table[RandomChoice[{0, 1}], {n + h - 1}];
A = Table[v[[i ;; i + n - 1]], {i, 1, h}];
ns = NullSpace[A, Modulus -> 3] /. {2 -> -1};
If[TestFullNullSpace[A, ns],
Print[v]];,
{1000}]

Salida de muestra para la ejecución anterior, después de unos segundos de cálculo:

{0,0,1,1,0,0,0,0,0,0,1,0,1,1,1,0,1,0,1,1,0,0,0,1,1,0}
{1,1,0,1,0,0,0,1,1,0,1,1,1,1,1,0,1,0,1,0,0,1,1,0,0,0}
{1,0,1,1,1,1,1,0,0,0,1,1,0,1,0,1,1,0,0,1,1,0,1,1,1,0}
{0,0,0,0,1,0,1,1,1,0,1,1,0,0,1,1,1,1,0,1,0,1,0,0,1,1}

Así que de 1000 vectores verificados, 4 fueron "buenos" (a menos que tenga un error). Desafortunadamente, para n=18 , corrí esto durante varios minutos y aún no encontré un vector "bueno". No sé si no existen o son extremadamente raras.

Answer 5

Puede haber una manera más rápida *

Lo que está buscando está relacionado con el concepto del núcleo o espacio nulo de una matriz .

En particular, para cada n + h-1 "longtuple" y dada n, construya una matriz h por n cuyas filas son las n subtítulos de la longtuple. En otras palabras, si su longuuple es [0,0,0,1,0,0] y n = 3, entonces su matriz es:

[[0 0 0]
[0 0 1]
[0 1 0]
[1 0 0]]

Llama a esta matriz A. Estás buscando un vector x tal que Ax = 0 , donde 0 es un vector de todos los 0s. Si existe tal x (que no es en sí misma todos los 0s) y se puede escalar para que contenga solo {-1, 0, 1}, entonces querrá tirar A y pasar al siguiente bloque largo.

No estoy muy seguro de cuál es la complejidad computacional (teóricamente más eficiente) de computar el kernel, pero parece ser del orden de O (h + n) ^ 3 o más, lo cual es, en cualquier caso, mucho mejor que O (3 ^ n). Consulte el enlace de Wikipedia anterior o Python (NumPy, SciPy), donde se encuentra el espacio nulo de una matriz para ver algunos ejemplos sobre cómo calcular el kernel.

De todos modos, una vez que identifiques el kernel, tendrás que hacer un trabajo adicional para averiguar si algún vector de la forma {-1, 0, 1} ^ n reside allí, pero no creo que sea tan grande. Una carga computacional.

* NB : En los comentarios, @vib señala que esto podría de hecho ser una gran carga computacional. No estoy seguro de cuál es el mejor algoritmo para averiguar si estos vectores se intersecan con el kernel. ¡Quizás no se pueda resolver en tiempo polinomial, en cuyo caso esta respuesta no acelera el problema original!

Código de ejemplo

Adaptar el código de la otra pregunta de desbordamiento de pila vinculada al ejemplo anterior que dio en los comentarios:

import scipy
from scipy import linalg, matrix
def null(A, eps=1e-15):
u, s, vh = scipy.linalg.svd(A)
null_mask = (s <= eps)
null_space = scipy.compress(null_mask, vh, axis=0)
return scipy.transpose(null_space)
A = matrix([[0,0,0,1],[0,0,1,0],[0,1,0,0],[0,0,0,0]])
print null(A)
#> [[-1.]
#> [ 0.]
#> [ 0.]
#> [ 0.]]

El código da un ejemplo (de hecho, el mismo ejemplo que usted dio) de una n-tupla que invalida [0, 0, 0, 1, 0, 0] como un "buen" conjunto de datos. Si el código devolvió [] , entonces presumiblemente no existe tal tupla n, y el tupla es "bueno". (Sin embargo, si el código devuelve algo, aún debe verificar la parte {-1, 0, 1}).

Pensamientos adicionales

Si existe tal x , ignorar la restricción {-1, 0, 1} por ahora, es equivalente a la pregunta de si la nulidad (dimensión del kernel) de A es mayor que 0. Esto sería equivalente a preguntar si el el rango de A es igual a n. Entonces, si encontró alguna manera de ser inteligente con la restricción {-1, 0, 1} y la dividió solo por la necesidad de calcular el rango de A , estoy seguro de que esto podría hacerse incluso más rápido.

Por cierto, parece muy probable que usted (o quien le haya dado este problema) ya sepa todo esto ... De lo contrario, ¿por qué habría llamado a la longitud del longtuple "n + h-1", si no hubiera Ya comencé con la matriz de altura h ...!