🔢 Linear Algebra Basics: The Mathematics of Data Science

Imagine you're a movie recommendation system. You have millions of users and movies, and you need to find patterns in who likes what. That's linear algebra! 🎬 It's the language that lets us talk about high-dimensional spaces, transformations, and the hidden structures in data. From Google's PageRank to neural networks, linear algebra is the mathematical foundation of modern data science.

Vectors and Matrices: The Building Blocks 🏗️

If data science were LEGO, vectors would be the individual bricks and matrices would be the assembled structures. Vectors represent single observations (like one customer's purchases), while matrices represent entire datasets (all customers' purchases). Let's build with these mathematical LEGOs!

The GPS Analogy 🗺️

Think of vectors as GPS coordinates. A 2D vector [latitude, longitude] pinpoints a location. A 3D vector adds altitude. In data science, we work with hundreds of dimensions - each feature is another coordinate in our data space. Matrices? They're like having multiple GPS coordinates at once - tracking an entire fleet of vehicles!

import numpy as np

# Linear Algebra Fundamentals in NumPy 🔢

def vectors_and_basic_operations():
    """
    Vectors: The atoms of linear algebra.
    Everything starts here!
    """
    print("🎯 Vectors and Basic Operations")
    print("=" * 60)
    
    # Creating vectors
    v1 = np.array([3, 4])  # 2D vector
    v2 = np.array([1, 2])
    v3 = np.array([1, 2, 3, 4, 5])  # 5D vector
    
    print("Vector creation:")
    print(f"  v1 = {v1}")
    print(f"  v2 = {v2}")
    print(f"  v3 = {v3}")
    
    # Vector operations
    print("\n📐 Basic Vector Operations:")
    
    # Addition and subtraction
    print(f"  v1 + v2 = {v1 + v2}")
    print(f"  v1 - v2 = {v1 - v2}")
    
    # Scalar multiplication
    print(f"  2 * v1 = {2 * v1}")
    
    # Dot product (inner product)
    dot_product = np.dot(v1, v2)
    print(f"  v1 · v2 = {dot_product}")
    
    # Alternative syntax
    dot_product_alt = v1 @ v2  # Python 3.5+
    print(f"  v1 @ v2 = {dot_product_alt}")
    
    # Vector magnitude (norm)
    magnitude_v1 = np.linalg.norm(v1)
    print(f"  ||v1|| = {magnitude_v1:.3f}")
    
    # Unit vector (normalization)
    unit_v1 = v1 / magnitude_v1
    print(f"  v1_unit = {unit_v1}")
    print(f"  ||v1_unit|| = {np.linalg.norm(unit_v1):.3f}")
    
    # Angle between vectors
    cos_angle = dot_product / (np.linalg.norm(v1) * np.linalg.norm(v2))
    angle_rad = np.arccos(cos_angle)
    angle_deg = np.degrees(angle_rad)
    print(f"\n  Angle between v1 and v2: {angle_deg:.2f}°")
    
    # Cross product (only for 3D vectors)
    a = np.array([1, 0, 0])
    b = np.array([0, 1, 0])
    cross_product = np.cross(a, b)
    print(f"\n  [1,0,0] × [0,1,0] = {cross_product}")

