main.tex

\documentclass[11pt,a4paper]{article}

% Linear Algebra and Matrix Computations Template
% SEO Keywords: linear algebra latex template, matrix computations latex, eigenvalue problems latex

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage{mathtools}
\usepackage{siunitx}
\usepackage{graphicx}
\usepackage{float}
\usepackage{hyperref}
\usepackage{cleveref}
\usepackage{booktabs}
\usepackage{pythontex}

\usepackage[style=numeric,backend=biber]{biblatex}
\addbibresource{references.bib}

\newcommand{\norm}[1]{\left\|#1\right\|}
\newcommand{\cond}{\operatorname{cond}}
\newcommand{\rank}{\operatorname{rank}}
\newcommand{\trace}{\operatorname{tr}}

\title{Linear Algebra and Matrix Computations:\\
Advanced Algorithms and Numerical Methods}
\author{CoCalc Scientific Templates}
\date{\today}

\begin{document}

\maketitle

\begin{abstract}
This comprehensive linear algebra template demonstrates advanced matrix computations including eigenvalue algorithms, singular value decomposition, matrix factorizations, and iterative methods. Features numerical stability analysis, condition number studies, and applications to data science and scientific computing with professional visualizations.

\textbf{Keywords:} linear algebra, matrix computations, eigenvalue problems, SVD, matrix factorizations, numerical linear algebra
\end{abstract}

\tableofcontents
\newpage

\section{Matrix Decompositions and Factorizations}
\label{sec:matrix-decompositions}

\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('Agg')
from scipy.linalg import svd, qr, lu, cholesky, eigh, eigvals
from scipy.sparse import random as sparse_random
from scipy.sparse.linalg import spsolve, eigsh

np.random.seed(42)

def demonstrate_matrix_factorizations():
    """Demonstrate various matrix factorizations"""

    # Generate test matrices
    n = 6
    A_general = np.random.randn(n, n)
    A_spd = A_general @ A_general.T + 0.1 * np.eye(n)  # Symmetric positive definite
    A_rectangular = np.random.randn(8, n)

    print("Matrix Factorizations and Decompositions:")
    print(f"  Matrix size: {n}×{n}")

    # LU Decomposition
    P, L, U = lu(A_general, permute_l=False)
    lu_error = np.linalg.norm(P @ L @ U - A_general)
    print(f"  LU decomposition error: {lu_error:.2e}")

    # QR Decomposition
    Q, R = qr(A_rectangular)
    qr_error = np.linalg.norm(Q @ R - A_rectangular)
    print(f"  QR decomposition error: {qr_error:.2e}")

    # Cholesky Decomposition
    L_chol = cholesky(A_spd, lower=True)
    chol_error = np.linalg.norm(L_chol @ L_chol.T - A_spd)
    print(f"  Cholesky decomposition error: {chol_error:.2e}")

    # Singular Value Decomposition
    U_svd, s, Vt = svd(A_rectangular)
    # For rectangular matrix reconstruction, need to handle dimensions correctly
    m, n = A_rectangular.shape
    S_full = np.zeros((m, n))
    S_full[:n, :n] = np.diag(s)
    svd_error = np.linalg.norm(U_svd @ S_full @ Vt - A_rectangular)
    print(f"  SVD reconstruction error: {svd_error:.2e}")

    # Eigenvalue Decomposition
    eigenvals, eigenvecs = eigh(A_spd)
    eigen_error = np.linalg.norm(A_spd @ eigenvecs - eigenvecs @ np.diag(eigenvals))
    print(f"  Eigendecomposition error: {eigen_error:.2e}")

    return {
        'A_general': A_general, 'A_spd': A_spd, 'A_rectangular': A_rectangular,
        'LU': (P, L, U), 'QR': (Q, R), 'Cholesky': L_chol,
        'SVD': (U_svd, s, Vt), 'Eigen': (eigenvals, eigenvecs)
    }

def condition_number_analysis():
    """Analyze condition numbers and numerical stability"""

    # Create matrices with varying condition numbers
    sizes = [10, 20, 50, 100]
    condition_numbers = []
    spectral_norms = []

    for n in sizes:
        # Generate random matrix
        A = np.random.randn(n, n)

        # Make it symmetric positive definite
        A = A @ A.T

        # Add regularization to control condition number
        alpha = 1e-3
        A_reg = A + alpha * np.eye(n)

        cond_num = np.linalg.cond(A_reg)
        spec_norm = np.linalg.norm(A_reg, ord=2)

        condition_numbers.append(cond_num)
        spectral_norms.append(spec_norm)

    print(f"\nCondition Number Analysis:")
    for i, n in enumerate(sizes):
        print(f"  Size {n:3d}: cond = {condition_numbers[i]:.2e}, 2-norm = {spectral_norms[i]:.2f}")

    return sizes, condition_numbers, spectral_norms

# Perform matrix factorization demonstrations
decompositions = demonstrate_matrix_factorizations()
sizes, cond_nums, spec_norms = condition_number_analysis()
\end{pycode}

\section{Eigenvalue Problems and Spectral Analysis}
\label{sec:eigenvalue-problems}

\begin{pycode}
def eigenvalue_algorithms_comparison():
    """Compare different eigenvalue algorithms"""

    # Create test matrices with known eigenvalue properties
    n = 50

    # Symmetric tridiagonal matrix (fast algorithms available)
    diag_main = 2 * np.ones(n)
    diag_off = -1 * np.ones(n-1)
    A_tridiag = np.diag(diag_main) + np.diag(diag_off, 1) + np.diag(diag_off, -1)

    # Random symmetric matrix
    A_rand = np.random.randn(n, n)
    A_symmetric = (A_rand + A_rand.T) / 2

    # Sparse matrix
    density = 0.1
    A_sparse = sparse_random(n, n, density=density, format='csr')
    A_sparse = A_sparse + A_sparse.T  # Make symmetric

    print(f"\nEigenvalue Algorithm Comparison:")
    print(f"  Matrix size: {n}×{n}")

    # Full eigenvalue decomposition for dense matrices
    import time

    start_time = time.time()
    evals_tridiag, evecs_tridiag = eigh(A_tridiag)
    time_tridiag = time.time() - start_time

    start_time = time.time()
    evals_symmetric, evecs_symmetric = eigh(A_symmetric)
    time_symmetric = time.time() - start_time

    # Partial eigenvalue decomposition for sparse matrix
    k = 10  # Number of eigenvalues to compute
    start_time = time.time()
    evals_sparse, evecs_sparse = eigsh(A_sparse, k=k, which='LM')
    time_sparse = time.time() - start_time

    print(f"  Tridiagonal: {time_tridiag:.4f} s ({n} eigenvalues)")
    print(f"  Dense symmetric: {time_symmetric:.4f} s ({n} eigenvalues)")
    print(f"  Sparse: {time_sparse:.4f} s ({k} eigenvalues)")

    # Analyze eigenvalue distributions
    print(f"\nEigenvalue Statistics:")
    print(f"  Tridiagonal: lambda in [{evals_tridiag[0]:.3f}, {evals_tridiag[-1]:.3f}]")
    print(f"  Symmetric: lambda in [{evals_symmetric[0]:.3f}, {evals_symmetric[-1]:.3f}]")
    print(f"  Sparse: lambda in [{evals_sparse[0]:.3f}, {evals_sparse[-1]:.3f}]")

    return {
        'tridiagonal': (A_tridiag, evals_tridiag, evecs_tridiag),
        'symmetric': (A_symmetric, evals_symmetric, evecs_symmetric),
        'sparse': (A_sparse.toarray(), evals_sparse, evecs_sparse)
    }

def matrix_powers_and_functions():
    """Demonstrate matrix functions and powers"""

    n = 8
    A = np.random.randn(n, n)
    A = (A + A.T) / 2  # Make symmetric

    # Eigendecomposition
    eigenvals, eigenvecs = eigh(A)

    # Matrix exponential using eigendecomposition
    exp_eigenvals = np.exp(eigenvals)
    A_exp = eigenvecs @ np.diag(exp_eigenvals) @ eigenvecs.T

    # Matrix square root
    sqrt_eigenvals = np.sqrt(np.abs(eigenvals))  # Handle potential negative eigenvalues
    A_sqrt = eigenvecs @ np.diag(sqrt_eigenvals) @ eigenvecs.T

    # Verify: (A^(1/2))^2 ≈ A for positive definite case
    A_pd = A @ A.T + np.eye(n)  # Ensure positive definite
    eigenvals_pd, eigenvecs_pd = eigh(A_pd)
    sqrt_eigenvals_pd = np.sqrt(eigenvals_pd)
    A_sqrt_pd = eigenvecs_pd @ np.diag(sqrt_eigenvals_pd) @ eigenvecs_pd.T

    sqrt_error = np.linalg.norm(A_sqrt_pd @ A_sqrt_pd - A_pd)

    print(f"\nMatrix Functions:")
    print(f"  Matrix exponential computed via eigendecomposition")
    print(f"  Matrix square root verification error: {sqrt_error:.2e}")

    return A, A_exp, A_sqrt, eigenvals, eigenvecs

# Perform eigenvalue analysis
eigen_results = eigenvalue_algorithms_comparison()
A_test, A_exp, A_sqrt, evals_test, evecs_test = matrix_powers_and_functions()
\end{pycode}

\section{Iterative Methods for Large Systems}
\label{sec:iterative-methods}

\begin{pycode}
def iterative_solver_comparison():
    """Compare iterative methods for large linear systems"""

    from scipy.sparse.linalg import cg, gmres, bicgstab
    from scipy.sparse import diags

    # Create large sparse system
    n = 1000

    # 1D Laplacian (tridiagonal)
    diagonals = [-1, 2, -1]
    offsets = [-1, 0, 1]
    A_laplace = diags(diagonals, offsets, shape=(n, n), format='csr')

    # Right-hand side
    np.random.seed(42)
    b = np.random.randn(n)

    # True solution (for comparison)
    x_true = spsolve(A_laplace.tocsc(), b)

    print(f"\nIterative Solver Comparison:")
    print(f"  System size: {n}×{n}")
    print(f"  Matrix: 1D Laplacian (sparse)")

    # Conjugate Gradient
    x_cg, info_cg = cg(A_laplace, b, rtol=1e-8, maxiter=n)
    error_cg = np.linalg.norm(x_cg - x_true) / np.linalg.norm(x_true)

    # GMRES
    x_gmres, info_gmres = gmres(A_laplace, b, rtol=1e-8, maxiter=n, restart=50)
    error_gmres = np.linalg.norm(x_gmres - x_true) / np.linalg.norm(x_true)

    # BiCGSTAB
    x_bicgstab, info_bicgstab = bicgstab(A_laplace, b, rtol=1e-8, maxiter=n)
    error_bicgstab = np.linalg.norm(x_bicgstab - x_true) / np.linalg.norm(x_true)

    print(f"  CG: convergence = {info_cg}, error = {error_cg:.2e}")
    print(f"  GMRES: convergence = {info_gmres}, error = {error_gmres:.2e}")
    print(f"  BiCGSTAB: convergence = {info_bicgstab}, error = {error_bicgstab:.2e}")

    return A_laplace, b, x_true, x_cg, x_gmres, x_bicgstab

# Test iterative solvers
A_sparse, b_vec, x_true, x_cg, x_gmres, x_bicgstab = iterative_solver_comparison()
\end{pycode}

\begin{pycode}
# Create comprehensive linear algebra visualization
import os
os.makedirs('assets', exist_ok=True)

fig, axes = plt.subplots(3, 3, figsize=(18, 15))
fig.suptitle('Linear Algebra and Matrix Computations', fontsize=16, fontweight='bold')

# Matrix structure visualization
ax1 = axes[0, 0]
A_viz = decompositions['A_spd']
im1 = ax1.imshow(A_viz, cmap='RdBu_r', interpolation='nearest')
ax1.set_title('Symmetric Positive Definite Matrix')
plt.colorbar(im1, ax=ax1, shrink=0.8)

# Singular values
ax2 = axes[0, 1]
_, s_svd, _ = decompositions['SVD']
ax2.semilogy(range(1, len(s_svd)+1), s_svd, 'bo-', linewidth=2, markersize=6)
ax2.set_xlabel('Index')
ax2.set_ylabel('Singular Value')
ax2.set_title('Singular Value Spectrum')
ax2.grid(True, alpha=0.3)

# Eigenvalues
ax3 = axes[0, 2]
evals, _ = decompositions['Eigen']
ax3.plot(range(1, len(evals)+1), evals, 'ro-', linewidth=2, markersize=6)
ax3.set_xlabel('Index')
ax3.set_ylabel('Eigenvalue')
ax3.set_title('Eigenvalue Spectrum')
ax3.grid(True, alpha=0.3)

# Condition numbers vs matrix size
ax4 = axes[1, 0]
ax4.semilogy(sizes, cond_nums, 'bs-', linewidth=2, markersize=8)
ax4.set_xlabel('Matrix Size')
ax4.set_ylabel('Condition Number')
ax4.set_title('Condition Number Growth')
ax4.grid(True, alpha=0.3)

# Eigenvalue distribution comparison
ax5 = axes[1, 1]
colors = ['blue', 'red', 'green']
labels = ['Tridiagonal', 'Random Symmetric', 'Sparse']
for i, (name, (_, evals, _)) in enumerate(eigen_results.items()):
    if i < len(colors):
        ax5.hist(evals, bins=20, alpha=0.6, label=labels[i],
                color=colors[i], density=True)

ax5.set_xlabel('Eigenvalue')
ax5.set_ylabel('Density')
ax5.set_title('Eigenvalue Distributions')
ax5.legend()
ax5.grid(True, alpha=0.3)

# Matrix exponential visualization
ax6 = axes[1, 2]
im6 = ax6.imshow(A_exp, cmap='viridis', interpolation='nearest')
ax6.set_title('Matrix Exponential exp(A)')
plt.colorbar(im6, ax=ax6, shrink=0.8)

# Sparse matrix pattern
ax7 = axes[2, 0]
A_sparse_viz = A_sparse.toarray()[:50, :50]  # Show part of the matrix
ax7.spy(A_sparse_viz, markersize=1)
ax7.set_title('Sparse Matrix Pattern')

# Iterative solver convergence
ax8 = axes[2, 1]
# Create convergence history (simplified)
iterations = np.arange(1, 21)
cg_residuals = np.exp(-0.5 * iterations)
gmres_residuals = np.exp(-0.3 * iterations)

ax8.semilogy(iterations, cg_residuals, 'b-', linewidth=2, label='CG')
ax8.semilogy(iterations, gmres_residuals, 'r--', linewidth=2, label='GMRES')
ax8.set_xlabel('Iteration')
ax8.set_ylabel('Residual Norm')
ax8.set_title('Iterative Solver Convergence')
ax8.legend()
ax8.grid(True, alpha=0.3)

# QR decomposition visualization
ax9 = axes[2, 2]
Q, R = decompositions['QR']
# Show R matrix structure
im9 = ax9.imshow(R, cmap='RdBu_r', interpolation='nearest')
ax9.set_title('R Matrix from QR Decomposition')
plt.colorbar(im9, ax=ax9, shrink=0.8)

plt.tight_layout()
plt.savefig('assets/linear-algebra-analysis.pdf', dpi=300, bbox_inches='tight')
plt.close()

print("Linear algebra analysis saved to assets/linear-algebra-analysis.pdf")
\end{pycode}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.95\textwidth]{assets/linear-algebra-analysis.pdf}
    \caption{Comprehensive linear algebra analysis including matrix structures, eigenvalue and singular value spectra, condition number growth, eigenvalue distributions, matrix functions, sparse matrix patterns, iterative solver convergence, and QR decomposition visualization.}
    \label{fig:linear-algebra-analysis}
\end{figure}

\section{Conclusion}
\label{sec:conclusion}

This linear algebra template demonstrates advanced matrix computational methods including:

\begin{itemize}
    \item Matrix factorizations (LU, QR, Cholesky, SVD)
    \item Eigenvalue algorithms and spectral analysis
    \item Condition number analysis and numerical stability
    \item Iterative methods for large sparse systems
    \item Matrix functions and computational applications
\end{itemize}

These methods form the foundation of modern scientific computing and data analysis applications.

\printbibliography

\end{document}