Re: [eigen] makeHouseholder()

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

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