Hessenberg decomposition of a matrix

geometry/r3x3_matrix_test.cpp

@@ -65,9 +65,9 @@ TEST_F(R3x3MatrixTest, QRDecomposition) {
   EXPECT_THAT(
       q,
       AlmostEquals(R3x3Matrix<double>(
-          {6.0 / 7.0, 3.0 / 7.0, 2.0 / 7.0},
-          {-30.0 / Sqrt(3577.0), 34.0 / Sqrt(3577.0), 39.0 / Sqrt(3577.0)},
-          {1.0 / Sqrt(73.0), -6.0 / Sqrt(73.0), 6.0 / Sqrt(73.0)}), 10));
+          {6.0 / 7.0, -30.0 / Sqrt(3577.0),  1.0 / Sqrt(73.0)},
+          {3.0 / 7.0,  34.0 / Sqrt(3577.0), -6.0 / Sqrt(73.0)},
+          {2.0 / 7.0,  39.0 / Sqrt(3577.0), 6.0 / Sqrt(73.0)}), 569));
   EXPECT_THAT(
       r,
       AlmostEquals(R3x3Matrix<Length>(

numerics/fixed_arrays.hpp

@@ -19,6 +19,8 @@ using namespace principia::numerics::_transposed_view;
 using namespace principia::quantities::_named_quantities;
 using namespace principia::quantities::_si;
 
+// TODO(phl): This should support the same operations as unbounded_arrays.hpp.
+
 template<typename Scalar, int rows, int columns>
 class FixedMatrix;
 
@@ -221,10 +223,10 @@ template<typename Scalar, int size>
 constexpr FixedVector<double, size> Normalize(
     FixedVector<Scalar, size> const& vector);
 
-template<typename LScalar, typename RScalar, int lsize, int rsize>
-constexpr FixedMatrix<Product<LScalar, RScalar>, lsize, rsize> SymmetricProduct(
-    FixedVector<LScalar, lsize> const& left,
-    FixedVector<RScalar, rsize> const& right);
+template<typename LScalar, typename RScalar, int size>
+constexpr FixedMatrix<Product<LScalar, RScalar>, size, size> SymmetricProduct(
+    FixedVector<LScalar, size> const& left,
+    FixedVector<RScalar, size> const& right);
 
 template<typename Scalar, int size>
 constexpr FixedMatrix<Square<Scalar>, size, size> SymmetricSquare(
@@ -284,23 +286,23 @@ constexpr FixedMatrix<Scalar, rows, columns>& operator-=(
 
 template<typename LScalar, typename RScalar, int size>
 constexpr FixedVector<Product<LScalar, RScalar>, size> operator*(
-    LScalar const left,
+    LScalar const& left,
     FixedVector<RScalar, size> const& right);
 
 template<typename LScalar, typename RScalar, int size>
 constexpr FixedVector<Product<LScalar, RScalar>, size> operator*(
     FixedVector<LScalar, size> const& left,
-    RScalar const right);
+    RScalar const& right);
 
 template<typename LScalar, typename RScalar, int rows, int columns>
 constexpr FixedMatrix<Product<LScalar, RScalar>, rows, columns>
-operator*(LScalar const left,
+operator*(LScalar const& left,
           FixedMatrix<RScalar, rows, columns> const& right);
 
 template<typename LScalar, typename RScalar, int rows, int columns>
 constexpr FixedMatrix<Product<LScalar, RScalar>, rows, columns>
 operator*(FixedMatrix<LScalar, rows, columns> const& left,
-          RScalar const right);
+          RScalar const& right);
 
 template<typename LScalar, typename RScalar, int size>
 constexpr FixedVector<Quotient<LScalar, RScalar>, size> operator/(
@@ -310,7 +312,27 @@ constexpr FixedVector<Quotient<LScalar, RScalar>, size> operator/(
 template<typename LScalar, typename RScalar, int rows, int columns>
 constexpr FixedMatrix<Quotient<LScalar, RScalar>, rows, columns>
 operator/(FixedMatrix<LScalar, rows, columns> const& left,
-          RScalar const right);
+          RScalar const& right);
+
+template<typename Scalar, int size>
+constexpr FixedVector<Scalar, size>& operator*=(
+    FixedVector<Scalar, size>& left,
+    double right);
+
+template<typename Scalar, int rows, int columns>
+constexpr FixedMatrix<Scalar, rows, columns>& operator*=(
+    FixedMatrix<Scalar, rows, columns>& left,
+    double right);
+
+template<typename Scalar, int size>
+constexpr FixedVector<Scalar, size>& operator/=(
+    FixedVector<Scalar, size>& left,
+    double right);
+
+template<typename Scalar, int rows, int columns>
+constexpr FixedMatrix<Scalar, rows, columns>& operator/=(
+    FixedMatrix<Scalar, rows, columns>& left,
+    double right);
 
 // Hilbert space and algebra.
 
@@ -341,6 +363,12 @@ constexpr FixedVector<Product<LScalar, RScalar>, rows> operator*(
     FixedMatrix<LScalar, rows, columns> const& left,
     FixedVector<RScalar, columns> const& right);
 
+// Use this operator to multiply a row vector with a matrix.
+template<typename LScalar, typename RScalar, int rows, int columns>
+constexpr FixedVector<Product<LScalar, RScalar>, columns> operator*(
+    TransposedView<FixedMatrix<LScalar, rows, columns>> const& left,
+    FixedVector<RScalar, rows> const& right);
+
 // Ouput.
 
 template<typename Scalar, int size>

numerics/fixed_arrays_body.hpp

@@ -83,7 +83,7 @@ Scalar FixedVector<Scalar, size_>::Norm() const {
 }
 
 template<typename Scalar, int size_>
-inline Square<Scalar> FixedVector<Scalar, size_>::Norm²() const {
+Square<Scalar> FixedVector<Scalar, size_>::Norm²() const {
   return DotProduct<Scalar, Scalar, size_>::Compute(data_, data_);
 }
 
@@ -397,12 +397,12 @@ constexpr FixedVector<double, size> Normalize(
   return vector / vector.Norm();
 }
 
-template<typename LScalar, typename RScalar, int lsize, int rsize>
-constexpr FixedMatrix<Product<LScalar, RScalar>, lsize, rsize> SymmetricProduct(
-    FixedVector<LScalar, lsize> const& left,
-    FixedVector<RScalar, rsize> const& right) {
-  FixedMatrix<Product<LScalar, RScalar>, lsize, rsize> result(uninitialized);
-  for (int i = 0; i < lsize; ++i) {
+template<typename LScalar, typename RScalar, int size>
+constexpr FixedMatrix<Product<LScalar, RScalar>, size, size> SymmetricProduct(
+    FixedVector<LScalar, size> const& left,
+    FixedVector<RScalar, size> const& right) {
+  FixedMatrix<Product<LScalar, RScalar>, size, size> result(uninitialized);
+  for (int i = 0; i < size; ++i) {
     for (int j = 0; j < i; ++j) {
       auto const r = 0.5 * (left[i] * right[j] + left[j] * right[i]);
       result(i, j) = r;
@@ -533,7 +533,7 @@ constexpr FixedMatrix<Scalar, rows, columns>& operator-=(
 
 template<typename LScalar, typename RScalar, int size>
 constexpr FixedVector<Product<LScalar, RScalar>, size> operator*(
-    LScalar const left,
+    LScalar const& left,
     FixedVector<RScalar, size> const& right) {
   std::array<Product<LScalar, RScalar>, size> result;
   for (int i = 0; i < size; ++i) {
@@ -545,7 +545,7 @@ constexpr FixedVector<Product<LScalar, RScalar>, size> operator*(
 template<typename LScalar, typename RScalar, int size>
 constexpr FixedVector<Product<LScalar, RScalar>, size> operator*(
     FixedVector<LScalar, size> const& left,
-    RScalar const right) {
+    RScalar const& right) {
   std::array<Product<LScalar, RScalar>, size> result;
   for (int i = 0; i < size; ++i) {
     result[i] = left[i] * right;
@@ -555,7 +555,7 @@ constexpr FixedVector<Product<LScalar, RScalar>, size> operator*(
 
 template<typename LScalar, typename RScalar, int rows, int columns>
 constexpr FixedMatrix<Product<LScalar, RScalar>, rows, columns> operator*(
-    LScalar const left,
+    LScalar const& left,
     FixedMatrix<RScalar, rows, columns> const& right) {
   FixedMatrix<Product<LScalar, RScalar>, rows, columns> result(
       uninitialized);
@@ -570,7 +570,7 @@ constexpr FixedMatrix<Product<LScalar, RScalar>, rows, columns> operator*(
 template<typename LScalar, typename RScalar, int rows, int columns>
 constexpr FixedMatrix<Product<LScalar, RScalar>, rows, columns> operator*(
     FixedMatrix<LScalar, rows, columns> const& left,
-    RScalar const right) {
+    RScalar const& right) {
   FixedMatrix<Product<LScalar, RScalar>, rows, columns> result(
       uninitialized);
   for (int i = 0; i < rows; ++i) {
@@ -595,7 +595,7 @@ constexpr FixedVector<Quotient<LScalar, RScalar>, size> operator/(
 template<typename LScalar, typename RScalar, int rows, int columns>
 constexpr FixedMatrix<Quotient<LScalar, RScalar>, rows, columns> operator/(
     FixedMatrix<LScalar, rows, columns> const& left,
-    RScalar const right) {
+    RScalar const& right) {
   FixedMatrix<Quotient<LScalar, RScalar>, rows, columns> result(
       uninitialized);
   for (int i = 0; i < rows; ++i) {
@@ -606,6 +606,32 @@ constexpr FixedMatrix<Quotient<LScalar, RScalar>, rows, columns> operator/(
   return result;
 }
 
+template<typename Scalar, int size>
+constexpr FixedVector<Scalar, size>& operator*=(FixedVector<Scalar, size>& left,
+                                                double const right) {
+  return left = left * right;
+}
+
+template<typename Scalar, int rows, int columns>
+constexpr FixedMatrix<Scalar, rows, columns>& operator*=(
+    FixedMatrix<Scalar, rows, columns>& left,
+    double const right) {
+  return left = left * right;
+}
+
+template<typename Scalar, int size>
+constexpr FixedVector<Scalar, size>& operator/=(FixedVector<Scalar, size>& left,
+                                                double const right) {
+  return left = left / right;
+}
+
+template<typename Scalar, int rows, int columns>
+constexpr FixedMatrix<Scalar, rows, columns>& operator/=(
+    FixedMatrix<Scalar, rows, columns>& left,
+    double const right) {
+  return left = left / right;
+}
+
 template<typename LScalar, typename RScalar, int size>
 constexpr Product<LScalar, RScalar> operator*(
     LScalar* const left,
@@ -664,6 +690,20 @@ constexpr FixedVector<Product<LScalar, RScalar>, rows> operator*(
   return FixedVector<Product<LScalar, RScalar>, rows>(std::move(result));
 }
 
+template<typename LScalar, typename RScalar, int rows, int columns>
+constexpr FixedVector<Product<LScalar, RScalar>, columns> operator*(
+    TransposedView<FixedMatrix<LScalar, rows, columns>> const& left,
+    FixedVector<RScalar, rows> const& right) {
+  std::array<Product<LScalar, RScalar>, columns> result{};
+  for (int j = 0; j < columns; ++j) {
+    auto& result_j = result[j];
+    for (int i = 0; i < rows; ++i) {
+      result_j += left.transpose(i, j) * right[i];
+    }
+  }
+  return FixedVector<Product<LScalar, RScalar>, columns>(std::move(result));
+}
+
 template<typename Scalar, int size>
 std::ostream& operator<<(std::ostream& out,
                          FixedVector<Scalar, size> const& vector) {

numerics/fixed_arrays_test.cpp

@@ -130,6 +130,8 @@ TEST_F(FixedArraysTest, Algebra) {
                                         14, -63,   5, -92})),
             m23_ * m34_);
   EXPECT_EQ(v3_, m34_ * v4_);
+  EXPECT_EQ((FixedVector<double, 4>({-486, -229, 333, 198})),
+            TransposedView{m34_} * v3_);  // NOLINT
 }
 
 TEST_F(FixedArraysTest, VectorIndexing) {

numerics/matrix_computations.hpp

@@ -29,6 +29,15 @@ struct ᵗRDRDecompositionGenerator;
 template<typename M, typename V>
 struct SubstitutionGenerator;
 
+// Declares:
+//   struct Result {
+//     (matrix) H;
+//     (matrix) U;
+//   }
+// TODO(phl): Add support for U.
+template<typename M>
+struct HessenbergDecompositionGenerator;
+
 // Declares:
 //   struct Result {
 //     ⟨matrix⟩ rotation;
@@ -81,6 +90,13 @@ typename SubstitutionGenerator<LowerTriangularMatrix, Vector>::Result
 ForwardSubstitution(LowerTriangularMatrix const& L,
                     Vector const& b);
 
+// If A is a square matrix, returns U and H so that A = ᵗU H U, where H is an
+// upper Hessenberg matrix.
+// TODO(phl): Add support for returning U.
+template<typename Matrix>
+typename HessenbergDecompositionGenerator<Matrix>::Result
+HessenbergDecomposition(Matrix const& A);
+
 // Returns the eigensystem of A, which must be symmetric.
 // As a safety measure we limit the number of iterations.  We prefer to exit
 // when the matrix is nearly diagonal, though.
@@ -95,7 +111,8 @@ template<typename Matrix, typename Vector>
 typename RayleighQuotientGenerator<Matrix, Vector>::Result
 RayleighQuotient(Matrix const& A, Vector const& x);
 
-// Returns the eigenvector closest to x and its eigenvalue.
+// Returns the eigenvector closest to x and its eigenvalue.  A must be
+// symmetric.
 template<typename Matrix, typename Vector>
 typename RayleighQuotientIterationGenerator<Matrix, Vector>::Result
 RayleighQuotientIteration(Matrix const& A, Vector const& x);
@@ -111,6 +128,7 @@ using internal::BackSubstitution;
 using internal::CholeskyDecomposition;
 using internal::ClassicalJacobi;
 using internal::ForwardSubstitution;
+using internal::HessenbergDecomposition;
 using internal::RayleighQuotient;
 using internal::RayleighQuotientIteration;
 using internal::Solve;