From 4815366e36bea3bc5f192ed99c6af2b591c82423 Mon Sep 17 00:00:00 2001
From: Nuwan Yapa <yapanuwan@gmail.com>
Date: Mon, 6 Jan 2025 15:16:53 -0500
Subject: [PATCH] More testing for Gram-Schmidt

---
 helper.jl    | 19 ++++++++++---------
 test/misc.jl | 32 +++++++++++++++++++++++++++-----
 2 files changed, 37 insertions(+), 14 deletions(-)

diff --git a/helper.jl b/helper.jl
index da9f150..8dd46a9 100644
--- a/helper.jl
+++ b/helper.jl
@@ -51,11 +51,9 @@ end
 
 "c-orthonormalization via Gram-Schmidt"
 function gram_schmidt(vecs, ws, precision=1e-10, verbose=false)
-    function proj(ui, vi)
-        num = sum(vecs[ui] .* ws .* vecs[vi])
-        den = sum(vecs[ui] .* ws .* vecs[ui])
-        return (num / den) .* vecs[ui]
-    end
+    c_product(i, j) = sum(vecs[i] .* ws .* vecs[j])
+    proj(i, j) = (c_product(i, j) / c_product(i, i)) .* vecs[i] # component of vec[i] in vec[j]
+
     mask = ones(Bool, size(vecs)) # vectors to retain
     for i in eachindex(vecs)
         for j in findall(mask[1 : (i - 1)])
@@ -73,9 +71,12 @@ function gram_schmidt(vecs, ws, precision=1e-10, verbose=false)
 end
 
 "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...)
     N = transpose(basis) * spdiagm(ws) * basis
-    abs_singular_values = eigvals(N) .|> abs .|> sqrt
-    return count(abs_singular_values .> threshold)
-end
\ No newline at end of file
+    singular_values = eigvals(N) .|> complex .|> sqrt
+    return singular_values
+end
diff --git a/test/misc.jl b/test/misc.jl
index c7126d8..19d47db 100644
--- a/test/misc.jl
+++ b/test/misc.jl
@@ -23,11 +23,7 @@ evals = eigvals(H_EC, N_EC)
 println("Eigenvalues with GEVP:")
 display(evals)
 
-push!(vecs, vecs[end]) # duplicate last basis vector to test orthogonalization
-ortho_vecs = gram_schmidt(vecs, ws)
-@assert length(ortho_vecs) == N "Basis dimensionality mismatch after Gram-Schmidt"
-ortho_basis = hcat(ortho_vecs...)
-
+ortho_basis = hcat(gram_schmidt(vecs, ws)...)
 N_ortho = transpose(ortho_basis) * ws_mat * ortho_basis
 println("Norm matrix after Gram-Schmidt:")
 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"
 
+############
+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)
+
 ######################
\ No newline at end of file