add Gram-Schmidt process #11

.gitignore

@@ -4,3 +4,5 @@ include/Eigen/*
 .idea/
 cmake-build*/
 build/*
+.vscode
+*~

include/Spectra/LinAlg/Orthogonalization.h

@@ -0,0 +1,68 @@
+// Copyright (C) 2020 Netherlands eScience Center <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#ifndef SPECTRA_ORTHOGONALIZATION_H
+#define SPECTRA_ORTHOGONALIZATION_H
+
+#include <Eigen/Core>
+
+namespace Spectra {
+
+/// Gram-Schmidt process to orthoganlize a given basis A. Each column correspond
+/// to a vector on the basis.
+/// The start index indicates from what vector to start
+/// the orthogonalization, 0 is the default.
+/// Warnings: Starting the normalization at index n imples that the previous
+/// n-1 vectors are already orthonormal.
+
+template <typename Scalar>
+class Orthogonalization
+{
+private:
+    using Index = Eigen::Index;
+    using Matrix = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
+    using Vector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
+    const Matrix& basis;
+    Index start_index = 0;
+
+public:
+    Orthogonalization(const Matrix& A, Index nstart) :
+        basis{A}, start_index{nstart} {}
+
+    Orthogonalization(const Matrix& A) :
+        basis{A} {}
+
+    Matrix orthonormalize()
+    {
+        Matrix Q = basis;
+        Index nstart = start_index;
+        if (start_index == 0)
+        {
+            Q.col(0).normalize();
+            nstart = 1;
+        }
+        for (Index j = nstart; j < basis.cols(); ++j)
+        {
+            Q.col(j) -= Q.leftCols(j) * (Q.leftCols(j).transpose() * basis.col(j));
+            Q.col(j).normalize();
+            // two consecutive Gram-schmidt iterations are enough
+            // http://stoppels.blog/posts/orthogonalization-performance
+            Q.col(j) -= Q.leftCols(j) * (Q.leftCols(j).transpose() * Q.col(j));
+            if (Q.col(j).norm() <= 1E-12 * basis.col(j).norm())
+            {
+                throw std::runtime_error(
+                    "There is a Linear dependencies in Gram-Schmidt."
+                    "Hint: try the QR method.");
+            }
+            Q.col(j).normalize();
+        }
+        return Q;
+    }
+};
+
+}  // namespace Spectra
+
+#endif  //SPECTRA_ORTHOGONALIZATION_H

test/CMakeLists.txt

@@ -9,6 +9,7 @@ list(APPEND test_target_sources
         GenEigs.cpp
         GenEigsRealShift.cpp
         GenEigsComplexShift.cpp
+        Orthogonalization.cpp
         SymGEigsCholesky.cpp
         SymGEigsRegInv.cpp
         SVD.cpp

test/Orthogonalization.cpp

@@ -0,0 +1,48 @@
+#include <Eigen/Core>
+#include <Spectra/LinAlg/Orthogonalization.h>
+
+using namespace Spectra;
+
+#define CATCH_CONFIG_MAIN
+#include "catch.hpp"
+
+using Eigen::MatrixXd;
+using Eigen::VectorXd;
+using Eigen::Index;
+
+template <typename Matrix>
+void check_orthogonality(const Matrix& basis)
+{
+    const double tol = 1e-12;
+    Matrix xs = basis.transpose() * basis;
+    INFO("The orthonormalized basis must fulfill that basis.T * basis = I");
+    REQUIRE(xs.isIdentity(tol));
+}
+
+TEST_CASE("complete orthonormalization", "[Gram-Schmidt]")
+{
+    std::srand(123);
+    const Index n = 100;
+
+    MatrixXd mat = MatrixXd::Random(n, n);
+    Orthogonalization<double> gs{mat};
+    check_orthogonality(gs.orthonormalize());
+}
+
+TEST_CASE("Partial orthonormalization", "[Gram-Schmidt]")
+{
+    std::srand(123);
+    const Index n = 100;
+
+    // Create a n x 20 orthonormal basis
+    MatrixXd mat = MatrixXd::Random(n, n - 20);
+    Orthogonalization<double> gs{mat};
+    mat.leftCols(n - 20) = gs.orthonormalize();
+
+    mat.conservativeResize(Eigen::NoChange, n);
+    mat.rightCols(20) = MatrixXd::Random(n, 20);
+
+    // Orthogonalize from 80 onwards
+    Orthogonalization<double> new_gs{mat, 80};
+    check_orthogonality(new_gs.orthonormalize());
+}