直方图算法 #

┌─────────────────────────────────────────────────────────────┐
│                    传统预排序方法                            │
├─────────────────────────────────────────────────────────────┤
│                                                             │
│  1. 对每个特征的所有值进行排序                               │
│  2. 遍历所有可能的分裂点                                     │
│  3. 时间复杂度: O(#data × #feature)                         │
│  4. 内存消耗: O(#data × #feature)                           │
│                                                             │
│  特征值: [1.2, 3.5, 2.1, 4.8, 0.5, 2.9]                     │
│  排序后: [0.5, 1.2, 2.1, 2.9, 3.5, 4.8]                     │
│  分裂点: 0.5, 1.2, 2.1, 2.9, 3.5 (5个候选)                  │
│                                                             │
└─────────────────────────────────────────────────────────────┘

┌─────────────────────────────────────────────────────────────┐
│                    直方图方法                                │
├─────────────────────────────────────────────────────────────┤
│                                                             │
│  1. 将连续值离散化到 k 个桶                                  │
│  2. 基于桶进行分裂点搜索                                     │
│  3. 时间复杂度: O(k × #feature)，k << #data                 │
│  4. 内存消耗: O(k × #feature)                               │
│                                                             │
│  特征值: [1.2, 3.5, 2.1, 4.8, 0.5, 2.9]                     │
│  桶边界: [0, 1, 2, 3, 4, 5]                                 │
│  桶索引: [ 1,  3,  2,  4,  0,  2]                           │
│  分裂点: 只需考虑 5 个桶边界                                 │
│                                                             │
└─────────────────────────────────────────────────────────────┘

import numpy as np

def build_histogram(feature_values, max_bin=255):
    """
    构建特征直方图
    
    Args:
        feature_values: 特征值数组
        max_bin: 最大桶数
    
    Returns:
        bin_upper_bounds: 桶的上边界
        bin_indices: 每个值所属的桶索引
    """
    feature_values = np.array(feature_values)
    
    min_val = feature_values.min()
    max_val = feature_values.max()
    
    bin_width = (max_val - min_val) / max_bin
    
    bin_upper_bounds = [min_val + (i + 1) * bin_width for i in range(max_bin)]
    bin_upper_bounds[-1] = max_val
    
    bin_indices = np.floor((feature_values - min_val) / bin_width).astype(int)
    bin_indices = np.clip(bin_indices, 0, max_bin - 1)
    
    return bin_upper_bounds, bin_indices

feature_values = np.array([1.2, 3.5, 2.1, 4.8, 0.5, 2.9, 3.1, 1.8])
bin_bounds, bin_indices = build_histogram(feature_values, max_bin=5)

print(f"特征值: {feature_values}")
print(f"桶边界: {bin_bounds}")
print(f"桶索引: {bin_indices}")

def compute_histogram_statistics(bin_indices, gradients, hessians, max_bin):
    """
    计算直方图统计量
    
    Args:
        bin_indices: 桶索引
        gradients: 梯度数组
        hessians: Hessian 数组
        max_bin: 最大桶数
    
    Returns:
        grad_hist: 梯度直方图
        hess_hist: Hessian 直方图
        count_hist: 计数直方图
    """
    grad_hist = np.zeros(max_bin)
    hess_hist = np.zeros(max_bin)
    count_hist = np.zeros(max_bin)
    
    for i, bin_idx in enumerate(bin_indices):
        grad_hist[bin_idx] += gradients[i]
        hess_hist[bin_idx] += hessians[i]
        count_hist[bin_idx] += 1
    
    return grad_hist, hess_hist, count_hist

gradients = np.array([0.1, -0.2, 0.3, -0.1, 0.2, -0.3, 0.1, 0.2])
hessians = np.array([0.5, 0.6, 0.4, 0.7, 0.5, 0.6, 0.4, 0.5])

grad_hist, hess_hist, count_hist = compute_histogram_statistics(
    bin_indices, gradients, hessians, max_bin=5
)

print(f"梯度直方图: {grad_hist}")
print(f"Hessian 直方图: {hess_hist}")
print(f"计数直方图: {count_hist}")

def find_best_split_histogram(grad_hist, hess_hist, count_hist, lambda_l2=0.1):
    """
    基于直方图找最佳分裂点
    
    Args:
        grad_hist: 梯度直方图
        hess_hist: Hessian 直方图
        count_hist: 计数直方图
        lambda_l2: L2 正则化参数
    
    Returns:
        best_split: 最佳分裂桶索引
        best_gain: 最佳增益
    """
    total_grad = np.sum(grad_hist)
    total_hess = np.sum(hess_hist)
    
    best_gain = -np.inf
    best_split = 0
    
    left_grad = 0
    left_hess = 0
    
    for i in range(len(grad_hist) - 1):
        left_grad += grad_hist[i]
        left_hess += hess_hist[i]
        
        right_grad = total_grad - left_grad
        right_hess = total_hess - left_hess
        
        if left_hess < 1e-10 or right_hess < 1e-10:
            continue
        
        left_score = (left_grad ** 2) / (left_hess + lambda_l2)
        right_score = (right_grad ** 2) / (right_hess + lambda_l2)
        total_score = (total_grad ** 2) / (total_hess + lambda_l2)
        
        gain = left_score + right_score - total_score
        
        if gain > best_gain:
            best_gain = gain
            best_split = i
    
    return best_split, best_gain

best_split, best_gain = find_best_split_histogram(grad_hist, hess_hist, count_hist)
print(f"最佳分裂桶: {best_split}")
print(f"最佳增益: {best_gain:.4f}")

import numpy as np

class HistogramDecisionTree:
    """基于直方图的决策树"""
    
    def __init__(self, max_bin=255, max_depth=3, min_samples=10, lambda_l2=0.1):
        self.max_bin = max_bin
        self.max_depth = max_depth
        self.min_samples = min_samples
        self.lambda_l2 = lambda_l2
        self.bin_bounds = {}
        self.tree = None
    
    def fit(self, X, y, gradients, hessians):
        """训练决策树"""
        n_features = X.shape[1]
        
        for feature in range(n_features):
            _, bin_bounds = build_histogram(X[:, feature], self.max_bin)
            self.bin_bounds[feature] = bin_bounds
        
        bin_matrix = np.zeros_like(X, dtype=int)
        for feature in range(n_features):
            bin_matrix[:, feature] = self._discretize(X[:, feature], feature)
        
        self.tree = self._build_tree(bin_matrix, gradients, hessians, depth=0)
    
    def _discretize(self, feature_values, feature):
        """离散化特征"""
        min_val = feature_values.min()
        max_val = feature_values.max()
        bin_width = (max_val - min_val) / self.max_bin
        
        bin_indices = np.floor((feature_values - min_val) / bin_width).astype(int)
        return np.clip(bin_indices, 0, self.max_bin - 1)
    
    def _build_tree(self, bin_matrix, gradients, hessians, depth):
        """递归构建树"""
        if depth >= self.max_depth or len(gradients) < self.min_samples:
            return {
                'leaf': True,
                'value': -np.sum(gradients) / (np.sum(hessians) + self.lambda_l2)
            }
        
        best_feature, best_split, best_gain = self._find_best_split(
            bin_matrix, gradients, hessians
        )
        
        if best_gain <= 0:
            return {
                'leaf': True,
                'value': -np.sum(gradients) / (np.sum(hessians) + self.lambda_l2)
            }
        
        left_mask = bin_matrix[:, best_feature] <= best_split
        right_mask = ~left_mask
        
        return {
            'leaf': False,
            'feature': best_feature,
            'split': best_split,
            'left': self._build_tree(
                bin_matrix[left_mask], 
                gradients[left_mask], 
                hessians[left_mask], 
                depth + 1
            ),
            'right': self._build_tree(
                bin_matrix[right_mask], 
                gradients[right_mask], 
                hessians[right_mask], 
                depth + 1
            )
        }
    
    def _find_best_split(self, bin_matrix, gradients, hessians):
        """找最佳分裂点"""
        best_gain = -np.inf
        best_feature = None
        best_split = None
        
        n_features = bin_matrix.shape[1]
        
        for feature in range(n_features):
            grad_hist, hess_hist, _ = compute_histogram_statistics(
                bin_matrix[:, feature], gradients, hessians, self.max_bin
            )
            
            split, gain = find_best_split_histogram(
                grad_hist, hess_hist, None, self.lambda_l2
            )
            
            if gain > best_gain:
                best_gain = gain
                best_feature = feature
                best_split = split
        
        return best_feature, best_split, best_gain
    
    def predict(self, X):
        """预测"""
        bin_matrix = np.zeros_like(X, dtype=int)
        for feature in range(X.shape[1]):
            bin_matrix[:, feature] = self._discretize(X[:, feature], feature)
        
        predictions = np.zeros(len(X))
        for i in range(len(X)):
            predictions[i] = self._predict_single(bin_matrix[i], self.tree)
        
        return predictions
    
    def _predict_single(self, bin_row, node):
        """单个样本预测"""
        if node['leaf']:
            return node['value']
        
        if bin_row[node['feature']] <= node['split']:
            return self._predict_single(bin_row, node['left'])
        else:
            return self._predict_single(bin_row, node['right'])

┌─────────────────────────────────────────────────────────────┐
│                    时间复杂度对比                            │
├─────────────────────────────────────────────────────────────┤
│                                                             │
│  传统预排序方法:                                             │
│  - 排序: O(#data × log(#data))                              │
│  - 分裂搜索: O(#data × #feature)                            │
│  - 总计: O(#data × #feature)                                │
│                                                             │
│  直方图方法:                                                 │
│  - 构建直方图: O(#data)                                      │
│  - 分裂搜索: O(k × #feature)，k << #data                    │
│  - 总计: O(k × #feature)                                    │
│                                                             │
│  加速比: #data / k (通常 100-1000 倍)                        │
│                                                             │
└─────────────────────────────────────────────────────────────┘

def compare_memory_usage(n_samples=100000, n_features=100, max_bin=255):
    """比较内存占用"""
    traditional_memory = n_samples * n_features * 4
    
    histogram_memory = max_bin * n_features * 4
    
    print(f"传统方法内存: {traditional_memory / 1024 / 1024:.2f} MB")
    print(f"直方图方法内存: {histogram_memory / 1024 / 1024:.2f} MB")
    print(f"内存减少: {(1 - histogram_memory / traditional_memory) * 100:.2f}%")

compare_memory_usage()

import time
from sklearn.tree import DecisionTreeRegressor

def benchmark_histogram_vs_traditional(n_samples=10000, n_features=50):
    """性能对比测试"""
    X = np.random.randn(n_samples, n_features)
    y = np.random.randn(n_samples)
    gradients = y - np.mean(y)
    hessians = np.ones(n_samples)
    
    start = time.time()
    tree = DecisionTreeRegressor(max_depth=5)
    tree.fit(X, gradients)
    traditional_time = time.time() - start
    
    start = time.time()
    hist_tree = HistogramDecisionTree(max_depth=5)
    hist_tree.fit(X, y, gradients, hessians)
    histogram_time = time.time() - start
    
    print(f"传统方法时间: {traditional_time:.4f}s")
    print(f"直方图方法时间: {histogram_time:.4f}s")
    print(f"加速比: {traditional_time / histogram_time:.2f}x")

benchmark_histogram_vs_traditional()

import matplotlib.pyplot as plt

def tune_max_bin(X, y):
    """调优 max_bin 参数"""
    max_bins = [15, 31, 63, 127, 255, 511]
    times = []
    mses = []
    
    gradients = y - np.mean(y)
    hessians = np.ones(len(y))
    
    for max_bin in max_bins:
        start = time.time()
        tree = HistogramDecisionTree(max_bin=max_bin, max_depth=5)
        tree.fit(X, y, gradients, hessians)
        elapsed = time.time() - start
        
        y_pred = tree.predict(X)
        mse = np.mean((y - y_pred) ** 2)
        
        times.append(elapsed)
        mses.append(mse)
    
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    
    axes[0].plot(max_bins, times, 'o-')
    axes[0].set_xlabel('max_bin')
    axes[0].set_ylabel('训练时间 (s)')
    axes[0].set_title('训练时间 vs max_bin')
    
    axes[1].plot(max_bins, mses, 'o-')
    axes[1].set_xlabel('max_bin')
    axes[1].set_ylabel('MSE')
    axes[1].set_title('MSE vs max_bin')
    
    plt.tight_layout()
    plt.show()

X = np.random.randn(5000, 20)
y = 2 * X[:, 0] + 3 * X[:, 1] + np.random.randn(5000) * 0.1
tune_max_bin(X, y)

✅ 训练速度快
   - 减少计算量
   - 适合大规模数据

✅ 内存占用低
   - 直方图存储
   - uint8 替代 float32

✅ 正则化效果
   - 离散化减少过拟合
   - 对噪声更鲁棒

✅ 支持类别特征
   - 自然处理类别特征
   - 无需独热编码

⚠️ 精度损失
   - 离散化丢失信息
   - 桶数太少影响精度

⚠️ 参数敏感
   - max_bin 需要调优
   - 不同数据最优值不同

⚠️ 不适合稀疏特征
   - 需要特殊处理
   - 可能增加计算量