-added code to compute only components of gradients/hessians and tests

-use in nystrom

kernel_exp_family/estimators/full/develop/gaussian_nystrom.py

@@ -50,9 +50,9 @@ def grad_naive(x, X, sigma, alpha, beta, inds):
 
     for ind, (a,i) in enumerate(ais):
         x_a = X[a]
-        xi_gradient_vec = gaussian_kernel_dx_i_dx_i_dx_j(x, x_a, sigma)
+        xi_gradient_mat = gaussian_kernel_dx_i_dx_i_dx_j(x, x_a, sigma)
         left_arg_hessian = gaussian_kernel_dx_i_dx_j(x, x_a, sigma)
-        xi_grad += xi_gradient_vec[i] / N
+        xi_grad += xi_gradient_mat[i] / N
         betasum_grad += beta[ind] * left_arg_hessian[i]
 
     return alpha * xi_grad + betasum_grad

kernel_exp_family/estimators/full/gaussian_nystrom.py

@@ -1,6 +1,9 @@
 from kernel_exp_family.estimators.estimator_oop import EstimatorBase
-from kernel_exp_family.estimators.full.develop.gaussian_nystrom import log_pdf_naive,\
-    nystrom_naive, grad_naive
+from kernel_exp_family.estimators.full.develop.gaussian_nystrom import nystrom_naive, grad_naive, ind_to_ai
+from kernel_exp_family.kernels.kernels import gaussian_kernel_dx_dx, gaussian_kernel_dx_component,\
+    gaussian_kernel_dx_dx_component, gaussian_kernel_dx_i_dx_i_dx_j,\
+    gaussian_kernel_dx_i_dx_j, gaussian_kernel_dx_i_dx_i_dx_j_component,\
+    gaussian_kernel_dx_i_dx_j_component
 from kernel_exp_family.tools.assertions import assert_array_shape
 import numpy as np
 
@@ -20,11 +23,39 @@ def fit(X, sigma, lmbda, inds):
     return alpha, beta
 
 def log_pdf(x, X, sigma, alpha, beta, inds):
-    return log_pdf_naive(x, X, sigma, alpha, beta, inds)
+    N, D = X.shape
+
+    xi = 0
+    betasum = 0
+
+    ais = [ind_to_ai(ind, D) for ind in range(len(inds))]
+
+    for ind, (a, i) in enumerate(ais):
+        gradient_x_xa_i = gaussian_kernel_dx_component(x, X[a, :], i, sigma)
+        xi_grad_i = gaussian_kernel_dx_dx_component(x, X[a,:], i, sigma)
+
+        xi += xi_grad_i / N
+        betasum += gradient_x_xa_i * beta[ind]
+
+    return np.float(alpha * xi + betasum)
 
 def grad(x, X, sigma, alpha, beta, inds):
-    return grad_naive(x, X, sigma, alpha, beta, inds)
+    N, D = X.shape
 
+    xi_grad = 0
+    betasum_grad = 0
+
+    ais = [ind_to_ai(ind, D) for ind in range(len(inds))]
+
+    for ind, (a,i) in enumerate(ais):
+        x_a = X[a]
+        xi_gradient_mat_component = gaussian_kernel_dx_i_dx_i_dx_j_component(x, x_a, i, sigma)
+        left_arg_hessian_component = gaussian_kernel_dx_i_dx_j_component(x, x_a, i, sigma)
+
+        xi_grad += xi_gradient_mat_component / N
+        betasum_grad += beta[ind] * left_arg_hessian_component
+
+    return alpha * xi_grad + betasum_grad
 
 class KernelExpFullNystromGaussian(EstimatorBase):
     def __init__(self, sigma, lmbda, D, N, m):
@@ -38,7 +69,7 @@ def __init__(self, sigma, lmbda, D, N, m):
         self.beta = np.zeros(m)
         self.X = np.zeros((0, D))
 
-        self.inds = np.sort(np.random.permutation(N*D)[:m])
+        self.inds = np.sort(np.random.permutation(N * D)[:m])
         self.m = m
 
     def fit(self, X):

kernel_exp_family/kernels/kernels.py

@@ -100,6 +100,19 @@ def gaussian_kernel_grad(x, Y, sigma=1.):
     G = (2.0 / sigma) * (k.T * differences)
     return G
 
+def gaussian_kernel_dx_component(x, y, ell, sigma=1.):
+    """
+    Ell'th partial derivative wrt left argument
+    """
+    assert(len(x.shape) == 1)
+    assert(len(y.shape) == 1)
+    assert(len(x) == len(y))
+
+    k = gaussian_kernel(np.atleast_2d(x), np.atleast_2d(y), sigma)[0,0]
+    difference = y[ell] - x[ell]
+    G = (2.0 / sigma) * (k * difference)
+    return G
+
 def gaussian_kernel_hessian(x, y, sigma=1.):
     assert(len(x.shape) == 1)
     assert(len(y.shape) == 1)
@@ -152,6 +165,14 @@ def gaussian_kernel_dx_dx(x, Y, sigma=1.):
     sq_differences = (Y - x)**2
     return k.T * (sq_differences*(2.0 / sigma)**2 - 2.0/sigma)
 
+def gaussian_kernel_dx_dx_component(x, y, ell, sigma=1.):
+    assert(len(x.shape) == 1)
+    assert(len(y.shape) == 1)
+    assert(len(x) == len(y))
+
+    k = gaussian_kernel(np.atleast_2d(x), np.atleast_2d(y), sigma)[0,0]
+    sq_difference = (y[ell] - x[ell])**2
+    return k.T * (sq_difference*(2.0 / sigma)**2 - 2.0/sigma)
 
 def gaussian_kernel_dx_dx_dy(x, y, sigma=1.):
     assert(len(x.shape) == 1)
