From f5e2a548e1259aad494ee577a0e3849cd0fb07c8 Mon Sep 17 00:00:00 2001 From: Marcus Edel Date: Sat, 18 Feb 2017 18:34:59 +0100 Subject: [PATCH 1/6] Add implementation of the randomized block krylov svd method. --- src/mlpack/methods/CMakeLists.txt | 1 + .../methods/block_krylov_svd/CMakeLists.txt | 15 ++ .../randomized_block_krylov_svd.cpp | 88 ++++++++++++ .../randomized_block_krylov_svd.hpp | 130 ++++++++++++++++++ 4 files changed, 234 insertions(+) create mode 100644 src/mlpack/methods/block_krylov_svd/CMakeLists.txt create mode 100644 src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp create mode 100644 src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.hpp diff --git a/src/mlpack/methods/CMakeLists.txt b/src/mlpack/methods/CMakeLists.txt index dde69de0b1..2218c378c1 100644 --- a/src/mlpack/methods/CMakeLists.txt +++ b/src/mlpack/methods/CMakeLists.txt @@ -21,6 +21,7 @@ set(DIRS approx_kfn amf ann + block_krylov_svd cf decision_stump det diff --git a/src/mlpack/methods/block_krylov_svd/CMakeLists.txt b/src/mlpack/methods/block_krylov_svd/CMakeLists.txt new file mode 100644 index 0000000000..6380befb28 --- /dev/null +++ b/src/mlpack/methods/block_krylov_svd/CMakeLists.txt @@ -0,0 +1,15 @@ +# Define the files we need to compile. +# Anything not in this list will not be compiled into mlpack. +set(SOURCES + randomized_block_krylov_svd.hpp + randomized_block_krylov_svd.cpp +) + +# Add directory name to sources. +set(DIR_SRCS) +foreach(file ${SOURCES}) + set(DIR_SRCS ${DIR_SRCS} ${CMAKE_CURRENT_SOURCE_DIR}/${file}) +endforeach() +# Append sources (with directory name) to list of all mlpack sources (used at +# the parent scope). +set(MLPACK_SRCS ${MLPACK_SRCS} ${DIR_SRCS} PARENT_SCOPE) diff --git a/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp new file mode 100644 index 0000000000..9009c5f642 --- /dev/null +++ b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp @@ -0,0 +1,88 @@ +/** + * @file randomized_block_krylov_svd.cpp + * @author Marcus Edel + * + * Implementation of the randomized SVD method. + * + * mlpack is free software; you may redistribute it and/or modify it under the + * terms of the 3-clause BSD license. You should have received a copy of the + * 3-clause BSD license along with mlpack. If not, see + * http://www.opensource.org/licenses/BSD-3-Clause for more information. + */ + +#include "randomized_block_krylov_svd.hpp" + +namespace mlpack { +namespace svd { + +RandomizedBlockKrylovSVD::RandomizedBlockKrylovSVD(const arma::mat& data, + arma::mat& u, + arma::vec& s, + arma::mat& v, + const size_t maxIterations, + const size_t rank, + const size_t blockSize) : + maxIterations(maxIterations), + blockSize(blockSize) +{ + if (rank == 0) + { + Apply(data, u, s, v, data.n_rows); + } + else + { + Apply(data, u, s, v, rank); + } +} + +RandomizedBlockKrylovSVD::RandomizedBlockKrylovSVD(const size_t maxIterations, + const size_t blockSize) : + maxIterations(maxIterations), + blockSize(blockSize) +{ + /* Nothing to do here */ +} + +void RandomizedBlockKrylovSVD::Apply(const arma::mat& data, + arma::mat& u, + arma::vec& s, + arma::mat& v, + const size_t rank) +{ + arma::mat Q, R, block; + + if (blockSize == 0) + { + blockSize = rank + 10; + } + + // Random block initialization. + arma::mat G = arma::randn(data.n_rows, blockSize); + + // Construct and orthonormalize Krlov subspace. + arma::mat K(data.n_rows, blockSize * (maxIterations + 1)); + arma::qr_econ(block, R, data * G); + + // Copy the temporary memory to the right place. + K.submat(0, 0, block.n_rows - 1, block.n_cols - 1) = block; + + for (size_t i = 0, b = block.n_cols; i < maxIterations; ++i, + b += block.n_cols) + { + arma::qr_econ(block, R, data * (data.t() * block)); + K.submat(0, b, block.n_rows - 1, b + block.n_cols - 1) = block; + } + + arma::qr_econ(Q, R, K); + + // Approximate eigenvalues and eigenvectors using Rayleigh–Ritz method. + arma::svd_econ(u, s, v, Q.t() * data); + + // Do economical singular value decomposition and compute only the + // approximations of the left singular vectors by using the centered data + // applied to Q. + u = Q * u; +} + +} // namespace svd +} // namespace mlpack diff --git a/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.hpp b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.hpp new file mode 100644 index 0000000000..964a654ea9 --- /dev/null +++ b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.hpp @@ -0,0 +1,130 @@ +/** + * @file randomized_block_krylov_svd.hpp + * @author Marcus Edel + * + * An implementation of the randomized SVD method. + * + * mlpack is free software; you may redistribute it and/or modify it under the + * terms of the 3-clause BSD license. You should have received a copy of the + * 3-clause BSD license along with mlpack. If not, see + * http://www.opensource.org/licenses/BSD-3-Clause for more information. + */ + +#ifndef MLPACK_METHODS_BLOCK_KRYLOV_SVD_RANDOMIZED_BLOCK_KRYLOV_SVD_HPP +#define MLPACK_METHODS_BLOCK_KRYLOV_SVD_RANDOMIZED_BLOCK_KRYLOV_SVD_HPP + +#include + +namespace mlpack { +namespace svd { + +/** + * Randomized block krylov SVD is a matrix factorization that is based on + * randomized matrix approximation techniques, developed in in + * "Randomized Block Krylov Methods for Stronger and Faster Approximate + * Singular Value Decomposition". + * + * For more information, see the following. + * + * @code + * @inproceedings{Musco2015, + * author = {Cameron Musco and Christopher Musco}, + * title = {Randomized Block Krylov Methods for Stronger and Faster + * Approximate Singular Value Decomposition}, + * booktitle = {Advances in Neural Information Processing Systems 28: Annual + * Conference on Neural Information Processing Systems 2015, + * December 7-12, 2015, Montreal, Quebec, Canada}, + * pages = {1396--1404}, + * year = {2015}, + * } + * @endcode + * + * @code + * + * An example of how to use the interface is shown below: + * + * @code + * arma::mat data; // Rating data in the form of coordinate list. + * + * const size_t rank = 20; // Rank used for the decomposition. + * + * // Make a RandomizedBlockKrylovSVD object. + * RandomizedBlockKrylovSVD bSVD(); + * + * arma::mat u, s, v; + * + * // Use the Apply() method to get a factorization. + * bSVD.Apply(data, u, s, v, rank); + * @endcode + */ +class RandomizedBlockKrylovSVD +{ + public: + /** + * Create object for the randomized block krylov SVD method. + * + * @param data Data matrix. + * @param u First unitary matrix. + * @param v Second unitary matrix. + * @param s Diagonal matrix of singular values. + * @param maxIterations Number of iterations for the power method + * (Default: 2). + * @param rank Rank of the approximation (Default: number of rows.) + * @param blockSize The block size, must be >= rank (Default: rank + 10). + */ + RandomizedBlockKrylovSVD(const arma::mat& data, + arma::mat& u, + arma::vec& s, + arma::mat& v, + const size_t maxIterations = 2, + const size_t rank = 0, + const size_t blockSize = 0); + + /** + * Create object for the randomized block krylov SVD method. + * + * @param maxIterations Number of iterations for the power method + * (Default: 2). + * @param blockSize The block size, must be >= rank (Default: rank + 10). + */ + RandomizedBlockKrylovSVD(const size_t maxIterations = 2, + const size_t blockSize = 0); + + /** + * Apply Principal Component Analysis to the provided data set using the + * randomized block krylov SVD. + * + * @param data Data matrix. + * @param u First unitary matrix. + * @param v Second unitary matrix. + * @param s Diagonal matrix of singular values. + * @param rank Rank of the approximation. + */ + void Apply(const arma::mat& data, + arma::mat& u, + arma::vec& s, + arma::mat& v, + const size_t rank); + + //! Get the number of iterations for the power method. + size_t MaxIterations() const { return maxIterations; } + //! Modify the number of iterations for the power method. + size_t& MaxIterations() { return maxIterations; } + + //! Get the block size. + size_t BlockSize() const { return blockSize; } + //! Modify the block size. + size_t& BlockSize() { return blockSize; } + + private: + //! Locally stored number of iterations for the power method. + size_t maxIterations; + + //! The block size value. + size_t blockSize; +}; + +} // namespace svd +} // namespace mlpack + +#endif From 4db6379cd8297932145aff7df9f649d779e7186f Mon Sep 17 00:00:00 2001 From: Marcus Edel Date: Sat, 18 Feb 2017 18:35:49 +0100 Subject: [PATCH 2/6] Add test cases for the randomized block krylov svd method. --- src/mlpack/tests/CMakeLists.txt | 1 + src/mlpack/tests/block_krylov_svd_test.cpp | 116 +++++++++++++++++++++ 2 files changed, 117 insertions(+) create mode 100644 src/mlpack/tests/block_krylov_svd_test.cpp diff --git a/src/mlpack/tests/CMakeLists.txt b/src/mlpack/tests/CMakeLists.txt index 514c453450..8e95e6069b 100644 --- a/src/mlpack/tests/CMakeLists.txt +++ b/src/mlpack/tests/CMakeLists.txt @@ -11,6 +11,7 @@ add_executable(mlpack_test armadillo_svd_test.cpp aug_lagrangian_test.cpp binarize_test.cpp + block_krylov_svd_test.cpp cf_test.cpp cli_test.cpp convolution_test.cpp diff --git a/src/mlpack/tests/block_krylov_svd_test.cpp b/src/mlpack/tests/block_krylov_svd_test.cpp new file mode 100644 index 0000000000..4ca9c6258d --- /dev/null +++ b/src/mlpack/tests/block_krylov_svd_test.cpp @@ -0,0 +1,116 @@ +/** + * @file block_krylov_svd_test.cpp + * @author Marcus Edel + * + * Test file for the Randomized Block Krylov SVD class. + * + * mlpack is free software; you may redistribute it and/or modify it under the + * terms of the 3-clause BSD license. You should have received a copy of the + * 3-clause BSD license along with mlpack. If not, see + * http://www.opensource.org/licenses/BSD-3-Clause for more information. + */ + +#include +#include + +#include +#include "test_tools.hpp" + +BOOST_AUTO_TEST_SUITE(BlockKrylovSVDTest); + +using namespace mlpack; + +// Generate a low rank matrix with bell-shaped singular values. +void CreateNoisyLowRankMatrix(arma::mat& data, + const size_t rows, + const size_t cols, + const size_t rank, + const double strength) +{ + arma::mat R, U, V; + const size_t n = std::min(rows, cols); + + arma::qr_econ(U, R, arma::randn(rows, n)); + arma::qr_econ(V, R, arma::randn(cols, n)); + + arma::vec ids = arma::linspace(0, n - 1, n); + + arma::vec lowRank = ((1 - strength) * + arma::exp(-1.0 * arma::pow((ids / rank), 2))); + arma::vec tail = strength * arma::exp(-0.1 * ids / rank); + + arma::mat s = arma::eye(n, n) * (lowRank + tail); + data = (U * s) * V.t(); +} + +/** + * The reconstruction and sigular value error of the obtained SVD should be + * small. + */ +BOOST_AUTO_TEST_CASE(RandomizedBlockKrylovSVDReconstructionError) +{ + arma::mat U = arma::randn(3, 20); + arma::mat V = arma::randn(10, 3); + + arma::mat R; + arma::qr_econ(U, R, U); + arma::qr_econ(V, R, V); + + arma::mat s = arma::diagmat(arma::vec("1 0.1 0.01")); + + arma::mat data = arma::trans(U * arma::diagmat(s) * V.t()); + + // Center the data into a temporary matrix. + arma::mat centeredData; + math::Center(data, centeredData); + + arma::mat U1, U2, V1, V2; + arma::vec s1, s2, s3; + + arma::svd_econ(U1, s1, V1, centeredData); + + svd::RandomizedBlockKrylovSVD rSVD(20, 10); + rSVD.Apply(centeredData, U2, s2, V2, 3); + + // Use the same amount of data for the compariosn (matrix rank). + s3 = s1.subvec(0, s2.n_elem - 1); + + // The sigular value error should be small. + double error = arma::norm(s2 - s3, "frob") / arma::norm(s2, "frob"); + BOOST_REQUIRE_SMALL(error, 1e-5); + + arma::mat reconstruct = U2 * arma::diagmat(s2) * V2.t(); + + // The relative reconstruction error should be small. + error = arma::norm(centeredData - reconstruct, "frob") / + arma::norm(centeredData, "frob"); + BOOST_REQUIRE_SMALL(error, 1e-5); +} + +/* + * Check if the method can handle noisy matrices. + */ +BOOST_AUTO_TEST_CASE(RandomizedBlockKrylovSVDNosiyLowRankTest) +{ + arma::mat data; + CreateNoisyLowRankMatrix(data, 100, 1000, 5, 1.0); + + const size_t rank = 5; + + arma::mat U1, U2, V1, V2; + arma::vec s1, s2, s3; + + arma::svd_econ(U1, s1, V1, data); + + svd::RandomizedBlockKrylovSVD rSVDA(data, U2, s2, V2, 1, rank, 5); + + double error = arma::max(arma::abs(s1.subvec(0, rank) - s2.subvec(0, rank))); + BOOST_REQUIRE_SMALL(error, 0.1); + + svd::RandomizedBlockKrylovSVD rSVDB(data, U2, s2, V2, 10, rank, 20); + + error = arma::max(arma::abs(s1.subvec(0, rank) - s2.subvec(0, rank))); + BOOST_REQUIRE_SMALL(error, 1e-3); +} + +BOOST_AUTO_TEST_SUITE_END(); From 397ef54120ceed8ad463e9db6d93f158e6a02eff Mon Sep 17 00:00:00 2001 From: Marcus Edel Date: Sun, 19 Feb 2017 17:53:02 +0100 Subject: [PATCH 3/6] Remove unused open @code block. --- .../methods/block_krylov_svd/randomized_block_krylov_svd.hpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.hpp b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.hpp index 964a654ea9..f7d77ab4c5 100644 --- a/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.hpp +++ b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.hpp @@ -39,8 +39,6 @@ namespace svd { * } * @endcode * - * @code - * * An example of how to use the interface is shown below: * * @code From e879e68a77fb0cbc2ed0a1817af77ab6da9a2e7b Mon Sep 17 00:00:00 2001 From: Marcus Edel Date: Wed, 5 Apr 2017 01:21:45 +0200 Subject: [PATCH 4/6] Use auxiliary memory to avoid copy in the iteration step. --- .../randomized_block_krylov_svd.cpp | 25 ++++++++++++------- .../randomized_block_krylov_svd.hpp | 2 +- 2 files changed, 17 insertions(+), 10 deletions(-) diff --git a/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp index 9009c5f642..41356e8e50 100644 --- a/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp +++ b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp @@ -2,7 +2,7 @@ * @file randomized_block_krylov_svd.cpp * @author Marcus Edel * - * Implementation of the randomized SVD method. + * Implementation of the randomized block krylov SVD method. * * mlpack is free software; you may redistribute it and/or modify it under the * terms of the 3-clause BSD license. You should have received a copy of the @@ -49,7 +49,7 @@ void RandomizedBlockKrylovSVD::Apply(const arma::mat& data, arma::mat& v, const size_t rank) { - arma::mat Q, R, block; + arma::mat Q, R, block, blockIteration; if (blockSize == 0) { @@ -61,16 +61,23 @@ void RandomizedBlockKrylovSVD::Apply(const arma::mat& data, // Construct and orthonormalize Krlov subspace. arma::mat K(data.n_rows, blockSize * (maxIterations + 1)); + + // Create a working matrix using data from writable auxiliary memory + // (K matrix). Doing so avoids an uncessary copy in upcoming step. + block = arma::mat(K.memptr(), data.n_rows, blockSize, false); arma::qr_econ(block, R, data * G); - // Copy the temporary memory to the right place. - K.submat(0, 0, block.n_rows - 1, block.n_cols - 1) = block; - - for (size_t i = 0, b = block.n_cols; i < maxIterations; ++i, - b += block.n_cols) + for (size_t blockOffset = block.n_elem; blockOffset < K.n_elem; + blockOffset += block.n_elem) { - arma::qr_econ(block, R, data * (data.t() * block)); - K.submat(0, b, block.n_rows - 1, b + block.n_cols - 1) = block; + // Temporary working matrix to store the result in the correct place. + blockIteration = arma::mat(K.memptr() + blockOffset, data.n_rows, + blockSize, false); + + arma::qr_econ(blockIteration, R, data * (data.t() * block)); + + // Update working matrix for the next iteration. + block = arma::mat(K.memptr() + blockOffset, data.n_rows, blockSize, false); } arma::qr_econ(Q, R, K); diff --git a/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.hpp b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.hpp index f7d77ab4c5..06ef8b4c2b 100644 --- a/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.hpp +++ b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.hpp @@ -2,7 +2,7 @@ * @file randomized_block_krylov_svd.hpp * @author Marcus Edel * - * An implementation of the randomized SVD method. + * An implementation of the randomized block krylov SVD method. * * mlpack is free software; you may redistribute it and/or modify it under the * terms of the 3-clause BSD license. You should have received a copy of the From 0a2e8cc34e3f84db4faf2a05dec50cf63bd0d6df Mon Sep 17 00:00:00 2001 From: Marcus Edel Date: Wed, 5 Apr 2017 01:24:12 +0200 Subject: [PATCH 5/6] Minor spelling fix. --- .../methods/block_krylov_svd/randomized_block_krylov_svd.cpp | 2 +- src/mlpack/tests/block_krylov_svd_test.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp index 41356e8e50..e0f689b63e 100644 --- a/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp +++ b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp @@ -59,7 +59,7 @@ void RandomizedBlockKrylovSVD::Apply(const arma::mat& data, // Random block initialization. arma::mat G = arma::randn(data.n_rows, blockSize); - // Construct and orthonormalize Krlov subspace. + // Construct and orthonormalize Krylov subspace. arma::mat K(data.n_rows, blockSize * (maxIterations + 1)); // Create a working matrix using data from writable auxiliary memory diff --git a/src/mlpack/tests/block_krylov_svd_test.cpp b/src/mlpack/tests/block_krylov_svd_test.cpp index 4ca9c6258d..5d7e6e2e28 100644 --- a/src/mlpack/tests/block_krylov_svd_test.cpp +++ b/src/mlpack/tests/block_krylov_svd_test.cpp @@ -90,7 +90,7 @@ BOOST_AUTO_TEST_CASE(RandomizedBlockKrylovSVDReconstructionError) /* * Check if the method can handle noisy matrices. */ -BOOST_AUTO_TEST_CASE(RandomizedBlockKrylovSVDNosiyLowRankTest) +BOOST_AUTO_TEST_CASE(RandomizedBlockKrylovSVDNoisyLowRankTest) { arma::mat data; CreateNoisyLowRankMatrix(data, 100, 1000, 5, 1.0); From 3026b5026ef1c71ae49e666ea55ecddef08634f8 Mon Sep 17 00:00:00 2001 From: Marcus Edel Date: Wed, 5 Apr 2017 23:50:41 +0200 Subject: [PATCH 6/6] The default setting of strict (strict = false use the auxiliary memory until a size change) in versions 5.600 and earlier is true, so to make sure it's always false we set it explicitly. --- .../methods/block_krylov_svd/randomized_block_krylov_svd.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp index e0f689b63e..702d706fd2 100644 --- a/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp +++ b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp @@ -64,7 +64,7 @@ void RandomizedBlockKrylovSVD::Apply(const arma::mat& data, // Create a working matrix using data from writable auxiliary memory // (K matrix). Doing so avoids an uncessary copy in upcoming step. - block = arma::mat(K.memptr(), data.n_rows, blockSize, false); + block = arma::mat(K.memptr(), data.n_rows, blockSize, false, false); arma::qr_econ(block, R, data * G); for (size_t blockOffset = block.n_elem; blockOffset < K.n_elem; @@ -77,7 +77,8 @@ void RandomizedBlockKrylovSVD::Apply(const arma::mat& data, arma::qr_econ(blockIteration, R, data * (data.t() * block)); // Update working matrix for the next iteration. - block = arma::mat(K.memptr() + blockOffset, data.n_rows, blockSize, false); + block = arma::mat(K.memptr() + blockOffset, data.n_rows, blockSize, false, + false); } arma::qr_econ(Q, R, K);