def matrix_operations():
    """
    Matrices: Where linear algebra gets powerful!
    """
    print("\n🔲 Matrix Operations")
    print("=" * 60)
    
    # Creating matrices
    A = np.array([[1, 2], 
                  [3, 4]])
    
    B = np.array([[5, 6], 
                  [7, 8]])
    
    print("Matrix creation:")
    print(f"A =\n{A}")
    print(f"\nB =\n{B}")
    
    # Basic operations
    print("\n📊 Basic Matrix Operations:")
    
    # Addition and subtraction
    print(f"A + B =\n{A + B}")
    print(f"\nA - B =\n{A - B}")
    
    # Scalar multiplication
    print(f"\n2 * A =\n{2 * A}")
    
    # Element-wise multiplication (Hadamard product)
    print(f"\nA * B (element-wise) =\n{A * B}")
    
    # Matrix multiplication
    print(f"\nA @ B (matrix multiplication) =\n{A @ B}")
    
    # Alternative syntax
    print(f"\nnp.matmul(A, B) =\n{np.matmul(A, B)}")
    
    # Transpose
    print(f"\nA.T (transpose) =\n{A.T}")
    
    # Matrix properties
    print("\n📈 Matrix Properties:")
    
    # Determinant
    det_A = np.linalg.det(A)
    print(f"  det(A) = {det_A:.3f}")
    
    # Trace (sum of diagonal)
    trace_A = np.trace(A)
    print(f"  trace(A) = {trace_A}")
    
    # Rank
    rank_A = np.linalg.matrix_rank(A)
    print(f"  rank(A) = {rank_A}")
    
    # Inverse
    A_inv = np.linalg.inv(A)
    print(f"\nA^(-1) =\n{A_inv}")
    
    # Verify: A @ A^(-1) = I
    identity = A @ A_inv
    print(f"\nA @ A^(-1) =\n{identity.round(10)}")  # Round to avoid floating point errors

# Run basic examples
vectors_and_basic_operations()
matrix_operations()

Matrix Multiplication: The Heart of Transformations ❤️

Matrix multiplication is like a recipe transformer. If you have a recipe (matrix A) and ingredients (matrix B), multiplication gives you the final dish. But here's the magic: the order matters! A×B ≠ B×A, just like adding eggs before or after flour gives different results!

import numpy as np

