| Re: [eigen] makeHouseholder() |
[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]
I guess the main cases (e.g. for the QR-decomposition use case) where this branch is/was taken are matrices which (partially) are already (almost) triangular or matrices which are near rank-deficient.
At least in the latter case abs2(c0) should be small as well.If we want to avoid applying an actual reflection in case the tail is (near) zero, we should probably also make a special case for tailSqNorm<=tol but numext::abs2(numext::imag(c0))>tol (i.e., also clear the essential part, but do a proper reflection based on c0). I was assuming the special branch was mostly to avoid divisions by near zeros in the main branch.
Btw:The Lapack implementation [1] seems to be using a stable norm implementation [2] and then always uses slapy2 or slapy3 [3] (which are essentially like std::hypot) for further computations of norms. It does not short-circuit for norm(t)^2==0 or similar, only if the actual stable norm is zero. On the other hand, we (always/usually?) normalize our matrices before calling makeHouseholder, so I guess we should only hit zero cases if we are actually calculating with only numerical noise left.
Overall, I guess the best solution would be to find test cases which demonstrate an advantage (or at least a difference) for either implementation.
Christoph
[1] https://github.com/Reference-LAPACK/lapack/blob/master/SRC/clarfgp.f
[2] https://github.com/Reference-LAPACK/lapack/blob/master/BLAS/SRC/scnrm2.f
[3] https://github.com/Reference-LAPACK/lapack/blob/master/SRC/slapy2.f
https://github.com/Reference-LAPACK/lapack/blob/master/SRC/slapy3.f
On 2016-09-17 23:19, Gael Guennebaud wrote:
Sorry for late reply, but I don't see how this change can increase stability. Recall that the goal here is to compute a reflection H such that H*v=[x 0 ... 0] with x real. The initial condition simply check whether v was already in the desired form, and if so we set H=I, and we are done. This identity transformation will also be skipped during later computation. In contrast, the new condition won't produce an identity transformation in this case and will thus lost some bit of precision during the useless computations. gael On Fri, Sep 16, 2016 at 2:26 PM, Christoph Hertzberg < chtz@xxxxxxxxxxxxxxxxxxxxxxxx> wrote:Pushed: https://bitbucket.org/eigen/eigen/commits/cc01b1eb81f23a41 On 2016-09-16 14:09, Christoph Hertzberg wrote:On 2016-09-12 20:03, Peter wrote:Dear All, is there are specific reason for if(tailSqNorm <= tol && numext::abs2(numext::imag(c0))<=tol) or could it be changed into if(tailSqNorm <= tol && numext::abs2(c0)<=tol)As discussed at the beginning of the week on IRC, I second that change. I'll be pushing this change (if I don't get any objections within the next minutes). We should probably try to make some unit tests which would fail without this change (preferably not requiring custom types). Cheers, Christoph within makeHouseholder() in Householder.C, [*] see below?I'm asking as I stumbled over the specific case of c0=1.0+0.0i and a 2x2 matrix, where tailSqNorm==0. The current implementation with the use the special case tau = RealScalar(0); beta = numext::real(c0); essential.setZero(); while if(tailSqNorm <= tol && numext::abs2(c0)<=tol) would select the normal path in that case, beta = sqrt(numext::abs2(c0) + tailSqNorm); if (numext::real(c0)>=RealScalar(0)) beta = -beta; essential = tail / (c0 - beta); tau = conj((beta - c0) / beta); which is fine to evaluate. Due to sign choice of beta, (c0 - beta) will even be further away from zero as beta itself, so numerically that shouldn't be a problem. I'm asking because I could pin-point, with the help of Christoph on IRC, this particular if-statement causing problems for my own custom type. Best regards, Peter [*] template<typename Derived> template<typename EssentialPart> void MatrixBase<Derived>::makeHouseholder( EssentialPart& essential, Scalar& tau, RealScalar& beta) const { using std::sqrt; using numext::conj; EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) VectorBlock<const Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1); RealScalar tailSqNorm = size()==1 ? RealScalar(0) : tail.squaredNorm(); Scalar c0 = coeff(0); const RealScalar tol = (std::numeric_limits<RealScalar>::min)(); if(tailSqNorm <= tol && numext::abs2(numext::imag(c0))<=tol) { tau = RealScalar(0); beta = numext::real(c0); essential.setZero(); } else { beta = sqrt(numext::abs2(c0) + tailSqNorm); if (numext::real(c0)>=RealScalar(0)) beta = -beta; essential = tail / (c0 - beta); tau = conj((beta - c0) / beta); } }-- Dipl. Inf., Dipl. Math. Christoph Hertzberg Universität Bremen FB 3 - Mathematik und Informatik AG Robotik Robert-Hooke-Straße 1 28359 Bremen, Germany Zentrale: +49 421 178 45-6611 Besuchsadresse der Nebengeschäftsstelle: Robert-Hooke-Straße 5 28359 Bremen, Germany Tel.: +49 421 178 45-4021 Empfang: +49 421 178 45-6600 Fax: +49 421 178 45-4150 E-Mail: chtz@xxxxxxxxxxxxxxxxxxxxxxxx Weitere Informationen: http://www.informatik.uni-bremen.de/robotik
-- Dipl. Inf., Dipl. Math. Christoph Hertzberg Universität Bremen FB 3 - Mathematik und Informatik AG Robotik Robert-Hooke-Straße 1 28359 Bremen, Germany Zentrale: +49 421 178 45-6611 Besuchsadresse der Nebengeschäftsstelle: Robert-Hooke-Straße 5 28359 Bremen, Germany Tel.: +49 421 178 45-4021 Empfang: +49 421 178 45-6600 Fax: +49 421 178 45-4150 E-Mail: chtz@xxxxxxxxxxxxxxxxxxxxxxxx Weitere Informationen: http://www.informatik.uni-bremen.de/robotik
| Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |