graph_algorithms.h 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344
  1. // Ceres Solver - A fast non-linear least squares minimizer
  2. // Copyright 2023 Google Inc. All rights reserved.
  3. // http://ceres-solver.org/
  4. //
  5. // Redistribution and use in source and binary forms, with or without
  6. // modification, are permitted provided that the following conditions are met:
  7. //
  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.
  16. //
  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.
  28. //
  29. // Author: sameeragarwal@google.com (Sameer Agarwal)
  30. //
  31. // Various algorithms that operate on undirected graphs.
  32. #ifndef CERES_INTERNAL_GRAPH_ALGORITHMS_H_
  33. #define CERES_INTERNAL_GRAPH_ALGORITHMS_H_
  34. #include <algorithm>
  35. #include <memory>
  36. #include <unordered_map>
  37. #include <unordered_set>
  38. #include <utility>
  39. #include <vector>
  40. #include "ceres/graph.h"
  41. #include "ceres/internal/export.h"
  42. #include "ceres/wall_time.h"
  43. #include "glog/logging.h"
  44. namespace ceres::internal {
  45. // Compare two vertices of a graph by their degrees, if the degrees
  46. // are equal then order them by their ids.
  47. template <typename Vertex>
  48. class CERES_NO_EXPORT VertexTotalOrdering {
  49. public:
  50. explicit VertexTotalOrdering(const Graph<Vertex>& graph) : graph_(graph) {}
  51. bool operator()(const Vertex& lhs, const Vertex& rhs) const {
  52. if (graph_.Neighbors(lhs).size() == graph_.Neighbors(rhs).size()) {
  53. return lhs < rhs;
  54. }
  55. return graph_.Neighbors(lhs).size() < graph_.Neighbors(rhs).size();
  56. }
  57. private:
  58. const Graph<Vertex>& graph_;
  59. };
  60. template <typename Vertex>
  61. class VertexDegreeLessThan {
  62. public:
  63. explicit VertexDegreeLessThan(const Graph<Vertex>& graph) : graph_(graph) {}
  64. bool operator()(const Vertex& lhs, const Vertex& rhs) const {
  65. return graph_.Neighbors(lhs).size() < graph_.Neighbors(rhs).size();
  66. }
  67. private:
  68. const Graph<Vertex>& graph_;
  69. };
  70. // Order the vertices of a graph using its (approximately) largest
  71. // independent set, where an independent set of a graph is a set of
  72. // vertices that have no edges connecting them. The maximum
  73. // independent set problem is NP-Hard, but there are effective
  74. // approximation algorithms available. The implementation here uses a
  75. // breadth first search that explores the vertices in order of
  76. // increasing degree. The same idea is used by Saad & Li in "MIQR: A
  77. // multilevel incomplete QR preconditioner for large sparse
  78. // least-squares problems", SIMAX, 2007.
  79. //
  80. // Given a undirected graph G(V,E), the algorithm is a greedy BFS
  81. // search where the vertices are explored in increasing order of their
  82. // degree. The output vector ordering contains elements of S in
  83. // increasing order of their degree, followed by elements of V - S in
  84. // increasing order of degree. The return value of the function is the
  85. // cardinality of S.
  86. template <typename Vertex>
  87. int IndependentSetOrdering(const Graph<Vertex>& graph,
  88. std::vector<Vertex>* ordering) {
  89. const std::unordered_set<Vertex>& vertices = graph.vertices();
  90. const int num_vertices = vertices.size();
  91. CHECK(ordering != nullptr);
  92. ordering->clear();
  93. ordering->reserve(num_vertices);
  94. // Colors for labeling the graph during the BFS.
  95. const char kWhite = 0;
  96. const char kGrey = 1;
  97. const char kBlack = 2;
  98. // Mark all vertices white.
  99. std::unordered_map<Vertex, char> vertex_color;
  100. std::vector<Vertex> vertex_queue;
  101. for (const Vertex& vertex : vertices) {
  102. vertex_color[vertex] = kWhite;
  103. vertex_queue.push_back(vertex);
  104. }
  105. std::sort(vertex_queue.begin(),
  106. vertex_queue.end(),
  107. VertexTotalOrdering<Vertex>(graph));
  108. // Iterate over vertex_queue. Pick the first white vertex, add it
  109. // to the independent set. Mark it black and its neighbors grey.
  110. for (const Vertex& vertex : vertex_queue) {
  111. if (vertex_color[vertex] != kWhite) {
  112. continue;
  113. }
  114. ordering->push_back(vertex);
  115. vertex_color[vertex] = kBlack;
  116. const std::unordered_set<Vertex>& neighbors = graph.Neighbors(vertex);
  117. for (const Vertex& neighbor : neighbors) {
  118. vertex_color[neighbor] = kGrey;
  119. }
  120. }
  121. int independent_set_size = ordering->size();
  122. // Iterate over the vertices and add all the grey vertices to the
  123. // ordering. At this stage there should only be black or grey
  124. // vertices in the graph.
  125. for (const Vertex& vertex : vertex_queue) {
  126. DCHECK(vertex_color[vertex] != kWhite);
  127. if (vertex_color[vertex] != kBlack) {
  128. ordering->push_back(vertex);
  129. }
  130. }
  131. CHECK_EQ(ordering->size(), num_vertices);
  132. return independent_set_size;
  133. }
  134. // Same as above with one important difference. The ordering parameter
  135. // is an input/output parameter which carries an initial ordering of
  136. // the vertices of the graph. The greedy independent set algorithm
  137. // starts by sorting the vertices in increasing order of their
  138. // degree. The input ordering is used to stabilize this sort, i.e., if
  139. // two vertices have the same degree then they are ordered in the same
  140. // order in which they occur in "ordering".
  141. //
  142. // This is useful in eliminating non-determinism from the Schur
  143. // ordering algorithm over all.
  144. template <typename Vertex>
  145. int StableIndependentSetOrdering(const Graph<Vertex>& graph,
  146. std::vector<Vertex>* ordering) {
  147. CHECK(ordering != nullptr);
  148. const std::unordered_set<Vertex>& vertices = graph.vertices();
  149. const int num_vertices = vertices.size();
  150. CHECK_EQ(vertices.size(), ordering->size());
  151. // Colors for labeling the graph during the BFS.
  152. const char kWhite = 0;
  153. const char kGrey = 1;
  154. const char kBlack = 2;
  155. std::vector<Vertex> vertex_queue(*ordering);
  156. std::stable_sort(vertex_queue.begin(),
  157. vertex_queue.end(),
  158. VertexDegreeLessThan<Vertex>(graph));
  159. // Mark all vertices white.
  160. std::unordered_map<Vertex, char> vertex_color;
  161. for (const Vertex& vertex : vertices) {
  162. vertex_color[vertex] = kWhite;
  163. }
  164. ordering->clear();
  165. ordering->reserve(num_vertices);
  166. // Iterate over vertex_queue. Pick the first white vertex, add it
  167. // to the independent set. Mark it black and its neighbors grey.
  168. for (int i = 0; i < vertex_queue.size(); ++i) {
  169. const Vertex& vertex = vertex_queue[i];
  170. if (vertex_color[vertex] != kWhite) {
  171. continue;
  172. }
  173. ordering->push_back(vertex);
  174. vertex_color[vertex] = kBlack;
  175. const std::unordered_set<Vertex>& neighbors = graph.Neighbors(vertex);
  176. for (const Vertex& neighbor : neighbors) {
  177. vertex_color[neighbor] = kGrey;
  178. }
  179. }
  180. int independent_set_size = ordering->size();
  181. // Iterate over the vertices and add all the grey vertices to the
  182. // ordering. At this stage there should only be black or grey
  183. // vertices in the graph.
  184. for (const Vertex& vertex : vertex_queue) {
  185. DCHECK(vertex_color[vertex] != kWhite);
  186. if (vertex_color[vertex] != kBlack) {
  187. ordering->push_back(vertex);
  188. }
  189. }
  190. CHECK_EQ(ordering->size(), num_vertices);
  191. return independent_set_size;
  192. }
  193. // Find the connected component for a vertex implemented using the
  194. // find and update operation for disjoint-set. Recursively traverse
  195. // the disjoint set structure till you reach a vertex whose connected
  196. // component has the same id as the vertex itself. Along the way
  197. // update the connected components of all the vertices. This updating
  198. // is what gives this data structure its efficiency.
  199. template <typename Vertex>
  200. Vertex FindConnectedComponent(const Vertex& vertex,
  201. std::unordered_map<Vertex, Vertex>* union_find) {
  202. auto it = union_find->find(vertex);
  203. DCHECK(it != union_find->end());
  204. if (it->second != vertex) {
  205. it->second = FindConnectedComponent(it->second, union_find);
  206. }
  207. return it->second;
  208. }
  209. // Compute a degree two constrained Maximum Spanning Tree/forest of
  210. // the input graph. Caller owns the result.
  211. //
  212. // Finding degree 2 spanning tree of a graph is not always
  213. // possible. For example a star graph, i.e. a graph with n-nodes
  214. // where one node is connected to the other n-1 nodes does not have
  215. // a any spanning trees of degree less than n-1.Even if such a tree
  216. // exists, finding such a tree is NP-Hard.
  217. // We get around both of these problems by using a greedy, degree
  218. // constrained variant of Kruskal's algorithm. We start with a graph
  219. // G_T with the same vertex set V as the input graph G(V,E) but an
  220. // empty edge set. We then iterate over the edges of G in decreasing
  221. // order of weight, adding them to G_T if doing so does not create a
  222. // cycle in G_T} and the degree of all the vertices in G_T remains
  223. // bounded by two. This O(|E|) algorithm results in a degree-2
  224. // spanning forest, or a collection of linear paths that span the
  225. // graph G.
  226. template <typename Vertex>
  227. std::unique_ptr<WeightedGraph<Vertex>> Degree2MaximumSpanningForest(
  228. const WeightedGraph<Vertex>& graph) {
  229. // Array of edges sorted in decreasing order of their weights.
  230. std::vector<std::pair<double, std::pair<Vertex, Vertex>>> weighted_edges;
  231. auto forest = std::make_unique<WeightedGraph<Vertex>>();
  232. // Disjoint-set to keep track of the connected components in the
  233. // maximum spanning tree.
  234. std::unordered_map<Vertex, Vertex> disjoint_set;
  235. // Sort of the edges in the graph in decreasing order of their
  236. // weight. Also add the vertices of the graph to the Maximum
  237. // Spanning Tree graph and set each vertex to be its own connected
  238. // component in the disjoint_set structure.
  239. const std::unordered_set<Vertex>& vertices = graph.vertices();
  240. for (const Vertex& vertex1 : vertices) {
  241. forest->AddVertex(vertex1, graph.VertexWeight(vertex1));
  242. disjoint_set[vertex1] = vertex1;
  243. const std::unordered_set<Vertex>& neighbors = graph.Neighbors(vertex1);
  244. for (const Vertex& vertex2 : neighbors) {
  245. if (vertex1 >= vertex2) {
  246. continue;
  247. }
  248. const double weight = graph.EdgeWeight(vertex1, vertex2);
  249. weighted_edges.push_back(
  250. std::make_pair(weight, std::make_pair(vertex1, vertex2)));
  251. }
  252. }
  253. // The elements of this vector, are pairs<edge_weight,
  254. // edge>. Sorting it using the reverse iterators gives us the edges
  255. // in decreasing order of edges.
  256. std::sort(weighted_edges.rbegin(), weighted_edges.rend());
  257. // Greedily add edges to the spanning tree/forest as long as they do
  258. // not violate the degree/cycle constraint.
  259. for (int i = 0; i < weighted_edges.size(); ++i) {
  260. const std::pair<Vertex, Vertex>& edge = weighted_edges[i].second;
  261. const Vertex vertex1 = edge.first;
  262. const Vertex vertex2 = edge.second;
  263. // Check if either of the vertices are of degree 2 already, in
  264. // which case adding this edge will violate the degree 2
  265. // constraint.
  266. if ((forest->Neighbors(vertex1).size() == 2) ||
  267. (forest->Neighbors(vertex2).size() == 2)) {
  268. continue;
  269. }
  270. // Find the id of the connected component to which the two
  271. // vertices belong to. If the id is the same, it means that the
  272. // two of them are already connected to each other via some other
  273. // vertex, and adding this edge will create a cycle.
  274. Vertex root1 = FindConnectedComponent(vertex1, &disjoint_set);
  275. Vertex root2 = FindConnectedComponent(vertex2, &disjoint_set);
  276. if (root1 == root2) {
  277. continue;
  278. }
  279. // This edge can be added, add an edge in either direction with
  280. // the same weight as the original graph.
  281. const double edge_weight = graph.EdgeWeight(vertex1, vertex2);
  282. forest->AddEdge(vertex1, vertex2, edge_weight);
  283. forest->AddEdge(vertex2, vertex1, edge_weight);
  284. // Connected the two connected components by updating the
  285. // disjoint_set structure. Always connect the connected component
  286. // with the greater index with the connected component with the
  287. // smaller index. This should ensure shallower trees, for quicker
  288. // lookup.
  289. if (root2 < root1) {
  290. std::swap(root1, root2);
  291. }
  292. disjoint_set[root2] = root1;
  293. }
  294. return forest;
  295. }
  296. } // namespace ceres::internal
  297. #endif // CERES_INTERNAL_GRAPH_ALGORITHMS_H_