From 856183f08069a87d0d8f2131cdff2296059b0d41 Mon Sep 17 00:00:00 2001 From: Garry Boyer Date: Wed, 6 Jun 2007 21:30:14 +0000 Subject: [PATCH] hi --- fastlib/u/garryb/nbr/affinity.cc | 440 +++++++++++++++++++----------- fastlib/u/garryb/nbr/affinity.tex | 89 ++++++ fastlib/u/garryb/nbr/cache.h | 8 + fastlib/u/garryb/nbr/dfs.h | 4 +- fastlib/u/garryb/nbr/spbounds.h | 18 ++ 5 files changed, 396 insertions(+), 163 deletions(-) create mode 100644 fastlib/u/garryb/nbr/affinity.tex diff --git a/fastlib/u/garryb/nbr/affinity.cc b/fastlib/u/garryb/nbr/affinity.cc index 3baaec3a8e..d505d75d30 100644 --- a/fastlib/u/garryb/nbr/affinity.cc +++ b/fastlib/u/garryb/nbr/affinity.cc @@ -5,6 +5,32 @@ #include "fastlib/fastlib.h" +/* + +As two-variable functions: + +\rho(i, k) = \sum_{j != i, j != k} max(0, S(j,k) - \alpha(j,k)) + +\alpha(i, k) = min(0, \max_{j != k} S(j,j) + \alpha(j,j) + \rho(i,j)) + +As one-variable rho: + +\rho(k) = \sum_{j != k} max(0, S(j,k) - \alpha(j,k)) + +\sum max(0, S(j,k) - \alpha(j,k)) - max(0, S(k,k) - alpha(k,k)) + +\alpha(i, k) = min(0, \max_{j != k} S(j,j) + \alpha(j,j) + + \rho(i) - max(0, S(i, j) - \alpha(i, j))) + +\alpha(i) = min(0, \max^2{j} + S(j,j) + \alpha(j,j) + \rho(j) - max(0, S(i, j) - \alpha(i, j))) + except when i = j in which case we don't need to do the second part + + S(j,j) + \alpha(j,j) + \rho(j) - S(i, j) + min(S(i,j), \alpha(i, j)) +*/ + +#define SIMILARITY_MAX 0 + struct AffinityCommon { /** The bounding type. Required by NBR. */ typedef SpHrectBound<2> Bound; @@ -50,8 +76,10 @@ struct AffinityCommon { struct Param { public: - /** The epsilon for approximation. */ + /** The epsilon for approximation for rho - ABSOLUTE error. */ double eps; + /** Epsilon per point. */ + double eps_per_point; /** The dimensionality of the data sets. */ index_t dim; /** Number of points */ @@ -63,6 +91,7 @@ struct AffinityCommon { OT_DEF(Param) { OT_MY_OBJECT(eps); + OT_MY_OBJECT(eps_per_point); OT_MY_OBJECT(dim); OT_MY_OBJECT(n_points); OT_MY_OBJECT(pref); @@ -71,16 +100,17 @@ struct AffinityCommon { public: void Copy(const Param& other) { - pref = other.pref; eps = other.eps; + eps_per_point = other.eps_per_point; dim = other.dim; n_points = other.n_points; + pref = other.pref; lambda = other.lambda; } void Init(datanode *module) { dim = -1; - eps = fx_param_double(module, "eps", 1.0e-2); + eps = fx_param_double(module, "eps", 0.0); pref = fx_param_double_req(module, "pref"); lambda = fx_param_double(module, "lambda", 0.8); } @@ -88,11 +118,17 @@ struct AffinityCommon { void BootstrapMonochromatic(CombinedPoint *point, index_t count) { dim = point->vec().length(); n_points = count; - // TODO: Realistic values - point->info().rho = 0; - point->info().alpha.max1 = 0; - point->info().alpha.max2 = 0; - point->info().alpha.max1_index = 0; + // NOTE: These values are manually assigned to be different later. + point->info().rho = DBL_NAN; + point->info().alpha.max1 = DBL_NAN; + point->info().alpha.max2 = DBL_NAN; + point->info().alpha.max1_index = -1; + eps_per_point = eps / n_points; + } + + void SetEpsilon(double eps_in) { + eps = eps_in; + eps_per_point = eps / n_points; } }; @@ -171,10 +207,16 @@ struct AffinityCommon { } return lo; } - - static double ErrorShare(const Param& param, - double abs_error_used, const CombinedNode& r_node) { - return (param.eps - abs_error_used) * r_node.count() / param.n_points; + static double SimilarityHi( + const Param& param, + const Vector& a, index_t a_index, const CombinedNode& b) { + double distsq = b.bound().MinDistanceSqToPoint(a); + double hi = Similarity(distsq); + if (a_index < b.end() && a_index >= b.begin() + && param.pref > hi) { + hi = param.pref; + } + return hi; } }; }; @@ -248,6 +290,7 @@ class AffinityAlpha { } void ApplyPostponed(const Param& param, const QPostponed& postponed, const QNode& q_node) {} + void StartReaccumulate(const Param& param, const QNode& q_node) { alpha.InitEmptySet(); } @@ -261,31 +304,6 @@ class AffinityAlpha { void FinishReaccumulate(const Param& param, const QNode& q_node) {} }; - - /* - - As two-variable functions: - - \rho(i, k) = \sum_{j != i, j != k} max(0, S(j,k) - \alpha(j,k)) - - \alpha(i, k) = min(0, \max_{j != k} S(j,j) + \alpha(j,j) + \rho(i,j)) - - As one-variable rho: - - \rho(k) = \sum_{j != k} max(0, S(j,k) - \alpha(j,k)) - - \sum max(0, S(j,k) - \alpha(j,k)) - max(0, S(k,k) - alpha(k,k)) - - \alpha(i, k) = min(0, \max_{j != k} S(j,j) + \alpha(j,j) - + \rho(i) - max(0, S(i, j) - \alpha(i, j))) - - \alpha(i) = min(0, \max^2{j} - S(j,j) + \alpha(j,j) + \rho(j) - max(0, S(i, j) - \alpha(i, j))) - except when i = j in which case we don't need to do the second part - - S(j,j) + \alpha(j,j) + \rho(j) - S(i, j) + min(S(i,j), \alpha(i, j)) - */ - struct PairVisitor { public: Alpha old_alpha; @@ -295,20 +313,33 @@ class AffinityAlpha { void Init(const Param& param) {} bool StartVisitingQueryPoint(const Param& param, - const QPoint& q, + const QPoint& q, index_t q_index, const RNode& r_node, const QMassResult& unapplied_mass_results, QResult* q_result, GlobalResult* global_result) { - // We could add a pruning rule here to speed things up quite a bit. + double alpha_hi; + alpha = q_result->alpha; old_alpha = q.info().alpha; - return true; + + if (unlikely(q_index >= r_node.begin() && q_index < r_node.end())) { + alpha_hi = r_node.stat().rho.hi + old_alpha.max1; + } else { + double sim = AffinityCommon::Helpers::Similarity( + r_node.bound().MinDistanceSqToPoint(q.vec())); + alpha_hi = min( + min(sim, old_alpha.max1) + r_node.stat().rho.hi, + sim); + } + + return (alpha_hi > alpha.max2); } void VisitPair(const Param& param, const QPoint& q, index_t q_index, const RPoint& r, index_t r_index) { double candidate_alpha; if (likely(q_index != r_index)) { - double sim = AffinityCommon::Helpers::Similarity(q.vec(), r.vec()); + double sim = AffinityCommon::Helpers::Similarity( + la::DistanceSqEuclidean(q.vec(), r.vec())); candidate_alpha = min( min(sim, old_alpha.get(r_index)) + r.info().rho, sim); @@ -326,7 +357,7 @@ class AffinityAlpha { } } void FinishVisitingQueryPoint(const Param& param, - const QPoint& q, + const QPoint& q, index_t q_index, const RNode& r_node, const QMassResult& unapplied_mass_results, QResult* q_result, GlobalResult* global_result) { q_result->alpha = alpha; @@ -339,22 +370,19 @@ class AffinityAlpha { const QNode& q_node, const RNode& r_node, Delta* delta, GlobalResult* global_result, QPostponed* q_postponed) { - //WALDO double sim_lo = AffinityCommon::Helpers::SimilarityLo( param, q_node, r_node); + delta->alpha.lo = min(min(q_node.stat().alpha.lo, sim_lo) + + r_node.stat().rho.lo, sim_lo); - delta->alpha.lo = min( - min(q_node.stat().alpha.lo, sim_lo) + r_node.stat().rho.lo, - sim_lo); if (q_node.begin() < r_node.end() && r_node.begin() < q_node.end()) { delta->alpha.hi = q_node.stat().alpha.hi + r_node.stat().rho.hi; } else { - double sim_hi = AffinityCommon::Helpers::SimilarityHi( - param, q_node, r_node); - delta->alpha.hi = min( - min(q_node.stat().alpha.hi, sim_hi) + r_node.stat().rho.hi, - sim_hi); + double sim_hi = AffinityCommon::Helpers::Similarity( + q_node.bound().MinDistanceSqToBound(r_node.bound())); + delta->alpha.hi = min(min(q_node.stat().alpha.hi, sim_hi) + + r_node.stat().rho.hi, sim_hi); } return true; @@ -363,11 +391,7 @@ class AffinityAlpha { const QNode& q_node, const RNode& r_node, const Delta& delta, const QMassResult& q_mass_result, const GlobalResult& global_result, QPostponed* q_postponed) { - if (delta.alpha.hi <= q_mass_result.alpha.lo) { - return false; - } else { - return true; - } + return (delta.alpha.hi > q_mass_result.alpha.lo); } static bool ConsiderQueryTermination(const Param& param, const QNode& q_node, @@ -396,14 +420,15 @@ class AffinityRho { typedef BlankGlobalResult GlobalResult; +#ifdef APPROX struct QPostponed { public: double d_rho; - double abs_error_used; + double error_buffer; OT_DEF(QPostponed) { OT_MY_OBJECT(d_rho); - OT_MY_OBJECT(abs_error_used); + OT_MY_OBJECT(error_buffer); } public: @@ -413,14 +438,17 @@ class AffinityRho { void Reset(const Param& param) { d_rho = 0; - abs_error_used = 0; + error_buffer = 0; } void ApplyPostponed(const Param& param, const QPostponed& other) { d_rho += other.d_rho; - abs_error_used += other.abs_error_used; + error_buffer += other.error_buffer; } }; +#else + typedef BlankQPostponed QPostponed; +#endif struct Delta { public: @@ -438,53 +466,57 @@ class AffinityRho { struct QResult { public: double rho; - double abs_error_used; + #ifdef APPROX + double error_buffer; + #endif OT_DEF(QResult) { OT_MY_OBJECT(rho); - OT_MY_OBJECT(abs_error_used); + #ifdef APPROX + OT_MY_OBJECT(error_buffer); + #endif } public: void Init(const Param& param) { rho = 0; - abs_error_used = 0; + #ifdef APPROX + error_buffer = 0; + #endif } void Postprocess(const Param& param, const QPoint& q, index_t q_index, const RNode& r_root) { - double self_responsibility = - param.pref - q.info().alpha.get(q_index); - - // Make sure we count ourselves regardless of sign. - if (self_responsibility < 0) { - rho += self_responsibility; - } + double self_responsibility = param.pref - q.info().alpha.get(q_index); + rho += self_responsibility; } void ApplyPostponed(const Param& param, const QPostponed& postponed, const QPoint& q) { + #ifdef APPRLX rho += postponed.d_rho; - abs_error_used += postponed.abs_error_used; + error_buffer += postponed.error_buffer; + #endif } }; +#ifdef APPROX struct QMassResult { public: SpRange rho; - double abs_error_used; + double error_buffer; OT_DEF(QMassResult) { OT_MY_OBJECT(rho); - OT_MY_OBJECT(abs_error_used); + OT_MY_OBJECT(error_buffer); } public: void Init(const Param& param) { rho.Init(0, 0); - abs_error_used = 0; + error_buffer = 0; } void ApplyMassResult(const Param& param, const QMassResult& mass_result) { rho += mass_result.rho; - abs_error_used += mass_result.abs_error_used; + error_buffer += mass_result.error_buffer; } void ApplyDelta(const Param& param, const Delta& delta) { rho += delta.d_rho; @@ -492,23 +524,26 @@ class AffinityRho { void ApplyPostponed(const Param& param, const QPostponed& postponed, const QNode& q_node) { rho += postponed.d_rho; - abs_error_used += postponed.abs_error_used; + error_buffer += postponed.error_buffer; } void StartReaccumulate(const Param& param, const QNode& q_node) { rho.InitEmptySet(); - abs_error_used = 0; + error_buffer = DBL_MAX; } void Accumulate(const Param& param, const QResult& result) { rho |= result.rho; - abs_error_used = max(abs_error_used, result.abs_error_used); + error_buffer = min(error_buffer, result.error_buffer); } void Accumulate(const Param& param, const QMassResult& result, index_t n_points) { rho |= result.rho; - abs_error_used = max(abs_error_used, result.abs_error_used); + error_buffer = min(error_buffer, result.error_buffer); } void FinishReaccumulate(const Param& param, const QNode& q_node) {} }; +#else + typedef BlankQMassResult QMassResult; +#endif struct PairVisitor { public: @@ -517,26 +552,36 @@ class AffinityRho { public: void Init(const Param& param) {} bool StartVisitingQueryPoint(const Param& param, - const QPoint& q, + const QPoint& q, index_t q_index, const RNode& r_node, const QMassResult& unapplied_mass_results, QResult* q_result, GlobalResult* global_result) { - rho = q_result->rho; - return true; + // do the point-node prune check + double sim_hi = AffinityCommon::Helpers::SimilarityHi( + param, q.vec(), q_index, r_node); + #ifdef APPROX + q_result->error_buffer += param.eps_per_point * r_node.count(); + #endif + if (sim_hi < r_node.stat().alpha.lo) { + return false; + } else { + rho = q_result->rho; + return true; + } } void VisitPair(const Param& param, const QPoint& q, index_t q_index, const RPoint& r, index_t r_index) { - double responsibility = - AffinityCommon::Helpers::Similarity( - param, q.vec(), q_index, r.vec(), r_index) - - r.info().alpha.get(q_index); - - if (responsibility > 0) { - rho += responsibility; + if (likely(q_index != r_index)) { + double responsibility = + AffinityCommon::Helpers::Similarity(q.vec(), r.vec()) + - r.info().alpha.get(q_index); + if (responsibility > 0) { + rho += responsibility; + } } } void FinishVisitingQueryPoint(const Param& param, - const QPoint& q, + const QPoint& q, index_t q_index, const RNode& r_node, const QMassResult& unapplied_mass_results, QResult* q_result, GlobalResult* global_result) { q_result->rho = rho; @@ -554,11 +599,6 @@ class AffinityRho { double sim_lo = AffinityCommon::Helpers::SimilarityLo( param, q_node, r_node); - // fprintf(stderr, "(%d,%d) alpha.lo,hi = (%f, %f)\n", - // r_node.begin(), r_node.end(), - // r_node.stat().alpha.hi, - // r_node.stat().alpha.lo - // ); delta->d_rho.lo = max(0.0, sim_lo - r_node.stat().alpha.hi) * r_node.count(); delta->d_rho.hi = max(0.0, sim_hi - r_node.stat().alpha.lo) @@ -570,21 +610,24 @@ class AffinityRho { const QNode& q_node, const RNode& r_node, const Delta& delta, const QMassResult& q_mass_result, const GlobalResult& global_result, QPostponed* q_postponed) { - /* - double abs_error = delta.d_rho.width() / 2; - double rel_error_hi = abs_error / q_mass_result.rho.lo; - - if (rel_error_hi < AffinityCommon::Helpers::ErrorShare( - param, q_mass_result.abs_error_used, r_node)) { - q_postponed->abs_error_used += abs_error; + #ifdef APPROX + double my_abs_error = delta.d_rho.width() * 0.5; + double allotted_error = param.eps_per_point * r_node.count(); + double allowed_error = allotted_error + q_mass_result.error_buffer; + + //if (param.eps != 0) + //fprintf(stderr, "%g %g %g %g\n", my_abs_error, param.eps_per_point, allotted_error, q_mass_result.error_buffer); + + if (my_abs_error <= allowed_error) { + q_postponed->error_buffer += allotted_error - my_abs_error; q_postponed->d_rho += delta.d_rho.mid(); return false; } else { return true; } - */ - + #else return true; + #endif } static bool ConsiderQueryTermination(const Param& param, const QNode& q_node, @@ -594,34 +637,87 @@ class AffinityRho { } static double Heuristic(const Param& param, const QNode& q_node, const RNode& r_node, const Delta& delta) { - // favor whatever brings our lower bound up the fastest - return -delta.d_rho.lo; + // TODO: If approximating, favor upper bound + return -delta.d_rho.hi; } }; }; -void FindExemplars(index_t dimensionality, index_t n_points, +struct Cluster { + Vector exemplar; + Vector centroid; + index_t count; +}; + +void FindExemplars(datanode *module, const AffinityCommon::Param& param, + index_t dimensionality, index_t n_points, CacheArray *data_points) { - ArrayList exemplars; + fx_timer_start(module, "exemplars"); + + ArrayList clusters; + ArrayList assignments; CacheReadIterator point(data_points, 0); - exemplars.Init(); + clusters.Init(); for (index_t point_i = 0; point_i < n_points; point_i++, point.Next()) { if (point->info().rho > 0) { - exemplars.AddBack()->Copy(point->vec()); + Cluster *cluster = clusters.AddBack(); + cluster->exemplar.Copy(point->vec()); + cluster->centroid.Init(dimensionality); + cluster->centroid.SetZero(); + cluster->count = 0; } } - ot::Print(exemplars); + // Run all nearest neighbors to assign clusters and find centroids + assignments.Init(n_points); + point.SetIndex(0); + for (index_t point_i = 0; point_i < n_points; point_i++, point.Next()) { + double best_distsq = DBL_MAX; + index_t best_k = -1; + + for (index_t k = 0; k < clusters.size(); k++) { + double distsq = la::DistanceSqEuclidean( + clusters[k].exemplar, point->vec()); + if (unlikely(distsq < best_distsq)) { + best_distsq = distsq; + best_k = k; + } + } + + assignments[point_i] = best_k; + clusters[best_k].count++; + la::AddTo(point->vec(), &clusters[best_k].centroid); + } + + // Divide centroids by size + for (index_t i = 0; i < clusters.size(); i++) { + la::Scale(1.0 / clusters[i].count, &clusters[i].centroid); + } + + // Calculate net similarity + double netsim = param.pref * clusters.size(); + + point.SetIndex(0); + for (index_t point_i = 0; point_i < n_points; point_i++, point.Next()) { + Cluster *cluster = &clusters[assignments[point_i]]; + double distsq = la::DistanceSqEuclidean(cluster->exemplar, point->vec()); + netsim -= distsq; + } + + fx_timer_start(module, "exemplars"); + fx_format_result(module, "netsim", "%f", netsim); + + // Make a matrix of exemplars so we can save it to file Matrix m; - m.Init(dimensionality, exemplars.size()); + m.Init(dimensionality, clusters.size()); - for (index_t i = 0; i < exemplars.size(); i++) { + for (index_t i = 0; i < clusters.size(); i++) { Vector dest; m.MakeColumnVector(i, &dest); - dest.CopyValues(exemplars[i]); + dest.CopyValues(clusters[i].exemplar); } data::Save("exemplars.txt", m); @@ -662,12 +758,10 @@ void FindCovariance(const Matrix& dataset) { Matrix u; // eigenvectors Matrix ui; // the inverse of eigenvectors - //cov.PrintDebug("cov"); la::EigenvectorsInit(cov, &d, &u); d.PrintDebug("covariance_eigenvectors"); la::TransposeInit(u, &ui); -// // for (index_t i = 0; i < d.length(); i++) { // d[i] = 1.0 / sqrt(d[i]); // } @@ -717,7 +811,7 @@ void AffinityMain(datanode *module, const char *gnp_name) { AffinityAlpha::QPoint default_point; default_point.vec().Init(data_matrix.n_rows()); param.BootstrapMonochromatic(&default_point, data_matrix.n_cols()); - data_points.Init(default_point, data_matrix.n_cols(), n_block_points); + data_points.Init(default_point, n_points, n_block_points); for (index_t i = 0; i < data_matrix.n_cols(); i++) { CacheWrite point(&data_points, i); point->vec().CopyValues(data_matrix.GetColumnPtr(i)); @@ -730,6 +824,19 @@ void AffinityMain(datanode *module, const char *gnp_name) { AffinityAlpha::QNode data_example_node; data_example_node.Init(dimensionality, param); + + // Initial conditions for alpha + for (index_t i = 0; i < n_points; i++) { + CacheWrite point(&data_points, i); + point->info().alpha.max1 = 0; + point->info().alpha.max2 = param.pref; + point->info().alpha.max1_index = i; + if (rand() % 16 == 0) { + point->info().rho = -param.pref / 2; + } else { + point->info().rho = 0; + } + } data_nodes.Init(data_example_node, 0, n_block_nodes); KdTreeMidpointBuilder @@ -742,37 +849,34 @@ void AffinityMain(datanode *module, const char *gnp_name) { // All the above is not any different for affinity than for anything // else. Now, time for the iteration. int n_iter = fx_param_int(module, "n_iter", 10000); - int stable_iterations = 0; + int iter; + int n_convergence_iter = fx_param_int(module, "n_convergence_iter", 8);; + int stable_iter = 0; index_t n_changed = n_points / 2; + index_t n_exemplars = 0; + timer *timer_rho = fx_timer(module, "all_rho"); + timer *timer_alpha = fx_timer(module, "all_alpha"); ArrayList is_exemplar; - ArrayList alpha_d; - ArrayList alpha_dd; is_exemplar.Init(n_points); - alpha_d.Init(n_points); - alpha_dd.Init(n_points); - fprintf(stderr, " --- killing alpha ---\n"); for (index_t i = 0; i < n_points; i++) { - CacheWrite point(&data_points, i); - point->info().alpha.max1 = 0; - point->info().alpha.max2 = param.pref; - point->info().alpha.max1_index = i; - alpha_d[i] = 0; - alpha_dd[i] = 0; - point->info().rho = 0; is_exemplar[i] = 0; } - for (int iter = 0; iter < n_iter && stable_iterations < 50; iter++) + for (iter = 0; iter < n_iter && stable_iter < n_convergence_iter; iter++) { double lambda = param.lambda; index_t n_alpha_changed = 0; - index_t n_exemplars = 0; double sum_alpha = 0; double sum_alpha2 = 0; double sum_rho = 0; + index_t unclassifieds; + n_exemplars = 0; + unclassifieds = 0; + + fx_timer_start(module, "all_alpha"); { TempCacheArray q_results_alpha; @@ -783,7 +887,7 @@ void AffinityMain(datanode *module, const char *gnp_name) { nbr_utils::ThreadedDualTreeSolver < AffinityAlpha, DualTreeDepthFirst >::Solve( - fx_submodule(module, "threads", "iter%d_alpha", iter), param, + fx_submodule(module, "threads", "iters/%d/alpha", iter), param, &data_points, &data_nodes, &data_points, &data_nodes, &q_results_alpha); @@ -795,26 +899,30 @@ void AffinityMain(datanode *module, const char *gnp_name) { n_alpha_changed++; } -/* - double max1 = damp(lambda, alpha_d[i], result->alpha.max1); - //double max1 = damp(lambda, point->info().alpha.max1, result->alpha.max1); - max1 = max(result->alpha.max2, max1); - alpha_d[i] = max1; - - point->info().alpha.max1 = max1; -*/ point->info().alpha = result->alpha; sum_alpha += result->alpha.max1; sum_alpha2 += result->alpha.max2; + + if (!is_exemplar[result->alpha.max1_index] && point->info().rho <= 0) { + unclassifieds++; + } } - + nbr_utils::StatFixer< AffinityAlpha::Param, AffinityAlpha::QPoint, AffinityAlpha::QNode> ::Fix(param, &data_points, &data_nodes); } + fx_timer_stop(module, "all_alpha"); + fprintf(stderr, + "\033[31m -------- %04d alpha: sum_alpha = %f, sum_alpha2 = %f, %"LI"d alphas changed, %"LI"d unclassifieds, %.3fsec cum. alpha\033[0m\n", + iter, sum_alpha, sum_alpha2, n_alpha_changed, unclassifieds, + timer_alpha->total.micros / 1.0e6); + + fx_timer_start(module, "all_rho"); { + double rho_squared_difference = 0; TempCacheArray q_results_rho; AffinityRho::QResult default_result_rho; @@ -823,7 +931,7 @@ void AffinityMain(datanode *module, const char *gnp_name) { data_points.n_block_elems()); nbr_utils::ThreadedDualTreeSolver< AffinityRho, DualTreeDepthFirst >::Solve( - fx_submodule(module, "threads", "iter%d_rho", iter), param, + fx_submodule(module, "threads", "iters/%d/rho", iter), param, &data_points, &data_nodes, &data_points, &data_nodes, &q_results_rho); @@ -835,52 +943,62 @@ void AffinityMain(datanode *module, const char *gnp_name) { double old_rho = point->info().rho; double new_rho = damp(lambda, old_rho, result->rho); - point->info().rho = new_rho; - sum_rho += new_rho; - if ((old_rho > 0) != (new_rho > 0)) { - point->info().rho *= math::Random(0.2, 1.8); + new_rho *= math::Random(0.4, 1.6); n_changed++; } + rho_squared_difference += math::Sqr(new_rho - old_rho); + sum_rho += new_rho; + if (new_rho > 0) { is_exemplar[i] = 1; n_exemplars++; } else { is_exemplar[i] = 0; } + + point->info().rho = new_rho; } if (n_changed == 0) { - stable_iterations++; + stable_iter++; } else { - stable_iterations = 0; + stable_iter = 0; } nbr_utils::StatFixer< AffinityAlpha::Param, AffinityAlpha::QPoint, AffinityAlpha::QNode> ::Fix(param, &data_points, &data_nodes); - } -/* - for (index_t i = 0; i < n_points; i++) { - CacheWrite point(&data_points, i); - double max1 = damp(lambda, alpha_dd[i], alpha_d[i]); - max1 = max(point->info().alpha.max2, max1); - alpha_dd[i] = max1; - point->info().alpha.max1 = max1; + param.SetEpsilon(sqrt(rho_squared_difference / n_points)); + //param.SetEpsilon(0); } -*/ - fprintf(stderr, "------------- iter %04d: %"LI"d rhos changed, (%"LI"d exemplars), %d alphas, sum_alpha = %f, sum_alpha2 = %f, sum_rho = %f\n", - iter, n_changed, n_exemplars, n_alpha_changed, sum_alpha, sum_alpha2, sum_rho); + fx_timer_stop(module, "all_rho"); + + fprintf(stderr, "\033[32m------------- %04d rho: %"LI"d rhos changed, (%"LI"d exemplars), sum_rho = %f, eps = %f, %.3fsec cum. rho\033[0m\n", + iter, n_changed, n_exemplars, + sum_rho, param.eps, + timer_rho->total.micros / 1.0e6); } - FindExemplars(dimensionality, n_points, &data_points); + fx_format_result(module, "n_iterations", "%d", iter); + fx_format_result(module, "n_exemplars", "%d", n_exemplars); + + // This will take too long if there are too many exemplars. + if (n_exemplars >= 10000) { + NONFATAL("Too many exemplars (%"LI"d >= 10000), NOT performing clustering.\n", + n_exemplars); + } else { + FindExemplars(module, param, dimensionality, n_points, &data_points); + } } int main(int argc, char *argv[]) { fx_init(argc, argv); + srand(time(NULL)); + AffinityMain(fx_root, "affinity"); fx_done(); diff --git a/fastlib/u/garryb/nbr/affinity.tex b/fastlib/u/garryb/nbr/affinity.tex new file mode 100644 index 0000000000..dc72180ec4 --- /dev/null +++ b/fastlib/u/garryb/nbr/affinity.tex @@ -0,0 +1,89 @@ +\documentclass[times, 10pt]{article} +\usepackage{times} +\usepackage{amsmath} + +\DeclareMathOperator*{\map}{map} + +\newcommand{\eqspace}{\!\!\!\!} +\newcommand{\cposij}{c^{+}_{ij}} + +\newcommand{\intersect}{\cap} + +\newcommand{\respo}[2]{R_{#1#2}} +\newcommand{\avail}[2]{A_{#1#2}} +\newcommand{\simil}[2]{S_{#1#2}} + +\newcommand{\vecrho}{\vec{\rho}} +\newcommand{\vecalpha}{\vec{\alpha}} +\newcommand{\frho}[1]{\rho_{#1}} +\newcommand{\falpha}[2]{\alpha_{#1#2}} + +\begin{document} + +Affinity propagation\cite{affinity} is a recent clustering technique that minimizes a metric the authors deem {\it net similarity}. +The similarity between points $x_i$ and $x_k$ is labeled $\simil{i}{j}$, which is often the negative squared Euclidean distance $-||x_i - x_j||^2$, with the special case that $\simil{i}{i}$ is a parameter $p$. +The goal is to find a set of exemplars that maximizes the sum of similarities from each exemplar to all its cluster members, including itself; the number of clusters is determined by the penalty $p$ imposes on cluster creation. +The computation involves iterative refinement of the responsibility message $\respo{i}{j}$ and availability message $\avail{i}{j}$, whose update steps are defined: +\begin{eqnarray*} + \respo{i}{j} \eqspace&\gets&\eqspace \simil{i}{j} - \max_{j' \neq j} (\avail{i}{j'} + \simil{i}{j'}) + \\ + \avail{i}{j} \eqspace&\gets&\eqspace \min \! \left(\! 0, \respo{j}{j} + \sum_{i' \neq i,j} max(0, \respo{i'}{j})\!\right) + \\ + \avail{j}{j} \eqspace&\gets&\eqspace \sum_{i' \neq j} \max(0, \respo{i'}{k}) +\end{eqnarray*} + +% Ryan, I am defining \alpha as the negative max, because it makes things look much prettier. + +\noindent Trivially, $\respo{}{}$ and $\avail{}{}$ can be stored as matrices. +However, it is not necessary to store the entirety of $\respo{}{}$ and $\avail{}{}$; instead, only a sum and a maximum is needed. +We define single-argument $\frho{j}$ and $\falpha{i}{}$ as the aggregate operations in use\footnote{ +Note $\falpha{i}{j}$ requires argument $j$ only for exclusion from the $\max$; computationally, one finds the first and second maximum, and chooses the first or second maximum depending on $j$.}: +\begin{eqnarray*} + \falpha{i}{j} \eqspace&\gets&\eqspace - \max_{j' \neq i,j} (\avail{i}{j'} + \simil{i}{j'}) + \\ + \frho{j} \eqspace&\gets&\eqspace \respo{j}{j} + \sum_{i' \neq j} \left( \max(0, \respo{i'}{j}) \right) +\end{eqnarray*} + +\noindent from which we can easily compute the original expressions: +\begin{eqnarray*} + \respo{i}{j} \eqspace&\gets&\eqspace \simil{i}{j} + \falpha{i}{j} + \\ + \avail{i}{j} \eqspace&\gets&\eqspace \min \left(0, \frho{j} - \max(0, \respo{i}{j}) \right) + \\ + \avail{j}{j} \eqspace&\gets&\eqspace \frho{j} - \respo{j}{j} +\end{eqnarray*} + +\noindent +After a bit of algebraic substitution, we can define $\falpha{i}{}$ and $\frho{j}$ in terms of each other: + +\begin{eqnarray*} + \vecrho \eqspace&\gets&\eqspace \map_{i} \sum_{j} \left( \cposij(\simil{i}{j} + \falpha{i}{j}) \right) + \\ + \vecalpha \eqspace&\gets&\eqspace \map_{j} \min^2_{i} \left( \cposij(\cposij(\simil{i}{j} + \falpha{i}{j}) - \frho{j}) - \simil{i}{j} \right) +\end{eqnarray*} + +\noindent where $\cposij(x)$ is $x$ when $i=j$, and $\max(x,0)$ otherwise. +We now have a problem nearly\footnote{To aid convergence, affinity propagation adds noise to the similarity matrix and dampens both $\respo{}{}$ and $\avail{}{}$ between time-steps. Since we store no matrices, we both dampen and add noise to $\vecrho$.} indentical to traditional affinity propagation. +However, we can show that this is a generalized $N$-body problem. + +{\bf Proposition}. +The computations of both $\vecrho$ and $\vecalpha$ are generalized $N$-body problems. + +{\bf Proof}. +Ryan, prove this please. + +We may now develop a dual-tree algorithm. +First, we build a tree on the data points, assuming they are in a metric space\footnote{Euclidean distance over real vectors is the example we implemented.}. +Notationally, we denote $Q'$ a subset or node of points that are being $\map$ed over, and $R'$ as a subset or node of points for the inner operator. + +For $\vecrho$, there is a simple intrinsic prune. +For subsets $Q'$ and $R'$, can show that $R'$ contributes nothing to the summation for all of $Q'$ and prune accordingly in the following case: + +$$\min{i \in Q', j \in R'} \left( \simil{i}{j} + \alpha{i}{j} \right) < 0$$ + +For $\vecalpha$, we use an extrinsic prune. +For each node of queries $Q'$, we keep track of the largest candidate second-minimum found. +Next, we provide a lower bound candidate minimum between $Q'$ and $R'$. + + +\end{document} diff --git a/fastlib/u/garryb/nbr/cache.h b/fastlib/u/garryb/nbr/cache.h index 49193575f6..c6babd501a 100644 --- a/fastlib/u/garryb/nbr/cache.h +++ b/fastlib/u/garryb/nbr/cache.h @@ -567,6 +567,14 @@ class CacheReadIterator { const Element & operator * () const { return *element_; } + + void SetIndex(index_t begin_index) { + DEBUG_ONLY(cache_->ReleaseBlock_(blockid_)); + blockid_ = begin_index >> cache_->n_block_elems_log_; + element_ = cache_->StartRead(begin_index); + unsigned int mask = cache_->n_block_elems_mask_; + left_ = (begin_index & mask) ^ mask; + } void Next() { element_ = mem::PointerAdd(element_, stride_); diff --git a/fastlib/u/garryb/nbr/dfs.h b/fastlib/u/garryb/nbr/dfs.h index 1c7b56f39e..19d1f9a553 100644 --- a/fastlib/u/garryb/nbr/dfs.h +++ b/fastlib/u/garryb/nbr/dfs.h @@ -356,7 +356,7 @@ void DualTreeDepthFirst::BaseCase_( q_result->ApplyPostponed(param_, q_node_mut->postponed, *q_point); - if (visitor.StartVisitingQueryPoint(param_, *q_point, *r_node, + if (visitor.StartVisitingQueryPoint(param_, *q_point, q_i, *r_node, exclusive_unvisited, q_result, &global_result_)) { CacheReadIterator r_iter(&r_points_, r_node->begin()); @@ -367,7 +367,7 @@ void DualTreeDepthFirst::BaseCase_( visitor.VisitPair(param_, *q_point, q_i, *r_point, r_i); } - visitor.FinishVisitingQueryPoint(param_, *q_point, *r_node, + visitor.FinishVisitingQueryPoint(param_, *q_point, q_i, *r_node, exclusive_unvisited, q_result, &global_result_); DEBUG_ONLY(n_naive_ += r_node->count()); diff --git a/fastlib/u/garryb/nbr/spbounds.h b/fastlib/u/garryb/nbr/spbounds.h index 3c55967380..880953602c 100644 --- a/fastlib/u/garryb/nbr/spbounds.h +++ b/fastlib/u/garryb/nbr/spbounds.h @@ -285,6 +285,24 @@ struct SpRange { hi = range.hi; } } + + void MaxWith(double v) { + if (unlikely(v > lo)) { + lo = v; + if (unlikely(v > hi)) { + hi = v; + } + } + } + + void MinWith(double v) { + if (unlikely(v < hi)) { + hi = v; + if (unlikely(v < lo)) { + lo = v; + } + } + } friend bool operator < (const SpRange& a, const SpRange& b) { return a.hi < b.lo;