Sparse sampling

Sparse sampling#

In this page, we describe how to infer IR expansion coefficients using the sparse-sampling techiniques.

using SparseIR
using Plots
gr() # USE GR backend
using OMEinsum
using LaTeXStrings

Setup#

We consider a semicircular spectral modeli (full bandwidth of 2):

\[ \rho(\omega) = \frac{2}{\pi}\sqrt{1-\omega^2}. \]

First, we compute the numerically exact expansion coefficients \(g_l\). Below, we plot the data for even \(l\).

function rho(omega)
    return abs.(omega) < 1 ? (2/π) .* sqrt.(1-omega^2) : 0.0
end

beta = 10000.0
wmax = 1.0
eps = 1e-15 # cutoff for SVD
basis = FiniteTempBasis(Fermionic(), beta, wmax, eps)

rhol = overlap(basis.v, rho)
gl = - basis.s .* rhol

ls = collect(0:length(basis)-1)
p = plot(marker=:x, yaxis=:log, ylabel=L"|g_l|", xlabel=L"l", ylims=(1e-15,10))
plot!(p, ls[1:2:end], abs.(gl[1:2:end]), marker=:x, yaxis=:log, ylabel=L"|g_l|", xlabel=L"l")

From sampling times#

We first create a TauSampling object for the default sampling times.

smpl_tau = TauSampling(basis)
println("sampling times: ", smpl_tau.sampling_points)
println("Condition number: ", SparseIR.cond(smpl_tau))

sampling times: 

[0.06472594109241392, 0.34174300046152517, 0.8430214467258779, 1.573743570455477, 2.5414665603407105, 3.7564568454939717, 5.2320383801574, 6.985039124735026, 9.036355822161779, 11.41165361439933, 14.14221593762588, 17.26595749197979, 20.828608291441064, 24.885070293481082, 29.500941448261653, 34.75419798499302, 40.73702793479739, 47.557820265421434, 55.343334912598685, 64.24110556933593, 74.42215209472747, 86.08409590177114, 99.45477716380213, 114.79647003097449, 132.410786798014, 152.64435851652047, 175.8953792526291, 202.62110258773336, 233.34637878568708, 268.673314868072, 309.29212188731793, 355.99317607327816, 409.6802519273501, 471.3847692357287, 542.2807084296078, 623.6995557013659, 717.1441939083745, 824.2999960143416, 947.0404299720092, 1087.4231696403858, 1247.6709717638912, 1430.1294543119952, 1637.191623041102, 1871.1771467179162, 2134.1542045122987, 2427.695275407337, 2752.5682814454813, 3108.383597604443, 3493.245653743826, 3903.489199194745, 4333.600190763889, 4776.409486578611, 5223.590513421389, 5666.399809236111, 6096.510800805255, 6506.754346256173, 6891.616402395557, 7247.431718554519, 7572.304724592664, 7865.845795487701, 8128.822853282084, 8362.808376958897, 8569.870545688005, 8752.329028236109, 8912.576830359614, 9052.959570027992, 9175.700003985658, 9282.855806091626, 9376.300444298635, 9457.719291570393, 9528.615230764271, 9590.31974807265, 9644.00682392672, 9690.707878112682, 9731.326685131928, 9766.653621214313, 9797.378897412267, 9824.104620747372, 9847.355641483478, 9867.589213201985, 9885.203529969025, 9900.545222836197, 9913.915904098229, 9925.577847905273, 9935.758894430664, 9944.656665087401, 9952.442179734579, 9959.262972065202, 9965.245802015006, 9970.499058551739, 9975.114929706519, 9979.171391708558, 9982.73404250802, 9985.857784062375, 9988.5883463856, 9990.963644177838, 9993.014960875265, 9994.767961619842, 9996.243543154505, 9997.45853343966, 9998.426256429544, 9999.156978553274, 9999.658256999539, 9999.935274058906]
Condition number: 51.83054439129436

