diff --git a/Project.toml b/Project.toml
index 39688f5..50d3d73 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "DiagonalArrays"
 uuid = "74fd4be6-21e2-4f6f-823a-4360d37c7a77"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
-version = "0.3.15"
+version = "0.3.16"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -11,11 +11,18 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MapBroadcast = "ebd9b9da-f48d-417c-9660-449667d60261"
 SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208"
 
+[weakdeps]
+MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
+
+[extensions]
+DiagonalArraysMatrixAlgebraKitExt = "MatrixAlgebraKit"
+
 [compat]
 ArrayLayouts = "1.10.4"
 DerivableInterfaces = "0.5.5"
 FillArrays = "1.13.0"
 LinearAlgebra = "1.10.0"
 MapBroadcast = "0.1.10"
+MatrixAlgebraKit = "0.2"
 SparseArraysBase = "0.7.2"
 julia = "1.10"
diff --git a/ext/DiagonalArraysMatrixAlgebraKitExt/DiagonalArraysMatrixAlgebraKitExt.jl b/ext/DiagonalArraysMatrixAlgebraKitExt/DiagonalArraysMatrixAlgebraKitExt.jl
new file mode 100644
index 0000000..fd90124
--- /dev/null
+++ b/ext/DiagonalArraysMatrixAlgebraKitExt/DiagonalArraysMatrixAlgebraKitExt.jl
@@ -0,0 +1,303 @@
+module DiagonalArraysMatrixAlgebraKitExt
+
+using DiagonalArrays:
+  AbstractDiagonalMatrix,
+  DeltaMatrix,
+  DiagonalMatrix,
+  ScaledDeltaMatrix,
+  δ,
+  diagview,
+  dual,
+  issquare
+using LinearAlgebra: LinearAlgebra, isdiag, ishermitian
+using MatrixAlgebraKit:
+  MatrixAlgebraKit,
+  AbstractAlgorithm,
+  check_input,
+  default_qr_algorithm,
+  eig_full,
+  eig_full!,
+  eig_vals,
+  eig_vals!,
+  eigh_full,
+  eigh_full!,
+  eigh_vals,
+  eigh_vals!,
+  left_null,
+  left_null!,
+  left_orth,
+  left_orth!,
+  left_polar,
+  left_polar!,
+  lq_compact,
+  lq_compact!,
+  lq_full,
+  lq_full!,
+  qr_compact,
+  qr_compact!,
+  qr_full,
+  qr_full!,
+  right_null,
+  right_null!,
+  right_orth,
+  right_orth!,
+  right_polar,
+  right_polar!,
+  svd_compact,
+  svd_compact!,
+  svd_full,
+  svd_full!,
+  svd_vals,
+  svd_vals!
+
+abstract type AbstractDiagonalAlgorithm <: AbstractAlgorithm end
+
+struct DeltaAlgorithm{KWargs<:NamedTuple} <: AbstractDiagonalAlgorithm
+  kwargs::KWargs
+end
+DeltaAlgorithm(; kwargs...) = DeltaAlgorithm((; kwargs...))
+
+struct ScaledDeltaAlgorithm{KWargs<:NamedTuple} <: AbstractDiagonalAlgorithm
+  kwargs::KWargs
+end
+ScaledDeltaAlgorithm(; kwargs...) = ScaledDeltaAlgorithm((; kwargs...))
+
+for f in [
+  :eig_full,
+  :eig_vals,
+  :eigh_full,
+  :eigh_vals,
+  :qr_compact,
+  :qr_full,
+  :left_null,
+  :left_orth,
+  :left_polar,
+  :lq_compact,
+  :lq_full,
+  :right_null,
+  :right_orth,
+  :right_polar,
+  :svd_compact,
+  :svd_full,
+  :svd_vals,
+]
+  @eval begin
+    MatrixAlgebraKit.copy_input(::typeof($f), a::AbstractDiagonalMatrix) = copy(a)
+  end
+end
+
+for f in [
+  :default_eig_algorithm,
+  :default_eigh_algorithm,
+  :default_lq_algorithm,
+  :default_qr_algorithm,
+  :default_polar_algorithm,
+  :default_svd_algorithm,
+]
+  @eval begin
+    function MatrixAlgebraKit.$f(::Type{<:DeltaMatrix}; kwargs...)
+      return DeltaAlgorithm(; kwargs...)
+    end
+    function MatrixAlgebraKit.$f(::Type{<:ScaledDeltaMatrix}; kwargs...)
+      return ScaledDeltaAlgorithm(; kwargs...)
+    end
+  end
+end
+
+for f in [
+  :eig_full!,
+  :eig_vals!,
+  :eigh_full!,
+  :eigh_vals!,
+  :left_null!,
+  :left_orth!,
+  :left_polar!,
+  :lq_compact!,
+  :lq_full!,
+  :qr_compact!,
+  :qr_full!,
+  :right_null!,
+  :right_orth!,
+  :right_polar!,
+  :svd_compact!,
+  :svd_full!,
+  :svd_vals!,
+]
+  for Alg in [:ScaledDeltaAlgorithm, :DeltaAlgorithm]
+    @eval begin
+      function MatrixAlgebraKit.initialize_output(::typeof($f), a, alg::$Alg)
+        return nothing
+      end
+    end
+  end
+end
+
+for f in [
+  :left_null!,
+  :left_orth!,
+  :left_polar!,
+  :lq_compact!,
+  :lq_full!,
+  :qr_compact!,
+  :qr_full!,
+  :right_null!,
+  :right_orth!,
+  :right_polar!,
+  :svd_compact!,
+  :svd_full!,
+  :svd_vals!,
+]
+  @eval begin
+    function MatrixAlgebraKit.check_input(::typeof($f), a, F, alg::DeltaAlgorithm)
+      @assert size(a, 1) == size(a, 2)
+      @assert isdiag(a)
+      @assert all(isone, diagview(a))
+      return nothing
+    end
+    function MatrixAlgebraKit.check_input(::typeof($f), a, F, alg::ScaledDeltaAlgorithm)
+      @assert size(a, 1) == size(a, 2)
+      @assert isdiag(a)
+      @assert allequal(diagview(a))
+      return nothing
+    end
+  end
+end
+for f in [:eig_full!, :eig_vals!, :eigh_full!, :eigh_vals!]
+  @eval begin
+    function MatrixAlgebraKit.check_input(::typeof($f), a, F, alg::DeltaAlgorithm)
+      @assert issquare(a)
+      @assert isdiag(a)
+      @assert all(isone, diagview(a))
+      return nothing
+    end
+    function MatrixAlgebraKit.check_input(::typeof($f), a, F, alg::ScaledDeltaAlgorithm)
+      @assert issquare(a)
+      @assert isdiag(a)
+      @assert allequal(diagview(a))
+      return nothing
+    end
+  end
+end
+
+# eig
+for Alg in [:DeltaAlgorithm, :ScaledDeltaAlgorithm]
+  @eval begin
+    function MatrixAlgebraKit.eig_full!(a, F, alg::$Alg)
+      check_input(eig_full!, a, F, alg)
+      d = complex(a)
+      v = δ(complex(eltype(a)), axes(a))
+      return (d, v)
+    end
+    function MatrixAlgebraKit.eigh_full!(a, F, alg::$Alg)
+      check_input(eigh_full!, a, F, alg)
+      ishermitian(a) || throw(ArgumentError("Matrix must be Hermitian"))
+      d = real(a)
+      v = δ(eltype(a), axes(a))
+      return (d, v)
+    end
+    function MatrixAlgebraKit.eig_vals!(a, F, alg::$Alg)
+      check_input(eig_vals!, a, F, alg)
+      return complex(diagview(a))
+    end
+    function MatrixAlgebraKit.eigh_vals!(a, F, alg::$Alg)
+      check_input(eigh_vals!, a, F, alg)
+      return real(diagview(a))
+    end
+  end
+end
+
+# svd
+for f in [:svd_compact!, :svd_full!]
+  @eval begin
+    function MatrixAlgebraKit.$f(a, F, alg::DeltaAlgorithm)
+      check_input($f, a, F, alg)
+      u = δ(eltype(a), (axes(a, 1), dual(axes(a, 1))))
+      s = real(a)
+      v = δ(eltype(a), (dual(axes(a, 2)), axes(a, 2)))
+      return (u, s, v)
+    end
+    function MatrixAlgebraKit.$f(a, F, alg::ScaledDeltaAlgorithm)
+      check_input($f, a, F, alg)
+      diagvalue = only(unique(diagview(a)))
+      u = δ(eltype(a), (axes(a, 1), dual(axes(a, 1))))
+      s = abs(diagvalue) * δ(Bool, axes(a))
+      # Sign is applied arbitarily to `v`, alternatively
+      # we could apply it to `u`.
+      v = sign(diagvalue) * δ(Bool, (dual(axes(a, 2)), axes(a, 2)))
+      return (u, s, v)
+    end
+  end
+end
+function MatrixAlgebraKit.svd_vals!(a, F, alg::DeltaAlgorithm)
+  check_input(svd_vals!, a, F, alg)
+  # Using `real` instead of `abs.` helps to preserve `Ones`.
+  return real(diagview(a))
+end
+function MatrixAlgebraKit.svd_vals!(a, F, alg::ScaledDeltaAlgorithm)
+  check_input(svd_vals!, a, F, alg)
+  return abs.(diagview(a))
+end
+
+# orth
+# left_orth is implicitly defined by defining backends like
+# qr_compact and left_polar.
+for f in [:left_polar!, :qr_compact!, :qr_full!]
+  @eval begin
+    function MatrixAlgebraKit.$f(a, F, alg::DeltaAlgorithm)
+      check_input($f, a, F, alg)
+      q = δ(eltype(a), (axes(a, 1), dual(axes(a, 1))))
+      r = copy(a)
+      return (q, r)
+    end
+    function MatrixAlgebraKit.$f(a, F, alg::ScaledDeltaAlgorithm)
+      check_input($f, a, F, alg)
+      diagvalue = only(unique(diagview(a)))
+      q = sign(diagvalue) * δ(Bool, (axes(a, 1), dual(axes(a, 1))))
+      # We're a bit pessimistic about the element type for type stability,
+      # since in the future we might provide the option to do non-positive QR.
+      r = eltype(a)(abs(diagvalue)) * δ(Bool, axes(a))
+      return (q, r)
+    end
+  end
+end
+# right_orth is implicitly defined by defining backends like
+# lq_compact and right_polar.
+for f in [:right_polar!, :lq_compact!, :lq_full!]
+  @eval begin
+    function MatrixAlgebraKit.$f(a, F, alg::DeltaAlgorithm)
+      check_input($f, a, F, alg)
+      l = copy(a)
+      q = δ(eltype(a), (dual(axes(a, 2)), axes(a, 2)))
+      return (l, q)
+    end
+    function MatrixAlgebraKit.$f(a, F, alg::ScaledDeltaAlgorithm)
+      check_input($f, a, F, alg)
+      diagvalue = only(unique(diagview(a)))
+      # We're a bit pessimistic about the element type for type stability,
+      # since in the future we might provide the option to do non-positive LQ.
+      l = eltype(a)(abs(diagvalue)) * δ(Bool, axes(a))
+      q = sign(diagvalue) * δ(Bool, (dual(axes(a, 2)), axes(a, 2)))
+      return (l, q)
+    end
+  end
+end
+
+# null
+for T in [:DeltaMatrix, :ScaledDeltaMatrix]
+  @eval begin
+    # TODO: Right now we can't overload `left_null!` on an algorithm,
+    # make a PR to MatrixAlgebraKit.jl to allow that.
+    function MatrixAlgebraKit.left_null!(a::$T, F)
+      check_input(left_null!, a, F, default_qr_algorithm(a))
+      return error("Not implemented.")
+    end
+    # TODO: Right now we can't overload `right_null!` on an algorithm,
+    # make a PR to MatrixAlgebraKit.jl to allow that.
+    function MatrixAlgebraKit.right_null!(a::$T, F)
+      check_input(right_null!, a, F, default_qr_algorithm(a))
+      return error("Not implemented.")
+    end
+  end
+end
+
+end
diff --git a/src/DiagonalArrays.jl b/src/DiagonalArrays.jl
index 51550ac..5140355 100644
--- a/src/DiagonalArrays.jl
+++ b/src/DiagonalArrays.jl
@@ -1,5 +1,6 @@
 module DiagonalArrays
 
