Implement min_blockdim keyword for blocksparse SVD (#923)

* Implement block_mindim keyword for QN SVD

* Add test and improve some docs and an export

* Change parameter name to min_blockdim

* Formatting

NDTensors/src/blocksparse/linearalgebra.jl

@@ -4,18 +4,21 @@ const DiagBlockSparseMatrix{ElT,StoreT,IndsT} = DiagBlockSparseTensor{ElT,2,Stor
 const DiagMatrix{ElT,StoreT,IndsT} = DiagTensor{ElT,2,StoreT,IndsT}
 
 function _truncated_blockdim(
-  S::DiagMatrix, docut::Float64; singular_values=false, truncate=true
+  S::DiagMatrix, docut::Float64; singular_values=false, truncate=true, min_blockdim=0
 )
-  !truncate && return diaglength(S)
+  full_dim = diaglength(S)
+  !truncate && return full_dim
+  min_blockdim = min(min_blockdim, full_dim)
   newdim = 0
   val = singular_values ? getdiagindex(S, newdim + 1)^2 : abs(getdiagindex(S, newdim + 1))
-  while newdim + 1 ≤ diaglength(S) && val > docut
+  while newdim + 1 ≤ full_dim && val > docut
     newdim += 1
-    if newdim + 1 ≤ diaglength(S)
+    if newdim + 1 ≤ full_dim
       val =
         singular_values ? getdiagindex(S, newdim + 1)^2 : abs(getdiagindex(S, newdim + 1))
     end
   end
+  (newdim >= min_blockdim) || (newdim = min_blockdim)
   return newdim
 end
 
@@ -31,7 +34,7 @@ computed from the dense svds of seperate blocks.
 """
 function LinearAlgebra.svd(T::BlockSparseMatrix{ElT}; kwargs...) where {ElT}
   alg::String = get(kwargs, :alg, "divide_and_conquer")
-
+  min_blockdim::Int = get(kwargs, :min_blockdim, 0)
   truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)
 
   #@timeit_debug timer "block sparse svd" begin
@@ -74,7 +77,9 @@ function LinearAlgebra.svd(T::BlockSparseMatrix{ElT}; kwargs...) where {ElT}
   if truncate
     truncerr, docut = truncate!(d; kwargs...)
     for n in 1:nnzblocks(T)
-      blockdim = _truncated_blockdim(Ss[n], docut; singular_values=true, truncate=truncate)
+      blockdim = _truncated_blockdim(
+        Ss[n], docut; min_blockdim, singular_values=true, truncate
+      )
       if blockdim == 0
         push!(dropblocks, n)
       else

src/decomp.jl

@@ -78,6 +78,7 @@ Utrunc2, Strunc2, Vtrunc2 = svd(A, i, k; cutoff=1e-10);
   - `"recursive"` - ITensor's custom svd. Very reliable, but may be slow if high precision is needed. To get an `svd` of a matrix `A`, an eigendecomposition of ``A^{\\dagger} A`` is used to compute `U` and then a `qr` of ``A^{\\dagger} U`` is used to compute `V`. This is performed recursively to compute small singular values.
 - `use_absolute_cutoff::Bool = false`: set if all probability weights below the `cutoff` value should be discarded, rather than the sum of discarded weights.
 - `use_relative_cutoff::Bool = true`: set if the singular values should be normalized for the sake of truncation.
+- `min_blockdim::Int = 0`: for SVD of block-sparse or QN ITensors, require that the number of singular values kept be greater than or equal to this value when possible
 
 See also: [`factorize`](@ref), [`eigen`](@ref)
 """

src/exports.jl

@@ -381,6 +381,7 @@ export
   val,
 
   # qn/qnindex.jl
+  blockdim,
   flux,
   hasqns,
   nblocks,

src/physics/sitetype.jl

@@ -693,7 +693,7 @@ siteind(tag::String, n; kwargs...) = siteind(SiteType(tag), n; kwargs...)
 # Special case of `siteind` where integer (dim) provided
 # instead of a tag string
 #siteind(d::Integer, n::Integer; kwargs...) = Index(d, "Site,n=$n")
-function siteind(d::Integer, n::Integer; addtags="", kwargs...) 
+function siteind(d::Integer, n::Integer; addtags="", kwargs...)
   return Index(d, "Site,n=$n, $addtags")
 end
 

src/qn/qnindex.jl

@@ -168,6 +168,18 @@ end
 
 dim(i::QNIndex) = dim(space(i))
 
+"""
+    nblocks(i::QNIndex)
+
+Returns the number of QN blocks, or subspaces, of the QNIndex `i`.
+
+### Example
+```
+julia> i = Index([QN("Sz",-1)=>2, QN("Sz",0)=>4, QN("Sz",1)=>2], "i")
+julia> nblocks(i)
+3
+```
+"""
 nblocks(i::QNIndex) = nblocks(space(i))
 # Define to be 1 for non-QN Index
 nblocks(i::Index) = 1

test/decomp.jl

@@ -119,6 +119,30 @@ using ITensors, LinearAlgebra, Test
 
     @test flux(F.Vt) == QN("Sz", 0)
   end
+
+  @testset "SVD block_mindim keyword" begin
+    i = Index(
+      [
+        QN("Sz", 4) => 1,
+        QN("Sz", 2) => 4,
+        QN("Sz", 0) => 6,
+        QN("Sz", -2) => 4,
+        QN("Sz", -4) => 1,
+      ],
+      "i",
+    )
+    j = sim(i)
+    X = randomITensor(QN("Sz", 0), i, j)
+
+    min_blockdim = 2
+    U, S, V = svd(X, i; cutoff=1E-1, min_blockdim)
+    u = commonind(S, U)
+
+    @test nblocks(u) == nblocks(i)
+    for b in 1:nblocks(u)
+      @test blockdim(u, b) == blockdim(i, b) || blockdim(u, b) >= min_blockdim
+    end
+  end
 end
 
 nothing