Matrix valued continuation & Hardy optimization
[1]:
import sys, os
import numpy as np
import scipy
from triqs.plot.mpl_interface import plt,oplot
from h5 import HDFArchive
from triqs.atom_diag import *
from triqs.gf import *
from triqs.operators import c, c_dag, n, dagger
from itertools import product
from triqs.lattice.tight_binding import TBLattice
from triqs.lattice.utils import k_space_path
from triqs_tprf.lattice import lattice_dyson_g0_wk, lattice_dyson_g0_fk
from triqs_Nevanlinna import Solver
from triqs_maxent import *
import seaborn as sns
import scienceplots
plt.style.use(['science','notebook'])
# sns.set_palette('muted')
Warning: could not identify MPI environment!
Starting serial run at: 2023-07-19 16:58:29.118244
Setup simple two orbital 2D Hubbard model on square lattice
[2]:
n_orb = 2
loc = np.zeros((n_orb,n_orb))
for i in range(n_orb):
for j in range(n_orb):
if i != 0 and i==j:
loc[i,j] = 0.2
if j > i or j < i:
loc[i,j] = -0.4
#nearest neighbor hopping
if n_orb == 2:
t = np.array([[-1.,0.1],[0.1,-0.4]])
else:
t = -1.0 * np.eye(n_orb)
tp = 0.1 * np.eye(n_orb) #next nearest neighbor hopping
hop= { (0,0) : loc,
(1,0) : t,
(-1,0) : t,
(0,1) : t,
(0,-1) : t,
(1,1) : tp,
(-1,-1): tp,
(1,-1) : tp,
(-1,1) : tp}
TB = TBLattice(units = [(1, 0, 0) , (0, 1, 0)], hoppings = hop, orbital_positions= [(0., 0., 0.)]*n_orb)
plot dispersion along high-symmetry lines
[3]:
n_pts = 101
G = np.array([ 0.00, 0.00, 0.00])
M = np.array([0.5, 0.5, 0.0])
X = np.array([0.5, 0.0, 0.0])
R = np.array([0.5, 0.5, 0.5])
paths = [(R, G), (G, X), (X, M), (M, G)]
kvecs, k = k_space_path(paths, num=n_pts, bz=TB.bz)
e_mat = TB.fourier(kvecs)
e_val = TB.dispersion(kvecs)
eps_k = {'k': k, 'K': np.concatenate(([0],k[n_pts-1::n_pts])), 'eval': e_val, 'emat' : e_mat}
[4]:
fig, ax = plt.subplots(1,1, figsize=(8,4), dpi=100, squeeze=False)
ax = ax.reshape(-1)
for band in range(eps_k['eval'].shape[1]):
ax[0].plot(eps_k['k'], eps_k['eval'][:,band].real, color='C0', zorder=2.5)
ax[0].axhline(y=0,zorder=2,color='gray',alpha=0.5,ls='--')
ax[0].set_xticks(eps_k['K'])
ax[0].set_xticklabels([r'R', '$\Gamma$', 'X', 'M', r'$\Gamma$'])
ax[0].set_xlim([eps_k['K'].min(), eps_k['K'].max()])
ax[0].set_ylabel(r'$\omega$ (eV)')
plt.show()
Setup lattice Green’s function and calculate \(G_{loc}\) on Matsubara axis
[5]:
k_grid = (100,100,1)
k_mesh = TB.get_kmesh(n_k = k_grid)
e_k = TB.fourier(k_mesh)
eps_k = TB.dispersion(k_mesh)
mu = 0.
beta = 100
n_iw = 1000
iw_mesh = MeshImFreq(beta=beta, S='Fermion', n_max=n_iw)
G0_iwk = lattice_dyson_g0_wk(mu=mu, e_k=e_k, mesh=iw_mesh)
calculate local Green’s function
[6]:
G_Γ_iw = G0_iwk[:, Idx(0,0,0)]
G_X_iw = G0_iwk[:, Idx(k_grid[1]-1,0,0)]
G_loc_iw = Gf(mesh=iw_mesh, target_shape=(n_orb,n_orb))
iw_arr = np.array(list(iw_mesh.values()))
for k in G0_iwk.mesh[1]:
G_loc_iw[:] += G0_iwk[:,k]
G_loc_iw = G_loc_iw / len(k_mesh)
[7]:
fig, ax = plt.subplots(n_orb, n_orb, figsize=(6*n_orb,3*n_orb), dpi=100, squeeze=False,sharex=True)
shp = G_Γ_iw.target_shape
for i in range(n_orb):
for j in range(n_orb):
ax[i,j].oplot(G_loc_iw[i,j].imag, label=f'[{i},{j}] sum')
ax[i,j].oplot(G_Γ_iw[i,j].imag, label=f'[{i},{j}]@$\Gamma$')
ax[i,j].oplot(G_X_iw[i,j].imag, label=f'[{i},{j}]@X')
ax[i,j].set_xlim(0,50)
if i == shp[0]-1:
ax[i,j].set_xlabel(r'$i\omega_n$')
else:
ax[i,j].set_xlabel(r'')
if j == 0:
ax[i,j].set_ylabel(r'Im $G (i\omega_n)$ (eV)')
else:
ax[i,j].set_ylabel(r'')
plt.tight_layout(pad=0.4, w_pad=0.1, h_pad=0.4)
plt.show()
The local Green’s function is off-diagonal in orbital space.
Calculate reference data directly on real frequency axis
[8]:
n_w = 1001
delta = 0.1
w_window = (-7, 7)
w_mesh = MeshReFreq(window=w_window, n_w=n_w)
G0_wk = lattice_dyson_g0_fk(mu=mu, e_k=e_k, mesh=w_mesh, delta=delta)
G_Γ_w = G0_wk[:, Idx(0,0,0)]
G_X_w = G0_wk[:, Idx(k_grid[1]-1,0,0)]
G_loc_w = Gf(mesh=w_mesh, target_shape=(n_orb,n_orb))
w_arr = np.array(list(w_mesh.values()))
for k in G0_wk.mesh[1]:
G_loc_w[:] += G0_wk[:,k]
G_loc_w = G_loc_w / len(k_mesh)
Analytic continuation with Pade, Nevanlinna, and MaxEnt
Pade
[9]:
G_loc_w_P = G_loc_w.copy()
G_loc_w_P.set_from_pade(G_loc_iw, n_points=n_iw, freq_offset=delta)
Nevanlinna & Caratheodory
[10]:
# setup Nevanlinna kernel solver
nevan_mesh = MeshReFreq(window=w_window, n_w=n_w//2)
solver = Solver(kernel='NEVANLINNA')
solver.solve(G_loc_iw)
G_loc_w_N = solver.evaluate(nevan_mesh, delta)
This is Nevanlinna analytical continuation. All off-diagonal elements will be ignored.
[11]:
# setup matrix valued Nevanlinna kernel solver
solver = Solver(kernel='CARATHEODORY')
solver.solve(G_loc_iw)
G_loc_w_CN = solver.evaluate(nevan_mesh, delta)
run MaxEnt
[12]:
tm = PoormanMaxEnt(use_hermiticity=True)
tm.set_G_iw(G_loc_iw)
tm.set_error(1.e-5)
tm.omega = HyperbolicOmegaMesh(omega_min=w_window[0], omega_max=w_window[1], n_points=200)
tm.alpha_mesh = LogAlphaMesh(alpha_min=1e-4, alpha_max=1e2, n_points=30)
result = tm.run()
k_diag = tm.maxent_diagonal.K
k_off_diag = tm.maxent_offdiagonal.K
print(result.analyzer_results[0][0]['LineFitAnalyzer']['alpha_index'])
print(result.analyzer_results[0][1]['LineFitAnalyzer']['alpha_index'])
print(result.analyzer_results[1][1]['LineFitAnalyzer']['alpha_index'])
w_mesh_arr = np.linspace(w_mesh.w_min, w_mesh.w_max, len(w_mesh))
G_loc_w_ME = np.zeros((n_orb,n_orb,len(w_mesh_arr)))
for i in range(n_orb):
for j in range(n_orb):
G_loc_w_ME[i,j,:] = np.interp(w_mesh_arr, np.array(result.omega), result.get_A_out('LineFitAnalyzer')[i,j])
Calculating diagonal elements.
Calling MaxEnt for element 0 0
2023-07-19 17:02:03.351843
MaxEnt run
TRIQS application maxent
Copyright(C) 2018 Gernot J. Kraberger
Copyright(C) 2018 Simons Foundation
Authors: Gernot J. Kraberger and Manuel Zingl
This program comes with ABSOLUTELY NO WARRANTY.
This is free software, and you are welcome to redistributeit under certain conditions; see file LICENSE.
Please cite this code and the appropriate original papers (see documentation).
Minimal chi2: 0.05023027078433784
scaling alpha by a factor 6001 (number of data points)
alpha[ 0] = 6.00100000e+05, chi2 = 1.68505715e+04, n_iter= 439
alpha[ 1] = 3.72672267e+05, chi2 = 9.36636420e+03, n_iter= 14
alpha[ 2] = 2.31435791e+05, chi2 = 5.17253363e+03, n_iter= 15
alpha[ 3] = 1.43725547e+05, chi2 = 2.85449082e+03, n_iter= 17
alpha[ 4] = 8.92560000e+04, chi2 = 1.57641824e+03, n_iter= 19
alpha[ 5] = 5.54294881e+04, chi2 = 8.66716052e+02, n_iter= 22
alpha[ 6] = 3.44226512e+04, chi2 = 4.72864090e+02, n_iter= 24
alpha[ 7] = 2.13770496e+04, chi2 = 2.57524710e+02, n_iter= 26
alpha[ 8] = 1.32755100e+04, chi2 = 1.41518184e+02, n_iter= 29
alpha[ 9] = 8.24431660e+03, chi2 = 7.89959548e+01, n_iter= 32
alpha[10] = 5.11986028e+03, chi2 = 4.48827884e+01, n_iter= 37
alpha[11] = 3.17951998e+03, chi2 = 2.60389633e+01, n_iter= 43
alpha[12] = 1.97453577e+03, chi2 = 1.54578466e+01, n_iter= 49
alpha[13] = 1.22622017e+03, chi2 = 9.31711843e+00, n_iter= 58
alpha[14] = 7.61503498e+02, chi2 = 5.61319442e+00, n_iter= 70
alpha[15] = 4.72906574e+02, chi2 = 3.34876762e+00, n_iter= 86
alpha[16] = 2.93682994e+02, chi2 = 1.98993106e+00, n_iter= 103
alpha[17] = 1.82382115e+02, chi2 = 1.19430780e+00, n_iter= 117
alpha[18] = 1.13262383e+02, chi2 = 7.30419326e-01, n_iter= 142
alpha[19] = 7.03378589e+01, chi2 = 4.56776182e-01, n_iter= 176
alpha[20] = 4.36810020e+01, chi2 = 2.95476130e-01, n_iter= 207
alpha[21] = 2.71266423e+01, chi2 = 2.02760019e-01, n_iter= 254
alpha[22] = 1.68461044e+01, chi2 = 1.51063747e-01, n_iter= 252
alpha[23] = 1.04617163e+01, chi2 = 1.22535756e-01, n_iter= 304
alpha[24] = 6.49690304e+00, chi2 = 1.06609542e-01, n_iter= 351
alpha[25] = 4.03468686e+00, chi2 = 9.76252631e-02, n_iter= 453
alpha[26] = 2.50560889e+00, chi2 = 9.25768858e-02, n_iter= 498
alpha[27] = 1.55602557e+00, chi2 = 8.97376449e-02, n_iter= 628
alpha[28] = 9.66318243e-01, chi2 = 8.80835577e-02, n_iter= 725
alpha[29] = 6.00100000e-01, chi2 = 8.70484247e-02, n_iter= 909
MaxEnt loop finished in 0:01:07.877541
Calling MaxEnt for element 1 1
Minimal chi2: 0.00854902316184783
scaling alpha by a factor 6001 (number of data points)
alpha[ 0] = 6.00100000e+05, chi2 = 4.21558814e+03, n_iter= 701
alpha[ 1] = 3.72672267e+05, chi2 = 2.04530448e+03, n_iter= 12
alpha[ 2] = 2.31435791e+05, chi2 = 1.03675545e+03, n_iter= 12
alpha[ 3] = 1.43725547e+05, chi2 = 5.49084178e+02, n_iter= 13
alpha[ 4] = 8.92560000e+04, chi2 = 3.05459136e+02, n_iter= 14
alpha[ 5] = 5.54294881e+04, chi2 = 1.80819227e+02, n_iter= 14
alpha[ 6] = 3.44226512e+04, chi2 = 1.14469855e+02, n_iter= 16
alpha[ 7] = 2.13770496e+04, chi2 = 7.64565053e+01, n_iter= 19
alpha[ 8] = 1.32755100e+04, chi2 = 5.23791300e+01, n_iter= 23
alpha[ 9] = 8.24431660e+03, chi2 = 3.56553240e+01, n_iter= 30
alpha[10] = 5.11986028e+03, chi2 = 2.35431227e+01, n_iter= 38
alpha[11] = 3.17951998e+03, chi2 = 1.49362893e+01, n_iter= 46
alpha[12] = 1.97453577e+03, chi2 = 9.16221686e+00, n_iter= 54
alpha[13] = 1.22622017e+03, chi2 = 5.53586336e+00, n_iter= 59
alpha[14] = 7.61503498e+02, chi2 = 3.35887712e+00, n_iter= 92
alpha[15] = 4.72906574e+02, chi2 = 2.05530424e+00, n_iter= 96
alpha[16] = 2.93682994e+02, chi2 = 1.24861624e+00, n_iter= 109
alpha[17] = 1.82382115e+02, chi2 = 7.35949338e-01, n_iter= 166
alpha[18] = 1.13262383e+02, chi2 = 4.15297466e-01, n_iter= 220
alpha[19] = 7.03378589e+01, chi2 = 2.25869801e-01, n_iter= 182
alpha[20] = 4.36810020e+01, chi2 = 1.21854134e-01, n_iter= 256
alpha[21] = 2.71266423e+01, chi2 = 6.82062067e-02, n_iter= 251
alpha[22] = 1.68461044e+01, chi2 = 4.14814836e-02, n_iter= 265
alpha[23] = 1.04617163e+01, chi2 = 2.82331641e-02, n_iter= 260
alpha[24] = 6.49690304e+00, chi2 = 2.15754451e-02, n_iter= 237
alpha[25] = 4.03468686e+00, chi2 = 1.81381916e-02, n_iter= 273
alpha[26] = 2.50560889e+00, chi2 = 1.62620878e-02, n_iter= 405
alpha[27] = 1.55602557e+00, chi2 = 1.51426742e-02, n_iter= 571
alpha[28] = 9.66318243e-01, chi2 = 1.44205067e-02, n_iter= 1000!
alpha[29] = 6.00100000e-01, chi2 = 1.39449518e-02, n_iter= 1000!
! ... The minimizer did not converge. Results might be wrong.
MaxEnt loop finished in 0:01:13.534457
Calculating off-diagonal elements using default model from diagonal solution
Calling MaxEnt for element 0 1
2023-07-19 17:04:31.955070
MaxEnt run
TRIQS application maxent
Copyright(C) 2018 Gernot J. Kraberger
Copyright(C) 2018 Simons Foundation
Authors: Gernot J. Kraberger and Manuel Zingl
This program comes with ABSOLUTELY NO WARRANTY.
This is free software, and you are welcome to redistributeit under certain conditions; see file LICENSE.
Please cite this code and the appropriate original papers (see documentation).
Minimal chi2: 0.04209151703611798
scaling alpha by a factor 6001 (number of data points)
alpha[ 0] = 6.00100000e+05, chi2 = 3.12438015e+02, n_iter= 5
alpha[ 1] = 3.72672267e+05, chi2 = 1.73027132e+02, n_iter= 4
alpha[ 2] = 2.31435791e+05, chi2 = 9.53652431e+01, n_iter= 4
alpha[ 3] = 1.43725547e+05, chi2 = 5.45231975e+01, n_iter= 4
alpha[ 4] = 8.92560000e+04, chi2 = 3.33402096e+01, n_iter= 4
alpha[ 5] = 5.54294881e+04, chi2 = 2.17268186e+01, n_iter= 4
alpha[ 6] = 3.44226512e+04, chi2 = 1.45184692e+01, n_iter= 5
alpha[ 7] = 2.13770496e+04, chi2 = 9.56991411e+00, n_iter= 7
alpha[ 8] = 1.32755100e+04, chi2 = 6.26274467e+00, n_iter= 7
alpha[ 9] = 8.24431660e+03, chi2 = 4.22753762e+00, n_iter= 10
alpha[10] = 5.11986028e+03, chi2 = 2.98557094e+00, n_iter= 11
alpha[11] = 3.17951998e+03, chi2 = 2.15156745e+00, n_iter= 12
alpha[12] = 1.97453577e+03, chi2 = 1.52273579e+00, n_iter= 15
alpha[13] = 1.22622017e+03, chi2 = 1.03073788e+00, n_iter= 19
alpha[14] = 7.61503498e+02, chi2 = 6.64792988e-01, n_iter= 22
alpha[15] = 4.72906574e+02, chi2 = 4.16473335e-01, n_iter= 26
alpha[16] = 2.93682994e+02, chi2 = 2.62424657e-01, n_iter= 31
alpha[17] = 1.82382115e+02, chi2 = 1.72797964e-01, n_iter= 36
alpha[18] = 1.13262383e+02, chi2 = 1.22405765e-01, n_iter= 43
alpha[19] = 7.03378589e+01, chi2 = 9.44432582e-02, n_iter= 48
alpha[20] = 4.36810020e+01, chi2 = 7.90202600e-02, n_iter= 53
alpha[21] = 2.71266423e+01, chi2 = 7.06342428e-02, n_iter= 59
alpha[22] = 1.68461044e+01, chi2 = 6.62091182e-02, n_iter= 66
alpha[23] = 1.04617163e+01, chi2 = 6.38902337e-02, n_iter= 73
alpha[24] = 6.49690304e+00, chi2 = 6.25851689e-02, n_iter= 74
alpha[25] = 4.03468686e+00, chi2 = 6.17447850e-02, n_iter= 89
alpha[26] = 2.50560889e+00, chi2 = 6.10999316e-02, n_iter= 164
alpha[27] = 1.55602557e+00, chi2 = 6.05099728e-02, n_iter= 206
alpha[28] = 9.66318243e-01, chi2 = 5.99675991e-02, n_iter= 237
alpha[29] = 6.00100000e-01, chi2 = 5.93482651e-02, n_iter= 297
MaxEnt loop finished in 0:00:21.167255
Element 1 0 not calculated, can be determined from hermiticity
22
19
24
Final plot for comparison
[13]:
fig, ax = plt.subplots(n_orb, n_orb, figsize=(8*n_orb,4*n_orb), dpi=100, squeeze=False,sharex=True)
shp = G_loc_w.target_shape
for i in range(n_orb):
for j in range(n_orb):
# plotting results
ax[i,j].oplot(-1*G_loc_w[i,j].imag, lw=5, label=f'[{i},{j}] ref')
ax[i,j].oplot(-1*G_loc_w_N[i,j].imag, label=f'[{i},{j}] Nev')
ax[i,j].oplot(-1*G_loc_w_CN[i,j].imag, label=f'[{i},{j}] Car')
ax[i,j].oplot(-1*G_loc_w_P[i,j].imag, label=f'[{i},{j}] Pade')
ax[i,j].plot(w_mesh_arr,np.pi*G_loc_w_ME[i,j], label=f'[{i},{j}] ME')
ax[i,j].set_xlim(w_window)
ax[i,j].legend()
if i == shp[0]-1:
ax[i,j].set_xlabel(r'$\omega$')
else:
ax[i,j].set_xlabel(r'')
if j == 0:
ax[i,j].set_ylabel(r'Im $G (\omega)$ (eV)')
else:
ax[i,j].set_ylabel(r'')
plt.tight_layout(pad=0.4, w_pad=0.1, h_pad=0.4)
plt.show()
fig, ax = plt.subplots(n_orb, n_orb, figsize=(8*n_orb,4*n_orb), dpi=100, squeeze=False,sharex=True)
for i in range(n_orb):
for j in range(n_orb):
# plotting results
ax[i,j].oplot(G_loc_w[i,j].real, lw=5, label=f'[{i},{j}] ref')
ax[i,j].oplot(G_loc_w_N[i,j].real, label=f'[{i},{j}] Nev')
ax[i,j].oplot(G_loc_w_CN[i,j].real, label=f'[{i},{j}] Car')
ax[i,j].oplot(G_loc_w_P[i,j].real, label=f'[{i},{j}] Pade')
ax[i,j].set_xlim(w_window)
if i == shp[0]-1:
ax[i,j].set_xlabel(r'$\omega$')
else:
ax[i,j].set_xlabel(r'')
if j == 0:
ax[i,j].set_ylabel(r'Re $G (\omega)$ (eV)')
else:
ax[i,j].set_ylabel(r'')
plt.tight_layout(pad=0.4, w_pad=0.1, h_pad=0.4)
plt.show()
Hardy optimization
Both Nevanlinna and Caratheodory continuation work well for finite systems. However, for infinite systems / continuous spectra the continuation will return a highly oscillating function. To address this issue an additional optimization of the continuation can be performed.
In the Nevanlinna package such a Hardy function optimization is implemented. As described in PRL 126 056402 (2021) the continued fraction form of the Nevanlinna factorization has a freedom in choosing \(\theta_{M+1}\), the last element in the factorization:
\[G(z) = - h^{-1}\left[\theta(z)\right] \nonumber
\theta(z) = \frac{a(z) \theta_{M+1}(z) +b(z)}{c(z)\theta_{M+1}(z) +d(z)},\]
where \(a(z)\), \(b(z)\), \(c(z)\) and \(d(z)\) are Nevanlinna factorization coefficients, and \(h^{-1}\) is the inverse Mobius transform.
To perform the optimization the user has to define target function for the optimization problem, that takes a retarded Green’s function as an input. The optimization is then performed using SciPy.
First we repeat the continuation using the normal Nevanlinna kernel of the [0,0] orbital of the Green’s function from above.
[14]:
# setup Nevanlinna kernel solver
nevan_mesh = MeshReFreq(window=w_window, n_w=n_w//2)
solver = Solver(kernel='NEVANLINNA')
solver.solve(G_loc_iw[0:1,0:1])
G_w_non_opt = solver.evaluate(nevan_mesh, delta)
This is Nevanlinna analytical continuation. All off-diagonal elements will be ignored.
Now we perform the optimization, by defining first a target function that optimizes the smoothness of the resulting Green’s function. We then perform the optimization for a few eta values:
[21]:
nevan_mesh_arr = np.array([x for x in nevan_mesh.values()])
Gw = G_loc_iw[0:1,0:1]
eta_list = [0.4, 0.2, 0.1]
def SmoothnessOptimization (Gw):
lmbda = 1e-4
diff = 0.0
# loop over all orbitals
for i in range(Gw.data.shape[1]):
# extract spectral function
Aw = -Gw.data[:,i,i].imag/np.pi
# compute estimate for second derivative
Aw2 = (Aw[:-2] - 2*Aw[1:-1] + Aw[2:])/(Gw.mesh.delta**2)
# evaluate normalization criteria weight
diff += np.abs(1 - scipy.integrate.simpson(Aw, nevan_mesh_arr))**2
# evaluate smoothness criteria weight
diff += lmbda * scipy.integrate.simpson(Aw2**2, nevan_mesh_arr[1:-1])
return diff
# perform Hardy optimization
G_w_opt_res = []
for eta in eta_list:
print(f'eta={eta}')
G_w_opt_res.append(solver.optimize(nevan_mesh, eta, SmoothnessOptimization, gtol=1e-4, maxiter=1000, verbose=True))
print('-------')
eta=0.4
Optimization terminated successfully.
optimizer performed 0 functions evaluations.
-------
eta=0.2
Optimization terminated successfully.
optimizer performed 99 functions evaluations.
-------
eta=0.1
Optimization terminated successfully.
optimizer performed 293 functions evaluations.
-------
Now we plot all results in one plot:
[22]:
fig, ax = plt.subplots(1,1, figsize=(8,4), dpi=150, squeeze=False)
ax = ax.reshape(-1)
ax[0].oplot(-1*G_loc_w[0,0].imag, lw=5, label=f'ref')
ax[0].oplot(-1*G_w_non_opt.imag, label=f'Nev')
for eta, G_w_opt in zip(eta_list, G_w_opt_res):
ax[0].oplot(-1*G_w_opt.imag, label=fr'$\eta=${eta}')
ax[0].set_xlim(w_window)
ax[0].set_xlabel(r'$\omega$')
ax[0].set_ylabel(r'Im $G (\omega)$ (eV)')
plt.legend()
plt.tight_layout(pad=0.4, w_pad=0.1, h_pad=0.4)
plt.show()
It becomes clear that the larger eta values will introduce more smoothness. However, this goes hand in hand with a lost in features. An optimal eta has to be choosen carefully.
[ ]: