自定义MLP模型中logit(expit(50))=inf引发NaN问题的解决方法

问题判断 #

你对问题的分析完全正确：这是Sigmoid激活函数饱和引发的数值不稳定问题，最终导致反向传播中出现Inf和NaN。当Sigmoid输入过大时输出趋近于1，此时用logit计算导数会得到np.inf，后续运算中正负Inf相加产生NaN，进而污染整个参数矩阵。

核心解决方法 #

1. 修正Sigmoid导数的计算逻辑 #

Sigmoid的导数不需要通过logit反推，其数学性质为：
σ'(x) = σ(x) * (1 - σ(x))
直接用激活后的输出值计算导数，完全避免logit(1)导致的Inf问题。

2. 优化权重初始化 #

当前的np.random.rand(...) *2 -1均匀初始化容易导致加权和数值过大，改用Xavier初始化（适配Sigmoid/Tanh）：

np.random.randn(number_of_nodes[i], number_of_nodes[i-1]) * np.sqrt(1 / number_of_nodes[i-1])

该初始化方式保证每层输入的方差稳定，减少激活函数饱和的概率。

3. 数值截断（可选） #

对进入Sigmoid的加权和进行截断，限制在[-10,10]范围内，因为expit(10)≈0.99995，已经足够接近1，不会触发logit的Inf：

weighted_sum = np.matmul(weights[j-1], neuron_values[j-1]) + biases[j-1]
weighted_sum = np.clip(weighted_sum, -10, 10)
neuron_values[j] = activate(weighted_sum, "Sigmoid")

4. 替换激活函数（可选） #

隐藏层改用ReLU/Leaky ReLU，这类激活函数不会出现饱和问题，从根源上避免数值不稳定。

代码中的关键问题修正 #

1. 修正deactivate函数 #

def deactivate(values, activation_type):
    if activation_type == "Sigmoid":
        # 正确的Sigmoid导数计算
        return values * (1 - values)
    if activation_type == "ReLU":
        return np.where(values > 0, 1, 0)
    if activation_type == "Softmax":
        # Softmax导数的正确形式（针对单个样本）
        diag = np.diagflat(values)
        outer = np.outer(values, values)
        return diag - outer
    if activation_type == "Tanh":
        return 1 - np.power(values, 2)

2. 权重初始化优化 #

for i in range(1, len(number_of_nodes)):
    # Xavier初始化
    weights.append(np.random.randn(number_of_nodes[i], number_of_nodes[i-1]) * np.sqrt(1 / number_of_nodes[i-1]))
    biases.append(np.zeros((number_of_nodes[i], 1)))  # 偏置初始化为0即可
    neuron_values.append(np.zeros((number_of_nodes[i], 1)))

3. 输入处理修正 #

在feed_forward中确保输入是列向量：

def feed_forward(inputs):
    # 将输入转为(784,1)的列向量
    neuron_values[0] = inputs.reshape(-1, 1)
    for j in range(1, len(number_of_nodes) - 1):
        weighted_sum = np.matmul(weights[j-1], neuron_values[j-1]) + biases[j-1]
        # 可选：数值截断
        weighted_sum = np.clip(weighted_sum, -10, 10)
        neuron_values[j] = activate(weighted_sum, "Sigmoid")
    # 输出层Softmax
    weighted_sum = np.matmul(weights[j], neuron_values[j]) + biases[j]
    neuron_values[j+1] = activate(weighted_sum, "Softmax")
    return neuron_values[-1]

4. 损失函数优化（可选） #

分类任务更适合交叉熵损失，替代当前的MSE损失，能提升收敛速度和数值稳定性：

# 新增交叉熵损失计算
def cross_entropy_loss(outputs, targets):
    # 防止log(0)，加极小值
    return -np.sum(targets * np.log(outputs + 1e-10))

# 反向传播中输出层误差修正（交叉熵+Softmax的简化导数）
errors[-1] = outputs - targets  # 交叉熵损失下的误差项，无需乘以2

完整修正后的代码 #

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import expit, softmax


LEARNING_RATE = 0.01
TEST_DATA_AMOUNT = 0.8
EPOCHS = 20
number_of_nodes = [28 * 28, 16, 16, 10]


df = pd.read_csv("5000_records_784I_1O_mnist.csv", index_col=0)
x_data = df.iloc[:, : 28 * 28]
y_data = df.iloc[:, 28 * 28 :]
normalized_x_data = x_data / 255
x_train = normalized_x_data[
    : round(TEST_DATA_AMOUNT * len(normalized_x_data))
].to_numpy()
y_train = y_data[: round(TEST_DATA_AMOUNT * len(y_data))].to_numpy()
x_test = normalized_x_data[
    round(TEST_DATA_AMOUNT * len(normalized_x_data)) :
].to_numpy()
y_test = y_data[round(TEST_DATA_AMOUNT * len(y_data)) :].to_numpy()

