ENH: line search for step size

Author: jasmainak

pyglmnet/pyglmnet.py

@@ -3,10 +3,13 @@
 from copy import deepcopy
 
 import numpy as np
+from scipy import optimize
 from scipy.special import expit
 from scipy.stats import norm
 
-from .utils import logger, set_log_level, power_iteration
+from functools import partial
+
+from .utils import logger, set_log_level
 from .base import BaseEstimator, is_classifier, check_version
 
 
@@ -192,17 +195,6 @@ def _grad_L2loss(distr, alpha, Tau, reg_lambda, X, y, eta, beta):
     return g
 
 
-def _learning_rate(distr, X, reg_lambda, alpha):
-    if distr == 'gaussian':
-        s = power_iteration(X.T.dot(X)) + reg_lambda * (1 - alpha)
-        return 0.99 / s
-    elif distr == 'binomial':
-        s = (np.linalg.norm(X.T.dot(X)) ** 2) / 4
-        return 0.99 / s
-    else:
-        return 1e-4
-
-
 def _gradhess_logloss_1d(distr, xk, y, z, eta):
     """
     Compute gradient (1st derivative)
@@ -380,8 +372,8 @@ class GLM(BaseEstimator):
         'cdfast' (Newton coordinate gradient descent).
         default: 'batch-gradient'
     learning_rate : float | 'auto'
-        learning rate for gradient descent. If "auto", it is 0.95 / L
-        where the differentiable part of the loss function is L-smooth.
+        learning rate for gradient descent. If "auto", backtracking line
+        search is performed using scipy.optimize.line_search.
         default: "auto"
     max_iter : int
         maximum iterations for the model.
@@ -627,12 +619,6 @@ def fit(self, X, y):
         self : instance of GLM
             The fitted model.
         """
-        if self.learning_rate == 'auto':
-            step_size = _learning_rate(self.distr, X,
-                                       self.reg_lambda, self.alpha)
-            print('Step size calculated as %f' % step_size)
-        else:
-            step_size = self.learning_rate
         np.random.RandomState(self.random_state)
 
         # checks for group
@@ -682,6 +668,17 @@ def fit(self, X, y):
                                     reg_lambda, X, y, self.eta,
                                     beta)
 
+                if self.learning_rate == 'auto':
+                    func = partial(_loss, self.distr, alpha, self.Tau,
+                                   reg_lambda, X, y, self.eta, self.group)
+                    fprime = partial(_grad_L2loss, self.distr, alpha, self.Tau,
+                                     reg_lambda, X, y, self.eta)
+                    step_size, _, _, _, _, _ = optimize.linesearch.line_search(
+                        func, fprime, beta, -grad, grad, c1=1e-4)
+                    if step_size is None:
+                        step_size = 1e-4
+                else:
+                    step_size = self.learning_rate
                 beta = beta - step_size * grad
             elif self.solver == 'cdfast':
                 beta, z = \
@@ -698,16 +695,15 @@ def fit(self, X, y):
             # Compute and save loss
             L.append(_loss(self.distr, alpha, self.Tau, reg_lambda,
                            X, y, self.eta, self.group, beta))
-            print(L[-1])
-            # if t > 1:
-            #     DL.append(L[-1] - L[-2])
-            #     if np.abs(DL[-1] / L[-1]) < tol:
-            #         msg = ('\tConverged. Loss function:'
-            #                ' {0:.2f}').format(L[-1])
-            #         logger.info(msg)
-            #         msg = ('\tdL/L: {0:.6f}\n'.format(DL[-1] / L[-1]))
-            #         logger.info(msg)
-            #         break
+            if t > 1:
+                DL.append(L[-1] - L[-2])
+                if np.abs(DL[-1] / L[-1]) < tol:
+                    msg = ('\tConverged. Loss function:'
+                           ' {0:.2f}').format(L[-1])
+                    logger.info(msg)
+                    msg = ('\tdL/L: {0:.6f}\n'.format(DL[-1] / L[-1]))
+                    logger.info(msg)
+                    break
 
         # Update the estimated variables
         self.beta0_ = beta[0]
@@ -906,8 +902,8 @@ class GLMCV(object):
         'cdfast' (Newton coordinate gradient descent).
         default: 'batch-gradient'
     learning_rate : float | 'auto'
-        learning rate for gradient descent. If "auto", it is 0.95 / L
-        where the differentiable part of the loss function is L-smooth.
+        learning rate for gradient descent. If "auto", backtracking line
+        search is performed using scipy.optimize.line_search.
         default: "auto"
     max_iter : int
         maximum iterations for the model.

pyglmnet/utils.py

@@ -4,7 +4,6 @@
 
 import numpy as np
 from copy import copy
-from scipy import linalg
 import logging
 
 
@@ -91,38 +90,6 @@ def tikhonov_from_prior(prior_cov, n_samples, threshold=0.0001):
     return Tau
 
 
-def power_iteration(A, max_iter=1000, tol=1e-7, random_state=None):
-    """Estimate dominant eigenvalue of matrix A.
-    Parameters
-    ----------
-    A : array, shape (n_points, n_points)
-        The matrix whose largest eigenvalue is to be found.
-    b_hat_0 : array, shape (n_points, )
-        init vector
-    Returns
-    -------
-    mu_hat : float
-        The largest eigenvalue
-    """
-    rng = np.random.RandomState(random_state)
-    b_hat = rng.rand((A.shape[1]))
-
-    Ab_hat = A.dot(b_hat)
-    mu_hat = np.nan
-    for ii in range(max_iter):
-        b_hat = A.dot(b_hat)
-        b_hat /= linalg.norm(b_hat)
-        Ab_hat = A.dot(b_hat)
-        mu_old = mu_hat
-        mu_hat = np.dot(b_hat, Ab_hat)
-        # note, we might exit the loop before b_hat converges
-        # since we care only about mu_hat converging
-        if (mu_hat - mu_old) / mu_old < tol:
-            break
-
-    return mu_hat
-
-
 def set_log_level(verbose):
     """Convenience function for setting the log level.
 

tests/test_pyglmnet.py

@@ -15,7 +15,6 @@
 
 from pyglmnet import (GLM, GLMCV, _grad_L2loss, _L2loss, simulate_glm,
                       _gradhess_logloss_1d, _loss, datasets)
-from pyglmnet.utils import power_iteration
 
 
 def test_gradients():
@@ -313,11 +312,3 @@ def test_cdfast():
 def test_fetch_datasets():
     """Test fetching datasets."""
     datasets.fetch_community_crime_data('/tmp/glm-tools')
-
-
-def test_power_iterations():
-    """Test power iteration."""
-    A = np.diag((1, 2, 3))
-    mu, b = np.linalg.eig(A)
-    mu_hat = power_iteration(A)
-    assert_allclose(mu_hat, mu.max())