From 053261de88f062f8870b5f550177c7991f8a7738 Mon Sep 17 00:00:00 2001
From: Hauke Heibel <hauke.heibel@gmail.com>
Date: Fri, 24 Sep 2010 16:28:20 +0200
Subject: [PATCH] Make the SVD's output unitary and improved unit tests.

---
 Eigen/src/SVD/SVD.h | 13 +++++++++++++
 test/svd.cpp        | 12 ++++++------
 2 files changed, 19 insertions(+), 6 deletions(-)

diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h
index 15e284dea..b8f0c6496 100644
--- a/Eigen/src/SVD/SVD.h
+++ b/Eigen/src/SVD/SVD.h
@@ -445,6 +445,19 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
   m_matU.leftCols(n) = A;
   m_matU.rightCols(m-n).setZero();
 
+  // Gram Schmidt orthogonalization to fill up U
+  for (int col = A.cols(); col < A.rows(); ++col)
+  {
+    typename MatrixUType::ColXpr colVec = m_matU.col(col);
+    colVec(col) = 1;
+    for (int prevCol = 0; prevCol < col; ++prevCol)
+    {
+      typename MatrixUType::ColXpr prevColVec = m_matU.col(prevCol);
+      colVec -= colVec.dot(prevColVec)*prevColVec;
+    }
+    m_matU.col(col) = colVec.normalized();
+  }
+
   m_isInitialized = true;
   return *this;
 }
diff --git a/test/svd.cpp b/test/svd.cpp
index 8ba510d9b..5aa766fee 100644
--- a/test/svd.cpp
+++ b/test/svd.cpp
@@ -47,9 +47,9 @@ template<typename MatrixType> void svd(const MatrixType& m)
 
     sigma.diagonal() = svd.singularValues();
     matU = svd.matrixU();
-    //VERIFY_IS_UNITARY(matU);
+    VERIFY_IS_UNITARY(matU);
     matV = svd.matrixV();
-    //VERIFY_IS_UNITARY(matV);
+    VERIFY_IS_UNITARY(matV);
     VERIFY_IS_APPROX(a, matU * sigma * matV.transpose());
   }
 
@@ -57,10 +57,10 @@ template<typename MatrixType> void svd(const MatrixType& m)
   if (rows>=cols)
   {
     SVD<MatrixType> svd(a);
-    Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> x = Matrix<Scalar, MatrixType::ColsAtCompileTime, 1>::Random(cols,1);
-    Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> b = a * x;
-    Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> result = svd.solve(b);
-    VERIFY_IS_APPROX(a * result, b);
+    Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> b = Matrix<Scalar, MatrixType::ColsAtCompileTime, 1>::Random(rows,1);
+    Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> x = svd.solve(b);
+    // evaluate normal equation which works also for least-squares solutions
+    VERIFY_IS_APPROX(a.adjoint()*a*x,a.adjoint()*b);
   }