如何在逻辑回归的 numpy 实现中避免 NaN？

Posted 2023-03-11

【中文标题】如何在逻辑回归的 numpy 实现中避免 NaN？

【英文标题】：How to avoid NaN in numpy implementation of logistic regression?

【发布时间】：2019-01-24 16:25:10

【问题描述】：

编辑：我已经取得了重大进展。我当前的问题是在我最后一次编辑之后写的，可以在没有上下文的情况下回答。

我目前在 Coursera 上关注 Andrew Ng 的 Machine Learning Course，并尝试在今天实现 logistic regression。

符号：

X 是一个(m x n)-矩阵，输入变量的向量为行（m 训练样本为n-1 变量，第一列的条目处处等于 1 以表示一个常数）。 y 是预期输出样本的对应向量（具有m 条目的列向量等于0 或1） theta 是模型系数的向量（带有n 条目的行向量）

对于输入行向量x，该模型将预测sigmoid(x * theta.T) 的概率，以获得积极的结果。

这是我的 Python3/numpy 实现：

import numpy as np

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

vec_sigmoid = np.vectorize(sigmoid)

def logistic_cost(X, y, theta):
    summands = np.multiply(y, np.log(vec_sigmoid(X*theta.T))) + np.multiply(1 - y, np.log(1 - vec_sigmoid(X*theta.T)))
    return - np.sum(summands) / len(y)


def gradient_descent(X, y, learning_rate, num_iterations):
    num_parameters = X.shape[1]                                 # dim theta
    theta = np.matrix([0.0 for i in range(num_parameters)])     # init theta
    cost = [0.0 for i in range(num_iterations)]

    for it in range(num_iterations):
        error = np.repeat(vec_sigmoid(X * theta.T) - y, num_parameters, axis=1)
        error_derivative = np.sum(np.multiply(error, X), axis=0)
        theta = theta - (learning_rate / len(y)) * error_derivative
        cost[it] = logistic_cost(X, y, theta)

    return theta, cost

这个实现似乎工作正常，但我在计算物流成本时遇到了问题。在某些时候，梯度下降算法会收敛到一个非常合适的 theta，然后会发生以下情况：

对于某些输入行 X_i 与预期结果 1 X * theta.T 将变为正数并具有良好的边距（例如 23.207）。这将导致sigmoid(X_i * theta) 变为完全 1.0000（这是因为我认为失去了精度）。这是一个很好的预测（因为预期结果等于1），但这会破坏后勤成本的计算，因为np.log(1 - vec_sigmoid(X*theta.T)) 将评估为NaN。这应该不是问题，因为该术语与1 - y = 0 相乘，但是一旦出现NaN 的值，整个计算就会中断（0 * NaN = NaN）。

我应该如何在矢量化实现中处理这个问题，因为np.multiply(1 - y, np.log(1 - vec_sigmoid(X*theta.T))) 是在X 的每一行中计算的（不仅是y = 0 的位置）？

输入示例：

X = np.matrix([[1. , 0. , 0. ],
               [1. , 1. , 0. ],
               [1. , 0. , 1. ],
               [1. , 0.5, 0.3],
               [1. , 1. , 0.2]])

y = np.matrix([[0],
               [1],
               [1],
               [0],
               [1]])

那么theta, _ = gradient_descent(X, y, 10000, 10000)（是的，在这种情况下我们可以设置学习率this很大）将theta设置为：

theta = np.matrix([[-3000.04008972,  3499.97995514,  4099.98797308]])

这将导致vec_sigmoid(X * theta.T) 成为真正好的预测：

np.matrix([[0.00000000e+00],      # 0
           [1.00000000e+00],      # 1
           [1.00000000e+00],      # 1
           [1.95334953e-09],      # nearly zero
           [1.00000000e+00]])     # 1

但 logistic_cost(X, y, theta) 的计算结果为 NaN。

编辑：

我想出了以下解决方案。我刚刚将logistic_cost 函数替换为：

def new_logistic_cost(X, y, theta):
    term1 = vec_sigmoid(X*theta.T)
    term1[y == 0] = 1
    term2 = 1 - vec_sigmoid(X*theta.T)
    term2[y == 1] = 1
    summands = np.multiply(y, np.log(term1)) + np.multiply(1 - y, np.log(term2))
    return - np.sum(summands) / len(y)

通过使用掩码，我只是在结果将乘以零的地方计算log(1)。现在log(0) 只会在梯度下降的错误实现中发生。

开放式问题：我怎样才能使这个解决方案更干净？是否有可能以更清洁的方式实现类似的效果？

【问题讨论】：

仅供参考：如果您不介意对 scipy 的依赖，可以将 vec_sigmoid 替换为 scipy.special.expit。

@WarrenWeckesser 我一般不介意，只是想通过自己实施一次来更好地理解算法。但这并不能解决我的问题，不是吗？

不，但它可能会更快，并且在函数的尾部更准确。我认为您可以使用scipy.special.xlog1py 来解决您的问题。即，将np.multiply(1 - y, np.log(1 - vec_sigmoid(X*theta.T))) 替换为xlog1py(1 - y, -expit(X*theta.T)))。我将此建议作为答案。

【参考方案1】：

如果您不介意使用 SciPy，可以从 scipy.special 导入 expit 和 xlog1py：

from scipy.special import expit, xlog1py

并替换表达式

np.multiply(1 - y, np.log(1 - vec_sigmoid(X*theta.T)))

与

xlog1py(1 - y, -expit(X*theta.T))

【参考方案2】：

我知道这是一个老问题，但我遇到了同样的问题，也许它可以在将来帮助其他人，我实际上通过在附加 X0 之前对数据实施规范化来解决它。

def normalize_data(X):
    mean = np.mean(X, axis=0)
    std = np.std(X, axis=0)
    return (X-mean) / std

之后一切正常！