Since the recent update of scipy library version 1.8.0 (February 5th), I'm encountering a bug in my usual pylearn-parsimony based routine. It happens when fitting a Conesta solver with a non-zero total variation parameter. The error occurred in the svds function of scipy.sparse.linalg, called by nipals.py#343. So I needed to downgrade scipy.
I am not sure if I should report the issue here.

I want to apply LinearRegressionL2SmoothedL1TV into my dataset:

import parsimony.estimators as estimators
import parsimony.algorithms as algorithms
import parsimony.functions.nesterov.l1tv as l1tv
from util import Dataset
import parsimony.algorithms.primaldual as primaldual
import numpy as np
X=pd.read_csv('/home/lemma/train.csv')
y=pd.read_csv('/home/lemma/y.csv')
y=np.log(y)
l1=0.1
l2=0.9
tv=1.0
mask=np.zeros(X.shape)
mask[1:(X.shape[0] - 1), 0:(X.shape[1] - 1)] = 1
A = l1tv.linear_operator_from_mask(mask, X.shape[1], penalty_start=0)
lr = estimators.LinearRegressionL2SmoothedL1TV(l2, l1, tv, A,
algorithm=primaldual.ExcessiveGapMethod(),
algorithm_params=dict(max_iter=1000),
mean=False)
res = lr.fit(X, y)
y_pred=res.predict(X)
ValueError: operands could not be broadcast together with shapes (14,1) (2448108,1) (14,1)

Not sure what's going on?
My train data shape is (188314,14) and y_train is (188314,)

The problem occurs when max_iter is very large ie.: 50000

1. The Outer FISTA max_iter is 50000,
2. The Outer FISTA set the max_iter of the inner FISTA to the same value "max_iter" = 50000
3. According to the theory the eps of the inner loop should decrease with the iteration of the outer loop (=1.0 / (float(i) ** (4.0 + consts.FLOAT_EPSILON)))

But the problem is that the inner loop always hit max_iter=50000, which make it VERY slow (hours) for a 50 x 50 x 1 x 300(samples) dataset

Tommy, Fouad, what should we do ?

1. Limiting the number of iteration in the inner loop ? This breaks the theory.

2. Limiting the eps in the inner loop: which also breaks the theory. If we chose solution I suggest to limit eps to a relevant numerical precision: consts.FLOAT_EPSILON

so NesterovFunction.prox(..., eps) will start with something like:
eps = max(consts.FLOAT_EPSILON, eps)

Do you have any other idea ? We will have to argue that choice in the paper. I think the solution 2
eps = max(consts.FLOAT_EPSILON, eps) is easily arguable, since it is useless of dig further than consts.FLOAT_EPSILON.

I have install pylearn-parsimony accordingly but when I run the example code for elastic net:

import parsimony.estimators as estimators
alpha = 0.1  # global penalty
l = 0.1  # l1 ratio (lasso)
enet_estimator = estimators.ElasticNet(l=l, alpha=alpha)
enet_estimator.fit(X, y)
print("Estimated beta error =", np.linalg.norm(enet_estimator.beta - beta))

It then returns:

/usr/local/lib/python3.6/dist-packages/pylearn_parsimony-0.3.1-py3.6.egg/parsimony/config.py:55: RuntimeWarning: Could not locate the config file.
Estimated beta error = 0.8064548347250173
/usr/local/lib/python3.6/dist-packages/pylearn_parsimony-0.3.1-py3.6.egg/parsimony/config.py:66: RuntimeWarning: Could not locate the config file.

Hello,

I've been trying to use your LogisticRegressionL1L2GL for a binary classification problem. I've done so using the scikit-learn Pipeline, in which the previous steps, before using the LogisticRegressionL1L2GL, are a countVectorizer and a tfidftransformer. To be able to utilize your LogisticRegressionL1L2GL in this setup I have to convert the csr_matrix format that the previous step outputs into an array to be used in the X attribute in the fit and predict methods, which I do. However, the results of the classification is always of one class. I've looked into the code and I've found that the probability given to each instance is 0.5 and as such the classification given is 1.0.

