Skip to content

Commit

Permalink
Merge pull request #3833 from pleroy/Hessenberg
Browse files Browse the repository at this point in the history 
Hessenberg decomposition of a matrix
  • Loading branch information
@pleroy
pleroy authored Jan 3, 2024
2 parents ff7f1aa + 7b9a381 commit 43b413c
Show file tree
Hide file tree
Showing 11 changed files with 750 additions and 77 deletions.
6 changes: 3 additions & 3 deletions geometry/r3x3_matrix_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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>(
Expand Down
46 changes: 37 additions & 9 deletions numerics/fixed_arrays.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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/(
Expand All @@ -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.

Expand Down Expand Up @@ -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>
Expand Down
64 changes: 52 additions & 12 deletions numerics/fixed_arrays_body.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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_);
}

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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) {
Expand All @@ -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;
Expand All @@ -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);
Expand All @@ -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) {
Expand All @@ -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) {
Expand All @@ -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,
Expand Down Expand Up @@ -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) {
Expand Down
2 changes: 2 additions & 0 deletions numerics/fixed_arrays_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down
20 changes: 19 additions & 1 deletion numerics/matrix_computations.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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.
Expand All @@ -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);
Expand All @@ -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;
Expand Down
Loading

0 comments on commit 43b413c

Please sign in to comment.