SageMath/Python使用科学计数法导致Ricci曲率判断异常问题

问题描述 #

我正在编写检测图是否为Ricci平坦（Ricci-flat）的程序，核心逻辑为：计算图中每条边的Ricci曲率，若所有边曲率均为0则判定图为Ricci平坦，反之则不是。受数值误差影响，实际为0的曲率有时会被计算为约10^{-16的极小值，因此程序设置判定规则：每条边Ricci曲率的绝对值小于10}-10即判定为0。

简化实现逻辑 #

# 输入图G，若为Ricci平坦返回True，否则返回False
def ricci_flat(G):
    for i in G.edges(): # 遍历图中所有边
        curvature = [计算边i的Ricci曲率]
        # 只要存在边曲率非零，图就不是Ricci平坦
        if abs(curvature) >= 10**(-10): return(False)
    # 所有边都满足曲率为0判定条件，返回True
    return(True)

异常现象 #

我发现函数的输出结果取决于判断语句中10^-10阈值的表示形式：

• 若使用10**(-10)或10^(-10)作为阈值，即使边的Ricci曲率实际为0（return前通过print(curvature)可验证），函数也总会在第一条边就错误返回False
• 若使用0.0000000001或1/(10**10)作为阈值，函数运行正常，仅当边曲率大于10^-10时返回False
  我可以用小数表示法临时解决问题，但想知道该异常出现的原因。

完整代码参考 #

检测函数代码

import numpy as np
from sage.all import *

def ricci_flat(G):   # 逐边计算图G的Ricci曲率
    # 若图为OLLY-Ricci平坦返回True，否则返回False
    # 提取图信息
    Adjacency=G.adjacency_matrix()
    Vertices=G.vertices()
    Degrees=G.degree()
    Edges=G.edges()
    Neighbors=[G.neighbors(i) for i in Vertices]
    n=len(Vertices)
    ############################################
    # 定义距离与邻居函数
    def d(x, y):
        return G.distance(x, y)
    # 定义alpha函数
    def alpha(a, b):
        return 1/(max(Degrees[a], Degrees[b])+1)
    # 定义概率分布函数
    def mu(a,b,z):
        if z==a: return alpha(a,b)
        elif z in Neighbors[a]: return (1-alpha(a,b))/Degrees[a]
        else: return 0
    ############################################
    # 计算OLLY Ricci平坦性
    for i,j,_ in Edges:
        # 设置优化方法
        p = MixedIntegerLinearProgram(maximization=False)
        A = p.new_variable(real=True) # A[x, y]为变换函数
        p.set_objective(sum(A[x, y]*d(x, y) for x in range(n) for y in range(n)))
        # 添加约束
        for x in range(n):
            p.add_constraint(sum(A[x, y] for y in range(n))==mu(i,j,x)) # 行和约束
            p.add_constraint(sum(A[y, x] for y in range(n))==mu(j,i,x)) # 列和约束
            for y in range(n): p.add_constraint(0<=A[x, y]) # 所有变量非负
        # 计算边的OLLY曲率
        curvature = 1-p.solve()
        # 调试时输出显示，曲率要么为0，要么为小于10^-15的极小值
        if abs(curvature) >= 10**(-10): return(False)
    # 所有边均满足Ricci平坦判定条件，返回True
    return(True)

测试用例构造代码

以下是已知的Ricci平坦图构造代码：若阈值使用1/(10**10)或0.0000000001，ricci_flat(A1A1)会正确返回True；若替换为10**(-10)或10^(-10)，则会错误返回False。

import numpy as np
from sage.all import *
n = 6
v = n**2
adj = np.zeros((v,v)) # 图的邻接矩阵
for i in range(v):
    if i<v-n:
        adj[i,i+n]=1
        adj[i+n,i]=1
    if i<n:
        adj[i,i+v-n]=1
        adj[i+v-n,i]=1
    if i%n!=0:
        adj[i,i-1]=1
        adj[i-1,i]=1
    else:
        adj[i,i+n-1]=1
        adj[i+n-1,i]=1
mat = matrix(adj).change_ring(ZZ)
A1A1 = Graph(mat)

补充说明 #

更换其他Ricci平坦图测试该问题也会复现，且在Sage shell中输入10**(-10) == 1/(10**10)会返回预期的True。 #

问题解答 #

核心原因 #

这个异常是Sage环境特有的类型系统特性导致的：

1. MixedIntegerLinearProgram.solve()返回的结果默认是符号表达式类型，而非直接的数值类型。当你用abs(curvature) >= 阈值做比较时：
   • 如果阈值是0.0000000001这类浮点数，Sage会自动将符号表达式转换为实际数值再做比较，逻辑正常。
   • 如果阈值是10**(-10)这类精确有理数类型，Sage不会触发自动数值转换，而是返回一个未求值的符号不等式。所有非布尔类型的对象在Python的if判断中都会被视为True，因此不管曲率实际值是多少，判断条件都会触发，直接返回False。
2. 你在Sage shell中测试10**(-10) == 1/(10**10)返回True，是因为两个都是有理数对象，直接做有理数相等比较，没有涉及符号表达式的转换逻辑。

验证方法 #

你可以在判断语句前加入以下打印代码验证类型：

print(type(curvature))
print(type(10**(-10)))
print(type(0.0000000001))
print(abs(curvature) >= 10**(-10)) # 这里会输出符号不等式而非布尔值

解决方案 #

有两种稳妥的修正方式：

• 方式1：在计算曲率后主动转换为数值类型，比如修改为curvature = float(1 - p.solve())
• 方式2：阈值统一使用浮点数表示，比如直接写1e-10

内容的提问来源于stack exchange，提问作者rtg142857