diff --git a/geometry/r3x3_matrix_test.cpp b/geometry/r3x3_matrix_test.cpp
index 3513d00146..266a4eda82 100644
--- a/geometry/r3x3_matrix_test.cpp
+++ b/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>(
diff --git a/numerics/fixed_arrays.hpp b/numerics/fixed_arrays.hpp
index 21885606a4..4eed4af7e4 100644
--- a/numerics/fixed_arrays.hpp
+++ b/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>
diff --git a/numerics/fixed_arrays_body.hpp b/numerics/fixed_arrays_body.hpp
index 27cc2d2224..f519c2a574 100644
--- a/numerics/fixed_arrays_body.hpp
+++ b/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) {
diff --git a/numerics/fixed_arrays_test.cpp b/numerics/fixed_arrays_test.cpp
index c84b1dabf0..92dc97a2d7 100644
--- a/numerics/fixed_arrays_test.cpp
+++ b/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) {
diff --git a/numerics/matrix_computations.hpp b/numerics/matrix_computations.hpp
index 73ba275b03..3b9da3ecfa 100644
--- a/numerics/matrix_computations.hpp
+++ b/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;
diff --git a/numerics/matrix_computations_body.hpp b/numerics/matrix_computations_body.hpp
index 6673538b57..84a77c2086 100644
--- a/numerics/matrix_computations_body.hpp
+++ b/numerics/matrix_computations_body.hpp
@@ -8,6 +8,7 @@
 
 #include "base/tags.hpp"
 #include "numerics/fixed_arrays.hpp"
+#include "numerics/transposed_view.hpp"
 #include "numerics/unbounded_arrays.hpp"
 #include "quantities/elementary_functions.hpp"
 #include "quantities/si.hpp"
@@ -19,10 +20,251 @@ namespace internal {
 
 using namespace principia::base::_tags;
 using namespace principia::numerics::_fixed_arrays;
+using namespace principia::numerics::_transposed_view;
 using namespace principia::numerics::_unbounded_arrays;
 using namespace principia::quantities::_elementary_functions;
 using namespace principia::quantities::_si;
 
+// TODO(phl): The view stuff should be (1) made completed, i.e., have all the
+// operations that exist for fixed/unbounded vectors/matrices; (2) moved to a
+// common place (probably together with TransposedView); (3) unified with
+// fixed/unbounded arrays so that we don't have to write each algorithm N times.
+
+// A view of a column of a matrix.
+// TODO(phl): Could we extract Scalar from Matrix?
+template<typename Scalar, typename Matrix>
+struct ColumnView {
+  Matrix& matrix;
+  int first_row;
+  int last_row;
+  int column;
+
+  Scalar Norm() const;
+  Square<Scalar> Norm²() const;
+  constexpr int size() const;
+
+  // Constructs an unbounded vector by copying data from the view.  Note that
+  // the result is unbounded even if the matrix being viewed is a FixedMatrix.
+  explicit operator UnboundedVector<Scalar>() const;
+
+  constexpr Scalar& operator[](int index);
+  constexpr Scalar const& operator[](int index) const;
+
+  ColumnView& operator/=(double right);
+};
+
+template<typename Scalar, typename Matrix>
+Square<Scalar> ColumnView<Scalar, Matrix>::Norm²() const {
+  Square<Scalar> result{};
+  for (int i = first_row; i <= last_row; ++i) {
+    result += Pow<2>(matrix(i, column));
+  }
+  return result;
+}
+
+template<typename Scalar, typename Matrix>
+Scalar ColumnView<Scalar, Matrix>::Norm() const {
+  return Sqrt(Norm²());
+}
+
+template<typename Scalar, typename Matrix>
+constexpr int ColumnView<Scalar, Matrix>::size() const {
+  return last_row - first_row + 1;
+}
+
+template<typename Scalar, typename Matrix>
+ColumnView<Scalar, Matrix>::operator UnboundedVector<Scalar>() const {
+  UnboundedVector<Scalar> result(size(), uninitialized);
+  for (int i = first_row; i <= last_row; ++i) {
+    result[i - first_row] = matrix(i, column);
+  }
+  return result;
+}
+
+template<typename Scalar, typename Matrix>
+constexpr Scalar& ColumnView<Scalar, Matrix>::operator[](int const index) {
+  DCHECK_LE(index, last_row - first_row);
+  return matrix(first_row + index, column);
+}
+
+template<typename Scalar, typename Matrix>
+constexpr Scalar const& ColumnView<Scalar, Matrix>::operator[](
+    int const index) const {
+  DCHECK_LE(index, last_row - first_row);
+  return matrix(first_row + index, column);
+}
+
+template<typename Scalar, typename Matrix>
+ColumnView<Scalar, Matrix>& ColumnView<Scalar, Matrix>::operator/=(
+    double const right) {
+  for (int i = first_row; i < last_row; ++i) {
+    matrix(i, column) /= right;
+  }
+}
+
+template<typename Scalar, typename Matrix>
+UnboundedVector<double> Normalize(ColumnView<Scalar, Matrix> const& view) {
+  return UnboundedVector<Scalar>(view) / view.Norm();
+}
+
+template<typename Scalar, typename Matrix>
+std::ostream& operator<<(std::ostream& out,
+                         ColumnView<Scalar, Matrix> const& view) {
+  std::stringstream s;
+  for (int i = 0; i < view.size(); ++i) {
+    s << (i == 0 ? "{" : "") << view[i]
+      << (i == view.size() - 1 ? "}" : ", ");
+  }
+  out << s.str();
+  return out;
+}
+
+template<typename Scalar, typename Matrix>
+struct BlockView {
+  Matrix& matrix;
+  int first_row;
+  int last_row;
+  int first_column;
+  int last_column;
+
+  constexpr int rows() const;
+  constexpr int columns() const;
+
+  constexpr Scalar& operator()(int row, int column);
+  constexpr Scalar const& operator()(int row, int column) const;
+
+  BlockView& operator-=(UnboundedMatrix<Scalar> const& right);
+};
+
+template<typename Scalar, typename Matrix>
+constexpr int BlockView<Scalar, Matrix>::rows() const {
+  return last_row - first_row + 1;
+}
+
+template<typename Scalar, typename Matrix>
+constexpr int BlockView<Scalar, Matrix>::columns() const {
+  return last_column - first_column + 1;
+}
+
+template<typename Scalar, typename Matrix>
+constexpr Scalar& BlockView<Scalar, Matrix>::operator()(int const row,
+                                                        int const column) {
+  DCHECK_LE(row, last_row - first_row);
+  DCHECK_LE(column, last_column - first_column);
+  return matrix(first_row + row, first_column + column);
+}
+
+template<typename Scalar, typename Matrix>
+constexpr Scalar const& BlockView<Scalar, Matrix>::operator()(
+    int const row,
+    int const column) const {
+  DCHECK_LE(row, last_row - first_row);
+  DCHECK_LE(column, last_column - first_column);
+  return matrix(first_row + row, first_column + column);
+}
+
+template<typename Scalar, typename Matrix>
+BlockView<Scalar, Matrix>& BlockView<Scalar, Matrix>::operator-=(
+    UnboundedMatrix<Scalar> const& right) {
+  CHECK_EQ(rows(), right.rows());
+  CHECK_EQ(columns(), right.columns());
+  for (int i = 0; i < right.rows(); ++i) {
+    for (int j = 0; j < right.columns(); ++j) {
+      matrix(first_row + i, first_column + j) -= right(i, j);
+    }
+  }
+  return *this;
+}
+
+
+template<typename LScalar, typename RScalar, typename Matrix>
+UnboundedVector<Product<LScalar, RScalar>> operator*(
+    BlockView<LScalar, Matrix> const& left,
+    UnboundedVector<RScalar> const& right) {
+  CHECK_EQ(left.columns(), right.size());
+  UnboundedVector<Product<LScalar, RScalar>> result(left.rows());
+  for (int i = 0; i < left.rows(); ++i) {
+    auto& result_i = result[i];
+    for (int j = 0; j < left.columns(); ++j) {
+      result_i += left(i, j) * right[j];
+    }
+  }
+  return result;
+}
+
+template<typename LScalar, typename RScalar, typename Matrix>
+UnboundedVector<Product<LScalar, RScalar>> operator*(
+    TransposedView<BlockView<LScalar, Matrix>> const& left,
+    UnboundedVector<RScalar> const& right) {
+  CHECK_EQ(left.transpose.rows(), right.size());
+  UnboundedVector<Product<LScalar, RScalar>> result(left.transpose.columns());
+  for (int j = 0; j < left.transpose.columns(); ++j) {
+    auto& result_j = result[j];
+    for (int i = 0; i < left.transpose.rows(); ++i) {
+      result_j += left.transpose(i, j) * right[i];
+    }
+  }
+  return result;
+}
+
+// As mentioned in [GV13] section 5.1.4, "It is critical to exploit structure
+// when applying [the Householder reflection] to a matrix".
+struct HouseholderReflection {
+  UnboundedVector<double> v;
+  double β;
+};
+
+template<typename Vector>
+HouseholderReflection ComputeHouseholderReflection(Vector const& x) {
+  int const m = x.size();
+  // In order to avoid issues with quantities, we start by normalizing x.  This
+  // implies that μ is 1.
+  auto const normalized_x = Normalize(x);
+  HouseholderReflection result{.v = normalized_x, .β = 0};
+  double const& x₁ = normalized_x[0];
+  double& v₁ = result.v[0];
+  auto x₂ₘ = normalized_x;
+  x₂ₘ[0] = 0;
+  auto const σ = x₂ₘ.Norm²();
+  v₁ = 1;
+  if (σ == 0) {
+    if (x₁ < 0) {
+      result.β = -2;
+    }
+  } else {
+    static constexpr double μ = 1;
+    if (x₁ <= 0) {
+      v₁ = x₁ - μ;
+    } else {
+      v₁ = -σ / (x₁ + μ);
+    }
+    double const v₁² = Pow<2>(v₁);
+    result.β = 2 * v₁² / (σ + v₁²);
+    result.v /= v₁;
+  }
+  return result;
+}
+
+// A becomes P A.
+template<typename Matrix>
+void Premultiply(HouseholderReflection const& P, Matrix& A) {
+  // We don't have a multiplication TransposedView<Vector> * Matrix because the
+  // ownership of the result is problematic.  Instead, we transpose twice.  That
+  // costs essentially nothing.
+  auto const ᵗAv = TransposedView{A} * P.v;  // NOLINT
+  auto const ᵗvA = TransposedView{ᵗAv};      // NOLINT
+  auto const βv = P.β * P.v;
+  A -= βv * ᵗvA;
+}
+
+// A becomes A P.
+template<typename Matrix>
+void PostMultiply(Matrix& A, HouseholderReflection const& P) {
+  auto const βᵗv = TransposedView{P.β * P.v};  // NOLINT
+  auto const Av = A * P.v;
+  A -= Av * βᵗv;
+}
+
 // This is J(p, q, θ) in [GV13] section 8.5.1.  This matrix is also called a
 // Givens rotation.  As mentioned in [GV13] section 5.1.9, "It is critical that
 // the special structure of a Givens rotation matrix be exploited".
@@ -145,6 +387,25 @@ struct SubstitutionGenerator<TriangularMatrix<LScalar, dimension>,
   static Result Uninitialized(TriangularMatrix<LScalar, dimension> const& m);
 };
 
+template<typename Scalar_>
+struct HessenbergDecompositionGenerator<UnboundedMatrix<Scalar_>> {
+  using Scalar = Scalar_;
+  using Reflector = UnboundedVector<double>;
+  struct Result {
+    UnboundedMatrix<Scalar> H;
+  };
+};
+
+template<typename Scalar_, int dimension>
+struct HessenbergDecompositionGenerator<
+    FixedMatrix<Scalar_, dimension, dimension>> {
+  using Scalar = Scalar_;
+  using Reflector = FixedVector<double, dimension>;
+  struct Result {
+    FixedMatrix<Scalar, dimension, dimension> H;
+  };
+};
+
 template<typename Scalar_>
 struct ClassicalJacobiGenerator<UnboundedMatrix<Scalar_>> {
   using Scalar = Scalar_;
@@ -433,6 +694,42 @@ ForwardSubstitution(LowerTriangularMatrix const& L,
   return x;
 }
 
+template<typename Matrix>
+typename HessenbergDecompositionGenerator<Matrix>::Result
+HessenbergDecomposition(Matrix const& A) {
+  using G = HessenbergDecompositionGenerator<Matrix>;
+  using Scalar = typename G::Scalar;
+  typename HessenbergDecompositionGenerator<Matrix>::Result result{.H = A};
+  auto& H = result.H;
+  int const n = A.rows();
+
+  // [GV13], Algorithm 7.4.2.
+  for (int k = 0; k < n - 2; ++k) {
+    auto const P = ComputeHouseholderReflection(
+        ColumnView<Scalar, Matrix>{.matrix = H,
+                                   .first_row = k + 1,
+                                   .last_row = n - 1,
+                                   .column = k});
+    {
+      auto block = BlockView<Scalar, Matrix>{.matrix = H,
+                                             .first_row = k + 1,
+                                             .last_row = n - 1,
+                                             .first_column = k,
+                                             .last_column = n - 1};
+      Premultiply(P, block);
+    }
+    {
+      auto block = BlockView<Scalar, Matrix>{.matrix = H,
+                                             .first_row = 0,
+                                             .last_row = n - 1,
+                                             .first_column = k + 1,
+                                             .last_column = n - 1};
+      PostMultiply(block, P);
+    }
+  }
+  return result;
+}
+
 template<typename Matrix>
 typename ClassicalJacobiGenerator<Matrix>::Result
 ClassicalJacobi(Matrix const& A,  int max_iterations, double ε) {
diff --git a/numerics/matrix_computations_test.cpp b/numerics/matrix_computations_test.cpp
index 8ecc372523..f2dbbaf0b8 100644
--- a/numerics/matrix_computations_test.cpp
+++ b/numerics/matrix_computations_test.cpp
@@ -8,6 +8,7 @@
 #include "numerics/unbounded_arrays.hpp"
 #include "quantities/elementary_functions.hpp"
 #include "testing_utilities/almost_equals.hpp"
+#include "testing_utilities/vanishes_before.hpp"
 
 namespace principia {
 namespace numerics {
@@ -17,6 +18,7 @@ using namespace principia::numerics::_matrix_computations;
 using namespace principia::numerics::_unbounded_arrays;
 using namespace principia::quantities::_elementary_functions;
 using namespace principia::testing_utilities::_almost_equals;
+using namespace principia::testing_utilities::_vanishes_before;
 
 template<typename T>
 class MatrixComputationsTest : public ::testing::Test {
@@ -102,6 +104,39 @@ TYPED_TEST(MatrixComputationsTest, ForwardSubstitution) {
   EXPECT_THAT(x4_actual, AlmostEquals(x4_expected, 0));
 }
 
+TYPED_TEST(MatrixComputationsTest, HessenbergForm) {
+  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});
+  Matrix const h4_expected({  1,
+                              1.4852968963237645012,
+                             -0.63177627362517020518,
+                              5.1375822981102002423,
+
+                             11.445523142259597039,
+                              7.7328244274809160305,
+                             10.840153034634636509,
+                             -4.4908563773743662120,
+
+                              0,
+                             10.020455740333542876,
+                             -3.9158079318984636532,
+                              0.14244486136993501485,
+
+                              0,
+                              0,
+                              2.4140754086886336638,
+                             -2.8170164955824523773});
+  auto h4_actual = HessenbergDecomposition(m4).H;
+  // This component should really use VanishesBefore, but we don't have a good
+  // way to do that.
+  EXPECT_THAT(h4_actual(3, 1), VanishesBefore(1, 24));
+  h4_actual(3, 1) = 0;
+  EXPECT_THAT(h4_actual, AlmostEquals(h4_expected, 14));
+}
+
 TYPED_TEST(MatrixComputationsTest, ClassicalJacobi) {
   using Vector = typename std::tuple_element<0, TypeParam>::type;
   using Matrix = typename std::tuple_element<3, TypeParam>::type;
diff --git a/numerics/unbounded_arrays.hpp b/numerics/unbounded_arrays.hpp
index e623e842c3..bf67267edf 100644
--- a/numerics/unbounded_arrays.hpp
+++ b/numerics/unbounded_arrays.hpp
@@ -14,6 +14,8 @@ namespace numerics {
 namespace _unbounded_arrays {
 namespace internal {
 
+// TODO(phl): This should support the same operations as fixed_arrays.hpp.
+
 using namespace principia::base::_tags;
 using namespace principia::numerics::_transposed_view;
 using namespace principia::quantities::_named_quantities;
@@ -50,6 +52,7 @@ class UnboundedVector final {
   void EraseToEnd(int begin_index);
 
   Scalar Norm() const;
+  Square<Scalar> Norm²() const;
 
   int size() const;
 
@@ -196,25 +199,110 @@ class UnboundedUpperTriangularMatrix final {
   friend class Row;
 };
 
-template<typename ScalarLeft, typename ScalarRight>
-UnboundedVector<Quotient<ScalarLeft, ScalarRight>> operator/(
-    UnboundedVector<ScalarLeft> const& left,
-    ScalarRight const& right);
-
-template<typename ScalarLeft, typename ScalarRight>
-Product<ScalarLeft, ScalarRight> operator*(
-    TransposedView<UnboundedVector<ScalarLeft>> const& left,
-    UnboundedVector<ScalarRight> const& right);
-
-template<typename ScalarLeft, typename ScalarRight>
-UnboundedMatrix<Product<ScalarLeft, ScalarRight>> operator*(
-    UnboundedMatrix<ScalarLeft> const& left,
-    UnboundedMatrix<ScalarRight> const& right);
-
-template<typename ScalarLeft, typename ScalarRight>
-UnboundedVector<Product<ScalarLeft, ScalarRight>> operator*(
-    UnboundedMatrix<ScalarLeft> const& left,
-    UnboundedVector<ScalarRight> const& right);
+template<typename Scalar>
+UnboundedVector<double> Normalize(UnboundedVector<Scalar> const& vector);
+
+// Additive groups.
+
+template<typename Scalar>
+constexpr UnboundedVector<Scalar>& operator+=(
+    UnboundedVector<Scalar>& left,
+    UnboundedVector<Scalar> const& right);
+
+template<typename Scalar>
+constexpr UnboundedMatrix<Scalar>& operator+=(
+    UnboundedMatrix<Scalar>& left,
+    UnboundedMatrix<Scalar> const& right);
+
+template<typename Scalar>
+constexpr UnboundedVector<Scalar>& operator-=(
+    UnboundedVector<Scalar>& left,
+    UnboundedVector<Scalar> const& right);
+
+template<typename Scalar>
+constexpr UnboundedMatrix<Scalar>& operator-=(
+    UnboundedMatrix<Scalar>& left,
+    UnboundedMatrix<Scalar> const& right);
+
+// Vector spaces.
+
+template<typename LScalar, typename RScalar>
+UnboundedVector<Product<LScalar, RScalar>> operator*(
+    LScalar const& left,
+    UnboundedVector<RScalar> const& right);
+
+template<typename LScalar, typename RScalar>
+UnboundedVector<Product<LScalar, RScalar>> operator*(
+    UnboundedVector<LScalar> const& left,
+    RScalar const& right);
+
+template<typename LScalar, typename RScalar>
+UnboundedMatrix<Product<LScalar, RScalar>>
+operator*(LScalar const& left,
+          UnboundedMatrix<RScalar> const& right);
+
+template<typename LScalar, typename RScalar>
+UnboundedMatrix<Product<LScalar, RScalar>>
+operator*(UnboundedMatrix<LScalar> const& left,
+          RScalar const& right);
+
+template<typename LScalar, typename RScalar>
+UnboundedVector<Quotient<LScalar, RScalar>> operator/(
+    UnboundedVector<LScalar> const& left,
+    RScalar const& right);
+
+template<typename LScalar, typename RScalar>
+constexpr UnboundedMatrix<Quotient<LScalar, RScalar>>
+operator/(UnboundedMatrix<LScalar> const& left,
+          RScalar const& right);
+
+template<typename Scalar>
+constexpr UnboundedVector<Scalar>& operator*=(
+    UnboundedVector<Scalar>& left,
+    double right);
+
+template<typename Scalar>
+constexpr UnboundedMatrix<Scalar>& operator*=(
+    UnboundedMatrix<Scalar>& left,
+    double right);
+
+template<typename Scalar>
+constexpr UnboundedVector<Scalar>& operator/=(
+    UnboundedVector<Scalar>& left,
+    double right);
+
+template<typename Scalar>
+constexpr UnboundedMatrix<Scalar>& operator/=(
+    UnboundedMatrix<Scalar>& left,
+    double right);
+
+// Hilbert space and algebra.
+
+template<typename LScalar, typename RScalar>
+Product<LScalar, RScalar> operator*(
+    TransposedView<UnboundedVector<LScalar>> const& left,
+    UnboundedVector<RScalar> const& right);
+
+template<typename LScalar, typename RScalar>
+UnboundedMatrix<Product<LScalar, RScalar>> operator*(
+    UnboundedVector<LScalar> const& left,
+    TransposedView<UnboundedVector<RScalar>> const& right);
+
+template<typename LScalar, typename RScalar>
+UnboundedMatrix<Product<LScalar, RScalar>> operator*(
+    UnboundedMatrix<LScalar> const& left,
+    UnboundedMatrix<RScalar> const& right);
+
+template<typename LScalar, typename RScalar>
+UnboundedVector<Product<LScalar, RScalar>> operator*(
+    UnboundedMatrix<LScalar> const& left,
+    UnboundedVector<RScalar> const& right);
+
+// Use this operator to multiply a row vector with a matrix.
+template<typename LScalar, typename RScalar>
+UnboundedVector<Product<LScalar, RScalar>> operator*(
+    TransposedView<UnboundedMatrix<LScalar>> const& left,
+    UnboundedVector<RScalar> const& right);
 
 template<typename Scalar>
 std::ostream& operator<<(std::ostream& out,
diff --git a/numerics/unbounded_arrays_body.hpp b/numerics/unbounded_arrays_body.hpp
index ee0c9c61d3..0e0f9f4943 100644
--- a/numerics/unbounded_arrays_body.hpp
+++ b/numerics/unbounded_arrays_body.hpp
@@ -63,11 +63,16 @@ void UnboundedVector<Scalar>::EraseToEnd(int const begin_index) {
 
 template<typename Scalar>
 Scalar UnboundedVector<Scalar>::Norm() const {
+  return Sqrt(Norm²());
+}
+
+template<typename Scalar>
+Square<Scalar> UnboundedVector<Scalar>::Norm²() const {
   Square<Scalar> norm²{};
   for (auto const c : data_) {
     norm² += c * c;
   }
-  return Sqrt(norm²);
+  return norm²;
 }
 
 template<typename Scalar>
@@ -468,36 +473,177 @@ UnboundedUpperTriangularMatrix<Scalar>::Transpose(
   return result;
 }
 
-template<typename ScalarLeft, typename ScalarRight>
-UnboundedVector<Quotient<ScalarLeft, ScalarRight>> operator/(
-    UnboundedVector<ScalarLeft> const& left,
-    ScalarRight const& right) {
-  UnboundedVector<Quotient<ScalarLeft, ScalarRight>> result(left.size(),
-                                                            uninitialized);
+template<typename Scalar>
+UnboundedVector<double> Normalize(UnboundedVector<Scalar> const& vector) {
+  return vector / vector.Norm();
+}
+
+template<typename Scalar>
+constexpr UnboundedVector<Scalar>& operator+=(
+    UnboundedVector<Scalar>& left,
+    UnboundedVector<Scalar> const& right) {
+  return left = left + right;
+}
+
+template<typename Scalar>
+constexpr UnboundedMatrix<Scalar>& operator+=(
+    UnboundedMatrix<Scalar>& left,
+    UnboundedMatrix<Scalar> const& right) {
+  return left = left + right;
+}
+
+template<typename Scalar>
+constexpr UnboundedVector<Scalar>& operator-=(
+    UnboundedVector<Scalar>& left,
+    UnboundedVector<Scalar> const& right) {
+  return left = left - right;
+}
+
+template<typename Scalar>
+constexpr UnboundedMatrix<Scalar>& operator-=(
+    UnboundedMatrix<Scalar>& left,
+    UnboundedMatrix<Scalar> const& right) {
+  return left = left - right;
+}
+
+template<typename LScalar, typename RScalar>
+UnboundedVector<Product<LScalar, RScalar>> operator*(
+    LScalar const& left,
+    UnboundedVector<RScalar> const& right) {
+  UnboundedVector<Product<LScalar, RScalar>> result(right.size(),
+                                                    uninitialized);
+  for (int i = 0; i < right.size(); ++i) {
+    result[i] = left * right[i];
+  }
+  return result;
+}
+
+template<typename LScalar, typename RScalar>
+UnboundedVector<Product<LScalar, RScalar>> operator*(
+    UnboundedVector<LScalar> const& left,
+    RScalar const& right) {
+  UnboundedVector<Product<LScalar, RScalar>> result(left.size(),
+                                                    uninitialized);
+  for (int i = 0; i < left.size(); ++i) {
+    result[i] = left[i] * right;
+  }
+  return result;
+}
+
+template<typename LScalar, typename RScalar>
+UnboundedMatrix<Product<LScalar, RScalar>> operator*(
+    LScalar const& left,
+    UnboundedMatrix<RScalar> const& right) {
+  UnboundedMatrix<Product<LScalar, RScalar>> result(right.rows(),
+                                                    right.columns(),
+                                                    uninitialized);
+  for (int i = 0; i < right.rows(); ++i) {
+    for (int j = 0; j < right.columns(); ++j) {
+      result(i, j) = left * right(i, j);
+    }
+  }
+  return result;
+}
+
+template<typename LScalar, typename RScalar>
+UnboundedMatrix<Product<LScalar, RScalar>> operator*(
+    UnboundedMatrix<LScalar> const& left,
+    RScalar const& right) {
+  UnboundedMatrix<Product<LScalar, RScalar>> result(left.rows(),
+                                                    left.columns(),
+                                                    uninitialized);
+  for (int i = 0; i < left.rows(); ++i) {
+    for (int j = 0; j < left.columns(); ++j) {
+      result(i, j) = left(i, j) * right;
+    }
+  }
+  return result;
+}
+
+template<typename LScalar, typename RScalar>
+UnboundedVector<Quotient<LScalar, RScalar>> operator/(
+    UnboundedVector<LScalar> const& left,
+    RScalar const& right) {
+  UnboundedVector<Quotient<LScalar, RScalar>> result(left.size(),
+                                                     uninitialized);
   for (int i = 0; i < left.size(); ++i) {
     result[i] = left[i] / right;
   }
   return result;
 }
 
-template<typename ScalarLeft, typename ScalarRight>
-Product<ScalarLeft, ScalarRight> operator*(
-    TransposedView<UnboundedVector<ScalarLeft>> const& left,
-    UnboundedVector<ScalarRight> const& right) {
+template<typename LScalar, typename RScalar>
+UnboundedMatrix<Quotient<LScalar, RScalar>> operator/(
+    UnboundedMatrix<LScalar> const& left,
+    RScalar const right) {
+  UnboundedMatrix<Quotient<LScalar, RScalar>> result(left.rows(),
+                                                     left.columns(),
+                                                     uninitialized);
+  for (int i = 0; i < left.rows(); ++i) {
+    for (int j = 0; j < left.columns(); ++j) {
+      result(i, j) = left(i, j) / right;
+    }
+  }
+  return result;
+}
+
+template<typename Scalar>
+constexpr UnboundedVector<Scalar>& operator*=(UnboundedVector<Scalar>& left,
+                                              double const right) {
+  return left = left * right;
+}
+
+template<typename Scalar>
+constexpr UnboundedMatrix<Scalar>& operator*=(UnboundedMatrix<Scalar>& left,
+                                              double const right) {
+  return left = left * right;
+}
+
+template<typename Scalar>
+constexpr UnboundedVector<Scalar>& operator/=(UnboundedVector<Scalar>& left,
+                                              double const right) {
+  return left = left / right;
+}
+
+template<typename Scalar>
+constexpr UnboundedMatrix<Scalar>& operator/=(UnboundedMatrix<Scalar>& left,
+                                              double const right) {
+  return left = left / right;
+}
+
+template<typename LScalar, typename RScalar>
+Product<LScalar, RScalar> operator*(
+    TransposedView<UnboundedVector<LScalar>> const& left,
+    UnboundedVector<RScalar> const& right) {
   CHECK_EQ(left.transpose.size(), right.size());
-  Product<ScalarLeft, ScalarRight> result{};
+  Product<LScalar, RScalar> result{};
   for (int i = 0; i < left.transpose.size(); ++i) {
     result += left.transpose[i] * right[i];
   }
   return result;
 }
 
-template<typename ScalarLeft, typename ScalarRight>
-UnboundedMatrix<Product<ScalarLeft, ScalarRight>> operator*(
-    UnboundedMatrix<ScalarLeft> const& left,
-    UnboundedMatrix<ScalarRight> const& right) {
+template<typename LScalar, typename RScalar>
+UnboundedMatrix<Product<LScalar, RScalar>> operator*(
+    UnboundedVector<LScalar> const& left,
+    TransposedView<UnboundedVector<RScalar>> const& right) {
+  UnboundedMatrix<Product<LScalar, RScalar>> result(left.size(),
+                                                    right.transpose.size(),
+                                                    uninitialized);
+  for (int i = 0; i < result.rows(); ++i) {
+    for (int j = 0; j < result.columns(); ++j) {
+      result(i, j) = left[i] * right.transpose[j];
+    }
+  }
+  return result;
+}
+
+template<typename LScalar, typename RScalar>
+UnboundedMatrix<Product<LScalar, RScalar>> operator*(
+    UnboundedMatrix<LScalar> const& left,
+    UnboundedMatrix<RScalar> const& right) {
   CHECK_EQ(left.columns(), right.rows());
-  UnboundedMatrix<Product<ScalarLeft, ScalarRight>> result(left.rows(),
+  UnboundedMatrix<Product<LScalar, RScalar>> result(left.rows(),
                                                            right.columns());
   for (int i = 0; i < left.rows(); ++i) {
     for (int j = 0; j < right.columns(); ++j) {
@@ -509,12 +655,12 @@ UnboundedMatrix<Product<ScalarLeft, ScalarRight>> operator*(
   return result;
 }
 
-template<typename ScalarLeft, typename ScalarRight>
-UnboundedVector<Product<ScalarLeft, ScalarRight>> operator*(
-    UnboundedMatrix<ScalarLeft> const& left,
-    UnboundedVector<ScalarRight> const& right) {
+template<typename LScalar, typename RScalar>
+UnboundedVector<Product<LScalar, RScalar>> operator*(
+    UnboundedMatrix<LScalar> const& left,
+    UnboundedVector<RScalar> const& right) {
   CHECK_EQ(left.columns(), right.size());
-  UnboundedVector<Product<ScalarLeft, ScalarRight>> result(left.rows());
+  UnboundedVector<Product<LScalar, RScalar>> result(left.rows());
   for (int i = 0; i < left.rows(); ++i) {
     auto& result_i = result[i];
     for (int j = 0; j < left.columns(); ++j) {
@@ -524,6 +670,21 @@ UnboundedVector<Product<ScalarLeft, ScalarRight>> operator*(
   return result;
 }
 
+template<typename LScalar, typename RScalar>
+UnboundedVector<Product<LScalar, RScalar>> operator*(
+    TransposedView<UnboundedMatrix<LScalar>> const& left,
+    UnboundedVector<RScalar> const& right) {
+  CHECK_EQ(left.rows(), right.size());
+  UnboundedVector<Product<LScalar, RScalar>> result(left.columns());
+  for (int j = 0; j < left.columns(); ++j) {
+    auto& result_j = result[j];
+    for (int i = 0; i < left.rows(); ++i) {
+      result_j += left(i, j) * right[i];
+    }
+  }
+  return result;
+}
+
 template<typename Scalar>
 std::ostream& operator<<(std::ostream& out,
                          UnboundedVector<Scalar> const& vector) {
diff --git a/numerics/unbounded_arrays_test.cpp b/numerics/unbounded_arrays_test.cpp
index 4f8350afd3..6811ed6c1d 100644
--- a/numerics/unbounded_arrays_test.cpp
+++ b/numerics/unbounded_arrays_test.cpp
@@ -108,6 +108,10 @@ TEST_F(UnboundedArraysTest, MultiplicationDivision) {
                  661611, 1070697, 1732308, 2803005}), m4_ * m4_);
 }
 
+TEST_F(UnboundedArraysTest, Algebra) {
+  EXPECT_EQ(3270, v3_.Transpose() * v3_);
+}
+
 TEST_F(UnboundedArraysTest, VectorIndexing) {
   EXPECT_EQ(31, v3_[1]);
   v3_[2] = -666;
diff --git a/testing_utilities/almost_equals_body.hpp b/testing_utilities/almost_equals_body.hpp
index 96d2f3170f..301a61b6e3 100644
--- a/testing_utilities/almost_equals_body.hpp
+++ b/testing_utilities/almost_equals_body.hpp
@@ -166,7 +166,7 @@ bool AlmostEqualsMatcher<T>::MatchAndExplain(
   int max_j = -1;
   for (int i = 0; i < 3; ++i) {
     for (int j = 0; j < 3; ++j) {
-      int const distance = NormalizedNaNULPDistance(
+      std::int64_t const distance = NormalizedNaNULPDistance(
           DoubleValue(actual(i, j)), DoubleValue(expected_(i, j)));
       if (distance > max_distance) {
         max_distance = distance;
@@ -338,8 +338,8 @@ bool AlmostEqualsMatcher<T>::MatchAndExplain(
   std::int64_t max_distance = -1;
   int max_i = -1;
   for (int i = 0; i < expected_.size(); ++i) {
-    int const distance = NormalizedNaNULPDistance(DoubleValue(actual[i]),
-                                                  DoubleValue(expected_[i]));
+    std::int64_t const distance = NormalizedNaNULPDistance(
+        DoubleValue(actual[i]), DoubleValue(expected_[i]));
     if (distance > max_distance) {
       max_distance = distance;
       max_i = i;
@@ -368,7 +368,7 @@ bool AlmostEqualsMatcher<T>::MatchAndExplain(
   int max_j = -1;
   for (int i = 0; i < expected_.rows(); ++i) {
     for (int j = 0; j < expected_.columns(); ++j) {
-      int const distance = NormalizedNaNULPDistance(
+      std::int64_t const distance = NormalizedNaNULPDistance(
           DoubleValue(actual(i, j)), DoubleValue(expected_(i, j)));
       if (distance > max_distance) {
         max_distance = distance;
@@ -400,7 +400,7 @@ bool AlmostEqualsMatcher<T>::MatchAndExplain(
   int max_j = -1;
   for (int i = 0; i < expected_.rows(); ++i) {
     for (int j = 0; j <= i; ++j) {
-      int const distance = NormalizedNaNULPDistance(
+      std::int64_t const distance = NormalizedNaNULPDistance(
           DoubleValue(actual(i, j)), DoubleValue(expected_(i, j)));
       if (distance > max_distance) {
         max_distance = distance;
@@ -432,7 +432,7 @@ bool AlmostEqualsMatcher<T>::MatchAndExplain(
   int max_j = -1;
   for (int i = 0; i < expected_.columns(); ++i) {
     for (int j = i; j < expected_.columns(); ++j) {
-      int const distance = NormalizedNaNULPDistance(
+      std::int64_t const distance = NormalizedNaNULPDistance(
           DoubleValue(actual(i, j)), DoubleValue(expected_(i, j)));
       if (distance > max_distance) {
         max_distance = distance;
@@ -466,8 +466,8 @@ bool AlmostEqualsMatcher<T>::MatchAndExplain(
   std::int64_t max_distance = -1;
   int max_i = -1;
   for (int i = 0; i < expected_.size(); ++i) {
-    int const distance = NormalizedNaNULPDistance(DoubleValue(actual[i]),
-                                                  DoubleValue(expected_[i]));
+    std::int64_t const distance = NormalizedNaNULPDistance(
+        DoubleValue(actual[i]), DoubleValue(expected_[i]));
     if (distance > max_distance) {
       max_distance = distance;
       max_i = i;
@@ -500,7 +500,7 @@ bool AlmostEqualsMatcher<T>::MatchAndExplain(
   int max_j = -1;
   for (int i = 0; i < expected_.rows(); ++i) {
     for (int j = 0; j < expected_.columns(); ++j) {
-      int const distance = NormalizedNaNULPDistance(
+      std::int64_t const distance = NormalizedNaNULPDistance(
           DoubleValue(actual(i, j)), DoubleValue(expected_(i, j)));
       if (distance > max_distance) {
         max_distance = distance;
@@ -536,7 +536,7 @@ bool AlmostEqualsMatcher<T>::MatchAndExplain(
   int max_j = -1;
   for (int i = 0; i < expected_.rows(); ++i) {
     for (int j = 0; j <= i; ++j) {
-      int const distance = NormalizedNaNULPDistance(
+      std::int64_t const distance = NormalizedNaNULPDistance(
           DoubleValue(actual(i, j)), DoubleValue(expected_(i, j)));
       if (distance > max_distance) {
         max_distance = distance;
@@ -573,7 +573,7 @@ bool AlmostEqualsMatcher<T>::MatchAndExplain(
   int max_j = -1;
   for (int i = 0; i < expected_.columns(); ++i) {
     for (int j = i; j < expected_.columns(); ++j) {
-      int const distance = NormalizedNaNULPDistance(
+      std::int64_t const distance = NormalizedNaNULPDistance(
           DoubleValue(actual(i, j)), DoubleValue(expected_(i, j)));
       if (distance > max_distance) {
         max_distance = distance;