From 1b0a8a2f702eafee776f64d77a3672c5d7a409f2 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Fri, 29 May 2026 13:43:14 -0500 Subject: [PATCH 1/4] can generate eigen systems --- include/tensorwrapper/generate/generate.hpp | 23 +++ .../generate/generate_eigen_system.hpp | 33 +++++ .../generate/generate_eigenvalues.hpp | 131 ++++++++++++++++++ .../tensorwrapper/generate/generate_utils.hpp | 62 +++++++++ .../generate/random_orthogonal_matrix.hpp | 25 ++++ include/tensorwrapper/tensorwrapper.hpp | 1 + .../generate/generate_eigen_system.cpp | 45 ++++++ .../generate/random_orthogonal_matrix.cpp | 49 +++++++ .../generate/generate_eigen_system.cpp | 127 +++++++++++++++++ .../generate/generate_eigenvalues.cpp | 129 +++++++++++++++++ .../tensorwrapper/generate/generate_utils.cpp | 35 +++++ .../generate/random_orthogonal_matrix.cpp | 58 ++++++++ 12 files changed, 718 insertions(+) create mode 100644 include/tensorwrapper/generate/generate.hpp create mode 100644 include/tensorwrapper/generate/generate_eigen_system.hpp create mode 100644 include/tensorwrapper/generate/generate_eigenvalues.hpp create mode 100644 include/tensorwrapper/generate/generate_utils.hpp create mode 100644 include/tensorwrapper/generate/random_orthogonal_matrix.hpp create mode 100644 src/tensorwrapper/generate/generate_eigen_system.cpp create mode 100644 src/tensorwrapper/generate/random_orthogonal_matrix.cpp create mode 100644 tests/cxx/unit_tests/tensorwrapper/generate/generate_eigen_system.cpp create mode 100644 tests/cxx/unit_tests/tensorwrapper/generate/generate_eigenvalues.cpp create mode 100644 tests/cxx/unit_tests/tensorwrapper/generate/generate_utils.cpp create mode 100644 tests/cxx/unit_tests/tensorwrapper/generate/random_orthogonal_matrix.cpp diff --git a/include/tensorwrapper/generate/generate.hpp b/include/tensorwrapper/generate/generate.hpp new file mode 100644 index 00000000..624e4f4c --- /dev/null +++ b/include/tensorwrapper/generate/generate.hpp @@ -0,0 +1,23 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include +#include +#include + +namespace tensorwrapper::generate {} diff --git a/include/tensorwrapper/generate/generate_eigen_system.hpp b/include/tensorwrapper/generate/generate_eigen_system.hpp new file mode 100644 index 00000000..9e288569 --- /dev/null +++ b/include/tensorwrapper/generate/generate_eigen_system.hpp @@ -0,0 +1,33 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include +#include + +namespace tensorwrapper::generate { + +struct EigenSystem { + std::size_t n = 0; + Tensor eigenvalues; + Tensor eigenvectors; + Tensor matrix; +}; + +EigenSystem generate_eigen_system(const SymmetricMatrixSpec& spec); + +} // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/generate/generate_eigenvalues.hpp b/include/tensorwrapper/generate/generate_eigenvalues.hpp new file mode 100644 index 00000000..9f1ef47d --- /dev/null +++ b/include/tensorwrapper/generate/generate_eigenvalues.hpp @@ -0,0 +1,131 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include +#include +#include +#include + +namespace tensorwrapper::generate { +namespace { +inline auto linear_spacing(std::size_t n, double lambda_min, + double lambda_max) { + const auto dx = (lambda_max - lambda_min) / (n - 1); + std::vector values(n); + for(std::size_t i = 0; i < n; ++i) { + values[i] = lambda_min + static_cast(i) * dx; + } + return values; +} + +inline auto logarithmic_spacing(std::size_t n, double lambda_min, + double lambda_max) { + const double log_min = std::log(lambda_min); + const double log_max = std::log(lambda_max); + const double dlog = (log_max - log_min) / (n - 1); + std::vector values(n); + for(std::size_t i = 0; i < n; ++i) { + values[i] = std::exp(log_min + static_cast(i) * dlog); + } + return values; +} + +inline auto clustered_spacing(std::size_t n, std::size_t n_clusters, + double lambda_min, double lambda_max, + double cluster_width, std::mt19937& gen) { + std::uniform_real_distribution jitter(-cluster_width, + cluster_width); + std::vector values(n); + std::vector cluster_centers(n_clusters); + if(n_clusters == 1) { + cluster_centers[0] = lambda_min; + } else { + const double dx = (lambda_max - lambda_min) / (n_clusters - 1); + for(std::size_t c = 0; c < n_clusters; ++c) { + cluster_centers[c] = lambda_min + static_cast(c) * dx; + } + } + + for(std::size_t i = 0; i < n; ++i) { + const auto cluster_id = i % n_clusters; + values[i] = cluster_centers[cluster_id] + jitter(gen); + } + return values; +} + +inline auto degenerate_spacing(std::size_t n, std::size_t n_clusters, + double lambda_min, double lambda_max) { + std::vector values(n); + if(n_clusters <= 1) { + std::fill(values.begin(), values.end(), lambda_min); + } else { + const auto n_plateaus = std::min(n_clusters, n); + const auto nm1 = std::max(1, n_plateaus - 1); + const double dx = (lambda_max - lambda_min) / static_cast(nm1); + for(std::size_t i = 0; i < n; ++i) { + const auto plateau = i % n_plateaus; + values[i] = lambda_min + static_cast(plateau) * dx; + } + } + return values; +} +} // namespace + +inline Tensor generate_eigenvalues(const SymmetricMatrixSpec& spec, + std::mt19937& gen) { + require_valid_n(spec.n); + const auto n = spec.n; + + const double lambda_min = spec.min_eigenvalue; + const double lambda_max = spec.min_eigenvalue * spec.condition_number; + + std::vector values; + + if(n == 1) return utilities::make_tensor({n}, values); + + switch(spec.spacing) { + case EigenvalueSpacing::Linear: { + values = linear_spacing(n, lambda_min, lambda_max); + break; + } + case EigenvalueSpacing::Logarithmic: { + values = logarithmic_spacing(n, lambda_min, lambda_max); + break; + } + case EigenvalueSpacing::Clustered: { + const auto n_clusters = std::max(1, spec.n_clusters); + const auto n_clusters_clamped = std::min(n_clusters, n); + values = clustered_spacing(n, n_clusters_clamped, lambda_min, + lambda_max, spec.cluster_width, gen); + break; + } + case EigenvalueSpacing::Degenerate: { + values = + degenerate_spacing(n, spec.n_clusters, lambda_min, lambda_max); + break; + } + default: { + throw std::invalid_argument("Invalid eigenvalue spacing"); + } + } + + std::sort(values.begin(), values.end()); + return utilities::make_tensor({n}, values); +} + +} // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/generate/generate_utils.hpp b/include/tensorwrapper/generate/generate_utils.hpp new file mode 100644 index 00000000..11192234 --- /dev/null +++ b/include/tensorwrapper/generate/generate_utils.hpp @@ -0,0 +1,62 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include +#include + +namespace tensorwrapper::generate { + +constexpr std::size_t kMaxMatrixDim = 10; + +enum class EigenvalueSpacing { + // delta = (lambda_max - lambda_min) / (n - 1). + Linear, + // delta = log(lambda_max / lambda_min) / (n - 1). + Logarithmic, + // clusters of width cluster_width are separated by + // (lambda_max - lambda_min) / (n_clusters - 1). + Clustered, + // same as clustered, but all eigenvalues in a cluster are the same. + Degenerate +}; + +struct SymmetricMatrixSpec { + std::size_t n = 2; + double condition_number = 10.0; + double min_eigenvalue = 1.0; + EigenvalueSpacing spacing = EigenvalueSpacing::Linear; + std::size_t n_clusters = 1; + double cluster_width = 1e-6; + std::uint64_t seed = 42; +}; + +inline std::mt19937 make_rng(std::uint64_t seed) { + if(seed == 0) { + std::random_device rd; + return std::mt19937(rd()); + } + return std::mt19937(static_cast(seed)); +} + +inline void require_valid_n(std::size_t n) { + if(n < 1 || n > kMaxMatrixDim) { + throw std::invalid_argument("n must be in [1, kMaxMatrixDim]"); + } +} + +} // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/generate/random_orthogonal_matrix.hpp b/include/tensorwrapper/generate/random_orthogonal_matrix.hpp new file mode 100644 index 00000000..85c8c366 --- /dev/null +++ b/include/tensorwrapper/generate/random_orthogonal_matrix.hpp @@ -0,0 +1,25 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include + +namespace tensorwrapper::generate { + +Tensor random_orthogonal_matrix(std::size_t n, std::mt19937& gen); + +} // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/tensorwrapper.hpp b/include/tensorwrapper/tensorwrapper.hpp index 861b0131..b6505f26 100644 --- a/include/tensorwrapper/tensorwrapper.hpp +++ b/include/tensorwrapper/tensorwrapper.hpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include diff --git a/src/tensorwrapper/generate/generate_eigen_system.cpp b/src/tensorwrapper/generate/generate_eigen_system.cpp new file mode 100644 index 00000000..c3b67d7f --- /dev/null +++ b/src/tensorwrapper/generate/generate_eigen_system.cpp @@ -0,0 +1,45 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include + +namespace tensorwrapper::generate { + +EigenSystem generate_eigen_system(const SymmetricMatrixSpec& spec) { + require_valid_n(spec.n); + auto gen = make_rng(spec.seed); + const auto n = spec.n; + + EigenSystem rv; + rv.n = n; + rv.eigenvalues = generate_eigenvalues(spec, gen); + rv.eigenvectors = random_orthogonal_matrix(n, gen); + + const auto D = utilities::diagonal_matrix(rv.eigenvalues); + + Tensor qd; + qd("i,k") = rv.eigenvectors("i,l") * D("l,k"); + + Tensor matrix; + matrix("i,j") = qd("i,k") * rv.eigenvectors("j,k"); + rv.matrix = std::move(matrix); + return rv; +} + +} // namespace tensorwrapper::generate diff --git a/src/tensorwrapper/generate/random_orthogonal_matrix.cpp b/src/tensorwrapper/generate/random_orthogonal_matrix.cpp new file mode 100644 index 00000000..3822c62f --- /dev/null +++ b/src/tensorwrapper/generate/random_orthogonal_matrix.cpp @@ -0,0 +1,49 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include + +namespace tensorwrapper::generate { + +Tensor random_orthogonal_matrix(std::size_t n, std::mt19937& gen) { + require_valid_n(n); + + Eigen::MatrixXd M(n, n); + std::normal_distribution dist(0.0, 1.0); + for(std::size_t i = 0; i < n; ++i) { + for(std::size_t j = 0; j < n; ++j) { + M(static_cast(i), static_cast(j)) = + dist(gen); + } + } + Eigen::HouseholderQR qr(M); + const Eigen::MatrixXd Q = qr.householderQ(); + + std::vector data(n * n); + for(std::size_t i = 0; i < n; ++i) { + for(std::size_t j = 0; j < n; ++j) { + data[i * n + j] = + Q(static_cast(i), static_cast(j)); + } + } + return utilities::make_tensor({n, n}, data); +} + +} // namespace tensorwrapper::generate diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigen_system.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigen_system.cpp new file mode 100644 index 00000000..3b72828b --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigen_system.cpp @@ -0,0 +1,127 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include +#include + +using namespace tensorwrapper; +using namespace tensorwrapper::buffer; +using namespace tensorwrapper::generate; +using namespace tensorwrapper::operations; +using namespace tensorwrapper::utilities; + +TEST_CASE("generate_eigen_system") { + SECTION("shapes") { + SymmetricMatrixSpec spec; + spec.n = 4; + spec.condition_number = 1e3; + spec.spacing = EigenvalueSpacing::Linear; + spec.seed = 11; + auto system = generate_eigen_system(spec); + + REQUIRE(system.n == 4); + REQUIRE( + make_contiguous(system.eigenvalues.buffer()).shape().extent(0) == 4); + REQUIRE( + make_contiguous(system.eigenvectors.buffer()).shape().extent(0) == 4); + REQUIRE( + make_contiguous(system.eigenvectors.buffer()).shape().extent(1) == 4); + REQUIRE(make_contiguous(system.matrix.buffer()).shape().extent(0) == 4); + REQUIRE(make_contiguous(system.matrix.buffer()).shape().extent(1) == 4); + } + + SECTION("symmetry") { + SymmetricMatrixSpec spec; + spec.n = 3; + spec.seed = 5; + auto system = generate_eigen_system(spec); + + Tensor sym; + sym("i,j") = system.matrix("i,j"); + Tensor tran; + tran("i,j") = system.matrix("j,i"); + REQUIRE(approximately_equal(sym, tran, 1e-12)); + } + + SECTION("orthonormal eigenvectors") { + SymmetricMatrixSpec spec; + spec.n = 3; + spec.seed = 7; + auto system = generate_eigen_system(spec); + + Tensor product; + product("i,k") = + system.eigenvectors("i,j") * system.eigenvectors("k,j"); + + auto ones = make_tensor({3}, std::vector{1, 1, 1}); + auto ident = diagonal_matrix(ones); + REQUIRE(approximately_equal(product, ident, 1e-12)); + } + + SECTION("exact eigenpair") { + SymmetricMatrixSpec spec; + spec.n = 4; + spec.seed = 13; + auto system = generate_eigen_system(spec); + + Tensor av; + av("i,k") = system.matrix("i,j") * system.eigenvectors("j,k"); + + const auto D = diagonal_matrix(system.eigenvalues); + Tensor scaled; + scaled("i,k") = system.eigenvectors("i,l") * D("l,k"); + REQUIRE(approximately_equal(av, scaled, 1e-12)); + } + + SECTION("deterministic for fixed seed") { + SymmetricMatrixSpec spec; + spec.n = 4; + spec.seed = 17; + auto system1 = generate_eigen_system(spec); + auto system2 = generate_eigen_system(spec); + REQUIRE(approximately_equal(system1.eigenvalues, system2.eigenvalues)); + REQUIRE( + approximately_equal(system1.eigenvectors, system2.eigenvectors)); + REQUIRE(approximately_equal(system1.matrix, system2.matrix)); + } + + SECTION("degenerate") { + SymmetricMatrixSpec spec; + spec.n = 2; + spec.condition_number = 1.0; + spec.min_eigenvalue = 1.0; + spec.spacing = EigenvalueSpacing::Degenerate; + spec.n_clusters = 1; + spec.seed = 7; + auto system = generate_eigen_system(spec); + + auto evals = make_tensor({2}, std::vector{1.0, 1.0}); + REQUIRE(approximately_equal(system.eigenvalues, evals)); + + auto ident = make_tensor({2, 2}, std::vector{1, 0, 0, 1}); + REQUIRE(approximately_equal(system.matrix, ident, 1e-12)); + } + + SECTION("invalid n throws") { + SymmetricMatrixSpec spec; + spec.n = 0; + REQUIRE_THROWS_AS(generate_eigen_system(spec), std::invalid_argument); + } +} diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigenvalues.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigenvalues.cpp new file mode 100644 index 00000000..98cc4e0f --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/generate/generate_eigenvalues.cpp @@ -0,0 +1,129 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace tensorwrapper; +using namespace tensorwrapper::buffer; +using namespace tensorwrapper::generate; +using namespace tensorwrapper::operations; +using namespace tensorwrapper::utilities; + +namespace { +double elem_as_double(const Contiguous::const_reference& elem) { + using wtf::fp::float_cast; + try { + return float_cast(elem); + } catch(const std::runtime_error&) { + return static_cast(float_cast(elem)); + } +} + +std::vector tensor_to_vector(const Tensor& t) { + auto buf = make_contiguous(t.buffer()); + std::vector rv(buf.shape().extent(0)); + for(std::size_t i = 0; i < rv.size(); ++i) { + rv[i] = elem_as_double(buf.get_elem({i})); + } + return rv; +} +} // namespace + +TEST_CASE("generate_eigenvalues") { + SECTION("linear spacing") { + SymmetricMatrixSpec spec; + spec.n = 4; + spec.min_eigenvalue = 1.0; + spec.condition_number = 10.0; + spec.spacing = EigenvalueSpacing::Linear; + auto gen = make_rng(1); + auto result = generate_eigenvalues(spec, gen); + auto corr = make_tensor({4}, std::vector{1, 4, 7, 10}); + REQUIRE(approximately_equal(result, corr)); + } + + SECTION("logarithmic spacing") { + SymmetricMatrixSpec spec; + spec.n = 3; + spec.min_eigenvalue = 1.0; + spec.condition_number = 100.0; + spec.spacing = EigenvalueSpacing::Logarithmic; + auto gen = make_rng(1); + auto values = tensor_to_vector(generate_eigenvalues(spec, gen)); + REQUIRE(values.size() == 3); + REQUIRE(values[0] == Catch::Approx(1.0)); + REQUIRE(values[1] == Catch::Approx(10.0)); + REQUIRE(values[2] == Catch::Approx(100.0)); + REQUIRE(std::is_sorted(values.begin(), values.end())); + } + + SECTION("degenerate spacing") { + SymmetricMatrixSpec spec; + spec.n = 3; + spec.min_eigenvalue = 2.5; + spec.condition_number = 10.0; + spec.spacing = EigenvalueSpacing::Degenerate; + spec.n_clusters = 1; + auto gen = make_rng(1); + auto values = tensor_to_vector(generate_eigenvalues(spec, gen)); + REQUIRE(values.size() == 3); + for(const auto v : values) { REQUIRE(v == Catch::Approx(2.5)); } + } + + SECTION("clustered spacing") { + SymmetricMatrixSpec spec; + spec.n = 6; + spec.min_eigenvalue = 1.0; + spec.condition_number = 100.0; + spec.spacing = EigenvalueSpacing::Clustered; + spec.n_clusters = 3; + spec.cluster_width = 1e-6; + spec.seed = 23; + auto gen = make_rng(spec.seed); + auto values = tensor_to_vector(generate_eigenvalues(spec, gen)); + REQUIRE(values.size() == 6); + REQUIRE(std::is_sorted(values.begin(), values.end())); + + const double lambda_max = spec.min_eigenvalue * spec.condition_number; + const double dx = (lambda_max - spec.min_eigenvalue) / 2.0; + std::vector centers(3); + for(std::size_t c = 0; c < 3; ++c) { + centers[c] = spec.min_eigenvalue + static_cast(c) * dx; + } + for(const auto v : values) { + const auto near_center = + std::any_of(centers.begin(), centers.end(), [&](double center) { + return std::abs(v - center) <= spec.cluster_width; + }); + REQUIRE(near_center); + } + } + + SECTION("invalid n throws") { + SymmetricMatrixSpec spec; + spec.n = 0; + auto gen = make_rng(1); + REQUIRE_THROWS_AS(generate_eigenvalues(spec, gen), + std::invalid_argument); + } +} diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/generate_utils.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/generate_utils.cpp new file mode 100644 index 00000000..04b6b98f --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/generate/generate_utils.cpp @@ -0,0 +1,35 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include + +using namespace tensorwrapper::generate; + +TEST_CASE("generate_utils") { + SECTION("make_rng is deterministic for fixed seed") { + auto gen1 = make_rng(42); + auto gen2 = make_rng(42); + REQUIRE(gen1() == gen2()); + } + + SECTION("require_valid_n throws for invalid n") { + REQUIRE_THROWS_AS(require_valid_n(0), std::invalid_argument); + REQUIRE_THROWS_AS(require_valid_n(11), std::invalid_argument); + REQUIRE_NOTHROW(require_valid_n(1)); + REQUIRE_NOTHROW(require_valid_n(10)); + } +} diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/random_orthogonal_matrix.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/random_orthogonal_matrix.cpp new file mode 100644 index 00000000..6761c082 --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/generate/random_orthogonal_matrix.cpp @@ -0,0 +1,58 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include +#include + +using namespace tensorwrapper; +using namespace tensorwrapper::buffer; +using namespace tensorwrapper::generate; +using namespace tensorwrapper::operations; +using namespace tensorwrapper::utilities; + +TEST_CASE("random_orthogonal_matrix") { + SECTION("shape") { + auto gen = make_rng(7); + auto Q = random_orthogonal_matrix(3, gen); + auto buf = make_contiguous(Q.buffer()); + REQUIRE(buf.shape().rank() == 2); + REQUIRE(buf.shape().extent(0) == 3); + REQUIRE(buf.shape().extent(1) == 3); + } + + SECTION("orthogonality") { + auto gen = make_rng(11); + auto Q = random_orthogonal_matrix(3, gen); + Tensor product; + product("i,k") = Q("i,j") * Q("k,j"); + + auto ones = make_tensor({3}, std::vector{1, 1, 1}); + auto ident = diagonal_matrix(ones); + REQUIRE(approximately_equal(product, ident, 1e-12)); + } + + SECTION("deterministic for fixed seed") { + auto gen1 = make_rng(99); + auto gen2 = make_rng(99); + auto Q1 = random_orthogonal_matrix(4, gen1); + auto Q2 = random_orthogonal_matrix(4, gen2); + REQUIRE(approximately_equal(Q1, Q2)); + } +} From 893ebd442b6b629f20b8f3a7a7a0bb19fdf1b6bb Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Fri, 29 May 2026 14:14:38 -0500 Subject: [PATCH 2/4] generate noisy tensor --- include/tensorwrapper/generate/add_noise.hpp | 41 +++++++++ include/tensorwrapper/generate/generate.hpp | 1 + src/tensorwrapper/generate/add_noise.cpp | 58 +++++++++++++ .../tensorwrapper/generate/add_noise.cpp | 87 +++++++++++++++++++ 4 files changed, 187 insertions(+) create mode 100644 include/tensorwrapper/generate/add_noise.hpp create mode 100644 src/tensorwrapper/generate/add_noise.cpp create mode 100644 tests/cxx/unit_tests/tensorwrapper/generate/add_noise.cpp diff --git a/include/tensorwrapper/generate/add_noise.hpp b/include/tensorwrapper/generate/add_noise.hpp new file mode 100644 index 00000000..50fc47a8 --- /dev/null +++ b/include/tensorwrapper/generate/add_noise.hpp @@ -0,0 +1,41 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include +#include + +namespace tensorwrapper::generate { + +/** @brief Adds clamped normal noise to each element of @p matrix. + * + * Draws `delta ~ Normal(0, t)` and clamps to `[-t, t]` before adding to each + * element. When @p t is zero the input is copied unchanged. + * + * @param[in] matrix The tensor to perturb. + * @param[in] t Non-negative noise scale (standard deviation and clamp bound). + * @param[in,out] gen Random number generator used for the normal draws. + * + * @return A new tensor with the same shape as @p matrix. + * + * @throw std::invalid_argument if @p t is negative. + */ +Tensor add_noise(const Tensor& matrix, double t, std::mt19937& gen); + +Tensor add_noise(const Tensor& matrix, double t, std::uint64_t seed = 42); + +} // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/generate/generate.hpp b/include/tensorwrapper/generate/generate.hpp index 624e4f4c..367436a3 100644 --- a/include/tensorwrapper/generate/generate.hpp +++ b/include/tensorwrapper/generate/generate.hpp @@ -15,6 +15,7 @@ */ #pragma once +#include #include #include #include diff --git a/src/tensorwrapper/generate/add_noise.cpp b/src/tensorwrapper/generate/add_noise.cpp new file mode 100644 index 00000000..f6cf6b08 --- /dev/null +++ b/src/tensorwrapper/generate/add_noise.cpp @@ -0,0 +1,58 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace tensorwrapper::generate { +namespace { + +std::vector shape_extents(const auto& shape) { + std::vector rv(shape.rank()); + for(std::size_t d = 0; d < shape.rank(); ++d) rv[d] = shape.extent(d); + return rv; +} + +} // namespace + +Tensor add_noise(const Tensor& matrix, double t, std::mt19937& gen) { + if(t < 0.0) { throw std::invalid_argument("t must be non-negative"); } + + auto& in_buf = buffer::make_contiguous(matrix.buffer()); + auto in_data = buffer::get_raw_data(in_buf); + std::vector data(in_data.begin(), in_data.end()); + + if(t > 0.0) { + std::normal_distribution dist(0.0, t); + for(auto& x : data) { x += std::clamp(dist(gen), -t, t); } + } + + return utilities::make_tensor(shape_extents(in_buf.shape()), data.begin(), + data.end()); +} + +Tensor add_noise(const Tensor& matrix, double t, std::uint64_t seed) { + auto gen = make_rng(seed); + return add_noise(matrix, t, gen); +} + +} // namespace tensorwrapper::generate diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/add_noise.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/add_noise.cpp new file mode 100644 index 00000000..9aaca48f --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/generate/add_noise.cpp @@ -0,0 +1,87 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include +#include + +using namespace tensorwrapper; +using namespace tensorwrapper::buffer; +using namespace tensorwrapper::generate; +using namespace tensorwrapper::operations; +using namespace tensorwrapper::utilities; + +namespace { +double elem_as_double(const Contiguous::const_reference& elem) { + using wtf::fp::float_cast; + try { + return float_cast(elem); + } catch(const std::runtime_error&) { + return static_cast(float_cast(elem)); + } +} +} // namespace + +TEST_CASE("add_noise") { + const auto matrix = make_tensor({2, 2}, std::vector{1, 2, 3, 4}); + + SECTION("t equals zero returns copy") { + auto gen = make_rng(1); + auto out = add_noise(matrix, 0.0, gen); + REQUIRE(approximately_equal(out, matrix)); + } + + SECTION("within t of input") { + const double t = 0.05; + auto gen = make_rng(17); + auto out = add_noise(matrix, t, gen); + + auto in_buf = make_contiguous(matrix.buffer()); + auto out_buf = make_contiguous(out.buffer()); + for(std::size_t i = 0; i < 2; ++i) { + for(std::size_t j = 0; j < 2; ++j) { + const auto x = elem_as_double(in_buf.get_elem({i, j})); + const auto y = elem_as_double(out_buf.get_elem({i, j})); + REQUIRE(std::abs(y - x) <= t); + } + } + } + + SECTION("deterministic for fixed seed") { + const double t = 0.01; + auto out1 = add_noise(matrix, t, 99); + auto out2 = add_noise(matrix, t, 99); + REQUIRE(approximately_equal(out1, out2)); + } + + SECTION("shape preserved") { + auto gen = make_rng(3); + auto out = add_noise(matrix, 0.01, gen); + auto in_shape = make_contiguous(matrix.buffer()).shape(); + auto out_shape = make_contiguous(out.buffer()).shape(); + REQUIRE(out_shape.rank() == in_shape.rank()); + REQUIRE(out_shape.extent(0) == in_shape.extent(0)); + REQUIRE(out_shape.extent(1) == in_shape.extent(1)); + } + + SECTION("invalid t throws") { + auto gen = make_rng(1); + REQUIRE_THROWS_AS(add_noise(matrix, -1.0, gen), std::invalid_argument); + } +} From a0ccb6969afe69d092a526a6c4236c4cf04a2e86 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Fri, 29 May 2026 14:19:02 -0500 Subject: [PATCH 3/4] adds identity matrix --- include/tensorwrapper/generate/generate.hpp | 1 + .../generate/identity_matrix.hpp | 33 ++++++++++++++ .../generate/identity_matrix.cpp | 44 +++++++++++++++++++ 3 files changed, 78 insertions(+) create mode 100644 include/tensorwrapper/generate/identity_matrix.hpp create mode 100644 tests/cxx/unit_tests/tensorwrapper/generate/identity_matrix.cpp diff --git a/include/tensorwrapper/generate/generate.hpp b/include/tensorwrapper/generate/generate.hpp index 367436a3..f92e96d5 100644 --- a/include/tensorwrapper/generate/generate.hpp +++ b/include/tensorwrapper/generate/generate.hpp @@ -19,6 +19,7 @@ #include #include #include +#include #include namespace tensorwrapper::generate {} diff --git a/include/tensorwrapper/generate/identity_matrix.hpp b/include/tensorwrapper/generate/identity_matrix.hpp new file mode 100644 index 00000000..fd386f09 --- /dev/null +++ b/include/tensorwrapper/generate/identity_matrix.hpp @@ -0,0 +1,33 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#pragma once +#include +#include +#include +#include +#include + +namespace tensorwrapper::generate { + +inline Tensor identity_matrix(std::size_t n) { + require_valid_n(n); + std::vector values(n, 1.0); + auto values_tensor = utilities::make_tensor({n}, values); + return utilities::diagonal_matrix(values_tensor); +} + +} // namespace tensorwrapper::generate diff --git a/tests/cxx/unit_tests/tensorwrapper/generate/identity_matrix.cpp b/tests/cxx/unit_tests/tensorwrapper/generate/identity_matrix.cpp new file mode 100644 index 00000000..c092a8de --- /dev/null +++ b/tests/cxx/unit_tests/tensorwrapper/generate/identity_matrix.cpp @@ -0,0 +1,44 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include + +using namespace tensorwrapper; +using namespace tensorwrapper::generate; +using namespace tensorwrapper::operations; +using namespace tensorwrapper::utilities; + +TEST_CASE("identity_matrix") { + SECTION("2 by 2") { + auto result = identity_matrix(2); + auto corr = make_tensor({2, 2}, std::vector{1, 0, 0, 1}); + REQUIRE(approximately_equal(result, corr)); + } + + SECTION("3 by 3") { + auto result = identity_matrix(3); + auto corr = + make_tensor({3, 3}, std::vector{1, 0, 0, 0, 1, 0, 0, 0, 1}); + REQUIRE(approximately_equal(result, corr)); + } + + SECTION("invalid n throws") { + REQUIRE_THROWS_AS(identity_matrix(0), std::invalid_argument); + } +} From 2ebe656973c197b38be735fa384aac10862cc821 Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Fri, 29 May 2026 22:19:00 -0500 Subject: [PATCH 4/4] adds documentation, refactors generate_eigenvalues --- include/tensorwrapper/generate/add_noise.hpp | 11 ++ include/tensorwrapper/generate/generate.hpp | 6 + .../generate/generate_eigen_system.hpp | 23 +++ .../generate/generate_eigenvalues.hpp | 128 ++----------- .../tensorwrapper/generate/generate_utils.hpp | 63 +++++-- .../generate/identity_matrix.hpp | 8 + .../generate/random_orthogonal_matrix.hpp | 12 ++ .../generate/generate_eigenvalues.cpp | 172 ++++++++++++++++++ 8 files changed, 304 insertions(+), 119 deletions(-) create mode 100644 src/tensorwrapper/generate/generate_eigenvalues.cpp diff --git a/include/tensorwrapper/generate/add_noise.hpp b/include/tensorwrapper/generate/add_noise.hpp index 50fc47a8..ee10a448 100644 --- a/include/tensorwrapper/generate/add_noise.hpp +++ b/include/tensorwrapper/generate/add_noise.hpp @@ -36,6 +36,17 @@ namespace tensorwrapper::generate { */ Tensor add_noise(const Tensor& matrix, double t, std::mt19937& gen); +/** @brief Overload of add_noise that creates its own RNG from @p seed. + * + * @param[in] matrix The tensor to perturb. + * @param[in] t Non-negative noise scale (standard deviation and clamp bound). + * @param[in] seed Seed for the internal random number generator. A value of + * zero selects a non-deterministic seed. + * + * @return A new tensor with the same shape as @p matrix. + * + * @throw std::invalid_argument if @p t is negative. + */ Tensor add_noise(const Tensor& matrix, double t, std::uint64_t seed = 42); } // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/generate/generate.hpp b/include/tensorwrapper/generate/generate.hpp index f92e96d5..7b672c65 100644 --- a/include/tensorwrapper/generate/generate.hpp +++ b/include/tensorwrapper/generate/generate.hpp @@ -22,4 +22,10 @@ #include #include +/** @brief Namespace for utilities that generate synthetic tensors and matrices. + * + * The functions in this namespace create reproducible test matrices, + * eigenvalue + * spectra, and related quantities for unit testing and benchmarking. + */ namespace tensorwrapper::generate {} diff --git a/include/tensorwrapper/generate/generate_eigen_system.hpp b/include/tensorwrapper/generate/generate_eigen_system.hpp index 9e288569..cbaeb653 100644 --- a/include/tensorwrapper/generate/generate_eigen_system.hpp +++ b/include/tensorwrapper/generate/generate_eigen_system.hpp @@ -21,13 +21,36 @@ namespace tensorwrapper::generate { +/** @brief A symmetric matrix together with its eigen-decomposition. + * + * The member @p matrix satisfies @f$M = Q D Q^T@f$ where @p eigenvectors is + * @f$Q@f$ and @p eigenvalues holds the diagonal entries of @f$D@f$. + */ struct EigenSystem { + /// Dimension of the square matrix. std::size_t n = 0; + /// Rank-1 tensor of length @p n containing the eigenvalues. Tensor eigenvalues; + /// Rank-2 tensor of shape `(n, n)` whose columns are the eigenvectors. Tensor eigenvectors; + /// Rank-2 tensor of shape `(n, n)` representing the symmetric matrix. Tensor matrix; }; +/** @brief Generates a reproducible symmetric matrix and its + * eigen-decomposition. + * + * Constructs @f$M = Q D Q^T@f$ where @f$D@f$ is built from eigenvalues + * generated according to @p spec and @f$Q@f$ is a random orthogonal matrix. + * + * @param[in] spec Parameters controlling the matrix dimension, spectrum, and + * random seed. + * + * @return An @ref EigenSystem populated with @p spec.n, the eigenvalues, + * eigenvectors, and the assembled matrix. + * + * @throw std::invalid_argument if @p spec.n is outside `[1, kMaxMatrixDim]`. + */ EigenSystem generate_eigen_system(const SymmetricMatrixSpec& spec); } // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/generate/generate_eigenvalues.hpp b/include/tensorwrapper/generate/generate_eigenvalues.hpp index 9f1ef47d..2ae3d46c 100644 --- a/include/tensorwrapper/generate/generate_eigenvalues.hpp +++ b/include/tensorwrapper/generate/generate_eigenvalues.hpp @@ -15,117 +15,29 @@ */ #pragma once -#include -#include +#include #include -#include -#include +#include namespace tensorwrapper::generate { -namespace { -inline auto linear_spacing(std::size_t n, double lambda_min, - double lambda_max) { - const auto dx = (lambda_max - lambda_min) / (n - 1); - std::vector values(n); - for(std::size_t i = 0; i < n; ++i) { - values[i] = lambda_min + static_cast(i) * dx; - } - return values; -} -inline auto logarithmic_spacing(std::size_t n, double lambda_min, - double lambda_max) { - const double log_min = std::log(lambda_min); - const double log_max = std::log(lambda_max); - const double dlog = (log_max - log_min) / (n - 1); - std::vector values(n); - for(std::size_t i = 0; i < n; ++i) { - values[i] = std::exp(log_min + static_cast(i) * dlog); - } - return values; -} - -inline auto clustered_spacing(std::size_t n, std::size_t n_clusters, - double lambda_min, double lambda_max, - double cluster_width, std::mt19937& gen) { - std::uniform_real_distribution jitter(-cluster_width, - cluster_width); - std::vector values(n); - std::vector cluster_centers(n_clusters); - if(n_clusters == 1) { - cluster_centers[0] = lambda_min; - } else { - const double dx = (lambda_max - lambda_min) / (n_clusters - 1); - for(std::size_t c = 0; c < n_clusters; ++c) { - cluster_centers[c] = lambda_min + static_cast(c) * dx; - } - } - - for(std::size_t i = 0; i < n; ++i) { - const auto cluster_id = i % n_clusters; - values[i] = cluster_centers[cluster_id] + jitter(gen); - } - return values; -} - -inline auto degenerate_spacing(std::size_t n, std::size_t n_clusters, - double lambda_min, double lambda_max) { - std::vector values(n); - if(n_clusters <= 1) { - std::fill(values.begin(), values.end(), lambda_min); - } else { - const auto n_plateaus = std::min(n_clusters, n); - const auto nm1 = std::max(1, n_plateaus - 1); - const double dx = (lambda_max - lambda_min) / static_cast(nm1); - for(std::size_t i = 0; i < n; ++i) { - const auto plateau = i % n_plateaus; - values[i] = lambda_min + static_cast(plateau) * dx; - } - } - return values; -} -} // namespace - -inline Tensor generate_eigenvalues(const SymmetricMatrixSpec& spec, - std::mt19937& gen) { - require_valid_n(spec.n); - const auto n = spec.n; - - const double lambda_min = spec.min_eigenvalue; - const double lambda_max = spec.min_eigenvalue * spec.condition_number; - - std::vector values; - - if(n == 1) return utilities::make_tensor({n}, values); - - switch(spec.spacing) { - case EigenvalueSpacing::Linear: { - values = linear_spacing(n, lambda_min, lambda_max); - break; - } - case EigenvalueSpacing::Logarithmic: { - values = logarithmic_spacing(n, lambda_min, lambda_max); - break; - } - case EigenvalueSpacing::Clustered: { - const auto n_clusters = std::max(1, spec.n_clusters); - const auto n_clusters_clamped = std::min(n_clusters, n); - values = clustered_spacing(n, n_clusters_clamped, lambda_min, - lambda_max, spec.cluster_width, gen); - break; - } - case EigenvalueSpacing::Degenerate: { - values = - degenerate_spacing(n, spec.n_clusters, lambda_min, lambda_max); - break; - } - default: { - throw std::invalid_argument("Invalid eigenvalue spacing"); - } - } - - std::sort(values.begin(), values.end()); - return utilities::make_tensor({n}, values); -} +/** @brief Generates a sorted rank-1 tensor of eigenvalues from @p spec. + * + * The spectrum spans + * @f$[\texttt{spec.min\_eigenvalue}, + * \texttt{spec.min\_eigenvalue} \times \texttt{spec.condition\_number}]@f$ + * using the spacing strategy given by @p spec.spacing. + * + * @param[in] spec Parameters controlling the eigenvalue distribution. + * @param[in,out] gen Random number generator used when @p spec.spacing + * requires random draws. + * + * @return A rank-1 tensor of length @p spec.n containing the eigenvalues in + * ascending order. + * + * @throw std::invalid_argument if @p spec.n is outside `[1, kMaxMatrixDim]` or + * if @p spec.spacing is invalid. + */ +Tensor generate_eigenvalues(const SymmetricMatrixSpec& spec, std::mt19937& gen); } // namespace tensorwrapper::generate diff --git a/include/tensorwrapper/generate/generate_utils.hpp b/include/tensorwrapper/generate/generate_utils.hpp index 11192234..15a2bfb2 100644 --- a/include/tensorwrapper/generate/generate_utils.hpp +++ b/include/tensorwrapper/generate/generate_utils.hpp @@ -21,30 +21,65 @@ namespace tensorwrapper::generate { +/// Maximum supported dimension for generated square matrices. constexpr std::size_t kMaxMatrixDim = 10; +/** @brief Specifies how eigenvalues are distributed between the endpoints. + * + * Each spacing mode fills the interval + * @f$[\lambda_{\min}, \lambda_{\max}]@f$ with @f$n@f$ eigenvalues, where + * @f$\lambda_{\max} = \lambda_{\min} \times@f$ the condition number. + */ enum class EigenvalueSpacing { - // delta = (lambda_max - lambda_min) / (n - 1). + /// Uniform spacing with @f$\Delta\lambda = (\lambda_{\max} - + /// \lambda_{\min}) / (n - 1)@f$. Linear, - // delta = log(lambda_max / lambda_min) / (n - 1). + /// Uniform spacing in log space with + /// @f$\Delta\log\lambda = \log(\lambda_{\max} / \lambda_{\min}) / (n - + /// 1)@f$. Logarithmic, - // clusters of width cluster_width are separated by - // (lambda_max - lambda_min) / (n_clusters - 1). + /// Eigenvalues are grouped into clusters of width @p cluster_width that are + /// separated by @f$(\lambda_{\max} - \lambda_{\min}) / (n_{\text{clusters}} + /// - 1)@f$. Clustered, - // same as clustered, but all eigenvalues in a cluster are the same. + /// Same cluster centers as @ref Clustered, but all eigenvalues in a cluster + /// are identical. Degenerate }; +/** @brief Parameters controlling the generation of a symmetric test matrix. + * + * The resulting matrix has eigenvalues in + * @f$[\texttt{min\_eigenvalue}, + * \texttt{min\_eigenvalue} \times \texttt{condition\_number}]@f$ with the + * spacing determined by @p spacing. + */ struct SymmetricMatrixSpec { - std::size_t n = 2; - double condition_number = 10.0; - double min_eigenvalue = 1.0; + /// Dimension of the square matrix. + std::size_t n = 2; + /// Ratio of the largest to smallest eigenvalue. + double condition_number = 10.0; + /// Smallest eigenvalue in the spectrum. + double min_eigenvalue = 1.0; + /// Strategy used to distribute eigenvalues between the endpoints. EigenvalueSpacing spacing = EigenvalueSpacing::Linear; - std::size_t n_clusters = 1; - double cluster_width = 1e-6; - std::uint64_t seed = 42; + /// Number of eigenvalue clusters when @p spacing is @ref Clustered or + /// @ref Degenerate. + std::size_t n_clusters = 1; + /// Half-width of each cluster when @p spacing is @ref Clustered. + double cluster_width = 1e-6; + /// Seed for random number generation. A value of zero selects a + /// non-deterministic seed. + std::uint64_t seed = 42; }; +/** @brief Creates a Mersenne Twister RNG from @p seed. + * + * @param[in] seed The seed value. When zero, a seed is drawn from + * `std::random_device`. + * + * @return A `std::mt19937` generator initialized with @p seed. + */ inline std::mt19937 make_rng(std::uint64_t seed) { if(seed == 0) { std::random_device rd; @@ -53,6 +88,12 @@ inline std::mt19937 make_rng(std::uint64_t seed) { return std::mt19937(static_cast(seed)); } +/** @brief Validates that @p n is an allowed matrix dimension. + * + * @param[in] n The dimension to validate. + * + * @throw std::invalid_argument if @p n is not in `[1, kMaxMatrixDim]`. + */ inline void require_valid_n(std::size_t n) { if(n < 1 || n > kMaxMatrixDim) { throw std::invalid_argument("n must be in [1, kMaxMatrixDim]"); diff --git a/include/tensorwrapper/generate/identity_matrix.hpp b/include/tensorwrapper/generate/identity_matrix.hpp index fd386f09..9e55f45e 100644 --- a/include/tensorwrapper/generate/identity_matrix.hpp +++ b/include/tensorwrapper/generate/identity_matrix.hpp @@ -23,6 +23,14 @@ namespace tensorwrapper::generate { +/** @brief Creates an @p n x @p n identity matrix. + * + * @param[in] n Matrix dimension. Must be in `[1, kMaxMatrixDim]`. + * + * @return A rank-2 tensor with ones on the diagonal and zeros elsewhere. + * + * @throw std::invalid_argument if @p n is outside the allowed range. + */ inline Tensor identity_matrix(std::size_t n) { require_valid_n(n); std::vector values(n, 1.0); diff --git a/include/tensorwrapper/generate/random_orthogonal_matrix.hpp b/include/tensorwrapper/generate/random_orthogonal_matrix.hpp index 85c8c366..0be2a3f9 100644 --- a/include/tensorwrapper/generate/random_orthogonal_matrix.hpp +++ b/include/tensorwrapper/generate/random_orthogonal_matrix.hpp @@ -20,6 +20,18 @@ namespace tensorwrapper::generate { +/** @brief Creates a random @p n x @p n orthogonal matrix. + * + * Draws entries from a standard normal distribution and applies Householder QR + * factorization to obtain an orthogonal matrix @f$Q@f$. + * + * @param[in] n Matrix dimension. Must be in `[1, kMaxMatrixDim]`. + * @param[in,out] gen Random number generator used for the normal draws. + * + * @return A rank-2 tensor whose columns form an orthonormal basis. + * + * @throw std::invalid_argument if @p n is outside the allowed range. + */ Tensor random_orthogonal_matrix(std::size_t n, std::mt19937& gen); } // namespace tensorwrapper::generate diff --git a/src/tensorwrapper/generate/generate_eigenvalues.cpp b/src/tensorwrapper/generate/generate_eigenvalues.cpp new file mode 100644 index 00000000..1a95eca1 --- /dev/null +++ b/src/tensorwrapper/generate/generate_eigenvalues.cpp @@ -0,0 +1,172 @@ +/* + * Copyright 2026 NWChemEx-Project + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include + +namespace tensorwrapper::generate { +namespace { + +/** @brief Generates @p n eigenvalues with uniform spacing in @f$\lambda@f$. + * + * @param[in] n Number of eigenvalues. + * @param[in] lambda_min Smallest eigenvalue. + * @param[in] lambda_max Largest eigenvalue. + * + * @return A vector of length @p n with values from @p lambda_min to + * @p lambda_max. + */ +auto linear_spacing(std::size_t n, double lambda_min, double lambda_max) { + const auto dx = (lambda_max - lambda_min) / (n - 1); + std::vector values(n); + for(std::size_t i = 0; i < n; ++i) { + values[i] = lambda_min + static_cast(i) * dx; + } + return values; +} + +/** @brief Generates @p n eigenvalues with uniform spacing in @f$\log\lambda@f$. + * + * @param[in] n Number of eigenvalues. + * @param[in] lambda_min Smallest eigenvalue. Must be positive. + * @param[in] lambda_max Largest eigenvalue. Must be positive. + * + * @return A vector of length @p n with values from @p lambda_min to + * @p lambda_max on a log scale. + */ +auto logarithmic_spacing(std::size_t n, double lambda_min, double lambda_max) { + const double log_min = std::log(lambda_min); + const double log_max = std::log(lambda_max); + const double dlog = (log_max - log_min) / (n - 1); + std::vector values(n); + for(std::size_t i = 0; i < n; ++i) { + values[i] = std::exp(log_min + static_cast(i) * dlog); + } + return values; +} + +/** @brief Generates @p n eigenvalues grouped into @p n_clusters jittered + * clusters. + * + * @param[in] n Number of eigenvalues. + * @param[in] n_clusters Number of cluster centers. + * @param[in] lambda_min Smallest cluster center. + * @param[in] lambda_max Largest cluster center. + * @param[in] cluster_width Half-width of the uniform jitter around each + * center. + * @param[in,out] gen Random number generator used for the jitter draws. + * + * @return A vector of length @p n with eigenvalues assigned cyclically to + * clusters. + */ +auto clustered_spacing(std::size_t n, std::size_t n_clusters, double lambda_min, + double lambda_max, double cluster_width, + std::mt19937& gen) { + std::uniform_real_distribution jitter(-cluster_width, + cluster_width); + std::vector values(n); + std::vector cluster_centers(n_clusters); + if(n_clusters == 1) { + cluster_centers[0] = lambda_min; + } else { + const double dx = (lambda_max - lambda_min) / (n_clusters - 1); + for(std::size_t c = 0; c < n_clusters; ++c) { + cluster_centers[c] = lambda_min + static_cast(c) * dx; + } + } + + for(std::size_t i = 0; i < n; ++i) { + const auto cluster_id = i % n_clusters; + values[i] = cluster_centers[cluster_id] + jitter(gen); + } + return values; +} + +/** @brief Generates @p n eigenvalues with exact degeneracies within each + * cluster. + * + * @param[in] n Number of eigenvalues. + * @param[in] n_clusters Number of distinct eigenvalue plateaus. + * @param[in] lambda_min Value of the smallest plateau. + * @param[in] lambda_max Value of the largest plateau. + * + * @return A vector of length @p n with eigenvalues assigned cyclically to + * @p n_clusters plateaus. + */ +auto degenerate_spacing(std::size_t n, std::size_t n_clusters, + double lambda_min, double lambda_max) { + std::vector values(n); + if(n_clusters <= 1) { + std::fill(values.begin(), values.end(), lambda_min); + } else { + const auto n_plateaus = std::min(n_clusters, n); + const auto nm1 = std::max(1, n_plateaus - 1); + const double dx = (lambda_max - lambda_min) / static_cast(nm1); + for(std::size_t i = 0; i < n; ++i) { + const auto plateau = i % n_plateaus; + values[i] = lambda_min + static_cast(plateau) * dx; + } + } + return values; +} +} // namespace + +Tensor generate_eigenvalues(const SymmetricMatrixSpec& spec, + std::mt19937& gen) { + require_valid_n(spec.n); + const auto n = spec.n; + + const double lambda_min = spec.min_eigenvalue; + const double lambda_max = spec.min_eigenvalue * spec.condition_number; + + std::vector values; + + if(n == 1) return utilities::make_tensor({n}, values); + + switch(spec.spacing) { + case EigenvalueSpacing::Linear: { + values = linear_spacing(n, lambda_min, lambda_max); + break; + } + case EigenvalueSpacing::Logarithmic: { + values = logarithmic_spacing(n, lambda_min, lambda_max); + break; + } + case EigenvalueSpacing::Clustered: { + const auto n_clusters = std::max(1, spec.n_clusters); + const auto n_clusters_clamped = std::min(n_clusters, n); + values = clustered_spacing(n, n_clusters_clamped, lambda_min, + lambda_max, spec.cluster_width, gen); + break; + } + case EigenvalueSpacing::Degenerate: { + values = + degenerate_spacing(n, spec.n_clusters, lambda_min, lambda_max); + break; + } + default: { + throw std::invalid_argument("Invalid eigenvalue spacing"); + } + } + + std::sort(values.begin(), values.end()); + return utilities::make_tensor({n}, values); +} + +} // namespace tensorwrapper::generate