diff --git a/helper.jl b/helper.jl
index c921493..fc7a71b 100644
--- a/helper.jl
+++ b/helper.jl
@@ -49,7 +49,8 @@ function exportCSV(file::String, data, labels=nothing)
     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)
     c_product(i, j) = sum(vecs[i] .* ws .* vecs[j])
     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)
     verbose && println("Target dimensionality = $target_dim")
 
-    for i in eachindex(vecs)
-        for j in (i + 1):length(vecs)
+    selected_vecs_i = Integer[]
+
+    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)
         end
     end
     
-    abs_norms = abs.(norm.(eachindex(vecs)))
-    verbose && println("Absolute norms = $(round.(abs_norms; sigdigits=1))")
-    selected_vecs_i = partialsortperm(abs_norms, 1:target_dim; rev=true)
+    verbose && println("Absolute norms of selected vectors = $(round.(selected_vecs_i .|> norm .|> abs; sigdigits=1))")
+    verbose && println("Absolute norms of dropped vectors = $(round.(setdiff(eachindex(vecs), selected_vecs_i) .|> norm .|> abs; sigdigits=1))")
 
+    # final normalization
     for i in selected_vecs_i
         vecs[i] ./= sqrt(norm(i))
     end
diff --git a/test/misc.jl b/test/misc.jl
index 7ff7d12..9329e2e 100644
--- a/test/misc.jl
+++ b/test/misc.jl
@@ -50,9 +50,6 @@ 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; verbose=true)
 ortho_basis = hcat(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:")
 display(evals_ortho)
 
+@assert isapprox(evals, evals_ortho; atol=1e-3) "Gram-Schmidt did not approximately preserve the eigenvalues"
+
 ######################
\ No newline at end of file