SVM Optimization

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

import numpy as np
import matplotlib.pyplot as plt
 
def analyze_convergence_theory():
    """
    Demonstrate key concepts in SMO convergence.
    """
    
    # Example: Dual objective function (2D case for visualization)
    # W(α) = α₁ + α₂ - 0.5 * (α₁²Q₁₁ + 2α₁α₂Q₁₂ + α₂²Q₂₂)
    
    # Kernel matrix (positive definite for this example)
    Q = np.array([[1.0, 0.3],
                  [0.3, 1.0]])
    
    def objective(alpha):
        """Dual objective W(α)."""
        return np.sum(alpha) - 0.5 * alpha @ Q @ alpha
    
    def gradient(alpha):
        """Gradient ∇W(α) = 1 - Qα."""
        return np.ones(2) - Q @ alpha
    
    # Simulate SMO optimization trajectory
    C = 2.0  # Upper bound
    alpha = np.array([0.0, 0.0])  # Start at origin
    
    trajectory = [alpha.copy()]
    objectives = [objective(alpha)]
    
    for iteration in range(50):
        # Find which variable to update (simplified)
        grad = gradient(alpha)
        
        # KKT check: is α at a bound where it shouldn't move?
        kkt_violated = False
        
        for i in range(2):
            if alpha[i] == 0 and grad[i] > 1e-6:
                kkt_violated = True
            elif alpha[i] == C and grad[i] < -1e-6:
                kkt_violated = True
            elif 0 < alpha[i] < C and abs(grad[i]) > 1e-6:
                kkt_violated = True
        
        if not kkt_violated:
            break
        
        # Simple coordinate ascent (one variable at a time)
        i = iteration % 2
        
        # Optimal step size along direction i
        step = grad[i] / Q[i, i]
        new_alpha_i = alpha[i] + step
        
        # Clip to [0, C]
        new_alpha_i = np.clip(new_alpha_i, 0, C)
        
        alpha[i] = new_alpha_i
        trajectory.append(alpha.copy())
        objectives.append(objective(alpha))
    
    trajectory = np.array(trajectory)
    objectives = np.array(objectives)
    
    # Compute optimal solution analytically
    # For unbounded case: α* = Q⁻¹ · 1
    alpha_opt_unbounded = np.linalg.solve(Q, np.ones(2))
    alpha_opt = np.clip(alpha_opt_unbounded, 0, C)
    W_opt = objective(alpha_opt)
    
    print("SMO Convergence Analysis")
    print("=" * 50)
    print(f"Kernel matrix Q:\n{Q}")
    print(f"\nOptimal solution: α* = {alpha_opt}")
    print(f"Optimal objective: W* = {W_opt:.4f}")
    print(f"\nTrajectory ({len(trajectory)} points):")
    for i, (a, w) in enumerate(zip(trajectory[:10], objectives[:10])):
        gap = W_opt - w
        print(f"  Iter {i:2d}: α = [{a[0]:.4f}, {a[1]:.4f}], W = {w:.4f}, gap = {gap:.4f}")
    
    return trajectory, objectives, W_opt
 
 
def convergence_rate_visualization():
    """
    Visualize convergence rate: linear vs sublinear.
    """
    iterations = np.arange(1, 101)
    
    # Linear convergence: error decays exponentially
    linear_rate = 0.9
    linear_error = (1 - linear_rate) ** iterations
    
    # Sublinear convergence: error decays as 1/k
    sublinear_error = 1 / iterations
    
    # Very slow convergence: error decays as 1/sqrt(k)
    slow_error = 1 / np.sqrt(iterations)
    
    plt.figure(figsize=(10, 6))
    plt.semilogy(iterations, linear_error, 'b-', linewidth=2, label='Linear: O(ρᵏ)')
    plt.semilogy(iterations, sublinear_error, 'g-', linewidth=2, label='Sublinear: O(1/k)')
    plt.semilogy(iterations, slow_error, 'r-', linewidth=2, label='Slow: O(1/√k)')
    
    plt.xlabel('Iteration k', fontsize=12)
    plt.ylabel('Optimality Gap (log scale)', fontsize=12)
    plt.title('Convergence Rate Comparison', fontsize=14)
    plt.legend(fontsize=11)
    plt.grid(True, alpha=0.3)
    plt.xlim([0, 100])
    plt.ylim([1e-6, 1])
    
    return plt.gcf()
 
 
# Run analysis
trajectory, objectives, W_opt = analyze_convergence_theory()

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

import numpy as np
 
