Computational Considerations in Linear Regression

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

import numpy as np
 
def demonstrate_catastrophic_cancellation():
    """
    Demonstrates how subtracting nearly equal numbers destroys precision.
    """
    # Case 1: sqrt(x+1) - sqrt(x) for large x
    x = 1e20
    
    # Naive computation (WRONG for large x)
    naive_result = np.sqrt(x + 1) - np.sqrt(x)
    
    # Algebraically equivalent but numerically stable form:
    # sqrt(x+1) - sqrt(x) = 1 / (sqrt(x+1) + sqrt(x))
    stable_result = 1.0 / (np.sqrt(x + 1) + np.sqrt(x))
    
    # True value (computed with arbitrary precision)
    true_value = 5e-11  # approximately
    
    print(f"x = {x:.0e}")
    print(f"Naive result:  {naive_result:.10e}")  # Shows 0.0!
    print(f"Stable result: {stable_result:.10e}")  # Correct answer
    print(f"Relative error (naive): {abs(naive_result - stable_result) / stable_result:.2%}")
    
    # Case 2: The classic quadratic formula failure
    # Solving x² - 10⁸x + 1 = 0
    # Roots are approximately 10⁸ and 10⁻⁸
    a, b, c = 1.0, -1e8, 1.0
    discriminant = b**2 - 4*a*c
    
    # Standard formula (fails for the small root)
    x1_standard = (-b + np.sqrt(discriminant)) / (2*a)
    x2_standard = (-b - np.sqrt(discriminant)) / (2*a)
    
    # Stable computation using alternative form for small root
    x1_stable = (-b + np.sqrt(discriminant)) / (2*a)
    x2_stable = c / (a * x1_stable)  # Uses Vieta's formula: x1 * x2 = c/a
    
    print(f"
Quadratic formula comparison:")
    print(f"Large root - Standard: {x1_standard:.10e}, Stable: {x1_stable:.10e}")
    print(f"Small root - Standard: {x2_standard:.10e}, Stable: {x2_stable:.10e}")
    print(f"True small root ≈ 1e-8")
 
demonstrate_catastrophic_cancellation()

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

import numpy as np
from numpy.linalg import cond, svd, lstsq
 
def analyze_condition_numbers():
    """
    Demonstrates condition number computation and its implications.
    """
    np.random.seed(42)
    
    # Create a design matrix with controlled condition number
    n, p = 100, 5
    
    # Well-conditioned case: X with orthonormal-ish columns
    X_good = np.random.randn(n, p)
    X_good, _ = np.linalg.qr(X_good)  # Make columns orthonormal
    X_good *= np.sqrt(n)  # Scale back
    
    # Ill-conditioned case: nearly collinear features
    X_bad = np.random.randn(n, p)
    X_bad[:, -1] = X_bad[:, 0] + 1e-8 * np.random.randn(n)  # Near-duplicate column
    
    print("=== Condition Number Analysis ===
")
    
    for name, X in [("Well-conditioned", X_good), ("Ill-conditioned", X_bad)]:
        # Compute singular values
        U, singular_values, Vt = svd(X, full_matrices=False)
        
        # Condition numbers
        kappa_X = cond(X)
        kappa_XtX = cond(X.T @ X)
        
        print(f"{name} Matrix:")
        print(f"  Singular values: {singular_values}")
        print(f"  κ(X) = {kappa_X:.2e}")
        print(f"  κ(XᵀX) = {kappa_XtX:.2e}")
        print(f"  κ(XᵀX) / κ(X)² = {kappa_XtX / kappa_X**2:.2f} (should be ≈1)")
        print(f"  Expected digits lost solving normal equations: ~{np.log10(kappa_XtX):.1f}")
        print()
        
        # Demonstrate actual error amplification
        # Create true coefficients and generate data
        beta_true = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
        y = X @ beta_true
        
        # Solve via normal equations (potentially unstable)
        try:
            beta_normal = np.linalg.solve(X.T @ X, X.T @ y)
        except np.linalg.LinAlgError:
            beta_normal = np.full(p, np.nan)
        
        # Solve via SVD (stable)
        beta_svd = lstsq(X, y, rcond=None)[0]
        
        error_normal = np.linalg.norm(beta_normal - beta_true)
        error_svd = np.linalg.norm(beta_svd - beta_true)
        
        print(f"  Reconstruction error (Normal Equations): {error_normal:.2e}")
        print(f"  Reconstruction error (SVD):              {error_svd:.2e}")
        print()
 
analyze_condition_numbers()

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

import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import cond, solve, lstsq
 
def demonstrate_normal_equations_failure():
    """
    Shows how normal equations fail for polynomial regression
    while QR/SVD methods remain stable.
    """
    np.random.seed(42)
    
    # Generate simple data
    n = 50
    x = np.linspace(0, 1, n)
    y_true = np.sin(2 * np.pi * x)
    y = y_true + 0.1 * np.random.randn(n)
    
    degrees = [2, 5, 8, 10, 12, 14, 16]
    
    print("Polynomial Regression: Normal Equations vs SVD")
    print("=" * 65)
    print(f"{'Degree':<8} {'κ(X)':<12} {'κ(XᵀX)':<15} {'MSE Normal':<12} {'MSE SVD':<12}")
    print("-" * 65)
    
    for d in degrees:
        # Build Vandermonde matrix
        X = np.vander(x, d + 1, increasing=True)
        
        # Condition numbers
        kappa_X = cond(X)
        XtX = X.T @ X
        
        # Check if XtX is numerically singular
        try:
            kappa_XtX = cond(XtX)
        except:
            kappa_XtX = np.inf
        
        # Normal equations (potentially fails)
        try:
            Xty = X.T @ y
            beta_normal = solve(XtX, Xty)
            y_pred_normal = X @ beta_normal
            mse_normal = np.mean((y_pred_normal - y)**2)
            mse_normal_str = f"{mse_normal:.2e}"
        except np.linalg.LinAlgError:
            mse_normal_str = "FAILED"
        
        # SVD-based solution (robust)
        beta_svd = lstsq(X, y, rcond=None)[0]
        y_pred_svd = X @ beta_svd
        mse_svd = np.mean((y_pred_svd - y)**2)
        
        print(f"{d:<8} {kappa_X:<12.2e} {kappa_XtX:<15.2e} {mse_normal_str:<12} {mse_svd:<12.2e}")
    
    print()
    print("Observation: As degree increases, κ(XᵀX) explodes exponentially.")
    print("Normal equations become unreliable around degree 10-12.")
    print("SVD remains stable regardless of conditioning.")
 
demonstrate_normal_equations_failure()

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

import numpy as np
from numpy.linalg import matrix_rank, cond, svd
 
def diagnose_conditioning(X, feature_names=None):
                                        """
    Comprehensive diagnostic for design matrix conditioning.
    
    Parameters:
    -----------
    X : ndarray of shape (n_samples, n_features)
        Design matrix to analyze
    feature_names : list of str, optional
        Names of features for interpretable output
        
    Returns:
    --------
    dict : Diagnostic information
    """
    n, p = X.shape
    
    if feature_names is None:
        feature_names = [f"x{i}" for i in range(p)]
    
    # Basic properties
    U, s, Vt = svd(X, full_matrices=False)
    rank = matrix_rank(X)
    kappa = cond(X)
    
    # Identify problematic directions
    # Small singular values indicate near-linear-dependence
    sv_threshold = s[0] * 1e-10  # Relative threshold
    problematic_idx = np.where(s < sv_threshold)[0]
    
    # Analyze which features contribute to ill-conditioning
    # Look at right singular vectors corresponding to small singular values
    V = Vt.T
    
    diagnostics = {
        "shape": (n, p),
        "rank": rank,
        "condition_number": kappa,
        "singular_values": s,
        "log10_condition": np.log10(kappa) if kappa < np.inf else np.inf,
    }
    
    print("=== Design Matrix Conditioning Diagnostics ===")
    print(f"Shape: {n} samples × {p} features")
    print(f"Numerical rank: {rank} (expecting {min(n, p)})")
    print(f"Condition number: {kappa:.2e}")
    print(f"Log₁₀(κ): {diagnostics['log10_condition']:.1f}")
    print()
    
    # Singular value analysis
    print("Singular Value Spectrum:")
    print(f"  Largest:  {s[0]:.4e}")
    print(f"  Smallest: {s[-1]:.4e}")
    print(f"  Ratio:    {s[0]/s[-1]:.4e}")
    print()
    
    # Check for scale disparity
    col_norms = np.linalg.norm(X, axis=0)
    scale_ratio = col_norms.max() / col_norms.min()
    print(f"Column Norm Disparity: {scale_ratio:.2e}")
    if scale_ratio > 100:
        print("  ⚠️ Consider standardizing features!")
    print()
    
    # Identify potential multicollinearity
    if len(problematic_idx) > 0:
        print("⚠️ Near-linear dependencies detected!")
        for idx in problematic_idx:
            # Which features contribute most to this null-space direction?
            v = V[:, idx]
            contributors = np.argsort(np.abs(v))[-3:][::-1]
            print(f"  σ_{idx} = {s[idx]:.2e}: involves {[feature_names[i] for i in contributors]}")
    else:
        print("✓ No near-exact linear dependencies detected.")
    
    # Guidance
    print()
    print("Recommendations:")
    if kappa > 1e12:
        print("  🔴 CRITICAL: Condition number too high for reliable double-precision.")
        print("     → Use regularization OR reformulate the problem")
    elif kappa > 1e8:
        print("  🟠 WARNING: Use QR or SVD, avoid normal equations entirely.")
    elif kappa > 1e4:
        print("  🟡 CAUTION: Use stable algorithms (QR preferred over normal equations).")
    else:
        print("  🟢 GOOD: Standard algorithms should work reliably.")
    
    return diagnostics
 
# Example usage
np.random.seed(42)
n, p = 100, 5
X = np.random.randn(n, p)
X[:, 4] = X[:, 0] + X[:, 1] + 1e-8 * np.random.randn(n)  # Introduce collinearity
diagnose_conditioning(X, ["A", "B", "C", "D", "A+B (collinear)"])