Why does orth use svd instead of QR ?
Hi, I wanted to know if there was a rationale for using svd to orthonormalize the columns of a matrix (in scipy.linalg). QR-based methods are likely to be much faster, and I thought this was the standard, numerically-stable method to orthonormalize a basis ? If the reason is to deal with rank-deficient matrices, maybe we could add an option to choose between them ? cheers, David
On Thu, Feb 4, 2010 at 8:45 PM, David Cournapeau <david@silveregg.co.jp>wrote:
Hi,
I wanted to know if there was a rationale for using svd to orthonormalize the columns of a matrix (in scipy.linalg). QR-based methods are likely to be much faster, and I thought this was the standard, numerically-stable method to orthonormalize a basis ? If the reason is to deal with rank-deficient matrices, maybe we could add an option to choose between them ?
QR with column rotation would deal with rank-deficient matrices and routines for that are available in LAPACK <http://netlib.org/lapack/lug/node42.html>. The SVD was probably used because it was available. The diagonal elements of the R matrix can somewhat take the place of the singular values when column rotation is used. Chuck
Charles R Harris wrote:
On Thu, Feb 4, 2010 at 8:45 PM, David Cournapeau <david@silveregg.co.jp <mailto:david@silveregg.co.jp>> wrote:
Hi,
I wanted to know if there was a rationale for using svd to orthonormalize the columns of a matrix (in scipy.linalg). QR-based methods are likely to be much faster, and I thought this was the standard, numerically-stable method to orthonormalize a basis ? If the reason is to deal with rank-deficient matrices, maybe we could add an option to choose between them ?
QR with column rotation would deal with rank-deficient matrices and routines for that are available in LAPACK <http://netlib.org/lapack/lug/node42.html>. The SVD was probably used because it was available. The diagonal elements of the R matrix can somewhat take the place of the singular values when column rotation is used.
So would be it ok to use this column-rotated QR in place of svd for every case in orth ? I would have to check that QR with column rotation is still significantly faster than svd, but I would surprised if if were not the case. QR has also the advantage of being implemented in PLASMA already contrary to eigen/svd solvers, cheers, David
On Fri, Feb 5, 2010 at 12:18 AM, David Cournapeau <david@silveregg.co.jp>wrote:
Charles R Harris wrote:
On Thu, Feb 4, 2010 at 8:45 PM, David Cournapeau <david@silveregg.co.jp <mailto:david@silveregg.co.jp>> wrote:
Hi,
I wanted to know if there was a rationale for using svd to orthonormalize the columns of a matrix (in scipy.linalg). QR-based methods are likely to be much faster, and I thought this was the standard, numerically-stable method to orthonormalize a basis ? If
the
reason is to deal with rank-deficient matrices, maybe we could add an option to choose between them ?
QR with column rotation would deal with rank-deficient matrices and routines for that are available in LAPACK <http://netlib.org/lapack/lug/node42.html>. The SVD was probably used because it was available. The diagonal elements of the R matrix can somewhat take the place of the singular values when column rotation is
used.
So would be it ok to use this column-rotated QR in place of svd for every case in orth ? I would have to check that QR with column rotation is still significantly faster than svd, but I would surprised if if were not the case. QR has also the advantage of being implemented in PLASMA already contrary to eigen/svd solvers,
I don't know how the two methods compare in practice. SVD algorithms generally use iterated QR reductions in their implementation, so QR reductions can't be worse numerically. But the SVD probably provides a better metric for rank determination. A google search turns up some literature on the subject that I can't access from home. Chuck
On Fri, Feb 5, 2010 at 12:47 AM, Charles R Harris <charlesr.harris@gmail.com
wrote:
On Fri, Feb 5, 2010 at 12:18 AM, David Cournapeau <david@silveregg.co.jp>wrote:
Charles R Harris wrote:
On Thu, Feb 4, 2010 at 8:45 PM, David Cournapeau <david@silveregg.co.jp <mailto:david@silveregg.co.jp>> wrote:
Hi,
I wanted to know if there was a rationale for using svd to orthonormalize the columns of a matrix (in scipy.linalg). QR-based methods are likely to be much faster, and I thought this was the standard, numerically-stable method to orthonormalize a basis ? If
the
reason is to deal with rank-deficient matrices, maybe we could add
an
option to choose between them ?
QR with column rotation would deal with rank-deficient matrices and routines for that are available in LAPACK <http://netlib.org/lapack/lug/node42.html>. The SVD was probably used because it was available. The diagonal elements of the R matrix can somewhat take the place of the singular values when column rotation is
used.
So would be it ok to use this column-rotated QR in place of svd for every case in orth ? I would have to check that QR with column rotation is still significantly faster than svd, but I would surprised if if were not the case. QR has also the advantage of being implemented in PLASMA already contrary to eigen/svd solvers,
I don't know how the two methods compare in practice. SVD algorithms generally use iterated QR reductions in their implementation, so QR reductions can't be worse numerically. But the SVD probably provides a better metric for rank determination. A google search turns up some literature on the subject that I can't access from home.
OK, here's a good reference<http://www.math.sjsu.edu/%7Efoster/rank/rank_revealing_s.pdf>. A quick look seems to indicate that the SVD is the way to go. Chuck
Charles R Harris wrote:
On Fri, Feb 5, 2010 at 12:47 AM, Charles R Harris <charlesr.harris@gmail.com <mailto:charlesr.harris@gmail.com>> wrote:
On Fri, Feb 5, 2010 at 12:18 AM, David Cournapeau <david@silveregg.co.jp <mailto:david@silveregg.co.jp>> wrote:
Charles R Harris wrote: > > > On Thu, Feb 4, 2010 at 8:45 PM, David Cournapeau <david@silveregg.co.jp <mailto:david@silveregg.co.jp> > <mailto:david@silveregg.co.jp <mailto:david@silveregg.co.jp>>> wrote: > > Hi, > > I wanted to know if there was a rationale for using svd to > orthonormalize the columns of a matrix (in scipy.linalg). QR-based > methods are likely to be much faster, and I thought this was the > standard, numerically-stable method to orthonormalize a basis ? If the > reason is to deal with rank-deficient matrices, maybe we could add an > option to choose between them ? > > > QR with column rotation would deal with rank-deficient matrices and > routines for that are available in LAPACK > <http://netlib.org/lapack/lug/node42.html>. The SVD was probably used > because it was available. The diagonal elements of the R matrix can > somewhat take the place of the singular values when column rotation is used.
So would be it ok to use this column-rotated QR in place of svd for every case in orth ? I would have to check that QR with column rotation is still significantly faster than svd, but I would surprised if if were not the case. QR has also the advantage of being implemented in PLASMA already contrary to eigen/svd solvers,
I don't know how the two methods compare in practice. SVD algorithms generally use iterated QR reductions in their implementation, so QR reductions can't be worse numerically. But the SVD probably provides a better metric for rank determination. A google search turns up some literature on the subject that I can't access from home.
OK, here's a good reference <http://www.math.sjsu.edu/%7Efoster/rank/rank_revealing_s.pdf>. A quick look seems to indicate that the SVD is the way to go.
AFAIK, SVD is indeed the way to go, but do we really need this for the orth function ? I am wrong to think that orthonormalizing a matrix of linearly independent vectors is the most common usage for orth ? The difference in terms of speed is really significant (for example, svd of a 2000x100 matrix takes ~1.9 second vs 0.1 s for QR). cheers, David
On Fri, Feb 5, 2010 at 1:29 AM, David Cournapeau <david@silveregg.co.jp>wrote:
Charles R Harris wrote:
On Fri, Feb 5, 2010 at 12:47 AM, Charles R Harris <charlesr.harris@gmail.com <mailto:charlesr.harris@gmail.com>> wrote:
On Fri, Feb 5, 2010 at 12:18 AM, David Cournapeau <david@silveregg.co.jp <mailto:david@silveregg.co.jp>> wrote:
Charles R Harris wrote: > > > On Thu, Feb 4, 2010 at 8:45 PM, David Cournapeau <david@silveregg.co.jp <mailto:david@silveregg.co.jp> > <mailto:david@silveregg.co.jp <mailto:david@silveregg.co.jp>>> wrote: > > Hi, > > I wanted to know if there was a rationale for using svd to > orthonormalize the columns of a matrix (in scipy.linalg). QR-based > methods are likely to be much faster, and I thought this was the > standard, numerically-stable method to orthonormalize a basis ? If the > reason is to deal with rank-deficient matrices, maybe we could add an > option to choose between them ? > > > QR with column rotation would deal with rank-deficient matrices and > routines for that are available in LAPACK > <http://netlib.org/lapack/lug/node42.html>. The SVD was probably used > because it was available. The diagonal elements of the R matrix can > somewhat take the place of the singular values when column rotation is used.
So would be it ok to use this column-rotated QR in place of svd
for
every case in orth ? I would have to check that QR with column rotation is still significantly faster than svd, but I would surprised if if were not the case. QR has also the advantage of being implemented in PLASMA already contrary to eigen/svd solvers,
I don't know how the two methods compare in practice. SVD algorithms generally use iterated QR reductions in their implementation, so QR reductions can't be worse numerically. But the SVD probably provides a better metric for rank determination. A google search turns up some literature on the subject that I can't access from home.
OK, here's a good reference <http://www.math.sjsu.edu/%7Efoster/rank/rank_revealing_s.pdf>. A quick look seems to indicate that the SVD is the way to go.
AFAIK, SVD is indeed the way to go, but do we really need this for the orth function ? I am wrong to think that orthonormalizing a matrix of linearly independent vectors is the most common usage for orth ? The difference in terms of speed is really significant (for example, svd of a 2000x100 matrix takes ~1.9 second vs 0.1 s for QR).
This looks a bit like sorting: quicksort is almost always fastest, but no quarantee. The other methods are safer but slower. Maybe the way to go is use a keyword to choose between methods. Chuck
Charles R Harris wrote:
Maybe the way to go is use a keyword to choose between methods.
'safe' vs 'fast' ? :) I wonder whether it would be ok to change the default to QR, though. I will look at the pivoted QR thing, cheers, David
On Fri, Feb 5, 2010 at 2:45 AM, David Cournapeau <david@silveregg.co.jp>wrote:
Charles R Harris wrote:
Maybe the way to go is use a keyword to choose between methods.
'safe' vs 'fast' ? :)
I wonder whether it would be ok to change the default to QR, though. I will look at the pivoted QR thing,
That sounds OK. The documentation should probably mention that the orthogonal columns produced by the two methods will be different, otherwise we will probably get some bug reports. Chuck
On Fri, Feb 5, 2010 at 9:53 AM, Charles R Harris <charlesr.harris@gmail.com>wrote:
On Fri, Feb 5, 2010 at 2:45 AM, David Cournapeau <david@silveregg.co.jp>wrote:
Charles R Harris wrote:
Maybe the way to go is use a keyword to choose between methods.
'safe' vs 'fast' ? :)
I wonder whether it would be ok to change the default to QR, though. I will look at the pivoted QR thing,
That sounds OK. The documentation should probably mention that the orthogonal columns produced by the two methods will be different, otherwise we will probably get some bug reports.
I don't know if you read the paper, but the fastest and safest QR factorizations seem to be the DGEQPX and DGEQPY algorithms in ACM Algorithm 782. It looks like one might need to get permission from the ACM to use the algorithms in a BSD licensed library. Submittal of an algorithm for publication in one of the ACM Transactions impliesthat unrestricted use of the algorithm within a computer is permissible.General permission to copy and distribute the algorithm without fee is granted provided that the copies are not made or distributed for direct commercial advantage. The ACM copyright notice and the title of the publication and its date appear, and notice is given that copying is by permission of the Association for Computing Machinery. To copy otherwise, or to republish, requires a fee and/or specific permission. Chuck
On Fri, Feb 5, 2010 at 1:12 AM, Charles R Harris <charlesr.harris@gmail.com>wrote:
On Fri, Feb 5, 2010 at 12:47 AM, Charles R Harris < charlesr.harris@gmail.com> wrote:
On Fri, Feb 5, 2010 at 12:18 AM, David Cournapeau <david@silveregg.co.jp>wrote:
Charles R Harris wrote:
On Thu, Feb 4, 2010 at 8:45 PM, David Cournapeau <
<mailto:david@silveregg.co.jp>> wrote:
Hi,
I wanted to know if there was a rationale for using svd to orthonormalize the columns of a matrix (in scipy.linalg). QR-based methods are likely to be much faster, and I thought this was the standard, numerically-stable method to orthonormalize a basis ? If
david@silveregg.co.jp the
reason is to deal with rank-deficient matrices, maybe we could add
an
option to choose between them ?
QR with column rotation would deal with rank-deficient matrices and routines for that are available in LAPACK <http://netlib.org/lapack/lug/node42.html>. The SVD was probably used because it was available. The diagonal elements of the R matrix can somewhat take the place of the singular values when column rotation is
used.
So would be it ok to use this column-rotated QR in place of svd for every case in orth ? I would have to check that QR with column rotation is still significantly faster than svd, but I would surprised if if were not the case. QR has also the advantage of being implemented in PLASMA already contrary to eigen/svd solvers,
I don't know how the two methods compare in practice. SVD algorithms generally use iterated QR reductions in their implementation, so QR reductions can't be worse numerically. But the SVD probably provides a better metric for rank determination. A google search turns up some literature on the subject that I can't access from home.
OK, here's a good reference<http://www.math.sjsu.edu/%7Efoster/rank/rank_revealing_s.pdf>. A quick look seems to indicate that the SVD is the way to go.
I take that back. The QR algorithms are faster, but the SVD is more robust. In practice the LAPACK QR algorithm with column pivoting works well for most things, but there are even faster versions. Chuck
On Fri, Feb 05, 2010 at 04:18:02PM +0900, David Cournapeau wrote:
So would be it ok to use this column-rotated QR in place of svd for every case in orth ? I would have to check that QR with column rotation is still significantly faster than svd, but I would surprised if if were not the case. QR has also the advantage of being implemented in PLASMA already contrary to eigen/svd solvers,
Out of curiosity, what's PLASMA, in this context? Gaël
Gael Varoquaux wrote:
On Fri, Feb 05, 2010 at 04:18:02PM +0900, David Cournapeau wrote:
So would be it ok to use this column-rotated QR in place of svd for every case in orth ? I would have to check that QR with column rotation is still significantly faster than svd, but I would surprised if if were not the case. QR has also the advantage of being implemented in PLASMA already contrary to eigen/svd solvers,
Out of curiosity, what's PLASMA, in this context?
http://icl.cs.utk.edu/projectsfiles/plasma/html/README.html """ The main purpose of PLASMA is to address the performance shortcomings of the LAPACK and ScaLAPACK libraries on multicore processors and multi-socket systems of multicore processors. PLASMA provides routines to solve dense general systems of linear equations, symmetric positive definite systems of linear equations and linear least squares problems, using LU, Cholesky, QR and LQ factorizations. Real arithmetic and complex arithmetic are supported in both single precision and double precision. """ It is under BSD license, and as a bonus point, may be compiled easily on windows (MS is one of the sponsor of the project). The main drawback I can see is that it requires a serial BLAS, and ATLAS cannot be switched dynamically between serial and parallel (you have to relink). I am hoping to provide a basic set of wrappers for scipy, cheers, David
On Fri, Feb 5, 2010 at 1:33 AM, David Cournapeau <david@silveregg.co.jp>wrote:
Gael Varoquaux wrote:
On Fri, Feb 05, 2010 at 04:18:02PM +0900, David Cournapeau wrote:
So would be it ok to use this column-rotated QR in place of svd for every case in orth ? I would have to check that QR with column rotation is still significantly faster than svd, but I would surprised if if were not the case. QR has also the advantage of being implemented in PLASMA already contrary to eigen/svd solvers,
Out of curiosity, what's PLASMA, in this context?
http://icl.cs.utk.edu/projectsfiles/plasma/html/README.html
""" The main purpose of PLASMA is to address the performance shortcomings of the LAPACK and ScaLAPACK libraries on multicore processors and multi-socket systems of multicore processors. PLASMA provides routines to solve dense general systems of linear equations, symmetric positive definite systems of linear equations and linear least squares problems, using LU, Cholesky, QR and LQ factorizations. Real arithmetic and complex arithmetic are supported in both single precision and double precision. """
It is under BSD license, and as a bonus point, may be compiled easily on windows (MS is one of the sponsor of the project). The main drawback I can see is that it requires a serial BLAS, and ATLAS cannot be switched dynamically between serial and parallel (you have to relink).
I am hoping to provide a basic set of wrappers for scipy,
Looks like the column pivoting QR algorithm isn't there. I'm not sure what orth is supposed to be used for, but if there is no danger of rank deficiency, then the usual QR algorithm should work just fine. Chuck
Charles R Harris <charlesr.harris <at> gmail.com> writes:
On Fri, Feb 5, 2010 at 1:33 AM, David Cournapeau <david <at> silveregg.co.jp>
Gael Varoquaux wrote:
On Fri, Feb 05, 2010 at 04:18:02PM +0900, David Cournapeau wrote:
So would be it ok to use this column-rotated QR in place of svd for every case in orth ? I would have to check that QR with column rotation is still significantly faster than svd, but I would surprised if if were not the case. QR has also the advantage of being implemented in PLASMA already contrary to eigen/svd solvers,
Out of curiosity, what's PLASMA, in this context? http://icl.cs.utk.edu/projectsfiles/plasma/html/README.html It is under BSD license, and as a bonus point, may be compiled easily on windows (MS is one of the sponsor of the project). The main drawback I can see is that it requires a serial BLAS, and ATLAS cannot be switched dynamically between serial and parallel (you have to relink). I am hoping to provide a basic set of wrappers for scipy,
Looks like the column pivoting QR algorithm isn't there. I'm not sure what orth is supposed to be used for, but if there is no danger of rank deficiency,
wrote: then the usual QR algorithm should work just fine.Chuck Indeed, there is no rank-revealing QR in PLASMA (I'm assuming this is what you're after). And there are no immediate plans for it. Piotr PS. I'm part of PLASMA team.
David Cournapeau <david <at> silveregg.co.jp> writes:
Gael Varoquaux wrote:
On Fri, Feb 05, 2010 at 04:18:02PM +0900, David Cournapeau wrote:
So would be it ok to use this column-rotated QR in place of svd for every case in orth ? I would have to check that QR with column rotation is still significantly faster than svd, but I would surprised if if were not the case. QR has also the advantage of being implemented in PLASMA already contrary to eigen/svd solvers,
Out of curiosity, what's PLASMA, in this context? http://icl.cs.utk.edu/projectsfiles/plasma/html/README.html It is under BSD license, and as a bonus point, may be compiled easily on windows (MS is one of the sponsor of the project). The main drawback I can see is that it requires a serial BLAS, and ATLAS cannot be switched dynamically between serial and parallel (you have to relink).
I am hoping to provide a basic set of wrappers for scipy,
We (PLASMA team) are looking into this problem. In principle it should be possible to just use parallel ATLAS and call internal functions of ATLAS to get at the serial interface. We might add this feature to PLASMA some time soon. Piotr
On Sat, Feb 6, 2010 at 1:05 PM, Piotr Luszczek <luszczek@eecs.utk.edu> wrote:
We (PLASMA team) are looking into this problem. In principle it should be possible to just use parallel ATLAS and call internal functions of ATLAS to get at the serial interface. We might add this feature to PLASMA some time soon.
Great, thank you for the information. cheers, David
participants (5)
-
Charles R Harris -
David Cournapeau -
David Cournapeau -
Gael Varoquaux -
Piotr Luszczek