[1]:
import sys
import numpy as np
import matplotlib.pyplot as plt

sys.path.insert(0, '../..')
import ana_cont.continuation as cont

4. Square lattice Hubbard model, \(U=2\), \(\beta=5\), half-filling

4.1. Analytic continuation of self-energy

[2]:
# load data
siw = np.load('data/square_U2_b5/siw.npy')
err = np.load('data/square_U2_b5/siw_err.npy')
beta = 5.
mu = 1.
[3]:
# define the problem
niw = 100  # number of Matsubara frequencies
wn = np.pi / beta * (2. * np.arange(niw) + 1.)  # fermionic Matsubara frequencies

w = np.linspace(-20., 20., num=501, endpoint=True)  # real-frequency grid
model = np.ones_like(w)  # flat default model
model /= np.trapz(model, w)  # normalization, not strictly necessary

# the problem is defined by: imaginary grid, real grid, data, kernel type
probl = cont.AnalyticContinuationProblem(im_axis=wn,  # imaginary grid
                                         re_axis=w,  # real grid
                                         im_data=1j * siw.imag[:niw],  # data
                                         kernel_mode='freq_fermionic')  # kernel type
[4]:
# most basic way to solve the problem
sol, _ = probl.solve(method='maxent_svd',  # Maxent solver that works with singular value decomp.
                     optimizer='newton',  # newton root finding for the optimization problem
                     alpha_determination='chi2kink',
                     model=model,  # default model
                     stdev=err[:niw],  # standard deviation of the data
                     alpha_start=1e12,  # largest alpha, starting value
                     interactive=True)
501 data points on real axis
200 data points on imaginary axis
31 significant singular values
Precomputation of coefficient matrices...
log10(alpha) = 12.00,   chi2 = 4.238e+06,   S = -5.061e-11,   nfev = 2,   norm = 1.000
log10(alpha) = 11.00,   chi2 = 4.236e+06,   S = -5.058e-09,   nfev = 8,   norm = 1.000
log10(alpha) = 10.00,   chi2 = 4.218e+06,   S = -5.030e-07,   nfev = 11,   norm = 1.001
log10(alpha) = 9.00,    chi2 = 4.045e+06,   S = -4.756e-05,   nfev = 14,   norm = 1.008
log10(alpha) = 8.00,    chi2 = 2.898e+06,   S = -2.878e-03,   nfev = 17,   norm = 1.062
log10(alpha) = 7.00,    chi2 = 1.088e+06,   S = -3.482e-02,   nfev = 18,   norm = 1.179
log10(alpha) = 6.00,    chi2 = 2.258e+05,   S = -1.751e-01,   nfev = 22,   norm = 1.133
log10(alpha) = 5.00,    chi2 = 2.230e+04,   S = -4.523e-01,   nfev = 55,   norm = 1.047
log10(alpha) = 4.00,    chi2 = 1.078e+03,   S = -6.940e-01,   nfev = 26,   norm = 1.012
log10(alpha) = 3.00,    chi2 = 8.521e+01,   S = -8.045e-01,   nfev = 32,   norm = 1.003
log10(alpha) = 2.00,    chi2 = 4.769e+01,   S = -8.460e-01,   nfev = 23,   norm = 1.001
log10(alpha) = 1.00,    chi2 = 4.528e+01,   S = -8.820e-01,   nfev = 36,   norm = 1.001
log10(alpha) = 0.00,    chi2 = 4.411e+01,   S = -1.084e+00,   nfev = 43,   norm = 1.000
log10(alpha) = -1.00,   chi2 = 4.371e+01,   S = -1.722e+00,   nfev = 48,   norm = 1.000
log10(alpha) = -2.00,   chi2 = 4.362e+01,   S = -2.824e+00,   nfev = 31,   norm = 1.000
log10(alpha) = -3.00,   chi2 = 4.362e+01,   S = -3.880e+00,   nfev = 39,   norm = 1.000
Fit parameters [1.58559603 5.0305896  4.91072294 1.09647514]
Optimal log alpha 2.6306894961847327
log10(alpha) = 2.63,    chi2 = 5.803e+01,   S = -8.241e-01,   nfev = 17,   norm = 1.002
_images/tutorial_5_1.png
[5]:
# load fitted self-energy
w_fit, sw_fit = np.load('data/square_U2_b5/sw_fit.npy').T

#plot results
fig, ax = plt.subplots(ncols=3, nrows=1, figsize=(15, 4))
ax[0].plot(w, sol.A_opt, label='Maxent solution')
ax[0].plot(w_fit, sw_fit, color='gray', ls='--', label='2-parameter fit')
ax[0].legend()
ax[0].set_xlim(-10., 10.)
ax[0].set_xlabel(r'$\omega$')
ax[1].plot(wn, siw.imag[:niw], marker='x', ls='None')
ax[1].plot(wn, sol.backtransform.imag)
ax[1].set_xlabel(r'$i\omega_n$')
ax[2].plot(wn, (siw[:niw] - sol.backtransform).imag)
ax[2].plot(wn, err[:niw], ls='--', color='gray')
ax[2].plot(wn, -err[:niw], ls='--', color='gray')
ax[2].set_xlabel(r'$i\omega_n$')
plt.tight_layout()
plt.show()
_images/tutorial_6_0.png
[6]:
# more sophisticated solution, using preblur

bw = [0., 0.5, 1., 1.5, 1.75]  # some values for the blur width
bt_list = []  # here we store the backtransform arrays
spec_list = []  # here we store the spectral functions
chi_list= []  # here we store the chi2 deviation

for i in range(len(bw)):  # loop over some values of the blur width
    sol, sol_all = probl.solve(method='maxent_svd',  # Maxent solver that works with singular value decomp.
                         optimizer='newton',  # newton root finding for the optimization problem
                         alpha_determination='chi2kink',
                         model=model,  # default model
                         stdev=err[:niw],  # standard deviation of the data
                         preblur=bw[i]>0.,  # switch on preblur if blur_width > 0.
                         blur_width=bw[i],  # 1. is good, 1.5 is limit
                         alpha_start=1e12,  # largest alpha, starting value
                         interactive=False,
                         verbose=False)
    bt_list.append(sol.backtransform)
    spec_list.append(sol.A_opt)
    chi_list.append(sol_all[-1].chi2)
[7]:
#plot results
fig, ax = plt.subplots(ncols=2, nrows=1, figsize=(15, 4))
for i in range(len(bw)):
    ax[0].plot(w, spec_list[i], label='b={}'.format(bw[i]))
ax[0].legend()
ax[0].set_xlim(-5., 5.)
for i in range(len(bw)):
    ax[1].plot(wn, (siw[:niw] - bt_list[i]).imag,
               label='b={}, min $\chi^2$={:.1f}'.format(bw[i], chi_list[i]))
