5 incorporate normalization factors

[arthur-lin1027]

Starting a PR so I can show what each normalization method does.

[arthur-lin1027]

Commit 3ad3c6c normalizes post-hoc and looks nice, but likely isn't correct because it only contains a normalization factor consisting only of the gto's gaussian width, not of the 3 axes of our anisotropic gaussian.

This results in a anisoap vector (at the end of doing the tensor product) that looks like this:

To obtain a normalization factor consisting of all three gaussian widths, I normalized within the pairwise expansion, as seen in this commit: cdb48a2

This results in a anisoap vector (at the end of doing the tensor product) that looks like this:

I'd love to hear your comments on how these look!

[arthur-lin1027]

The normalization code is ready for review. It is pretty much the same as 3ad3c6c

• The normalization constant is taken from equation 10 in this paper, but what's giving me some doubts is the fact that $\sigman$ is coupled to $n$ and $n{max}$ in the paper, while I believe no such coupling occurs in our code, i.e. $\sigma_n$ is a tunable hyperparameter. I'm also unsure what $Nn$ does, it certainly scales the function down, but I'm not sure if it ensures $\int{r=0}^{\infty} r^n e^{-r^2/2\sigman^2} dr= 1$, $\int{r=0}^{\infty} (r^n e^{-r^2/2\sigma_n^2})^2 dr = 1$, or if the integral is over $R^3$, etc.

  • Maybe these considerations don't actually matter too much.

• I've included scratch code here that reads in the benzene files, and compares the normalized vs unnormalized representations.

• Am also open to feedback on what kinds of tests I should include. I manually checked that contract(normalize(pairwise_features)) == normalize(contract(pairwise_features)), should I write that as a separate test? I see that there aren't individual tests for pairwise_ellip_expansion and contract_pairwise_feat so I wasn't sure if I should include one for the normalization function.

[arthur-lin1027]

I've completed my first iteration of orthonormalization. Orthonormalization functions are now all in the radial_basis.py function. The orthonormalized results perform well! Based off some rudimentary tests the orthonormalized results generally outperform the unnormalized results, and seems to be less prone to overfitting.

• First, super important to check the math. I'm pretty confident in the math but it's good to get a second pair of eyes. I've taken a slightly different approach to librascal by calculating overlap with $n{eff}$ and $\sigma{eff}$ but it should be equivalent (and to me, more intuitive).
• I'm super welcome to any API changes, e.g. maybe these new methods that I added should not be part of the class since they are largely static...?
• If any docstrings are confusing then let me know as well

[rosecers]

I could go either way — if you keep the checks, they should have coinciding tests.

On 28 Jun 2023, at 17:35, arthur-lin1027 @.***> wrote:

@arthur-lin1027 commented on this pull request.

In anisoap/representations/radial_basis.py https://github.com/cersonsky-lab/anisoap/pull/6#discussion_r1245891014:

• if not np.allclose(matrix, matrix.T):
• raise ValueError("Matrix is not hermitian")
• eva, eve = np.linalg.eigh(matrix)
• if (eva < 0).all():
• raise ValueError(
• "Matrix is not positive semidefinite. Check that a valid gram matrix is passed."
• ) Hmm in its current usage these checks will not be triggered. However, these are general purpose util functions that could theoretically be (mis)used elsewhere.

So, should it be the function's role to verify that the matrix is symmetric positive semidefinite, or should such a verification be done outside the function?

— Reply to this email directly, view it on GitHub https://github.com/cersonsky-lab/anisoap/pull/6#discussion_r1245891014, or unsubscribe https://github.com/notifications/unsubscribe-auth/ALKVP3WANDINOYTJ35O2BD3XNSWUXANCNFSM6AAAAAAY6NPPWE. You are receiving this because your review was requested.

[arthur-lin1027]

I don't know how to test for orthonormality... the numeric integration isn't great LOL. But I've been able to test for everything else.

[arthur-lin1027]

Unfortunately this code does not pass when we are dealing with large number of lmax (I'm working with lmax=15, which is what is needed for the gay berne ellipsoids). For whatever reason, our gram matrix is no longer positive semidefinite and has negative eigenvalues, which I find surprising -- I'm suspecting this is a numerical stability issue, as when I perform the eigendecomposition on this matrix, and perform eve @ eva @ eve.T, I do not obtain my original matrix

[arthur-lin1027]

With this latest commit I made on June 30, I first normalize the values in the tensor block, then perform orthonormalization using an overlap matrix consisting of the overlaps of normalized GTOs.

Previously, I was using an overlap matrix consisting of the overlaps of nonnormalized GTOs, which led to numerical stability issues. The code should now be ready to merge in pending a final review.