From 4c76130dc18b1e3431d9d640c013df049537377a Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Tue, 14 Oct 2025 11:42:34 -0400
Subject: [PATCH 01/10] remove hermitian checks from yalapack wrappers

---
 src/yalapack.jl | 22 +---------------------
 1 file changed, 1 insertion(+), 21 deletions(-)

diff --git a/src/yalapack.jl b/src/yalapack.jl
index e5e9f35c..96a6240a 100644
--- a/src/yalapack.jl
+++ b/src/yalapack.jl
@@ -10,7 +10,7 @@ module YALAPACK # Yet another lapack wrapper
 
 using LinearAlgebra: BlasFloat, BlasReal, BlasComplex, BlasInt, Char, LAPACK,
     LAPACKException, SingularException, PosDefException, checksquare, chkstride1,
-    require_one_based_indexing, triu!, issymmetric, ishermitian, isposdef, adjoint!
+    require_one_based_indexing, triu!, isposdef, adjoint!
 
 using LinearAlgebra.BLAS: @blasfunc, libblastrampoline
 using LinearAlgebra.LAPACK: chkfinite, chktrans, chkside, chkuplofinite, chklapackerror
@@ -989,11 +989,6 @@ for (heev, heevx, heevr, heevd, hegvd, elty, relty) in
             require_one_based_indexing(A, V, W)
             chkstride1(A, V, W)
             n = checksquare(A)
-            if $elty <: Real
-                issymmetric(A) || throw(ArgumentError("A must be symmetric"))
-            else
-                ishermitian(A) || throw(ArgumentError("A must be Hermitian"))
-            end
             chkuplofinite(A, uplo)
             n == length(W) || throw(DimensionMismatch("length mismatch between A and W"))
             if length(V) == 0
@@ -1063,11 +1058,6 @@ for (heev, heevx, heevr, heevd, hegvd, elty, relty) in
             require_one_based_indexing(A, V, W)
             chkstride1(A, V, W)
             n = checksquare(A)
-            if $elty <: Real
-                issymmetric(A) || throw(ArgumentError("A must be symmetric"))
-            else
-                ishermitian(A) || throw(ArgumentError("A must be Hermitian"))
-            end
             chkuplofinite(A, uplo)
             if haskey(kwargs, :irange)
                 il = first(kwargs[:irange])
@@ -1175,11 +1165,6 @@ for (heev, heevx, heevr, heevd, hegvd, elty, relty) in
             require_one_based_indexing(A, V, W)
             chkstride1(A, V, W)
             n = checksquare(A)
-            if $elty <: Real
-                issymmetric(A) || throw(ArgumentError("A must be symmetric"))
-            else
-                ishermitian(A) || throw(ArgumentError("A must be Hermitian"))
-            end
             chkuplofinite(A, uplo)
             if haskey(kwargs, :irange)
                 il = first(irange)
@@ -1294,11 +1279,6 @@ for (heev, heevx, heevr, heevd, hegvd, elty, relty) in
             require_one_based_indexing(A, V, W)
             chkstride1(A, V, W)
             n = checksquare(A)
-            if $elty <: Real
-                issymmetric(A) || throw(ArgumentError("A must be symmetric"))
-            else
-                ishermitian(A) || throw(ArgumentError("A must be Hermitian"))
-            end
             uplo = 'U' # shouldn't matter but 'U' seems slightly faster than 'L'
             chkuplofinite(A, uplo)
             n == length(W) || throw(DimensionMismatch("length mismatch between A and W"))

From 6509c68e1b3d1c06f36d2874ff7d11e1ca9abf93 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Tue, 14 Oct 2025 13:48:47 -0400
Subject: [PATCH 02/10] accept keywords in all eigh solvers

---
 src/yalapack.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/yalapack.jl b/src/yalapack.jl
index 96a6240a..5bc87073 100644
--- a/src/yalapack.jl
+++ b/src/yalapack.jl
@@ -984,7 +984,8 @@ for (heev, heevx, heevr, heevd, hegvd, elty, relty) in
                 A::AbstractMatrix{$elty},
                 W::AbstractVector{$relty} = similar(A, $relty, size(A, 1)),
                 V::AbstractMatrix{$elty} = A;
-                uplo::AbstractChar = 'U'
+                uplo::AbstractChar = 'U',
+                kwargs...
             ) # shouldn't matter but 'U' seems slightly faster than 'L'
             require_one_based_indexing(A, V, W)
             chkstride1(A, V, W)
