-
-
Notifications
You must be signed in to change notification settings - Fork 10.8k
BUG: Fix pinv for stacked matrices #8827
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
Conversation
numpy/linalg/tests/test_linalg.py
Outdated
@@ -423,25 +423,6 @@ def test_generalized_empty_herm_cases(self): | |||
exclude={'none'}) | |||
|
|||
|
|||
def dot_generalized(a, b): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ahead of its time...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Except not, as it turns out. This doesn't broadcast, but instead tries to work out which are vectors and which are matrices. Guess we're keeping it, for now
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Funny, I had exactly the same thought -- couldn't you use matmul
instead -- but then realised it was actually different (in the actual expression, I might have used einsum
instead, but it is not really important at all).
@@ -226,6 +227,22 @@ def _assertNoEmpty2d(*arrays): | |||
if _isEmpty2d(a): | |||
raise LinAlgError("Arrays cannot be empty") | |||
|
|||
def transpose(a): | |||
""" | |||
Transpose each matrix in a stack of matrices. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
TODO: expose this publicly? (#7495)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Well, there is at least one implementation that could (eventually) be abandoned if this is exposed: https://github.com/astropy/astropy/blob/master/astropy/coordinates/matrix_utilities.py#L32
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Any reason for calling the method and not the free function?
Also, i'm increasingly considering using astropy for completely non-astrophysics related things in future - your units sounds more actively-developed than Pint
, and rotation matrices are something I end up writing too often.
Regarding the function before the one you link - #8719
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
The method is the real thing; the function just turns the input into an array if needed and then calls the method. https://github.com/numpy/numpy/blob/master/numpy/core/fromnumeric.py#L458
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yep, I know :). Guess people are unlikely to pass around lists of lists
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
p.s. In the astropy case this is only used internally in cases where I know the inputs are arrays; if we exposed this publically in numpy, I think using the function makes more sense.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We need something like this, but I think it should be an array method as well. Maybe something like a.M(atrix)T(ranspose
. When the *.T
method was proposed there was discussion of what it should do, transposing matrices in a stack was one of ideas. I could also see a vector transpose, a.VT
which just appends 1 to the dimensions. These ideas need discussion on the list.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Should I add a leading underscore to the name for now?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@eric-wieser That would help distinguish it from np.transpose
and avoid some confusion. Long term, it probably needs a different name.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
At any rate, this isn't going to appear in the docs at the moment, so isn't really public API yet.
I'll send something around the mailing list requesting a name at some point. I think we can start off with the conservative assert ndim >= 2
semantics initially, and put off the decision in #7495 and #9530
numpy/linalg/linalg.py
Outdated
@@ -2008,13 +2020,13 @@ def lstsq(a, b, rcond=-1): | |||
resids = array([sum((ravel(bstar)[n:])**2)], | |||
dtype=result_real_t) | |||
else: | |||
x = array(transpose(bstar)[:n,:], dtype=result_t, copy=True) | |||
x = array(nd_transpose(bstar)[:n,:], dtype=result_t, copy=True) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If you're changing it anyway, maybe just do bstar.T
? So much clearer...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Of course, since bstar
is 2-D, you can also just leave it as it was...
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I mean, I plan to throw this all out in #8720 anyway...
The only reason I changed this line is because I wanted to use the name transpose
for something else.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Done
8e02ce0
to
ad36b66
Compare
doc/release/1.13.0-notes.rst
Outdated
@@ -159,6 +159,10 @@ being iterated over. | |||
For consistency with ``ndarray`` and ``broadcast``, ``d.ndim`` is a shorthand | |||
for ``len(d.shape)``. | |||
|
|||
``matrix_rank`` and ``pinv`` in ``np.linalg`` now works for stacked arrays | |||
-------------------------------------------------------------------------- |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
works -> work
numpy/linalg/linalg.py
Outdated
matrices was added in Numpy 1.13.0. | ||
rcond : (...) arr_like of float | ||
Cutoff for small singular values. | ||
Singular values smaller (in modulus) than |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
How about "less than" rather than smaller. The singular values are always non-negative reals.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Even for complex matrices.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Note that I haven't actually changed these comments, just adjusted their indentation.
The singular values are always non-negative reals.
Is this the case even under numerical error?
Either way, the description right now matches the implementation. If taking the modulus is not required, then we should remove it from the doc and implementation at the same time (and in another PR). Edit: I'm wrong
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Puzzling -- as @charris notes, singular values are by definition positive real, so reading something about a modulus would seem only confusing.
☔ The latest upstream changes (presumably #8886) made this pull request unmergeable. Please resolve the merge conflicts. |
☔ The latest upstream changes (presumably #9050) made this pull request unmergeable. Please resolve the merge conflicts. |
Updated to refer to 1.14 not 1.13 in release notes |
8b29605
to
bffaf0d
Compare
@eric-wieser Want to rebase this? I'm trying to get all the linalg fixes/enhancements in. |
Yes, but I can't until I regain access to a git installation in two weeks - so maybe someone else should. |
Conflict is just in the imports |
7b6f25d
to
219d2a5
Compare
@charris: Ok, rebased. Ready to go in? |
I'm putting this in as is. Thanks Eric. If you want to change the name of |
Fixes #8826
Duplicate of #4678