[mlpack-svn] r10701 - mlpack/trunk/src/mlpack/methods/lars
fastlab-svn at coffeetalk-1.cc.gatech.edu
fastlab-svn at coffeetalk-1.cc.gatech.edu
Sun Dec 11 14:12:53 EST 2011
Author: niche
Date: 2011-12-11 14:12:53 -0500 (Sun, 11 Dec 2011)
New Revision: 10701
Modified:
mlpack/trunk/src/mlpack/methods/lars/lars.cpp
mlpack/trunk/src/mlpack/methods/lars/lars.hpp
Log:
many style fixes, made some functions private, commented out some unused functions
Modified: mlpack/trunk/src/mlpack/methods/lars/lars.cpp
===================================================================
--- mlpack/trunk/src/mlpack/methods/lars/lars.cpp 2011-12-10 21:33:21 UTC (rev 10700)
+++ mlpack/trunk/src/mlpack/methods/lars/lars.cpp 2011-12-11 19:12:53 UTC (rev 10701)
@@ -7,9 +7,12 @@
#include "lars.hpp"
-//using namespace std; // we are explicit with std:: to avoid confusing std::vector with arma::vec
-using namespace arma; // we use this too often to explicitly use arma:: everywhere
+// we are explicit with std:: to avoid confusing std::vector with arma::vec
+// we use arma namespace too often to explicitly use arma:: everywhere
+//using namespace std;
+using namespace arma;
+
namespace mlpack {
namespace lars {
@@ -58,9 +61,13 @@
void LARS::ComputeGram()
{
if (elasticNet)
+ {
matGram = trans(matX) * matX + lambda2 * eye(matX.n_cols, matX.n_cols);
+ }
else
+ {
matGram = trans(matX) * matX;
+ }
}
@@ -69,7 +76,7 @@
matXTy = trans(matX) * y;
}
-
+ /*
void LARS::UpdateX(const std::vector<int>& colInds, const mat& matNewCols)
{
for (u32 i = 0; i < colInds.size(); i++)
@@ -109,21 +116,28 @@
i != colInds.end(); ++i)
matXTy(*i) = dot(matX.col(*i), y);
}
+ */
+/*
void LARS::PrintGram()
{
matGram.print("Gram matrix");
}
+*/
+ /*
void LARS::SetY(const vec& y)
{
this->y = y;
}
+ */
+ /*
void LARS::PrintY()
{
y.print();
}
+ */
const std::vector<u32> LARS::ActiveSet()
{
@@ -140,17 +154,21 @@
return lambdaPath;
}
+ /*
void LARS::SetDesiredLambda(double lambda1)
{
this->lambda1 = lambda1;
}
+ */
void LARS::DoLARS()
{
// compute Gram matrix, XtY, and initialize active set varibles
ComputeXty();
if (!useCholesky && matGram.is_empty())
+ {
ComputeGram();
+ }
// set up active set variables
nActive = 0;
@@ -216,7 +234,9 @@
{
vec newGramCol = vec(nActive);
for (u32 i = 0; i < nActive; i++)
+ {
newGramCol[i] = dot(matX.col(activeSet[i]), matX.col(changeInd));
+ }
CholeskyInsert(matX.col(changeInd), newGramCol);
}
@@ -228,7 +248,9 @@
// compute signs of correlations
vec s = vec(nActive);
for (u32 i = 0; i < nActive; i++)
+ {
s(i) = corr(activeSet[i]) / fabs(corr(activeSet[i]));
+ }
// compute "equiangular" direction in parameter space (betaDirection)
/* We use quotes because in the case of non-unit norm variables,
@@ -249,8 +271,8 @@
* = Solve(R % S, Solve(R^T, s)
* = s % Solve(R, Solve(R^T, s))
*/
- unnormalizedBetaDirection = solve(trimatu(utriCholFactor),
- solve(trimatl(trans(utriCholFactor)), s));
+ unnormalizedBetaDirection = solve(trimatu(matUtriCholFactor),
+ solve(trimatl(trans(matUtriCholFactor)), s));
normalization = 1.0 / sqrt(dot(s, unnormalizedBetaDirection));
betaDirection = normalization * unnormalizedBetaDirection;
@@ -266,9 +288,9 @@
}
}
- mat S = s * ones<mat>(1, nActive);
+ mat matS = s * ones<mat>(1, nActive);
unnormalizedBetaDirection =
- solve(matGramActive % trans(S) % S, ones<mat>(nActive, 1));
+ solve(matGramActive % trans(matS) % matS, ones<mat>(nActive, 1));
normalization = 1.0 / sqrt(sum(unnormalizedBetaDirection));
betaDirection = normalization * unnormalizedBetaDirection % s;
}
@@ -294,9 +316,13 @@
double val1 = (maxCorr - corr(ind)) / (normalization - dirCorr);
double val2 = (maxCorr + corr(ind)) / (normalization + dirCorr);
if ((val1 > 0) && (val1 < gamma))
- gamma = val1;
+ {
+ gamma = val1;
+ }
if((val2 > 0) && (val2 < gamma))
- gamma = val2;
+ {
+ gamma = val2;
+ }
}
}
@@ -389,11 +415,14 @@
beta = BetaPath().back();
}
-void LARS::GetCholFactor(mat& matR)
+void LARS::GetCholFactor(mat& matUtriCholFactor)
{
- matR = utriCholFactor;
+ matUtriCholFactor = this->matUtriCholFactor;
}
+
+ // private functions
+
void LARS::Deactivate(u32 activeVarInd)
{
nActive--;
@@ -413,7 +442,9 @@
{
yHatDirection.fill(0);
for(u32 i = 0; i < nActive; i++)
+ {
yHatDirection += betaDirection(i) * matX.col(activeSet[i]);
+ }
}
void LARS::InterpolateBeta()
@@ -434,13 +465,17 @@
void LARS::CholeskyInsert(const vec& newX, const mat& X)
{
- if (utriCholFactor.n_rows == 0)
+ if (matUtriCholFactor.n_rows == 0)
{
- utriCholFactor = mat(1, 1);
+ matUtriCholFactor = mat(1, 1);
if (elasticNet)
- utriCholFactor(0, 0) = sqrt(dot(newX, newX) + lambda2);
+ {
+ matUtriCholFactor(0, 0) = sqrt(dot(newX, newX) + lambda2);
+ }
else
- utriCholFactor(0, 0) = norm(newX, 2);
+ {
+ matUtriCholFactor(0, 0) = norm(newX, 2);
+ }
}
else
{
@@ -450,15 +485,19 @@
}
void LARS::CholeskyInsert(const vec& newX, const vec& newGramCol) {
- int n = utriCholFactor.n_rows;
+ int n = matUtriCholFactor.n_rows;
if (n == 0)
{
- utriCholFactor = mat(1, 1);
+ matUtriCholFactor = mat(1, 1);
if (elasticNet)
- utriCholFactor(0, 0) = sqrt(dot(newX, newX) + lambda2);
+ {
+ matUtriCholFactor(0, 0) = sqrt(dot(newX, newX) + lambda2);
+ }
else
- utriCholFactor(0, 0) = norm(newX, 2);
+ {
+ matUtriCholFactor(0, 0) = norm(newX, 2);
+ }
}
else
{
@@ -466,41 +505,45 @@
double sqNormNewX;
if (elasticNet)
+ {
sqNormNewX = dot(newX, newX) + lambda2;
+ }
else
+ {
sqNormNewX = dot(newX, newX);
+ }
- vec utriCholFactork = solve(trimatl(trans(utriCholFactor)),
+ vec matUtriCholFactork = solve(trimatl(trans(matUtriCholFactor)),
newGramCol);
- matNewR(span(0, n - 1), span(0, n - 1)) = utriCholFactor;
- matNewR(span(0, n - 1), n) = utriCholFactork;
+ matNewR(span(0, n - 1), span(0, n - 1)) = matUtriCholFactor;
+ matNewR(span(0, n - 1), n) = matUtriCholFactork;
matNewR(n, span(0, n - 1)).fill(0.0);
- matNewR(n, n) = sqrt(sqNormNewX - dot(utriCholFactork, utriCholFactork));
+ matNewR(n, n) = sqrt(sqNormNewX - dot(matUtriCholFactork, matUtriCholFactork));
- utriCholFactor = matNewR;
+ matUtriCholFactor = matNewR;
}
}
-void LARS::GivensRotate(const vec& x, vec& rotatedX, mat& G)
+void LARS::GivensRotate(const vec& x, vec& rotatedX, mat& matG)
{
if (x(1) == 0)
{
- G = eye(2, 2);
+ matG = eye(2, 2);
rotatedX = x;
}
else
{
double r = norm(x, 2);
- G = mat(2, 2);
+ matG = mat(2, 2);
double scaledX1 = x(0) / r;
double scaledX2 = x(1) / r;
- G(0, 0) = scaledX1;
- G(1, 0) = -scaledX2;
- G(0, 1) = scaledX2;
- G(1, 1) = scaledX1;
+ matG(0, 0) = scaledX1;
+ matG(1, 0) = -scaledX2;
+ matG(0, 1) = scaledX2;
+ matG(1, 1) = scaledX1;
rotatedX = vec(2);
rotatedX(0) = r;
@@ -510,30 +553,30 @@
void LARS::CholeskyDelete(u32 colToKill)
{
- u32 n = utriCholFactor.n_rows;
+ u32 n = matUtriCholFactor.n_rows;
if (colToKill == (n - 1))
{
- utriCholFactor = utriCholFactor(span(0, n - 2), span(0, n - 2));
+ matUtriCholFactor = matUtriCholFactor(span(0, n - 2), span(0, n - 2));
}
else
{
- utriCholFactor.shed_col(colToKill); // remove column colToKill
+ matUtriCholFactor.shed_col(colToKill); // remove column colToKill
n--;
for(u32 k = colToKill; k < n; k++)
{
- mat G;
+ mat matG;
vec rotatedVec;
- GivensRotate(utriCholFactor(span(k, k + 1), k), rotatedVec, G);
- utriCholFactor(span(k, k + 1), k) = rotatedVec;
+ GivensRotate(matUtriCholFactor(span(k, k + 1), k), rotatedVec, matG);
+ matUtriCholFactor(span(k, k + 1), k) = rotatedVec;
if (k < n - 1)
{
- utriCholFactor(span(k, k + 1), span(k + 1, n - 1)) = G *
- utriCholFactor(span(k, k + 1), span(k + 1, n - 1));
+ matUtriCholFactor(span(k, k + 1), span(k + 1, n - 1)) =
+ matG * matUtriCholFactor(span(k, k + 1), span(k + 1, n - 1));
}
}
- utriCholFactor.shed_row(n);
+ matUtriCholFactor.shed_row(n);
}
}
Modified: mlpack/trunk/src/mlpack/methods/lars/lars.hpp
===================================================================
--- mlpack/trunk/src/mlpack/methods/lars/lars.hpp 2011-12-10 21:33:21 UTC (rev 10700)
+++ mlpack/trunk/src/mlpack/methods/lars/lars.hpp 2011-12-11 19:12:53 UTC (rev 10701)
@@ -19,32 +19,51 @@
// beta is the estimator
// yHat is the prediction from the current estimator
+/**
+ * An implementation of LARS, a stage-wise homotopy-based algorithm for
+ * l1 regularized linear regression (LASSO) and l1+l2 regularized linear
+ * regression (Elastic Net).
+ * Let X be a matrix where each row is a point and each column is a dimension,
+ * and let y be a vector of targets.
+ * The Elastic Net problem is to solve
+ * min_beta ||X beta - y||_2^2 + lambda_1 ||beta||_1 + 0.5 lambda_2 ||beta||_2^2
+ * If lambda_1 > 0, lambda_2 = 0, the problem is the LASSO.
+ * If lambda_1 > 0, lambda_2 > 0, the problem is the Elastic Net.
+ * If lambda_1 = 0, lambda_2 > 0, the problem is Ridge Regression.
+ * If lambda_1 = 0, lambda_2 = 0, the problem is unregularized linear regression.
+ * Note: This algorithm is not recommended for use (in terms of efficiency)
+ * when lambda_1 = 0.
+ *
+ * For more details, see the following papers:
+ *
+ *
+ * @article{efron2004least,
+ * title={Least angle regression},
+ * author={Efron, B. and Hastie, T. and Johnstone, I. and Tibshirani, R.},
+ * journal={The Annals of statistics},
+ * volume={32},
+ * number={2},
+ * pages={407--499},
+ * year={2004},
+ * publisher={Institute of Mathematical Statistics}
+ * }
+ *
+ *
+ * @article{zou2005regularization,
+ * title={Regularization and variable selection via the elastic net},
+ * author={Zou, H. and Hastie, T.},
+ * journal={Journal of the Royal Statistical Society Series B},
+ * volume={67},
+ * number={2},
+ * pages={301--320},
+ * year={2005},
+ * publisher={Royal Statistical Society}
+ * }
+ */
class LARS {
- private:
- arma::mat matX;
- arma::vec y;
- arma::vec matXTy;
- arma::mat matGram;
- //! Upper triangular cholesky factor; initially 0x0 arma::matrix.
- arma::mat utriCholFactor;
-
- bool useCholesky;
-
- bool lasso;
- double lambda1;
-
- bool elasticNet;
- double lambda2;
-
- std::vector<arma::vec> betaPath;
- std::vector<double> lambdaPath;
-
- arma::u32 nActive;
- std::vector<arma::u32> activeSet;
- std::vector<bool> isActive;
-
public:
+
LARS(const arma::mat& matX,
const arma::vec& y,
const bool useCholesky);
@@ -68,32 +87,59 @@
void ComputeXty();
- void UpdateX(const std::vector<int>& colInds, const arma::mat& matNewCols);
+ //void UpdateX(const std::vector<int>& colInds, const arma::mat& matNewCols);
- void UpdateGram(const std::vector<int>& colInds);
+ //void UpdateGram(const std::vector<int>& colInds);
- void UpdateXty(const std::vector<int>& colInds);
+ //void UpdateXty(const std::vector<int>& colInds);
- void PrintGram();
+ //void PrintGram();
- void SetY(const arma::vec& y);
+ //void SetY(const arma::vec& y);
- void PrintY();
+ //void PrintY();
const std::vector<arma::u32> ActiveSet();
const std::vector<arma::vec> BetaPath();
const std::vector<double> LambdaPath();
+
+ //void SetDesiredLambda(double lambda1);
- void SetDesiredLambda(double lambda1);
-
void DoLARS();
void Solution(arma::vec& beta);
- void GetCholFactor(arma::mat& matR);
+ void GetCholFactor(arma::mat& matUtriCholFactor);
+
+
+private:
+ arma::mat matX;
+ arma::vec y;
+
+ arma::vec matXTy;
+ arma::mat matGram;
+ //! Upper triangular cholesky factor; initially 0x0 arma::matrix.
+ arma::mat matUtriCholFactor;
+
+ bool useCholesky;
+
+ bool lasso;
+ double lambda1;
+
+ bool elasticNet;
+ double lambda2;
+
+ std::vector<arma::vec> betaPath;
+ std::vector<double> lambdaPath;
+
+ arma::u32 nActive;
+ std::vector<arma::u32> activeSet;
+ std::vector<bool> isActive;
+
+
void Deactivate(arma::u32 activeVarInd);
void Activate(arma::u32 varInd);
@@ -110,6 +156,8 @@
void GivensRotate(const arma::vec& x, arma::vec& rotatedX, arma::mat& G);
void CholeskyDelete(arma::u32 colToKill);
+
+
};
}; // namespace lars
More information about the mlpack-svn
mailing list