Ignore:
Timestamp:
05/04/12 18:07:12 (2 years ago)
Author:
Lan Zagar <lan.zagar@…>
Branch:
default
Message:

Reimplemented lasso. Breaks compatibility.

It now uses a proximal gradient method for optimization instead of using scipy.optimize (see #1118).
The formulation is slightly different so there are new parameters (mainly lasso_lambda instead of t/s).
Improved some other things as well.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • Orange/regression/lasso.py

    r10580 r10859  
     1from numpy import dot, std, array, zeros, maximum, sqrt, sign, log, abs, \ 
     2                  ascontiguousarray, random as rnd 
     3from scipy.linalg import norm, eigh 
     4 
    15import Orange 
    2 import numpy 
    3  
    4 from Orange.regression import base 
    5  
    66from Orange.utils import deprecated_members, deprecated_keywords 
    77 
    8 def center(X): 
    9     """Centers the data, i.e. subtracts the column means. 
    10     Returns the centered data and the mean. 
    11  
    12     :param X: the data arry 
    13     :type table: :class:`numpy.array` 
     8 
     9def get_bootstrap_sample(data): 
     10    """Generate a bootstrap sample of a given data set. 
     11 
     12    :param data: the original data sample 
     13    :type data: :class:`Orange.data.Table` 
    1414    """ 
    15     mu = X.mean(axis=0) 
    16     return X - mu, mu 
    17  
    18 def standardize(X): 
    19     """Standardizes the data, i.e. subtracts the column means and divide by  
    20     standard deviation. 
    21     Returns the centered data, the mean, and deviations. 
    22  
    23     :param X: the data arry 
    24     :type table: :class:`numpy.array` 
    25     """ 
    26     mu = numpy.mean(X, axis=0) 
    27     std = numpy.std(X, axis=0) 
    28     return (X - mu) / std, mu, std 
    29  
    30 def get_bootstrap_sample(table): 
    31     """Generates boostrap sample from an Orange Example Table 
    32     and stores it in a new :class:`Orange.data.Table` object 
    33  
    34     :param table: the original data sample 
    35     :type table: :class:`Orange.data.Table` 
    36     """ 
    37     n = len(table) 
    38     bootTable = Orange.data.Table(table.domain) 
    39     for i in range(n): 
    40         id = numpy.random.randint(0, n) 
    41         bootTable.append(table[id]) 
    42     return bootTable 
    43  
    44 def permute_responses(table): 
    45     """ Permutes values of the class (response) variable. 
     15    n = len(data) 
     16    bootstrap = Orange.data.Table(data.domain) 
     17    for id in rnd.randint(0, n, n): 
     18        bootstrap.append(data[id]) 
     19    return bootstrap 
     20 
     21def permute_responses(data): 
     22    """Permute values of the class (response) variable. 
    4623    The independence between independent variables and the response 
    4724    is obtained but the distribution of the response variable is kept. 
    4825 
    49     :param table: the original data sample 
    50     :type table: :class:`Orange.data.Table` 
     26    :param data: Original data. 
     27    :type data: :class:`Orange.data.Table` 
    5128    """ 
    52     n = len(table) 
    53     perm = numpy.random.permutation(n) 
    54     permTable = Orange.data.Table(table.domain, table) 
    55     for i, ins in enumerate(table): 
    56         permTable[i].set_class(table[perm[i]].get_class()) 
    57     return permTable 
    58  
    59 class LassoRegressionLearner(base.BaseRegressionLearner): 
    60     """Fits the lasso regression model, i.e. learns the regression parameters 
    61     The class is derived from 
    62     :class:`Orange.regression.base.BaseRegressionLearner` 
    63     which is used for preprocessing the data (continuization and imputation) 
    64     before fitting the regression parameters 
    65  
     29    n = len(data) 
     30    perm = rnd.permutation(n) 
     31    perm_data = Orange.data.Table(data.domain, data) 
     32    for i, ins in enumerate(data): 
     33        perm_data[i].set_class(data[perm[i]].get_class()) 
     34    return perm_data 
     35 
     36class LassoRegressionLearner(Orange.regression.base.BaseRegressionLearner): 
     37    """Fits the lasso regression model using FISTA 
     38    (Fast Iterative Shrinkage-Thresholding Algorithm). 
    6639    """ 
    6740 
    68     def __init__(self, name='lasso regression', t=1, s=None, tol=0.001, \ 
    69                  n_boot=100, n_perm=100, imputer=None, continuizer=None): 
    70         """ 
    71         :param name: name of the linear model, default 'lasso regression' 
    72         :type name: string 
     41    def __init__(self, lasso_lambda=0.1, max_iter=20000, eps=1e-6, 
     42                 n_boot=0, n_perm=0, imputer=None, continuizer=None, 
     43                 name='Lasso'): 
     44        """ 
     45        :param lasso_lambda: Regularization parameter. 
     46        :type lasso_lambda: float 
     47 
     48        :param max_iter: Maximum number of iterations for 
     49                         the optimization method. 
     50        :type max_iter: int 
     51 
     52        :param eps: Stop optimization when improvements are lower than eps. 
     53        :type eps: float 
     54 
     55        :param n_boot: Number of bootstrap samples used for non-parametric 
     56                       estimation of standard errors. 
     57        :type n_boot: int 
     58 
     59        :param n_perm: Number of permuations used for non-parametric 
     60                       estimation of p-values. 
     61        :type n_perm: int 
     62 
     63        :param name: Learner name. 
     64        :type name: str 
    7365         
    74         :param t: tuning parameter, upper bound for the L1-norm of the 
    75             regression coefficients 
    76         :type t: float 
    77          
    78         :param s: An alternative way to specify the tuning parameter ``t``. 
    79             Here ``t`` is taken to be t = s * sum(abs(B)) where B are the 
    80             coefficients of an ordinary least square linear fit. ``t`` parameter is ignored if ``s`` is specified (by default it 
    81             is None). 
    82         :type s: float 
    83          
    84         :param tol: tolerance parameter, regression coefficients 
    85             (absoulute value) under tol are set to 0, 
    86             default=0.001 
    87         :type tol: float 
    88          
    89         :param n_boot: number of bootstrap samples used for non-parametric 
    90             estimation of standard errors 
    91         :type n_boot: int 
    92          
    93         :param n_perm: number of permuations used for non-parametric 
    94             estimation of p-values 
    95         :type n_perm: int 
    96          
    97         """ 
    98  
    99         self.name = name 
    100         self.t = t 
    101         self.s = s 
    102         self.tol = tol 
     66        """ 
     67        self.lasso_lambda = lasso_lambda 
     68        self.max_iter = max_iter 
     69        self.eps = eps 
    10370        self.n_boot = n_boot 
    10471        self.n_perm = n_perm 
    10572        self.set_imputer(imputer=imputer) 
    10673        self.set_continuizer(continuizer=continuizer) 
    107  
    108  
    109     def __call__(self, table, weight=None): 
    110         """ 
    111         :param table: data instances. 
    112         :type table: :class:`Orange.data.Table` 
    113         :param weight: the weights for instances. Default: None, i.e. 
    114             all data instances are eqaully important in fitting 
    115             the regression parameters 
    116         :type weight: None or list of Orange.feature.Continuous 
    117             which stores weights for instances 
     74        self.name = name 
     75 
     76    def get_lipschitz(self, X): 
     77        """Return the Lipschitz constant of :math:`\\nabla f`, 
     78        where :math:`f(w) = \\frac{1}{2}||Xw-y||^2`. 
     79        """ 
     80        n, m = X.shape 
     81        if n > m: 
     82            X = ascontiguousarray(X.T) 
     83        k = min(n, m) - 1 
     84        eigvals = eigh(dot(X, X.T), eigvals_only=True, eigvals=(k, k)) 
     85        return eigvals[-1] 
     86 
     87    def fista(self, X, y, l, lipschitz, w_init=None): 
     88        """Fast Iterative Shrinkage-Thresholding Algorithm (FISTA).""" 
     89        z = w_old = zeros(X.shape[1]) if w_init is None else w_init 
     90        t_old, obj_old = 1, 1e400 
     91        XT = ascontiguousarray(X.T) 
     92        for i in range(self.max_iter): 
     93            z -= 1. / lipschitz * dot(XT, dot(X, z) - y) 
     94            w = maximum(0, abs(z) - l / lipschitz) * sign(z) 
     95            t = (1 + sqrt(1 + 4 * t_old**2)) / 2 
     96            z = w + (t_old - 1) / t * (w - w_old) 
     97            obj = ((y - dot(X, w))**2).sum() + l * norm(w, 1)  
     98            if abs(obj_old - obj) / obj < self.eps: 
     99                stop += 1 
     100                if obj < obj_old and stop > log(i + 1): 
     101                    break 
     102            else: 
     103                stop = 0 
     104            w_old, t_old = w, t 
     105            obj_old = obj 
     106        return w 
     107 
     108    def __call__(self, data, weight=None): 
     109        """ 
     110        :param data: Training data. 
     111        :type data: :class:`Orange.data.Table` 
     112        :param weight: Weights for instances. Not implemented yet. 
    118113         
    119114        """ 
    120115        # dicrete values are continuized         
    121         table = self.continuize_table(table) 
     116        data = self.continuize_table(data) 
    122117        # missing values are imputed 
    123         table = self.impute_table(table) 
    124  
    125         domain = table.domain 
    126         X, y, w = table.to_numpy() 
    127         n, m = numpy.shape(X) 
    128  
    129         X, mu_x, sigma_x = standardize(X) 
    130         y, coef0 = center(y) 
    131  
    132         t = self.t 
    133  
    134         if self.s is not None: 
    135             beta_full, rss, _, _ = numpy.linalg.lstsq(X, y) 
    136             t = self.s * numpy.sum(numpy.abs(beta_full)) 
    137             print "t =", t 
    138  
    139         import scipy.optimize 
    140  
    141         # objective function to be minimized 
    142         objective = lambda beta: numpy.linalg.norm(y - numpy.dot(X, beta)) 
    143         # initial guess for the regression parameters 
    144         beta_init = numpy.random.random(m) 
    145         # constraints for the regression coefficients 
    146         cnstr = lambda beta: t - numpy.sum(numpy.abs(beta)) 
    147         # optimal solution 
    148         coefficients = scipy.optimize.fmin_cobyla(objective, beta_init, cnstr, iprint=0) 
    149  
    150         # set small coefficients to 0 
    151         def set_2_0(c): return c if abs(c) > self.tol else 0 
    152         coefficients = numpy.array(map(set_2_0, coefficients)) 
    153         coefficients /= sigma_x 
     118        data = self.impute_table(data) 
     119        domain = data.domain 
     120        # prepare numpy matrices 
     121        X, y, _ = data.to_numpy() 
     122        n, m = X.shape 
     123        coefficients = zeros(m) 
     124        std_errors = array([float('nan')] * m) 
     125        p_vals = array([float('nan')] * m) 
     126        # standardize y 
     127        coef0, sigma_y = y.mean(), y.std() + 1e-6 
     128        y = (y - coef0) / sigma_y 
     129        # standardize X and remove constant vars 
     130        mu_x = X.mean(axis=0) 
     131        X -= mu_x 
     132        sigma_x = X.std(axis=0) 
     133        nz = sigma_x != 0 
     134        X = ascontiguousarray(X[:, nz]) 
     135        sigma_x = sigma_x[nz] 
     136        X /= sigma_x 
     137        m = sum(nz) 
     138 
     139        # run optimization method 
     140        lipschitz = self.get_lipschitz(X) 
     141        l = 0.5 * self.lasso_lambda * n / m 
     142        coefficients[nz] = self.fista(X, y, l, lipschitz) 
     143        coefficients[nz] *= sigma_y / sigma_x 
     144 
     145        d = dict(self.__dict__) 
     146        d.update({'n_boot': 0, 'n_perm': 0}) 
    154147 
    155148        # bootstrap estimator of standard error of the coefficient estimators 
    156         # assumption: fixed t 
     149        # assumption: fixed regularization parameter 
    157150        if self.n_boot > 0: 
    158151            coeff_b = [] # bootstrapped coefficients 
    159152            for i in range(self.n_boot): 
    160                 tmp_table = get_bootstrap_sample(table) 
    161                 l = LassoRegressionLearner(t=t, n_boot=0, n_perm=0) 
    162                 c = l(tmp_table) 
     153                tmp_data = get_bootstrap_sample(data) 
     154                l = LassoRegressionLearner(**d) 
     155                c = l(tmp_data) 
    163156                coeff_b.append(c.coefficients) 
    164             std_errors_fixed_t = numpy.std(coeff_b, axis=0) 
    165         else: 
    166             std_errors_fixed_t = [float("nan")] * m 
    167  
    168         # permutation test to obtain the significance of the regression 
    169         #coefficients 
     157            std_errors[nz] = std(coeff_b, axis=0) 
     158 
     159        # permutation test to obtain the significance of 
     160        # the regression coefficients 
    170161        if self.n_perm > 0: 
    171162            coeff_p = [] 
    172163            for i in range(self.n_perm): 
    173                 tmp_table = permute_responses(table) 
    174                 l = LassoRegressionLearner(t=t, n_boot=0, n_perm=0) 
    175                 c = l(tmp_table) 
     164                tmp_data = permute_responses(data) 
     165                l = LassoRegressionLearner(**d) 
     166                c = l(tmp_data) 
    176167                coeff_p.append(c.coefficients) 
    177             p_vals = \ 
    178                    numpy.sum(abs(numpy.array(coeff_p)) > \ 
    179                              abs(numpy.array(coefficients)), \ 
    180                              axis=0) / float(self.n_perm) 
    181         else: 
    182             p_vals = [float("nan")] * m 
     168            p_vals[nz] = (abs(coeff_p) > abs(coefficients)).sum(axis=0) 
     169            p_vals[nz] /= float(self.n_perm) 
    183170 
    184171        # dictionary of regression coefficients with standard errors 
    185172        # and p-values 
    186         dict_model = {} 
     173        model = {} 
    187174        for i, var in enumerate(domain.attributes): 
    188             dict_model[var.name] = (coefficients[i], std_errors_fixed_t[i], p_vals[i]) 
     175            model[var.name] = (coefficients[i], std_errors[i], p_vals[i]) 
    189176 
    190177        return LassoRegression(domain=domain, class_var=domain.class_var, 
    191                                coef0=coef0, coefficients=coefficients, 
    192                                std_errors_fixed_t=std_errors_fixed_t, 
    193                                p_vals=p_vals, 
    194                                dict_model=dict_model, 
    195                                mu_x=mu_x) 
     178            coef0=coef0, coefficients=coefficients, std_errors=std_errors, 
     179            p_vals=p_vals, model=model, mu_x=mu_x) 
    196180 
    197181deprecated_members({"nBoot": "n_boot", 
     
    201185 
    202186class LassoRegression(Orange.classification.Classifier): 
    203     """Lasso regression predicts value of the response variable 
     187    """Lasso regression predicts the value of the response variable 
    204188    based on the values of independent variables. 
    205189 
     
    210194    .. attribute:: coefficients 
    211195 
    212         Regression coefficients, sotred in list.  
    213  
    214     .. attribute:: std_errors_fixed_t 
    215  
    216         Standard errors of the coefficient estimator for the fixed 
    217         tuning parameter t. The standard errors are estimated using 
    218         bootstrapping method. 
     196        Regression coefficients.  
     197 
     198    .. attribute:: std_errors 
     199 
     200        Standard errors of coefficient estimates for a fixed 
     201        regularization parameter. The standard errors are estimated 
     202        using the bootstrapping method. 
    219203 
    220204    .. attribute:: p_vals 
    221205 
    222         List of p-values for the null hypothesis that the regression 
    223         coefficients equal 0 based on non-parametric permutation test. 
    224  
    225     .. attribute:: dict_model 
    226  
    227         Statistical properties of the model stored in dictionary: 
     206        List of p-values for the null hypotheses that the regression 
     207        coefficients equal 0 based on a non-parametric permutation test. 
     208 
     209    .. attribute:: model 
     210 
     211        Dictionary with the statistical properties of the model: 
    228212        Keys - names of the independent variables 
    229213        Values - tuples (coefficient, standard error, p-value)  
     
    231215    .. attribute:: mu_x 
    232216 
    233         Sample mean of the all independent variables.     
     217        Sample mean of independent variables.     
    234218 
    235219    """ 
    236220    def __init__(self, domain=None, class_var=None, coef0=None, 
    237                  coefficients=None, std_errors_fixed_t=None, p_vals=None, 
    238                  dict_model=None, mu_x=None): 
     221                 coefficients=None, std_errors=None, p_vals=None, 
     222                 model=None, mu_x=None): 
    239223        self.domain = domain 
    240224        self.class_var = class_var 
    241225        self.coef0 = coef0 
    242226        self.coefficients = coefficients 
    243         self.std_errors_fixed_t = std_errors_fixed_t 
     227        self.std_errors = std_errors 
    244228        self.p_vals = p_vals 
    245         self.dict_model = dict_model 
     229        self.model = model 
    246230        self.mu_x = mu_x 
     231 
     232    def _miss_2_0(self, x): 
     233        return x if x != '?' else 0 
    247234 
    248235    @deprecated_keywords({"resultType": "result_type"}) 
    249236    def __call__(self, instance, result_type=Orange.core.GetValue): 
    250237        """ 
    251         :param instance: data instance for which the value of the response 
    252             variable will be predicted 
    253         :type instance:  
     238        :param instance: Data instance for which the value of the response 
     239                         variable will be predicted. 
     240        :type instance: :obj:`Orange.data.Instance` 
    254241        """ 
    255242        ins = Orange.data.Instance(self.domain, instance) 
    256         if "?" in ins: # missing value -> corresponding coefficient omitted 
    257             def miss_2_0(x): return x if x != "?" else 0 
    258             ins = map(miss_2_0, ins) 
    259             ins = numpy.array(ins)[:-1] - self.mu_x 
     243        if '?' in ins: # missing value -> corresponding coefficient omitted 
     244            ins = map(self._miss_2_0, ins) 
     245            ins = array(ins)[:-1] - self.mu_x 
    260246        else: 
    261             ins = numpy.array(ins.native())[:-1] - self.mu_x 
    262  
    263         y_hat = numpy.dot(self.coefficients, ins) + self.coef0 
     247            ins = array(ins.native())[:-1] - self.mu_x 
     248 
     249        y_hat = dot(self.coefficients, ins) + self.coef0 
    264250        y_hat = self.class_var(y_hat) 
    265251        dist = Orange.statistics.distribution.Continuous(self.class_var) 
    266         dist[y_hat] = 1.0 
     252        dist[y_hat] = 1. 
    267253        if result_type == Orange.core.GetValue: 
    268254            return y_hat 
     
    274260    @deprecated_keywords({"skipZero": "skip_zero"}) 
    275261    def to_string(self, skip_zero=True): 
    276         """Pretty-prints Lasso regression model, 
     262        """Pretty-prints a lasso regression model, 
    277263        i.e. estimated regression coefficients with standard errors 
    278         and significances. Standard errors are obtained using bootstrapping 
    279         method and significances by the permuation test 
    280  
    281         :param skip_zero: if True variables with estimated coefficient equal to 0 
    282             are omitted 
    283         :type skip_zero: boolean 
    284         """ 
    285  
    286         from string import join 
     264        and significances. Standard errors are obtained using the 
     265        bootstrapping method and significances by a permuation test. 
     266 
     267        :param skip_zero: If True, variables with estimated coefficient 
     268                          equal to 0 are omitted. 
     269        :type skip_zero: bool 
     270        """ 
    287271        labels = ('Variable', 'Coeff Est', 'Std Error', 'p') 
    288         lines = [join(['%10s' % l for l in labels], ' ')] 
    289  
    290         fmt = "%10s " + join(["%10.3f"] * 3, " ") + " %5s" 
    291         fmt1 = "%10s %10.3f" 
     272        lines = [' '.join(['%10s' % l for l in labels])] 
     273 
     274        fmt = '%10s ' + ' '.join(['%10.3f'] * 3) + ' %5s' 
     275        fmt1 = '%10s %10.3f' 
    292276 
    293277        def get_star(p): 
    294             if p < 0.001: return  "*"*3 
    295             elif p < 0.01: return "*"*2 
    296             elif p < 0.05: return "*" 
    297             elif p < 0.1: return  "." 
    298             else: return " " 
     278            if p < 0.001: return  '*' * 3 
     279            elif p < 0.01: return '*' * 2 
     280            elif p < 0.05: return '*' 
     281            elif p < 0.1: return  '.' 
     282            else: return ' ' 
    299283 
    300284        stars = get_star(self.p_vals[0]) 
     
    307291            stars = get_star(self.p_vals[i]) 
    308292            lines.append(fmt % (self.domain.attributes[i].name, 
    309                          self.coefficients[i], self.std_errors_fixed_t[i], 
     293                         self.coefficients[i], self.std_errors[i], 
    310294                         self.p_vals[i], stars)) 
    311         lines.append("Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1 empty 1") 
    312         lines.append("\n") 
     295        lines.append('Signif. codes:  0 *** 0.001 ** 0.01 * 0.05 . 0.1 empty 1\n') 
    313296        if skip_zero: 
    314297            k = len(skipped) 
    315298            if k == 0: 
    316                 lines.append("All variables have non-zero regression coefficients. ") 
     299                lines.append('All variables have non-zero regression coefficients.') 
    317300            else: 
    318                 suff = "s" if k > 1 else "" 
    319                 lines.append("For %d variable%s the regression coefficient equals 0: " \ 
    320                       % (k, suff)) 
    321                 for var in skipped: 
    322                     lines.append(var) 
    323         return "\n".join(lines) 
     301                lines.append('For %d variable%s the regression coefficient equals 0:' 
     302                             % (k, 's' if k > 1 else '')) 
     303                lines.append(', '.join(var for var in skipped)) 
     304        return '\n'.join(lines) 
    324305 
    325306    def __str__(self): 
     
    327308 
    328309deprecated_members({"muX": "mu_x", 
    329                     "stdErrorsFixedT": "std_errors_fixed_t", 
     310                    "stdErrorsFixedT": "std_errors", 
    330311                    "pVals": "p_vals", 
    331                     "dictModel": "dict_model"}, 
     312                    "dictModel": "model"}, 
    332313                   wrap_methods=["__init__"], 
    333314                   in_place=True)(LassoRegression) 
    334315 
    335316if __name__ == "__main__": 
    336  
    337     import Orange 
    338  
    339     table = Orange.data.Table("housing.tab") 
    340  
    341     c = LassoRegressionLearner(table, t=len(table.domain)) 
     317    data = Orange.data.Table('housing') 
     318    c = LassoRegressionLearner(data) 
    342319    print c 
Note: See TracChangeset for help on using the changeset viewer.