$\Phi$-SO：物理符号优化

物理符号优化（$\Phi$-SO）—— 一个专为物理学设计的符号优化工具包。

源代码：WassimTenachi/PhySO
文档：physo.readthedocs.io

最新动态 ✨

2025-08 ：📦 现可通过 pip install physo 和 conda 安装！
2025-07 ：🐍 支持 Python 3.12 及最新版本的 NumPy、PyTorch 和 SymPy。
2024-06 ：📚 全面更新文档。
2024-05 ：🔬 Class SR：多数据集符号回归。
2024-02 ：🎯 考虑不确定性的拟合。
2023-08 ：⚡ 加速维度分析。
2023-03 ：🌟 PhySO 初始发布（专注于物理学的符号回归）。

亮点

$\Phi$-SO 的符号回归模块利用深度强化学习，从函数形式空间中搜索并推导出能够拟合数据点的解析物理定律。

physo 能够利用：

• 物理单位约束，通过维度分析缩小搜索空间（[Tenachi 等, 2023]）

• 类约束，寻找一个能够准确拟合多个数据集的单一解析函数形式——每个数据集可能由自己独特的拟合参数集所支配（[Tenachi 等, 2024]）

$\Phi$-SO 恢复阻尼谐振子方程的过程：

https://github.com/WassimTenachi/PhySO/assets/63928316/655b0eea-70ba-4975-8a80-00553a6e2786

在来自 SRBench 的标准费曼基准测试中的表现，该基准包含《费曼物理学讲义》中的 120 个表达式，并与流行的符号回归软件进行了对比。

$\Phi$-SO 在噪声存在时（超过 0.1%）达到了最先进的性能，即使在高达 10% 的显著噪声下也能保持稳健的表现：

安装

该软件包已在以下平台上测试通过：

• Linux
• OSX（ARM 和 Intel）
• Windows

如果您在安装过程中遇到问题，建议您尝试从源码安装（源码安装指南）。如果问题仍然存在，请参阅常见问题解答（FAQ），或在 GitHub 仓库 中提交问题。

使用 pip 安装

从 PyPI 安装 physo：

pip install physo

使用 conda 安装

使用 conda 安装 physo：

conda install -c conda-forge physo

快速入门（SR）

本教程将展示如何使用 physo 进行符号回归（SR）。本教程的参考笔记本可在以下链接找到：📙sr_quick_start.ipynb。

准备工作

导入必要的库：

# 外部库
import numpy as np
import matplotlib.pyplot as plt
import torch

导入 physo：

# 内部代码
import physo
import physo.learn.monitoring as monitoring

为了确保结果的可重复性，建议固定随机种子：

# 随机种子
seed = 0
np.random.seed(seed)
torch.manual_seed(seed)

创建合成数据集

创建一个玩具合成数据集：

# 创建玩具合成数据
z = np.random.uniform(-10, 10, 50)
v = np.random.uniform(-10, 10, 50)
X = np.stack((z, v), axis=0)
y = 1.234*9.807*z + 1.234*v**2

需要注意的是，自由常数的搜索默认从 1 开始。因此，在使用默认超参数时，强烈建议将数据归一化到数量级为 1 左右。

---

DA 补充说明：
$\Phi$-SO 可以利用 DA（维度分析）来提高符号回归的效率。可以考虑 $X=(z,v)$ 的物理单位：$z$ 是长度，其维度为 $L^{1}, T^{0}, M^{0}$；$v$ 是速度，其维度为 $L^{1}, T^{-1}, M^{0}$；而 $y=E$ 是能量，其维度为 $L^{2}, T^{-2}, M^{1}$。
如果您处理的问题并非物理相关，且所有变量和常数均为无量纲，则无需指定任何 xx_units 参数（或将所有变量/常数的单位设置为 [0,0]），此时 physo 将执行无量纲的符号回归任务。

---

绘制数据集：

n_dim = X.shape[0]
fig, ax = plt.subplots(n_dim, 1, figsize=(10,5))
for i in range (n_dim):
    curr_ax = ax if n_dim==1 else ax[i]
    curr_ax.plot(X[i], y, 'k.',)
    curr_ax.set_xlabel("X[%i]"%(i))
    curr_ax.set_ylabel("y")
plt.show()

SR 配置

需要注意的是，physo 的 SR 能力高度依赖于超参数，因此建议根据具体问题调整超参数以进行科学研究。目前可用的超参数预设配置总结如下：

配置	推荐用例	速度	效果	备注
config0	演示	★★★	★	轻量且快速的配置。
config1	带有量纲分析 $^$ 的 SR；带有量纲分析 $^$ 的分类 SR	★	★★★	用于费曼基准和 MW 流基准的配置。
config2	SR；分类 SR	★★	★★	用于分类基准的配置。

