Skip to content
This repository

Pinv enhancements #289

Closed
wants to merge 7 commits into from

5 participants

Jake Vanderplas Pauli Virtanen Charles Harris Josef Perktold Alexandre Gramfort
Jake Vanderplas
Collaborator

This is a PR based on some of @vene's work in scikit-learn: see scikit-learn/scikit-learn#1015

There are two enhancements:

  • scipy.linalg.pinv2 is sped up by only computing necessary singular values, and not allocating the entire psigma matrix

  • scipy.linalg.pinvh added: this uses eigh to significantly speed up the pseudo-inverse computation in the symmetric/hermitian case

Here are some simple benchmarks:

old pinv2 vs new pinv2

old version:

In [1]: import numpy as np
In [2]: from scipy.linalg import pinv2
In [3]: np.random.seed(0)
In [4]: X = np.random.random((500, 100))
In [5]: %timeit pinv2(X)
10 loops, best of 3: 172 ms per loop

new version:

...
In [5]: %timeit pinv2(X)
10 loops, best of 3: 38.1 ms per loop

pinvh vs new pinv2:

In [1]: import numpy as np
In [2]: from scipy.linalg import pinv2, pinvh
In [3]: np.random.seed(0)
In [4]: X = np.random.random((500, 400))
In [5]: X = np.dot(X, X.T) # make symmetric positive semi-definite
In [6]: %timeit pinv2(X)
1 loops, best of 3: 736 ms per loop
In [7]: %timeit pinvh(X)
1 loops, best of 3: 320 ms per loop
Charles Harris charris commented on the diff August 16, 2012
scipy/linalg/basic.py
((8 lines not shown))
587 587
     if rcond is not None:
588 588
         cond = rcond
589 589
     if cond in [None,-1]:
590  
-        eps = np.finfo(np.float).eps
591  
-        feps = np.finfo(np.single).eps
592  
-        _array_precision = {'f': 0, 'd': 1, 'F': 0, 'D': 1}
593  
-        cond = {0: feps*1e3, 1: eps*1e6}[_array_precision[t]]
594  
-    m, n = a.shape
595  
-    cutoff = cond*np.maximum.reduce(s)
596  
-    psigma = np.zeros((m, n), t)
597  
-    for i in range(len(s)):
598  
-        if s[i] > cutoff:
599  
-            psigma[i,i] = 1.0/np.conjugate(s[i])
600  
-    #XXX: use lapack/blas routines for dot
601  
-    return np.transpose(np.conjugate(np.dot(np.dot(u,psigma),vh)))
  590
+        t = u.dtype.char.lower()
  591
+        factor = {'f': 1E3, 'd': 1E6}
2
Charles Harris Collaborator
charris added a note August 16, 2012

I wonder if the updated version in np.matrix_rank would be a better choice for the default condition? See the discussion at http://thread.gmane.org/gmane.comp.python.numeric.general/50396/focus=50912

Jake Vanderplas Collaborator
jakevdp added a note August 16, 2012

Interesting discussion: here I just duplicated the previous behavior. Any change would be a potential compatibility issue: before making that decision we should probably discuss it on the mailing list.

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

I just had another thought: we should add a lower keyword to pinvh which is passed to eigh.

scipy/linalg/basic.py
((53 lines not shown))
  620
+
  621
+    Returns
  622
+    -------
  623
+    B : array, shape (N, N)
  624
+
  625
+    Raises
  626
+    ------
  627
+    LinAlgError
  628
+        If eigenvalue does not converge
  629
+
  630
+    Examples
  631
+    --------
  632
+    >>> from numpy import *
  633
+    >>> a = random.randn(9, 6)
  634
+    >>> a = np.dot(a, a.T)
  635
+    >>> B = pinv_symmetric(a)
2
Alexandre Gramfort
agramfort added a note August 17, 2012

docstring is not up to date right?

Jake Vanderplas Collaborator
jakevdp added a note August 17, 2012

right - good catch

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

@pv - I think this is ready for merge. @charris's comment about matrix rank would be interesting to explore, but I think it's an enhancement that would affect more than just this PR.

scipy/linalg/basic.py
((70 lines not shown))
  637
+    True
  638
+    >>> allclose(B, dot(B, dot(a, B)))
  639
