Add truncerr truncation strategy (#59)

Author: mtfishman

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 @@ export eigen,
   svd,
   svd!,
   svdvals,
-  svdvals!
+  svdvals!,
+  truncerr
 
 using LinearAlgebra: LinearAlgebra
 using MatrixAlgebraKit
@@ -133,4 +134,43 @@ for (factorize, orth_f) in ((:factorize, :(MatrixAlgebra.orth)), (:factorize!, :
   end
 end
 
+using MatrixAlgebraKit: MatrixAlgebraKit, TruncationStrategy
+
+struct TruncationError{T<:Real} <: TruncationStrategy
+  atol::T
+  rtol::T
+  p::Int
+end
+
+"""
+    truncerr(; atol::Real=0, rtol::Real=0, p::Int=2)
+
+Create a truncation strategy for truncating such that the error in the factorization
+is smaller than `max(atol, rtol * norm)`, where the error is determined using the `p`-norm.
+"""
+function truncerr(; atol::Real=0, rtol::Real=0, p::Int=2)
+  return TruncationError(promote(atol, rtol)..., p)
+end
+
+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)
+  ϵᵖ = max(strategy.atol ^ strategy.p, strategy.rtol ^ strategy.p * normᵖ)
+  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