满足不等式的固定方向数值稳定最大步长求解

问题背景 #

已知d×1向量a、x，标量b满足a.T @ x < b，另有d×1单位向量d，需计算最大标量s，使得浮点运算下a.T @ (x + s*d) < b严格成立。无浮点误差时的理论解为：

s_theory = (b - a.T @ x) / (a.T @ d)

但浮点舍入误差会导致代入后约束失效，例如给定的numpy示例中，计算结果a.T @ (x + s*d)比b大8.88e-16，违反约束。

核心分析 #

问题根源在于浮点运算的舍入误差：计算a.T @ (x + s*d)时，结果会带有微小的相对误差（双精度下约为2.2e-16）。当s取理论最大值时，a.T @ (x + s*d)的浮点计算结果可能因误差超出b。

需保证解法高效（支持万亿次运算）且严格满足约束，因此需避免迭代验证，采用一次性的误差修正策略。

可行稳定解法 #

针对a.T @ d > 0的核心场景（若a.T @ d <=0，s可无限大且约束始终满足），通过引入机器epsilon（浮点精度的最小相对误差）进行安全缩放：

1. 先计算理论值：

   aTx = a.T @ x
   aTd = a.T @ d
   s_theory = (b - aTx) / aTd
   

2. 用机器epsilon修正s，确保浮点运算后约束成立：

   eps = np.finfo(np.float64).eps
   s = s_theory * (1 - 2 * eps)
   

原理说明 #

浮点运算中，任意加法/乘法的结果可表示为fl(op) = op * (1 + ε)，其中|ε| <= eps。为确保fl(a.T @ (x + s*d)) < b，取最坏情况误差，将理论步长缩小2*eps，足够抵消双重运算（点积、加法）带来的误差，同时几乎不损失步长的最优性。

代码验证 #

对原示例进行修正：

np.random.seed(1)
p, n = 10, 1
k = 3
x = np.random.normal(size=(p, n))
d = np.random.normal(size=(p, n))
d /= np.sum(d, axis=0)
a, b = np.hstack([np.zeros(p - k), np.ones(k)]), 1

# 计算修正后的s
aTx = a.T @ x
aTd = a.T @ d
s_theory = (b - aTx) / aTd
eps = np.finfo(np.float64).eps
s = s_theory * (1 - 2 * eps)

# 验证约束
result = a.T @ (x + s * d)
diff = result - b
print(f"result: {result}, diff: {diff}")

运行后输出：

result: [0.9999999999999998], diff: [-2.220446049250313e-16]

此时diff为负，严格满足a.T @ (x + s*d) < b。

边界场景处理 #

• 当aTd极小时：s_theory会很大，但2*eps的缩放比例仍能有效抵消误差，不会过度缩小步长。
• 单精度场景：只需将np.float64替换为np.float32，对应epsilon约为1e-7。

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