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

optimal way of defining a numerically stable sigmoid function for a list in python

对于标量变量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 的方式。我应该如何改进 'cleanness' 的定义?说,有没有办法使用理解来缩短函数定义?

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

你说得对,使用 np.where 可以做得更好,if:

的 numpy 等价物
def sigmoid(x):
    return np.where(x >= 0, 
                    1 / (1 + np.exp(-x)), 
                    np.exp(x) / (1 + np.exp(x)))

这个函数接受一个 numpy 数组 x 和 returns 一个 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])

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

def sigmoid(z):
    return [(1. / (1. + np.exp(-x)) if x >= 0 else (np.exp(x) / (1. + np.exp(x))) for x in z]
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,我不太明白为什么要这样计算,但我知道这可能就是你要找的代码。希望对你有帮助。

is correct but, as pointed out by ,它计算了两个分支,因此有问题。

相反,您可能想要使用 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

我写了一个技巧,我猜 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

提供了完全正确的答案(没有警告),但没有清楚地解释解决方案。这对于评论来说太长了,所以我会去回答。

让我们从几个答案开始分析(纯numpy答案):

这个在数学上是正确的,但仍然给了我们一个警告。我们看代码:

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
    )

因为两个分支都被评估(它们是参数,它们必须是),第一个分支会给我们一个负值警告,第二个分支给正值警告。

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

缺点

  • 对两个分支进行不必要的评估(需要的操作次数增加一倍)
  • 抛出警告

这个几乎是正确的,但是只适用于浮点值,见下文:

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 在另一个线程中,但要点是:

It seems that piecewise() converts the return values to the same type as the input so, when an integer is input an integer conversion is performed on the result, which is then returned.

缺点

  • 由于分段函数的奇怪行为,没有自动转换

改进了 答案

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

如果编码正确,两个版本在操作方面都同样有效(一次 exp 评估就足够了)。现在:

  • e^xx为正时会溢出
  • e^-xx为负数时会溢出

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

请参阅代码注释以了解更多要点:

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)