$^*$ DA = 量纲分析

鼓励用户编辑配置文件（可在：physo/config/ 中找到）。默认使用 config0，但为了科学研究，建议遵循上述推荐配置。

---

关于量纲分析的补充说明:

1. 在前几十次迭代中，神经网络通常仍在学习量纲分析的规则，导致大多数候选表达式被丢弃而未被训练，从而有效批次大小显著减小（通常缩小 10 倍），使得评估过程的计算开销大大降低。因此，建议通过使用较大的批次大小来弥补这一现象，以便为神经网络提供足够的学习信息。

---

日志与可视化设置：

save_path_training_curves = 'demo_curves.png'
save_path_log             = 'demo.log'

run_logger     = lambda : monitoring.RunLogger(save_path = save_path_log,
                                                do_save = True)

run_visualiser = lambda : monitoring.RunVisualiser (epoch_refresh_rate = 1,
                                           save_path = save_path_training_curves,
                                           do_show   = False,
                                           do_prints = True,
                                           do_save   = True, )

运行 SR

给定变量数据 $(x_0,..., x_n)$（此处为 $(z, v)$），根变量 $y$（此处为 $E$）以及自由常数和固定常数，您可以通过以下命令运行 SR 任务以恢复 $f$。

---

关于量纲分析的补充说明:
在此我们允许使用一个维度为 $L^{0}, T^{0}, M^{0}$（即无量纲）的固定常数 $1$，以及维度分别为 $L^{0}, T^{0}, M^{1}$ 的自由常数 $m$ 和 $L^{1}, T^{-2}, M^{0}$ 的自由常数 $g$。需要注意的是，这里的单位向量长度为 3（例如：[1, 0, 0]），因为本示例中的变量仅涉及长度、时间与质量三个物理量。然而，只要在 $X$、$y$ 和常数之间保持一致，单位向量的长度可以是任意不超过 7 的值，从而允许用户表示任何单位（包括长度、时间、质量、温度、电流、光通量或物质的量等）。此外，无论单位的顺序如何，都可以进行量纲分析，因此用户可以采用任何约定（如 [长度、质量、时间] 或 [质量、时间、长度] 等），只要在 $X$、$y$ 和常数之间保持一致即可。

---

# 运行 SR 任务
expression, logs = physo.SR(X, y,
                            # 给出变量名称（用于显示）
                            X_names = [ "z"       , "v"        ],
                            # 相应的物理单位（若不相关可忽略或传入零向量）
                            X_units = [ [1, 0, 0] , [1, -1, 0] ],
                            # 给出根变量名称（用于显示）
                            y_name  = "E",
                            y_units = [2, -2, 1],
                            # 固定常数
                            fixed_consts       = [ 1.      ],
                            fixed_consts_units = [ [0,0,0] ],
                            # 共享常数名称（用于显示）
                            free_consts_names = [ "m"       , "g"        ],
                            free_consts_units = [ [0, 0, 1] , [1, -2, 0] ],
                            # 可用于构建 f 的符号运算
                            op_names = ["mul", "add", "sub", "div", "inv", "n2", "sqrt", "neg", "exp", "log", "sin", "cos"],
                            get_run_logger     = run_logger,
                            get_run_visualiser = run_visualiser,
                            # 运行配置
                            run_config = physo.config.config0.config0,
                            # 并行模式（仅在从 Python 脚本运行时可用，笔记本中不可用）
                            parallel_mode = False,
                            # 迭代次数
                            epochs = 20
)

检查找到的最佳表达式

获取最佳表达式：

在准确度方面找到的最佳表达式会返回在 expression 变量中：

best_expr = expression
print(best_expr.get_infix_pretty())

>>> 
                     2           
    -g⋅m⋅z + -v⋅v⋅sin (1.0)⋅1.0⋅m

该表达式也可以稍后从日志文件中加载：

import physo
from physo.benchmark.utils import symbolic_utils as su
import sympy

# 加载帕累托前沿表达式
pareto_expressions = physo.read_pareto_pkl("demo_curves_pareto.pkl")

# 最准确的表达式是帕累托前沿中的最后一个：
best_expr = pareto_expressions[-1]
print(best_expr.get_infix_pretty())

显示：

该表达式可以转换为……
一个 SymPy 表达式：

best_expr.get_infix_sympy()

>>> -g*m*z - v*v*sin(1.0)**2*1.0*m

一个 SymPy 表达式（已计算自由常数的值）：

best_expr.get_infix_sympy(evaluate_consts=True)[0]

>>> 1.74275713004454*v**2*sin(1.0)**2 + 12.1018380702846*z

