main.pytxcode

=>PYTHONTEX#py#default#default#0#code#####136#
import numpy as np
import matplotlib.pyplot as plt

# Set random seed for reproducibility
np.random.seed(42)

def generate_test_matrix(n, condition_number):
    """Generate a symmetric positive definite test matrix with a given condition number."""
    lambda_min = 1.0
    lambda_max = condition_number * lambda_min
    eigenvals = np.logspace(np.log10(lambda_min), np.log10(lambda_max), n)
    Q, _ = np.linalg.qr(np.random.randn(n, n))
    A = Q @ np.diag(eigenvals) @ Q.T
    return A, eigenvals

def conjugate_gradient(A, b, x0, tol=1e-8, max_iter=None):
    """Classical conjugate gradient method."""
    n = len(b)
    if max_iter is None:
        max_iter = n
    x = x0.copy()
    r = b - A @ x
    p = r.copy()
    residuals = [np.linalg.norm(r)]
    for k in range(max_iter):
        Ap = A @ p
        alpha = (r.T @ r) / (p.T @ Ap)
        x = x + alpha * p
        r_new = r - alpha * Ap
        residual_norm = np.linalg.norm(r_new)
        residuals.append(residual_norm)
        if residual_norm < tol:
            break
        beta = (r_new.T @ r_new) / (r.T @ r)
        p = r_new + beta * p
        r = r_new
    return x, residuals, k+1

def accelerated_cg(A, b, x0, tol=1e-8, max_iter=None):
    """Accelerated CG with a simple adaptive momentum term."""
    n = len(b)
    if max_iter is None:
        max_iter = n
    x = x0.copy()
    r = b - A @ x
    p = r.copy()
    r_old = r.copy()
    residuals = [np.linalg.norm(r)]
    for k in range(max_iter):
        Ap = A @ p
        alpha = (r.T @ r) / (p.T @ Ap)
        x = x + alpha * p
        r_new = r - alpha * Ap
        residual_norm = np.linalg.norm(r_new)
        residuals.append(residual_norm)
        if residual_norm < tol:
            break
        beta_cg = (r_new.T @ r_new) / (r.T @ r)
        if k > 0:
            accel_term = 0.1 * (r_new.T @ r_old) / (np.linalg.norm(r_old)**2 + 1e-12)
            beta_acc = max(0.0, beta_cg + accel_term)
        else:
            beta_acc = beta_cg
        p = r_new + beta_acc * p
        r_old = r.copy()
        r = r_new
    return x, residuals, k+1

# Problem set
condition_numbers = [10, 100, 1000]
matrix_size = 50
problems = []
for kappa in condition_numbers:
    A, eigenvals = generate_test_matrix(matrix_size, kappa)
    x_true = np.random.randn(matrix_size)
    b = A @ x_true
    x0 = np.zeros(matrix_size)
    problems.append({
        'A': A,
        'b': b,
        'x0': x0,
        'x_true': x_true,
        'kappa': kappa,
        'eigenvals': eigenvals
    })

# Run solvers and collect results
results = []
for prob in problems:
    A = prob['A']; b = prob['b']; x0 = prob['x0']; x_true = prob['x_true']
    x_cg, res_cg, it_cg = conjugate_gradient(A, b, x0, tol=1e-8, max_iter=matrix_size)
    x_acc, res_acc, it_acc = accelerated_cg(A, b, x0, tol=1e-8, max_iter=matrix_size)
    results.append({
        'kappa': prob['kappa'],
        'residuals_cg': res_cg,
        'residuals_acc': res_acc,
        'iters_cg': it_cg,
        'iters_acc': it_acc,
        'err_cg': np.linalg.norm(x_cg - x_true),
        'err_acc': np.linalg.norm(x_acc - x_true),
        'eigenvals': prob['eigenvals']
    })

print(f"Generated {len(problems)} problems and built {len(results)} result entries.")
=>PYTHONTEX#py#default#default#1#i#####243#
len(problems)
=>PYTHONTEX#py#default#default#2#i#####243#
matrix_size
=>PYTHONTEX#py#default#default#3#i#####243#
matrix_size
=>PYTHONTEX#py#default#default#4#i#####243#
min(condition_numbers)
=>PYTHONTEX#py#default#default#5#i#####243#
max(condition_numbers)
=>PYTHONTEX#py#default#default#6#code#####247#
import numpy as np
import matplotlib.pyplot as plt

# Build a 2x2 summary figure and save it
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
colors = list(plt.cm.tab10.colors)

# (0,0) Convergence curves
ax = axes[0, 0]
for i, result in enumerate(results):
    kappa = result['kappa']
    ax.semilogy(result['residuals_cg'], label=f'CG κ={kappa}', color=colors[i%len(colors)], linestyle='-')
    ax.semilogy(result['residuals_acc'], label=f'AACG κ={kappa}', color=colors[i%len(colors)], linestyle='--')
ax.set_xlabel('Iteration')
ax.set_ylabel('Residual norm')
ax.set_title('Convergence')
ax.grid(True, which='both', ls=':')
ax.legend()

