python – Symbolic Regression on domain of matrices (using summation operator)

I am attempting to get a python SR library to derive a formula similar to the Lennard-Jones potential, which features nested sigma operators (requiring for a 2D array of inputs for a single output):

Lennard-Jones potential format

My input and output are formatted as the following for 2000 configurations of a 40-atom system:

Note that “r” in the above equation represents one of the entries in the 40×40 matrix below

X: (ndarray(shape=(40,40))

len(X) = 2000

y: (float64)

len(y) = 2000

My question is whether there exist any python SR libraries which support this input format and feature summation operators. The one I’m using currently ( only features basic arithmetic operators on scalar values.

Alternatively, is there a way I could reformat this data to conform to existing libraries’ requirements?