Chapter 11: Matrix Diagonalization Theory

Basic Concepts of Matrix Diagonalization

Definition

For an $n \times n$ matrix $A$, if there exists an invertible matrix $P$ and diagonal matrix $D$ such that: $P^{-1}AP = D$

Then matrix $A$ is called diagonalizable, and $P$ is called the diagonalizing matrix.

Equivalent Conditions for Diagonalization

Matrix $A$ is diagonalizable if and only if:

1. $A$ has $n$ linearly independent eigenvectors
2. For each eigenvalue, its geometric multiplicity equals its algebraic multiplicity
3. $A$ is similar to some diagonal matrix

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eig, inv
import matplotlib.patches as patches

def diagonalization_demo():
    """Demonstrate matrix diagonalization process"""

    # Diagonalizable matrix example
    A = np.array([[4, 1],
                  [2, 3]], dtype=float)

    print("Original matrix A:")
    print(A)

    # Calculate eigenvalues and eigenvectors
    eigenvals, eigenvecs = eig(A)

    print(f"\nEigenvalues: {eigenvals}")
    print(f"Eigenvector matrix P:\n{eigenvecs}")

    # Construct diagonal matrix
    D = np.diag(eigenvals)
    P = eigenvecs
    P_inv = inv(P)

    print(f"\nDiagonal matrix D:\n{D}")

    # Verify diagonalization
    result = P_inv @ A @ P
    print(f"\nP^(-1)AP =\n{result}")
    print(f"Equal to D: {np.allclose(result, D)}")

    # Verify inverse transformation
    reconstructed_A = P @ D @ P_inv
    print(f"\nPDP^(-1) =\n{reconstructed_A}")
    print(f"Equal to original matrix A: {np.allclose(reconstructed_A, A)}")

    return A, P, D

A, P, D = diagonalization_demo()

Properties of Similar Matrices

Definition of Similar Matrices

If there exists an invertible matrix $P$ such that $B = P^{-1}AP$, then matrices $A$ and $B$ are called similar, denoted as $A \sim B$.

Important Properties of Similar Matrices

Similar matrices have the same:

1. Eigenvalues
2. Determinant
3. Trace (sum of diagonal elements)
4. Rank
5. Characteristic polynomial

def similarity_properties_demo():
    """Demonstrate properties of similar matrices"""

    # Original matrix
    A = np.array([[2, 1, 0],
                  [1, 2, 1],
                  [0, 1, 2]], dtype=float)

    # Random invertible transformation matrix
    np.random.seed(42)
    P = np.random.randn(3, 3)
    while np.abs(np.linalg.det(P)) < 0.1:  # Ensure P is invertible
        P = np.random.randn(3, 3)

    P_inv = inv(P)

    # Calculate similar matrix
    B = P_inv @ A @ P

    print("Original matrix A:")
    print(A)
    print(f"\nTransformation matrix P:\n{P}")
    print(f"\nSimilar matrix B = P^(-1)AP:\n{B}")

    # Verify properties of similar matrices
    properties = {
        "Eigenvalues": (np.sort(eig(A)[0]), np.sort(eig(B)[0])),
        "Determinant": (np.linalg.det(A), np.linalg.det(B)),
        "Trace": (np.trace(A), np.trace(B)),
        "Rank": (np.linalg.matrix_rank(A), np.linalg.matrix_rank(B))
    }

    print(f"\nVerifying properties of similar matrices:")
    print("=" * 50)
    for prop_name, (val_A, val_B) in properties.items():
        if prop_name == "Eigenvalues":
            equal = np.allclose(val_A, val_B)
            print(f"{prop_name}:")
            print(f"  A: {val_A}")
            print(f"  B: {val_B}")
        else:
            equal = np.isclose(val_A, val_B)
            print(f"{prop_name}: A={val_A:.4f}, B={val_B:.4f}")
        print(f"  Equal: {equal}")
        print("-" * 30)

similarity_properties_demo()

Diagonalization Algorithm

Standard Diagonalization Steps

1. Calculate eigenvalues: solve $\det(A - \lambda I) = 0$
2. For each eigenvalue, find corresponding eigenvectors
3. Check if there are enough linearly independent eigenvectors
4. Construct matrix $P = [v_1, v_2, \ldots, v_n]$
5. Verify $P^{-1}AP = D$

