diff --git a/internal/ceres/graph_algorithms.h b/internal/ceres/graph_algorithms.h index 2e6eec0e6..f38a13fdf 100644 --- a/internal/ceres/graph_algorithms.h +++ b/internal/ceres/graph_algorithms.h @@ -45,9 +45,9 @@ namespace internal { // Compare two vertices of a graph by their degrees. template -class VertexDegreeLessThan { +class VertexTotalOrdering { public: - explicit VertexDegreeLessThan(const Graph& graph) + explicit VertexTotalOrdering(const Graph& graph) : graph_(graph) {} bool operator()(const Vertex& lhs, const Vertex& rhs) const { @@ -61,6 +61,20 @@ class VertexDegreeLessThan { const Graph& graph_; }; +template +class VertexDegreeLessThan { + public: + explicit VertexDegreeLessThan(const Graph& graph) + : graph_(graph) {} + + bool operator()(const Vertex& lhs, const Vertex& rhs) const { + return graph_.Neighbors(lhs).size() < graph_.Neighbors(rhs).size(); + } + + private: + const Graph& graph_; +}; + // Order the vertices of a graph using its (approximately) largest // independent set, where an independent set of a graph is a set of // vertices that have no edges connecting them. The maximum @@ -104,7 +118,7 @@ int IndependentSetOrdering(const Graph& graph, sort(vertex_queue.begin(), vertex_queue.end(), - VertexDegreeLessThan(graph)); + VertexTotalOrdering(graph)); // Iterate over vertex_queue. Pick the first white vertex, add it // to the independent set. Mark it black and its neighbors grey. @@ -143,6 +157,81 @@ int IndependentSetOrdering(const Graph& graph, return independent_set_size; } +// Same as above with one important difference. The ordering parameter +// is an input/output parameter which carries an initial ordering of +// the vertices of the graph. The greedy independent set algorithm +// starts by sorting the vertices in increasing order of their +// degree. The input ordering is used to stabilize this sort, i.e., if +// two vertices have the same degree then they are ordered in the same +// order in which they occur in "ordering". +// +// This is useful in eliminating non-determinism from the Schur +// ordering algorithm over all. +template +int StableIndependentSetOrdering(const Graph& graph, + vector* ordering) { + CHECK_NOTNULL(ordering); + const HashSet& vertices = graph.vertices(); + const int num_vertices = vertices.size(); + CHECK_EQ(vertices.size(), ordering->size()); + + // Colors for labeling the graph during the BFS. + const char kWhite = 0; + const char kGrey = 1; + const char kBlack = 2; + + vector vertex_queue(*ordering); + + stable_sort(vertex_queue.begin(), vertex_queue.end(), + VertexDegreeLessThan(graph)); + + // Mark all vertices white. + HashMap vertex_color; + for (typename HashSet::const_iterator it = vertices.begin(); + it != vertices.end(); + ++it) { + vertex_color[*it] = kWhite; + } + + ordering->clear(); + ordering->reserve(num_vertices); + // Iterate over vertex_queue. Pick the first white vertex, add it + // to the independent set. Mark it black and its neighbors grey. + for (int i = 0; i < vertex_queue.size(); ++i) { + const Vertex& vertex = vertex_queue[i]; + if (vertex_color[vertex] != kWhite) { + continue; + } + + ordering->push_back(vertex); + vertex_color[vertex] = kBlack; + const HashSet& neighbors = graph.Neighbors(vertex); + for (typename HashSet::const_iterator it = neighbors.begin(); + it != neighbors.end(); + ++it) { + vertex_color[*it] = kGrey; + } + } + + int independent_set_size = ordering->size(); + + // Iterate over the vertices and add all the grey vertices to the + // ordering. At this stage there should only be black or grey + // vertices in the graph. + for (typename vector::const_iterator it = vertex_queue.begin(); + it != vertex_queue.end(); + ++it) { + const Vertex vertex = *it; + DCHECK(vertex_color[vertex] != kWhite); + if (vertex_color[vertex] != kBlack) { + ordering->push_back(vertex); + } + } + + CHECK_EQ(ordering->size(), num_vertices); + return independent_set_size; +} + // Find the connected component for a vertex implemented using the // find and update operation for disjoint-set. Recursively traverse // the disjoint set structure till you reach a vertex whose connected diff --git a/internal/ceres/graph_algorithms_test.cc b/internal/ceres/graph_algorithms_test.cc index 78a0452cc..7c244766b 100644 --- a/internal/ceres/graph_algorithms_test.cc +++ b/internal/ceres/graph_algorithms_test.cc @@ -165,7 +165,7 @@ TEST(Degree2MaximumSpanningForest, StarGraph) { } } -TEST(VertexDegreeLessThan, TotalOrdering) { +TEST(VertexTotalOrdering, TotalOrdering) { Graph graph; graph.AddVertex(0); graph.AddVertex(1); @@ -178,7 +178,7 @@ TEST(VertexDegreeLessThan, TotalOrdering) { // 0,1 and 2 have degree 1 and 3 has degree 2. graph.AddEdge(0, 1, 1.0); graph.AddEdge(2, 3, 1.0); - VertexDegreeLessThan less_than(graph); + VertexTotalOrdering less_than(graph); for (int i = 0; i < 4; ++i) { EXPECT_FALSE(less_than(i, i)) << "Failing vertex: " << i; @@ -196,5 +196,49 @@ TEST(VertexDegreeLessThan, TotalOrdering) { } } + +TEST(StableIndependentSet, BreakTies) { + Graph graph; + graph.AddVertex(0); + graph.AddVertex(1); + graph.AddVertex(2); + graph.AddVertex(3); + + graph.AddEdge(0, 1); + graph.AddEdge(0, 2); + graph.AddEdge(0, 3); + graph.AddEdge(1, 2); + graph.AddEdge(1, 3); + graph.AddEdge(2, 3); + + // Since this is a completely connected graph, the independent set + // contains exactly one vertex. StableIndependentSetOrdering + // guarantees that it will always be the first vertex in the + // ordering vector. + { + vector ordering; + ordering.push_back(0); + ordering.push_back(1); + ordering.push_back(2); + ordering.push_back(3); + const int independent_set_size = + StableIndependentSetOrdering(graph, &ordering); + EXPECT_EQ(independent_set_size, 1); + EXPECT_EQ(ordering[0], 0); + } + + { + vector ordering; + ordering.push_back(1); + ordering.push_back(0); + ordering.push_back(2); + ordering.push_back(3); + const int independent_set_size = + StableIndependentSetOrdering(graph, &ordering); + EXPECT_EQ(independent_set_size, 1); + EXPECT_EQ(ordering[0], 1); + } + +} } // namespace internal } // namespace ceres diff --git a/internal/ceres/parameter_block_ordering.cc b/internal/ceres/parameter_block_ordering.cc index e8f626f8e..190715bee 100644 --- a/internal/ceres/parameter_block_ordering.cc +++ b/internal/ceres/parameter_block_ordering.cc @@ -42,6 +42,32 @@ namespace ceres { namespace internal { +int ComputeStableSchurOrdering(const Program& program, + vector* ordering) { + CHECK_NOTNULL(ordering)->clear(); + + scoped_ptr > graph(CreateHessianGraph(program)); + const vector& parameter_blocks = program.parameter_blocks(); + const HashSet& vertices = graph->vertices(); + for (int i = 0; i < parameter_blocks.size(); ++i) { + if (vertices.count(parameter_blocks[i]) > 0) { + ordering->push_back(parameter_blocks[i]); + } + } + + int independent_set_size = StableIndependentSetOrdering(*graph, ordering); + + // Add the excluded blocks to back of the ordering vector. + for (int i = 0; i < parameter_blocks.size(); ++i) { + ParameterBlock* parameter_block = parameter_blocks[i]; + if (parameter_block->IsConstant()) { + ordering->push_back(parameter_block); + } + } + + return independent_set_size; +} + int ComputeSchurOrdering(const Program& program, vector* ordering) { CHECK_NOTNULL(ordering)->clear(); diff --git a/internal/ceres/parameter_block_ordering.h b/internal/ceres/parameter_block_ordering.h index a5277a44c..4675cb8dc 100644 --- a/internal/ceres/parameter_block_ordering.h +++ b/internal/ceres/parameter_block_ordering.h @@ -58,6 +58,12 @@ class ParameterBlock; int ComputeSchurOrdering(const Program& program, vector* ordering); +// Same as above, except that ties while computing the independent set +// ordering are resolved in favour of the order in which the parameter +// blocks occur in the program. +int ComputeStableSchurOrdering(const Program& program, + vector* ordering); + // Use an approximate independent set ordering to decompose the // parameter blocks of a problem in a sequence of independent // sets. The ordering covers all the non-constant parameter blocks in diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc index ddcbb3736..56993c874 100644 --- a/internal/ceres/solver_impl.cc +++ b/internal/ceres/solver_impl.cc @@ -1451,8 +1451,8 @@ bool SolverImpl::ReorderProgramForSchurTypeLinearSolver( // this means that the user wishes for Ceres to identify the // e_blocks, which we do by computing a maximal independent set. vector schur_ordering; - const int num_eliminate_blocks = ComputeSchurOrdering(*program, - &schur_ordering); + const int num_eliminate_blocks = ComputeStableSchurOrdering(*program, + &schur_ordering); CHECK_EQ(schur_ordering.size(), program->NumParameterBlocks()) << "Congratulations, you found a Ceres bug! Please report this error "