From 0a8927850bff2022e6d063d52d0d92224f3dad83 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/12] 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 b33aff15a02e47f55ce7377ac029a01a6451646e 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/12] Add `ishermitian` and `isantihermitian`

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

diff --git a/src/common/matrixproperties.jl b/src/common/matrixproperties.jl
index 61ff73e9..4a8698d3 100644
--- a/src/common/matrixproperties.jl
+++ b/src/common/matrixproperties.jl
@@ -144,3 +144,67 @@ function _ishermitian_exact_offdiag(Al, Au, ::Val{anti}) where {anti}
     end
     return true
 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 72ed30eafdcd3cc8b7cdfeb0c955840daa066127 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 03/12] blocked (anti)hermitian

---
 src/common/initialization.jl     | 50 -----------------
 src/implementations/hermitian.jl | 95 ++++++++++++++++++++++++++++++++
 src/interface/hermitian.jl       | 45 +++++++++++++++
 3 files changed, 140 insertions(+), 50 deletions(-)
 create mode 100644 src/implementations/hermitian.jl
 create mode 100644 src/interface/hermitian.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

From 4750e944ec47e449a256912fc9c08097ef4e1157 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 04/12] rename to project_

---
 src/implementations/hermitian.jl | 95 --------------------------------
 src/interface/hermitian.jl       | 45 ---------------
 2 files changed, 140 deletions(-)
 delete mode 100644 src/implementations/hermitian.jl
 delete mode 100644 src/interface/hermitian.jl

diff --git a/src/implementations/hermitian.jl b/src/implementations/hermitian.jl
deleted file mode 100644
index 048c0768..00000000
--- a/src/implementations/hermitian.jl
+++ /dev/null
@@ -1,95 +0,0 @@
-# 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
deleted file mode 100644
index 704fb974..00000000
--- a/src/interface/hermitian.jl
+++ /dev/null
@@ -1,45 +0,0 @@
-@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

From 7de2c634228362e4adcab72a9a0f28d824753468 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 05/12] also change includes

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

diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl
index 69ee4171..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

From 9a9f7670d4a59686c2ed7da50a6fd3ee6604cd4a 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 06/12] add blocked ishermitian

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

