Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

added LinearOperator support in safe_sparse_dot() #16463

Open
wants to merge 13 commits into
base: main
Choose a base branch
from

Conversation

PavelStishenko
Copy link

Reference Issues/PRs

Fixes #16461.

@PavelStishenko
Copy link
Author

Please, help! I can not see details of the codecov/patch check. It shows 500 server error.

@@ -148,7 +149,10 @@ def safe_sparse_dot(a, b, dense_output=False):
else:
ret = np.dot(a, b)
else:
ret = a @ b
if isinstance(b, LinearOperator) and not isinstance(a, LinearOperator):
ret = (b.T @ a.T).T
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you explain this line? a @ b does not work? Do we actually need to handle this case? I'm OK with the rest of changes, but if we can avoid this special case it would be better IMO.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you explain this line? a @ b does not work? Do we actually need to handle this case? I'm OK with the rest of changes, but if we can avoid this special case it would be better IMO.

Yes, a @ b doesn't work if b is a LinearOperator and a is not. Because the first argument of the @ operator defines actual function that do the job. LinearOperator can be multiplied by anything, but ndarray or sparse matrices know nothing about LinearOperator and do not know how to proceed.

import numpy as np
from scipy.sparse.linalg import aslinearoperator
A = np.random.random((4,4))
op = aslinearoperator(A)

print (op @ A) # works fine
print (A @ op) # ValueError: matmul: Input operand 1 does not have enough dimensions (has 0, gufunc core with signature (n?,k),(k,m?)->(n?,m?) requires 1)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And so this is necessary for randomized_svd?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And so this is necessary for randomized_svd?

Definitely. LinearOperator can incapsulate anything that behave like a matrix via matvec() implementation. Obvious examples are SVD decomposition, einsum() routine, any tensor network with two free legs.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not return b@a then, it should work since here b and a are 1d vectors?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not return b@a then, it should work since here b and a are 1d vectors?

They can be 2D.

@PavelStishenko
Copy link
Author

Any help with codecov/patch please. It still shows 500-th server error

Copy link
Member

@rth rth left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please add a unit test, testing the behaviour of safe_sparse_dot(array, LinearOperator) and safe_sparse_dot(LinearOperator, array) to sklearn/utils/tests/test_extmath.py.

That would also fix the coverage error.

Also please add an entry to the change log at doc/whats_new/v0.23.rst. Like the other entries there, please reference this pull request with :pr: and credit yourself (and other contributors if applicable) with :user:.

pvstishenko added 2 commits February 21, 2020 17:16
…_dot(LinearOperator, array) adjusted docstrings and added whats_new entry
Copy link
Member

@rth rth left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Currently the added test fails in some CI with,

AttributeError: 'MatrixLinearOperator' object has no attribute '_transpose'

Maybe there is a minimal scipy version required? If so, that is fine, but it should be indicated in the what's new, with a comment in the code and the test skipped using,

pytest.importorskip("docutils", minversion="???")

@@ -341,6 +341,10 @@ Changelog
pandas sparse DataFrame.
:pr:`16021` by :user:`Rushabh Vasani <rushabh-v>`.

- |Feature| :func:`utils.extmath.randomized_svd` now accepts
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

and utils.extmath.safe_dot_product

@@ -341,7 +341,8 @@ Changelog
pandas sparse DataFrame.
:pr:`16021` by :user:`Rushabh Vasani <rushabh-v>`.

- |Feature| :func:`utils.extmath.randomized_svd` now accepts
- |Feature| :func:`utils.extmath.randomized_svd` and
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Tests are still failing. Also please add a test for randomized_svd on a LinearOperator.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fixed. The test added

Copy link
Member

@rth rth left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Minor comment otherwise LGTM. Thanks @Mazay0 !

sklearn/utils/tests/test_extmath.py Outdated Show resolved Hide resolved
… with naming convention and added comment about motivation behind this class
@PavelStishenko
Copy link
Author

@glemaitre review please.

@rth rth requested a review from NicolasHug April 8, 2020 11:36
Copy link
Member

@NicolasHug NicolasHug left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the PR @Mazay0 . Made a few comments.

@rth I can't say I'm sold on the need to support LinearOperator. The changes seem minimal here, but who knows how they might complicate things in the future? I'm not saying we shouldn't do this, but I think we need to be super careful when introducing support for a new data structure. WDYT?

@@ -148,7 +149,10 @@ def safe_sparse_dot(a, b, dense_output=False):
else:
ret = np.dot(a, b)
else:
ret = a @ b
if isinstance(b, LinearOperator) and not isinstance(a, LinearOperator):
ret = (b.T @ a.T).T
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not return b@a then, it should work since here b and a are 1d vectors?

sklearn/utils/tests/test_extmath.py Show resolved Hide resolved
@rth
Copy link
Member

rth commented Apr 8, 2020

