Update scaling and initialization.

Regularization scaling is now ON by default. I think this is
sensible, because it keeps the choice independent of data split.

Adagrad seems very sensitive to the initial norm of P, so I changed
the init to have unit variance rather than 0.01.
Makes benchmark more reasonable but norms are still weird.
Finnicky tests (fm warm starts) had to be updated, but most
things behave well.

Author: vene

benchmarks/bench_20newsgroups.py

@@ -16,42 +16,43 @@
 
 
 estimators = {
-    'fm-2': FactorizationMachineClassifier(beta=0.0001,
-                                           n_components=30,
+    'fm-2': FactorizationMachineClassifier(beta=1e-15,
+                                           scale_regularization=True,
+                                           n_components=12,
                                            fit_linear=False,
                                            fit_lower=None,
                                            degree=2,
                                            random_state=0,
-                                           max_iter=20,
+                                           max_iter=100,
                                            tol=1e-10),
 
-    'polynet-2': PolynomialNetworkClassifier(beta=0.0001,
-                                             n_components=15,
+    'polynet-2': PolynomialNetworkClassifier(beta=1e-15,
+                                             n_components=6,
                                              degree=2,
                                              fit_lower=None,
-                                             max_iter=20,
+                                             max_iter=100,
                                              random_state=0,
                                              tol=1e-10)
 }
 
 estimators['fm-3'] = clone(estimators['fm-2']).set_params(degree=3)
-estimators['fm-2-ada'] = clone(estimators['fm-2']).set_params(
-    solver='adagrad')
-estimators['fm-3-ada'] = clone(estimators['fm-3']).set_params(
-    solver='adagrad')
+estimators['fm-2-ada'] = clone(estimators['fm-2']).set_params(solver='adagrad')
+estimators['fm-3-ada'] = clone(estimators['fm-3']).set_params(solver='adagrad')
 estimators['polynet-3'] = (clone(estimators['polynet-2'])
-                           .set_params(degree=3, n_components=10))
+                           .set_params(degree=3, n_components=4))
 
 if __name__ == '__main__':
     data_train = fetch_20newsgroups_vectorized(subset="train")
     data_test = fetch_20newsgroups_vectorized(subset="test")
-    X_train_csc = sp.csc_matrix(data_train.data)
-    X_test_csc = sp.csc_matrix(data_test.data)
-    X_train_csr = sp.csr_matrix(data_train.data)
-    X_test_csr = sp.csr_matrix(data_test.data)
+    train_mask = np.in1d(data_train.target, [0, 15])  # atheism vs christianity
+    test_mask = np.in1d(data_test.target, [0, 15])  # atheism vs christianity
+    X_train_csc = sp.csc_matrix(data_train.data[train_mask])
+    X_test_csc = sp.csc_matrix(data_test.data[test_mask])
+    X_train_csr = sp.csr_matrix(data_train.data[train_mask])
+    X_test_csr = sp.csr_matrix(data_test.data[test_mask])
 
-    y_train = data_train.target == 0  # atheism vs rest
-    y_test = data_test.target == 0
+    y_train = data_train.target[train_mask] == 0  # atheism is positive
+    y_test = data_test.target[test_mask] == 0
 
     print("20 newsgroups")
     print("=============")
@@ -60,6 +61,7 @@
     print("X_train density = {0}"
           "".format(X_train_csr.nnz / np.product(X_train_csr.shape)))
     print("y_train {0}".format(y_train.shape))
+    print("Training class ratio: {0}".format(y_train.mean()))
     print("X_test {0}".format(X_test_csr.shape))
     print("X_test.dtype = {0}".format(X_test_csr.dtype))
     print("y_test {0}".format(y_test.shape))
@@ -84,6 +86,10 @@
         accuracy[name] = accuracy_score(y_test, y_pred)
         f1[name] = f1_score(y_test, y_pred)
         print("done")
+        try:
+            print("||P|| =", np.linalg.norm(clf.P_))
+        except:
+            pass
 
     print("Classification performance:")
     print("===========================")

polylearn/factorization_machine.py

@@ -32,9 +32,10 @@ class _BaseFactorizationMachine(six.with_metaclass(ABCMeta, _BasePoly)):
     @abstractmethod
     def __init__(self, degree=2, loss='squared', n_components=2, alpha=1,
                  beta=1, tol=1e-6, fit_lower='explicit', fit_linear=True,
-                 learning_rate=0.001, solver='cd', warm_start=False,
-                 init_lambdas='ones', max_iter=10000, verbose=False,
-                 callback=None, n_calls=100, random_state=None):
+                 learning_rate=0.001, scale_regularization=True,
+                 solver='cd', warm_start=False, init_lambdas='ones',
+                 max_iter=10000, verbose=False, callback=None, n_calls=100,
+                 random_state=None):
         self.degree = degree
         self.loss = loss
         self.n_components = n_components