diff --git a/src/common/matrixproperties.jl b/src/common/matrixproperties.jl
index 4a8698d3..36856c07 100644
--- a/src/common/matrixproperties.jl
+++ b/src/common/matrixproperties.jl
@@ -151,30 +151,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
 
 """
@@ -183,28 +171,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...)
-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
+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
+
+# 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
+    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
-    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_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 074e56f6a7be96ef770ae1245a7302d634f5f7ab Mon Sep 17 00:00:00 2001
From: Jutho Haegeman <jutho.haegeman@ugent.be>
Date: Mon, 6 Oct 2025 14:54:24 +0200
Subject: [PATCH 07/12] add PolarNewton

---
 src/MatrixAlgebraKit.jl      |   1 +
 src/implementations/polar.jl | 149 ++++++++++++++++++++++++++++++++---
 src/interface/polar.jl       |   9 +++
 test/polar.jl                |  15 ++--
 4 files changed, 154 insertions(+), 20 deletions(-)

diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl
index 276523c9..99e3bf1a 100644
--- a/src/MatrixAlgebraKit.jl
+++ b/src/MatrixAlgebraKit.jl
@@ -34,6 +34,7 @@ export LAPACK_HouseholderQR, LAPACK_HouseholderLQ, LAPACK_Simple, LAPACK_Expert,
     LAPACK_QRIteration, LAPACK_Bisection, LAPACK_MultipleRelativelyRobustRepresentations,
     LAPACK_DivideAndConquer, LAPACK_Jacobi
 export LQViaTransposedQR
+export PolarViaSVD, PolarNewton
 export DiagonalAlgorithm
 export NativeBlocked
 export CUSOLVER_Simple, CUSOLVER_HouseholderQR, CUSOLVER_QRIteration, CUSOLVER_SVDPolar,
diff --git a/src/implementations/polar.jl b/src/implementations/polar.jl
index 1128ccda..abc902e4 100644
--- a/src/implementations/polar.jl
+++ b/src/implementations/polar.jl
@@ -11,7 +11,7 @@ function check_input(::typeof(left_polar!), A::AbstractMatrix, WP, ::AbstractAlg
     @assert W isa AbstractMatrix && P isa AbstractMatrix
     @check_size(W, (m, n))
     @check_scalar(W, A)
-    @check_size(P, (n, n))
+    isempty(P) || @check_size(P, (n, n))
     @check_scalar(P, A)
     return nothing
 end
@@ -21,7 +21,7 @@ function check_input(::typeof(right_polar!), A::AbstractMatrix, PWᴴ, ::Abstrac
     n >= m ||
         throw(ArgumentError("input matrix needs at least as many columns as rows"))
     @assert P isa AbstractMatrix && Wᴴ isa AbstractMatrix
-    @check_size(P, (m, m))
+    isempty(P) || @check_size(P, (m, m))
     @check_scalar(P, A)
     @check_size(Wᴴ, (m, n))
     @check_scalar(Wᴴ, A)
@@ -43,16 +43,18 @@ function initialize_output(::typeof(right_polar!), A::AbstractMatrix, ::Abstract
     return (P, Wᴴ)
 end
 
-# Implementation
-# --------------
+# Implementation via SVD
+# -----------------------
 function left_polar!(A::AbstractMatrix, WP, alg::PolarViaSVD)
     check_input(left_polar!, A, WP, alg)
     U, S, Vᴴ = svd_compact!(A, alg.svdalg)
     W, P = WP
     W = mul!(W, U, Vᴴ)
-    S .= sqrt.(S)
-    SsqrtVᴴ = lmul!(S, Vᴴ)
-    P = mul!(P, SsqrtVᴴ', SsqrtVᴴ)
+    if !isempty(P)
+        S .= sqrt.(S)
+        SsqrtVᴴ = lmul!(S, Vᴴ)
+        P = mul!(P, SsqrtVᴴ', SsqrtVᴴ)
+    end
     return (W, P)
 end
 function right_polar!(A::AbstractMatrix, PWᴴ, alg::PolarViaSVD)
@@ -60,8 +62,135 @@ function right_polar!(A::AbstractMatrix, PWᴴ, alg::PolarViaSVD)
     U, S, Vᴴ = svd_compact!(A, alg.svdalg)
     P, Wᴴ = PWᴴ
     Wᴴ = mul!(Wᴴ, U, Vᴴ)
-    S .= sqrt.(S)
-    USsqrt = rmul!(U, S)
-    P = mul!(P, USsqrt, USsqrt')
+    if !isempty(P)
+        S .= sqrt.(S)
+        USsqrt = rmul!(U, S)
+        P = mul!(P, USsqrt, USsqrt')
+    end
     return (P, Wᴴ)
 end
+
+# Implementation via Newton
+# --------------------------
+function left_polar!(A::AbstractMatrix, WP, alg::PolarNewton)
+    check_input(left_polar!, A, WP, alg)
+    W, P = WP
+    if isempty(P)
+        W = _left_polarnewton!(A, W, P; alg.kwargs...)
+        return W, P
+    else
+        W = _left_polarnewton!(copy(A), W, P; alg.kwargs...)
+        # we still need `A` to compute `P`
+        P = project_hermitian!(mul!(P, W', A))
+        return W, P
+    end
+end
+
+function right_polar!(A::AbstractMatrix, PWᴴ, alg::PolarNewton)
+    check_input(right_polar!, A, PWᴴ, alg)
+    P, Wᴴ = PWᴴ
+    if isempty(P)
+        Wᴴ = _right_polarnewton!(A, Wᴴ, P; alg.kwargs...)
+        return P, Wᴴ
+    else
+        Wᴴ = _right_polarnewton!(copy(A), Wᴴ, P; alg.kwargs...)
+        # we still need `A` to compute `P`
+        P = project_hermitian!(mul!(P, A, Wᴴ'))
+        return P, Wᴴ
+    end
+end
+
+# these methods only compute W and destroy A in the process
+function _left_polarnewton!(A::AbstractMatrix, W, P = similar(A, (0, 0)); tol = defaulttol(A), maxiter = 10)
+    m, n = size(A) # we must have m >= n
+    Rᴴinv = isempty(P) ? similar(P, (n, n)) : P # use P as workspace when available
+    if m > n # initial QR
+        Q, R = qr_compact!(A)
+        Rc = view(A, 1:n, 1:n)
+        copy!(Rc, R)
+        Rᴴinv = ldiv!(UpperTriangular(Rc)', one!(Rᴴinv))
+    else # m == n
+        R = A
+        Rc = view(W, 1:n, 1:n)
+        copy!(Rc, R)
+        Rᴴinv = ldiv!(lu!(Rc)', one!(Rᴴinv))
+    end
+    γ = sqrt(norm(Rᴴinv) / norm(R)) # scaling factor
+    rmul!(R, γ)
+    rmul!(Rᴴinv, 1 / γ)
+    R, Rᴴinv = _avgdiff!(R, Rᴴinv)
+    copy!(Rc, R)
+    i = 1
+    conv = norm(Rᴴinv, Inf)
+    while i < maxiter && conv > tol
+        Rᴴinv = ldiv!(lu!(Rc)', one!(Rᴴinv))
+        γ = sqrt(norm(Rᴴinv) / norm(R)) # scaling factor
+        rmul!(R, γ)
+        rmul!(Rᴴinv, 1 / γ)
+        R, Rᴴinv = _avgdiff!(R, Rᴴinv)
+        copy!(Rc, R)
+        conv = norm(Rᴴinv, Inf)
+        i += 1
+    end
+    if conv > tol
+        @warn "`left_polar!` via Newton iteration did not converge within $maxiter iterations (final residual: $conv)"
+    end
+    if m > n
+        return mul!(W, Q, Rc)
+    end
+    return W
+end
+
+function _right_polarnewton!(A::AbstractMatrix, Wᴴ, P = similar(A, (0, 0)); tol = defaulttol(A), maxiter = 10)
+    m, n = size(A) # we must have m <= n
+    Lᴴinv = isempty(P) ? similar(P, (m, m)) : P # use P as workspace when available
+    if m < n # initial QR
+        L, Q = lq_compact!(A)
+        Lc = view(A, 1:m, 1:m)
+        copy!(Lc, L)
+        Lᴴinv = ldiv!(LowerTriangular(Lc)', one!(Lᴴinv))
+    else # m == n
+        L = A
+        Lc = view(Wᴴ, 1:m, 1:m)
+        copy!(Lc, L)
+        Lᴴinv = ldiv!(lu!(Lc)', one!(Lᴴinv))
+    end
+    γ = sqrt(norm(Lᴴinv) / norm(L)) # scaling factor
+    rmul!(L, γ)
+    rmul!(Lᴴinv, 1 / γ)
+    L, Lᴴinv = _avgdiff!(L, Lᴴinv)
+    copy!(Lc, L)
+    i = 1
+    conv = norm(Lᴴinv, Inf)
+    while i < maxiter && conv > tol
+        Lᴴinv = ldiv!(lu!(Lc)', one!(Lᴴinv))
+        γ = sqrt(norm(Lᴴinv) / norm(L)) # scaling factor
+        rmul!(L, γ)
+        rmul!(Lᴴinv, 1 / γ)
+        L, Lᴴinv = _avgdiff!(L, Lᴴinv)
+        copy!(Lc, L)
+        conv = norm(Lᴴinv, Inf)
+        i += 1
+    end
+    if conv > tol
+        @warn "`right_polar!` via Newton iteration did not converge within $maxiter iterations (final residual: $conv)"
+    end
+    if m < n
+        return mul!(Wᴴ, Lc, Q)
+    end
+    return Wᴴ
+end
+
+# in place computation of the average and difference of two arrays
+function _avgdiff!(A::AbstractArray, B::AbstractArray)
+    axes(A) == axes(B) || throw(DimensionMismatch())
+    @simd for I in eachindex(A, B)
+        @inbounds begin
+            a = A[I]
+            b = B[I]
+            A[I] = (a + b) / 2
+            B[I] = b - a
+        end
+    end
+    return A, B
+end
diff --git a/src/interface/polar.jl b/src/interface/polar.jl
index dee14d23..0f407b4a 100644
--- a/src/interface/polar.jl
+++ b/src/interface/polar.jl
@@ -49,6 +49,15 @@ struct PolarViaSVD{SVDAlg} <: AbstractAlgorithm
     svdalg::SVDAlg
 end
 
+"""
+    PolarNewton(; maxiter = 10, tol = defaulttol(A))
+    
+Algorithm for computing the polar decomposition of a matrix `A` via
+scaled Newton iteration, with a maximum of `maxiter` iterations and
+until convergence up to tolerance `tol`.
+"""
+@algdef PolarNewton
+
 # Algorithm selection
 # -------------------
 default_polar_algorithm(A; kwargs...) = default_polar_algorithm(typeof(A); kwargs...)
diff --git a/test/polar.jl b/test/polar.jl
index f60153b0..ad5228c0 100644
--- a/test/polar.jl
+++ b/test/polar.jl
@@ -3,7 +3,6 @@ using Test
 using TestExtras
 using StableRNGs
 using LinearAlgebra: LinearAlgebra, I, isposdef
-using MatrixAlgebraKit: PolarViaSVD
 
 @testset "left_polar! for T = $T" for T in (Float32, Float64, ComplexF32, ComplexF64)
     rng = StableRNG(123)
@@ -11,14 +10,11 @@ using MatrixAlgebraKit: PolarViaSVD
     @testset "size ($m, $n)" for n in (37, m)
         k = min(m, n)
         if LinearAlgebra.LAPACK.version() < v"3.12.0"
-            algs = PolarViaSVD.(
-                (LAPACK_DivideAndConquer(), LAPACK_QRIteration(), LAPACK_Bisection())
-            )
+            svdalgs = (LAPACK_DivideAndConquer(), LAPACK_QRIteration(), LAPACK_Bisection())
         else
-            algs = PolarViaSVD.(
-                (LAPACK_DivideAndConquer(), LAPACK_QRIteration(), LAPACK_Bisection(), LAPACK_Jacobi())
-            )
+            svdalgs = (LAPACK_DivideAndConquer(), LAPACK_QRIteration(), LAPACK_Bisection(), LAPACK_Jacobi())
         end
+        algs = (PolarViaSVD.(svdalgs)..., PolarNewton())
         @testset "algorithm $alg" for alg in algs
             A = randn(rng, T, m, n)
 
@@ -45,9 +41,8 @@ end
     n = 54
     @testset "size ($m, $n)" for m in (37, n)
         k = min(m, n)
-        algs = PolarViaSVD.(
-            (LAPACK_DivideAndConquer(), LAPACK_QRIteration(), LAPACK_Bisection())
-        )
+        svdalgs = (LAPACK_DivideAndConquer(), LAPACK_QRIteration(), LAPACK_Bisection())
+        algs = (PolarViaSVD.(svdalgs)..., PolarNewton())
         @testset "algorithm $alg" for alg in algs
             A = randn(rng, T, m, n)
 

From 85a6a41570f3babaa737cb1026ef4370b023eba8 Mon Sep 17 00:00:00 2001
From: Jutho Haegeman <jutho.haegeman@ugent.be>
Date: Tue, 7 Oct 2025 12:55:35 +0200
Subject: [PATCH 08/12] add more tests and project_isometric

---
 src/MatrixAlgebraKit.jl            |  4 ++--
 src/implementations/projections.jl | 35 +++++++++++++++++++++++-----
 src/interface/projections.jl       | 37 +++++++++++++++++++++++++-----
 test/polar.jl                      | 18 +++++++++++++++
 test/projections.jl                | 37 +++++++++++++++++++++++++++++-
 5 files changed, 116 insertions(+), 15 deletions(-)

diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl
index 99e3bf1a..72b011ca 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 project_hermitian, project_antihermitian
-export project_hermitian!, project_antihermitian!
+export project_hermitian, project_antihermitian, project_isometric
+export project_hermitian!, project_antihermitian!, project_isometric!
 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/projections.jl b/src/implementations/projections.jl
index c905180b..b60cfd89 100644
--- a/src/implementations/projections.jl
+++ b/src/implementations/projections.jl
@@ -5,6 +5,8 @@ function copy_input(::typeof(project_hermitian), A::AbstractMatrix)
 end
 copy_input(::typeof(project_antihermitian), A) = copy_input(project_hermitian, A)
 
+copy_input(::typeof(project_isometric), A) = copy_input(left_polar, A)
+
 function check_input(::typeof(project_hermitian!), A::AbstractMatrix, B::AbstractMatrix, ::AbstractAlgorithm)
     LinearAlgebra.checksquare(A)
     n = Base.require_one_based_indexing(A)
@@ -18,6 +20,16 @@ function check_input(::typeof(project_antihermitian!), A::AbstractMatrix, B::Abs
     return nothing
 end
 
+function check_input(::typeof(project_isometric!), A::AbstractMatrix, W::AbstractMatrix, ::AbstractAlgorithm)
+    m, n = size(A)
+    m >= n ||
+        throw(ArgumentError("input matrix needs at least as many rows as columns"))
+    @assert W isa AbstractMatrix
+    @check_size(W, (m, n))
+    @check_scalar(W, A)
+    return nothing
+end
+
 # Outputs
 # -------
 function initialize_output(::typeof(project_hermitian!), A::AbstractMatrix, ::NativeBlocked)
@@ -27,15 +39,26 @@ function initialize_output(::typeof(project_antihermitian!), A::AbstractMatrix,
     return A
 end
 
+function initialize_output(::typeof(project_isometric!), A::AbstractMatrix, ::AbstractAlgorithm)
+    return similar(A)
+end
+
 # Implementation
 # --------------
-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...)
+function project_hermitian!(A::AbstractMatrix, Aₕ, alg::NativeBlocked)
+    check_input(project_hermitian!, A, Aₕ, alg)
+    return project_hermitian_native!(A, Aₕ, Val(false); alg.kwargs...)
 end
-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...)
+function project_antihermitian!(A::AbstractMatrix, Aₐ, alg::NativeBlocked)
+    check_input(project_antihermitian!, A, Aₐ, alg)
+    return project_hermitian_native!(A, Aₐ, Val(true); alg.kwargs...)
+end
+
+function project_isometric!(A::AbstractMatrix, W, alg::AbstractAlgorithm)
+    check_input(project_isometric!, A, W, alg)
+    noP = similar(W, (0, 0))
+    W, _ = left_polar!(A, (W, noP), alg)
+    return W
 end
 
 function project_hermitian_native!(A::AbstractMatrix, B::AbstractMatrix, anti::Val; blocksize = 32)
diff --git a/src/interface/projections.jl b/src/interface/projections.jl
index dbf62b92..9011426f 100644
--- a/src/interface/projections.jl
+++ b/src/interface/projections.jl
@@ -1,11 +1,13 @@
 @doc """
     project_hermitian(A; kwargs...)
     project_hermitian(A, alg)
