1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127

import numpy as np
from typing import Tuple, Optional
 
def exact_greedy_split(
    X: np.ndarray,
    g: np.ndarray,
    h: np.ndarray,
    lambda_: float = 1.0,
    gamma: float = 0.0,
    min_child_weight: float = 1.0
) -> Tuple[Optional[int], Optional[float], float]:
    """
    Exact greedy algorithm to find the best split.
    
    Parameters:
    -----------
    X : array of shape (n_samples, n_features)
        Feature matrix
    g : array of shape (n_samples,)
        First-order gradients
    h : array of shape (n_samples,)
        Second-order gradients (Hessians)
    lambda_ : float
        L2 regularization parameter
    gamma : float
        Minimum split gain threshold
    min_child_weight : float
        Minimum sum of Hessians in each child
        
    Returns:
    --------
    best_feature : int or None
        Index of best feature (None if no valid split)
    best_threshold : float or None
        Threshold value for split
    best_gain : float
        Gain of the best split
    """
    n_samples, n_features = X.shape
    G_total = np.sum(g)
    H_total = np.sum(h)
    
    best_gain = -np.inf
    best_feature = None
    best_threshold = None
    
    # Score of the current node (before split)
    score_before = (G_total ** 2) / (H_total + lambda_)
    
    # For each feature
    for feature_idx in range(n_features):
        # Sort samples by this feature
        sorted_indices = np.argsort(X[:, feature_idx])
        
        # Running sums for left child
        G_left = 0.0
        H_left = 0.0
        
        # Scan through sorted samples
        for i in range(n_samples - 1):
            idx = sorted_indices[i]
            G_left += g[idx]
            H_left += h[idx]
            G_right = G_total - G_left
            H_right = H_total - H_left
            
            # Check minimum child weight constraint
            if H_left < min_child_weight or H_right < min_child_weight:
                continue
            
            # Skip if this value equals the next (no split point between)
            if X[sorted_indices[i], feature_idx] == X[sorted_indices[i+1], feature_idx]:
                continue
            
            # Compute gain
            score_left = (G_left ** 2) / (H_left + lambda_)
            score_right = (G_right ** 2) / (H_right + lambda_)
            gain = 0.5 * (score_left + score_right - score_before) - gamma
            
            if gain > best_gain:
                best_gain = gain
                best_feature = feature_idx
                # Threshold is midpoint between current and next value
                best_threshold = (X[sorted_indices[i], feature_idx] + 
                                  X[sorted_indices[i+1], feature_idx]) / 2
    
    if best_gain <= 0:
        return None, None, best_gain
    
    return best_feature, best_threshold, best_gain
 
 
# Demonstration
np.random.seed(42)
 
# Generate sample data
n_samples = 1000
n_features = 5
X = np.random.randn(n_samples, n_features)
y = 2 * X[:, 0] - 3 * X[:, 1] + np.random.randn(n_samples) * 0.5  # Linear relationship
 
# Simulate being at some point in boosting
y_pred = np.zeros(n_samples)  # Initial prediction
g = y_pred - y  # Gradient for MSE
h = np.ones(n_samples)  # Hessian for MSE
 
print("Exact Greedy Split Finding Algorithm")
print("=" * 55)
print(f"Dataset: {n_samples} samples, {n_features} features")
print(f"Regularization: λ=1.0, γ=0.0")
print()
 
# Find best split
best_feat, best_thresh, best_gain = exact_greedy_split(X, g, h)
 
print(f"Best split found:")
print(f"  Feature index: {best_feat}")
print(f"  Threshold: {best_thresh:.4f}")
print(f"  Gain: {best_gain:.4f}")
print()
 
# Show complexity
print("Computational Complexity:")
print(f"  - Sorting: O({n_features} × {n_samples} × log({n_samples}))")
print(f"           = O({n_features * n_samples * int(np.log2(n_samples))})")
print(f"  - Scanning: O({n_features} × {n_samples})")
print(f"  Total operations: ~{n_features * n_samples * (int(np.log2(n_samples)) + 1)}")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128

import numpy as np
from typing import List, Tuple
 
def compute_quantile_candidates(X: np.ndarray, n_bins: int = 256) -> List[np.ndarray]:
    """
    Compute quantile-based candidate split points for each feature.
    
    Parameters:
    -----------
    X : array of shape (n_samples, n_features)
    n_bins : int
        Number of bins (candidates = n_bins - 1)
        
    Returns:
    --------
    candidates : list of arrays
        Candidate threshold values for each feature
    """
    n_samples, n_features = X.shape
    candidates = []
    
    for feature_idx in range(n_features):
        feature_values = X[:, feature_idx]
        
        # Compute percentiles
        percentiles = np.linspace(0, 100, n_bins + 1)[1:-1]  # Exclude 0 and 100
        thresholds = np.percentile(feature_values, percentiles)
        
        # Remove duplicates
        thresholds = np.unique(thresholds)
        candidates.append(thresholds)
    
    return candidates
 
 