class ConvergenceMonitor:
    """
    Monitor convergence during SMO training.
    
    Provides multiple stopping criteria and diagnostic information.
    """
    
    def __init__(self, tol=1e-3, max_iter=10000, 
                 relative_tol=1e-6, patience=100):
        """
        Parameters:
            tol: KKT tolerance (optimality gap threshold)
            max_iter: Maximum number of iterations
            relative_tol: Relative objective change threshold
            patience: Iterations of no improvement before early stopping
        """
        self.tol = tol
        self.max_iter = max_iter
        self.relative_tol = relative_tol
        self.patience = patience
        
        # History for diagnostics
        self.gaps = []
        self.objectives = []
        self.sv_counts = []
        
        self.iteration = 0
        self.no_improvement_count = 0
        self.best_objective = -np.inf
    
    def check_convergence(self, alpha, y, grad, C, objective=None):
        """
        Check if optimization should stop.
        
        Parameters:
            alpha: Current Lagrange multipliers
            y: Labels
            grad: Gradient G = -y * ∇W
            C: Regularization parameter
            objective: Current dual objective (optional)
        
        Returns:
            (converged: bool, reason: str)
        """
        self.iteration += 1
        
        # Compute optimality gap (KKT-based)
        gap = self._compute_optimality_gap(alpha, y, grad, C)
        self.gaps.append(gap)
        
        # Track objective if provided
        if objective is not None:
            self.objectives.append(objective)
            
            if objective > self.best_objective + self.relative_tol * abs(self.best_objective):
                self.best_objective = objective
                self.no_improvement_count = 0
            else:
                self.no_improvement_count += 1
        
        # Track support vector count
        sv_count = np.sum(alpha > 0)
        self.sv_counts.append(sv_count)
        
        # Check stopping criteria
        
        # 1. KKT tolerance
        if gap < self.tol:
            return True, f"Converged: optimality gap {gap:.2e} < tolerance {self.tol}"
        
        # 2. Maximum iterations
        if self.iteration >= self.max_iter:
            return True, f"Max iterations ({self.max_iter}) reached, gap = {gap:.2e}"
        
        # 3. No improvement for too long
        if self.no_improvement_count >= self.patience:
            return True, f"Early stopping: no improvement for {self.patience} iterations"
        
        # 4. Relative objective change (if enough history)
        if len(self.objectives) > 100:
            old_obj = self.objectives[-100]
            new_obj = self.objectives[-1]
            rel_change = abs(new_obj - old_obj) / (abs(new_obj) + 1e-10)
            
            if rel_change < self.relative_tol:
                return True, f"Converged: relative change {rel_change:.2e} < {self.relative_tol}"
        
        return False, ""
    
    def _compute_optimality_gap(self, alpha, y, grad, C):
        """Compute the optimality gap."""
        # I_up: indices where α can increase
        I_up = np.where(
            ((alpha < C) & (y == 1)) | 
            ((alpha > 0) & (y == -1))
        )[0]
        
        # I_low: indices where α can decrease
        I_low = np.where(
            ((alpha < C) & (y == -1)) | 
            ((alpha > 0) & (y == 1))
        )[0]
        
        if len(I_up) == 0 or len(I_low) == 0:
            return 0.0
        
        return np.max(grad[I_up]) - np.min(grad[I_low])
    
    def summary(self):
        """Return convergence summary statistics."""
        return {
            'iterations': self.iteration,
            'final_gap': self.gaps[-1] if self.gaps else None,
            'final_objective': self.objectives[-1] if self.objectives else None,
            'final_sv_count': self.sv_counts[-1] if self.sv_counts else None,
            'gap_reduction': self.gaps[0] / self.gaps[-1] if len(self.gaps) > 1 else 1,
        }
    
    def plot_convergence(self):
        """Generate convergence plots."""
        import matplotlib.pyplot as plt
        
        fig, axes = plt.subplots(2, 2, figsize=(12, 10))
        
        # Optimality gap (log scale)
        ax = axes[0, 0]
        ax.semilogy(self.gaps, 'b-', linewidth=1)
        ax.axhline(y=self.tol, color='r', linestyle='--', label=f'tolerance={self.tol}')
        ax.set_xlabel('Iteration')
        ax.set_ylabel('Optimality Gap')
        ax.set_title('Convergence of Optimality Gap')
        ax.legend()
        ax.grid(True, alpha=0.3)
        
        # Objective value
        if self.objectives:
            ax = axes[0, 1]
            ax.plot(self.objectives, 'g-', linewidth=1)
            ax.set_xlabel('Iteration')
            ax.set_ylabel('Dual Objective W(α)')
            ax.set_title('Objective Function Progress')
            ax.grid(True, alpha=0.3)
        
        # Support vector count
        ax = axes[1, 0]
        ax.plot(self.sv_counts, 'm-', linewidth=1)
        ax.set_xlabel('Iteration')
        ax.set_ylabel('Number of Support Vectors')
        ax.set_title('Support Vector Discovery')
        ax.grid(True, alpha=0.3)
        
        # Gap per iteration (rate estimate)
        if len(self.gaps) > 10:
            ax = axes[1, 1]
            gaps = np.array(self.gaps)
            gap_ratio = gaps[1:] / gaps[:-1]
            ax.plot(gap_ratio, 'c-', linewidth=0.5, alpha=0.7)
            ax.axhline(y=1.0, color='r', linestyle='--')
            ax.set_xlabel('Iteration')
            ax.set_ylabel('Gap Ratio (gap[k]/gap[k-1])')
            ax.set_title('Convergence Rate Estimate')
            ax.set_ylim([0, 1.5])
            ax.grid(True, alpha=0.3)
        
        plt.tight_layout()
        return fig
 
 
