main.pytxcode

=>PYTHONTEX#py#default#default#0#code#####47#
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()
=>PYTHONTEX#py#default#default#1#code#####144#
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()
=>PYTHONTEX#py#default#default#2#code#####241#
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()
=>PYTHONTEX#py#default#default#3#code#####289#
# 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")
=>PYTHONTEX:SETTINGS#
version=0.18
outputdir=pythontex-files-main
workingdir=.
workingdirset=false
gobble=none
rerun=default
hashdependencies=default
makestderr=false
stderrfilename=full
keeptemps=none
pyfuture=default
pyconfuture=none
pygments=true
pygglobal=:GLOBAL||
fvextfile=-1
pyconbanner=none
pyconfilename=stdin
depythontex=false
pygfamily=py|python3|
pygfamily=pycon|pycon|
pygfamily=sympy|python3|
pygfamily=sympycon|pycon|
pygfamily=pylab|python3|
pygfamily=pylabcon|pycon|