class MatrixMultiplicationMaster:
    """
    Master the art of matrix multiplication.
    The most important operation in linear algebra!
    """
    
    def __init__(self):
        np.random.seed(42)
    
    def multiplication_basics(self):
        """Understanding matrix multiplication step by step"""
        print("🎯 Matrix Multiplication Deep Dive")
        print("=" * 60)
        
        # Simple example
        A = np.array([[1, 2, 3],
                      [4, 5, 6]])  # Shape: (2, 3)
        
        B = np.array([[7, 8],
                      [9, 10],
                      [11, 12]])  # Shape: (3, 2)
        
        print("Matrix A (2×3):")
        print(A)
        print("\nMatrix B (3×2):")
        print(B)
        
        # Manual calculation for understanding
        print("\n📐 Manual Calculation:")
        C = np.zeros((2, 2))
        
        for i in range(2):
            for j in range(2):
                for k in range(3):
                    C[i, j] += A[i, k] * B[k, j]
                print(f"  C[{i},{j}] = {' + '.join([f'{A[i,k]}×{B[k,j]}' for k in range(3)])} = {C[i,j]}")
        
        print(f"\nResult C (2×2):\n{C}")
        
        # NumPy's efficient way
        C_numpy = A @ B
        print(f"\nNumPy result (A @ B):\n{C_numpy}")
        print(f"Results match: {np.allclose(C, C_numpy)}")
        
        # Dimension rules
        print("\n📏 Dimension Rules:")
        print(f"  A shape: {A.shape}")
        print(f"  B shape: {B.shape}")
        print(f"  C shape: {C.shape}")
        print(f"  Rule: (m×n) @ (n×p) = (m×p)")
        print(f"  Inner dimensions must match: 3 == 3 ✓")
    
    def real_world_applications(self):
        """Real-world applications of matrix multiplication"""
        print("\n🌍 Real-World Applications")
        print("=" * 60)
        
        # 1. Linear Transformations
        print("1️⃣ Image Rotation:")
        
        # Rotation matrix for 45 degrees
        theta = np.pi / 4  # 45 degrees
        rotation_matrix = np.array([[np.cos(theta), -np.sin(theta)],
                                   [np.sin(theta), np.cos(theta)]])
        
        # Points of a square
        square = np.array([[0, 1, 1, 0],   # x-coordinates
                          [0, 0, 1, 1]])  # y-coordinates
        
        # Apply rotation
        rotated_square = rotation_matrix @ square
        
        print(f"  Original square corners:")
        for i in range(4):
            print(f"    Point {i+1}: ({square[0,i]:.1f}, {square[1,i]:.1f})")
        
        print(f"  Rotated square corners (45°):")
        for i in range(4):
            print(f"    Point {i+1}: ({rotated_square[0,i]:.2f}, {rotated_square[1,i]:.2f})")
        
        # 2. Neural Network Forward Pass
        print("\n2️⃣ Neural Network Layer:")
        
        # Input batch (4 samples, 3 features)
        X = np.random.randn(4, 3)
        
        # Weights (3 inputs, 2 neurons)
        W = np.random.randn(3, 2)
        
        # Bias
        b = np.array([0.5, -0.5])
        
        # Forward pass: Y = XW + b
        Y = X @ W + b
        
        # Apply activation (ReLU)
        Y_activated = np.maximum(0, Y)
        
        print(f"  Input shape: {X.shape} (batch_size, features)")
        print(f"  Weight shape: {W.shape} (features, neurons)")
        print(f"  Output shape: {Y.shape} (batch_size, neurons)")
        print(f"  Sample output:\n{Y_activated[:2].round(3)}")
        
        # 3. PageRank (Simplified)
        print("\n3️⃣ PageRank Algorithm (Simplified):")
        
        # Adjacency matrix (who links to whom)
        # Pages: A, B, C, D
        links = np.array([[0, 1, 1, 0],  # A links to B, C
                         [1, 0, 0, 1],  # B links to A, D
                         [0, 0, 0, 1],  # C links to D
                         [1, 1, 0, 0]]) # D links to A, B
        
        # Normalize columns (transition matrix)
        column_sums = links.sum(axis=0)
        column_sums[column_sums == 0] = 1  # Avoid division by zero
        transition_matrix = links / column_sums
        
        # Initial page ranks (equal)
        page_rank = np.ones(4) / 4
        
        # Power iteration
        damping = 0.85
        for iteration in range(5):
            page_rank = damping * (transition_matrix @ page_rank) + (1 - damping) / 4
        
        print(f"  Page ranks after 5 iterations:")
        pages = ['A', 'B', 'C', 'D']
        for i, page in enumerate(pages):
            print(f"    Page {page}: {page_rank[i]:.3f}")
        
        print(f"  Most important page: {pages[np.argmax(page_rank)]}")
    
    def matrix_chain_multiplication(self):
        """Optimizing matrix chain multiplication"""
        print("\n⚡ Matrix Chain Multiplication")
        print("=" * 60)
        
        # Three matrices to multiply
        A = np.random.randn(10, 100)
        B = np.random.randn(100, 5)
        C = np.random.randn(5, 50)
        
        print("Matrix dimensions:")
        print(f"  A: {A.shape}")
        print(f"  B: {B.shape}")
        print(f"  C: {C.shape}")
        
        # Two ways to compute A @ B @ C
        import time
        
        # Method 1: (A @ B) @ C
        start = time.perf_counter()
        result1 = (A @ B) @ C
        time1 = time.perf_counter() - start
        
        # Method 2: A @ (B @ C)
        start = time.perf_counter()
        result2 = A @ (B @ C)
        time2 = time.perf_counter() - start
        
        print(f"\n⏱️ Performance Comparison:")
        print(f"  (A @ B) @ C: {time1*1000:.3f} ms")
        print(f"  A @ (B @ C): {time2*1000:.3f} ms")
        print(f"  Speedup: {time1/time2:.2f}x")
        
        # Operations count
        ops1 = 10*100*5 + 10*5*50  # (A@B) then (@C)
        ops2 = 100*5*50 + 10*100*50  # (B@C) then (A@)
        
        print(f"\n🔢 Operations Count:")
        print(f"  (A @ B) @ C: {ops1:,} multiplications")
        print(f"  A @ (B @ C): {ops2:,} multiplications")
        print(f"  Ratio: {ops2/ops1:.2f}x more operations")
        
        print(f"\n💡 Lesson: Parentheses matter in matrix multiplication!")

