apaam/mlpack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Put the code in order, and ready to play with the original formulation

Browse Source
This commit is contained in:
vasiloglou
2008-05-22 23:05:00 +00:00
parent 0cd9faa025
commit 4976c5e848
4 changed files with 376 additions and 111 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib2/contrib/nvasil/convex_nmf/main.cc
+5 -5
View File
@@ -5,14 +5,14 @@ int main(int argc, char *argv[]) {
fx_module *fx_root;
fx_root=fx_init(argc, argv, NULL);
fx_module *nmf_module=fx_submodule(fx_root, "/engine");
NmfEngine engine;
NmfEngine<BigSdpNmfObjectiveMinVar> engine;
fx_set_param_str(nmf_module, "data_file",
"/net/hg200/nvasil/dataset/orl_faces/orl_test_faces_100.csv");
fx_set_param_int(nmf_module, "sdp_rank",1 );
"/net/hg200/nvasil/dataset/orl_faces/orl_test_faces_474.csv");
fx_set_param_int(nmf_module, "sdp_rank", 4);
fx_set_param_int(nmf_module, "new_dim", 20);
fx_set_param_double(nmf_module, "l_bfgs/sigma", 0.010);
fx_set_param_double(nmf_module, "l_bfgs/sigma", 0.1);
fx_set_param_double(nmf_module, "l_bfgs/norm_grad_tolerance", 10);
fx_set_param_double(nmf_module, "l_bfgs/desired_feasibility", 1);
fx_set_param_double(nmf_module, "l_bfgs/desired_feasibility", 5);
fx_set_param_double(nmf_module, "l_bfgs/feasibility_tolerance", 0.01);
fx_set_param_int(nmf_module, "l_bfgs/mem_bfgs", 5);
engine.Init(nmf_module);
fastlib2/contrib/nvasil/convex_nmf/nmf_engine.h
+47 -35
View File
@@ -4,6 +4,7 @@
#include "../l_bfgs/l_bfgs.h"
#include "nmf_objectives.h"
template<typename NmfObjective>
class NmfEngine {
public:
void Init(fx_module *module) {
@@ -59,44 +60,55 @@ class NmfEngine {
FATAL("Svd failed soething is wrong...\n");
}
bool negative_flag=false;
index_t positives=0;
index_t negatives=0;
index_t zeros=0;
for (index_t i=0; i<vt_mat.n_cols(); i++) {
if (vt_mat.get(0, i)<0) {
negative_flag=true;
negatives++;
}
if (vt_mat.get(0,i)>0) {
positives++;
}
if (negative_flag==true && vt_mat.get(0, i)>0) {
NONFATAL("Method failed, first eigenvector has positive and negative elements");
break;
}
}
if (negative_flag==true) {
la::Scale(-1, &vt_mat);
}
opt_function_.Project(&vt_mat);
for(index_t i=0; i<num_of_rows_; i++) {
for(index_t j=0; j<new_dim_; j++) {
w_mat_.set(i, j, s[0]*vt_mat.get(0, i*new_dim_+j));
}
}
index_t offset_h=num_of_rows_*new_dim_;
for(index_t i=0; i<num_of_columns_; i++) {
for(index_t j=0; j<new_dim_; j++) {
h_mat_.set(j, i , s[0]*vt_mat.get(0, offset_h+i*new_dim_+j ));
}
}
// now compute reconstruction error
Matrix v_rec;
la::MulInit(w_mat_, h_mat_, &v_rec);
double error=0;
double v_sum=0;
for(index_t i=0; i<values_.size(); i++) {
index_t r=rows_[i];
index_t c=columns_[i];
error+=fabs(v_rec.get(r, c)-values_[i]);
v_sum+=values_[i];
}
data::Save("result.csv", result);
NOTIFY("Reconstruction error: %lg%%\n", error*100/v_sum);
if (vt_mat.get(0, i)==0) {
zeros++;
}
}
data::Save("result.csv", result);
data::Save("vt_mat.csv", vt_mat);
NOTIFY("We found %i positives, %i negatives and %i zeros",
positives, negatives, zeros);
if (negatives>=positives) {
la::Scale(-1, &vt_mat);
}
OptUtils::NonNegativeProjection(&vt_mat);
for(index_t i=0; i<num_of_rows_; i++) {
for(index_t j=0; j<new_dim_; j++) {
w_mat_.set(i, j, s[0]*vt_mat.get(0, i*new_dim_+j));
}
}
index_t offset_h=num_of_rows_*new_dim_;
for(index_t i=0; i<num_of_columns_; i++) {
for(index_t j=0; j<new_dim_; j++) {
h_mat_.set(j, i , s[0]*vt_mat.get(0, offset_h+i*new_dim_+j));
}
}
// now compute reconstruction error
Matrix v_rec;
la::MulInit(w_mat_, h_mat_, &v_rec);
double error=0;
double v_sum=0;
for(index_t i=0; i<values_.size(); i++) {
index_t r=rows_[i];
index_t c=columns_[i];
error+=fabs(v_rec.get(r, c)-values_[i]);
v_sum+=values_[i];
}
NOTIFY("Reconstruction error: %lg%%\n", error*100/v_sum);
}
void GetW(Matrix *w_mat) {
w_mat->Copy(w_mat_);
}
@@ -106,8 +118,8 @@ class NmfEngine {
private:
fx_module *module_;
LBfgs<BigSdpNmfObjective> engine_;
BigSdpNmfObjective opt_function_;
LBfgs<NmfObjective> engine_;
NmfObjective opt_function_;
ArrayList<index_t> rows_;
ArrayList<index_t> columns_;
ArrayList<double> values_;
fastlib2/contrib/nvasil/convex_nmf/nmf_objectives.h
+65 -1
View File
@@ -3,7 +3,7 @@
#include "fastlib/fastlib.h"
#include "../l_bfgs/optimization_utils.h"
class BigSdpNmfObjective {
class BigSdpNmfObjectiveMaxVar {
public:
void Init(fx_module *module,
ArrayList<index_t> &rows,
@@ -33,7 +33,71 @@ class BigSdpNmfObjective {
ArrayList<index_t> columns_;
ArrayList<double> values_;
Vector eq_lagrange_mult_;
};
class BigSdpNmfObjectiveMinVar {
public:
void Init(fx_module *module,
ArrayList<index_t> &rows,
ArrayList<index_t> &columns,
ArrayList<double> &values);
void Destruct();
void ComputeGradient(Matrix &coordinates, Matrix *gradient);
void ComputeObjective(Matrix &coordinates, double *objective);
void ComputeFeasibilityError(Matrix &coordinates, double *error);
double ComputeLagrangian(Matrix &coordinates);
void UpdateLagrangeMult(Matrix &coordinates);
void Project(Matrix *coordinates);
void set_sigma(double sigma);
void GiveInitMatrix(Matrix *init_data);
bool IsDiverging(double objective);
private:
fx_module *module_;
double sigma_;
index_t num_of_columns_;
index_t num_of_rows_;
index_t new_dim_;
index_t rank_;
index_t offset_h_;
index_t offset_h_mat_; // we need this to know wher h matrix starts
ArrayList<index_t> rows_;
ArrayList<index_t> columns_;
ArrayList<double> values_;
Vector eq_lagrange_mult_;
};
class BigSdpNmfObjectiveMinVarIneq {
public:
void Init(fx_module *module,
ArrayList<index_t> &rows,
ArrayList<index_t> &columns,
ArrayList<double> &values);
void Destruct();
void ComputeGradient(Matrix &coordinates, Matrix *gradient);
void ComputeObjective(Matrix &coordinates, double *objective);
void ComputeFeasibilityError(Matrix &coordinates, double *error);
double ComputeLagrangian(Matrix &coordinates);
void UpdateLagrangeMult(Matrix &coordinates);
void Project(Matrix *coordinates);
void set_sigma(double sigma);
void GiveInitMatrix(Matrix *init_data);
bool IsDiverging(double objective);
private:
fx_module *module_;
double sigma_;
index_t num_of_columns_;
index_t num_of_rows_;
index_t new_dim_;
index_t rank_;
index_t offset_h_;
index_t offset_h_mat_; // we need this to know wher h matrix starts
ArrayList<index_t> rows_;
ArrayList<index_t> columns_;
ArrayList<double> values_;
Vector ineqeq_lagrange_mult_;
};
#include "nmf_objectives_impl.h"
fastlib2/contrib/nvasil/convex_nmf/nmf_objectives_impl.h
+259 -70
View File
@@ -1,33 +1,33 @@
void BigSdpNmfObjective::Init(fx_module *module,
void BigSdpNmfObjectiveMaxVar::Init(fx_module *module,
ArrayList<index_t> &rows,
ArrayList<index_t> &columns,
ArrayList<double> &values) {
module_=module;
rows_.InitCopy(rows);
new_dim_=fx_param_int(module_, "new_dim", 5);
rank_=fx_param_int(module_, "rank", 3);
rows_.InitCopy(rows);
columns_.InitCopy(columns);
values_.InitCopy(values);
num_of_rows_=*std::max_element(rows_.begin(), rows_.end())+1;
num_of_columns_=*std::max_element(columns_.begin(), columns_.end())+1;
eq_lagrange_mult_.Init(values_.size());
for(index_t i=0; i<eq_lagrange_mult_.length(); i++) {
eq_lagrange_mult_[i]=math::Random(-0.0, 1.0);
}
rank_=fx_param_int(module_, "rank", 3);
new_dim_=fx_param_int(module_, "new_dim", 5);
eq_lagrange_mult_.SetAll(0);
offset_h_ = num_of_rows_*new_dim_;
}
void BigSdpNmfObjective::Destruct() {
void BigSdpNmfObjectiveMaxVar::Destruct() {
rows_.Renew();
columns_.Renew();
values_.Renew();
eq_lagrange_mult_.Destruct();
}
void BigSdpNmfObjective::ComputeGradient(Matrix &coordinates, Matrix *gradient) {
void BigSdpNmfObjectiveMaxVar::ComputeGradient(Matrix &coordinates, Matrix *gradient) {
// gradient of the variance
gradient->CopyValues(coordinates);
la::Scale(2.0, gradient);
la::Scale(-2.0, gradient);
// Gradient on the equalities and inequalities
for(index_t i=0; i<values_.size(); i++) {
index_t w =rows_[i]*new_dim_;
index_t h =offset_h_+columns_[i]*new_dim_;
@@ -38,19 +38,18 @@ void BigSdpNmfObjective::ComputeGradient(Matrix &coordinates, Matrix *gradient)
coordinates.GetColumnPtr(w), gradient->GetColumnPtr(h));
la::AddExpert(rank_*new_dim_, -eq_lagrange_mult_[i]+2*sigma_*diff,
coordinates.GetColumnPtr(h), gradient->GetColumnPtr(w));
}
}
}
void BigSdpNmfObjective::ComputeObjective(Matrix &coordinates, double *objective) {
void BigSdpNmfObjectiveMaxVar::ComputeObjective(Matrix &coordinates, double *objective) {
*objective=0;
for(index_t i=0; i<coordinates.n_cols(); i++) {
*objective+=la::Dot(coordinates.n_rows(),
for(index_t i=0; i<coordinates.n_cols(); i++) {
*objective-=la::Dot(coordinates.n_rows(),
coordinates.GetColumnPtr(i), coordinates.GetColumnPtr(i));
}
}
void BigSdpNmfObjective::ComputeFeasibilityError(Matrix &coordinates, double *error) {
void BigSdpNmfObjectiveMaxVar::ComputeFeasibilityError(Matrix &coordinates, double *error) {
*error=0;
for(index_t i=0; i<values_.size(); i++) {
index_t w =rows_[i]*new_dim_;
@@ -59,14 +58,14 @@ void BigSdpNmfObjective::ComputeFeasibilityError(Matrix &coordinates, double *er
coordinates.GetColumnPtr(w),
coordinates.GetColumnPtr(h))-values_[i];
*error+=diff*diff;
}
}
}
double BigSdpNmfObjective::ComputeLagrangian(Matrix &coordinates) {
double BigSdpNmfObjectiveMaxVar::ComputeLagrangian(Matrix &coordinates) {
double lagrangian=0;
ComputeObjective(coordinates, &lagrangian);
for(index_t i=0; i<values_.size(); i++) {
// equalities
index_t w =rows_[i]*new_dim_;
index_t h =offset_h_+columns_[i]*new_dim_;
double diff=la::Dot(rank_*new_dim_,
@@ -77,70 +76,27 @@ double BigSdpNmfObjective::ComputeLagrangian(Matrix &coordinates) {
return lagrangian;
}
void BigSdpNmfObjective::UpdateLagrangeMult(Matrix &coordinates) {
void BigSdpNmfObjectiveMaxVar::UpdateLagrangeMult(Matrix &coordinates) {
for(index_t i=0; i<values_.size(); i++) {
index_t w =rows_[i]*new_dim_;
index_t h =offset_h_+columns_[i]*new_dim_;
// equalities
double diff=la::Dot(rank_*new_dim_,
coordinates.GetColumnPtr(w),
coordinates.GetColumnPtr(h))-values_[i];
eq_lagrange_mult_[i]-=sigma_*diff;
}
}
}
void BigSdpNmfObjective::Project(Matrix *coordinates) {
for(index_t i=0; i<coordinates->n_cols(); i++) {
for(index_t j=0; j< coordinates->n_rows(); j++) {
if (coordinates->get(j, i)<0.0) {
coordinates->set(j, i, 0.0);
}
}
}
/* Matrix u_mat, vt_mat;
Vector s;
success_t success=la::SVDInit(*coordinates, &s, &u_mat, &vt_mat);
if (success==SUCCESS_FAIL) {
FATAL("Svd failed...\n");
}
index_t positives=0;
index_t negatives=0;
for(index_t i=0; i<vt_mat.n_cols(); i++) {
if (vt_mat.get(0, i)<0) {
negatives++;
} else {
positives++;
}
}
if (positives>=negatives) {
for(index_t i=0; i<vt_mat.n_cols(); i++) {
if (vt_mat.get(0, i)<0) {
vt_mat.set(0, i, vt_mat.get(0, i));
}
}
} else {
for(index_t i=0; i<vt_mat.n_cols(); i++) {
if (vt_mat.get(0, i)>0) {
vt_mat.set(0, i, vt_mat.get(0, i));
}
}
}
Matrix s_mat;
s_mat.Init(s.length(), s.length());
s_mat.SetDiagonal(s);
Matrix tmp;
coordinates->Destruct();
la::MulInit(u_mat, s_mat, &tmp);
la::MulInit(tmp, vt_mat, coordinates);
*/
void BigSdpNmfObjectiveMaxVar::Project(Matrix *coordinates) {
OptUtils::NonNegativeProjection(coordinates);
}
void BigSdpNmfObjective::set_sigma(double sigma) {
void BigSdpNmfObjectiveMaxVar::set_sigma(double sigma) {
sigma_=sigma;
}
void BigSdpNmfObjective::GiveInitMatrix(Matrix *init_data) {
void BigSdpNmfObjectiveMaxVar::GiveInitMatrix(Matrix *init_data) {
init_data->Init(rank_, (num_of_rows_+num_of_columns_)*new_dim_);
for(index_t i=0; i<init_data->n_rows(); i++) {
for(index_t j=0; j<init_data->n_cols(); j++) {
@@ -150,7 +106,240 @@ void BigSdpNmfObjective::GiveInitMatrix(Matrix *init_data) {
}
bool BigSdpNmfObjective::IsDiverging(double objective) {
bool BigSdpNmfObjectiveMaxVar::IsDiverging(double objective) {
return false;
}
///////////////////////////////////////////////////////////////////////////
void BigSdpNmfObjectiveMinVar::Init(fx_module *module,
ArrayList<index_t> &rows,
ArrayList<index_t> &columns,
ArrayList<double> &values) {
module_=module;
new_dim_=fx_param_int(module_, "new_dim", 5);
rank_=fx_param_int(module_, "rank", 3);
rows_.InitCopy(rows);
columns_.InitCopy(columns);
values_.InitCopy(values);
num_of_rows_=*std::max_element(rows_.begin(), rows_.end())+1;
num_of_columns_=*std::max_element(columns_.begin(), columns_.end())+1;
eq_lagrange_mult_.Init(values_.size());
eq_lagrange_mult_.SetAll(0);
offset_h_ = num_of_rows_*new_dim_;
}
void BigSdpNmfObjectiveMinVar::Destruct() {
rows_.Renew();
columns_.Renew();
values_.Renew();
eq_lagrange_mult_.Destruct();
}
void BigSdpNmfObjectiveMinVar::ComputeGradient(Matrix &coordinates,
Matrix *gradient) {
// gradient of the variance
gradient->CopyValues(coordinates);
la::Scale(2.0, gradient);
// Gradient on the equalities and inequalities
for(index_t i=0; i<values_.size(); i++) {
index_t w =rows_[i]*new_dim_;
index_t h =offset_h_+columns_[i]*new_dim_;
double diff=la::Dot(rank_*new_dim_,
coordinates.GetColumnPtr(w),
coordinates.GetColumnPtr(h))-values_[i];
la::AddExpert(rank_*new_dim_, -eq_lagrange_mult_[i]+2*sigma_*diff,
coordinates.GetColumnPtr(w), gradient->GetColumnPtr(h));
la::AddExpert(rank_*new_dim_, -eq_lagrange_mult_[i]+2*sigma_*diff,
coordinates.GetColumnPtr(h), gradient->GetColumnPtr(w));
}
}
void BigSdpNmfObjectiveMinVar::ComputeObjective(Matrix &coordinates,
double *objective) {
*objective=0;
for(index_t i=0; i<coordinates.n_cols(); i++) {
*objective+=la::Dot(coordinates.n_rows(),
coordinates.GetColumnPtr(i), coordinates.GetColumnPtr(i));
}
}
void BigSdpNmfObjectiveMinVar::ComputeFeasibilityError(Matrix &coordinates,
double *error) {
*error=0;
for(index_t i=0; i<values_.size(); i++) {
index_t w =rows_[i]*new_dim_;
index_t h =offset_h_+columns_[i]*new_dim_;
double diff=la::Dot(rank_*new_dim_,
coordinates.GetColumnPtr(w),
coordinates.GetColumnPtr(h))-values_[i];
*error+=diff*diff;
}
}
double BigSdpNmfObjectiveMinVar::ComputeLagrangian(Matrix &coordinates) {
double lagrangian=0;
ComputeObjective(coordinates, &lagrangian);
for(index_t i=0; i<values_.size(); i++) {
// equalities
index_t w =rows_[i]*new_dim_;
index_t h =offset_h_+columns_[i]*new_dim_;
double diff=la::Dot(rank_*new_dim_,
coordinates.GetColumnPtr(w),
coordinates.GetColumnPtr(h))-values_[i];
lagrangian+=(sigma_*diff-eq_lagrange_mult_[i])*diff;
}
return lagrangian;
}
void BigSdpNmfObjectiveMinVar::UpdateLagrangeMult(Matrix &coordinates) {
for(index_t i=0; i<values_.size(); i++) {
index_t w =rows_[i]*new_dim_;
index_t h =offset_h_+columns_[i]*new_dim_;
// equalities
double diff=la::Dot(rank_*new_dim_,
coordinates.GetColumnPtr(w),
coordinates.GetColumnPtr(h))-values_[i];
eq_lagrange_mult_[i]-=sigma_*diff;
}
}
void BigSdpNmfObjectiveMinVar::Project(Matrix *coordinates) {
OptUtils::NonNegativeProjection(coordinates);
}
void BigSdpNmfObjectiveMinVar::set_sigma(double sigma) {
sigma_=sigma;
}
void BigSdpNmfObjectiveMinVar::GiveInitMatrix(Matrix *init_data) {
init_data->Init(rank_, (num_of_rows_+num_of_columns_)*new_dim_);
for(index_t i=0; i<init_data->n_rows(); i++) {
for(index_t j=0; j<init_data->n_cols(); j++) {
init_data->set(i, j, math::Random());
}
}
}
bool BigSdpNmfObjectiveMinVar::IsDiverging(double objective) {
return false;
}
///////////////////////////////////////////////////////////////////////////////
void BigSdpNmfObjectiveMinVarIneq::Init(fx_module *module,
ArrayList<index_t> &rows,
ArrayList<index_t> &columns,
ArrayList<double> &values) {
module_=module;
new_dim_=fx_param_int(module_, "new_dim", 5);
rank_=fx_param_int(module_, "rank", 3);
rows_.InitCopy(rows);
columns_.InitCopy(columns);
values_.InitCopy(values);
num_of_rows_=*std::max_element(rows_.begin(), rows_.end())+1;
num_of_columns_=*std::max_element(columns_.begin(), columns_.end())+1;
eq_lagrange_mult_.Init(values_.size());
eq_lagrange_mult_.SetAll(0);
offset_h_ = num_of_rows_*new_dim_;
}
void BigSdpNmfObjectiveMinVar::Destruct() {
rows_.Renew();
columns_.Renew();
values_.Renew();
eq_lagrange_mult_.Destruct();
}
void BigSdpNmfObjectiveMinVarIneq::ComputeGradient(Matrix &coordinates,
Matrix *gradient) {
// gradient of the variance
gradient->CopyValues(coordinates);
la::Scale(2.0, gradient);
// Gradient on the equalities and inequalities
for(index_t i=0; i<values_.size(); i++) {
index_t w =rows_[i]*new_dim_;
index_t h =offset_h_+columns_[i]*new_dim_;
double diff=la::Dot(rank_*new_dim_,
coordinates.GetColumnPtr(w),
coordinates.GetColumnPtr(h))-values_[i];
la::AddExpert(rank_*new_dim_, -eq_lagrange_mult_[i]+2*sigma_*diff,
coordinates.GetColumnPtr(w), gradient->GetColumnPtr(h));
la::AddExpert(rank_*new_dim_, -eq_lagrange_mult_[i]+2*sigma_*diff,
coordinates.GetColumnPtr(h), gradient->GetColumnPtr(w));
}
}
void BigSdpNmfObjectiveMinVarIneq::ComputeObjective(Matrix &coordinates,
double *objective) {
*objective=0;
for(index_t i=0; i<coordinates.n_cols(); i++) {
*objective+=la::Dot(coordinates.n_rows(),
coordinates.GetColumnPtr(i), coordinates.GetColumnPtr(i));
}
}
void BigSdpNmfObjectiveMinVarIneq::ComputeFeasibilityError(Matrix &coordinates,
double *error) {
*error=0;
for(index_t i=0; i<values_.size(); i++) {
index_t w =rows_[i]*new_dim_;
index_t h =offset_h_+columns_[i]*new_dim_;
double diff=la::Dot(rank_*new_dim_,
coordinates.GetColumnPtr(w),
coordinates.GetColumnPtr(h))-values_[i];
*error+=diff*diff;
}
}
double BigSdpNmfObjectiveMinVarIneq::ComputeLagrangian(Matrix &coordinates) {
double lagrangian=0;
ComputeObjective(coordinates, &lagrangian);
for(index_t i=0; i<values_.size(); i++) {
// equalities
index_t w =rows_[i]*new_dim_;
index_t h =offset_h_+columns_[i]*new_dim_;
double diff=la::Dot(rank_*new_dim_,
coordinates.GetColumnPtr(w),
coordinates.GetColumnPtr(h))-values_[i];
lagrangian+=(sigma_*diff-eq_lagrange_mult_[i])*diff;
}
return lagrangian;
}
void BigSdpNmfObjectiveMinVarIneq::UpdateLagrangeMult(Matrix &coordinates) {
for(index_t i=0; i<values_.size(); i++) {
index_t w =rows_[i]*new_dim_;
index_t h =offset_h_+columns_[i]*new_dim_;
// equalities
double diff=la::Dot(rank_*new_dim_,
coordinates.GetColumnPtr(w),
coordinates.GetColumnPtr(h))-values_[i];
eq_lagrange_mult_[i]-=sigma_*diff;
}
}
void BigSdpNmfObjectiveMinVarIneq::Project(Matrix *coordinates) {
OptUtils::NonNegativeProjection(coordinates);
}
void BigSdpNmfObjectiveMinVarIneq::set_sigma(double sigma) {
sigma_=sigma;
}
void BigSdpNmfObjectiveMinVarIneq::GiveInitMatrix(Matrix *init_data) {
init_data->Init(rank_, (num_of_rows_+num_of_columns_)*new_dim_);
for(index_t i=0; i<init_data->n_rows(); i++) {
for(index_t j=0; j<init_data->n_cols(); j++) {
init_data->set(i, j, math::Random());
}
}
}
bool BigSdpNmfObjectiveMinVarIneq::IsDiverging(double objective) {
return false;
}