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

Improved stability for Gram-Schmidt

This commit is contained in:
Nuwan Yapa 2025-01-07 18:15:50 -05:00
parent 22c9c1eaf1
commit 99739a011c
2 changed files with 13 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
helper.jl
View File

@ -49,7 +49,8 @@ function exportCSV(file::String, data, labels=nothing)
end end
end end
"In-place c-orthonormalization via (modified) Gram-Schmidt. Most significant vectors determined by the singular values (compared to the threshold) are returned." "In-place c-orthogonalization via (modified) Gram-Schmidt. Only significant vectors are returned (c-normalized).
The number of significant vectors to return are determined by the original singular values (compared to the threshold) "
function gram_schmidt!(vecs, ws, threshold=1e-5; verbose=false) function gram_schmidt!(vecs, ws, threshold=1e-5; verbose=false)
c_product(i, j) = sum(vecs[i] .* ws .* vecs[j]) c_product(i, j) = sum(vecs[i] .* ws .* vecs[j])
norm(i) = c_product(i, i) norm(i) = c_product(i, i)
@ -64,16 +65,20 @@ function gram_schmidt!(vecs, ws, threshold=1e-5; verbose=false)
target_dim = c_rank(vecs, ws, threshold) target_dim = c_rank(vecs, ws, threshold)
verbose && println("Target dimensionality = $target_dim") verbose && println("Target dimensionality = $target_dim")
for i in eachindex(vecs) selected_vecs_i = Integer[]
for j in (i + 1):length(vecs)
while length(selected_vecs_i) < target_dim
i = argmax(i -> abs(norm(i)), setdiff(eachindex(vecs), selected_vecs_i)) # find the largest vector from the remaining
push!(selected_vecs_i, i)
for j in setdiff(eachindex(vecs), selected_vecs_i)
vecs[j] .-= proj(i, j) vecs[j] .-= proj(i, j)
end end
end end
abs_norms = abs.(norm.(eachindex(vecs))) verbose && println("Absolute norms of selected vectors = $(round.(selected_vecs_i .|> norm .|> abs; sigdigits=1))")
verbose && println("Absolute norms = $(round.(abs_norms; sigdigits=1))") verbose && println("Absolute norms of dropped vectors = $(round.(setdiff(eachindex(vecs), selected_vecs_i) .|> norm .|> abs; sigdigits=1))")
selected_vecs_i = partialsortperm(abs_norms, 1:target_dim; rev=true)
# final normalization
for i in selected_vecs_i for i in selected_vecs_i
vecs[i] ./= sqrt(norm(i)) vecs[i] ./= sqrt(norm(i))
end end

5
test/misc.jl
View File

@ -50,9 +50,6 @@ end
println("Dimensionality before Gram-Schmidt = $(length(vecs))") 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; verbose=true) ortho_vecs = gram_schmidt!(vecs, ws; verbose=true)
ortho_basis = hcat(ortho_vecs...) ortho_basis = hcat(ortho_vecs...)
println("Dimensionality after Gram-Schmidt = $(length(ortho_vecs))") println("Dimensionality after Gram-Schmidt = $(length(ortho_vecs))")
@ -63,4 +60,6 @@ evals_ortho = eigvals(H_EC_ortho)
println("Eigenvalues after Gram-Schmidt:") println("Eigenvalues after Gram-Schmidt:")
display(evals_ortho) display(evals_ortho)
@assert isapprox(evals, evals_ortho; atol=1e-3) "Gram-Schmidt did not approximately preserve the eigenvalues"
###################### ######################