为python中的列表定义数值稳定的sigmoid函数的最佳方法

Posted 2023-03-12

【中文标题】为python中的列表定义数值稳定的sigmoid函数的最佳方法

【英文标题】：optimal way of defining a numerically stable sigmoid function for a list in python

【发布时间】：2019-01-29 06:45:59

【问题描述】：

对于一个标量变量x，我们知道如何在python中写下一个数值稳定的sigmoid函数：

def sigmoid(x):
    if x >= 0:
        return 1. / ( 1. + np.exp(-x) )
    else:
        return exp(x) / ( 1. + np.exp(x) )

对于一个标量列表，比如z = [x_1, x_2, x_3, ...]，假设我们事先不知道每个x_i 的符号，我们可以推广上述定义并尝试：

def sigmoid(z):
    result = []
    for x in z:
        if x >= 0:
            result.append(1. / ( 1. + np.exp(-x) ) )
        else:
            result.append( exp(x) / ( 1. + np.exp(x) ) )
    return result

这似乎有效。但是，我觉得这可能不是最 Pythonic 的方式。我应该如何改进“清洁度”的定义？话说，有没有办法使用推导来缩短函数定义？

如果有人问过这个问题，我很抱歉，因为我在 SO 上找不到类似的问题。非常感谢您的宝贵时间和帮助！

【问题讨论】：

它曾经像scipy.special.expit(x) 一样简单，但在2016 年是someone changed it。

不过，我不确定它的重要性 - 我可以找到一些案例，其中所谓的更稳定的版本大约更准确 1 ULP，但我找不到可靠的分析来证明它是更好，the closest thing to a source 我发现分析溢出的行为错误，得出的错误结论是直接的 1/(1+exp(-x)) 实现将返回 x=-100 的无穷大。

看起来“稳定”版本的优势可能是一个神话，或者最多只是一个微小的精度优势。

@user2357112supportsMonica Infinity 将返回 ~-710 与幼稚实现（numpy，double 精度）

【参考方案1】：

你是对的，你可以通过使用np.where 做得更好，相当于if：

def sigmoid(x):
    return np.where(x >= 0, 
                    1 / (1 + np.exp(-x)), 
                    np.exp(x) / (1 + np.exp(x)))

这个函数接受一个 numpy 数组 x 并返回一个 numpy 数组：

data = np.arange(-5,5)
sigmoid(data)
#array([0.00669285, 0.01798621, 0.04742587, 0.11920292, 0.26894142,
#       0.5       , 0.73105858, 0.88079708, 0.95257413, 0.98201379])

【讨论】：

感谢您的快速回复。这是一个可爱的解决方案。请问是否有可能摆脱元素明智的 for 循环，我们需要一次检查 x 的符号？

我不关注你。建议的解决方案中没有元素循环。

好的，当我将 python 列表 z 传递给您的定义时，我犯了一个错误，因此引发了异常。所以我虽然我仍然需要原始定义中的 for 循环。传递 np.array(z) 有效。我会将问题标记为已解决。非常感谢！

如果要传递一个列表，只需在传递之前将其转换为数组即可（使用np.array(mylist)）。

似乎np.where 评估两个分支，然后选择它需要的分支，这会导致误导性溢出警告。类似sigmoid(np.array(-300, np.float32))

【参考方案2】：

def sigmoid(x):
    """
    A numerically stable version of the logistic sigmoid function.
    """
    pos_mask = (x >= 0)
    neg_mask = (x < 0)
    z = np.zeros_like(x)
    z[pos_mask] = np.exp(-x[pos_mask])
    z[neg_mask] = np.exp(x[neg_mask])
    top = np.ones_like(x)
    top[neg_mask] = z[neg_mask]
    return top / (1 + z)

这段代码来自cs231n的assignment3，我不太明白为什么要这样计算，但我知道这可能是你要找的代码。希望对您有所帮助。

【参考方案3】：

@hao peng 提供了完全正确的答案（没有警告），但没有清楚地解释解决方案。评论太长了，所以我会去找答案。

我们先分析几个答案（纯numpy答案而已）：

@DYZ accepted answer

这个在数学上是正确的，但仍然给我们一个警告。我们来看代码：

def sigmoid(x):
    return np.where(
            x >= 0, # condition
            1 / (1 + np.exp(-x)), # For positive values
            np.exp(x) / (1 + np.exp(x)) # For negative values
    )

由于两个分支都被评估（它们是参数，它们必须是），第一个分支会给我们一个负值警告，第二个是正值。

虽然会发出警告，但溢出的结果不会被合并，因此结果是正确的。

缺点

对两个分支都进行了不必要的评估（操作次数是需要的两倍） 发出警告

@ynn answer

这个几乎是正确的，BUT 只适用于浮点值，见下文：

def sigmoid(x):
    return np.piecewise(
        x,
        [x > 0],
        [lambda i: 1 / (1 + np.exp(-i)), lambda i: np.exp(i) / (1 + np.exp(i))],
    )


sigmoid(np.array([0.0, 1.0]))  # [0.5 0.73105858] correct
sigmoid(np.array([0, 1]))  # [0, 0] incorrect

为什么？ @mhawke 在另一个线程中提供了更长的答案，但要点是：

似乎piecewise() 将返回值转换为相同类型 作为输入，当输入整数时，整数转换是 对结果执行，然后返回。