-    project_hermitian!(A; kwargs...)
-    project_hermitian!(A, alg)
+    project_hermitian!(A, [Aₕ]; kwargs...)
+    project_hermitian!(A, [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`.
+For real matrices this corresponds to the symmetric part of `A`. In the bang method,
+the output storage can be provided via the optional argument `Aₕ`; by default it is
+equal to `A` and so the input matrix `A` is replaced by its hermitian projection.
 
 See also [`project_antihermitian`](@ref).
 """
@@ -14,16 +16,36 @@ See also [`project_antihermitian`](@ref).
 @doc """
     project_antihermitian(A; kwargs...)
     project_antihermitian(A, alg)
-    project_antihermitian!(A; kwargs...)
-    project_antihermitian!(A, alg)
+    project_antihermitian!(A, [Aₐ]; kwargs...)
+    project_antihermitian!(A, [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`.
+For real matrices this corresponds to the antisymmetric part of `A`. In the bang method,
+the output storage can be provided via the optional argument `Aₐ``; by default it is
+equal to `A` and so the input matrix `A` is replaced by its antihermitian projection.
 
 See also [`project_hermitian`](@ref).
 """
 @functiondef project_antihermitian
 
+@doc """
+    project_isometric(A; kwargs...)
+    project_isometric(A, alg)
+    project_isometric!(A, [W]; kwargs...)
+    project_isometric!(A, [W], alg)
+
+Compute the projection of `A` onto the manifold of isometric matrices, i.e. matrices
+satisfying `A' * A ≈ I`. This projection is computed via the polar decomposition, i.e.
+`W` corresponds to the first return value of `left_polar!`, but avoids computing the
+positive definite factor explicitly.
+
+!!! note
+    The bang method `project_isometric!` optionally accepts the output structure and
+    possibly destroys the input matrix `A`. Always use the return value of the function
+    as it may not always be possible to use the provided `W` as output.
+"""
+@functiondef project_isometric
+
 """
 NativeBlocked(; blocksize = 32)
 
@@ -43,3 +65,6 @@ for f in (:project_hermitian!, :project_antihermitian!)
         return default_hermitian_algorithm(A; kwargs...)
     end
 end
+
+default_algorithm(::typeof(project_isometric!), ::Type{A}; kwargs...) where {A} =
+    default_polar_algorithm(A; kwargs...)
diff --git a/test/polar.jl b/test/polar.jl
index ad5228c0..9e4ee341 100644
--- a/test/polar.jl
+++ b/test/polar.jl
@@ -32,6 +32,15 @@ using LinearAlgebra: LinearAlgebra, I, isposdef
             @test W * P ≈ A
             @test isisometry(W)
             @test isposdef(P)
+
+            noP = similar(P, (0, 0))
+            W2, P2 = @constinferred left_polar!(copy!(Ac, A), (W, noP), alg)
+            @test P2 === noP
+            @test W2 === W
+            @test isisometry(W)
+            P = W' * A # compute P explicitly to verify W correctness
+            @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P))
+            @test isposdef(project_hermitian!(P))
         end
     end
 end
