LeastSquares.h

Go to the documentation of this file.

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
//
// Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.

#ifndef EIGEN2_LEASTSQUARES_H
#define EIGEN2_LEASTSQUARES_H

template<typename VectorType>
void linearRegression(int numPoints,
                      VectorType **points,
                      VectorType *result,
                      int funcOfOthers )
{
  typedef typename VectorType::Scalar Scalar;
  typedef Hyperplane<Scalar, VectorType::SizeAtCompileTime> HyperplaneType;
  const int size = points[0]->size();
  result->resize(size);
  HyperplaneType h(size);
  fitHyperplane(numPoints, points, &h);
  for(int i = 0; i < funcOfOthers; i++)
    result->coeffRef(i) = - h.coeffs()[i] / h.coeffs()[funcOfOthers];
  for(int i = funcOfOthers; i < size; i++)
    result->coeffRef(i) = - h.coeffs()[i+1] / h.coeffs()[funcOfOthers];
}

template<typename VectorType, typename HyperplaneType>
void fitHyperplane(int numPoints,
                   VectorType **points,
                   HyperplaneType *result,
                   typename NumTraits<typename VectorType::Scalar>::Real* soundness = 0)
{
  typedef typename VectorType::Scalar Scalar;
  typedef Matrix<Scalar,VectorType::SizeAtCompileTime,VectorType::SizeAtCompileTime> CovMatrixType;
  EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType)
  ei_assert(numPoints >= 1);
  int size = points[0]->size();
  ei_assert(size+1 == result->coeffs().size());

  // compute the mean of the data
  VectorType mean = VectorType::Zero(size);
  for(int i = 0; i < numPoints; ++i)
    mean += *(points[i]);
  mean /= numPoints;

  // compute the covariance matrix
  CovMatrixType covMat = CovMatrixType::Zero(size, size);
  VectorType remean = VectorType::Zero(size);
  for(int i = 0; i < numPoints; ++i)
  {
    VectorType diff = (*(points[i]) - mean).conjugate();
    covMat += diff * diff.adjoint();
  }

  // now we just have to pick the eigen vector with smallest eigen value
  SelfAdjointEigenSolver<CovMatrixType> eig(covMat);
  result->normal() = eig.eigenvectors().col(0);
  if (soundness)
    *soundness = eig.eigenvalues().coeff(0)/eig.eigenvalues().coeff(1);

  // let's compute the constant coefficient such that the
  // plane pass trough the mean point:
  result->offset() = - (result->normal().cwise()* mean).sum();
}


#endif // EIGEN2_LEASTSQUARES_H