def diagonalization_algorithm():
    """Complete matrix diagonalization algorithm"""

    def diagonalize_matrix(A):
        """
        Matrix diagonalization algorithm
        Returns: (is_diagonalizable, P, D)
        """
        n = A.shape[0]

        # Step 1: Calculate eigenvalues and eigenvectors
        eigenvals, eigenvecs = eig(A)

        # Step 2: Check linear independence of eigenvectors
        rank_eigenvecs = np.linalg.matrix_rank(eigenvecs)

        if rank_eigenvecs == n:
            # Diagonalizable
            P = eigenvecs
            D = np.diag(eigenvals)
            return True, P, D
        else:
            # Not diagonalizable
            return False, None, None

    # Test matrix 1: Diagonalizable matrix
    print("Test matrix 1: Diagonalizable matrix")
    A1 = np.array([[3, 1],
                   [0, 2]], dtype=float)

    print(f"Matrix A1:\n{A1}")

    is_diag1, P1, D1 = diagonalize_matrix(A1)

    if is_diag1:
        print("Matrix is diagonalizable!")
        print(f"Diagonalizing matrix P:\n{P1}")
        print(f"Diagonal matrix D:\n{D1}")

        # Verification
        verification = inv(P1) @ A1 @ P1
        print(f"Verification P^(-1)AP:\n{verification}")
        print(f"Error: {np.linalg.norm(verification - D1):.2e}")
    else:
        print("Matrix is not diagonalizable")

    print("\n" + "="*50 + "\n")

    # Test matrix 2: Non-diagonalizable matrix (repeated eigenvalue with insufficient geometric multiplicity)
    print("Test matrix 2: Non-diagonalizable matrix")
    A2 = np.array([[2, 1],
                   [0, 2]], dtype=float)

    print(f"Matrix A2:\n{A2}")

    is_diag2, P2, D2 = diagonalize_matrix(A2)

    if is_diag2:
        print("Matrix is diagonalizable!")
        print(f"Diagonalizing matrix P:\n{P2}")
        print(f"Diagonal matrix D:\n{D2}")
    else:
        print("Matrix is not diagonalizable")

        # Analyze reason
        eigenvals2, eigenvecs2 = eig(A2)
        print(f"Eigenvalues: {eigenvals2}")
        print(f"Rank of eigenvector matrix: {np.linalg.matrix_rank(eigenvecs2)}")
        print("Reason: Geometric multiplicity of repeated eigenvalue less than algebraic multiplicity")

    print("\n" + "="*50 + "\n")

    # Test matrix 3: Symmetric matrix (always diagonalizable)
    print("Test matrix 3: Symmetric matrix")
    A3 = np.array([[1, 2, 3],
                   [2, 4, 5],
                   [3, 5, 6]], dtype=float)

    print(f"Matrix A3:\n{A3}")
    print(f"Is symmetric: {np.allclose(A3, A3.T)}")

    is_diag3, P3, D3 = diagonalize_matrix(A3)

    if is_diag3:
        print("Symmetric matrix is diagonalizable!")
        print(f"Eigenvalues: {np.diag(D3)}")

        # Verify orthogonality of eigenvectors (for symmetric matrices)
        dot_products = []
        for i in range(3):
            for j in range(i+1, 3):
                dot_prod = np.dot(P3[:, i], P3[:, j])
                dot_products.append(dot_prod)
                print(f"v{i+1} · v{j+1} = {dot_prod:.2e}")

        print(f"Eigenvectors are orthogonal: {np.allclose(dot_products, 0, atol=1e-10)}")

diagonalization_algorithm()

Geometric Meaning of Diagonalization

Coordinate Transformation Perspective

Diagonalization essentially finds an eigenvector basis, where the linear transformation appears as simple coordinate scaling.

