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

the bug in the gradient was found now it is working but it seems to me that the minimum has a lot of points

Browse Source
This commit is contained in:
vasiloglou
2008-05-27 06:14:27 +00:00
parent ce471f9889
commit ef2fdfb8c3
4 changed files with 97 additions and 51 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fastlib2/contrib/nvasil/convex_nmf/main.cc
+9 -8
View File
@@ -24,15 +24,16 @@ int main(int argc, char *argv[]) {
fx_root=fx_init(argc, argv, NULL);
fx_module *nmf_module=fx_submodule(fx_root, "/engine");
SdpNmfEngine<SmallSdpNmf> engine;
fx_set_param_str(nmf_module, "data_file", "1.csv");
fx_set_param_str(nmf_module, "data_file", "5.csv");
// "/net/hg200/nvasil/dataset/orl_faces/orl_test_faces_100.csv");
fx_set_param_int(nmf_module, "new_dim", 1);
fx_set_param_int(nmf_module, "new_dim", 2);
fx_set_param_double(nmf_module, "l_bfgs/sigma", 1);
fx_set_param_double(nmf_module, "l_bfgs/gamma", 1.7);
fx_set_param_double(nmf_module, "l_bfgs/wolfe_sigma2", 0.9999);
fx_set_param_double(nmf_module, "l_bfgs/gamma", 2);
fx_set_param_double(nmf_module, "l_bfgs/wolfe_sigma2", 0.99);
fx_set_param_double(nmf_module, "l_bfgs/wolfe_sigma1", 1e-4);
fx_set_param_double(nmf_module, "l_bfgs/use_default_termination", false);
fx_set_param_int(nmf_module, "l_bfgs/mem_bfgs", 10);
fx_set_param_double(nmf_module, "optfun/desired_duality_gap", 1e-3);
fx_set_param_int(nmf_module, "l_bfgs/mem_bfgs", 4);
fx_set_param_double(nmf_module, "optfun/desired_duality_gap", 1e-4);
fx_set_param_double(nmf_module, "optfun/gradient_tolerance", 1e-2);
engine.Init(nmf_module);
engine.ComputeNmf();
@@ -40,6 +41,6 @@ int main(int argc, char *argv[]) {
Matrix h_mat;
engine.GetW(&w_mat);
engine.GetH(&h_mat);
data::Save("W.csv", w_mat);
data::Save("H.csv", h_mat);
data::Save("w.csv", w_mat);
data::Save("h.csv", h_mat);
}
fastlib2/contrib/nvasil/convex_nmf/sdp_nmf_engine.h
+35 -13
View File
@@ -25,6 +25,14 @@
template<typename SdpNmfObjective>
class SdpNmfEngine {
public:
SdpNmfEngine() {
engine_=NULL;
}
~SdpNmfEngine() {
if (engine_!=NULL) {
delete engine_;
}
}
void Init(fx_module *module) {
module_=module;
std::string data_file=fx_param_str_req(module_, "data_file");
@@ -40,21 +48,36 @@ class SdpNmfEngine {
fx_module *opt_function_module=fx_submodule(module_, "optfun");
fx_set_param_int(opt_function_module, "new_dim", new_dim_);
opt_function_.Init(opt_function_module, rows_, columns_, values_);
fx_module *l_bfgs_module=fx_submodule(module_, "l_bfgs");
Matrix init_data;
opt_function_.GiveInitMatrix(&init_data);
fx_set_param_int(l_bfgs_module, "num_of_points", init_data.n_cols());
fx_set_param_int(l_bfgs_module, "new_dimension", init_data.n_rows());
engine_.Init(&opt_function_, l_bfgs_module);
engine_.set_coordinates(init_data);
l_bfgs_module_=fx_submodule(module_, "l_bfgs");
}
void Destruct() {
if (engine_!=NULL) {
delete engine_;
engine_ = NULL;
}
SdpNmfObjective opt_function_;
rows_.Renew();
columns_.Renew();
values_.Renew();
w_mat_.Destruct();
h_mat_.Destruct();
};
void ComputeNmf() {
engine_.ComputeLocalOptimumBFGS();
Matrix result;
engine_.GetResults(&result);
opt_function_.GiveInitMatrix(&result);
fx_set_param_int(l_bfgs_module_, "num_of_points", result.n_cols());
fx_set_param_int(l_bfgs_module_, "new_dimension", result.n_rows());
for(double sigma=1.0; sigma<=1e4; sigma*=4) {
engine_=new LBfgs<SdpNmfObjective>();
engine_->Init(&opt_function_, l_bfgs_module_);
engine_->set_coordinates(result);
result.Destruct();
engine_->set_sigma(sigma);
engine_->ComputeLocalOptimumBFGS();
engine_->GetResults(&result);
}
w_mat_.Init(new_dim_, num_of_rows_);
h_mat_.Init(new_dim_, num_of_columns_);
w_mat_.CopyColumnFromMat(0, 0, num_of_rows_, result);
@@ -72,8 +95,6 @@ class SdpNmfEngine {
error+=fabs(v_rec.get(r, c)-values_[i]);
v_sum+=values_[i];
}
data::Save("w.csv", w_mat_);
data::Save("h.csv", h_mat_);
NOTIFY("Reconstruction error: %lg%%\n", error*100/v_sum);
}
@@ -86,7 +107,8 @@ class SdpNmfEngine {
private:
fx_module *module_;
LBfgs<SdpNmfObjective> engine_;
fx_module *l_bfgs_module_;
LBfgs<SdpNmfObjective> *engine_;
SdpNmfObjective opt_function_;
ArrayList<index_t> rows_;
ArrayList<index_t> columns_;
fastlib2/contrib/nvasil/convex_nmf/sdp_objectives_impl.h
+41 -22
View File
@@ -41,32 +41,32 @@ void SmallSdpNmf::Init(fx_module *module,
// W
for(index_t i=0; i<num_of_rows_; i++) {
for(index_t j=0; j<new_dim_; j++) {
objective_factor_.set(j, i, -1);
objective_factor_.set(j, i, -1.0);
}
}
// H
for(index_t i=offset_h_; i<offset_h_+num_of_columns_; i++) {
for(index_t j=0; j<new_dim_; j++) {
objective_factor_.set(j, i, -1);
objective_factor_.set(j, i, -1.0);
}
}
// tw
for(index_t i=offset_tw_; i<offset_tw_+num_of_rows_; i++) {
for(index_t j=0; j<new_dim_; j++) {
objective_factor_.set(j, i, 1);
objective_factor_.set(j, i, 1.0);
}
}
// th
for(index_t i=offset_th_; i<offset_th_+num_of_rows_; i++) {
for(index_t j=0; j<new_dim_; j++) {
objective_factor_.set(j, i, 1);
objective_factor_.set(j, i, 1.0);
}
}
// v
for(index_t i=offset_v_; i<offset_v_+values_.size(); i++) {
for(index_t j=0; j<new_dim_; j++) {
objective_factor_.set(j, i, 1);
objective_factor_.set(j, i, 1.0);
}
}
@@ -84,6 +84,7 @@ void SmallSdpNmf::ComputeGradient(Matrix &coordinates, Matrix *gradient) {
diff+=coordinates.get(j, v_i);
}
diff-=values_[i];
//NOTIFY("diff:%lg", diff);
for(index_t j=0; j<new_dim_; j++) {
gradient->set(j, v_i, gradient->get(j, v_i)-1.0/diff);
}
@@ -110,16 +111,18 @@ void SmallSdpNmf::ComputeGradient(Matrix &coordinates, Matrix *gradient) {
double t1_minus_ww=(t1-w*w);
double t2_minus_hh=(t2-h*h);
double wh_minus_v=(w*h-v);
double determinant=t1_minus_ww*t2_minus_hh-math::Pow<2,1>(wh_minus_v);
NOTIFY("determinant:%lg", determinant);
DEBUG_ERR_MSG_IF(determinant==0.0, "Determinant equal to zero");
if (determinant==0.0) {
double determinant=t1_minus_ww*t2_minus_hh
-math::Pow<2,1>(wh_minus_v);
//NOTIFY("determinant:%lg", determinant);
if (determinant<=0.0) {
gradient->set(j, w_i, -DBL_MAX);
gradient->set(j, h_i, -DBL_MAX);
gradient->set(j, t1_i, -DBL_MAX);
gradient->set(j, t2_i, -DBL_MAX);
gradient->set(j, v_i, -DBL_MAX);
NONFATAL("Determinant is zero");
NONFATAL("Determinant is less or equal to zero");
return;
} else {
double dw=(-2*w*(t2_minus_hh)-2*h*(wh_minus_v))/determinant;
double dh=(-2*h*(t1_minus_ww)-2*w*(wh_minus_v))/determinant;
@@ -129,7 +132,7 @@ void SmallSdpNmf::ComputeGradient(Matrix &coordinates, Matrix *gradient) {
gradient->set(j, w_i, gradient->get(j, w_i)-dw);
gradient->set(j, h_i, gradient->get(j, h_i)-dh);
gradient->set(j, t1_i, gradient->get(j, t1_i)-dt1);
gradient->set(j, t2_i, gradient->get(j, t1_i)-dt2);
gradient->set(j, t2_i, gradient->get(j, t2_i)-dt2);
gradient->set(j, v_i, gradient->get(j, v_i)-dv);
}
}
@@ -144,6 +147,10 @@ void SmallSdpNmf::ComputeObjective(Matrix &coordinates, double *objective) {
void SmallSdpNmf::ComputeFeasibilityError(Matrix &coordinates, double *error) {
// return duality gap instead
*error=number_of_cones_/sigma_;
double lagrangian=ComputeLagrangian(coordinates);
double objective;
ComputeObjective(coordinates, &objective);
NOTIFY("sum of log barriers:%lg ", lagrangian-sigma_*objective);
}
double SmallSdpNmf::ComputeLagrangian(Matrix &coordinates) {
@@ -152,7 +159,9 @@ double SmallSdpNmf::ComputeLagrangian(Matrix &coordinates) {
ComputeObjective(coordinates, &lagrangian);
lagrangian*=sigma_;
// from the LP cones
double temp_prod=1.0;
for(index_t i=0; i<values_.size(); i++) {
index_t v_i=offset_v_+i;
double diff=0;
for(index_t j=0; j<new_dim_; j++) {
@@ -162,13 +171,12 @@ double SmallSdpNmf::ComputeLagrangian(Matrix &coordinates) {
if unlikely(diff<0) {
return DBL_MAX;
}
// DEBUG_ERR_MSG_IF(diff<=0, "LP cone is invalid, you are "
// " out of the feasible region, constraint %i, diff %lg",
// i, diff);
lagrangian-=log(diff+1e-100);
}
temp_prod*=diff;
}
lagrangian-=log(temp_prod);
// from the SDP cones
// determinant=(t1-w*w)*(t2-h*h)-math::Pow<2,1>(w*h-v);
temp_prod=1.0;
for(index_t i=0; i<values_.size(); i++) {
index_t w_i=rows_[i];
index_t h_i=offset_h_+columns_[i];
@@ -183,16 +191,25 @@ double SmallSdpNmf::ComputeLagrangian(Matrix &coordinates) {
double v=coordinates.get(j, v_i);
double t1_minus_ww=(t1-w*w);
double t2_minus_hh=(t2-h*h);
// we need this, because during Wolfe step
// we might accidently cross the barrier
// and never realize it since the determinant
// will stay positive
if (t1_minus_ww<0 || t2_minus_hh<0) {
return DBL_MAX;
}
double wh_minus_v=(w*h-v);
double determinant=t1_minus_ww*t2_minus_hh-math::Pow<2,1>(wh_minus_v);
if (unlikely(determinant<0)) {
double determinant=t1_minus_ww*t2_minus_hh
-math::Pow<2,1>(wh_minus_v);
if (unlikely(determinant<=0)) {
return DBL_MAX;
}
// DEBUG_ERR_MSG_IF(determinant<=0, "SDP cone is invalid, you are "
// " out of the feasible region");
lagrangian-=log(determinant+1e-100);
temp_prod*=determinant;
}
}
lagrangian-=log(temp_prod);
return lagrangian;
}
@@ -201,7 +218,7 @@ void SmallSdpNmf::UpdateLagrangeMult(Matrix &coordinates) {
}
void SmallSdpNmf::Project(Matrix *coordinates) {
OptUtils::NonNegativeProjection(coordinates);
//OptUtils::NonNegativeProjection(coordinates);
}
void SmallSdpNmf::set_sigma(double sigma) {
@@ -213,7 +230,7 @@ void SmallSdpNmf::GiveInitMatrix(Matrix *init_data) {
2*(num_of_rows_+num_of_columns_)+values_.size());
for(index_t i=0; i<num_of_rows_+num_of_columns_; i++) {
for(index_t j=0; j<new_dim_; j++) {
init_data->set(j, i, math::Random(0.0, 1.0));
init_data->set(j, i, 0*math::Random(0.0, 1.0));
}
}
for(index_t i=0; i<values_.size(); i++) {
@@ -244,7 +261,7 @@ bool SmallSdpNmf::IsDiverging(double objective) {
bool SmallSdpNmf::IsOptimizationOver(Matrix &coordinates,
Matrix &gradient, double step) {
double norm_gradient = la::Dot(gradient.n_elements(),
/* double norm_gradient = la::Dot(gradient.n_elements(),
gradient.ptr(), gradient.ptr());
// one of our barriers is zero
if (norm_gradient>=DBL_MAX) {
@@ -255,6 +272,8 @@ bool SmallSdpNmf::IsOptimizationOver(Matrix &coordinates,
} else {
return false;
}
*/
return true;
}
bool SmallSdpNmf::IsIntermediateStepOver(Matrix &coordinates,
fastlib2/contrib/nvasil/l_bfgs/l_bfgs_impl.h
+12 -8
View File
@@ -40,7 +40,7 @@ void LBfgs<OptimizedFunction>::Init(OptimizedFunction *optimized_function,
}
DEBUG_ASSERT(wolfe_sigma1_>0);
DEBUG_ASSERT(wolfe_sigma1_<1);
DEBUG_ASSERT(wolfe_sigma2_>0);
DEBUG_ASSERT(wolfe_sigma2_>wolfe_sigma1_);
DEBUG_ASSERT(wolfe_sigma2_<1);
if (unlikely(wolfe_sigma1_>=wolfe_sigma2_)) {
FATAL("Wolfe sigma1 %lg should be less than sigma2 %lg",
@@ -48,7 +48,7 @@ void LBfgs<OptimizedFunction>::Init(OptimizedFunction *optimized_function,
}
DEBUG_ASSERT(wolfe_sigma1_>0);
DEBUG_ASSERT(wolfe_sigma1_<1);
DEBUG_ASSERT(wolfe_sigma2_>0);
DEBUG_ASSERT(wolfe_sigma2_>wolfe_sigma1_);
DEBUG_ASSERT(wolfe_sigma2_<1);
norm_grad_tolerance_ = fx_param_double(module_, "norm_grad_tolerance", 0.1);
wolfe_beta_ = fx_param_double(module_, "wolfe_beta", 0.8);
@@ -108,7 +108,7 @@ void LBfgs<OptimizedFunction>::ComputeLocalOptimumBFGS() {
success_t success=ComputeBFGS_(&step_, gradient_, i);
if (success==SUCCESS_FAIL) {
NOTIFY("LBFGS failed to find a direction, continuing with gradient descent\n");
ComputeWolfeStep_(&step_, gradient_);
//ComputeWolfeStep_(&step_, gradient_);
}
optimized_function_->ComputeGradient(coordinates_, &gradient_);
UpdateBFGS_();
@@ -118,7 +118,7 @@ void LBfgs<OptimizedFunction>::ComputeLocalOptimumBFGS() {
if (silent_==false) {
ReportProgressFile_();
}
if (use_default_termination_== true) {
/* if (use_default_termination_== true) {
optimized_function_->ComputeFeasibilityError(coordinates_,
&feasibility_error);
if (feasibility_error < desired_feasibility_) {
@@ -132,6 +132,7 @@ void LBfgs<OptimizedFunction>::ComputeLocalOptimumBFGS() {
return;
}
}
*/
}
NOTIFY("Now starting optimizing with BFGS...\n");
@@ -145,8 +146,11 @@ void LBfgs<OptimizedFunction>::ComputeLocalOptimumBFGS() {
gradient_.ptr(), gradient_.ptr());
num_of_iterations_++;
if (success_bfgs==SUCCESS_FAIL){
NOTIFY("LBFGS failed to find a direction, continuing with gradient descent\n");
ComputeWolfeStep_(&step_, gradient_);
NOTIFY("LBFGS failed to find a direction, continuing with gradient descent\n");
// if (ComputeWolfeStep_(&step_, gradient_)==SUCCESS_FAIL) {
// NONFATAL("Gradient descent failed too");
// }
break;
}
if (silent_==false) {
ReportProgressFile_();
@@ -357,10 +361,10 @@ success_t LBfgs<OptimizedFunction>::ComputeBFGS_(double *step, Matrix &grad, ind
double y_y = la::Dot(num_of_points_*
new_dimension_, y_bfgs_[index_bfgs_].ptr(),
y_bfgs_[index_bfgs_].ptr());
if (unlikely(y_y<1e-10)){
if (unlikely(y_y<1e-40)){
NONFATAL("Gradient differences close to singular...norm=%lg\n", y_y);
}
double norm_scale=s_y/(y_y+1e-10);
double norm_scale=s_y/(y_y+1e-40);
la::Scale(norm_scale, &temp_direction);
Matrix scaled_s;
double beta;