ax[1].plot(wn, err[:niw], ls='--', color='gray')
ax[1].plot(wn, -err[:niw], ls='--', color='gray')
ax[1].legend()
plt.tight_layout()
plt.show()
_images/tutorial_8_0.png
[8]:
# take b=1., which is list element #2 (in zero-based indexing)
sw = cont.GreensFunction(spectrum=spec_list[2], wgrid=w, kind='fermionic').kkt() + 1.
[9]:
plt.plot(w, sw.real, label='real part')
plt.plot(w, sw.imag, label='imaginary part')
plt.legend()
plt.show()
_images/tutorial_10_0.png

4.2. Calculate spectrum from self-energy and DOS

In cases where the self-energy is proportional to a unit matrix in orbital space, we can calculate the spectral function from the non-interacting DOS \(A_0(\omega)\) and the self-energy in real frequencies as

\[A(\omega) = -\frac{1}{\pi} \mathrm{Im}\; \int_\text{band} d\omega' \frac{A_0(\omega')}{\omega - \omega' + \mu - \Sigma(\omega)}\]
[10]:
from scipy.special import ellipk
w_dos = np.linspace(-5., 5., num=1000, endpoint=True)
dos = np.nan_to_num(ellipk(np.sqrt(1. - (w_dos / 4.)**2)**2) / (2. * np.pi**2))
<ipython-input-10-eed5e07aad46>:3: RuntimeWarning: invalid value encountered in sqrt
  dos = np.nan_to_num(ellipk(np.sqrt(1. - (w_dos / 4.)**2)**2) / (2. * np.pi**2))
[11]:
g_from_s = np.trapz(dos[None, :]
                    / (w[:, None] + mu - w_dos[None, :] - sw[:, None]),
                    w_dos, axis=1)
[12]:
plt.plot(w, -g_from_s.imag / np.pi, color='red', lw=3., label='spectral function')
plt.plot(w_dos, dos, color='gray', ls='--', label='DOS')
plt.xlim(-8., 8.)
plt.legend()
plt.show()
_images/tutorial_14_0.png

4.3. Direct analytic continuation of the Green’s function

[13]:
giw = np.load('data/square_U2_b5/giw.npy')
err_g = np.ones_like(giw, dtype=np.float) * 1e-4
[14]:
niw_g = 100
iw_g = np.pi / beta * (2. * np.arange(niw_g) + 1.)

w_g = np.linspace(-12., 12., num=501, endpoint=True)
model = np.ones_like(w_g)
model /= np.trapz(model, w_g)  # normalization, not strictly necessary

probl = cont.AnalyticContinuationProblem(im_axis=iw_g,
                                         re_axis=w_g,
                                         im_data=1j * giw.imag[:niw_g],
                                         kernel_mode='freq_fermionic')
[15]:
sol, _ = probl.solve(method='maxent_svd',
                     optimizer='newton',
                     alpha_determination='chi2kink',
                     model=model,
                     stdev=err_g[:niw_g],
                     interactive=True,
                     alpha_start=1e11,
                     alpha_end=1e-4)
501 data points on real axis
200 data points on imaginary axis
27 significant singular values
Precomputation of coefficient matrices...
log10(alpha) = 11.00,   chi2 = 2.113e+07,   S = -1.867e-08,   nfev = 6,   norm = 1.000
log10(alpha) = 10.00,   chi2 = 2.106e+07,   S = -1.866e-06,   nfev = 12,   norm = 1.001
log10(alpha) = 9.00,    chi2 = 2.039e+07,   S = -1.857e-04,   nfev = 15,   norm = 1.011
log10(alpha) = 8.00,    chi2 = 1.432e+07,   S = -1.636e-02,   nfev = 18,   norm = 1.102
log10(alpha) = 7.00,    chi2 = 1.698e+06,   S = -2.150e-01,   nfev = 18,   norm = 1.281
log10(alpha) = 6.00,    chi2 = 2.101e+05,   S = -4.029e-01,   nfev = 21,   norm = 1.173
log10(alpha) = 5.00,    chi2 = 1.556e+04,   S = -6.481e-01,   nfev = 22,   norm = 1.056
log10(alpha) = 4.00,    chi2 = 9.667e+02,   S = -8.289e-01,   nfev = 38,   norm = 1.013
log10(alpha) = 3.00,    chi2 = 2.968e+01,   S = -9.280e-01,   nfev = 35,   norm = 1.003
log10(alpha) = 2.00,    chi2 = 1.042e+00,   S = -9.596e-01,   nfev = 20,   norm = 1.000
log10(alpha) = 1.00,    chi2 = 1.699e-02,   S = -9.698e-01,   nfev = 27,   norm = 1.000
log10(alpha) = 0.00,    chi2 = 7.127e-04,   S = -9.713e-01,   nfev = 21,   norm = 1.000
log10(alpha) = -1.00,   chi2 = 2.308e-04,   S = -9.720e-01,   nfev = 19,   norm = 1.000
log10(alpha) = -2.00,   chi2 = 1.796e-04,   S = -9.726e-01,   nfev = 31,   norm = 1.000
log10(alpha) = -3.00,   chi2 = 1.649e-04,   S = -9.754e-01,   nfev = 39,   norm = 1.000
log10(alpha) = -4.00,   chi2 = 1.446e-04,   S = -1.019e+00,   nfev = 29,   norm = 1.000
Fit parameters [-4.28238488 11.83879519  3.22272695  0.5777884 ]
Optimal log alpha -1.1041165907794408
log10(alpha) = -1.10,   chi2 = 2.173e-04,   S = -9.721e-01,   nfev = 11,   norm = 1.000
_images/tutorial_18_1.png
[16]:
fig, ax = plt.subplots(ncols=3, nrows=1, figsize=(15, 4))
ax[0].plot(w_g, sol.A_opt)
ax[1].plot(iw_g, giw.imag[:niw_g])
ax[1].plot(iw_g, sol.backtransform.imag)
ax[2].plot(iw_g, (giw[:niw_g] - sol.backtransform).imag)
plt.tight_layout()
plt.show()
_images/tutorial_19_0.png
[17]:
plt.plot(w, -g_from_s.imag / np.pi, color='red', lw=3., label=r'spectral function from $\Sigma$')
plt.plot(w_g, sol.A_opt, color='green', alpha=1., lw=3., label=r'spectral function from $G$')
plt.plot(w_dos, dos, color='gray', ls='--', label='DOS')
plt.xlim(-8., 8.)
plt.legend()
plt.show()
_images/tutorial_20_0.png

5. Strontium vanadate

[18]:
siw = np.load('data/svo/siw.npy')
err = np.load('data/svo/siw_err.npy')
hartree = np.load('data/svo/hartree.npy')
beta = np.load('data/svo/beta.npy')
[19]:
niw = 300
wn = np.pi / beta * (2. * np.arange(niw) + 1.)

# here we want a non-equispaced grid to cover a larger interval
w = 30 * np.tan(np.linspace(-np.pi/2.1, np.pi/2.1, num=501, endpoint=True)) / np.tan(np.pi / 2.1)
model = np.exp(-(w / 20)**6)  # smooth cut-off at |w| > 20
model /= np.trapz(model, w) / 4.77  # normalization, not strictly necessary

probl = cont.AnalyticContinuationProblem(im_axis=wn,
                                         re_axis=w,
                                         im_data=siw[:niw]-hartree,
                                         kernel_mode='freq_fermionic')
[20]:
sol, _ = probl.solve(method='maxent_svd',
                     optimizer='newton',
                     alpha_determination='chi2kink',
                     model=model,
                     stdev=err[:niw],
                     preblur=True,
                     blur_width=0.3,
                     interactive=True,
                     alpha_start=1e12,
                     alpha_end=1e-1)
501 data points on real axis
600 data points on imaginary axis
51 significant singular values
Precomputation of coefficient matrices...
log10(alpha) = 12.00,   chi2 = 5.181e+07,   S = -3.514e-10,   nfev = 3,   norm = 4.770
log10(alpha) = 11.00,   chi2 = 5.179e+07,   S = -3.513e-08,   nfev = 9,   norm = 4.770
log10(alpha) = 10.00,   chi2 = 5.167e+07,   S = -3.507e-06,   nfev = 12,   norm = 4.772
log10(alpha) = 9.00,    chi2 = 5.042e+07,   S = -3.444e-04,   nfev = 16,   norm = 4.787
log10(alpha) = 8.00,    chi2 = 3.975e+07,   S = -2.854e-02,   nfev = 21,   norm = 4.937
log10(alpha) = 7.00,    chi2 = 9.856e+06,   S = -5.861e-01,   nfev = 53,   norm = 5.640
log10(alpha) = 6.00,    chi2 = 1.282e+06,   S = -1.744e+00,   nfev = 28,   norm = 5.621
log10(alpha) = 5.00,    chi2 = 9.610e+04,   S = -3.252e+00,   nfev = 31,   norm = 5.080
log10(alpha) = 4.00,    chi2 = 6.853e+03,   S = -4.306e+00,   nfev = 27,   norm = 4.853
log10(alpha) = 3.00,    chi2 = 5.785e+02,   S = -5.065e+00,   nfev = 45,   norm = 4.782
log10(alpha) = 2.00,    chi2 = 2.517e+02,   S = -5.413e+00,   nfev = 60,   norm = 4.767
log10(alpha) = 1.00,    chi2 = 2.330e+02,   S = -5.760e+00,   nfev = 64,   norm = 4.766
log10(alpha) = 0.00,    chi2 = 2.174e+02,   S = -8.239e+00,   nfev = 50,   norm = 4.769
log10(alpha) = -1.00,   chi2 = 2.132e+02,   S = -1.370e+01,   nfev = 56,   norm = 4.770
Fit parameters [2.25041422 5.51210904 5.05097152 0.98324766]
Optimal log alpha 2.50837711444609
log10(alpha) = 2.51,    chi2 = 3.032e+02,   S = -5.284e+00,   nfev = 36,   norm = 4.771
_images/tutorial_24_1.png
[21]:
fig, ax = plt.subplots(ncols=3, nrows=1, figsize=(15, 4))
ax[0].plot(w, sol.A_opt, label='spectrum')
ax[0].plot(w, model * sol.norm, color='gray', label='model')
ax[0].legend()
ax[0].set_xlim(-5., 15.)
ax[0].grid()
ax[1].plot(wn, siw.imag[:niw])
ax[1].plot(wn, sol.backtransform.imag)
ax[2].plot(wn, (siw[:niw] - sol.backtransform).imag)
ax[2].plot(wn, err[:niw], ls='--', color='gray')
ax[2].plot(wn, -err[:niw], ls='--', color='gray')
plt.tight_layout()
plt.show()
_images/tutorial_25_0.png

6. Analytical continuation of a matrix-valued function (Kagome)

As an example for matrix-valued analytical continuation we take a D\(\Gamma\)A self-energy for the kagome Hubbard model at \(U=6t\), \(\beta=3/t\), half-filling, for a certain \(\mathbf{k}\)-point.

There, the self-energy is a 3$:nbsphinx-math:times`$3-matrix, with non-zero offdiagonal elements. In this case we are lucky and the projection of the self-energy on the eigenstates of the tight-binding Hamiltonian for this :math:mathbf{k}`-point is diagonal. This means that here we have two ways of doing the analytic continuation: 1. Via projections: * Project self-energy on eigenbasis of tb-Hamiltonian * Continue only the diagonal elements (standard problem) * Rotate back to original basis 2. Direct * Continue the diagonal elements in original basis (standard problem) * construct better default model for offdiag elements * continue the offdiag elements (pos-neg problem)

In the end we can compare both results.

[22]:
# load data
siw = np.load('data/kagome_offdiag/siw.npy')  # DGA self-energy at certain k-point
states = np.load('data/kagome_offdiag/states.npy')  # eigenstates of tb-Hamiltonian at this k-point
hartree = 3.  # for U=6 and half-filling
beta = 3.
niw = siw.shape[0]
wn = np.pi / beta * (2. * np.arange(niw) + 1.)
[23]:
# project self-energy to Hamiltonian eigenbasis
siw_projected = states.T.conj() @ siw @ states

# plot original and projected matrix elements
fig, ax = plt.subplots(ncols=3, nrows=3, figsize=(15, 10))
for i in [0, 1, 2]:
    for j in [0, 1, 2]:
        ax[i, j].plot(wn, siw[:, i, j].real, color='green', ls='--')
        ax[i, j].plot(wn, siw[:, i, j].imag, color='green', ls='-', label='original')
        ax[i, j].plot(wn, siw_projected[:, i, j].real, color='orange', ls='--')
        ax[i, j].plot(wn, siw_projected[:, i, j].imag, color='orange', ls='-', label='projected')
ax[0, 0].legend()
plt.show()
_images/tutorial_28_0.png
[24]:
# define real frequencies, default model and error for continuation
nw = 401
w = np.linspace(-20., 20., num=nw, endpoint=True)
model = np.ones_like(w)
model /= np.trapz(model, w) / 9.  # for U=6 and half-filling, the norm is 9.
wn = np.pi / beta * (2. * np.arange(niw) + 1.)
err = np.ones_like(wn) * 2e-3
sw_diag = np.zeros((nw, 3, 3), dtype=complex)

for i in [0, 1, 2]:  # we have 3 eigenmodes
    # specify the problem
    probl = cont.AnalyticContinuationProblem(im_axis=wn,
                                             re_axis=w,
                                             im_data=siw_projected[:, i, i] - hartree,  # subtract Hartree term
                                             kernel_mode='freq_fermionic')

    # solve the problem
    sol, _ = probl.solve(method='maxent_svd',
                         optimizer='newton',
                         alpha_determination='chi2kink',
                         model=model,
                         stdev=err,
                         interactive=True)

    # get full self-energy from Kramers Kronig
    sw_diag[:, i, i] = cont.GreensFunction(spectrum=sol.A_opt, wgrid=w, kind='fermionic').kkt() + hartree

    # plot
    fig, ax = plt.subplots(ncols=3, nrows=1, figsize=(15, 4))
    ax[0].plot(w, sol.A_opt)
    ax[1].plot(wn, siw_projected[:, i, i].real - hartree, marker='x', ls='None')
    ax[1].plot(wn, siw_projected[:, i, i].imag, marker='+', ls='None')
    ax[1].plot(wn, sol.backtransform.real)
    ax[1].plot(wn, sol.backtransform.imag)
    ax[2].plot(wn, (siw_projected[:, i, i] - hartree - sol.backtransform).real)
    ax[2].plot(wn, (siw_projected[:, i, i] - hartree - sol.backtransform).imag)
    plt.show()

# go back to original basis
sw_rotated_back = states @ sw_diag @ states.T.conj()
401 data points on real axis
80 data points on imaginary axis
26 significant singular values
Precomputation of coefficient matrices...
log10(alpha) = 9.00,    chi2 = 6.783e+05,   S = -4.278e-08,   nfev = 4,   norm = 9.001
log10(alpha) = 8.00,    chi2 = 6.767e+05,   S = -4.272e-06,   nfev = 11,   norm = 9.005
log10(alpha) = 7.00,    chi2 = 6.615e+05,   S = -4.214e-04,   nfev = 14,   norm = 9.050
log10(alpha) = 6.00,    chi2 = 5.274e+05,   S = -3.608e-02,   nfev = 17,   norm = 9.466
log10(alpha) = 5.00,    chi2 = 1.281e+05,   S = -7.818e-01,   nfev = 19,   norm = 11.124
log10(alpha) = 4.00,    chi2 = 2.349e+04,   S = -2.224e+00,   nfev = 20,   norm = 11.403
log10(alpha) = 3.00,    chi2 = 1.953e+03,   S = -5.132e+00,   nfev = 39,   norm = 10.332
log10(alpha) = 2.00,    chi2 = 1.346e+02,   S = -7.279e+00,   nfev = 37,   norm = 9.823
log10(alpha) = 1.00,    chi2 = 2.966e+01,   S = -8.397e+00,   nfev = 35,   norm = 9.684
log10(alpha) = 0.00,    chi2 = 2.604e+01,   S = -8.865e+00,   nfev = 46,   norm = 9.663
log10(alpha) = -1.00,   chi2 = 2.551e+01,   S = -9.576e+00,   nfev = 26,   norm = 9.665
log10(alpha) = -2.00,   chi2 = 2.540e+01,   S = -1.143e+01,   nfev = 28,   norm = 9.666
log10(alpha) = -3.00,   chi2 = 2.535e+01,   S = -1.939e+01,   nfev = 31,   norm = 9.667
Fit parameters [1.34029571 4.51440873 3.34634067 1.12696789]
Optimal log alpha 1.1279988482793266
log10(alpha) = 1.13,    chi2 = 3.146e+01,   S = -8.320e+00,   nfev = 14,   norm = 9.691
_images/tutorial_29_1.png
_images/tutorial_29_2.png
401 data points on real axis
80 data points on imaginary axis
26 significant singular values
Precomputation of coefficient matrices...
log10(alpha) = 9.00,    chi2 = 2.587e+05,   S = -1.459e-08,   nfev = 4,   norm = 9.000
log10(alpha) = 8.00,    chi2 = 2.581e+05,   S = -1.455e-06,   nfev = 10,   norm = 9.004
log10(alpha) = 7.00,    chi2 = 2.530e+05,   S = -1.423e-04,   nfev = 14,   norm = 9.040
log10(alpha) = 6.00,    chi2 = 2.101e+05,   S = -1.138e-02,   nfev = 17,   norm = 9.364
log10(alpha) = 5.00,    chi2 = 8.735e+04,   S = -2.497e-01,   nfev = 18,   norm = 10.733
log10(alpha) = 4.00,    chi2 = 2.374e+04,   S = -1.311e+00,   nfev = 20,   norm = 11.318
log10(alpha) = 3.00,    chi2 = 2.221e+03,   S = -4.292e+00,   nfev = 35,   norm = 10.354
log10(alpha) = 2.00,    chi2 = 1.513e+02,   S = -6.728e+00,   nfev = 37,   norm = 9.831
log10(alpha) = 1.00,    chi2 = 2.791e+01,   S = -8.065e+00,   nfev = 45,   norm = 9.683
log10(alpha) = 0.00,    chi2 = 2.456e+01,   S = -8.398e+00,   nfev = 23,   norm = 9.657
log10(alpha) = -1.00,   chi2 = 2.444e+01,   S = -8.599e+00,   nfev = 52,   norm = 9.655
log10(alpha) = -2.00,   chi2 = 2.438e+01,   S = -9.493e+00,   nfev = 26,   norm = 9.657
log10(alpha) = -3.00,   chi2 = 2.436e+01,   S = -1.312e+01,   nfev = 29,   norm = 9.657
Fit parameters [1.32887755 4.07280426 3.0878616  1.23959631]
Optimal log alpha 1.0710759853132283
log10(alpha) = 1.07,    chi2 = 2.901e+01,   S = -8.015e+00,   nfev = 12,   norm = 9.688
_images/tutorial_29_4.png
_images/tutorial_29_5.png
401 data points on real axis
80 data points on imaginary axis
26 significant singular values
Precomputation of coefficient matrices...
log10(alpha) = 9.00,    chi2 = 2.592e+05,   S = -1.385e-08,   nfev = 4,   norm = 9.000
log10(alpha) = 8.00,    chi2 = 2.587e+05,   S = -1.382e-06,   nfev = 10,   norm = 9.004
log10(alpha) = 7.00,    chi2 = 2.538e+05,   S = -1.352e-04,   nfev = 13,   norm = 9.039
log10(alpha) = 6.00,    chi2 = 2.129e+05,   S = -1.087e-02,   nfev = 17,   norm = 9.350
log10(alpha) = 5.00,    chi2 = 9.153e+04,   S = -2.497e-01,   nfev = 18,   norm = 10.698
log10(alpha) = 4.00,    chi2 = 2.489e+04,   S = -1.360e+00,   nfev = 20,   norm = 11.318
log10(alpha) = 3.00,    chi2 = 2.364e+03,   S = -4.505e+00,   nfev = 36,   norm = 10.354
log10(alpha) = 2.00,    chi2 = 1.469e+02,   S = -7.057e+00,   nfev = 36,   norm = 9.829
log10(alpha) = 1.00,    chi2 = 2.922e+01,   S = -8.340e+00,   nfev = 46,   norm = 9.684
log10(alpha) = 0.00,    chi2 = 2.578e+01,   S = -8.695e+00,   nfev = 36,   norm = 9.658
log10(alpha) = -1.00,   chi2 = 2.558e+01,   S = -9.015e+00,   nfev = 39,   norm = 9.657
log10(alpha) = -2.00,   chi2 = 2.549e+01,   S = -1.046e+01,   nfev = 27,   norm = 9.659
log10(alpha) = -3.00,   chi2 = 2.546e+01,   S = -1.500e+01,   nfev = 30,   norm = 9.659
Fit parameters [1.35107894 4.04894961 3.08105855 1.26503099]
Optimal log alpha 1.1048223759601463
log10(alpha) = 1.10,    chi2 = 3.094e+01,   S = -8.265e+00,   nfev = 13,   norm = 9.690
_images/tutorial_29_7.png
_images/tutorial_29_8.png
[25]:
sw_full = np.zeros_like(sw_diag)


# first loop only over the diagonal elements
for i in [0, 1, 2]:
    print('Diagonal element {}'.format(i))

    # specify the problem
    probl = cont.AnalyticContinuationProblem(im_axis=wn,
                                             re_axis=w,
                                             im_data=siw[:, i, i] - hartree,  # subtract Hartree term
                                             kernel_mode='freq_fermionic')

    # solve the problem
    sol, _ = probl.solve(method='maxent_svd',
                         optimizer='newton',
                         alpha_determination='chi2kink',
                         model=model,
                         stdev=err,
                         interactive=True)

    # construct full self-energy by Kramers Kronig
    sw_full[:, i, i] = cont.GreensFunction(spectrum=sol.A_opt, wgrid=w, kind='fermionic').kkt() + hartree

    # plot
    fig, ax = plt.subplots(ncols=3, nrows=1, figsize=(15, 4))
    ax[0].plot(w, sol.A_opt)
    ax[1].plot(wn, siw[:, i, i].real - hartree, marker='x', ls='None')
    ax[1].plot(wn, siw[:, i, i].imag, marker='+', ls='None')
    ax[1].plot(wn, sol.backtransform.real)
    ax[1].plot(wn, sol.backtransform.imag)
    ax[2].plot(wn, (siw[:, i, i] - hartree - sol.backtransform).real)
    ax[2].plot(wn, (siw[:, i, i] - hartree - sol.backtransform).imag)
    plt.show()

# now loop over the offdiagonal elements
for i in [0, 1, 2]:
    for j in [0, 1, 2]:
        if i == j: continue
        print('Offdiagonal element {} {}'.format(i, j))

        # construct a good default model
        offdiag_model = np.sqrt(sw_full.imag[:, i, i] * sw_full.imag[:, j, j])  # "poorman's method" Kraberger et al.

        # specify the problem
        probl = cont.AnalyticContinuationProblem(im_axis=wn,
                                                 re_axis=w,
                                                 im_data=siw[:, i, j],  # no Hartree term here!
                                                 kernel_mode='freq_fermionic')

        # solve the problem
        sol, _ = probl.solve(method='maxent_svd',
                             optimizer='newton',
                             alpha_determination='chi2kink',
                             model=offdiag_model,  # switch on offdiag mode
                             stdev=err,
                             offdiag=True,
                             alpha_start=1e10,
                             alpha_end=1e-6,
                             fit_position=1.5,  # smaller value to prevent overfitting
                             interactive=True)

        # construct full self-energy by Kramers Kronig
        sw_full[:, i, j] = cont.GreensFunction(spectrum=sol.A_opt, wgrid=w, kind='fermionic').kkt()

        # plot
        fig, ax = plt.subplots(ncols=3, nrows=1, figsize=(15, 4))
        ax[0].plot(w, sol.A_opt)
        ax[1].plot(wn, siw[:, i, j].real, marker='x', ls='None')
        ax[1].plot(wn, siw[:, i, j].imag, marker='+', ls='None')
        ax[1].plot(wn, sol.backtransform.real)
        ax[1].plot(wn, sol.backtransform.imag)
        ax[2].plot(wn, (siw[:, i, j] - sol.backtransform).real)
        ax[2].plot(wn, (siw[:, i, j] - sol.backtransform).imag)
        plt.show()
Diagonal element 0
401 data points on real axis
80 data points on imaginary axis
26 significant singular values
Precomputation of coefficient matrices...
log10(alpha) = 9.00,    chi2 = 3.003e+05,   S = -1.797e-08,   nfev = 4,   norm = 9.000
log10(alpha) = 8.00,    chi2 = 2.996e+05,   S = -1.793e-06,   nfev = 10,   norm = 9.004
log10(alpha) = 7.00,    chi2 = 2.933e+05,   S = -1.755e-04,   nfev = 14,   norm = 9.043
log10(alpha) = 6.00,    chi2 = 2.400e+05,   S = -1.417e-02,   nfev = 17,   norm = 9.389
log10(alpha) = 5.00,    chi2 = 8.949e+04,   S = -3.016e-01,   nfev = 18,   norm = 10.814
log10(alpha) = 4.00,    chi2 = 2.324e+04,   S = -1.381e+00,   nfev = 20,   norm = 11.331
log10(alpha) = 3.00,    chi2 = 2.095e+03,   S = -4.284e+00,   nfev = 44,   norm = 10.348
log10(alpha) = 2.00,    chi2 = 1.421e+02,   S = -6.596e+00,   nfev = 38,   norm = 9.830
log10(alpha) = 1.00,    chi2 = 2.607e+01,   S = -7.838e+00,   nfev = 31,   norm = 9.684
log10(alpha) = 0.00,    chi2 = 2.308e+01,   S = -8.137e+00,   nfev = 28,   norm = 9.657
log10(alpha) = -1.00,   chi2 = 2.298e+01,   S = -8.278e+00,   nfev = 48,   norm = 9.653
log10(alpha) = -2.00,   chi2 = 2.293e+01,   S = -9.278e+00,   nfev = 48,   norm = 9.653
log10(alpha) = -3.00,   chi2 = 2.290e+01,   S = -1.440e+01,   nfev = 29,   norm = 9.652
Fit parameters [1.29842143 4.17127903 3.12916432 1.21176069]
Optimal log alpha 1.0660506905566396
log10(alpha) = 1.07,    chi2 = 2.697e+01,   S = -7.796e+00,   nfev = 12,   norm = 9.688
_images/tutorial_30_1.png
_images/tutorial_30_2.png
Diagonal element 1
401 data points on real axis
80 data points on imaginary axis
26 significant singular values
Precomputation of coefficient matrices...
log10(alpha) = 9.00,    chi2 = 3.010e+05,   S = -1.795e-08,   nfev = 4,   norm = 9.000
log10(alpha) = 8.00,    chi2 = 3.004e+05,   S = -1.791e-06,   nfev = 10,   norm = 9.004
log10(alpha) = 7.00,    chi2 = 2.940e+05,   S = -1.753e-04,   nfev = 14,   norm = 9.043
log10(alpha) = 6.00,    chi2 = 2.408e+05,   S = -1.416e-02,   nfev = 17,   norm = 9.389
log10(alpha) = 5.00,    chi2 = 9.008e+04,   S = -3.022e-01,   nfev = 18,   norm = 10.814
log10(alpha) = 4.00,    chi2 = 2.340e+04,   S = -1.389e+00,   nfev = 20,   norm = 11.332
log10(alpha) = 3.00,    chi2 = 2.099e+03,   S = -4.314e+00,   nfev = 44,   norm = 10.347
log10(alpha) = 2.00,    chi2 = 1.402e+02,   S = -6.623e+00,   nfev = 38,   norm = 9.830
log10(alpha) = 1.00,    chi2 = 2.604e+01,   S = -7.844e+00,   nfev = 31,   norm = 9.684
log10(alpha) = 0.00,    chi2 = 2.308e+01,   S = -8.141e+00,   nfev = 28,   norm = 9.657
log10(alpha) = -1.00,   chi2 = 2.298e+01,   S = -8.284e+00,   nfev = 47,   norm = 9.653
log10(alpha) = -2.00,   chi2 = 2.293e+01,   S = -9.296e+00,   nfev = 46,   norm = 9.653
log10(alpha) = -3.00,   chi2 = 2.290e+01,   S = -1.446e+01,   nfev = 29,   norm = 9.652
Fit parameters [1.29883645 4.17134996 3.12935829 1.2159353 ]
Optimal log alpha 1.0733278533586619
log10(alpha) = 1.07,    chi2 = 2.704e+01,   S = -7.798e+00,   nfev = 12,   norm = 9.689
_images/tutorial_30_4.png
_images/tutorial_30_5.png
Diagonal element 2
401 data points on real axis
80 data points on imaginary axis
26 significant singular values
Precomputation of coefficient matrices...
log10(alpha) = 9.00,    chi2 = 2.994e+05,   S = -1.801e-08,   nfev = 4,   norm = 9.000
log10(alpha) = 8.00,    chi2 = 2.987e+05,   S = -1.797e-06,   nfev = 10,   norm = 9.004
log10(alpha) = 7.00,    chi2 = 2.924e+05,   S = -1.759e-04,   nfev = 14,   norm = 9.043
log10(alpha) = 6.00,    chi2 = 2.390e+05,   S = -1.420e-02,   nfev = 17,   norm = 9.390
log10(alpha) = 5.00,    chi2 = 8.866e+04,   S = -3.010e-01,   nfev = 18,   norm = 10.815
log10(alpha) = 4.00,    chi2 = 2.303e+04,   S = -1.370e+00,   nfev = 20,   norm = 11.330
log10(alpha) = 3.00,    chi2 = 2.089e+03,   S = -4.243e+00,   nfev = 45,   norm = 10.349
log10(alpha) = 2.00,    chi2 = 1.450e+02,   S = -6.559e+00,   nfev = 38,   norm = 9.831
log10(alpha) = 1.00,    chi2 = 2.612e+01,   S = -7.832e+00,   nfev = 35,   norm = 9.684
log10(alpha) = 0.00,    chi2 = 2.307e+01,   S = -8.135e+00,   nfev = 35,   norm = 9.657
log10(alpha) = -1.00,   chi2 = 2.298e+01,   S = -8.274e+00,   nfev = 55,   norm = 9.653
log10(alpha) = -2.00,   chi2 = 2.293e+01,   S = -9.260e+00,   nfev = 49,   norm = 9.653
log10(alpha) = -3.00,   chi2 = 2.290e+01,   S = -1.433e+01,   nfev = 29,   norm = 9.652
Fit parameters [1.29780736 4.17145334 3.12884244 1.205656  ]
Optimal log alpha 1.0552824850877567
log10(alpha) = 1.06,    chi2 = 2.687e+01,   S = -7.797e+00,   nfev = 12,   norm = 9.687
_images/tutorial_30_7.png
_images/tutorial_30_8.png
Offdiagonal element 0 1
401 data points on real axis
80 data points on imaginary axis
26 significant singular values
Precomputation of coefficient matrices...
log10(alpha) = 10.00,   chi2 = 6.525e+04,   S = -4.689e-10,   nfev = 1,   norm = 0.000
log10(alpha) = 9.00,    chi2 = 6.508e+04,   S = -4.676e-08,   nfev = 3,   norm = 0.000
log10(alpha) = 8.00,    chi2 = 6.344e+04,   S = -4.552e-06,   nfev = 9,   norm = 0.003
log10(alpha) = 7.00,    chi2 = 5.003e+04,   S = -3.542e-04,   nfev = 13,   norm = 0.018
log10(alpha) = 6.00,    chi2 = 1.220e+04,   S = -7.678e-03,   nfev = 19,   norm = 0.009
log10(alpha) = 5.00,    chi2 = 9.045e+02,   S = -2.131e-02,   nfev = 21,   norm = -0.079
log10(alpha) = 4.00,    chi2 = 5.157e+01,   S = -3.166e-02,   nfev = 16,   norm = -0.080
log10(alpha) = 3.00,    chi2 = 7.075e+00,   S = -3.665e-02,   nfev = 25,   norm = -0.048
log10(alpha) = 2.00,    chi2 = 2.438e+00,   S = -4.344e-02,   nfev = 16,   norm = -0.029
log10(alpha) = 1.00,    chi2 = 1.600e+00,   S = -5.622e-02,   nfev = 26,   norm = -0.020
log10(alpha) = 0.00,    chi2 = 1.265e+00,   S = -1.187e-01,   nfev = 18,   norm = -0.015
log10(alpha) = -1.00,   chi2 = 9.880e-01,   S = -6.166e-01,   nfev = 20,   norm = -0.010
log10(alpha) = -2.00,   chi2 = 6.846e-01,   S = -6.943e+00,   nfev = 72,   norm = -0.007
log10(alpha) = -3.00,   chi2 = 4.620e-01,   S = -5.035e+01,   nfev = 74,   norm = -0.005
log10(alpha) = -4.00,   chi2 = 2.742e-01,   S = -3.289e+02,   nfev = 88,   norm = -0.005
log10(alpha) = -5.00,   chi2 = 1.697e-01,   S = -2.243e+03,   nfev = 127,   norm = -0.007
log10(alpha) = -6.00,   chi2 = 1.462e-01,   S = -6.135e+03,   nfev = 62,   norm = -0.011
Fit parameters [-0.39097276  5.43464373  4.36226097  0.80648082]
Optimal log alpha 2.502328332228192
log10(alpha) = 2.50,    chi2 = 3.918e+00,   S = -3.943e-02,   nfev = 12,   norm = -0.036
_images/tutorial_30_10.png
_images/tutorial_30_11.png
Offdiagonal element 0 2
401 data points on real axis
80 data points on imaginary axis
26 significant singular values
Precomputation of coefficient matrices...
log10(alpha) = 10.00,   chi2 = 3.497e+04,   S = -2.028e-10,   nfev = 1,   norm = -0.000
log10(alpha) = 9.00,    chi2 = 3.490e+04,   S = -2.023e-08,   nfev = 2,   norm = -0.001
log10(alpha) = 8.00,    chi2 = 3.418e+04,   S = -1.978e-06,   nfev = 9,   norm = -0.006
log10(alpha) = 7.00,    chi2 = 2.817e+04,   S = -1.601e-04,   nfev = 12,   norm = -0.053
log10(alpha) = 6.00,    chi2 = 8.684e+03,   S = -4.093e-03,   nfev = 16,   norm = -0.266
log10(alpha) = 5.00,    chi2 = 8.224e+02,   S = -1.443e-02,   nfev = 16,   norm = -0.337
log10(alpha) = 4.00,    chi2 = 8.786e+01,   S = -2.348e-02,   nfev = 17,   norm = -0.215
log10(alpha) = 3.00,    chi2 = 1.631e+01,   S = -3.355e-02,   nfev = 20,   norm = -0.126
log10(alpha) = 2.00,    chi2 = 2.891e+00,   S = -5.109e-02,   nfev = 23,   norm = -0.076
log10(alpha) = 1.00,    chi2 = 1.253e+00,   S = -7.461e-02,   nfev = 26,   norm = -0.053
log10(alpha) = 0.00,    chi2 = 7.296e-01,   S = -1.677e-01,   nfev = 18,   norm = -0.040
log10(alpha) = -1.00,   chi2 = 4.765e-01,   S = -6.223e-01,   nfev = 19,   norm = -0.029
log10(alpha) = -2.00,   chi2 = 2.909e-01,   S = -3.844e+00,   nfev = 72,   norm = -0.020
log10(alpha) = -3.00,   chi2 = 1.945e-01,   S = -2.190e+01,   nfev = 85,   norm = -0.012
log10(alpha) = -4.00,   chi2 = 1.311e-01,   S = -1.214e+02,   nfev = 40,   norm = -0.007
log10(alpha) = -5.00,   chi2 = 9.567e-02,   S = -8.655e+02,   nfev = 88,   norm = -0.006
log10(alpha) = -6.00,   chi2 = 7.966e-02,   S = -3.558e+03,   nfev = 41,   norm = -0.011
Fit parameters [-0.84341047  5.72437401  3.83747472  0.61814642]
Optimal log alpha 1.410865161600237
log10(alpha) = 1.41,    chi2 = 1.628e+00,   S = -6.305e-02,   nfev = 12,   norm = -0.060
_images/tutorial_30_13.png
_images/tutorial_30_14.png
Offdiagonal element 1 0
401 data points on real axis
80 data points on imaginary axis
26 significant singular values
Precomputation of coefficient matrices...
log10(alpha) = 10.00,   chi2 = 6.523e+04,   S = -4.929e-10,   nfev = 1,   norm = 0.000
log10(alpha) = 9.00,    chi2 = 6.505e+04,   S = -4.915e-08,   nfev = 3,   norm = 0.001
log10(alpha) = 8.00,    chi2 = 6.332e+04,   S = -4.772e-06,   nfev = 9,   norm = 0.009
log10(alpha) = 7.00,    chi2 = 4.951e+04,   S = -3.635e-04,   nfev = 13,   norm = 0.076
log10(alpha) = 6.00,    chi2 = 1.316e+04,   S = -7.272e-03,   nfev = 15,   norm = 0.263
log10(alpha) = 5.00,    chi2 = 1.551e+03,   S = -2.290e-02,   nfev = 15,   norm = 0.203
log10(alpha) = 4.00,    chi2 = 6.559e+01,   S = -4.067e-02,   nfev = 16,   norm = 0.087
log10(alpha) = 3.00,    chi2 = 6.448e+00,   S = -4.680e-02,   nfev = 17,   norm = 0.048
log10(alpha) = 2.00,    chi2 = 2.318e+00,   S = -5.281e-02,   nfev = 16,   norm = 0.029
log10(alpha) = 1.00,    chi2 = 1.561e+00,   S = -6.424e-02,   nfev = 26,   norm = 0.021
log10(alpha) = 0.00,    chi2 = 1.270e+00,   S = -1.196e-01,   nfev = 24,   norm = 0.016
log10(alpha) = -1.00,   chi2 = 9.898e-01,   S = -6.360e-01,   nfev = 38,   norm = 0.012
log10(alpha) = -2.00,   chi2 = 6.704e-01,   S = -7.233e+00,   nfev = 102,   norm = 0.007
log10(alpha) = -3.00,   chi2 = 4.413e-01,   S = -5.117e+01,   nfev = 61,   norm = 0.001
log10(alpha) = -4.00,   chi2 = 2.615e-01,   S = -3.186e+02,   nfev = 83,   norm = -0.002
log10(alpha) = -5.00,   chi2 = 1.667e-01,   S = -2.055e+03,   nfev = 122,   norm = -0.005
log10(alpha) = -6.00,   chi2 = 1.488e-01,   S = -4.818e+03,   nfev = 52,   norm = -0.006
Fit parameters [-0.39264746  5.3892644   4.25688119  0.84949537]
Optimal log alpha 2.4911270061596706
log10(alpha) = 2.49,    chi2 = 3.578e+00,   S = -4.934e-02,   nfev = 12,   norm = 0.037
_images/tutorial_30_16.png
_images/tutorial_30_17.png
Offdiagonal element 1 2
401 data points on real axis
80 data points on imaginary axis
26 significant singular values
Precomputation of coefficient matrices...
log10(alpha) = 10.00,   chi2 = 4.747e+04,   S = -2.987e-10,   nfev = 1,   norm = -0.000
log10(alpha) = 9.00,    chi2 = 4.736e+04,   S = -2.980e-08,   nfev = 3,   norm = -0.000
log10(alpha) = 8.00,    chi2 = 4.631e+04,   S = -2.910e-06,   nfev = 9,   norm = -0.004
log10(alpha) = 7.00,    chi2 = 3.755e+04,   S = -2.328e-04,   nfev = 13,   norm = -0.037
log10(alpha) = 6.00,    chi2 = 1.036e+04,   S = -5.645e-03,   nfev = 16,   norm = -0.211
log10(alpha) = 5.00,    chi2 = 8.102e+02,   S = -1.752e-02,   nfev = 15,   norm = -0.296
log10(alpha) = 4.00,    chi2 = 8.474e+01,   S = -2.655e-02,   nfev = 17,   norm = -0.199
log10(alpha) = 3.00,    chi2 = 1.522e+01,   S = -3.605e-02,   nfev = 27,   norm = -0.116
log10(alpha) = 2.00,    chi2 = 3.075e+00,   S = -5.217e-02,   nfev = 17,   norm = -0.070
log10(alpha) = 1.00,    chi2 = 1.498e+00,   S = -7.510e-02,   nfev = 26,   norm = -0.049
log10(alpha) = 0.00,    chi2 = 9.631e-01,   S = -1.713e-01,   nfev = 18,   norm = -0.036
log10(alpha) = -1.00,   chi2 = 6.728e-01,   S = -6.903e-01,   nfev = 27,   norm = -0.027
log10(alpha) = -2.00,   chi2 = 4.379e-01,   S = -5.034e+00,   nfev = 71,   norm = -0.017
log10(alpha) = -3.00,   chi2 = 2.997e-01,   S = -3.111e+01,   nfev = 85,   norm = -0.010
log10(alpha) = -4.00,   chi2 = 1.906e-01,   S = -1.983e+02,   nfev = 36,   norm = -0.007
log10(alpha) = -5.00,   chi2 = 1.283e-01,   S = -1.436e+03,   nfev = 108,   norm = -0.007
log10(alpha) = -6.00,   chi2 = 1.073e-01,   S = -4.961e+03,   nfev = 64,   norm = -0.012
Fit parameters [-0.65366492  5.66409807  4.03328656  0.64482825]
Optimal log alpha 1.7070857589775366
log10(alpha) = 1.71,    chi2 = 2.324e+00,   S = -5.737e-02,   nfev = 20,   norm = -0.062
_images/tutorial_30_19.png
_images/tutorial_30_20.png
Offdiagonal element 2 0
401 data points on real axis
80 data points on imaginary axis
26 significant singular values
Precomputation of coefficient matrices...
log10(alpha) = 10.00,   chi2 = 3.499e+04,   S = -2.410e-10,   nfev = 1,   norm = 0.000
log10(alpha) = 9.00,    chi2 = 3.490e+04,   S = -2.402e-08,   nfev = 2,   norm = 0.001
log10(alpha) = 8.00,    chi2 = 3.406e+04,   S = -2.331e-06,   nfev = 9,   norm = 0.012
log10(alpha) = 7.00,    chi2 = 2.734e+04,   S = -1.765e-04,   nfev = 12,   norm = 0.101
log10(alpha) = 6.00,    chi2 = 9.609e+03,   S = -3.584e-03,   nfev = 15,   norm = 0.402
log10(alpha) = 5.00,    chi2 = 1.559e+03,   S = -1.607e-02,   nfev = 16,   norm = 0.399
log10(alpha) = 4.00,    chi2 = 1.045e+02,   S = -3.359e-02,   nfev = 19,   norm = 0.220
log10(alpha) = 3.00,    chi2 = 1.491e+01,   S = -4.495e-02,   nfev = 19,   norm = 0.126
log10(alpha) = 2.00,    chi2 = 2.747e+00,   S = -6.085e-02,   nfev = 24,   norm = 0.076
log10(alpha) = 1.00,    chi2 = 1.197e+00,   S = -8.279e-02,   nfev = 33,   norm = 0.053
log10(alpha) = 0.00,    chi2 = 7.387e-01,   S = -1.655e-01,   nfev = 25,   norm = 0.040
log10(alpha) = -1.00,   chi2 = 4.813e-01,   S = -6.434e-01,   nfev = 49,   norm = 0.030
log10(alpha) = -2.00,   chi2 = 2.774e-01,   S = -4.180e+00,   nfev = 34,   norm = 0.020
log10(alpha) = -3.00,   chi2 = 1.689e-01,   S = -2.425e+01,   nfev = 70,   norm = 0.011
log10(alpha) = -4.00,   chi2 = 1.132e-01,   S = -1.098e+02,   nfev = 73,   norm = 0.004
log10(alpha) = -5.00,   chi2 = 8.941e-02,   S = -5.790e+02,   nfev = 62,   norm = 0.000
log10(alpha) = -6.00,   chi2 = 8.205e-02,   S = -1.574e+03,   nfev = 35,   norm = 0.002
Fit parameters [-0.85180161  5.67618507  3.72214498  0.65382993]
Optimal log alpha 1.4279704254349954
log10(alpha) = 1.43,    chi2 = 1.560e+00,   S = -7.192e-02,   nfev = 12,   norm = 0.061
_images/tutorial_30_22.png
_images/tutorial_30_23.png
Offdiagonal element 2 1
401 data points on real axis
80 data points on imaginary axis
26 significant singular values
Precomputation of coefficient matrices...
log10(alpha) = 10.00,   chi2 = 4.751e+04,   S = -3.436e-10,   nfev = 1,   norm = 0.000
log10(alpha) = 9.00,    chi2 = 4.738e+04,   S = -3.426e-08,   nfev = 2,   norm = 0.001
log10(alpha) = 8.00,    chi2 = 4.618e+04,   S = -3.323e-06,   nfev = 9,   norm = 0.012
log10(alpha) = 7.00,    chi2 = 3.663e+04,   S = -2.511e-04,   nfev = 12,   norm = 0.104
log10(alpha) = 6.00,    chi2 = 1.183e+04,   S = -4.974e-03,   nfev = 15,   norm = 0.403
log10(alpha) = 5.00,    chi2 = 1.842e+03,   S = -2.002e-02,   nfev = 18,   norm = 0.385
log10(alpha) = 4.00,    chi2 = 1.071e+02,   S = -4.088e-02,   nfev = 19,   norm = 0.204
log10(alpha) = 3.00,    chi2 = 1.381e+01,   S = -5.215e-02,   nfev = 19,   norm = 0.116
log10(alpha) = 2.00,    chi2 = 2.879e+00,   S = -6.659e-02,   nfev = 17,   norm = 0.070
log10(alpha) = 1.00,    chi2 = 1.428e+00,   S = -8.733e-02,   nfev = 40,   norm = 0.049
log10(alpha) = 0.00,    chi2 = 9.735e-01,   S = -1.707e-01,   nfev = 25,   norm = 0.037
log10(alpha) = -1.00,   chi2 = 6.780e-01,   S = -7.193e-01,   nfev = 39,   norm = 0.028
log10(alpha) = -2.00,   chi2 = 4.174e-01,   S = -5.514e+00,   nfev = 70,   norm = 0.017
log10(alpha) = -3.00,   chi2 = 2.640e-01,   S = -3.380e+01,   nfev = 48,   norm = 0.008
log10(alpha) = -4.00,   chi2 = 1.680e-01,   S = -1.817e+02,   nfev = 101,   norm = 0.002
log10(alpha) = -5.00,   chi2 = 1.218e-01,   S = -1.083e+03,   nfev = 91,   norm = -0.001
log10(alpha) = -6.00,   chi2 = 1.110e-01,   S = -2.586e+03,   nfev = 55,   norm = 0.000
Fit parameters [-0.65452172  5.58782077  3.88284786  0.69667046]
Optimal log alpha 1.7297495431654122
log10(alpha) = 1.73,    chi2 = 2.224e+00,   S = -7.102e-02,   nfev = 25,   norm = 0.063
_images/tutorial_30_25.png
_images/tutorial_30_26.png
[26]:
# now plot all matrix elements from both methods and compare
fig, ax = plt.subplots(ncols=3, nrows=3, figsize=(15, 10))
for i in [0, 1, 2]:
    for j in [0, 1, 2]:
        ax[i, j].plot(w, sw_full[:, i, j].imag, color='green', ls='-', label='direct')
        ax[i, j].plot(w, sw_rotated_back[:, i, j].imag, color='orange', ls='-', label='from projection')
        ax[i, j].set_xlim(-10., 10)
        ax[i, j].axhline(0., color='gray', lw=1)
ax[0, 0].legend()
plt.show()
_images/tutorial_31_0.png
[ ]: