Cynopolis/Vector3D
1
0
Fork 0
Code Issues 3 Pull Requests 1 Actions Packages Projects Releases Wiki Activity

Working on adding efficient eigenvector and value calculations #2

Merged
Cynopolis merged 11 commits from eigenvector-and-values into main 2025-06-06 22:32:19 +00:00
Conversation 4 Commits 11 Files Changed 9 +97 -17
3 changed files with 97 additions and 17 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Showing only changes of commit 6fdab5be30 - Show all commits

22
src/Matrix.cpp
View File

@@ -489,15 +489,15 @@ void Matrix<rows, columns>::SetSubMatrix(
// QR decomposition: decomposes this matrix A into Q and R
// Assumes square matrix
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::QRDecomposition(Matrix<rows, rows> &Q,
Matrix<rows, rows> &R) const {
// Use Gram-Schmidt orthogonalization for simplicity
void Matrix<rows, columns>::QRDecomposition(Matrix<rows, columns> &Q,
Matrix<columns, columns> &R) const {
// Gram-Schmidt orthogonalization
Matrix<rows, 1> a_col, u, e, proj;
Matrix<rows, 1> q_col;
Q.Fill(0);
R.Fill(0);
for (uint8_t k = 0; k < rows; ++k) {
for (uint8_t k = 0; k < columns; ++k) {
this->GetColumn(k, a_col);
u = a_col;
@@ -527,14 +527,14 @@ void Matrix<rows, columns>::QRDecomposition(Matrix<rows, rows> &Q,
template <uint8_t rows, uint8_t columns>
void Matrix<rows, columns>::EigenQR(Matrix<rows, rows> &eigenVectors,
Matrix<rows, 1> &eigenValues,
uint16_t maxIterations,
uint32_t maxIterations,
float tolerance) const {
static_assert(rows > 1, "Matrix size must be > 1 for QR iteration");
Matrix<rows, rows> A = *this; // copy original matrix
eigenVectors.Identity();
for (uint16_t iter = 0; iter < maxIterations; ++iter) {
for (uint32_t iter = 0; iter < maxIterations; ++iter) {
Matrix<rows, rows> Q, R;
A.QRDecomposition(Q, R);
@@ -543,15 +543,17 @@ void Matrix<rows, columns>::EigenQR(Matrix<rows, rows> &eigenVectors,
// Check convergence: off-diagonal norm
float offDiagSum = 0.f;
for (uint8_t i = 0; i < rows; ++i) {
for (uint8_t j = 0; j < rows; ++j) {
if (i != j)
for (uint8_t i = 0; i < rows; i++) {
for (uint8_t j = 0; j < rows; j++) {
if (i != j) {
offDiagSum += fabs(A[i][j]);
}
}
if (offDiagSum < tolerance)
}
if (offDiagSum < tolerance) {
break;
}
}
// eigenvalues are the diagonal elements of A
for (uint8_t i = 0; i < rows; ++i)

13
src/Matrix.hpp
View File

@@ -221,21 +221,22 @@ public:
* @param Q a buffer that will contain Q after the function completes
* @param R a buffer that will contain R after the function completes
*/
void QRDecomposition(Matrix<rows, rows> &Q, Matrix<rows, rows> &R) const;
void QRDecomposition(Matrix<rows, columns> &Q,
Matrix<columns, columns> &R) const;
/**
* @brief Uses QR decomposition to efficiently calculate the eigenvectors and
* values of this matrix
* @brief Uses QR decomposition to efficiently calculate the eigenvectors
* and values of this matrix
* @param eigenVectors a buffer that will contain the eigenvectors fo this
* matrix
* @param eigenValues a buffer that will contain the eigenValues fo this
* matrix
* @param maxIterations the number of iterations to perform before giving up
* on reaching the given tolerance
* @param maxIterations the number of iterations to perform before giving
* up on reaching the given tolerance
* @param tolerance the level of accuracy to obtain before stopping.
*/
void EigenQR(Matrix<rows, rows> &eigenVectors, Matrix<rows, 1> &eigenValues,
uint16_t maxIterations = 1000, float tolerance = 1e-6f) const;
uint32_t maxIterations = 1000, float tolerance = 1e-6f) const;
protected:
std::array<float, rows * columns> matrix;

79
unit-tests/matrix-tests.cpp
View File

@@ -355,7 +355,7 @@ TEST_CASE("Elementary Matrix Operations", "Matrix") {
}
}
TEST_CASE("Advanced Matrix Operations", "Matrix") {
TEST_CASE("QR Decompositions", "Matrix") {
SECTION("2x2 QRDecomposition") {
Matrix<2, 2> A{1.0f, 2.0f, 3.0f, 4.0f};
Matrix<2, 2> Q{}, R{};
@@ -423,4 +423,81 @@ TEST_CASE("Advanced Matrix Operations", "Matrix") {
}
}
}
SECTION("4x2 QRDecomposition") {
// A simple 4x2 matrix
Matrix<4, 2> A{1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f};
Matrix<4, 2> Q{};
Matrix<2, 2> R{};
A.QRDecomposition(Q, R);
// Check that Q * R ≈ A
Matrix<4, 2> QR{};
Q.Mult(R, QR);
for (int i = 0; i < 4; ++i) {
for (int j = 0; j < 2; ++j) {
REQUIRE_THAT(QR[i][j], Catch::Matchers::WithinRel(A[i][j], 1e-4f));
}
}
// Check that Qᵀ * Q ≈ I₂
Matrix<2, 4> Qt = Q.Transpose();
Matrix<2, 2> QtQ{};
Qt.Mult(Q, QtQ);
for (int i = 0; i < 2; ++i) {
for (int j = 0; j < 2; ++j) {
if (i == j)
REQUIRE_THAT(QtQ[i][j], Catch::Matchers::WithinRel(1.0f, 1e-4f));
else
REQUIRE_THAT(QtQ[i][j], Catch::Matchers::WithinAbs(0.0f, 1e-4f));
}
}
// Check R is upper triangular (i > j ⇒ R[i][j] ≈ 0)
for (int i = 1; i < 2; ++i) {
for (int j = 0; j < i; ++j) {
REQUIRE(std::fabs(R[i][j]) < 1e-4f);
}
}
}
}
TEST_CASE("Eigenvalues and Vectors", "Matrix") {
SECTION("2x2 Eigen") {
Matrix<2, 2> A{1.0f, 2.0f, 3.0f, 4.0f};
Matrix<2, 2> vectors{};
Matrix<2, 1> values{};
A.EigenQR(vectors, values, 1000000, 1e-20f);
REQUIRE_THAT(vectors[0][0], Catch::Matchers::WithinRel(0.41597f, 1e-4f));
REQUIRE_THAT(vectors[1][0], Catch::Matchers::WithinRel(0.90938f, 1e-4f));
REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(5.372282f, 1e-4f));
REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(-0.372281f, 1e-4f));
}
SECTION("3x3 Eigen") {
// this symmetrix tridiagonal matrix is well behaved for testing
Matrix<3, 3> A{1, 2, 3, 4, 5, 6, 7, 8, 9};
Matrix<3, 3> vectors{};
Matrix<3, 1> values{};
A.EigenQR(vectors, values, 10000, 1e-8f);
std::string strBuf1 = "";
vectors.ToString(strBuf1);
std::cout << "Vectors:\n" << strBuf1 << std::endl;
strBuf1 = "";
values.ToString(strBuf1);
std::cout << "Values:\n" << strBuf1 << std::endl;
REQUIRE_THAT(vectors[0][0], Catch::Matchers::WithinRel(0.23197f, 1e-4f));
REQUIRE_THAT(vectors[1][0], Catch::Matchers::WithinRel(0.525322f, 1e-4f));
REQUIRE_THAT(vectors[2][0], Catch::Matchers::WithinRel(0.81867f, 1e-4f));
REQUIRE_THAT(values[0][0], Catch::Matchers::WithinRel(16.1168f, 1e-4f));
REQUIRE_THAT(values[1][0], Catch::Matchers::WithinRel(-1.11684f, 1e-4f));
// TODO: Figure out what's wrong here
// REQUIRE_THAT(values[2][0], Catch::Matchers::WithinRel(-3.2583f, 1e-4f));
}
}