# Example usage
def demonstrate_convergence_monitoring():
    """
    Simulate SMO training with convergence monitoring.
    """
    np.random.seed(42)
    n = 100
    
    monitor = ConvergenceMonitor(tol=1e-3, max_iter=500)
    
    # Simulate convergence: gap decreases exponentially then slows
    for i in range(500):
        # Simulate decreasing gap
        if i < 200:
            gap = 10 * (0.97 ** i)  # Fast initial convergence
        else:
            gap = 0.01 + 0.005 * np.random.rand()  # Slow tail
        
        # Create fake state consistent with this gap
        alpha = np.random.rand(n) * 0.5
        y = np.random.choice([-1, 1], n)
        grad = np.random.randn(n) * gap
        C = 1.0
        objective = -np.random.rand() * gap  # Larger gap = more negative
        
        converged, reason = monitor.check_convergence(alpha, y, grad, C, objective)
        
        if converged:
            print(f"Stopped at iteration {i}: {reason}")
            break
    
    print(f"\nConvergence Summary:")
    for key, value in monitor.summary().items():
        print(f"  {key}: {value}")
 
demonstrate_convergence_monitoring()

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

import numpy as np
 
class ConvergenceDiagnostics:
    """
    Tools for diagnosing SMO convergence issues.
    """
    
    @staticmethod
    def diagnose_kernel_conditioning(K, threshold=1e6):
        """
        Check kernel matrix conditioning.
        
        Parameters:
            K: Kernel matrix (n×n)
            threshold: Condition number above which we flag issues
        
        Returns:
            Dictionary with diagnostic information
        """
        # Compute eigenvalues
        eigenvalues = np.linalg.eigvalsh(K)
        
        lambda_max = np.max(eigenvalues)
        lambda_min = np.max([np.min(eigenvalues), 1e-10])  # Avoid division by zero
        condition_number = lambda_max / lambda_min
        
        # Check for numerical issues
        n_negative = np.sum(eigenvalues < 0)
        n_near_zero = np.sum(np.abs(eigenvalues) < 1e-10)
        
        is_well_conditioned = condition_number < threshold
        is_positive_semidefinite = n_negative == 0
        
        diagnosis = {
            'condition_number': condition_number,
            'lambda_max': lambda_max,
            'lambda_min': lambda_min,
            'n_negative_eigenvalues': n_negative,
            'n_near_zero_eigenvalues': n_near_zero,
            'well_conditioned': is_well_conditioned,
            'positive_semidefinite': is_positive_semidefinite,
        }
        
        if not is_positive_semidefinite:
            diagnosis['recommendation'] = "Kernel has negative eigenvalues. Use a valid Mercer kernel."
        elif not is_well_conditioned:
            diagnosis['recommendation'] = f"Condition number {condition_number:.2e} is high. Try adjusting kernel parameters."
        else:
            diagnosis['recommendation'] = "Kernel matrix looks healthy."
        
        return diagnosis
    
    @staticmethod
    def diagnose_support_vectors(alpha, C, threshold_ratio=0.9):
        """
        Analyze support vector distribution.
        
        Parameters:
            alpha: Learned Lagrange multipliers
            C: Regularization parameter
            threshold_ratio: Fraction at which to warn about too many SVs
        
        Returns:
            Dictionary with diagnostic information
        """
        n = len(alpha)
        
        n_sv = np.sum(alpha > 0)
        n_bound = np.sum(np.isclose(alpha, C))
        n_free = np.sum((alpha > 0) & (alpha < C - 1e-8))
        n_zero = np.sum(alpha == 0)
        
        sv_ratio = n_sv / n
        bound_ratio = n_bound / n_sv if n_sv > 0 else 0
        
        diagnosis = {
            'n_samples': n,
            'n_support_vectors': n_sv,
            'n_bound_svs': n_bound,
            'n_free_svs': n_free,
            'n_non_svs': n_zero,
            'sv_ratio': sv_ratio,
            'bound_sv_ratio': bound_ratio,
        }
        
        if sv_ratio > threshold_ratio:
            diagnosis['issue'] = 'too_many_svs'
            diagnosis['recommendation'] = f"SV ratio {sv_ratio:.1%} is high. Consider reducing C or changing kernel."
        elif bound_ratio > 0.9 and n_free < 10:
            diagnosis['issue'] = 'almost_all_bound'
            diagnosis['recommendation'] = "Almost all SVs at bounds. C may be too large or problem is degenerate."
        elif n_sv == 0:
            diagnosis['issue'] = 'no_svs'
            diagnosis['recommendation'] = "No support vectors! Check data and parameters."
        else:
            diagnosis['issue'] = None
            diagnosis['recommendation'] = "Support vector distribution looks reasonable."
        
        return diagnosis
    
    @staticmethod
    def diagnose_convergence_trajectory(gaps, window=100):
        """
        Analyze convergence trajectory to detect issues.
        
        Parameters:
            gaps: List of optimality gaps over iterations
            window: Window for trend analysis
        
        Returns:
            Dictionary with diagnostic information
        """
        gaps = np.array(gaps)
        
        if len(gaps) < window:
            return {'issue': 'insufficient_data', 'recommendation': 'Need more iterations for diagnosis.'}
        
        # Compute rate of decrease
        early_gaps = gaps[:window]
        late_gaps = gaps[-window:]
        
        early_decrease = (early_gaps[0] - early_gaps[-1]) / early_gaps[0]
        late_decrease = (late_gaps[0] - late_gaps[-1]) / late_gaps[0]
        
        # Check for oscillation
        late_std = np.std(late_gaps)
        late_mean = np.mean(late_gaps)
        oscillation_ratio = late_std / late_mean if late_mean > 0 else 0
        
        # Check for plateau
        is_plateau = late_decrease < 0.01
        is_oscillating = oscillation_ratio > 0.1
        
        # Estimate convergence rate (geometric)
        if len(gaps) > 2:
            ratios = gaps[1:] / gaps[:-1]
            avg_ratio = np.mean(ratios[-window:])
        else:
            avg_ratio = 1.0
        
        diagnosis = {
            'early_decrease': early_decrease,
            'late_decrease': late_decrease,
            'oscillation_ratio': oscillation_ratio,
            'avg_convergence_ratio': avg_ratio,
            'is_plateau': is_plateau,
            'is_oscillating': is_oscillating,
        }
        
        if is_oscillating:
            diagnosis['issue'] = 'oscillation'
            diagnosis['recommendation'] = "Gap is oscillating. Possible numerical issues or tolerance too tight."
        elif is_plateau and late_gaps[-1] > 1e-2:
            diagnosis['issue'] = 'premature_plateau'
            diagnosis['recommendation'] = "Converged to suboptimal solution. Check kernel and parameters."
        elif avg_ratio > 0.999:
            diagnosis['issue'] = 'very_slow'
            diagnosis['recommendation'] = "Very slow convergence. Consider early stopping or parameter tuning."
        else:
            diagnosis['issue'] = None
            diagnosis['recommendation'] = "Convergence trajectory looks normal."
        
        return diagnosis
 
 
