| Re: [eigen] Householder.h: ::min() and operator<= |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
Dear Christoph,
Am 08.09.2016 um 18:32 schrieb Christoph Hertzberg:
On 2016-09-08 18:04, Peter wrote:
But that's asking for trouble. Wasn't the idea of
NumTraits<TpFloat>::epsilon() to provide this information?
NumTraits::epsilon() is too large in this case -- as you pointed out, this is to prevent underflows when calculating beta.
The related bug is here:
http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1017
Interesting.
I think that the default to not rotate at all, could lead to my problems.
However, to be sure I'll write a Jacobi-diagonalization myself. Maybe I messed up my definitions.
I guess it would be slightly better, if we provided something like a NumTraits::smallest(), ourselves.
I agree. The current approach can't work if numerics_limits<CustomType> is not defined.
I looked into <limits> of gcc4.8, and there is a template
template<typename _Tp>
struct numeric_limits : public __numeric_limits_base
{
/** The minimum finite value, or for floating types with
denormalization, the minimum positive normalized value. */
static _GLIBCXX_CONSTEXPR _Tp
min() _GLIBCXX_USE_NOEXCEPT { return _Tp(); }
[...]
which explains the behaviour. If there is no specialization for a custom type, the min() just calls the default constructor.
Quiet a tricky pitfall.
line 80: if(tailSqNorm <= tol && numext::abs2(numext::imag(c0))<=tol)
Not sure exactly what you are doing with the "tuple of arithmetic types". Is it some kind of interval algebra?
Actually not. I'd love to have all the algorithms working for interval types.
But that's probably calling for a complete rewrite. But you are right, for an interval type
one would be in trouble here.
In that specific case operator<= should be interpreted as "could be smaller or equal".
Well an interval [-1,1] wouldn't be small, but still contains a zero.
Would it help for you, if we replaced the condition by this:
!(tailSqNorm > tol) && !(.... > tol)
I don't think so, but I'll do my homework first (checking that with a homebrew Jacobi method everything is fine).
Yes, it is a bit difficult to make our decompositions robust to arbitrary self-defined types.
Sure, I'm even amazed on how far you got.
We gave some opportunity for customization for pivoting LU decompositions, though (i.e., if your use case is somewhat reasonable, we may consider giving some customization opportunity here as well).
o.k., it will take me some time to understand the details. But I'll be back.
Out of curiosity, is there an advantage of using householder for a 2x2
matrix, as Eigen::SelfAdjointEigenSolver does,
instead of just performing a Jacobi rotation?
Usually, you should use SelfAdjointEigenSolver::computeDirect for 2x2 (and 3x3) matrices. We could make runtime checks for n==2 and n==3, but especially the 3x3 method might be a tiny bit less accurate than the Householder based
method.
Mmmh, I tried SelfAdjointEigenSolver::computeDirect even for fixed size 2x2 and 3x3 matrices.
It always calles Householder. I've added printing a debug message to line 79.
Best regards,
Peter