diff --git a/src/mlpack/methods/CMakeLists.txt b/src/mlpack/methods/CMakeLists.txt
index b7bedde853..e0ceb91c56 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
   dbscan
   decision_stump
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..702d706fd2
--- /dev/null
+++ b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.cpp
@@ -0,0 +1,96 @@
+/**
+ * @file randomized_block_krylov_svd.cpp
+ * @author Marcus Edel
+ *
+ * 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
+ * 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, blockIteration;
+
+  if (blockSize == 0)
+  {
+    blockSize = rank + 10;
+  }
+
+  // Random block initialization.
+  arma::mat G = arma::randn(data.n_rows, blockSize);
+
+  // Construct and orthonormalize Krylov 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, false);
+  arma::qr_econ(block, R, data * G);
+
+  for (size_t blockOffset = block.n_elem; blockOffset < K.n_elem;
+      blockOffset += block.n_elem)
+  {
+    // 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,
+        false);
+  }
+
+  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..06ef8b4c2b
--- /dev/null
+++ b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.hpp
@@ -0,0 +1,128 @@
+/**
+ * @file randomized_block_krylov_svd.hpp
+ * @author Marcus Edel
+ *
+ * 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
+ * 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 <mlpack/prereqs.hpp>
+
+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
+ *
+ * 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
diff --git a/src/mlpack/tests/CMakeLists.txt b/src/mlpack/tests/CMakeLists.txt
index b93077b25e..475764580e 100644
--- a/src/mlpack/tests/CMakeLists.txt
+++ b/src/mlpack/tests/CMakeLists.txt
@@ -12,6 +12,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..5d7e6e2e28
--- /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 <mlpack/core.hpp>
+#include <mlpack/methods/block_krylov_svd/randomized_block_krylov_svd.hpp>
+
+#include <boost/test/unit_test.hpp>
+#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<arma::mat>(rows, n));
+  arma::qr_econ(V, R, arma::randn<arma::mat>(cols, n));
+
+  arma::vec ids = arma::linspace<arma::vec>(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<arma::mat>(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<arma::mat>(3, 20);
+  arma::mat V = arma::randn<arma::mat>(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(RandomizedBlockKrylovSVDNoisyLowRankTest)
+{
+  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();