1 // Ceres Solver - A fast non-linear least squares minimizer
2 // Copyright 2015 Google Inc. All rights reserved.
3 // http://ceres-solver.org/
5 // Redistribution and use in source and binary forms, with or without
6 // modification, are permitted provided that the following conditions are met:
8 // * Redistributions of source code must retain the above copyright notice,
9 // this list of conditions and the following disclaimer.
10 // * Redistributions in binary form must reproduce the above copyright notice,
11 // this list of conditions and the following disclaimer in the documentation
12 // and/or other materials provided with the distribution.
13 // * Neither the name of Google Inc. nor the names of its contributors may be
14 // used to endorse or promote products derived from this software without
15 // specific prior written permission.
17 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
21 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 // POSSIBILITY OF SUCH DAMAGE.
29 // Author: sameeragarwal@google.com (Sameer Agarwal)
31 // Various algorithms that operate on undirected graphs.
33 #ifndef CERES_INTERNAL_GRAPH_ALGORITHMS_H_
34 #define CERES_INTERNAL_GRAPH_ALGORITHMS_H_
39 #include "ceres/collections_port.h"
40 #include "ceres/graph.h"
41 #include "ceres/wall_time.h"
42 #include "glog/logging.h"
47 // Compare two vertices of a graph by their degrees, if the degrees
48 // are equal then order them by their ids.
49 template <typename Vertex>
50 class VertexTotalOrdering {
52 explicit VertexTotalOrdering(const Graph<Vertex>& graph)
55 bool operator()(const Vertex& lhs, const Vertex& rhs) const {
56 if (graph_.Neighbors(lhs).size() == graph_.Neighbors(rhs).size()) {
59 return graph_.Neighbors(lhs).size() < graph_.Neighbors(rhs).size();
63 const Graph<Vertex>& graph_;
66 template <typename Vertex>
67 class VertexDegreeLessThan {
69 explicit VertexDegreeLessThan(const Graph<Vertex>& graph)
72 bool operator()(const Vertex& lhs, const Vertex& rhs) const {
73 return graph_.Neighbors(lhs).size() < graph_.Neighbors(rhs).size();
77 const Graph<Vertex>& graph_;
80 // Order the vertices of a graph using its (approximately) largest
81 // independent set, where an independent set of a graph is a set of
82 // vertices that have no edges connecting them. The maximum
83 // independent set problem is NP-Hard, but there are effective
84 // approximation algorithms available. The implementation here uses a
85 // breadth first search that explores the vertices in order of
86 // increasing degree. The same idea is used by Saad & Li in "MIQR: A
87 // multilevel incomplete QR preconditioner for large sparse
88 // least-squares problems", SIMAX, 2007.
90 // Given a undirected graph G(V,E), the algorithm is a greedy BFS
91 // search where the vertices are explored in increasing order of their
92 // degree. The output vector ordering contains elements of S in
93 // increasing order of their degree, followed by elements of V - S in
94 // increasing order of degree. The return value of the function is the
96 template <typename Vertex>
97 int IndependentSetOrdering(const Graph<Vertex>& graph,
98 std::vector<Vertex>* ordering) {
99 const HashSet<Vertex>& vertices = graph.vertices();
100 const int num_vertices = vertices.size();
102 CHECK_NOTNULL(ordering);
104 ordering->reserve(num_vertices);
106 // Colors for labeling the graph during the BFS.
107 const char kWhite = 0;
108 const char kGrey = 1;
109 const char kBlack = 2;
111 // Mark all vertices white.
112 HashMap<Vertex, char> vertex_color;
113 std::vector<Vertex> vertex_queue;
114 for (typename HashSet<Vertex>::const_iterator it = vertices.begin();
115 it != vertices.end();
117 vertex_color[*it] = kWhite;
118 vertex_queue.push_back(*it);
122 std::sort(vertex_queue.begin(), vertex_queue.end(),
123 VertexTotalOrdering<Vertex>(graph));
125 // Iterate over vertex_queue. Pick the first white vertex, add it
126 // to the independent set. Mark it black and its neighbors grey.
127 for (int i = 0; i < vertex_queue.size(); ++i) {
128 const Vertex& vertex = vertex_queue[i];
129 if (vertex_color[vertex] != kWhite) {
133 ordering->push_back(vertex);
134 vertex_color[vertex] = kBlack;
135 const HashSet<Vertex>& neighbors = graph.Neighbors(vertex);
136 for (typename HashSet<Vertex>::const_iterator it = neighbors.begin();
137 it != neighbors.end();
139 vertex_color[*it] = kGrey;
143 int independent_set_size = ordering->size();
145 // Iterate over the vertices and add all the grey vertices to the
146 // ordering. At this stage there should only be black or grey
147 // vertices in the graph.
148 for (typename std::vector<Vertex>::const_iterator it = vertex_queue.begin();
149 it != vertex_queue.end();
151 const Vertex vertex = *it;
152 DCHECK(vertex_color[vertex] != kWhite);
153 if (vertex_color[vertex] != kBlack) {
154 ordering->push_back(vertex);
158 CHECK_EQ(ordering->size(), num_vertices);
159 return independent_set_size;
162 // Same as above with one important difference. The ordering parameter
163 // is an input/output parameter which carries an initial ordering of
164 // the vertices of the graph. The greedy independent set algorithm
165 // starts by sorting the vertices in increasing order of their
166 // degree. The input ordering is used to stabilize this sort, i.e., if
167 // two vertices have the same degree then they are ordered in the same
168 // order in which they occur in "ordering".
170 // This is useful in eliminating non-determinism from the Schur
171 // ordering algorithm over all.
172 template <typename Vertex>
173 int StableIndependentSetOrdering(const Graph<Vertex>& graph,
174 std::vector<Vertex>* ordering) {
175 CHECK_NOTNULL(ordering);
176 const HashSet<Vertex>& vertices = graph.vertices();
177 const int num_vertices = vertices.size();
178 CHECK_EQ(vertices.size(), ordering->size());
180 // Colors for labeling the graph during the BFS.
181 const char kWhite = 0;
182 const char kGrey = 1;
183 const char kBlack = 2;
185 std::vector<Vertex> vertex_queue(*ordering);
187 std::stable_sort(vertex_queue.begin(), vertex_queue.end(),
188 VertexDegreeLessThan<Vertex>(graph));
190 // Mark all vertices white.
191 HashMap<Vertex, char> vertex_color;
192 for (typename HashSet<Vertex>::const_iterator it = vertices.begin();
193 it != vertices.end();
195 vertex_color[*it] = kWhite;
199 ordering->reserve(num_vertices);
200 // Iterate over vertex_queue. Pick the first white vertex, add it
201 // to the independent set. Mark it black and its neighbors grey.
202 for (int i = 0; i < vertex_queue.size(); ++i) {
203 const Vertex& vertex = vertex_queue[i];
204 if (vertex_color[vertex] != kWhite) {
208 ordering->push_back(vertex);
209 vertex_color[vertex] = kBlack;
210 const HashSet<Vertex>& neighbors = graph.Neighbors(vertex);
211 for (typename HashSet<Vertex>::const_iterator it = neighbors.begin();
212 it != neighbors.end();
214 vertex_color[*it] = kGrey;
218 int independent_set_size = ordering->size();
220 // Iterate over the vertices and add all the grey vertices to the
221 // ordering. At this stage there should only be black or grey
222 // vertices in the graph.
223 for (typename std::vector<Vertex>::const_iterator it = vertex_queue.begin();
224 it != vertex_queue.end();
226 const Vertex vertex = *it;
227 DCHECK(vertex_color[vertex] != kWhite);
228 if (vertex_color[vertex] != kBlack) {
229 ordering->push_back(vertex);
233 CHECK_EQ(ordering->size(), num_vertices);
234 return independent_set_size;
237 // Find the connected component for a vertex implemented using the
238 // find and update operation for disjoint-set. Recursively traverse
239 // the disjoint set structure till you reach a vertex whose connected
240 // component has the same id as the vertex itself. Along the way
241 // update the connected components of all the vertices. This updating
242 // is what gives this data structure its efficiency.
243 template <typename Vertex>
244 Vertex FindConnectedComponent(const Vertex& vertex,
245 HashMap<Vertex, Vertex>* union_find) {
246 typename HashMap<Vertex, Vertex>::iterator it = union_find->find(vertex);
247 DCHECK(it != union_find->end());
248 if (it->second != vertex) {
249 it->second = FindConnectedComponent(it->second, union_find);
255 // Compute a degree two constrained Maximum Spanning Tree/forest of
256 // the input graph. Caller owns the result.
258 // Finding degree 2 spanning tree of a graph is not always
259 // possible. For example a star graph, i.e. a graph with n-nodes
260 // where one node is connected to the other n-1 nodes does not have
261 // a any spanning trees of degree less than n-1.Even if such a tree
262 // exists, finding such a tree is NP-Hard.
264 // We get around both of these problems by using a greedy, degree
265 // constrained variant of Kruskal's algorithm. We start with a graph
266 // G_T with the same vertex set V as the input graph G(V,E) but an
267 // empty edge set. We then iterate over the edges of G in decreasing
268 // order of weight, adding them to G_T if doing so does not create a
269 // cycle in G_T} and the degree of all the vertices in G_T remains
270 // bounded by two. This O(|E|) algorithm results in a degree-2
271 // spanning forest, or a collection of linear paths that span the
273 template <typename Vertex>
274 WeightedGraph<Vertex>*
275 Degree2MaximumSpanningForest(const WeightedGraph<Vertex>& graph) {
276 // Array of edges sorted in decreasing order of their weights.
277 std::vector<std::pair<double, std::pair<Vertex, Vertex> > > weighted_edges;
278 WeightedGraph<Vertex>* forest = new WeightedGraph<Vertex>();
280 // Disjoint-set to keep track of the connected components in the
281 // maximum spanning tree.
282 HashMap<Vertex, Vertex> disjoint_set;
284 // Sort of the edges in the graph in decreasing order of their
285 // weight. Also add the vertices of the graph to the Maximum
286 // Spanning Tree graph and set each vertex to be its own connected
287 // component in the disjoint_set structure.
288 const HashSet<Vertex>& vertices = graph.vertices();
289 for (typename HashSet<Vertex>::const_iterator it = vertices.begin();
290 it != vertices.end();
292 const Vertex vertex1 = *it;
293 forest->AddVertex(vertex1, graph.VertexWeight(vertex1));
294 disjoint_set[vertex1] = vertex1;
296 const HashSet<Vertex>& neighbors = graph.Neighbors(vertex1);
297 for (typename HashSet<Vertex>::const_iterator it2 = neighbors.begin();
298 it2 != neighbors.end();
300 const Vertex vertex2 = *it2;
301 if (vertex1 >= vertex2) {
304 const double weight = graph.EdgeWeight(vertex1, vertex2);
305 weighted_edges.push_back(
306 std::make_pair(weight, std::make_pair(vertex1, vertex2)));
310 // The elements of this vector, are pairs<edge_weight,
311 // edge>. Sorting it using the reverse iterators gives us the edges
312 // in decreasing order of edges.
313 std::sort(weighted_edges.rbegin(), weighted_edges.rend());
315 // Greedily add edges to the spanning tree/forest as long as they do
316 // not violate the degree/cycle constraint.
317 for (int i =0; i < weighted_edges.size(); ++i) {
318 const std::pair<Vertex, Vertex>& edge = weighted_edges[i].second;
319 const Vertex vertex1 = edge.first;
320 const Vertex vertex2 = edge.second;
322 // Check if either of the vertices are of degree 2 already, in
323 // which case adding this edge will violate the degree 2
325 if ((forest->Neighbors(vertex1).size() == 2) ||
326 (forest->Neighbors(vertex2).size() == 2)) {
330 // Find the id of the connected component to which the two
331 // vertices belong to. If the id is the same, it means that the
332 // two of them are already connected to each other via some other
333 // vertex, and adding this edge will create a cycle.
334 Vertex root1 = FindConnectedComponent(vertex1, &disjoint_set);
335 Vertex root2 = FindConnectedComponent(vertex2, &disjoint_set);
337 if (root1 == root2) {
341 // This edge can be added, add an edge in either direction with
342 // the same weight as the original graph.
343 const double edge_weight = graph.EdgeWeight(vertex1, vertex2);
344 forest->AddEdge(vertex1, vertex2, edge_weight);
345 forest->AddEdge(vertex2, vertex1, edge_weight);
347 // Connected the two connected components by updating the
348 // disjoint_set structure. Always connect the connected component
349 // with the greater index with the connected component with the
350 // smaller index. This should ensure shallower trees, for quicker
353 std::swap(root1, root2);
356 disjoint_set[root2] = root1;
361 } // namespace internal
364 #endif // CERES_INTERNAL_GRAPH_ALGORITHMS_H_