13_talbot_spectral_validation.cpp¶
High-precision validation of Talbot contour inversion against Krylov Arnoldi expv using num:: toolkit.
#include <chrono>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <random>
#include <vector>
#include <numerics.hpp>
#include <string_view>
int main(int argc, char **argv) {
const bool plot = argc > 1 && std::string_view(argv[1]) == "--plot";
using namespace num;
std::cout << "========================================================================\n";
std::cout << " NUMERICS: TALBOT CONTOUR VS KRYLOV ARNOILDI EXPV VALIDATION\n";
std::cout << "========================================================================\n\n";
const idx N = 100;
std::mt19937_64 rng(42);
std::uniform_real_distribution<double> weight_dist(0.5, 2.5);
std::uniform_real_distribution<double> prob_dist(0.0, 1.0);
// 1. Construct a connected Markov jump generator Q with row sums = 0
mat<real> Q(N, N, 0.0);
for (idx i = 0; i < N; ++i) {
for (idx j = i + 1; j < N; ++j) {
if (prob_dist(rng) < 0.10 || j == i + 1) { // ensure connected chain
double rate_ij = weight_dist(rng);
double rate_ji = weight_dist(rng);
Q(i, j) = rate_ij;
Q(j, i) = rate_ji;
}
}
}
for (idx j = 0; j < N; ++j) {
double col_rate = 0.0;
for (idx i = 0; i < N; ++i) {
if (i != j) col_rate += Q(i, j);
}
Q(j, j) = -col_rate;
}
vec<real> p0(N, 0.0);
p0[0] = 1.0; // Initial state at node 0
// 2. Precompute Hessenberg decomposition of Q once in O(N^3)
hessenberg_resolvent hess_solver(Q);
const double t_eval = 1.0;
std::cout << "Evaluating Markov diffusion at t = " << t_eval << " (N = " << N << " states)...\n";
// 3. High-precision Ground Truth via Krylov Arnoldi expv
operators::dense_op Q_op(Q);
auto t_krylov_start = std::chrono::high_resolution_clock::now();
vec<real> p_exact = expv(t_eval, Q_op, p0, 50, 1e-15);
auto t_krylov_end = std::chrono::high_resolution_clock::now();
double krylov_ms = std::chrono::duration<double, std::milli>(t_krylov_end - t_krylov_start).count();
std::cout << "Krylov-Arnoldi ground truth computed in " << std::fixed << std::setprecision(2) << krylov_ms << " ms.\n\n";
// 4. Parameter sweep over Talbot quadrature modes M
const std::vector<idx> modes = {4, 6, 8, 10, 12, 14, 16, 20, 24, 28, 32};
std::vector<double> modes_dbl, err_inf_list, err_l1_list;
std::cout << std::left << std::setw(12) << "Nodes (M)"
<< std::setw(18) << "Walltime (ms)"
<< std::setw(18) << "L_inf Error"
<< std::setw(18) << "L_1 Error" << "\n";
std::cout << std::string(66, '-') << "\n";
for (idx M : modes) {
auto t0 = std::chrono::high_resolution_clock::now();
std::vector<cplx> density(N, cplx(0.0, 0.0));
for (const auto &[shift, weight] : talbot_contour(t_eval, M)) {
auto sol = solve(num::shift(hess_solver, shift), p0);
for (idx i = 0; i < N; ++i) {
density[i] += weight * sol[i];
}
}
vec<real> p_talbot(N, 0.0);
for (idx i = 0; i < N; ++i) p_talbot[i] = std::max(0.0, density[i].real());
clip_and_normalize_nonnegative(p_talbot);
auto t1 = std::chrono::high_resolution_clock::now();
double ms = std::chrono::duration<double, std::milli>(t1 - t0).count();
double l_inf = 0.0, l_1 = 0.0;
for (idx i = 0; i < N; ++i) {
double diff = std::abs(p_talbot[i] - p_exact[i]);
l_inf = std::max(l_inf, diff);
l_1 += diff;
}
modes_dbl.push_back(static_cast<double>(M));
err_inf_list.push_back(std::max(1e-16, l_inf));
err_l1_list.push_back(std::max(1e-16, l_1));
std::cout << std::left << std::setw(12) << M
<< std::setw(18) << std::fixed << std::setprecision(3) << ms
<< std::setw(18) << std::scientific << std::setprecision(3) << l_inf
<< std::setw(18) << std::scientific << std::setprecision(3) << l_1 << "\n";
}
// 5. Convergence plot, with --plot
if (plot) {
plt::plot(modes_dbl, err_inf_list, "L_inf Error", "linespoints");
plt::plot(modes_dbl, err_l1_list, "L_1 Error", "linespoints");
plt::title("13 Talbot Spectral Convergence vs Modes M (t = 1.0)");
plt::xlabel("Talbot Quadrature Nodes M");
plt::ylabel("Absolute Error vs Arnoldi expv");
plt::semilogy();
plt::legend();
plt::show();
}
return 0;
}
Possible output:
========================================================================
NUMERICS: TALBOT CONTOUR VS KRYLOV ARNOILDI EXPV VALIDATION
========================================================================
Evaluating Markov diffusion at t = 1 (N = 100 states)...
Krylov-Arnoldi ground truth computed in 0.47 ms.
Nodes (M) Walltime (ms) L_inf Error L_1 Error
------------------------------------------------------------------
4 0.142 2.773e-04 1.117e-03
6 0.187 3.059e-05 7.035e-05
8 0.251 1.327e-06 6.272e-06
10 0.310 1.384e-07 5.747e-07
12 0.380 6.870e-09 4.174e-08
14 0.437 4.816e-10 2.211e-09
16 0.480 2.409e-11 1.619e-10
20 0.605 1.367e-13 1.141e-12
24 0.721 3.197e-15 3.620e-14
28 0.842 4.384e-15 3.555e-14
32 0.972 3.069e-14 6.691e-14