Add truncerr truncation strategy

Project.toml

@@ -1,7 +1,7 @@
 name = "TensorAlgebra"
 uuid = "68bd88dc-f39d-4e12-b2ca-f046b68fcc6a"
 authors = ["ITensor developers <[email protected]> and contributors"]
-version = "0.3.3"
+version = "0.3.4"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

docs/src/reference.md

@@ -1,5 +1,5 @@
 # Reference
 
 ```@autodocs
-Modules = [TensorAlgebra]
+Modules = [TensorAlgebra, TensorAlgebra.MatrixAlgebra]
 ```

src/MatrixAlgebra.jl

@@ -17,7 +17,8 @@
   svd,
   svd!,
   svdvals,
-  svdvals!
+  svdvals!,
+  truncerr
 
 using LinearAlgebra: LinearAlgebra
 using MatrixAlgebraKit
@@ -133,4 +134,44 @@
   end
 end
 
+using MatrixAlgebraKit: MatrixAlgebraKit, TruncationStrategy
+
+struct TruncationError{T<:Real} <: TruncationStrategy
+  ϵ::T
+  p::Int
+  relative::Bool
+end
+
+"""
+    truncerr(epsilon::Real, p::Int=2; relative=true)
+
+Create a truncation strategy for truncating such that the error in the factorization
+is smaller than `epsilon`, where the error is determined using the `p`-norm.
+
+The keyword argument `relative` specifies if the error `epsilon` is a relative error
+or not.
+"""
+truncerr(epsilon::Real, p::Int=2; relative=true) = TruncationError(epsilon, p, relative)
+
+function MatrixAlgebraKit.findtruncated(values::AbstractVector, strategy::TruncationError)
+  Base.require_one_based_indexing(values)
+  issorted(values; rev=true) || error("Not sorted.")
+  # norm(values, p) ^ p
+  normᵖ = sum(Base.Fix2(^, strategy.p) ∘ abs, values)
+  ϵᵖ = strategy.relative ? strategy.ϵ ^ strategy.p * normᵖ : strategy.ϵ ^ strategy.p
+  if ϵᵖ ≥ normᵖ
+    return Base.OneTo(0)
+  end
+  truncerrᵖ = zero(real(eltype(values)))
+  rank = length(values)
+  for i in reverse(eachindex(values))
+    truncerrᵖ += abs(values[i]) ^ strategy.p
+    if truncerrᵖ ≥ ϵᵖ
+      rank = i
+      break
+    end
+  end
+  return Base.OneTo(rank)
+end
+
 end

test/test_exports.jl

@@ -45,6 +45,7 @@ using TensorAlgebra: TensorAlgebra
     :svd!,
     :svdvals,
     :svdvals!,
+    :truncerr,
   ]
   @test issetequal(names(TensorAlgebra.MatrixAlgebra), exports)
 end