# Explore matrix multiplication
mm = MatrixMultiplicationMaster()
mm.multiplication_basics()
mm.real_world_applications()
mm.matrix_chain_multiplication()

Eigenvalues and Eigenvectors: Finding Hidden Patterns 🔍

Eigenvectors are like the "natural directions" of a transformation - directions that don't change, only stretch or shrink. Eigenvalues tell us how much stretching happens. This concept powers Principal Component Analysis (PCA), Google's PageRank, and facial recognition systems!

import numpy as np

class EigenExplorer:
    """
    Exploring eigenvalues and eigenvectors.
    The key to understanding data's hidden structure!
    """
    
    def __init__(self):
        np.random.seed(42)
    
    def eigen_basics(self):
        """Understanding eigenvalues and eigenvectors"""
        print("🔍 Eigenvalues and Eigenvectors")
        print("=" * 60)
        
        # Simple 2x2 matrix
        A = np.array([[3, 1],
                      [1, 3]])
        
        print(f"Matrix A:\n{A}")
        
        # Compute eigenvalues and eigenvectors
        eigenvalues, eigenvectors = np.linalg.eig(A)
        
        print(f"\n📊 Eigenvalues: {eigenvalues}")
        print(f"\n📐 Eigenvectors:")
        for i in range(len(eigenvalues)):
            print(f"  λ{i+1} = {eigenvalues[i]:.3f}")
            print(f"  v{i+1} = {eigenvectors[:, i]}")
        
        # Verification: A @ v = λ * v
        print("\n✅ Verification (A @ v = λ × v):")
        for i in range(len(eigenvalues)):
            v = eigenvectors[:, i]
            λ = eigenvalues[i]
            
            Av = A @ v
            λv = λ * v
            
            print(f"\n  Eigenvector {i+1}:")
            print(f"    A @ v = {Av}")
            print(f"    λ × v = {λv}")
            print(f"    Match? {np.allclose(Av, λv)}")
        
        # Geometric interpretation
        print("\n🎨 Geometric Interpretation:")
        print("  • Eigenvectors = directions that don't change")
        print("  • Eigenvalues = scaling factors in those directions")
        print(f"  • In this case: stretch by {eigenvalues[0]:.1f} in one direction,")
        print(f"    and by {eigenvalues[1]:.1f} in the perpendicular direction")
    
    def pca_with_eigen(self):
        """PCA using eigendecomposition"""
        print("\n📈 Principal Component Analysis (PCA)")
        print("=" * 60)
        
        # Generate correlated 2D data
        mean = [0, 0]
        cov = [[1, 0.8], 
               [0.8, 1]]
        
        data = np.random.multivariate_normal(mean, cov, 1000)
        
        print(f"Dataset shape: {data.shape}")
        print(f"Features: 2, Samples: 1000")
        
        # Step 1: Center the data
        data_centered = data - data.mean(axis=0)
        
        # Step 2: Compute covariance matrix
        cov_matrix = np.cov(data_centered.T)
        print(f"\nCovariance matrix:\n{cov_matrix}")
        
        # Step 3: Eigendecomposition
        eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)
        
        # Sort by eigenvalue (descending)
        idx = eigenvalues.argsort()[::-1]
        eigenvalues = eigenvalues[idx]
        eigenvectors = eigenvectors[:, idx]
        
        print(f"\nEigenvalues (variance explained):")
        total_var = eigenvalues.sum()
        for i, λ in enumerate(eigenvalues):
            var_explained = λ / total_var * 100
            print(f"  PC{i+1}: λ={λ:.3f} ({var_explained:.1f}% variance)")
        
        print(f"\nPrincipal Components (eigenvectors):")
        for i in range(len(eigenvalues)):
            print(f"  PC{i+1}: {eigenvectors[:, i]}")
        
        # Step 4: Transform data
        data_pca = data_centered @ eigenvectors
        
        print(f"\nTransformed data shape: {data_pca.shape}")
        print(f"First principal component captures {eigenvalues[0]/total_var*100:.1f}% of variance!")
        
        # Dimensionality reduction
        # Keep only first component
        data_reduced = data_pca[:, :1] @ eigenvectors[:, :1].T
        
        print(f"\n🔽 Dimensionality Reduction:")
        print(f"  Original dimensions: 2")
        print(f"  Reduced dimensions: 1")
        print(f"  Information retained: {eigenvalues[0]/total_var*100:.1f}%")
    
    def power_iteration(self):
        """Finding dominant eigenvalue using power iteration"""
        print("\n🔄 Power Iteration Method")
        print("=" * 60)
        
        # Create a matrix
        A = np.array([[2, -1, 0],
                      [-1, 2, -1],
                      [0, -1, 2]])
        
        print(f"Matrix A:\n{A}")
        
        # Power iteration
        v = np.random.randn(3)  # Random starting vector
        v = v / np.linalg.norm(v)  # Normalize
        
        print("\nPower Iteration:")
        for iteration in range(10):
            # Multiply by A
            v_new = A @ v
            
            # Estimate eigenvalue (Rayleigh quotient)
            eigenvalue_estimate = (v_new @ v) / (v @ v)
            
            # Normalize
            v = v_new / np.linalg.norm(v_new)
            
            if iteration % 3 == 0:
                print(f"  Iteration {iteration}: λ ≈ {eigenvalue_estimate:.6f}")
        
        print(f"\nFinal estimate:")
        print(f"  Dominant eigenvalue: {eigenvalue_estimate:.6f}")
        print(f"  Corresponding eigenvector: {v}")
        
        # Compare with NumPy
        eigenvalues, eigenvectors = np.linalg.eig(A)
        max_idx = np.argmax(np.abs(eigenvalues))
        
        print(f"\nNumPy result:")
        print(f"  Dominant eigenvalue: {eigenvalues[max_idx]:.6f}")
        print(f"  Match? {np.abs(eigenvalue_estimate - eigenvalues[max_idx]) < 0.001}")
    
    def spectral_clustering_demo(self):
        """Using eigenvalues for clustering"""
        print("\n🎯 Spectral Clustering with Eigenvalues")
        print("=" * 60)
        
        # Create a simple graph (adjacency matrix)
        # Two clusters with weak connection
        adjacency = np.array([
            [0, 1, 1, 0, 0, 0],  # Cluster 1
            [1, 0, 1, 0, 0, 0],
            [1, 1, 0, 0.1, 0, 0],  # Weak link
            [0, 0, 0.1, 0, 1, 1],  # Cluster 2
            [0, 0, 0, 1, 0, 1],
            [0, 0, 0, 1, 1, 0]
        ])
        
        print("Graph adjacency matrix (connection strengths):")
        print(adjacency)
        
        # Compute degree matrix
        degree = np.diag(adjacency.sum(axis=1))
        
        # Laplacian matrix
        laplacian = degree - adjacency
        
        # Eigendecomposition of Laplacian
        eigenvalues, eigenvectors = np.linalg.eig(laplacian)
        
        # Sort by eigenvalue
        idx = eigenvalues.argsort()
        eigenvalues = eigenvalues[idx]
        eigenvectors = eigenvectors[:, idx]
        
        print(f"\n📊 Laplacian Eigenvalues:")
        for i, λ in enumerate(eigenvalues[:4]):
            print(f"  λ{i}: {λ:.3f}")
        
        print(f"\n🔍 Fiedler Vector (2nd smallest eigenvector):")
        fiedler = eigenvectors[:, 1]
        print(f"  {fiedler}")
        
        # Cluster based on sign of Fiedler vector
        cluster1 = np.where(fiedler < 0)[0]
        cluster2 = np.where(fiedler >= 0)[0]
        
        print(f"\n✨ Clustering Result:")
        print(f"  Cluster 1: Nodes {cluster1 + 1}")
        print(f"  Cluster 2: Nodes {cluster2 + 1}")
        print(f"  Perfect separation of the two groups!")

