From 58e086b8c81a72daf92c7b4d379713ea9d9fcab3 Mon Sep 17 00:00:00 2001 From: Kolja Brix Date: Mon, 23 Aug 2021 16:00:05 +0000 Subject: [PATCH] Add random matrix generation via SVD --- test/CMakeLists.txt | 1 + test/main.h | 172 +++++++++++++++++++++++++++++++++++++++-- test/random_matrix.cpp | 136 ++++++++++++++++++++++++++++++++ 3 files changed, 301 insertions(+), 8 deletions(-) create mode 100644 test/random_matrix.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5136f82aa..2d6af34ae 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -285,6 +285,7 @@ ei_add_test(array_of_string) ei_add_test(num_dimensions) ei_add_test(stl_iterators) ei_add_test(blasutil) +ei_add_test(random_matrix) if(EIGEN_TEST_CXX11) ei_add_test(initializer_list_construction) ei_add_test(diagonal_matrix_variadic_ctor) diff --git a/test/main.h b/test/main.h index 786673dea..01667388b 100644 --- a/test/main.h +++ b/test/main.h @@ -357,7 +357,7 @@ namespace Eigen #endif // EIGEN_NO_ASSERTION_CHECKING #define EIGEN_INTERNAL_DEBUGGING -#include // required for createRandomPIMatrixOfRank +#include // required for createRandomPIMatrixOfRank and generateRandomMatrixSvs inline void verify_impl(bool condition, const char *testname, const char *file, int line, const char *condition_as_string) { @@ -403,6 +403,36 @@ inline void verify_impl(bool condition, const char *testname, const char *file, } while (0) + namespace Eigen { + +// Forward declarations to avoid ICC warnings +template +bool test_is_equal(const T& actual, const U& expected, bool expect_equal=true); + +template +void createRandomPIMatrixOfRank(Index desired_rank, Index rows, Index cols, MatrixType& m); + +template +void randomPermutationVector(PermutationVectorType& v, Index size); + +template +MatrixType generateRandomUnitaryMatrix(const Index dim); + +template +void generateRandomMatrixSvs(const RealScalarVectorType &svs, const Index rows, const Index cols, MatrixType& M); + +template +VectorType setupRandomSvs(const Index dim, const RealScalar max); + +template +VectorType setupRangeSvs(const Index dim, const RealScalar min, const RealScalar max); + +} // end namespace Eigen + +// Forward declaration to avoid ICC warnings +template std::string type_name(); + + namespace Eigen { template @@ -625,10 +655,6 @@ inline bool test_isUnitary(const MatrixBase& m) return m.isUnitary(test_precision::Scalar>()); } -// Forward declaration to avoid ICC warning -template -bool test_is_equal(const T& actual, const U& expected, bool expect_equal=true); - template bool test_is_equal(const T& actual, const U& expected, bool expect_equal) { @@ -683,7 +709,7 @@ void createRandomPIMatrixOfRank(Index desired_rank, Index rows, Index cols, Matr MatrixType d = MatrixType::Identity(rows,cols); MatrixBType b = MatrixBType::Random(cols,cols); - // set the diagonal such that only desired_rank non-zero entries reamain + // set the diagonal such that only desired_rank non-zero entries remain const Index diag_size = (std::min)(d.rows(),d.cols()); if(diag_size != desired_rank) d.diagonal().segment(desired_rank, diag_size-desired_rank) = VectorType::Zero(diag_size-desired_rank); @@ -719,6 +745,138 @@ void randomPermutationVector(PermutationVectorType& v, Index size) } } +/** + * Generate a random unitary matrix of prescribed dimension. + * + * The algorithm is using a random Householder sequence to produce + * a random unitary matrix. + * + * @tparam MatrixType type of matrix to generate + * @param dim row and column dimension of the requested square matrix + * @return random unitary matrix + */ +template +MatrixType generateRandomUnitaryMatrix(const Index dim) +{ + typedef typename internal::traits::Scalar Scalar; + typedef Matrix VectorType; + + MatrixType v = MatrixType::Identity(dim, dim); + VectorType h = VectorType::Zero(dim); + for (Index i = 0; i < dim; ++i) + { + v.col(i).tail(dim - i - 1) = VectorType::Random(dim - i - 1); + h(i) = 2 / v.col(i).tail(dim - i).squaredNorm(); + } + + const Eigen::HouseholderSequence HSeq(v, h); + return MatrixType(HSeq); +} + +/** + * Generation of random matrix with prescribed singular values. + * + * We generate random matrices with given singular values by setting up + * a singular value decomposition. By choosing the number of zeros as + * singular values we can specify the rank of the matrix. + * Moreover, we also control its spectral norm, which is the largest + * singular value, as well as its condition number with respect to the + * l2-norm, which is the quotient of the largest and smallest singular + * value. + * + * Reference: For details on the method see e.g. Section 8.1 (pp. 62 f) in + * + * C. C. Paige, M. A. Saunders, + * LSQR: An algorithm for sparse linear equations and sparse least squares. + * ACM Transactions on Mathematical Software 8(1), pp. 43-71, 1982. + * https://web.stanford.edu/group/SOL/software/lsqr/lsqr-toms82a.pdf + * + * and also the LSQR webpage https://web.stanford.edu/group/SOL/software/lsqr/. + * + * @tparam MatrixType matrix type to generate + * @tparam RealScalarVectorType vector type with real entries used for singular values + * @param svs vector of desired singular values + * @param rows row dimension of requested random matrix + * @param cols column dimension of requested random matrix + * @param M generated matrix with prescribed singular values + */ +template +void generateRandomMatrixSvs(const RealScalarVectorType &svs, const Index rows, const Index cols, MatrixType& M) +{ + enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime }; + typedef typename internal::traits::Scalar Scalar; + typedef Matrix MatrixAType; + typedef Matrix MatrixBType; + + const Index min_dim = (std::min)(rows, cols); + + const MatrixAType U = generateRandomUnitaryMatrix(rows); + const MatrixBType V = generateRandomUnitaryMatrix(cols); + + M = U.block(0, 0, rows, min_dim) * svs.asDiagonal() * V.block(0, 0, cols, min_dim).transpose(); +} + +/** + * Setup a vector of random singular values with prescribed upper limit. + * For use with generateRandomMatrixSvs(). + * + * Singular values are non-negative real values. By convention (to be consistent with + * singular value decomposition) we sort them in decreasing order. + * + * This strategy produces random singular values in the range [0, max], in particular + * the singular values can be zero or arbitrarily close to zero. + * + * @tparam VectorType vector type with real entries used for singular values + * @tparam RealScalar data type used for real entry + * @param dim number of singular values to generate + * @param max upper bound for singular values + * @return vector of singular values + */ +template +VectorType setupRandomSvs(const Index dim, const RealScalar max) +{ + VectorType svs = max / RealScalar(2) * (VectorType::Random(dim) + VectorType::Ones(dim)); + std::sort(svs.begin(), svs.end(), std::greater()); + return svs; +} + +/** + * Setup a vector of random singular values with prescribed range. + * For use with generateRandomMatrixSvs(). + * + * Singular values are non-negative real values. By convention (to be consistent with + * singular value decomposition) we sort them in decreasing order. + * + * For dim > 1 this strategy generates a vector with largest entry max, smallest entry + * min, and remaining entries in the range [min, max]. For dim == 1 the only entry is + * min. + * + * @tparam VectorType vector type with real entries used for singular values + * @tparam RealScalar data type used for real entry + * @param dim number of singular values to generate + * @param min smallest singular value to use + * @param max largest singular value to use + * @return vector of singular values + */ +template +VectorType setupRangeSvs(const Index dim, const RealScalar min, const RealScalar max) +{ + VectorType svs = VectorType::Random(dim); + if(dim == 0) + return svs; + if(dim == 1) + { + svs(0) = min; + return svs; + } + std::sort(svs.begin(), svs.end(), std::greater()); + + // scale to range [min, max] + const RealScalar c_min = svs(dim - 1), c_max = svs(0); + svs = (svs - VectorType::Constant(dim, c_min)) / (c_max - c_min); + return min * (VectorType::Ones(dim) - svs) + max * svs; +} + /** * Check if number is "not a number" (NaN). * @@ -764,8 +922,6 @@ template<> struct GetDifferentType { typedef float type; }; template struct GetDifferentType > { typedef std::complex::type> type; }; -// Forward declaration to avoid ICC warning -template std::string type_name(); template std::string type_name() { return "other"; } template<> std::string type_name() { return "float"; } template<> std::string type_name() { return "double"; } diff --git a/test/random_matrix.cpp b/test/random_matrix.cpp new file mode 100644 index 000000000..fb877dea4 --- /dev/null +++ b/test/random_matrix.cpp @@ -0,0 +1,136 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2021 Kolja Brix +// +// 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 + + +template +void check_generateRandomUnitaryMatrix(const Index dim) +{ + const MatrixType Q = generateRandomUnitaryMatrix(dim); + + // validate dimensions + VERIFY_IS_EQUAL(Q.rows(), dim); + VERIFY_IS_EQUAL(Q.cols(), dim); + + VERIFY_IS_UNITARY(Q); +} + +template +void check_setupRandomSvs(const Index dim, const RealScalarType max) +{ + const VectorType v = setupRandomSvs(dim, max); + + // validate dimensions + VERIFY_IS_EQUAL(v.size(), dim); + + // check entries + for(Index i = 0; i < v.size(); ++i) + VERIFY_GE(v(i), 0); + for(Index i = 0; i < v.size()-1; ++i) + VERIFY_GE(v(i), v(i+1)); +} + +template +void check_setupRangeSvs(const Index dim, const RealScalarType min, const RealScalarType max) +{ + const VectorType v = setupRangeSvs(dim, min, max); + + // validate dimensions + VERIFY_IS_EQUAL(v.size(), dim); + + // check entries + if(dim == 1) { + VERIFY_IS_APPROX(v(0), min); + } else { + VERIFY_IS_APPROX(v(0), max); + VERIFY_IS_APPROX(v(dim-1), min); + } + for(Index i = 0; i < v.size()-1; ++i) + VERIFY_GE(v(i), v(i+1)); +} + +template +void check_generateRandomMatrixSvs(const Index rows, const Index cols, const Index diag_size, + const RealScalar min_svs, const RealScalar max_svs) +{ + RealVectorType svs = setupRangeSvs(diag_size, min_svs, max_svs); + + MatrixType M; + generateRandomMatrixSvs(svs, rows, cols, M); + + // validate dimensions + VERIFY_IS_EQUAL(M.rows(), rows); + VERIFY_IS_EQUAL(M.cols(), cols); + VERIFY_IS_EQUAL(svs.size(), diag_size); + + // validate singular values + Eigen::JacobiSVD SVD(M); + VERIFY_IS_APPROX(svs, SVD.singularValues()); +} + +template +void check_random_matrix(const MatrixType &m) +{ + enum { + Rows = MatrixType::RowsAtCompileTime, + Cols = MatrixType::ColsAtCompileTime, + DiagSize = EIGEN_SIZE_MIN_PREFER_DYNAMIC(Rows, Cols) + }; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + typedef Matrix RealVectorType; + + const Index rows = m.rows(), cols = m.cols(); + const Index diag_size = (std::min)(rows, cols); + const RealScalar min_svs = 1.0, max_svs = 1000.0; + + // check generation of unitary random matrices + typedef Matrix MatrixAType; + typedef Matrix MatrixBType; + check_generateRandomUnitaryMatrix(rows); + check_generateRandomUnitaryMatrix(cols); + + // test generators for singular values + check_setupRandomSvs(diag_size, max_svs); + check_setupRangeSvs(diag_size, min_svs, max_svs); + + // check generation of random matrices + check_generateRandomMatrixSvs(rows, cols, diag_size, min_svs, max_svs); +} + +EIGEN_DECLARE_TEST(random_matrix) +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1(check_random_matrix(Matrix())); + CALL_SUBTEST_2(check_random_matrix(Matrix())); + CALL_SUBTEST_3(check_random_matrix(Matrix())); + CALL_SUBTEST_4(check_random_matrix(Matrix())); + + CALL_SUBTEST_5(check_random_matrix(Matrix())); + CALL_SUBTEST_6(check_random_matrix(Matrix())); + CALL_SUBTEST_7(check_random_matrix(Matrix())); + CALL_SUBTEST_8(check_random_matrix(Matrix())); + + CALL_SUBTEST_9(check_random_matrix(Matrix, 12, 12>())); + CALL_SUBTEST_10(check_random_matrix(Matrix, 7, 14>())); + CALL_SUBTEST_11(check_random_matrix(Matrix, 15, 11>())); + CALL_SUBTEST_12(check_random_matrix(Matrix, 6, 9>())); + + CALL_SUBTEST_13(check_random_matrix( + MatrixXf(internal::random(1, EIGEN_TEST_MAX_SIZE), internal::random(1, EIGEN_TEST_MAX_SIZE)))); + CALL_SUBTEST_14(check_random_matrix( + MatrixXd(internal::random(1, EIGEN_TEST_MAX_SIZE), internal::random(1, EIGEN_TEST_MAX_SIZE)))); + CALL_SUBTEST_15(check_random_matrix( + MatrixXcf(internal::random(1, EIGEN_TEST_MAX_SIZE), internal::random(1, EIGEN_TEST_MAX_SIZE)))); + CALL_SUBTEST_16(check_random_matrix( + MatrixXcd(internal::random(1, EIGEN_TEST_MAX_SIZE), internal::random(1, EIGEN_TEST_MAX_SIZE)))); + } +}