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