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#endifnamespace CGAL {// replacement for C99 functionsinline intfeclearexcept(int exceptions) {  _clearfp();  return 0;}inline intfetestexcept(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