Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Merge pull request #492 from pluskid/rename-LARS
Rename LARS to LeastAngleRegression
  • Loading branch information
Soeren Sonnenburg committed May 2, 2012
2 parents 800be9b + 9935912 commit d587b61
Show file tree
Hide file tree
Showing 5 changed files with 46 additions and 95 deletions.
26 changes: 13 additions & 13 deletions examples/undocumented/python_modular/graphical/regression_lars.py
Expand Up @@ -4,12 +4,12 @@
import matplotlib.pyplot as plt

from shogun.Features import Labels, RealFeatures
from shogun.Regression import LARS, LinearRidgeRegression, LeastSquaresRegression
from shogun.Regression import LeastAngleRegression, LinearRidgeRegression, LeastSquaresRegression

# we compare LASSO with ordinary least-squares (OLE)
# in the ideal case, the MSE of OLE should coincide
# with LASSO at the end of the path
#
#
# if OLE is unstable, we may use RidgeRegression (with
# a small regularization coefficient) to simulate OLE
use_ridge = False
Expand Down Expand Up @@ -42,9 +42,9 @@
y -= np.mean(y)

# train LASSO
lars = LARS()
lars.set_labels(Labels(y))
lars.train(RealFeatures(X.T))
LeastAngleRegression = LeastAngleRegression()
LeastAngleRegression.set_labels(Labels(y))
LeastAngleRegression.train(RealFeatures(X.T))

# train ordinary LSR
if use_ridge:
Expand All @@ -56,24 +56,24 @@
lsr.train(RealFeatures(X.T))

# gather LASSO path
path = np.zeros((p, lars.get_path_size()))
path = np.zeros((p, LeastAngleRegression.get_path_size()))
for i in xrange(path.shape[1]):
path[:,i] = lars.get_w(i)
path[:,i] = LeastAngleRegression.get_w(i)

# apply on training data
mse_train = np.zeros(lars.get_path_size())
mse_train = np.zeros(LeastAngleRegression.get_path_size())
for i in xrange(mse_train.shape[0]):
lars.switch_w(i)
ypred = lars.apply(RealFeatures(X.T)).get_labels()
LeastAngleRegression.switch_w(i)
ypred = LeastAngleRegression.apply(RealFeatures(X.T)).get_labels()
mse_train[i] = np.dot(ypred - y, ypred - y) / y.shape[0]
ypred = lsr.apply(RealFeatures(X.T)).get_labels()
mse_train_lsr = np.dot(ypred - y, ypred - y) / y.shape[0]

# apply on test data
mse_test = np.zeros(lars.get_path_size())
mse_test = np.zeros(LeastAngleRegression.get_path_size())
for i in xrange(mse_test.shape[0]):
lars.switch_w(i)
ypred = lars.apply(RealFeatures(Xtest.T)).get_labels()
LeastAngleRegression.switch_w(i)
ypred = LeastAngleRegression.apply(RealFeatures(Xtest.T)).get_labels()
mse_test[i] = np.dot(ypred - ytest, ypred - ytest) / ytest.shape[0]
ypred = lsr.apply(RealFeatures(Xtest.T)).get_labels()
mse_test_lsr = np.dot(ypred - ytest, ypred - ytest) / ytest.shape[0]
Expand Down
4 changes: 2 additions & 2 deletions src/interfaces/modular/Regression.i
Expand Up @@ -14,7 +14,7 @@
%rename(LinearRidgeRegression) CLinearRidgeRegression;
%rename(LeastSquaresRegression) CLeastSquaresRegression;
%rename(GaussianProcessRegression) CGaussianProcessRegression;
%rename(LARS) CLARS;
%rename(LeastAngleRegression) CLeastAngleRegression;
%rename(LibSVR) CLibSVR;
%rename(MKL) CMKL;
%rename(MKLRegression) CMKLRegression;
Expand All @@ -28,7 +28,7 @@
%include <shogun/regression/LinearRidgeRegression.h>
%include <shogun/regression/LeastSquaresRegression.h>
%include <shogun/regression/GaussianProcessRegression.h>
%include <shogun/regression/LARS.h>
%include <shogun/regression/LeastAngleRegression.h>
%include <shogun/regression/svr/LibSVR.h>
%include <shogun/classifier/mkl/MKL.h>
%include <shogun/regression/svr/MKLRegression.h>
Expand Down
2 changes: 1 addition & 1 deletion src/interfaces/modular/Regression_includes.i
Expand Up @@ -6,7 +6,7 @@
#include <shogun/regression/LinearRidgeRegression.h>
#include <shogun/regression/LeastSquaresRegression.h>
#include <shogun/regression/GaussianProcessRegression.h>
#include <shogun/regression/LARS.h>
#include <shogun/regression/LeastAngleRegression.h>
#include <shogun/classifier/svm/SVM.h>
#include <shogun/classifier/svm/LibSVM.h>
#include <shogun/regression/svr/LibSVR.h>
Expand Down
Expand Up @@ -19,7 +19,7 @@
#include <shogun/features/SimpleFeatures.h>
#include <shogun/mathematics/Math.h>
#include <shogun/mathematics/lapack.h>
#include <shogun/regression/LARS.h>
#include <shogun/regression/LeastAngleRegression.h>

