CatBoost

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
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255

import numpy as np
from typing import List, Tuple
 
class SymmetricTree:
    """
    A symmetric (oblivious) decision tree where all nodes at the same
    depth use the identical split.
    
    For depth D, the tree is fully described by:
    - D feature indices: which feature to split at each level
    - D thresholds: the split threshold at each level
    - 2^D leaf values: the prediction at each leaf
    """
    
    def __init__(self, depth: int):
        self.depth = depth
        self.n_leaves = 2 ** depth
        
        # Tree structure: splits at each level
        self.feature_indices = np.zeros(depth, dtype=np.int32)
        self.thresholds = np.zeros(depth, dtype=np.float32)
        
        # Leaf values
        self.leaf_values = np.zeros(self.n_leaves, dtype=np.float32)
    
    def get_leaf_index(self, x: np.ndarray) -> int:
        """
        Compute which leaf a sample reaches using bit manipulation.
        
        The leaf index is a D-bit number where bit k is 1 if 
        x[feature_k] >= threshold_k.
        
        This is O(D) but can be vectorized across samples.
        """
        leaf_index = 0
        for level in range(self.depth):
            feature_idx = self.feature_indices[level]
            threshold = self.thresholds[level]
            
            # Set bit if going right (feature >= threshold)
            if x[feature_idx] >= threshold:
                leaf_index |= (1 << level)
        
        return leaf_index
    
    def get_leaf_indices_vectorized(self, X: np.ndarray) -> np.ndarray:
        """
        Compute leaf indices for all samples using vectorized operations.
        
        For N samples and depth D, this is O(N * D) but highly parallelizable.
        """
        n_samples = X.shape[0]
        leaf_indices = np.zeros(n_samples, dtype=np.int32)
        
        for level in range(self.depth):
            feature_idx = self.feature_indices[level]
            threshold = self.thresholds[level]
            
            # Vectorized comparison
            goes_right = X[:, feature_idx] >= threshold
            
            # Update leaf indices for samples going right
            leaf_indices += goes_right.astype(np.int32) << level
        
        return leaf_indices
    
    def predict(self, X: np.ndarray) -> np.ndarray:
        """Predict for all samples. O(N * D) but vectorized."""
        leaf_indices = self.get_leaf_indices_vectorized(X)
        return self.leaf_values[leaf_indices]
 
 