@@ -60,6 +69,15 @@ end
             @test P * Wᴴ ≈ A
             @test isisometry(Wᴴ; side = :right)
             @test isposdef(P)
+
+            noP = similar(P, (0, 0))
+            P2, Wᴴ2 = @constinferred right_polar!(copy!(Ac, A), (noP, Wᴴ), alg)
+            @test P2 === noP
+            @test Wᴴ2 === Wᴴ
+            @test isisometry(Wᴴ; side = :right)
+            P = A * Wᴴ' # compute P explicitly to verify W correctness
+            @test ishermitian(P; rtol = MatrixAlgebraKit.defaulttol(P))
+            @test isposdef(project_hermitian!(P))
         end
     end
 end
diff --git a/test/projections.jl b/test/projections.jl
index 2eb2e28e..5698a702 100644
--- a/test/projections.jl
+++ b/test/projections.jl
@@ -2,7 +2,7 @@ using MatrixAlgebraKit
 using Test
 using TestExtras
 using StableRNGs
-using LinearAlgebra: Diagonal
+using LinearAlgebra: LinearAlgebra, Diagonal, norm
 
 const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64)
 
@@ -44,3 +44,38 @@ const BLASFloats = (Float32, Float64, ComplexF32, ComplexF64)
         @test Ba ≈ Aa
     end
 end
