Implement complex lasso in progress ......

Author: cxzheng

sklearn/linear_model/cd_fast.pyx

@@ -874,3 +874,95 @@ def enet_coordinate_descent_multi_task(floating[::1, :] W, floating l1_reg,
                     break
 
     return np.asarray(W), gap, tol, n_iter + 1
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.cdivision(True)
+def enet_coordinate_descent_complex(floating[::1, :] W, floating l1_reg,
+                                    floating l2_reg,
+                                    floating[::1, :] Xr,
+                                    floating[::1, :] Xi,
+                                    floating[::1, :] Y,
+                                    int max_iter, floating tol, object rng,
+                                    bint random=0):
+    """Cython version of the coordinate descent algorithm
+        for Elastic-Net mult-task regression in complex domain
+
+    """
+
+    # fused types version of BLAS functions
+    if floating is float:
+        dtype = np.float32
+        gemv = sgemv
+        dot = sdot
+        copy = scopy
+    else:
+        dtype = np.float64
+        gemv = dgemv
+        dot = ddot
+        copy = dcopy
+
+    # get the data information into easy vars
+    cdef unsigned int n_samples = Xr.shape[0]
+    cdef unsigned int n_features = Xr.shape[1]
+
+    # initial value of the residuals
+    cdef floating[::1] Rr = np.empty(n_samples, dtype=dtype)
+    cdef floating[::1] Ri = np.empty(n_samples, dtype=dtype)
+
+    cdef floating[:] w_ii = np.zeros(2, dtype=dtype)
+    cdef unsigned int ii
+    cdef unsigned int jj
+    cdef unsigned int n_iter = 0
+    cdef unsigned int f_iter
+    cdef UINT32_t rand_r_state_seed = rng.randint(0, RAND_R_MAX)
+    cdef UINT32_t* rand_r_state = &rand_r_state_seed
+
+    cdef floating* W_ptr = &W[0, 0]
+    cdef floating* Y_ptr = &Y[0, 0]
+
+    if l1_reg == 0:
+        warnings.warn("Coordinate descent with l1_reg=0 may lead to unexpected"
+            " results and is discouraged.")
+
+    with nogil:
+        # Compute Rr and Ri: real and imaginary parts of the residual
+        # real part: Yr - np.dot(Xr, Wr) + np.dot(Xi, Wi)
+        copy(n_samples, Y_ptr, 1, &Rr[0], 1)
+        gemv(CblasColMajor, CblasNoTrans,
+             n_samples, n_features, -1.0, &Xr[0, 0], n_samples,
+             W_ptr, 2, 1.0, &Rr[0], 1)
+        gemv(CblasColMajor, CblasNoTrans,
+             n_samples, n_features, 1.0, &Xi[0, 0], n_samples,
+             W_ptr + 1, 2, 1.0, &Rr[0], 1)
+
+        # imaginary part:
+        # real part: Yr - np.dot(Xr, Wi) - np.dot(Xi, Wr)
+        copy(n_samples, Y_ptr + n_samples, 1, &Ri[0], 1)
+        gemv(CblasColMajor, CblasNoTrans,
+             n_samples, n_features, -1.0, &Xr[0, 0], n_samples,
+             W_ptr + 1, 2, 1.0, &Ri[0], 1)
+        gemv(CblasColMajor, CblasNoTrans,
+             n_samples, n_features, -1.0, &Xi[0, 0], n_samples,
+             W_ptr, 2, 1.0, &Ri[0], 1)
+
+        # tol = tol * linalg.norm(Y, ord='fro') ** 2
+        tol = tol * dot(n_samples * 2, Y_ptr, 1, Y_ptr, 1)
+
+        for n_iter in range(max_iter):
+            for f_iter in range(n_features):  # Loop over coordinates
+                # select a coordinate
+                if random:
+                    ii = rand_int(n_features, rand_r_state)
+                else:
+                    ii = f_iter
+
+                # w_ii = W[:, ii] # Store previous value
+                w_ii[0] = W[0, ii]
+                w_ii[1] = W[1, ii]
+
+                # if np.sum(w_ii ** 2) != 0.0:  # can do better
+                if w_ii[0] != 0.0 or w_ii[1] != 0.0:
+                    # Remove contributions of w_ii from R
+
+                # prepare for the soft-thresholding

sklearn/linear_model/coordinate_descent.py

@@ -2260,3 +2260,87 @@ def __init__(self, eps=1e-3, n_alphas=100, alphas=None, fit_intercept=True,
             max_iter=max_iter, tol=tol, copy_X=copy_X,
             cv=cv, verbose=verbose, n_jobs=n_jobs, random_state=random_state,
             selection=selection)
+
+###############################################################################
+# Complex ElasticNet and Lasso models 
+
+class ComplexElasticNet(ElasticNet):
+    """ElasticNet model in complex domain trained with L1/L2 mixed-norm as regularizer
+    """
+
+    def __init__(self, alpha=1.0, l1_ratio=0.5, max_iter=1000, 
+                 tol=1e-4, warm_start=False, random_state=None, 
+                 selection='cyclic'):
+        self.l1_ratio = l1_ratio
+        self.alpha = alpha
+        self.max_iter = max_iter
+        self.tol = tol
+        self.warm_start = warm_start
+        self.random_state = random_state
+        self.selection = selection
+
+    def fit(self, Xr, Xi, Y):
+        """Fit MultiTaskElasticNet model with coordinate descent
+
+        Parameters
+        -----------
+        Xr : ndarray, shape (n_samples, n_features)
+             real part of the X matrix
+        Xi : ndarray, shape (n_samples, n_features)
+             imaginary part of the X matrix
+        Y  : ndarray, shape (n_samples, 2)
+             Target. Will be cast to X's dtype if necessary
+             Y[:,0] is the real part and Y[:,1] is the imaginary part
+        """
+        Xr = check_array(Xr, dtype=[np.float64, np.float32], order='F', copy=False)
+        Xi = check_array(Xi, dtype=[np.float64, np.float32], order='F', copy=False)
+        Y = check_array(Y, dtype=X.dtype.type, order='F', copy=False)
+
+        n_samples, n_features = X.shape
+
+        if n_samples != Y.shape[0]:
+            raise ValueError("X and Y have inconsistent dimensions (%d != %d)"
+                             % (n_samples, y.shape[0]))
+        if Y.shape[1] != 2:
+            raise ValueError("Y should have %d columns, but two are expected (real and imaginary part)"
+                             % y.shape[1])
+
+        if not self.warm_start or self.coef_ is None:
+            self.coef_ = np.zeros((2, n_features), dtype=X.dtype.type, order='F')
+
+        l1_reg = self.alpha * self.l1_ratio * n_samples
+        l2_reg = self.alpha * (1.0 - self.l1_ratio) * n_samples
+
+        self.coef_ = np.asfortranarray(self.coef_)  # coef contiguous in memory
+
+        if self.selection not in ['random', 'cyclic']:
+            raise ValueError("selection should be either random or cyclic.")
+        random = (self.selection == 'random')
+
+        # now call coordinate descient
+        self.coef_, self.dual_gap_, self.eps_, self.n_iter_ = \
+            cd_fast.enet_coordinate_descent_complex(
+                self.coef_, l1_reg, l2_reg, Xr, Xi, Y, self.max_iter, self.tol,
+                check_random_state(self.random_state), random)
+
+        if self.dual_gap_ > self.eps_:
+            warnings.warn('Objective did not converge, you might want'
+                          ' to increase the number of iterations',
+                          ConvergenceWarning)
+
+        # return self for chaining fit and predict calls
+        return self
+
+class ComplexLasso(ComplexElasticNet):
+
+    def __init__(self, alpha=1.0, max_iter=1000, tol=1e-4, 
+                 warm_start=False, random_state=None, selection='cyclic'):
+        self.l1_ratio = 1.0
+        self.alpha = alpha
+        self.max_iter = max_iter
+        self.tol = tol
+        self.warm_start = warm_start
+        self.random_state = random_state
+        self.selection = selection
+
+###############################################################################