def geometric_interpretation():
    """Geometric interpretation of diagonalization"""

    # Define diagonalizable matrix
    A = np.array([[3, 1],
                  [1, 2]], dtype=float)

    # Diagonalize
    eigenvals, eigenvecs = eig(A)
    P = eigenvecs
    D = np.diag(eigenvals)

    print("Original matrix A:")
    print(A)
    print(f"Eigenvalues: {eigenvals}")

    # Create test vectors
    test_vectors = np.array([[1, 0, 1, 1],
                            [0, 1, 1, -1]])

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

    # First row: Transformation in standard basis
    ax1, ax2, ax3 = axes[0]

    # Original vectors (standard basis)
    for i in range(test_vectors.shape[1]):
        ax1.arrow(0, 0, test_vectors[0, i], test_vectors[1, i],
                 head_width=0.1, head_length=0.1, fc='blue', ec='blue',
                 alpha=0.7, label=f'v{i+1}' if i < 2 else None)

    # Draw standard basis
    ax1.arrow(0, 0, 1, 0, head_width=0.05, head_length=0.05,
             fc='red', ec='red', linewidth=2, label='e1')
    ax1.arrow(0, 0, 0, 1, head_width=0.05, head_length=0.05,
             fc='green', ec='green', linewidth=2, label='e2')

    ax1.set_xlim(-1.5, 2)
    ax1.set_ylim(-1.5, 2)
    ax1.grid(True, alpha=0.3)
    ax1.set_aspect('equal')
    ax1.legend()
    ax1.set_title('Original Vectors (Standard Basis)')

    # Transformed vectors (standard basis)
    transformed_vectors = A @ test_vectors

    for i in range(test_vectors.shape[1]):
        ax2.arrow(0, 0, test_vectors[0, i], test_vectors[1, i],
                 head_width=0.1, head_length=0.1, fc='blue', ec='blue',
                 alpha=0.3, linestyle='--')
        ax2.arrow(0, 0, transformed_vectors[0, i], transformed_vectors[1, i],
                 head_width=0.1, head_length=0.1, fc='red', ec='red',
                 alpha=0.7)

    # Draw eigenvectors
    for i, (val, vec) in enumerate(zip(eigenvals, eigenvecs.T)):
        length = 1.5
        ax2.arrow(0, 0, vec[0]*length, vec[1]*length,
                 head_width=0.1, head_length=0.1, fc=f'C{i+2}', ec=f'C{i+2}',
                 linewidth=3, label=f'Eigenvector{i+1}')

    ax2.set_xlim(-1, 4)
    ax2.set_ylim(-1, 3)
    ax2.grid(True, alpha=0.3)
    ax2.set_aspect('equal')
    ax2.legend()
    ax2.set_title('Transformed (Standard Basis)')

    # Representation in eigenvector basis
    P_inv = inv(P)
    coords_in_eigenbasis = P_inv @ test_vectors

    for i in range(test_vectors.shape[1]):
        ax3.arrow(0, 0, coords_in_eigenbasis[0, i], coords_in_eigenbasis[1, i],
                 head_width=0.1, head_length=0.1, fc='blue', ec='blue',
                 alpha=0.7)

    # Draw eigenvector basis (unit vectors in eigenvector coordinate system)
    ax3.arrow(0, 0, 1, 0, head_width=0.05, head_length=0.05,
             fc='purple', ec='purple', linewidth=2, label='Eigenbasis1')
    ax3.arrow(0, 0, 0, 1, head_width=0.05, head_length=0.05,
             fc='orange', ec='orange', linewidth=2, label='Eigenbasis2')

    ax3.set_xlim(-2, 2)
    ax3.set_ylim(-2, 2)
    ax3.grid(True, alpha=0.3)
    ax3.set_aspect('equal')
    ax3.legend()
    ax3.set_title('Eigenvector Coordinate System')

    # Second row: Transformation in eigenvector basis
    ax4, ax5, ax6 = axes[1]

    # Transformation in eigenvector basis (diagonal matrix action)
    transformed_coords = D @ coords_in_eigenbasis

    for i in range(test_vectors.shape[1]):
        ax4.arrow(0, 0, coords_in_eigenbasis[0, i], coords_in_eigenbasis[1, i],
                 head_width=0.1, head_length=0.1, fc='blue', ec='blue',
                 alpha=0.3, linestyle='--')
        ax4.arrow(0, 0, transformed_coords[0, i], transformed_coords[1, i],
                 head_width=0.1, head_length=0.1, fc='red', ec='red',
                 alpha=0.7)

    ax4.arrow(0, 0, 1, 0, head_width=0.05, head_length=0.05,
             fc='purple', ec='purple', linewidth=2)
    ax4.arrow(0, 0, 0, 1, head_width=0.05, head_length=0.05,
             fc='orange', ec='orange', linewidth=2)

    ax4.set_xlim(-2, 4)
    ax4.set_ylim(-2, 3)
    ax4.grid(True, alpha=0.3)
    ax4.set_aspect('equal')
    ax4.set_title(f'Transformation in Eigenbasis\\nD = diag({eigenvals[0]:.2f}, {eigenvals[1]:.2f})')

    # Display diagonalization decomposition
    ax5.text(0.1, 0.8, r'$A = PDP^{-1}$', fontsize=16, transform=ax5.transAxes)
    ax5.text(0.1, 0.6, r'$P = $ Eigenvector matrix', fontsize=12, transform=ax5.transAxes)
    ax5.text(0.1, 0.4, r'$D = $ Diagonal matrix', fontsize=12, transform=ax5.transAxes)
    ax5.text(0.1, 0.2, r'$P^{-1} = $ Coordinate transformation', fontsize=12, transform=ax5.transAxes)

    ax5.text(0.1, 0.1, f'P = \n{P}', fontsize=10, transform=ax5.transAxes,
             fontfamily='monospace')

    ax5.set_xlim(0, 1)
    ax5.set_ylim(0, 1)
    ax5.axis('off')
    ax5.set_title('Diagonalization Decomposition')

    # Transformation process diagram
    ax6.text(0.5, 0.9, 'Three Steps of Diagonalization:', fontsize=14, ha='center',
             transform=ax6.transAxes, weight='bold')

    steps = [
        '1. Standard basis → Eigenbasis',
        '   (multiply by P⁻¹)',
        '2. Scaling in eigenbasis',
        '   (multiply by D)',
        '3. Eigenbasis → Standard basis',
        '   (multiply by P)'
    ]

    for i, step in enumerate(steps):
        y_pos = 0.75 - i * 0.1
        color = 'blue' if i % 2 == 0 else 'gray'
        weight = 'bold' if i % 2 == 0 else 'normal'
        ax6.text(0.1, y_pos, step, fontsize=11, transform=ax6.transAxes,
                color=color, weight=weight)

    ax6.set_xlim(0, 1)
    ax6.set_ylim(0, 1)
    ax6.axis('off')
    ax6.set_title('Transformation Steps')

    plt.tight_layout()
    plt.show()