+    True
  640
+
  641
+    """
  642
+    a = np.asarray_chkfinite(a)
  643
+    s, u = decomp.eigh(a, lower=lower)
  644
+    
  645
+    if rcond is not None:
  646
+        cond = rcond
  647
+    if cond in [None, -1]:
  648
+        t = u.dtype.char.lower()
  649
+        factor = {'f': 1E3, 'd': 1E6}
  650
+        cond = factor[t] * np.finfo(t).eps
  651
+    
  652
+    # unlike svd case, eigh can lead to negative eigenvalues
1
Pauli Virtanen Owner
pv added a note August 25, 2012

More illuminating comment: "For hermitian matrices, singular values equal abs(eigenvalues)"

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Pauli Virtanen
Owner
pv commented August 25, 2012

@jakevdp: looks good to me. Changing the default cutoff is maybe out of scope for this PR. It seems we don't have tests for single-precision data types?

Charles Harris
Collaborator

It might be worth bringing over some of Matthew's documentation for np.matrix_rank, I think it is quite good. But I agree that changing the current default choices for the cutoff is independent of this PR.

Josef Perktold
Collaborator

a suggestion: can we add a keyword argument to optionally return the rank that was used, or eigen or singular values?

looking briefly through the scikit learn PR
"Yes... though I don't think it's possible to do this quickly using the current pinv or pinv2. We should add a flag that optionally returns the number of singular values used in the computation."

In statsmodels, we have currently the problem that we are doing two svd, one for the pinv and one to get the rank. The only option I saw so far was to copy pinv directly into statsmodels with the additional return.

Jake Vanderplas
Collaborator

@josef-pkt: good suggestion on the flag to return the rank. It's a few lines of code that add a lot of flexibility to the use of the routines.

I think I've addressed all the comments (short of the re-working of singular value cutoff selection). This should be ready to merge.

Pauli Virtanen pv referenced this pull request from a commit October 21, 2012
Pauli Virtanen ENH: linalg: Merge PR #289 Pinv enhancements
This is a PR based on some of @vene's work in scikit-learn: see
scikit-learn/scikit-learn#1015

There are three enhancements:

- scipy.linalg.pinv2 is sped up by only computing necessary singular
  values, and not allocating the entire psigma matrix

- scipy.linalg.pinvh added: this uses eigh to significantly speed up the
  pseudo-inverse computation in the symmetric/hermitian case

- add return_rank keyword to obtain the effective rank from all pinv*
  routines.
5adeb12
Pauli Virtanen
Owner
pv commented October 21, 2012

Thanks, merged.

Pauli Virtanen pv closed this October 21, 2012
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
This page is out of date. Refresh to see the latest.
1  scipy/linalg/__init__.py
@@ -31,6 +31,7 @@
31 31
    lstsq - Solve a linear least-squares problem
32 32
    pinv - Pseudo-inverse (Moore-Penrose) using lstsq
33 33
    pinv2 - Pseudo-inverse using svd
  34
+   pinvh - Pseudo-inverse of hermitian matrix
34 35
    kron - Kronecker product of two arrays
35 36
    tril - Construct a lower-triangular matrix from a given matrix
36 37
    triu - Construct an upper-triangular matrix from a given matrix
123  scipy/linalg/basic.py
@@ -2,9 +2,10 @@
2 2
 # Author: Pearu Peterson, March 2002
3 3
 #
4 4
 # w/ additions by Travis Oliphant, March 2002
  5
+#              and Jake Vanderplas, August 2012
5 6
 
6 7
 __all__ = ['solve', 'solve_triangular', 'solveh_banded', 'solve_banded',
7  
-            'inv', 'det', 'lstsq', 'pinv', 'pinv2']
  8
+            'inv', 'det', 'lstsq', 'pinv', 'pinv2', 'pinvh']
8 9
 
9 10
 import numpy as np
10 11
 
@@ -13,7 +14,7 @@
13 14
 from misc import LinAlgError, _datacopied
14 15
 from scipy.linalg import calc_lwork
15 16
 from decomp_schur import schur
16  
-import decomp_svd
  17
+import decomp, decomp_svd
17 18
 
18 19
 
19 20
 # Linear equations
@@ -500,7 +501,7 @@ def lstsq(a, b, cond=None, overwrite_a=False, overwrite_b=False):
500 501
     return x, resids, rank, s
501 502
 
502 503
 
503  
-def pinv(a, cond=None, rcond=None):
  504
+def pinv(a, cond=None, rcond=None, return_rank=False):
504 505
     """
