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

import numpy as np
from concurrent.futures import ThreadPoolExecutor
from typing import Tuple, List
import time
 
def find_best_split_one_feature(
    X: np.ndarray,
    g: np.ndarray,
    h: np.ndarray,
    feature_idx: int,
    lambda_: float = 1.0,
    gamma: float = 0.0
) -> Tuple[int, float, float]:
    """Find best split for a single feature (runs in parallel)."""
    n = len(g)
    G_total = np.sum(g)
    H_total = np.sum(h)
    score_parent = (G_total ** 2) / (H_total + lambda_)
    
    sorted_idx = np.argsort(X[:, feature_idx])
    G_left = 0.0
    H_left = 0.0
    best_gain = -np.inf
    best_threshold = None
    
    for i in range(n - 1):
        idx = sorted_idx[i]
        G_left += g[idx]
        H_left += h[idx]
        
        if X[sorted_idx[i], feature_idx] == X[sorted_idx[i+1], feature_idx]:
            continue
            
        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
            best_threshold = (X[sorted_idx[i], feature_idx] + 
                            X[sorted_idx[i+1], feature_idx]) / 2
    
    return feature_idx, best_threshold, best_gain
 
 
def parallel_split_finding(
    X: np.ndarray,
    g: np.ndarray,
    h: np.ndarray,
    n_threads: int = 4
) -> Tuple[int, float, float]:
    """Find best split using parallel feature evaluation."""
    n_features = X.shape[1]
    
    with ThreadPoolExecutor(max_workers=n_threads) as executor:
        futures = [
            executor.submit(find_best_split_one_feature, X, g, h, j)
            for j in range(n_features)
        ]
        results = [f.result() for f in futures]
    
    # Find global best across all features
    best = max(results, key=lambda x: x[2])
    return best
 
 
# Demonstration
np.random.seed(42)
n_samples = 50000
n_features = 100
 
X = np.random.randn(n_samples, n_features)
y = X[:, 0] + X[:, 1] - X[:, 2] + np.random.randn(n_samples) * 0.5
g = -y  # gradient for MSE
h = np.ones(n_samples)
 
print("Parallel Split Finding Demonstration")
print("=" * 60)
print(f"Data: {n_samples:,} samples × {n_features} features")
print()
 
# Sequential
start = time.time()
sequential_result = None
for j in range(n_features):
    result = find_best_split_one_feature(X, g, h, j)
    if sequential_result is None or result[2] > sequential_result[2]:
        sequential_result = result
sequential_time = time.time() - start
 
# Parallel with different thread counts
print(f"Sequential time: {sequential_time:.3f}s")
print()
 
for n_threads in [2, 4, 8]:
    start = time.time()
    parallel_result = parallel_split_finding(X, g, h, n_threads)
    parallel_time = time.time() - start
    speedup = sequential_time / parallel_time
    print(f"{n_threads} threads: {parallel_time:.3f}s (speedup: {speedup:.1f}×)")
 
print()
print(f"Best split: feature {sequential_result[0]}, "
      f"threshold {sequential_result[1]:.4f}, gain {sequential_result[2]:.4f}")

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

import numpy as np
import time
 
class ColumnBlock:
    """
    XGBoost-style column block for efficient split finding.
    Pre-sorts data once, then enables fast linear scans.
    """
    
    def __init__(self, X: np.ndarray):
        """
        Initialize column block by pre-sorting all features.
        
        Parameters:
        -----------
        X : array of shape (n_samples, n_features)
        """
        self.n_samples, self.n_features = X.shape
        
        # Pre-sort each feature (one-time cost)
        self.sorted_indices = []
        self.sorted_values = []
        
        for j in range(self.n_features):
            sorted_idx = np.argsort(X[:, j])
            self.sorted_indices.append(sorted_idx)
            self.sorted_values.append(X[sorted_idx, j])
    
    def find_split_feature(
        self,
        feature_idx: int,
        g: np.ndarray,
        h: np.ndarray,
        lambda_: float = 1.0,
        gamma: float = 0.0
    ) -> tuple:
        """
        Find best split for one feature using pre-sorted indices.
        """
        G_total = np.sum(g)
        H_total = np.sum(h)
        score_parent = (G_total ** 2) / (H_total + lambda_)
        
        sorted_idx = self.sorted_indices[feature_idx]
        sorted_vals = self.sorted_values[feature_idx]
        
        G_left = 0.0
        H_left = 0.0
        best_gain = -np.inf
        best_threshold = None
        
        # Linear scan through pre-sorted samples
        for i in range(self.n_samples - 1):
            sample_idx = sorted_idx[i]
            G_left += g[sample_idx]
            H_left += h[sample_idx]
            
            # Skip if same value
            if sorted_vals[i] == sorted_vals[i + 1]:
                continue
            
            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
                best_threshold = (sorted_vals[i] + sorted_vals[i + 1]) / 2
        
        return best_threshold, best_gain
 
 
