ENH Improve Ridge's conjugate gradient descent

Reorder terms in the algorithm so that X'X is not explicitly computed and
add LinearOperator as accepted types. This improves precedent code in two ways:

    * Speed: I get speedups ranging from 2.5x to more than 10x.
      Benchmark suite can be seen here: https://gist.github.com/1322711
      where I print the time spent for various sizes and densities.

    * Interface. By allowing a LinearOperator type, we can use the
      conjugate gradient method almost for free on dense matrices with
      special structure, pytables, etc., about anything that
      implements a fast matrix-vector product.

Author: Fabian Pedregosa

sklearn/linear_model/ridge.py

@@ -2,14 +2,19 @@
 Ridge regression
 """
 
-# Author:   Mathieu Blondel <[email protected]>
-#           Reuben Fletcher-Costin <[email protected]>
+# Author: Mathieu Blondel <[email protected]>
+#         Reuben Fletcher-Costin <[email protected]>
+#         Fabian Pedregosa <[email protected]>
 # License: Simplified BSD
 
 
 from abc import ABCMeta, abstractmethod
 import warnings
+
 import numpy as np
+from scipy import linalg
+from scipy import sparse
+from scipy.sparse import linalg as sp_linalg
 
 from .base import LinearClassifierMixin, LinearModel
 from ..base import RegressorMixin
@@ -19,44 +24,12 @@
 from ..grid_search import GridSearchCV
 
 
-def _solve(A, b, solver, tol):
-    # helper method for ridge_regression, A is symmetric positive
-
-    if solver == 'auto':
-        if hasattr(A, 'todense'):
-            solver = 'sparse_cg'
-        else:
-            solver = 'dense_cholesky'
-
-    if solver == 'sparse_cg':
-        if b.ndim < 2:
-            from scipy.sparse import linalg as sp_linalg
-            sol, error = sp_linalg.cg(A, b, tol=tol)
-            if error:
-                raise ValueError("Failed with error code %d" % error)
-            return sol
-        else:
-            # sparse_cg cannot handle a 2-d b.
-            sol = []
-            for j in range(b.shape[1]):
-                sol.append(_solve(A, b[:, j], solver="sparse_cg", tol=tol))
-            return np.array(sol).T
-
-    elif solver == 'dense_cholesky':
-        from scipy import linalg
-        if hasattr(A, 'todense'):
-            A = A.todense()
-        return linalg.solve(A, b, sym_pos=True, overwrite_a=True)
-    else:
-        raise NotImplementedError('Solver %s not implemented' % solver)
-
-
 def ridge_regression(X, y, alpha, sample_weight=1.0, solver='auto', tol=1e-3):
     """Solve the ridge equation by the method of normal equations.
 
     Parameters
     ----------
-    X : {array-like, sparse matrix}, shape = [n_samples, n_features]
+    X : {array-like, sparse matrix, LinearOperator}, shape = [n_samples, n_features]
         Training data
 
     y : array-like, shape = [n_samples] or [n_samples, n_responses]
@@ -86,26 +59,54 @@ def ridge_regression(X, y, alpha, sample_weight=1.0, solver='auto', tol=1e-3):
     """
 
     n_samples, n_features = X.shape
-    is_sparse = False
-
-    if hasattr(X, 'todense'):  # lazy import of scipy.sparse
-        from scipy import sparse
-        is_sparse = sparse.issparse(X)
 
-    if is_sparse:
-        if n_features > n_samples or \
-           isinstance(sample_weight, np.ndarray) or \
-           sample_weight != 1.0:
+    if solver == 'auto':
+        # cholesky if it's a dense array and cg in
+        # any other case
+        if hasattr(X, '__array__'):
+            solver = 'dense_cholesky'
+        else:
+            solver = 'sparse_cg'
 
-            I = sparse.lil_matrix((n_samples, n_samples))
-            I.setdiag(np.ones(n_samples) * alpha * sample_weight)
-            c = _solve(X * X.T + I, y, solver, tol)
-            coef = X.T * c
+    if solver == 'sparse_cg':
+        # gradient descent
+        X1 = sp_linalg.aslinearoperator(X)
+        if y.ndim == 1:
+            y1 = np.reshape(y, (-1, 1))
         else:
-            I = sparse.lil_matrix((n_features, n_features))
-            I.setdiag(np.ones(n_features) * alpha)
-            coef = _solve(X.T * X + I, X.T * y, solver, tol)
+            y1 = y
+        coefs = np.empty((y1.shape[1], n_features))
+
+        for i in range(y1.shape[1]):
+            y_column = y1[:, i]
+            if n_features > n_samples:
+                # kernel ridge
+                # w = X.T * inv(X X^t + alpha*Id) y
+                def mv(x):
+                    return X1.matvec(X1.rmatvec(x)) + alpha * x
+                C = sp_linalg.LinearOperator(
+                    (n_samples, n_samples), matvec=mv, dtype=X.dtype)
+                coef, info = sp_linalg.cg(C, y_column, tol=tol)
+                coefs[i] = X1.rmatvec(coef)
+            else:
+                # ridge
+                # w = inv(X^t X + alpha*Id) * X.T y
+                def mv(x):
+                    return X1.rmatvec(X1.matvec(x)) + alpha * x
+                y_column = X1.rmatvec(y_column)
+                C = sp_linalg.LinearOperator(
+                    (n_features, n_features), matvec=mv, dtype=X.dtype)
+                coefs[i], info = sp_linalg.cg(C, y_column, tol=tol)
+            if info != 0:
+                raise ValueError("Failed with error code %d" % info)
+
+        if y.ndim == 1:
+            return np.ravel(coefs)
+        return coefs
     else:
+        # normal equations (cholesky) method
+        if hasattr(X, 'todense'):
+            X = X.todense()
         if n_features > n_samples or \
            isinstance(sample_weight, np.ndarray) or \
            sample_weight != 1.0:
@@ -114,13 +115,13 @@ def ridge_regression(X, y, alpha, sample_weight=1.0, solver='auto', tol=1e-3):
             # w = X.T * inv(X X^t + alpha*Id) y
             A = np.dot(X, X.T)
             A.flat[::n_samples + 1] += alpha * sample_weight
-            coef = np.dot(X.T, _solve(A, y, solver, tol))
+            coef = np.dot(X.T, linalg.solve(A, y, sym_pos=True, overwrite_a=True))
         else:
             # ridge
             # w = inv(X^t X + alpha*Id) * X.T y
             A = np.dot(X.T, X)
             A.flat[::n_features + 1] += alpha
-            coef = _solve(A, np.dot(X.T, y), solver, tol)
+            coef = linalg.solve(A, np.dot(X.T, y), sym_pos=True, overwrite_a=True)
 
     return coef.T
 
@@ -381,7 +382,6 @@ def __init__(self, alphas=[0.1, 1.0, 10.0], fit_intercept=True,
     def _pre_compute(self, X, y):
         # even if X is very sparse, K is usually very dense
         K = safe_sparse_dot(X, X.T, dense_output=True)
-        from scipy import linalg
         v, Q = linalg.eigh(K)
         QT_y = np.dot(Q.T, y)
         return v, Q, QT_y
@@ -418,7 +418,6 @@ def _values(self, alpha, y, v, Q, QT_y):
         return y - (c / G_diag), c
 
     def _pre_compute_svd(self, X, y):
-        from scipy import sparse
         if sparse.issparse(X) and hasattr(X, 'toarray'):
             X = X.toarray()
         U, s, _ = np.linalg.svd(X, full_matrices=0)