calc_cn.py

from sklearn.datasets import load_svmlight_file
from numpy.random import random
from numpy.random import uniform
import numpy as np

K=40
data = load_svmlight_file('a9a.tr')
X = data[0].todense()
y = data[1]
U = uniform(0, 0.1/(K**0.5), size=(X.shape[1], K))
V = uniform(0, 0.1/(K**0.5), size=(X.shape[1], K))

H_U = np.zeros((X.shape[1]*K, X.shape[1]*K))
w = random((X.shape[1],1))
D = np.zeros(X.shape[0])
for i in range(X.shape[0]):
    print i
    y_i = y[i]
    x_i = X[i, :]
    q_i = np.dot(x_i, V)
    p_i = np.dot(x_i, V)
    val = np.inner(w.T, x_i) + np.inner(q_i, p_i)
    sigma = 1/(1+np.exp(-y_i*val))[0][0]
    D[i] = sigma*(1-sigma)
    q_i_x_i = np.outer(q_i, x_i).flatten()
    H_U += D[i] * np.outer(q_i_x_i, q_i_x_i)

H = np.dot(np.dot(X.T, np.diagflat(D)), X)
np.save('H_U', H_U)
np.save('H', H)

H = np.load('H.npy')
H_U = np.load('H_U.npy')
H += np.identity(H.shape[0])
H_U += np.identity(H_U.shape[0])


H_diag = np.diag(H)
H_U_diag = np.diag(H_U)

H_diag = H_diag ** -0.5
H_U_diag = H_U_diag ** -0.5

P_H = np.diagflat(H_diag)
P_H_U = np.diagflat(H_U_diag)



from numpy.linalg import eigh
w,v = eigh(H)
print abs(w[-1]/w[0])
w,v = eigh(H_U)
print abs(w[-1]/w[0])
H_precon = np.dot(np.dot(P_H, H), P_H.T)
w,v = eigh(H_precon)
print abs(w[-1]/w[0])
H_U_precon = np.dot(np.dot(P_H_U, H_U), P_H_U.T)
w,v = eigh(H_U_precon)
print abs(w[-1]/w[0])