class OrderedSymmetricTreeBuilder:
    """
    Build a symmetric tree with ordered boosting.
    
    Key insight: tree STRUCTURE is computed once, but leaf VALUES
    are computed incrementally using ordered target statistics.
    """
    
    def __init__(self, depth: int, prior: float = 0.0, prior_weight: float = 1.0):
        self.depth = depth
        self.prior = prior
        self.prior_weight = prior_weight
    
    def build_tree(
        self, 
        X: np.ndarray, 
        y: np.ndarray, 
        permutation: np.ndarray
    ) -> SymmetricTree:
        """
        Build a symmetric tree with ordered residual computation.
        
        Parameters:
        -----------
        X : array of shape (n_samples, n_features)
        y : array of shape (n_samples,) - target values
        permutation : array of shape (n_samples,) - ordering for ordered boosting
        """
        tree = SymmetricTree(self.depth)
        
        # Step 1: Find optimal splits for tree structure
        # (This uses all data - structure doesn't cause leakage)
        tree.feature_indices, tree.thresholds = self._find_splits(X, y)
        
        # Step 2: Compute ordered residuals
        residuals = self._compute_ordered_residuals(X, y, tree, permutation)
        
        # Step 3: Fit leaf values to residuals
        tree.leaf_values = self._fit_leaf_values(X, residuals, tree)
        
        return tree
    
    def _find_splits(
        self, 
        X: np.ndarray, 
        y: np.ndarray
    ) -> Tuple[np.ndarray, np.ndarray]:
        """Find optimal splits at each level using greedy optimization."""
        n_samples, n_features = X.shape
        
        feature_indices = np.zeros(self.depth, dtype=np.int32)
        thresholds = np.zeros(self.depth, dtype=np.float32)
        
        # Current partition indicator for each sample
        partition = np.zeros(n_samples, dtype=np.int32)
        
        for level in range(self.depth):
            best_gain = -np.inf
            best_feature = 0
            best_threshold = 0.0
            
            # Try each feature and threshold
            for feature in range(n_features):
                values = X[:, feature]
                unique_vals = np.unique(values)
                
                for thresh in unique_vals[:-1]:  # Candidate thresholds
                    gain = self._compute_split_gain(
                        y, partition, values, thresh, level
                    )
                    
                    if gain > best_gain:
                        best_gain = gain
                        best_feature = feature
                        best_threshold = thresh
            
            feature_indices[level] = best_feature
            thresholds[level] = best_threshold
            
            # Update partition for next level
            goes_right = X[:, best_feature] >= best_threshold
            partition = partition * 2 + goes_right.astype(np.int32)
        
        return feature_indices, thresholds
    
    def _compute_split_gain(
        self,
        y: np.ndarray,
        partition: np.ndarray,
        feature_values: np.ndarray,
        threshold: float,
        level: int
    ) -> float:
        """
        Compute gain from a split, summed across all nodes at this level.
        
        In symmetric trees, the same split is applied to all 2^level nodes.
        """
        n_partitions = 2 ** level
        total_gain = 0.0
        
        for p in range(n_partitions):
            # Samples in this partition
            mask = partition == p
            if mask.sum() == 0:
                continue
            
            y_part = y[mask]
            vals_part = feature_values[mask]
            
            # Split this partition
            left_mask = vals_part < threshold
            right_mask = ~left_mask
            
            if left_mask.sum() == 0 or right_mask.sum() == 0:
                continue
            
            # Variance reduction gain
            var_before = np.var(y_part) * len(y_part)
            var_left = np.var(y_part[left_mask]) * left_mask.sum()
            var_right = np.var(y_part[right_mask]) * right_mask.sum()
            
            total_gain += var_before - var_left - var_right
        
        return total_gain
    
    def _compute_ordered_residuals(
        self,
        X: np.ndarray,
        y: np.ndarray,
        tree: SymmetricTree,
        permutation: np.ndarray
    ) -> np.ndarray:
        """
        Compute residuals using ordered leaf values.
        
        For each sample, prediction uses only preceding samples' statistics.
        """
        n_samples = len(y)
        residuals = np.zeros(n_samples)
        
        # Running statistics per leaf
        leaf_sums = np.zeros(tree.n_leaves)
        leaf_counts = np.zeros(tree.n_leaves)
        
        for pos in range(n_samples):
            sample_idx = permutation[pos]
            leaf_idx = tree.get_leaf_index(X[sample_idx])
            
            # Prediction using ordered statistics
            if leaf_counts[leaf_idx] > 0:
                pred = (
                    (leaf_sums[leaf_idx] + self.prior_weight * self.prior) /
                    (leaf_counts[leaf_idx] + self.prior_weight)
                )
            else:
                pred = self.prior
            
            # Ordered residual (no leakage!)
            residuals[sample_idx] = y[sample_idx] - pred
            
            # Update statistics for subsequent samples
            leaf_sums[leaf_idx] += y[sample_idx]
            leaf_counts[leaf_idx] += 1
        
        return residuals
    
    def _fit_leaf_values(
        self,
        X: np.ndarray,
        residuals: np.ndarray,
        tree: SymmetricTree
    ) -> np.ndarray:
        """Fit final leaf values as mean residual per leaf."""
        leaf_indices = tree.get_leaf_indices_vectorized(X)
        leaf_values = np.zeros(tree.n_leaves)
        
        for leaf in range(tree.n_leaves):
            mask = leaf_indices == leaf
            if mask.sum() > 0:
                leaf_values[leaf] = residuals[mask].mean()
        
        return leaf_values

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

import numpy as np
import time
from catboost import CatBoostClassifier, Pool
 
