Chapter 14: Practical Applications of Linear Algebra

Principal Component Analysis (PCA)

Theoretical Foundation

Principal component analysis projects data onto a low-dimensional space through linear transformation while preserving as much variance as possible. The core idea is to perform eigenvalue decomposition on the covariance matrix.

Mathematical Principles

Given data matrix $X \in \mathbb{R}^{n \times p}$, the covariance matrix is: $C = \frac{1}{n-1}X^TX$

Principal components are eigenvectors of $C$, sorted by eigenvalue magnitude.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris, make_blobs
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import seaborn as sns

def pca_comprehensive_demo():
    """Comprehensive demonstration of principal component analysis"""

    print("Principal Component Analysis (PCA) Application Demonstration")
    print("=" * 50)

    # Generate example data
    np.random.seed(42)
    n_samples = 300

    # Create three correlated features
    X1 = np.random.randn(n_samples)
    X2 = 0.8 * X1 + 0.6 * np.random.randn(n_samples)
    X3 = -0.5 * X1 + 0.3 * X2 + 0.7 * np.random.randn(n_samples)
    X4 = 0.1 * np.random.randn(n_samples)  # Noise feature

    X = np.column_stack([X1, X2, X3, X4])
    feature_names = ['Feature 1', 'Feature 2', 'Feature 3', 'Feature 4']

    print(f"Data Shape: {X.shape}")
    print(f"Feature Statistics:")
    for i, name in enumerate(feature_names):
        print(f"  {name}: Mean={np.mean(X[:, i]):.3f}, Std={np.std(X[:, i]):.3f}")

    # Standardize data
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)

    print(f"\nStandardized Data Statistics:")
    for i, name in enumerate(feature_names):
        print(f"  {name}: Mean={np.mean(X_scaled[:, i]):.3f}, Std={np.std(X_scaled[:, i]):.3f}")

    # Manual PCA implementation
    print(f"\nManual PCA Implementation:")

    # Compute covariance matrix
    cov_matrix = np.cov(X_scaled.T)
    print(f"Covariance Matrix:\n{cov_matrix}")

    # Eigenvalue decomposition
    eigenvals, eigenvecs = np.linalg.eig(cov_matrix)

    # Sort
    idx = np.argsort(eigenvals)[::-1]
    eigenvals = eigenvals[idx]
    eigenvecs = eigenvecs[:, idx]

    print(f"\nEigenvalues: {eigenvals}")
    print(f"Variance Explained Ratio: {eigenvals / np.sum(eigenvals)}")
    print(f"Cumulative Variance Explained: {np.cumsum(eigenvals / np.sum(eigenvals))}")

    # Principal component transformation
    X_pca_manual = X_scaled @ eigenvecs

    # Verify with sklearn PCA
    pca_sklearn = PCA()
    X_pca_sklearn = pca_sklearn.fit_transform(X_scaled)

    print(f"\nVerification with sklearn:")
    print(f"Eigenvalues Match: {np.allclose(eigenvals, pca_sklearn.explained_variance_)}")
    print(f"Principal Components Match: {np.allclose(np.abs(X_pca_manual), np.abs(X_pca_sklearn))}")

    # Visualization
    fig, axes = plt.subplots(2, 3, figsize=(18, 12))

    # Original data correlation
    ax1 = axes[0, 0]
    correlation_matrix = np.corrcoef(X_scaled.T)
    im = ax1.imshow(correlation_matrix, cmap='RdBu', vmin=-1, vmax=1)
    ax1.set_xticks(range(4))
    ax1.set_yticks(range(4))
    ax1.set_xticklabels(feature_names, rotation=45)
    ax1.set_yticklabels(feature_names)
    ax1.set_title('Original Feature Correlation Matrix')

    # Add numerical annotations
    for i in range(4):
        for j in range(4):
            ax1.text(j, i, f'{correlation_matrix[i,j]:.2f}',
                    ha='center', va='center', color='white' if abs(correlation_matrix[i,j]) > 0.5 else 'black')

    plt.colorbar(im, ax=ax1)

    # Eigenvalues and variance explained
    ax2 = axes[0, 1]
    x_pos = range(1, 5)
    bars = ax2.bar(x_pos, eigenvals, alpha=0.7, color='skyblue')
    ax2.set_xlabel('Principal Component')
    ax2.set_ylabel('Eigenvalue (Variance)')
    ax2.set_title('Variance Contribution of Each Principal Component')
    ax2.set_xticks(x_pos)
    ax2.set_xticklabels([f'PC{i}' for i in x_pos])

    # Add variance explained ratio labels
    for bar, val, ratio in zip(bars, eigenvals, eigenvals/np.sum(eigenvals)):
        ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.02,
                f'{ratio:.1%}', ha='center', va='bottom')

    ax2.grid(True, alpha=0.3)

    # Cumulative variance explained
    ax3 = axes[0, 2]
    cumsum_ratio = np.cumsum(eigenvals / np.sum(eigenvals))
    ax3.plot(x_pos, cumsum_ratio, 'o-', linewidth=2, markersize=8, color='red')
    ax3.axhline(y=0.95, color='green', linestyle='--', label='95% Threshold')
    ax3.set_xlabel('Number of Principal Components')
    ax3.set_ylabel('Cumulative Variance Explained Ratio')
    ax3.set_title('Cumulative Variance Explained Ratio')
    ax3.set_xticks(x_pos)
    ax3.set_xticklabels([f'PC{i}' for i in x_pos])
    ax3.grid(True, alpha=0.3)
    ax3.legend()

    # Data distribution in first two principal components
    ax4 = axes[1, 0]
    scatter = ax4.scatter(X_pca_manual[:, 0], X_pca_manual[:, 1],
                         alpha=0.6, c=X_scaled[:, 0], cmap='viridis')
    ax4.set_xlabel(f'PC1 ({eigenvals[0]/np.sum(eigenvals)*100:.1f}%)')
    ax4.set_ylabel(f'PC2 ({eigenvals[1]/np.sum(eigenvals)*100:.1f}%)')
    ax4.set_title('Data Distribution in First Two Principal Components')
    ax4.grid(True, alpha=0.3)
    plt.colorbar(scatter, ax=ax4, label='Feature 1 value')

    # Principal component loadings plot
    ax5 = axes[1, 1]
    for i, name in enumerate(feature_names):
        ax5.arrow(0, 0, eigenvecs[i, 0]*2, eigenvecs[i, 1]*2,
                 head_width=0.1, head_length=0.1, fc=f'C{i}', ec=f'C{i}',
                 linewidth=2, label=name)

    ax5.set_xlim(-1, 1)
    ax5.set_ylim(-1, 1)
    ax5.set_xlabel('PC1 Loading')
    ax5.set_ylabel('PC2 Loading')
    ax5.set_title('Principal Component Loadings Plot')
    ax5.legend()
    ax5.grid(True, alpha=0.3)
    ax5.set_aspect('equal')

    # Reconstruction error with dimensionality reduction
    ax6 = axes[1, 2]
    reconstruction_errors = []

    for n_components in range(1, 5):
        # Reconstruct using first n principal components
        X_reconstructed = X_pca_manual[:, :n_components] @ eigenvecs[:, :n_components].T
        error = np.mean(np.sum((X_scaled - X_reconstructed)**2, axis=1))
        reconstruction_errors.append(error)

    ax6.plot(range(1, 5), reconstruction_errors, 'o-', linewidth=2, markersize=8)
    ax6.set_xlabel('Number of Principal Components Used')
    ax6.set_ylabel('Mean Reconstruction Error')
    ax6.set_title('Dimensionality Reduction Reconstruction Error')
    ax6.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()

    return X_scaled, X_pca_manual, eigenvals, eigenvecs