# Demonstration
np.random.seed(42)
n_samples = 100000
n_features = 50
 
X = np.random.randn(n_samples, n_features)
g = np.random.randn(n_samples)
h = np.ones(n_samples)
 
print("Column Block Pre-Sorting Efficiency")
print("=" * 60)
 
# Without pre-sorting (sort each time)
def naive_split_finding(X, g, h, feature_idx):
    sorted_idx = np.argsort(X[:, feature_idx])  # Sort every time!
    G_left = 0.0
    for i in range(len(g) - 1):
        G_left += g[sorted_idx[i]]
    return G_left
 
# Time naive approach (re-sort each call)
start = time.time()
for iteration in range(10):  # Simulate 10 boosting iterations
    for j in range(n_features):
        _ = naive_split_finding(X, g, h, j)
naive_time = time.time() - start
 
# Time column block approach (pre-sort once)
start = time.time()
block = ColumnBlock(X)  # One-time pre-sort
presort_time = time.time() - start
 
start = time.time()
for iteration in range(10):
    for j in range(n_features):
        _ = block.find_split_feature(j, g, h)
block_time = time.time() - start
 
print(f"Naive (re-sort each time): {naive_time:.2f}s")
print(f"Column block pre-sort:     {presort_time:.2f}s (one-time)")
print(f"Column block iterations:   {block_time:.2f}s")
print(f"Column block total:        {presort_time + block_time:.2f}s")
print()
print(f"Speedup: {naive_time / (presort_time + block_time):.1f}×")
print()
print("With more iterations, the column block advantage grows!")

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

import xgboost as xgb
import numpy as np
 
# Note: This demonstrates the API. Actual out-of-core training 
# requires data in external memory format.
 
print("XGBoost Out-of-Core Training")
print("=" * 60)
 
# Method 1: External Memory DMatrix
# Data is kept on disk and streamed during training
 
external_memory_example = '''
# Create data file in LibSVM format (can be very large)
# data.txt:
# 0 0:1.0 2:3.5 10:2.1
# 1 1:2.0 5:1.5
# ...
 
# Load with external memory caching
dtrain = xgb.DMatrix('data.txt#cache_prefix')
# The #cache_prefix creates cache files for efficient access
 
# Or use a path with ?format specification
# dtrain = xgb.DMatrix('data.csv?format=csv#cache')
'''
 
print("External Memory API:")
print(external_memory_example)
print()
 
# Method 2: Iterative Loading with DataIter
data_iter_example = '''
import xgboost as xgb
import numpy as np
 
class DataIterator(xgb.DataIter):
    """Custom iterator for external data."""
    
    def __init__(self, file_paths):
        self.file_paths = file_paths
        self.current_idx = 0
        super().__init__()
    
    def next(self, input_data):
        """Load next batch of data."""
        if self.current_idx >= len(self.file_paths):
            return 0  # No more data
        
        # Load from disk (your actual loading logic)
        X, y = load_data_from_file(self.file_paths[self.current_idx])
        input_data(data=X, label=y)
        self.current_idx += 1
        return 1  # More data available
    
    def reset(self):
        """Reset to beginning of data."""
        self.current_idx = 0
 
# Usage
files = ['data_part1.npz', 'data_part2.npz', 'data_part3.npz']
iterator = DataIterator(files)
dtrain = xgb.DMatrix(iterator)
'''
 
print("Data Iterator API (for custom data loading):")
print(data_iter_example)
print()
 