def approximate_split(
    X: np.ndarray,
    g: np.ndarray,
    h: np.ndarray,
    candidates: List[np.ndarray],
    lambda_: float = 1.0,
    gamma: float = 0.0
) -> Tuple[int, float, float]:
    """
    Approximate split finding using pre-computed quantile candidates.
    """
    n_samples, n_features = X.shape
    G_total = np.sum(g)
    H_total = np.sum(h)
    
    best_gain = -np.inf
    best_feature = -1
    best_threshold = 0.0
    
    score_parent = (G_total ** 2) / (H_total + lambda_)
    
    for feature_idx in range(n_features):
        feature_candidates = candidates[feature_idx]
        feature_values = X[:, feature_idx]
        
        for threshold in feature_candidates:
            # Split samples
            left_mask = feature_values <= threshold
            
            G_left = np.sum(g[left_mask])
            H_left = np.sum(h[left_mask])
            G_right = G_total - G_left
            H_right = H_total - H_left
            
            if H_left < 1.0 or H_right < 1.0:
                continue
            
            # Compute gain
            score_left = (G_left ** 2) / (H_left + lambda_)
            score_right = (G_right ** 2) / (H_right + lambda_)
            gain = 0.5 * (score_left + score_right - score_parent) - gamma
            
            if gain > best_gain:
                best_gain = gain
                best_feature = feature_idx
                best_threshold = threshold
    
    return best_feature, best_threshold, best_gain
 
 
# Demonstration: Compare exact vs approximate
np.random.seed(42)
n_samples = 10000
n_features = 10
X = np.random.randn(n_samples, n_features)
y = np.sin(X[:, 0]) + 0.5 * X[:, 1] ** 2 + np.random.randn(n_samples) * 0.1
 
g = -y  # Simplified gradients
h = np.ones(n_samples)
 
print("Approximate vs Exact Split Finding")
print("=" * 60)
print(f"Dataset: {n_samples} samples, {n_features} features")
print()
 
# Exact (for comparison)
import time
start = time.time()
exact_feat, exact_thresh, exact_gain = exact_greedy_split(X, g, h)
exact_time = time.time() - start
print(f"Exact Greedy:")
print(f"  Best feature: {exact_feat}, threshold: {exact_thresh:.4f}")
print(f"  Gain: {exact_gain:.4f}")
print(f"  Time: {exact_time:.4f}s")
print()
 
# Approximate with different numbers of bins
for n_bins in [16, 64, 256]:
    candidates = compute_quantile_candidates(X, n_bins)
    start = time.time()
    approx_feat, approx_thresh, approx_gain = approximate_split(
        X, g, h, candidates
    )
    approx_time = time.time() - start
    
    total_candidates = sum(len(c) for c in candidates)
    gain_ratio = approx_gain / exact_gain if exact_gain != 0 else 1.0
    
    print(f"Approximate ({n_bins} bins, {total_candidates} total candidates):")
    print(f"  Best feature: {approx_feat}, threshold: {approx_thresh:.4f}")
    print(f"  Gain: {approx_gain:.4f} ({gain_ratio:.1%} of exact)")
    print(f"  Time: {approx_time:.4f}s ({exact_time/approx_time:.1f}x speedup)")
    print()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100

import numpy as np
 
def weighted_quantile(values: np.ndarray, weights: np.ndarray, 
                      percentiles: np.ndarray) -> np.ndarray:
    """
    Compute weighted quantiles.
    
    Parameters:
    -----------
    values : array of shape (n,)
        Feature values
    weights : array of shape (n,)
        Weights (Hessians in XGBoost context)
    percentiles : array
        Percentiles to compute (0-100 scale)
        
    Returns:
    --------
    quantiles : array
        Weighted quantile values
    """
    # Sort by values
    sorted_idx = np.argsort(values)
    sorted_values = values[sorted_idx]
    sorted_weights = weights[sorted_idx]
    
    # Cumulative weight sum
    cumsum = np.cumsum(sorted_weights)
    total = cumsum[-1]
    
    # Normalize to 0-100 scale
    cumsum_pct = (cumsum / total) * 100
    
    # Find quantiles
    quantiles = np.interp(percentiles, cumsum_pct, sorted_values)
    
    return quantiles
 
 