def pca_real_world_example():
    """PCA application on real data"""

    print(f"\nReal Data Example: Iris Dataset")
    print("=" * 40)

    # Load iris dataset
    iris = load_iris()
    X = iris.data
    y = iris.target
    feature_names = iris.feature_names
    target_names = iris.target_names

    print(f"Data Shape: {X.shape}")
    print(f"Feature Names: {feature_names}")
    print(f"Classes: {target_names}")

    # Standardize
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)

    # PCA
    pca = PCA()
    X_pca = pca.fit_transform(X_scaled)

    print(f"\nVariance Explained Ratio for Each Principal Component:")
    for i, ratio in enumerate(pca.explained_variance_ratio_):
        print(f"PC{i+1}: {ratio:.3f} ({ratio*100:.1f}%)")

    print(f"Cumulative Variance Explained: {np.cumsum(pca.explained_variance_ratio_)}")

    # Visualization
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))

    # First two features of original data
    ax1 = axes[0, 0]
    colors = ['red', 'green', 'blue']
    for i, (color, target) in enumerate(zip(colors, target_names)):
        mask = y == i
        ax1.scatter(X_scaled[mask, 0], X_scaled[mask, 1],
                   c=color, label=target, alpha=0.7)

    ax1.set_xlabel(f'{feature_names[0]} (Standardized)')
    ax1.set_ylabel(f'{feature_names[1]} (Standardized)')
    ax1.set_title('Original Feature Space')
    ax1.legend()
    ax1.grid(True, alpha=0.3)

    # First two principal components after PCA
    ax2 = axes[0, 1]
    for i, (color, target) in enumerate(zip(colors, target_names)):
        mask = y == i
        ax2.scatter(X_pca[mask, 0], X_pca[mask, 1],
                   c=color, label=target, alpha=0.7)

    ax2.set_xlabel(f'PC1 ({pca.explained_variance_ratio_[0]*100:.1f}%)')
    ax2.set_ylabel(f'PC2 ({pca.explained_variance_ratio_[1]*100:.1f}%)')
    ax2.set_title('Principal Component Space')
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    # Principal component loadings
    ax3 = axes[1, 0]
    loadings = pca.components_[:2].T * np.sqrt(pca.explained_variance_[:2])

    for i, feature in enumerate(feature_names):
        ax3.arrow(0, 0, loadings[i, 0], loadings[i, 1],
                 head_width=0.1, head_length=0.1, fc=f'C{i}', ec=f'C{i}',
                 linewidth=2)
        ax3.text(loadings[i, 0]*1.1, loadings[i, 1]*1.1, feature,
                fontsize=10, ha='center', va='center')

    ax3.set_xlim(-3, 3)
    ax3.set_ylim(-3, 3)
    ax3.set_xlabel('PC1 Loading')
    ax3.set_ylabel('PC2 Loading')
    ax3.set_title('Feature Loadings Plot')
    ax3.grid(True, alpha=0.3)
    ax3.axhline(0, color='black', linewidth=0.5)
    ax3.axvline(0, color='black', linewidth=0.5)

    # Variance explained ratio
    ax4 = axes[1, 1]
    x_pos = range(1, 5)
    bars = ax4.bar(x_pos, pca.explained_variance_ratio_, alpha=0.7)
    ax4.set_xlabel('Principal Component')
    ax4.set_ylabel('Variance Explained Ratio')
    ax4.set_title('Variance Contribution of Each Principal Component')
    ax4.set_xticks(x_pos)

    # Cumulative line
    ax4_twin = ax4.twinx()
    ax4_twin.plot(x_pos, np.cumsum(pca.explained_variance_ratio_),
                 'ro-', linewidth=2, label='Cumulative')
    ax4_twin.set_ylabel('Cumulative Variance Explained Ratio')
    ax4_twin.legend()

    plt.tight_layout()
    plt.show()

