Skip to content

Commit

Permalink
Remove exp(), pow(), and sqrt() methods
Browse files Browse the repository at this point in the history
  • Loading branch information
@cm-jones
cm-jones committed Jul 5, 2024
1 parent 4ecbe5e commit 3d82433
Show file tree
Hide file tree
Showing 4 changed files with 49 additions and 194 deletions.
98 changes: 49 additions & 49 deletions benchmarks/benchmark_matrix.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@ using namespace cramer;
// Helper function to create a random matrix
template <typename T>
Matrix<T> createRandomMatrix(size_t rows, size_t cols) {
static std::random_device rd;
static std::random_device rand;
static std::mt19937 gen(rd());
static std::uniform_real_distribution<> dis(-1.0, 1.0);

Expand All @@ -28,20 +28,20 @@ Matrix<T> createRandomMatrix(size_t rows, size_t cols) {
// Benchmark Matrix Constructor
static void BM_MatrixConstructor(benchmark::State& state) {
const size_t size = state.range(0);
for (auto _ : state) {
Matrix<double> m(size, size);
benchmark::DoNotOptimize(m);
for (auto iter : state) {
Matrix<double> matrix(size, size);
benchmark::DoNotOptimize(matrix);
}
}
BENCHMARK(BM_MatrixConstructor)->Range(8, 1024);

// Benchmark Matrix Addition
static void BM_MatrixAddition(benchmark::State& state) {
const size_t size = state.range(0);
auto m1 = createRandomMatrix<double>(size, size);
auto m2 = createRandomMatrix<double>(size, size);
for (auto _ : state) {
auto result = m1 + m2;
auto matrix1 = createRandomMatrix<double>(size, size);
auto matrix2 = createRandomMatrix<double>(size, size);
for (auto iter : state) {
auto result = matrix1 + matrix2;
benchmark::DoNotOptimize(result);
}
}
Expand All @@ -50,10 +50,10 @@ BENCHMARK(BM_MatrixAddition)->Range(8, 1024);
// Benchmark Matrix Subtraction
static void BM_MatrixSubtraction(benchmark::State& state) {
const size_t size = state.range(0);
auto m1 = createRandomMatrix<double>(size, size);
auto m2 = createRandomMatrix<double>(size, size);
for (auto _ : state) {
auto result = m1 - m2;
auto matrix1 = createRandomMatrix<double>(size, size);
auto matrix2 = createRandomMatrix<double>(size, size);
for (auto iter : state) {
auto result = matrix1 - matrix2;
benchmark::DoNotOptimize(result);
}
}
Expand All @@ -62,10 +62,10 @@ BENCHMARK(BM_MatrixSubtraction)->Range(8, 1024);
// Benchmark Matrix Multiplication
static void BM_MatrixMultiplication(benchmark::State& state) {
const size_t size = state.range(0);
auto m1 = createRandomMatrix<double>(size, size);
auto m2 = createRandomMatrix<double>(size, size);
for (auto _ : state) {
auto result = m1 * m2;
auto matrix1 = createRandomMatrix<double>(size, size);
auto matrix2 = createRandomMatrix<double>(size, size);
for (auto iter : state) {
auto result = matrix1 * matrix2;
benchmark::DoNotOptimize(result);
}
}
Expand All @@ -74,10 +74,10 @@ BENCHMARK(BM_MatrixMultiplication)->Range(8, 256);
// Benchmark Matrix Scalar Multiplication
static void BM_MatrixScalarMultiplication(benchmark::State& state) {
const size_t size = state.range(0);
auto m = createRandomMatrix<double>(size, size);
auto matrix = createRandomMatrix<double>(size, size);
double scalar = 2.0;
for (auto _ : state) {
auto result = m * scalar;
for (auto iter : state) {
auto result = matrix * 2.0;
benchmark::DoNotOptimize(result);
}
}
Expand All @@ -86,9 +86,9 @@ BENCHMARK(BM_MatrixScalarMultiplication)->Range(8, 1024);
// Benchmark Matrix Transpose
static void BM_MatrixTranspose(benchmark::State& state) {
const size_t size = state.range(0);
auto m = createRandomMatrix<double>(size, size);
for (auto _ : state) {
auto result = m.transpose();
auto matrix = createRandomMatrix<double>(size, size);
for (auto iter : state) {
auto result = matrix.transpose();
benchmark::DoNotOptimize(result);
}
}
Expand All @@ -97,9 +97,9 @@ BENCHMARK(BM_MatrixTranspose)->Range(8, 1024);
// Benchmark Matrix Determinant
static void BM_MatrixDeterminant(benchmark::State& state) {
const size_t size = state.range(0);
auto m = createRandomMatrix<double>(size, size);
for (auto _ : state) {
auto result = m.det();
auto matrix = createRandomMatrix<double>(size, size);
for (auto iter : state) {
auto result = matrix.det();
benchmark::DoNotOptimize(result);
}
}
Expand All @@ -108,9 +108,9 @@ BENCHMARK(BM_MatrixDeterminant)->Range(2, 128);
// Benchmark Matrix Inverse
static void BM_MatrixInverse(benchmark::State& state) {
const size_t size = state.range(0);
auto m = createRandomMatrix<double>(size, size);
for (auto _ : state) {
auto result = m.inverse();
auto matrix = createRandomMatrix<double>(size, size);
for (auto iter : state) {
auto result = matrix.inverse();
benchmark::DoNotOptimize(result);
}
}
Expand All @@ -119,9 +119,9 @@ BENCHMARK(BM_MatrixInverse)->Range(2, 128);
// Benchmark Matrix LU Decomposition
static void BM_MatrixLU(benchmark::State& state) {
const size_t size = state.range(0);
auto m = createRandomMatrix<double>(size, size);
for (auto _ : state) {
auto [L, U] = m.lu();
auto matrix = createRandomMatrix<double>(size, size);
for (auto iter : state) {
auto [L, U] = matrix.lu();
benchmark::DoNotOptimize(L);
benchmark::DoNotOptimize(U);
}
Expand All @@ -131,9 +131,9 @@ BENCHMARK(BM_MatrixLU)->Range(2, 128);
// Benchmark Matrix QR Decomposition
static void BM_MatrixQR(benchmark::State& state) {
const size_t size = state.range(0);
auto m = createRandomMatrix<double>(size, size);
for (auto _ : state) {
auto [Q, R] = m.qr();
auto matrix = createRandomMatrix<double>(size, size);
for (auto iter : state) {
auto [Q, R] = matrix.qr();
benchmark::DoNotOptimize(Q);
benchmark::DoNotOptimize(R);
}
Expand All @@ -143,9 +143,9 @@ BENCHMARK(BM_MatrixQR)->Range(2, 128);
// Benchmark Matrix SVD
static void BM_MatrixSVD(benchmark::State& state) {
const size_t size = state.range(0);
auto m = createRandomMatrix<double>(size, size);
for (auto _ : state) {
auto [U, S, V] = m.svd();
auto matrix = createRandomMatrix<double>(size, size);
for (auto iter : state) {
auto [U, S, V] = matrix.svd();
benchmark::DoNotOptimize(U);
benchmark::DoNotOptimize(S);
benchmark::DoNotOptimize(V);
Expand All @@ -156,9 +156,9 @@ BENCHMARK(BM_MatrixSVD)->Range(2, 64);
// Benchmark Matrix Eigenvalues
static void BM_MatrixEigenvalues(benchmark::State& state) {
const size_t size = state.range(0);
auto m = createRandomMatrix<double>(size, size);
for (auto _ : state) {
auto result = m.eigenvalues();
auto matrix = createRandomMatrix<double>(size, size);
for (auto iter : state) {
auto result = matrix.eigenvalues();
benchmark::DoNotOptimize(result);
}
}
Expand All @@ -167,9 +167,9 @@ BENCHMARK(BM_MatrixEigenvalues)->Range(2, 64);
// Benchmark Matrix Eigenvectors
static void BM_MatrixEigenvectors(benchmark::State& state) {
const size_t size = state.range(0);
auto m = createRandomMatrix<double>(size, size);
for (auto _ : state) {
auto result = m.eigenvectors();
auto matrix = createRandomMatrix<double>(size, size);
for (auto iter : state) {
auto result = matrix.eigenvectors();
benchmark::DoNotOptimize(result);
}
}
Expand All @@ -178,9 +178,9 @@ BENCHMARK(BM_MatrixEigenvectors)->Range(2, 64);
// Benchmark Matrix Rank
static void BM_MatrixRank(benchmark::State& state) {
const size_t size = state.range(0);
auto m = createRandomMatrix<double>(size, size);
for (auto _ : state) {
auto result = m.rank();
auto matrix = createRandomMatrix<double>(size, size);
for (auto iter : state) {
auto result = matrix.rank();
benchmark::DoNotOptimize(result);
}
}
Expand All @@ -189,9 +189,9 @@ BENCHMARK(BM_MatrixRank)->Range(2, 128);
// Benchmark Matrix Trace
static void BM_MatrixTrace(benchmark::State& state) {
const size_t size = state.range(0);
auto m = createRandomMatrix<double>(size, size);
for (auto _ : state) {
auto result = m.trace();
auto matrix = createRandomMatrix<double>(size, size);
for (auto iter : state) {
auto result = matrix.trace();
benchmark::DoNotOptimize(result);
}
}
Expand Down
19 changes: 0 additions & 19 deletions include/matrix.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -264,25 +264,6 @@ class Matrix {
*/
Matrix<T> conjugate() const;

/**
* @brief Calculates the matrix exponential.
* @return The exponential of the matrix.
*/
Matrix<T> exp() const;

/**
* @brief Calculates the matrix power.
* @param n The power to raise the matrix to.
* @return The matrix raised to the power n.
*/
Matrix<T> pow(int n) const;

/**
* @brief Calculates the square root of the matrix.
* @return The square root of the matrix.
*/
Matrix<T> sqrt() const;

/**
* @brief Performs LU decomposition of the matrix.
* @return A pair of matrices (L, U) where L is lower triangular and U is
Expand Down
99 changes: 0 additions & 99 deletions src/matrix.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -621,63 +621,6 @@ Matrix<T> Matrix<T>::exp() const {
return F;
}

template <typename T>
Matrix<T> Matrix<T>::pow(int n) const {
if (!is_square()) {
throw std::invalid_argument(
"Matrix must be square to calculate its integer power.");
}

if (n == 0) {
return identity(rows);
}

Matrix<T> result = identity(rows);
Matrix<T> base = *this;
bool negative = (n < 0);
n = std::abs(n);

while (n > 0) {
if (n % 2 == 1) {
result = result * base; // Use matrix multiplication instead of *=
}
base = base * base; // Use matrix multiplication instead of *=
n /= 2;
}

return negative ? result.inverse() : result;
}

template <typename T>
Matrix<T> Matrix<T>::sqrt() const {
if (!is_square()) {
throw std::invalid_argument(
"Matrix must be square to calculate its square root.");
}

// This is a placeholder implementation using Denman-Beavers iteration
// For a more robust implementation, consider using Schur decomposition or
// other methods
Matrix<T> X = *this;
Matrix<T> Y = identity(rows);
const int max_iterations = 20;
const T tolerance = std::numeric_limits<T>::epsilon();

for (int i = 0; i < max_iterations; ++i) {
Matrix<T> X_next = (X + Y.inverse()) * static_cast<T>(0.5);
Matrix<T> Y_next = (Y + X.inverse()) * static_cast<T>(0.5);

if (std::abs((X_next - X).max_norm()) < std::abs(tolerance)) {
return X_next;
}

X = X_next;
Y = Y_next;
}

throw std::runtime_error("Square root iteration did not converge.");
}

// Matrix decompositions

template <typename T>
Expand Down Expand Up @@ -717,48 +660,6 @@ std::pair<Matrix<T>, Matrix<T>> Matrix<T>::lu() const {
return std::make_pair(L, U);
}

template <typename T>
Matrix<T> Matrix<T>::exp() const {
if (!is_square()) {
throw std::invalid_argument("Matrix must be square to calculate the exponential.");
}

const int q = 6;
const int p = 1;
T norm = max_norm();
int s = std::max(0, static_cast<int>(std::ceil(std::log2(norm))));

Matrix<T> A = *this * (1.0 / std::pow(2.0, s));
Matrix<T> X = identity(rows) + A;
Matrix<T> cX = identity(rows) - A;
Matrix<T> A2 = A * A;

for (int k = 2; k <= q; ++k) {
T c = 1.0;
for (int j = 1; j <= std::min(k, p); ++j) {
c *= static_cast<T>(k - j + 1) / static_cast<T>(j * (2 * k - j + 1));
}
Matrix<T> Y = c * A2;
X += Y;
if (k % 2 == 0) {
cX += Y;
} else {
cX -= Y;
}
if (k < q) {
A2 *= A;
}
}

Matrix<T> E = X * cX.inverse();

for (int k = 0; k < s; ++k) {
E = E * E;
}

return E;
}

template <typename T>
std::tuple<Matrix<T>, Matrix<T>, Matrix<T>> Matrix<T>::svd() const {
const size_t m = rows;
Expand Down
27 changes: 0 additions & 27 deletions tests/test_matrix.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -334,33 +334,6 @@ TEST_F(MatrixTest, Conjugate) {
EXPECT_EQ(result(1, 1), std::complex<double>(4, -4));
}

TEST_F(MatrixTest, Exp) {
Matrix<double> mat(2, 2, {0, 1, -1, 0});
Matrix<double> result = mat.exp();
EXPECT_NEAR(result(0, 0), std::cos(1), 1e-9);
EXPECT_NEAR(result(0, 1), std::sin(1), 1e-9);
EXPECT_NEAR(result(1, 0), -std::sin(1), 1e-9);
EXPECT_NEAR(result(1, 1), std::cos(1), 1e-9);
}

TEST_F(MatrixTest, Pow) {
Matrix<double> mat(2, 2, {1, 1, 1, 0});
Matrix<double> result = mat.pow(3);
EXPECT_EQ(result(0, 0), 3);
EXPECT_EQ(result(0, 1), 2);
EXPECT_EQ(result(1, 0), 2);
EXPECT_EQ(result(1, 1), 1);
}

TEST_F(MatrixTest, Sqrt) {
Matrix<double> mat(2, 2, {4, 0, 0, 9});
Matrix<double> result = mat.sqrt();
EXPECT_NEAR(result(0, 0), 2, 1e-9);
EXPECT_NEAR(result(0, 1), 0, 1e-9);
EXPECT_NEAR(result(1, 0), 0, 1e-9);
EXPECT_NEAR(result(1, 1), 3, 1e-9);
}

TEST_F(MatrixTest, LU) {
Matrix<double> mat(3, 3, {1, 2, 3, 4, 5, 6, 7, 8, 9});
auto [lower, upper] = mat.lu();
Expand Down

0 comments on commit 3d82433

Please sign in to comment.