From 67cb4b9c5efd9de14e90a19aacefdd2d11e66aad Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Wed, 1 Oct 2025 14:39:03 -0400
Subject: [PATCH 01/13] Add `hermitianpart!` and `antihermitianpart!`

---
 src/common/initialization.jl | 50 ++++++++++++++++++++++++++++++++++++
 1 file changed, 50 insertions(+)

diff --git a/src/common/initialization.jl b/src/common/initialization.jl
index cdabcf43..5903da5d 100644
--- a/src/common/initialization.jl
+++ b/src/common/initialization.jl
@@ -30,3 +30,53 @@ function lowertriangular!(A::AbstractMatrix)
     end
     return A
 end
+
+@doc """
+    hermitianpart(A)
+    hermitianpart!(A)
+
+In-place implementation of the Hermitian part of a (square) matrix `A`, defined as `(A + A') / 2`.
+For real matrices this is also called the symmetric part of `A`.
+
+See also [`antihermitianpart`](@ref).
+""" hermitianpart, hermitianpart!
+
+hermitianpart(A) = hermitianpart!(copy(A))
+function hermitianpart!(A::AbstractMatrix)
+    Base.require_one_based_indexing(A)
+    n = LinearAlgebra.checksquare(A)
+    @inbounds for j in 1:n
+        A[j, j] = real(A[j, j])
+        for i in 1:(j - 1)
+            val = (A[i, j] + adjoint(A[j, i])) / 2
+            A[i, j] = val
+            A[j, i] = adjoint(val)
+        end
+    end
+    return A
+end
+
+@doc """
+    antihermitianpart(A)
+    antihermitianpart!(A)
+
+In-place implementation of the anti-Hermitian part of a (square) matrix `A`, defined as `(A - A') / 2`.
+For real matrices this is also called the anti-symmetric part of `A`.
+
+See also [`hermitianpart`](@ref).
+""" antihermitianpart, antihermitianpart!
+
+antihermitianpart(A) = antihermitianpart!(copy(A))
+function antihermitianpart!(A::AbstractMatrix)
+    Base.require_one_based_indexing(A)
+    n = LinearAlgebra.checksquare(A)
+    @inbounds for j in 1:n
+        A[j, j] = imag(A[j, j]) * im
+        for i in 1:(j - 1)
+            val = (A[i, j] - adjoint(A[j, i])) / 2
+            A[i, j] = val
+            A[j, i] = -adjoint(val)
+        end
+    end
+    return A
+end

From 3e3eb130cd421424bb7a22ca58d70575f001acb5 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Wed, 1 Oct 2025 15:03:46 -0400
Subject: [PATCH 02/13] Add `ishermitian` and `isantihermitian`

---
 src/MatrixAlgebraKit.jl        |  2 +-
 src/common/matrixproperties.jl | 64 ++++++++++++++++++++++++++++++++++
 2 files changed, 65 insertions(+), 1 deletion(-)

diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl
index 7dc44b72..f0f319bb 100644
--- a/src/MatrixAlgebraKit.jl
+++ b/src/MatrixAlgebraKit.jl
@@ -4,7 +4,7 @@ using LinearAlgebra: LinearAlgebra
 using LinearAlgebra: norm # TODO: eleminate if we use VectorInterface.jl?
 using LinearAlgebra: mul!, rmul!, lmul!, adjoint!, rdiv!, ldiv!
 using LinearAlgebra: sylvester
-using LinearAlgebra: isposdef, ishermitian, issymmetric
+using LinearAlgebra: isposdef, issymmetric
 using LinearAlgebra: Diagonal, diag, diagind, isdiag
 using LinearAlgebra: UpperTriangular, LowerTriangular
 using LinearAlgebra: BlasFloat, BlasReal, BlasComplex, BlasInt
