lisyarus/psemek
1
0
Fork 1
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Add QR eigenvalue algorithm and tests (non-symmetric tests need some more checking)

Browse source
This commit is contained in:
Nikita Lisitsa 2023-12-28 17:16:12 +03:00
parent 5eef1e13f5
commit 2364c15120
2 changed files with 87 additions and 1 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

21
libs/geom/include/psemek/geom/qr.hpp
View file

@ -9,9 +9,10 @@ namespace psemek::geom
struct qr_result struct qr_result
{ {
matrix<T, N, N> q; // orthogonal matrix<T, N, N> q; // orthogonal
matrix<T, N, N> r; // upper-triangular matrix<T, N, N> r; // upper-triangular (for decomposition) or block upper-triangular (for eigenvalues)
}; };
// returns {q, r} such that m = q * r
template <typename T, std::size_t N> template <typename T, std::size_t N>
qr_result<T, N> qr_decomposition(matrix<T, N, N> const & m) qr_result<T, N> qr_decomposition(matrix<T, N, N> const & m)
{ {
@ -59,4 +60,22 @@ namespace psemek::geom
return {q, r}; return {q, r};
} }
// returns {q, r} such that m = q * r * q^(-1)
// NB: q^(-1) = q^T
template <typename T, std::size_t N>
qr_result<T, N> qr_eigenvalues(matrix<T, N, N> const & m, std::size_t const iterations)
{
matrix<T, N, N> r = m;
matrix<T, N, N> q = matrix<T, N, N>::identity();
for (std::size_t iteration = 0; iteration < iterations; ++iteration)
{
auto qr = qr_decomposition(r);
r = qr.r * qr.q;
q = q * qr.q;
}
return {q, r};
}
} }

67
libs/geom/tests/matrix.cpp
View file

@ -217,3 +217,70 @@ test_case(geom_matrix_qr)
expect_small(result.r[j][i], 1e-6); expect_small(result.r[j][i], 1e-6);
} }
} }
test_case(geom_matrix_qr__eig_symmetric)
{
random::generator rng;
random::normal_distribution<float> normal;
for (int iteration = 0; iteration < 64; ++iteration)
{
matrix<float, 8, 8> m;
for (std::size_t i = 0; i < m.rows(); ++i)
for (std::size_t j = 0; j < m.columns(); ++j)
m[i][j] = normal(rng);
m = transpose(m) * m;
auto result = qr_eigenvalues(m, 1024);
expect_small(frobenius_norm(m * result.q - result.q * result.r), 1e-3);
expect_small(frobenius_norm(matrix<float, 8, 8>::identity() - result.q * transpose(result.q)), 1e-4);
for (std::size_t i = 0; i < m.rows(); ++i)
for (std::size_t j = i + 1; j < m.columns(); ++j)
expect_small(result.r[j][i], 1e-3);
std::string test;
for (std::size_t i = 0; i < m.columns(); ++i)
{
auto v = column(result.q, i);
expect_small(length(m * v - result.r[i][i] * v), 1e-3);
}
}
}
test_case(geom_matrix_qr__eig_general)
{
random::generator rng;
random::normal_distribution<float> normal;
for (int iteration = 0; iteration < 64; ++iteration)
{
matrix<float, 8, 8> m;
for (std::size_t i = 0; i < m.rows(); ++i)
for (std::size_t j = 0; j < m.columns(); ++j)
m[i][j] = normal(rng);
auto result = qr_eigenvalues(m, 4096);
expect_small(frobenius_norm(m * result.q - result.q * result.r), 1e-2);
expect_small(frobenius_norm(matrix<float, 8, 8>::identity() - result.q * transpose(result.q)), 1e-3);
for (std::size_t i = 0; i < m.rows(); ++i)
for (std::size_t j = i + 2; j < m.columns(); ++j)
expect_small(result.r[j][i], 1e-2);
for (std::size_t i = 0; i < m.columns(); ++i)
{
// Don't check complex eigenvalues
if (i > 0 && std::abs(result.r[i][i - 1]) > 1e-3f) continue;
if (i + 1 < m.rows() && std::abs(result.r[i + 1][i]) > 1e-3f) continue;
// TODO: properly check eigenvalues
// auto v = column(result.q, i);
// expect_small(length(m * v - result.r[i][i] * v), 1e-3);
}
}
}