using namespace shogun;
using namespace std;
Expand All @@ -31,7 +31,7 @@ static vector<float64_t> make_vector(int32_t size, float64_t val)
return result;
}

static void plane_rot(float64_t x0, float64_t x1,
static void plane_rot(float64_t x0, float64_t x1,
float64_t &y0, float64_t &y1, SGMatrix<float64_t> &G)
{
G.zero();
Expand All @@ -57,7 +57,7 @@ static void plane_rot(float64_t x0, float64_t x1,
}
}

static void find_max_abs(const vector<float64_t> &vec, const vector<bool> &ignore,
static void find_max_abs(const vector<float64_t> &vec, const vector<bool> &ignore,
int32_t &imax, float64_t& vmax)
{
imax = -1;
Expand All @@ -74,20 +74,20 @@ static void find_max_abs(const vector<float64_t> &vec, const vector<bool> &ignor
}
}

CLARS::CLARS(bool lasso) :
CLeastAngleRegression::CLeastAngleRegression(bool lasso) :
CLinearMachine(), m_lasso(lasso),
m_max_nonz(0), m_max_l1_norm(0)
{
SG_ADD(&m_max_nonz, "max_nonz", "Max number of non-zero variables", MS_AVAILABLE);
SG_ADD(&m_max_l1_norm, "max_l1_norm", "Max l1-norm of estimator", MS_AVAILABLE);
}

CLARS::~CLARS()
CLeastAngleRegression::~CLeastAngleRegression()
{

}

bool CLARS::train_machine(CFeatures* data)
bool CLeastAngleRegression::train_machine(CFeatures* data)
{
if (!m_labels)
SG_ERROR("No labels set\n");
Expand Down Expand Up @@ -133,13 +133,13 @@ bool CLARS::train_machine(CFeatures* data)
for (int32_t j=0; j < n_fea; ++j)
X(i,j) = Xorig(j,i);
Xorig.free_matrix();

// beta is the estimator
vector<float64_t> beta = make_vector(n_fea, 0);

vector<float64_t> Xy = make_vector(n_fea, 0);
// Xy = X' * y
cblas_dgemv(CblasColMajor, CblasTrans, n_vec, n_fea, 1, X.matrix, n_vec,
cblas_dgemv(CblasColMajor, CblasTrans, n_vec, n_fea, 1, X.matrix, n_vec,
y.vector, 1, 0, &Xy[0], 1);

// mu is the prediction
Expand All @@ -161,14 +161,14 @@ bool CLARS::train_machine(CFeatures* data)
m_beta_idx.push_back(0);

//========================================
// main loop
// main loop
//========================================
int32_t nloop=0;
while (m_num_active < n_fea && max_corr > CMath::MACHINE_EPSILON && !stop_cond)
{
// corr = X' * (y-mu) = - X'*mu + Xy
copy(Xy.begin(), Xy.end(), corr.begin());
cblas_dgemv(CblasColMajor, CblasTrans, n_vec, n_fea, -1,
cblas_dgemv(CblasColMajor, CblasTrans, n_vec, n_fea, -1,
X.matrix, n_vec, &mu[0], 1, 1, &corr[0], 1);

// corr_sign = sign(corr)
Expand All @@ -193,7 +193,7 @@ bool CLARS::train_machine(CFeatures* data)

// GA1 = R\(R'\corr_sign_a)
vector<float64_t> GA1(corr_sign_a);
cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans, CblasNonUnit,
cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans, CblasNonUnit,
m_num_active, 1, 1, R.matrix, m_num_active, &GA1[0], m_num_active);
cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit,
m_num_active, 1, 1, R.matrix, m_num_active, &GA1[0], m_num_active);
Expand All @@ -213,12 +213,12 @@ bool CLARS::train_machine(CFeatures* data)
for (int32_t i=0; i < m_num_active; ++i)
{
// u += wA[i] * X[:,m_active_set[i]]
cblas_daxpy(n_vec, wA[i],
cblas_daxpy(n_vec, wA[i],
X.get_column_vector(m_active_set[i]), 1, &u[0], 1);
}

// step size
float64_t gamma = max_corr / AA;
float64_t gamma = max_corr / AA;
if (m_num_active < n_fea)
{
for (int32_t i=0; i < n_fea; ++i)
Expand All @@ -227,7 +227,7 @@ bool CLARS::train_machine(CFeatures* data)
continue;

// correlation between X[:,i] and u
float64_t dir_corr = cblas_ddot(n_vec,
float64_t dir_corr = cblas_ddot(n_vec,
X.get_column_vector(i), 1, &u[0], 1);

float64_t tmp1 = (max_corr-corr[i])/(AA-dir_corr);
Expand Down Expand Up @@ -327,21 +327,7 @@ bool CLARS::train_machine(CFeatures* data)
return true;
}

bool CLARS::load(FILE* srcfile)
{
SG_SET_LOCALE_C;
SG_RESET_LOCALE;
return false;
}

bool CLARS::save(FILE* dstfile)
{
SG_SET_LOCALE_C;
SG_RESET_LOCALE;
return false;
}

void CLARS::cholesky_insert(const SGMatrix<float64_t> &X, SGMatrix<float64_t> &R, int32_t i_max_corr)
void CLeastAngleRegression::cholesky_insert(const SGMatrix<float64_t> &X, SGMatrix<float64_t> &R, int32_t i_max_corr)
{
// diag_k = X[:,i_max_corr]' * X[:,i_max_corr]
float64_t diag_k = cblas_ddot(X.num_rows, X.get_column_vector(i_max_corr), 1,
Expand Down Expand Up @@ -369,10 +355,10 @@ void CLARS::cholesky_insert(const SGMatrix<float64_t> &X, SGMatrix<float64_t> &R

// R' * R_k = (X' * X)_k = col_k, solving to get R_k
vector<float64_t> R_k(col_k);
cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, m_num_active, 1,
cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, m_num_active, 1,
1, R.matrix, m_num_active, &R_k[0], m_num_active);

float64_t R_kk = CMath::sqrt(diag_k -
float64_t R_kk = CMath::sqrt(diag_k -
cblas_ddot(m_num_active, &R_k[0], 1, &R_k[0], 1));

// new_R = [R R_k; zeros(...) R_kk]
Expand All @@ -397,7 +383,7 @@ void CLARS::cholesky_insert(const SGMatrix<float64_t> &X, SGMatrix<float64_t> &R
}
}

void CLARS::cholesky_delete(SGMatrix<float64_t> &R, int32_t i_kick)
void CLeastAngleRegression::cholesky_delete(SGMatrix<float64_t> &R, int32_t i_kick)
{
if (i_kick != m_num_active-1)
{
Expand Down
Expand Up @@ -8,23 +8,23 @@
* Copyright (C) 2012 Chiyuan Zhang
*/

#ifndef _LARS_H__
#define _LARS_H__
#ifndef LEASTANGLEREGRESSION_H__
#define LEASTANGLEREGRESSION_H__

#include <shogun/lib/config.h>

#ifdef HAVE_LAPACK
#include <vector>
#include <shogun/machine/LinearMachine.h>

namespace shogun
namespace shogun
{
class CFeatures;

/** @brief Class for Least Angle Regression, can be used to solve LASSO.
*
* LASSO is basically L1 regulairzed least square regression
*
*
* \f[
* \min \|X^T\beta - y\|^2 + \lambda\|\beta\|_{1}
* \f]
Expand All @@ -48,9 +48,9 @@ class CFeatures;
*
* There is a correspondence between the regularization coefficient lambda
* and the hard constraint constant C. The latter form is easier to control
* by explicitly constraining the l1-norm of the estimator. In this
* by explicitly constraining the l1-norm of the estimator. In this
* implementation, we provide support for the latter form, moreover, we
* allow explicit control of the number of non-zero variables.
* allow explicit control of the number of non-zero variables.
*
* When no constraints is provided, the full path is generated.
*
Expand All @@ -69,17 +69,17 @@ class CFeatures;
* }
* @endcode
*/
class CLARS: public CLinearMachine
class CLeastAngleRegression: public CLinearMachine
{
public:
/** default constructor
*
* @param lasso - when true, it runs the LASSO, when false, it runs LARS
* */
CLARS(bool lasso=true);
CLeastAngleRegression(bool lasso=true);

/** default destructor */
virtual ~CLARS();
virtual ~CLeastAngleRegression();

/** set max number of non-zero variables for early stopping
*
Expand Down Expand Up @@ -142,28 +142,6 @@ class CLARS: public CLinearMachine
return m_beta_idx.size();
}

/** get w
*
* @param num_var number of non-zero coefficients
* @param dst_w store w in this argument
* @param dst_dims dimension of w
*
* **Note** the returned memory references to some internal structures. The
* pointer will become invalid if train is called *again*. So make a copy
* if you want to call train multiple times.
*
* @see switch_w
*/
void get_w(int32_t num_var, float64_t*& dst_w, int32_t& dst_dims)
{
if (w_dim <= 0)
SG_ERROR("cannot get estimator before training");
if (size_t(num_var) >= m_beta_idx.size() || num_var < 0)
SG_ERROR("cannot get an estimator of %d non-zero coefficients", num_var);
dst_dims=w_dim;
dst_w=&m_beta_path[m_beta_idx[num_var]][0];
}

/** get w
*
* @param num_var number of non-zero coefficients
Expand All @@ -176,25 +154,12 @@ class CLARS: public CLinearMachine
SGVector<float64_t> get_w(int32_t num_var)
{
SGVector<float64_t> vec;
get_w(num_var, vec.vector, vec.vlen);
vec.vector = &m_beta_path[m_beta_idx[num_var]][0];
vec.vlen = w_dim;
vec.do_free = false;
return vec;
}

/** load regression from file
*
* @param srcfile file to load from
* @return if loading was successful
*/
virtual bool load(FILE* srcfile);

/** save regression to file
*
* @param dstfile file to save to
* @return if saving was successful
*/
virtual bool save(FILE* dstfile);

/** get classifier type
*
* @return classifier type LinearRidgeRegression
Expand Down Expand Up @@ -244,4 +209,4 @@ class CLARS: public CLinearMachine
} // namespace shogun

#endif // HAVE_LAPACK
#endif // _LARS_H__
#endif // LEASTANGLEREGRESSION_H__

0 comments on commit d587b61

Please sign in to comment.