# Explore eigenvalues and eigenvectors
eigen = EigenExplorer()
eigen.eigen_basics()
eigen.pca_with_eigen()
eigen.power_iteration()
eigen.spectral_clustering_demo()

Matrix Decompositions: Breaking Down Complexity 🧩

Matrix decompositions are like taking apart a complex machine to understand its components. Just as you might disassemble a clock to see the gears, springs, and wheels, we decompose matrices to reveal their fundamental structures.

import numpy as np

class MatrixDecompositions:
    """
    Matrix decompositions: The power tools of linear algebra.
    Breaking down complexity into manageable pieces!
    """
    
    def __init__(self):
        np.random.seed(42)
    
    def lu_decomposition(self):
        """LU Decomposition for solving linear systems"""
        print("📊 LU Decomposition")
        print("=" * 60)
        
        # Matrix to decompose
        A = np.array([[2, 1, 1],
                      [4, -6, 0],
                      [-2, 7, 2]], dtype=float)
        
        print(f"Original matrix A:\n{A}")
        
        # LU decomposition using scipy
        from scipy.linalg import lu
        
        P, L, U = lu(A)
        
        print(f"\nP (Permutation matrix):\n{P}")
        print(f"\nL (Lower triangular):\n{L}")
        print(f"\nU (Upper triangular):\n{U}")
        
        # Verify: P @ A = L @ U
        reconstructed = L @ U
        original_permuted = P @ A
        
        print(f"\nVerification:")
        print(f"  P @ A =\n{original_permuted}")
        print(f"  L @ U =\n{reconstructed}")
        print(f"  Match? {np.allclose(original_permuted, reconstructed)}")
        
        # Application: Solving Ax = b
        b = np.array([1, 2, 3])
        print(f"\n🎯 Solving Ax = b where b = {b}")
        
        # Forward substitution: Ly = Pb
        # Back substitution: Ux = y
        x = np.linalg.solve(A, b)
        
        print(f"  Solution x = {x}")
        print(f"  Verification: A @ x = {A @ x}")
        print(f"  Equals b? {np.allclose(A @ x, b)}")
    
    def qr_decomposition(self):
        """QR Decomposition for orthogonalization"""
        print("\n🎯 QR Decomposition")
        print("=" * 60)
        
        # Create a matrix
        A = np.array([[1, 2],
                      [3, 4],
                      [5, 6]])
        
        print(f"Matrix A (3×2):\n{A}")
        
        # QR decomposition
        Q, R = np.linalg.qr(A)
        
        print(f"\nQ (Orthogonal, 3×3):\n{Q}")
        print(f"\nR (Upper triangular, 3×2):\n{R}")
        
        # Properties of Q
        print(f"\n📐 Properties of Q:")
        print(f"  Q^T @ Q =\n{(Q.T @ Q).round(10)}")
        print(f"  Q is orthogonal? {np.allclose(Q.T @ Q, np.eye(3))}")
        
        # Reconstruction
        reconstructed = Q @ R
        print(f"\nReconstruction Q @ R:\n{reconstructed}")
        print(f"Match with A? {np.allclose(reconstructed, A)}")
        
        # Application: Least squares
        print(f"\n📈 Application: Least Squares Fitting")
        
        # Overdetermined system (more equations than unknowns)
        b = np.array([1, 2, 4])
        
        # Solve using QR
        # A @ x ≈ b  =>  Q @ R @ x = b  =>  R @ x = Q^T @ b
        x = np.linalg.solve(R[:2, :], (Q.T @ b)[:2])
        
        print(f"  Least squares solution: x = {x}")
        print(f"  Residual: ||Ax - b|| = {np.linalg.norm(A @ x - b):.4f}")
    
    def svd_decomposition(self):
        """SVD: The most important decomposition"""
        print("\n⭐ Singular Value Decomposition (SVD)")
        print("=" * 60)
        
        # Create a data matrix (e.g., user-movie ratings)
        # 5 users, 4 movies
        ratings = np.array([[5, 4, 0, 1],
                           [4, 5, 0, 1],
                           [1, 1, 5, 4],
                           [1, 0, 5, 4],
                           [0, 1, 4, 5]], dtype=float)
        
        print("User-Movie Rating Matrix:")
        print("(5 users × 4 movies)")
        print(ratings)
        
        # SVD decomposition
        U, s, Vt = np.linalg.svd(ratings)
        
        print(f"\nSVD Components:")
        print(f"  U shape: {U.shape} (user features)")
        print(f"  Σ values: {s.round(3)} (importance)")
        print(f"  V^T shape: {Vt.shape} (movie features)")
        
        # Interpretation
        print(f"\n📊 Singular Values (importance of each component):")
        total_variance = np.sum(s**2)
        for i, σ in enumerate(s):
            variance_explained = (σ**2 / total_variance) * 100
            print(f"  Component {i+1}: σ={σ:.3f} ({variance_explained:.1f}% variance)")
        
        # Low-rank approximation
        k = 2  # Keep only 2 components
        
        # Truncated SVD
        U_k = U[:, :k]
        s_k = s[:k]
        Vt_k = Vt[:k, :]
        
        # Reconstruct
        Σ_k = np.diag(s_k)
        ratings_approx = U_k @ Σ_k @ Vt_k
        
        print(f"\n🔽 Rank-{k} Approximation:")
        print(ratings_approx.round(2))
        
        print(f"\nApproximation quality:")
        print(f"  Frobenius norm error: {np.linalg.norm(ratings - ratings_approx):.3f}")
        print(f"  Compression ratio: {ratings.size}/{U_k.size + s_k.size + Vt_k.size:.0f} = "
              f"{ratings.size/(U_k.size + s_k.size + Vt_k.size):.1f}x")
        
        # Recommendation system insight
        print(f"\n💡 Recommendation Insights:")
        print(f"  User 1 and 2 are similar (both like movies 1 & 2)")
        print(f"  User 3, 4, 5 form another group (like movies 3 & 4)")
        print(f"  The rank-2 approximation captures these two taste groups!")
    
    def cholesky_decomposition(self):
        """Cholesky decomposition for positive definite matrices"""
        print("\n🔺 Cholesky Decomposition")
        print("=" * 60)
        
        # Create a positive definite matrix (covariance matrix)
        A_temp = np.random.randn(3, 3)
        A = A_temp @ A_temp.T  # Ensure positive definite
        
        print(f"Positive definite matrix A:\n{A}")
        
        # Cholesky decomposition: A = L @ L^T
        L = np.linalg.cholesky(A)
        
        print(f"\nL (Lower triangular):\n{L}")
        
        # Verify
        reconstructed = L @ L.T
        print(f"\nL @ L^T:\n{reconstructed}")
        print(f"Match with A? {np.allclose(reconstructed, A)}")
        
        # Application: Efficient random sampling
        print(f"\n🎲 Application: Multivariate Normal Sampling")
        
        # Generate samples from N(0, A)
        n_samples = 5
        z = np.random.randn(n_samples, 3)  # Standard normal
        samples = z @ L.T  # Transform to N(0, A)
        
        print(f"  Generated {n_samples} samples from N(0, A)")
        print(f"  Sample covariance:\n{np.cov(samples.T).round(2)}")
        print(f"  Should approximate A")

