1 // Ceres Solver - A fast non-linear least squares minimizer
2 // Copyright 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: moll.markus@arcor.de (Markus Moll)
30 //         sameeragarwal@google.com (Sameer Agarwal)
31 
32 #ifndef CERES_INTERNAL_POLYNOMIAL_SOLVER_H_
33 #define CERES_INTERNAL_POLYNOMIAL_SOLVER_H_
34 
35 #include <vector>
36 #include "ceres/internal/eigen.h"
37 #include "ceres/internal/port.h"
38 
39 namespace ceres {
40 namespace internal {
41 
42 // All polynomials are assumed to be the form
43 //
44 //   sum_{i=0}^N polynomial(i) x^{N-i}.
45 //
46 // and are given by a vector of coefficients of size N + 1.
47 
48 // Evaluate the polynomial at x using the Horner scheme.
EvaluatePolynomial(const Vector & polynomial,double x)49 inline double EvaluatePolynomial(const Vector& polynomial, double x) {
50   double v = 0.0;
51   for (int i = 0; i < polynomial.size(); ++i) {
52     v = v * x + polynomial(i);
53   }
54   return v;
55 }
56 
57 // Use the companion matrix eigenvalues to determine the roots of the
58 // polynomial.
59 //
60 // This function returns true on success, false otherwise.
61 // Failure indicates that the polynomial is invalid (of size 0) or
62 // that the eigenvalues of the companion matrix could not be computed.
63 // On failure, a more detailed message will be written to LOG(ERROR).
64 // If real is not NULL, the real parts of the roots will be returned in it.
65 // Likewise, if imaginary is not NULL, imaginary parts will be returned in it.
66 bool FindPolynomialRoots(const Vector& polynomial,
67                          Vector* real,
68                          Vector* imaginary);
69 
70 // Return the derivative of the given polynomial. It is assumed that
71 // the input polynomial is at least of degree zero.
72 Vector DifferentiatePolynomial(const Vector& polynomial);
73 
74 // Find the minimum value of the polynomial in the interval [x_min,
75 // x_max]. The minimum is obtained by computing all the roots of the
76 // derivative of the input polynomial. All real roots within the
77 // interval [x_min, x_max] are considered as well as the end points
78 // x_min and x_max. Since polynomials are differentiable functions,
79 // this ensures that the true minimum is found.
80 void MinimizePolynomial(const Vector& polynomial,
81                         double x_min,
82                         double x_max,
83                         double* optimal_x,
84                         double* optimal_value);
85 
86 // Structure for storing sample values of a function.
87 //
88 // Clients can use this struct to communicate the value of the
89 // function and or its gradient at a given point x.
90 struct FunctionSample {
FunctionSampleFunctionSample91   FunctionSample()
92       : x(0.0),
93         value(0.0),
94         value_is_valid(false),
95         gradient(0.0),
96         gradient_is_valid(false) {
97   }
98   string ToDebugString() const;
99 
100   double x;
101   double value;      // value = f(x)
102   bool value_is_valid;
103   double gradient;   // gradient = f'(x)
104   bool gradient_is_valid;
105 };
106 
107 // Given a set of function value and/or gradient samples, find a
108 // polynomial whose value and gradients are exactly equal to the ones
109 // in samples.
110 //
111 // Generally speaking,
112 //
113 // degree = # values + # gradients - 1
114 //
115 // Of course its possible to sample a polynomial any number of times,
116 // in which case, generally speaking the spurious higher order
117 // coefficients will be zero.
118 Vector FindInterpolatingPolynomial(const vector<FunctionSample>& samples);
119 
120 // Interpolate the function described by samples with a polynomial,
121 // and minimize it on the interval [x_min, x_max]. Depending on the
122 // input samples, it is possible that the interpolation or the root
123 // finding algorithms may fail due to numerical difficulties. But the
124 // function is guaranteed to return its best guess of an answer, by
125 // considering the samples and the end points as possible solutions.
126 void MinimizeInterpolatingPolynomial(const vector<FunctionSample>& samples,
127                                      double x_min,
128                                      double x_max,
129                                      double* optimal_x,
130                                      double* optimal_value);
131 
132 }  // namespace internal
133 }  // namespace ceres
134 
135 #endif  // CERES_INTERNAL_POLYNOMIAL_SOLVER_H_
136