+
+@testset "project_isometric! for T = $T" for T in BLASFloats
+    rng = StableRNG(123)
+    m = 54
+    @testset "size ($m, $n)" for n in (37, m)
+        k = min(m, n)
+        if LinearAlgebra.LAPACK.version() < v"3.12.0"
+            svdalgs = (LAPACK_DivideAndConquer(), LAPACK_QRIteration(), LAPACK_Bisection())
+        else
+            svdalgs = (LAPACK_DivideAndConquer(), LAPACK_QRIteration(), LAPACK_Bisection(), LAPACK_Jacobi())
+        end
+        algs = (PolarViaSVD.(svdalgs)..., PolarNewton())
+        @testset "algorithm $alg" for alg in algs
+            A = randn(rng, T, m, n)
+            W = project_isometric(A, alg)
+            @test isisometry(W)
+            W2 = project_isometric(W, alg)
+            @test W2 ≈ W # stability of the projection
+            @test W * (W' * A) ≈ A
+
+            Ac = similar(A)
+            W2 = @constinferred project_isometric!(copy!(Ac, A), W, alg)
+            @test W2 === W
+            @test isisometry(W)
+
+            # test that W is closer to A then any other isometry
+            for k in 1:10
+                δA = randn(rng, T, m, n)
+                W = project_isometric(A, alg)
+                W2 = project_isometric(A + δA / 100, alg)
+                @test norm(A - W2) > norm(A - W)
+            end
+        end
+    end
+end

From 257a237cc44be3bcd49df362c1a24d2eafa6cbbe Mon Sep 17 00:00:00 2001
From: Jutho Haegeman <jutho.haegeman@ugent.be>
Date: Tue, 7 Oct 2025 14:39:47 +0200
Subject: [PATCH 09/12] fix/restore matrixproperties

---
 src/common/matrixproperties.jl | 76 ----------------------------------
 1 file changed, 76 deletions(-)

diff --git a/src/common/matrixproperties.jl b/src/common/matrixproperties.jl
index 36856c07..61ff73e9 100644
--- a/src/common/matrixproperties.jl
+++ b/src/common/matrixproperties.jl
@@ -144,79 +144,3 @@ function _ishermitian_exact_offdiag(Al, Au, ::Val{anti}) where {anti}
     end
     return true
 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, norm = LinearAlgebra.norm, kwargs...)