geometric_interpretation()

Application of Diagonalization: Matrix Powers

Fast Computation of Matrix Powers

If $A = PDP^{-1}$, then: $A^n = PD^nP^{-1}$

This greatly simplifies the computation of matrix powers.

def matrix_power_application():
    """Application of matrix diagonalization in computing matrix powers"""

    # Define matrix
    A = np.array([[2, 1],
                  [1, 2]], dtype=float)

    print("Original matrix A:")
    print(A)

    # Diagonalize
    eigenvals, eigenvecs = eig(A)
    P = eigenvecs
    P_inv = inv(P)
    D = np.diag(eigenvals)

    print(f"\nDiagonalization: A = PDP^(-1)")
    print(f"P = \n{P}")
    print(f"D = \n{D}")

    # Calculate different powers
    powers = [2, 5, 10, 20]

    print(f"\nMatrix power calculation comparison:")
    print("="*60)

    for n in powers:
        # Method 1: Direct calculation (inefficient)
        A_power_direct = np.linalg.matrix_power(A, n)

        # Method 2: Using diagonalization (efficient)
        D_power = np.diag(eigenvals**n)
        A_power_diag = P @ D_power @ P_inv

        print(f"\nA^{n}:")
        print("Direct calculation:")
        print(A_power_direct)
        print("Diagonalization calculation:")
        print(A_power_diag)
        print(f"Error: {np.linalg.norm(A_power_direct - A_power_diag):.2e}")

    # Visualize eigenvalue power growth
    n_values = np.arange(1, 11)
    eigenval_powers = np.array([eigenvals**n for n in n_values])

    plt.figure(figsize=(12, 5))

    plt.subplot(1, 2, 1)
    for i, eigenval in enumerate(eigenvals):
        plt.plot(n_values, eigenval_powers[:, i], 'o-',
                label=f'λ{i+1} = {eigenval:.2f}', linewidth=2)

    plt.xlabel('n')
    plt.ylabel('λⁿ')
    plt.title('Powers of Eigenvalues')
    plt.legend()
    plt.grid(True, alpha=0.3)

    plt.subplot(1, 2, 2)
    for i, eigenval in enumerate(eigenvals):
        plt.semilogy(n_values, eigenval_powers[:, i], 'o-',
                    label=f'λ{i+1} = {eigenval:.2f}', linewidth=2)

    plt.xlabel('n')
    plt.ylabel('λⁿ (log scale)')
    plt.title('Powers of Eigenvalues (Log Scale)')
    plt.legend()
    plt.grid(True, alpha=0.3)

    plt.tight_layout()
    plt.show()

    # Fibonacci sequence application
    print(f"\nApplication case: Fibonacci sequence")
    print("Fibonacci recurrence can be written in matrix form:")
    print("[ F(n+1) ]   [ 1  1 ] [ F(n)   ]")
    print("[ F(n)   ] = [ 1  0 ] [ F(n-1) ]")

    fib_matrix = np.array([[1, 1],
                          [1, 0]], dtype=float)

    # Diagonalize Fibonacci matrix
    fib_eigenvals, fib_eigenvecs = eig(fib_matrix)
    fib_P = fib_eigenvecs
    fib_P_inv = inv(fib_P)

    print(f"\nEigenvalues of Fibonacci matrix: {fib_eigenvals}")
    print("These are the golden ratio φ and -1/φ!")

    phi = (1 + np.sqrt(5)) / 2
    print(f"Golden ratio φ = {phi:.6f}")
    print(f"Eigenvalue1 = {fib_eigenvals[0]:.6f}")
    print(f"Eigenvalue2 = {fib_eigenvals[1]:.6f}")

    # Calculate nth Fibonacci number
    def fibonacci_formula(n):
        """Calculate Fibonacci number using eigenvalue formula"""
        phi = (1 + np.sqrt(5)) / 2
        psi = (1 - np.sqrt(5)) / 2
        return int((phi**n - psi**n) / np.sqrt(5))

    print(f"\nFirst 15 Fibonacci numbers:")
    for i in range(1, 16):
        fib_n = fibonacci_formula(i)
        print(f"F({i}) = {fib_n}")

