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.06472594109352414, 0.34174300046152517, 0.8430214467258779, 1.5737435704565872, 2.5414665603407105, 3.756456845495082, 5.2320383801574, 6.985039124735026, 9.036355822161779, 11.41165361439933, 14.14221593762588, 17.26595749197979, 20.82860829144051, 24.88507029347886, 29.500941448261653, 34.75419798499302, 40.737027934795165, 47.557820265420325, 55.34333491259758, 64.24110556933593, 74.42215209472747, 86.08409590177058, 99.45477716380103, 114.79647003097226, 132.41078679801288, 152.64435851651993, 175.89537925262744, 202.62110258773114, 233.34637878568597, 268.6733148680698, 309.29212188731793, 355.9931760732793, 409.68025192735126, 471.3847692357265, 542.2807084296055, 623.6995557013647, 717.1441939083734, 824.2999960143404, 947.0404299720082, 1087.4231696403847, 1247.6709717638923, 1430.1294543119957, 1637.1916230411027, 1871.177146717915, 2134.1542045122987, 2427.6952754073354, 2752.56828144548, 3108.383597604443, 3493.245653743827, 3903.4891991947457, 4333.600190763889, 4776.409486578611, 5223.590513421389, 5666.399809236111, 6096.510800805255, 6506.7543462561725, 6891.616402395557, 7247.43171855452, 7572.304724592665, 7865.845795487701, 8128.822853282085, 8362.808376958897, 8569.870545688003, 8752.329028236109, 8912.576830359616, 9052.959570027992, 9175.70000398566, 9282.855806091628, 9376.300444298635, 9457.719291570394, 9528.615230764273, 9590.319748072649, 9644.00682392672, 9690.707878112682, 9731.32668513193, 9766.653621214315, 9797.378897412269, 9824.104620747372, 9847.35564148348, 9867.589213201987, 9885.203529969027, 9900.5452228362, 9913.915904098229, 9925.577847905273, 9935.758894430664, 9944.656665087403, 9952.44217973458, 9959.262972065204, 9965.245802015006, 9970.499058551739, 9975.11492970652, 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.83054439106973

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.6509087273407376
  2.092072236239416e-16
 -0.5874748878888251
 -2.8156356914357664e-16
 -0.26969535847287046
 -1.362855857315814e-16
 -0.08358548387824515
 -2.3272182830160755e-16
 -0.016527848807111285
 -1.2367292954584764e-16
  0.0010656038792232932
  2.98448505041147e-17
  0.0034739976272104515
  ⋮
 -3.2105426469124623e-16
  1.1007397811945174e-16
 -4.490409110140299e-16
 -4.288930746586756e-17
 -4.752011152050265e-16
 -3.115278682145597e-17
 -1.8032138732119027e-16
 -4.802673910674388e-18
  1.4516746547912205e-16
  3.9514263117750846e-17
  4.7193642813299473e-17
  1.0668059985106709e-17

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

The condition number is around 210, about four times as large as 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(\mathrm{i}\nu)")

# Fit G(iν) on the sampling frequencies
gl_reconst_from_matsu = fit(smpl_matsu, giv_smpl)

104-element Vector{ComplexF64}:
     -0.6509087273407425 - 2.9939686890181265e-16im
 -1.0408340855860843e-17 - 2.45029690981724e-17im
     -0.5874748878888317 - 4.837710093630321e-16im
   5.941427905220564e-17 - 8.803721640582296e-17im
     -0.2696953584728777 - 5.660619542546819e-16im
  4.3747557659790104e-17 - 3.8380756905986857e-16im
    -0.08358548387825197 - 5.8223011915187195e-16im
   4.683753385137379e-17 - 5.750608322863116e-16im
   -0.016527848807118536 - 5.5348520905385e-16im
  1.0343288725511712e-16 - 8.222589276130066e-16im
    0.001065603879215595 - 5.399326818977812e-16im
  2.1033522146218786e-17 - 8.427774537272947e-16im
   0.0034739976272030217 - 5.710492842481152e-16im
                         ⋮
  -5.385074430849074e-16 - 2.2184819916038672e-16im
  4.0454767369920255e-16 - 3.8723231914198676e-15im
   -3.79630955476076e-16 - 1.884174498585948e-16im
   3.553603513381371e-16 - 3.925646218320351e-15im
  -3.465552420193592e-16 - 1.5003708999930196e-16im
   4.197232814183331e-16 - 4.534115763816936e-15im
 -1.2156181596555393e-16 - 9.272778708269676e-17im
   3.288421882880508e-16 - 3.3964989194365222e-15im
  -4.487574727065868e-16 - 4.626491310925028e-17im
   5.085984020061702e-16 - 4.548972182390017e-15im
   3.391980350214681e-18 - 1.4794911025394037e-17im
  4.2710581916812243e-16 - 3.748328590226898e-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")