# (0,1) Iteration counts
ax = axes[0, 1]
xpos = np.arange(len(results))
width = 0.35
iters_cg = [r['iters_cg'] for r in results]
iters_acc = [r['iters_acc'] for r in results]
labels = [f'κ={r["kappa"]}' for r in results]
ax.bar(xpos - width/2, iters_cg, width, label='CG')
ax.bar(xpos + width/2, iters_acc, width, label='AACG')
ax.set_xticks(xpos)
ax.set_xticklabels(labels)
ax.set_ylabel('Iterations')
ax.set_title('Iteration counts to tolerance')
ax.grid(axis='y', ls=':', alpha=0.5)
ax.legend()

# (1,0) Final solution errors
ax = axes[1, 0]
err_cg = [r['err_cg'] for r in results]
err_acc = [r['err_acc'] for r in results]
ax.bar(xpos - width/2, err_cg, width, label='CG')
ax.bar(xpos + width/2, err_acc, width, label='AACG')
ax.set_xticks(xpos)
ax.set_xticklabels(labels)
ax.set_ylabel('||x - x*||_2')
ax.set_yscale('log')
ax.set_title('Final solution error')
ax.grid(axis='y', which='both', ls=':', alpha=0.5)
ax.legend()

# (1,1) Example eigenvalue spectra
ax = axes[1, 1]
for i, r in enumerate(results):
    ev = np.sort(r['eigenvals'])
    ax.semilogy(ev, marker='.', linestyle='none', label=f'κ={r["kappa"]}', color=colors[i%len(colors)])
ax.set_xlabel('Index')
ax.set_ylabel('Eigenvalue')
ax.set_title('Eigenvalue spectra')
ax.grid(True, which='both', ls=':')
ax.legend()

fig.tight_layout()
fig.savefig('conv_summary.pdf', bbox_inches='tight')
plt.close(fig)
=>PYTHONTEX#py#default#default#7#i#####320#
f"{avg_improvement:.1f}"
=>PYTHONTEX#py#default#default#8#code#####324#
# Create comprehensive convergence plots
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Plot 1: Convergence curves for different condition numbers
ax = axes[0, 0]
colors = ['blue', 'red', 'green']
for i, result in enumerate(results):
    kappa = result['kappa']
    iterations_cg = range(len(result['residuals_cg']))
    iterations_acc = range(len(result['residuals_acc']))

    ax.semilogy(iterations_cg, result['residuals_cg'],
                color=colors[i], linestyle='-', alpha=0.7,
                label=f'CG κ={kappa}')
    ax.semilogy(iterations_acc, result['residuals_acc'],
                color=colors[i], linestyle='--', linewidth=2,
                label=f'ACC κ={kappa}')

ax.set_xlabel('Iteration')
ax.set_ylabel('Residual Norm')
ax.set_title('Convergence Comparison')
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)

# Plot 2: Iteration count vs condition number
ax = axes[0, 1]
kappas = [r['kappa'] for r in results]
iters_cg = [r['iter_cg'] for r in results]
iters_acc = [r['iter_acc'] for r in results]

ax.semilogx(kappas, iters_cg, 'bo-',
            label='Classical CG', linewidth=2)
ax.semilogx(kappas, iters_acc, 'ro-',
            label='Accelerated CG', linewidth=2)
ax.set_xlabel('Condition Number')
ax.set_ylabel('Iterations to Convergence')
ax.set_title('Iteration Count vs Condition Number')
ax.legend()
ax.grid(True, alpha=0.3)

# Plot 3: Performance improvement
ax = axes[1, 0]
improvements = [r['improvement'] for r in results]
ax.bar(range(len(kappas)), improvements, color='green', alpha=0.7)
ax.set_xlabel('Problem Index')
ax.set_ylabel('Improvement (%)')
ax.set_title('Iteration Reduction by Problem')
ax.set_xticks(range(len(kappas)))
ax.set_xticklabels([f'κ={k}' for k in kappas])
ax.grid(True, alpha=0.3)

# Plot 4: Theoretical vs observed convergence
ax = axes[1, 1]
for i, result in enumerate(results):
    kappa = result['kappa']
    # Theoretical CG bound
    theoretical_rate = (np.sqrt(kappa) - 1) / (np.sqrt(kappa) + 1)

    # Empirical convergence rate (from residuals)
    residuals = np.array(result['residuals_cg'][1:])  # Skip initial residual
    if len(residuals) > 1:
        empirical_rate = np.mean(residuals[1:] / residuals[:-1])
    else:
        empirical_rate = 0

    ax.scatter(theoretical_rate, empirical_rate, s=100,
               label=f'κ={kappa}',
               color=colors[i])

# Add diagonal line
rate_range = np.linspace(0, 1, 100)
ax.plot(rate_range, rate_range, 'k--', alpha=0.5,
        label='Theory = Practice')
