yapanuwan
/
BergEC-jl
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

More testing for Gram-Schmidt

This commit is contained in:
Nuwan Yapa 2025-01-06 15:16:53 -05:00
parent f51c0f32e6
commit 4815366e36
2 changed files with 37 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
helper.jl
View File

@ -51,11 +51,9 @@ end
"c-orthonormalization via Gram-Schmidt" "c-orthonormalization via Gram-Schmidt"
function gram_schmidt(vecs, ws, precision=1e-10, verbose=false) function gram_schmidt(vecs, ws, precision=1e-10, verbose=false)
function proj(ui, vi) c_product(i, j) = sum(vecs[i] .* ws .* vecs[j])
num = sum(vecs[ui] .* ws .* vecs[vi]) proj(i, j) = (c_product(i, j) / c_product(i, i)) .* vecs[i] # component of vec[i] in vec[j]
den = sum(vecs[ui] .* ws .* vecs[ui])
return (num / den) .* vecs[ui]
end
mask = ones(Bool, size(vecs)) # vectors to retain mask = ones(Bool, size(vecs)) # vectors to retain
for i in eachindex(vecs) for i in eachindex(vecs)
for j in findall(mask[1 : (i - 1)]) for j in findall(mask[1 : (i - 1)])
@ -73,9 +71,12 @@ function gram_schmidt(vecs, ws, precision=1e-10, verbose=false)
end end
"Rank of a basis set (pre-normalized) under c-product" "Rank of a basis set (pre-normalized) under c-product"
function c_rank(vecs, ws, threshold=1e-8) c_rank(vecs, ws, threshold=1e-8) = count(abs.(c_singular_values(vecs, ws)) .> threshold)
"Singular values of a basis set (pre-normalized) under c-product"
function c_singular_values(vecs, ws)
basis = hcat(vecs...) basis = hcat(vecs...)
N = transpose(basis) * spdiagm(ws) * basis N = transpose(basis) * spdiagm(ws) * basis
abs_singular_values = eigvals(N) .|> abs .|> sqrt singular_values = eigvals(N) .|> complex .|> sqrt
return count(abs_singular_values .> threshold) return singular_values
end end

32
test/misc.jl
View File

@ -23,11 +23,7 @@ evals = eigvals(H_EC, N_EC)
println("Eigenvalues with GEVP:") println("Eigenvalues with GEVP:")
display(evals) display(evals)
push!(vecs, vecs[end]) # duplicate last basis vector to test orthogonalization ortho_basis = hcat(gram_schmidt(vecs, ws)...)
ortho_vecs = gram_schmidt(vecs, ws)
@assert length(ortho_vecs) == N "Basis dimensionality mismatch after Gram-Schmidt"
ortho_basis = hcat(ortho_vecs...)
N_ortho = transpose(ortho_basis) * ws_mat * ortho_basis N_ortho = transpose(ortho_basis) * ws_mat * ortho_basis
println("Norm matrix after Gram-Schmidt:") println("Norm matrix after Gram-Schmidt:")
display(round.(N_ortho; digits=2)) # should be ≈I display(round.(N_ortho; digits=2)) # should be ≈I
@ -41,4 +37,30 @@ display(evals_ortho)
@assert evals ≈ evals_ortho "Gram-Schmidt did not preserve the eigenvalues" @assert evals ≈ evals_ortho "Gram-Schmidt did not preserve the eigenvalues"
############
println("\nRepeat with a redundant basis\n")
println("Original dimensionality = $(length(vecs))")
for pow in 6:9
noise = rand(n) ./ 10^pow
new_vec = vecs[1] + noise
push!(vecs, new_vec)
end
println("Dimensionality before Gram-Schmidt = $(length(vecs))")
println("Absolute singular values before Gram-Schmidt:")
display(abs.(c_singular_values(vecs, ws)))
ortho_vecs = gram_schmidt(vecs, ws)
ortho_basis = hcat(ortho_vecs...)
println("Dimensionality after Gram-Schmidt = $(length(ortho_vecs))")
H_EC_ortho = transpose(ortho_basis) * ws_mat * H * ws_mat * ortho_basis
evals_ortho = eigvals(H_EC_ortho)
println("Eigenvalues after Gram-Schmidt:")
display(evals_ortho)
###################### ######################