I can't say I'm sold on the need to support LinearOperator. The changes seem minimal here, but who knows how they might complicate things in the future? I'm not saying we shouldn't do this, but I think we need to be super careful when introducing support for a new data structure. WDYT?

@NicolasHug I generally share your feeling. However the changes are indeed minimal and LinearOperator is part of of scipy.sparse so it's not like we are adding support for a new library. Of course, I'm not saying that we would support it everywhere, not will we provide much official support. But if it makes life easier for some contributors, they are willing make a PR to fix it, and changes are minimal, then why not.

@PavelStishenko
Copy link
Author

Thanks for the PR @Mazay0 . Made a few comments.

@rth I can't say I'm sold on the need to support LinearOperator. The changes seem minimal here, but who knows how they might complicate things in the future? I'm not saying we shouldn't do this, but I think we need to be super careful when introducing support for a new data structure. WDYT?

LinearOperator is supposed to mimic matrix behavior. Progress in its development goes in this direction, so in future in should be even simpler to support it. It is especially true given that sklearn supports scipy's sparse matrices and anyway shouldn't rely on internal data representation.
Key feature of Randomized-SVD is that it is a matrix-free method and support of LinearOperator is the most natural way to employ this feature.

@rth
Copy link
Member

rth commented Apr 19, 2020

@NicolasHug Would be good to make a decision on this one (in one direction or another). Has your position changed since last review?

@NicolasHug
Copy link
Member

Has your position changed since last review?

Not really honestly. Maybe other @scikit-learn/core-devs should chime in

Copy link
Member

@adrinjalali adrinjalali left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Where would this apply? only when people directly pass a LinearOperator to these functions? what are the consequences of the changed returned value for the user in places where these functions are used internally?

@@ -131,7 +132,7 @@ def safe_sparse_dot(a, b, dense_output=False):
dot_product : array or sparse matrix
sparse if ``a`` and ``b`` are sparse and ``dense_output=False``.
"""
if a.ndim > 2 or b.ndim > 2:
if len(a.shape) > 2 or len(b.shape) > 2:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I vaguely remember this potentially being expensive compared to a.ndim. Is that not a concern here?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Probably it is more expensive, but the difference is negligible compared with subsequent operations.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I vaguely remember this potentially being expensive compared to a.ndim. Is that not a concern here?

scipy folks provisionally agreed to accept PR with added ndim attribute to the LinearOperator: scipy/scipy#11908 (comment)
So this row can be reverted to the original state.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

SciPy folks merged PR with ndim attribute for LinearOperator: scipy/scipy#11915 So in future this row can be reverted to ndim again.

@@ -131,7 +132,7 @@ def safe_sparse_dot(a, b, dense_output=False):
dot_product : array or sparse matrix
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This PR also changes the output type.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes indeed. I will fix the docstring and will add an appropriate test, if you at least in principle agree that this PR can be accepted.

@rth
Copy link
Member

rth commented Apr 20, 2020

Where would this apply? only when people directly pass a LinearOperator to these functions?

Yes, only if you pass LinearOperator as input.

what are the consequences of the changed returned value for the user in places where these functions are used internally?

Right now it would have no consequence internally since it's safe_dot_product is never used with LinearOperator. We do use LinearOperator is other parts of the code base. However it would allow users to use LinearOperator as input to randomized_svd (cf parent issue) in their own code. The proposed changes seem to me like a small price for that. Of course we are not promising any official support of it.

Several scipy function, particularly for matrix decomposition (svds etc) support LinearOperator as input alongside sparse matrices https://docs.scipy.org/doc/scipy/reference/sparse.linalg.html#matrix-factorizations

@PavelStishenko
Copy link
Author

Where would this apply? only when people directly pass a LinearOperator to these functions? what are the consequences of the changed returned value for the user in places where these functions are used internally?

Only when people pass some objects derived from LinearOperator. This is nice if you want to encapsulate internal data structure in some custom way, other than scipy.sparse classes. Returned value will not change for cases where numpy or scipy sparse arrays are used. LinearOperators could not be used before, so no code will be broken by this PR.

@adrinjalali
Copy link
Member

Thanks for following up the ndim issue on scipy @Mazay0 .

I checked the LinearOperator doc, and it seems it's kind of a separate thing compared to __array_function__. I was hoping that when we move to supporting either __array_function__ or __array_module__, then LinearOperator would automatically be supported, but it seems that's not the case.

I guess the other question is, why is LinearOperator not implementing __array_function__?

@jnothman
Copy link
Member

We should probably have discussed this at the monthly meeting. Labelling as needs decision.

@jnothman jnothman added the Needs Decision Requires decision label Apr 27, 2020
Copy link
Member

@jnothman jnothman left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unless we explicitly document some estimators or equivalent functions as accepting a LinearOperator, I find this support in a utility a bit obscure.

Base automatically changed from master to main January 22, 2021 10:52
@glemaitre glemaitre removed their request for review April 10, 2021 16:58
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Add support of LinearOperator in randomized_svd
5 participants