diff options
author | Rasmus Munk Larsen <rmlarsen@google.com> | 2016-02-06 16:32:00 -0800 |
---|---|---|
committer | Rasmus Munk Larsen <rmlarsen@google.com> | 2016-02-06 16:32:00 -0800 |
commit | d904c8ac8f5963b0c6f2528f4ee1823350277713 (patch) | |
tree | 58a49b2ab896ec432af0bb744b165202f5b8196a /test/qr_colpivoting.cpp | |
parent | 093f2b3c01518697fad88f9cfe91c576cd8e6f6e (diff) |
Implement complete orthogonal decomposition in Eigen.
Diffstat (limited to 'test/qr_colpivoting.cpp')
-rw-r--r-- | test/qr_colpivoting.cpp | 77 |
1 files changed, 77 insertions, 0 deletions
diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp index 9c989823e..648250af6 100644 --- a/test/qr_colpivoting.cpp +++ b/test/qr_colpivoting.cpp @@ -10,6 +10,83 @@ #include "main.h" #include <Eigen/QR> +#include <Eigen/SVD> + +template <typename MatrixType> +void cod() { + typedef typename MatrixType::Index Index; + + 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); + + MatrixType exact_solution = MatrixType::Random(cols, cols2); + MatrixType rhs = matrix * exact_solution; + MatrixType cod_solution = cod.solve(rhs); + VERIFY_IS_APPROX(rhs, matrix * cod_solution); + + // Verify that we get the same minimum-norm solution as the SVD. + JacobiSVD<MatrixType> svd(matrix, ComputeThinU | ComputeThinV); + MatrixType svd_solution = svd.solve(rhs); + VERIFY_IS_APPROX(cod_solution, svd_solution); +} + +template <typename MatrixType, int Cols2> +void cod_fixedsize() { + enum { + Rows = MatrixType::RowsAtCompileTime, + Cols = MatrixType::ColsAtCompileTime + }; + typedef typename MatrixType::Scalar Scalar; + int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols)) - 1); + Matrix<Scalar, Rows, Cols> matrix; + createRandomPIMatrixOfRank(rank, Rows, Cols, matrix); + CompleteOrthogonalDecomposition<Matrix<Scalar, Rows, Cols> > 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())); + + 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); + VERIFY_IS_APPROX(rhs, matrix * cod_solution); + + // Verify that we get the same minimum-norm solution as the SVD. + JacobiSVD<MatrixType> svd(matrix, ComputeFullU | ComputeFullV); + Matrix<Scalar, Cols, Cols2> svd_solution = svd.solve(rhs); + VERIFY_IS_APPROX(cod_solution, svd_solution); +} template<typename MatrixType> void qr() { |