1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. Eigen itself is part of the KDE project. 3 // 4 // Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1 (at) 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> 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 VectorType remean = VectorType::Zero(size); 151 for(int i = 0; i < numPoints; ++i) 152 { 153 VectorType diff = (*(points[i]) - mean).conjugate(); 154 covMat += diff * diff.adjoint(); 155 } 156 157 // now we just have to pick the eigen vector with smallest eigen value 158 SelfAdjointEigenSolver<CovMatrixType> eig(covMat); 159 result->normal() = eig.eigenvectors().col(0); 160 if (soundness) 161 *soundness = eig.eigenvalues().coeff(0)/eig.eigenvalues().coeff(1); 162 163 // let's compute the constant coefficient such that the 164 // plane pass trough the mean point: 165 result->offset() = - (result->normal().cwise()* mean).sum(); 166 } 167 168 } // end namespace Eigen 169 170 #endif // EIGEN2_LEASTSQUARES_H 171