# Execute PCA demonstration
X_scaled, X_pca, eigenvals, eigenvecs = pca_comprehensive_demo()
pca_real_world_example()

Linear Algebra in Graph Theory

Adjacency Matrix and Laplacian Matrix

Graph structure can be represented using matrices, and graph properties are closely related to matrix eigenvalues.

PageRank Algorithm

Google’s PageRank algorithm is essentially an eigenvector problem.

import networkx as nx
from scipy.sparse import csc_matrix
from scipy.sparse.linalg import eigs

def graph_theory_applications():
    """Linear algebra applications in graph theory"""

    print("Linear Algebra Applications in Graph Theory")
    print("=" * 50)

    # Create example graph
    G = nx.Graph()
    edges = [(1, 2), (1, 3), (2, 3), (2, 4), (3, 4), (3, 5), (4, 5), (4, 6), (5, 6)]
    G.add_edges_from(edges)

    print(f"Graph Basic Information:")
    print(f"Number of Nodes: {G.number_of_nodes()}")
    print(f"Number of Edges: {G.number_of_edges()}")

    # Adjacency matrix
    A = nx.adjacency_matrix(G).toarray()
    print(f"\nAdjacency Matrix:")
    print(A)

    # Degree matrix
    degrees = dict(G.degree())
    D = np.diag([degrees[i] for i in range(1, 7)])
    print(f"\nDegree Matrix:")
    print(D)

    # Laplacian matrix
    L = D - A
    print(f"\nLaplacian Matrix (L = D - A):")
    print(L)

    # Compute eigenvalues of Laplacian matrix
    eigenvals_L, eigenvecs_L = np.linalg.eig(L)
    eigenvals_L = np.real(eigenvals_L)
    idx = np.argsort(eigenvals_L)
    eigenvals_L = eigenvals_L[idx]
    eigenvecs_L = eigenvecs_L[:, idx]

    print(f"\nEigenvalues of Laplacian Matrix:")
    print(eigenvals_L)

    # Second smallest eigenvalue (algebraic connectivity)
    algebraic_connectivity = eigenvals_L[1]
    print(f"Algebraic Connectivity (Second Smallest Eigenvalue): {algebraic_connectivity:.6f}")

    # Graph connectivity
    num_components = np.sum(eigenvals_L < 1e-10)
    print(f"Number of Connected Components: {num_components}")

    # Visualization
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))

    # Graph visualization
    ax1 = axes[0, 0]
    pos = nx.spring_layout(G, seed=42)
    nx.draw(G, pos, ax=ax1, with_labels=True, node_color='lightblue',
            node_size=500, font_size=10, font_weight='bold')
    ax1.set_title('Original Graph')

    # Adjacency matrix heatmap
    ax2 = axes[0, 1]
    im2 = ax2.imshow(A, cmap='Blues')
    ax2.set_title('Adjacency Matrix')
    ax2.set_xticks(range(6))
    ax2.set_yticks(range(6))
    ax2.set_xticklabels(range(1, 7))
    ax2.set_yticklabels(range(1, 7))
    plt.colorbar(im2, ax=ax2)

    # Laplacian matrix heatmap
    ax3 = axes[0, 2]
    im3 = ax3.imshow(L, cmap='RdBu')
    ax3.set_title('Laplacian Matrix')
    ax3.set_xticks(range(6))
    ax3.set_yticks(range(6))
    ax3.set_xticklabels(range(1, 7))
    ax3.set_yticklabels(range(1, 7))
    plt.colorbar(im3, ax=ax3)

    # Eigenvalue spectrum
    ax4 = axes[1, 0]
    ax4.bar(range(len(eigenvals_L)), eigenvals_L, alpha=0.7)
    ax4.set_xlabel('Eigenvalue Index')
    ax4.set_ylabel('Eigenvalue')
    ax4.set_title('Laplacian Matrix Eigenvalue Spectrum')
    ax4.grid(True, alpha=0.3)

    # Fiedler vector (second smallest eigenvector)
    ax5 = axes[1, 1]
    fiedler_vector = eigenvecs_L[:, 1]
    colors = ['red' if x < 0 else 'blue' for x in fiedler_vector]

    pos_fixed = {i+1: pos[i+1] for i in range(6)}
    node_colors = ['red' if fiedler_vector[i] < 0 else 'blue' for i in range(6)]

    nx.draw(G, pos_fixed, ax=ax5, with_labels=True, node_color=node_colors,
            node_size=500, font_size=10, font_weight='bold')
    ax5.set_title('Fiedler Vector Partition')

    # Eigenvector analysis
    ax6 = axes[1, 2]
    for i in range(min(3, len(eigenvals_L))):
        ax6.plot(eigenvecs_L[:, i], 'o-', label=f'Eigenvector {i+1} (λ={eigenvals_L[i]:.3f})')

    ax6.set_xlabel('Node')
    ax6.set_ylabel('Eigenvector Value')
    ax6.set_title('First Three Eigenvectors')
    ax6.legend()
    ax6.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()

