how to calculate the QR decomposition of a matrix

In many papers(such as [1], [2]), we need to solve the polar decomposition of a matrix which decompose a matrix A into a rotation matrix R and a symmetric matrix S, s.t. A = RS

Here, we will be using Eigen library to solve this problem contained in paper [2].

In Eigen library, there is no method can solve polar decomposition directly. Thus, we must find another ways to do this. Fortunately, up to now, I have already found two ways to deal with this problem. So let's cut into our theme:

First, I would like to introduce something about  how to obtain matrix R, S from the polar decomposition, please look at the following picture:

Then CGAL source code is presented:

    bool polar_eigen(const Matrix& A, Matrix& R)
    {
      if(A.determinant() < 0)
      { return false; }

      typedef Matrix::Scalar Scalar;

      const Scalar th = std::sqrt(Eigen::NumTraits<Scalar>::dummy_precision());

      Eigen::SelfAdjointEigenSolver<Matrix> eig;
      CGAL::feclearexcept(FE_UNDERFLOW);
      eig.computeDirect(A.transpose()*A);//solve eigen values and eigen vectors
      if(CGAL::fetestexcept(FE_UNDERFLOW) || eig.eigenvalues()(0)/eig.eigenvalues()(2)<th) // to make sure the matrix is semi-positive
      { return false; }

      Vector S = eig.eigenvalues().cwiseSqrt();
      R = A  * eig.eigenvectors() * S.asDiagonal().inverse()
        * eig.eigenvectors().transpose();

      if(std::abs(R.squaredNorm()-3.) > th || R.determinant() < 0)// because R is an orthonormal matrix the squared norms of which column vectors should be 1
      { return false; }

      R.transposeInPlace(); // the optimal rotation matrix should be transpose of decomposition result, only in paper[2], generally we don't need tranpose it
      return true;
    }

feclearexcept and fetestexcept are used to detect some excetion  which are defined in the source code below:

// Copyright (c) 2010-2011  GeometryFactory Sarl (France)
//
// This file is part of CGAL (www.cgal.org); 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.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// $URL$
// $Id$
//
//
// Author(s)     : Laurent Rineau

#include <cfloat>

#ifndef FE_INVALID
#  define FE_INEXACT     _EM_INEXACT
#  define FE_UNDERFLOW   _EM_UNDERFLOW
#  define FE_OVERFLOW    _EM_OVERFLOW
#  define FE_DIVBYZERO   _EM_ZERODIVIDE
#  define FE_INVALID     _EM_INVALID
#endif

namespace CGAL {

// replacement for C99 functions

inline int
feclearexcept(int exceptions) {
  _clearfp();
  return 0;
}

inline int
fetestexcept(int exceptions)
{
#if defined(_M_IX86) && _M_IX86_FP > 0
  // On x86/x64 processors, when SSE is used.
  unsigned int i1;
  unsigned int i2;
  _statusfp2(&i1, &i2);
  return (i1 & exceptions) | (i2 & exceptions);
#else
  // On x86 processors without SSE, or on other processors supported by MSVC
  return _statusfp() & exceptions;
#endif
}

} // end namespace CGAL

the second way is to use single value decomposition which written in CGAL like below:

  /// Computes the closest rotation to `m` and places it into `R`
  void compute_close_rotation(const Matrix& m, Matrix& R)
  {
    Eigen::JacobiSVD<Eigen::Matrix3d> solver;
    solver.compute( m, Eigen::ComputeFullU | Eigen::ComputeFullV );

    const Matrix& u = solver.matrixU(); const Matrix& v = solver.matrixV();
    R = v * u.transpose();

    if( R.determinant() < 0 ) {
      Matrix u_copy = u;
      u_copy.col(2) *= -1;        // singular values sorted ascendingly
      R = v * u_copy.transpose(); // re-extract rotation matrix
    }
  }

Here R is denoted as the way in the paper [2] needs. Normally R should be expressed as  R = U * V.transpose();

[1] Meshless Deformations Based on Shape Matching

[2] As-Rigid-As-Possible surface Modeling