Could I be doing something wrong, and if so, is there any insight you could provide?

Thanks in advance

Problem lambda max of the l1tv linear operator is defined twice:

• in l1tv.lambda_max()
• in NesterovFunction.prox.F.step()

Moreover the lambda_max is calculated again for each new inner inexactFista loop (in NesterovFunction.prox.F.step)

Proposition : in NesterovFunction.prox.F.step call l1tv.lambda_max (of A'A) passing the function to F similarly to proj. Then multiply by t^2. Since the step is changing. See p47 of ols tv paper.

=> lambda_max would be computed only once over the whole Fista-fista loops, since it is stored.
=> code to compute lambda_max is factorized

Tommy: does is sound right to you ?

Documentation should be verified and completed.

• To check: the Estimators section (some errors)
  • GL
  • LogisticRegression
• To be completed: the Algorithms section
  • Precise what are the default algorithms for the estimators
  • FISTA (use estimator with FISTA)

I realized that we could make an impressive gain of computation time by vectorizing the computation to simultaneously solve several minimization problems. Instead of having on l1, l2 tv scalars values we could provide vectors of length k to simultaneously solve k minimization problems. If we provide l1 = [0.1, 0.9], l2=[0.9, 0.1], tv = [0.2, 0.8], we would solve two problems: one with l1, l2, tv = 0.1, 0.9, 0.2 and one with l1, l2, tv = 0.9, 0.1, 0.8. We would simultaneously estimate k beta vectors and model.predict(X) will then provide k prediction. This would be the only consequence on the global API.

First, this is useful for multinomial classification problems.

Second, such framework should considerably speed up the computation over a grid of parameters. The rationale is the following, most of the computation is spent in the np.dot(X, beta). My hypothesis was that the time-consuming procedure is to load X data int the L3 cache of the processor. I realized that by observing a considerable slowdown of the computation when using several cores of the same processor. Using the 8 cores of the same processor slow down the minimization by a factor of 5! It is a very poor scaling! So if the same X data is used for several problems (beta is p x k array instead of being a p x 1) we could simultaneously solve several problems with a minimum overhead. Test it with the following code:

import time
import matplotlib.pylab as plt

n_features = int(550**2)
n_samples = 200

X = np.random.randn(n_samples, n_features)
beta = np.random.randn(n_features, 1)
y = np.dot(X, beta)

sizes = np.arange(1, 101)
elapsed = np.zeros(sizes.shape[0])


for s in sizes:
    beta = np.random.randn(X.shape[1], s)
    t_ = time.clock()

    # regression loss pattern
    grad = np.dot(X.T, np.dot(X, beta) - y)
    elapsed[s-1] = time.clock() - t_

# plt.plot(sizes, elapsed)
#plt.plot([1, sizes[-1]], [elapsed[0], elapsed[0]*sizes[-1]])

plt.plot(sizes, elapsed / elapsed[0])
plt.xlabel("beta.shape[1] (in X * beta)")
plt.ylabel("Ratio of CPU time: time(X beta) / time(X beta.shape[1]==1) ")

I observed an acceleration factor of 10 (anaconda with MKL single thread). Solving 50 problems beta.shape = [p, 50] is only 5 times slower than solving a single problem.

Most of the code should smoothly move to a vectorized version since beta is a [p x 1] vector and can be easly extend to a [p x k]. See an example of FISTA:

z = betanew + ((i - 2.0) / (i + 1.0)) * (betanew - betaold)
step = function.step(z)

betaold = betanew
betanew = function.prox(z - step * function.grad(z),
                                    step,
                                    eps=1.0 / (float(i) ** (4.0 + consts.FLOAT_EPSILON)),
                                    max_iter=self.max_iter)

Most of change will ocure when testing for convergences:
gap_mu + mu * gM < self.eps will become np.min(gap_mu + mu * gM) < self.eps

I suggest to start with FISTA and elasticnet, and evaluated the acceleration.

After doing a git clone and then running python setup.py install I get many errors. Should I be able to install this on Windows? Trace is below.

D:\Source\pylearn-parsimony>python setup.py install
C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\site-packages\setuptools-20.7.0-py2.7.egg\setuptools\dist.py:294: UserWarning: The version specified ('0
x') is an invalid version, this may not work as expected with newer versions of setuptools, pip, and PyPI. Please see PEP 440 for more details.
running install
running bdist_egg
running egg_info
creating pylearn_parsimony.egg-info
writing requirements to ./pylearn_parsimony.egg-info\requires.txt
writing ./pylearn_parsimony.egg-info\PKG-INFO
writing top-level names to ./pylearn_parsimony.egg-info\top_level.txt
writing dependency_links to ./pylearn_parsimony.egg-info\dependency_links.txt
writing manifest file './pylearn_parsimony.egg-info\SOURCES.txt'
Traceback (most recent call last):
File "setup.py", line 84, in
setup(**params)
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\distutils\core.py", line 151, in setup
dist.run_commands()
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\distutils\dist.py", line 953, in run_commands
self.run_command(cmd)
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\distutils\dist.py", line 972, in run_command
cmd_obj.run()
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\site-packages\setuptools-20.7.0-py2.7.egg\setuptools\command\install.py", line 67, in run
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\site-packages\setuptools-20.7.0-py2.7.egg\setuptools\command\install.py", line 109, in do_egg_in
ll
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\distutils\cmd.py", line 326, in run_command
self.distribution.run_command(command)
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\distutils\dist.py", line 972, in run_command
cmd_obj.run()
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\site-packages\setuptools-20.7.0-py2.7.egg\setuptools\command\bdist_egg.py", line 152, in run
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\distutils\cmd.py", line 326, in run_command
self.distribution.run_command(command)
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\distutils\dist.py", line 972, in run_command
cmd_obj.run()
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\site-packages\setuptools-20.7.0-py2.7.egg\setuptools\command\egg_info.py", line 186, in run
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\site-packages\setuptools-20.7.0-py2.7.egg\setuptools\command\egg_info.py", line 209, in find_sou
s
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\site-packages\setuptools-20.7.0-py2.7.egg\setuptools\command\egg_info.py", line 293, in run
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\site-packages\setuptools-20.7.0-py2.7.egg\setuptools\command\egg_info.py", line 322, in add_defa
s
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\site-packages\setuptools-20.7.0-py2.7.egg\setuptools\command\sdist.py", line 120, in add_default
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\distutils\cmd.py", line 312, in get_finalized_command
cmd_obj.ensure_finalized()
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\distutils\cmd.py", line 109, in ensure_finalized
self.finalize_options()
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\site-packages\setuptools-20.7.0-py2.7.egg\setuptools\command\build_py.py", line 33, in finalize_
ions
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\distutils\command\build_py.py", line 56, in finalize_options
self.package_dir[name] = convert_path(path)
File "C:\Users\jennl\AppData\Local\Continuum\Anaconda2\lib\distutils\util.py", line 126, in convert_path
raise ValueError, "path '%s' cannot end with '/'" % pathname
ValueError: path './' cannot end with '/'

When running the block
k = 0.0 # l2 ridge regression coefficient
l = 0.0 # l1 lasso coefficient
g = 1.0 # tv coefficient
A, n_compacts = tv.linear_operator_from_shape(shape) # Memory allocation for TV
olstv = estimators.LinearRegressionL1L2TV(k, l, g, A, mu=0.0001,
algorithm=algorithms.explicit.FISTA(),
algorithm_params=dict(max_iter=1000))

you end up with the following error in the call to LinearRegressionL1L2TV:
AttributeError: 'module' object has no attribute 'explicit'

Apparently, algorithms no longer has an 'explici't attribute

Please fix it!

1. Harmonize linear operator construction helper: some helper return A, some return A, n_compact. Proposition all helper return only A

2. linear operator is actually a simple list. We would like to:

• store information like the Lipschitz constant (that could be computed only once and reused) and eventually the n_compacts
• Provide serialization to disk. This would also also a lot of computation time

I propose to create an object LinearOperator which is a list (inherits from a list) this way it will not impact the code. I propose to add, at least, 3 methods:

• L() return the Lipschitz constant, could be None
• load(filename): npz format not pickle !
• save(filename)

Question should we move from helper to statics constructors ?

Currently this:

    def set_params(self, **kwargs):
        """Set the given input parameters in the estimator.
        """
        for k in kwargs:
            self.__setattr__(k, kwargs[k])

Should be this:

    def set_params(self, **kwargs):
        """Set the given input parameters in the estimator.
        """
        for k in kwargs:
            self.__setattr__(k, kwargs[k])
        return self

I actually think you can probably just delete the set_params method as there is no need to overwrite it (but guessing I'm wrong given the sklearn developers involved in parsimony!)

#43 fixes few failures, but some remain:

========================================================================= short test summary info =========================================================================
FAILED tests/test_linear_regression.py::TestLinearRegression::test_large - ValueError: not enough values to unpack (expected 2, got 1)
FAILED tests/test_linear_regression_large.py::test_weights_vs_sklearn - AttributeError: 'Ridge' object has no attribute 'coef_'
FAILED tests/test_logistic_regression.py::TestLogisticRegression::test_large - ValueError: not enough values to unpack (expected 2, got 1)
FAILED tests/test_logistic_regression_large.py::test_weights_vs_sklearn - AttributeError: 'LogisticRegression' object has no attribute 'coef_'
FAILED tests/test_multiblock.py::TestMultiblock::test_multiblock_fista - AttributeError: 'MultiblockFunctionWrapper' object has no attribute 'at_point'
FAILED tests/test_svd.py::TestSVD::test_fast_sparse_svd - ValueError: not enough values to unpack (expected 2, got 1)
FAILED tests/test_svd.py::TestSVD::test_svd_adjacency_laplacian_coincide - AssertionError: 44.91093483101495 != 44.91093483101497

The "has no attribute coef_" hints at sklearn estimators not being fitted.

Also, the tests take 5 minutes to run. Launching solvers on smaller data may make them faster and thus improve the feedback loop when running the tests.

Some properties of pylearn-parsimony Estimators are incompatible with sklearn Estimators. For instance in sklearn all parameters of the constructor have a default value (this allow to default construct the Estimator). There are many guidelines about the Estimator interface here (see chapters "APIs of scikit-learn objects" and "Rolling your own estimator"). Respecting this interface could help to mix parsimony and sklearn Estimators.

There are some tests in sklearn to ensure that all Estimators respect some assumptions of the API. One option to ensure the compatibility would be to use the same tests.

In my case, the group size in group lasso in algorithm is not unified and for different group it may should use different penalty (like $\sqrt(p_j)$ where $p_j$ is the degree of freedom of each group).

However, it seems that in your algorithm the penalty coefficients for groups is always the same, which may lead to some inconsistency.

The duality gap may infinite when the l2 penalty coefficient is null. Thus, we need to use different dual variable to keep it finite while maintaining the other key properties, this leads to a different expression.

While trying to use this module with a recent version of Python and numpy, scipy, and understand the error underlying desanou/multiscale_glasso#1, I got the following error:

Module 'numpy' has no attribute 'float'.np.floatwas a deprecated alias for the builtinfloat. To avoid this error in existing code, use floatby itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, usenp.float64` here. The aliases was originally deprecated in NumPy 1.20; for more details and guidance see the original release note at:

https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations

It seems that parsimony use at least once np.float :

Do you plan to publish a release correcting this?

Thanks

80% of the computation time is spent in the np.dot(X, beta) of the loss and it is done twice in the gradient and in the gap inside the same iteration. We could cache this result. Memory cache solution such as joblib exists (the Memory class) but it provides complex and unnecessary features:, ie, storing on disk, long term memorizing that will slow down the performances.
We only need a "short" time memory with only need a one step back cache. Here are my proposed specifications:

• One step back cache mechanisms: only the last result is memorized avoiding the computation with a second call with the same parameters.
• We assume that all parameters can be converted into numpy array (this is for sake of simplicity). Such constraint can be released latter using a more complex mechanism.
• Equality of parameters are tested based on the values. The same memory chunk that has been modified will fail and different memory chunks with same values will success.
• The check of some parameters can be avoided, typically, the constant large dataset or the linear operator.

The pros arguments are

• Simplicity
• The gain in execution time, that should me measured, but I assume between 30% and 40% of gain if 80% is spent in np.dot(X, beta).
• Negligible impact on the code of Parsimony. The cache object will be instanciated in the constructor of the function to be minimized, (eg, LinearRegressionL1L2TV, etc. ). The cache instance will be given to the loss. Then, only a few lines where np.dot(X, beta) will be touched.

Bellow an example of implementation

import numpy as np
import hashlib

class CacheOneStepBack:
    def __init__(self):
        self.memory = dict()

    def cache(self, func, *args, checkargs=None):
        """Cache result of func(*args). This is one step back memory ie, only the
        last result is memorized avoiding the computation with a second call with
        the same parameters. Parameters a compared based of their values. Assume
        args can be converted into numpy arrays.

        Parameters
        ----------
        func : Python function

        *args: Arguments of the python function

        checkargs: Optional, list of index of argument to be checked. Default is
            None, ie, all parameters are checked. Use this to avoid checking of
            aparameters that are known to be constant, ie, large data matrix.
            This should be used carefully because it raises risks of errors.

        """
        if checkargs is None:
            checkargs = list(range(len(args)))
        key = str(func)
        hashcodes = [hashlib.sha1(np.asarray(args[i])).hexdigest()
            if i in checkargs else None for i in range(len(args))]

        if not key in self.memory:
            self.memory[key] = [hashcodes, func(*args)]
        else:
            #same = [CACHE_ONE_STEP_BACK[key][0][i] == hashs[i] for i in range(len(hashs))]
            same = [self.memory[key][0][i] == hashcodes[i] for i in range(len(args))]
            if not np.all(same):
                self.memory[key] = [hashcodes, func(*args)]
        return self.memory[key][1]


memory = CacheOneStepBack()

# Test with floats and integers
# All arguments are checked, no risk of confusion.
def add(a, b, c):
    print("Call add", a, b, c)
    return a + b + c

memory.cache(add, 1.0, 2.0, 3) # first computation
memory.cache(add, 1.0, 2.0, 3) # Use cached results
memory.cache(add, 1.0, 2.0, 3.0)  # recompute hash 3.0 != hash 3
memory.cache(add, 1., 2., 4) # re-compute

# Test with numpy arrays
X = np.random.rand(150, int(5e5))
beta = np.random.rand(int(5e5), 1)
betaref = beta
betacpy = beta.copy()


# Here we avoid checking X
%time a1 = memory.cache(np.dot, X, beta, checkargs=[1]) # first computation (~454 ms)
%time a2 = memory.cache(np.dot, X, betaref, checkargs=[1]) # use cache, same data chunck (~12.2 ms)
%time a3 = memory.cache(np.dot, X, betacpy, checkargs=[1]) # use cache (~12.5 ms)
assert np.all(a1 == a2) and np.all(a1 == a3)
beta[0] = 33
%time a1 = memory.cache(np.dot, X, beta, checkargs=[1]) # re-compute (~454 ms)
%time a2 = memory.cache(np.dot, X, betaref, checkargs=[1]) # use cache, same data chunck (~12.2 ms)
%time a2 = memory.cache(np.dot, X, betacpy, checkargs=[1]) # recompute (~387 ms)

This is a question about the simulation parameter in the classes available in proximal.py.

In practice, what happen when simulation = True ? Is there any possibility to return the estimated regression vector for each iteration of the optimization alogorihms ?

When tying to add a unit test for the static Conesta, I noticed that using the gap is a very good idea to test whether the algorithm converged correctly. This is especially useful when testing the logistic regression model since in this case no explicit solution is known. So, I added Info.gap to the satic_Conesta 's INFO_PROVIDED in the proximal file.

I was studying the code of L1 regularization and I noticed that you use a norm with the maximum of the columns when it is applied to a matrix. I didn't understand why this, and I think that the subgradient and proximal is not suited to this type of norm.

Sorry if I said something wrong, I just didn't understand completely.

Thanks

It is currently not possible to use the scikit-learn GridSearchCV module with estimators from Parsimony because the method get_params is not compatible:

import parsimony.estimators as estimators
from sklearn.model_selection import GridSearchCV

rr = estimators.RidgeRegression(l=1.0)
param = {'ridgeregression__l': [0.2 * k for k in range(6)]}
grid_search = GridSearchCV(estimator=rr, param_grid=param, cv=5)

grid_search.fit(Xtr, ytr)

Traceback (most recent call last):

  File "<ipython-input-51-4f6ebd276276>", line 8, in <module>
    grid_search.fit(Xtr, ytr)

  File "/home/python-environments/env/lib/python2.7/site-packages/sklearn/model_selection/_search.py", line 625, in fit
    base_estimator = clone(self.estimator)

  File "/home/python-environments/env/lib/python2.7/site-packages/sklearn/base.py", line 60, in clone
    new_object_params = estimator.get_params(deep=False)

TypeError: get_params() got an unexpected keyword argument 'deep'

This could be fixed adding the parameter deep and as in scikit-learn:

https://github.com/scikit-learn/scikit-learn/blob/f0ab589f1541b1ca4570177d93fd7979613497e3/sklearn/base.py#L244-L248

how can I use meshes for LinearRegressionL1L2TV/ LinearRegressionL1L2GL?
when I use linear_operator from_meshes and give its output as A for LinearRegressionL1L2TV/ LinearRegressionL1L2GL I get Error: "that shapes of A don't match beta"
How should I use it? or if there algorithms which supposed to be used with meshes?

For example use Travis.

Actually, by default, FISTA is using the old stopping criteria. We (with Fouad) suggest to change that: If the function inherit from DualFunction (implement the gap method) FISTA should use is because the duality gap provides explicit upper bound of the error.

It's not clear whether x ** 0.5 is faster than math.sqrt(x):
Which is faster in Python: x**.5 or math.sqrt(x)?

I find x ** 0.5 less readable so I suggest getting back to math.sqrt(x).

The following two parameters of LinearRegressionL1L2TV:
A : Numpy or (usually) scipy.sparse array. The linear operator for the smoothed total variation Nesterov function. A must be given.
mu : Non-negative float. The regularization constant for the smoothing.

I don't know why it works. How does it work?
Could anyone show some papers about the two parameters?

Thanks!

To avoid the re-computation of lambda max, while searching a grid of parameters we would like to provide pre-computed values of lambda max for both the data (X) and the linear operator (A).

First solution: provide lambda_max in the constructors of Estimators and Functions. The drawback is the large modifications of the API.

Second solution: provide helpers function, ie
precompute(X) and precompute(A) that computes the lambda max and eventually other quantities, and store it on the numpy objects X and A. This way, the lambda max will be propagated inside the loss and tv functions. There is no modification of the API, just a simple addition of code in the constructors of losses (eg. RidgeLogisticRegression, etc.) and penalties (eg. TotalVariation) to check the existence of such lambda max and to store it. Potentially, the precompute() function could add a decorator such lambda_max() to X and A numpy object instead of a simple additional attribute.

@duchesnay Is it possible to supply the Lasso regularization parameter as a vector instead of a scalar?
Thanks