一个 LaTeX 字符串：

best_expr.get_infix_latex()

>>> '\\frac{m \\left(- 1000000000000000 g z - 708073418273571 v^{2}\\right)}{1000000000000000}'

一个 LaTeX 字符串（已计算自由常数的值）：

sympy.latex(best_expr.get_infix_sympy(evaluate_consts=True))

>>> '\\mathtt{\\text{[1.74275713004454*v**2*sin(1.0)**2 + 12.1018380702846*z]}}'

获取自由常数值：

best_expr.free_consts

>>> FreeConstantsTable
     -> 类常数（['g' 'm']）：(1, 2)
     -> 特殊常数（[]）：(1, 0, 1)

best_expr.free_consts.class_values

>>> tensor([[ 6.9441, -1.7428]], dtype=torch.float64)

检查精确的符号恢复

# 转换为 SymPy
best_expr = best_expr.get_infix_sympy(evaluate_consts=True)

best_expr = best_expr[0]

# 打印简化后的最佳表达式，并对常数进行四舍五入
print("best_expr : ", su.clean_sympy_expr(best_expr, round_decimal = 4))

# 目标表达式为：
target_expr = sympy.parse_expr("1.234*9.807*z + 1.234*v**2")
print("target_expr : ", su.clean_sympy_expr(target_expr, round_decimal = 4))

# 检查等价性
print("\n检查等价性：")
is_equivalent, log = su.compare_expression(
                        trial_expr  = best_expr,
                        target_expr = target_expr,
                        handle_trigo            = True,
                        prevent_zero_frac       = True,
                        prevent_inf_equivalence = True,
                        verbose                 = True,
)
print("是否等价：", is_equivalent)

>>> best_expr :  1.234*v**2 + 12.1018*z
    target_expr :  1.234*v**2 + 12.1018*z
    
    检查等价性：
      -> 评估 1.234*v**2 + 12.101838*z（目标）是否与 1.74275713004454*v**2*sin(1.0)**2 + 12.1018380702846*z（试验）等价
       -> 简化后的表达式：1.23*v**2 + 12.1*z
       -> 符号误差：0
       -> 符号分数：1
       -> 三角函数符号误差：0
       -> 三角函数符号分数：1
       -> 是否等价：是
    是否等价：是

文档

更多文档可在 physo.readthedocs.io 查阅。

符号回归 快速入门指南：这里

分类符号回归 快速入门指南：这里

引用本工作

基于强化学习和量纲分析的符号回归