matrix_power_application()

Introduction to Jordan Canonical Form

Canonical Form for Non-Diagonalizable Matrices

Not all matrices are diagonalizable, but every matrix is similar to a Jordan canonical form matrix.

Structure of Jordan Blocks

A Jordan block has the following form:

\lambda & 1 & 0 & \cdots & 0 \\ 0 & \lambda & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & \lambda & 1 \\ 0 & 0 & 0 & 0 & \lambda \end{pmatrix}$$ ```python def jordan_form_introduction(): """Introduction to Jordan canonical form""" # Non-diagonalizable matrix example A = np.array([[2, 1, 0], [0, 2, 1], [0, 0, 2]], dtype=float) print("Non-diagonalizable matrix A:") print(A) # Calculate eigenvalues eigenvals, eigenvecs = eig(A) print(f"\nEigenvalues: {eigenvals}") print(f"Eigenvalue multiplicity: 3 (algebraic multiplicity)") # Calculate dimension of eigenspace lambda_val = eigenvals[0] # All eigenvalues are 2 null_space_matrix = A - lambda_val * np.eye(3) rank_null = 3 - np.linalg.matrix_rank(null_space_matrix) print(f"Eigenspace dimension: {rank_null} (geometric multiplicity)") print(f"\nBecause geometric multiplicity < algebraic multiplicity, matrix is not diagonalizable") # Construct Jordan canonical form print(f"\nJordan canonical form (3×3 Jordan block):") J = np.array([[2, 1, 0], [0, 2, 1], [0, 0, 2]]) print(J) print(f"\nMatrix A is already in Jordan canonical form!") # Demonstrate properties of Jordan blocks print(f"\nPowers of Jordan block:") for n in [2, 3, 4, 5]: J_power = np.linalg.matrix_power(J, n) print(f"\nJ^{n} =") print(J_power) # Visualize Jordan block structure fig, axes = plt.subplots(1, 4, figsize=(16, 4)) jordan_blocks = [ (np.array([[2]]), "1×1 Jordan block"), (np.array([[2, 1], [0, 2]]), "2×2 Jordan block"), (np.array([[2, 1, 0], [0, 2, 1], [0, 0, 2]]), "3×3 Jordan block"), (np.array([[2, 1, 0, 0], [0, 2, 1, 0], [0, 0, 2, 1], [0, 0, 0, 2]]), "4×4 Jordan block") ] for i, (block, title) in enumerate(jordan_blocks): ax = axes[i] im = ax.imshow(block, cmap='RdYlBu', aspect='equal') # Add numerical labels for row in range(block.shape[0]): for col in range(block.shape[1]): value = block[row, col] if value != 0: ax.text(col, row, f'{value:.0f}', ha='center', va='center', fontsize=12, weight='bold') ax.set_title(title) ax.set_xticks([]) ax.set_yticks([]) # Add borders for row in range(block.shape[0]): for col in range(block.shape[1]): rect = patches.Rectangle((col-0.5, row-0.5), 1, 1, linewidth=1, edgecolor='black', facecolor='none') ax.add_patch(rect) plt.tight_layout() plt.show() # General structure of Jordan canonical form print(f"\nGeneral structure of Jordan canonical form:") print("J = diag(J₁, J₂, ..., Jₖ)") print("where each Jᵢ is a Jordan block") print("\nImportant properties:") print("1. Each Jordan block corresponds to an eigenvalue") print("2. Size of Jordan block determines geometric multiplicity of eigenvalue") print("3. Every matrix is similar to a unique Jordan canonical form") jordan_form_introduction() ``` ## Diagonalization of Symmetric Matrices ### Spectral Theorem Symmetric matrices can always be **orthogonally diagonalized**, i.e., there exists an orthogonal matrix $Q$ such that: $$Q^TAQ = D$$ ```python def symmetric_matrix_diagonalization(): """Orthogonal diagonalization of symmetric matrices""" # Construct symmetric matrix A = np.array([[3, 1, 2], [1, 4, 0], [2, 0, 5]], dtype=float) print("Symmetric matrix A:") print(A) print(f"Is symmetric: {np.allclose(A, A.T)}") # Calculate eigenvalues and eigenvectors eigenvals, eigenvecs = eig(A) # Sort by eigenvalue magnitude idx = np.argsort(eigenvals)[::-1] eigenvals = eigenvals[idx] eigenvecs = eigenvecs[:, idx] print(f"\nEigenvalues: {eigenvals}") print("All eigenvalues are real:", np.all(np.isreal(eigenvals))) # Check orthogonality of eigenvectors print(f"\nChecking orthogonality of eigenvectors:") for i in range(3): for j in range(i+1, 3): dot_product = np.dot(eigenvecs[:, i], eigenvecs[:, j]) print(f"v{i+1} · v{j+1} = {dot_product:.2e}") # Normalize eigenvectors to get orthogonal matrix Q = eigenvecs / np.linalg.norm(eigenvecs, axis=0) D = np.diag(eigenvals) print(f"\nOrthogonal matrix Q:") print(Q) # Verify orthogonality Q^T Q = I QTQ = Q.T @ Q print(f"\nQ^T Q =") print(QTQ) print(f"Is identity matrix: {np.allclose(QTQ, np.eye(3))}") # Verify diagonalization Q^T A Q = D QTAQ = Q.T @ A @ Q print(f"\nQ^T A Q =") print(QTAQ) print(f"Diagonal matrix D =") print(D) print(f"Diagonalization successful: {np.allclose(QTAQ, D)}") # Visualize quadratic form of symmetric matrix fig = plt.figure(figsize=(15, 5)) # 1. Contour lines in original coordinate system ax1 = fig.add_subplot(131) x = np.linspace(-2, 2, 100) y = np.linspace(-2, 2, 100) X, Y = np.meshgrid(x, y) # Quadratic form x^T A x Z = A[0,0]*X**2 + A[1,1]*Y**2 + 2*A[0,1]*X*Y contours1 = ax1.contour(X, Y, Z, levels=[1, 4, 9, 16], colors='blue') ax1.clabel(contours1, inline=True, fontsize=8) # Draw eigenvectors for i, vec in enumerate(eigenvecs.T): ax1.arrow(0, 0, vec[0], vec[1], head_width=0.1, head_length=0.1, fc=f'C{i}', ec=f'C{i}', linewidth=3, label=f'Eigenvector{i+1}') ax1.set_xlim(-2, 2) ax1.set_ylim(-2, 2) ax1.set_aspect('equal') ax1.grid(True, alpha=0.3) ax1.legend() ax1.set_title('Original Coordinate System') # 2. Principal axis coordinate system ax2 = fig.add_subplot(132) # In principal axis coordinate system, quadratic form becomes diagonal u = np.linspace(-2, 2, 100) v = np.linspace(-2, 2, 100) U, V = np.meshgrid(u, v) Z_diag = eigenvals[0]*U**2 + eigenvals[1]*V**2 contours2 = ax2.contour(U, V, Z_diag, levels=[1, 4, 9, 16], colors='red') ax2.clabel(contours2, inline=True, fontsize=8) # Coordinate axes ax2.arrow(0, 0, 1.5, 0, head_width=0.1, head_length=0.1, fc='red', ec='red', linewidth=3, label='Principal axis 1') ax2.arrow(0, 0, 0, 1.5, head_width=0.1, head_length=0.1, fc='green', ec='green', linewidth=3, label='Principal axis 2') ax2.set_xlim(-2, 2) ax2.set_ylim(-2, 2) ax2.set_aspect('equal') ax2.grid(True, alpha=0.3) ax2.legend() ax2.set_title('Principal Axis Coordinate System') # 3. Eigenvalue spectrum ax3 = fig.add_subplot(133) bars = ax3.bar(range(1, 4), eigenvals, color=['red', 'green', 'blue'], alpha=0.7) ax3.set_xlabel('Eigenvalue index') ax3.set_ylabel('Eigenvalue magnitude') ax3.set_title('Eigenvalue Spectrum') ax3.grid(True, alpha=0.3) # Add numerical labels for i, (bar, val) in enumerate(zip(bars, eigenvals)): ax3.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.1, f'{val:.2f}', ha='center', va='bottom') plt.tight_layout() plt.show() symmetric_matrix_diagonalization() ``` ## Application: Diagonalization of Quadratic Forms ### Standardization of Quadratic Forms Through orthogonal transformation, a quadratic form can be reduced to standard form: $$f(x_1, x_2, \ldots, x_n) = \lambda_1 y_1^2 + \lambda_2 y_2^2 + \cdots + \lambda_n y_n^2$$ ```python def quadratic_form_diagonalization(): """Application of quadratic form diagonalization""" # Define quadratic form matrix A = np.array([[2, 1, 0], [1, 2, 1], [0, 1, 2]], dtype=float) print("Quadratic form matrix A:") print(A) # Corresponding quadratic form print("\nCorresponding quadratic form:") print("f(x₁,x₂,x₃) = 2x₁² + 2x₂² + 2x₃² + 2x₁x₂ + 2x₂x₃") # Orthogonal diagonalization eigenvals, eigenvecs = eig(A) # Sort idx = np.argsort(eigenvals)[::-1] eigenvals = eigenvals[idx] eigenvecs = eigenvecs[:, idx] # Orthogonalize (already orthogonal, just normalize) Q = eigenvecs / np.linalg.norm(eigenvecs, axis=0) print(f"\nEigenvalues: {eigenvals}") print(f"Orthogonal transformation matrix Q:\n{Q}") # Standard form print(f"\nStandard form after orthogonal transformation x = Qy:") print(f"f = {eigenvals[0]:.3f}y₁² + {eigenvals[1]:.3f}y₂² + {eigenvals[2]:.3f}y₃²") # Determine sign of quadratic form positive_eigenvals = np.sum(eigenvals > 0) negative_eigenvals = np.sum(eigenvals < 0) zero_eigenvals = np.sum(np.abs(eigenvals) < 1e-10) print(f"\nSign characteristics of quadratic form:") print(f"Number of positive eigenvalues: {positive_eigenvals}") print(f"Number of negative eigenvalues: {negative_eigenvals}") print(f"Number of zero eigenvalues: {zero_eigenvals}") if negative_eigenvals == 0 and zero_eigenvals == 0: quadratic_type = "Positive definite" elif positive_eigenvals == 0 and zero_eigenvals == 0: quadratic_type = "Negative definite" elif negative_eigenvals == 0: quadratic_type = "Positive semidefinite" elif positive_eigenvals == 0: quadratic_type = "Negative semidefinite" else: quadratic_type = "Indefinite" print(f"Quadratic form type: {quadratic_type}") # Visualize different quadratic forms fig, axes = plt.subplots(2, 3, figsize=(15, 10)) quadratic_examples = [ (np.array([[1, 0], [0, 1]]), "Positive definite", [1, 4, 9]), (np.array([[-1, 0], [0, -1]]), "Negative definite", [1, 4, 9]), (np.array([[1, 0], [0, 0]]), "Positive semidefinite", [0, 1, 4]), (np.array([[1, 0], [0, -1]]), "Indefinite", [-4, -1, 1, 4]), (np.array([[2, 1], [1, 1]]), "Positive definite (non-diagonal)", [1, 4, 9]), (np.array([[1, 2], [2, -1]]), "Indefinite (non-diagonal)", [-4, -1, 1, 4]) ] for i, (matrix, title, levels) in enumerate(quadratic_examples): ax = axes[i//3, i%3] x = np.linspace(-3, 3, 100) y = np.linspace(-3, 3, 100) X, Y = np.meshgrid(x, y) Z = matrix[0,0]*X**2 + matrix[1,1]*Y**2 + 2*matrix[0,1]*X*Y if "Negative definite" in title: contours = ax.contour(X, Y, Z, levels=[-l for l in levels], colors='red') else: contours = ax.contour(X, Y, Z, levels=levels, colors='blue') ax.clabel(contours, inline=True, fontsize=8) # Draw eigenvectors evals, evecs = eig(matrix) for j, vec in enumerate(evecs.T): if np.abs(evals[j]) > 1e-10: # Non-zero eigenvalue length = 1.5 if evals[j] > 0 else 1.0 ax.arrow(0, 0, vec[0]*length, vec[1]*length, head_width=0.1, head_length=0.1, fc=f'C{j}', ec=f'C{j}', linewidth=2, alpha=0.7) ax.set_xlim(-3, 3) ax.set_ylim(-3, 3) ax.set_aspect('equal') ax.grid(True, alpha=0.3) ax.set_title(title) plt.tight_layout() plt.show() quadratic_form_diagonalization() ``` ## Chapter Summary ### Core Concept Summary | Concept | Definition | Condition | Application | |---------|------------|-----------|-------------| | Matrix diagonalization | $P^{-1}AP = D$ | n linearly independent eigenvectors | Matrix powers, differential equations | | Similar matrices | $B = P^{-1}AP$ | Invertible matrix P exists | Invariant preservation | | Jordan canonical form | Quasi-diagonal form | All matrices have one | Analysis of non-diagonalizable matrices | | Orthogonal diagonalization | $Q^TAQ = D$ | Symmetric matrix | Quadratic forms, principal component analysis | ### Diagonalization Decision Flow ```mermaid graph TD A[Given matrix A] --> B{Calculate eigenvalues} B --> C{Check repeated eigenvalues} C -->|No repeated eigenvalues| D[Diagonalizable] C -->|Has repeated eigenvalues| E{Geometric multiplicity = Algebraic multiplicity?} E -->|Yes| D E -->|No| F[Not diagonalizable] D --> G[Construct P and D] F --> H[Consider Jordan canonical form] I{Matrix symmetric?} --> J[Orthogonally diagonalizable] I -->|No| K[General diagonalization] ``` ### Important Application Domains ```python def applications_summary(): """Summary of important applications of diagonalization""" applications = { "Matrix power computation": { "Formula": "A^n = PD^nP^(-1)", "Advantage": "Complexity reduced from O(n³) to O(n)", "Application": "Fibonacci sequence, Markov chains" }, "Differential equation systems": { "Formula": "x' = Ax has solution x(t) = Pe^(Dt)P^(-1)x₀", "Advantage": "Decouples variables", "Application": "Dynamical systems, vibration analysis" }, "Quadratic form analysis": { "Formula": "x^TAx = y^TDy (via orthogonal transformation)", "Advantage": "Eliminates cross terms", "Application": "Optimization problems, ellipse analysis" }, "Principal component analysis": { "Formula": "Eigenvectors of covariance matrix are principal components", "Advantage": "Dimensionality reduction, decorrelation", "Application": "Data analysis, machine learning" }, "Graph theory applications": { "Formula": "Eigenvalues of adjacency matrix reflect graph properties", "Advantage": "Spectral graph theory", "Application": "Network analysis, clustering" } } print("Important applications of matrix diagonalization:") print("=" * 60) for app_name, details in applications.items(): print(f"\n{app_name}:") for key, value in details.items(): print(f" {key}: {value}") print("-" * 40) applications_summary() ``` Matrix diagonalization is one of the most important theoretical achievements in linear algebra. It not only provides a powerful tool for understanding the essence of linear transformations but also plays a critical role in scientific computing, engineering applications, and data analysis. Mastering diagonalization theory lays a solid foundation for subsequent study of advanced topics such as inner product spaces and quadratic forms, as well as for applying linear algebra methods in practical problems.