Unitriangular Gram-Schmidt decomposition

documentation/bibliography.bib

@@ -1650,6 +1650,14 @@ @book{Higham2002
   title     = {Accuracy and Stability of Numerical Algorithms},
 }
 
+@book{HoffsteinPipherSilverman2014,
+  author    = {Hoffstein, Jeffrey and Pipher, Jill and Silverman, Joseph H.},
+  publisher = {Springer},
+  date      = {2014},
+  series    = {Undergraduate Texts in Mathematics},
+  title     = {An Introduction to Mathematical Cryptography},
+}
+
 @book{Householder1970,
   author    = {Householder, Alston Scott},
   publisher = {McGraw-Hill},

numerics/matrix_computations.hpp

@@ -37,6 +37,14 @@ struct SubstitutionGenerator;
 template<typename M>
 struct GramSchmidtGenerator;
 
+// Declares:
+//   struct Result {
+//     ⟨upper triangular matrix⟩ R;
+//     ⟨matrix⟩ Q;
+//   };
+template<typename M>
+struct UnitriangularGramSchmidtGenerator;
+
 // Declares:
 //   struct Result {
 //     (matrix) H;
@@ -107,9 +115,17 @@ typename SubstitutionGenerator<LowerTriangularMatrix, Vector>::Result
 ForwardSubstitution(LowerTriangularMatrix const& L,
                     Vector const& b);
 
+// Returns Q and R such that A = Q R where R is upper triangular and Q
+// orthogonal.
 template<typename Matrix>
 typename GramSchmidtGenerator<Matrix>::Result
-ClassicalGramSchmidt(Matrix const& L);
+ClassicalGramSchmidt(Matrix const& A);
+
+// Returns Q and R such that A = Q R where R is upper unitriangular and Q
+// has orthogonal columns.
+template<typename Matrix>
+typename UnitriangularGramSchmidtGenerator<Matrix>::Result
+UnitriangularGramSchmidt(Matrix const& A);
 
 // If A is a square matrix, returns U and H so that A = ᵗU H U, where H is an
 // upper Hessenberg matrix.
@@ -161,6 +177,7 @@ using internal::RayleighQuotient;
 using internal::RayleighQuotientIteration;
 using internal::RealSchurDecomposition;
 using internal::Solve;
+using internal::UnitriangularGramSchmidt;
 using internal::ᵗRDRDecomposition;
 
 }  // namespace _matrix_computations

numerics/matrix_computations_body.hpp

@@ -308,6 +308,28 @@ struct GramSchmidtGenerator<FixedMatrix<Scalar, dimension, dimension>> {
       FixedMatrix<Scalar, dimension, dimension> const& m);
 };
 