-    if iszero(atol) && iszero(rtol)
-        return ishermitian_exact(A; kwargs...)
-    else
-        return 2 * norm(project_antihermitian(A; kwargs...)) ≤ max(atol, rtol * norm(A))
-    end
-end
-function ishermitian_exact(A)
-    return A == A'
-end
-function ishermitian_exact(A::AbstractMatrix; kwargs...)
-    return _ishermitian_exact(A, Val(false); kwargs...)
-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, 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
-
-# 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
-    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
-    end
-    return true
-end

From 0bea6ed5e5532a1dbf4988afc2d38781a054827e Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Tue, 7 Oct 2025 12:59:41 -0400
Subject: [PATCH 10/12] move polar alg definition

---
 src/interface/decompositions.jl | 22 ++++++++++++++++++++++
 src/interface/polar.jl          | 19 -------------------
 2 files changed, 22 insertions(+), 19 deletions(-)

diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl
index e0ce2d52..543bbf35 100644
--- a/src/interface/decompositions.jl
+++ b/src/interface/decompositions.jl
@@ -117,6 +117,28 @@ const LAPACK_SVDAlgorithm = Union{
     LAPACK_Jacobi,
 }
 
+# =========================
+# Polar decompositions
+# =========================
+"""
+    PolarViaSVD(svdalg)
+
+Algorithm for computing the polar decomposition of a matrix `A` via the singular value
+decomposition (SVD) of `A`. The `svdalg` argument specifies the SVD algorithm to use.
+"""
+struct PolarViaSVD{SVDAlg} <: AbstractAlgorithm
+    svdalg::SVDAlg
+end
+
+"""
+    PolarNewton(; maxiter = 10, tol = defaulttol(A))
+
+Algorithm for computing the polar decomposition of a matrix `A` via
+scaled Newton iteration, with a maximum of `maxiter` iterations and
+until convergence up to tolerance `tol`.
+"""
+@algdef PolarNewton
+
 # =========================
 # DIAGONAL ALGORITHMS
 # =========================