505 506
     Compute the (Moore-Penrose) pseudo-inverse of a matrix.
506 507
 
@@ -515,11 +516,15 @@ def pinv(a, cond=None, rcond=None):
515 516
         Cutoff for 'small' singular values in the least-squares solver.
516 517
         Singular values smaller than ``rcond * largest_singular_value``
517 518
         are considered zero.
  519
+    return_rank : bool, optional
  520
+        if True, return the effective rank of the matrix
518 521
 
519 522
     Returns
520 523
     -------
521 524
     B : array, shape (N, M)
522 525
         The pseudo-inverse of matrix `a`.
  526
+    rank : int
  527
+        The effective rank of the matrix.  Returned if return_rank == True
523 528
 
524 529
     Raises
525 530
     ------
@@ -540,10 +545,16 @@ def pinv(a, cond=None, rcond=None):
540 545
     b = np.identity(a.shape[0], dtype=a.dtype)
541 546
     if rcond is not None:
542 547
         cond = rcond
543  
-    return lstsq(a, b, cond=cond)[0]
  548
+
  549
+    x, resids, rank, s = lstsq(a, b, cond=cond)
  550
+
  551
+    if return_rank:
  552
+        return x, rank
  553
+    else:
  554
+        return x
544 555
 
545 556
 
546  
-def pinv2(a, cond=None, rcond=None):
  557
+def pinv2(a, cond=None, rcond=None, return_rank=False):
547 558
     """
548 559
     Compute the (Moore-Penrose) pseudo-inverse of a matrix.
549 560
 
@@ -560,11 +571,15 @@ def pinv2(a, cond=None, rcond=None):
560 571
         Singular values smaller than ``rcond*largest_singular_value``
561 572
         are considered zero.
562 573
         If None or -1, suitable machine precision is used.
  574
+    return_rank : bool, optional
  575
+        if True, return the effective rank of the matrix
563 576
 
564 577
     Returns
565 578
     -------
566 579
     B : array, shape (N, M)
567 580
         The pseudo-inverse of matrix `a`.
  581
+    rank : int
  582
+        The effective rank of the matrix.  Returned if return_rank == True
568 583
 
569 584
     Raises
570 585
     ------
@@ -582,21 +597,91 @@ def pinv2(a, cond=None, rcond=None):
582 597
 
583 598
     """
584 599
     a = np.asarray_chkfinite(a)
585  
-    u, s, vh = decomp_svd.svd(a)
586  
-    t = u.dtype.char
  600
+    u, s, vh = decomp_svd.svd(a, full_matrices=False)
  601
+    
587 602
     if rcond is not None:
588 603
         cond = rcond
589 604
     if cond in [None,-1]:
590  
-        eps = np.finfo(np.float).eps
591  
-        feps = np.finfo(np.single).eps
592  
-        _array_precision = {'f': 0, 'd': 1, 'F': 0, 'D': 1}
593  
-        cond = {0: feps*1e3, 1: eps*1e6}[_array_precision[t]]
594  
-    m, n = a.shape
595  
-    cutoff = cond*np.maximum.reduce(s)
596  
-    psigma = np.zeros((m, n), t)
597  
-    for i in range(len(s)):
598  
-        if s[i] > cutoff:
599  
-            psigma[i,i] = 1.0/np.conjugate(s[i])
600  
-    #XXX: use lapack/blas routines for dot
601  
-    return np.transpose(np.conjugate(np.dot(np.dot(u,psigma),vh)))
  605
+        t = u.dtype.char.lower()
  606
+        factor = {'f': 1E3, 'd': 1E6}
  607
+        cond = factor[t] * np.finfo(t).eps
  608
+    
  609
+    above_cutoff = (s > cond * np.max(s))
  610
+    psigma_diag = np.zeros_like(s)
  611
+    psigma_diag[above_cutoff] = 1.0 / s[above_cutoff]
  612
+ 
  613
