如何为自定义Python符号回归工具确定必要算子？

我用Python开发了一个符号回归工具，支持指定一元（unary）和二元（binary）算子。当未指定算子时，希望能自动确定潜在的必要算子。尝试过机器学习和相关性分析方法，但针对2*log(x1)+3*x2-x1*x2+5和exp(x1)这两个函数定义，所得效果并不理想。曾考虑一种思路：先基于相关性最高的项进行回归，移除主成分后再次迭代，但认为该思路对cos(log)这类复合函数可能无效。

实现代码 #

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestRegressor
from sklearn.preprocessing import FunctionTransformer
from sklearn.model_selection import train_test_split
import operator
import sympy

np.random.seed(42)
n = 1000
x1 = np.random.uniform(1, 10, size = n)
x2 = np.random.uniform(1, 10, size = n)
y = 2 * np.log(x1) + 3 * x2 - x1 * x2 + 5
#y = np.exp(x1)

unary_operators = {"neg": (lambda x: sympy.sympify("neg(" + str(x) + ")"), operator.neg),
                   "abs": (sympy.Abs, operator.abs),
                   "inv_": (lambda x: sympy.sympify("inv_(" + str(x) + ")"), lambda x: 1 / x),
                   "sqrt": (sympy.sqrt, np.sqrt),
                   "cos": (sympy.cos, np.cos),
                   "sin": (sympy.sin, np.sin),
                   "exp": (sympy.tan, np.tan),
                   "log": (sympy.log, np.log),
                   "exp": (sympy.exp, np.exp),
                   "sinh": (sympy.sinh, np.sinh),
                   "cosh": (sympy.cosh, np.cosh),
                   "floor": (lambda x: sympy.sympify("floor(" + str(x) + ")"), np.ceil),
                   "ceil": (lambda x: sympy.sympify("ceil(" + str(x) + ")"), np.floor)}

binary_operators = {"+": (operator.add, operator.add),
                    "-": (operator.sub, operator.sub),
                    "*": (operator.mul, operator.mul),
                    "/": (operator.truediv, operator.truediv),
                    "//": (operator.floordiv, operator.floordiv),
                    "%": (operator.mod, operator.mod),
                    "**": (sympy.Pow, operator.pow)}

symmetric_binary_operators = {"+": True, "-": True, "*": False, "conv": False}
symbols = sympy.symbols("x1 x2")
X = [x1, x2]

X_raw = pd.DataFrame()
        
for i in range(0, len(symbols)):
    X_raw[str(symbols[i])] = X[i]

feature_dict = {}
base_vars = [str(x) for x in symbols]
ops = {}

for k, v in unary_operators.items():
    sym_op, num_op = v

    for var in symbols:
        feature_dict[str(sym_op(var))] = num_op(X_raw[str(var)])
        ops[str(sym_op(var))] = k

for k, v in binary_operators.items():
    sym_op, num_op = v
    
    indices1 = list(range(0, len(symbols)))

    for i1 in indices1:
        indices2 = list(range(0, len(symbols)))

        if (k in list(symmetric_binary_operators.keys())):
            indices2 = list(range(i1 + 1 if v else i1, len(base_vars)))
                        
        for i2 in indices2:
            feature_dict[str(sym_op(symbols[i1], symbols[i2]))] = num_op(X_raw[str(symbols[i1])], X_raw[str(symbols[i2])])
            ops[str(sym_op(symbols[i1], symbols[i2]))] = k

X_feat = pd.DataFrame(feature_dict)
X_feat = X_feat.replace([np.inf, -np.inf, np.nan], 0)
X_feat = X_feat.loc[:, (X_feat.abs().max() < 1e6)]

correlations = X_feat.apply(lambda col: np.corrcoef(col, y)[0, 1])

# 按绝对值降序排序
correlations_sorted = correlations.abs().sort_values(ascending = False)

print("Top 10 相关性最高的特征（Pearson）:\n")
for name in correlations_sorted.head(10).index:
    print(f"{name:30}  相关系数: {correlations[name]:+.4f}")

from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_feat_scaled_array = scaler.fit_transform(X_feat)
X_feat = pd.DataFrame(X_feat_scaled_array, columns=X_feat.columns)

from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X_feat, y, test_size = 0.2, random_state = 0)

model = RandomForestRegressor(n_estimators = 200, random_state = 0)#, max_depth = 10)
model.fit(X_train, y_train)

importances = model.feature_importances_

features_sorted = sorted(zip(X_feat.columns, importances), key = lambda x: x[1], reverse = True)

imp = 0

for i in range(0, len(features_sorted)):
    features_sorted[i] = list(features_sorted[i])
    x = features_sorted[i]
    imp_ = x[1]
    x.append(abs(x[1] - imp))
    imp = imp_

print("\nTop 10 重要性最高的特征:\n")
for name, score, r in features_sorted[:10]:
    print(f"{name:30}  重要性: {score:.4f} {r:.4f}")

min_score = 0.015
min_r = 0.05
un_ops = set()
bin_ops = set()

for name, score, r in features_sorted:
    #if (score > min_score):
    if (r > min_r):
        if (ops[name] in unary_operators):
            un_ops.add(ops[name])
        elif (ops[name] in binary_operators):
            bin_ops.add(ops[name])  

print("\n识别出的一元算子:", un_ops)
print("识别出的二元算子:", bin_ops)

运行输出 #

Top 10 相关性最高的特征（Pearson）:

x1*x2                           相关系数: -0.9480
x1 + x2                         相关系数: -0.8527
neg(x1)                         相关系数: +0.8298
Abs(x1)                         相关系数: -0.8298
ceil(x1)                        相关系数: -0.8244
floor(x1)                       相关系数: -0.8244
sqrt(x1)                        相关系数: -0.8170
log(x1)                         相关系数: -0.7880
inv_(x1)                        相关系数: +0.6826
sinh(x1)                        相关系数: -0.6305

Top 10 重要性最高的特征:

x1*x2                           重要性: 0.7583 0.7583
x1 + x2                         重要性: 0.0850 0.6734
log(x1)                         重要性: 0.0192 0.0658
inv_(x1)                        重要性: 0.0170 0.0022
sinh(x1)                        重要性: 0.0164 0.0007
neg(x1)                         重要性: 0.0162 0.0002
sqrt(x1)                        重要性: 0.0159 0.0003
exp(x1)                         重要性: 0.0159 0.0001
cosh(x1)                        重要性: 0.0151 0.0008
Abs(x1)                         重要性: 0.0137 0.0014

识别出的一元算子: {'log'}
识别出的二元算子: {'*', '+'}

内容的提问来源于stack exchange，提问作者P'tit Ju