@@ -1274,7 +1275,8 @@ for (heev, heevx, heevr, heevd, hegvd, elty, relty) in
                 A::AbstractMatrix{$elty},
                 W::AbstractVector{$relty} = similar(A, $relty, size(A, 1)),
                 V::AbstractMatrix{$elty} = A;
-                uplo::AbstractChar = 'U'
+                uplo::AbstractChar = 'U',
+                kwargs...
             ) # shouldn't matter but 'U' seems slightly faster than 'L'
             require_one_based_indexing(A, V, W)
             chkstride1(A, V, W)

From 475b171a5f321433c9207c4a032179d9f44d3b1c Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Tue, 14 Oct 2025 13:51:33 -0400
Subject: [PATCH 03/10] rework hermitian checks

---
 src/common/defaults.jl      | 10 ++++++++++
 src/implementations/eigh.jl | 38 +++++++++++++++++++++++--------------
 2 files changed, 34 insertions(+), 14 deletions(-)

diff --git a/src/common/defaults.jl b/src/common/defaults.jl
index 2f1e9e50..c275ee87 100644
--- a/src/common/defaults.jl
+++ b/src/common/defaults.jl
@@ -19,3 +19,13 @@ function default_pullback_gaugetol(a)
     n = norm(a, Inf)
     return eps(eltype(n))^(3 / 4) * max(n, one(n))
 end
+
+"""
+    default_hermitian_tol(A)
+
+Default tolerance for deciding to warn if the provided `A` is not hermitian.
+"""
+function default_hermitian_tol(A)
+    n = norm(A, Inf)
+    return eps(eltype(n))^(3 / 4) * max(n, one(n))
+end
diff --git a/src/implementations/eigh.jl b/src/implementations/eigh.jl
index 7c2a2010..7f14a06e 100644
--- a/src/implementations/eigh.jl
+++ b/src/implementations/eigh.jl
@@ -8,10 +8,20 @@ copy_input(::typeof(eigh_trunc), A) = copy_input(eigh_full, A)
 
 copy_input(::typeof(eigh_full), A::Diagonal) = copy(A)
 
-function check_input(::typeof(eigh_full!), A::AbstractMatrix, DV, ::AbstractAlgorithm)
+_hermitian_tol(A, alg) = default_hermitian_tol(A)
+_hermitian_tol(A, alg::Algorithm) = get(alg.kwargs, :hermitian_tol, default_hermitian_tol(A))
+function check_hermitian(A, alg, context::Symbol)
     m, n = size(A)
     m == n || throw(DimensionMismatch("square input matrix expected"))
+    ishermitian(A; atol = _hermitian_tol(A, alg)) ||
+        throw(DomainError(A, "`eigh_$(context)!(A)` was called on a non-hermitian input matrix `A`. Try `eig_$(context)!(A)` or `eigh_$(context)(project_hermitian(A))` instead."))
+    return nothing
+end
+
+function check_input(::typeof(eigh_full!), A::AbstractMatrix, DV, alg::AbstractAlgorithm)
+    check_hermitian(A, alg, :full)
     D, V = DV
+    m = size(A, 1)
     @assert D isa Diagonal && V isa AbstractMatrix
     @check_size(D, (m, m))
     @check_scalar(D, A, real)