# Memory requirements analysis
print("Memory Requirements Analysis")
print("-" * 60)
 
def estimate_memory(n_samples, n_features, sparse=False, density=1.0):
    """Estimate memory requirements for XGBoost training."""
    bytes_per_value = 4  # float32
    
    if sparse:
        data_bytes = n_samples * n_features * density * (4 + 4)  # value + index
    else:
        data_bytes = n_samples * n_features * bytes_per_value
    
    # Gradients and Hessians
    gradient_bytes = n_samples * 4 * 2  # g and h
    
    # Sorted indices (for exact split finding)
    sorted_indices_bytes = n_samples * n_features * 4  # int32
    
    # Histogram (for histogram-based)
    histogram_bytes = 256 * n_features * 4 * 2  # bins × features × (G, H)
    
    return {
        'data': data_bytes,
        'gradients': gradient_bytes,
        'sorted_indices': sorted_indices_bytes,
        'histograms': histogram_bytes,
        'total_exact': data_bytes + gradient_bytes + sorted_indices_bytes,
        'total_hist': data_bytes + gradient_bytes + histogram_bytes
    }
 
# Example calculations
scenarios = [
    (1_000_000, 100, False, 1.0, "1M × 100 (dense)"),
    (10_000_000, 500, False, 1.0, "10M × 500 (dense)"),
    (10_000_000, 500, True, 0.01, "10M × 500 (1% sparse)"),
]
 