@@ -44,6 +45,7 @@ def __init__(self, degree=2, loss='squared', n_components=2, alpha=1,
         self.fit_lower = fit_lower
         self.fit_linear = fit_linear
         self.learning_rate = learning_rate
+        self.scale_regularization = scale_regularization
         self.solver = solver
         self.warm_start = warm_start
         self.init_lambdas = init_lambdas
@@ -82,10 +84,20 @@ def fit(self, X, y):
 
         X, y = self._check_X_y(X, y)
         X = self._augment(X)
-        n_features = X.shape[1]  # augmented
+        n_samples, n_features = X.shape  # augmented
         rng = check_random_state(self.random_state)
         loss_obj = self._get_loss(self.loss)
 
+        # Scale regularization params to make losses equivalent.
+        if self.scale_regularization and self.solver == 'cd':
+            alpha = 0.5 * self.alpha * n_samples
+            beta = 0.5 * self.beta * n_samples
+        elif not self.scale_regularization and self.solver == 'adagrad':
+            alpha = self.alpha / 0.5 * n_samples
+            beta = self.beta / 0.5 * n_samples
+        else:
+            alpha, beta = self.alpha, self.beta
+
         if not (self.warm_start and hasattr(self, 'w_')):
             self.w_ = np.zeros(n_features, dtype=np.double)
 
@@ -95,7 +107,7 @@ def fit(self, X, y):
             n_orders = 1
 
         if not (self.warm_start and hasattr(self, 'P_')):
-            self.P_ = 0.01 * rng.randn(n_orders, self.n_components, n_features)
+            self.P_ = rng.randn(n_orders, self.n_components, n_features)
         if 'ada' in self.solver:
             # ensure each slice P[0], P[1]... is in F-order
             self.P_ = np.transpose(self.P_, [1, 2, 0])
@@ -125,7 +137,7 @@ def fit(self, X, y):
 
             converged = _cd_direct_ho(self.P_, self.w_, dataset, X_col_norms,
                                       y, y_pred, self.lams_, self.degree,
-                                      self.alpha, self.beta, self.fit_linear,
+                                      alpha, beta, self.fit_linear,
                                       self.fit_lower == 'explicit', loss_obj,
                                       self.max_iter, self.tol, self.verbose)
             if not converged:
@@ -141,9 +153,9 @@ def fit(self, X, y):
 
             dataset = get_dataset(X, order="c")
             _fast_fm_adagrad(self, self.w_, self.P_[0], dataset, y,
-                             self.degree, self.alpha, self.beta,
-                             self.fit_linear, loss_obj, self.max_iter,
-                             self.learning_rate, self.callback, self.n_calls)
+                             self.degree, alpha, beta, self.fit_linear,
+                             loss_obj, self.max_iter, self.learning_rate,
+                             self.callback, self.n_calls)
         return self
 
     def _get_output(self, X):
