Home | History | Annotate | Download | only in ceres
      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