@@ -201,6 +222,21 @@ def gaussian_kernel_dx_i_dx_j(x, y, sigma=1.):
 
     return term1 - term2
 
+def gaussian_kernel_dx_i_dx_j_component(x, y, ell, sigma=1.):
+    assert(len(x.shape) == 1)
+    assert(len(y.shape) == 1)
+    d = x.size
+
+    pairwise_dist = np.outer(y-x, y-x)
+
+    x_2d = x[np.newaxis,:]
+    y_2d = y[np.newaxis,:]
+
+    k = gaussian_kernel(x_2d, y_2d, sigma)
+    term1 = k*pairwise_dist * (2.0/sigma)**2
+    term2 = k*np.eye(d) * (2.0/sigma)
+
+    return (term1 - term2)[ell]
 
 def gaussian_kernel_dx_i_dx_i_dx_j(x, y, sigma=1.):
     """ Matrix of \frac{\partial k}{\partial x_i^2 \partial x_j}"""
@@ -221,6 +257,23 @@ def gaussian_kernel_dx_i_dx_i_dx_j(x, y, sigma=1.):
 
     return term1 - term2 - term3
 
+def gaussian_kernel_dx_i_dx_i_dx_j_component(x, y, ell, sigma=1.):
+    assert(len(x.shape) == 1), x
+    assert(len(y.shape) == 1)
+    d = x.size
+
+    pairwise_dist_squared_i = np.outer((y-x)**2, y-x)
+    row_repeated_distances = np.tile(y-x, [d,1])
+
+    x_2d = x[np.newaxis,:]
+    y_2d = y[np.newaxis,:]
+    k = gaussian_kernel(x_2d, y_2d, sigma)
+
+    term1 = k*pairwise_dist_squared_i * (2.0/sigma)**3
+    term2 = k*row_repeated_distances * (2.0/sigma)**2
+    term3 = term2*2*np.eye(d)
+
+    return (term1 - term2 - term3)[ell]
 
 def rff_sample_basis(D, m, sigma):
     # rbf sampler is parametrised in gamma, which is at the same time

tests/kernels/test_kernels.py

@@ -12,8 +12,9 @@
     rff_feature_map_grad2, rff_feature_map_grad_single, rff_sample_basis, \
     gaussian_kernel_hessian, gaussian_kernel_hessians, gaussian_kernel_dx_dx_dy, \
     gaussian_kernel_dx_dx, gaussian_kernel_dx_dx_dy_dy, gaussian_kernel_dx_i_dx_j, \
-    gaussian_kernel_dx_i_dx_i_dx_j
-
+    gaussian_kernel_dx_i_dx_i_dx_j, gaussian_kernel_dx_component,\
+    gaussian_kernel_dx_dx_component, gaussian_kernel_dx_i_dx_i_dx_j_component,\
+    gaussian_kernel_dx_i_dx_j_component
 import numpy as np
 
 
@@ -406,3 +407,47 @@ def test_gaussian_kernel_dx_i_dx_i_dx_j_equals_SE_dx_i_dx_i_dx_j():
     reference = SE_dx_i_dx_i_dx_j(x.reshape(-1, 1), y.reshape(-1, 1), l=np.sqrt(sigma/2.0))
 
     assert_close(implementation, reference)
+
+def test_gaussian_kernel_dx_component_equals_grad():
+    D = 4
+    x = np.random.randn(D)
+    y = np.random.randn(D)
+    sigma = 0.5
+
+    grad = gaussian_kernel_grad(x,np.atleast_2d(y), sigma)[0]
+    for i in range(D):
+        dxi = gaussian_kernel_dx_component(x, y, i, sigma)
+        assert_allclose(grad[i], dxi)
+
+def test_gaussian_kernel_dx_dx_component_equals_gaussian_kernel_dx_dx():
+    D = 4
+    x = np.random.randn(D)
+    y = np.random.randn(D)
+    sigma = 0.5
+
+    dx_dx = gaussian_kernel_dx_dx(x, np.atleast_2d(y), sigma)[0]
+    for i in range(D):
+        dxi = gaussian_kernel_dx_dx_component(x, y, i, sigma)
+        assert_allclose(dx_dx[i], dxi)
+
+def test_gaussian_kernel_dx_i_dx_i_dx_j_component_equals_gaussian_kernel_dx_i_dx_i_dx_j():
+    D = 4
+    x = np.random.randn(D)
+    y = np.random.randn(D)
+    sigma = 0.5
+
+    dx_i_dx_i_dx_j = gaussian_kernel_dx_i_dx_i_dx_j(x, y, sigma)
+    for i in range(D):
+        a = gaussian_kernel_dx_i_dx_i_dx_j_component(x, y, i, sigma)
+        assert_allclose(dx_i_dx_i_dx_j[i], a)
+
+def test_gaussian_kernel_dx_i_dx_j_component_equals_gaussian_kernel_dx_i_dx_j():
+    D = 4
+    x = np.random.randn(D)
+    y = np.random.randn(D)
+    sigma = 0.5
+
+    dx_i_dx_j = gaussian_kernel_dx_i_dx_j(x, y, sigma)
+    for i in range(D):
+        a = gaussian_kernel_dx_i_dx_j_component(x, y, i, sigma)
+        assert_allclose(dx_i_dx_j[i], a)