def pagerank_demo():
    """PageRank algorithm demonstration"""

    print(f"\nPageRank Algorithm Demonstration")
    print("=" * 30)

    # Create directed graph
    G = nx.DiGraph()
    edges = [(1, 2), (1, 3), (2, 3), (3, 1), (3, 4), (4, 3), (4, 5), (5, 4), (5, 6), (6, 5)]
    G.add_edges_from(edges)

    print(f"Directed Graph Information:")
    print(f"Number of Nodes: {G.number_of_nodes()}")
    print(f"Number of Edges: {G.number_of_edges()}")

    # Construct transition matrix
    A = nx.adjacency_matrix(G).toarray().astype(float)

    # Column normalization (each column sums to 1)
    out_degrees = np.sum(A, axis=0)
    for i in range(len(out_degrees)):
        if out_degrees[i] > 0:
            A[:, i] = A[:, i] / out_degrees[i]

    # PageRank matrix (add random jump)
    damping_factor = 0.85
    n = G.number_of_nodes()
    M = damping_factor * A + (1 - damping_factor) / n * np.ones((n, n))

    print(f"\nPageRank Transition Matrix:")
    print(M)

    # Compute dominant eigenvector (corresponding to eigenvalue 1)
    eigenvals, eigenvecs = np.linalg.eig(M)

    # Find eigenvalue closest to 1
    idx = np.argmin(np.abs(eigenvals - 1))
    pagerank_vector = np.real(eigenvecs[:, idx])
    pagerank_vector = pagerank_vector / np.sum(pagerank_vector)  # Normalize

    print(f"\nManually Computed PageRank Values:")
    for i in range(n):
        print(f"Node {i+1}: {pagerank_vector[i]:.6f}")

    # Verify with NetworkX
    pagerank_nx = nx.pagerank(G, alpha=damping_factor)

    print(f"\nNetworkX Computed PageRank Values:")
    for node in sorted(pagerank_nx.keys()):
        print(f"Node {node}: {pagerank_nx[node]:.6f}")

    # Power method for solving PageRank
    def power_method_pagerank(M, max_iterations=100, tolerance=1e-6):
        n = M.shape[0]
        x = np.ones(n) / n  # Initialize

        for i in range(max_iterations):
            x_new = M @ x

            if np.linalg.norm(x_new - x) < tolerance:
                break

            x = x_new

        return x, i+1

    pr_power, iterations = power_method_pagerank(M)
    print(f"\nPower Method Computed PageRank Values ({iterations} iterations):")
    for i in range(n):
        print(f"Node {i+1}: {pr_power[i]:.6f}")

    # Visualization
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))

    # Directed graph
    ax1 = axes[0]
    pos = nx.spring_layout(G, seed=42)
    nx.draw(G, pos, ax=ax1, with_labels=True, node_color='lightgreen',
            node_size=1000, font_size=12, font_weight='bold',
            arrows=True, arrowsize=20, arrowstyle='->')
    ax1.set_title('Directed Graph Structure')

    # PageRank value comparison
    ax2 = axes[1]
    nodes = list(range(1, n+1))
    width = 0.25

    ax2.bar([x - width for x in nodes], pagerank_vector, width,
            label='Manual Calculation', alpha=0.7)
    ax2.bar(nodes, [pagerank_nx[i] for i in nodes], width,
            label='NetworkX', alpha=0.7)
    ax2.bar([x + width for x in nodes], pr_power, width,
            label='Power Method', alpha=0.7)

    ax2.set_xlabel('Node')
    ax2.set_ylabel('PageRank Value')
    ax2.set_title('PageRank Value Comparison')
    ax2.set_xticks(nodes)
    ax2.legend()
    ax2.grid(True, alpha=0.3)

    # Graph with node sizes adjusted by PageRank values
    ax3 = axes[2]
    node_sizes = [pr_power[i-1] * 5000 for i in G.nodes()]  # Magnify for display
    nx.draw(G, pos, ax=ax3, with_labels=True,
            node_size=node_sizes, font_size=12, font_weight='bold',
            arrows=True, arrowsize=20, arrowstyle='->',
            node_color='lightcoral')
    ax3.set_title('Node Size Reflects PageRank Value')

    plt.tight_layout()
    plt.show()

# Execute graph theory demonstration
graph_theory_applications()
pagerank_demo()

Markov Chains

Transition Matrix

State transitions of a Markov chain can be represented by a matrix, and the steady-state distribution is the dominant eigenvector of the transition matrix.

Steady-State Analysis

The long-term behavior of the system is determined by the eigenvalues and eigenvectors of the transition matrix.