NUMBER_OF_RECORDS = x_train.shape[0]

neuron_values = [np.zeros((number_of_nodes[0], 1))]
weights = []
biases = []
for i in range(1, len(number_of_nodes)):
    # Xavier初始化
    weights.append(np.random.randn(number_of_nodes[i], number_of_nodes[i-1]) * np.sqrt(1 / number_of_nodes[i-1]))
    biases.append(np.zeros((number_of_nodes[i], 1)))
    neuron_values.append(np.zeros((number_of_nodes[i], 1)))

predicted_outputs = []
network_errors = []


def activate(values, activation_type):
    if activation_type == "Sigmoid":
        return expit(values)
    if activation_type == "ReLU":
        return np.maximum(values, 0)
    if activation_type == "Softmax":
        return softmax(values)
    if activation_type == "Tanh":
        return np.tanh(values)


def deactivate(values, activation_type):
    if activation_type == "Sigmoid":
        # 正确的Sigmoid导数
        return values * (1 - values)
    if activation_type == "ReLU":
        return np.where(values > 0, 1, 0)
    if activation_type == "Softmax":
        # Softmax导数矩阵（单个样本）
        diag = np.diagflat(values)
        outer = np.outer(values, values)
        return diag - outer
    if activation_type == "Tanh":
        return 1 - np.power(values, 2)


def one_hot(value):
    one_hot_value = np.zeros((10, 1))
    one_hot_value[value] = 1
    return one_hot_value


def feed_forward(inputs):
    # 确保输入是列向量
    neuron_values[0] = inputs.reshape(-1, 1)
    for j in range(1, len(number_of_nodes) - 1):
        weighted_sum = np.matmul(weights[j-1], neuron_values[j-1]) + biases[j-1]
        # 数值截断，防止Sigmoid饱和
        weighted_sum = np.clip(weighted_sum, -10, 10)
        neuron_values[j] = activate(weighted_sum, "Sigmoid")
    # 输出层
    weighted_sum = np.matmul(weights[j], neuron_values[j]) + biases[j]
    neuron_values[j+1] = activate(weighted_sum, "Softmax")
    return neuron_values[-1]


def back_propagate(outputs, targets, learning_rate):
    errors = [[] for _ in range(1, len(number_of_nodes))]
    gradients = [[] for _ in range(1, len(number_of_nodes))]
    
    # 输出层：交叉熵损失下的误差项（替代MSE的2*(outputs-targets)）
    errors[-1] = outputs - targets
    # Softmax导数与误差的乘积，这里简化为直接相乘（因为单个样本时矩阵乘法等价于点乘）
    gradients[-1] = learning_rate * np.matmul(errors[-1].T, deactivate(outputs, "Softmax")).T
    biases[-1] -= gradients[-1]
    weights[-1] -= np.matmul(gradients[-1], neuron_values[-2].T)
    
    # 隐藏层反向传播
    for j in range(len(number_of_nodes) - 3, -1, -1):
        errors[j] = np.matmul(weights[j+1].T, errors[j+1])
        gradients[j] = learning_rate * errors[j] * deactivate(neuron_values[j+1], "Sigmoid")
        biases[j] -= gradients[j]
        weights[j] -= np.matmul(gradients[j], neuron_values[j].T)
    
    return outputs - targets


def guess(inputs):
    return feed_forward(inputs)


def train(inputs, target, epochs, learning_rate):
    for epoch in range(epochs):
        p = np.random.permutation(len(inputs))
        shuffled_inputs = inputs[p]
        shuffled_targets = target[p]
        epoch_total_error = 0
        for j, set_of_inputs in enumerate(shuffled_inputs):
            epoch_input = set_of_inputs
            epoch_target = shuffled_targets[j]
            outputs = feed_forward(epoch_input)
            error = back_propagate(outputs, one_hot(epoch_target), learning_rate)
            sum_of_squared_errors = np.sum(np.power(error, 2))
            epoch_total_error += sum_of_squared_errors
            network_errors.append(sum_of_squared_errors)
        # 可选：学习率衰减
        # learning_rate *= 0.99
        print(
            f"Epoch #{epoch+1}: Average squared error: {epoch_total_error/len(shuffled_inputs):.4f}, Last Target: {epoch_target}, Last Predict: {outputs.argmax()}"
        )


train(x_train, y_train, EPOCHS, LEARNING_RATE)
plt.figure()
plt.plot(network_errors)
plt.grid()
plt.title("Iteration Number vs Error")
plt.xlabel("Iteration Number")
plt.ylabel("Squared Error")
plt.show()
print("Trained weights:")
for layer_weights in weights:
    print(layer_weights.shape)
    print(layer_weights)
print("Trained biases:")
for layer_biases in biases:
    print(layer_biases.shape)
    print(layer_biases)

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