Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP

Loading…

Minor improvement of pinv2 and pinvh #393

Merged
merged 2 commits into from

2 participants

@akhmerov

Avoid matrix operations on array elements that are multiplied by 0
in the process of calculation.

@akhmerov akhmerov Improve pinv2 and pinvh.
Avoid matrix operations on array elements that are multiplied by 0
in the process of calculation.
7e4298d
scipy/linalg/basic.py
((9 lines not shown))
- B = np.transpose(np.conjugate(np.dot(u * psigma_diag, vh)))
+ B = np.transpose(np.conjugate(np.dot(u[:, : rank] * psigma_diag, vh[: rank])))
@jakevdp Collaborator
jakevdp added a note

PEP8: line is too long

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@jakevdp
Collaborator

Looks great - I was thinking about doing something along these lines after I added pinvh, but never got to it.
Tests all pass for me; other than my one pep8 comment above, I think this is ready to merge.

@akhmerov

Thanks, and my bad about PEP8.

@jakevdp
Collaborator

Great - I'll let this sit for a few days in case others have comments, but I'm +1 for merge. Thanks for the contribution!

@jakevdp
Collaborator

Merging, thanks!

@jakevdp jakevdp merged commit a0c0a2a into from
@akhmerov

@jakevdp Is there anything that needs to be done for this to propagate to numpy.linalg?

@jakevdp
Collaborator

numpy.linalg is a separate package - you might consider copying this code and submitting a numpy PR. Note that numpy.linalg only contains one pinv function, which is implemented via SVD similar to scipy's pinv2.

@jakevdp
Collaborator

I should add that I'm not sure if the numpy devs would want functionality added to np.linalg (e.g. alternative forms of pinv may not be a good fit, but you could ask) but and improved implementation of the current function will certainly be appreciated there.

@ClemensFMN ClemensFMN referenced this pull request from a commit
Commit has since been removed from the repository and is no longer available.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Commits on Jan 1, 2013
  1. @akhmerov

    Improve pinv2 and pinvh.

    akhmerov authored
    Avoid matrix operations on array elements that are multiplied by 0
    in the process of calculation.
  2. @akhmerov

    Fix PEP8

    akhmerov authored
This page is out of date. Refresh to see the latest.
Showing with 13 additions and 12 deletions.
  1. +13 −12 scipy/linalg/basic.py
View
25 scipy/linalg/basic.py
@@ -659,24 +659,25 @@ def pinv2(a, cond=None, rcond=None, return_rank=False, check_finite=True):
factor = {'f': 1E3, 'd': 1E6}
cond = factor[t] * np.finfo(t).eps
- above_cutoff = (s > cond * np.max(s))
- psigma_diag = np.zeros_like(s)
- psigma_diag[above_cutoff] = 1.0 / s[above_cutoff]
+ rank = np.sum(s > cond * np.max(s))
+ psigma_diag = 1.0 / s[: rank]
- B = np.transpose(np.conjugate(np.dot(u * psigma_diag, vh)))
+ B = np.transpose(np.conjugate(np.dot(u[:, : rank] *
+ psigma_diag, vh[: rank])))
if return_rank:
- return B, np.sum(above_cutoff)
+ return B, rank
else:
return B
def pinvh(a, cond=None, rcond=None, lower=True, return_rank=False,
check_finite=True):
- """Compute the (Moore-Penrose) pseudo-inverse of a hermetian matrix.
+ """Compute the (Moore-Penrose) pseudo-inverse of a Hermitian matrix.
- Calculate a generalized inverse of a symmetric matrix using its
- eigenvalue decomposition and including all 'large' eigenvalues.
+ Calculate a generalized inverse of a Hermitian or real symmetric matrix
+ using its eigenvalue decomposition and including all eigenvalues with
+ 'large' absolute value.
Parameters
----------
@@ -735,14 +736,14 @@ def pinvh(a, cond=None, rcond=None, lower=True, return_rank=False,
factor = {'f': 1E3, 'd': 1E6}
cond = factor[t] * np.finfo(t).eps
- # For hermitian matrices, singular values equal abs(eigenvalues)
+ # For Hermitian matrices, singular values equal abs(eigenvalues)
above_cutoff = (abs(s) > cond * np.max(abs(s)))
- psigma_diag = np.zeros_like(s)
- psigma_diag[above_cutoff] = 1.0 / s[above_cutoff]
+ psigma_diag = 1.0 / s[above_cutoff]
+ u = u[:, above_cutoff]
B = np.dot(u * psigma_diag, np.conjugate(u).T)
if return_rank:
- return B, np.sum(above_cutoff)
+ return B, len(psigma_diag)
else:
return B
Something went wrong with that request. Please try again.