diff --git a/EC.jl b/EC.jl index 01797ed..61cadb6 100644 --- a/EC.jl +++ b/EC.jl @@ -23,9 +23,11 @@ end If a list is provided for ref_eval, they are used as reference values for picking the closest eigenvalues at each sampling point. If a single number is provided for ref_eval, it is used as a reference for the first point, and the previous eigenvalue is used as the reference for each successive point. If pseudo_inv_rtol > 0, the GEVP is avoided using Moore-Penrose psuedoinverse, using this value as the relative tolerance for dropping redundant vectors. -If orthonormalize_threshold > 0, Gram-Schmidt orthonormalization is performed, overriding Moore-Penrose, using this value as the threshold for dropping redundant vectors. +If orthonormalize_threshold > 0, Gram-Schmidt orthonormalization is performed, using this value as the threshold for dropping redundant vectors. If both are = 0, GEVP is solved instead." -function train!(EC::affine_EC, c_vals; ref_eval=-10.0, pseudo_inv_rtol=1e-6, orthonormalize_threshold=0, CAEC=false, verbose=true, tol=1e-5) +function train!(EC::affine_EC, c_vals; ref_eval=-10.0, pseudo_inv_rtol=1e-6, gram_schmidt_threshold=0, CAEC=false, verbose=true, tol=1e-5) + @assert pseudo_inv_rtol == 0 || gram_schmidt_threshold == 0 "Cannot perform both psuedoinverse and Gram-Schmidt" + training_vecs = Vector{ComplexF64}[] for c in c_vals @@ -50,23 +52,21 @@ function train!(EC::affine_EC, c_vals; ref_eval=-10.0, pseudo_inv_rtol=1e-6, ort CAEC && append!(training_vecs, conj.(training_vecs)) - if orthonormalize_threshold ≠ 0 - training_vecs = gram_schmidt!(training_vecs, EC.weights, orthonormalize_threshold; verbose=verbose) + if gram_schmidt_threshold ≠ 0 + training_vecs = gram_schmidt!(training_vecs, EC.weights, gram_schmidt_threshold; verbose=verbose) end EC_basis = hcat(training_vecs...) weights_mat = spdiagm(EC.weights) EC.H0_EC = transpose(EC_basis) * weights_mat * EC.H0 * EC_basis EC.H1_EC = transpose(EC_basis) * weights_mat * EC.H1 * EC_basis + EC.N_EC = transpose(EC_basis) * weights_mat * EC_basis - if orthonormalize_threshold == 0 - EC.N_EC = transpose(EC_basis) * weights_mat * EC_basis - if pseudo_inv_rtol > 0 - inv_N_EC = pinv(EC.N_EC; rtol=pseudo_inv_rtol) - EC.H0_EC = inv_N_EC * EC.H0_EC - EC.H1_EC = inv_N_EC * EC.H1_EC - EC.N_EC = nothing - end + if pseudo_inv_rtol > 0 + inv_N_EC = pinv(EC.N_EC; rtol=pseudo_inv_rtol) + EC.H0_EC = inv_N_EC * EC.H0_EC + EC.H1_EC = inv_N_EC * EC.H1_EC + EC.N_EC = nothing end EC.trained = true