# Explore decompositions
decomp = MatrixDecompositions()
decomp.lu_decomposition()
decomp.qr_decomposition()
decomp.svd_decomposition()
decomp.cholesky_decomposition()

Practical Applications: Linear Algebra in Action 🚀

import numpy as np

# Complete Practical Example: Image Compression with SVD

def compress_image_with_svd():
    """
    Compress an image using SVD - real-world linear algebra!
    """
    print("🖼️ Image Compression with SVD")
    print("=" * 60)
    
    # Create a synthetic "image" (normally you'd load a real image)
    # Let's create a pattern
    size = 50
    x = np.linspace(-5, 5, size)
    y = np.linspace(-5, 5, size)
    X, Y = np.meshgrid(x, y)
    
    # Create an interesting pattern
    image = np.sin(np.sqrt(X**2 + Y**2)) + np.cos(2*X) * np.sin(2*Y)
    
    # Normalize to 0-255 range
    image = ((image - image.min()) / (image.max() - image.min()) * 255).astype(np.uint8)
    
    print(f"Original image: {size}×{size} pixels")
    print(f"Original size: {image.size} bytes")
    
    # Perform SVD
    U, s, Vt = np.linalg.svd(image, full_matrices=False)
    
    print(f"\nSingular values (first 10): {s[:10].round(1)}")
    
    # Compress by keeping only k components
    compression_levels = [5, 10, 20, 30]
    
    for k in compression_levels:
        # Reconstruct with k components
        compressed = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]
        
        # Calculate compression ratio
        original_size = image.size
        compressed_size = k * (U.shape[0] + Vt.shape[1] + 1)
        ratio = original_size / compressed_size
        
        # Calculate error
        mse = np.mean((image - compressed)**2)
        psnr = 10 * np.log10(255**2 / mse) if mse > 0 else float('inf')
        
        print(f"\nRank-{k} approximation:")
        print(f"  Compression ratio: {ratio:.1f}x")
        print(f"  PSNR: {psnr:.1f} dB")
        print(f"  Storage: {compressed_size} values (vs {original_size})")

