[SciPy-Dev] Proposal to merge linalg.pinv and linalg.pinv2 and deprecate pinv2

[Stefan van der Walt]

On Sun, Mar 21, 2021, at 16:17, Robert Kern wrote:
> On Sun, Mar 21, 2021 at 6:53 PM Robert Kern <robert.kern at gmail.com> wrote:
>> On Sun, Mar 21, 2021 at 6:00 PM Ralf Gommers <ralf.gommers at gmail.com> wrote:
>>> 
>>> Do you happen to know the history of how we ended up with pinv2?
>> 
>> I suspect that when `pinv2()` was added, the `lstsq()` call underlying `pinv()` was not SVD-based. The precise LAPACK driver has changed over the years. We might have started with the QR-based driver.
> 
> It's going to be very hard to tell definitively because I think the history got lost in the SVN->git conversion due to some directory renames that happened in the early days. The pinv/pinv2 split seems to have been very early, though, so it may have dated from the original Multipack library (the source tarballs of which are also linkrotted away). 

I think you're right, this was pre-SVN.  I was looking at the following commit from https://github.com/scipy/scipy-svn

commit c6ef539392f31bda0b56541a1c8fdd61a0c0e6eb (HEAD)
Author: pearu <pearu at d6536bca-fef9-0310-8506-e4c0a848fbcf>
Date:   Sun Apr 7 15:03:50 2002 +0000

    Replacing linalg with linalg2: linalg->linalg/linalg1 and linalg2->linalg


There, the linalg/basic.py file is added, and inside it both pinv and pinv2 already exist:

def pinv(a, cond=-1):
    """Compute generalized inverse of A using least-squares solver.
    """
    a = asarray(a)
    t = a.typecode()
    b = scipy.identity(a.shape[0],t)
    return lstsq(a, b, cond=cond)[0]

def pinv2(a, cond=-1):
    """Compute the generalized inverse of A using svd.
    """
    a = asarray(a)
    u, s, vh = decomp.svd(a)
    m = u.shape[1]
    n = vh.shape[0]
    t = u.typecode()
    if cond is -1 or cond is None:
        cond = {0: feps*1e3, 1: eps*1e6}[_array_precision[t]]
    cutoff = cond*scipy_base.maximum.reduce(s)
    for i in range(min(n,m)):
        if s[i] > cutoff:
            s[i] = 1.0/s[i]
        else:
            s[i] = 0.0
    return dot(tran(conj(vh)),tran(conj(u))*s[:,NewAxis])


I have not been able to find a copy of multipack-0.7.tar.gz

Stéfan
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mail.python.org/pipermail/scipy-dev/attachments/20210322/41c964ff/attachment.html>