def benchmark_inference_speed():
    """
    Benchmark CatBoost inference speed across batch sizes.
    """
    # Train a model
    np.random.seed(42)
    n_train = 10000
    n_features = 50
    
    X_train = np.random.randn(n_train, n_features)
    y_train = (X_train[:, 0] + X_train[:, 1] > 0).astype(int)
    
    model = CatBoostClassifier(
        iterations=100,
        depth=6,
        learning_rate=0.1,
        random_seed=42,
        verbose=0,
        # Enable for fastest inference
        task_type='CPU',
    )
    model.fit(X_train, y_train)
    
    # Benchmark different batch sizes
    batch_sizes = [1, 10, 100, 1000, 10000]
    
    print("Inference Speed Benchmark")
    print("=" * 60)
    
    for batch_size in batch_sizes:
        # Generate test data
        X_test = np.random.randn(batch_size, n_features)
        
        # Warm up
        _ = model.predict(X_test)
        
        # Time multiple iterations
        n_iterations = max(1000 // batch_size, 10)
        
        start = time.perf_counter()
        for _ in range(n_iterations):
            _ = model.predict(X_test)
        elapsed = time.perf_counter() - start
        
        total_predictions = batch_size * n_iterations
        throughput = total_predictions / elapsed
        latency_us = 1e6 * elapsed / total_predictions
        
        print(f"Batch size {batch_size:5d}: "
              f"{throughput:>10,.0f} pred/sec, "
              f"{latency_us:>6.2f} μs/sample")
    
    print("\nNote: Actual speeds depend on CPU, tree depth, and tree count.")
 
 
def compare_prediction_modes():
    """
    Compare different CatBoost prediction modes for speed vs accuracy.
    """
    np.random.seed(42)
    X_train = np.random.randn(5000, 30)
    y_train = (X_train[:, 0] + X_train[:, 1] > 0).astype(int)
    X_test = np.random.randn(10000, 30)
    
    model = CatBoostClassifier(
        iterations=200,
        depth=6,
        random_seed=42,
        verbose=0,
    )
    model.fit(X_train, y_train)
    
    print("\nPrediction Mode Comparison")
    print("=" * 60)
    
    # Standard prediction
    start = time.perf_counter()
    preds_standard = model.predict_proba(X_test)
    time_standard = time.perf_counter() - start
    
    # Virtual ensembles (uncertainty quantification, slower)
    start = time.perf_counter()
    preds_ve = model.virtual_ensembles_predict(
        X_test,
        prediction_type='TotalUncertainty',
        virtual_ensembles_count=10
    )
    time_ve = time.perf_counter() - start
    
    print(f"Standard prediction:       {1000*time_standard:>6.2f} ms")
    print(f"Virtual ensembles (10):    {1000*time_ve:>6.2f} ms")
    print(f"Slowdown for uncertainty:  {time_ve/time_standard:.1f}x")
 
 
# Run benchmarks
if __name__ == "__main__":
    benchmark_inference_speed()
    compare_prediction_modes()

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

from catboost import CatBoostClassifier
import numpy as np
import pandas as pd
from sklearn.model_selection import cross_val_score
 
def tune_depth_vs_iterations(X, y, depth_range=(3, 10)):
    """
    Find optimal depth and iteration balance.
    
    Strategy: For each depth, find optimal iterations via early stopping,
    then compare final performance.
    """
    results = []
    
    for depth in range(depth_range[0], depth_range[1] + 1):
        # Use early stopping to find optimal iterations
        model = CatBoostClassifier(
            iterations=2000,  # Max iterations
            learning_rate=0.05,
            depth=depth,
            early_stopping_rounds=50,
            random_seed=42,
            verbose=0,
        )
        
        # Cross-validation with early stopping
        cv_scores = []
        for train_idx, val_idx in KFold(5, shuffle=True, random_state=42).split(X):
            X_train, X_val = X[train_idx], X[val_idx]
            y_train, y_val = y[train_idx], y[val_idx]
            
            model.fit(X_train, y_train, eval_set=(X_val, y_val))
            cv_scores.append(model.best_score_['validation']['Logloss'])
        
        avg_score = np.mean(cv_scores)
        avg_iters = model.best_iteration_
        
        results.append({
            'depth': depth,
            'avg_logloss': avg_score,
            'avg_iterations': avg_iters,
            'total_leaves': 2 ** depth * avg_iters,  # Model complexity proxy
        })
        
        print(f"Depth {depth}: Logloss={avg_score:.4f}, Iters={avg_iters}")
    
    df = pd.DataFrame(results)
    best = df.loc[df['avg_logloss'].idxmin()]
    print(f"\nBest: Depth={int(best['depth'])}, "
          f"Logloss={best['avg_logloss']:.4f}")
    
    return df
 
 
def demonstrate_depth_overfitting():
    """
    Show how excessive depth leads to overfitting.
    """
    np.random.seed(42)
    n_samples = 1000
    n_features = 20
    
    # Simple linear relationship with noise
    X = np.random.randn(n_samples, n_features)
    y = (X[:, 0] + 0.5 * X[:, 1] + np.random.randn(n_samples) * 0.5 > 0).astype(int)
    
    from sklearn.model_selection import train_test_split
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, test_size=0.3, random_state=42
    )
    
    print("Overfitting Demonstration: Simple Data with Various Depths")
    print("=" * 65)
    
    for depth in [3, 6, 10, 14]:
        model = CatBoostClassifier(
            iterations=200,
            learning_rate=0.1,
            depth=depth,
            random_seed=42,
            verbose=0,
        )
        model.fit(X_train, y_train)
        
        train_acc = (model.predict(X_train) == y_train).mean()
        test_acc = (model.predict(X_test) == y_test).mean()
        gap = train_acc - test_acc
        
        print(f"Depth {depth:2d}: Train={train_acc:.3f}, "
              f"Test={test_acc:.3f}, Gap={gap:.3f}")
    
    print("\nNote: Higher depth = larger gap = more overfitting")
 
from sklearn.model_selection import KFold
 
# Example usage
if __name__ == "__main__":
    demonstrate_depth_overfitting()