+    B = np.transpose(np.conjugate(np.dot(u * psigma_diag, vh)))
  614
+
  615
+    if return_rank:
  616
+        return B, np.sum(above_cutoff)
  617
+    else:
  618
+        return B
  619
+
  620
+
  621
+def pinvh(a, cond=None, rcond=None, lower=True, return_rank=False):
  622
+    """Compute the (Moore-Penrose) pseudo-inverse of a hermetian matrix.
  623
+
  624
+    Calculate a generalized inverse of a symmetric matrix using its
  625
+    eigenvalue decomposition and including all 'large' eigenvalues.
  626
+
  627
+    Parameters
  628
+    ----------
  629
+    a : array, shape (N, N)
  630
+        Real symmetric or complex hermetian matrix to be pseudo-inverted
  631
+    cond, rcond : float or None
  632
+        Cutoff for 'small' eigenvalues.
  633
+        Singular values smaller than rcond * largest_eigenvalue are considered
  634
+        zero.
  635
+
  636
+        If None or -1, suitable machine precision is used.
  637
+    lower : boolean
  638
+        Whether the pertinent array data is taken from the lower or upper
  639
+        triangle of a. (Default: lower)
  640
+    return_rank : bool, optional
  641
+        if True, return the effective rank of the matrix
  642
+
  643
+    Returns
  644
+    -------
  645
+    B : array, shape (N, N)
  646
+        The pseudo-inverse of matrix `a`.
  647
+    rank : int
  648
+        The effective rank of the matrix.  Returned if return_rank == True
  649
+
  650
+    Raises
  651
+    ------
  652
+    LinAlgError
  653
+        If eigenvalue does not converge
  654
+
  655
+    Examples
  656
+    --------
  657
+    >>> from numpy import *
  658
+    >>> a = random.randn(9, 6)
  659
+    >>> a = np.dot(a, a.T)
  660
+    >>> B = pinvh(a)
  661
+    >>> allclose(a, dot(a, dot(B, a)))
  662
+    True
  663
+    >>> allclose(B, dot(B, dot(a, B)))
  664
+    True
602 665
 
  666
+    """
  667
+    a = np.asarray_chkfinite(a)
  668
+    s, u = decomp.eigh(a, lower=lower)
  669
+    
  670
+    if rcond is not None:
  671
+        cond = rcond
  672
+    if cond in [None, -1]:
  673
+        t = u.dtype.char.lower()
  674
+        factor = {'f': 1E3, 'd': 1E6}
  675
+        cond = factor[t] * np.finfo(t).eps
  676
+    
  677
+    # For hermitian matrices, singular values equal abs(eigenvalues)
  678
+    above_cutoff = (abs(s) > cond * np.max(abs(s)))
  679
+    psigma_diag = np.zeros_like(s)
  680
+    psigma_diag[above_cutoff] = 1.0 / s[above_cutoff]
  681
+
  682
+    B = np.dot(u * psigma_diag, np.conjugate(u).T)
  683
+
  684
+    if return_rank:
  685
+        return B, np.sum(above_cutoff)
  686
+    else:
  687
+        return B
56  scipy/linalg/tests/test_basic.py
@@ -28,7 +28,7 @@
28 28
     assert_equal, assert_almost_equal, assert_array_almost_equal, assert_, \
29 29
     assert_allclose
30 30
 
31  
-from scipy.linalg import solve, inv, det, lstsq, pinv, pinv2, norm,\
  31
+from scipy.linalg import solve, inv, det, lstsq, pinv, pinv2, pinvh, norm,\
32 32
         solve_banded, solveh_banded, solve_triangular
33 33
 
34 34
 from scipy.linalg._testutils import assert_no_overwrite
@@ -534,30 +534,66 @@ def test_random_complex_overdet(self):
534 534
 
535 535
 class TestPinv(TestCase):
536 536
 
537  
-    def test_simple(self):
538  
-        a=array([[1,2,3],[4,5,6.],[7,8,10]])
  537
+    def test_simple_real(self):
  538
+        a = array([[1, 2, 3], [4, 5, 6], [7, 8, 10]], dtype=float)
539 539
         a_pinv = pinv(a)
540  
-        assert_array_almost_equal(dot(a,a_pinv),[[1,0,0],[0,1,0],[0,0,1]])
  540