+template<typename Scalar>
+struct UnitriangularGramSchmidtGenerator<UnboundedMatrix<Scalar>> {
+  struct Result {
+    UnboundedMatrix<Scalar> Q;
+    UnboundedUpperTriangularMatrix<double> R;
+  };
+  using QVector = UnboundedVector<Scalar>;
+  static Result Uninitialized(UnboundedMatrix<Scalar> const& m);
+};
+
+template<typename Scalar, int dimension>
+struct UnitriangularGramSchmidtGenerator<
+    FixedMatrix<Scalar, dimension, dimension>> {
+  struct Result {
+    FixedMatrix<Scalar, dimension, dimension> Q;
+    FixedUpperTriangularMatrix<double, dimension> R;
+  };
+  using QVector = FixedVector<Scalar, dimension>;
+  static Result Uninitialized(
+      FixedMatrix<Scalar, dimension, dimension> const& m);
+};
+
 template<typename Scalar>
 struct HessenbergDecompositionGenerator<UnboundedMatrix<Scalar>> {
   struct Result {
@@ -471,6 +493,23 @@ Uninitialized(FixedMatrix<Scalar, dimension, dimension> const& m) -> Result {
       .R = FixedUpperTriangularMatrix<Scalar, dimension>(uninitialized)};
 }
 
+template<typename Scalar>
+auto UnitriangularGramSchmidtGenerator<UnboundedMatrix<Scalar>>::
+Uninitialized(UnboundedMatrix<Scalar> const& m) -> Result {
+  return Result{
+      .Q = UnboundedMatrix<Scalar>(m.rows(), m.columns(), uninitialized),
+      .R = UnboundedUpperTriangularMatrix<double>(m.columns(), uninitialized)};
+}
+
+template<typename Scalar, int dimension>
+auto UnitriangularGramSchmidtGenerator<
+    FixedMatrix<Scalar, dimension, dimension>>::
+Uninitialized(FixedMatrix<Scalar, dimension, dimension> const& m) -> Result {
+  return Result{
+      .Q = FixedMatrix<Scalar, dimension, dimension>(uninitialized),
+      .R = FixedUpperTriangularMatrix<double, dimension>(uninitialized)};
+}
+
 template<typename Scalar>
 auto ClassicalJacobiGenerator<UnboundedMatrix<Scalar>>::Identity(
     UnboundedMatrix<Scalar> const& m) -> Rotation {
@@ -680,6 +719,50 @@ typename GramSchmidtGenerator<Matrix>::Result ClassicalGramSchmidt(
   return result;
 }
 
+template<typename Matrix>
+typename UnitriangularGramSchmidtGenerator<Matrix>::Result
+UnitriangularGramSchmidt(Matrix const& A) {
+  using G = UnitriangularGramSchmidtGenerator<Matrix>;
+  auto result = G::Uninitialized(A);
+  auto& Q = result.Q;
+  auto& R = result.R;
+  int const n = A.rows();
+
+  // The algorithm is derived from [HPS14], Theorem 7.13, but we keep the code
+  // similar to the one for classical Gram-Schmidt, i.e., adopt the conventions
+  // from [Hig02], Algorithm 19.11.
+  // The correspondances are as follows:
+  //   vⱼ  ≘ aⱼ
+  //   v⭑ⱼ ≘ qⱼ
+  //   μᵢⱼ ≘ Rⱼᵢ
+  for (int j = 0; j < n; ++j) {
+    auto const aⱼ = ColumnView<Matrix const>{.matrix = A,
+                                             .first_row = 0,
+                                             .last_row = n - 1,
+                                             .column = j};
+    for (int i = 0; i < j; ++i) {
+      auto const qᵢ = ColumnView{.matrix = Q,
+                                 .first_row = 0,
+                                 .last_row = n - 1,
+                                 .column = i};
+      R(i, j) = TransposedView{.transpose = qᵢ} * aⱼ / qᵢ.Norm²();  // NOLINT
+    }
+    auto qⱼ = ColumnView{.matrix = Q,  // NOLINT
+                         .first_row = 0,
+                         .last_row = n - 1,
+                         .column = j};
+    qⱼ = aⱼ;
+    for (int k = 0; k < j; ++k) {
+      qⱼ -= R(k, j) * typename G::QVector(ColumnView{.matrix = Q,
+                                                     .first_row = 0,
+                                                     .last_row = n - 1,
+                                                     .column = k});
+    }
+    R(j, j) = 1;
+  }
+  return result;
+}
+
 template<typename Matrix>
 typename HessenbergDecompositionGenerator<Matrix>::Result
 HessenbergDecomposition(Matrix const& A) {

numerics/matrix_computations_test.cpp

@@ -5,6 +5,7 @@
 #include "gmock/gmock.h"
 #include "gtest/gtest.h"
 #include "numerics/fixed_arrays.hpp"
+#include "numerics/matrix_views.hpp"
 #include "numerics/transposed_view.hpp"
 #include "numerics/unbounded_arrays.hpp"
 #include "quantities/elementary_functions.hpp"
@@ -19,6 +20,7 @@ namespace numerics {
 
 using namespace principia::numerics::_fixed_arrays;
 using namespace principia::numerics::_matrix_computations;
+using namespace principia::numerics::_matrix_views;
 using namespace principia::numerics::_transposed_view;
 using namespace principia::numerics::_unbounded_arrays;
 using namespace principia::quantities::_elementary_functions;
@@ -137,6 +139,40 @@ TYPED_TEST(MatrixComputationsTest, ClassicalGramSchmidt) {
   }
 }
 
+TYPED_TEST(MatrixComputationsTest, UnitriangularGramSchmidt) {
+  using Matrix = typename std::tuple_element<3, TypeParam>::type;
+  Matrix const m4({1, 2, 3, -4,
+                   5, 6, 7, 8,
+                   9, 8, -7, 6,
+                   5, 4, 3, 2});
+  auto const qr = UnitriangularGramSchmidt(m4);
+
+  // Check that the decomposition is correct.
+  auto const near_m4 = qr.Q * Matrix(qr.R);
+  EXPECT_THAT(near_m4, AlmostEquals(m4, 1));
+
+  // Check that R is unitriangular.
+  for (int i = 0; i < qr.R.rows(); ++i) {
+    EXPECT_THAT(qr.R(i, i), AlmostEquals(1, 0));
+  }
+
+  // Check that the columns of Q are approximately orthogonal.
+  for (int c1 = 0; c1 < qr.Q.columns(); ++c1) {
+    auto const column_c1 = ColumnView{.matrix = qr.Q,
+                                      .first_row = 0,
+                                      .last_row = qr.Q.rows() - 1,
+                                      .column = c1};
+    for (int c2 = c1 + 1; c2 < qr.Q.columns(); ++c2) {
+      auto const column_c2 = ColumnView{.matrix = qr.Q,
+                                        .first_row = 0,
+                                        .last_row = qr.Q.rows() - 1,
+                                        .column = c2};
+      EXPECT_THAT(TransposedView{.transpose = column_c1} * column_c2,  // NOLINT
+                  VanishesBefore(1, 24, 176));
+    }
+  }
+}
+
 TYPED_TEST(MatrixComputationsTest, HessenbergForm) {
   using Matrix = typename std::tuple_element<3, TypeParam>::type;
   Matrix const m4({1, 2, 3, -4,