+include("dual.jl")
 include("diaginterface/diaginterface.jl")
 include("diaginterface/diagindex.jl")
 include("diaginterface/diagindices.jl")
diff --git a/src/abstractdiagonalarray/abstractdiagonalarray.jl b/src/abstractdiagonalarray/abstractdiagonalarray.jl
index 3c0bdd6..1c23777 100644
--- a/src/abstractdiagonalarray/abstractdiagonalarray.jl
+++ b/src/abstractdiagonalarray/abstractdiagonalarray.jl
@@ -1,3 +1,21 @@
 using SparseArraysBase: AbstractSparseArray
 
 abstract type AbstractDiagonalArray{T,N} <: AbstractSparseArray{T,N} end
+const AbstractDiagonalMatrix{T} = AbstractDiagonalArray{T,2}
+const AbstractDiagonalVector{T} = AbstractDiagonalArray{T,1}
+
+using LinearAlgebra: LinearAlgebra, ishermitian, isposdef, issymmetric
+LinearAlgebra.ishermitian(a::AbstractDiagonalMatrix{<:Real}) = issquare(a)
+function LinearAlgebra.ishermitian(a::AbstractDiagonalMatrix{<:Number})
+  return issquare(a) && isreal(diagview(a))
+end
+function LinearAlgebra.ishermitian(a::AbstractDiagonalMatrix)
+  return issquare(a) && all(ishermitian, diagview(a))
+end
+LinearAlgebra.issymmetric(a::AbstractDiagonalMatrix{<:Number}) = issquare(a)
+function LinearAlgebra.issymmetric(a::AbstractDiagonalMatrix)
+  return issquare(a) && all(issymmetric, diagview(a))
+end
+function LinearAlgebra.isposdef(a::AbstractDiagonalMatrix)
+  return issquare(a) && all(isposdef, diagview(a))
+end
diff --git a/src/diagonalarray/diagonalarray.jl b/src/diagonalarray/diagonalarray.jl
index 96e5d67..47ce8b2 100644
--- a/src/diagonalarray/diagonalarray.jl
+++ b/src/diagonalarray/diagonalarray.jl
@@ -170,9 +170,13 @@ function Base.similar(a::DiagonalArray, unstored::Unstored)
   return DiagonalArray(undef, unstored)
 end
 
