Vectorized solver

examples/plot_logistic_regression_vectorized.py

@@ -0,0 +1,137 @@
+# -*- coding: utf-8 -*-
+"""
+Vectorizer solver for linear problem
+====================================
+
+Parsimony authorizes parallel problem solving leading to important (~5) acceleration
+factor.
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+import time
+import parsimony.datasets as datasets
+import parsimony.estimators as estimators
+import parsimony.functions.nesterov.tv as nesterov_tv
+from parsimony.utils.penalties import l1_max_logistic_loss
+import parsimony.utils as utils
+
+###############################################################################
+# Fetch dice5 dataset
+dataset_name = "%s_%s_%ix%ix%i_%i_dataset_v%s.npz" % \
+         tuple(["dice5", "classif", 50, 50, 1, 500, '0.3.1'])
+_, data  = datasets.utils.download_dataset(dataset_name)
+
+X, y, beta3d, = data['X'], data['y'], data['beta3d']
+
+###############################################################################
+# Solve in parallel many Enet-TV problems
+# ---------------------------------------
+#
+# Empirically set the global penalty, based on maximum l1 penaly
+
+alpha = l1_max_logistic_loss(X, y)
+
+###############################################################################
+# Penalization parameters are now vectors of equal length
+
+l1 = alpha * np.array([0.5, 0.5, 0.5])
+l2 = alpha * np.array([0.5, 0.5, 0.5])
+tv = alpha * np.array([0.01, 0.2, 0.8])
+max_iter = 1000
+
+###############################################################################
+# Build linear operator and fit the model:
+A = nesterov_tv.linear_operator_from_shape(beta3d.shape, calc_lambda_max=True)
+enettv = estimators.LogisticRegressionL1L2TV(
+               l1=l1,
+               l2=l2,
+               tv=tv,
+               A = A,
+               algorithm_params=dict(max_iter=max_iter))
+
+enettv.fit(X, y)
+
+###############################################################################
+# Plot coeffitients maps
+
+plt.clf()
+plot = plt.subplot(221)
+utils.plots.map2d(beta3d, plot, title="beta star")
+
+for i in range(len(l1)):
+    print(i)
+    plot = plt.subplot(222+i)
+    utils.plots.map2d(enettv.beta[:, i].reshape(beta3d.shape), plot, #limits=[-0.01, 0.01],
+               title="enettv (l1:%.3f, l2:%.3f, tv:%.3f)" % (l1[i], l2[i], tv[i]))
+
+plt.tight_layout()
+
+###############################################################################
+# Evaluate the accélération factor
+# --------------------------------
+#
+# ranges of alphas and l1 ratios
+
+tvs = alpha * np.array([0.2, 0.4, 0.6, 0.8])
+l1s = np.array([0.1, 0.5])
+#l1s = np.array([0.1, 0.5, .9])
+l2s = 1 - l1
+l1s *= alpha
+l2s *= alpha
+
+# Cartesian product (all possible combinations)
+import itertools
+params = np.array([param for param in itertools.product(l1s, l2s, tvs)])
+#print(params)
+print(params.shape)
+
+step_size = 2
+sizes = np.arange(1, min(100, params.shape[0]+1), step_size)
+max_iter = 1000
+elapsed = list()
+
+for s in sizes:
+    enettv = estimators.LogisticRegressionL1L2TV(
+               l1=params[:s, 0],
+               l2=params[:s, 1],
+               tv=params[:s, 2],
+               A = A,
+               algorithm_params=dict(max_iter=max_iter))
+
+
+    t_ = time.clock()
+    yte_pred_enettv = enettv.fit(X, y)
+    delta_time = time.clock() - t_
+    print("Vectorized fit of %i problems took %.2f sec" % (s, delta_time))
+    elapsed.append(delta_time)
+
+elapsed = np.array(elapsed)
+
+plt.clf()
+plt.subplot(311)
+plt.plot(sizes, elapsed,  'b-', label="Achieved")
+plt.plot([sizes[0], sizes[-1]], [elapsed[0],  sizes[-1]*elapsed[0]], 'g-', label="Expected")
+plt.legend()
+plt.xlabel("Nb of problems (k) solved in parallel")
+plt.ylabel("CPU time (S)")
+plt.grid()
+
+plt.subplot(312)
+plt.plot(sizes, elapsed / elapsed[0])
+plt.xlabel("Nb of problems (k) solved in parallel")
+plt.ylabel("CPU time ratio: size k / size 1")
+plt.grid()
+
+plt.subplot(313)
+plt.plot(sizes, (elapsed[0] * sizes) / elapsed)
+plt.xlabel("Nb of problems (k) solved in parallel")
+plt.ylabel("Acceleration factor")
+plt.grid()
+
+plt.tight_layout()
+
+accel = elapsed[0] / (elapsed[-1] / sizes[-1])
+
+print("Solving %i pbs in parallel is %.2f time longer. Acceleration factor is %.2f" % \
+      (sizes[-1], elapsed[-1] / elapsed[0], accel))
+

parsimony/algorithms/proximal.py

@@ -365,25 +365,26 @@ def run(self, function, beta):
                                    max_iter=self.max_iter)
 
                 # TODO: Warn if G_new < -consts.TOLERANCE.
-                gap = abs(gap)  # May happen close to machine epsilon.
+                gap = np.abs(gap)  # May happen close to machine epsilon.
                 if self.info_requested(Info.gap):
                     gap_.append(gap)
 
                 if not self.simulation:
                     if self.info_requested(Info.verbose):
                         print("FISTA ite:%i, gap:%g" % (i, gap))
-                    if gap < self.eps:
+                    if np.all(gap < self.eps):
                         if self.info_requested(Info.converged):
                             self.info_set(Info.converged, True)
 
                         break
             else:
                 if not self.simulation:
-                    eps_cur = maths.norm(betanew - z)
+                    eps_cur = maths.norm(betanew - z, axis=0)
                     if self.info_requested(Info.verbose):
                         print("FISTA ite: %i, eps_cur:%g" % (i, eps_cur))
-                    if step > 0.0:
-                        if (1.0 / step) * eps_cur < self.eps \
+                    if np.all(step > 0.0):
+                        converged = (1.0 / step) * eps_cur < self.eps
+                        if np.all(converged) \
                                 and i >= self.min_iter:
 
                             if self.info_requested(Info.converged):
@@ -392,7 +393,8 @@ def run(self, function, beta):
                             break
 
                     else:  # TODO: Fix this!
-                        if maths.norm(betanew - z) < self.eps \
+                        converged = maths.norm(betanew - z, axis=0) < self.eps
+                        if np.all(converged) \
                                 and i >= self.min_iter:
 
                             if self.info_requested(Info.converged):
@@ -522,30 +524,32 @@ def run(self, function, beta):
         # Compute current gap, precision eps (gap decreased by tau) and mu.
         function.set_mu(consts.TOLERANCE)
         gap = function.gap(beta, eps=self.eps, max_iter=self.max_iter)
-        eps = self.tau * abs(gap)
+        eps = self.tau * np.abs(gap)
         # Warning below if gap < -consts.TOLERANCE: See Special case 1
         gM = function.eps_max(1.0)
         loop = True
 
         # Special case 1: gap is very small: stopping criterion satisfied
-        if gap < self.eps:  # "- mu * gM" has been removed since mu == 0
+        if np.any(gap < self.eps):  # "- mu * gM" has been removed since mu == 0
             warnings.warn(
                 "Stopping criterion satisfied before the first iteration."
                 " Either beta_start a the solution (given eps)."
                 " If beta_start is null the problem might be over-penalized. "
                 " Then try smaller penalization.")
+
+        if np.all(gap < self.eps):
             loop = False
 
         # Special case 2: gap infinite or NaN => eps is not finite or NaN
         # => mu is NaN etc. Force eps to a large value, to force some FISTA
-        # iteration to getbetter starting point
-        if not np.isfinite(eps):
-            eps = self.eps_max
+        # iteration to get better starting point
+        eps[~np.isfinite(eps)] = self.eps_max
+        # if not np.isfinite(eps):
+        #    eps = self.eps_max
 
         if loop:  # mu is useless if loop is False
             mu = function.mu_opt(eps)
             function.set_mu(mu)
-        #gM = function.eps_max(1.0)
 
         # Initialise info variables. Info variables have the suffix "_".
         if self.info_requested(Info.time):
@@ -563,10 +567,9 @@ def run(self, function, beta):
 
         i = 0  # Iteration counter.
         while loop:
-            converged = False
 
             # Current precision.
-            eps_mu = max(eps, self.eps) - mu * gM
+            eps_mu = np.maximum(eps, self.eps) - mu * gM
 
             # Set current parameters to algorithm.
             algorithm.set_params(eps=eps_mu,
@@ -596,16 +599,15 @@ def run(self, function, beta):
 
             # Obtain the gap from the last FISTA run. May be small and negative
             # close to machine epsilon.
-            gap_mu = abs(algorithm.info_get(Info.gap)[-1])
+            gap_mu = np.abs(algorithm.info_get(Info.gap)[-1])
             # TODO: Warn if gap_mu < -consts.TOLERANCE.
 
             if not self.simulation:
-                if gap_mu + mu * gM < self.eps:
+                converged = gap_mu + mu * gM < self.eps
+                if np.all(converged):
 
                     if self.info_requested(Info.converged):
-                        self.info_set(Info.converged, True)
-
-                    converged = True
+                        self.info_set(Info.converged, converged)
 
             if self.callback is not None:
                 self.callback(locals())
@@ -614,16 +616,16 @@ def run(self, function, beta):
                       "eps_mu: %g" % (i, gap_mu, eps, mu, eps_mu))
 
             # Stopping criteria.
-            if (converged or self.num_iter >= self.max_iter) \
+            if (np.all(converged) or self.num_iter >= self.max_iter) \
                     and self.num_iter >= self.min_iter:
                 break
 
             # Update the precision eps.
 #            eps = self.tau * (gap_mu + mu * gM)
-            eps = max(self.eps, self.tau * (gap_mu + mu * gM))
+            eps = np.maximum(self.eps, self.tau * (gap_mu + mu * gM))
             # Compute and update mu.
 #            mu = max(self.mu_min, min(function.mu_opt(eps), mu))
-            mu = min(function.mu_opt(eps), mu)
+            mu = np.minimum(function.mu_opt(eps), mu)
             function.set_mu(mu)
 
             i = i + 1

parsimony/datasets/simulate/grad.py

@@ -446,6 +446,7 @@ def _Nesterov_grad(beta, A, rng=RandomUniform(-1, 1), grad_norm=grad_l2):
     grad_Ab = 0
     for i in range(len(A)):
         Ai = A[i]
+        print("grad.py, _Nesterov_grad", Ai.shape, beta.shape)
         Ab = Ai.dot(beta)
         grad_Ab += Ai.T.dot(grad_norm(Ab, rng))
 

parsimony/datasets/utils.py

@@ -13,6 +13,39 @@
 from scipy import ndimage
 #import matplotlib.pyplot as plt
 
+def download_dataset(dataset):
+    """Download dataset.
+
+    Parameters
+    ----------
+    dataset: string
+
+    Return
+    ------
+    dataset_filename, dataset: path_to_archive, numpy_archive
+
+    Examples
+    --------
+
+    >>> dataset = "%s_%s_%ix%ix%i_%i_dataset_v%s.npz" % \
+    ...         tuple(["dice5", "classif", 50, 50, 1, 500, '0.3.1'])
+    >>> archiv, filename = download_dataset(dataset)
+    """
+    import os
+    import urllib
+    import tempfile
+    # base_ftp_url = "ftp://ftp.cea.fr/pub/dsv/anatomist/parsimony/%s"
+    base_ftp_url = "ftp://ftp.cea.fr/pub/unati/share/parsimony/datasets/%s"
+    tmp_dir = tempfile.gettempdir()
+    # dataset
+    dataset_url = base_ftp_url % dataset
+    dataset_filename = os.path.join(tmp_dir, os.path.basename(dataset_url))
+    if not os.path.exists(dataset_filename):
+        print("Download dataset from: %s => %s" % (dataset_url, dataset_filename))
+        urllib.request.urlretrieve(dataset_url, dataset_filename)
+    data = np.load(dataset_filename)
+    return(dataset_filename, data)
+
 
 def corr_to_coef(v_x, v_e, cov_xe, cor):
     """In a linear model y = bx + e. Calculate b such cor(bx + e, x) = cor.