def compare_weighted_unweighted_quantiles():
    """
    Demonstrate the difference between weighted and unweighted quantiles.
    """
    np.random.seed(42)
    
    # Simulate binary classification scenario
    n = 1000
    
    # Feature values (predictive of class)
    class_0 = np.random.randn(n // 2) - 1  # Class 0 samples
    class_1 = np.random.randn(n // 2) + 1  # Class 1 samples
    X = np.concatenate([class_0, class_1])
    
    # Simulate predictions from previous round
    # Samples far from boundary have extreme predictions
    logits = X * 0.5  # Correlated with feature
    probs = 1 / (1 + np.exp(-logits))
    
    # Hessians for logistic loss: h = p(1-p)
    h = probs * (1 - probs)
    
    print("Weighted vs Unweighted Quantile Comparison")
    print("=" * 60)
    print(f"Scenario: Binary classification with {n} samples")
    print()
    
    # Compute quantile candidates
    percentiles = np.array([10, 25, 50, 75, 90])
    
    # Unweighted (standard) quantiles
    unweighted = np.percentile(X, percentiles)
    
    # Weighted quantiles (by Hessian)
    weighted = weighted_quantile(X, h, percentiles)
    
    print("Quantile candidates (split thresholds):")
    print("-" * 60)
    print(f"{'Percentile':>12} {'Unweighted':>15} {'Weighted (by h)':>18}")
    print("-" * 60)
    for p, uw, w in zip(percentiles, unweighted, weighted):
        print(f"{p:>12}% {uw:>15.4f} {w:>18.4f}")
    
    print()
    print("Analysis:")
    print("-" * 60)
    
    # Show Hessian distribution
    boundary_mask = np.abs(X) < 1  # Near decision boundary
    h_boundary = np.mean(h[boundary_mask])
    h_extreme = np.mean(h[~boundary_mask])
    
    print(f"Mean Hessian near boundary (|X| < 1): {h_boundary:.4f}")
    print(f"Mean Hessian away from boundary: {h_extreme:.4f}")
    print(f"Ratio: {h_boundary/h_extreme:.2f}x")
    print()
    print("Weighted quantiles shift toward the boundary region where")
    print("samples have higher Hessians (more uncertain predictions).")
    print("This places more candidate split points where they matter most!")
 
compare_weighted_unweighted_quantiles()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183

import numpy as np
from typing import Tuple, List
 
class HistogramSplitFinder:
    """
    Histogram-based split finding for gradient boosting.
    """
    
    def __init__(self, n_bins: int = 256):
        self.n_bins = n_bins
        self.bin_edges = None  # Per-feature bin edges
        
    def fit_bin_edges(self, X: np.ndarray):
        """
        Compute bin edges for each feature (done once before training).
        """
        n_features = X.shape[1]
        self.bin_edges = []
        
        for j in range(n_features):
            # Compute percentile-based bin edges
            edges = np.percentile(
                X[:, j], 
                np.linspace(0, 100, self.n_bins + 1)
            )
            # Ensure uniqueness
            edges = np.unique(edges)
            self.bin_edges.append(edges)
    
    def bin_data(self, X: np.ndarray) -> np.ndarray:
        """
        Convert continuous features to bin indices.
        """
        n_samples, n_features = X.shape
        binned = np.zeros((n_samples, n_features), dtype=np.uint8)
        
        for j in range(n_features):
            binned[:, j] = np.digitize(X[:, j], self.bin_edges[j][1:-1])
        
        return binned
    
    def build_histogram(self, binned_data: np.ndarray, g: np.ndarray, 
                        h: np.ndarray, feature_idx: int) -> Tuple[np.ndarray, np.ndarray]:
        """
        Build histogram for a single feature.
        Returns (G_per_bin, H_per_bin).
        """
        n_bins = len(self.bin_edges[feature_idx]) - 1
        G_hist = np.zeros(n_bins)
        H_hist = np.zeros(n_bins)
        
        # Aggregate gradients into bins
        bins = binned_data[:, feature_idx]
        np.add.at(G_hist, bins.astype(int), g)
        np.add.at(H_hist, bins.astype(int), h)
        
        return G_hist, H_hist
    
    def find_best_split_from_histogram(
        self, 
        G_hist: np.ndarray, 
        H_hist: np.ndarray,
        feature_idx: int,
        lambda_: float = 1.0,
        gamma: float = 0.0
    ) -> Tuple[float, float]:
        """
        Find best split by scanning histogram.
        """
        G_total = np.sum(G_hist)
        H_total = np.sum(H_hist)
        score_parent = (G_total ** 2) / (H_total + lambda_)
        
        best_gain = -np.inf
        best_threshold = None
        
        G_left = 0.0
        H_left = 0.0
        
        bin_edges = self.bin_edges[feature_idx]
        
        for bin_idx in range(len(G_hist) - 1):
            G_left += G_hist[bin_idx]
            H_left += H_hist[bin_idx]
            G_right = G_total - G_left
            H_right = H_total - H_left
            
            if H_left < 1.0 or H_right < 1.0:
                continue
            
            score_left = (G_left ** 2) / (H_left + lambda_)
            score_right = (G_right ** 2) / (H_right + lambda_)
            gain = 0.5 * (score_left + score_right - score_parent) - gamma
            
            if gain > best_gain:
                best_gain = gain
                # Threshold is the right edge of current bin
                best_threshold = bin_edges[bin_idx + 1]
        
        return best_threshold, best_gain
    
    def find_best_split(
        self,
        X_binned: np.ndarray,
        g: np.ndarray,
        h: np.ndarray,
        lambda_: float = 1.0,
        gamma: float = 0.0
    ) -> Tuple[int, float, float]:
        """
        Find the best split across all features using histograms.
        """
        n_features = X_binned.shape[1]
        
        best_gain = -np.inf
        best_feature = -1
        best_threshold = None
        
        for j in range(n_features):
            # Build histogram for this feature
            G_hist, H_hist = self.build_histogram(X_binned, g, h, j)
            
            # Find best split from histogram
            threshold, gain = self.find_best_split_from_histogram(
                G_hist, H_hist, j, lambda_, gamma
            )
            
            if gain > best_gain:
                best_gain = gain
                best_feature = j
                best_threshold = threshold
        
        return best_feature, best_threshold, best_gain
 
 
# Demonstration
np.random.seed(42)
n_samples = 100000  # Large dataset
n_features = 50
 
X = np.random.randn(n_samples, n_features)
y = np.sin(X[:, 0] * 2) + 0.5 * X[:, 1] + np.random.randn(n_samples) * 0.1
g = -y
h = np.ones(n_samples)
 
print("Histogram-Based Split Finding")
print("=" * 60)
print(f"Dataset: {n_samples:,} samples, {n_features} features")
print()
 
# Initialize and fit histogram finder
hist_finder = HistogramSplitFinder(n_bins=256)
hist_finder.fit_bin_edges(X)
 
import time
 
# Bin the data (done once)
start = time.time()
X_binned = hist_finder.bin_data(X)
bin_time = time.time() - start
print(f"Binning time: {bin_time:.4f}s (one-time preprocessing)")
 
# Find split using histogram
start = time.time()
feature, threshold, gain = hist_finder.find_best_split(X_binned, g, h)
split_time = time.time() - start
 
print(f"Split finding time: {split_time:.4f}s")
print(f"Best split: feature {feature}, threshold {threshold:.4f}, gain {gain:.4f}")
print()
 
# Compare with exact (on smaller subset for fairness)
print("Comparison with Exact Greedy (on 10K sample subset):")
X_small = X[:10000]
g_small = g[:10000]
h_small = h[:10000]
 
start = time.time()
exact_feat, exact_thresh, exact_gain = exact_greedy_split(X_small, g_small, h_small)
exact_time = time.time() - start
print(f"  Exact on 10K samples: {exact_time:.4f}s")
print(f"  Histogram on 100K samples: {split_time:.4f}s")
print(f"  Histogram processes 10x more data in similar time!")

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

import xgboost as xgb
import numpy as np
import time
 
# Generate sample data
np.random.seed(42)
n_samples = 50000
n_features = 100
 
X = np.random.randn(n_samples, n_features)
y = (X[:, 0] + X[:, 1] > 0).astype(int)  # Binary classification
 
dtrain = xgb.DMatrix(X, label=y)
 
# Common parameters
base_params = {
    'objective': 'binary:logistic',
    'max_depth': 6,
    'learning_rate': 0.1,
    'n_jobs': -1,
    'verbosity': 0
}
 
print("XGBoost tree_method Comparison")
print("=" * 60)
print(f"Dataset: {n_samples:,} samples, {n_features} features")
print()
 
# Compare different tree methods
for tree_method in ['exact', 'approx', 'hist']:
    params = {**base_params, 'tree_method': tree_method}
    
    start = time.time()
    model = xgb.train(params, dtrain, num_boost_round=100)
    elapsed = time.time() - start
    
    print(f"tree_method='{tree_method}':")
    print(f"  Training time: {elapsed:.2f}s")
    print()
 
# The 'hist' method is usually 5-20x faster than 'exact'
# with virtually identical prediction quality
 
print("Practical Recommendations:")
print("-" * 60)
print("1. Use 'hist' (default) for most cases")
print("2. Use 'gpu_hist' if you have a GPU and large data")
print("3. Use 'exact' only for small data or debugging")
print("4. Use 'approx' rarely - 'hist' is usually better")