缺点

由于分段函数的奇怪行为，无法自动转换

改进了@hao peng答案

稳定 sigmoid 的想法来自以下事实：

如果编码正确，这两个版本在操作方面的效率相同（一个 exp 评估就足够了）。现在：

当x 为正时，e^x 会溢出 当x 为负数时，e^-x 会溢出

因此我们必须在x 等于零时进行分支。使用numpy 的掩码，我们可以通过特定的 sigmoid 实现仅转换数组的正数或负数部分。

更多点见代码cmets：

def _positive_sigmoid(x):
    return 1 / (1 + np.exp(-x))


def _negative_sigmoid(x):
    # Cache exp so you won't have to calculate it twice
    exp = np.exp(x)
    return exp / (exp + 1)


def sigmoid(x):
    positive = x >= 0
    # Boolean array inversion is faster than another comparison
    negative = ~positive

    # empty contains junk hence will be faster to allocate
    # Zeros has to zero-out the array after allocation, no need for that
    result = np.empty_like(x)
    result[positive] = _positive_sigmoid(x[positive])
    result[negative] = _negative_sigmoid(x[negative])

    return result

时间测量

结果（来自ynn 的 50 次案例测试）：

289.5070939064026 #DYZ
222.49267292022705 #ynn
230.81086134910583 #this

确实分段似乎更快（不确定原因，也许掩码和额外的掩码操作使其更慢）。

使用了以下代码：

import time

import numpy as np


def _positive_sigmoid(x):
    return 1 / (1 + np.exp(-x))


def _negative_sigmoid(x):
    # Cache exp so you won't have to calculate it twice
    exp = np.exp(x)
    return exp / (exp + 1)


def sigmoid(x):
    positive = x >= 0
    # Boolean array inversion is faster than another comparison
    negative = ~positive

    # empty contains juke hence will be faster to allocate than zeros
    result = np.empty_like(x)
    result[positive] = _positive_sigmoid(x[positive])
    result[negative] = _negative_sigmoid(x[negative])

    return result


N = int(1e4)
x = np.random.uniform(size=(N, N))

start: float = time.time()
for _ in range(50):
    y1 = np.where(x > 0, 1 / (1 + np.exp(-x)), np.exp(x) / (1 + np.exp(x)))
    y1 += 1
end: float = time.time()
print(end - start)

start: float = time.time()
for _ in range(50):
    y2 = np.piecewise(
        x,
        [x > 0],
        [lambda i: 1 / (1 + np.exp(-i)), lambda i: np.exp(i) / (1 + np.exp(i))],
    )
    y2 += 1
end: float = time.time()
print(end - start)

start: float = time.time()
for _ in range(50):
    y2 = sigmoid(x)
    y2 += 1
end: float = time.time()
print(end - start)

【讨论】：

这是一个很棒的解释。感谢您抽出宝贵时间！

【参考方案4】：

The accepted answer 是正确的，但正如this comment 所指出的，它计算了两个分支，因此存在问题。

相反，您可能想要使用np.piecewise()。这更快、更有意义（np.where不是旨在定义分段函数）并且不会因进入两个分支而导致误导性警告。

基准测试

源代码

import numpy as np
import time

N: int = int(1e+4)

np.random.seed(0)

x: np.ndarray = np.random.random((N, N))
x *= 1e+3

start: float = time.time()
y1 = np.where(x > 0, 1 / (1 + np.exp(-x)), np.exp(x) / (1 + np.exp(x)))
end: float = time.time()
print()
print(end - start)

start: float = time.time()
y2 = np.piecewise(x, [x > 0], [lambda i: 1 / (1 + np.exp(-i)), lambda i: np.exp(i) / (1 + np.exp(i))])
end: float = time.time()
print(end - start)

assert (np.array_equal(y1, y2))

结果

np.piecewise() 静音，速度快两倍！

test.py:12: RuntimeWarning: overflow encountered in exp
  y1 = np.where(x > 0, 1 / (1 + np.exp(-x)), np.exp(x) / (1 + np.exp(x)))
test.py:12: RuntimeWarning: invalid value encountered in true_divide
  y1 = np.where(x > 0, 1 / (1 + np.exp(-x)), np.exp(x) / (1 + np.exp(x)))

6.32736349105835
3.138420343399048

【讨论】：

参见my answer 进行比较以及为什么您的解决方案不适用于边缘情况（np.array 是整数）。

【参考方案5】：

您的代码的另一种替代方法如下：

def sigmoid(z):
    return [(1. / (1. + np.exp(-x)) if x >= 0 else (np.exp(x) / (1. + np.exp(x))) for x in z]

【讨论】：

使用numpy函数而不使用numpy向量化真的没有意义。

这是一个非常好的和干净的选择。非常感谢！

@DYZ 虽然我同意您的解决方案更短且更重要，但最好使用列表推导而不是 for 循环。我将编辑我的答案，以便每个人都可以看到这是问题中发布的代码的替代方案，而不是您的解决方案。

【参考方案6】：

我写了一个技巧，我猜 np.where 或 torch.where 都是以同样的方式实现来处理二进制条件的：

def sigmoid(x, max_v=1.0):    
    sign = (torch.sign(x) + 3)//3
    x = torch.abs(x)
    res = max_v/(1 + torch.exp(-x))
    res = res * sign + (1 - sign) * (max_v - res)
    return res