Hi All, I'm thinking about adding some functionality to lstsq because I find myself doing the same fixes over and over. List follows. 1. Add weights so data points can be weighted. 2. Use column scaling so condition numbers make more sense. 3. Compute covariance approximation? Unfortunately, the last will require using svd since there no linear least squares routines in LAPACK that also compute the covariance, at least that google knows about. Thoughts? Chuck
On Mon, Jul 19, 2010 at 5:50 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
Hi All,
I'm thinking about adding some functionality to lstsq because I find myself doing the same fixes over and over. List follows.
Add weights so data points can be weighted. Use column scaling so condition numbers make more sense. Compute covariance approximation?
Unfortunately, the last will require using svd since there no linear least squares routines in LAPACK that also compute the covariance, at least that google knows about.
Thoughts?
Chuck
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Maybe make 2 functions--one which implements 1 and 2, and another which implements 3? I think weights is an excellent idea! --Josh
On Mon, Jul 19, 2010 at 6:53 PM, Joshua Holbrook <josh.holbrook@gmail.com> wrote:
On Mon, Jul 19, 2010 at 5:50 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
Hi All,
I'm thinking about adding some functionality to lstsq because I find myself doing the same fixes over and over. List follows.
Add weights so data points can be weighted. Use column scaling so condition numbers make more sense. Compute covariance approximation?
Unfortunately, the last will require using svd since there no linear least squares routines in LAPACK that also compute the covariance, at least that google knows about.
Thoughts?
Maybe make 2 functions--one which implements 1 and 2, and another which implements 3? I think weights is an excellent idea!
I'd like a lstsq that did less, like not calculate the sum of squared residuals. That's useful in tight loops. So I also think having two lstsq makes sense. One as basic as possible; one with bells. How does scipy's lstsq fit into all this?
On Mon, Jul 19, 2010 at 9:02 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 5:50 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
Hi All,
I'm thinking about adding some functionality to lstsq because I find myself doing the same fixes over and over. List follows.
Add weights so data points can be weighted. Use column scaling so condition numbers make more sense. Compute covariance approximation?
Unfortunately, the last will require using svd since there no linear least squares routines in LAPACK that also compute the covariance, at least
On Mon, Jul 19, 2010 at 6:53 PM, Joshua Holbrook <josh.holbrook@gmail.com> wrote: that
google knows about.
Thoughts?
Maybe make 2 functions--one which implements 1 and 2, and another which implements 3? I think weights is an excellent idea!
I'd like a lstsq that did less, like not calculate the sum of squared residuals. That's useful in tight loops. So I also think having two lstsq makes sense. One as basic as possible; one with bells. How does scipy's lstsq fit into all this? _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On Mon, Jul 19, 2010 at 9:02 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 5:50 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
Hi All,
I'm thinking about adding some functionality to lstsq because I find myself doing the same fixes over and over. List follows.
Add weights so data points can be weighted. Use column scaling so condition numbers make more sense. Compute covariance approximation?
Unfortunately, the last will require using svd since there no linear least squares routines in LAPACK that also compute the covariance, at least
On Mon, Jul 19, 2010 at 6:53 PM, Joshua Holbrook <josh.holbrook@gmail.com> wrote: that
google knows about.
Thoughts?
Maybe make 2 functions--one which implements 1 and 2, and another which implements 3? I think weights is an excellent idea!
I'd like a lstsq that did less, like not calculate the sum of squared residuals. That's useful in tight loops. So I also think having two lstsq makes sense. One as basic as possible; one with bells. How does scipy's lstsq fit into all this?
I think the computation of the residues is cheap in lstsq. The algolrithm used starts by reducing the design matrix to bidiagonal from and reduces the rhs at the same time. In other words an mxn problem becomes a (n+1)xn problem. That's why the summed square of residuals is available but not the individual residuals, after the reduction there is only one residual and its square it the residue. Chuck
On Mon, Jul 19, 2010 at 8:27 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:02 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 6:53 PM, Joshua Holbrook <josh.holbrook@gmail.com> wrote:
On Mon, Jul 19, 2010 at 5:50 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
Hi All,
I'm thinking about adding some functionality to lstsq because I find myself doing the same fixes over and over. List follows.
Add weights so data points can be weighted. Use column scaling so condition numbers make more sense. Compute covariance approximation?
Unfortunately, the last will require using svd since there no linear least squares routines in LAPACK that also compute the covariance, at least that google knows about.
Thoughts?
Maybe make 2 functions--one which implements 1 and 2, and another which implements 3? I think weights is an excellent idea!
I'd like a lstsq that did less, like not calculate the sum of squared residuals. That's useful in tight loops. So I also think having two lstsq makes sense. One as basic as possible; one with bells. How does scipy's lstsq fit into all this?
I think the computation of the residues is cheap in lstsq. The algolrithm used starts by reducing the design matrix to bidiagonal from and reduces the rhs at the same time. In other words an mxn problem becomes a (n+1)xn problem. That's why the summed square of residuals is available but not the individual residuals, after the reduction there is only one residual and its square it the residue.
That does sound good. But it must take some time. There's indexing, array creation, if statement, summing: if results['rank'] == n and m > n: resids = sum((transpose(bstar)[n:,:])**2, axis=0).astype(result_t) Here are the timings after removing the sum of squared residuals:
x = np.random.rand(1000,10) y = np.random.rand(1000) timeit np.linalg.lstsq(x, y) 1000 loops, best of 3: 369 us per loop timeit np.linalg.lstsq2(x, y) 1000 loops, best of 3: 344 us per loop
x = np.random.rand(10,2) y = np.random.rand(10) timeit np.linalg.lstsq(x, y) 10000 loops, best of 3: 102 us per loop timeit np.linalg.lstsq2(x, y) 10000 loops, best of 3: 77 us per loop
On Mon, Jul 19, 2010 at 9:40 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 8:27 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:02 PM, Keith Goodman <kwgoodman@gmail.com>
On Mon, Jul 19, 2010 at 6:53 PM, Joshua Holbrook <josh.holbrook@gmail.com> wrote:
On Mon, Jul 19, 2010 at 5:50 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
Hi All,
I'm thinking about adding some functionality to lstsq because I find myself doing the same fixes over and over. List follows.
Add weights so data points can be weighted. Use column scaling so condition numbers make more sense. Compute covariance approximation?
Unfortunately, the last will require using svd since there no linear least squares routines in LAPACK that also compute the covariance, at least that google knows about.
Thoughts?
Maybe make 2 functions--one which implements 1 and 2, and another which implements 3? I think weights is an excellent idea!
I'd like a lstsq that did less, like not calculate the sum of squared residuals. That's useful in tight loops. So I also think having two lstsq makes sense. One as basic as possible; one with bells. How does scipy's lstsq fit into all this?
I think the computation of the residues is cheap in lstsq. The algolrithm used starts by reducing the design matrix to bidiagonal from and reduces
rhs at the same time. In other words an mxn problem becomes a (n+1)xn problem. That's why the summed square of residuals is available but not
wrote: the the
individual residuals, after the reduction there is only one residual and its square it the residue.
That does sound good. But it must take some time. There's indexing, array creation, if statement, summing:
if results['rank'] == n and m > n: resids = sum((transpose(bstar)[n:,:])**2, axis=0).astype(result_t)
Here are the timings after removing the sum of squared residuals:
x = np.random.rand(1000,10) y = np.random.rand(1000) timeit np.linalg.lstsq(x, y) 1000 loops, best of 3: 369 us per loop timeit np.linalg.lstsq2(x, y) 1000 loops, best of 3: 344 us per loop
x = np.random.rand(10,2) y = np.random.rand(10) timeit np.linalg.lstsq(x, y) 10000 loops, best of 3: 102 us per loop timeit np.linalg.lstsq2(x, y) 10000 loops, best of 3: 77 us per loop _
Now that I've looked through the driver program I see that you are right. Also, all the info needed for the covariance is almost available in the LAPACK driver program. Hmm, it seems that maybe the best thing to do here is to pull out the lapack_lite driver program, modify it, and make it a standalone python module that links to either the lapack_lite programs or the ATLAS versions. But that is more work than just doing things in python :-\ It does have the added advantage that all the work arrays can be handled internally. Chuck
On Tue, Jul 20, 2010 at 1:08 AM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:40 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 8:27 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:02 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 6:53 PM, Joshua Holbrook <josh.holbrook@gmail.com> wrote:
On Mon, Jul 19, 2010 at 5:50 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
Hi All,
I'm thinking about adding some functionality to lstsq because I find myself doing the same fixes over and over. List follows.
Add weights so data points can be weighted. Use column scaling so condition numbers make more sense. Compute covariance approximation?
Unfortunately, the last will require using svd since there no linear least squares routines in LAPACK that also compute the covariance, at least that google knows about.
Thoughts?
Maybe make 2 functions--one which implements 1 and 2, and another which implements 3? I think weights is an excellent idea!
I'd like a lstsq that did less, like not calculate the sum of squared residuals. That's useful in tight loops. So I also think having two lstsq makes sense. One as basic as possible; one with bells. How does scipy's lstsq fit into all this?
I think the computation of the residues is cheap in lstsq. The algolrithm used starts by reducing the design matrix to bidiagonal from and reduces the rhs at the same time. In other words an mxn problem becomes a (n+1)xn problem. That's why the summed square of residuals is available but not the individual residuals, after the reduction there is only one residual and its square it the residue.
That does sound good. But it must take some time. There's indexing, array creation, if statement, summing:
if results['rank'] == n and m > n: resids = sum((transpose(bstar)[n:,:])**2, axis=0).astype(result_t)
Here are the timings after removing the sum of squared residuals:
x = np.random.rand(1000,10) y = np.random.rand(1000) timeit np.linalg.lstsq(x, y) 1000 loops, best of 3: 369 us per loop timeit np.linalg.lstsq2(x, y) 1000 loops, best of 3: 344 us per loop
x = np.random.rand(10,2) y = np.random.rand(10) timeit np.linalg.lstsq(x, y) 10000 loops, best of 3: 102 us per loop timeit np.linalg.lstsq2(x, y) 10000 loops, best of 3: 77 us per loop _
Now that I've looked through the driver program I see that you are right. Also, all the info needed for the covariance is almost available in the LAPACK driver program. Hmm, it seems that maybe the best thing to do here is to pull out the lapack_lite driver program, modify it, and make it a standalone python module that links to either the lapack_lite programs or the ATLAS versions. But that is more work than just doing things in python :-\ It does have the added advantage that all the work arrays can be handled internally.
IIUC, I tried this with statsmodels earlier this summer, but found the result to be slower than our current implementation (and surely slower than lstsq since we have a bit more Python overhead around the LAPACK function calls). I will have to dig through the revisions and make sure about this. But for the record, we have weighted least squares (WLS) and you can get the covariance matrix. I think you can do column-scaling with our GLS (generalized least squares class). I haven't seen the terms row-scaling and column-scaling used much in econometrics books do you have a good reference offhand? Skipper
On Tue, Jul 20, 2010 at 8:32 AM, Skipper Seabold <jsseabold@gmail.com>wrote:
On Tue, Jul 20, 2010 at 1:08 AM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:40 PM, Keith Goodman <kwgoodman@gmail.com>
On Mon, Jul 19, 2010 at 8:27 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:02 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 6:53 PM, Joshua Holbrook <josh.holbrook@gmail.com> wrote:
On Mon, Jul 19, 2010 at 5:50 PM, Charles R Harris <charlesr.harris@gmail.com> wrote: > Hi All, > > I'm thinking about adding some functionality to lstsq because I
find
> myself > doing the same fixes over and over. List follows. > > Add weights so data points can be weighted. > Use column scaling so condition numbers make more sense. > Compute covariance approximation? > > Unfortunately, the last will require using svd since there no
> least > squares routines in LAPACK that also compute the covariance, at > least > that > google knows about. > > Thoughts?
Maybe make 2 functions--one which implements 1 and 2, and another which implements 3? I think weights is an excellent idea!
I'd like a lstsq that did less, like not calculate the sum of squared residuals. That's useful in tight loops. So I also think having two lstsq makes sense. One as basic as possible; one with bells. How does scipy's lstsq fit into all this?
I think the computation of the residues is cheap in lstsq. The algolrithm used starts by reducing the design matrix to bidiagonal from and reduces the rhs at the same time. In other words an mxn problem becomes a (n+1)xn problem. That's why the summed square of residuals is available but not the individual residuals, after the reduction there is only one residual and its square it the residue.
That does sound good. But it must take some time. There's indexing, array creation, if statement, summing:
if results['rank'] == n and m > n: resids = sum((transpose(bstar)[n:,:])**2, axis=0).astype(result_t)
Here are the timings after removing the sum of squared residuals:
x = np.random.rand(1000,10) y = np.random.rand(1000) timeit np.linalg.lstsq(x, y) 1000 loops, best of 3: 369 us per loop timeit np.linalg.lstsq2(x, y) 1000 loops, best of 3: 344 us per loop
x = np.random.rand(10,2) y = np.random.rand(10) timeit np.linalg.lstsq(x, y) 10000 loops, best of 3: 102 us per loop timeit np.linalg.lstsq2(x, y) 10000 loops, best of 3: 77 us per loop _
Now that I've looked through the driver program I see that you are right. Also, all the info needed for the covariance is almost available in the LAPACK driver program. Hmm, it seems that maybe the best thing to do here is to pull out the lapack_lite driver program, modify it, and make it a standalone python module that links to either the lapack_lite programs or the ATLAS versions. But that is more work than just doing things in
wrote: linear python
:-\ It does have the added advantage that all the work arrays can be handled internally.
IIUC, I tried this with statsmodels earlier this summer, but found the result to be slower than our current implementation (and surely slower than lstsq since we have a bit more Python overhead around the LAPACK function calls). I will have to dig through the revisions and make sure about this.
You modified the C version of dgelsd in lapack_lite? That is the driver I'm talking about.
But for the record, we have weighted least squares (WLS) and you can get the covariance matrix. I think you can do column-scaling with our GLS (generalized least squares class). I haven't seen the terms row-scaling and column-scaling used much in econometrics books do you have a good reference offhand?
The dgelsd driver scales the matrix, but only enough to keep the numbers in the working range. Column scaling affects the singular values, so if you are using them to track the condition number of the matrix it can effect the outcome. The deeper thing that is going on is that some of the usefulness of the svd is based on the assumption that the euclidean distance between parameter vectors means something. I just tend to normalize the columns to unit length. It is a bit like using z-scores in statistics except the column norm is normalized instead of the column variance. I'm sure there are references out there but a quick google doesn't show anything obvious. Normalization seems to be a bit neglected in the usual classroom notes. Chuck
On Tue, Jul 20, 2010 at 11:23 AM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Tue, Jul 20, 2010 at 8:32 AM, Skipper Seabold <jsseabold@gmail.com> wrote:
On Tue, Jul 20, 2010 at 1:08 AM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:40 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 8:27 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:02 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 6:53 PM, Joshua Holbrook <josh.holbrook@gmail.com> wrote: > On Mon, Jul 19, 2010 at 5:50 PM, Charles R Harris > <charlesr.harris@gmail.com> wrote: >> Hi All, >> >> I'm thinking about adding some functionality to lstsq because I >> find >> myself >> doing the same fixes over and over. List follows. >> >> Add weights so data points can be weighted. >> Use column scaling so condition numbers make more sense. >> Compute covariance approximation? >> >> Unfortunately, the last will require using svd since there no >> linear >> least >> squares routines in LAPACK that also compute the covariance, at >> least >> that >> google knows about. >> >> Thoughts? > > Maybe make 2 functions--one which implements 1 and 2, and another > which implements 3? I think weights is an excellent idea!
I'd like a lstsq that did less, like not calculate the sum of squared residuals. That's useful in tight loops. So I also think having two lstsq makes sense. One as basic as possible; one with bells. How does scipy's lstsq fit into all this?
I think the computation of the residues is cheap in lstsq. The algolrithm used starts by reducing the design matrix to bidiagonal from and reduces the rhs at the same time. In other words an mxn problem becomes a (n+1)xn problem. That's why the summed square of residuals is available but not the individual residuals, after the reduction there is only one residual and its square it the residue.
That does sound good. But it must take some time. There's indexing, array creation, if statement, summing:
if results['rank'] == n and m > n: resids = sum((transpose(bstar)[n:,:])**2, axis=0).astype(result_t)
Here are the timings after removing the sum of squared residuals:
x = np.random.rand(1000,10) y = np.random.rand(1000) timeit np.linalg.lstsq(x, y) 1000 loops, best of 3: 369 us per loop timeit np.linalg.lstsq2(x, y) 1000 loops, best of 3: 344 us per loop
x = np.random.rand(10,2) y = np.random.rand(10) timeit np.linalg.lstsq(x, y) 10000 loops, best of 3: 102 us per loop timeit np.linalg.lstsq2(x, y) 10000 loops, best of 3: 77 us per loop _
Now that I've looked through the driver program I see that you are right. Also, all the info needed for the covariance is almost available in the LAPACK driver program. Hmm, it seems that maybe the best thing to do here is to pull out the lapack_lite driver program, modify it, and make it a standalone python module that links to either the lapack_lite programs or the ATLAS versions. But that is more work than just doing things in python :-\ It does have the added advantage that all the work arrays can be handled internally.
IIUC, I tried this with statsmodels earlier this summer, but found the result to be slower than our current implementation (and surely slower than lstsq since we have a bit more Python overhead around the LAPACK function calls). I will have to dig through the revisions and make sure about this.
You modified the C version of dgelsd in lapack_lite? That is the driver I'm talking about.
Ah, then, no, I didn't understand correctly, but I see it in the source now.
But for the record, we have weighted least squares (WLS) and you can get the covariance matrix. I think you can do column-scaling with our GLS (generalized least squares class). I haven't seen the terms row-scaling and column-scaling used much in econometrics books do you have a good reference offhand?
The dgelsd driver scales the matrix, but only enough to keep the numbers in the working range. Column scaling affects the singular values, so if you are using them to track the condition number of the matrix it can effect the outcome. The deeper thing that is going on is that some of the usefulness of the svd is based on the assumption that the euclidean distance between parameter vectors means something. I just tend to normalize the columns to unit length. It is a bit like using z-scores in statistics except the column norm is normalized instead of the column variance. I'm sure there are references out there but a quick google doesn't show anything obvious. Normalization seems to be a bit neglected in the usual classroom notes.
That makes sense. Thanks, Skipper
On Mon, Jul 19, 2010 at 10:08 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:40 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 8:27 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:02 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 6:53 PM, Joshua Holbrook <josh.holbrook@gmail.com> wrote:
On Mon, Jul 19, 2010 at 5:50 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
Hi All,
I'm thinking about adding some functionality to lstsq because I find myself doing the same fixes over and over. List follows.
Add weights so data points can be weighted. Use column scaling so condition numbers make more sense. Compute covariance approximation?
Unfortunately, the last will require using svd since there no linear least squares routines in LAPACK that also compute the covariance, at least that google knows about.
Thoughts?
Maybe make 2 functions--one which implements 1 and 2, and another which implements 3? I think weights is an excellent idea!
I'd like a lstsq that did less, like not calculate the sum of squared residuals. That's useful in tight loops. So I also think having two lstsq makes sense. One as basic as possible; one with bells. How does scipy's lstsq fit into all this?
I think the computation of the residues is cheap in lstsq. The algolrithm used starts by reducing the design matrix to bidiagonal from and reduces the rhs at the same time. In other words an mxn problem becomes a (n+1)xn problem. That's why the summed square of residuals is available but not the individual residuals, after the reduction there is only one residual and its square it the residue.
That does sound good. But it must take some time. There's indexing, array creation, if statement, summing:
if results['rank'] == n and m > n: resids = sum((transpose(bstar)[n:,:])**2, axis=0).astype(result_t)
Here are the timings after removing the sum of squared residuals:
x = np.random.rand(1000,10) y = np.random.rand(1000) timeit np.linalg.lstsq(x, y) 1000 loops, best of 3: 369 us per loop timeit np.linalg.lstsq2(x, y) 1000 loops, best of 3: 344 us per loop
x = np.random.rand(10,2) y = np.random.rand(10) timeit np.linalg.lstsq(x, y) 10000 loops, best of 3: 102 us per loop timeit np.linalg.lstsq2(x, y) 10000 loops, best of 3: 77 us per loop _
Now that I've looked through the driver program I see that you are right. Also, all the info needed for the covariance is almost available in the LAPACK driver program. Hmm, it seems that maybe the best thing to do here is to pull out the lapack_lite driver program, modify it, and make it a standalone python module that links to either the lapack_lite programs or the ATLAS versions. But that is more work than just doing things in python :-\ It does have the added advantage that all the work arrays can be handled internally.
While taking a quick look at lstsq to see what I could cut, this line caught my eye: bstar[:b.shape[0],:n_rhs] = b.copy() I thought that array.__setitem__ never gives a view of the right hand side. If that's so then the copy is not needed. That would save some time since b can be large. All unit test pass when I remove the copy, but you know how that goes... I also noticed that "import math" is done inside the lstsq function. Why is that?
On Tue, Jul 20, 2010 at 6:35 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 10:08 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:40 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 8:27 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:02 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 6:53 PM, Joshua Holbrook <josh.holbrook@gmail.com> wrote:
On Mon, Jul 19, 2010 at 5:50 PM, Charles R Harris <charlesr.harris@gmail.com> wrote: > Hi All, > > I'm thinking about adding some functionality to lstsq because I find > myself > doing the same fixes over and over. List follows. > > Add weights so data points can be weighted. > Use column scaling so condition numbers make more sense. > Compute covariance approximation? > > Unfortunately, the last will require using svd since there no linear > least > squares routines in LAPACK that also compute the covariance, at > least > that > google knows about. > > Thoughts?
Maybe make 2 functions--one which implements 1 and 2, and another which implements 3? I think weights is an excellent idea!
I'd like a lstsq that did less, like not calculate the sum of squared residuals. That's useful in tight loops. So I also think having two lstsq makes sense. One as basic as possible; one with bells. How does scipy's lstsq fit into all this?
I think the computation of the residues is cheap in lstsq. The algolrithm used starts by reducing the design matrix to bidiagonal from and reduces the rhs at the same time. In other words an mxn problem becomes a (n+1)xn problem. That's why the summed square of residuals is available but not the individual residuals, after the reduction there is only one residual and its square it the residue.
That does sound good. But it must take some time. There's indexing, array creation, if statement, summing:
if results['rank'] == n and m > n: resids = sum((transpose(bstar)[n:,:])**2, axis=0).astype(result_t)
Here are the timings after removing the sum of squared residuals:
x = np.random.rand(1000,10) y = np.random.rand(1000) timeit np.linalg.lstsq(x, y) 1000 loops, best of 3: 369 us per loop timeit np.linalg.lstsq2(x, y) 1000 loops, best of 3: 344 us per loop
x = np.random.rand(10,2) y = np.random.rand(10) timeit np.linalg.lstsq(x, y) 10000 loops, best of 3: 102 us per loop timeit np.linalg.lstsq2(x, y) 10000 loops, best of 3: 77 us per loop _
Now that I've looked through the driver program I see that you are right. Also, all the info needed for the covariance is almost available in the LAPACK driver program. Hmm, it seems that maybe the best thing to do here is to pull out the lapack_lite driver program, modify it, and make it a standalone python module that links to either the lapack_lite programs or the ATLAS versions. But that is more work than just doing things in python :-\ It does have the added advantage that all the work arrays can be handled internally.
While taking a quick look at lstsq to see what I could cut, this line caught my eye:
bstar[:b.shape[0],:n_rhs] = b.copy()
Even if bstar contained a view of b, the next line in lstsq would take care of it: a, bstar = _fastCopyAndTranspose(t, a, bstar)
I thought that array.__setitem__ never gives a view of the right hand side. If that's so then the copy is not needed. That would save some time since b can be large. All unit test pass when I remove the copy, but you know how that goes...
I also noticed that "import math" is done inside the lstsq function. Why is that?
On Tue, Jul 20, 2010 at 10:12 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Tue, Jul 20, 2010 at 6:35 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 10:08 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:40 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 8:27 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:02 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 6:53 PM, Joshua Holbrook <josh.holbrook@gmail.com> wrote: > On Mon, Jul 19, 2010 at 5:50 PM, Charles R Harris > <charlesr.harris@gmail.com> wrote: >> Hi All, >> >> I'm thinking about adding some functionality to lstsq because I find >> myself >> doing the same fixes over and over. List follows. >> >> Add weights so data points can be weighted. >> Use column scaling so condition numbers make more sense. >> Compute covariance approximation? >> >> Unfortunately, the last will require using svd since there no linear >> least >> squares routines in LAPACK that also compute the covariance, at >> least >> that >> google knows about. >> >> Thoughts? > > Maybe make 2 functions--one which implements 1 and 2, and another > which implements 3? I think weights is an excellent idea!
I'd like a lstsq that did less, like not calculate the sum of squared residuals. That's useful in tight loops. So I also think having two lstsq makes sense. One as basic as possible; one with bells. How does scipy's lstsq fit into all this?
I think the computation of the residues is cheap in lstsq. The algolrithm used starts by reducing the design matrix to bidiagonal from and reduces the rhs at the same time. In other words an mxn problem becomes a (n+1)xn problem. That's why the summed square of residuals is available but not the individual residuals, after the reduction there is only one residual and its square it the residue.
That does sound good. But it must take some time. There's indexing, array creation, if statement, summing:
if results['rank'] == n and m > n: resids = sum((transpose(bstar)[n:,:])**2, axis=0).astype(result_t)
Here are the timings after removing the sum of squared residuals:
x = np.random.rand(1000,10) y = np.random.rand(1000) timeit np.linalg.lstsq(x, y) 1000 loops, best of 3: 369 us per loop timeit np.linalg.lstsq2(x, y) 1000 loops, best of 3: 344 us per loop
x = np.random.rand(10,2) y = np.random.rand(10) timeit np.linalg.lstsq(x, y) 10000 loops, best of 3: 102 us per loop timeit np.linalg.lstsq2(x, y) 10000 loops, best of 3: 77 us per loop _
Now that I've looked through the driver program I see that you are right. Also, all the info needed for the covariance is almost available in the LAPACK driver program. Hmm, it seems that maybe the best thing to do here is to pull out the lapack_lite driver program, modify it, and make it a standalone python module that links to either the lapack_lite programs or the ATLAS versions. But that is more work than just doing things in python :-\ It does have the added advantage that all the work arrays can be handled internally.
While taking a quick look at lstsq to see what I could cut, this line caught my eye:
bstar[:b.shape[0],:n_rhs] = b.copy()
Even if bstar contained a view of b, the next line in lstsq would take care of it:
a, bstar = _fastCopyAndTranspose(t, a, bstar)
If a view is never taken, then what does the copy method actually do? A redundant copy? The b matrix to *gelsd is in/out, so if the original b is not going to be written to then taking out the redundant copy seems like a good idea. Out of curiosity, is there an explicit way to check if these share memory? import numpy as np b = np.empty((1000,2)) c = np.arange(1000)[:,None] b[:,1:] = c Then what? Skipper
I thought that array.__setitem__ never gives a view of the right hand side. If that's so then the copy is not needed. That would save some time since b can be large. All unit test pass when I remove the copy, but you know how that goes...
I also noticed that "import math" is done inside the lstsq function. Why is that?
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On Tue, Jul 20, 2010 at 7:16 PM, Skipper Seabold <jsseabold@gmail.com> wrote:
On Tue, Jul 20, 2010 at 10:12 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Tue, Jul 20, 2010 at 6:35 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 10:08 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:40 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Mon, Jul 19, 2010 at 8:27 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Mon, Jul 19, 2010 at 9:02 PM, Keith Goodman <kwgoodman@gmail.com> wrote: > > On Mon, Jul 19, 2010 at 6:53 PM, Joshua Holbrook > <josh.holbrook@gmail.com> wrote: > > On Mon, Jul 19, 2010 at 5:50 PM, Charles R Harris > > <charlesr.harris@gmail.com> wrote: > >> Hi All, > >> > >> I'm thinking about adding some functionality to lstsq because I find > >> myself > >> doing the same fixes over and over. List follows. > >> > >> Add weights so data points can be weighted. > >> Use column scaling so condition numbers make more sense. > >> Compute covariance approximation? > >> > >> Unfortunately, the last will require using svd since there no linear > >> least > >> squares routines in LAPACK that also compute the covariance, at > >> least > >> that > >> google knows about. > >> > >> Thoughts? > > > > Maybe make 2 functions--one which implements 1 and 2, and another > > which implements 3? I think weights is an excellent idea! > > I'd like a lstsq that did less, like not calculate the sum of squared > residuals. That's useful in tight loops. So I also think having two > lstsq makes sense. One as basic as possible; one with bells. How does > scipy's lstsq fit into all this?
I think the computation of the residues is cheap in lstsq. The algolrithm used starts by reducing the design matrix to bidiagonal from and reduces the rhs at the same time. In other words an mxn problem becomes a (n+1)xn problem. That's why the summed square of residuals is available but not the individual residuals, after the reduction there is only one residual and its square it the residue.
That does sound good. But it must take some time. There's indexing, array creation, if statement, summing:
if results['rank'] == n and m > n: resids = sum((transpose(bstar)[n:,:])**2, axis=0).astype(result_t)
Here are the timings after removing the sum of squared residuals:
> x = np.random.rand(1000,10) > y = np.random.rand(1000) > timeit np.linalg.lstsq(x, y) 1000 loops, best of 3: 369 us per loop > timeit np.linalg.lstsq2(x, y) 1000 loops, best of 3: 344 us per loop > > x = np.random.rand(10,2) > y = np.random.rand(10) > timeit np.linalg.lstsq(x, y) 10000 loops, best of 3: 102 us per loop > timeit np.linalg.lstsq2(x, y) 10000 loops, best of 3: 77 us per loop _
Now that I've looked through the driver program I see that you are right. Also, all the info needed for the covariance is almost available in the LAPACK driver program. Hmm, it seems that maybe the best thing to do here is to pull out the lapack_lite driver program, modify it, and make it a standalone python module that links to either the lapack_lite programs or the ATLAS versions. But that is more work than just doing things in python :-\ It does have the added advantage that all the work arrays can be handled internally.
While taking a quick look at lstsq to see what I could cut, this line caught my eye:
bstar[:b.shape[0],:n_rhs] = b.copy()
Even if bstar contained a view of b, the next line in lstsq would take care of it:
a, bstar = _fastCopyAndTranspose(t, a, bstar)
If a view is never taken, then what does the copy method actually do? A redundant copy? The b matrix to *gelsd is in/out, so if the original b is not going to be written to then taking out the redundant copy seems like a good idea.
Good point. Looks like we can get rid of 2 copies! I didn't get rid of the second copy but I did cut things down just to see what the timing was like. I also threw out the ability to handle complex numbers (only saves an if iscomplex statement):
x = np.random.rand(10,2) y = np.random.rand(10)
timeit np.linalg.lstsq(x, y) 10000 loops, best of 3: 99.6 us per loop timeit np.linalg.lstsq3(x, y) 10000 loops, best of 3: 61 us per loop
Here's what it looks like (any other tweaks possible?): import math def lstsq3(a, b, rcond=-1): a, _ = _makearray(a) b, wrap = _makearray(b) is_1d = len(b.shape) == 1 if is_1d: b = b[:, newaxis] _assertRank2(a, b) m = a.shape[0] n = a.shape[1] n_rhs = b.shape[1] ldb = max(n, m) if m != b.shape[0]: raise LinAlgError, 'Incompatible dimensions' t, result_t = _commonType(a, b) real_t = _linalgRealType(t) bstar = zeros((ldb, n_rhs), t) bstar[:b.shape[0],:n_rhs] = b a, bstar = _fastCopyAndTranspose(t, a, bstar) s = zeros((min(m, n),), real_t) nlvl = max( 0, int( math.log( float(min(m, n))/2. ) ) + 1 ) iwork = zeros((3*min(m, n)*nlvl+11*min(m, n),), fortran_int) lapack_routine = lapack_lite.dgelsd lwork = 1 work = zeros((lwork,), t) results = lapack_routine(m, n, n_rhs, a, m, bstar, ldb, s, rcond, 0, work, -1, iwork, 0) lwork = int(work[0]) work = zeros((lwork,), t) results = lapack_routine(m, n, n_rhs, a, m, bstar, ldb, s, rcond, 0, work, lwork, iwork, 0) if results['info'] > 0: raise LinAlgError, 'SVD did not converge in Linear Least Squares' if is_1d: x = array(ravel(bstar)[:n], dtype=result_t, copy=True) else: x = array(transpose(bstar)[:n,:], dtype=result_t, copy=True) return wrap(x)
Out of curiosity, is there an explicit way to check if these share memory?
import numpy as np b = np.empty((1000,2)) c = np.arange(1000)[:,None] b[:,1:] = c
Then what?
Skipper
I thought that array.__setitem__ never gives a view of the right hand side. If that's so then the copy is not needed. That would save some time since b can be large. All unit test pass when I remove the copy, but you know how that goes...
I also noticed that "import math" is done inside the lstsq function. Why is that?
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On Tue, Jul 20, 2010 at 10:24 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
Good point. Looks like we can get rid of 2 copies! I didn't get rid of the second copy but I did cut things down just to see what the timing was like. I also threw out the ability to handle complex numbers (only saves an if iscomplex statement):
Just for timing purposes though right?
On Wed, Jul 21, 2010 at 4:35 PM, Skipper Seabold <jsseabold@gmail.com> wrote:
On Tue, Jul 20, 2010 at 10:24 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
Good point. Looks like we can get rid of 2 copies! I didn't get rid of the second copy but I did cut things down just to see what the timing was like. I also threw out the ability to handle complex numbers (only saves an if iscomplex statement):
Just for timing purposes though right?
Yes. Can someone confirm that the copy in np.linalg.lstsq bstar[:b.shape[0],:n_rhs] = b.copy() is not needed? I'm assuming that ndarray.__setitem__ never gives a view of the right-hand side.
On Wed, Jul 21, 2010 at 5:44 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
On Wed, Jul 21, 2010 at 4:35 PM, Skipper Seabold <jsseabold@gmail.com> wrote:
On Tue, Jul 20, 2010 at 10:24 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
Good point. Looks like we can get rid of 2 copies! I didn't get rid of the second copy but I did cut things down just to see what the timing was like. I also threw out the ability to handle complex numbers (only saves an if iscomplex statement):
Just for timing purposes though right?
Yes.
Can someone confirm that the copy in np.linalg.lstsq
bstar[:b.shape[0],:n_rhs] = b.copy()
is not needed? I'm assuming that ndarray.__setitem__ never gives a view of the right-hand side. ________
bstar doesn't share memory with b, so there is no need for a copy. The whole first part of lstsq could use a cleanup, there are a lot of other things that could be made cleaner and some other bits that don't look right. Anyone want to take a shot at it? Chuck
On Wed, Jul 21, 2010 at 8:24 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Jul 21, 2010 at 5:44 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
Can someone confirm that the copy in np.linalg.lstsq
bstar[:b.shape[0],:n_rhs] = b.copy()
is not needed? I'm assuming that ndarray.__setitem__ never gives a view of the right-hand side. ________
bstar doesn't share memory with b, so there is no need for a copy. The whole first part of lstsq could use a cleanup, there are a lot of other things that could be made cleaner and some other bits that don't look right. Anyone want to take a shot at it?
Chuck
Some minor tweaks suggestions:
b = np.random.rand(1000)
timeit b[:,newaxis] # actual 1000000 loops, best of 3: 1.03 us per loop timeit b.reshape(-1,1) # alternative 1000000 loops, best of 3: 445 ns per loop
timeit len(b.shape) # actual 10000000 loops, best of 3: 145 ns per loop timeit b.ndim # alternative 10000000 loops, best of 3: 74.1 ns per loop
On 2010-07-20, at 10:16 PM, Skipper Seabold wrote:
Out of curiosity, is there an explicit way to check if these share memory?
You could do the exact calculations (I think) but this isn't actually implemented in NumPy, though np.may_share_memory is a conservative test for it that will err on the side of false positives. David
On Wed, Jul 21, 2010 at 7:25 PM, David Warde-Farley <dwf@cs.toronto.edu> wrote:
On 2010-07-20, at 10:16 PM, Skipper Seabold wrote:
Out of curiosity, is there an explicit way to check if these share memory?
You could do the exact calculations (I think) but this isn't actually implemented in NumPy, though np.may_share_memory is a conservative test for it that will err on the side of false positives.
Ah, that's good to know. In [5]: np.may_share_memory(b,c) Out[5]: False In [6]: a = b.view() In [7]: np.may_share_memory(a,b) Out[7]: True In [8]: np.may_share_memory(a,c) Out[8]: False
participants (5)
-
Charles R Harris -
David Warde-Farley -
Joshua Holbrook -
Keith Goodman -
Skipper Seabold