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