@@ -212,9 +224,17 @@ class FactorizationMachineRegressor(_BaseFactorizationMachine,
         coordinate descent. If False, the model can still capture linear
         effects if ``fit_lower == 'augment'``.
 
-    learning_rate: double, default: 0.001
+    learning_rate : double, default: 0.001
         Learning rate for 'adagrad' solver. Ignored by other solvers.
 
+    scale_regularization : boolean, default: True
+        Whether to adjust regularization according to the number of samples.
+        This helps if, after tuning regularization, the model will be retrained
+        on more data.
+
+        If set, the loss optimized is mean_i(l_i) + 0.5 || params || ^2
+        If not set, the loss becomes sum_i(l_i) + || params || ^ 2
+
     solver : {'cd'|'adagrad'}, default: 'cd'
         - 'cd': Uses a coordinate descent solver. Currently limited to
         degree=3.
@@ -292,14 +312,15 @@ class FactorizationMachineRegressor(_BaseFactorizationMachine,
     """
     def __init__(self, degree=2, n_components=2, alpha=1, beta=1, tol=1e-6,
                  fit_lower='explicit', fit_linear=True, learning_rate=0.001,
-                 solver='cd', warm_start=False, init_lambdas='ones',
-                 max_iter=10000, verbose=False, callback=None, n_calls=100,
-                 random_state=None):
+                 scale_regularization=True, solver='cd', warm_start=False,
+                 init_lambdas='ones', max_iter=10000, verbose=False,
+                 callback=None, n_calls=100, random_state=None):
 
         super(FactorizationMachineRegressor, self).__init__(
             degree, 'squared', n_components, alpha, beta, tol, fit_lower,
-            fit_linear, learning_rate, solver, warm_start, init_lambdas,
-            max_iter, verbose, callback, n_calls, random_state)
+            fit_linear, learning_rate, scale_regularization, solver,
+            warm_start, init_lambdas, max_iter, verbose, callback, n_calls,
+            random_state)
 
 
 class FactorizationMachineClassifier(_BaseFactorizationMachine,
@@ -355,9 +376,17 @@ class FactorizationMachineClassifier(_BaseFactorizationMachine,
         coordinate descent. If False, the model can still capture linear
         effects if ``fit_lower == 'augment'``.
 
-    learning_rate: double, default: 0.001
+    learning_rate : double, default: 0.001
         Learning rate for 'adagrad' solver. Ignored by other solvers.
 
+    scale_regularization : boolean, default: True
+        Whether to adjust regularization according to the number of samples.
+        This helps if, after tuning regularization, the model will be retrained
+        on more data.
+
+        If set, the loss optimized is mean_i(l_i) + 0.5 || params || ^2
+        If not set, the loss becomes sum_i(l_i) + || params || ^ 2
+
     solver : {'cd'|'adagrad'}, default: 'cd'
         - 'cd': Uses a coordinate descent solver. Currently limited to
         degree=3.
@@ -436,11 +465,12 @@ class FactorizationMachineClassifier(_BaseFactorizationMachine,
 
     def __init__(self, degree=2, loss='squared_hinge', n_components=2, alpha=1,
                  beta=1, tol=1e-6, fit_lower='explicit', fit_linear=True,
-                 learning_rate=0.001, solver='cd', warm_start=False,
-                 init_lambdas='ones', max_iter=10000, verbose=False,
-                 callback=None, n_calls=100, random_state=None):
+                 learning_rate=0.001, scale_regularization=True, solver='cd',
+                 warm_start=False, init_lambdas='ones', max_iter=10000,
+                 verbose=False, callback=None, n_calls=100, random_state=None):
 
         super(FactorizationMachineClassifier, self).__init__(
             degree, loss, n_components, alpha, beta, tol, fit_lower,
-            fit_linear, learning_rate, solver, warm_start, init_lambdas,
-            max_iter, verbose, callback, n_calls, random_state)
+            fit_linear, learning_rate, scale_regularization, solver,
+            warm_start, init_lambdas, max_iter, verbose, callback, n_calls,
+            random_state)

polylearn/tests/test_adagrad.py

@@ -12,13 +12,17 @@
 from .test_kernels import dumb_anova_grad
 
 
-def sg_adagrad_slow(P, X, y, degree, beta, max_iter, learning_rate):
+def sg_adagrad_slow(P, X, y, degree, beta, max_iter, learning_rate,
+                    scale_regularization=True):
 
     n_samples = X.shape[0]
     n_components = P.shape[0]
 
     grad_norms = np.zeros_like(P)
 
+    if not scale_regularization:
+        beta /= 0.5 * n_samples
+
     for it in range(max_iter):
 
         for i in range(n_samples):
@@ -92,13 +96,14 @@ def test_adagrad_decrease():
 def check_adagrad_fit(degree):
     y = _poly_predict(X, P, lams, kernel="anova", degree=degree)
 
-    est = FactorizationMachineRegressor(degree=degree, n_components=3,
+    est = FactorizationMachineRegressor(degree=degree, n_components=5,
                                         fit_linear=True, fit_lower=None,
                                         solver='adagrad',
                                         init_lambdas='ones',
-                                        max_iter=30000,
-                                        learning_rate=0.1,
-                                        beta=1e-8,
+                                        max_iter=2000,
+                                        learning_rate=0.25,
+                                        alpha=1e-10,
+                                        beta=1e-10,
                                         random_state=0)
 
     est.fit(X, y)
@@ -116,7 +121,7 @@ def test_adagrad_fit():
 
 def check_adagrad_same_as_slow(degree, sparse):
 
-    beta = 0.00001
+    beta = 1e-5
     lr = 0.01
 
     if sparse:
@@ -128,7 +133,7 @@ def check_adagrad_same_as_slow(degree, sparse):
 
     y = _poly_predict(X, P, lams, kernel="anova", degree=degree)
 
-    P_fast = 0.01 * np.random.RandomState(42).randn(1, P.shape[0], P.shape[1])
+    P_fast = np.random.RandomState(42).randn(1, P.shape[0], P.shape[1])
     P_slow = P_fast[0].copy()
 
     reg = FactorizationMachineRegressor(degree=degree, n_components=P.shape[0],

polylearn/tests/test_factorization_machine.py

@@ -26,7 +26,7 @@ def cd_direct_slow(X, y, lams=None, degree=2, n_components=5, beta=1.,
     n_samples, n_features = X.shape
 
     rng = check_random_state(random_state)
-    P = 0.01 * rng.randn(n_components, n_features)
+    P = rng.randn(n_components, n_features)
     if lams is None:
         lams = np.ones(n_components)
 
@@ -222,6 +222,7 @@ def check_same_as_slow(degree):
     y = _poly_predict(X, P, lams, kernel="anova", degree=degree)
 
     reg = FactorizationMachineRegressor(degree=degree, n_components=5,
+                                        scale_regularization=False,
                                         fit_lower=None, fit_linear=False,
                                         beta=1, warm_start=False, tol=1e-3,
                                         max_iter=5, random_state=0)
@@ -303,21 +304,20 @@ def check_warm_start(degree):
     X_train, X_test = X[:10], X[10:]
     y_train, y_test = noisy_y[:10], noisy_y[10:]
 
-    beta_low = 0.5
-    beta = 0.1
-    beta_hi = 1
+    beta_low = 0.001
+    beta = 0.002
+    beta_hi = 0.003
     ref = FactorizationMachineRegressor(degree=degree, n_components=5,
                                         fit_linear=False, fit_lower=None,
-                                        beta=beta, max_iter=20000,
-                                        random_state=0)
+                                        beta=beta, random_state=0)
     ref.fit(X_train, y_train)
     y_pred_ref = ref.predict(X_test)
 
     # (a) starting from lower beta, increasing and refitting
     from_low = FactorizationMachineRegressor(degree=degree, n_components=5,
                                              fit_lower=None, fit_linear=False,
                                              beta=beta_low, warm_start=True,
-                                             random_state=0)
+                                             max_iter=20000, random_state=0)
     from_low.fit(X_train, y_train)
     from_low.set_params(beta=beta)
     from_low.fit(X_train, y_train)
@@ -339,4 +339,4 @@ def check_warm_start(degree):
 
 def test_warm_start():
     yield check_warm_start, 2
-    yield check_warm_start, 3
+    yield check_warm_start, 3