@ARTICLE{PhySO_RL_DA,
       author = {{Tenachi}, Wassim and {Ibata}, Rodrigo and {Diakogiannis}, Foivos I.},
        title = "{利用单位约束引导的深度符号回归：迈向物理定律的自动化发现}",
      journal = {ApJ},
         year = 2023,
        month = dec,
       volume = {959},
       number = {2},
          eid = {99},
        pages = {99},
          doi = {10.3847/1538-4357/ad014c},
archivePrefix = {arXiv},
       eprint = {2303.03192},
 primaryClass = {astro-ph.IM},
       adsurl = {https://ui.adsabs.harvard.edu/abs/2023ApJ...959...99T},
      adsnote = {由 SAO/NASA 天体物理数据系统提供}
}

分类符号回归

@ARTICLE{PhySO_ClassSR,
       author = {{Tenachi}, Wassim and {Ibata}, Rodrigo and {Fran{\c{c}}ois}, Thibaut L. and {Diakogiannis}, Foivos I.},
        title = "{分类符号回归：必须适应所有情况}",
      journal = {天体物理学杂志快报},
         year = {2024},
        month = {jul},
       volume = {969},
       number = {2},
          eid = {arXiv:2312.01816},
        pages = {L26},
          doi = {10.3847/2041-8213/ad5970},
archivePrefix = {arXiv},
       eprint = {2312.01816},
 primaryClass = {cs.LG},
       adsurl = {https://ui.adsabs.harvard.edu/abs/2023arXiv231201816T},
      adsnote = {由 SAO/NASA 天体物理数据系统提供}
}

PhySO 快速上手指南

PhySO ($\Phi$-SO) 是一个专为物理学设计的符号优化包。它利用深度强化学习，在考虑物理单位约束（量纲分析）和类别约束的情况下，从数据点中推断出解析形式的物理定律。

环境准备

系统要求

PhySO 已在以下操作系统上经过测试：

• Linux
• macOS (支持 ARM 和 Intel 芯片)
• Windows

前置依赖

• Python: 支持 Python 3.12 及更高版本。
• 核心库: 自动安装 NumPy, PyTorch, SymPy 等依赖。
• 建议: 为了获得最佳科学计算性能，建议使用虚拟环境（如 venv 或 conda）。

注意：由于涉及深度学习训练，推荐使用具备 CUDA 支持的 GPU 环境以加速运算，但 CPU 亦可运行。

安装步骤

你可以选择通过 pip 或 conda 进行安装。国内用户若遇到下载速度慢的问题，可配置国内镜像源。

方式一：使用 pip 安装

标准安装：

pip install physo

使用国内镜像加速安装（推荐）：

pip install physo -i https://pypi.tuna.tsinghua.edu.cn/simple

方式二：使用 conda 安装

标准安装：

conda install -c conda-forge physo

使用国内镜像加速安装（推荐）：

conda install -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge physo

基本使用

以下是一个最简单的符号回归（Symbolic Regression, SR）示例，演示如何从合成数据中恢复物理公式。

1. 导入库与设置种子

import numpy as np
import torch
import physo
import physo.learn.monitoring as monitoring

# 设置随机种子以保证结果可复现
seed = 0
np.random.seed(seed)
torch.manual_seed(seed)

2. 构建合成数据集

我们构造一个符合物理规律 $E = m \cdot g \cdot z + m \cdot v^2$ 的数据集。 注意：默认情况下自由常数的搜索范围在 1 附近，因此建议将数据归一化到数量级为 1 的范围。

# 生成随机数据
z = np.random.uniform(-10, 10, 50)
v = np.random.uniform(-10, 10, 50)
X = np.stack((z, v), axis=0)

# 目标公式：y = 1.234 * 9.807 * z + 1.234 * v^2
# 这里隐含了常数 m=1.234, g=9.807
y = 1.234*9.807*z + 1.234*v**2

3. 配置日志与可视化

save_path_training_curves = 'demo_curves.png'
save_path_log             = 'demo.log'

run_logger     = lambda : monitoring.RunLogger(save_path = save_path_log, do_save = True)
run_visualiser = lambda : monitoring.RunVisualiser(epoch_refresh_rate = 1,
                                           save_path = save_path_training_curves,
                                           do_show   = False,
                                           do_prints = True,
                                           do_save   = True)

4. 运行符号回归

在此示例中，我们利用**量纲分析（Dimensional Analysis）**来辅助搜索：

• 输入变量 $z$ (长度) 和 $v$ (速度)。
• 输出变量 $y$ (能量)。
• 定义固定常数 $1$ (无量纲) 和自由常数 $m$ (质量), $g$ (加速度)。

expression, logs = physo.SR(X, y,
                            # 变量名称
                            X_names = ["z", "v"],
                            # 变量物理单位 [长度 L, 时间 T, 质量 M]
                            # z: 长度 [1, 0, 0], v: 速度 [1, -1, 0]
                            X_units = [[1, 0, 0], [1, -1, 0]],
                            
                            # 根变量名称
                            y_name  = "E",
                            # 根变量单位 (能量): [2, -2, 1]
                            y_units = [2, -2, 1],
                            
                            # 固定常数及其单位 (无量纲)
                            fixed_consts       = [1.],
                            fixed_consts_units = [[0, 0, 0]],
                            
                            # 自由常数名称及其单位
                            # m: 质量 [0, 0, 1], g: 加速度 [1, -2, 0]
                            free_consts_names = ["m", "g"],
                            free_consts_units = [[0, 0, 1], [1, -2, 0]],
                            
                            # 可用的数学操作符
                            op_names = ["mul", "add", "sub", "div", "inv", "n2", "sqrt", "neg", "exp", "log", "sin", "cos"],
                            
                            get_run_logger     = run_logger,
                            get_run_visualiser = run_visualiser,
                            
                            # 运行配置：科研用途建议使用 config1 或 config2，此处演示使用轻量级的 config0
                            run_config = physo.config.config0.config0,
                            
                            parallel_mode = False,
                            epochs = 20  # 迭代次数
)

5. 查看结果

获取并打印找到的最佳表达式：

best_expr = expression
print(best_expr.get_infix_pretty())

输出示例：

                     2           
    -g⋅m⋅z + -v⋅v⋅sin (1.0)⋅1.0⋅m

你也可以将其转换为标准的 SymPy 表达式并计算出具体的常数值：

# 获取带具体数值的 SymPy 表达式
sympy_expr = best_expr.get_infix_sympy(evaluate_consts=True)[0]
print(sympy_expr)

输出示例：

1.74275713004454*v**2*sin(1.0)**2 + 12.1018380702846*z

(注：由于迭代次数较少且存在随机性，具体数值可能与理论值 $m \approx 1.234, g \approx 9.807$ 略有偏差，增加 epochs 或使用 config1 可提高精度)