[MRG] Uniformize initialization for all algorithms

metric_learn/_util.py

@@ -1,10 +1,16 @@
-import warnings
 import numpy as np
 import six
 from numpy.linalg import LinAlgError
+from sklearn.datasets import make_spd_matrix
+from sklearn.decomposition import PCA
 from sklearn.utils import check_array
-from sklearn.utils.validation import check_X_y
-from metric_learn.exceptions import PreprocessorError
+from sklearn.utils.validation import check_X_y, check_random_state
+from .exceptions import PreprocessorError
+from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
+from sklearn.utils.multiclass import type_of_target
+from scipy.linalg import pinvh
+import sys
+import time
 
 # hack around lack of axis kwarg in older numpy versions
 try:
@@ -405,3 +411,262 @@ def validate_vector(u, dtype=None):
   if u.ndim > 1:
     raise ValueError("Input vector should be 1-D.")
   return u
+
+
+def _initialize_transformer(num_dims, input, y=None, init='auto',
+                            verbose=False, random_state=None,
+                            has_classes=True):
+  """Returns the initial transformer to be used depending on the arguments.
+
+  Parameters
+  ----------
+  num_dims : int
+    The number of components to take. (Note: it should have been checked
+    before, meaning it should not be None and it should be a value in
+    [1, X.shape[1]])
+
+  input : array-like
+    The input samples (can be tuples or regular samples).
+
+  y : array-like or None
+    The input labels (or not if there are no labels).
+
+  init : string or numpy array, optional (default='auto')
+      Initialization of the linear transformation. Possible options are
+      'auto', 'pca', 'lda', 'identity', 'random', and a numpy array of shape
+      (n_features_a, n_features_b).
+
+      'auto'
+          Depending on ``num_dims``, the most reasonable initialization will
+          be chosen. If ``num_dims <= n_classes`` we use 'lda' (if possible,
+          see the description of 'lda' init), as it uses labels information.
+          If not, but ``num_dims < min(n_features, n_samples)``, we use
+          'pca', as it projects data in meaningful directions (those of
+          higher variance). Otherwise, we just use 'identity'.
+
+      'pca'
+          ``num_dims`` principal components of the inputs passed
+          to :meth:`fit` will be used to initialize the transformation.
+          (See `sklearn.decomposition.PCA`)
+
+      'lda'
+          ``min(num_dims, n_classes)`` most discriminative
+          components of the inputs passed to :meth:`fit` will be used to
+          initialize the transformation. (If ``num_dims > n_classes``,
+          the rest of the components will be zero.) (See
+          `sklearn.discriminant_analysis.LinearDiscriminantAnalysis`).
+          This initialization is possible only if `has_classes == True`.
+
+      'identity'
+          If ``num_dims`` is strictly smaller than the
+          dimensionality of the inputs passed to :meth:`fit`, the identity
+          matrix will be truncated to the first ``num_dims`` rows.
+
+      'random'
+          The initial transformation will be a random array of shape
+          `(num_dims, n_features)`. Each value is sampled from the
+          standard normal distribution.
+
+      numpy array
+          n_features_b must match the dimensionality of the inputs passed to
+          :meth:`fit` and n_features_a must be less than or equal to that.
+          If ``num_dims`` is not None, n_features_a must match it.
+
+  verbose : bool
+    Whether to print the details of the initialization or not.
+
+  random_state : int or `numpy.RandomState` or None, optional (default=None)
+    A pseudo random number generator object or a seed for it if int. If
+    ``init='random'``, ``random_state`` is used to initialize the random
+    transformation. If ``init='pca'``, ``random_state`` is passed as an
+    argument to PCA when initializing the transformation.
+
+  has_classes : bool (default=True)
+    Whether the labels are in fact classes. If true, this will allow to use
+    the 'lda' initialization.
+
+  Returns
+  -------
+  init_transformer : `numpy.ndarray`
+    The initial transformer to use.
+  """
+  # if we are doing a regression we cannot use lda:
+  n_features = input.shape[-1]
+  authorized_inits = ['auto', 'pca', 'identity', 'random']
+  if has_classes:
+    authorized_inits.append('lda')
+
+  if isinstance(init, np.ndarray):
+    init = check_array(init)
+
+    # Assert that init.shape[1] = X.shape[1]
+    if init.shape[1] != n_features:
+      raise ValueError('The input dimensionality ({}) of the given '
+                       'linear transformation `init` must match the '
+                       'dimensionality of the given inputs `X` ({}).'
+                       .format(init.shape[1], n_features))
+
+    # Assert that init.shape[0] <= init.shape[1]
+    if init.shape[0] > init.shape[1]:
+      raise ValueError('The output dimensionality ({}) of the given '
+                       'linear transformation `init` cannot be '
+                       'greater than its input dimensionality ({}).'
+                       .format(init.shape[0], init.shape[1]))
+
+    # Assert that self.num_dims = init.shape[0]
+    if num_dims != init.shape[0]:
+      raise ValueError('The preferred dimensionality of the '
+                       'projected space `num_dims` ({}) does'
+                       ' not match the output dimensionality of '
+                       'the given linear transformation '
+                       '`init` ({})!'
+                       .format(num_dims,
+                               init.shape[0]))
+  elif init in authorized_inits:
+    pass
+  else:
+    raise ValueError(
+        "`init` must be '{}' "
+        "or a numpy array of shape (num_dims, n_features)."
+        .format("', '".join(authorized_inits)))
+
+  random_state = check_random_state(random_state)
+  transformation = init
+  if isinstance(init, np.ndarray):
+    pass
+  else:
+    n_samples = input.shape[0]
+    if init == 'auto':
+      if has_classes:
+        n_classes = len(np.unique(y))
+      if (has_classes and num_dims <= min(n_features, n_classes - 1)):
+        init = 'lda'
+      elif num_dims < min(n_features, n_samples):
+        init = 'pca'
+      else:
+        init = 'identity'
+    if init == 'identity':
+      transformation = np.eye(num_dims, input.shape[-1])
+    elif init == 'random':
+      transformation = random_state.randn(num_dims,
+                                          input.shape[-1])
+    elif init in {'pca', 'lda'}:
+      init_time = time.time()
+      if init == 'pca':
+        pca = PCA(n_components=num_dims,
+                  random_state=random_state)
+        if verbose:
+          print('Finding principal components... ')
+          sys.stdout.flush()
+        pca.fit(input)
+        transformation = pca.components_
+      elif init == 'lda':
+        lda = LinearDiscriminantAnalysis(n_components=num_dims)
+        if verbose:
+          print('Finding most discriminative components... ')
+          sys.stdout.flush()
+        lda.fit(input, y)
+        transformation = lda.scalings_.T[:num_dims]
+      if verbose:
+        print('done in {:5.2f}s'.format(time.time() - init_time))
+  return transformation
+
+
+def _initialize_metric_mahalanobis(input, init='identity', random_state=None,
+                                   return_inverse=False):
+  """Returns the initial mahalanobis matrix to be used depending on the
+  arguments.
+
+  Parameters
+  ----------
+  input : array-like
+    The input samples (can be tuples or regular samples).
+
+  init : string or numpy array, optional (default='identity')
+         Initialization of the linear transformation. Possible options are
+         'identity', 'covariance', 'random', and a numpy array of shape
+         (n_features, n_features).
+
+         'identity'
+            An identity matrix of shape (n_features, n_features).
+
+         'covariance'
+            The inverse covariance matrix.
+
+         'random'
+             The initial transformation will be a random SPD matrix of shape
+             `(n_features, n_features)`, generated using
+             `sklearn.datasets.make_spd_matrix`.
+
+         numpy array
+             An SPD matrix of shape (n_features, n_features), that will
+             be used as such to initialize the metric.
+
+  random_state : int or `numpy.RandomState` or None, optional (default=None)
+    A pseudo random number generator object or a seed for it if int. If
+    ``init='random'``, ``random_state`` is used to initialize the random
+    matrix. If ``init='pca'``, ``random_state`` is passed as an
+    argument to PCA when initializing the matrix.
+
+  return_inverse : bool, optional (default=False)
+    Whether to return the inverse of the matrix initializing the metric. This
+    can be sometimes useful.
+
+  Returns
+  -------
+  M, or (M, M_inv) : `numpy.ndarray`
+    The initial matrix to use M, and its inverse if `return_inverse=True`.
+  """
+  n_features = input.shape[-1]
+  if isinstance(init, np.ndarray):
+    init = check_array(init)  # TODO: do we want to copy the array ?
+    # see how they do it in scikit-learn for instance
+
+    # Assert that init.shape[1] = n_features
+    if (init.shape) != (n_features,) * 2:
+      raise ValueError('The input dimensionality {} of the given '
+                       'mahalanobis matrix `init` must match the '
+                       'dimensionality of the given inputs ({}).'
+                       .format(init.shape, n_features))
+
+    # Assert that the matrix is symmetric
+    if not np.allclose(init, init.T):
+      raise ValueError("The initialization matrix should be semi-definite "
+                       "positive (SPD). It is not, since it appears not to be "
+                       "symmetric.")
+
+  elif init in ['identity', 'covariance', 'random']:
+    pass
+  else:
+    raise ValueError(
+        "`init` must be 'identity', 'covariance', 'random' "
+        "or a numpy array of shape (n_features, n_features).")
+
+  random_state = check_random_state(random_state)
+  M = init
+  if isinstance(init, np.ndarray):
+    if return_inverse:
+      M_inv = pinvh(M)
+  else:
+    if init == 'identity':
+      M = np.eye(n_features, n_features)
+      if return_inverse:
+        M_inv = M.copy()
+    if init == 'covariance':
+      if input.ndim == 3:
+        # if the input are tuples, we need to form an X by deduplication
+        X = np.vstack({tuple(row) for row in input.reshape(-1, n_features)})
+      else:
+        X = input
+      M_inv = np.atleast_2d(np.cov(X, rowvar=False))
+      # TODO: check atleast_2d necessary
+      M = pinvh(M_inv)
+    elif init == 'random':
+      # we need to create a random symmetric matrix
+      M = make_spd_matrix(n_features, random_state=random_state)
+      if return_inverse:
+        M_inv = pinvh(M)
+  if return_inverse:
+    return (M, M_inv)
+  else:
+    return M

metric_learn/base_metric.py

@@ -288,7 +288,7 @@ def get_mahalanobis_matrix(self):
 
     Returns
     -------
-    M : `numpy.ndarray`, shape=(n_components, n_features)
+    M : `numpy.ndarray`, shape=(num_dims, n_features)
       The copy of the learned Mahalanobis matrix.
     """
     return self.transformer_.T.dot(self.transformer_)

metric_learn/covariance.py

@@ -21,7 +21,7 @@ class Covariance(MahalanobisMixin, TransformerMixin):
 
   Attributes
   ----------
-  transformer_ : `numpy.ndarray`, shape=(num_dims, n_features)
+  transformer_ : `numpy.ndarray`, shape=(n_features, n_features)
       The linear transformation ``L`` deduced from the learned Mahalanobis
       metric (See function `transformer_from_metric`.)
   """