diff --git a/src/interface/polar.jl b/src/interface/polar.jl
index 0f407b4a..38461093 100644
--- a/src/interface/polar.jl
+++ b/src/interface/polar.jl
@@ -39,25 +39,6 @@ See also [`left_polar(!)`](@ref left_polar).
 """
 @functiondef right_polar
 
-"""
-    PolarViaSVD(svdalg)
-    
-Algorithm for computing the polar decomposition of a matrix `A` via the singular value
-decomposition (SVD) of `A`. The `svdalg` argument specifies the SVD algorithm to use.
-"""
-struct PolarViaSVD{SVDAlg} <: AbstractAlgorithm
-    svdalg::SVDAlg
-end
-
-"""
-    PolarNewton(; maxiter = 10, tol = defaulttol(A))
-    
-Algorithm for computing the polar decomposition of a matrix `A` via
-scaled Newton iteration, with a maximum of `maxiter` iterations and
-until convergence up to tolerance `tol`.
-"""
-@algdef PolarNewton
-
 # Algorithm selection
 # -------------------
 default_polar_algorithm(A; kwargs...) = default_polar_algorithm(typeof(A); kwargs...)

From 923c1a96eedad2aa94bab36b93ad41d0882a05cd Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Tue, 7 Oct 2025 13:01:01 -0400
Subject: [PATCH 11/12] rename `svdalg` to `svd_alg`

---
 src/implementations/polar.jl    | 4 ++--
 src/interface/decompositions.jl | 6 +++---
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/implementations/polar.jl b/src/implementations/polar.jl
index abc902e4..9a767441 100644
--- a/src/implementations/polar.jl
+++ b/src/implementations/polar.jl
@@ -47,7 +47,7 @@ end
 # -----------------------
 function left_polar!(A::AbstractMatrix, WP, alg::PolarViaSVD)
     check_input(left_polar!, A, WP, alg)
-    U, S, Vᴴ = svd_compact!(A, alg.svdalg)
+    U, S, Vᴴ = svd_compact!(A, alg.svd_alg)
     W, P = WP
     W = mul!(W, U, Vᴴ)
     if !isempty(P)
@@ -59,7 +59,7 @@ function left_polar!(A::AbstractMatrix, WP, alg::PolarViaSVD)
 end
 function right_polar!(A::AbstractMatrix, PWᴴ, alg::PolarViaSVD)
     check_input(right_polar!, A, PWᴴ, alg)
-    U, S, Vᴴ = svd_compact!(A, alg.svdalg)
+    U, S, Vᴴ = svd_compact!(A, alg.svd_alg)
     P, Wᴴ = PWᴴ
     Wᴴ = mul!(Wᴴ, U, Vᴴ)
     if !isempty(P)
diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl
index 543bbf35..bdda6612 100644
--- a/src/interface/decompositions.jl
+++ b/src/interface/decompositions.jl
@@ -121,13 +121,13 @@ const LAPACK_SVDAlgorithm = Union{
 # Polar decompositions
 # =========================
 """
