python log - logistic 回归中遇到的除以零

python divide by zero encountered in log - logistic regression

我正在尝试实现区分 k 不同 类 的多类逻辑回归分类器。

这是我的代码。

import numpy as np
from scipy.special import expit


def cost(X,y,theta,regTerm):
    (m,n) = X.shape
    J = (np.dot(-(y.T),np.log(expit(np.dot(X,theta))))-np.dot((np.ones((m,1))-y).T,np.log(np.ones((m,1)) - (expit(np.dot(X,theta))).reshape((m,1))))) / m + (regTerm / (2 * m)) * np.linalg.norm(theta[1:])
    return J

def gradient(X,y,theta,regTerm):
    (m,n) = X.shape
    grad = np.dot(((expit(np.dot(X,theta))).reshape(m,1) - y).T,X)/m + (np.concatenate(([0],theta[1:].T),axis=0)).reshape(1,n)
    return np.asarray(grad)

def train(X,y,regTerm,learnRate,epsilon,k):
    (m,n) = X.shape
    theta = np.zeros((k,n))
    for i in range(0,k):
        previousCost = 0;
        currentCost = cost(X,y,theta[i,:],regTerm)
        while(np.abs(currentCost-previousCost) > epsilon):
            print(theta[i,:])
            theta[i,:] = theta[i,:] - learnRate*gradient(X,y,theta[i,:],regTerm)
            print(theta[i,:])
            previousCost = currentCost
            currentCost = cost(X,y,theta[i,:],regTerm)
    return theta

trX = np.load('trX.npy')
trY = np.load('trY.npy')
theta = train(trX,trY,2,0.1,0.1,4)

我可以验证成本和梯度正在返回正确维度的值(成本 returns 一个标量,梯度 returns 一个 1 × n 行向量),但我得到了错误

RuntimeWarning: divide by zero encountered in log
  J = (np.dot(-(y.T),np.log(expit(np.dot(X,theta))))-np.dot((np.ones((m,1))-y).T,np.log(np.ones((m,1)) - (expit(np.dot(X,theta))).reshape((m,1))))) / m + (regTerm / (2 * m)) * np.linalg.norm(theta[1:])

为什么会发生这种情况,我该如何避免这种情况?

我猜你的数据中有负值。您不能记录负数。

import numpy as np
np.log(2)
> 0.69314718055994529
np.log(-2)
> nan

如果是这种情况,有很多不同的方法可以帮助您转换数据。

您可以通过适当地使用广播、用于向量点积的运算符 * 和用于矩阵乘法的运算符 @ 来清理公式——并按照评论中的建议将其分解.

这是你的成本函数:

def cost(X, y, theta, regTerm):
    m = X.shape[0]  # or y.shape, or even p.shape after the next line, number of training set
    p = expit(X @ theta)
    log_loss = -np.average(y*np.log(p) + (1-y)*np.log(1-p))
    J = log_loss + regTerm * np.linalg.norm(theta[1:]) / (2*m)
    return J

您可以按照相同的方式清理渐变函数。

顺便问一下,你确定要np.linalg.norm(theta[1:])。如果您尝试进行 L2 正则化,则该术语应为 np.linalg.norm(theta[1:]) ** 2.

此处正确的解决方案是在 log 函数的参数中添加一些小的 epsilon。对我有用的是

epsilon = 1e-5    

def cost(X, y, theta):
    m = X.shape[0]
    yp = expit(X @ theta)
    cost = - np.average(y * np.log(yp + epsilon) + (1 - y) * np.log(1 - yp + epsilon))
    return cost
def cost(X, y, theta):
    yp = expit(X @ theta)
    cost = - np.average(y * np.log(yp) + (1 - y) * np.log(1 - yp))
    return cost

警告来自 np.log(yp)yp==0np.log(1 - yp)yp==1。一种选择是过滤掉这些值,而不是将它们传递给 np.log。另一种选择是添加一个小常量以防止值恰好为 0(如上面的评论之一所建议)

原因:

这是因为在某些情况下,每当 y[i] 等于 1 时,Sigmoid 函数的值 (theta) 也等于 1.

代价函数:

J = (np.dot(-(y.T),np.log(expit(np.dot(X,theta))))-np.dot((np.ones((m,1))-y).T,np.log(np.ones((m,1)) - (expit(np.dot(X,theta))).reshape((m,1))))) / m + (regTerm / (2 * m)) * np.linalg.norm(theta[1:])

现在,考虑上面代码片段中的以下部分:

np.log(np.ones((m,1)) - (expit(np.dot(X,theta))).reshape((m,1)))

在这里,当 theta 的值为 1 时,您正在执行 (1 - theta)。因此,这将有效地变成未定义的 log (1 - 1) = log (0)。

将epsilon值[这是一个微型值]添加到日志值中,这样它就不会成为问题了。 但是我不确定它是否会给出准确的结果。