[mlpack-svn] r12376 - mlpack/trunk/src/mlpack/methods/sparse_coding
fastlab-svn at coffeetalk-1.cc.gatech.edu
fastlab-svn at coffeetalk-1.cc.gatech.edu
Fri Apr 13 18:29:38 EDT 2012
Author: rcurtin
Date: 2012-04-13 18:29:38 -0400 (Fri, 13 Apr 2012)
New Revision: 12376
Modified:
mlpack/trunk/src/mlpack/methods/sparse_coding/sparse_coding.cpp
mlpack/trunk/src/mlpack/methods/sparse_coding/sparse_coding.hpp
mlpack/trunk/src/mlpack/methods/sparse_coding/sparse_coding_main.cpp
Log:
Pass one: comment grammar Nazi and fix code style issues.
Modified: mlpack/trunk/src/mlpack/methods/sparse_coding/sparse_coding.cpp
===================================================================
--- mlpack/trunk/src/mlpack/methods/sparse_coding/sparse_coding.cpp 2012-04-13 21:46:25 UTC (rev 12375)
+++ mlpack/trunk/src/mlpack/methods/sparse_coding/sparse_coding.cpp 2012-04-13 22:29:38 UTC (rev 12376)
@@ -2,128 +2,133 @@
* @file sparse_coding.cpp
* @author Nishant Mehta
*
- * Implementation of Sparse Coding with Dictionary Learning using l1 (LASSO) or
- * l1+l2 (Elastic Net) regularization
+ * Implementation of Sparse Coding with Dictionary Learning using l1 (LASSO) or
+ * l1+l2 (Elastic Net) regularization.
*/
-
#include "sparse_coding.hpp"
using namespace std;
using namespace arma;
+using namespace mlpack;
using namespace mlpack::regression;
+using namespace mlpack::sparse_coding;
+// TODO: parameterizable; options to methods?
#define OBJ_TOL 1e-2 // 1E-9
#define NEWTON_TOL 1e-6 // 1E-9
-namespace mlpack {
-namespace sparse_coding {
+SparseCoding::SparseCoding(const mat& matX,
+ uword nAtoms,
+ double lambda1,
+ double lambda2) :
+ nDims(matX.n_rows),
+ nAtoms(nAtoms),
+ nPoints(matX.n_cols),
+ matX(matX),
+ matZ(mat(nAtoms, nPoints)),
+ lambda1(lambda1),
+ lambda2(lambda2)
+{ /* Nothing left to do. */ }
-
-SparseCoding::SparseCoding(const mat& matX, uword nAtoms, double lambda1, double lambda2) :
- nDims(matX.n_rows),
- nAtoms(nAtoms),
- nPoints(matX.n_cols),
- matX(matX),
- matZ(mat(nAtoms, nPoints)),
- lambda1(lambda1),
- lambda2(lambda2)
-{ /* nothing left to do */ }
-
-
-void SparseCoding::SetData(const mat& matX) {
- this -> matX = matX;
+void SparseCoding::SetData(const mat& matX)
+{
+ this->matX = matX;
}
-
-void SparseCoding::SetDictionary(const mat& matD) {
- this -> matD = matD;
+void SparseCoding::SetDictionary(const mat& matD)
+{
+ this->matD = matD;
}
-
-void SparseCoding::InitDictionary() {
+void SparseCoding::InitDictionary()
+{
DataDependentRandomInitDictionary();
}
-
-void SparseCoding::LoadDictionary(const char* dictionaryFilename) {
+void SparseCoding::LoadDictionary(const char* dictionaryFilename)
+{
matD.load(dictionaryFilename);
}
-// always a not good decision!
+// Always a not good decision!
void SparseCoding::RandomInitDictionary() {
matD = randn(nDims, nAtoms);
- for(uword j = 0; j < nAtoms; j++) {
+
+ for (uword j = 0; j < nAtoms; ++j)
matD.col(j) /= norm(matD.col(j), 2);
- }
}
-// the sensible heuristic
-void SparseCoding::DataDependentRandomInitDictionary() {
+// The sensible heuristic.
+void SparseCoding::DataDependentRandomInitDictionary()
+{
matD = mat(nDims, nAtoms);
- for(uword j = 0; j < nAtoms; j++) {
+ for (uword j = 0; j < nAtoms; ++j)
+ {
vec vecD_j = matD.unsafe_col(j);
RandomAtom(vecD_j);
}
}
-
-void SparseCoding::RandomAtom(vec& atom) {
+void SparseCoding::RandomAtom(vec& atom)
+{
atom.zeros();
const uword nSeedAtoms = 3;
- for(uword i = 0; i < nSeedAtoms; i++) {
- atom += matX.col(rand() % nPoints);
- }
+ for(uword i = 0; i < nSeedAtoms; i++)
+ atom += matX.col(rand() % nPoints);
+
atom /= norm(atom, 2);
}
+void SparseCoding::DoSparseCoding(uword nIterations)
+{
+ bool converged = false;
+ double lastObjVal = DBL_MAX;
-void SparseCoding::DoSparseCoding(uword nIterations) {
+ Log::Info << "Initial Coding Step." << endl;
- bool converged = false;
- double lastObjVal = 1e99;
-
- Log::Info << "Initial Coding Step" << endl;
OptimizeCode();
uvec adjacencies = find(matZ);
- Log::Info << "\tSparsity level: "
- << 100.0 * ((double)(adjacencies.n_elem))
- / ((double)(nAtoms * nPoints))
- << "%\n";
- Log::Info << "\tObjective value: " << Objective() << endl;
-
- for(uword t = 1; t <= nIterations && !converged; t++) {
- Log::Info << "Iteration " << t << " of " << nIterations << endl;
- Log::Info << "Dictionary Step\n";
+ Log::Info << " Sparsity level: "
+ << 100.0 * ((double) (adjacencies.n_elem)) / ((double) (nAtoms * nPoints))
+ << "%" << endl;
+ Log::Info << " Objective value: " << Objective() << "." << endl;
+
+ for (uword t = 1; t <= nIterations && !converged; ++t)
+ {
+ Log::Info << "Iteration " << t << " of " << nIterations << "." << endl;
+
+ Log::Info << "Performing dictionary step... ";
OptimizeDictionary(adjacencies);
- //ProjectDictionary(); // is this necessary? solutions to OptimizeDictionary should be feasible
- Log::Info << "\tObjective value: " << Objective() << endl;
+ Log::Info << "objective value: " << Objective() << "." << endl;
- Log::Info << "Coding Step" << endl;
+ Log::Info << "Performing coding step..." << endl;
OptimizeCode();
adjacencies = find(matZ);
- Log::Info << "\tSparsity level: "
- << 100.0 * ((double)(adjacencies.n_elem))
- / ((double)(nAtoms * nPoints))
- << "%\n";
+ Log::Info << " Sparsity level: "
+ << 100.0 *
+ ((double) (adjacencies.n_elem)) / ((double) (nAtoms * nPoints))
+ << "%" << endl;
+
double curObjVal = Objective();
- Log::Info << "\tObjective value: " << curObjVal << endl;
-
+ Log::Info << " Objective value: " << curObjVal << "." << endl;
+
double objValImprov = lastObjVal - curObjVal;
- Log::Info << "\t\t\t\t\tImprovement: " << std::scientific
- << objValImprov << endl;
- if(objValImprov < OBJ_TOL) {
+ Log::Info << " Improvement: " << scientific << objValImprov << "." << endl;
+
+ if (objValImprov < OBJ_TOL)
+ {
converged = true;
- Log::Info << "Converged within tolerance\n";
+ Log::Info << "Converged within tolerance " << OBJ_TOL << ".\n";
}
-
+
lastObjVal = curObjVal;
}
}
-
-void SparseCoding::OptimizeCode() {
- // when using Cholesky version of LARS, this is correct even if lambda2 > 0
+void SparseCoding::OptimizeCode()
+{
+ // When using Cholesky version of LARS, this is correct even if lambda2 > 0.
mat matGram = trans(matD) * matD;
// mat matGram;
// if(lambda2 > 0) {
@@ -132,123 +137,139 @@
// else {
// matGram = trans(matD) * matD;
// }
-
- for(uword i = 0; i < nPoints; i++) {
+
+ for (uword i = 0; i < nPoints; i++)
+ {
// report progress
- if((i % 100) == 0) {
- Log::Debug << "\t" << i << endl;
- }
-
+ if ((i % 100) == 0)
+ Log::Debug << "Optimization at point " << i << "." << endl;
+
bool useCholesky = true;
LARS* lars;
- if(lambda2 > 0) {
+ if(lambda2 > 0)
lars = new LARS(useCholesky, lambda1, lambda2);
- }
- else {
+ else
lars = new LARS(useCholesky, lambda1);
- }
- lars -> SetGramMem(matGram.memptr(), matGram.n_rows);
- lars -> DoLARS(matD, matX.unsafe_col(i));
-
+
+ lars->SetGramMem(matGram.memptr(), matGram.n_rows);
+ lars->DoLARS(matD, matX.unsafe_col(i));
+
vec beta;
- lars -> Solution(beta);
+ lars->Solution(beta);
matZ.col(i) = beta;
+
delete lars;
}
}
-void SparseCoding::OptimizeDictionary(uvec adjacencies) {
- // count number of atomic neighbors for each point x^i
+void SparseCoding::OptimizeDictionary(uvec adjacencies)
+{
+ // Count the number of atomic neighbors for each point x^i.
uvec neighborCounts = zeros<uvec>(nPoints, 1);
- if(adjacencies.n_elem > 0) {
- // this gets the column index
+
+ if (adjacencies.n_elem > 0)
+ {
+ // This gets the column index.
uword curPointInd = (uword)(adjacencies(0) / nAtoms);
uword curCount = 1;
- for(uword l = 1; l < adjacencies.n_elem; l++) {
- if((uword)(adjacencies(l) / nAtoms) == curPointInd) {
- curCount++;
+
+ for (uword l = 1; l < adjacencies.n_elem; ++l)
+ {
+ if ((uword)(adjacencies(l) / nAtoms) == curPointInd)
+ {
+ curCount++;
}
- else {
- neighborCounts(curPointInd) = curCount;
- curPointInd = (uword)(adjacencies(l) / nAtoms);
- curCount = 1;
+ else
+ {
+ neighborCounts(curPointInd) = curCount;
+ curPointInd = (uword)(adjacencies(l) / nAtoms);
+ curCount = 1;
}
}
+
neighborCounts(curPointInd) = curCount;
}
-
- // handle the case of inactive atoms (atoms not used in the given coding)
+
+ // Handle the case of inactive atoms (atoms not used in the given coding).
std::vector<uword> inactiveAtoms;
std::vector<uword> activeAtoms;
activeAtoms.reserve(nAtoms);
- for(uword j = 0; j < nAtoms; j++) {
- if(accu(matZ.row(j) != 0) == 0) {
+
+ for (uword j = 0; j < nAtoms; j++)
+ {
+ if (accu(matZ.row(j) != 0) == 0)
inactiveAtoms.push_back(j);
- }
- else {
+ else
activeAtoms.push_back(j);
- }
}
+
uword nActiveAtoms = activeAtoms.size();
uword nInactiveAtoms = inactiveAtoms.size();
- // efficient construction of Z restricted to active atoms
+ // Efficient construction of Z restricted to active atoms.
mat matActiveZ;
- if(inactiveAtoms.empty()) {
+ if (inactiveAtoms.empty())
+ {
matActiveZ = matZ;
}
- else {
+ else
+ {
uvec inactiveAtomsVec = conv_to< uvec >::from(inactiveAtoms);
RemoveRows(matZ, inactiveAtomsVec, matActiveZ);
}
-
+
uvec atomReverseLookup = uvec(nAtoms);
- for(uword i = 0; i < nActiveAtoms; i++) {
+ for(uword i = 0; i < nActiveAtoms; i++)
+ {
atomReverseLookup(activeAtoms[i]) = i;
}
- if(nInactiveAtoms > 0) {
- Log::Info << "There are " << nInactiveAtoms << " inactive atoms. They will be re-initialized randomly.\n";
+ if (nInactiveAtoms > 0)
+ {
+ Log::Info << "There are " << nInactiveAtoms
+ << " inactive atoms. They will be re-initialized randomly.\n";
}
-
-
- Log::Debug << "Solving Dual via Newton's Method\n";
-
+
+ Log::Debug << "Solving Dual via Newton's Method.\n";
+
mat matDEstimate;
- // solve using Newton's method in the dual - note that the final dot multiplication with inv(A) seems to be unavoidable. Although more expensive, the code written this way (we use solve()) should be more numerically stable than just using inv(A) for everything.
+ // Solve using Newton's method in the dual - note that the final dot
+ // multiplication with inv(A) seems to be unavoidable. Although more
+ // expensive, the code written this way (we use solve()) should be more
+ // numerically stable than just using inv(A) for everything.
vec dualVars = zeros<vec>(nActiveAtoms);
+
//vec dualVars = 1e-14 * ones<vec>(nActiveAtoms);
- //vec dualVars = 10.0 * randu(nActiveAtoms, 1); // method used by feature sign code - fails miserably here. perhaps the MATLAB optimizer fmincon does something clever?
- /*vec dualVars = diagvec(solve(matD, matX * trans(matZ)) - matZ * trans(matZ));
- for(uword i = 0; i < dualVars.n_elem; i++) {
- if(dualVars(i) < 0) {
- dualVars(i) = 0;
- }
- }
- */
- //dualVars.print("dual vars");
+ // Method used by feature sign code - fails miserably here. Perhaps the
+ // MATLAB optimizer fmincon does something clever?
+ //vec dualVars = 10.0 * randu(nActiveAtoms, 1);
+
+ //vec dualVars = diagvec(solve(matD, matX * trans(matZ))
+ // - matZ * trans(matZ));
+ //for (uword i = 0; i < dualVars.n_elem; i++)
+ // if (dualVars(i) < 0)
+ // dualVars(i) = 0;
+
bool converged = false;
mat matZXT = matActiveZ * trans(matX);
mat matZZT = matActiveZ * trans(matActiveZ);
- for(uword t = 1; !converged; t++) {
+
+ for (uword t = 1; !converged; ++t)
+ {
mat A = matZZT + diagmat(dualVars);
-
+
mat matAInvZXT = solve(A, matZXT);
-
- vec gradient = -( sum(square(matAInvZXT), 1) - ones<vec>(nActiveAtoms) );
-
- mat hessian =
- -( -2 * (matAInvZXT * trans(matAInvZXT)) % inv(A) );
-
- //printf("solving for dual variable update...");
+
+ vec gradient = -(sum(square(matAInvZXT), 1) - ones<vec>(nActiveAtoms));
+
+ mat hessian = -(-2 * (matAInvZXT * trans(matAInvZXT)) % inv(A));
+
vec searchDirection = -solve(hessian, gradient);
//vec searchDirection = -gradient;
-
-
- // BEGIN ARMIJO LINE SEARCH
+ // Armijo line search.
const double c = 1e-4;
double alpha = 1.0;
const double rho = 0.9;
@@ -256,164 +277,166 @@
/*
{
- double sumDualVars = sum(dualVars);
- double fOld =
- -( -trace(trans(matZXT) * matAInvZXT) - sumDualVars );
- printf("fOld = %f\t", fOld);
- double fNew =
- -( -trace(trans(matZXT) * solve(matZZT + diagmat(dualVars + alpha * searchDirection), matZXT))
- - (sumDualVars + alpha * sum(searchDirection)) );
- printf("fNew = %f\n", fNew);
+ double sumDualVars = sum(dualVars);
+ double fOld = -(-trace(trans(matZXT) * matAInvZXT) - sumDualVars);
+ Log::Debug << "fOld = " << fOld << "." << endl;
+ double fNew =
+ -(-trace(trans(matZXT) * solve(matZZT +
+ diagmat(dualVars + alpha * searchDirection), matZXT))
+ - (sumDualVars + alpha * sum(searchDirection)) );
+ Log::Debug << "fNew = " << fNew << "." << endl;
}
*/
-
+
double improvement;
- while(true) {
- // objective
+ while (true)
+ {
+ // Calculate objective.
double sumDualVars = sum(dualVars);
- double fOld =
- -( -trace(trans(matZXT) * matAInvZXT) - sumDualVars );
- double fNew =
- -( -trace(trans(matZXT) * solve(matZZT + diagmat(dualVars + alpha * searchDirection), matZXT))
- - (sumDualVars + alpha * sum(searchDirection)) );
+ double fOld = -(-trace(trans(matZXT) * matAInvZXT) - sumDualVars);
+ double fNew = -(-trace(trans(matZXT) * solve(matZZT +
+ diagmat(dualVars + alpha * searchDirection), matZXT)) -
+ (sumDualVars + alpha * sum(searchDirection)));
- // printf("alpha = %e\n", alpha);
- // printf("norm of gradient = %e\n", norm(gradient, 2));
- // printf("sufficientDecrease = %e\n", sufficientDecrease);
- // printf("fNew - fOld - sufficientDecrease = %e\n",
- // fNew - fOld - alpha * sufficientDecrease);
- if(fNew <= fOld + alpha * sufficientDecrease) {
- searchDirection = alpha * searchDirection;
- improvement = fOld - fNew;
- break;
+ if (fNew <= fOld + alpha * sufficientDecrease)
+ {
+ searchDirection = alpha * searchDirection;
+ improvement = fOld - fNew;
+ break;
}
+
alpha *= rho;
}
- // END ARMIJO LINE SEARCH
-
+
+ // End of Armijo line search code.
+
dualVars += searchDirection;
- //printf("\n");
double normGradient = norm(gradient, 2);
- Log::Debug << "Newton Iteration " << t << ":" << endl;
- Log::Debug << "\tnorm of gradient = " << std::scientific << normGradient << endl;
- Log::Debug << "\timprovement = " << std::scientific << improvement << endl;
+ Log::Debug << "Newton Method iteration " << t << ":" << endl;
+ Log::Debug << " Gradient norm: " << std::scientific << normGradient
+ << "." << endl;
+ Log::Debug << " Improvement: " << std::scientific << improvement << ".\n";
- // if(normGradient < NEWTON_TOL) {
- // converged = true;
- // }
- if(improvement < NEWTON_TOL) {
+ if (improvement < NEWTON_TOL)
converged = true;
- }
}
- //dualVars.print("dual solution");
- if(inactiveAtoms.empty()) {
+
+ if (inactiveAtoms.empty())
+ {
matDEstimate = trans(solve(matZZT + diagmat(dualVars), matZXT));
}
- else {
+ else
+ {
mat matDActiveEstimate = trans(solve(matZZT + diagmat(dualVars), matZXT));
matDEstimate = zeros(nDims, nAtoms);
- for(uword i = 0; i < nActiveAtoms; i++) {
+
+ for (uword i = 0; i < nActiveAtoms; i++)
matDEstimate.col(activeAtoms[i]) = matDActiveEstimate.col(i);
- }
- for(uword i = 0; i < nInactiveAtoms; i++) {
+
+ for (uword i = 0; i < nInactiveAtoms; i++)
+ {
vec vecmatDi = matDEstimate.unsafe_col(inactiveAtoms[i]);
RandomAtom(vecmatDi);
}
}
+
matD = matDEstimate;
}
-
-void SparseCoding::ProjectDictionary() {
- for(uword j = 0; j < nAtoms; j++) {
+void SparseCoding::ProjectDictionary()
+{
+ for (uword j = 0; j < nAtoms; j++)
+ {
double normD_j = norm(matD.col(j), 2);
- if(normD_j > 1) {
- if(normD_j - 1.0 > 1e-9) {
- Log::Warn << "Norm Exceeded 1 by " << std::scientific << normD_j - 1.0
- << "\n\tShrinking...\n";
- matD.col(j) /= normD_j;
- }
- // no need to normalize if the dictionary wasn't that infeasible
- //matD.col(j) /= normD_j;
+ if ((normD_j > 1) && (normD_j - 1.0 > 1e-9))
+ {
+ Log::Warn << "Norm exceeded 1 by " << std::scientific << normD_j - 1.0
+ << ". Shrinking...\n";
+ matD.col(j) /= normD_j;
}
}
}
-
-double SparseCoding::Objective() {
+double SparseCoding::Objective()
+{
double l11NormZ = sum(sum(abs(matZ)));
double froNormResidual = norm(matX - matD * matZ, "fro");
- if(lambda2 > 0) {
+
+ if (lambda2 > 0)
+ {
double froNormZ = norm(matZ, "fro");
- return
- 0.5 * (froNormResidual * froNormResidual + lambda2 * froNormZ * froNormZ)
- + lambda1 * l11NormZ;
+ return 0.5 *
+ (froNormResidual * froNormResidual + lambda2 * froNormZ * froNormZ) +
+ lambda1 * l11NormZ;
}
- else {
+ else
+ {
return 0.5 * froNormResidual * froNormResidual + lambda1 * l11NormZ;
}
}
-
-void SparseCoding::PrintDictionary() {
- matD.print("Dictionary");
+void SparseCoding::PrintDictionary()
+{
+ Log::Info << "Dictionary: " << endl << matD;
}
-
-void SparseCoding::PrintCoding() {
- matZ.print("Coding matrix");
+void SparseCoding::PrintCoding()
+{
+ Log::Info << "Coding matrix: " << endl << matZ;
}
-
-void RemoveRows(const mat& X, uvec rows_to_remove, mat& X_mod) {
-
+void mlpack::sparse_coding::RemoveRows(const mat& X,
+ uvec rows_to_remove,
+ mat& X_mod)
+{
uword n_cols = X.n_cols;
uword n_rows = X.n_rows;
uword n_to_remove = rows_to_remove.n_elem;
uword n_to_keep = n_rows - n_to_remove;
-
- if(n_to_remove == 0) {
+
+ if (n_to_remove == 0)
+ {
X_mod = X;
}
- else {
+ else
+ {
X_mod.set_size(n_to_keep, n_cols);
uword cur_row = 0;
uword remove_ind = 0;
- // first, check 0 to first row to remove
- if(rows_to_remove(0) > 0) {
- // note that this implies that n_rows > 1
+ // First, check 0 to first row to remove.
+ if (rows_to_remove(0) > 0)
+ {
+ // Note that this implies that n_rows > 1.
uword height = rows_to_remove(0);
X_mod(span(cur_row, cur_row + height - 1), span::all) =
- X(span(0, rows_to_remove(0) - 1), span::all);
+ X(span(0, rows_to_remove(0) - 1), span::all);
cur_row += height;
}
- // now, check i'th row to remove to (i + 1)'th row to remove, until i = penultimate row
- while(remove_ind < n_to_remove - 1) {
- uword height =
- rows_to_remove[remove_ind + 1]
- - rows_to_remove[remove_ind]
- - 1;
- if(height > 0) {
- X_mod(span(cur_row, cur_row + height - 1),
- span::all) =
- X(span(rows_to_remove[remove_ind] + 1,
- rows_to_remove[remove_ind + 1] - 1),
- span::all);
- cur_row += height;
+ // Now, check i'th row to remove to (i + 1)'th row to remove, until i is the
+ // penultimate row.
+ while (remove_ind < n_to_remove - 1)
+ {
+ uword height = rows_to_remove[remove_ind + 1] - rows_to_remove[remove_ind]
+ - 1;
+
+ if (height > 0)
+ {
+ X_mod(span(cur_row, cur_row + height - 1), span::all) =
+ X(span(rows_to_remove[remove_ind] + 1,
+ rows_to_remove[remove_ind + 1] - 1), span::all);
+ cur_row += height;
}
+
remove_ind++;
}
- // now that i is last row to remove, check last row to remove to last row
- if(rows_to_remove[remove_ind] < n_rows - 1) {
- X_mod(span(cur_row, n_to_keep - 1),
- span::all) =
- X(span(rows_to_remove[remove_ind] + 1, n_rows - 1),
- span::all);
+
+ // Now that i is the last row to remove, check last row to remove to last
+ // row.
+ if (rows_to_remove[remove_ind] < n_rows - 1)
+ {
+ X_mod(span(cur_row, n_to_keep - 1), span::all) =
+ X(span(rows_to_remove[remove_ind] + 1, n_rows - 1), span::all);
}
}
}
-
-
-}; // namespace sparse_coding
-}; // namespace mlpack
Modified: mlpack/trunk/src/mlpack/methods/sparse_coding/sparse_coding.hpp
===================================================================
--- mlpack/trunk/src/mlpack/methods/sparse_coding/sparse_coding.hpp 2012-04-13 21:46:25 UTC (rev 12375)
+++ mlpack/trunk/src/mlpack/methods/sparse_coding/sparse_coding.hpp 2012-04-13 22:29:38 UTC (rev 12376)
@@ -1,15 +1,13 @@
-/**
+/**
* @file sparse_coding.hpp
* @author Nishant Mehta
*
- * Definition of the SparseCoding class, which performs l1 (LASSO) or
- * l1+l2 (Elastic Net)-regularized sparse coding with dictionary learning
+ * Definition of the SparseCoding class, which performs L1 (LASSO) or
+ * L1+L2 (Elastic Net)-regularized sparse coding with dictionary learning
*/
-
#ifndef __MLPACK_METHODS_SPARSE_CODING_SPARSE_CODING_HPP
#define __MLPACK_METHODS_SPARSE_CODING_SPARSE_CODING_HPP
-//#include <armadillo>
#include <mlpack/core.hpp>
#include <mlpack/methods/lars/lars.hpp>
@@ -17,13 +15,14 @@
namespace sparse_coding {
/**
- * An implementation of Sparse Coding with Dictionary Learning that achieves
- * sparsity via an l1-norm regularizer on the codes (LASSO) or an (l1+l2)-norm
+ * An implementation of Sparse Coding with Dictionary Learning that achieves
+ * sparsity via an l1-norm regularizer on the codes (LASSO) or an (l1+l2)-norm
* regularizer on the codes (the Elastic Net).
- * Let d be the number of dimensions in the original space, m the number of
- * training points, and k the number of atoms in the dictionary (the dimension
- * of the learned feature space). The training data X is a d-by-m matrix where
- * each column is a point and each row is a dimension. The dictionary D is a
+ *
+ * Let d be the number of dimensions in the original space, m the number of
+ * training points, and k the number of atoms in the dictionary (the dimension
+ * of the learned feature space). The training data X is a d-by-m matrix where
+ * each column is a point and each row is a dimension. The dictionary D is a
* d-by-k matrix, and the sparse codes matrix Z is a k-by-m matrix.
* This program seeks to minimize the objective:
* min_{D,Z} 0.5 ||X - D Z||_{Fro}^2\ + lambda_1 sum_{i=1}^m ||Z_i||_1
@@ -32,10 +31,10 @@
* where typically lambda_1 > 0 and lambda_2 = 0.
*
* This problem is solved by an algorithm that alternates between a dictionary
- * learning step and a sparse coding step. The dictionary learning step updates
- * the dictionary D using a Newton method based on the Lagrange dual (see the
- * paper below for details). The sparse coding step involves solving a large
- * number of sparse linear regression problems; this can be done efficiently
+ * learning step and a sparse coding step. The dictionary learning step updates
+ * the dictionary D using a Newton method based on the Lagrange dual (see the
+ * paper below for details). The sparse coding step involves solving a large
+ * number of sparse linear regression problems; this can be done efficiently
* using LARS, an algorithm that can solve the LASSO or the Elastic Net (papers below).
*
* Here are those papers:
@@ -79,15 +78,10 @@
* }
* @endcode
*/
-class SparseCoding {
+class SparseCoding
+{
public:
- // void Init(double* memX, uword nDims, uword nPoints,
- // uword nAtoms, double lambda1);
-
- //void SetDictionary(double* memD);
-
-
/**
* Set the parameters to SparseCoding. lambda2 defaults to 0.
*
@@ -96,93 +90,89 @@
* @param lambda1 Regularization parameter for l1-norm penalty
* @param lambda2 Regularization parameter for l2-norm penalty
*/
- SparseCoding(const arma::mat& matX, arma::uword nAtoms, double lambda1, double lambda2 = 0);
-
+ SparseCoding(const arma::mat& matX,
+ arma::uword nAtoms,
+ double lambda1,
+ double lambda2 = 0);
/**
- * Initialize dictionary somehow
+ * Initialize dictionary.
*/
void InitDictionary();
-
- /**
+
+ /**
* Load dictionary from a file
- *
+ *
* @param dictionaryFilename Filename containing dictionary
*/
void LoadDictionary(const char* dictionaryFilename);
-
+
/**
- * Initialize dictionary by drawing k vectors uniformly at random from the
- * unit sphere
+ * Initialize dictionary by drawing k vectors uniformly at random from the
+ * unit sphere.
*/
void RandomInitDictionary();
/**
* Initialize dictionary by initializing each atom to a normalized mixture of
- * a small number of randomly selected points in X
+ * a small number of randomly selected points in X.
*/
void DataDependentRandomInitDictionary();
/**
* Initialize an atom to a normalized mixture of a small number of randomly
- * selected points in X
+ * selected points in X.
*
* @param atom The atom to initialize
*/
void RandomAtom(arma::vec& atom);
-
// core algorithm functions
/**
- * Run Sparse Coding with Dictionary Learning
+ * Run Sparse Coding with Dictionary Learning.
*
- * @param nIterations Maximum number of iterations to run algorithm
+ * @param nIterations Maximum number of iterations to run algorithm.
*/
void DoSparseCoding(arma::uword nIterations);
/**
- * Sparse code each point via LARS
+ * Sparse code each point via LARS.
*/
void OptimizeCode();
-
- /**
+
+ /**
* Learn dictionary via Newton method based on Lagrange dual
*
- * @param adjacencies Indices of entries (unrolled column by column) of
- * the coding matrix Z that are non-zero (the adjacency matrix for the
- * bipartite graph of points and atoms)
+ * @param adjacencies Indices of entries (unrolled column by column) of
+ * the coding matrix Z that are non-zero (the adjacency matrix for the
+ * bipartite graph of points and atoms).
*/
void OptimizeDictionary(arma::uvec adjacencies);
/**
- * Project each atom of the dictionary onto the unit ball
+ * Project each atom of the dictionary onto the unit ball.
*/
void ProjectDictionary();
/**
- * Compute objective function
+ * Compute objective function.
*/
double Objective();
// accessors, modifiers, printers
- //! Modifier for matX
+ //! Modifier for matX.
void SetData(const arma::mat& matX);
- //! Modifier for matD
+ //! Modifier for matD.
void SetDictionary(const arma::mat& matD);
-
- //! Accessor for matD
- const arma::mat& MatD() {
- return matD;
- }
-
+
+ //! Accessor for matD.
+ const arma::mat& MatD() { return matD; }
//! Accessor for matZ
- const arma::mat& MatZ() {
- return matZ;
- }
+ const arma::mat& MatZ() { return matZ; }
// Print the dictionary matD
void PrintDictionary();
@@ -191,30 +181,33 @@
void PrintCoding();
private:
+ //! Number of dimensions.
arma::uword nDims;
+ //! Number of atoms.
arma::uword nAtoms;
+ //! Number of points.
arma::uword nPoints;
- // data (columns are points)
+ //! data (columns are points).
arma::mat matX;
- // dictionary (columns are atoms)
- arma::mat matD;
-
- // sparse codes (columns are points)
- arma::mat matZ;
-
- // l1 regularization term
- double lambda1;
-
- // l2 regularization term
- double lambda2;
-
-};
+ //! dictionary (columns are atoms).
+ arma::mat matD;
-void RemoveRows(const arma::mat& X, arma::uvec rows_to_remove, arma::mat& X_mod);
+ //! Sparse codes (columns are points).
+ arma::mat matZ;
+ //! l1 regularization term.
+ double lambda1;
+ //! l2 regularization term.
+ double lambda2;
+};
+
+void RemoveRows(const arma::mat& X,
+ arma::uvec rows_to_remove,
+ arma::mat& X_mod);
+
}; // namespace sparse_coding
}; // namespace mlpack
Modified: mlpack/trunk/src/mlpack/methods/sparse_coding/sparse_coding_main.cpp
===================================================================
--- mlpack/trunk/src/mlpack/methods/sparse_coding/sparse_coding_main.cpp 2012-04-13 21:46:25 UTC (rev 12375)
+++ mlpack/trunk/src/mlpack/methods/sparse_coding/sparse_coding_main.cpp 2012-04-13 22:29:38 UTC (rev 12376)
@@ -1,109 +1,117 @@
-/** @file sparse_coding_main.cpp
- * @author Nishant Mehta
+/**
+ * @file sparse_coding_main.cpp
+ * @author Nishant Mehta
*
- * Executable for Sparse Coding
+ * Executable for Sparse Coding.
*/
-
#include <mlpack/core.hpp>
#include "sparse_coding.hpp"
PROGRAM_INFO("Sparse Coding", "An implementation of l1-norm and l1+l2-norm "
- "regularized Sparse Coding with Dictionary Learning");
+ "regularized Sparse Coding with Dictionary Learning");
-PARAM_DOUBLE_REQ("lambda1", "sparse coding l1-norm regularization parameter.", "l");
-PARAM_DOUBLE("lambda2", "sparse coding l2-norm regularization parameter.", "", 0);
+PARAM_DOUBLE_REQ("lambda1", "Sparse coding l1-norm regularization parameter.",
+ "l");
+PARAM_DOUBLE("lambda2", "Sparse coding l2-norm regularization parameter.", "L",
+ 0);
-PARAM_INT_REQ("n_atoms", "number of atoms in dictionary.", "k");
+PARAM_INT_REQ("n_atoms", "Number of atoms in the dictionary.", "k");
-PARAM_INT_REQ("n_iterations", "number of iterations for sparse coding.", "");
+PARAM_INT_REQ("n_iterations", "Number of iterations for sparse coding.", "n");
-PARAM_STRING_REQ("data", "path to the input data.", "");
-PARAM_STRING("initial_dictionary", "Filename for initial dictionary.", "", "");
-PARAM_STRING("results_dir", "Directory for results.", "", "");
+PARAM_STRING_REQ("input_file", "Filename of the input data.", "i");
+PARAM_STRING("initial_dictionary", "Filename for optional initial dictionary.", "d", "");
+PARAM_STRING("results_dir", "Directory to store results in.", "r", "");
+PARAM_INT("seed", "Random seed. If 0, 'std::time(NULL) is used.", "s", 0);
-
-
using namespace arma;
using namespace std;
using namespace mlpack;
using namespace mlpack::sparse_coding;
-int main(int argc, char* argv[]) {
+int main(int argc, char* argv[])
+{
CLI::ParseCommandLine(argc, argv);
-
- std::srand(time(NULL));
-
+
+ if (CLI::GetParam<int>("seed") != 0)
+ math::RandomSeed((size_t) CLI::GetParam<int>("seed"));
+ else
+ math::RandomSeed((size_t) std::time(NULL));
+
double lambda1 = CLI::GetParam<double>("lambda1");
double lambda2 = CLI::GetParam<double>("lambda2");
-
- // if using fx-run, one could just leave results_dir blank
+
const char* resultsDir = CLI::GetParam<string>("results_dir").c_str();
-
const char* dataFullpath = CLI::GetParam<string>("data").c_str();
+ const char* initialDictionaryFullpath =
+ CLI::GetParam<string>("initial_dictionary").c_str();
- const char* initialDictionaryFullpath =
- CLI::GetParam<string>("initial_dictionary").c_str();
-
size_t nIterations = CLI::GetParam<int>("n_iterations");
-
+
size_t nAtoms = CLI::GetParam<int>("n_atoms");
-
+
mat matX;
data::Load(dataFullpath, matX);
-
+
uword nPoints = matX.n_cols;
- printf("Loaded %d points in %d dimensions\n", nPoints, matX.n_rows);
+ Log::Info << "Loaded " << nPoints << " points in " << matX.n_rows <<
+ " dimensions." << endl;
- // normalize each point since these are images
- for(uword i = 0; i < nPoints; i++) {
+ // Normalize each point since these are images.
+ for (uword i = 0; i < nPoints; ++i)
matX.col(i) /= norm(matX.col(i), 2);
- }
-
- // run Sparse Coding
+
+ // Run the sparse coding algorithm.
SparseCoding sc(matX, nAtoms, lambda1, lambda2);
-
- if(strlen(initialDictionaryFullpath) == 0) {
+
+ if (strlen(initialDictionaryFullpath) == 0)
+ {
sc.DataDependentRandomInitDictionary();
}
- else {
+ else
+ {
mat matInitialD;
data::Load(initialDictionaryFullpath, matInitialD);
- if(matInitialD.n_cols != nAtoms) {
- Log::Fatal << "The specified initial dictionary to load has "
- << matInitialD.n_cols << " atoms, but the learned dictionary "
- << "was specified to have " << nAtoms << " atoms!\n";
- return EXIT_FAILURE;
+
+ if (matInitialD.n_cols != nAtoms)
+ {
+ Log::Fatal << "The specified initial dictionary to load has "
+ << matInitialD.n_cols << " atoms, but the learned dictionary "
+ << "was specified to have " << nAtoms << " atoms!" << endl;
}
- if(matInitialD.n_rows != matX.n_rows) {
+
+ if (matInitialD.n_rows != matX.n_rows)
+ {
Log::Fatal << "The specified initial dictionary to load has "
- << matInitialD.n_rows << " dimensions, but the specified data "
- << "has " << matX.n_rows << " dimensions!\n";
- return EXIT_FAILURE;
+ << matInitialD.n_rows << " dimensions, but the specified data "
+ << "has " << matX.n_rows << " dimensions!" << endl;
}
+
sc.SetDictionary(matInitialD);
}
-
-
+
Timer::Start("sparse_coding");
sc.DoSparseCoding(nIterations);
- Timer::Stop("sparse_coding");
-
+ Timer::Stop("sparse_coding");
+
mat learnedD = sc.MatD();
- mat learnedZ = sc.MatZ();
+ mat learnedZ = sc.MatZ();
- if(strlen(resultsDir) == 0) {
+ if (strlen(resultsDir) == 0)
+ {
data::Save("D.csv", learnedD);
data::Save("Z.csv", learnedZ);
}
- else {
- char* dataFullpath = (char*) malloc(320 * sizeof(char));
+ else
+ {
+ stringstream datapath;
+ datapath << resultsDir << "/D.csv";
- sprintf(dataFullpath, "%s/D.csv", resultsDir);
- data::Save(dataFullpath, learnedD);
-
- sprintf(dataFullpath, "%s/Z.csv", resultsDir);
- data::Save(dataFullpath, learnedZ);
-
- free(dataFullpath);
+ data::Save(datapath.str(), learnedD);
+
+ datapath.clear();
+ datapath << resultsDir << "/Z.csv";
+
+ data::Save(datapath.str(), learnedZ);
}
}
More information about the mlpack-svn
mailing list