diff --git a/HISTORY.md b/HISTORY.md
index 231b5e7a5f..e2e90effb5 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -34,6 +34,9 @@
* Fix a bug for the stddev and mean in `RandNormal()` #(3651).
+ * Allow PCA to take different matrix types (#3677).
+
+>>>>>>> master
### mlpack 4.3.0
###### 2023-11-27
* Fix include ordering issue for `LinearRegression` (#3541).
diff --git a/doc/css/gfm-mod.css b/doc/css/gfm-mod.css
index 5cd333d7a9..29518d3836 100644
--- a/doc/css/gfm-mod.css
+++ b/doc/css/gfm-mod.css
@@ -1062,6 +1062,7 @@ div#sidebar {
top: 5px;
min-width: 200px;
font-size: 90%;
+ max-width: 15.1515%;
}
div#sidebar ul {
diff --git a/doc/index.md b/doc/index.md
index 83b2b2a624..d9ebf87571 100644
--- a/doc/index.md
+++ b/doc/index.md
@@ -120,7 +120,7 @@ Prepare data for machine learning algorithms.
Transform data from one space to another.
-
+ * [`PCA`](user/methods/pca.md): principal components analysis
### Modeling utilities
diff --git a/doc/sidebar.html b/doc/sidebar.html
index 63f668b5bf..478ae85ff7 100644
--- a/doc/sidebar.html
+++ b/doc/sidebar.html
@@ -172,9 +172,20 @@ when the sidebar is built for each page.
-
- Transformations
-
+
+
+
+ Transformations
+
+
+
+
diff --git a/doc/user/core.md b/doc/user/core.md
index e3ec6a9122..4fd8c02da8 100644
--- a/doc/user/core.md
+++ b/doc/user/core.md
@@ -15,6 +15,9 @@ classes, each of which are documented on this page.
mlpack provides a number of additional mathematical utility classes and
functions on top of Armadillo.
+ * [Aliases](#aliases): utilities to create and manage aliases (`MakeAlias()`,
+ `ClearAlias()`, `UnwrapAlias()`).
+
* [`Range`](#range): simple mathematical range (i.e. `[0, 3]`)
* [`ColumnCovariance()`](#columncovariance): compute covariance of
@@ -31,9 +34,6 @@ functions on top of Armadillo.
* [Logarithmic utilities](#logarithmic-utilities): `LogAdd()`, `AccuLog()`,
`LogSumExp()`, `LogSumExpT()`.
-
-
* [`MultiplyCube2Cube()`](#multiplycube2cube): multiply each slice in a cube by each slice in another cube
* [`MultiplyMat2Cube()`](#multiplymat2cube): multiply a matrix by each slice in a cube
* [`MultiplyCube2Mat()`](#multiplycube2mat): multiply each slice in a cube by a matrix
@@ -47,6 +47,80 @@ functions on top of Armadillo.
---
+### Aliases
+
+Aliases are matrix, vector, or cube objects that share memory with another
+matrix, vector, or cube. They are often used internally inside of mlpack to
+avoid copies.
+
+***Important caveats about aliases***:
+
+ - An alias represents the same memory block as the input. As such, changes to
+ the alias object will also be reflected in the original object.
+
+ - The `MakeAlias()` function is not guaranteed to return an alias; it only
+ returns an alias *if possible*, and makes a copy otherwise.
+
+ - If `mat` goes out of scope or is destructed, then `a` ***becomes invalid***.
+ _You are responsible for ensuring an invalid alias is not used!_
+
+---
+
+ * `MakeAlias(a, mat, rows, cols, strict=true)`
+ - Make `a` into an alias of `mat` with the given size.
+ - If `strict` is `true`, the size of `a` cannot be changed.
+ - `mat` and `a` should have the same matrix type (e.g. `arma::mat`,
+ `arma::fmat`, `arma::sp_mat`).
+ - If an alias cannot be created, the matrix will be copied. Sparse types
+ cannot have aliases and will be copied.
+
+ * `MakeAlias(a, cube, rows, cols, slices, strict=true)`
+ - Make `a` into an alias of `cube` with the given size.
+ - If `strict` is `true`, the size of `a` cannot be changed.
+ - `cube` and `a` should have the same matrix type (e.g. `arma::cube`,
+ `arma::fcube`).
+ - If an alias cannot be created, the matrix will be copied.
+
+ * `MakeAlias(a, memptr, rows, cols, strict=true)`
+ - Make `a` into an alias of the memory block starting at `memptr` of size
+ `rows` by `cols`.
+ - The memory at `memptr` should be arranged in a [column-major
+ ordering](matrices.md#representing-data-in-mlpack).
+ - If `strict` is `true`, the size of `a` cannot be changed.
+ - `a` should be a dense matrix type (e.g. `arma::mat`, `arma::fmat`), and
+ `memptr` should be a non-const pointer of the matrix's element type (e.g.
+ `double*`, `float*`).
+
+ * `MakeAlias(a, memptr, rows, cols, slices, strict=true)`
+ - Make `a` into an alias of the memory block starting at `memptr` of size
+ `rows` by `cols` by `slices`.
+ - The memory at `memptr` should be arranged in a [column-major
+ ordering](matrices.md#representing-data-in-mlpack).
+ - If `strict` is `true`, the size of `a` cannot be changed.
+ - `a` should be a cube type (e.g. `arma::cube`, `arma::fcube`), and `memptr`
+ should be a non-const pointer of the matrix's element type (e.g. `double*`,
+ `float*`).
+
+---
+
+ * `ClearAlias(a)`
+ - If `a` is an alias, reset `a` to an empty matrix, without modifying the
+ aliased memory. `a` is no longer an alias after this call.
+
+---
+
+ * `UnwrapAlias(a, in)`
+ - If `in` is a matrix type (e.g. `arma::mat`), make `a` into an alias of
+ `in`.
+ - If `in` is not a matrix type, but instead, e.g., an Armadillo expression,
+ fill `a` with the results of the evaluated expression `in`.
+ - This can be used in place of, e.g., `a = in`, to avoid a copy when
+ possible.
+ - `a` should be a matrix type that matches the type of the expression or
+ matrix `in`.
+
+---
+
### `Range`
The `Range` class represents a simple mathematical range (i.e. `[0, 3]`),
diff --git a/doc/user/load_save.md b/doc/user/load_save.md
index 9785c968e9..c7fa95627c 100644
--- a/doc/user/load_save.md
+++ b/doc/user/load_save.md
@@ -157,7 +157,7 @@ With a `data::DatasetInfo` object, categorical data can be loaded:
contents of the file:
- `.csv`, `.tsv`, or `.txt` for CSV/TSV (tab-separated)/ASCII
(space-separated)
- - `.arff` for [ARFF](https://www.cs.waikato.ac.nz/~ml/weka/arff.html)
+ - `.arff` for [ARFF](https://ml.cms.waikato.ac.nz/weka/arff.html)
* `matrix` is an `arma::mat&`, `arma::Mat&`, or similar (e.g., a
reference to an Armadillo object that data will be loaded into or saved
@@ -723,7 +723,7 @@ The format of mixed categorical data is detected automatically based on the
file extension and inspecting the file contents:
- `.csv`, `.txt`, or `.tsv` indicates CSV/TSV/ASCII format
- - `.arff` indicates [ARFF](https://www.cs.waikato.ac.nz/~ml/weka/arff.html)
+ - `.arff` indicates [ARFF](https://ml.cms.waikato.ac.nz/weka/arff.html)
---
diff --git a/doc/user/methods/pca.md b/doc/user/methods/pca.md
new file mode 100644
index 0000000000..ccf65048d2
--- /dev/null
+++ b/doc/user/methods/pca.md
@@ -0,0 +1,326 @@
+## `PCA`
+
+The `PCA` class implements principal components analysis (PCA), a standard
+machine learning data preparation technique. PCA can be used to reduce the
+number of dimensions in a dataset, or to preserve a certain percentage of the
+variance of a dataset.
+
+By default, `PCA` uses the full exact singular value decomposition (SVD), but
+supports the use of other more efficient decompositions, including approximate
+singular value decompositions.
+
+#### Simple usage example:
+
+```c++
+// Use PCA to reduce the number of dimensions to 5 on uniform random data.
+
+// This dataset is uniform random in 10 dimensions.
+// Replace with a data::Load() call or similar for a real application.
+arma::mat dataset(10, 1000, arma::fill::randu); // 1000 points.
+
+mlpack::PCA pca; // Step 1: create PCA object.
+pca.Apply(dataset, 5); // Step 2: reduce data dimension to 5.
+
+// Print some information about the modified dataset.
+std::cout << "The transformed data matrix has size " << dataset.n_rows /* 5 */
+ << " x " << dataset.n_cols << "." << std::endl;
+```
+More examples...
+
+#### Quick links:
+
+ * [Constructors](#constructors): create `PCA` objects.
+ * [`Apply()`](#applying-transformations): apply PCA transformation to data.
+ * [Examples](#simple-examples) of simple usage and links to detailed example
+ projects.
+ * [Template parameters](#advanced-functionality-different-decomposition-strategies)
+ for using different decomposition strategies.
+
+#### See also:
+
+
+
+ * [mlpack preprocessing utilities](../../index.md#preprocessing-utilities)
+ * [mlpack transformations](../../index.md#transformations)
+ * [Principal component analysis on Wikipedia](https://en.wikipedia.org/wiki/Principal_component_analysis)
+
+### Constructors
+
+ * `pca = PCA(scaleData=false)`
+ - Construct a `PCA` object.
+ - If `scaleData` is `true`, then all dimensions will have variance scaled to
+ 1 before applying PCA.
+ - The `scaleData` parameter can be inspected with `pca.ScaleData()`, and also
+ set; `pca.ScaleData() = true` will enable data variance scaling.
+
+---
+
+ * `pca = PCA(scaleData, decompositionPolicy)`
+ - Construct a `PCA` object with a custom decomposition policy.
+ - See the documentation for using
+ [different decomposition strategies](#advanced-functionality-different-decomposition-strategies).
+
+### Applying Transformations
+
+ * `pca.Apply(data, transformedData)`
+ * `pca.Apply(data, transformedData, eigVal)`
+ * `pca.Apply(data, transformedData, eigVal, eigVec)`
+ - Transform the
+ [column-major matrix](../matrices.md#representing-data-in-mlpack) `data`
+ using PCA, storing the result in `transformedData`.
+ - `data` should be a floating-point matrix (e.g. `arma::mat`, `arma::fmat`,
+ `arma::sp_mat`, etc.) or an expression that evaluates to one.
+ - `transformedData` should be a dense floating-point matrix (e.g.,
+ `arma::mat`, `arma::sp_mat`).
+ - The size of `transformedData` will be the same as the size of `data`.
+ - Dimensions in `transformedData` will be ordered decreasing in variance;
+ that is, the first row of `transformedData` will correspond to the
+ dimension with maximum variance.
+ - Optionally, eigenvalues and eigenvectors of the covariance matrix can be
+ returned:
+ * If specified, `eigVal` should be a dense floating-point vector (e.g.
+ `arma::vec`, `arma::fvec`, etc.) and will be filled with the eigenvalues
+ of `transformedData`.
+ * If specified, `eigvec` should be a dense floating-point matrix (e.g.
+ `arma::mat`, `arma::fmat`, etc.) and will be filled with the eigenvectors
+ of `transformedData`.
+
+---
+
+ * `double varRetained = pca.Apply(data, transformedData, newDimension)`
+ - Use PCA to reduce the number of dimensions in the
+ [column-major matrix](../matrices.md#representing-data-in-mlpack) `data`
+ to `newDimension`, storing the result in `transformedData`.
+ - `data` should be a floating-point matrix (e.g. `arma::mat`,
+ `arma::fmat`, `arma::sp_mat`, etc.) or an expression that evaluates to
+ one.
+ - `transformedData` should be a dense floating-point matrix with the same
+ element type as `data` (e.g. `arma::mat`, `arma::fmat`).
+ - `transformedData` will have `newDimension` rows after the transformation.
+ - Returns a `double` indicating the percentage of variance retained (between
+ `0.0` and `1.0`).
+
+---
+
+ * `double varRetained = pca.Apply(data, transformedData, varianceToKeep)`
+ - Use PCA to retain the dimensions of the
+ [column-major matrix](../matrices.md#representing-data-in-mlpack) `data`
+ that capture a factor of `varianceToKeep` of the data variance.
+ - `data` should be a floating-point matrix (e.g. `arma::mat`, `arma::fmat`,
+ `arma::sp_mat`, etc.) or an expression that evaluates to one.
+ - `transformedData` should be a dense floating-point matrix with the same
+ element type as `data` (e.g. `arma::mat`, `arma::fmat`).
+ - `transformedData` will have `newDimension` rows after the transformation.
+ - `varianceToKeep` should be a floating-point value between `0.0` and `1.0`.
+ If `1.0`, all of the data variance is retained, and this is equivalent to
+ the first version of `Apply()` (above).
+ - Returns a `double` indicating the percentage of variance actually retained
+ (between `0.0` and `1.0`).
+
+---
+
+ * `double varRetained = pca.Apply(data, newDimension)`
+ * `double varRetained = pca.Apply(data, varianceToKeep)`
+ - In-place versions of the two `Apply()` functions above.
+ - Equivalent to `pca.Apply(data, data, newDimension)` or
+ `pca.Apply(data, data, varianceToKeep)`.
+ - `data` should be a dense floating-point matrix (e.g. `arma::mat`,
+ `arma::fmat`, etc.).
+
+---
+
+### Simple Examples
+
+See also the [simple usage example](#simple-usage-example) for a trivial usage
+of the `PCA` class.
+
+---
+
+Apply PCA to a dataset, keeping dimensions that capture 90% of the data
+variance.
+
+```c++
+// See https://datasets.mlpack.org/satellite.train.csv.
+arma::mat data;
+mlpack::data::Load("satellite.train.csv", data, true);
+const size_t origDim = data.n_rows;
+
+mlpack::PCA pca;
+
+// Keep 90% of the data variance.
+pca.Apply(data, 0.9);
+
+std::cout << "PCA kept " << data.n_rows << " of " << origDim << " dimensions "
+ << "to capture 90\% of the data variance." << std::endl;
+```
+
+---
+
+Apply PCA to a 32-bit floating point dataset with dimension scaling, keeping all
+dimensions, and printing the 5 largest eigenvalues of the covariance matrix of
+the transformed data.
+
+```c++
+// See https://datasets.mlpack.org/iris.csv.
+arma::fmat data;
+mlpack::data::Load("iris.csv", data, true);
+
+mlpack::PCA pca(true /* scale data when transforming */);
+
+arma::fvec eigval;
+arma::fmat transformedData;
+
+pca.Apply(data, transformedData, eigval);
+
+std::cout << "First point, before PCA: " << data.col(0).t();
+std::cout << "First point, after PCA: " << transformedData.col(0).t();
+std::cout << std::endl;
+
+// Now print the top 5 eigenvalues.
+for (size_t i = 0; i < 5; ++i)
+ std::cout << "Eigenvalue " << i << ": " << eigval[i] << "." << std::endl;
+```
+
+---
+
+Apply PCA to a random sparse dataset, to reduce the dimensionality to a
+20-dimensional dense dataset.
+
+```c++
+arma::sp_mat data;
+// This dataset has 10k points in 1k dimensions, with 1% density.
+data.sprandn(1000, 10000, 0.01);
+
+mlpack::PCA pca(true /* scale data when transforming */);
+
+arma::mat transformedData;
+const double varianceRetained = pca.Apply(data, transformedData, 20);
+
+std::cout << "First point, before PCA: " << data.col(0).t();
+std::cout << "First point, after PCA: " << transformedData.col(0).t();
+
+// Note that for random uniform data, this won't capture very much of the
+// variance! It would be much more for a real, structured dataset.
+std::cout << "50 dimensions captured " << (100.0 * varianceRetained) << "\% of "
+ << "the data variance." << std::endl;
+```
+
+---
+
+### Advanced Functionality: Different Decomposition Strategies
+
+By default, `PCA` uses the full exact singular value decomposition (SVD) to
+transform data. However, for very large datasets, it may be faster to use
+alternative strategies, some of which may be approximate. The `PCA` class has
+one template parameter that allows different decomposition strategies to be
+used. The full signature of the class is:
+
+```
+PCA
+```
+
+`DecompositionPolicy` specifies the strategy to be used to compute the singular
+values and vectors of a data matrix.
+
+Several decomposition policies are already implemented and ready for drop-in
+usage:
+
+ * `ExactSVDPolicy` _(default)_: use Armadillo's `svd()` and `svd_econ()`
+ functions to compute the SVD
+ * `RandomizedSVDPCAPolicy`: use the randomized SVD algorithm to compute the SVD
+
+ * `RandomizedBlockKrylovSVDPolicy`: use the randomized Block Krylov SVD
+ algorithm to compute the SVD
+ * `QUICSVDPolicy`: use the tree-based `QUIC-SVD` algorithm to compute the SVD
+
+
+The simple example program below uses all four decomposition types on the same
+MNIST data, timing how long each decomposition takes.
+
+```c++
+arma::mat data;
+// See https://datasets.mlpack.org/mnist.train.csv.
+mlpack::data::Load("mnist.train.csv", data, true);
+
+arma::mat output1, output2, output3, output4;
+
+mlpack::PCA pca1;
+mlpack::PCA pca2;
+mlpack::PCA pca3;
+mlpack::PCA pca4;
+
+// Compute decompositions on all four, timing each one.
+arma::wall_clock c;
+
+c.tic();
+pca1.Apply(data, output1);
+const double pca1Time = c.toc();
+
+c.tic();
+pca2.Apply(data, output2);
+const double pca2Time = c.toc();
+
+c.tic();
+pca3.Apply(data, output3);
+const double pca3Time = c.toc();
+
+c.tic();
+pca4.Apply(data, output4);
+const double pca4Time = c.toc();
+
+std::cout << "PCA computation times for " << data.n_rows << " x " << data.n_cols
+ << " data:" << std::endl;
+std::cout << " - ExactSVDPolicy: " << pca1Time << "s."
+ << std::endl;
+std::cout << " - RandomizedSVDPCAPolicy: " << pca2Time << "s."
+ << std::endl;
+std::cout << " - RandomizedBlockKrylovSVDPolicy: " << pca3Time << "s."
+ << std::endl;
+std::cout << " - QUICSVDPolicy: " << pca4Time << "s."
+ << std::endl;
+```
+
+---
+
+#### Custom decomposition policies
+
+Instead of using the predefined classes above, it is also possible to implement
+fully custom functionality via a new decomposition policy. Any new
+decomposition policy must implement one method:
+
+```c++
+class CustomDecompositionPolicy
+{
+ public:
+ // Given input data `data` and `centeredData`, compute the singular value
+ // decomposition of the data, and then project the data onto the first `rank`
+ // singular vectors.
+ //
+ // * `data` is the input matrix. It is not guaranteed to be centered or
+ // scaled.
+ // * `centeredData` is the centered (and possibly scaled) version of the
+ // input matrix (e.g. the mean of each dimension is 0).
+ // * `transformedData` should be overwritten with the centered data's
+ // projection onto the singular vectors.
+ // * `svals` and `svecs` should be filled with the singular values and
+ // vectors of the centered data.
+ // * `rank` specifies the number of singular values/vectors to keep, and the
+ // dimension of `transformedData` should be equivalent to `rank`. `rank`
+ // will be at most equal to `data.n_rows`.
+ //
+ // * `InMatType` is a dense floating-point matrix type, but may be a subview
+ // or expression.
+ // * `MatType` is the type of matrix used to represent data, and will be a
+ // dense floating-point matrix type (e.g. `arma::mat`, `arma::fmat`,
+ // etc.).
+ // * `VecType` is the corresponding vector type to `MatType` (e.g., a
+ // `MatType` of `arma::mat` would mean a `VecType` of `arma::vec`, etc.).
+ template
+ static void Apply(const InMatType& data,
+ const MatType& centeredData,
+ MatType& transformedData,
+ VecType& svals,
+ MatType& svecs,
+ const size_t rank);
+};
+```
diff --git a/scripts/build-docs.sh b/scripts/build-docs.sh
index 49b6360164..104e59310c 100755
--- a/scripts/build-docs.sh
+++ b/scripts/build-docs.sh
@@ -320,7 +320,7 @@ do
echo "Checking links in $f...";
# To run checklink we have to strip out some perl stderr warnings...
- checklink -qs --follow-file-links --suppress-broken 405 "$f" 2>&1 |
+ checklink -qs --follow-file-links --suppress-broken 405 --suppress-broken 301 "$f" 2>&1 |
grep -v 'Use of uninitialized value' > checklink_out;
if [ -s checklink_out ];
then
diff --git a/src/mlpack/core/math/make_alias.hpp b/src/mlpack/core/math/make_alias.hpp
index f73578e6ed..93bcd0138a 100644
--- a/src/mlpack/core/math/make_alias.hpp
+++ b/src/mlpack/core/math/make_alias.hpp
@@ -16,89 +16,69 @@
namespace mlpack {
/**
- * Make an alias of a dense cube. If strict is true, then the alias cannot be
- * resized or pointed at new memory.
+ * Reconstruct `m` as an alias around the memory `newMem`, with size `numRows` x
+ * `numCols`.
*/
-template
-arma::Cube MakeAlias(arma::Cube& input,
- const bool strict = true)
+template
+void MakeAlias(MatType& m,
+ typename MatType::elem_type* newMem,
+ const size_t numRows,
+ const size_t numCols,
+ const bool strict = true,
+ const typename std::enable_if_t::value>* = 0)
{
- // Use the advanced constructor.
- return arma::Cube(input.memptr(), input.n_rows, input.n_cols,
- input.n_slices, false, strict);
+ // We use placement new to reinitialize the object, since the copy and move
+ // assignment operators in Armadillo will end up copying memory instead of
+ // making an alias.
+ m.~MatType();
+ new (&m) MatType(newMem, numRows, numCols, false, strict);
}
/**
- * Make an alias of a dense matrix. If strict is true, then the alias cannot be
- * resized or pointed at new memory.
+ * Reconstruct `c` as an alias around the memory` newMem`, with size `numRows` x
+ * `numCols` x `numSlices`.
*/
-template
-arma::Mat MakeAlias(arma::Mat& input,
- const bool strict = true)
+template
+void MakeAlias(CubeType& c,
+ typename CubeType::elem_type* newMem,
+ const size_t numRows,
+ const size_t numCols,
+ const size_t numSlices,
+ const bool strict = true,
+ const typename std::enable_if_t::value>* = 0)
{
- // Use the advanced constructor.
- return arma::Mat(input.memptr(), input.n_rows, input.n_cols, false,
- strict);
+ // We use placement new to reinitialize the object, since the copy and move
+ // assignment operators in Armadillo will end up copying memory instead of
+ // making an alias.
+ c.~CubeType();
+ new (&c) CubeType(newMem, numRows, numCols, numSlices, false, strict);
}
/**
- * Make an alias of a dense row. If strict is true, then the alias cannot be
- * resized or pointed at new memory.
+ * Make `m` an alias of `in`, using the given size.
*/
-template
-arma::Row MakeAlias(arma::Row& input,
- const bool strict = true)
+template
+void MakeAlias(arma::Mat& m,
+ const arma::Mat& in,
+ const size_t numRows,
+ const size_t numCols,
+ const bool strict = true)
{
- // Use the advanced constructor.
- return arma::Row(input.memptr(), input.n_elem, false, strict);
+ MakeAlias(m, (eT*) in.memptr(), numRows, numCols, strict);
}
/**
- * Make an alias of a dense column. If strict is true, then the alias cannot be
- * resized or pointed at new memory.
+ * Make `m` an alias of `in`, using the given size.
*/
-template
-arma::Col MakeAlias(arma::Col& input,
- const bool strict = true)
+template
+void MakeAlias(arma::SpMat& m,
+ const arma::SpMat& in,
+ const size_t /* numRows */,
+ const size_t /* numCols */,
+ const bool /* strict */)
{
- // Use the advanced constructor.
- return arma::Col(input.memptr(), input.n_elem, false, strict);
-}
-
-/**
- * Make a copy of a sparse matrix (an alias is not possible). The strict
- * parameter is ignored.
- */
-template
-arma::SpMat MakeAlias(const arma::SpMat& input,
- const bool /* strict */ = true)
-{
- // Make a copy...
- return arma::SpMat(input);
-}
-
-/**
- * Make a copy of a sparse row (an alias is not possible). The strict
- * parameter is ignored.
- */
-template
-arma::SpRow MakeAlias(const arma::SpRow& input,
- const bool /* strict */ = true)
-{
- // Make a copy...
- return arma::SpRow(input);
-}
-
-/**
- * Make a copy of a sparse column (an alias is not possible). The strict
- * parameter is ignored.
- */
-template
-arma::SpCol MakeAlias(const arma::SpCol& input,
- const bool /* strict */ = true)
-{
- // Make a copy...
- return arma::SpCol(input);
+ // We can't make aliases of sparse objects, so just copy it.
+ m = in;
}
/**
@@ -113,16 +93,15 @@ void ClearAlias(arma::Mat& mat)
}
/**
- * Clear an alias for a sparse matrix. This does nothing because no sparse
- * matrices can have aliases.
+ * Clear an alias so that no data is overwritten. This resets the matrix if it
+ * is an alias (and does nothing otherwise).
*/
template
void ClearAlias(arma::SpMat& /* mat */)
{
- // Nothing to do.
+ // We cannot make aliases of sparse matrices, so, nothing to do.
}
-
} // namespace mlpack
#endif
diff --git a/src/mlpack/core/math/math.hpp b/src/mlpack/core/math/math.hpp
index 66fcfbdd59..3e254863fa 100644
--- a/src/mlpack/core/math/math.hpp
+++ b/src/mlpack/core/math/math.hpp
@@ -25,5 +25,6 @@
#include "range.hpp"
#include "shuffle_data.hpp"
#include "trigamma.hpp"
+#include "unwrap_alias.hpp"
#endif
diff --git a/src/mlpack/core/math/unwrap_alias.hpp b/src/mlpack/core/math/unwrap_alias.hpp
new file mode 100644
index 0000000000..7afaf30d96
--- /dev/null
+++ b/src/mlpack/core/math/unwrap_alias.hpp
@@ -0,0 +1,38 @@
+/**
+ * @file core/math/unwrap_alias.hpp
+ * @author Ryan Curtin
+ *
+ * Make an alias of a matrix if possible, or unwrap it if an expression is
+ * given.
+ *
+ * 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_CORE_MATH_UNWRAP_ALIAS_HPP
+#define MLPACK_CORE_MATH_UNWRAP_ALIAS_HPP
+
+namespace mlpack {
+
+/**
+ * If `in` is an expression, unwrap it into `m`. If `in` is a matrix, then
+ * create `m` as an alias of `in`.
+ * `numCols`.
+ */
+template
+void UnwrapAlias(MatType& m, const MatType& in)
+{
+ MakeAlias(m, in, in.n_rows, in.n_cols);
+}
+
+template
+void UnwrapAlias(MatType& m,
+ const InMatType& in)
+{
+ m = in;
+}
+
+} // namespace mlpack
+
+#endif
diff --git a/src/mlpack/core/tree/cosine_tree/cosine_tree.hpp b/src/mlpack/core/tree/cosine_tree/cosine_tree.hpp
index 5966688e37..cb5b52fe25 100644
--- a/src/mlpack/core/tree/cosine_tree/cosine_tree.hpp
+++ b/src/mlpack/core/tree/cosine_tree/cosine_tree.hpp
@@ -20,14 +20,20 @@ namespace mlpack {
// Predeclare classes for CosineNodeQueue typedef.
class CompareCosineNode;
+
+template
class CosineTree;
// CosineNodeQueue typedef.
-typedef std::vector CosineNodeQueue;
+template
+using CosineNodeQueue = std::vector*>;
+template
class CosineTree
{
public:
+ typedef typename GetDenseColType::type VecType;
+
/**
* CosineTree constructor for the root node of the tree. It initializes the
* necessary variables required for splitting of the node, and building the
@@ -36,7 +42,7 @@ class CosineTree
*
* @param dataset Matrix for which cosine tree is constructed.
*/
- CosineTree(const arma::mat& dataset);
+ CosineTree(const MatType& dataset);
/**
* CosineTree constructor for nodes other than the root node of the tree. It
@@ -63,7 +69,7 @@ class CosineTree
* @param epsilon Error tolerance fraction for calculated subspace.
* @param delta Cumulative probability for Monte Carlo error lower bound.
*/
- CosineTree(const arma::mat& dataset,
+ CosineTree(const MatType& dataset,
const double epsilon,
const double delta);
@@ -110,10 +116,10 @@ class CosineTree
* @param newBasisVector Orthonormalized centroid of the node.
* @param addBasisVector Address to additional basis vector.
*/
- void ModifiedGramSchmidt(CosineNodeQueue& treeQueue,
- arma::vec& centroid,
- arma::vec& newBasisVector,
- arma::vec* addBasisVector = NULL);
+ void ModifiedGramSchmidt(CosineNodeQueue& treeQueue,
+ VecType& centroid,
+ VecType& newBasisVector,
+ VecType* addBasisVector = NULL);
/**
* Estimates the squared error of the projection of the input node's matrix
@@ -128,16 +134,16 @@ class CosineTree
* @param addBasisVector2 Address to second additional basis vector.
*/
double MonteCarloError(CosineTree* node,
- CosineNodeQueue& treeQueue,
- arma::vec* addBasisVector1 = NULL,
- arma::vec* addBasisVector2 = NULL);
+ CosineNodeQueue& treeQueue,
+ VecType* addBasisVector1 = NULL,
+ VecType* addBasisVector2 = NULL);
/**
* Constructs the final basis matrix, after the cosine tree construction.
*
* @param treeQueue Priority queue of cosine nodes.
*/
- void ConstructBasis(CosineNodeQueue& treeQueue);
+ void ConstructBasis(CosineNodeQueue& treeQueue);
/**
* This function splits the cosine node into two children based on the cosines
@@ -153,7 +159,8 @@ class CosineTree
* randomly generated values in the range [0, 1].
*/
void ColumnSamplesLS(std::vector& sampledIndices,
- arma::vec& probabilities, size_t numSamples);
+ VecType& probabilities,
+ size_t numSamples);
/**
* Sample a point from the Length-Squared distribution of the cosine node. The
@@ -175,7 +182,9 @@ class CosineTree
* @param start Starting index of the distribution interval to search in.
* @param end Ending index of the distribution interval to search in.
*/
- size_t BinarySearch(arma::vec& cDistribution, double value, size_t start,
+ size_t BinarySearch(VecType& cDistribution,
+ double value,
+ size_t start,
size_t end);
/**
@@ -185,7 +194,7 @@ class CosineTree
*
* @param cosines Vector to store the cosine values in.
*/
- void CalculateCosines(arma::vec& cosines);
+ void CalculateCosines(VecType& cosines);
/**
* Calculate centroid of the columns present in the node. The calculated
@@ -194,10 +203,10 @@ class CosineTree
void CalculateCentroid();
//! Returns the basis of the constructed subspace.
- void GetFinalBasis(arma::mat& finalBasis) { finalBasis = basis; }
+ void GetFinalBasis(MatType& finalBasis) { finalBasis = basis; }
//! Get pointer to the dataset matrix.
- const arma::mat& GetDataset() const { return *dataset; }
+ const MatType& GetDataset() const { return *dataset; }
//! Get the indices of columns in the node.
std::vector& VectorIndices() { return indices; }
@@ -208,13 +217,13 @@ class CosineTree
double L2Error() const { return l2Error; }
//! Get pointer to the centroid vector.
- arma::vec& Centroid() { return centroid; }
+ VecType& Centroid() { return centroid; }
//! Set the basis vector of the node.
- void BasisVector(arma::vec& bVector) { this->basisVector = bVector; }
+ void BasisVector(VecType& bVector) { this->basisVector = bVector; }
//! Get the basis vector of the node.
- arma::vec& BasisVector() { return basisVector; }
+ VecType& BasisVector() { return basisVector; }
//! Get pointer to the parent node.
CosineTree* Parent() const { return parent; }
@@ -242,11 +251,11 @@ class CosineTree
private:
//! Matrix for which cosine tree is constructed.
- const arma::mat* dataset;
+ const MatType* dataset;
//! Cumulative probability for Monte Carlo error lower bound.
double delta;
//! Subspace basis of the input dataset.
- arma::mat basis;
+ MatType basis;
//! Parent of the node.
CosineTree* parent;
//! Left child of the node.
@@ -256,11 +265,11 @@ class CosineTree
//! Indices of columns of input matrix in the node.
std::vector indices;
//! L2-norm squared of columns in the node.
- arma::vec l2NormsSquared;
+ VecType l2NormsSquared;
//! Centroid of columns of input matrix in the node.
- arma::vec centroid;
+ VecType centroid;
//! Orthonormalized basis vector of the node.
- arma::vec basisVector;
+ VecType basisVector;
//! Index of split point of cosine node.
size_t splitPointIndex;
//! Number of columns of input matrix in the node.
@@ -277,7 +286,9 @@ class CompareCosineNode
{
public:
// Comparison function for construction of priority queue.
- bool operator() (const CosineTree* a, const CosineTree* b) const
+ template
+ bool operator() (const CosineTree* a,
+ const CosineTree* b) const
{
return a->L2Error() < b->L2Error();
}
diff --git a/src/mlpack/core/tree/cosine_tree/cosine_tree_impl.hpp b/src/mlpack/core/tree/cosine_tree/cosine_tree_impl.hpp
index dcc036465f..7ef0393399 100644
--- a/src/mlpack/core/tree/cosine_tree/cosine_tree_impl.hpp
+++ b/src/mlpack/core/tree/cosine_tree/cosine_tree_impl.hpp
@@ -16,7 +16,8 @@
namespace mlpack {
-inline CosineTree::CosineTree(const arma::mat& dataset) :
+template
+inline CosineTree::CosineTree(const MatType& dataset) :
dataset(&dataset),
parent(NULL),
left(NULL),
@@ -32,7 +33,7 @@ inline CosineTree::CosineTree(const arma::mat& dataset) :
for (size_t i = 0; i < numColumns; ++i)
{
indices[i] = i;
- double l2Norm = arma::norm(dataset.col(i), 2);
+ double l2Norm = (double) arma::norm(dataset.col(i), 2);
l2NormsSquared(i) = l2Norm * l2Norm;
}
@@ -45,8 +46,9 @@ inline CosineTree::CosineTree(const arma::mat& dataset) :
splitPointIndex = ColumnSampleLS();
}
-inline CosineTree::CosineTree(CosineTree& parentNode,
- const std::vector& subIndices) :
+template
+inline CosineTree::CosineTree(CosineTree& parentNode,
+ const std::vector& subIndices) :
dataset(&parentNode.GetDataset()),
parent(&parentNode),
left(NULL),
@@ -74,9 +76,10 @@ inline CosineTree::CosineTree(CosineTree& parentNode,
splitPointIndex = ColumnSampleLS();
}
-inline CosineTree::CosineTree(const arma::mat& dataset,
- const double epsilon,
- const double delta) :
+template
+inline CosineTree::CosineTree(const MatType& dataset,
+ const double epsilon,
+ const double delta) :
dataset(&dataset),
delta(delta),
left(NULL),
@@ -84,15 +87,16 @@ inline CosineTree::CosineTree(const arma::mat& dataset,
localDataset(false)
{
// Declare the cosine tree priority queue.
- CosineNodeQueue treeQueue;
+ CosineNodeQueue treeQueue;
CompareCosineNode comp;
// Define root node of the tree and add it to the queue.
CosineTree root(dataset);
- arma::vec tempVector = arma::zeros(dataset.n_rows);
+ VecType tempVector = arma::zeros(dataset.n_rows);
root.L2Error(-1.0); // We don't know what the error is.
root.BasisVector(tempVector);
- treeQueue.push_back(&root); // treeQueue is empty now, so we don't need to call std::push_heap here.
+ treeQueue.push_back(&root);
+ // treeQueue is empty now, so we don't need to call std::push_heap here.
// Initialize Monte Carlo error estimate for comparison.
double monteCarloError = root.FrobNormSquared();
@@ -128,7 +132,7 @@ inline CosineTree::CosineTree(const arma::mat& dataset,
currentRight = currentNode->Right();
// Calculate basis vectors of left and right children.
- arma::vec lBasisVector, rBasisVector;
+ VecType lBasisVector, rBasisVector;
ModifiedGramSchmidt(treeQueue, currentLeft->Centroid(), lBasisVector);
ModifiedGramSchmidt(treeQueue, currentRight->Centroid(), rBasisVector,
@@ -157,9 +161,10 @@ inline CosineTree::CosineTree(const arma::mat& dataset,
}
//! Copy the given tree.
-inline CosineTree::CosineTree(const CosineTree& other) :
+template
+inline CosineTree::CosineTree(const CosineTree& other) :
// Copy matrix, but only if we are the root.
- dataset((other.parent == NULL) ? new arma::mat(*other.dataset) : NULL),
+ dataset((other.parent == NULL) ? new MatType(*other.dataset) : NULL),
delta(other.delta),
parent(NULL),
left(NULL),
@@ -210,7 +215,9 @@ inline CosineTree::CosineTree(const CosineTree& other) :
}
//! Copy assignment operator: copy the given other tree.
-inline CosineTree& CosineTree::operator=(const CosineTree& other)
+template
+inline CosineTree& CosineTree::operator=(
+ const CosineTree& other)
{
// Return if it's the same tree.
if (this == &other)
@@ -224,7 +231,7 @@ inline CosineTree& CosineTree::operator=(const CosineTree& other)
delete right;
// Performing a deep copy of the dataset.
- dataset = (other.parent == NULL) ? new arma::mat(*other.dataset) : NULL;
+ dataset = (other.parent == NULL) ? new MatType(*other.dataset) : NULL;
delta = other.delta;
parent = other.Parent();
@@ -278,7 +285,8 @@ inline CosineTree& CosineTree::operator=(const CosineTree& other)
}
//! Move the given tree.
-inline CosineTree::CosineTree(CosineTree&& other) :
+template
+inline CosineTree::CosineTree(CosineTree&& other) :
dataset(other.dataset),
delta(std::move(other.delta)),
parent(other.parent),
@@ -313,7 +321,8 @@ inline CosineTree::CosineTree(CosineTree&& other) :
}
//! Move assignment operator: take ownership of the given tree.
-inline CosineTree& CosineTree::operator=(CosineTree&& other)
+template
+inline CosineTree& CosineTree::operator=(CosineTree&& other)
{
// Return if it's the same tree.
if (this == &other)
@@ -360,7 +369,8 @@ inline CosineTree& CosineTree::operator=(CosineTree&& other)
return *this;
}
-inline CosineTree::~CosineTree()
+template
+inline CosineTree::~CosineTree()
{
if (localDataset)
delete dataset;
@@ -370,17 +380,19 @@ inline CosineTree::~CosineTree()
delete right;
}
-inline void CosineTree::ModifiedGramSchmidt(CosineNodeQueue& treeQueue,
- arma::vec& centroid,
- arma::vec& newBasisVector,
- arma::vec* addBasisVector)
+template
+inline void CosineTree::ModifiedGramSchmidt(
+ CosineNodeQueue& treeQueue,
+ typename CosineTree::VecType& centroid,
+ typename CosineTree::VecType& newBasisVector,
+ typename CosineTree::VecType* addBasisVector)
{
// Set new basis vector to centroid.
newBasisVector = centroid;
// Variables for iterating throught the priority queue.
CosineTree *currentNode;
- CosineNodeQueue::const_iterator i = treeQueue.cbegin();
+ typename CosineNodeQueue::const_iterator i = treeQueue.cbegin();
// For every vector in the current basis, remove its projection from the
// centroid.
@@ -388,14 +400,14 @@ inline void CosineTree::ModifiedGramSchmidt(CosineNodeQueue& treeQueue,
{
currentNode = *i;
- double projection = dot(currentNode->BasisVector(), centroid);
+ double projection = (double) dot(currentNode->BasisVector(), centroid);
newBasisVector -= projection * currentNode->BasisVector();
}
// If additional basis vector is passed, take it into account.
if (addBasisVector)
{
- double projection = dot(*addBasisVector, centroid);
+ double projection = (double) dot(*addBasisVector, centroid);
newBasisVector -= *addBasisVector * projection;
}
@@ -404,13 +416,15 @@ inline void CosineTree::ModifiedGramSchmidt(CosineNodeQueue& treeQueue,
newBasisVector /= arma::norm(newBasisVector, 2);
}
-inline double CosineTree::MonteCarloError(CosineTree* node,
- CosineNodeQueue& treeQueue,
- arma::vec* addBasisVector1,
- arma::vec* addBasisVector2)
+template
+inline double CosineTree::MonteCarloError(
+ CosineTree* node,
+ CosineNodeQueue& treeQueue,
+ typename CosineTree::VecType* addBasisVector1,
+ typename CosineTree::VecType* addBasisVector2)
{
std::vector sampledIndices;
- arma::vec probabilities;
+ VecType probabilities;
// Sample O(log m) points from the input node's distribution.
// 'm' is the number of columns present in the node.
@@ -418,10 +432,10 @@ inline double CosineTree::MonteCarloError(CosineTree* node,
node->ColumnSamplesLS(sampledIndices, probabilities, numSamples);
// Get pointer to the original dataset.
- const arma::mat& dataset = node->GetDataset();
+ const MatType& dataset = node->GetDataset();
// Initialize weighted projection magnitudes as zeros.
- arma::vec weightedMagnitudes;
+ VecType weightedMagnitudes;
weightedMagnitudes.zeros(numSamples);
// Set size of projection vector, depending on whether additional basis
@@ -436,11 +450,11 @@ inline double CosineTree::MonteCarloError(CosineTree* node,
for (size_t i = 0; i < numSamples; ++i)
{
// Initialize projection as a vector of zeros.
- arma::vec projection;
+ VecType projection;
projection.zeros(projectionSize);
CosineTree *currentNode;
- CosineNodeQueue::const_iterator j = treeQueue.cbegin();
+ typename CosineNodeQueue::const_iterator j = treeQueue.cbegin();
size_t k = 0;
// Compute the projection of the sampled vector onto the existing subspace.
@@ -488,14 +502,16 @@ inline double CosineTree::MonteCarloError(CosineTree* node,
return (node->FrobNormSquared() - lowerBound);
}
-inline void CosineTree::ConstructBasis(CosineNodeQueue& treeQueue)
+template
+inline void CosineTree::ConstructBasis(
+ CosineNodeQueue& treeQueue)
{
// Initialize basis as matrix of zeros.
basis.zeros(dataset->n_rows, treeQueue.size());
// Variables for iterating through the priority queue.
CosineTree *currentNode;
- CosineNodeQueue::const_iterator i = treeQueue.cbegin();
+ typename CosineNodeQueue::const_iterator i = treeQueue.cbegin();
// Transfer basis vectors from the queue to the basis matrix.
size_t j = 0;
@@ -506,7 +522,8 @@ inline void CosineTree::ConstructBasis(CosineNodeQueue& treeQueue)
}
}
-inline void CosineTree::CosineNodeSplit()
+template
+inline void CosineTree::CosineNodeSplit()
{
// If less than two points, splitting does not make sense---there is nothing
// to split.
@@ -514,7 +531,7 @@ inline void CosineTree::CosineNodeSplit()
return;
// Calculate cosines with respect to the splitting point.
- arma::vec cosines;
+ VecType cosines;
CalculateCosines(cosines);
// Compute maximum and minimum cosine values.
@@ -543,12 +560,14 @@ inline void CosineTree::CosineNodeSplit()
right = new CosineTree(*this, rightIndices);
}
-inline void CosineTree::ColumnSamplesLS(std::vector& sampledIndices,
- arma::vec& probabilities,
- size_t numSamples)
+template
+inline void CosineTree::ColumnSamplesLS(
+ std::vector& sampledIndices,
+ typename CosineTree::VecType& probabilities,
+ size_t numSamples)
{
// Initialize the cumulative distribution vector size.
- arma::vec cDistribution;
+ VecType cDistribution;
cDistribution.zeros(numColumns + 1);
// Calculate cumulative length-squared distribution for the node.
@@ -575,7 +594,8 @@ inline void CosineTree::ColumnSamplesLS(std::vector& sampledIndices,
}
}
-inline size_t CosineTree::ColumnSampleLS()
+template
+inline size_t CosineTree::ColumnSampleLS()
{
// If only one element is present, there can only be one sample.
if (numColumns < 2)
@@ -584,7 +604,7 @@ inline size_t CosineTree::ColumnSampleLS()
}
// Initialize the cumulative distribution vector size.
- arma::vec cDistribution;
+ VecType cDistribution;
cDistribution.zeros(numColumns + 1);
// Calculate cumulative length-squared distribution for the node.
@@ -602,10 +622,12 @@ inline size_t CosineTree::ColumnSampleLS()
return BinarySearch(cDistribution, randValue, start, end);
}
-inline size_t CosineTree::BinarySearch(arma::vec& cDistribution,
- double value,
- size_t start,
- size_t end)
+template
+inline size_t CosineTree::BinarySearch(
+ typename CosineTree::VecType& cDistribution,
+ double value,
+ size_t start,
+ size_t end)
{
size_t pivot = (start + end) / 2;
@@ -630,7 +652,9 @@ inline size_t CosineTree::BinarySearch(arma::vec& cDistribution,
}
}
-inline void CosineTree::CalculateCosines(arma::vec& cosines)
+template
+inline void CosineTree::CalculateCosines(
+ typename CosineTree::VecType& cosines)
{
// Initialize cosine vector as a vector of zeros.
cosines.zeros(numColumns);
@@ -652,7 +676,8 @@ inline void CosineTree::CalculateCosines(arma::vec& cosines)
}
}
-inline void CosineTree::CalculateCentroid()
+template
+inline void CosineTree::CalculateCentroid()
{
// Initialize centroid as vector of zeros.
centroid.zeros(dataset->n_rows);
diff --git a/src/mlpack/core/util/arma_traits.hpp b/src/mlpack/core/util/arma_traits.hpp
index 923ea3da72..2dd360761b 100644
--- a/src/mlpack/core/util/arma_traits.hpp
+++ b/src/mlpack/core/util/arma_traits.hpp
@@ -254,4 +254,49 @@ struct GetSparseMatType>
typedef arma::SpMat type;
};
+// Get whether or not the given type is a base matrix type (e.g. not an
+// expression).
+
+template
+struct IsBaseMatType
+{
+ constexpr static bool value = false;
+};
+
+template
+struct IsBaseMatType>
+{
+ constexpr static bool value = true;
+};
+
+template
+struct IsBaseMatType>
+{
+ constexpr static bool value = true;
+};
+
+template
+struct IsBaseMatType>
+{
+ constexpr static bool value = true;
+};
+
+template
+struct IsBaseMatType>
+{
+ constexpr static bool value = true;
+};
+
+template
+struct IsBaseMatType>
+{
+ constexpr static bool value = true;
+};
+
+template
+struct IsBaseMatType>
+{
+ constexpr static bool value = true;
+};
+
#endif
diff --git a/src/mlpack/methods/ann/ann.hpp b/src/mlpack/methods/ann/ann.hpp
index 9c20ec46a3..e2e2a468d3 100644
--- a/src/mlpack/methods/ann/ann.hpp
+++ b/src/mlpack/methods/ann/ann.hpp
@@ -17,7 +17,6 @@
#define MLPACK_METHODS_ANN_ANN_HPP
#include "forward_decls.hpp"
-#include "make_alias.hpp"
#include "activation_functions/activation_functions.hpp"
#include "augmented/augmented.hpp"
diff --git a/src/mlpack/methods/ann/ffn_impl.hpp b/src/mlpack/methods/ann/ffn_impl.hpp
index aaf7545606..f95eec3ea1 100644
--- a/src/mlpack/methods/ann/ffn_impl.hpp
+++ b/src/mlpack/methods/ann/ffn_impl.hpp
@@ -15,8 +15,6 @@
// In case it hasn't been included yet.
#include "ffn.hpp"
-#include "make_alias.hpp"
-
namespace mlpack {
template& InputDimensions() const { return inputDimensions; }
diff --git a/src/mlpack/methods/ann/layer/multi_layer.hpp b/src/mlpack/methods/ann/layer/multi_layer.hpp
index 991bafa31b..bad1ef6038 100644
--- a/src/mlpack/methods/ann/layer/multi_layer.hpp
+++ b/src/mlpack/methods/ann/layer/multi_layer.hpp
@@ -12,7 +12,6 @@
#ifndef MLPACK_METHODS_ANN_LAYER_MULTI_LAYER_HPP
#define MLPACK_METHODS_ANN_LAYER_MULTI_LAYER_HPP
-#include "../make_alias.hpp"
#include "layer.hpp"
namespace mlpack {
diff --git a/src/mlpack/methods/ann/make_alias.hpp b/src/mlpack/methods/ann/make_alias.hpp
deleted file mode 100644
index a0de66a004..0000000000
--- a/src/mlpack/methods/ann/make_alias.hpp
+++ /dev/null
@@ -1,58 +0,0 @@
-/**
- * @file make_alias.hpp
- * @author Ryan Curtin
- *
- * Implementation of `MakeAlias()`, a utility function. This is meant to be
- * used in `SetWeights()` calls in various layers, to wrap internal weight
- * objects as aliases around the given memory pointers.
- *
- * 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_ANN_MAKE_ALIAS_HPP
-#define MLPACK_METHODS_ANN_MAKE_ALIAS_HPP
-
-#include
-
-namespace mlpack {
-
-/**
- * Reconstruct `m` as an alias around the memory `newMem`, with size `numRows` x
- * `numCols`.
- */
-template
-void MakeAlias(MatType& m,
- typename MatType::elem_type* newMem,
- const size_t numRows,
- const size_t numCols)
-{
- // We use placement new to reinitialize the object, since the copy and move
- // assignment operators in Armadillo will end up copying memory instead of
- // making an alias.
- m.~MatType();
- new (&m) MatType(newMem, numRows, numCols, false, true);
-}
-
-/**
- * Reconstruct `c` as an alias around the memory` newMem`, with size `numRows` x
- * `numCols` x `numSlices`.
- */
-template
-void MakeAlias(CubeType& c,
- typename CubeType::elem_type* newMem,
- const size_t numRows,
- const size_t numCols,
- const size_t numSlices)
-{
- // We use placement new to reinitialize the object, since the copy and move
- // assignment operators in Armadillo will end up copying memory instead of
- // making an alias.
- c.~CubeType();
- new (&c) CubeType(newMem, numRows, numCols, numSlices, false, true);
-}
-
-} // namespace mlpack
-
-#endif
diff --git a/src/mlpack/methods/bias_svd/bias_svd_function_impl.hpp b/src/mlpack/methods/bias_svd/bias_svd_function_impl.hpp
index 729b3f3c98..735b2ec66a 100644
--- a/src/mlpack/methods/bias_svd/bias_svd_function_impl.hpp
+++ b/src/mlpack/methods/bias_svd/bias_svd_function_impl.hpp
@@ -19,13 +19,14 @@
namespace mlpack {
template
-BiasSVDFunction::BiasSVDFunction(const MatType& data,
+BiasSVDFunction::BiasSVDFunction(const MatType& dataIn,
const size_t rank,
const double lambda) :
- data(MakeAlias(const_cast(data), false)),
rank(rank),
lambda(lambda)
{
+ MakeAlias(data, dataIn, dataIn.n_rows, dataIn.n_cols, false);
+
// Number of users and items in the data.
numUsers = max(data.row(0)) + 1;
numItems = max(data.row(1)) + 1;
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 add30fb990..ad32145c06 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
@@ -69,10 +69,11 @@ class RandomizedBlockKrylovSVD
* @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,
+ template
+ RandomizedBlockKrylovSVD(const InMatType& data,
+ MatType& u,
+ VecType& s,
+ MatType& v,
const size_t maxIterations = 2,
const size_t rank = 0,
const size_t blockSize = 0);
@@ -97,10 +98,11 @@ class RandomizedBlockKrylovSVD
* @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,
+ template
+ void Apply(const InMatType& data,
+ MatType& u,
+ VecType& s,
+ MatType& v,
const size_t rank);
//! Get the number of iterations for the power method.
diff --git a/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd_impl.hpp b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd_impl.hpp
index 982321ecbf..567dd2ccbf 100644
--- a/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd_impl.hpp
+++ b/src/mlpack/methods/block_krylov_svd/randomized_block_krylov_svd_impl.hpp
@@ -16,11 +16,12 @@
namespace mlpack {
+template
inline RandomizedBlockKrylovSVD::RandomizedBlockKrylovSVD(
- const arma::mat& data,
- arma::mat& u,
- arma::vec& s,
- arma::mat& v,
+ const InMatType& data,
+ MatType& u,
+ VecType& s,
+ MatType& v,
const size_t maxIterations,
const size_t rank,
const size_t blockSize) :
@@ -46,42 +47,46 @@ inline RandomizedBlockKrylovSVD::RandomizedBlockKrylovSVD(
/* Nothing to do here */
}
-inline void RandomizedBlockKrylovSVD::Apply(const arma::mat& data,
- arma::mat& u,
- arma::vec& s,
- arma::mat& v,
+template
+inline void RandomizedBlockKrylovSVD::Apply(const InMatType& data,
+ MatType& u,
+ VecType& s,
+ MatType& v,
const size_t rank)
{
- arma::mat Q, R, block, blockIteration;
+ MatType Q, R, block, blockIteration;
if (blockSize == 0)
{
- blockSize = rank + 10;
+ // The block size cannot be greater than the number of points in the
+ // dataset or the dimensionality of the dataset.
+ blockSize = std::min((size_t) data.n_rows, std::min((size_t) data.n_cols,
+ rank + 10));
}
// Random block initialization.
- arma::mat G = arma::randn(data.n_cols, blockSize);
+ MatType G = arma::randn(data.n_cols, blockSize);
// Construct and orthonormalize Krylov subspace.
- arma::mat K(data.n_rows, blockSize * (maxIterations + 1));
+ MatType 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);
+ // (K matrix). Doing so avoids an unnecessary copy in upcoming step.
+ MakeAlias(block, K.memptr(), data.n_rows, blockSize, 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, block.n_rows,
- block.n_cols, false, false);
+ MakeAlias(blockIteration, K.memptr() + blockOffset, block.n_rows,
+ block.n_cols, false);
arma::qr_econ(blockIteration, R, data * (data.t() * block));
// Update working matrix for the next iteration.
- block = arma::mat(K.memptr() + blockOffset, block.n_rows, block.n_cols,
- false, false);
+ MakeAlias(block, K.memptr() + blockOffset, block.n_rows, block.n_cols,
+ false);
}
arma::qr_econ(Q, R, K);
diff --git a/src/mlpack/methods/cf/decomposition_policies/quic_svd_method.hpp b/src/mlpack/methods/cf/decomposition_policies/quic_svd_method.hpp
index 775a7c8803..a8b6cdbe39 100644
--- a/src/mlpack/methods/cf/decomposition_policies/quic_svd_method.hpp
+++ b/src/mlpack/methods/cf/decomposition_policies/quic_svd_method.hpp
@@ -74,7 +74,7 @@ class QUIC_SVDPolicy
arma::mat data(cleanedData);
// Do singular value decomposition using the quic SVD algorithm.
- QUIC_SVD quicsvd;
+ QUIC_SVD<> quicsvd;
quicsvd.Apply(data, w, h, sigma);
// Sigma matrix is multiplied to w.
diff --git a/src/mlpack/methods/linear_svm/linear_svm_function_impl.hpp b/src/mlpack/methods/linear_svm/linear_svm_function_impl.hpp
index beea27d0d1..86bffc0dc4 100644
--- a/src/mlpack/methods/linear_svm/linear_svm_function_impl.hpp
+++ b/src/mlpack/methods/linear_svm/linear_svm_function_impl.hpp
@@ -24,18 +24,19 @@ namespace mlpack {
template
LinearSVMFunction::LinearSVMFunction(
- const MatType& dataset,
+ const MatType& datasetIn,
const arma::Row& labels,
const size_t numClasses,
const double lambda,
const double delta,
const bool fitIntercept) :
- dataset(MakeAlias(const_cast(dataset), false)),
numClasses(numClasses),
lambda(lambda),
delta(delta),
fitIntercept(fitIntercept)
{
+ MakeAlias(dataset, datasetIn, datasetIn.n_rows, datasetIn.n_cols, false);
+
InitializeWeights(initialPoint, dataset.n_rows, numClasses, fitIntercept);
initialPoint *= 0.005;
diff --git a/src/mlpack/methods/lmnn/lmnn_function_impl.hpp b/src/mlpack/methods/lmnn/lmnn_function_impl.hpp
index 73a9baa497..1197ec49f0 100644
--- a/src/mlpack/methods/lmnn/lmnn_function_impl.hpp
+++ b/src/mlpack/methods/lmnn/lmnn_function_impl.hpp
@@ -19,23 +19,24 @@
namespace mlpack {
template
-LMNNFunction::LMNNFunction(const arma::mat& dataset,
- const arma::Row& labels,
+LMNNFunction::LMNNFunction(const arma::mat& datasetIn,
+ const arma::Row& labelsIn,
size_t k,
double regularization,
size_t range,
MetricType metric) :
- dataset(MakeAlias(const_cast(dataset), false)),
- labels(MakeAlias(const_cast&>(labels), false)),
k(k),
metric(metric),
regularization(regularization),
iteration(0),
range(range),
- constraint(dataset, labels, k),
- points(dataset.n_cols),
+ constraint(datasetIn, labelsIn, k),
+ points(datasetIn.n_cols),
impBounds(false)
{
+ MakeAlias(dataset, datasetIn, datasetIn.n_rows, datasetIn.n_cols, false);
+ MakeAlias(labels, labelsIn, labelsIn.n_rows, labelsIn.n_cols, false);
+
// Initialize the initial learning point.
initialPoint.eye(dataset.n_rows, dataset.n_rows);
// Initialize transformed dataset to base dataset.
diff --git a/src/mlpack/methods/logistic_regression/logistic_regression_function_impl.hpp b/src/mlpack/methods/logistic_regression/logistic_regression_function_impl.hpp
index 4e861ac287..3b7c2282b5 100644
--- a/src/mlpack/methods/logistic_regression/logistic_regression_function_impl.hpp
+++ b/src/mlpack/methods/logistic_regression/logistic_regression_function_impl.hpp
@@ -21,15 +21,17 @@ namespace mlpack {
template
LogisticRegressionFunction::LogisticRegressionFunction(
- const MatType& predictors,
- const arma::Row& responses,
+ const MatType& predictorsIn,
+ const arma::Row& responsesIn,
const double lambda) :
- // We promise to be well-behaved... the elements won't be modified.
- predictors(MakeAlias(const_cast(predictors), false)),
- responses(MakeAlias(const_cast&>(responses),
- false)),
lambda(lambda)
{
+ // We promise to be well-behaved... the elements won't be modified.
+ MakeAlias(this->predictors, predictorsIn, predictorsIn.n_rows,
+ predictorsIn.n_cols, false);
+ MakeAlias(this->responses, responsesIn, responsesIn.n_rows,
+ responsesIn.n_cols, false);
+
// Sanity check.
if (responses.n_elem != predictors.n_cols)
{
diff --git a/src/mlpack/methods/nca/nca_softmax_error_function_impl.hpp b/src/mlpack/methods/nca/nca_softmax_error_function_impl.hpp
index c83f3afb7f..1ce40fb853 100644
--- a/src/mlpack/methods/nca/nca_softmax_error_function_impl.hpp
+++ b/src/mlpack/methods/nca/nca_softmax_error_function_impl.hpp
@@ -22,14 +22,15 @@ namespace mlpack {
// Initialize with the given kernel.
template
SoftmaxErrorFunction::SoftmaxErrorFunction(
- const arma::mat& dataset,
- const arma::Row& labels,
+ const arma::mat& datasetIn,
+ const arma::Row& labelsIn,
MetricType metric) :
- dataset(MakeAlias(const_cast(dataset), false)),
- labels(MakeAlias(const_cast&>(labels), false)),
metric(metric),
precalculated(false)
-{ /* nothing to do */ }
+{
+ MakeAlias(dataset, datasetIn, datasetIn.n_rows, datasetIn.n_cols, false);
+ MakeAlias(labels, labelsIn, labelsIn.n_rows, labelsIn.n_cols, false);
+}
//! Shuffle the dataset.
template
diff --git a/src/mlpack/methods/pca/decomposition_policies/exact_svd_method.hpp b/src/mlpack/methods/pca/decomposition_policies/exact_svd_method.hpp
index 4f42d23eac..cd37fee403 100644
--- a/src/mlpack/methods/pca/decomposition_policies/exact_svd_method.hpp
+++ b/src/mlpack/methods/pca/decomposition_policies/exact_svd_method.hpp
@@ -37,19 +37,20 @@ class ExactSVDPolicy
* @param eigvec Matrix to put eigenvectors (loadings) into.
* @param * (rank) Rank of the decomposition.
*/
- void Apply(const arma::mat& data,
- const arma::mat& centeredData,
- arma::mat& transformedData,
- arma::vec& eigVal,
- arma::mat& eigvec,
+ template
+ void Apply(const InMatType& /* data */,
+ const MatType& centeredData,
+ MatType& transformedData,
+ VecType& eigVal,
+ MatType& eigvec,
const size_t /* rank */)
{
// This matrix will store the right singular vectors; we do not need them.
- arma::mat v;
+ MatType v;
// Do singular value decomposition. Use the economical singular value
// decomposition if the columns are much larger than the rows.
- if (data.n_rows < data.n_cols)
+ if (centeredData.n_rows < centeredData.n_cols)
{
// Do economical singular value decomposition and compute only the left
// singular vectors.
@@ -63,7 +64,7 @@ class ExactSVDPolicy
// Now we must square the singular values to get the eigenvalues.
// In addition we must divide by the number of points, because the
// covariance matrix is X * X' / (N - 1).
- eigVal %= eigVal / (data.n_cols - 1);
+ eigVal %= eigVal / (centeredData.n_cols - 1);
// Project the samples to the principals.
transformedData = trans(eigvec) * centeredData;
diff --git a/src/mlpack/methods/pca/decomposition_policies/quic_svd_method.hpp b/src/mlpack/methods/pca/decomposition_policies/quic_svd_method.hpp
index 6a84ccbef6..240cbe3e1c 100644
--- a/src/mlpack/methods/pca/decomposition_policies/quic_svd_method.hpp
+++ b/src/mlpack/methods/pca/decomposition_policies/quic_svd_method.hpp
@@ -49,23 +49,24 @@ class QUICSVDPolicy
* @param eigvec Matrix to put eigenvectors (loadings) into.
* @param * (rank) Rank of the decomposition.
*/
- void Apply(const arma::mat& data,
- const arma::mat& centeredData,
- arma::mat& transformedData,
- arma::vec& eigVal,
- arma::mat& eigvec,
+ template
+ void Apply(const InMatType& /* data */,
+ const MatType& centeredData,
+ MatType& transformedData,
+ VecType& eigVal,
+ MatType& eigvec,
const size_t /* rank */)
{
// This matrix will store the right singular vectors; we do not need them.
- arma::mat v, sigma;
+ MatType v, sigma;
// Do singular value decomposition using the QUIC-SVD algorithm.
- QUIC_SVD quicsvd(centeredData, eigvec, v, sigma, epsilon, delta);
+ QUIC_SVD quicsvd(centeredData, eigvec, v, sigma, epsilon, delta);
// Now we must square the singular values to get the eigenvalues.
// In addition we must divide by the number of points, because the
// covariance matrix is X * X' / (N - 1).
- eigVal = pow(arma::diagvec(sigma), 2) / (data.n_cols - 1);
+ eigVal = pow(arma::diagvec(sigma), 2) / (centeredData.n_cols - 1);
// Project the samples to the principals.
transformedData = trans(eigvec) * centeredData;
diff --git a/src/mlpack/methods/pca/decomposition_policies/randomized_block_krylov_method.hpp b/src/mlpack/methods/pca/decomposition_policies/randomized_block_krylov_method.hpp
index 7836b72c8a..64aad7b4b8 100644
--- a/src/mlpack/methods/pca/decomposition_policies/randomized_block_krylov_method.hpp
+++ b/src/mlpack/methods/pca/decomposition_policies/randomized_block_krylov_method.hpp
@@ -52,15 +52,16 @@ class RandomizedBlockKrylovSVDPolicy
* @param eigvec Matrix to put eigenvectors (loadings) into.
* @param rank Rank of the decomposition.
*/
- void Apply(const arma::mat& data,
- const arma::mat& centeredData,
- arma::mat& transformedData,
- arma::vec& eigVal,
- arma::mat& eigvec,
+ template
+ void Apply(const InMatType& /* data */,
+ const MatType& centeredData,
+ MatType& transformedData,
+ VecType& eigVal,
+ MatType& eigvec,
const size_t rank)
{
// This matrix will store the right singular vectors; we do not need them.
- arma::mat v;
+ MatType v;
// Do singular value decomposition using the randomized block krylov SVD
// algorithm.
@@ -70,7 +71,7 @@ class RandomizedBlockKrylovSVDPolicy
// Now we must square the singular values to get the eigenvalues.
// In addition we must divide by the number of points, because the
// covariance matrix is X * X' / (N - 1).
- eigVal %= eigVal / (data.n_cols - 1);
+ eigVal %= eigVal / (centeredData.n_cols - 1);
// Project the samples to the principals.
transformedData = trans(eigvec) * centeredData;
diff --git a/src/mlpack/methods/pca/decomposition_policies/randomized_svd_method.hpp b/src/mlpack/methods/pca/decomposition_policies/randomized_svd_method.hpp
index ceb2fab0cd..cec751717e 100644
--- a/src/mlpack/methods/pca/decomposition_policies/randomized_svd_method.hpp
+++ b/src/mlpack/methods/pca/decomposition_policies/randomized_svd_method.hpp
@@ -53,15 +53,16 @@ class RandomizedSVDPCAPolicy
* @param eigvec Matrix to put eigenvectors (loadings) into.
* @param rank Rank of the decomposition.
*/
- void Apply(const arma::mat& data,
- const arma::mat& centeredData,
- arma::mat& transformedData,
- arma::vec& eigVal,
- arma::mat& eigvec,
+ template
+ void Apply(const InMatType& data,
+ const MatType& centeredData,
+ MatType& transformedData,
+ VecType& eigVal,
+ MatType& eigvec,
const size_t rank)
{
// This matrix will store the right singular vectors; we do not need them.
- arma::mat v;
+ MatType v;
// Do singular value decomposition using the randomized SVD algorithm.
RandomizedSVD rsvd(iteratedPower, maxIterations);
@@ -70,7 +71,7 @@ class RandomizedSVDPCAPolicy
// Now we must square the singular values to get the eigenvalues.
// In addition we must divide by the number of points, because the
// covariance matrix is X * X' / (N - 1).
- eigVal %= eigVal / (data.n_cols - 1);
+ eigVal %= eigVal / (centeredData.n_cols - 1);
// Project the samples to the principals.
transformedData = trans(eigvec) * centeredData;
diff --git a/src/mlpack/methods/pca/pca.hpp b/src/mlpack/methods/pca/pca.hpp
index 295b5a188c..bf92350d45 100644
--- a/src/mlpack/methods/pca/pca.hpp
+++ b/src/mlpack/methods/pca/pca.hpp
@@ -52,10 +52,13 @@ class PCA
* @param eigVal Vector to put eigenvalues into.
* @param eigvec Matrix to put eigenvectors (loadings) into.
*/
- void Apply(const arma::mat& data,
- arma::mat& transformedData,
- arma::vec& eigVal,
- arma::mat& eigvec);
+ template
+ void Apply(const MatType& data,
+ OutMatType& transformedData,
+ VecType& eigVal,
+ OutMatType& eigvec);
/**
* Apply Principal Component Analysis to the provided data set. It is safe
@@ -65,17 +68,21 @@ class PCA
* @param transformedData Matrix to store results of PCA in.
* @param eigVal Vector to put eigenvalues into.
*/
- void Apply(const arma::mat& data,
- arma::mat& transformedData,
- arma::vec& eigVal);
+ template
+ void Apply(const MatType& data,
+ OutMatType& transformedData,
+ VecType& eigVal);
/**
* Apply Principal Component Analysis to the provided data set. It is safe
* to pass the same matrix reference for both data and transformedData.
* @param data Data matrix.
* @param transformedData Matrix to store results of PCA in.
*/
- void Apply(const arma::mat& data,
- arma::mat& transformedData);
+ template
+ void Apply(const MatType& data,
+ OutMatType& transformedData);
/**
* Use PCA for dimensionality reduction on the given dataset. This will save
@@ -88,14 +95,43 @@ class PCA
* @param newDimension New dimension of the data.
* @return Amount of the variance of the data retained (between 0 and 1).
*/
- double Apply(arma::mat& data, const size_t newDimension);
+ template
+ double Apply(MatType& data, const size_t newDimension);
//! This overload is here to make sure int gets casted right to size_t.
- inline double Apply(arma::mat& data, const int newDimension)
+ template
+ inline double Apply(MatType& data, const int newDimension)
{
return Apply(data, size_t(newDimension));
}
+ /**
+ * Use PCA for dimensionality reduction on the given dataset. This will save
+ * the newDimension largest principal components of the data and remove the
+ * rest, storing the result in `transformedData`. The return value is the
+ * amount of variance of the data that is retained; this is a value between 0
+ * and 1. For instance, a value of 0.9 indicates that 90% of the variance
+ * present in the data was retained.
+ *
+ * @param data Data matrix.
+ * @param transformedData Output matrix to store transformed data in.
+ * @param newDimension New dimension of the data.
+ * @return Amount of the variance of the data retained (between 0 and 1).
+ */
+ template
+ double Apply(const MatType& data,
+ OutMatType& transformedData,
+ const size_t newDimension);
+
+ //! This overload is here to make sure int gets casted right to size_t.
+ template
+ inline double Apply(const MatType& data,
+ OutMatType& transformedData,
+ const int newDimension)
+ {
+ return Apply(data, transformedData, size_t(newDimension));
+ }
+
/**
* Use PCA for dimensionality reduction on the given dataset. This will save
* as many dimensions as necessary to retain at least the given amount of
@@ -111,7 +147,28 @@ class PCA
* between 0 and 1.
* @return Actual amount of variance retained (between 0 and 1).
*/
- double Apply(arma::mat& data, const double varRetained);
+ template
+ double Apply(MatType& data, const double varRetained);
+
+ /**
+ * Use PCA for dimensionality reduction on the given dataset. This will save
+ * as many dimensions as necessary to retain at least the given amount of
+ * variance (specified by parameter varRetained). The amount should be
+ * between 0 and 1; if the amount is 0, then only 1 dimension will be
+ * retained. If the amount is 1, then all dimensions will be retained.
+ *
+ * The method returns the actual amount of variance retained, which will
+ * always be greater than or equal to the varRetained parameter.
+ *
+ * @param data Data matrix.
+ * @param varRetained Lower bound on amount of variance to retain; should be
+ * between 0 and 1.
+ * @return Actual amount of variance retained (between 0 and 1).
+ */
+ template
+ double Apply(const MatType& data,
+ OutMatType& transformedData,
+ const double varRetained);
//! Get whether or not this PCA object will scale (by standard deviation)
//! the data when PCA is performed.
@@ -122,14 +179,15 @@ class PCA
private:
//! Scaling the data is when we reduce the variance of each dimension to 1.
- void ScaleData(arma::mat& centeredData)
+ template
+ void ScaleData(MatType& centeredData)
{
if (scaleData)
{
// Scaling the data is when we reduce the variance of each dimension
// to 1. We do this by dividing each dimension by its standard
// deviation.
- arma::vec stdDev = arma::stddev(
+ arma::Col stdDev = arma::stddev(
centeredData, 0, 1 /* for each dimension */);
// If there are any zeroes, make them very small.
diff --git a/src/mlpack/methods/pca/pca_impl.hpp b/src/mlpack/methods/pca/pca_impl.hpp
index 05db905443..865c968300 100644
--- a/src/mlpack/methods/pca/pca_impl.hpp
+++ b/src/mlpack/methods/pca/pca_impl.hpp
@@ -12,7 +12,6 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
-
#ifndef MLPACK_METHODS_PCA_PCA_IMPL_HPP
#define MLPACK_METHODS_PCA_PCA_IMPL_HPP
@@ -37,19 +36,31 @@ PCA::PCA(
* @param eigvec - PCA Loadings/Coeffs/EigenVectors
*/
template
-void PCA::Apply(const arma::mat& data,
- arma::mat& transformedData,
- arma::vec& eigVal,
- arma::mat& eigvec)
+template
+void PCA::Apply(const MatType& data,
+ OutMatType& transformedData,
+ VecType& eigVal,
+ OutMatType& eigvec)
{
+ // Sanity checks on input types.
+ static_assert(IsBaseMatType::value,
+ "PCA::Apply(): transformedData must be a matrix type!");
+ static_assert(IsBaseMatType::value,
+ "PCA::Apply(): eigVal must be a vector type!");
+ static_assert(std::is_same::value,
+ "PCA::Apply(): data and transformedData must have the same element "
+ "types!");
+
// Center the data into a temporary matrix.
- arma::mat centeredData = data.each_col() - arma::mean(data, 1);
+ OutMatType centeredData = arma::conv_to::from(data);
+ centeredData.each_col() -= arma::mean(centeredData, 1);
// Scale the data if the user asked for it.
ScaleData(centeredData);
decomposition.Apply(data, centeredData, transformedData, eigVal, eigvec,
- data.n_rows);
+ centeredData.n_rows);
}
/**
@@ -60,11 +71,22 @@ void PCA::Apply(const arma::mat& data,
* @param eigVal - contains eigen values in a column vector
*/
template
-void PCA::Apply(const arma::mat& data,
- arma::mat& transformedData,
- arma::vec& eigVal)
+template
+void PCA::Apply(const MatType& data,
+ OutMatType& transformedData,
+ VecType& eigVal)
{
- arma::mat eigvec;
+ // Sanity checks on input types.
+ static_assert(IsBaseMatType::value,
+ "PCA::Apply(): transformedData must be a matrix type!");
+ static_assert(IsBaseMatType::value,
+ "PCA::Apply(): eigVal must be a vector type!");
+ static_assert(std::is_same::value,
+ "PCA::Apply(): data and transformedData must have the same element "
+ "types!");
+
+ OutMatType eigvec;
Apply(data, transformedData, eigVal, eigvec);
}
@@ -75,11 +97,24 @@ void PCA::Apply(const arma::mat& data,
* @param transformedData Data with PCA applied.
*/
template
-void PCA::Apply(const arma::mat& data,
- arma::mat& transformedData)
+template
+void PCA::Apply(const MatType& data,
+ OutMatType& transformedData)
{
- arma::mat eigvec;
- arma::vec eigVal;
+ // Sanity checks on input types.
+ static_assert(IsBaseMatType::value,
+ "PCA::Apply(): transformedData must be a matrix type!");
+ static_assert(std::is_same::value,
+ "PCA::Apply(): data and transformedData must have the same element "
+ "types!");
+
+ // It's possible a user didn't pass in a matrix but instead an expression, but
+ // we need a type that we can store.
+ typedef typename GetDenseColType::type BaseColType;
+
+ OutMatType eigvec;
+ BaseColType eigVal;
Apply(data, transformedData, eigVal, eigvec);
}
@@ -95,32 +130,58 @@ void PCA::Apply(const arma::mat& data,
* @return Amount of the variance of the data retained (between 0 and 1).
*/
template
-double PCA::Apply(arma::mat& data,
+template
+double PCA::Apply(MatType& data,
+ const size_t newDimension)
+{
+ return Apply(data, data, newDimension);
+}
+
+template
+template
+double PCA::Apply(const MatType& data,
+ OutMatType& transformedData,
const size_t newDimension)
{
// Parameter validation.
if (newDimension == 0)
- Log::Fatal << "PCA::Apply(): newDimension (" << newDimension << ") cannot "
- << "be zero!" << std::endl;
- if (newDimension > data.n_rows)
- Log::Fatal << "PCA::Apply(): newDimension (" << newDimension << ") cannot "
- << "be greater than the existing dimensionality of the data ("
- << data.n_rows << ")!" << std::endl;
+ {
+ std::ostringstream oss;
+ oss << "PCA::Apply(): newDimension (" << newDimension << ") cannot be "
+ << "zero!";
+ throw std::invalid_argument(oss.str());
+ }
- arma::mat eigvec;
- arma::vec eigVal;
+ typedef typename GetDenseMatType::type BaseMatType;
+ typedef typename GetDenseColType::type BaseColType;
+
+ BaseMatType eigvec;
+ BaseColType eigVal;
// Center the data into a temporary matrix.
- arma::mat centeredData = data.each_col() - arma::mean(data, 1);
+ BaseMatType centeredData = arma::conv_to::from(data);
+ centeredData.each_col() -= arma::mean(centeredData, 1);
+
+ // This check cannot happen until here, as `data` may not have a .n_rows
+ // member if it is an expression.
+ if (newDimension > centeredData.n_rows)
+ {
+ std::ostringstream oss;
+ oss << "PCA::Apply(): newDimension (" << newDimension << ") cannot "
+ << "be greater than the existing dimensionality of the data ("
+ << data.n_rows << ")!";
+ throw std::invalid_argument(oss.str());
+ }
// Scale the data if the user ask for.
ScaleData(centeredData);
- decomposition.Apply(data, centeredData, data, eigVal, eigvec, newDimension);
+ decomposition.Apply(data, centeredData, transformedData, eigVal, eigvec,
+ newDimension);
if (newDimension < eigvec.n_rows)
// Drop unnecessary rows.
- data.shed_rows(newDimension, data.n_rows - 1);
+ transformedData.shed_rows(newDimension, data.n_rows - 1);
// The svd method returns only non-zero eigenvalues so we have to calculate
// the right dimension before calculating the amount of variance retained.
@@ -141,21 +202,43 @@ double PCA::Apply(arma::mat& data,
* always be greater than or equal to the varRetained parameter.
*/
template
-double PCA::Apply(arma::mat& data,
+template
+double PCA::Apply(MatType& data,
+ const double varRetained)
+{
+ return Apply(data, data, varRetained);
+}
+
+template
+template
+double PCA::Apply(const MatType& data,
+ OutMatType& transformedData,
const double varRetained)
{
// Parameter validation.
if (varRetained < 0)
- Log::Fatal << "PCA::Apply(): varRetained (" << varRetained << ") must be "
- << "greater than or equal to 0." << std::endl;
- if (varRetained > 1)
- Log::Fatal << "PCA::Apply(): varRetained (" << varRetained << ") should be "
- << "less than or equal to 1." << std::endl;
+ {
+ std::ostringstream oss;
+ oss << "PCA::Apply(): varRetained (" << varRetained << ") must be greater "
+ << "than or equal to 0.";
+ throw std::invalid_argument(oss.str());
+ }
+ else if (varRetained > 1)
+ {
+ std::ostringstream oss;
+ oss << "PCA::Apply(): varRetained (" << varRetained << ") should be less "
+ << "than or equal to 1.";
+ throw std::invalid_argument(oss.str());
+ }
- arma::mat eigvec;
- arma::vec eigVal;
+ typedef typename GetDenseMatType::type BaseMatType;
+ typedef typename GetDenseColType::type BaseColType;
- Apply(data, data, eigVal, eigvec);
+ BaseMatType eigvec;
+ BaseColType eigVal;
+ BaseMatType out;
+
+ Apply(data, transformedData, eigVal, eigvec);
// Calculate the dimension we should keep.
size_t newDimension = 0;
@@ -169,7 +252,7 @@ double PCA::Apply(arma::mat& data,
// varSum is the actual variance we will retain.
if (newDimension < eigVal.n_elem)
- data.shed_rows(newDimension, data.n_rows - 1);
+ transformedData.shed_rows(newDimension, transformedData.n_rows - 1);
return varSum;
}
diff --git a/src/mlpack/methods/quic_svd/quic_svd.hpp b/src/mlpack/methods/quic_svd/quic_svd.hpp
index d3c8d976a5..159fe04f1e 100644
--- a/src/mlpack/methods/quic_svd/quic_svd.hpp
+++ b/src/mlpack/methods/quic_svd/quic_svd.hpp
@@ -51,6 +51,7 @@ namespace mlpack {
* qSVD.Apply(data, u, v, sigma, epsilon, delta);
* @endcode
*/
+template
class QUIC_SVD
{
public:
@@ -64,10 +65,10 @@ class QUIC_SVD
* @param epsilon Error tolerance fraction for calculated subspace.
* @param delta Cumulative probability for Monte Carlo error lower bound.
*/
- QUIC_SVD(const arma::mat& dataset,
- arma::mat& u,
- arma::mat& v,
- arma::mat& sigma,
+ QUIC_SVD(const MatType& dataset,
+ MatType& u,
+ MatType& v,
+ MatType& sigma,
const double epsilon = 0.03,
const double delta = 0.1);
@@ -93,10 +94,10 @@ class QUIC_SVD
* @param epsilon Error tolerance fraction for calculated subspace.
* @param delta Cumulative probability for Monte Carlo error lower bound.
*/
- void Apply(const arma::mat& dataset,
- arma::mat& u,
- arma::mat& v,
- arma::mat& sigma,
+ void Apply(const MatType& dataset,
+ MatType& u,
+ MatType& v,
+ MatType& sigma,
const double epsilon = 0.03,
const double delta = 0.1);
@@ -108,14 +109,14 @@ class QUIC_SVD
* @param v Second unitary matrix.
* @param sigma Diagonal matrix of singular values.
*/
- void ExtractSVD(const arma::mat& dataset,
- arma::mat& u,
- arma::mat& v,
- arma::mat& sigma);
+ void ExtractSVD(const MatType& dataset,
+ MatType& u,
+ MatType& v,
+ MatType& sigma);
private:
//! Subspace basis of the input dataset.
- arma::mat basis;
+ MatType basis;
};
} // namespace mlpack
diff --git a/src/mlpack/methods/quic_svd/quic_svd_impl.hpp b/src/mlpack/methods/quic_svd/quic_svd_impl.hpp
index 4938c62892..f2a0d00124 100644
--- a/src/mlpack/methods/quic_svd/quic_svd_impl.hpp
+++ b/src/mlpack/methods/quic_svd/quic_svd_impl.hpp
@@ -17,39 +17,42 @@
namespace mlpack {
-inline QUIC_SVD::QUIC_SVD(
- const arma::mat& dataset,
- arma::mat& u,
- arma::mat& v,
- arma::mat& sigma,
+template
+inline QUIC_SVD::QUIC_SVD(
+ const MatType& dataset,
+ MatType& u,
+ MatType& v,
+ MatType& sigma,
const double epsilon,
const double delta)
{
Apply(dataset, u, v, sigma, epsilon, delta);
}
-inline QUIC_SVD::QUIC_SVD(
+template
+inline QUIC_SVD::QUIC_SVD(
const double /* epsilon */,
const double /* delta */)
{
/* Nothing to do here */
}
-inline void QUIC_SVD::Apply(
- const arma::mat& dataset,
- arma::mat& u,
- arma::mat& v,
- arma::mat& sigma,
+template
+inline void QUIC_SVD::Apply(
+ const MatType& dataset,
+ MatType& u,
+ MatType& v,
+ MatType& sigma,
const double epsilon,
const double delta)
{
// Since columns are sample in the implementation, the matrix is transposed if
// necessary for maximum speedup.
- CosineTree* ctree;
+ CosineTree* ctree;
if (dataset.n_cols > dataset.n_rows)
- ctree = new CosineTree(dataset, epsilon, delta);
+ ctree = new CosineTree(dataset, epsilon, delta);
else
- ctree = new CosineTree(dataset.t(), epsilon, delta);
+ ctree = new CosineTree(dataset.t(), epsilon, delta);
// Get subspace basis by creating the cosine tree.
ctree->GetFinalBasis(basis);
@@ -62,24 +65,25 @@ inline void QUIC_SVD::Apply(
ExtractSVD(dataset, u, v, sigma);
}
-inline void QUIC_SVD::ExtractSVD(const arma::mat& dataset,
- arma::mat& u,
- arma::mat& v,
- arma::mat& sigma)
+template
+inline void QUIC_SVD::ExtractSVD(const MatType& dataset,
+ MatType& u,
+ MatType& v,
+ MatType& sigma)
{
// Calculate A * V_hat, necessary for further calculations.
- arma::mat projectedMat;
+ MatType projectedMat;
if (dataset.n_cols > dataset.n_rows)
projectedMat = dataset.t() * basis;
else
projectedMat = dataset * basis;
// Calculate the squared projected matrix.
- arma::mat projectedMatSquared = projectedMat.t() * projectedMat;
+ MatType projectedMatSquared = projectedMat.t() * projectedMat;
// Calculate the SVD of the above matrix.
- arma::mat uBar, vBar;
- arma::vec sigmaBar;
+ MatType uBar, vBar;
+ arma::Col sigmaBar;
arma::svd(uBar, sigmaBar, vBar, projectedMatSquared);
// Calculate the approximate SVD of the original matrix, using the SVD of the
@@ -92,7 +96,7 @@ inline void QUIC_SVD::ExtractSVD(const arma::mat& dataset,
// the transposed matrix is not passed.
if (dataset.n_cols > dataset.n_rows)
{
- arma::mat tempMat = u;
+ MatType tempMat = u;
u = v;
v = tempMat;
}
diff --git a/src/mlpack/methods/randomized_svd/randomized_svd.hpp b/src/mlpack/methods/randomized_svd/randomized_svd.hpp
index 652fe3032a..80a6782293 100644
--- a/src/mlpack/methods/randomized_svd/randomized_svd.hpp
+++ b/src/mlpack/methods/randomized_svd/randomized_svd.hpp
@@ -80,10 +80,11 @@ class RandomizedSVD
* @param eps The eps coefficient to avoid division by zero (numerical
* stability).
*/
- RandomizedSVD(const arma::mat& data,
- arma::mat& u,
- arma::vec& s,
- arma::mat& v,
+ template
+ RandomizedSVD(const InMatType& data,
+ MatType& u,
+ VecType& s,
+ MatType& v,
const size_t iteratedPower = 0,
const size_t maxIterations = 2,
const size_t rank = 0,
@@ -113,10 +114,11 @@ class RandomizedSVD
* @param s Diagonal "Sigma" matrix of singular values.
* @param rank Rank of the approximation.
*/
- void Apply(const arma::sp_mat& data,
- arma::mat& u,
- arma::vec& s,
- arma::mat& v,
+ template
+ void Apply(const arma::SpMat& data,
+ MatType& u,
+ VecType& s,
+ MatType& v,
const size_t rank);
/**
@@ -129,10 +131,11 @@ class RandomizedSVD
* @param s Diagonal "Sigma" 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,
+ template