[mlpack-svn] r14241 - mlpack/branches/mlpack-1.x/src/mlpack/methods/radical
fastlab-svn at coffeetalk-1.cc.gatech.edu
fastlab-svn at coffeetalk-1.cc.gatech.edu
Fri Feb 8 15:38:36 EST 2013
Author: rcurtin
Date: 2013-02-08 15:38:36 -0500 (Fri, 08 Feb 2013)
New Revision: 14241
Modified:
mlpack/branches/mlpack-1.x/src/mlpack/methods/radical/
mlpack/branches/mlpack-1.x/src/mlpack/methods/radical/radical.cpp
mlpack/branches/mlpack-1.x/src/mlpack/methods/radical/radical.hpp
mlpack/branches/mlpack-1.x/src/mlpack/methods/radical/radical_main.cpp
Log:
Style changes for RADICAL.
Property changes on: mlpack/branches/mlpack-1.x/src/mlpack/methods/radical
___________________________________________________________________
Added: svn:mergeinfo
+ /mlpack/trunk/src/mlpack/methods/radical:13981-14240
Modified: mlpack/branches/mlpack-1.x/src/mlpack/methods/radical/radical.cpp
===================================================================
--- mlpack/branches/mlpack-1.x/src/mlpack/methods/radical/radical.cpp 2013-02-08 20:38:18 UTC (rev 14240)
+++ mlpack/branches/mlpack-1.x/src/mlpack/methods/radical/radical.cpp 2013-02-08 20:38:36 UTC (rev 14241)
@@ -14,80 +14,80 @@
// Set the parameters to RADICAL.
Radical::Radical(const double noiseStdDev,
- const size_t nReplicates,
- const size_t nAngles,
- const size_t nSweeps,
+ const size_t replicates,
+ const size_t angles,
+ const size_t sweeps,
const size_t m) :
noiseStdDev(noiseStdDev),
- nReplicates(nReplicates),
- nAngles(nAngles),
- nSweeps(nSweeps),
+ replicates(replicates),
+ angles(angles),
+ sweeps(sweeps),
m(m)
{
// Nothing to do here.
}
-void Radical::CopyAndPerturb(mat& matXNew, const mat& matX)
+void Radical::CopyAndPerturb(mat& xNew, const mat& x) const
{
- matXNew = repmat(matX, nReplicates, 1) + noiseStdDev *
- randn(nReplicates * matX.n_rows, matX.n_cols);
+ Timer::Start("radical_copy_and_perturb");
+ xNew = repmat(x, replicates, 1) + noiseStdDev * randn(replicates * x.n_rows,
+ x.n_cols);
+ Timer::Stop("radical_copy_and_perturb");
}
-double Radical::Vasicek(vec& z)
+double Radical::Vasicek(vec& z) const
{
z = sort(z);
- // Apparently slower
+ // Apparently slower.
/*
- vec logs = log(z.subvec(m, nPoints - 1) - z.subvec(0, nPoints - 1 - m));
+ vec logs = log(z.subvec(m, z.n_elem - 1) - z.subvec(0, z.n_elem - 1 - m));
//vec val = sum(log(z.subvec(m, nPoints - 1) - z.subvec(0, nPoints - 1 - m)));
return (double) sum(logs);
*/
- // Apparently faster
+ // Apparently faster.
double sum = 0;
uword range = z.n_elem - m;
for (uword i = 0; i < range; i++)
{
sum += log(z(i + m) - z(i));
}
+
return sum;
}
double Radical::DoRadical2D(const mat& matX)
{
- mat matXMod;
+ CopyAndPerturb(perturbed, matX);
- CopyAndPerturb(matXMod, matX);
+ mat::fixed<2, 2> matJacobi;
- mat matJacobi(2,2);
- mat candidateY;
+ vec values(angles);
- vec thetas = linspace<vec>(0, nAngles - 1, nAngles) /
- ((double) nAngles) * M_PI / 2;
- vec values(nAngles);
+ for (size_t i = 0; i < angles; i++)
+ {
+ const double theta = (i / (double) angles) * M_PI / 2.0;
+ const double cosTheta = cos(theta);
+ const double sinTheta = sin(theta);
- for (size_t i = 0; i < nAngles; i++)
- {
- double cosTheta = cos(thetas(i));
- double sinTheta = sin(thetas(i));
matJacobi(0, 0) = cosTheta;
matJacobi(1, 0) = -sinTheta;
matJacobi(0, 1) = sinTheta;
matJacobi(1, 1) = cosTheta;
- candidateY = matXMod * matJacobi;
- vec candidateY1 = candidateY.unsafe_col(0);
- vec candidateY2 = candidateY.unsafe_col(1);
+ candidate = perturbed * matJacobi;
+ vec candidateY1 = candidate.unsafe_col(0);
+ vec candidateY2 = candidate.unsafe_col(1);
values(i) = Vasicek(candidateY1) + Vasicek(candidateY2);
}
uword indOpt;
values.min(indOpt); // we ignore the return value; we don't care about it
- return thetas(indOpt);
+ return (indOpt / (double) angles) * M_PI / 2.0;
}
@@ -96,8 +96,10 @@
// matX is nPoints by nDims (although less intuitive than columns being
// points, and although this is the transpose of the ICA literature, this
// choice is for computational efficiency when repeatedly generating
- // two-dimensional coordinate projections for Radical2D.
+ // two-dimensional coordinate projections for Radical2D).
+ Timer::Start("radical_transpose_data");
mat matX = trans(matXT);
+ Timer::Stop("radical_transpose_data");
// If m was not specified, initialize m as recommended in
// (Learned-Miller and Fisher, 2003).
@@ -107,56 +109,71 @@
const size_t nDims = matX.n_cols;
const size_t nPoints = matX.n_rows;
+ Timer::Start("radical_whiten_data");
mat matXWhitened;
mat matWhitening;
WhitenFeatureMajorMatrix(matX, matY, matWhitening);
- // matY is now the whitened form of matX
+ Timer::Stop("radical_whiten_data");
+ // matY is now the whitened form of matX.
// In the RADICAL code, they do not copy and perturb initially, although the
- // paper does. we follow the code as it should match their reported results
+ // paper does. We follow the code as it should match their reported results
// and likely does a better job bouncing out of local optima.
//GeneratePerturbedX(X, X);
// Initialize the unmixing matrix to the whitening matrix.
+ Timer::Start("radical_do_radical");
matW = matWhitening;
mat matYSubspace(nPoints, 2);
- mat matEye = eye(nDims, nDims);
+ mat matJ = eye(nDims, nDims);
- for (size_t sweepNum = 0; sweepNum < nSweeps; sweepNum++)
+ for (size_t sweepNum = 0; sweepNum < sweeps; sweepNum++)
{
+ Log::Info << "RADICAL: sweep " << sweepNum << "." << std::endl;
+
for (size_t i = 0; i < nDims - 1; i++)
{
for (size_t j = i + 1; j < nDims; j++)
{
+ Log::Debug << "RADICAL 2D on dimensions " << i << " and " << j << "."
+ << std::endl;
+
matYSubspace.col(0) = matY.col(i);
matYSubspace.col(1) = matY.col(j);
- double thetaOpt = DoRadical2D(matYSubspace);
- mat matJ = matEye;
- double cosThetaOpt = cos(thetaOpt);
- double sinThetaOpt = sin(thetaOpt);
+
+ const double thetaOpt = DoRadical2D(matYSubspace);
+
+ const double cosThetaOpt = cos(thetaOpt);
+ const double sinThetaOpt = sin(thetaOpt);
+
+ // Set elements of transformation matrix.
matJ(i, i) = cosThetaOpt;
matJ(j, i) = -sinThetaOpt;
matJ(i, j) = sinThetaOpt;
matJ(j, j) = cosThetaOpt;
- matW = matW * matJ;
- matY = matX * matW; // to avoid any issues of mismatch between matW
- // and matY, do not use matY = matY * matJ,
- // even though it may be much more efficient
+
+ matY *= matJ;
+
+ // Unset elements of transformation matrix, so matJ = eye() again.
+ matJ(i, i) = 1;
+ matJ(j, i) = 0;
+ matJ(i, j) = 0;
+ matJ(j, j) = 1;
}
}
}
+ Timer::Stop("radical_do_radical");
- // the final transposes provide W and Y in the typical form from the ICA
- // literature
- //matW = trans(matWhitening * matW);
- //matY = trans(matY);
+ // The final transposes provide W and Y in the typical form from the ICA
+ // literature.
+ Timer::Start("radical_transpose_data");
matW = trans(matW);
matY = trans(matY);
+ Timer::Stop("radical_transpose_data");
}
-
void mlpack::radical::WhitenFeatureMajorMatrix(const mat& matX,
mat& matXWhitened,
mat& matWhitening)
Modified: mlpack/branches/mlpack-1.x/src/mlpack/methods/radical/radical.hpp
===================================================================
--- mlpack/branches/mlpack-1.x/src/mlpack/methods/radical/radical.hpp 2013-02-08 20:38:18 UTC (rev 14240)
+++ mlpack/branches/mlpack-1.x/src/mlpack/methods/radical/radical.hpp 2013-02-08 20:38:36 UTC (rev 14241)
@@ -43,28 +43,28 @@
*
* @param noiseStdDev Standard deviation of the Gaussian noise added to the
* replicates of the data points during Radical2D
- * @param nReplicates Number of Gaussian-perturbed replicates to use
- * (per point) in Radical2D
- * @param nAngles Number of angles to consider in brute-force search during
+ * @param replicates Number of Gaussian-perturbed replicates to use (per
+ * point) in Radical2D
+ * @param angles Number of angles to consider in brute-force search during
* Radical2D
- * @param nSweeps Number of sweeps
- * Each sweep calls Radical2D once for each pair of dimensions
+ * @param sweeps Number of sweeps. Each sweep calls Radical2D once for each
+ * pair of dimensions
* @param m The variable m from Vasicek's m-spacing estimator of entropy.
*/
Radical(const double noiseStdDev = 0.175,
- const size_t nReplicates = 30,
- const size_t nAngles = 150,
- const size_t nSweeps = 0,
+ const size_t replicates = 30,
+ const size_t angles = 150,
+ const size_t sweeps = 0,
const size_t m = 0);
/**
* Run RADICAL.
*
* @param matX Input data into the algorithm - a matrix where each column is
- * a point and each row is a dimension
+ * a point and each row is a dimension.
* @param matY Estimated independent components - a matrix where each column
- * is a point and each row is an estimated independent component
- * @param matW Estimated unmixing matrix, where matY = matW * matX
+ * is a point and each row is an estimated independent component.
+ * @param matW Estimated unmixing matrix, where matY = matW * matX.
*/
void DoRadical(const arma::mat& matX, arma::mat& matY, arma::mat& matW);
@@ -74,44 +74,60 @@
*
* @param x Empirical sample (one-dimensional) over which to estimate entropy.
*/
- double Vasicek(arma::vec& x);
+ double Vasicek(arma::vec& x) const;
/**
- * Make nReplicates copies of each data point and perturb data with Gaussian
+ * Make replicates of each data point (the number of replicates is set in
+ * either the constructor or with Replicates()) and perturb data with Gaussian
* noise with standard deviation noiseStdDev.
*/
- void CopyAndPerturb(arma::mat& matXNew, const arma::mat& matX);
+ void CopyAndPerturb(arma::mat& xNew, const arma::mat& x) const;
//! Two-dimensional version of RADICAL.
double DoRadical2D(const arma::mat& matX);
+ //! Get the standard deviation of the additive Gaussian noise.
+ double NoiseStdDev() const { return noiseStdDev; }
+ //! Modify the standard deviation of the additive Gaussian noise.
+ double& NoiseStdDev() { return noiseStdDev; }
+
+ //! Get the number of Gaussian-perturbed replicates used per point.
+ size_t Replicates() const { return replicates; }
+ //! Modify the number of Gaussian-perturbed replicates used per point.
+ size_t& Replicates() { return replicates; }
+
+ //! Get the number of angles considered during brute-force search.
+ size_t Angles() const { return angles; }
+ //! Modify the number of angles considered during brute-force search.
+ size_t& Angles() { return angles; }
+
+ //! Get the number of sweeps.
+ size_t Sweeps() const { return sweeps; }
+ //! Modify the number of sweeps.
+ size_t& Sweeps() { return sweeps; }
+
private:
- /**
- * standard deviation of the Gaussian noise added to the replicates of
- * the data points during Radical2D
- */
+ //! Standard deviation of the Gaussian noise added to the replicates of
+ //! the data points during Radical2D.
double noiseStdDev;
- /**
- * Number of Gaussian-perturbed replicates to use (per point) in Radical2D
- */
- size_t nReplicates;
+ //! Number of Gaussian-perturbed replicates to use (per point) in Radical2D.
+ size_t replicates;
- /**
- * Number of angles to consider in brute-force search during Radical2D
- */
- size_t nAngles;
+ //! Number of angles to consider in brute-force search during Radical2D.
+ size_t angles;
- /**
- * Number of sweeps
- * - Each sweep calls Radical2D once for each pair of dimensions
- */
- size_t nSweeps;
+ //! Number of sweeps; each sweep calls Radical2D once for each pair of
+ //! dimensions.
+ size_t sweeps;
- /**
- * m to use for Vasicek's m-spacing estimator of entropy
- */
+ //! Value of m to use for Vasicek's m-spacing estimator of entropy.
size_t m;
+
+ //! Internal matrix, held as member variable to prevent memory reallocations.
+ arma::mat perturbed;
+ //! Internal matrix, held as member variable to prevent memory reallocations.
+ arma::mat candidate;
};
void WhitenFeatureMajorMatrix(const arma::mat& matX,
Modified: mlpack/branches/mlpack-1.x/src/mlpack/methods/radical/radical_main.cpp
===================================================================
--- mlpack/branches/mlpack-1.x/src/mlpack/methods/radical/radical_main.cpp 2013-02-08 20:38:18 UTC (rev 14240)
+++ mlpack/branches/mlpack-1.x/src/mlpack/methods/radical/radical_main.cpp 2013-02-08 20:38:36 UTC (rev 14241)
@@ -11,7 +11,8 @@
PROGRAM_INFO("RADICAL", "An implementation of RADICAL, a method for independent"
"component analysis (ICA). Assuming that we have an input matrix X, the"
"goal is to find a square unmixing matrix W such that Y = W * X and the "
- "dimensions of Y are independent components.");
+ "dimensions of Y are independent components. If the algorithm is running"
+ "particularly slowly, try reducing the number of replicates.");
PARAM_STRING_REQ("input_file", "Input dataset filename for ICA.", "i");
@@ -26,9 +27,11 @@
"(per point) in Radical2D.", "r", 30);
PARAM_INT("angles", "Number of angles to consider in brute-force search "
"during Radical2D.", "a", 150);
-PARAM_INT("sweeps", "Number of sweeps (each sweep calls Radical2D once for "
- "each pair of dimensions", "S", 0);
+PARAM_INT("sweeps", "Number of sweeps; each sweep calls Radical2D once for "
+ "each pair of dimensions.", "S", 0);
PARAM_INT("seed", "Random seed. If 0, 'std::time(NULL)' is used.", "s", 0);
+PARAM_FLAG("objective", "If set, an estimate of the final objective function "
+ "is printed.", "O");
using namespace mlpack;
using namespace mlpack::radical;
@@ -69,22 +72,28 @@
mat matW;
rad.DoRadical(matX, matY, matW);
- // save results
+ // Save results.
const string matYFilename = CLI::GetParam<string>("output_ic");
data::Save(matYFilename, matY);
const string matWFilename = CLI::GetParam<string>("output_unmixing");
data::Save(matWFilename, matW);
- /*
- // compute and print objective
- mat matYT = trans(matY);
- double valEst = 0;
- for(size_t i = 0; i < matYT.n_cols; i++)
+ if (CLI::HasParam("objective"))
{
- vec Yi = vec(matYT.col(i));
- valEst += rad.Vasicek(Yi);
+ // Compute and print objective.
+ mat matYT = trans(matY);
+ double valEst = 0;
+ for (size_t i = 0; i < matYT.n_cols; i++)
+ {
+ vec y = vec(matYT.col(i));
+ valEst += rad.Vasicek(y);
+ }
+
+ // Force output even if --verbose is not given.
+ const bool ignoring = Log::Info.ignoreInput;
+ Log::Info.ignoreInput = false;
+ Log::Info << "Objective (estimate): " << valEst << "." << endl;
+ Log::Info.ignoreInput = ignoring;
}
- Log::Info << "Objective (estimate): " << valEst << "." << endl;
- */
}
More information about the mlpack-svn
mailing list