| // This file is part of Eigen, a lightweight C++ template library |
| // for linear algebra. |
| // |
| // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> |
| // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> |
| // |
| // This Source Code Form is subject to the terms of the Mozilla |
| // Public License v. 2.0. If a copy of the MPL was not distributed |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| |
| #include "main.h" |
| #include <Eigen/QR> |
| #include <Eigen/SVD> |
| #include "solverbase.h" |
| |
| template <typename MatrixType> |
| void cod() { |
| Index rows = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE); |
| Index cols = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE); |
| Index cols2 = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE); |
| Index rank = internal::random<Index>(1, (std::min)(rows, cols) - 1); |
| |
| typedef typename MatrixType::Scalar Scalar; |
| typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType; |
| MatrixType matrix; |
| createRandomPIMatrixOfRank(rank, rows, cols, matrix); |
| CompleteOrthogonalDecomposition<MatrixType> cod(matrix); |
| VERIFY(rank == cod.rank()); |
| VERIFY(cols - cod.rank() == cod.dimensionOfKernel()); |
| VERIFY(!cod.isInjective()); |
| VERIFY(!cod.isInvertible()); |
| VERIFY(!cod.isSurjective()); |
| |
| MatrixQType q = cod.householderQ(); |
| VERIFY_IS_UNITARY(q); |
| |
| MatrixType z = cod.matrixZ(); |
| VERIFY_IS_UNITARY(z); |
| |
| MatrixType t; |
| t.setZero(rows, cols); |
| t.topLeftCorner(rank, rank) = cod.matrixT().topLeftCorner(rank, rank).template triangularView<Upper>(); |
| |
| MatrixType c = q * t * z * cod.colsPermutation().inverse(); |
| VERIFY_IS_APPROX(matrix, c); |
| |
| check_solverbase<MatrixType, MatrixType>(matrix, cod, rows, cols, cols2); |
| |
| // Verify that we get the same minimum-norm solution as the SVD. |
| MatrixType exact_solution = MatrixType::Random(cols, cols2); |
| MatrixType rhs = matrix * exact_solution; |
| MatrixType cod_solution = cod.solve(rhs); |
| JacobiSVD<MatrixType, ComputeThinU | ComputeThinV> svd(matrix); |
| MatrixType svd_solution = svd.solve(rhs); |
| VERIFY_IS_APPROX(cod_solution, svd_solution); |
| |
| MatrixType pinv = cod.pseudoInverse(); |
| VERIFY_IS_APPROX(cod_solution, pinv * rhs); |
| } |
| |
| template <typename MatrixType, int Cols2> |
| void cod_fixedsize() { |
| enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime }; |
| typedef typename MatrixType::Scalar Scalar; |
| typedef CompleteOrthogonalDecomposition<Matrix<Scalar, Rows, Cols> > COD; |
| int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols)) - 1); |
| Matrix<Scalar, Rows, Cols> matrix; |
| createRandomPIMatrixOfRank(rank, Rows, Cols, matrix); |
| COD cod(matrix); |
| VERIFY(rank == cod.rank()); |
| VERIFY(Cols - cod.rank() == cod.dimensionOfKernel()); |
| VERIFY(cod.isInjective() == (rank == Rows)); |
| VERIFY(cod.isSurjective() == (rank == Cols)); |
| VERIFY(cod.isInvertible() == (cod.isInjective() && cod.isSurjective())); |
| |
| check_solverbase<Matrix<Scalar, Cols, Cols2>, Matrix<Scalar, Rows, Cols2> >(matrix, cod, Rows, Cols, Cols2); |
| |
| // Verify that we get the same minimum-norm solution as the SVD. |
| Matrix<Scalar, Cols, Cols2> exact_solution; |
| exact_solution.setRandom(Cols, Cols2); |
| Matrix<Scalar, Rows, Cols2> rhs = matrix * exact_solution; |
| Matrix<Scalar, Cols, Cols2> cod_solution = cod.solve(rhs); |
| JacobiSVD<MatrixType, ComputeFullU | ComputeFullV> svd(matrix); |
| Matrix<Scalar, Cols, Cols2> svd_solution = svd.solve(rhs); |
| VERIFY_IS_APPROX(cod_solution, svd_solution); |
| |
| typename Inverse<COD>::PlainObject pinv = cod.pseudoInverse(); |
| VERIFY_IS_APPROX(cod_solution, pinv * rhs); |
| } |
| |
| template <typename MatrixType> |
| void qr() { |
| using std::sqrt; |
| |
| Index rows = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE), cols = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE), |
| cols2 = internal::random<Index>(2, EIGEN_TEST_MAX_SIZE); |
| Index rank = internal::random<Index>(1, (std::min)(rows, cols) - 1); |
| |
| typedef typename MatrixType::Scalar Scalar; |
| typedef typename MatrixType::RealScalar RealScalar; |
| typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType; |
| MatrixType m1; |
| createRandomPIMatrixOfRank(rank, rows, cols, m1); |
| ColPivHouseholderQR<MatrixType> qr(m1); |
| VERIFY_IS_EQUAL(rank, qr.rank()); |
| VERIFY_IS_EQUAL(cols - qr.rank(), qr.dimensionOfKernel()); |
| VERIFY(!qr.isInjective()); |
| VERIFY(!qr.isInvertible()); |
| VERIFY(!qr.isSurjective()); |
| |
| MatrixQType q = qr.householderQ(); |
| VERIFY_IS_UNITARY(q); |
| |
| MatrixType r = qr.matrixQR().template triangularView<Upper>(); |
| MatrixType c = q * r * qr.colsPermutation().inverse(); |
| VERIFY_IS_APPROX(m1, c); |
| |
| // Verify that the absolute value of the diagonal elements in R are |
| // non-increasing until they reach the singularity threshold. |
| RealScalar threshold = sqrt(RealScalar(rows)) * numext::abs(r(0, 0)) * NumTraits<Scalar>::epsilon(); |
| for (Index i = 0; i < (std::min)(rows, cols) - 1; ++i) { |
| RealScalar x = numext::abs(r(i, i)); |
| RealScalar y = numext::abs(r(i + 1, i + 1)); |
| if (x < threshold && y < threshold) continue; |
| if (!test_isApproxOrLessThan(y, x)) { |
| for (Index j = 0; j < (std::min)(rows, cols); ++j) { |
| std::cout << "i = " << j << ", |r_ii| = " << numext::abs(r(j, j)) << std::endl; |
| } |
| std::cout << "Failure at i=" << i << ", rank=" << rank << ", threshold=" << threshold << std::endl; |
| } |
| VERIFY_IS_APPROX_OR_LESS_THAN(y, x); |
| } |
| |
| check_solverbase<MatrixType, MatrixType>(m1, qr, rows, cols, cols2); |
| |
| { |
| MatrixType m2, m3; |
| Index size = rows; |
| do { |
| m1 = MatrixType::Random(size, size); |
| qr.compute(m1); |
| } while (!qr.isInvertible()); |
| MatrixType m1_inv = qr.inverse(); |
| m3 = m1 * MatrixType::Random(size, cols2); |
| m2 = qr.solve(m3); |
| VERIFY_IS_APPROX(m2, m1_inv * m3); |
| } |
| } |
| |
| template <typename MatrixType, int Cols2> |
| void qr_fixedsize() { |
| using std::abs; |
| using std::sqrt; |
| enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime }; |
| typedef typename MatrixType::Scalar Scalar; |
| typedef typename MatrixType::RealScalar RealScalar; |
| int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols)) - 1); |
| Matrix<Scalar, Rows, Cols> m1; |
| createRandomPIMatrixOfRank(rank, Rows, Cols, m1); |
| ColPivHouseholderQR<Matrix<Scalar, Rows, Cols> > qr(m1); |
| VERIFY_IS_EQUAL(rank, qr.rank()); |
| VERIFY_IS_EQUAL(Cols - qr.rank(), qr.dimensionOfKernel()); |
| VERIFY_IS_EQUAL(qr.isInjective(), (rank == Rows)); |
| VERIFY_IS_EQUAL(qr.isSurjective(), (rank == Cols)); |
| VERIFY_IS_EQUAL(qr.isInvertible(), (qr.isInjective() && qr.isSurjective())); |
| |
| Matrix<Scalar, Rows, Cols> r = qr.matrixQR().template triangularView<Upper>(); |
| Matrix<Scalar, Rows, Cols> c = qr.householderQ() * r * qr.colsPermutation().inverse(); |
| VERIFY_IS_APPROX(m1, c); |
| |
| check_solverbase<Matrix<Scalar, Cols, Cols2>, Matrix<Scalar, Rows, Cols2> >(m1, qr, Rows, Cols, Cols2); |
| |
| // Verify that the absolute value of the diagonal elements in R are |
| // non-increasing until they reache the singularity threshold. |
| RealScalar threshold = sqrt(RealScalar(Rows)) * (std::abs)(r(0, 0)) * NumTraits<Scalar>::epsilon(); |
| for (Index i = 0; i < (std::min)(int(Rows), int(Cols)) - 1; ++i) { |
| RealScalar x = numext::abs(r(i, i)); |
| RealScalar y = numext::abs(r(i + 1, i + 1)); |
| if (x < threshold && y < threshold) continue; |
| if (!test_isApproxOrLessThan(y, x)) { |
| for (Index j = 0; j < (std::min)(int(Rows), int(Cols)); ++j) { |
| std::cout << "i = " << j << ", |r_ii| = " << numext::abs(r(j, j)) << std::endl; |
| } |
| std::cout << "Failure at i=" << i << ", rank=" << rank << ", threshold=" << threshold << std::endl; |
| } |
| VERIFY_IS_APPROX_OR_LESS_THAN(y, x); |
| } |
| } |
| |
| // This test is meant to verify that pivots are chosen such that |
| // even for a graded matrix, the diagonal of R falls of roughly |
| // monotonically until it reaches the threshold for singularity. |
| // We use the so-called Kahan matrix, which is a famous counter-example |
| // for rank-revealing QR. See |
| // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf |
| // page 3 for more detail. |
| template <typename MatrixType> |
| void qr_kahan_matrix() { |
| using std::abs; |
| using std::sqrt; |
| typedef typename MatrixType::Scalar Scalar; |
| typedef typename MatrixType::RealScalar RealScalar; |
| |
| Index rows = 300, cols = rows; |
| |
| MatrixType m1; |
| m1.setZero(rows, cols); |
| RealScalar s = std::pow(NumTraits<RealScalar>::epsilon(), 1.0 / rows); |
| RealScalar c = std::sqrt(1 - s * s); |
| RealScalar pow_s_i(1.0); // pow(s,i) |
| for (Index i = 0; i < rows; ++i) { |
| m1(i, i) = pow_s_i; |
| m1.row(i).tail(rows - i - 1) = -pow_s_i * c * MatrixType::Ones(1, rows - i - 1); |
| pow_s_i *= s; |
| } |
| m1 = (m1 + m1.transpose()).eval(); |
| ColPivHouseholderQR<MatrixType> qr(m1); |
| MatrixType r = qr.matrixQR().template triangularView<Upper>(); |
| |
| RealScalar threshold = std::sqrt(RealScalar(rows)) * numext::abs(r(0, 0)) * NumTraits<Scalar>::epsilon(); |
| for (Index i = 0; i < (std::min)(rows, cols) - 1; ++i) { |
| RealScalar x = numext::abs(r(i, i)); |
| RealScalar y = numext::abs(r(i + 1, i + 1)); |
| if (x < threshold && y < threshold) continue; |
| if (!test_isApproxOrLessThan(y, x)) { |
| for (Index j = 0; j < (std::min)(rows, cols); ++j) { |
| std::cout << "i = " << j << ", |r_ii| = " << numext::abs(r(j, j)) << std::endl; |
| } |
| std::cout << "Failure at i=" << i << ", rank=" << qr.rank() << ", threshold=" << threshold << std::endl; |
| } |
| VERIFY_IS_APPROX_OR_LESS_THAN(y, x); |
| } |
| } |
| |
| template <typename MatrixType> |
| void qr_invertible() { |
| using std::abs; |
| using std::log; |
| typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; |
| typedef typename MatrixType::Scalar Scalar; |
| |
| int size = internal::random<int>(10, 50); |
| |
| MatrixType m1(size, size), m2(size, size), m3(size, size); |
| m1 = MatrixType::Random(size, size); |
| |
| if (internal::is_same<RealScalar, float>::value) { |
| // let's build a matrix more stable to inverse |
| MatrixType a = MatrixType::Random(size, size * 2); |
| m1 += a * a.adjoint(); |
| } |
| |
| ColPivHouseholderQR<MatrixType> qr(m1); |
| |
| check_solverbase<MatrixType, MatrixType>(m1, qr, size, size, size); |
| |
| // now construct a matrix with prescribed determinant |
| m1.setZero(); |
| for (int i = 0; i < size; i++) m1(i, i) = internal::random<Scalar>(); |
| Scalar det = m1.diagonal().prod(); |
| RealScalar absdet = abs(det); |
| m3 = qr.householderQ(); // get a unitary |
| m1 = m3 * m1 * m3.adjoint(); |
| qr.compute(m1); |
| VERIFY_IS_APPROX(det, qr.determinant()); |
| VERIFY_IS_APPROX(absdet, qr.absDeterminant()); |
| VERIFY_IS_APPROX(log(absdet), qr.logAbsDeterminant()); |
| } |
| |
| template <typename MatrixType> |
| void qr_verify_assert() { |
| MatrixType tmp; |
| |
| ColPivHouseholderQR<MatrixType> qr; |
| VERIFY_RAISES_ASSERT(qr.matrixQR()) |
| VERIFY_RAISES_ASSERT(qr.solve(tmp)) |
| VERIFY_RAISES_ASSERT(qr.transpose().solve(tmp)) |
| VERIFY_RAISES_ASSERT(qr.adjoint().solve(tmp)) |
| VERIFY_RAISES_ASSERT(qr.householderQ()) |
| VERIFY_RAISES_ASSERT(qr.dimensionOfKernel()) |
| VERIFY_RAISES_ASSERT(qr.isInjective()) |
| VERIFY_RAISES_ASSERT(qr.isSurjective()) |
| VERIFY_RAISES_ASSERT(qr.isInvertible()) |
| VERIFY_RAISES_ASSERT(qr.inverse()) |
| VERIFY_RAISES_ASSERT(qr.determinant()) |
| VERIFY_RAISES_ASSERT(qr.absDeterminant()) |
| VERIFY_RAISES_ASSERT(qr.logAbsDeterminant()) |
| } |
| |
| template <typename MatrixType> |
| void cod_verify_assert() { |
| MatrixType tmp; |
| |
| CompleteOrthogonalDecomposition<MatrixType> cod; |
| VERIFY_RAISES_ASSERT(cod.matrixQTZ()) |
| VERIFY_RAISES_ASSERT(cod.solve(tmp)) |
| VERIFY_RAISES_ASSERT(cod.transpose().solve(tmp)) |
| VERIFY_RAISES_ASSERT(cod.adjoint().solve(tmp)) |
| VERIFY_RAISES_ASSERT(cod.householderQ()) |
| VERIFY_RAISES_ASSERT(cod.dimensionOfKernel()) |
| VERIFY_RAISES_ASSERT(cod.isInjective()) |
| VERIFY_RAISES_ASSERT(cod.isSurjective()) |
| VERIFY_RAISES_ASSERT(cod.isInvertible()) |
| VERIFY_RAISES_ASSERT(cod.pseudoInverse()) |
| VERIFY_RAISES_ASSERT(cod.determinant()) |
| VERIFY_RAISES_ASSERT(cod.absDeterminant()) |
| VERIFY_RAISES_ASSERT(cod.logAbsDeterminant()) |
| } |
| |
| EIGEN_DECLARE_TEST(qr_colpivoting) { |
| for (int i = 0; i < g_repeat; i++) { |
| CALL_SUBTEST_1(qr<MatrixXf>()); |
| CALL_SUBTEST_2(qr<MatrixXd>()); |
| CALL_SUBTEST_3(qr<MatrixXcd>()); |
| CALL_SUBTEST_4((qr_fixedsize<Matrix<float, 3, 5>, 4>())); |
| CALL_SUBTEST_5((qr_fixedsize<Matrix<double, 6, 2>, 3>())); |
| CALL_SUBTEST_5((qr_fixedsize<Matrix<double, 1, 1>, 1>())); |
| } |
| |
| for (int i = 0; i < g_repeat; i++) { |
| CALL_SUBTEST_1(cod<MatrixXf>()); |
| CALL_SUBTEST_2(cod<MatrixXd>()); |
| CALL_SUBTEST_3(cod<MatrixXcd>()); |
| CALL_SUBTEST_4((cod_fixedsize<Matrix<float, 3, 5>, 4>())); |
| CALL_SUBTEST_5((cod_fixedsize<Matrix<double, 6, 2>, 3>())); |
| CALL_SUBTEST_5((cod_fixedsize<Matrix<double, 1, 1>, 1>())); |
| } |
| |
| for (int i = 0; i < g_repeat; i++) { |
| CALL_SUBTEST_1(qr_invertible<MatrixXf>()); |
| CALL_SUBTEST_2(qr_invertible<MatrixXd>()); |
| CALL_SUBTEST_6(qr_invertible<MatrixXcf>()); |
| CALL_SUBTEST_3(qr_invertible<MatrixXcd>()); |
| } |
| |
| CALL_SUBTEST_7(qr_verify_assert<Matrix3f>()); |
| CALL_SUBTEST_8(qr_verify_assert<Matrix3d>()); |
| CALL_SUBTEST_1(qr_verify_assert<MatrixXf>()); |
| CALL_SUBTEST_2(qr_verify_assert<MatrixXd>()); |
| CALL_SUBTEST_6(qr_verify_assert<MatrixXcf>()); |
| CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>()); |
| |
| CALL_SUBTEST_7(cod_verify_assert<Matrix3f>()); |
| CALL_SUBTEST_8(cod_verify_assert<Matrix3d>()); |
| CALL_SUBTEST_1(cod_verify_assert<MatrixXf>()); |
| CALL_SUBTEST_2(cod_verify_assert<MatrixXd>()); |
| CALL_SUBTEST_6(cod_verify_assert<MatrixXcf>()); |
| CALL_SUBTEST_3(cod_verify_assert<MatrixXcd>()); |
| |
| // Test problem size constructors |
| CALL_SUBTEST_9(ColPivHouseholderQR<MatrixXf>(10, 20)); |
| |
| CALL_SUBTEST_1(qr_kahan_matrix<MatrixXf>()); |
| CALL_SUBTEST_2(qr_kahan_matrix<MatrixXd>()); |
| } |