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

ENH: orthogonal procrustes solver #3809

Merged
merged 8 commits into from Aug 18, 2014

Conversation

argriffing
Copy link
Contributor

This introduces a new feature rather than bug fixing or maintenance, so it's now open for the bike-shedding!

@argriffing
Copy link
Contributor Author

See the discussion at #3786.

@argriffing argriffing mentioned this pull request Jul 18, 2014
A = np.asarray_chkfinite(A)
B = np.asarray_chkfinite(B)
else:
A = np.asarray(A)
Copy link
Contributor

Choose a reason for hiding this comment

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

A = np.ascontiguousarray(A)

@sturlamolden
Copy link
Contributor

Hm, it seems np.dot always returns C order arrays. This will be better:

A = np.asanyarray(A)
B = np.asanyarray(B)
...
u, w, vt = svd(np.dot(B.T,A).T)

@sturlamolden
Copy link
Contributor

or
u, w, vt = svd(B.T.dot(A).T)

@argriffing
Copy link
Contributor Author

Is it cool to use asanyarray? I don't know how picky I should be about supporting {ndarray, np.matrix, sparse} correctly.

@sturlamolden
Copy link
Contributor

The input to dot should be an ndarray, but it does not matter if it C or Fortran contiguous. But then the input to svd should be Fortran contiguous. Which means that those transpositions are actually preventing a transposition and a temporary copy inside the f2py layer... It's an f2py nastyness that transpose prevents transpose.

@argriffing
Copy link
Contributor Author

Just to be clear, the concern about array order and contiguity in the procrustes implementation is for efficiency not correctness, right?

@sturlamolden
Copy link
Contributor

It is mostly for saving memory

@coveralls
Copy link

Coverage Status

Changes Unknown when pulling bcc46d8 on argriffing:linalg-orthogonal-procrustes into * on scipy:master*.

@sturlamolden
Copy link
Contributor

But if you're in doubt, correctness always prewails.

@argriffing
Copy link
Contributor Author

By the way, there is supposedly a better way to compute the svd of a matrix product -- Gene Golub et al. "Computing the SVD of a general matrix product/quotient."

@coveralls
Copy link

Coverage Status

Changes Unknown when pulling bcc46d8 on argriffing:linalg-orthogonal-procrustes into * on scipy:master*.

@coveralls
Copy link

Coverage Status

Changes Unknown when pulling 94cbac0 on argriffing:linalg-orthogonal-procrustes into * on scipy:master*.

@argriffing
Copy link
Contributor Author

Added more tests and tried to improve memory efficiency according to the suggestions.

@coveralls
Copy link

Coverage Status

Changes Unknown when pulling 94cbac0 on argriffing:linalg-orthogonal-procrustes into * on scipy:master*.

@argriffing
Copy link
Contributor Author

To be more useful for applications, maybe the sum of singular values of the intermediate matrix should be returned?

@coveralls
Copy link

Coverage Status

Changes Unknown when pulling 6602cf8 on argriffing:linalg-orthogonal-procrustes into * on scipy:master*.

@pv pv added the PR label Jul 21, 2014
@pv pv removed the PR label Aug 13, 2014
Whether to check that the input matrices contain only finite numbers.
Disabling may give a performance gain, but may result in problems
(crashes, non-termination) if the inputs do contain infinities or NaNs.

Copy link
Contributor

Choose a reason for hiding this comment

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

Should there be a Raises section here?

@argriffing
Copy link
Contributor Author

Should there be a Raises section here?

I added this section explaining under what conditions a ValueError will be raised.

@coveralls
Copy link

Coverage Status

Changes Unknown when pulling 671ff41 on argriffing:linalg-orthogonal-procrustes into * on scipy:master*.

@ElDeveloper
Copy link
Contributor

@argriffing, thanks!

👍

u, w, vt = svd(B.T.dot(A).T)
R = u.dot(vt)
# Always return R, and maybe return a scaling factor.
if compute_scale:
Copy link
Member

Choose a reason for hiding this comment

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

I'm not an immediate user of this feature, so take this with a spoonful of salt: are you sure you really want to have the number of returns to depend on the input? In my experience it often makes it harder to reason about the results and leads to a proliferation of if compute_scale clauses at call sites. In this case always returning the scale does not seem to cost much, what is the benefit of having a variable number of returns?

@argriffing
Copy link
Contributor Author

are you sure you really want to have the number of returns to depend on the input?

This is now removed.

@coveralls
Copy link

Coverage Status

Changes Unknown when pulling 1195439 on argriffing:linalg-orthogonal-procrustes into * on scipy:master*.

@ev-br
Copy link
Member

ev-br commented Aug 18, 2014

LGTM.
It's been open for a month, which I guess qualifies for a long enough time for review. Merging, thanks @argriffing

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement A new feature or improvement scipy.linalg
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

8 participants