def markov_chain_analysis():
    """Matrix analysis of Markov chains"""

    print("Markov Chain Matrix Analysis")
    print("=" * 50)

    # Weather forecasting model
    # States: 0=Sunny, 1=Cloudy, 2=Rainy
    weather_states = ['Sunny', 'Cloudy', 'Rainy']

    # Transition probability matrix P[i,j] = probability of transitioning from state i to state j
    P = np.array([[0.7, 0.2, 0.1],    # Sunny -> Sunny/Cloudy/Rainy
                  [0.3, 0.4, 0.3],    # Cloudy -> Sunny/Cloudy/Rainy
                  [0.2, 0.3, 0.5]])   # Rainy -> Sunny/Cloudy/Rainy

    print("Transition Probability Matrix:")
    print("     ", "  ".join(f"{state:>6}" for state in weather_states))
    for i, from_state in enumerate(weather_states):
        print(f"{from_state:>6}: {P[i]}")

    # Verify probability matrix (each row sums to 1)
    row_sums = np.sum(P, axis=1)
    print(f"\nRow Probability Sums: {row_sums}")
    print(f"Is Valid Probability Matrix: {np.allclose(row_sums, 1.0)}")

    # Compute n-step transition probabilities
    print(f"\nMulti-step Transition Probabilities:")
    for n in [2, 5, 10, 20]:
        P_n = np.linalg.matrix_power(P, n)
        print(f"\n{n}-step Transition Probability Matrix:")
        print(P_n)

    # Steady-state distribution (dominant eigenvector)
    eigenvals, eigenvecs = np.linalg.eig(P.T)  # Note the transpose

    # Find eigenvector corresponding to eigenvalue 1
    idx = np.argmin(np.abs(eigenvals - 1.0))
    steady_state = np.real(eigenvecs[:, idx])
    steady_state = steady_state / np.sum(steady_state)  # Normalize to probability

    print(f"\nSteady-State Distribution:")
    for i, (state, prob) in enumerate(zip(weather_states, steady_state)):
        print(f"{state}: {prob:.6f} ({prob*100:.2f}%)")

    # Verify steady-state property
    steady_check = P.T @ steady_state
    print(f"\nSteady-State Verification (Pᵀπ = π):")
    print(f"Error: {np.linalg.norm(steady_check - steady_state):.2e}")

    # Simulate Markov chain
    def simulate_markov_chain(P, initial_state, n_steps):
        """Simulate Markov chain"""
        states = [initial_state]
        current_state = initial_state

        for _ in range(n_steps):
            # Choose next state according to transition probabilities
            current_state = np.random.choice(3, p=P[current_state])
            states.append(current_state)

        return np.array(states)

    # Run multiple simulations
    np.random.seed(42)
    n_simulations = 1000
    n_steps = 100

    all_final_states = []
    for _ in range(n_simulations):
        # Random initial state
        initial = np.random.choice(3)
        trajectory = simulate_markov_chain(P, initial, n_steps)
        all_final_states.append(trajectory[-1])

    # Compute empirical steady-state distribution
    empirical_steady = np.bincount(all_final_states) / len(all_final_states)

    print(f"\nEmpirical Steady-State Distribution ({n_simulations} simulations):")
    for i, (state, prob) in enumerate(zip(weather_states, empirical_steady)):
        print(f"{state}: {prob:.6f} ({prob*100:.2f}%)")

    print(f"\nTheoretical vs Empirical Steady-State Distribution Error: {np.linalg.norm(steady_state - empirical_steady):.6f}")

    # Visualization
    fig, axes = plt.subplots(2, 3, figsize=(18, 12))

    # Transition probability matrix heatmap
    ax1 = axes[0, 0]
    im1 = ax1.imshow(P, cmap='Blues', vmin=0, vmax=1)
    ax1.set_xticks(range(3))
    ax1.set_yticks(range(3))
    ax1.set_xticklabels(weather_states, rotation=45)
    ax1.set_yticklabels(weather_states)
    ax1.set_title('Transition Probability Matrix')

    # Add numerical annotations
    for i in range(3):
        for j in range(3):
            ax1.text(j, i, f'{P[i,j]:.2f}', ha='center', va='center',
                    color='white' if P[i,j] > 0.5 else 'black')

    plt.colorbar(im1, ax=ax1)

    # Eigenvalues
    ax2 = axes[0, 1]
    eigenvals_magnitude = np.abs(eigenvals)
    ax2.bar(range(len(eigenvals)), eigenvals_magnitude, alpha=0.7)
    ax2.set_xlabel('Eigenvalue Index')
    ax2.set_ylabel('Eigenvalue Magnitude')
    ax2.set_title('Transition Matrix Eigenvalues')
    ax2.grid(True, alpha=0.3)

    # Steady-state distribution comparison
    ax3 = axes[0, 2]
    x_pos = range(3)
    width = 0.35

    ax3.bar([x - width/2 for x in x_pos], steady_state, width,
            label='Theoretical Steady-State', alpha=0.7)
    ax3.bar([x + width/2 for x in x_pos], empirical_steady, width,
            label='Empirical Steady-State', alpha=0.7)

    ax3.set_xlabel('Weather State')
    ax3.set_ylabel('Probability')
    ax3.set_title('Steady-State Distribution Comparison')
    ax3.set_xticks(x_pos)
    ax3.set_xticklabels(weather_states)
    ax3.legend()
    ax3.grid(True, alpha=0.3)

    # Convergence process
    ax4 = axes[1, 0]
    initial_dist = np.array([1, 0, 0])  # Start from sunny
    distributions = [initial_dist]

    for n in range(1, 21):
        dist = initial_dist @ np.linalg.matrix_power(P, n)
        distributions.append(dist)

    distributions = np.array(distributions)

    for i, state in enumerate(weather_states):
        ax4.plot(range(21), distributions[:, i], 'o-', label=state)

    ax4.axhline(y=steady_state[0], color='C0', linestyle='--', alpha=0.7)
    ax4.axhline(y=steady_state[1], color='C1', linestyle='--', alpha=0.7)
    ax4.axhline(y=steady_state[2], color='C2', linestyle='--', alpha=0.7)

    ax4.set_xlabel('Time Step')
    ax4.set_ylabel('Probability')
    ax4.set_title('Distribution Convergence Process')
    ax4.legend()
    ax4.grid(True, alpha=0.3)

    # Single simulation trajectory
    ax5 = axes[1, 1]
    trajectory = simulate_markov_chain(P, 0, 50)
    ax5.plot(trajectory, 'o-', markersize=4)
    ax5.set_xlabel('Time')
    ax5.set_ylabel('State')
    ax5.set_title('Markov Chain Trajectory Example')
    ax5.set_yticks(range(3))
    ax5.set_yticklabels(weather_states)
    ax5.grid(True, alpha=0.3)

    # State occupancy time statistics
    ax6 = axes[1, 2]
    long_trajectory = simulate_markov_chain(P, 0, 10000)
    state_counts = np.bincount(long_trajectory[1000:])  # Remove first 1000 steps
    state_frequencies = state_counts / len(long_trajectory[1000:])

    ax6.bar(range(3), state_frequencies, alpha=0.7, label='Long-term Simulation')
    ax6.bar(range(3), steady_state, alpha=0.5, label='Theoretical Steady-State', width=0.5)

    ax6.set_xlabel('Weather State')
    ax6.set_ylabel('Frequency/Probability')
    ax6.set_title('Long-term State Distribution')
    ax6.set_xticks(range(3))
    ax6.set_xticklabels(weather_states)
    ax6.legend()
    ax6.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()

    return P, steady_state

