neurospin / pylearn-parsimony

Structured Machine Learning in Python
BSD 3-Clause "New" or "Revised" License
45 stars 16 forks source link

'v0' must have shape {shape} #35

Closed desanou closed 1 year ago

desanou commented 2 years ago

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.

DimitriPapadopoulos commented 2 years ago

Could be a SciPy issue. This is not the SciPy site.

Also, you should post an example that reproduces the problem, if at all possible, and at least the complete error messages.

DimitriPapadopoulos commented 2 years ago

The 'v0' must have shape {shape} error message originates in SciPy:

scipy/sparse/linalg/_eigen/_svds.py

        shape = (A.shape[0],) if solver == 'propack' else (min(A.shape),)
        if v0.shape != shape:
            message = "`v0` must have shape {shape}."
            raise ValueError(message)
DimitriPapadopoulos commented 2 years ago

Note that the https://github.com/scipy/scipy/pull/15566 ticket I have opened against SciPy does not address this issue.

desanou commented 2 years ago

Okay, thanks for your feedback. I have a reproducible example but it is not minimal. I'll look into the different libraries a bit more and get back to you if the problem is parsimony related.

I'll close the issue.

desanou commented 2 years ago

@DimitriPapadopoulos I reopen the issue with a minimal example that can produce it. I used the same examples provided by the parsimony-pylearn quick start section.

import numpy as np
np.random.seed(42)
shape = (1, 4, 4)  # Three-dimension matrix
num_samples = 10  # The number of samples
num_ft = shape[0] * shape[1] * shape[2] # The number of features per sample

Randomly generate X

X = np.random.rand(num_samples, num_ft)
beta = np.random.rand(num_ft, 1) # Define beta
# Add noise to y
y = np.dot(X, beta) + 0.001 * np.random.rand(num_samples, 1)
X_train = X[0:6, :]
y_train = y[0:6]
X_test = X[6:10, :]
y_test = y[6:10]
import parsimony.estimators as estimators
import parsimony.algorithms as algorithms
import parsimony.functions.nesterov.tv as tv
l = 0.0  # l1 lasso coefficient
k = 0.0  # l2 ridge regression coefficient
g = 1.0  # tv coefficient
A = tv.linear_operator_from_shape(shape)  # Memory allocation for TV
olstv = estimators.LinearRegressionL1L2TV(l, k, g, A, mu=0.0001,
                                         algorithm=algorithms.proximal.FISTA(),
                                         algorithm_params=dict(max_iter=1000))

res = olstv.fit(X_train, y_train)
print("Estimated beta error = ", np.linalg.norm(olstv.beta - beta))
print("Prediction error = ", np.linalg.norm(olstv.predict(X_test) - y_test))

The code has been compiled with pylearn=parsimony 0.3.1 python 3.8 and scipy last version 1.9.0. I got the error: ValueError: v0 must have shape (16,).

A traceback is given below.

---------------------------------------------------------------------------
**ValueError**                                Traceback (most recent call last)
Input In [3], in <cell line: 1>()
----> 1 res = olstv.fit(X_train, y_train)

File ** .\parsimony\estimators.py:888**, in LinearRegressionL1L2TV.fit(self, X, y, beta)
--> 888 self.beta = self.algorithm.run(function, beta)

File **.\parsimony\algorithms\bases.py:101**, in force_reset.<locals>.wrapper(self, function, *args, **kwargs)
--> 101 return f(self, function, *args, **kwargs)

File **.\parsimony\algorithms\bases.py:81**, in check_compatibility.<locals>.wrapper(self, function, *args, **kwargs)
---> 81     return f(self, function, *args, **kwargs)

File **.\parsimony\algorithms\proximal.py:401**, in FISTA.run(self, function, beta)
--> 401 step = function.step(z)

File **.\parsimony\functions\combinedfunctions.py:1044**, in LinearRegressionL1L2TV.step(self, x, **kwargs)
-> 1044     return 1.0 / self.L()

File **.\parsimony\functions\combinedfunctions.py:769**, in LinearRegressionL1L2TV.L(self)
--> 769          + self.tv.L()

File **.\parsimony\functions\nesterov\tv.py:156**, in TotalVariation.L(self)
--> 156 lmaxA = self.lambda_max()

File **.\parsimony\functions\nesterov\tv.py:189**, in TotalVariation.lambda_max(self)
--> 189 v = RankOneSparseSVD().run(A)  # , max_iter=max_iter)

File **.\parsimony\algorithms\nipals.py:343**, in RankOneSparseSVD.run(self, X, start_vector)
--> 343     [_, _, v] = sparse_linalg.svds(X, k=1, v0=v0,

File **.\scipy\sparse\linalg\_eigen\_svds.py:252**, in svds(A, k, ncv, tol, which, v0, maxiter, return_singular_vectors, solver, random_state, options)
--> 252 args = _iv(A, k, ncv, tol, which, v0, maxiter, return_singular_vectors,

File **.\scipy\sparse\linalg\_eigen\_svds.py:84**, in _iv(A, k, ncv, tol, which, v0, maxiter, return_singular, solver, random_state)
---> 84         raise ValueError(message)

**ValueError**: `v0` must have shape (16,).

Using the scipy version 1.7.1 and the relevant numpy version does not raised an error. The error might be due to update of scipy > 1.7.1.

DimitriPapadopoulos commented 2 years ago

The error message had been introduced in SciPy 1.8.0rc1 by https://github.com/scipy/scipy/pull/14433.

v0 is initialised here, it would be great if you could print the value of X.shape: https://github.com/neurospin/pylearn-parsimony/blob/969ea182573dddf8cd70e93846613713598e3afa/parsimony/algorithms/nipals.py#L151

desanou commented 2 years ago

X.shape (10, 16)

DimitriPapadopoulos commented 2 years ago

According to the documentation of scipy.sparse.linalg.svds, the default solver, which we use here as we do not pass a solver argument, is ‘arpack’:

solver{‘arpack’, ‘propack’, ‘lobpcg’}, optional The solver used. ‘arpack’, ‘lobpcg’, and ‘propack’ are supported. Default: ‘arpack’.

According to the documentation of svds(solver=’arpack’):

In the descriptions below, let M, N = A.shape.

v0ndarray, optional The starting vector for iteration: an (approximate) left singular vector if N > M and a right singular vector otherwise. Must be of length min(M, N). Default: random

Because in our case we have M = 10 and N = 16. I would expect the length of v0 to be 10, and hence its shape (10,). ⟹ Could you perhaps print v0.shape?

Then, I don't understand why the input validation function _iv() of svds() raises this error where v0 must have shape (16,).