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.83054439129438

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.6509087273407412
  2.42518238293115e-16
 -0.5874748878888223
  1.2025834971937654e-16
 -0.2696953584728711
 -1.5092433054177123e-17
 -0.08358548387824541
  3.2134735952933646e-16
 -0.016527848807110883
  1.594175299038901e-16
  0.0010656038792224388
 -1.6729408928040784e-16
  0.0034739976272102234
  ⋮
 -8.326672684688674e-17
 -1.0248526914537637e-16
  5.551115123125783e-16
 -2.196828653098495e-16
 -6.661338147750939e-16
 -1.883475167008871e-16
  5.273559366969494e-16
 -1.0703267140182949e-16
 -2.0816681711721685e-16
  4.448092937343472e-17
  3.3306690738754696e-16
  3.512195646010355e-18

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.18566918573356

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.6509087273407417 + 8.40318513949398e-16im
  5.0415401020575956e-17 - 1.6263032587282567e-17im
     -0.5874748878888223 + 1.3788341128584403e-15im
   4.672911363412524e-17 - 8.047490625273657e-17im
     -0.2696953584728738 + 1.6246769554695284e-15im
   -2.64924800846833e-16 - 9.768661574094395e-17im
     -0.0835854838782459 + 1.6622445607461511e-15im
  -3.855422925358454e-16 - 2.09576279941448e-16im
   -0.016527848807112624 + 1.6400726263188226e-15im
  -4.675079767757495e-16 - 2.4549047690503034e-16im
   0.0010656038792221515 + 1.7054500173196985e-15im
  -3.877106968808164e-16 - 2.8492833092919057e-16im
   0.0034739976272091136 + 1.655143036516371e-15im
                         ⋮
 -4.1777766659082467e-16 + 1.921459896307347e-16im
   -3.89107636852318e-15 - 2.3112087529826044e-15im
  -4.138786746729396e-16 + 5.99713521075356e-17im
  -3.927040732336552e-15 - 2.297957981313995e-15im
  -3.892348664949515e-16 + 8.63229693723196e-17im
  -4.316483145064638e-15 - 2.585893107152371e-15im
 -2.6984791929085567e-16 + 5.69921078962174e-17im
 -3.3552126997279296e-15 - 1.9950522859780056e-15im
  -4.638747728779288e-16 + 6.645110672039275e-17im
  -4.628380016760369e-15 - 2.674310700908515e-15im
 -1.8928031907527884e-16 + 9.052785224128706e-18im
  -3.695589757462681e-15 - 2.1599269328678685e-15im

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")