title: Confused about Eigen QR decomposition
body: <p>I am confused about Eigen's QR decomposition. My understanding is that the matrix <code>Q</code> is stored implicitly as a sequence of Householder transformations, and that the matrix <code>R</code> is stored as an upper triangular matrix, and that the diagonal of <code>R</code> contains the eigenvalues of <code>A</code> (at least up to phase, which is all I care about).</p>
  
  <p>However, I wrote the following program which computes the eigenvalues of a matrix <code>A</code> via two different methods, one using the <code>Eigen::EigenSolver</code>, and the other using <code>QR</code>. I know that my <code>QR</code> method is returning the wrong results, and that the <code>EigenSolver</code> method is returning the correct results.</p>
  
  <p>What am I misunderstanding here?</p>
  
  <pre class="lang-cpp prettyprint-override"><code>#include &lt;iostream&gt;
  #include &lt;algorithm&gt;
  #include &lt;Eigen/Dense&gt;
  
  int main()
  {
      using Real = long double;
      long n = 2;
      Eigen::Matrix&lt;Real, Eigen::Dynamic, Eigen::Dynamic&gt; A(n,n);
  
      for(long i = 0; i &lt; n; ++i) {
          for (long j = 0; j &lt; n; ++j) {
              A(i,j) = Real(1)/Real(i+j+1);
          }
      }
  
      auto QR = A.householderQr();
  
      auto Rdiag = QR.matrixQR().diagonal().cwiseAbs();
      auto [min, max] = std::minmax_element(Rdiag.begin(), Rdiag.end());
      std::cout &lt;&lt; "\u03BA\u2082(A) = " &lt;&lt; (*max)/(*min) &lt;&lt; "\n";
  
      std::cout &lt;&lt; "\u2016A\u2016\u2082 via QR = " &lt;&lt; (*max) &lt;&lt; "\n";
      std::cout &lt;&lt; "Diagonal of R =\n" &lt;&lt; Rdiag &lt;&lt; "\n";
  
  
      // dblcheck:
  
      Eigen::SelfAdjointEigenSolver&lt;decltype(A)&gt; eigensolver(A);
      if (eigensolver.info() != Eigen::Success) {
          std::cout &lt;&lt; "Something went wrong.\n";
          return 1;
      }
      auto absolute_eigs = eigensolver.eigenvalues().cwiseAbs();
      auto [min1, max1] = std::minmax_element(absolute_eigs.begin(), absolute_eigs.end());
      std::cout &lt;&lt; "\u03BA\u2082(A) via eigensolver = " &lt;&lt; (*max1)/(*min1) &lt;&lt; "\n";
      std::cout &lt;&lt; "\u2016A\u2016\u2082 via eigensolver = " &lt;&lt; (*max1) &lt;&lt; "\n";
      std::cout &lt;&lt; "The absolute eigenvalues of A via eigensolver are:\n" &lt;&lt; absolute_eigs &lt;&lt; "\n";
  }
  </code></pre>
  
  <p>Output:</p>
  
  <pre class="lang-sh prettyprint-override"><code>κ₂(A) = 15
  ‖A‖₂ via QR = 1.11803
  Diagonal of R =
    1.11803
  0.0745356
  κ₂(A) via eigensolver = 19.2815
  ‖A‖₂ via eigensolver = 1.26759
  The absolute eigenvalues of A via eigensolver are:
  0.0657415
    1.26759
  </code></pre>
  
  <p>Other info:</p>
  
  <ul>
  <li>Cloned eigen from the head:</li>
  </ul>
  
  <pre class="lang-sh prettyprint-override"><code>$ hg log | more
  changeset:   11993:20cbc6576426
  tag:         tip
  date:        Tue May 07 16:44:55 2019 -0700
  summary:     Fix AVX512 &amp; GCC 6.3 compilation
  </code></pre>
  
  <ul>
  <li><p>Occurs when compiled with g++-8, g++-9, and Apple Clang, with and without <code>-ffast-math</code>. I obtain the same wrong result with <code>Eigen::FullPivHouseholderQR</code>.</p></li>
  <li><p>I also looked into the source <code>HouseholderQR.h</code>, and they compute the determinant via <code>m_qr.diagonal().prod()</code>, which makes me feel somewhat more confident that I'm using the API correctly. Taking the product of the eigenvalues from the EigenSolver returns the same values as <code>QR.absDeterminant()</code>.</p></li>
  <li><p>The following code snippet <em>does not return the original matrix A</em>:</p></li>
  </ul>
  
  <pre class="lang-cpp prettyprint-override"><code>Eigen::Matrix&lt;Real, Eigen::Dynamic, Eigen::Dynamic&gt; R = QR.matrixQR();
  Eigen::Matrix&lt;Real, Eigen::Dynamic, Eigen::Dynamic&gt; Q = QR.householderQ();
  std::cout &lt;&lt; "Q*R = \n" &lt;&lt; Q*R &lt;&lt; "\n";
  </code></pre>
  
  <ul>
  <li><p>I checked that <code>Q</code> has all the requisite properties: Q^-1 = Q^T, Q^TQ = I, and |det(Q)| = 1.</p></li>
  <li><p>I've also verified that <code>QR.householderQ().transpose()*QR.matrixQR()</code> is not equal <code>A</code>; although one column is correct and another is wrong. </p></li>
  </ul>