-    PolarViaSVD(svdalg)
+    PolarViaSVD(svd_alg)
 
 Algorithm for computing the polar decomposition of a matrix `A` via the singular value
-decomposition (SVD) of `A`. The `svdalg` argument specifies the SVD algorithm to use.
+decomposition (SVD) of `A`. The `svd_alg` argument specifies the SVD algorithm to use.
 """
 struct PolarViaSVD{SVDAlg} <: AbstractAlgorithm
-    svdalg::SVDAlg
+    svd_alg::SVDAlg
 end
 
 """

From d2252c0ce6e30d4ad7649149762b97d2ed6e7e23 Mon Sep 17 00:00:00 2001
From: Lukas Devos <ldevos98@gmail.com>
Date: Tue, 7 Oct 2025 13:01:33 -0400
Subject: [PATCH 12/12] remove `opnorm` from imports

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

diff --git a/src/MatrixAlgebraKit.jl b/src/MatrixAlgebraKit.jl
index 72b011ca..bb6434ab 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, lu!
-using LinearAlgebra: isposdef, issymmetric, opnorm
+using LinearAlgebra: isposdef, issymmetric
 using LinearAlgebra: Diagonal, diag, diagind, isdiag
 using LinearAlgebra: UpperTriangular, LowerTriangular
 using LinearAlgebra: BlasFloat, BlasReal, BlasComplex, BlasInt