1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN2_LEASTSQUARES_H
11 #define EIGEN2_LEASTSQUARES_H
12 
13 namespace Eigen {
14 
15 /** \ingroup LeastSquares_Module
16   *
17   * \leastsquares_module
18   *
19   * For a set of points, this function tries to express
20   * one of the coords as a linear (affine) function of the other coords.
21   *
22   * This is best explained by an example. This function works in full
23   * generality, for points in a space of arbitrary dimension, and also over
24   * the complex numbers, but for this example we will work in dimension 3
25   * over the real numbers (doubles).
26   *
27   * So let us work with the following set of 5 points given by their
28   * \f$(x,y,z)\f$ coordinates:
29   * @code
30     Vector3d points[5];
31     points[0] = Vector3d( 3.02, 6.89, -4.32 );
32     points[1] = Vector3d( 2.01, 5.39, -3.79 );
33     points[2] = Vector3d( 2.41, 6.01, -4.01 );
34     points[3] = Vector3d( 2.09, 5.55, -3.86 );
35     points[4] = Vector3d( 2.58, 6.32, -4.10 );
36   * @endcode
37   * Suppose that we want to express the second coordinate (\f$y\f$) as a linear
38   * expression in \f$x\f$ and \f$z\f$, that is,
39   * \f[ y=ax+bz+c \f]
40   * for some constants \f$a,b,c\f$. Thus, we want to find the best possible
41   * constants \f$a,b,c\f$ so that the plane of equation \f$y=ax+bz+c\f$ fits
42   * best the five above points. To do that, call this function as follows:
43   * @code
44     Vector3d coeffs; // will store the coefficients a, b, c
45     linearRegression(
46       5,
47       &points,
48       &coeffs,
49       1 // the coord to express as a function of
50         // the other ones. 0 means x, 1 means y, 2 means z.
51     );
52   * @endcode
53   * Now the vector \a coeffs is approximately
54   * \f$( 0.495 ,  -1.927 ,  -2.906 )\f$.
55   * Thus, we get \f$a=0.495, b = -1.927, c = -2.906\f$. Let us check for
56   * instance how near points[0] is from the plane of equation \f$y=ax+bz+c\f$.
57   * Looking at the coords of points[0], we see that:
58   * \f[ax+bz+c = 0.495 * 3.02 + (-1.927) * (-4.32) + (-2.906) = 6.91.\f]
59   * On the other hand, we have \f$y=6.89\f$. We see that the values
60   * \f$6.91\f$ and \f$6.89\f$
61   * are near, so points[0] is very near the plane of equation \f$y=ax+bz+c\f$.
62   *
63   * Let's now describe precisely the parameters:
64   * @param numPoints the number of points
65   * @param points the array of pointers to the points on which to perform the linear regression
66   * @param result pointer to the vector in which to store the result.
67                   This vector must be of the same type and size as the
68                   data points. The meaning of its coords is as follows.
69                   For brevity, let \f$n=Size\f$,
70                   \f$r_i=result[i]\f$,
71                   and \f$f=funcOfOthers\f$. Denote by
72                   \f$x_0,\ldots,x_{n-1}\f$
73                   the n coordinates in the n-dimensional space.
74                   Then the resulting equation is:
75                   \f[ x_f = r_0 x_0 + \cdots + r_{f-1}x_{f-1}
76                    + r_{f+1}x_{f+1} + \cdots + r_{n-1}x_{n-1} + r_n. \f]
77   * @param funcOfOthers Determines which coord to express as a function of the
78                         others. Coords are numbered starting from 0, so that a
79                         value of 0 means \f$x\f$, 1 means \f$y\f$,
80                         2 means \f$z\f$, ...
81   *
82   * \sa fitHyperplane()
83   */
84 template<typename VectorType>
linearRegression(int numPoints,VectorType ** points,VectorType * result,int funcOfOthers)85 void linearRegression(int numPoints,
86                       VectorType **points,
87                       VectorType *result,
88                       int funcOfOthers )
89 {
90   typedef typename VectorType::Scalar Scalar;
91   typedef Hyperplane<Scalar, VectorType::SizeAtCompileTime> HyperplaneType;
92   const int size = points[0]->size();
93   result->resize(size);
94   HyperplaneType h(size);
95   fitHyperplane(numPoints, points, &h);
96   for(int i = 0; i < funcOfOthers; i++)
97     result->coeffRef(i) = - h.coeffs()[i] / h.coeffs()[funcOfOthers];
98   for(int i = funcOfOthers; i < size; i++)
99     result->coeffRef(i) = - h.coeffs()[i+1] / h.coeffs()[funcOfOthers];
100 }
101 
102 /** \ingroup LeastSquares_Module
103   *
104   * \leastsquares_module
105   *
106   * This function is quite similar to linearRegression(), so we refer to the
107   * documentation of this function and only list here the differences.
108   *
109   * The main difference from linearRegression() is that this function doesn't
110   * take a \a funcOfOthers argument. Instead, it finds a general equation
111   * of the form
112   * \f[ r_0 x_0 + \cdots + r_{n-1}x_{n-1} + r_n = 0, \f]
113   * where \f$n=Size\f$, \f$r_i=retCoefficients[i]\f$, and we denote by
114   * \f$x_0,\ldots,x_{n-1}\f$ the n coordinates in the n-dimensional space.
115   *
116   * Thus, the vector \a retCoefficients has size \f$n+1\f$, which is another
117   * difference from linearRegression().
118   *
119   * In practice, this function performs an hyper-plane fit in a total least square sense
120   * via the following steps:
121   *  1 - center the data to the mean
122   *  2 - compute the covariance matrix
123   *  3 - pick the eigenvector corresponding to the smallest eigenvalue of the covariance matrix
124   * The ratio of the smallest eigenvalue and the second one gives us a hint about the relevance
125   * of the solution. This value is optionally returned in \a soundness.
126   *
127   * \sa linearRegression()
128   */
129 template<typename VectorType, typename HyperplaneType>
130 void fitHyperplane(int numPoints,
131                    VectorType **points,
132                    HyperplaneType *result,
133                    typename NumTraits<typename VectorType::Scalar>::Real* soundness = 0)
134 {
135   typedef typename VectorType::Scalar Scalar;
136   typedef Matrix<Scalar,VectorType::SizeAtCompileTime,VectorType::SizeAtCompileTime> CovMatrixType;
137   EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType)
138   ei_assert(numPoints >= 1);
139   int size = points[0]->size();
140   ei_assert(size+1 == result->coeffs().size());
141 
142   // compute the mean of the data
143   VectorType mean = VectorType::Zero(size);
144   for(int i = 0; i < numPoints; ++i)
145     mean += *(points[i]);
146   mean /= numPoints;
147 
148   // compute the covariance matrix
149   CovMatrixType covMat = CovMatrixType::Zero(size, size);
150   for(int i = 0; i < numPoints; ++i)
151   {
152     VectorType diff = (*(points[i]) - mean).conjugate();
153     covMat += diff * diff.adjoint();
154   }
155 
156   // now we just have to pick the eigen vector with smallest eigen value
157   SelfAdjointEigenSolver<CovMatrixType> eig(covMat);
158   result->normal() = eig.eigenvectors().col(0);
159   if (soundness)
160     *soundness = eig.eigenvalues().coeff(0)/eig.eigenvalues().coeff(1);
161 
162   // let's compute the constant coefficient such that the
163   // plane pass trough the mean point:
164   result->offset() = - (result->normal().cwise()* mean).sum();
165 }
166 
167 } // end namespace Eigen
168 
169 #endif // EIGEN2_LEASTSQUARES_H
170