Merge branch 'master' into coot_14

Signed-off-by: Omar Shrit <omar@avontech.fr>
This commit is contained in:
Omar Shrit
2024-04-10 12:25:25 +02:00
45 changed files with 1273 additions and 464 deletions
+3
View File
@@ -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).
+1
View File
@@ -1062,6 +1062,7 @@ div#sidebar {
top: 5px;
min-width: 200px;
font-size: 90%;
max-width: 15.1515%;
}
div#sidebar ul {
+1 -1
View File
@@ -120,7 +120,7 @@ Prepare data for machine learning algorithms.
Transform data from one space to another.
<!-- TODO: add some -->
* [`PCA`](user/methods/pca.md): principal components analysis
### Modeling utilities
+14 -3
View File
@@ -172,9 +172,20 @@ when the sidebar is built for each page.
<!-- Transformations -->
<li>
<a href="LINKROOTindex.html#transformations">
Transformations
</a>
<details>
<summary>
<a href="LINKROOTindex.html#transformations">
Transformations
</a>
</summary>
<ul>
<li>
<a href="LINKROOTuser/methods/pca.html">
<code>PCA</code>
</a>
</li>
</ul>
</details>
</li>
<!-- Modeling utilities -->
+77 -3
View File
@@ -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()`.
<!-- TODO: do something with MakeAlias(); but it needs to be refactored first
-->
* [`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]`),
+2 -2
View File
@@ -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<size_t>&`, 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)
---
+326
View File
@@ -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;
```
<p style="text-align: center; font-size: 85%"><a href="#simple-examples">More examples...</a></p>
#### 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:
<!-- TODO: add link -->
<!-- * [`RADICAL`](radical.md): independent components analysis -->
* [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>
```
`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
<!-- TODO: add link to documentation! -->
* `RandomizedBlockKrylovSVDPolicy`: use the randomized Block Krylov SVD
algorithm to compute the SVD <!-- TODO: add link to documentation! -->
* `QUICSVDPolicy`: use the tree-based `QUIC-SVD` algorithm to compute the SVD
<!-- TODO: add link to documentation -->
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<mlpack::ExactSVDPolicy> pca1;
mlpack::PCA<mlpack::RandomizedSVDPCAPolicy> pca2;
mlpack::PCA<mlpack::RandomizedBlockKrylovSVDPolicy> pca3;
mlpack::PCA<mlpack::QUICSVDPolicy> 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<typename MatType, typename MatType, typename VecType>
static void Apply(const InMatType& data,
const MatType& centeredData,
MatType& transformedData,
VecType& svals,
MatType& svecs,
const size_t rank);
};
```
+1 -1
View File
@@ -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
+49 -70
View File
@@ -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<typename ElemType>
arma::Cube<ElemType> MakeAlias(arma::Cube<ElemType>& input,
const bool strict = true)
template<typename MatType>
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<!IsCube<MatType>::value>* = 0)
{
// Use the advanced constructor.
return arma::Cube<ElemType>(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<typename ElemType>
arma::Mat<ElemType> MakeAlias(arma::Mat<ElemType>& input,
const bool strict = true)
template<typename CubeType>
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<IsCube<CubeType>::value>* = 0)
{
// Use the advanced constructor.
return arma::Mat<ElemType>(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<typename ElemType>
arma::Row<ElemType> MakeAlias(arma::Row<ElemType>& input,
const bool strict = true)
template<typename eT>
void MakeAlias(arma::Mat<eT>& m,
const arma::Mat<eT>& in,
const size_t numRows,
const size_t numCols,
const bool strict = true)
{
// Use the advanced constructor.
return arma::Row<ElemType>(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<typename ElemType>
arma::Col<ElemType> MakeAlias(arma::Col<ElemType>& input,
const bool strict = true)
template<typename eT>
void MakeAlias(arma::SpMat<eT>& m,
const arma::SpMat<eT>& in,
const size_t /* numRows */,
const size_t /* numCols */,
const bool /* strict */)
{
// Use the advanced constructor.
return arma::Col<ElemType>(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<typename ElemType>
arma::SpMat<ElemType> MakeAlias(const arma::SpMat<ElemType>& input,
const bool /* strict */ = true)
{
// Make a copy...
return arma::SpMat<ElemType>(input);
}
/**
* Make a copy of a sparse row (an alias is not possible). The strict
* parameter is ignored.
*/
template<typename ElemType>
arma::SpRow<ElemType> MakeAlias(const arma::SpRow<ElemType>& input,
const bool /* strict */ = true)
{
// Make a copy...
return arma::SpRow<ElemType>(input);
}
/**
* Make a copy of a sparse column (an alias is not possible). The strict
* parameter is ignored.
*/
template<typename ElemType>
arma::SpCol<ElemType> MakeAlias(const arma::SpCol<ElemType>& input,
const bool /* strict */ = true)
{
// Make a copy...
return arma::SpCol<ElemType>(input);
// We can't make aliases of sparse objects, so just copy it.
m = in;
}
/**
@@ -113,16 +93,15 @@ void ClearAlias(arma::Mat<ElemType>& 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<typename ElemType>
void ClearAlias(arma::SpMat<ElemType>& /* mat */)
{
// Nothing to do.
// We cannot make aliases of sparse matrices, so, nothing to do.
}
} // namespace mlpack
#endif
+1
View File
@@ -25,5 +25,6 @@
#include "range.hpp"
#include "shuffle_data.hpp"
#include "trigamma.hpp"
#include "unwrap_alias.hpp"
#endif
+38
View File
@@ -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<typename MatType>
void UnwrapAlias(MatType& m, const MatType& in)
{
MakeAlias(m, in, in.n_rows, in.n_cols);
}
template<typename MatType, typename InMatType>
void UnwrapAlias(MatType& m,
const InMatType& in)
{
m = in;
}
} // namespace mlpack
#endif
@@ -20,14 +20,20 @@ namespace mlpack {
// Predeclare classes for CosineNodeQueue typedef.
class CompareCosineNode;
template<typename MatType>
class CosineTree;
// CosineNodeQueue typedef.
typedef std::vector<CosineTree*> CosineNodeQueue;
template<typename MatType = arma::mat>
using CosineNodeQueue = std::vector<CosineTree<MatType>*>;
template<typename MatType = arma::mat>
class CosineTree
{
public:
typedef typename GetDenseColType<MatType>::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<MatType>& 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<MatType>& 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<MatType>& 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<size_t>& 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<size_t>& 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<size_t> 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<typename MatType>
bool operator() (const CosineTree<MatType>* a,
const CosineTree<MatType>* b) const
{
return a->L2Error() < b->L2Error();
}
@@ -16,7 +16,8 @@
namespace mlpack {
inline CosineTree::CosineTree(const arma::mat& dataset) :
template<typename MatType>
inline CosineTree<MatType>::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<size_t>& subIndices) :
template<typename MatType>
inline CosineTree<MatType>::CosineTree(CosineTree& parentNode,
const std::vector<size_t>& 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<typename MatType>
inline CosineTree<MatType>::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<MatType> 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<VecType>(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<typename MatType>
inline CosineTree<MatType>::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<typename MatType>
inline CosineTree<MatType>& CosineTree<MatType>::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<typename MatType>
inline CosineTree<MatType>::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<typename MatType>
inline CosineTree<MatType>& CosineTree<MatType>::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<typename MatType>
inline CosineTree<MatType>::~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<typename MatType>
inline void CosineTree<MatType>::ModifiedGramSchmidt(
CosineNodeQueue<MatType>& treeQueue,
typename CosineTree<MatType>::VecType& centroid,
typename CosineTree<MatType>::VecType& newBasisVector,
typename CosineTree<MatType>::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<MatType>::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<typename MatType>
inline double CosineTree<MatType>::MonteCarloError(
CosineTree* node,
CosineNodeQueue<MatType>& treeQueue,
typename CosineTree<MatType>::VecType* addBasisVector1,
typename CosineTree<MatType>::VecType* addBasisVector2)
{
std::vector<size_t> 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<MatType>::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<typename MatType>
inline void CosineTree<MatType>::ConstructBasis(
CosineNodeQueue<MatType>& 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<MatType>::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<typename MatType>
inline void CosineTree<MatType>::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<size_t>& sampledIndices,
arma::vec& probabilities,
size_t numSamples)
template<typename MatType>
inline void CosineTree<MatType>::ColumnSamplesLS(
std::vector<size_t>& sampledIndices,
typename CosineTree<MatType>::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<size_t>& sampledIndices,
}
}
inline size_t CosineTree::ColumnSampleLS()
template<typename MatType>
inline size_t CosineTree<MatType>::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<typename MatType>
inline size_t CosineTree<MatType>::BinarySearch(
typename CosineTree<MatType>::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<typename MatType>
inline void CosineTree<MatType>::CalculateCosines(
typename CosineTree<MatType>::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<typename MatType>
inline void CosineTree<MatType>::CalculateCentroid()
{
// Initialize centroid as vector of zeros.
centroid.zeros(dataset->n_rows);
+45
View File
@@ -254,4 +254,49 @@ struct GetSparseMatType<arma::SpMat<eT>>
typedef arma::SpMat<eT> type;
};
// Get whether or not the given type is a base matrix type (e.g. not an
// expression).
template<typename MatType>
struct IsBaseMatType
{
constexpr static bool value = false;
};
template<typename eT>
struct IsBaseMatType<arma::Mat<eT>>
{
constexpr static bool value = true;
};
template<typename eT>
struct IsBaseMatType<arma::Col<eT>>
{
constexpr static bool value = true;
};
template<typename eT>
struct IsBaseMatType<arma::Row<eT>>
{
constexpr static bool value = true;
};
template<typename eT>
struct IsBaseMatType<arma::SpMat<eT>>
{
constexpr static bool value = true;
};
template<typename eT>
struct IsBaseMatType<arma::SpCol<eT>>
{
constexpr static bool value = true;
};
template<typename eT>
struct IsBaseMatType<arma::SpRow<eT>>
{
constexpr static bool value = true;
};
#endif
-1
View File
@@ -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"
-2
View File
@@ -15,8 +15,6 @@
// In case it hasn't been included yet.
#include "ffn.hpp"
#include "make_alias.hpp"
namespace mlpack {
template<typename OutputLayerType,
@@ -12,7 +12,6 @@
#ifndef MLPACK_METHODS_ANN_LAYER_ADD_MERGE_HPP
#define MLPACK_METHODS_ANN_LAYER_ADD_MERGE_HPP
#include "../make_alias.hpp"
#include "multi_layer.hpp"
namespace mlpack {
+1 -1
View File
@@ -220,7 +220,7 @@ class Layer
//! Get the layer loss. Overload this if the layer should add any extra loss
//! to the loss function when computing the objective. (TODO: better comment)
virtual double Loss() { return 0; }
virtual double Loss() const { return 0; }
//! Get the input dimensions.
const std::vector<size_t>& InputDimensions() const { return inputDimensions; }
@@ -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 {
-58
View File
@@ -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 <mlpack/prereqs.hpp>
namespace mlpack {
/**
* Reconstruct `m` as an alias around the memory `newMem`, with size `numRows` x
* `numCols`.
*/
template<typename MatType>
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<typename CubeType>
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
@@ -19,13 +19,14 @@
namespace mlpack {
template <typename MatType>
BiasSVDFunction<MatType>::BiasSVDFunction(const MatType& data,
BiasSVDFunction<MatType>::BiasSVDFunction(const MatType& dataIn,
const size_t rank,
const double lambda) :
data(MakeAlias(const_cast<MatType&>(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;
@@ -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<typename InMatType, typename MatType, typename VecType>
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<typename InMatType, typename MatType, typename VecType>
void Apply(const InMatType& data,
MatType& u,
VecType& s,
MatType& v,
const size_t rank);
//! Get the number of iterations for the power method.
@@ -16,11 +16,12 @@
namespace mlpack {
template<typename InMatType, typename MatType, typename VecType>
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<typename InMatType, typename MatType, typename VecType>
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<MatType>(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);
@@ -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.
@@ -24,18 +24,19 @@ namespace mlpack {
template<typename MatType, typename ParametersType>
LinearSVMFunction<MatType, ParametersType>::LinearSVMFunction(
const MatType& dataset,
const MatType& datasetIn,
const arma::Row<size_t>& labels,
const size_t numClasses,
const double lambda,
const double delta,
const bool fitIntercept) :
dataset(MakeAlias(const_cast<MatType&>(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;
@@ -19,23 +19,24 @@
namespace mlpack {
template<typename MetricType>
LMNNFunction<MetricType>::LMNNFunction(const arma::mat& dataset,
const arma::Row<size_t>& labels,
LMNNFunction<MetricType>::LMNNFunction(const arma::mat& datasetIn,
const arma::Row<size_t>& labelsIn,
size_t k,
double regularization,
size_t range,
MetricType metric) :
dataset(MakeAlias(const_cast<arma::mat&>(dataset), false)),
labels(MakeAlias(const_cast<arma::Row<size_t>&>(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.
@@ -21,15 +21,17 @@ namespace mlpack {
template<typename MatType>
LogisticRegressionFunction<MatType>::LogisticRegressionFunction(
const MatType& predictors,
const arma::Row<size_t>& responses,
const MatType& predictorsIn,
const arma::Row<size_t>& responsesIn,
const double lambda) :
// We promise to be well-behaved... the elements won't be modified.
predictors(MakeAlias(const_cast<MatType&>(predictors), false)),
responses(MakeAlias(const_cast<arma::Row<size_t>&>(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)
{
@@ -22,14 +22,15 @@ namespace mlpack {
// Initialize with the given kernel.
template<typename MetricType>
SoftmaxErrorFunction<MetricType>::SoftmaxErrorFunction(
const arma::mat& dataset,
const arma::Row<size_t>& labels,
const arma::mat& datasetIn,
const arma::Row<size_t>& labelsIn,
MetricType metric) :
dataset(MakeAlias(const_cast<arma::mat&>(dataset), false)),
labels(MakeAlias(const_cast<arma::Row<size_t>&>(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<typename MetricType>
@@ -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<typename InMatType, typename MatType, typename VecType>
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;
@@ -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<typename InMatType, typename MatType, typename VecType>
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<MatType> 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;
@@ -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<typename InMatType, typename MatType, typename VecType>
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;
@@ -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<typename InMatType, typename MatType, typename VecType>
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;
+72 -14
View File
@@ -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<typename MatType = arma::mat,
typename OutMatType = arma::mat,
typename VecType = arma::vec>
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<typename MatType = arma::mat,
typename OutMatType = arma::mat,
typename VecType = arma::vec>
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<typename MatType = arma::mat, typename OutMatType = arma::mat>
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<typename MatType = arma::mat>
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<typename MatType = arma::mat>
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<typename MatType = arma::mat, typename OutMatType = arma::mat>
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<typename MatType = arma::mat, typename OutMatType = arma::mat>
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<typename MatType>
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<typename MatType = arma::mat, typename OutMatType = arma::mat>
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<typename MatType>
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<typename MatType::elem_type> stdDev = arma::stddev(
centeredData, 0, 1 /* for each dimension */);
// If there are any zeroes, make them very small.
+120 -37
View File
@@ -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<DecompositionPolicy>::PCA(
* @param eigvec - PCA Loadings/Coeffs/EigenVectors
*/
template<typename DecompositionPolicy>
void PCA<DecompositionPolicy>::Apply(const arma::mat& data,
arma::mat& transformedData,
arma::vec& eigVal,
arma::mat& eigvec)
template<typename MatType, typename OutMatType, typename VecType>
void PCA<DecompositionPolicy>::Apply(const MatType& data,
OutMatType& transformedData,
VecType& eigVal,
OutMatType& eigvec)
{
// Sanity checks on input types.
static_assert(IsBaseMatType<OutMatType>::value,
"PCA::Apply(): transformedData must be a matrix type!");
static_assert(IsBaseMatType<VecType>::value,
"PCA::Apply(): eigVal must be a vector type!");
static_assert(std::is_same<typename MatType::elem_type,
typename OutMatType::elem_type>::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<OutMatType>::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<DecompositionPolicy>::Apply(const arma::mat& data,
* @param eigVal - contains eigen values in a column vector
*/
template<typename DecompositionPolicy>
void PCA<DecompositionPolicy>::Apply(const arma::mat& data,
arma::mat& transformedData,
arma::vec& eigVal)
template<typename MatType, typename OutMatType, typename VecType>
void PCA<DecompositionPolicy>::Apply(const MatType& data,
OutMatType& transformedData,
VecType& eigVal)
{
arma::mat eigvec;
// Sanity checks on input types.
static_assert(IsBaseMatType<OutMatType>::value,
"PCA::Apply(): transformedData must be a matrix type!");
static_assert(IsBaseMatType<VecType>::value,
"PCA::Apply(): eigVal must be a vector type!");
static_assert(std::is_same<typename MatType::elem_type,
typename OutMatType::elem_type>::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<DecompositionPolicy>::Apply(const arma::mat& data,
* @param transformedData Data with PCA applied.
*/
template<typename DecompositionPolicy>
void PCA<DecompositionPolicy>::Apply(const arma::mat& data,
arma::mat& transformedData)
template<typename MatType, typename OutMatType>
void PCA<DecompositionPolicy>::Apply(const MatType& data,
OutMatType& transformedData)
{
arma::mat eigvec;
arma::vec eigVal;
// Sanity checks on input types.
static_assert(IsBaseMatType<OutMatType>::value,
"PCA::Apply(): transformedData must be a matrix type!");
static_assert(std::is_same<typename MatType::elem_type,
typename OutMatType::elem_type>::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<MatType>::type BaseColType;
OutMatType eigvec;
BaseColType eigVal;
Apply(data, transformedData, eigVal, eigvec);
}
@@ -95,32 +130,58 @@ void PCA<DecompositionPolicy>::Apply(const arma::mat& data,
* @return Amount of the variance of the data retained (between 0 and 1).
*/
template<typename DecompositionPolicy>
double PCA<DecompositionPolicy>::Apply(arma::mat& data,
template<typename MatType>
double PCA<DecompositionPolicy>::Apply(MatType& data,
const size_t newDimension)
{
return Apply(data, data, newDimension);
}
template<typename DecompositionPolicy>
template<typename MatType, typename OutMatType>
double PCA<DecompositionPolicy>::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<MatType>::type BaseMatType;
typedef typename GetDenseColType<MatType>::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<OutMatType>::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<DecompositionPolicy>::Apply(arma::mat& data,
* always be greater than or equal to the varRetained parameter.
*/
template<typename DecompositionPolicy>
double PCA<DecompositionPolicy>::Apply(arma::mat& data,
template<typename MatType>
double PCA<DecompositionPolicy>::Apply(MatType& data,
const double varRetained)
{
return Apply(data, data, varRetained);
}
template<typename DecompositionPolicy>
template<typename MatType, typename OutMatType>
double PCA<DecompositionPolicy>::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<MatType>::type BaseMatType;
typedef typename GetDenseColType<MatType>::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<DecompositionPolicy>::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;
}
+14 -13
View File
@@ -51,6 +51,7 @@ namespace mlpack {
* qSVD.Apply(data, u, v, sigma, epsilon, delta);
* @endcode
*/
template<typename MatType = arma::mat>
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
+27 -23
View File
@@ -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<typename MatType>
inline QUIC_SVD<MatType>::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<typename MatType>
inline QUIC_SVD<MatType>::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<typename MatType>
inline void QUIC_SVD<MatType>::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<MatType>* ctree;
if (dataset.n_cols > dataset.n_rows)
ctree = new CosineTree(dataset, epsilon, delta);
ctree = new CosineTree<MatType>(dataset, epsilon, delta);
else
ctree = new CosineTree(dataset.t(), epsilon, delta);
ctree = new CosineTree<MatType>(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<typename MatType>
inline void QUIC_SVD<MatType>::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<typename MatType::elem_type> 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;
}
@@ -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<typename InMatType, typename MatType, typename VecType>
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<typename eT, typename MatType, typename VecType>
void Apply(const arma::SpMat<eT>& 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<typename InMatType, typename MatType, typename VecType>
void Apply(const InMatType& data,
MatType& u,
VecType& s,
MatType& v,
const size_t rank);
/**
@@ -146,13 +149,16 @@ class RandomizedSVD
* @param rank Rank of the approximation.
* @param rowMean Centered mean value matrix.
*/
template<typename MatType>
void Apply(const MatType& data,
arma::mat& u,
arma::vec& s,
arma::mat& v,
template<typename InMatType,
typename MatType,
typename VecType,
typename MeanType>
void Apply(const InMatType& data,
MatType& u,
VecType& s,
MatType& v,
const size_t rank,
MatType rowMean);
const MeanType& rowMean);
//! Get the size of the normalized power iterations.
size_t IteratedPower() const { return iteratedPower; }
@@ -17,11 +17,12 @@
namespace mlpack {
template<typename InMatType, typename MatType, typename VecType>
inline RandomizedSVD::RandomizedSVD(
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 iteratedPower,
const size_t maxIterations,
const size_t rank,
@@ -51,42 +52,50 @@ inline RandomizedSVD::RandomizedSVD(
/* Nothing to do here */
}
inline void RandomizedSVD::Apply(const arma::sp_mat& data,
arma::mat& u,
arma::vec& s,
arma::mat& v,
template<typename eT, typename MatType, typename VecType>
inline void RandomizedSVD::Apply(const arma::SpMat<eT>& data,
MatType& u,
VecType& s,
MatType& v,
const size_t rank)
{
// Center the data into a temporary matrix for sparse matrix.
arma::sp_mat rowMean = sum(data, 1) / data.n_cols;
arma::SpMat<eT> rowMean = sum(data, 1) / data.n_cols;
Apply(data, u, s, v, rank, rowMean);
}
inline void RandomizedSVD::Apply(const arma::mat& data,
arma::mat& u,
arma::vec& s,
arma::mat& v,
template<typename InMatType, typename MatType, typename VecType>
inline void RandomizedSVD::Apply(const InMatType& dataIn,
MatType& u,
VecType& s,
MatType& v,
const size_t rank)
{
// Center the data into a temporary matrix.
arma::mat rowMean = sum(data, 1) / data.n_cols + eps;
MatType data;
UnwrapAlias(data, dataIn);
MatType rowMean = sum(data, 1) / data.n_cols + eps;
Apply(data, u, s, v, rank, rowMean);
}
template<typename MatType>
inline void RandomizedSVD::Apply(const MatType& data,
arma::mat& u,
arma::vec& s,
arma::mat& v,
template<typename InMatType,
typename MatType,
typename VecType,
typename MeanType>
inline void RandomizedSVD::Apply(const InMatType& data,
MatType& u,
VecType& s,
MatType& v,
const size_t rank,
MatType rowMean)
const MeanType& rowMean)
{
if (iteratedPower == 0)
iteratedPower = rank + 2;
arma::mat R, Q, Qdata;
MatType R, Q, Qdata;
// Apply the centered data matrix to a random matrix, obtaining Q.
if (data.n_cols >= data.n_rows)
@@ -97,14 +106,14 @@ inline void RandomizedSVD::Apply(const MatType& data,
else
{
R.randn(data.n_cols, iteratedPower);
Q = (data * R) - (rowMean * (ones(1, data.n_cols) * R));
Q = (data * R) - (rowMean * (ones<MatType>(1, data.n_cols) * R));
}
// Form a matrix Q whose columns constitute a
// well-conditioned basis for the columns of the earlier Q.
if (maxIterations == 0)
{
arma::qr_econ(Q, v, Q);
arma::qr_econ(Q, v, Q);
}
else
{
@@ -116,7 +125,7 @@ inline void RandomizedSVD::Apply(const MatType& data,
{
if (data.n_cols >= data.n_rows)
{
Q = (data * Q) - rowMean * (ones(1, data.n_cols) * Q);
Q = (data * Q) - rowMean * (ones<MatType>(1, data.n_cols) * Q);
arma::lu(Q, v, Q);
Q = (data.t() * Q) - repmat(rowMean.t() * Q, data.n_cols, 1);
}
@@ -124,7 +133,7 @@ inline void RandomizedSVD::Apply(const MatType& data,
{
Q = (data.t() * Q) - repmat(rowMean.t() * Q, data.n_cols, 1);
arma::lu(Q, v, Q);
Q = (data * Q) - (rowMean * (ones(1, data.n_cols) * Q));
Q = (data * Q) - (rowMean * (ones<MatType>(1, data.n_cols) * Q));
}
// Computing the LU decomposition is more efficient than computing the QR
@@ -146,7 +155,7 @@ inline void RandomizedSVD::Apply(const MatType& data,
// applied to Q.
if (data.n_cols >= data.n_rows)
{
Qdata = (data * Q) - rowMean * (ones(1, data.n_cols) * Q);
Qdata = (data * Q) - rowMean * (ones<MatType>(1, data.n_cols) * Q);
arma::svd_econ(u, s, v, Qdata);
v = Q * v;
}
@@ -18,13 +18,14 @@
namespace mlpack {
template <typename MatType>
RegularizedSVDFunction<MatType>::RegularizedSVDFunction(const MatType& data,
RegularizedSVDFunction<MatType>::RegularizedSVDFunction(const MatType& dataIn,
const size_t rank,
const double lambda) :
data(MakeAlias(const_cast<MatType&>(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;
@@ -18,16 +18,17 @@ namespace mlpack {
template<typename MatType>
inline SoftmaxRegressionFunction<MatType>::SoftmaxRegressionFunction(
const MatType& data,
const MatType& dataIn,
const arma::Row<size_t>& labels,
const size_t numClasses,
const double lambda,
const bool fitIntercept) :
data(MakeAlias(const_cast<MatType&>(data), false)),
numClasses(numClasses),
lambda(lambda),
fitIntercept(fitIntercept)
{
MakeAlias(data, dataIn, dataIn.n_rows, dataIn.n_cols, false);
// Initialize the parameters to suitable values.
initialPoint = InitializeWeights();
@@ -20,15 +20,16 @@ namespace mlpack {
template <typename MatType>
SVDPlusPlusFunction<MatType>::SVDPlusPlusFunction(
const MatType& data,
const MatType& dataIn,
const arma::sp_mat& implicitData,
const size_t rank,
const double lambda) :
data(MakeAlias(const_cast<MatType&>(data), false)),
implicitData(implicitData),
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;
+25 -25
View File
@@ -34,7 +34,7 @@ TEST_CASE("CosineTreeNoSplit", "[CosineTreeTest]")
// Make a cosine tree, with the generated dataset and the defined constants.
// Note that the value of epsilon is one.
CosineTree ctree(data, epsilon, delta);
CosineTree<> ctree(data, epsilon, delta);
arma::mat basis;
ctree.GetFinalBasis(basis);
@@ -57,17 +57,17 @@ TEST_CASE("CosineNodeCosineSplit", "[CosineTreeTest]")
// Make a random dataset and the root object.
arma::mat data = arma::randu(numRows, numCols);
CosineTree root(data);
CosineTree<> root(data);
// Stack for depth first search of the tree.
std::vector<CosineTree*> nodeStack;
std::vector<CosineTree<>*> nodeStack;
nodeStack.push_back(&root);
// While stack is not empty.
while (nodeStack.size())
{
// Pop a node from the stack and split it.
CosineTree *currentNode, *currentLeft, *currentRight;
CosineTree<> *currentNode, *currentLeft, *currentRight;
currentNode = nodeStack.back();
currentNode->CosineNodeSplit();
nodeStack.pop_back();
@@ -173,14 +173,14 @@ TEST_CASE("CosineTreeModifiedGramSchmidt", "[CosineTreeTest]")
// Declare a queue and a dummy CosineTree object.
CompareCosineNode comp;
CosineNodeQueue basisQueue;
CosineTree dummyTree(data, epsilon, delta);
CosineNodeQueue<> basisQueue;
CosineTree<> dummyTree(data, epsilon, delta);
for (size_t i = 0; i < numCols; ++i)
{
// Make a new CosineNode object.
CosineTree* basisNode;
basisNode = new CosineTree(data);
CosineTree<>* basisNode;
basisNode = new CosineTree<>(data);
// Use the columns of the dataset as random centroids.
arma::vec centroid = data.col(i);
@@ -190,8 +190,8 @@ TEST_CASE("CosineTreeModifiedGramSchmidt", "[CosineTreeTest]")
dummyTree.ModifiedGramSchmidt(basisQueue, centroid, newBasisVector);
// Check if the obtained vector is orthonormal to the basis vectors.
CosineNodeQueue::const_iterator j = basisQueue.cbegin();
CosineTree* currentNode;
CosineNodeQueue<>::const_iterator j = basisQueue.cbegin();
CosineTree<>* currentNode;
for (; j != basisQueue.cend(); ++j)
{
@@ -210,7 +210,7 @@ TEST_CASE("CosineTreeModifiedGramSchmidt", "[CosineTreeTest]")
// Deallocate memory given to the objects.
for (size_t i = 0; i < numCols; ++i)
{
CosineTree* currentNode;
CosineTree<>* currentNode;
currentNode = basisQueue.front();
std::pop_heap(basisQueue.begin(), basisQueue.end(), comp);
basisQueue.pop_back();
@@ -236,17 +236,17 @@ TEST_CASE("CopyConstructorAndOperatorCosineTreeTest", "[CosineTreeTest]")
arma::mat* data = new arma::mat(numRows, numCols, arma::fill::randu);
// Make a cosine tree, with the generated dataset.
CosineTree* ctree1 = new CosineTree(*data);
CosineTree<>* ctree1 = new CosineTree<>(*data);
// Stacks for depth first search of the tree.
std::vector<CosineTree*> nodeStack1, nodeStack2, nodeStack3;
std::vector<CosineTree<>*> nodeStack1, nodeStack2, nodeStack3;
nodeStack1.push_back(ctree1);
// While stack is not empty.
while (nodeStack1.size())
{
// Pop a node from the stack and split it.
CosineTree *currentNode1, *currentLeft1, *currentRight1;
CosineTree<> *currentNode1, *currentLeft1, *currentRight1;
currentNode1 = nodeStack1.back();
currentNode1->CosineNodeSplit();
@@ -268,8 +268,8 @@ TEST_CASE("CopyConstructorAndOperatorCosineTreeTest", "[CosineTreeTest]")
}
// Copy constructor and operator.
CosineTree ctree2(*ctree1);
CosineTree ctree3 = *ctree1;
CosineTree<> ctree2(*ctree1);
CosineTree<> ctree3 = *ctree1;
delete ctree1;
delete data;
@@ -281,8 +281,8 @@ TEST_CASE("CopyConstructorAndOperatorCosineTreeTest", "[CosineTreeTest]")
while (nodeStack2.size() && nodeStack3.size())
{
// Pop a node from the stack and split it.
CosineTree *currentNode2, *currentLeft2, *currentRight2;
CosineTree *currentNode3, *currentLeft3, *currentRight3;
CosineTree<> *currentNode2, *currentLeft2, *currentRight2;
CosineTree<> *currentNode3, *currentLeft3, *currentRight3;
currentNode2 = nodeStack2.back();
nodeStack2.pop_back();
@@ -337,17 +337,17 @@ TEST_CASE("MoveConstructorAndOperatorCosineTreeTest", "[CosineTreeTest]")
arma::mat data = arma::randu(numRows, numCols);
// Make a cosine tree, with the generated dataset.
CosineTree ctree1(data);
CosineTree<> ctree1(data);
// Stacks for depth first search of the tree.
std::vector<CosineTree*> nodeStack1, nodeStack2, nodeStack3;
std::vector<CosineTree<>*> nodeStack1, nodeStack2, nodeStack3;
nodeStack1.push_back(&ctree1);
// While stack is not empty.
while (nodeStack1.size())
{
// Pop a node from the stack and split it.
CosineTree *currentNode1, *currentLeft1, *currentRight1;
CosineTree<> *currentNode1, *currentLeft1, *currentRight1;
currentNode1 = nodeStack1.back();
currentNode1->CosineNodeSplit();
@@ -369,7 +369,7 @@ TEST_CASE("MoveConstructorAndOperatorCosineTreeTest", "[CosineTreeTest]")
}
// Move constructor.
CosineTree ctree2(std::move(ctree1));
CosineTree<> ctree2(std::move(ctree1));
nodeStack2.push_back(&ctree2);
@@ -377,7 +377,7 @@ TEST_CASE("MoveConstructorAndOperatorCosineTreeTest", "[CosineTreeTest]")
while (nodeStack2.size())
{
// Pop a node from the stack and split it.
CosineTree *currentNode2, *currentLeft2, *currentRight2;
CosineTree<> *currentNode2, *currentLeft2, *currentRight2;
currentNode2 = nodeStack2.back();
nodeStack2.pop_back();
@@ -398,7 +398,7 @@ TEST_CASE("MoveConstructorAndOperatorCosineTreeTest", "[CosineTreeTest]")
}
// Move operator.
CosineTree ctree3 = std::move(ctree2);
CosineTree<> ctree3 = std::move(ctree2);
nodeStack3.push_back(&ctree3);
@@ -406,7 +406,7 @@ TEST_CASE("MoveConstructorAndOperatorCosineTreeTest", "[CosineTreeTest]")
while (nodeStack3.size())
{
// Pop a node from the stack and split it.
CosineTree *currentNode3, *currentLeft3, *currentRight3;
CosineTree<> *currentNode3, *currentLeft3, *currentRight3;
currentNode3 = nodeStack3.back();
nodeStack3.pop_back();
+177
View File
@@ -330,3 +330,180 @@ TEST_CASE("PCAScalingTest", "[PCATest]")
// The eigenvalues should sum to three.
REQUIRE(accu(eigval) == Approx(3.0).epsilon(0.001));
}
/**
* Test PCA on a subview of a matrix with different decomposition strategies.
*/
TEMPLATE_TEST_CASE("PCASubviewTest", "[PCATest]", ExactSVDPolicy,
RandomizedSVDPCAPolicy, RandomizedBlockKrylovSVDPolicy, QUICSVDPolicy)
{
typedef TestType DecompositionPolicy;
// Generate an artifical dataset in 10 dimensions.
arma::mat data(3, 5000);
arma::vec mean("1.0 3.0 -12.0");
arma::mat cov("1.0 0.9 0.0;"
"0.9 1.0 0.0;"
"0.0 0.0 12.0");
GaussianDistribution g(mean, cov);
for (size_t i = 0; i < 5000; ++i)
data.col(i) = g.Random();
// Compute PCA on the first 2000 points.
arma::mat transData1, transData2, transData3, eigvec;
arma::vec eigval1, eigval2;
PCA<DecompositionPolicy> p;
p.Apply(data.cols(0, 1999), transData1);
p.Apply(data.cols(0, 1999), transData2, eigval1);
p.Apply(data.cols(0, 1999), transData3, eigval2, eigvec);
// Only check for deterministic policies.
if (std::is_same<DecompositionPolicy, ExactSVDPolicy>::value)
{
arma::mat trueTransData, trueEigvec;
arma::vec trueEigval;
arma::mat dataSub = data.cols(0, 1999);
p.Apply(dataSub, trueTransData, trueEigval, trueEigvec);
REQUIRE(arma::approx_equal(transData1, trueTransData, "both", 1e-5, 1e-5));
REQUIRE(arma::approx_equal(transData2, trueTransData, "both", 1e-5, 1e-5));
REQUIRE(arma::approx_equal(transData3, trueTransData, "both", 1e-5, 1e-5));
REQUIRE(arma::approx_equal(eigval1, trueEigval, "both", 1e-5, 1e-5));
REQUIRE(arma::approx_equal(eigval2, trueEigval, "both", 1e-5, 1e-5));
REQUIRE(arma::approx_equal(eigvec, trueEigvec, "both", 1e-5, 1e-5));
}
}
/**
* Test PCA on an input expression.
*/
TEMPLATE_TEST_CASE("PCAExpressionTest", "[PCATest]", ExactSVDPolicy,
RandomizedSVDPCAPolicy, RandomizedBlockKrylovSVDPolicy, QUICSVDPolicy)
{
typedef TestType DecompositionPolicy;
// Generate an artifical dataset in 10 dimensions.
arma::mat data(3, 5000);
arma::vec mean("1.0 3.0 -12.0");
arma::mat cov("1.0 0.9 0.0;"
"0.9 1.0 0.0;"
"0.0 0.0 12.0");
GaussianDistribution g(mean, cov);
for (size_t i = 0; i < 5000; ++i)
data.col(i) = g.Random();
// Compute PCA on an expression involving the input matrix.
arma::mat transData1, transData2, transData3, eigvec;
arma::vec eigval1, eigval2;
PCA<DecompositionPolicy> p;
p.Apply(2 * data + 1, transData1);
p.Apply(2 * data + 1, transData2, eigval1);
p.Apply(2 * data + 1, transData3, eigval2, eigvec);
// Only check for deterministic policies.
if (std::is_same<DecompositionPolicy, ExactSVDPolicy>::value)
{
arma::mat trueTransData, trueEigvec;
arma::vec trueEigval;
arma::mat dataSub = 2 * data + 1;
p.Apply(dataSub, trueTransData, trueEigval, trueEigvec);
REQUIRE(arma::approx_equal(transData1, trueTransData, "both", 1e-5, 1e-5));
REQUIRE(arma::approx_equal(transData2, trueTransData, "both", 1e-5, 1e-5));
REQUIRE(arma::approx_equal(transData3, trueTransData, "both", 1e-5, 1e-5));
REQUIRE(arma::approx_equal(eigval1, trueEigval, "both", 1e-5, 1e-5));
REQUIRE(arma::approx_equal(eigval2, trueEigval, "both", 1e-5, 1e-5));
REQUIRE(arma::approx_equal(eigvec, trueEigvec, "both", 1e-5, 1e-5));
}
}
/**
* Test PCA on 32-bit data.
*/
TEMPLATE_TEST_CASE("PCAFloatTest", "[PCATest]", ExactSVDPolicy,
RandomizedSVDPCAPolicy, RandomizedBlockKrylovSVDPolicy, QUICSVDPolicy)
{
typedef TestType DecompositionPolicy;
// Generate an artifical dataset in 10 dimensions.
arma::fmat data(3, 5000);
arma::vec mean("1.0 3.0 -12.0");
arma::mat cov("1.0 0.9 0.0;"
"0.9 1.0 0.0;"
"0.0 0.0 12.0");
GaussianDistribution g(mean, cov);
for (size_t i = 0; i < 5000; ++i)
data.col(i) = arma::conv_to<arma::fvec>::from(g.Random());
// Compute PCA on the floating-point data.
arma::fmat coeff, coeff1, score, score1;
arma::fvec eigVal, eigVal1;
PCA<DecompositionPolicy> p;
p.Apply(data, score1, eigVal1, coeff1);
princomp(coeff, score, eigVal, trans(data));
// Verify the PCA results based on the eigenvalues. We don't check for
// QUIC-SVD, since that method has a lot of noise.
if (!std::is_same<DecompositionPolicy, QUICSVDPolicy>::value)
{
for (size_t i = 0; i < eigVal.n_elem; ++i)
{
if (eigVal[i] == 0.0)
REQUIRE(eigVal1[i] == Approx(0.0).margin(1e-5));
else
REQUIRE(eigVal[i] == Approx(eigVal1[i]).epsilon(1e-3));
}
}
}
/**
* Test that we can convert sparse input matrices to dense output matrices.
* We check the general shape of the output, but not the exact details---those
* are handled in other tests.
*/
TEMPLATE_TEST_CASE("PCASparseToDenseTest", "[PCATest]", float, double)
{
typedef arma::Mat<TestType> MatType;
typedef arma::SpMat<TestType> SpMatType;
SpMatType dataset;
dataset.sprandu(1000, 50000, 0.01);
MatType transformedDataset1, transformedDataset2;
PCA<> p;
const double varRetained1 = p.Apply(dataset, transformedDataset1, 5);
const double varRetained2 = p.Apply(dataset, transformedDataset2, 0.6);
REQUIRE(transformedDataset1.n_cols == dataset.n_cols);
REQUIRE(transformedDataset2.n_cols == dataset.n_cols);
REQUIRE(varRetained1 >= 0.0);
REQUIRE(varRetained1 <= 1.0);
REQUIRE(varRetained2 >= 0.0);
REQUIRE(varRetained2 <= 1.0);
REQUIRE(transformedDataset1.n_rows == 5);
REQUIRE(transformedDataset2.n_rows <= dataset.n_rows);
// Ensure we get basically the same as if we had done it to a dense matrix.
MatType denseData1(dataset);
MatType denseData2(dataset);
const double varRetained3 = p.Apply(denseData1, 5);
const double varRetained4 = p.Apply(denseData2, 0.6);
REQUIRE(varRetained1 == Approx(varRetained3));
REQUIRE(varRetained2 == Approx(varRetained4));
REQUIRE(denseData2.n_rows == transformedDataset2.n_rows);
}
+2 -2
View File
@@ -192,7 +192,7 @@ TEST_CASE("OUNoiseTest", "[PolicyGradientTest]")
// Test the reset function.
ouNoise.reset();
arma::colvec state = ouNoise.sample();
REQUIRE(state.n_elem == size);
REQUIRE((int) state.n_elem == size);
// Verify that the sample is not equal to the reset state.
arma::colvec sample = ouNoise.sample();
@@ -213,7 +213,7 @@ TEST_CASE("GaussianNoiseTest", "[PolicyGradientTest]")
// Test the sample function.
arma::colvec noise = gaussianNoise.sample();
REQUIRE(noise.n_elem == size);
REQUIRE((int) noise.n_elem == size);
// Verify that the noise vector has values drawn from a
// Gaussian distribution with the specified mean and standard deviation.
+3 -3
View File
@@ -34,7 +34,7 @@ TEST_CASE("QUICSVDReconstructionError", "[QUICSVDTest]")
{
// Obtain the SVD using default parameters.
arma::mat u, v, sigma;
QUIC_SVD quicsvd(dataset, u, v, sigma);
QUIC_SVD<> quicsvd(dataset, u, v, sigma);
// Reconstruct the matrix using the SVD.
arma::mat reconstruct;
@@ -71,7 +71,7 @@ TEST_CASE("QUICSVDSingularValueError", "[QUICSVDTest]")
// Obtain the SVD using default parameters.
arma::svd_econ(U1, s1, V1, data);
QUIC_SVD quicsvd(data, U1, V1, s2);
QUIC_SVD<> quicsvd(data, U1, V1, s2);
s3 = arma::diagvec(s2);
s1 = s1.subvec(0, s3.n_elem - 1);
@@ -87,5 +87,5 @@ TEST_CASE("QUICSVDSameDimensionTest", "[QUICSVDTest]")
// Obtain the SVD using default parameters.
arma::mat u, v, sigma;
QUIC_SVD quicsvd(dataset, u, v, sigma);
QUIC_SVD<> quicsvd(dataset, u, v, sigma);
}