-# This definition is helpful for immutable diagonals
+# These definitions are helpful for immutable diagonals
 # such as FillArrays.
-Base.copy(a::DiagonalArray) = DiagonalArray(copy(diagview(a)), axes(a))
+for f in [:complex, :copy, :imag, :real]
+  @eval begin
+    Base.$f(a::DiagonalArray) = DiagonalArray($f(diagview(a)), axes(a))
+  end
+end
 
 # DiagonalArrays interface.
 diagview(a::DiagonalArray) = a.diag
diff --git a/src/diagonalarray/diagonalmatrix.jl b/src/diagonalarray/diagonalmatrix.jl
index 2fff225..4eb6fab 100644
--- a/src/diagonalarray/diagonalmatrix.jl
+++ b/src/diagonalarray/diagonalmatrix.jl
@@ -8,7 +8,7 @@ using LinearAlgebra: LinearAlgebra
 
 function mul_diagviews(a1, a2)
   # TODO: Compare that duals are equal, or define a function to overload.
-  axes(a1, 2) == axes(a2, 1) || throw(
+  dual(axes(a1, 2)) == axes(a2, 1) || throw(
     DimensionMismatch(
       lazy"Incompatible dimensions for multiplication: $(axes(a1)) and $(axes(a2))"
     ),
diff --git a/src/dual.jl b/src/dual.jl
new file mode 100644
index 0000000..9780a3b
--- /dev/null
+++ b/src/dual.jl
@@ -0,0 +1,3 @@
+# TODO: Define `TensorProducts.dual`.
+dual(x) = x
+issquare(a::AbstractMatrix) = (axes(a, 1) == dual(axes(a, 2)))
diff --git a/test/Project.toml b/test/Project.toml
index c8c35c0..97d324d 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,8 +4,10 @@ DerivableInterfaces = "6c5e35bf-e59e-4898-b73c-732dcc4ba65f"
 DiagonalArrays = "74fd4be6-21e2-4f6f-823a-4360d37c7a77"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SparseArraysBase = "0d5efcca-f356-4864-8770-e1ed8d78f208"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 TensorAlgebra = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
@@ -16,8 +18,10 @@ DerivableInterfaces = "0.5"
 DiagonalArrays = "0.3"
 FillArrays = "1"
 LinearAlgebra = "1"
+MatrixAlgebraKit = "0.2.5"
 SafeTestsets = "0.1"
 SparseArraysBase = "0.7"
+StableRNGs = "1"
 Suppressor = "0.2"
 TensorAlgebra = "0.3.10"
 Test = "1"
diff --git a/test/test_basics.jl b/test/test_basics.jl
index dbc5f16..4bc24cb 100644
--- a/test/test_basics.jl
+++ b/test/test_basics.jl
@@ -17,7 +17,7 @@ using DiagonalArrays:
   diagview
 using FillArrays: Fill, Ones
 using SparseArraysBase: SparseArrayDOK, sparsezeros, storedlength
-using LinearAlgebra: Diagonal, mul!
+using LinearAlgebra: Diagonal, mul!, ishermitian, isposdef, issymmetric
 
 @testset "Test DiagonalArrays" begin
   @testset "DiagonalArray (eltype=$elt)" for elt in (
@@ -135,6 +135,28 @@ using LinearAlgebra: Diagonal, mul!
       c = DiagonalArray{elt}(undef, (2, 3))
       @test_broken c .= a .+ 2
     end
+    @testset "LinearAlgebra matrix properties" begin
+      @test ishermitian(DiagonalMatrix([1, 2]))
+      @test !ishermitian(DiagonalMatrix([1, 2], (2, 3)))
+      @test !ishermitian(DiagonalMatrix([1 + 1im, 2 + 2im]))
+      @test ishermitian(DiagonalMatrix([ones(2, 2), ones(3, 3)]))
+      @test !ishermitian(DiagonalMatrix([randn(2, 2), randn(3, 3)]))
+
+      @test issymmetric(DiagonalMatrix([1, 2]))
+      @test !issymmetric(DiagonalMatrix([1, 2], (2, 3)))
+      @test issymmetric(DiagonalMatrix([1 + 1im, 2 + 2im]))
+      @test issymmetric(DiagonalMatrix([ones(2, 2), ones(3, 3)]))
+      @test !issymmetric(DiagonalMatrix([randn(2, 2), randn(3, 3)]))
+      @test !issymmetric(DiagonalMatrix([randn(2, 2), randn(2, 3)]))
+
+      @test isposdef(DiagonalMatrix([1, 2]))
+      @test !isposdef(DiagonalMatrix([1, -2]))
+      @test !isposdef(DiagonalMatrix([1, 2], (2, 3)))
+      @test !isposdef(DiagonalMatrix([1 + 1im, 2 + 2im]))
+      @test isposdef(DiagonalMatrix([[1 0; 0 1], [2 0; 0 2]]))
+      @test !isposdef(DiagonalMatrix([randn(2, 2), randn(3, 3)]))
+      @test !isposdef(DiagonalMatrix([randn(2, 2), randn(2, 3)]))
+    end
     @testset "Matrix multiplication" begin
       a1 = DiagonalArray{elt}(undef, (2, 3))
       a1[1, 1] = 11
diff --git a/test/test_matrixalgebrakit.jl b/test/test_matrixalgebrakit.jl
new file mode 100644
index 0000000..e281a83
--- /dev/null
+++ b/test/test_matrixalgebrakit.jl
@@ -0,0 +1,213 @@
+using Test: @test, @testset
+using LinearAlgebra: Diagonal
+using DiagonalArrays: DiagonalArrays, DeltaMatrix, ScaledDeltaMatrix, δ, dual
+using FillArrays: Ones
+using MatrixAlgebraKit:
+  eig_full,
+  eig_vals,
+  eigh_full,
+  eigh_vals,
+  left_orth,
+  left_polar,
+  lq_compact,
+  lq_full,
+  qr_compact,
+  qr_full,
+  right_orth,
+  right_polar,
+  svd_compact,
+  svd_full,
+  svd_vals
+using StableRNGs: StableRNG
+
+struct SU2 <: AbstractUnitRange{Int}
+  j::Int
+  isdual::Bool
+end
+SU2(j::Int) = SU2(j, false)
+Base.:(==)(s1::SU2, s2::SU2) = ((s1.j == s2.j) && (s1.isdual == s2.isdual))
+Base.first(s::SU2) = 1
+Base.last(s::SU2) = 2 * s.j + 1
+DiagonalArrays.dual(s::SU2) = SU2(s.j, !s.isdual)
+
+@testset "MatrixAlgebraKitExt" begin
+  @testset "DeltaMatrix factorizations (eltype=$elt)" for elt in (
+    Float32, Float64, ComplexF32, ComplexF64
+  )
+    @testset "SVD" begin
+      for f in (svd_compact, svd_full)
+        ax = SU2(2)
+        a = δ(elt, (ax, ax))
+        u, s, v = f(a)
+        @test u * s * v ≡ a
+        @test u ≡ δ(elt, (ax, dual(ax)))
+        @test s ≡ δ(real(elt), (ax, ax))
+        @test v ≡ δ(elt, (dual(ax), ax))
+      end
+    end
+    @testset "SVD values" begin
+      ax = SU2(2)
+      a = δ(elt, (ax, ax))
+      s = svd_vals(a)
+      @test s ≡ Ones(real(elt), length(ax))
+    end
+    @testset "left orth" begin
+      for f in (left_orth, left_polar, qr_compact, qr_full)
+        ax = SU2(2)
+        a = δ(elt, (ax, ax))
+        q, r = f(a)
+        @test q * r ≡ a
+        @test q ≡ δ(elt, (ax, dual(ax)))
+        @test r ≡ δ(elt, (ax, ax))
+      end
+    end
+    @testset "right orth" begin
+      for f in (lq_compact, lq_full, right_orth, right_polar)
+        ax = SU2(2)
+        a = δ(elt, (ax, ax))
+        l, q = f(a)
+        @test l * q ≡ a
+        @test l ≡ δ(elt, (ax, ax))
+        @test q ≡ δ(elt, (dual(ax), ax))
+      end
+    end
+    @testset "Eigendecomposition" begin
+      ax = SU2(2)
+      a = δ(elt, (dual(ax), ax))
+      d, v = eig_full(a)
+      @test a * v ≡ v * d
+      @test d ≡ δ(complex(elt), (dual(ax), ax))
+      @test v ≡ δ(complex(elt), (dual(ax), ax))
+    end
+    @testset "Hermitian eigendecomposition" begin
+      ax = SU2(2)
+      a = δ(elt, (dual(ax), ax))
+      d, v = eigh_full(a)
+      @test a * v ≡ v * d
+      @test d ≡ δ(real(elt), (dual(ax), ax))
+      @test v ≡ δ(elt, (dual(ax), ax))
+    end
+    @testset "Eigenvalues" begin
+      ax = SU2(2)
+      a = δ(elt, (dual(ax), ax))
+      d = eig_vals(a)
+      @test d ≡ Ones{complex(elt)}(length(ax))
+    end
+    @testset "Hermitian eigenvalues" begin
+      ax = SU2(2)
+      a = δ(elt, (dual(ax), ax))
+      d = eigh_vals(a)
+      @test d ≡ Ones{real(elt)}(length(ax))
+    end
+    @testset "left null" begin
+      ax = SU2(2)
+      a = δ(elt, (ax, ax))
+      @test_broken left_null(a)
+    end
+    @testset "right null" begin
+      ax = SU2(2)
+      a = δ(elt, (ax, ax))
+      @test_broken right_null(a)
+    end
+  end
+
+  @testset "ScaledDeltaMatrix factorizations (eltype=$elt)" for elt in (
+    Float32, Float64, ComplexF32, ComplexF64
+  )
+    @testset "SVD" begin
+      for f in (svd_compact, svd_full)
+        ax = SU2(2)
+        rng = StableRNG(1234)
+        scale = randn(rng, elt)
+        a = scale * δ(elt, (ax, ax))
+        u, s, v = f(a)
+        @test u * s * v ≡ a
+        @test u ≡ δ(elt, (ax, dual(ax)))
+        @test s ≡ abs(scale) * δ(real(elt), (ax, ax))
+        @test v ≡ sign(scale) * δ(elt, (dual(ax), ax))
+      end
+    end
+    @testset "SVD values" begin
+      ax = SU2(2)
+      rng = StableRNG(1234)
+      scale = randn(rng, elt)
+      a = scale * δ(elt, (ax, ax))
+      s = svd_vals(a)
+      @test s ≡ abs(scale) * Ones(real(elt), length(ax))
+    end
+    @testset "left orth" begin
+      for f in (left_orth, left_polar, qr_compact, qr_full)
+        ax = SU2(2)
+        rng = StableRNG(1234)
+        scale = randn(rng, elt)
+        a = scale * δ(elt, (ax, ax))
+        q, r = f(a)
+        @test q * r ≡ a
+        @test q ≡ sign(scale) * δ(elt, (ax, dual(ax)))
+        @test r ≡ abs(scale) * δ(elt, (ax, ax))
+      end
+    end
+    @testset "right orth" begin
+      for f in (lq_compact, lq_full, right_orth, right_polar)
+        ax = SU2(2)
+        rng = StableRNG(1234)
+        scale = randn(rng, elt)
+        a = scale * δ(elt, (ax, ax))
+        l, q = f(a)
+        @test l * q ≡ a
+        @test l ≡ abs(scale) * δ(elt, (ax, ax))
+        @test q ≡ sign(scale) * δ(elt, (dual(ax), ax))
+      end
+    end
+    @testset "Eigendecomposition" begin
+      ax = SU2(2)
+      rng = StableRNG(1234)
+      scale = randn(rng, elt)
+      a = scale * δ(elt, (dual(ax), ax))
+      d, v = eig_full(a)
+      @test a * v ≡ v * d
+      @test d ≡ scale * δ(complex(elt), (dual(ax), ax))
+      @test v ≡ δ(complex(elt), (dual(ax), ax))
+    end
+    @testset "Hermitian eigendecomposition" begin
+      ax = SU2(2)
+      rng = StableRNG(1234)
+      scale = randn(rng, real(elt))
+      a = scale * δ(elt, (dual(ax), ax))
+      d, v = eigh_full(a)
+      @test a * v ≡ v * d
+      @test d ≡ scale * δ(real(elt), (dual(ax), ax))
+      @test v ≡ δ(elt, (dual(ax), ax))
+    end
+    @testset "Eigenvalues" begin
+      ax = SU2(2)
+      rng = StableRNG(1234)
+      scale = randn(rng, elt)
+      a = scale * δ(elt, (dual(ax), ax))
+      d = eig_vals(a)
+      @test d ≡ scale * Ones{complex(elt)}(length(ax))
+    end
+    @testset "Hermitian eigenvalues" begin
+      ax = SU2(2)
+      rng = StableRNG(1234)
+      scale = randn(rng, real(elt))
+      a = scale * δ(elt, (dual(ax), ax))
+      d = eigh_vals(a)
+      @test d ≡ scale * Ones{real(elt)}(length(ax))
+    end
+    @testset "left null" begin
+      ax = SU2(2)
+      rng = StableRNG(1234)
+      scale = randn(rng, real(elt))
+      a = scale * δ(elt, (ax, ax))
+      @test_broken left_null(a)
+    end
+    @testset "right null" begin
+      ax = SU2(2)
+      rng = StableRNG(1234)
+      scale = randn(rng, real(elt))
+      a = scale * δ(elt, (ax, ax))
+      @test_broken right_null(a)
+    end
+  end
+end