# Demonstration
def run_diagnostics_demo():
    """Demonstrate convergence diagnostics."""
    np.random.seed(42)
    
    print("=== Convergence Diagnostics Demo ===\n")
    
    # 1. Kernel conditioning
    n = 50
    X = np.random.randn(n, 10)
    
    # Well-conditioned kernel (moderate gamma)
    gamma = 0.1
    K_good = np.exp(-gamma * np.sum((X[:, None, :] - X[None, :, :]) ** 2, axis=2))
    
    print("1. Well-conditioned RBF kernel (γ=0.1):")
    diag = ConvergenceDiagnostics.diagnose_kernel_conditioning(K_good)
    for k, v in diag.items():
        print(f"   {k}: {v}")
    
    # Ill-conditioned kernel (extreme gamma)
    gamma = 100
    K_bad = np.exp(-gamma * np.sum((X[:, None, :] - X[None, :, :]) ** 2, axis=2))
    
    print("\n2. Ill-conditioned RBF kernel (γ=100):")
    diag = ConvergenceDiagnostics.diagnose_kernel_conditioning(K_bad)
    for k, v in diag.items():
        print(f"   {k}: {v}")
    
    # 2. Support vector analysis
    print("\n3. Support vector analysis:")
    
    # Case: too many SVs
    alpha_bad = np.random.rand(100) * 0.8
    diag = ConvergenceDiagnostics.diagnose_support_vectors(alpha_bad, C=1.0)
    print(f"   High SV case: {diag['recommendation']}")
    
    # Case: healthy
    alpha_good = np.zeros(100)
    alpha_good[10:25] = 0.5  # 15 free SVs
    alpha_good[25:30] = 1.0  # 5 bound SVs
    diag = ConvergenceDiagnostics.diagnose_support_vectors(alpha_good, C=1.0)
    print(f"   Healthy case: {diag['recommendation']}")
 
run_diagnostics_demo()