ax.set_xlabel('Theoretical Convergence Rate')
ax.set_ylabel('Empirical Convergence Rate')
ax.set_title('Theory vs Practice')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('convergence_analysis.pdf', dpi=300,
            bbox_inches='tight')
plt.savefig('convergence_analysis.png', dpi=200,
            bbox_inches='tight')
print("Generated comprehensive convergence analysis plots")
=>PYTHONTEX#py#default#default#9#code#####423#
# Analyze spectral properties of test matrices
print("Spectral Analysis of Test Matrices:")
print("=" * 40)

for i, prob in enumerate(problems):
    kappa = prob['kappa']
    eigenvals = prob['eigenvals']

    # Compute spectral statistics
    lambda_min = np.min(eigenvals)
    lambda_max = np.max(eigenvals)
    spectral_gap = np.min(np.diff(np.sort(eigenvals)))

    print(f"\nProblem {i+1} (kappa = {kappa}):")
    print(f"  lambda min = {lambda_min:.2e}")
    print(f"  lambda max = {lambda_max:.2e}")
    print(f"  Smallest spectral gap = {spectral_gap:.2e}")
    print(f"  CG iterations: {results[i]['iter_cg']}")
    print(f"  ACC iterations: {results[i]['iter_acc']}")
    print(f"  Improvement: {results[i]['improvement']:.1f}%")

# Create eigenvalue distribution plot
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

for i, prob in enumerate(problems):
    ax = axes[i]
    eigenvals = prob['eigenvals']
    kappa = prob['kappa']

    ax.hist(eigenvals, bins=20, alpha=0.7,
            color=['blue', 'red', 'green'][i])
    ax.set_xlabel('Eigenvalue')
    ax.set_ylabel('Frequency')
    ax.set_title(f'Eigenvalue Distribution (kappa = {kappa})')
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('eigenvalue_distributions.pdf', dpi=300,
            bbox_inches='tight')
plt.savefig('eigenvalue_distributions.png', dpi=200,
            bbox_inches='tight')
plt.close()
print("\nGenerated eigenvalue distribution analysis")
=>PYTHONTEX#py#default#default#10#i#####469#
f"{results[-1]['improvement']:.1f}"
=>PYTHONTEX#py#default#default#11#code#####500#
# Analyze computational overhead of acceleration
print("Computational Overhead Analysis:")
print("=" * 40)

total_time_cg = sum(r['time_cg'] for r in results)
total_time_acc = sum(r['time_acc'] for r in results)
time_overhead = (total_time_acc - total_time_cg) / total_time_cg * 100

avg_time_per_iter_cg = total_time_cg / \
                      sum(r['iter_cg'] for r in results)
avg_time_per_iter_acc = total_time_acc / \
                       sum(r['iter_acc'] for r in results)

print(f"Total CG time: {total_time_cg:.4f}s")
print(f"Total ACC time: {total_time_acc:.4f}s")
print(f"Time overhead: {time_overhead:.1f}%")
print(f"Avg time per CG iteration: "
      f"{avg_time_per_iter_cg*1000:.2f}ms")
print(f"Avg time per ACC iteration: "
      f"{avg_time_per_iter_acc*1000:.2f}ms")

# Overall efficiency metric
iteration_reduction = avg_improvement
time_increase = time_overhead
net_efficiency = iteration_reduction - time_increase

print(f"\nNet efficiency gain: {net_efficiency:.1f}%")
print(f"(Iteration reduction: {iteration_reduction:.1f}%, "
      f"Time overhead: {time_overhead:.1f}%)")
=>PYTHONTEX#py#default#default#12#i#####532#
f"{time_overhead:.1f}"
=>PYTHONTEX#py#default#default#13#i#####532#
f"{net_efficiency:.1f}"
=>PYTHONTEX#py#default#default#14#code#####562#
# Analyze convergence detection strategies
print("Convergence Detection Analysis:")
print("=" * 35)

for result in results:
    kappa = result['kappa']
    residuals = np.array(result['residuals_cg'])

    # Different stopping criteria
    rel_residual = residuals / residuals[0]
    abs_residual = residuals

    # Find when different criteria would stop
    tol_abs = 1e-10
    tol_rel = 1e-8

    stop_abs = np.where(abs_residual < tol_abs)[0]
    stop_rel = np.where(rel_residual < tol_rel)[0]

    iter_abs = stop_abs[0] if len(stop_abs) > 0 else len(residuals)
    iter_rel = stop_rel[0] if len(stop_rel) > 0 else len(residuals)

    print(f"kappa = {kappa}:")
    print(f"  Absolute tol ({tol_abs}): {iter_abs} iterations")
    print(f"  Relative tol ({tol_rel}): {iter_rel} iterations")
    print(f"  Actual convergence: {result['iter_cg']} iterations")
=>PYTHONTEX#py#default#default#15#i#####596#
f"{avg_improvement:.1f}"
=>PYTHONTEX#py#default#default#16#i#####597#
f"{time_overhead:.1f}"
=>PYTHONTEX#py#default#default#17#i#####597#
f"{net_efficiency:.1f}"
=>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|