diff options
author | Benoit Jacob <jacob.benoit.1@gmail.com> | 2009-12-09 12:43:25 -0500 |
---|---|---|
committer | Benoit Jacob <jacob.benoit.1@gmail.com> | 2009-12-09 12:43:25 -0500 |
commit | d2e44f263631981d9e547caafe36b1de5ba785f9 (patch) | |
tree | 574d7aff0554739bea2cede8dc706cd5cd8c7c4a | |
parent | f0315295e9ae2fd8afdc05d3e5b790b4660ffc58 (diff) |
* 4x4 inverse: revert to cofactors method
* inverse tests: use createRandomMatrixOfRank, use more strict precision
* tests: createRandomMatrixOfRank: support 1x1 matrices
* determinant: nest the xpr
* Minor: add comment
-rw-r--r-- | Eigen/src/Core/Minor.h | 5 | ||||
-rw-r--r-- | Eigen/src/LU/Determinant.h | 4 | ||||
-rw-r--r-- | Eigen/src/LU/Inverse.h | 101 | ||||
-rw-r--r-- | test/inverse.cpp | 11 | ||||
-rw-r--r-- | test/main.h | 15 | ||||
-rw-r--r-- | test/prec_inverse_4x4.cpp | 28 |
6 files changed, 43 insertions, 121 deletions
diff --git a/Eigen/src/Core/Minor.h b/Eigen/src/Core/Minor.h index ab058b187..eb44f3760 100644 --- a/Eigen/src/Core/Minor.h +++ b/Eigen/src/Core/Minor.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com> // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -54,7 +54,8 @@ struct ei_traits<Minor<MatrixType> > MaxColsAtCompileTime = (MatrixType::MaxColsAtCompileTime != Dynamic) ? int(MatrixType::MaxColsAtCompileTime) - 1 : Dynamic, Flags = _MatrixTypeNested::Flags & HereditaryBits, - CoeffReadCost = _MatrixTypeNested::CoeffReadCost + CoeffReadCost = _MatrixTypeNested::CoeffReadCost // minor is used typically on tiny matrices, + // where loops are unrolled and the 'if' evaluates at compile time }; }; diff --git a/Eigen/src/LU/Determinant.h b/Eigen/src/LU/Determinant.h index 8870d9f20..27ad6abe9 100644 --- a/Eigen/src/LU/Determinant.h +++ b/Eigen/src/LU/Determinant.h @@ -118,7 +118,9 @@ template<typename Derived> inline typename ei_traits<Derived>::Scalar MatrixBase<Derived>::determinant() const { assert(rows() == cols()); - return ei_determinant_impl<Derived>::run(derived()); + typedef typename ei_nested<Derived,RowsAtCompileTime>::type Nested; + Nested nested(derived()); + return ei_determinant_impl<typename ei_cleantype<Nested>::type>::run(nested); } #endif // EIGEN_DETERMINANT_H diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h index 9d5e86845..8afbfda96 100644 --- a/Eigen/src/LU/Inverse.h +++ b/Eigen/src/LU/Inverse.h @@ -183,92 +183,27 @@ struct ei_compute_inverse_and_det_with_check<MatrixType, ResultType, 3> ****************************/ template<typename MatrixType, typename ResultType> -void ei_compute_inverse_size4_helper(const MatrixType& matrix, ResultType& result) -{ - /* Let's split M into four 2x2 blocks: - * (P Q) - * (R S) - * If P is invertible, with inverse denoted by P_inverse, and if - * (S - R*P_inverse*Q) is also invertible, then the inverse of M is - * (P' Q') - * (R' S') - * where - * S' = (S - R*P_inverse*Q)^(-1) - * P' = P1 + (P1*Q) * S' *(R*P_inverse) - * Q' = -(P_inverse*Q) * S' - * R' = -S' * (R*P_inverse) - */ - typedef Block<ResultType,2,2> XprBlock22; - typedef typename MatrixBase<XprBlock22>::PlainMatrixType Block22; - Block22 P_inverse; - ei_compute_inverse<XprBlock22, Block22>::run(matrix.template block<2,2>(0,0), P_inverse); - const Block22 Q = matrix.template block<2,2>(0,2); - const Block22 P_inverse_times_Q = P_inverse * Q; - const XprBlock22 R = matrix.template block<2,2>(2,0); - const Block22 R_times_P_inverse = R * P_inverse; - const Block22 R_times_P_inverse_times_Q = R_times_P_inverse * Q; - const XprBlock22 S = matrix.template block<2,2>(2,2); - const Block22 X = S - R_times_P_inverse_times_Q; - Block22 Y; - ei_compute_inverse<Block22, Block22>::run(X, Y); - result.template block<2,2>(2,2) = Y; - result.template block<2,2>(2,0) = - Y * R_times_P_inverse; - const Block22 Z = P_inverse_times_Q * Y; - result.template block<2,2>(0,2) = - Z; - result.template block<2,2>(0,0) = P_inverse + Z * R_times_P_inverse; -} - -template<typename MatrixType, typename ResultType> struct ei_compute_inverse<MatrixType, ResultType, 4> { - static inline void run(const MatrixType& _matrix, ResultType& result) + static inline void run(const MatrixType& matrix, ResultType& result) { - typedef typename ResultType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - - // we will do row permutations on the matrix. This copy should have negligible cost. - // if not, consider working in-place on the matrix (const-cast it, but then undo the permutations - // to nevertheless honor constness) - typename MatrixType::PlainMatrixType matrix(_matrix); - - // let's extract from the 2 first colums a 2x2 block whose determinant is as big as possible. - int good_row0, good_row1, good_i; - Matrix<RealScalar,6,1> absdet; - - // any 2x2 block with determinant above this threshold will be considered good enough. - // The magic value 1e-1 here comes from experimentation. The bigger it is, the higher the precision, - // the slower the computation. This value 1e-1 gives precision almost as good as the brutal cofactors - // algorithm, both in average and in worst-case precision. - RealScalar d = (matrix.col(0).squaredNorm()+matrix.col(1).squaredNorm()) * RealScalar(1e-1); - #define ei_inv_size4_helper_macro(i,row0,row1) \ - absdet[i] = ei_abs(matrix.coeff(row0,0)*matrix.coeff(row1,1) \ - - matrix.coeff(row0,1)*matrix.coeff(row1,0)); \ - if(absdet[i] > d) { good_row0=row0; good_row1=row1; goto good; } - ei_inv_size4_helper_macro(0,0,1) - ei_inv_size4_helper_macro(1,0,2) - ei_inv_size4_helper_macro(2,0,3) - ei_inv_size4_helper_macro(3,1,2) - ei_inv_size4_helper_macro(4,1,3) - ei_inv_size4_helper_macro(5,2,3) - - // no 2x2 block has determinant bigger than the threshold. So just take the one that - // has the biggest determinant - absdet.maxCoeff(&good_i); - good_row0 = good_i <= 2 ? 0 : good_i <= 4 ? 1 : 2; - good_row1 = good_i <= 2 ? good_i+1 : good_i <= 4 ? good_i-1 : 3; - - // now good_row0 and good_row1 are correctly set - good: - - // do row permutations to move this 2x2 block to the top - matrix.row(0).swap(matrix.row(good_row0)); - matrix.row(1).swap(matrix.row(good_row1)); - // now applying our helper function is numerically stable - ei_compute_inverse_size4_helper(matrix, result); - // Since we did row permutations on the original matrix, we need to do column permutations - // in the reverse order on the inverse - result.col(1).swap(result.col(good_row1)); - result.col(0).swap(result.col(good_row0)); + result.coeffRef(0,0) = matrix.minor(0,0).determinant(); + result.coeffRef(1,0) = -matrix.minor(0,1).determinant(); + result.coeffRef(2,0) = matrix.minor(0,2).determinant(); + result.coeffRef(3,0) = -matrix.minor(0,3).determinant(); + result.coeffRef(0,2) = matrix.minor(2,0).determinant(); + result.coeffRef(1,2) = -matrix.minor(2,1).determinant(); + result.coeffRef(2,2) = matrix.minor(2,2).determinant(); + result.coeffRef(3,2) = -matrix.minor(2,3).determinant(); + result.coeffRef(0,1) = -matrix.minor(1,0).determinant(); + result.coeffRef(1,1) = matrix.minor(1,1).determinant(); + result.coeffRef(2,1) = -matrix.minor(1,2).determinant(); + result.coeffRef(3,1) = matrix.minor(1,3).determinant(); + result.coeffRef(0,3) = -matrix.minor(3,0).determinant(); + result.coeffRef(1,3) = matrix.minor(3,1).determinant(); + result.coeffRef(2,3) = -matrix.minor(3,2).determinant(); + result.coeffRef(3,3) = matrix.minor(3,3).determinant(); + result /= (matrix.col(0).cwise()*result.row(0).transpose()).sum(); } }; diff --git a/test/inverse.cpp b/test/inverse.cpp index 59b791507..b80e139e0 100644 --- a/test/inverse.cpp +++ b/test/inverse.cpp @@ -38,18 +38,11 @@ template<typename MatrixType> void inverse(const MatrixType& m) typedef typename NumTraits<Scalar>::Real RealScalar; typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType; - MatrixType m1 = MatrixType::Random(rows, cols), + MatrixType m1(rows, cols), m2(rows, cols), mzero = MatrixType::Zero(rows, cols), identity = MatrixType::Identity(rows, rows); - - if (ei_is_same_type<RealScalar,float>::ret) - { - // let's build a more stable to inverse matrix - MatrixType a = MatrixType::Random(rows,cols); - m1 += m1 * m1.adjoint() + a * a.adjoint(); - } - + createRandomMatrixOfRank(rows,rows,rows,m1); m2 = m1.inverse(); VERIFY_IS_APPROX(m1, m2.inverse() ); diff --git a/test/main.h b/test/main.h index 8bc71b3ae..06c17a9ae 100644 --- a/test/main.h +++ b/test/main.h @@ -353,13 +353,26 @@ void createRandomMatrixOfRank(int desired_rank, int rows, int cols, MatrixType& typedef Matrix<Scalar, Rows, Rows> MatrixAType; typedef Matrix<Scalar, Cols, Cols> MatrixBType; + if(desired_rank == 0) + { + m.setZero(rows,cols); + return; + } + + if(desired_rank == 1) + { + m = VectorType::Random(rows) * VectorType::Random(cols).transpose(); + return; + } + MatrixAType a = MatrixAType::Random(rows,rows); MatrixType d = MatrixType::Identity(rows,cols); MatrixBType b = MatrixBType::Random(cols,cols); // set the diagonal such that only desired_rank non-zero entries reamain const int diag_size = std::min(d.rows(),d.cols()); - d.diagonal().segment(desired_rank, diag_size-desired_rank) = VectorType::Zero(diag_size-desired_rank); + if(diag_size != desired_rank) + d.diagonal().segment(desired_rank, diag_size-desired_rank) = VectorType::Zero(diag_size-desired_rank); HouseholderQR<MatrixAType> qra(a); HouseholderQR<MatrixBType> qrb(b); diff --git a/test/prec_inverse_4x4.cpp b/test/prec_inverse_4x4.cpp index 613535346..83b1a8a71 100644 --- a/test/prec_inverse_4x4.cpp +++ b/test/prec_inverse_4x4.cpp @@ -26,28 +26,6 @@ #include <Eigen/LU> #include <algorithm> -Matrix4f inverse(const Matrix4f& m) -{ - Matrix4f r; - r(0,0) = m.minor(0,0).determinant(); - r(1,0) = -m.minor(0,1).determinant(); - r(2,0) = m.minor(0,2).determinant(); - r(3,0) = -m.minor(0,3).determinant(); - r(0,2) = m.minor(2,0).determinant(); - r(1,2) = -m.minor(2,1).determinant(); - r(2,2) = m.minor(2,2).determinant(); - r(3,2) = -m.minor(2,3).determinant(); - r(0,1) = -m.minor(1,0).determinant(); - r(1,1) = m.minor(1,1).determinant(); - r(2,1) = -m.minor(1,2).determinant(); - r(3,1) = m.minor(1,3).determinant(); - r(0,3) = -m.minor(3,0).determinant(); - r(1,3) = m.minor(3,1).determinant(); - r(2,3) = -m.minor(3,2).determinant(); - r(3,3) = m.minor(3,3).determinant(); - return r / (m(0,0)*r(0,0) + m(1,0)*r(0,1) + m(2,0)*r(0,2) + m(3,0)*r(0,3)); -} - template<typename MatrixType> void inverse_permutation_4x4() { typedef typename MatrixType::Scalar Scalar; @@ -79,7 +57,7 @@ template<typename MatrixType> void inverse_general_4x4(int repeat) do { m = MatrixType::Random(); absdet = ei_abs(m.determinant()); - } while(absdet < 10 * epsilon<Scalar>()); + } while(absdet < epsilon<Scalar>()); MatrixType inv = m.inverse(); double error = double( (m*inv-MatrixType::Identity()).norm() * absdet / epsilon<Scalar>() ); error_sum += error; @@ -89,8 +67,8 @@ template<typename MatrixType> void inverse_general_4x4(int repeat) double error_avg = error_sum / repeat; EIGEN_DEBUG_VAR(error_avg); EIGEN_DEBUG_VAR(error_max); - VERIFY(error_avg < (NumTraits<Scalar>::IsComplex ? 8.4 : 1.4) ); - VERIFY(error_max < (NumTraits<Scalar>::IsComplex ? 160.0 : 75.) ); + VERIFY(error_avg < (NumTraits<Scalar>::IsComplex ? 8.0 : 1.0)); + VERIFY(error_max < (NumTraits<Scalar>::IsComplex ? 64.0 : 16.0)); } void test_prec_inverse_4x4() |