From 03318f611134da67336f80c47aa45160304a7a17 Mon Sep 17 00:00:00 2001 From: Arthur Lin Date: Thu, 11 Jan 2024 18:36:41 -0600 Subject: [PATCH 1/4] filtered out small evals and modified correctness check --- anisoap/representations/radial_basis.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/anisoap/representations/radial_basis.py b/anisoap/representations/radial_basis.py index 6f84f73..c3ad58d 100644 --- a/anisoap/representations/radial_basis.py +++ b/anisoap/representations/radial_basis.py @@ -22,12 +22,17 @@ def inverse_matrix_sqrt(matrix: np.array): if not np.allclose(matrix, matrix.conjugate().T): raise ValueError("Matrix is not hermitian") eva, eve = np.linalg.eigh(matrix) + eve = eve[:, eva > 1e-8] + eva = eva[eva > 1e-8] - if (eva < 0).any(): - raise ValueError( - "Matrix is not positive semidefinite. Check that a valid gram matrix is passed." - ) - return eve @ np.diag(1 / np.sqrt(eva)) @ eve.T + result = eve @ np.diag(1 / np.sqrt(eva)) @ eve.T + + # Do quick test to make sure inverse matrix sqrt succeeded (should succeed for overlap matrices up to order 100) + # If it doesn't, matrix likely wasn't a gram matrix to start with. + matrix2 = np.linalg.pinv(result @ result) # this should be very close to the original matrix + if np.linalg.norm(matrix-matrix2) > 1.0e-3: + raise ValueError(f"Incurred Numerical Imprecision {np.linalg.norm(matrix-matrix2)= :.3f}") + return result def gto_square_norm(n, sigma): From a9e8799806c1cf9f16e5f3712e70d07afba7ed10 Mon Sep 17 00:00:00 2001 From: Arthur Lin Date: Thu, 11 Jan 2024 18:40:45 -0600 Subject: [PATCH 2/4] filtered out small evals and modified correctness check --- anisoap/representations/radial_basis.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/anisoap/representations/radial_basis.py b/anisoap/representations/radial_basis.py index c3ad58d..14f858c 100644 --- a/anisoap/representations/radial_basis.py +++ b/anisoap/representations/radial_basis.py @@ -27,11 +27,13 @@ def inverse_matrix_sqrt(matrix: np.array): result = eve @ np.diag(1 / np.sqrt(eva)) @ eve.T - # Do quick test to make sure inverse matrix sqrt succeeded (should succeed for overlap matrices up to order 100) - # If it doesn't, matrix likely wasn't a gram matrix to start with. - matrix2 = np.linalg.pinv(result @ result) # this should be very close to the original matrix - if np.linalg.norm(matrix-matrix2) > 1.0e-3: - raise ValueError(f"Incurred Numerical Imprecision {np.linalg.norm(matrix-matrix2)= :.3f}") + # Do quick test to make sure inverse square of the inverse matrix sqrt succeeded + # This should succed for most cases (e.g. GTO orders up to 100), if not the matrix likely wasn't a gram matrix to start with. + matrix2 = np.linalg.pinv(result @ result) + if np.linalg.norm(matrix - matrix2) > 1.0e-3: + raise ValueError( + f"Incurred Numerical Imprecision {np.linalg.norm(matrix-matrix2)= :.3f}" + ) return result From 860095912cf7d420f581e053eb0e9160de0b24ca Mon Sep 17 00:00:00 2001 From: "Rose K. Cersonsky" <47536110+rosecers@users.noreply.github.com> Date: Fri, 12 Jan 2024 11:58:13 -0600 Subject: [PATCH 3/4] Moved rcond and tolerance to definable params --- anisoap/representations/radial_basis.py | 30 +++++++++++++++++++------ 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/anisoap/representations/radial_basis.py b/anisoap/representations/radial_basis.py index 14f858c..4a824f2 100644 --- a/anisoap/representations/radial_basis.py +++ b/anisoap/representations/radial_basis.py @@ -1,12 +1,11 @@ import warnings import numpy as np -import scipy.linalg from metatensor import TensorMap from scipy.special import gamma -def inverse_matrix_sqrt(matrix: np.array): +def inverse_matrix_sqrt(matrix: np.array, rcond=1e-8, tol=1e-3): """ Returns the inverse matrix square root. The inverse square root of the overlap matrix (or slices of the overlap matrix) yields the @@ -15,6 +14,13 @@ def inverse_matrix_sqrt(matrix: np.array): matrix: np.array Symmetric square matrix to find the inverse square root of + rcond: float + Lower bound for eigenvalues for inverse square root + + tol: float + Tolerance for differences between original matrix and reconstruction via + inverse square root + Returns: inverse_sqrt_matrix: S^{-1/2} @@ -22,15 +28,15 @@ def inverse_matrix_sqrt(matrix: np.array): if not np.allclose(matrix, matrix.conjugate().T): raise ValueError("Matrix is not hermitian") eva, eve = np.linalg.eigh(matrix) - eve = eve[:, eva > 1e-8] - eva = eva[eva > 1e-8] + eve = eve[:, eva > rcond] + eva = eva[eva > rcond] result = eve @ np.diag(1 / np.sqrt(eva)) @ eve.T # Do quick test to make sure inverse square of the inverse matrix sqrt succeeded # This should succed for most cases (e.g. GTO orders up to 100), if not the matrix likely wasn't a gram matrix to start with. matrix2 = np.linalg.pinv(result @ result) - if np.linalg.norm(matrix - matrix2) > 1.0e-3: + if np.linalg.norm(matrix - matrix2) > tol: raise ValueError( f"Incurred Numerical Imprecision {np.linalg.norm(matrix-matrix2)= :.3f}" ) @@ -106,13 +112,23 @@ class RadialBasis: """ def __init__( - self, radial_basis, max_angular, cutoff_radius, max_radial=None, **hypers + self, + radial_basis, + max_angular, + cutoff_radius, + max_radial=None, + rcond=1e-8, + tol=1e-3, + **hypers, ): # Store all inputs into internal variables self.radial_basis = radial_basis self.max_angular = max_angular self.cutoff_radius = cutoff_radius self.hypers = hypers + self.rcond = rcond + self.tol = tol + if self.radial_basis not in ["monomial", "gto"]: raise ValueError(f"{self.radial_basis} is not an implemented basis.") @@ -247,7 +263,7 @@ def orthonormalize_basis(self, features: TensorMap): gto_overlap_matrix_slice = self.overlap_matrix[l_2n_arr, :][:, l_2n_arr] orthonormalization_matrix = inverse_matrix_sqrt( - gto_overlap_matrix_slice + gto_overlap_matrix_slice, self.rcond, self.tol ) block.values[:, :, neighbor_mask] = np.einsum( "ijk,kl->ijl", From 9c43999122e20a8473c9c20a502dab05b6bb70d8 Mon Sep 17 00:00:00 2001 From: arthur-lin1027 <35580059+arthur-lin1027@users.noreply.github.com> Date: Fri, 19 Jan 2024 10:40:07 -0600 Subject: [PATCH 4/4] Update tests.yml Removed 3.7 from python versions to test. --- .github/workflows/tests.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index a12d1a7..dfa27ed 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -11,7 +11,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: ["3.7", "3.8", "3.9", "3.10"] + python-version: ["3.8", "3.9", "3.10"] steps: - uses: actions/checkout@v2