Issue/64/mp_sparse_encode

spmimage/decomposition/__init__.py

@@ -1,7 +1,8 @@
 from .ksvd import KSVD
-from .dict_learning import sparse_encode_with_mask
+from .dict_learning import sparse_encode_with_mask, sparse_encode
 
 __all__ = [
     'KSVD',
     'sparse_encode_with_mask'
+    'sparse_encode'
 ]

spmimage/decomposition/dict_learning.py

@@ -1,6 +1,9 @@
 import numpy as np
-
-from sklearn.decomposition import sparse_encode
+import sklearn
+from sklearn.externals.joblib import Parallel, delayed, effective_n_jobs
+from sklearn.utils import (check_array, check_random_state, gen_even_slices,
+                     gen_batches)
+from ..linear_model import matching_pursuit
 
 
 def sparse_encode_with_mask(X, dictionary, mask, **kwargs):
@@ -42,7 +45,100 @@ def sparse_encode_with_mask(X, dictionary, mask, **kwargs):
     """
     code = np.zeros((X.shape[0], dictionary.shape[0]))
     for idx in range(X.shape[0]):
-        code[idx, :] = sparse_encode(X[idx, :][mask[idx, :] == 1].reshape(1, -1),
-                                     dictionary[:, mask[idx, :] == 1],
-                                     **kwargs)
+        code[idx, :] = sparse_encode(
+            X[idx, :][mask[idx, :] == 1].reshape(1, -1),
+            dictionary[:, mask[idx, :] == 1], **kwargs)
+    return code
+
+def sparse_encode(X, dictionary, gram=None, cov=None, algorithm='lasso_lars',
+                  n_nonzero_coefs=None, alpha=None, copy_cov=True, init=None,
+                  max_iter=1000, n_jobs=1, check_input=True, verbose=0):
+    """Sparse coding
+    Each row of the result is the solution to a sparse coding problem.
+    The goal is to find a sparse array `code` such that::
+        X ~= code * dictionary
+    Read more in the :ref:`User Guide <SparseCoder>`.
+    Parameters
+    ----------
+        X : array of shape (n_samples, n_features)
+            Data matrix
+        dictionary : array of shape (n_components, n_features)
+            The dictionary matrix against which to solve the sparse coding of
+            the data. Some of the algorithms assume normalized rows for meaningful
+            output (particularly for 'omp' and 'mp').
+        gram : array, shape=(n_components, n_components)
+            Precomputed Gram matrix, dictionary * dictionary'
+        cov : array, shape=(n_components, n_samples)
+            Precomputed covariance, dictionary' * X
+        algorithm : {'lasso_lars', 'lasso_cd', 'lars', 'omp', 'threshold', 'mp'}
+            lars: uses the least angle regression method (linear_model.lars_path)
+            lasso_lars: uses Lars to compute the Lasso solution
+            lasso_cd: uses the coordinate descent method to compute the
+            Lasso solution (linear_model.Lasso). lasso_lars will be faster if
+            the estimated components are sparse.
+            omp: uses orthogonal matching pursuit to estimate the sparse solution
+            threshold: squashes to zero all coefficients less than alpha from
+            the projection dictionary * X'
+            mp: uses matching pursuit to estimate the sparse solution
+        n_nonzero_coefs : int, 0.1 * n_features by default
+            Number of nonzero coefficients to target in each column of the
+            solution. This is only used by `algorithm='lars'` and `algorithm='omp'`
+            and is overridden by `alpha` in the `omp` case.
+        alpha : float, 1. by default
+            If `algorithm='lasso_lars'` or `algorithm='lasso_cd'`, `alpha` is the
+            penalty applied to the L1 norm.
+            If `algorithm='threshold'`, `alpha` is the absolute value of the
+            threshold below which coefficients will be squashed to zero.
+            If `algorithm='omp'`, `alpha` is the tolerance parameter: the value of
+            the reconstruction error targeted. In this case, it overrides
+            `n_nonzero_coefs`.
+        copy_cov : boolean, optional
+            Whether to copy the precomputed covariance matrix; if False, it may be
+            overwritten.
+        init : array of shape (n_samples, n_components)
+            Initialization value of the sparse codes. Only used if
+            `algorithm='lasso_cd'`.
+        max_iter : int, 1000 by default
+            Maximum number of iterations to perform if `algorithm='lasso_cd'`.
+        n_jobs : int or None, optional (default=None)
+            Number of parallel jobs to run.
+            ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context.
+            ``-1`` means using all processors. See :term:`Glossary <n_jobs>`
+            for more details.
+        check_input : boolean, optional
+            If False, the input arrays X and dictionary will not be checked.
+        verbose : int, optional
+            Controls the verbosity; the higher, the more messages. Defaults to 0.
+    Returns
+    -------
+        code : array of shape (n_samples, n_components)
+            The sparse codes
+    """
+    if algorithm != 'mp':
+        return sklearn.decomposition.sparse_encode(
+            X, dictionary, gram, cov, algorithm, n_nonzero_coefs, alpha, copy_cov,
+            init, max_iter, n_jobs, check_input, verbose)
+    elif check_input:
+        dictionary = check_array(dictionary)
+        X = check_array(X)
+
+    n_samples, n_features = X.shape
+    n_components = dictionary.shape[0]
+
+    if n_nonzero_coefs is None:
+        n_nonzero_coefs = min(max(n_features / 10, 1), n_components)
+
+    if n_jobs == 1:
+        code = matching_pursuit(dictionary=dictionary, signal=X, n_nonzero_coefs=n_nonzero_coefs)
+        return code
+
+    # Enter parallel code block
+    code = np.empty((n_samples, n_components))
+    slices = list(gen_even_slices(n_samples, effective_n_jobs(n_jobs)))
+
+    code_views = Parallel(n_jobs=n_jobs, verbose=verbose)(
+        delayed(matching_pursuit)(dictionary=dictionary, signal=X[this_slice], n_nonzero_coefs=n_nonzero_coefs) 
+        for this_slice in slices)
+    for this_slice, this_view in zip(slices, code_views):
+        code[this_slice] = this_view
     return code

spmimage/linear_model/__init__.py

@@ -1,8 +1,10 @@
 from .admm import LassoADMM, FusedLassoADMM, TrendFilteringADMM, QuadraticTrendFilteringADMM
+from .mp import matching_pursuit
 
 __all__ = [
     'LassoADMM',
     'FusedLassoADMM',
     'TrendFilteringADMM',
     'QuadraticTrendFilteringADMM',
+    'matching_pursuit'
 ]

spmimage/linear_model/mp.py

@@ -0,0 +1,63 @@
+"""Matching pursuit algorithms
+"""
+
+import numpy as np
+from sklearn.utils import check_array
+from sklearn.preprocessing import normalize
+
+def matching_pursuit(dictionary, signal, n_nonzero_coefs=None, copy_dictionary=True, copy_signal=True, tol=None):
+    """Matching Pursuit (MP)
+
+    Solves n_targets Matching Pursuit problems.
+    The goal is to find a sparse array `coefs` such that::
+        signal ~= coefs * dictionary
+    Parameters
+    ----------
+    dictionary : array, shape (n_components, n_features)
+        Input data. Columns are assumed to have unit norm.
+
+    signal : array, shape (n_samples, n_features)
+        Input targets.
+
+    n_nonzero_coefs : int
+        Targeted number of non-zero elements
+
+    copy_X : bool, optional
+        Whether the input data must be copied by the algorithm. A false
+        value is only helpful if it is already Fortran-ordered, otherwise a
+        copy is made anyway.
+
+    copy_y : bool, optional
+        Whether the covariance vector Xy must be copied by the algorithm.
+        If False, it may be overwritten.
+
+    tol : float
+        Maximum norm of the residual.
+    """
+    normalized = normalize(dictionary)
+    if abs(np.linalg.norm(normalized - dictionary, 'fro')) > len(dictionary) * 1e-8:
+        raise ValueError("All columns of the dictionary must have unit norm.")
+    if tol is not None and tol < 0:
+        raise ValueError("Epsilon cannot be negative")
+    dictionary = check_array(dictionary, order='F', copy=copy_dictionary)
+    if copy_signal:
+        signal = signal.copy()
+
+    n_active = 0
+
+    n_samples = signal.shape[0]
+    n_components = dictionary.shape[0]
+    n_features = dictionary.shape[1]
+    coefs = np.zeros((n_samples, n_components))
+    while(True):
+        inners = dictionary.dot(signal.T)
+        max_ids = np.argmax(np.abs(inners), axis=0)
+        for i, max_id in enumerate(max_ids):
+            coefs[i, max_id] += inners[max_id, i]
+            signal[i] -= inners[max_id, i] * dictionary[max_id]
+        n_active += 1
+        if tol is not None and np.linalg.norm(signal) <= tol:
+                break
+        if n_active == n_features or n_active == n_nonzero_coefs:
+            break
+    return coefs

tests/test_decomposition_dict_learning.py

@@ -1,7 +1,9 @@
 import unittest
 
+from sklearn.preprocessing import normalize
 from sklearn.decomposition import sparse_encode
 from spmimage.decomposition import sparse_encode_with_mask
+import spmimage
 
 import numpy as np
 
@@ -41,7 +43,49 @@ def test_sparse_encode_with_mask(self):
 
         # check error of learning
         # print(np.linalg.norm(mask*(X-W.dot(A0)), 'fro'))
-        self.assertTrue(np.linalg.norm(mask * (X - W.dot(A0)), 'fro') < 50)
+        self.assertTrue(np.linalg.norm(mask * (X - W.dot(A0)), 'fro') < np.linalg.norm(mask * X, 'fro'))
+
+    def test_sparse_encode(self):
+        k0 = 3
+        n_samples = 64
+        n_features = 32
+        n_components = 10
+
+        A0, X = generate_dictionary_and_samples(n_samples, n_features, n_components, k0)
+
+        W1 = sparse_encode(X, A0, algorithm='omp', n_nonzero_coefs=k0)
+        W2 = spmimage.decomposition.sparse_encode(X, A0, algorithm='omp', n_nonzero_coefs=k0)
+
+        # check if W1 and W2 is almost same
+        self.assertTrue(abs(np.linalg.norm(X - W1.dot(A0), 'fro') - np.linalg.norm(X - W2.dot(A0), 'fro')) < 1e-8)
+
+    def test_sparse_encode_mp(self):
+        k0 = 3
+        n_samples = 64
+        n_features = 32
+        n_components = 10
+
+        A0, X = generate_dictionary_and_samples(n_samples, n_features, n_components, k0)
+
+        W = spmimage.decomposition.sparse_encode(X, A0, algorithm='mp', n_nonzero_coefs=k0)
+
+        # check error of learning
+        # print(np.linalg.norm(X - W.dot(A0), 'fro'))
+        self.assertTrue(np.linalg.norm(X - W.dot(A0), 'fro') < np.linalg.norm(X))
+
+    def test_parallel_sparse_encode_mp(self):
+        k0 = 3
+        n_samples = 64
+        n_features = 32
+        n_components = 10
+
+        A0, X = generate_dictionary_and_samples(n_samples, n_features, n_components, k0)
+
+        W = spmimage.decomposition.sparse_encode(X, A0, algorithm='mp', n_nonzero_coefs=k0, n_jobs=2)
+
+        # check error of learning
+        # print(np.linalg.norm(X - W.dot(A0), 'fro'))
+        self.assertTrue(np.linalg.norm(X - W.dot(A0), 'fro') < np.linalg.norm(X))
 
 
 if __name__ == '__main__':

tests/test_linear_model_mp.py

@@ -0,0 +1,22 @@
+import unittest
+
+from spmimage.linear_model import matching_pursuit
+import numpy as np
+from sklearn.preprocessing import normalize
+from numpy.testing import assert_array_almost_equal
+
+
+class TestMatchingPursuit(unittest.TestCase):
+    def setUp(self):
+        np.random.seed()
+
+    def test_matching_pursuit(self):
+        X = np.array([[1., 2.], [3., 4.], [5., 6.]])
+        beta = np.array([[0., 1., 0.], [3., 0., 0.]])
+        y = beta.dot(X)
+        with self.assertRaises(ValueError):
+            predict = matching_pursuit(dictionary=X, signal=y, tol=1e-6)
+        X = normalize(X)
+        y = beta.dot(X)
+        predict = matching_pursuit(dictionary=X, signal=y, tol=1e-6)
+        assert_array_almost_equal(beta, predict)

tests/utils.py

@@ -1,6 +1,7 @@
 from typing import Tuple
 
 import numpy as np
+from sklearn.preprocessing import normalize
 
 
 def generate_dictionary_and_samples(n_samples: int, n_features: int, n_components: int, n_nonzero_coefs: int) \
@@ -14,4 +15,5 @@ def generate_dictionary_and_samples(n_samples: int, n_features: int, n_component
         # select n_nonzero_coefs components from dictionary
         X[i, :] = np.dot(np.random.randn(n_nonzero_coefs),
                          A0[np.random.permutation(range(n_components))[:n_nonzero_coefs], :])
-    return A0, X
+    A0 = normalize(A0)
+    return A0, X # A0 is normalized