@@ -19,19 +29,19 @@ function check_input(::typeof(eigh_full!), A::AbstractMatrix, DV, ::AbstractAlgo
     @check_scalar(V, A)
     return nothing
 end
-function check_input(::typeof(eigh_vals!), A::AbstractMatrix, D, ::AbstractAlgorithm)
-    m, n = size(A)
-    m == n || throw(DimensionMismatch("square input matrix expected"))
+function check_input(::typeof(eigh_vals!), A::AbstractMatrix, D, alg::AbstractAlgorithm)
+    check_hermitian(A, alg, :vals)
+    m = size(A, 1)
     @assert D isa AbstractVector
-    @check_size(D, (n,))
+    @check_size(D, (m,))
     @check_scalar(D, A, real)
     return nothing
 end
 
-function check_input(::typeof(eigh_full!), A::AbstractMatrix, DV, ::DiagonalAlgorithm)
-    m, n = size(A)
-    @assert m == n && isdiag(A)
-    @assert (eltype(A) <: Real && issymmetric(A)) || ishermitian(A)
+function check_input(::typeof(eigh_full!), A::AbstractMatrix, DV, alg::DiagonalAlgorithm)
+    check_hermitian(A, alg, :full)
+    @assert isdiag(A)
+    m = size(A, 1)
     D, V = DV
     @assert D isa Diagonal && V isa Diagonal
     @check_size(D, (m, m))
@@ -40,12 +50,12 @@ function check_input(::typeof(eigh_full!), A::AbstractMatrix, DV, ::DiagonalAlgo
     @check_scalar(V, A)
     return nothing
 end
-function check_input(::typeof(eigh_vals!), A::AbstractMatrix, D, ::DiagonalAlgorithm)
-    m, n = size(A)
-    @assert m == n && isdiag(A)
-    @assert (eltype(A) <: Real && issymmetric(A)) || ishermitian(A)
+function check_input(::typeof(eigh_vals!), A::AbstractMatrix, D, alg::DiagonalAlgorithm)
+    check_hermitian(A, alg, :vals)
+    @assert isdiag(A)
+    m = size(A, 1)
     @assert D isa AbstractVector
-    @check_size(D, (n,))
+    @check_size(D, (m,))
     @check_scalar(D, A, real)
     return nothing
 end

From 4b59a7bc6e38269586f053f0ceb41c903999c702 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Tue, 14 Oct 2025 14:26:06 -0400
Subject: [PATCH 04/10] non-allocating `ishermitian_approx`

---
 src/common/matrixproperties.jl | 71 +++++++++++++++++++++++++++++-----
 1 file changed, 62 insertions(+), 9 deletions(-)

diff --git a/src/common/matrixproperties.jl b/src/common/matrixproperties.jl
index 224a16cd..c4afb368 100644
--- a/src/common/matrixproperties.jl
+++ b/src/common/matrixproperties.jl
@@ -69,19 +69,20 @@ is_right_isometric(A; kwargs...) = is_left_isometric(A'; kwargs...)
 Test whether a linear map is Hermitian, i.e. `A = A'`.
 The `isapprox_kwargs` can be used to control the tolerances of the equality.
 """
-function ishermitian(A; atol::Real = 0, rtol::Real = 0, norm = LinearAlgebra.norm, kwargs...)
+function ishermitian(A; atol::Real = 0, rtol::Real = 0, kwargs...)
     if iszero(atol) && iszero(rtol)
         return ishermitian_exact(A; kwargs...)
     else
-        return 2 * norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
+        return ishermitian_approx(A; atol, rtol, kwargs...)
     end
 end
-function ishermitian_exact(A)
-    return A == A'
-end
-function ishermitian_exact(A::StridedMatrix; kwargs...)
-    return strided_ishermitian_exact(A, Val(false); kwargs...)
+
+ishermitian_exact(A) = A == A'
+ishermitian_exact(A::StridedMatrix; kwargs...) = strided_ishermitian_exact(A, Val(false); kwargs...)
+function ishermitian_approx(A; atol, rtol, kwargs...)
+    return 2 * norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
 end
+ishermitian_approx(A::StridedMatrix; kwargs...) = strided_ishermitian_approx(A, Val(false); kwargs...)
 
 """
     isantihermitian(A; isapprox_kwargs...)
@@ -89,11 +90,11 @@ end
 Test whether a linear map is anti-Hermitian, i.e. `A = -A'`.
 The `isapprox_kwargs` can be used to control the tolerances of the equality.
 """
-function isantihermitian(A; atol::Real = 0, rtol::Real = 0, norm = LinearAlgebra.norm, kwargs...)
+function isantihermitian(A; atol::Real = 0, rtol::Real = 0, kwargs...)
     if iszero(atol) && iszero(rtol)
         return isantihermitian_exact(A; kwargs...)
     else
-        return 2 * norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
+        return isantihermitian_approx(A; atol, rtol, kwargs...)
     end
 end
 function isantihermitian_exact(A)
@@ -102,6 +103,10 @@ end
 function isantihermitian_exact(A::StridedMatrix; kwargs...)
     return strided_ishermitian_exact(A, Val(true); kwargs...)
 end
+function isantihermitian_approx(A; atol, rtol, kwargs...)
+    return 2 * norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
+end
+isantihermitian_approx(A::StridedMatrix; kwargs...) = strided_ishermitian_approx(A, Val(true); kwargs...)
 
 # blocked implementation of exact checks for strided matrices
 # -----------------------------------------------------------
@@ -139,3 +144,51 @@ function _ishermitian_exact_offdiag(Al, Au, ::Val{anti}) where {anti}
     end
     return true
 end
+
+
+function strided_ishermitian_approx(
+        A::AbstractMatrix, anti::Val;
+        blocksize = 32, atol::Real = default_hermitian_tol(A), rtol::Real = 0
+    )
+    n = size(A, 1)
+    ϵ = abs2(zero(eltype(A)))
+    ϵmax = oftype(ϵ, rtol > 0 ? max(atol, rtol * norm(A)) : atol)^2
+    for j in 1:blocksize:n
+        jb = min(blocksize, n - j + 1)
+        ϵ += _ishermitian_approx_diag(view(A, j:(j + jb - 1), j:(j + jb - 1)), anti)
+        ϵ < ϵmax || return false
+        for i in 1:blocksize:(j - 1)
+            ib = blocksize
+            ϵ += _ishermitian_approx_offdiag(
+                view(A, i:(i + ib - 1), j:(j + jb - 1)),
+                view(A, j:(j + jb - 1), i:(i + ib - 1)),
+                anti
+            )
+            ϵ < ϵmax || return false
+        end
+    end
+    return true
+end
+
+function _ishermitian_approx_diag(A, ::Val{anti}) where {anti}
+    n = size(A, 1)
+    ϵ = abs2(zero(eltype(A)))
+    @inbounds for j in 1:n
+        @simd for i in 1:j
+            val = anti ? (A[i, j] + adjoint(A[j, i])) : (A[i, j] - adjoint(A[j, i]))
+            ϵ += abs2(val)
+        end
+    end
+    return ϵ
+end
+function _ishermitian_approx_offdiag(Al, Au, ::Val{anti}) where {anti}
+    m, n = size(Al) # == reverse(size(Al))
+    ϵ = abs2(zero(eltype(Al)))
+    @inbounds for j in 1:n
+        @simd for i in 1:m
+            val = anti ? (Al[i, j] + adjoint(Au[j, i])) : (Al[i, j] - adjoint(Au[j, i]))
+            ϵ += abs2(val)
+        end
+    end
+    return ϵ
+end

From a658c42c7a412e8dce08d57b12b5688cd2616995 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Tue, 14 Oct 2025 19:24:40 -0400
Subject: [PATCH 05/10] code suggestions

---
 src/common/matrixproperties.jl | 24 ++++++++++++------------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/common/matrixproperties.jl b/src/common/matrixproperties.jl
index c4afb368..c7511d8d 100644
--- a/src/common/matrixproperties.jl
+++ b/src/common/matrixproperties.jl
@@ -151,20 +151,20 @@ function strided_ishermitian_approx(
         blocksize = 32, atol::Real = default_hermitian_tol(A), rtol::Real = 0
     )
     n = size(A, 1)
-    ϵ = abs2(zero(eltype(A)))
-    ϵmax = oftype(ϵ, rtol > 0 ? max(atol, rtol * norm(A)) : atol)^2
+    ϵ² = abs2(zero(eltype(A)))
+    ϵ²max = oftype(ϵ², rtol > 0 ? max(atol, rtol * norm(A)) : atol)^2
     for j in 1:blocksize:n
         jb = min(blocksize, n - j + 1)
-        ϵ += _ishermitian_approx_diag(view(A, j:(j + jb - 1), j:(j + jb - 1)), anti)
-        ϵ < ϵmax || return false
+        ϵ² += _ishermitian_approx_diag(view(A, j:(j + jb - 1), j:(j + jb - 1)), anti)
+        ϵ² < ϵ²max || return false
         for i in 1:blocksize:(j - 1)
             ib = blocksize
-            ϵ += _ishermitian_approx_offdiag(
+            ϵ² += _ishermitian_approx_offdiag(
                 view(A, i:(i + ib - 1), j:(j + jb - 1)),
                 view(A, j:(j + jb - 1), i:(i + ib - 1)),
                 anti
             )
-            ϵ < ϵmax || return false
+            ϵ² < ϵ²max || return false
         end
     end
     return true
@@ -172,23 +172,23 @@ end
 
 function _ishermitian_approx_diag(A, ::Val{anti}) where {anti}
     n = size(A, 1)
-    ϵ = abs2(zero(eltype(A)))
+    ϵ² = abs2(zero(eltype(A)))
     @inbounds for j in 1:n
         @simd for i in 1:j
             val = anti ? (A[i, j] + adjoint(A[j, i])) : (A[i, j] - adjoint(A[j, i]))
-            ϵ += abs2(val)
+            ϵ² += abs2(val) * (1 + Int(i < j))
         end
     end
-    return ϵ
+    return ϵ²
 end
 function _ishermitian_approx_offdiag(Al, Au, ::Val{anti}) where {anti}
     m, n = size(Al) # == reverse(size(Al))
-    ϵ = abs2(zero(eltype(Al)))
+    ϵ² = abs2(zero(eltype(Al)))
     @inbounds for j in 1:n
         @simd for i in 1:m
             val = anti ? (Al[i, j] + adjoint(Au[j, i])) : (Al[i, j] - adjoint(Au[j, i]))
-            ϵ += abs2(val)
+            ϵ² += abs2(val)
         end
     end
-    return ϵ
+    return 2ϵ²
 end

From 1e206bc70af7e123e141a3c3ae661be1963a7575 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Wed, 15 Oct 2025 08:44:12 -0400
Subject: [PATCH 06/10] refactor `check_hermitian`

---
 src/implementations/eigh.jl | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/src/implementations/eigh.jl b/src/implementations/eigh.jl
index 7f14a06e..0cfa6db0 100644
--- a/src/implementations/eigh.jl
+++ b/src/implementations/eigh.jl
@@ -8,18 +8,18 @@ copy_input(::typeof(eigh_trunc), A) = copy_input(eigh_full, A)
 
 copy_input(::typeof(eigh_full), A::Diagonal) = copy(A)
 
-_hermitian_tol(A, alg) = default_hermitian_tol(A)
-_hermitian_tol(A, alg::Algorithm) = get(alg.kwargs, :hermitian_tol, default_hermitian_tol(A))
-function check_hermitian(A, alg, context::Symbol)
+check_hermitian(A, ::AbstractAlgorithm) = check_hermitian(A)
+check_hermitian(A, alg::Algorithm) = check_hermitian(A; atol = get(alg.kwargs, :hermitian_tol, default_hermitian_tol(A)))
+function check_hermitian(A; atol::Real = default_hermitian_tol(A), rtol::Real = 0)
     m, n = size(A)
     m == n || throw(DimensionMismatch("square input matrix expected"))
-    ishermitian(A; atol = _hermitian_tol(A, alg)) ||
-        throw(DomainError(A, "`eigh_$(context)!(A)` was called on a non-hermitian input matrix `A`. Try `eig_$(context)!(A)` or `eigh_$(context)(project_hermitian(A))` instead."))
+    ishermitian(A; atol, rtol) ||
+        throw(DomainError(A, "Hermitian matrix was expected. Use `project_hermitian` to project onto the nearest hermitian matrix."))
     return nothing
 end
 
 function check_input(::typeof(eigh_full!), A::AbstractMatrix, DV, alg::AbstractAlgorithm)
-    check_hermitian(A, alg, :full)
+    check_hermitian(A, alg)
     D, V = DV
     m = size(A, 1)
     @assert D isa Diagonal && V isa AbstractMatrix
@@ -30,7 +30,7 @@ function check_input(::typeof(eigh_full!), A::AbstractMatrix, DV, alg::AbstractA
     return nothing
 end
 function check_input(::typeof(eigh_vals!), A::AbstractMatrix, D, alg::AbstractAlgorithm)
-    check_hermitian(A, alg, :vals)
+    check_hermitian(A, alg)
     m = size(A, 1)
     @assert D isa AbstractVector
     @check_size(D, (m,))
@@ -39,7 +39,7 @@ function check_input(::typeof(eigh_vals!), A::AbstractMatrix, D, alg::AbstractAl
 end
 
 function check_input(::typeof(eigh_full!), A::AbstractMatrix, DV, alg::DiagonalAlgorithm)
-    check_hermitian(A, alg, :full)
+    check_hermitian(A, alg)
     @assert isdiag(A)
     m = size(A, 1)
     D, V = DV
@@ -51,7 +51,7 @@ function check_input(::typeof(eigh_full!), A::AbstractMatrix, DV, alg::DiagonalA
     return nothing
 end
 function check_input(::typeof(eigh_vals!), A::AbstractMatrix, D, alg::DiagonalAlgorithm)
-    check_hermitian(A, alg, :vals)
+    check_hermitian(A, alg)
     @assert isdiag(A)
     m = size(A, 1)
     @assert D isa AbstractVector

From d76918d716479272d69d6ff0268d12398a268fa1 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Wed, 15 Oct 2025 09:37:16 -0400
Subject: [PATCH 07/10] rework tests

---
 src/common/matrixproperties.jl     | 12 ++++++------
 src/implementations/projections.jl |  8 +++++---
 test/projections.jl                | 16 +++++++++++++++-
 3 files changed, 26 insertions(+), 10 deletions(-)

diff --git a/src/common/matrixproperties.jl b/src/common/matrixproperties.jl
index c7511d8d..925f92e6 100644
--- a/src/common/matrixproperties.jl
+++ b/src/common/matrixproperties.jl
@@ -80,7 +80,7 @@ end
 ishermitian_exact(A) = A == A'
 ishermitian_exact(A::StridedMatrix; kwargs...) = strided_ishermitian_exact(A, Val(false); kwargs...)
 function ishermitian_approx(A; atol, rtol, kwargs...)
-    return 2 * norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
+    return norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
 end
 ishermitian_approx(A::StridedMatrix; kwargs...) = strided_ishermitian_approx(A, Val(false); kwargs...)
 
@@ -104,7 +104,7 @@ function isantihermitian_exact(A::StridedMatrix; kwargs...)
     return strided_ishermitian_exact(A, Val(true); kwargs...)
 end
 function isantihermitian_approx(A; atol, rtol, kwargs...)
-    return 2 * norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
+    return norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
 end
 isantihermitian_approx(A::StridedMatrix; kwargs...) = strided_ishermitian_approx(A, Val(true); kwargs...)
 
@@ -159,7 +159,7 @@ function strided_ishermitian_approx(
         ϵ² < ϵ²max || return false
         for i in 1:blocksize:(j - 1)
             ib = blocksize
-            ϵ² += _ishermitian_approx_offdiag(
+            ϵ² += 2 * _ishermitian_approx_offdiag(
                 view(A, i:(i + ib - 1), j:(j + jb - 1)),
                 view(A, j:(j + jb - 1), i:(i + ib - 1)),
                 anti
@@ -175,7 +175,7 @@ function _ishermitian_approx_diag(A, ::Val{anti}) where {anti}
     ϵ² = abs2(zero(eltype(A)))
     @inbounds for j in 1:n
         @simd for i in 1:j
-            val = anti ? (A[i, j] + adjoint(A[j, i])) : (A[i, j] - adjoint(A[j, i]))
+            val = _project_hermitian(A[i, j], A[j, i], !anti)
             ϵ² += abs2(val) * (1 + Int(i < j))
         end
     end
@@ -186,9 +186,9 @@ function _ishermitian_approx_offdiag(Al, Au, ::Val{anti}) where {anti}
     ϵ² = abs2(zero(eltype(Al)))
     @inbounds for j in 1:n
         @simd for i in 1:m
-            val = anti ? (Al[i, j] + adjoint(Au[j, i])) : (Al[i, j] - adjoint(Au[j, i]))
+            val = _project_hermitian(Al[i, j], Au[j, i], !anti)
             ϵ² += abs2(val)
         end
     end
-    return 2ϵ²
+    return ϵ²
 end
diff --git a/src/implementations/projections.jl b/src/implementations/projections.jl
index b60cfd89..8baf09a5 100644
--- a/src/implementations/projections.jl
+++ b/src/implementations/projections.jl
@@ -85,14 +85,16 @@ function project_hermitian_native!(A::AbstractMatrix, B::AbstractMatrix, anti::V
     return B
 end
 
+@inline function _project_hermitian(Aij::Number, Aji::Number, anti::Bool)
+    return anti ? (Aij - Aji') / 2 : (Aij + Aji') / 2
+end
 function _project_hermitian_offdiag!(
         Au::AbstractMatrix, Al::AbstractMatrix, Bu::AbstractMatrix, Bl::AbstractMatrix, ::Val{anti}
     ) where {anti}
-
     m, n = size(Au) # == reverse(size(Au))
     return @inbounds for j in 1:n
         @simd for i in 1:m
-            val = anti ? (Au[i, j] - adjoint(Al[j, i])) / 2 : (Au[i, j] + adjoint(Al[j, i])) / 2
+            val = _project_hermitian(Au[i, j], Al[j, i], anti)
             Bu[i, j] = val
             aval = adjoint(val)
             Bl[j, i] = anti ? -aval : aval
@@ -104,7 +106,7 @@ function _project_hermitian_diag!(A::AbstractMatrix, B::AbstractMatrix, ::Val{an
     n = size(A, 1)
     @inbounds for j in 1:n
         @simd for i in 1:(j - 1)
-            val = anti ? (A[i, j] - adjoint(A[j, i])) / 2 : (A[i, j] + adjoint(A[j, i])) / 2
+            val = _project_hermitian(A[i, j], A[j, i], anti)
             B[i, j] = val
             aval = adjoint(val)
             B[j, i] = anti ? -aval : aval
diff --git a/test/projections.jl b/test/projections.jl
index 684b3bac..3825d91e 100644
--- a/test/projections.jl
+++ b/test/projections.jl
@@ -2,7 +2,7 @@ using MatrixAlgebraKit
 using Test
 using TestExtras
 using StableRNGs
-using LinearAlgebra: LinearAlgebra, Diagonal, norm
+using LinearAlgebra: LinearAlgebra, Diagonal, norm, normalize!
 
 const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64)
 
@@ -43,6 +43,20 @@ const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64)
         @test isantihermitian(Ba)
         @test Ba ≈ Aa
     end
+
+    # test approximate error calculation
+    A = normalize!(randn(rng, T, m, m))
+    Ah = project_hermitian(A)
+    Ah_approx = Ah + noisefactor * Aa
+    ϵ = norm(project_antihermitian(Ah_approx))
+    @test !ishermitian(Ah_approx; atol = (999 // 1000) * ϵ)
+    @test ishermitian(Ah_approx; atol = (1001 // 1000) * ϵ)
+
+    Aa = project_antihermitian(A)
+    Aa_approx = Aa + noisefactor * Ah
+    ϵ = norm(project_hermitian(Aa_approx))
+    @test !isantihermitian(Aa_approx; atol = (999 // 1000) * ϵ)
+    @test isantihermitian(Aa_approx; atol = (1001 // 1000) * ϵ)
 end
 
 @testset "project_isometric! for T = $T" for T in BLASFloats

From b889afeb7709e32c11f4b4be4b5ce30583154ccb Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Wed, 15 Oct 2025 16:22:12 -0400
Subject: [PATCH 08/10] allocating path for GPU

---
 .../MatrixAlgebraKitAMDGPUExt.jl               | 13 ++++++++++---
 .../MatrixAlgebraKitCUDAExt.jl                 | 18 +++++++++++++-----
 2 files changed, 23 insertions(+), 8 deletions(-)

diff --git a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl
index e1d24de9..e8614461 100644
--- a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl
+++ b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl
@@ -128,10 +128,17 @@ function MatrixAlgebraKit._project_hermitian_diag!(A::StridedROCMatrix, B::Strid
 end
 
 MatrixAlgebraKit.ishermitian_exact(A::StridedROCMatrix) = all(A .== adjoint(A))
-MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T} = all(A.diag .== adjoint(A.diag))
+MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T} =
+    all(A.diag .== adjoint(A.diag))
+MatrixAlgebraKit.ishermitian_approx(A::StridedROCMatrix; kwargs...) =
+    @invoke MatrixAlgebraKit.ishermitian_approx(A::Any; kwargs...)
 
-MatrixAlgebraKit.isantihermitian_exact(A::StridedROCMatrix) = all(A .== -adjoint(A))
-MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T} = all(A.diag .== -adjoint(A.diag))
+MatrixAlgebraKit.isantihermitian_exact(A::StridedROCMatrix) =
+    all(A .== -adjoint(A))
+MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedROCVector{T}}) where {T} =
+    all(A.diag .== -adjoint(A.diag))
+MatrixAlgebraKit.isantihermitian_approx(A::StridedROCMatrix; kwargs...) =
+    @invoke MatrixAlgebraKit.isantihermitian_approx(A::Any; kwargs...)
 
 function MatrixAlgebraKit._avgdiff!(A::StridedROCMatrix, B::StridedROCMatrix)
     axes(A) == axes(B) || throw(DimensionMismatch())
diff --git a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl
index 4a2c77e3..0bcb61cb 100644
--- a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl
+++ b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl
@@ -134,11 +134,19 @@ function MatrixAlgebraKit._project_hermitian_diag!(A::StridedCuMatrix, B::Stride
     return nothing
 end
 
-MatrixAlgebraKit.ishermitian_exact(A::StridedCuMatrix) = all(A .== adjoint(A))
-MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T} = all(A.diag .== adjoint(A.diag))
-
-MatrixAlgebraKit.isantihermitian_exact(A::StridedCuMatrix) = all(A .== -adjoint(A))
-MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T} = all(A.diag .== -adjoint(A.diag))
+MatrixAlgebraKit.ishermitian_exact(A::StridedCuMatrix) =
+    all(A .== adjoint(A))
+MatrixAlgebraKit.ishermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T} =
+    all(A.diag .== adjoint(A.diag))
+MatrixAlgebraKit.ishermitian_approx(A::StridedCuMatrix; kwargs...) =
+    @invoke MatrixAlgebraKit.ishermitian_approx(A::Any; kwargs...)
+
+MatrixAlgebraKit.isantihermitian_exact(A::StridedCuMatrix) =
+    all(A .== -adjoint(A))
+MatrixAlgebraKit.isantihermitian_exact(A::Diagonal{T, <:StridedCuVector{T}}) where {T} =
+    all(A.diag .== -adjoint(A.diag))
+MatrixAlgebraKit.isantihermitian_approx(A::StridedCuMatrix; kwargs...) =
+    @invoke MatrixAlgebraKit.isantihermitian_approx(A::Any; kwargs...)
 
 function MatrixAlgebraKit._avgdiff!(A::StridedCuMatrix, B::StridedCuMatrix)
     axes(A) == axes(B) || throw(DimensionMismatch())

From be4344d54664442710192b8f34e1f662aa1f681c Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Wed, 15 Oct 2025 17:58:37 -0400
Subject: [PATCH 09/10] small test fix

---
 test/projections.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/test/projections.jl b/test/projections.jl
index 3825d91e..56bc1a04 100644
--- a/test/projections.jl
+++ b/test/projections.jl
@@ -47,12 +47,13 @@ const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64)
     # test approximate error calculation
     A = normalize!(randn(rng, T, m, m))
     Ah = project_hermitian(A)
+    Aa = project_antihermitian(A)
+
     Ah_approx = Ah + noisefactor * Aa
     ϵ = norm(project_antihermitian(Ah_approx))
     @test !ishermitian(Ah_approx; atol = (999 // 1000) * ϵ)
     @test ishermitian(Ah_approx; atol = (1001 // 1000) * ϵ)
 
-    Aa = project_antihermitian(A)
     Aa_approx = Aa + noisefactor * Ah
     ϵ = norm(project_hermitian(Aa_approx))
     @test !isantihermitian(Aa_approx; atol = (999 // 1000) * ϵ)

From cd5cd6ed9e7b0622aaf04f84f9915f896c271d2d Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Wed, 15 Oct 2025 18:00:56 -0400
Subject: [PATCH 10/10] add missing imports

---
 ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl | 2 +-
 ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl     | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl
index e8614461..8f495f46 100644
--- a/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl
+++ b/ext/MatrixAlgebraKitAMDGPUExt/MatrixAlgebraKitAMDGPUExt.jl
@@ -4,7 +4,7 @@ using MatrixAlgebraKit
 using MatrixAlgebraKit: @algdef, Algorithm, check_input
 using MatrixAlgebraKit: one!, zero!, uppertriangular!, lowertriangular!
 using MatrixAlgebraKit: diagview, sign_safe
-using MatrixAlgebraKit: LQViaTransposedQR, TruncationByValue
+using MatrixAlgebraKit: LQViaTransposedQR, TruncationByValue, AbstractAlgorithm
 using MatrixAlgebraKit: default_qr_algorithm, default_lq_algorithm, default_svd_algorithm, default_eigh_algorithm
 import MatrixAlgebraKit: _gpu_geqrf!, _gpu_ungqr!, _gpu_unmqr!, _gpu_gesvd!, _gpu_Xgesvdp!, _gpu_gesvdj!
 import MatrixAlgebraKit: _gpu_heevj!, _gpu_heevd!, _gpu_heev!, _gpu_heevx!
diff --git a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl
index 0bcb61cb..68d25c36 100644
--- a/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl
+++ b/ext/MatrixAlgebraKitCUDAExt/MatrixAlgebraKitCUDAExt.jl
@@ -4,7 +4,7 @@ using MatrixAlgebraKit
 using MatrixAlgebraKit: @algdef, Algorithm, check_input
 using MatrixAlgebraKit: one!, zero!, uppertriangular!, lowertriangular!
 using MatrixAlgebraKit: diagview, sign_safe
-using MatrixAlgebraKit: LQViaTransposedQR, TruncationByValue
+using MatrixAlgebraKit: LQViaTransposedQR, TruncationByValue, AbstractAlgorithm
 using MatrixAlgebraKit: default_qr_algorithm, default_lq_algorithm, default_svd_algorithm, default_eig_algorithm, default_eigh_algorithm
 import MatrixAlgebraKit: _gpu_geqrf!, _gpu_ungqr!, _gpu_unmqr!, _gpu_gesvd!, _gpu_Xgesvdp!, _gpu_Xgesvdr!, _gpu_gesvdj!, _gpu_geev!
 import MatrixAlgebraKit: _gpu_heevj!, _gpu_heevd!