The condition number is around 50, indicating that 1–2 significant digits may be lost in a fit from the sampling times. Let us fit from the sampling times!

# Evaluate G(τ) on the sampling times
gtau_smpl = evaluate(smpl_tau, gl)

plot(smpl_tau.sampling_points, gtau_smpl, marker=:x, xlabel=L"\tau", ylabel=L"G(\tau)")

# Fit G(τ) on the sampling times
gl_reconst_from_tau = fit(smpl_tau, gtau_smpl)

104-element Vector{Float64}:
 -0.6509087273407415
 -3.3183447444929196e-16
 -0.5874748878888212
 -7.977356297241e-17
 -0.2696953584728733
  7.322769235602877e-17
 -0.08358548387824463
  8.525894833882885e-17
 -0.01652784880711202
 -2.0939077972599932e-16
  0.0010656038792228421
  3.264126165539172e-17
  0.0034739976272101097
  ⋮
  1.6653345369377348e-16
 -2.9095581738185217e-18
  2.220446049250313e-16
 -3.294767582406096e-17
 -2.498001805406602e-16
 -2.8855511425109667e-16
  1.942890293094024e-16
 -3.70395649373044e-16
  2.7755575615628914e-17
  1.922470975713683e-16
 -1.8041124150158794e-16
  2.0812160673365715e-16

From sampling frequencies#

We create a MatsubaraSampling object for the default sampling frequencies.

smpl_matsu = MatsubaraSampling(basis)
println("sampling frequencies: ", smpl_matsu.sampling_points)
println("Condition number: ", SparseIR.cond(smpl_matsu))

sampling frequencies: FermionicFreq

[FermionicFreq(-45241), FermionicFreq(-14981), FermionicFreq(-8871), FermionicFreq(-6213), FermionicFreq(-4709), FermionicFreq(-3733), FermionicFreq(-3041), FermionicFreq(-2523), FermionicFreq(-2121), FermionicFreq(-1797), FermionicFreq(-1533), FermionicFreq(-1315), FermionicFreq(-1131), FermionicFreq(-975), FermionicFreq(-843), FermionicFreq(-729), FermionicFreq(-631), FermionicFreq(-547), FermionicFreq(-475), FermionicFreq(-411), FermionicFreq(-357), FermionicFreq(-309), FermionicFreq(-269), FermionicFreq(-233), FermionicFreq(-203), FermionicFreq(-177), FermionicFreq(-153), FermionicFreq(-133), FermionicFreq(-115), FermionicFreq(-101), FermionicFreq(-87), FermionicFreq(-75), FermionicFreq(-65), FermionicFreq(-57), FermionicFreq(-49), FermionicFreq(-43), FermionicFreq(-37), FermionicFreq(-33), FermionicFreq(-29), FermionicFreq(-25), FermionicFreq(-23), FermionicFreq(-21), FermionicFreq(-19), FermionicFreq(-17), FermionicFreq(-15), FermionicFreq(-13), FermionicFreq(-11), FermionicFreq(-9), FermionicFreq(-7), FermionicFreq(-5), FermionicFreq(-3), FermionicFreq(-1), FermionicFreq(1), FermionicFreq(3), FermionicFreq(5), FermionicFreq(7), FermionicFreq(9), FermionicFreq(11), FermionicFreq(13), FermionicFreq(15), FermionicFreq(17), FermionicFreq(19), FermionicFreq(21), FermionicFreq(23), FermionicFreq(25), FermionicFreq(29), FermionicFreq(33), FermionicFreq(37), FermionicFreq(43), FermionicFreq(49), FermionicFreq(57), FermionicFreq(65), FermionicFreq(75), FermionicFreq(87), FermionicFreq(101), FermionicFreq(115), FermionicFreq(133), FermionicFreq(153), FermionicFreq(177), FermionicFreq(203), FermionicFreq(233), FermionicFreq(269), FermionicFreq(309), FermionicFreq(357), FermionicFreq(411), FermionicFreq(475), FermionicFreq(547), FermionicFreq(631), FermionicFreq(729), FermionicFreq(843), FermionicFreq(975), FermionicFreq(1131), FermionicFreq(1315), FermionicFreq(1533), FermionicFreq(1797), FermionicFreq(2121), FermionicFreq(2523), FermionicFreq(3041), FermionicFreq(3733), FermionicFreq(4709), FermionicFreq(6213), FermionicFreq(8871), FermionicFreq(14981), FermionicFreq(45241)]
Condition number: 213.18566918573023

