[mlpack-svn] r16154 - mlpack/trunk/src/mlpack/methods/lars
fastlab-svn at coffeetalk-1.cc.gatech.edu
fastlab-svn at coffeetalk-1.cc.gatech.edu
Wed Jan 22 11:36:07 EST 2014
Author: rcurtin
Date: Wed Jan 22 11:36:06 2014
New Revision: 16154
Log:
Fix from Michael Fox for #242.
Modified:
mlpack/trunk/src/mlpack/methods/lars/lars.cpp
mlpack/trunk/src/mlpack/methods/lars/lars.hpp
Modified: mlpack/trunk/src/mlpack/methods/lars/lars.cpp
==============================================================================
--- mlpack/trunk/src/mlpack/methods/lars/lars.cpp (original)
+++ mlpack/trunk/src/mlpack/methods/lars/lars.cpp Wed Jan 22 11:36:06 2014
@@ -57,6 +57,9 @@
// (all dimensions are inactive).
isActive.resize(dataRef.n_cols, false);
+ // Set up ignores set variables. Initialized empty.
+ isIgnored.resize(dataRef.n_cols, false);
+
// Initialize yHat and beta.
beta = arma::zeros(dataRef.n_cols);
arma::vec yHat = arma::zeros(dataRef.n_rows);
@@ -100,13 +103,14 @@
}
// Main loop.
- while ((activeSet.size() < dataRef.n_cols) && (maxCorr > tolerance))
+ while (((activeSet.size() + ignoreSet.size()) < dataRef.n_cols) &&
+ (maxCorr > tolerance))
{
// Compute the maximum correlation among inactive dimensions.
maxCorr = 0;
for (size_t i = 0; i < dataRef.n_cols; i++)
{
- if ((!isActive[i]) && (fabs(corr(i)) > maxCorr))
+ if ((!isActive[i]) && (!isIgnored[i]) && (fabs(corr(i)) > maxCorr))
{
maxCorr = fabs(corr(i));
changeInd = i;
@@ -146,22 +150,41 @@
arma::vec betaDirection;
if (useCholesky)
{
- /**
- * Note that:
- * R^T R % S^T % S = (R % S)^T (R % S)
- * Now, for 1 the ones vector:
- * inv( (R % S)^T (R % S) ) 1
- * = inv(R % S) inv((R % S)^T) 1
- * = inv(R % S) Solve((R % S)^T, 1)
- * = inv(R % S) Solve(R^T, s)
- * = Solve(R % S, Solve(R^T, s)
- * = s % Solve(R, Solve(R^T, s))
- */
- unnormalizedBetaDirection = solve(trimatu(matUtriCholFactor),
- solve(trimatl(trans(matUtriCholFactor)), s));
+ // Check for singularity.
+ const double lastUtriElement = matUtriCholFactor(
+ matUtriCholFactor.n_cols - 1, matUtriCholFactor.n_rows - 1);
+ if (std::abs(lastUtriElement) > tolerance)
+ {
+ // Ok, no singularity.
+ /**
+ * Note that:
+ * R^T R % S^T % S = (R % S)^T (R % S)
+ * Now, for 1 the ones vector:
+ * inv( (R % S)^T (R % S) ) 1
+ * = inv(R % S) inv((R % S)^T) 1
+ * = inv(R % S) Solve((R % S)^T, 1)
+ * = inv(R % S) Solve(R^T, s)
+ * = Solve(R % S, Solve(R^T, s)
+ * = s % Solve(R, Solve(R^T, s))
+ */
+ unnormalizedBetaDirection = solve(trimatu(matUtriCholFactor),
+ solve(trimatl(trans(matUtriCholFactor)), s));
- normalization = 1.0 / sqrt(dot(s, unnormalizedBetaDirection));
- betaDirection = normalization * unnormalizedBetaDirection;
+ normalization = 1.0 / sqrt(dot(s, unnormalizedBetaDirection));
+ betaDirection = normalization * unnormalizedBetaDirection;
+ }
+ else
+ {
+ // Singularity, so remove variable from active set, add to ignores set,
+ // and look for new variable to add.
+ Log::Warn << "Encountered singularity when adding variable "
+ << changeInd << " to active set; permanently removing."
+ << std::endl;
+ Deactivate(activeSet.size() - 1);
+ Ignore(changeInd);
+ CholeskyDelete(matUtriCholFactor.n_rows - 1);
+ continue;
+ }
}
else
{
@@ -170,11 +193,28 @@
for (size_t j = 0; j < activeSet.size(); j++)
matGramActive(i, j) = matGram(activeSet[i], activeSet[j]);
+ // Check for singularity.
arma::mat matS = s * arma::ones<arma::mat>(1, activeSet.size());
- unnormalizedBetaDirection = solve(matGramActive % trans(matS) % matS,
+ const bool solvedOk = solve(unnormalizedBetaDirection,
+ matGramActive % trans(matS) % matS,
arma::ones<arma::mat>(activeSet.size(), 1));
- normalization = 1.0 / sqrt(sum(unnormalizedBetaDirection));
- betaDirection = normalization * unnormalizedBetaDirection % s;
+ if (solvedOk)
+ {
+ // Ok, no singularity.
+ normalization = 1.0 / sqrt(sum(unnormalizedBetaDirection));
+ betaDirection = normalization * unnormalizedBetaDirection % s;
+ }
+ else
+ {
+ // Singularity, so remove variable from active set, add to ignores set,
+ // and look for new variable to add.
+ Deactivate(activeSet.size() - 1);
+ Ignore(changeInd);
+ Log::Warn << "Encountered singularity when adding variable "
+ << changeInd << " to active set; permanently removing."
+ << std::endl;
+ continue;
+ }
}
// compute "equiangular" direction in output space
@@ -183,12 +223,12 @@
double gamma = maxCorr / normalization;
// If not all variables are active.
- if (activeSet.size() < dataRef.n_cols)
+ if ((activeSet.size() + ignoreSet.size()) < dataRef.n_cols)
{
// Compute correlations with direction.
for (size_t ind = 0; ind < dataRef.n_cols; ind++)
{
- if (isActive[ind])
+ if (isActive[ind] || isIgnored[ind])
continue;
double dirCorr = dot(dataRef.col(ind), yHatDirection);
@@ -295,6 +335,12 @@
activeSet.push_back(varInd);
}
+void LARS::Ignore(const size_t varInd)
+{
+ isIgnored[varInd] = true;
+ ignoreSet.push_back(varInd);
+}
+
void LARS::ComputeYHatDirection(const arma::mat& matX,
const arma::vec& betaDirection,
arma::vec& yHatDirection)
Modified: mlpack/trunk/src/mlpack/methods/lars/lars.hpp
==============================================================================
--- mlpack/trunk/src/mlpack/methods/lars/lars.hpp (original)
+++ mlpack/trunk/src/mlpack/methods/lars/lars.hpp Wed Jan 22 11:36:06 2014
@@ -189,6 +189,14 @@
//! Active set membership indicator (for each dimension).
std::vector<bool> isActive;
+ // Set of variables that are ignored (if any).
+
+ //! Set of ignored variables (for dimensions in span{active set dimensions}).
+ std::vector<size_t> ignoreSet;
+
+ //! Membership indicator for set of ignored variables.
+ std::vector<bool> isIgnored;
+
/**
* Remove activeVarInd'th element from active set.
*
@@ -203,6 +211,13 @@
*/
void Activate(const size_t varInd);
+ /**
+ * Add dimension varInd to ignores set (never removed).
+ *
+ * @param varInd Dimension to add to ignores set.
+ */
+ void Ignore(const size_t varInd);
+
// compute "equiangular" direction in output space
void ComputeYHatDirection(const arma::mat& matX,
const arma::vec& betaDirection,
More information about the mlpack-svn
mailing list