+        assert_array_almost_equal(dot(a,a_pinv), np.eye(3))
541 541
         a_pinv = pinv2(a)
542  
-        assert_array_almost_equal(dot(a,a_pinv),[[1,0,0],[0,1,0],[0,0,1]])
543  
-    def test_simple_0det(self):
544  
-        a=array([[1,2,3],[4,5,6.],[7,8,9]])
  542
+        assert_array_almost_equal(dot(a,a_pinv), np.eye(3))
  543
+
  544
+    def test_simple_complex(self):
  545
+        a = (array([[1, 2, 3], [4, 5, 6], [7, 8, 10]], dtype=float)
  546
+             + 1j * array([[10, 8, 7], [6, 5, 4], [3, 2, 1]], dtype=float))
  547
+        a_pinv = pinv(a)
  548
+        assert_array_almost_equal(dot(a, a_pinv), np.eye(3))
  549
+        a_pinv = pinv2(a)
  550
+        assert_array_almost_equal(dot(a, a_pinv), np.eye(3))
  551
+
  552
+    def test_simple_singular(self):
  553
+        a = array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=float)
545 554
         a_pinv = pinv(a)
546 555
         a_pinv2 = pinv2(a)
547 556
         assert_array_almost_equal(a_pinv,a_pinv2)
548 557
 
549 558
     def test_simple_cols(self):
550  
-        a=array([[1,2,3],[4,5,6.]])
  559
+        a = array([[1, 2, 3], [4, 5, 6]], dtype=float)
551 560
         a_pinv = pinv(a)
552 561
         a_pinv2 = pinv2(a)
553 562
         assert_array_almost_equal(a_pinv,a_pinv2)
554 563
 
555 564
     def test_simple_rows(self):
556  
-        a=array([[1,2],[3,4],[5,6]])
  565
+        a = array([[1, 2], [3, 4], [5, 6]], dtype=float)
557 566
         a_pinv = pinv(a)
558 567
         a_pinv2 = pinv2(a)
559 568
         assert_array_almost_equal(a_pinv,a_pinv2)
560 569
 
  570
+
  571
+class TestPinvSymmetric(TestCase):
  572
+
  573
+    def test_simple_real(self):
  574
+        a = array([[1, 2, 3], [4, 5, 6], [7, 8, 10]], dtype=float)
  575
+        a = np.dot(a, a.T)
  576
+        a_pinv = pinvh(a)
  577
+        assert_array_almost_equal(np.dot(a, a_pinv), np.eye(3))
  578
+
  579
+    def test_nonpositive(self):
  580
+        a = array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=float)
  581
+        a = np.dot(a, a.T)
  582
+        u, s, vt = np.linalg.svd(a)
  583
+        s[0] *= -1
  584
+        a = np.dot(u * s, vt)  # a is now symmetric non-positive and singular
  585
+        a_pinv = pinv2(a)
  586
+        a_pinvh = pinvh(a)
  587
+        assert_array_almost_equal(a_pinv, a_pinvh)
  588
+
  589
+    def test_simple_complex(self):
  590
+        a = (array([[1, 2, 3], [4, 5, 6], [7, 8, 10]], dtype=float)
  591
+             + 1j * array([[10, 8, 7], [6, 5, 4], [3, 2, 1]], dtype=float))
  592
+        a = np.dot(a, a.conj().T)
  593
+        a_pinv = pinvh(a)
  594
+        assert_array_almost_equal(np.dot(a, a_pinv), np.eye(3))
  595
+
  596
+
561 597
 class TestNorm(object):
562 598
 
563 599
     def test_types(self):
@@ -608,6 +644,8 @@ def test_pinv(self):
608 644
         assert_no_overwrite(pinv, [(3,3)])
609 645
     def test_pinv2(self):
610 646
         assert_no_overwrite(pinv2, [(3,3)])
  647
+    def test_pinvh(self):
  648
+        assert_no_overwrite(pinvh, [(3,3)])
611 649
 
612 650
 if __name__ == "__main__":
613 651
     run_module_suite()
Commit_comment_tip

Tip: You can add notes to lines in a file. Hover to the left of a line to make a note

Something went wrong with that request. Please try again.