diff --git a/src/common/matrixproperties.jl b/src/common/matrixproperties.jl
index d9541cab..9d0c3771 100644
--- a/src/common/matrixproperties.jl
+++ b/src/common/matrixproperties.jl
@@ -57,3 +57,67 @@ See also [`isisometry`](@ref) and [`is_left_isometry`](@ref).
 function is_right_isometry(A::AbstractMatrix; isapprox_kwargs...)
     return isapprox(A * A', LinearAlgebra.I; isapprox_kwargs...)
 end
+
+"""
+    ishermitian(A; isapprox_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, kwargs...)
+    return (atol == rtol == 0) ? (A == A') : isapprox(A, A'; atol, rtol, kwargs...)
+end
+function ishermitian(A::AbstractMatrix; atol::Real = 0, rtol::Real = 0, norm = LinearAlgebra.norm, kwargs...)
+    Base.require_one_based_indexing(A)
+    m, n = size(A)
+    m == n || return false
+    if atol == rtol == 0
+        for j in 1:n
+            for i in 1:j
+                A[i, j] == adjoint(A[j, i]) || return false
+            end
+        end
+    elseif norm === LinearAlgebra.norm
+        atol = max(atol, rtol * norm(A))
+        for j in 1:n
+            for i in 1:j
+                isapprox(A[i, j], adjoint(A[j, i]); atol, abs, kwargs...)  || return false
+            end
+        end
+    else
+        return isapprox(A, A'; atol, rtol, norm, kwargs...)
+    end
+    return true
+end
+
+"""
+    isantihermitian(A; isapprox_kwargs...)
+
+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, kwargs...)
+    return (atol == 0 & rtol == 0) ? (A == -A') : isapprox(A, -A'; atol, rtol, kwargs...)
+end
+function isantihermitian(A::AbstractMatrix; atol::Real = 0, rtol::Real = 0, norm = LinearAlgebra.norm, kwargs...)
+    Base.require_one_based_indexing(A)
+    m, n = size(A)
+    m == n || return false
+    if atol == rtol == 0
+        @inbounds for j in 1:n
+            for i in 1:j
+                A[i, j] == -adjoint(A[j, i]) || return false
+            end
+        end
+    elseif norm === LinearAlgebra.norm
+        atol = max(atol, rtol * norm(A))
+        @inbounds for j in 1:n
+            for i in 1:j
+                isapprox(A[i, j], -adjoint(A[j, i]); atol, abs, kwargs...)  || return false
+            end
+        end
+    else
+        return isapprox(A, -A'; atol, rtol, norm, kwargs...)
+    end
+    return true
+end

From 935ff795c692792bb4c59ad98e4e4e3180915988 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Wed, 1 Oct 2025 15:31:41 -0400
Subject: [PATCH 03/13] Start using in pullbacks

---
 src/pullbacks/eigh.jl  | 6 +++---
 src/pullbacks/lq.jl    | 4 ++--
 src/pullbacks/polar.jl | 4 ++--
 src/pullbacks/qr.jl    | 4 ++--
 src/pullbacks/svd.jl   | 8 ++++----
 5 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/src/pullbacks/eigh.jl b/src/pullbacks/eigh.jl
index 243d476b..cb000a06 100644
--- a/src/pullbacks/eigh.jl
+++ b/src/pullbacks/eigh.jl
@@ -42,7 +42,7 @@ function eigh_pullback!(
         indV = axes(V, 2)[ind]
         length(indV) == pV || throw(DimensionMismatch())
         mul!(view(VᴴΔV, :, indV), V', ΔV)
-        aVᴴΔV = rmul!(VᴴΔV - VᴴΔV', 1 / 2)
+        aVᴴΔV = antihermitianpart(VᴴΔV) # can't use in-place or recycling doesn't work
 
         mask = abs.(D' .- D) .< degeneracy_atol
         Δgauge = norm(view(aVᴴΔV, mask))
@@ -58,7 +58,7 @@ function eigh_pullback!(
             length(indD) == pD || throw(DimensionMismatch())
             view(diagview(aVᴴΔV), indD) .+= real.(ΔDvec)
         end
-        # recylce VdΔV space
+        # recycle VdΔV space
         ΔA = mul!(ΔA, mul!(VᴴΔV, V, aVᴴΔV), V', 1, 1)
     elseif !iszerotangent(ΔDmat)
         ΔDvec = diagview(ΔDmat)
@@ -112,7 +112,7 @@ function eigh_trunc_pullback!(
     if !iszerotangent(ΔV)
         (n, p) == size(ΔV) || throw(DimensionMismatch())
         VᴴΔV = V' * ΔV
-        aVᴴΔV = rmul!(VᴴΔV - VᴴΔV', 1 / 2)
+        aVᴴΔV = antihermitianpart!(VᴴΔV)
 
         mask = abs.(D' .- D) .< degeneracy_atol
         Δgauge = norm(view(aVᴴΔV, mask))
diff --git a/src/pullbacks/lq.jl b/src/pullbacks/lq.jl
index 5b4db9f3..b130a6d7 100644
--- a/src/pullbacks/lq.jl
+++ b/src/pullbacks/lq.jl
@@ -118,8 +118,8 @@ function lq_null_pullback!(
         gauge_atol::Real = tol
     )
     if !iszerotangent(ΔNᴴ) && size(Nᴴ, 1) > 0
-        NᴴΔN = Nᴴ * ΔNᴴ'
-        Δgauge = norm((NᴴΔN .- NᴴΔN') ./ 2)
+        aNᴴΔN = antihermitianpart!(Nᴴ * ΔNᴴ')
+        Δgauge = norm(aNᴴΔN)
         Δgauge < tol ||
             @warn "`lq_null` cotangent sensitive to gauge choice: (|Δgauge| = $Δgauge)"
         L, Q = lq_compact(A; positive = true) # should we be able to provide algorithm here?
diff --git a/src/pullbacks/polar.jl b/src/pullbacks/polar.jl
index 522b3201..bc5c3877 100644
--- a/src/pullbacks/polar.jl
+++ b/src/pullbacks/polar.jl
@@ -11,7 +11,7 @@ function left_polar_pullback!(ΔA::AbstractMatrix, A, WP, ΔWP)
     # Extract and check the cotangents
     ΔW, ΔP = ΔWP
     if !iszerotangent(ΔP)
-        ΔP = (ΔP + ΔP') / 2
+        ΔP = hermitianpart(ΔP)
     end
     M = zero(P)
     !iszerotangent(ΔW) && mul!(M, W', ΔW, 1, 1)
@@ -41,7 +41,7 @@ function right_polar_pullback!(ΔA::AbstractMatrix, A, PWᴴ, ΔPWᴴ)
     # Extract and check the cotangents
     ΔP, ΔWᴴ = ΔPWᴴ
     if !iszerotangent(ΔP)
-        ΔP = (ΔP + ΔP') / 2
+        ΔP = hermitianpart(ΔP)
     end
     M = zero(P)
     !iszerotangent(ΔWᴴ) && mul!(M, ΔWᴴ, Wᴴ', 1, 1)
diff --git a/src/pullbacks/qr.jl b/src/pullbacks/qr.jl
index e1e8ac6c..518dcd27 100644
--- a/src/pullbacks/qr.jl
+++ b/src/pullbacks/qr.jl
@@ -117,8 +117,8 @@ function qr_null_pullback!(
         gauge_atol::Real = tol
     )
     if !iszerotangent(ΔN) && size(N, 2) > 0
-        NᴴΔN = N' * ΔN
-        Δgauge = norm((NᴴΔN .- NᴴΔN') ./ 2)
+        aNᴴΔN = antihermitianpart!(N' * ΔN)
+        Δgauge = norm(aNᴴΔN)
         Δgauge < tol ||
             @warn "`qr_null` cotangent sensitive to gauge choice: (|Δgauge| = $Δgauge)"
 
diff --git a/src/pullbacks/svd.jl b/src/pullbacks/svd.jl
index c5f6e98b..40c61884 100644
--- a/src/pullbacks/svd.jl
+++ b/src/pullbacks/svd.jl
@@ -69,8 +69,8 @@ function svd_pullback!(
     end
 
     # Project onto antihermitian part; hermitian part outside of Grassmann tangent space
-    aUΔU = rmul!(UΔU - UΔU', 1 / 2)
-    aVΔV = rmul!(VΔV - VΔV', 1 / 2)
+    aUΔU = antihermitianpart!(UΔU)
+    aVΔV = antihermitianpart!(VΔV)
 
     # check whether cotangents arise from gauge-invariance objective function
     mask = abs.(Sr' .- Sr) .< degeneracy_atol
@@ -159,8 +159,8 @@ function svd_trunc_pullback!(
     end
 
     # Project onto antihermitian part; hermitian part outside of Grassmann tangent space
-    aUΔU = rmul!(UΔU - UΔU', 1 / 2)
-    aVΔV = rmul!(VΔV - VΔV', 1 / 2)
+    aUΔU = antihermitianpart!(UΔU)
+    aVΔV = antihermitianpart!(VΔV)
 
     # check whether cotangents arise from gauge-invariance objective function
     mask = abs.(S' .- S) .< degeneracy_atol

From f0044ae9117bed781222684f45f6248e75bf345b Mon Sep 17 00:00:00 2001
From: Jutho Haegeman <jutho.haegeman@ugent.be>
Date: Sun, 5 Oct 2025 09:23:55 +0200
Subject: [PATCH 04/13] blocked (anti)hermitian

---
 src/MatrixAlgebraKit.jl          |  7 ++-
 src/common/initialization.jl     | 50 -----------------
 src/implementations/hermitian.jl | 95 ++++++++++++++++++++++++++++++++
 src/interface/hermitian.jl       | 45 +++++++++++++++
 test/projections.jl              | 39 +++++++++++++
 test/runtests.jl                 |  3 +
 6 files changed, 188 insertions(+), 51 deletions(-)
 create mode 100644 src/implementations/hermitian.jl
 create mode 100644 src/interface/hermitian.jl
 create mode 100644 test/projections.jl

diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl
index f0f319bb..af66d7c6 100644
--- a/src/MatrixAlgebraKit.jl
+++ b/src/MatrixAlgebraKit.jl
@@ -9,8 +9,10 @@ using LinearAlgebra: Diagonal, diag, diagind, isdiag
 using LinearAlgebra: UpperTriangular, LowerTriangular
 using LinearAlgebra: BlasFloat, BlasReal, BlasComplex, BlasInt
 
-export isisometry, isunitary
+export isisometry, isunitary, ishermitian, isantihermitian
 
+export hermitianpart, antihermitianpart
+export hermitianpart!, antihermitianpart!
 export qr_compact, qr_full, qr_null, lq_compact, lq_full, lq_null
 export qr_compact!, qr_full!, qr_null!, lq_compact!, lq_full!, lq_null!
 export svd_compact, svd_full, svd_vals, svd_trunc
@@ -33,6 +35,7 @@ export LAPACK_HouseholderQR, LAPACK_HouseholderLQ, LAPACK_Simple, LAPACK_Expert,
     LAPACK_DivideAndConquer, LAPACK_Jacobi
 export LQViaTransposedQR
 export DiagonalAlgorithm
+export NativeBlocked
 export CUSOLVER_Simple, CUSOLVER_HouseholderQR, CUSOLVER_QRIteration, CUSOLVER_SVDPolar,
     CUSOLVER_Jacobi, CUSOLVER_Randomized, CUSOLVER_DivideAndConquer
 export ROCSOLVER_HouseholderQR, ROCSOLVER_QRIteration, ROCSOLVER_Jacobi,
@@ -74,6 +77,7 @@ include("common/gauge.jl")
 
 include("yalapack.jl")
 include("algorithms.jl")
+include("interface/hermitian.jl")
 include("interface/decompositions.jl")
 include("interface/truncation.jl")
 include("interface/qr.jl")
@@ -86,6 +90,7 @@ include("interface/schur.jl")
 include("interface/polar.jl")
 include("interface/orthnull.jl")
 
+include("implementations/hermitian.jl")
 include("implementations/truncation.jl")
 include("implementations/qr.jl")
 include("implementations/lq.jl")
diff --git a/src/common/initialization.jl b/src/common/initialization.jl
index 5903da5d..cdabcf43 100644
--- a/src/common/initialization.jl
+++ b/src/common/initialization.jl
@@ -30,53 +30,3 @@ function lowertriangular!(A::AbstractMatrix)
     end
     return A
 end
-
-@doc """
-    hermitianpart(A)
-    hermitianpart!(A)
-
-In-place implementation of the Hermitian part of a (square) matrix `A`, defined as `(A + A') / 2`.
-For real matrices this is also called the symmetric part of `A`.
-
-See also [`antihermitianpart`](@ref).
-""" hermitianpart, hermitianpart!
-
-hermitianpart(A) = hermitianpart!(copy(A))
-function hermitianpart!(A::AbstractMatrix)
-    Base.require_one_based_indexing(A)
-    n = LinearAlgebra.checksquare(A)
-    @inbounds for j in 1:n
-        A[j, j] = real(A[j, j])
-        for i in 1:(j - 1)
-            val = (A[i, j] + adjoint(A[j, i])) / 2
-            A[i, j] = val
-            A[j, i] = adjoint(val)
-        end
-    end
-    return A
-end
-
-@doc """
-    antihermitianpart(A)
-    antihermitianpart!(A)
-
-In-place implementation of the anti-Hermitian part of a (square) matrix `A`, defined as `(A - A') / 2`.
-For real matrices this is also called the anti-symmetric part of `A`.
-
-See also [`hermitianpart`](@ref).
-""" antihermitianpart, antihermitianpart!
-
-antihermitianpart(A) = antihermitianpart!(copy(A))
-function antihermitianpart!(A::AbstractMatrix)
-    Base.require_one_based_indexing(A)
-    n = LinearAlgebra.checksquare(A)
-    @inbounds for j in 1:n
-        A[j, j] = imag(A[j, j]) * im
-        for i in 1:(j - 1)
-            val = (A[i, j] - adjoint(A[j, i])) / 2
-            A[i, j] = val
-            A[j, i] = -adjoint(val)
-        end
-    end
-    return A
-end
diff --git a/src/implementations/hermitian.jl b/src/implementations/hermitian.jl
new file mode 100644
index 00000000..048c0768
--- /dev/null
+++ b/src/implementations/hermitian.jl
@@ -0,0 +1,95 @@
+# Inputs
+# ------
+function copy_input(::typeof(hermitianpart), A::AbstractMatrix)
+    return copy!(similar(A, float(eltype(A))), A)
+end
+copy_input(::typeof(antihermitianpart), A) = copy_input(hermitianpart, A)
+
+function check_input(::typeof(hermitianpart!), A::AbstractMatrix, B::AbstractMatrix, ::AbstractAlgorithm)
+    LinearAlgebra.checksquare(A)
+    n = Base.require_one_based_indexing(A)
+    B === A || @check_size(B, (n, n))
+    return nothing
+end
+function check_input(::typeof(antihermitianpart!), A::AbstractMatrix, B::AbstractMatrix, ::AbstractAlgorithm)
+    LinearAlgebra.checksquare(A)
+    n = Base.require_one_based_indexing(A)
+    B === A || @check_size(B, (n, n))
+    return nothing
+end
+
+# Outputs
+# -------
+function initialize_output(::typeof(hermitianpart!), A::AbstractMatrix, ::NativeBlocked)
+    return A
+end
+function initialize_output(::typeof(antihermitianpart!), A::AbstractMatrix, ::NativeBlocked)
+    return A
+end
+
+# Implementation
+# --------------
+function hermitianpart!(A::AbstractMatrix, B, alg::NativeBlocked)
+    check_input(hermitianpart!, A, B, alg)
+    return hermitianpart_native!(A, B, Val(false); alg.kwargs...)
+end
+function antihermitianpart!(A::AbstractMatrix, B, alg::NativeBlocked)
+    check_input(antihermitianpart!, A, B, alg)
+    return hermitianpart_native!(A, B, Val(true); alg.kwargs...)
+end
+
+function hermitianpart_native!(A::AbstractMatrix, B::AbstractMatrix, anti::Val; blocksize = 32)
+    n = size(A, 1)
+    for j in 1:blocksize:n
+        for i in 1:blocksize:(j - 1)
+            jb = min(blocksize, n - j + 1)
+            ib = blocksize
+            _hermitianpart_offdiag!(
+                view(A, i:(i + ib - 1), j:(j + jb - 1)),
+                view(A, j:(j + jb - 1), i:(i + ib - 1)),
+                view(B, i:(i + ib - 1), j:(j + jb - 1)),
+                view(B, j:(j + jb - 1), i:(i + ib - 1)),
+                anti
+            )
+        end
+        jb = min(blocksize, n - j + 1)
+        _hermitianpart_diag!(
+            view(A, j:(j + jb - 1), j:(j + jb - 1)),
+            view(B, j:(j + jb - 1), j:(j + jb - 1)),
+            anti
+        )
+    end
+    return B
+end
+
+function _hermitianpart_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
+            Bu[i, j] = val
+            aval = adjoint(val)
+            Bl[j, i] = anti ? -aval : aval
+        end
+    end
+    return nothing
+end
+function _hermitianpart_diag!(A::AbstractMatrix, B::AbstractMatrix, ::Val{anti}) where {anti}
+    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
+            B[i, j] = val
+            aval = adjoint(val)
+            B[j, i] = anti ? -aval : aval
+        end
+        B[j, j] = anti ? _imimag(A[j, j]) : real(A[j, j])
+    end
+    return nothing
+end
+
+_imimag(x::Real) = zero(x)
+_imimag(x::Complex) = im * imag(x)
diff --git a/src/interface/hermitian.jl b/src/interface/hermitian.jl
new file mode 100644
index 00000000..704fb974
--- /dev/null
+++ b/src/interface/hermitian.jl
@@ -0,0 +1,45 @@
+@doc """
+    hermitianpart(A; kwargs...)
+    hermitianpart(A, alg)
+    hermitianpart!(A; kwargs...)
+    hermitianpart!(A, alg)
+
+Compute the hermitian part of a (square) matrix `A`, defined as `(A + A') / 2`.
+For real matrices this corresponds to the symmetric part of `A`.
+
+See also [`antihermitianpart`](@ref).
+"""
+@functiondef hermitianpart
+
+@doc """
+    antihermitianpart(A; kwargs...)
+    antihermitianpart(A, alg)
+    antihermitianpart!(A; kwargs...)
+    antihermitianpart!(A, alg)
+
+Compute the anti-hermitian part of a (square) matrix `A`, defined as `(A - A') / 2`.
+For real matrices this corresponds to the antisymmetric part of `A`.
+
+See also [`hermitianpart`](@ref).
+"""
+@functiondef antihermitianpart
+
+"""
+NativeBlocked(; blocksize = 32)
+
+Algorithm type to denote a native blocked algorithm with given `blocksize` for computing
+the hermitian or anti-hermitian part of a matrix.
+"""
+@algdef NativeBlocked
+# TODO: multithreaded? numthreads keyword?
+
+default_hermitian_algorithm(A; kwargs...) = default_hermitian_algorithm(typeof(A); kwargs...)
+function default_hermitian_algorithm(::Type{A}; kwargs...) where {A <: AbstractMatrix}
+    return NativeBlocked(; kwargs...)
+end
+
+for f in (:hermitianpart!, :antihermitianpart!)
+    @eval function default_algorithm(::typeof($f), ::Type{A}; kwargs...) where {A}
+        return default_hermitian_algorithm(A; kwargs...)
+    end
+end
diff --git a/test/projections.jl b/test/projections.jl
new file mode 100644
index 00000000..7d95719d
--- /dev/null
+++ b/test/projections.jl
@@ -0,0 +1,39 @@
+using MatrixAlgebraKit
+using Test
+using TestExtras
+using StableRNGs
+using LinearAlgebra: Diagonal
+
+const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64)
+
+@testset "(anti)hermitianpart! for T = $T" for T in BLASFloats
+    rng = StableRNG(123)
+    m = 54
+    for alg in (NativeBlocked(blocksize = 16), NativeBlocked(blocksize = 32), NativeBlocked(blocksize = 64))
+        A = randn(rng, T, m, m)
+        Ah = (A + A') / 2
+        Aa = (A - A') / 2
+        Ac = copy(A)
+
+        Bh = hermitianpart(A, alg)
+        @test ishermitian(Bh)
+        @test Bh ≈ Ah
+        @test A == Ac
+
+        Ba = antihermitianpart(A, alg)
+        @test isantihermitian(Ba)
+        @test Ba ≈ Aa
+        @test A == Ac
+
+        Bh = hermitianpart!(Ac, alg)
+        @test Bh === Ac
+        @test ishermitian(Bh)
+        @test Bh ≈ Ah
+
+        copy!(Ac, A)
+        Ba = antihermitianpart!(Ac, alg)
+        @test Ba === Ac
+        @test isantihermitian(Ba)
+        @test Ba ≈ Aa
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 207215ec..10f82214 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,6 +7,9 @@ if !is_buildkite
     @safetestset "Algorithms" begin
         include("algorithms.jl")
     end
+    @safetestset "Projections" begin
+        include("projections.jl")
+    end
     @safetestset "Truncate" begin
         include("truncate.jl")
     end

From 90db13fc29394a591d9cc27fc647cc3cf40e492c Mon Sep 17 00:00:00 2001
From: Jutho Haegeman <jutho.haegeman@ugent.be>
Date: Sun, 5 Oct 2025 21:50:02 +0200
Subject: [PATCH 05/13] rename to project_

---
 src/MatrixAlgebraKit.jl                       |  4 +--
 .../{hermitian.jl => projections.jl}          | 34 +++++++++----------
 .../{hermitian.jl => projections.jl}          | 26 +++++++-------
 src/pullbacks/eigh.jl                         |  4 +--
 src/pullbacks/lq.jl                           |  2 +-
 src/pullbacks/qr.jl                           |  2 +-
 src/pullbacks/svd.jl                          |  8 ++---
 test/projections.jl                           | 10 +++---
 8 files changed, 45 insertions(+), 45 deletions(-)
 rename src/implementations/{hermitian.jl => projections.jl} (63%)
 rename src/interface/{hermitian.jl => projections.jl} (66%)

diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl
index af66d7c6..8b09f1ec 100644
--- a/src/MatrixAlgebraKit.jl
+++ b/src/MatrixAlgebraKit.jl
@@ -11,8 +11,8 @@ using LinearAlgebra: BlasFloat, BlasReal, BlasComplex, BlasInt
 
 export isisometry, isunitary, ishermitian, isantihermitian
 
-export hermitianpart, antihermitianpart
-export hermitianpart!, antihermitianpart!
+export project_hermitian, project_antihermitian
+export project_hermitian!, project_antihermitian!
 export qr_compact, qr_full, qr_null, lq_compact, lq_full, lq_null
 export qr_compact!, qr_full!, qr_null!, lq_compact!, lq_full!, lq_null!
 export svd_compact, svd_full, svd_vals, svd_trunc
diff --git a/src/implementations/hermitian.jl b/src/implementations/projections.jl
similarity index 63%
rename from src/implementations/hermitian.jl
rename to src/implementations/projections.jl
index 048c0768..c905180b 100644
--- a/src/implementations/hermitian.jl
+++ b/src/implementations/projections.jl
@@ -1,17 +1,17 @@
 # Inputs
 # ------
-function copy_input(::typeof(hermitianpart), A::AbstractMatrix)
+function copy_input(::typeof(project_hermitian), A::AbstractMatrix)
     return copy!(similar(A, float(eltype(A))), A)
 end
-copy_input(::typeof(antihermitianpart), A) = copy_input(hermitianpart, A)
+copy_input(::typeof(project_antihermitian), A) = copy_input(project_hermitian, A)
 
-function check_input(::typeof(hermitianpart!), A::AbstractMatrix, B::AbstractMatrix, ::AbstractAlgorithm)
+function check_input(::typeof(project_hermitian!), A::AbstractMatrix, B::AbstractMatrix, ::AbstractAlgorithm)
     LinearAlgebra.checksquare(A)
     n = Base.require_one_based_indexing(A)
     B === A || @check_size(B, (n, n))
     return nothing
 end
-function check_input(::typeof(antihermitianpart!), A::AbstractMatrix, B::AbstractMatrix, ::AbstractAlgorithm)
+function check_input(::typeof(project_antihermitian!), A::AbstractMatrix, B::AbstractMatrix, ::AbstractAlgorithm)
     LinearAlgebra.checksquare(A)
     n = Base.require_one_based_indexing(A)
     B === A || @check_size(B, (n, n))
@@ -20,31 +20,31 @@ end
 
 # Outputs
 # -------
-function initialize_output(::typeof(hermitianpart!), A::AbstractMatrix, ::NativeBlocked)
+function initialize_output(::typeof(project_hermitian!), A::AbstractMatrix, ::NativeBlocked)
     return A
 end
-function initialize_output(::typeof(antihermitianpart!), A::AbstractMatrix, ::NativeBlocked)
+function initialize_output(::typeof(project_antihermitian!), A::AbstractMatrix, ::NativeBlocked)
     return A
 end
 
 # Implementation
 # --------------
-function hermitianpart!(A::AbstractMatrix, B, alg::NativeBlocked)
-    check_input(hermitianpart!, A, B, alg)
-    return hermitianpart_native!(A, B, Val(false); alg.kwargs...)
+function project_hermitian!(A::AbstractMatrix, B, alg::NativeBlocked)
+    check_input(project_hermitian!, A, B, alg)
+    return project_hermitian_native!(A, B, Val(false); alg.kwargs...)
 end
-function antihermitianpart!(A::AbstractMatrix, B, alg::NativeBlocked)
-    check_input(antihermitianpart!, A, B, alg)
-    return hermitianpart_native!(A, B, Val(true); alg.kwargs...)
+function project_antihermitian!(A::AbstractMatrix, B, alg::NativeBlocked)
+    check_input(project_antihermitian!, A, B, alg)
+    return project_hermitian_native!(A, B, Val(true); alg.kwargs...)
 end
 
-function hermitianpart_native!(A::AbstractMatrix, B::AbstractMatrix, anti::Val; blocksize = 32)
+function project_hermitian_native!(A::AbstractMatrix, B::AbstractMatrix, anti::Val; blocksize = 32)
     n = size(A, 1)
     for j in 1:blocksize:n
         for i in 1:blocksize:(j - 1)
             jb = min(blocksize, n - j + 1)
             ib = blocksize
-            _hermitianpart_offdiag!(
+            _project_hermitian_offdiag!(
                 view(A, i:(i + ib - 1), j:(j + jb - 1)),
                 view(A, j:(j + jb - 1), i:(i + ib - 1)),
                 view(B, i:(i + ib - 1), j:(j + jb - 1)),
@@ -53,7 +53,7 @@ function hermitianpart_native!(A::AbstractMatrix, B::AbstractMatrix, anti::Val;
             )
         end
         jb = min(blocksize, n - j + 1)
-        _hermitianpart_diag!(
+        _project_hermitian_diag!(
             view(A, j:(j + jb - 1), j:(j + jb - 1)),
             view(B, j:(j + jb - 1), j:(j + jb - 1)),
             anti
@@ -62,7 +62,7 @@ function hermitianpart_native!(A::AbstractMatrix, B::AbstractMatrix, anti::Val;
     return B
 end
 
-function _hermitianpart_offdiag!(
+function _project_hermitian_offdiag!(
         Au::AbstractMatrix, Al::AbstractMatrix, Bu::AbstractMatrix, Bl::AbstractMatrix, ::Val{anti}
     ) where {anti}
 
@@ -77,7 +77,7 @@ function _hermitianpart_offdiag!(
     end
     return nothing
 end
-function _hermitianpart_diag!(A::AbstractMatrix, B::AbstractMatrix, ::Val{anti}) where {anti}
+function _project_hermitian_diag!(A::AbstractMatrix, B::AbstractMatrix, ::Val{anti}) where {anti}
     n = size(A, 1)
     @inbounds for j in 1:n
         @simd for i in 1:(j - 1)
diff --git a/src/interface/hermitian.jl b/src/interface/projections.jl
similarity index 66%
rename from src/interface/hermitian.jl
rename to src/interface/projections.jl
index 704fb974..dbf62b92 100644
--- a/src/interface/hermitian.jl
+++ b/src/interface/projections.jl
@@ -1,28 +1,28 @@
 @doc """
-    hermitianpart(A; kwargs...)
-    hermitianpart(A, alg)
-    hermitianpart!(A; kwargs...)
-    hermitianpart!(A, alg)
+    project_hermitian(A; kwargs...)
+    project_hermitian(A, alg)
+    project_hermitian!(A; kwargs...)
+    project_hermitian!(A, alg)
 
 Compute the hermitian part of a (square) matrix `A`, defined as `(A + A') / 2`.
 For real matrices this corresponds to the symmetric part of `A`.
 
-See also [`antihermitianpart`](@ref).
+See also [`project_antihermitian`](@ref).
 """
-@functiondef hermitianpart
+@functiondef project_hermitian
 
 @doc """
-    antihermitianpart(A; kwargs...)
-    antihermitianpart(A, alg)
-    antihermitianpart!(A; kwargs...)
-    antihermitianpart!(A, alg)
+    project_antihermitian(A; kwargs...)
+    project_antihermitian(A, alg)
+    project_antihermitian!(A; kwargs...)
+    project_antihermitian!(A, alg)
 
 Compute the anti-hermitian part of a (square) matrix `A`, defined as `(A - A') / 2`.
 For real matrices this corresponds to the antisymmetric part of `A`.
 
-See also [`hermitianpart`](@ref).
+See also [`project_hermitian`](@ref).
 """
-@functiondef antihermitianpart
+@functiondef project_antihermitian
 
 """
 NativeBlocked(; blocksize = 32)
@@ -38,7 +38,7 @@ function default_hermitian_algorithm(::Type{A}; kwargs...) where {A <: AbstractM
     return NativeBlocked(; kwargs...)
 end
 
-for f in (:hermitianpart!, :antihermitianpart!)
+for f in (:project_hermitian!, :project_antihermitian!)
     @eval function default_algorithm(::typeof($f), ::Type{A}; kwargs...) where {A}
         return default_hermitian_algorithm(A; kwargs...)
     end
diff --git a/src/pullbacks/eigh.jl b/src/pullbacks/eigh.jl
index cb000a06..b15f912e 100644
--- a/src/pullbacks/eigh.jl
+++ b/src/pullbacks/eigh.jl
@@ -42,7 +42,7 @@ function eigh_pullback!(
         indV = axes(V, 2)[ind]
         length(indV) == pV || throw(DimensionMismatch())
         mul!(view(VᴴΔV, :, indV), V', ΔV)
-        aVᴴΔV = antihermitianpart(VᴴΔV) # can't use in-place or recycling doesn't work
+        aVᴴΔV = project_antihermitian(VᴴΔV) # can't use in-place or recycling doesn't work
 
         mask = abs.(D' .- D) .< degeneracy_atol
         Δgauge = norm(view(aVᴴΔV, mask))
@@ -112,7 +112,7 @@ function eigh_trunc_pullback!(
     if !iszerotangent(ΔV)
         (n, p) == size(ΔV) || throw(DimensionMismatch())
         VᴴΔV = V' * ΔV
-        aVᴴΔV = antihermitianpart!(VᴴΔV)
+        aVᴴΔV = project_antihermitian!(VᴴΔV)
 
         mask = abs.(D' .- D) .< degeneracy_atol
         Δgauge = norm(view(aVᴴΔV, mask))
diff --git a/src/pullbacks/lq.jl b/src/pullbacks/lq.jl
index b130a6d7..05e23c38 100644
--- a/src/pullbacks/lq.jl
+++ b/src/pullbacks/lq.jl
@@ -118,7 +118,7 @@ function lq_null_pullback!(
         gauge_atol::Real = tol
     )
     if !iszerotangent(ΔNᴴ) && size(Nᴴ, 1) > 0
-        aNᴴΔN = antihermitianpart!(Nᴴ * ΔNᴴ')
+        aNᴴΔN = project_antihermitian!(Nᴴ * ΔNᴴ')
         Δgauge = norm(aNᴴΔN)
         Δgauge < tol ||
             @warn "`lq_null` cotangent sensitive to gauge choice: (|Δgauge| = $Δgauge)"
diff --git a/src/pullbacks/qr.jl b/src/pullbacks/qr.jl
index 518dcd27..bc49d6af 100644
--- a/src/pullbacks/qr.jl
+++ b/src/pullbacks/qr.jl
@@ -117,7 +117,7 @@ function qr_null_pullback!(
         gauge_atol::Real = tol
     )
     if !iszerotangent(ΔN) && size(N, 2) > 0
-        aNᴴΔN = antihermitianpart!(N' * ΔN)
+        aNᴴΔN = project_antihermitian!(N' * ΔN)
         Δgauge = norm(aNᴴΔN)
         Δgauge < tol ||
             @warn "`qr_null` cotangent sensitive to gauge choice: (|Δgauge| = $Δgauge)"
diff --git a/src/pullbacks/svd.jl b/src/pullbacks/svd.jl
index 40c61884..effceaa3 100644
--- a/src/pullbacks/svd.jl
+++ b/src/pullbacks/svd.jl
@@ -69,8 +69,8 @@ function svd_pullback!(
     end
 
     # Project onto antihermitian part; hermitian part outside of Grassmann tangent space
-    aUΔU = antihermitianpart!(UΔU)
-    aVΔV = antihermitianpart!(VΔV)
+    aUΔU = project_antihermitian!(UΔU)
+    aVΔV = project_antihermitian!(VΔV)
 
     # check whether cotangents arise from gauge-invariance objective function
     mask = abs.(Sr' .- Sr) .< degeneracy_atol
@@ -159,8 +159,8 @@ function svd_trunc_pullback!(
     end
 
     # Project onto antihermitian part; hermitian part outside of Grassmann tangent space
-    aUΔU = antihermitianpart!(UΔU)
-    aVΔV = antihermitianpart!(VΔV)
+    aUΔU = project_antihermitian!(UΔU)
+    aVΔV = project_antihermitian!(VΔV)
 
     # check whether cotangents arise from gauge-invariance objective function
     mask = abs.(S' .- S) .< degeneracy_atol
diff --git a/test/projections.jl b/test/projections.jl
index 7d95719d..089965f7 100644
--- a/test/projections.jl
+++ b/test/projections.jl
@@ -6,7 +6,7 @@ using LinearAlgebra: Diagonal
 
 const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64)
 
-@testset "(anti)hermitianpart! for T = $T" for T in BLASFloats
+@testset "project_(anti)hermitian! for T = $T" for T in BLASFloats
     rng = StableRNG(123)
     m = 54
     for alg in (NativeBlocked(blocksize = 16), NativeBlocked(blocksize = 32), NativeBlocked(blocksize = 64))
@@ -15,23 +15,23 @@ const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64)
         Aa = (A - A') / 2
         Ac = copy(A)
 
-        Bh = hermitianpart(A, alg)
+        Bh = project_hermitian(A, alg)
         @test ishermitian(Bh)
         @test Bh ≈ Ah
         @test A == Ac
 
-        Ba = antihermitianpart(A, alg)
+        Ba = project_antihermitian(A, alg)
         @test isantihermitian(Ba)
         @test Ba ≈ Aa
         @test A == Ac
 
-        Bh = hermitianpart!(Ac, alg)
+        Bh = project_hermitian!(Ac, alg)
         @test Bh === Ac
         @test ishermitian(Bh)
         @test Bh ≈ Ah
 
         copy!(Ac, A)
-        Ba = antihermitianpart!(Ac, alg)
+        Ba = project_antihermitian!(Ac, alg)
         @test Ba === Ac
         @test isantihermitian(Ba)
         @test Ba ≈ Aa

From 42a55650f3e04695dba7b15cb4c853db32cf7a76 Mon Sep 17 00:00:00 2001
From: Jutho Haegeman <jutho.haegeman@ugent.be>
Date: Mon, 6 Oct 2025 09:16:01 +0200
Subject: [PATCH 06/13] also change includes

---
 src/MatrixAlgebraKit.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl
index 8b09f1ec..276523c9 100644
--- a/src/MatrixAlgebraKit.jl
+++ b/src/MatrixAlgebraKit.jl
@@ -3,8 +3,8 @@ module MatrixAlgebraKit
 using LinearAlgebra: LinearAlgebra
 using LinearAlgebra: norm # TODO: eleminate if we use VectorInterface.jl?
 using LinearAlgebra: mul!, rmul!, lmul!, adjoint!, rdiv!, ldiv!
-using LinearAlgebra: sylvester
-using LinearAlgebra: isposdef, issymmetric
+using LinearAlgebra: sylvester, lu!
+using LinearAlgebra: isposdef, issymmetric, opnorm
 using LinearAlgebra: Diagonal, diag, diagind, isdiag
 using LinearAlgebra: UpperTriangular, LowerTriangular
 using LinearAlgebra: BlasFloat, BlasReal, BlasComplex, BlasInt
@@ -77,7 +77,7 @@ include("common/gauge.jl")
 
 include("yalapack.jl")
 include("algorithms.jl")
-include("interface/hermitian.jl")
+include("interface/projections.jl")
 include("interface/decompositions.jl")
 include("interface/truncation.jl")
 include("interface/qr.jl")
@@ -90,7 +90,7 @@ include("interface/schur.jl")
 include("interface/polar.jl")
 include("interface/orthnull.jl")
 
-include("implementations/hermitian.jl")
+include("implementations/projections.jl")
 include("implementations/truncation.jl")
 include("implementations/qr.jl")
 include("implementations/lq.jl")

From e22af665fac75bf0818995a975659f51c902b9bc Mon Sep 17 00:00:00 2001
From: Jutho Haegeman <jutho.haegeman@ugent.be>
Date: Mon, 6 Oct 2025 09:49:21 +0200
Subject: [PATCH 07/13] fix polar pullback

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

diff --git a/src/pullbacks/polar.jl b/src/pullbacks/polar.jl
index bc5c3877..fabc2c2e 100644
--- a/src/pullbacks/polar.jl
+++ b/src/pullbacks/polar.jl
@@ -11,7 +11,7 @@ function left_polar_pullback!(ΔA::AbstractMatrix, A, WP, ΔWP)
     # Extract and check the cotangents
     ΔW, ΔP = ΔWP
     if !iszerotangent(ΔP)
-        ΔP = hermitianpart(ΔP)
+        ΔP = project_hermitian(ΔP)
     end
     M = zero(P)
     !iszerotangent(ΔW) && mul!(M, W', ΔW, 1, 1)
@@ -41,7 +41,7 @@ function right_polar_pullback!(ΔA::AbstractMatrix, A, PWᴴ, ΔPWᴴ)
     # Extract and check the cotangents
     ΔP, ΔWᴴ = ΔPWᴴ
     if !iszerotangent(ΔP)
-        ΔP = hermitianpart(ΔP)
+        ΔP = project_hermitian(ΔP)
     end
     M = zero(P)
     !iszerotangent(ΔWᴴ) && mul!(M, ΔWᴴ, Wᴴ', 1, 1)

From 39241965859fc8be97f4fbf64d0cd145f1a33066 Mon Sep 17 00:00:00 2001
From: Jutho Haegeman <jutho.haegeman@ugent.be>
Date: Mon, 6 Oct 2025 10:45:39 +0200
Subject: [PATCH 08/13] add blocked ishermitian

---
 src/common/matrixproperties.jl | 95 +++++++++++++++++++---------------
 1 file changed, 54 insertions(+), 41 deletions(-)

diff --git a/src/common/matrixproperties.jl b/src/common/matrixproperties.jl
index 9d0c3771..4741d620 100644
--- a/src/common/matrixproperties.jl
+++ b/src/common/matrixproperties.jl
@@ -31,6 +31,7 @@ function isunitary(A; isapprox_kwargs...)
     return is_left_isometry(A; isapprox_kwargs...) &&
         is_right_isometry(A; isapprox_kwargs...)
 end
+# TODO: for matrices, isunitary corresponds to is_left_isometry + square.
 
 @doc """
     is_left_isometry(A; isapprox_kwargs...) -> Bool
@@ -64,30 +65,18 @@ end
 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, kwargs...)
-    return (atol == rtol == 0) ? (A == A') : isapprox(A, A'; atol, rtol, kwargs...)
-end
-function ishermitian(A::AbstractMatrix; atol::Real = 0, rtol::Real = 0, norm = LinearAlgebra.norm, kwargs...)
-    Base.require_one_based_indexing(A)
-    m, n = size(A)
-    m == n || return false
-    if atol == rtol == 0
-        for j in 1:n
-            for i in 1:j
-                A[i, j] == adjoint(A[j, i]) || return false
-            end
-        end
-    elseif norm === LinearAlgebra.norm
-        atol = max(atol, rtol * norm(A))
-        for j in 1:n
-            for i in 1:j
-                isapprox(A[i, j], adjoint(A[j, i]); atol, abs, kwargs...)  || return false
-            end
-        end
+function ishermitian(A; atol::Real = 0, rtol::Real = 0, norm = LinearAlgebra.norm, kwargs...)
+    if iszero(atol) && iszero(rtol)
+        return ishermitian_exact(A; kwargs...)
     else
-        return isapprox(A, A'; atol, rtol, norm, kwargs...)
+        return 2 * norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
     end
-    return true
+end
+function ishermitian_exact(A)
+    return A == A'
+end
+function ishermitian_exact(A::AbstractMatrix; kwargs...)
+    return _ishermitian_exact(A, Val(false); kwargs...)
 end
 
 """
@@ -96,28 +85,52 @@ 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, kwargs...)
-    return (atol == 0 & rtol == 0) ? (A == -A') : isapprox(A, -A'; atol, rtol, kwargs...)
+function isantihermitian(A; atol::Real = 0, rtol::Real = 0, norm = LinearAlgebra.norm, kwargs...)
+    if iszero(atol) && iszero(rtol)
+        return isantihermitian_exact(A; kwargs...)
+    else
+        return 2 * norm(project_hermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
+    end
+end
+function isantihermitian_exact(A)
+    return A == -A'
+end
+function isantihermitian_exact(A::AbstractMatrix; kwargs...)
+    return _ishermitian_exact(A, Val(true); kwargs...)
 end
-function isantihermitian(A::AbstractMatrix; atol::Real = 0, rtol::Real = 0, norm = LinearAlgebra.norm, kwargs...)
-    Base.require_one_based_indexing(A)
-    m, n = size(A)
-    m == n || return false
-    if atol == rtol == 0
-        @inbounds for j in 1:n
-            for i in 1:j
-                A[i, j] == -adjoint(A[j, i]) || return false
-            end
+
+# block implementation of exact checks
+function _ishermitian_exact(A::AbstractMatrix, anti::Val; blocksize = 32)
+    n = size(A, 1)
+    for j in 1:blocksize:n
+        jb = min(blocksize, n - j + 1)
+        _ishermitian_exact_diag(view(A, j:(j + jb - 1), j:(j + jb - 1)), anti) || return false
+        for i in 1:blocksize:(j - 1)
+            ib = blocksize
+            _ishermitian_exact_offdiag(
+                view(A, i:(i + ib - 1), j:(j + jb - 1)),
+                view(A, j:(j + jb - 1), i:(i + ib - 1)),
+                anti
+            ) || return false
         end
-    elseif norm === LinearAlgebra.norm
-        atol = max(atol, rtol * norm(A))
-        @inbounds for j in 1:n
-            for i in 1:j
-                isapprox(A[i, j], -adjoint(A[j, i]); atol, abs, kwargs...)  || return false
-            end
+    end
+    return true
+end
+function _ishermitian_exact_diag(A, ::Val{anti}) where {anti}
+    n = size(A, 1)
+    @inbounds for j in 1:n
+        @simd for i in 1:j
+            A[i, j] == (anti ? -adjoint(A[j, i]) : adjoint(A[j, i])) || return false
+        end
+    end
+    return true
+end
+function _ishermitian_exact_offdiag(Al, Au, ::Val{anti}) where {anti}
+    m, n = size(Al) # == reverse(size(Al))
+    @inbounds for j in 1:n
+        @simd for i in 1:m
+            Al[i, j] == (anti ? -adjoint(Au[j, i]) : adjoint(Au[j, i])) || return false
         end
-    else
-        return isapprox(A, -A'; atol, rtol, norm, kwargs...)
     end
     return true
 end

From f4ea67cf6e6fcb2e3c8ce7cded634a39584883c6 Mon Sep 17 00:00:00 2001
From: Jutho Haegeman <jutho.haegeman@ugent.be>
Date: Mon, 6 Oct 2025 16:49:31 +0200
Subject: [PATCH 09/13] improve matrixproperties

---
 src/common/matrixproperties.jl | 26 ++++++++++++++++++--------
 1 file changed, 18 insertions(+), 8 deletions(-)

diff --git a/src/common/matrixproperties.jl b/src/common/matrixproperties.jl
index 4741d620..94871f0a 100644
--- a/src/common/matrixproperties.jl
+++ b/src/common/matrixproperties.jl
@@ -31,7 +31,10 @@ function isunitary(A; isapprox_kwargs...)
     return is_left_isometry(A; isapprox_kwargs...) &&
         is_right_isometry(A; isapprox_kwargs...)
 end
-# TODO: for matrices, isunitary corresponds to is_left_isometry + square.
+function isunitary(A::AbstractMatrix; isapprox_kwargs...)
+    size(A, 1) == size(A, 2) || return false
+    return is_left_isometry(A; isapprox_kwargs...)
+end
 
 @doc """
     is_left_isometry(A; isapprox_kwargs...) -> Bool
@@ -42,8 +45,11 @@ The `isapprox_kwargs` can be used to control the tolerances of the equality.
 See also [`isisometry`](@ref) and [`is_right_isometry`](@ref).
 """ is_left_isometry
 
-function is_left_isometry(A::AbstractMatrix; isapprox_kwargs...)
-    return isapprox(A' * A, LinearAlgebra.I; isapprox_kwargs...)
+function is_left_isometry(A::AbstractMatrix; atol::Real = 0, rtol::Real = defaulttol(A), norm = LinearAlgebra.norm)
+    P = A' * A
+    nP = norm(P) # isapprox would use `rtol * max(norm(P), norm(I))`
+    diagview(P) .-= 1
+    return norm(P) <= max(atol, rtol * nP) # assume that the norm of I is `sqrt(n)`
 end
 
 @doc """
@@ -55,8 +61,11 @@ The `isapprox_kwargs` can be used to control the tolerances of the equality.
 See also [`isisometry`](@ref) and [`is_left_isometry`](@ref).
 """ is_right_isometry
 
-function is_right_isometry(A::AbstractMatrix; isapprox_kwargs...)
-    return isapprox(A * A', LinearAlgebra.I; isapprox_kwargs...)
+function is_right_isometry(A::AbstractMatrix; atol::Real = 0, rtol::Real = defaulttol(A), norm = LinearAlgebra.norm)
+    P = A * A'
+    nP = norm(P) # isapprox would use `rtol * max(norm(P), norm(I))`
+    diagview(P) .-= 1
+    return norm(P) <= max(atol, rtol * nP) # assume that the norm of I is `sqrt(n)`
 end
 
 """
@@ -75,7 +84,7 @@ end
 function ishermitian_exact(A)
     return A == A'
 end
-function ishermitian_exact(A::AbstractMatrix; kwargs...)
+function ishermitian_exact(A::StridedMatrix; kwargs...)
     return _ishermitian_exact(A, Val(false); kwargs...)
 end
 
@@ -95,11 +104,12 @@ end
 function isantihermitian_exact(A)
     return A == -A'
 end
-function isantihermitian_exact(A::AbstractMatrix; kwargs...)
+function isantihermitian_exact(A::StridedMatrix; kwargs...)
     return _ishermitian_exact(A, Val(true); kwargs...)
 end
 
-# block implementation of exact checks
+# blocked implementation of exact checks for strided matrices
+# -----------------------------------------------------------
 function _ishermitian_exact(A::AbstractMatrix, anti::Val; blocksize = 32)
     n = size(A, 1)
     for j in 1:blocksize:n

From 9e9df93d9cc293a4edcd0d40bcf21b97f4c1df90 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Mon, 6 Oct 2025 11:34:56 -0400
Subject: [PATCH 10/13] increase coverage

---
 test/projections.jl | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/test/projections.jl b/test/projections.jl
index 089965f7..2eb2e28e 100644
--- a/test/projections.jl
+++ b/test/projections.jl
@@ -9,6 +9,7 @@ const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64)
 @testset "project_(anti)hermitian! for T = $T" for T in BLASFloats
     rng = StableRNG(123)
     m = 54
+    noisefactor = eps(real(T))^(3 / 4)
     for alg in (NativeBlocked(blocksize = 16), NativeBlocked(blocksize = 32), NativeBlocked(blocksize = 64))
         A = randn(rng, T, m, m)
         Ah = (A + A') / 2
@@ -19,11 +20,17 @@ const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64)
         @test ishermitian(Bh)
         @test Bh ≈ Ah
         @test A == Ac
+        Bh_approx = Bh + noisefactor * Aa
+        @test !ishermitian(Bh_approx)
+        @test ishermitian(Bh_approx; rtol = 10 * noisefactor)
 
         Ba = project_antihermitian(A, alg)
         @test isantihermitian(Ba)
         @test Ba ≈ Aa
         @test A == Ac
+        Ba_approx = Ba + noisefactor * Ah
+        @test !isantihermitian(Ba_approx)
+        @test isantihermitian(Ba_approx; rtol = 10 * noisefactor)
 
         Bh = project_hermitian!(Ac, alg)
         @test Bh === Ac

From 9f7617a530d998098c473582cefa88a08ebb3302 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Mon, 6 Oct 2025 18:16:52 -0400
Subject: [PATCH 11/13] Update src/MatrixAlgebraKit.jl

Co-authored-by: Jutho <Jutho@users.noreply.github.com>
---
 src/MatrixAlgebraKit.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl
index 276523c9..050f7008 100644
--- a/src/MatrixAlgebraKit.jl
+++ b/src/MatrixAlgebraKit.jl
@@ -3,8 +3,8 @@ module MatrixAlgebraKit
 using LinearAlgebra: LinearAlgebra
 using LinearAlgebra: norm # TODO: eleminate if we use VectorInterface.jl?
 using LinearAlgebra: mul!, rmul!, lmul!, adjoint!, rdiv!, ldiv!
-using LinearAlgebra: sylvester, lu!
-using LinearAlgebra: isposdef, issymmetric, opnorm
+using LinearAlgebra: sylvester
+using LinearAlgebra: isposdef
 using LinearAlgebra: Diagonal, diag, diagind, isdiag
 using LinearAlgebra: UpperTriangular, LowerTriangular
 using LinearAlgebra: BlasFloat, BlasReal, BlasComplex, BlasInt

From bb077a127c4be6bf11c32b5ed4c94c6f76d21671 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Mon, 6 Oct 2025 18:58:29 -0400
Subject: [PATCH 12/13] small rename to strided_ishermitian_exact

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

diff --git a/src/common/matrixproperties.jl b/src/common/matrixproperties.jl
index 94871f0a..61ff73e9 100644
--- a/src/common/matrixproperties.jl
+++ b/src/common/matrixproperties.jl
@@ -85,7 +85,7 @@ function ishermitian_exact(A)
     return A == A'
 end
 function ishermitian_exact(A::StridedMatrix; kwargs...)
-    return _ishermitian_exact(A, Val(false); kwargs...)
+    return strided_ishermitian_exact(A, Val(false); kwargs...)
 end
 
 """
@@ -105,12 +105,12 @@ function isantihermitian_exact(A)
     return A == -A'
 end
 function isantihermitian_exact(A::StridedMatrix; kwargs...)
-    return _ishermitian_exact(A, Val(true); kwargs...)
+    return strided_ishermitian_exact(A, Val(true); kwargs...)
 end
 
 # blocked implementation of exact checks for strided matrices
 # -----------------------------------------------------------
-function _ishermitian_exact(A::AbstractMatrix, anti::Val; blocksize = 32)
+function strided_ishermitian_exact(A::AbstractMatrix, anti::Val; blocksize = 32)
     n = size(A, 1)
     for j in 1:blocksize:n
         jb = min(blocksize, n - j + 1)

From 2af93e899bcb64a76e417d7272af7e2f6c8c1d3b Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Mon, 6 Oct 2025 18:59:41 -0400
Subject: [PATCH 13/13] add back missing import

---
 src/MatrixAlgebraKit.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl
index 050f7008..69ee4171 100644
--- a/src/MatrixAlgebraKit.jl
+++ b/src/MatrixAlgebraKit.jl
@@ -4,7 +4,7 @@ using LinearAlgebra: LinearAlgebra
 using LinearAlgebra: norm # TODO: eleminate if we use VectorInterface.jl?
 using LinearAlgebra: mul!, rmul!, lmul!, adjoint!, rdiv!, ldiv!
 using LinearAlgebra: sylvester
-using LinearAlgebra: isposdef
+using LinearAlgebra: isposdef, issymmetric
 using LinearAlgebra: Diagonal, diag, diagind, isdiag
 using LinearAlgebra: UpperTriangular, LowerTriangular
 using LinearAlgebra: BlasFloat, BlasReal, BlasComplex, BlasInt