| Re: [eigen] Specializing max_coeff_visitor for some number types |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/eigen Archives
]
On Wed, 20 Aug 2014, Marc Glisse wrote:
On Wed, 20 Aug 2014, Christoph Hertzberg wrote:
On 19.08.2014 15:53, Christoph Hertzberg wrote:
On 17.08.2014 19:34, Marc Glisse wrote:
Any comment on this patch?
http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2014/06/msg00012.html
Hi Marc,
sorry I forgot about this thread, so thanks for the reminder.
As no further objections arose to your patch, I'll commit it soon.
One thing that is missing is that only LU decompositions are considered by
your patch.
Actually, a clean solution should somehow share the pivoting code for all
pivoting decompositions. At least for QR and LU the code is very similar at
the moment. Any thoughts?
The codes are similar, but not the same. QR has even more threshold,
precision, isMuchSmallerThan and stuff. Merging those codes requires a better
global understanding of what they do than I have (finding the line that says
"max" and replacing it with "best" is about as far as I looked at the
pivoting code...).
Introducing the score in QR shouldn't be much harder than in the full LU, as
long as we basically ignore the effect on biggest_pivot, the case of non-zero
epsilon, etc.
Maybe I should mention in the doc that when we specialize scores, non-zero
values of epsilon are not supported and biggest_pivot -type functions may
give unexpected results. It may be possible to define a sensible interaction
between scores and non-0 thresholds, but that requires a lot more thought.
So no, I don't particularly have any idea.
(and I am currently only using LU, somehow every time I tried to use
something else it failed because it required sqrt or some other thing not
available for rationals)
I think adding some documentation somewhere here would be nice:
http://eigen.tuxfamily.org/dox-devel/TopicCustomizingEigen.html#CustomScalarType
I'll write something as soon as the patch is in, so I know what to document.
An example with rationals (minimize the complexity of the pivot) seems easier
than intervals.
I added an example to #CustomScalarType and used the score in
FullPivHouseholderQR. I feel that doing more requires someone who actually
has an idea what the code is supposed to be doing. Refactoring pivots is
almost orthogonal to what I am doing...
Is the attached pair of patches ok?
--
Marc Glisse
# HG changeset patch
# User Marc Glisse <marc.glisse@xxxxxxxx>
# Date 1403387390 -7200
# Sat Jun 21 23:49:50 2014 +0200
# Node ID c24d264ec62355f527ef4f7160bc5c8dedb23c10
# Parent 9c5d2b7f6e5d13ec18766c2f4c3d2887944805f7
New scoring functor to select the pivot.
diff --git a/Eigen/src/Core/functors/UnaryFunctors.h b/Eigen/src/Core/functors/UnaryFunctors.h
--- a/Eigen/src/Core/functors/UnaryFunctors.h
+++ b/Eigen/src/Core/functors/UnaryFunctors.h
@@ -51,16 +51,44 @@ struct functor_traits<scalar_abs_op<Scal
{
enum {
Cost = NumTraits<Scalar>::AddCost,
PacketAccess = packet_traits<Scalar>::HasAbs
};
};
/** \internal
+ * \brief Template functor to compute the score of a scalar, to chose a pivot
+ *
+ * \sa class CwiseUnaryOp
+ */
+template<typename Scalar> struct scalar_score_coeff_op : scalar_abs_op<Scalar>
+{
+ typedef void Score_is_abs;
+};
+template<typename Scalar>
+struct functor_traits<scalar_score_coeff_op<Scalar> > : functor_traits<scalar_abs_op<Scalar> > {};
+
+/* Avoid recomputing abs when we know the score and they are the same. Not a true Eigen functor. */
+template<typename Scalar, typename=void> struct abs_knowing_score
+{
+ EIGEN_EMPTY_STRUCT_CTOR(abs_knowing_score)
+ typedef typename NumTraits<Scalar>::Real result_type;
+ template<typename Score>
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const result_type operator() (const Scalar& a, const Score&) const { using std::abs; return abs(a); }
+};
+template<typename Scalar> struct abs_knowing_score<Scalar, typename scalar_score_coeff_op<Scalar>::Score_is_abs>
+{
+ EIGEN_EMPTY_STRUCT_CTOR(abs_knowing_score)
+ typedef typename NumTraits<Scalar>::Real result_type;
+ template<typename Scal>
+ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const result_type operator() (const Scal&, const result_type& a) const { return a; }
+};
+
+/** \internal
* \brief Template functor to compute the squared absolute value of a scalar
*
* \sa class CwiseUnaryOp, Cwise::abs2
*/
template<typename Scalar> struct scalar_abs2_op {
EIGEN_EMPTY_STRUCT_CTOR(scalar_abs2_op)
typedef typename NumTraits<Scalar>::Real result_type;
EIGEN_DEVICE_FUNC
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -438,37 +438,40 @@ FullPivLU<MatrixType>& FullPivLU<MatrixT
m_maxpivot = RealScalar(0);
for(Index k = 0; k < size; ++k)
{
// First, we need to find the pivot.
// biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
Index row_of_biggest_in_corner, col_of_biggest_in_corner;
- RealScalar biggest_in_corner;
+ typedef internal::scalar_score_coeff_op<Scalar> Scoring;
+ typedef typename Scoring::result_type Score;
+ Score biggest_in_corner;
biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
- .cwiseAbs()
+ .unaryExpr(Scoring())
.maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
col_of_biggest_in_corner += k; // need to add k to them.
- if(biggest_in_corner==RealScalar(0))
+ if(biggest_in_corner==Score(0))
{
// before exiting, make sure to initialize the still uninitialized transpositions
// in a sane state without destroying what we already have.
m_nonzero_pivots = k;
for(Index i = k; i < size; ++i)
{
m_rowsTranspositions.coeffRef(i) = i;
m_colsTranspositions.coeffRef(i) = i;
}
break;
}
- if(biggest_in_corner > m_maxpivot) m_maxpivot = biggest_in_corner;
+ RealScalar abs_pivot = internal::abs_knowing_score<Scalar>()(m_lu(row_of_biggest_in_corner, col_of_biggest_in_corner), biggest_in_corner);
+ if(abs_pivot > m_maxpivot) m_maxpivot = abs_pivot;
// Now that we've found the pivot, we need to apply the row/col swaps to
// bring it to the location (k,k).
m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
if(k != row_of_biggest_in_corner) {
m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
--- a/Eigen/src/LU/PartialPivLU.h
+++ b/Eigen/src/LU/PartialPivLU.h
@@ -233,34 +233,36 @@ struct partial_lu_impl
* vector \a row_transpositions which must have a size equal to the number
* of columns of the matrix \a lu, and an integer \a nb_transpositions
* which returns the actual number of transpositions.
*
* \returns The index of the first pivot which is exactly zero if any, or a negative number otherwise.
*/
static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
{
+ typedef scalar_score_coeff_op<Scalar> Scoring;
+ typedef typename Scoring::result_type Score;
const Index rows = lu.rows();
const Index cols = lu.cols();
const Index size = (std::min)(rows,cols);
nb_transpositions = 0;
Index first_zero_pivot = -1;
for(Index k = 0; k < size; ++k)
{
Index rrows = rows-k-1;
Index rcols = cols-k-1;
Index row_of_biggest_in_col;
- RealScalar biggest_in_corner
- = lu.col(k).tail(rows-k).cwiseAbs().maxCoeff(&row_of_biggest_in_col);
+ Score biggest_in_corner
+ = lu.col(k).tail(rows-k).unaryExpr(Scoring()).maxCoeff(&row_of_biggest_in_col);
row_of_biggest_in_col += k;
row_transpositions[k] = PivIndex(row_of_biggest_in_col);
- if(biggest_in_corner != RealScalar(0))
+ if(biggest_in_corner != Score(0))
{
if(k != row_of_biggest_in_col)
{
lu.row(k).swap(lu.row(row_of_biggest_in_col));
++nb_transpositions;
}
// FIXME shall we introduce a safe quotient expression in cas 1/lu.coeff(k,k)
# HG changeset patch
# User Marc Glisse <marc.glisse@xxxxxxxx>
# Date 1418596029 -3600
# Sun Dec 14 23:27:09 2014 +0100
# Node ID 90ad87c399a50206a8cd170392b9611fd6ddcffa
# Parent c24d264ec62355f527ef4f7160bc5c8dedb23c10
Document the scoring functor and use it in QR.
diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h
--- a/Eigen/src/QR/FullPivHouseholderQR.h
+++ b/Eigen/src/QR/FullPivHouseholderQR.h
@@ -426,23 +426,25 @@ FullPivHouseholderQR<MatrixType>& FullPi
RealScalar biggest(0);
m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
m_maxpivot = RealScalar(0);
for (Index k = 0; k < size; ++k)
{
Index row_of_biggest_in_corner, col_of_biggest_in_corner;
- RealScalar biggest_in_corner;
+ typedef internal::scalar_score_coeff_op<Scalar> Scoring;
+ typedef typename Scoring::result_type Score;
- biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k)
- .cwiseAbs()
- .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
+ Score score = m_qr.bottomRightCorner(rows-k, cols-k)
+ .unaryExpr(Scoring())
+ .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
row_of_biggest_in_corner += k;
col_of_biggest_in_corner += k;
+ RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score);
if(k==0) biggest = biggest_in_corner;
// if the corner is negligible, then we have less than full rank, and we can finish early
if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
{
m_nonzero_pivots = k;
for(Index i = k; i < size; i++)
{
diff --git a/doc/CustomizingEigen.dox b/doc/CustomizingEigen.dox
--- a/doc/CustomizingEigen.dox
+++ b/doc/CustomizingEigen.dox
@@ -152,14 +152,75 @@ inline adouble imag(const adouble&) {
inline adouble abs(const adouble& x) { return fabs(x); }
inline adouble abs2(const adouble& x) { return x*x; }
}
#endif // ADOLCSUPPORT_H
\endcode
+This other example adds support for the \c mpq_class type from <a href="https://gmplib.org/";>GMP</a>. It shows in particular how to change the way Eigen picks the best pivot during LU factorization. It selects the coefficient with the highest score, where the score is by default the absolute value of a number, but we can define a different score, for instance to prefer pivots with a more compact representation (this is an example, not a recommendation). Note that the scores should always be non-negative and only zero is allowed to have a score of zero. Also, this can interact badly with thresholds for inexact scalar types.
+
+\code
+#include <gmpxx.h>
+#include <Eigen/Core>
+#include <boost/operators.hpp>
+
+namespace Eigen {
+ template<class> struct NumTraits;
+ template<> struct NumTraits<mpq_class>
+ {
+ typedef mpq_class Real;
+ typedef mpq_class NonInteger;
+ typedef mpq_class Nested;
+
+ static inline Real epsilon() { return 0; }
+ static inline Real dummy_precision() { return 0; }
+
+ enum {
+ IsInteger = 0,
+ IsSigned = 1,
+ IsComplex = 0,
+ RequireInitialization = 1,
+ ReadCost = 6,
+ AddCost = 150,
+ MulCost = 100
+ };
+ };
+
+ namespace internal {
+ template<>
+ struct significant_decimals_impl<mpq_class>
+ {
+ // Infinite precision when printing
+ static inline int run() { return 0; }
+ };
+
+ template<> struct scalar_score_coeff_op<mpq_class> {
+ struct result_type : boost::totally_ordered1<result_type> {
+ std::size_t len;
+ result_type(int i = 0) : len(i) {} // Eigen uses Score(0) and Score()
+ result_type(mpq_class const& q) :
+ len(mpz_size(q.get_num_mpz_t())+
+ mpz_size(q.get_den_mpz_t())-1) {}
+ friend bool operator<(result_type x, result_type y) {
+ // 0 is the worst possible pivot
+ if (x.len == 0) return y.len > 0;
+ if (y.len == 0) return false;
+ // Prefer a pivot with a small representation
+ return x.len > y.len;
+ }
+ friend bool operator==(result_type x, result_type y) {
+ // Only used to test if the score is 0
+ return x.len == y.len;
+ }
+ };
+ result_type operator()(mpq_class const& x) const { return x; }
+ };
+ }
+}
+\endcode
\sa \ref TopicPreprocessorDirectives
*/
}