add dynamic programming algorithm for high order ANOVA

Author: vene

polylearn/kernels.py

@@ -6,6 +6,9 @@
 from scipy.sparse import issparse
 
 import numpy as np
+from lightning.impl.dataset_fast import get_dataset
+
+from .kernels_fast import _fast_anova_kernel_batch, _fast_anova_grad
 
 
 def safe_power(X, degree=2):
@@ -65,7 +68,7 @@ def homogeneous_kernel(X, P, degree=2):
     return polynomial_kernel(X, P, degree=degree, gamma=1, coef0=0)
 
 
-def anova_kernel(X, P, degree=2):
+def anova_kernel(X, P, degree=2, method='auto'):
     """ANOVA kernel between X and P::
 
         K_A(x, p) = sum_i1>i2>...>id x_i1 p_i1 x_i2 p_i2 ... x_id p_id
@@ -81,11 +84,25 @@ def anova_kernel(X, P, degree=2):
 
     degree : int, default 2
 
+    method : string, default: 'auto'
+        - 'dp' :  dynamic programming recursion
+        - 'auto':  vectorized formula for degree 2 or 3, revert to 'dp' for
+            higher degrees.
+
     Returns
     -------
     Gram matrix : array of shape (n_samples_1, n_samples_2)
     """
-    if degree == 2:
+
+    if degree > 3 or method == 'dp':
+        n_samples = X.shape[0]
+        n_components = P.shape[0]
+        ds = get_dataset(X, 'c')
+
+        K = np.empty((n_samples, n_components))
+        _fast_anova_kernel_batch(ds, P, degree, K)
+
+    elif degree == 2:
         K = homogeneous_kernel(X, P, degree=2)
         K -= _D(X, P, degree=2)
         K /= 2
@@ -95,11 +112,20 @@ def anova_kernel(X, P, degree=2):
         K += 2 * _D(X, P, degree=3)
         K /= 6
     else:
-        raise NotImplementedError("ANOVA kernel for degree >= 4 not yet "
-                                  "implemented efficiently.")
+        raise ValueError("Unsupported parameters. Degree must be > 1.")
+
     return K
 
 
+def anova_grad(X, i, P, degree=2):
+    """Computes the ANOVA gradient of the i-th row of X, wrt to P"""
+
+    ds = get_dataset(X, 'c')
+    grad = np.empty_like(P)
+    _fast_anova_grad(ds, i, P, degree, grad)
+    return grad
+
+
 def _poly_predict(X, P, lams, kernel, degree=2):
     if kernel == "anova":
         K = anova_kernel(X, P, degree)