The condition number is slightly larger than that for the sampling times.

# Evaluate G(iv) on the sampling frequencies
giv_smpl = evaluate(smpl_matsu, gl)

# `value` function evaluate the actual values of Matsubara frequencies
plot(SparseIR.value.(smpl_matsu.ωn, beta), imag.(giv_smpl), marker=:x, xlabel=L"\nu", ylabel=L"\mathrm{Im} G(i\nu)")

# Fit G(τ) on the sampling times
gl_reconst_from_matsu = fit(smpl_matsu, giv_smpl)

104-element Vector{ComplexF64}:
     -0.6509087273407425 + 2.773486906791006e-15im
   7.018853814128034e-17 - 1.452830911130576e-17im
     -0.5874748878888228 + 4.542156581410772e-15im
   1.592692991381206e-16 - 1.452830911130576e-17im
     -0.2696953584728744 + 5.3399125399256064e-15im
  -7.383416794626285e-17 + 2.3635607360183997e-17im
    -0.08358548387824642 + 5.432720245890366e-15im
   4.412702842016003e-17 + 7.364443256607789e-17im
   -0.016527848807113325 + 5.54460991009087e-15im
   -9.12762703961234e-17 + 1.0169816377914032e-16im
    0.001065603879221143 + 5.5721486452720015e-15im
  1.4994516045474526e-16 + 1.140038584368508e-16im
    0.003473997627208071 + 5.5957842526321855e-15im
                         ⋮
  -5.806322179666563e-16 + 5.773556562986603e-16im
  -2.197431277335797e-15 - 7.90671057306644e-16im
  -5.599263718395718e-16 + 3.9831447735436636e-16im
 -2.1603282563958535e-15 - 6.952653160041986e-16im
 -3.7160004088853344e-16 + 3.2660716943398965e-16im
 -2.4435663066090185e-15 - 8.239861601413561e-16im
  -3.114725997994091e-16 + 2.1061628419944358e-16im
  -1.950409131183608e-15 - 5.906662860651165e-16im
 -1.3543461794222622e-16 + 1.1544757811897022e-16im
 -2.4661072020534442e-15 - 1.0319310777874495e-15im
 -3.3629483993451163e-17 + 3.9522357451370146e-17im
  -1.947528672874772e-15 - 6.880480873012983e-16im

Comparison with exact results#

We now compare the reconstructed expansion coefficients with the exact one. For clarity, we plot only the data for even \(l\).

p = plot(xlabel=L"l", ylabel=L"g_l", ylims=(1e-17, 10), yaxis=:log)
plot!(p, ls[1:2:end], abs.(gl[1:2:end]), marker=:none, label="Exact")
plot!(p, ls[1:2:end], abs.(gl_reconst_from_tau[1:2:end]), marker=:x, label="from sampling times")
plot!(p, ls[1:2:end], abs.(gl_reconst_from_matsu[1:2:end]), marker=:+, label="from sampling frequencies")

We saw a perfect match! Let us plot the differences from the exact one.

p = plot(xlabel=L"L", ylabel=L"Error in $g_l$", ylims=(1e-18, 10), yaxis=:log)
plot!(p, ls[1:2:end], abs.((gl_reconst_from_tau-gl)[1:2:end]), marker=:x, label="from sampling times")
plot!(p, ls[1:2:end], abs.((gl_reconst_from_matsu-gl)[1:2:end]), marker=:+, label="from sampling frequencies")
Contents