# Execute Markov chain demonstration
P, steady_state = markov_chain_analysis()

Linear Algebra in Computer Graphics

3D Transformations

Transformations such as rotation, scaling, and translation can all be represented by matrices.

Projection Transformations

Projection from 3D to 2D is an important application of linear transformations.

from mpl_toolkits.mplot3d.art3d import Poly3DCollection

def computer_graphics_demo():
    """Linear algebra applications in computer graphics"""

    print("Linear Algebra in Computer Graphics")
    print("=" * 50)

    # Define 3D cube vertices
    cube_vertices = np.array([
        [0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0],  # Bottom face
        [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]   # Top face
    ])

    print("Cube Vertex Coordinates:")
    for i, vertex in enumerate(cube_vertices):
        print(f"Vertex {i}: {vertex}")

    # Define cube faces (each face represented by 4 vertex indices)
    cube_faces = [
        [0, 1, 2, 3],  # Bottom face
        [4, 5, 6, 7],  # Top face
        [0, 1, 5, 4],  # Front face
        [2, 3, 7, 6],  # Back face
        [0, 3, 7, 4],  # Left face
        [1, 2, 6, 5]   # Right face
    ]

    # Transformation matrix definitions
    def rotation_x(angle):
        """Rotation around X-axis"""
        c, s = np.cos(angle), np.sin(angle)
        return np.array([[1, 0, 0], [0, c, -s], [0, s, c]])

    def rotation_y(angle):
        """Rotation around Y-axis"""
        c, s = np.cos(angle), np.sin(angle)
        return np.array([[c, 0, s], [0, 1, 0], [-s, 0, c]])

    def rotation_z(angle):
        """Rotation around Z-axis"""
        c, s = np.cos(angle), np.sin(angle)
        return np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]])

    def scaling(sx, sy, sz):
        """Scaling transformation"""
        return np.array([[sx, 0, 0], [0, sy, 0], [0, 0, sz]])

    def translation(tx, ty, tz):
        """Translation transformation (homogeneous coordinates)"""
        return np.array([[1, 0, 0, tx],
                        [0, 1, 0, ty],
                        [0, 0, 1, tz],
                        [0, 0, 0, 1]])

    def perspective_projection(d):
        """Perspective projection"""
        return np.array([[1, 0, 0, 0],
                        [0, 1, 0, 0],
                        [0, 0, 1, 0],
                        [0, 0, 1/d, 0]])

    # Apply various transformations
    print(f"\nApplying Transformations:")

    # 1. Rotation transformation
    angle = np.pi / 6  # 30 degrees
    R_x = rotation_x(angle)
    R_y = rotation_y(angle)
    R_z = rotation_z(angle)

    print(f"Rotation around X-axis by 30 degrees matrix:")
    print(R_x)

    # Combined rotation
    R_combined = R_z @ R_y @ R_x
    rotated_cube = (R_combined @ cube_vertices.T).T

    # 2. Scaling transformation
    S = scaling(1.5, 0.8, 1.2)
    scaled_cube = (S @ cube_vertices.T).T

    # 3. Homogeneous coordinate transformation (rotation + translation)
    # Convert 3D points to homogeneous coordinates
    vertices_homogeneous = np.column_stack([cube_vertices, np.ones(8)])

    # Combined transformation: rotate then translate
    T = translation(2, 1, 0.5)
    R_homo = np.eye(4)
    R_homo[:3, :3] = R_combined

    transform_combined = T @ R_homo
    transformed_vertices = (transform_combined @ vertices_homogeneous.T).T
    transformed_cube = transformed_vertices[:, :3]  # Remove homogeneous coordinate

    print(f"\nCombined Transformation Matrix (Rotate then Translate):")
    print(transform_combined)

    # Visualization
    fig = plt.figure(figsize=(20, 15))

    def draw_cube(ax, vertices, title, color='blue', alpha=0.7):
        """Draw cube"""
        # Create collection of faces
        faces = []
        for face in cube_faces:
            face_vertices = vertices[face]
            faces.append(face_vertices)

        # Draw cube faces
        poly3d = [[vertices[j] for j in face] for face in cube_faces]
        ax.add_collection3d(Poly3DCollection(poly3d, alpha=alpha,
                                           facecolor=color, edgecolor='black'))

        # Set coordinate axes
        all_coords = vertices.flatten()
        ax.set_xlim([all_coords.min()-0.5, all_coords.max()+0.5])
        ax.set_ylim([all_coords.min()-0.5, all_coords.max()+0.5])
        ax.set_zlim([all_coords.min()-0.5, all_coords.max()+0.5])

        ax.set_xlabel('X')
        ax.set_ylabel('Y')
        ax.set_zlabel('Z')
        ax.set_title(title)

    # Original cube
    ax1 = fig.add_subplot(2, 4, 1, projection='3d')
    draw_cube(ax1, cube_vertices, 'Original Cube', 'lightblue')

    # Rotated cube
    ax2 = fig.add_subplot(2, 4, 2, projection='3d')
    draw_cube(ax2, rotated_cube, 'Rotation Transformation', 'lightgreen')

    # Scaled cube
    ax3 = fig.add_subplot(2, 4, 3, projection='3d')
    draw_cube(ax3, scaled_cube, 'Scaling Transformation', 'lightcoral')

    # Combined transformation cube
    ax4 = fig.add_subplot(2, 4, 4, projection='3d')
    draw_cube(ax4, transformed_cube, 'Combined Transformation', 'lightyellow')

    # Projection transformation demonstration
    ax5 = fig.add_subplot(2, 4, 5)

    # Orthographic projection onto XY plane
    projected_xy = cube_vertices[:, :2]

    # Draw projected figure
    for face in cube_faces:
        face_2d = projected_xy[face]
        # Close polygon
        face_2d_closed = np.vstack([face_2d, face_2d[0]])
        ax5.plot(face_2d_closed[:, 0], face_2d_closed[:, 1], 'b-', alpha=0.7)

    ax5.scatter(projected_xy[:, 0], projected_xy[:, 1], c='red', s=50)
    ax5.set_xlabel('X')
    ax5.set_ylabel('Y')
    ax5.set_title('XY Plane Orthographic Projection')
    ax5.grid(True, alpha=0.3)
    ax5.set_aspect('equal')

    # Perspective projection
    ax6 = fig.add_subplot(2, 4, 6)

    # Simple perspective projection: z coordinate as depth information
    d = 3  # Viewing distance
    perspective_x = cube_vertices[:, 0] / (1 + cube_vertices[:, 2]/d)
    perspective_y = cube_vertices[:, 1] / (1 + cube_vertices[:, 2]/d)

    for face in cube_faces:
        face_x = perspective_x[face]
        face_y = perspective_y[face]
        # Close polygon
        face_x = np.append(face_x, face_x[0])
        face_y = np.append(face_y, face_y[0])
        ax6.plot(face_x, face_y, 'g-', alpha=0.7)

    ax6.scatter(perspective_x, perspective_y, c='red', s=50)
    ax6.set_xlabel('X')
    ax6.set_ylabel('Y')
    ax6.set_title('Perspective Projection')
    ax6.grid(True, alpha=0.3)
    ax6.set_aspect('equal')

    # Transformation matrix visualization
    ax7 = fig.add_subplot(2, 4, 7)

    matrices = [R_x, R_y, R_z, S]
    labels = ['Rot_X', 'Rot_Y', 'Rot_Z', 'Scale']

    for i, (matrix, label) in enumerate(zip(matrices, labels)):
        # Display matrix determinant or condition number
        det_val = np.linalg.det(matrix)
        ax7.bar(i, det_val, alpha=0.7, label=f'{label}\ndet={det_val:.3f}')

    ax7.set_xlabel('Transformation Type')
    ax7.set_ylabel('Determinant')
    ax7.set_title('Transformation Matrix Determinant')
    ax7.set_xticks(range(4))
    ax7.set_xticklabels(labels)
    ax7.legend()
    ax7.grid(True, alpha=0.3)

    # Transformation sequence
    ax8 = fig.add_subplot(2, 4, 8)

    # Show effects of transformation steps
    steps = ['Original', 'Rotate X', 'Rotate Y', 'Rotate Z', 'Scale']
    step_vertices = [
        cube_vertices,
        (R_x @ cube_vertices.T).T,
        (R_y @ R_x @ cube_vertices.T).T,
        (R_z @ R_y @ R_x @ cube_vertices.T).T,
        (S @ R_z @ R_y @ R_x @ cube_vertices.T).T
    ]

    # Compute center displacement after each step
    centers = [np.mean(vertices, axis=0) for vertices in step_vertices]
    center_distances = [np.linalg.norm(center - centers[0]) for center in centers]

    ax8.plot(range(len(steps)), center_distances, 'o-', linewidth=2, markersize=8)
    ax8.set_xlabel('Transformation Step')
    ax8.set_ylabel('Center Displacement Distance')
    ax8.set_title('Center Displacement in Transformation Sequence')
    ax8.set_xticks(range(len(steps)))
    ax8.set_xticklabels(steps, rotation=45)
    ax8.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()

    # Animation effect demonstration (static display of multiple frames)
    print(f"\nCreating Rotation Animation Frames:")

    fig_anim, axes_anim = plt.subplots(2, 4, figsize=(16, 8),
                                      subplot_kw={'projection': '3d'})
    axes_anim = axes_anim.flatten()

    angles = np.linspace(0, 2*np.pi, 8)

    for i, angle in enumerate(angles):
        R_anim = rotation_y(angle)
        rotated_vertices = (R_anim @ cube_vertices.T).T

        draw_cube(axes_anim[i], rotated_vertices,
                 f'Rotation {angle*180/np.pi:.0f}°', 'lightblue')

    plt.tight_layout()
    plt.show()

    return cube_vertices, transform_combined

