Cálculo de la suma $\frac{1}{2} \sum x^T \Sigma x$ para todos $x \in \{0,1\}^n$

Question

Nota: la ecuación dentro de la suma está relacionada con las Máquinas de Boltzmann / Redes de Hopfield, la función de energía de estas funciones son similares. Para más información, por ejemplo, sobre cómo obtener el estimador de máxima verosimilitud para $\Sigma$ podrías echar un vistazo al libro de David MacKay, sección 43 MacKay lo llama $W$ )

http://www.inference.phy.cam.ac.uk/itprnn/ps/521.526.pdf

Pregunta original:

Intento calcular la suma

$$ \frac{1}{2} \sum_{\forall x} x^T \Sigma x $$

donde $\Sigma$ es una matriz simétrica definida positiva, $x$ es un vector columna con valores binarios. Puedo asumir la normalización de $\Sigma$ como otro paso (para que todas las filas sumen 1), si eso me da convergencia. $\Sigma$ es el resultado del cálculo de la coocurrencia $A^TA$ sobre datos de muestra que sólo contienen valores binarios, 0 ó 1.

He realizado algunos experimentos numéricos con muestras generadas aleatoriamente,

from sklearn.preprocessing import normalize
import numpy as np
np.random.seed(0)
A = np.random.randint(0,2,(100,4))
cooc = A.T.dot(A).astype(float)
cc = normalize(cooc, norm='l1', axis=1) 

Cuando aplico la fórmula, es decir, recorro todos los valores posibles de $x \in \{0,1\}^n$ y calcula la suma,

import itertools, numpy.linalg as lin
lst = np.array(list(itertools.product([0, 1], repeat=A.shape[1])))
s = 0
for x in lst: 
    s += np.dot(np.dot(x.T,cc), x) / 2
print s

Obtengo 11,15 en datos con 4 dimensiones. Para n=5 obtengo alrededor de 26,6, en cualquier muestra que genere con estas dimensiones. A partir de esto he llegado a la conclusión de que este número está directamente relacionado con la dimensión con la que estoy trabajando y es una constante, por lo que me preguntaba si podría calcularse como un límite de alguna manera.

Nota:

Necesito este número como una constante de normalización, planeo usar $p(x;\Sigma) = \frac{1}{2C} x^T \Sigma x$ como función de masa de probabilidad. He aquí cómo he llegado a esta pregunta.

Intentaba capturar la frecuencia de cada variable individual y la dependencia entre todas las variables de muestras binarias multivariantes. Utilicé un simple cálculo de coocurrencia en la muestra,

import numpy as np
A = np.array([\
[0.,1.,1.,0],
[1.,1.,0, 0],
[1.,1.,1.,0],
[0, 1.,1.,1.],
[0, 0, 1.,0]
])
c = A.T.dot(A).astype(float)
print c 

Resultado

[[ 2.  2.  1.  0.]
 [ 2.  4.  3.  1.]
 [ 1.  3.  4.  1.]
 [ 0.  1.  1.  1.]]

Ahora, para cualquier nuevo punto de datos $x$ si quisiera calcular una "puntuación", ej:

x = np.array([[0,1,1,0]])
print np.dot(np.dot(x.T,c), x) / 2

me daría 7. La fórmula básicamente escoge los números 4,3,3,4 en el bloque central de la matriz cooc, y los suma, que era lo que yo quería porque el nuevo punto de datos $x=[0,1,1,0]$ tiene activadas las variables binarias 2 y 3 (es 1) por lo que me interesa la dependencia entre $x_2$ y $x_3$ así como la frecuencia de las variables por sí mismas.

Una vez que tuve esta puntuación, empecé a preguntarme si podría convertir esta función en una PMF, de ahí la necesidad de la constante de normalización y la necesidad de integrar la función para todos los valores posibles de $x$ .

He jugado con la idea de dividir la suma por $x^Tx$ lo que hace que la ecuación se parezca al Cociente de Rayleigh,

$$ = \frac{1}{2} \sum_{\forall x} \frac{x^T \Sigma x }{x^Tx} $$

entonces, si supusiera "x=todos los valores propios" en lugar de "x=todos los valores posibles", tal vez sumando todos los valores propios obtendría algo. Pero la suma debe ser para todas las x. ¿Representar todas las x usando los vectores propios como base quizás..?

Answer 1

Sea $X_{1}, \cdots, X_{n}$ sean variables aleatorias i.i.d. Bernoulli. Entonces con $\Sigma = (\sigma_{ij})$ ,

\begin{align*} \frac{1}{2} \sum_{\mathrm{x} \in \{0, 1\}^{n}} \mathrm{x}^{T} \Sigma \mathrm{x} &= 2^{n-1}\Bbb{E} \sum_{i,j} \sigma_{ij}X_{i}X_{j} \\ &= 2^{n-1}\sum_{i,j} \sigma_{ij} \Bbb{E} (X_{i}X_{j}) \\ &= 2^{n-1}\bigg( \sum_{i} \sigma_{ii} \Bbb{E} (X_{i}) + \sum_{i \neq j} \sigma_{ij} \Bbb{E} (X_{i}X_{j}) \bigg) \\ &= 2^{n-3}\bigg( \operatorname{tr} \Sigma + \sum_{i, j} \sigma_{ij} \bigg). \end{align*}

Asumiendo el paso de normalización para que $\Sigma$ es una matriz estocástica derecha, entonces la suma $\sum_{i, j} \sigma_{ij}$ se reduce a $n$ y por lo tanto

$$ \frac{1}{2} \sum_{\mathrm{x} \in \{0, 1\}^{n}} \mathrm{x}^{T} \Sigma \mathrm{x} = 2^{n-3}(n + \operatorname{tr}\Sigma). $$

La expectativa de esta cantidad sobre todas las posibles $\Sigma$ depende de la distribución de $A$ pero siempre tenemos $1 \leq \operatorname{tr}\Sigma \leq n$ por lo que el valor $(n+1) 2^{n-3}$ se aproxima bien a este valor hasta la magnitud.

Answer 2

Basándonos en la respuesta de @sos404, la ecuación final es

$$ p(x;\Sigma,n) = \frac{ \mathrm{x}^{T} \Sigma \mathrm{x}}{2^{n-2}(n + tr \ \Sigma)} $$

Algunas comprobaciones numéricas

from sklearn.preprocessing import normalize
import itertools
n = 4
A = np.random.randint(0,2,(200,n))
c = A.T.dot(A).astype(float)
cc = normalize(c, norm='l1', axis=1) 
lst = np.array(list(itertools.product([0, 1], repeat=n)))
res = [np.dot(np.dot(x.T,cc), x) / \
         ((2**(n-2)) * (n+np.sum(np.diag(cc)))) for x in lst]
print np.array(res).sum()

Independientemente del valor elegido para D, el resultado es siempre 1,0.

En itertools.product es la llamada a la función que permite iterar sobre todos los valores posibles de $x \in \{0,1\}^n$ .