for n, d, sparse, density, desc in scenarios:
    mem = estimate_memory(n, d, sparse, density)
    print(f"
{desc}:")
    print(f"  Data:       {mem['data'] / 1e9:.2f} GB")
    print(f"  Exact mode: {mem['total_exact'] / 1e9:.2f} GB")
    print(f"  Hist mode:  {mem['total_hist'] / 1e9:.2f} GB")

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

import xgboost as xgb
import numpy as np
import time
 
# GPU Training Setup
print("XGBoost GPU Training")
print("=" * 60)
 
# Generate sample data
np.random.seed(42)
n_samples = 500000
n_features = 100
X = np.random.randn(n_samples, n_features).astype(np.float32)
y = (X[:, 0] + X[:, 1] > 0).astype(np.float32)
 
print(f"Dataset: {n_samples:,} samples × {n_features} features")
print()
 
# GPU parameters
gpu_params = {
    'tree_method': 'gpu_hist',  # GPU histogram method
    'device': 'cuda',            # Use CUDA device
    'max_depth': 8,
    'learning_rate': 0.1,
    'objective': 'binary:logistic',
    'eval_metric': 'logloss',
    'n_estimators': 100,
}
 
# CPU parameters for comparison
cpu_params = {
    'tree_method': 'hist',  # CPU histogram method
    'n_jobs': -1,            # Use all CPU cores
    'max_depth': 8,
    'learning_rate': 0.1,
    'objective': 'binary:logistic',
    'eval_metric': 'logloss',
    'n_estimators': 100,
}
 
# Note: Actual GPU execution requires CUDA-enabled XGBoost
print("GPU Configuration Example:")
print("-" * 60)
for k, v in gpu_params.items():
    print(f"  {k}: {v}")
print()
 
# Additional GPU parameters
additional_gpu_params = """
Additional GPU-specific parameters:
-----------------------------------
gpu_id: int (default 0)
    GPU device ordinal (for multi-GPU)
    
max_bin: int (default 256)
    Number of histogram bins
    Lower = faster but less precise
    
deterministic_histogram: bool (default True)
    Ensure reproducible results
    Set to False for slight speedup
    
predictor: str
    'gpu_predictor' for GPU-based prediction
    Note: With tree_method='gpu_hist', this is automatic
    
sampling_method: str
    'gradient_based' can use GPU for GOSS-like sampling
"""
print(additional_gpu_params)
 
# Speedup expectations
print("
Expected Speedups (GPU vs CPU):")
print("-" * 60)
speedup_data = [
    ("Small data (<100K samples)", "1-2×", "Kernel launch overhead"),
    ("Medium data (100K-1M)", "5-10×", "Good GPU utilization"),
    ("Large data (>1M samples)", "10-20×", "Full GPU parallelism"),
    ("Very high dimensional", "3-8×", "Memory bandwidth limited"),
]
print(f"{'Scenario':<30} {'Speedup':<10} {'Notes'}")
print("-" * 60)
for scenario, speedup, notes in speedup_data:
    print(f"{scenario:<30} {speedup:<10} {notes}")
 
# Memory considerations
print("
" + "=" * 60)
print("GPU Memory Considerations:")
print("-" * 60)
 
def estimate_gpu_memory(n_samples, n_features, max_bin=256):
    """Estimate GPU memory requirements."""
    # Data
    data_bytes = n_samples * n_features * 4  # float32
    
    # Histogram buffers
    # Double-buffered: current + sibling
    hist_bytes = 2 * max_bin * n_features * 2 * 4  # (G, H) float32
    
    # Row indices, positions
    indices_bytes = n_samples * 4 * 2
    
    total = data_bytes + hist_bytes + indices_bytes
    return total
 
for n in [100_000, 1_000_000, 10_000_000]:
    mem = estimate_gpu_memory(n, 100)
    print(f"{n/1e6:.1f}M samples × 100 features: ~{mem/1e9:.2f} GB GPU memory")

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

import xgboost as xgb
 
# Method 1: Dask-based Distribution
dask_example = '''
import dask.dataframe as dd
import dask.array as da
from xgboost import dask as dxgb
from dask.distributed import Client
 
# Start Dask cluster
client = Client(n_workers=4)  # Or connect to existing cluster
 
# Load data as Dask DataFrame
df = dd.read_parquet("s3://bucket/large_data/")
X = df[feature_columns]
y = df[target_column]
 
# Create DaskDMatrix
dtrain = dxgb.DaskDMatrix(client, X, y)
 
# Train
params = {
    'tree_method': 'hist',  # or 'gpu_hist' if workers have GPUs
    'max_depth': 6,
    'learning_rate': 0.1,
    'objective': 'binary:logistic',
}
 
output = dxgb.train(
    client,
    params,
    dtrain,
    num_boost_round=100,
    evals=[(dtrain, 'train')]
)
 
model = output['booster']
'''
 
print("XGBoost Distributed Training")
print("=" * 60)
print("
Method 1: Dask-based Distribution")
print("-" * 60)
print(dask_example)
 
# Method 2: Spark-based Distribution
spark_example = '''
from xgboost.spark import SparkXGBClassifier
 
# Data is a Spark DataFrame
df = spark.read.parquet("hdfs://path/to/data")
 
# Configure and train
classifier = SparkXGBClassifier(
    features_col="features",
    label_col="label",
    num_workers=10,
    max_depth=6,
    n_estimators=100,
    use_gpu=True,  # If workers have GPUs
)
 
model = classifier.fit(df)
predictions = model.transform(test_df)
'''
 
print("
Method 2: Spark-based Distribution")
print("-" * 60)
print(spark_example)
 
# Method 3: Ray-based Distribution
ray_example = '''
from xgboost_ray import RayXGBClassifier
import ray
 
ray.init(address="auto")  # Connect to Ray cluster
 
# RayDMatrix for distributed data loading
from xgboost_ray import RayDMatrix
dtrain = RayDMatrix(
    data="s3://bucket/train.csv",
    label="target",
    filetype=RayFileType.CSV
)
 
params = {
    "max_depth": 8,
    "n_estimators": 100,
    "objective": "binary:logistic",
}
 
classifier = RayXGBClassifier(
    **params,
    ray_params=RayParams(
        num_actors=4,
        gpus_per_actor=1  # If using GPUs
    )
)
 
classifier.fit(dtrain)
'''
 
print("
Method 3: Ray-based Distribution")
print("-" * 60)
print(ray_example)
 
# Scaling considerations
print("
" + "=" * 60)
print("Distributed Training Considerations:")
print("-" * 60)
print("""
1. Data Partitioning:
   - Data should be evenly distributed across workers
   - Avoid data skew which creates stragglers
   
2. Communication Overhead:
   - Histogram reduction: O(features × bins × workers)
   - Becomes significant with many workers (>100)
   
3. Fault Tolerance:
   - Use checkpointing for long-running jobs
   - Dask/Ray handle worker failures automatically
   
4. Resource Allocation:
   - Memory per worker: ~2-4× data partition size
   - Network: 1 Gbps minimum, 10 Gbps recommended
   
5. When to Use Distributed:
   - Data doesn't fit in single machine memory
   - Training time on single machine is prohibitive
   - Need to utilize existing cluster infrastructure
""")

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

 
Performance Optimization Decision Tree
=======================================
 
Start here:
│
├── Data size < 100K samples?
│   └── Use tree_method='hist', n_jobs=-1 (CPU is fine)
│
├── Data size 100K - 10M samples?
│   │
│   ├── Have GPU with 8GB+ VRAM?
│   │   └── Use tree_method='gpu_hist'
│   │
│   └── No GPU?
│       └── Use tree_method='hist', n_jobs=-1
│
├── Data size > 10M samples?
│   │
│   ├── Fits in RAM of single machine?
│   │   ├── Have multi-GPU? → Multi-GPU with Dask-CUDA
│   │   └── Single GPU → tree_method='gpu_hist'
│   │
│   └── Doesn't fit in RAM?
│       ├── Have cluster? → Distributed (Dask/Spark/Ray)
│       └── Single machine → Out-of-core training
│
└── Very sparse data (>90% zeros)?
    └── Always use sparse matrices + tree_method='hist'
 

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

import xgboost as xgb
import numpy as np
 
def get_optimized_params(
    n_samples: int,
    n_features: int,
    sparsity: float = 0.0,
    has_gpu: bool = False,
    memory_gb: float = 16.0
) -> dict:
    """
    Get optimized XGBoost parameters based on data characteristics.
    
    Parameters:
    -----------
    n_samples : int
        Number of training samples
    n_features : int
        Number of features
    sparsity : float
        Fraction of zeros (0-1)
    has_gpu : bool
        Whether GPU is available
    memory_gb : float
        Available RAM in GB
        
    Returns:
    --------
    dict : Optimized XGBoost parameters
    """
    params = {
        'objective': 'binary:logistic',
        'eval_metric': 'logloss',
        'verbosity': 1,
    }
    
    # Estimate memory requirements
    dense_memory_gb = n_samples * n_features * 8 / 1e9
    effective_memory_gb = dense_memory_gb * (1 - sparsity)
    
    # Tree method selection
    if n_samples < 100_000:
        # Small data: exact or hist both fine
        params['tree_method'] = 'hist'
    elif has_gpu and effective_memory_gb < 14:  # Fits in GPU memory
        params['tree_method'] = 'gpu_hist'
        params['device'] = 'cuda'
    else:
        params['tree_method'] = 'hist'
    
    # Thread configuration
    if params['tree_method'] != 'gpu_hist':
        params['n_jobs'] = -1
    
    # Histogram bins
    if sparsity > 0.9:
        params['max_bin'] = 128  # Fewer bins for very sparse data
    else:
        params['max_bin'] = 256  # Default
    
    # Subsample for large data
    if n_samples > 1_000_000:
        params['subsample'] = 0.8
        params['colsample_bytree'] = 0.8
    
    # Memory constraints
    if effective_memory_gb > memory_gb * 0.6:
        # Memory tight: reduce features per tree
        params['colsample_bytree'] = min(0.5, params.get('colsample_bytree', 1.0))
        params['max_bin'] = 64
    
    return params
 
 
# Example usage
print("Optimized XGBoost Configuration Generator")
print("=" * 60)
 
scenarios = [
    (50_000, 100, 0.0, False, "Small dense data"),
    (1_000_000, 200, 0.0, True, "Medium dense data with GPU"),
    (5_000_000, 500, 0.95, False, "Large sparse data"),
    (10_000_000, 100, 0.0, True, "Large dense data with GPU"),
]
 
for n, d, sparsity, gpu, desc in scenarios:
    print(f"
{desc}:")
    print(f"  {n:,} samples × {d} features, {sparsity:.0%} sparse, GPU={gpu}")
    params = get_optimized_params(n, d, sparsity, gpu)
    for k, v in params.items():
        if k not in ['objective', 'eval_metric', 'verbosity']:
            print(f"  {k}: {v}")