# Execute computer graphics demonstration
cube_vertices, transform_matrix = computer_graphics_demo()

Chapter Summary

Summary of Application Areas

Application Area	Core Concept	Key Techniques	Practical Value
Data Science	Eigenvalue Decomposition	PCA, SVD	Dimensionality Reduction, Denoising
Graph Theory	Adjacency Matrix	Spectral Graph Theory	Network Analysis
Markov Chains	Transition Matrix	Steady-State Analysis	Predictive Modeling
Computer Graphics	Transformation Matrices	3D Transformations	Rendering, Animation
Machine Learning	Optimization Theory	Gradient Methods	Model Training

Mathematical Toolbox

def linear_algebra_toolbox_summary():
    """Linear algebra toolbox summary"""

    print("Practical Linear Algebra Toolbox")
    print("=" * 50)

    tools = {
        "Matrix Decompositions": {
            "Eigenvalue Decomposition": "A = QΛQ⁻¹, used for PCA, steady-state analysis",
            "Singular Value Decomposition": "A = UΣVᵀ, used for dimensionality reduction, compression",
            "QR Decomposition": "A = QR, used for least squares, orthogonalization",
            "Cholesky Decomposition": "A = LLᵀ, used for solving positive definite matrices"
        },

        "Numerical Methods": {
            "Power Method": "Find dominant eigenvalue, used for PageRank",
            "Jacobi Method": "Find all eigenvalues, used for symmetric matrices",
            "Gram-Schmidt": "Orthogonalization, used to construct orthonormal basis",
            "Least Squares": "Overdetermined system solving, used for data fitting"
        },

        "Application Algorithms": {
            "Principal Component Analysis": "Data dimensionality reduction and visualization",
            "Markov Chains": "State transition and prediction",
            "Graph Algorithms": "Network analysis and ranking",
            "Graphics Transformations": "3D modeling and rendering"
        }
    }

    for category, items in tools.items():
        print(f"\n{category}:")
        print("-" * 30)
        for tool, description in items.items():
            print(f"• {tool:15} : {description}")

    # Decision tree for choosing appropriate tools
    print(f"\nTool Selection Guide:")
    print("=" * 30)

    decision_tree = """
    Data Analysis Problem?
    ├─ Dimensionality Reduction → PCA (Eigenvalue Decomposition)
    ├─ Data Compression → SVD
    └─ Regression Analysis → Least Squares

    Network Problem?
    ├─ Node Importance → PageRank (Dominant Eigenvector)
    ├─ Community Detection → Spectral Clustering (Laplacian Matrix)
    └─ Connectivity Analysis → Algebraic Connectivity

    Dynamic System?
    ├─ Steady-State Prediction → Markov Chain (Transition Matrix)
    ├─ Vibration Analysis → Eigenfrequencies (Eigenvalues)
    └─ Stability → Sign of Real Part of Eigenvalues

    Graphics Problem?
    ├─ 3D Transformation → Homogeneous Coordinate Matrices
    ├─ Projection → Projection Matrices
    └─ Animation → Interpolation Transformations
    """

    print(decision_tree)

linear_algebra_toolbox_summary()

Cross-Disciplinary Impact

Linear algebra, as a fundamental branch of mathematics, plays an increasingly important role in modern science and technology:

1. Artificial Intelligence: Deep learning is fundamentally matrix operations
2. Quantum Computing: Quantum states are represented by vectors, quantum gates by matrices
3. Economics: Input-output models, game theory
4. Bioinformatics: Gene sequence analysis, protein structure prediction
5. Signal Processing: Fourier transforms, filter design

Linear algebra is not just a mathematical tool, but a way of thinking for understanding and solving complex problems. Through the abstraction of matrices and vectors, we can handle high-dimensional data, complex networks, and dynamic systems, making linear algebra an indispensable foundation for modern science and technology.