From b03cba9111184794e433f34805227863586331c0 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 27 May 2025 18:27:09 -0400
Subject: [PATCH 1/3] Implement blockwise polar decomposition

---
 src/BlockSparseArrays.jl    |  1 +
 src/factorizations/polar.jl | 53 +++++++++++++++++++++++++++++++++++++
 test/test_factorizations.jl | 22 +++++++++++++++
 3 files changed, 76 insertions(+)
 create mode 100644 src/factorizations/polar.jl

diff --git a/src/BlockSparseArrays.jl b/src/BlockSparseArrays.jl
index 40b55492..bb3bf9e0 100644
--- a/src/BlockSparseArrays.jl
+++ b/src/BlockSparseArrays.jl
@@ -48,5 +48,6 @@ include("factorizations/svd.jl")
 include("factorizations/truncation.jl")
 include("factorizations/qr.jl")
 include("factorizations/lq.jl")
+include("factorizations/polar.jl")
 
 end
diff --git a/src/factorizations/polar.jl b/src/factorizations/polar.jl
new file mode 100644
index 00000000..9192f011
--- /dev/null
+++ b/src/factorizations/polar.jl
@@ -0,0 +1,53 @@
+using MatrixAlgebraKit:
+  MatrixAlgebraKit,
+  PolarViaSVD,
+  check_input,
+  default_algorithm,
+  left_polar!,
+  right_polar!,
+  svd_compact!
+
+function MatrixAlgebraKit.check_input(::typeof(left_polar!), A::AbstractBlockSparseMatrix)
+  @views for I in eachblockstoredindex(A)
+    m, n = size(A[I])
+    m >= n ||
+      throw(ArgumentError("each input matrix block needs at least as many rows as columns"))
+  end
+  return nothing
+end
+function MatrixAlgebraKit.check_input(::typeof(right_polar!), A::AbstractBlockSparseMatrix)
+  @views for I in eachblockstoredindex(A)
+    m, n = size(A[I])
+    m <= n ||
+      throw(ArgumentError("each input matrix block needs at least as many columns as rows"))
+  end
+  return nothing
+end
+
+function MatrixAlgebraKit.left_polar!(A::AbstractBlockSparseMatrix, alg::PolarViaSVD)
+  check_input(left_polar!, A)
+  # TODO: Use more in-place operations here, avoid `copy`.
+  U, S, Vᴴ = svd_compact!(A, alg.svdalg)
+  W = U * Vᴴ
+  P = copy(Vᴴ') * S * Vᴴ
+  return (W, P)
+end
+function MatrixAlgebraKit.right_polar!(A::AbstractBlockSparseMatrix, alg::PolarViaSVD)
+  check_input(right_polar!, A)
+  # TODO: Use more in-place operations here, avoid `copy`.
+  U, S, Vᴴ = svd_compact!(A, alg.svdalg)
+  Wᴴ = U * Vᴴ
+  P = U * S * copy(U')
+  return (P, Wᴴ)
+end
+
+function MatrixAlgebraKit.default_algorithm(
+  ::typeof(left_polar!), a::AbstractBlockSparseMatrix; kwargs...
+)
+  return PolarViaSVD(default_algorithm(svd_compact!, a; kwargs...))
+end
+function MatrixAlgebraKit.default_algorithm(
+  ::typeof(right_polar!), a::AbstractBlockSparseMatrix; kwargs...
+)
+  return PolarViaSVD(default_algorithm(svd_compact!, a; kwargs...))
+end
diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl
index e528ca69..70639203 100644
--- a/test/test_factorizations.jl
+++ b/test/test_factorizations.jl
@@ -1,10 +1,12 @@
 using BlockArrays: Block, BlockedMatrix, BlockedVector, blocks, mortar
 using BlockSparseArrays: BlockSparseArray, BlockDiagonal, eachblockstoredindex
 using MatrixAlgebraKit:
+  left_polar,
   lq_compact,
   lq_full,
   qr_compact,
   qr_full,
+  right_polar,
   svd_compact,
   svd_full,
   svd_trunc,
@@ -215,3 +217,23 @@ end
     @test A ≈ L * Q
   end
 end
+
+@testset "left_polar (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
+  A = BlockSparseArray{T}(undef, ([3, 4], [2, 3]))
+  A[Block(1, 1)] = randn(T, 3, 2)
+  A[Block(2, 2)] = randn(T, 4, 3)
+
+  U, C = left_polar(A)
+  @test U * C ≈ A
+  @test Matrix(U'U) ≈ LinearAlgebra.I
+end
+
+@testset "right_polar (T=$T)" for T in (Float32, Float64, ComplexF32, ComplexF64)
+  A = BlockSparseArray{T}(undef, ([2, 3], [3, 4]))
+  A[Block(1, 1)] = randn(T, 2, 3)
+  A[Block(2, 2)] = randn(T, 3, 4)
+
+  C, U = right_polar(A)
+  @test C * U ≈ A
+  @test Matrix(U * U') ≈ LinearAlgebra.I
+end

From 49fbe814c243330045075fa08b6a50c1b35d23c2 Mon Sep 17 00:00:00 2001
From: mtfishman <mfishman@flatironinstitute.org>
Date: Tue, 27 May 2025 18:29:50 -0400
Subject: [PATCH 2/3] Bump version

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 7f3b7dfc..8a88983a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "BlockSparseArrays"
 uuid = "2c9a651f-6452-4ace-a6ac-809f4280fbb4"
 authors = ["ITensor developers <support@itensor.org> and contributors"]
-version = "0.6.3"
+version = "0.6.4"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"

From f22612cbae7141999482a63b3ccc321242f724eb Mon Sep 17 00:00:00 2001
From: Matt Fishman <mtfishman@users.noreply.github.com>
Date: Tue, 27 May 2025 19:03:35 -0400
Subject: [PATCH 3/3] Comments

---
 src/factorizations/polar.jl | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/factorizations/polar.jl b/src/factorizations/polar.jl
index 9192f011..9b9c2831 100644
--- a/src/factorizations/polar.jl
+++ b/src/factorizations/polar.jl
@@ -29,6 +29,9 @@ function MatrixAlgebraKit.left_polar!(A::AbstractBlockSparseMatrix, alg::PolarVi
   # TODO: Use more in-place operations here, avoid `copy`.
   U, S, Vᴴ = svd_compact!(A, alg.svdalg)
   W = U * Vᴴ
+  # TODO: `copy` is required for now because of:
+  # https://github.com/ITensor/BlockSparseArrays.jl/issues/24
+  # Remove when that is fixed.
   P = copy(Vᴴ') * S * Vᴴ
   return (W, P)
 end
@@ -37,6 +40,9 @@ function MatrixAlgebraKit.right_polar!(A::AbstractBlockSparseMatrix, alg::PolarV
   # TODO: Use more in-place operations here, avoid `copy`.
   U, S, Vᴴ = svd_compact!(A, alg.svdalg)
   Wᴴ = U * Vᴴ
+  # TODO: `copy` is required for now because of:
+  # https://github.com/ITensor/BlockSparseArrays.jl/issues/24
+  # Remove when that is fixed.
   P = U * S * copy(U')
   return (P, Wᴴ)
 end