compress_image_with_svd()

Summary: Your Linear Algebra Toolkit 🛠️

# Quick Linear Algebra Reference Card
import numpy as np

# BASIC OPERATIONS
# ----------------
# Vectors
v1 = np.array([1, 2, 3])
v2 = np.array([4, 5, 6])
dot_product = v1 @ v2  # or np.dot(v1, v2)
magnitude = np.linalg.norm(v1)

# Matrices
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
matrix_mult = A @ B  # or np.matmul(A, B)
transpose = A.T
inverse = np.linalg.inv(A)
determinant = np.linalg.det(A)

# SOLVING SYSTEMS
# ---------------
# Ax = b
x = np.linalg.solve(A, b)  # Direct solve
x = np.linalg.lstsq(A, b, rcond=None)[0]  # Least squares

# EIGENPROBLEMS
# -------------
eigenvalues, eigenvectors = np.linalg.eig(A)
# For symmetric matrices (faster, more stable)
eigenvalues, eigenvectors = np.linalg.eigh(A)

# DECOMPOSITIONS
# --------------
# SVD
U, s, Vt = np.linalg.svd(A)
# QR
Q, R = np.linalg.qr(A)
# Cholesky (for positive definite)
L = np.linalg.cholesky(A)

# COMMON PATTERNS
# ---------------
# PCA
cov_matrix = np.cov(data.T)
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)
principal_components = data @ eigenvectors

# Distance matrix
diff = points[:, np.newaxis, :] - points[np.newaxis, :, :]
distances = np.linalg.norm(diff, axis=2)

# Regularized regression
w = np.linalg.solve(X.T @ X + lambda * np.eye(n), X.T @ y)

print("🔢 Linear algebra: The foundation of data science!")