[NDTensors] [ITensors] Implement optional rank reduction for QR/RQ/QL/LQ decompositions

NDTensors/src/linearalgebra/linearalgebra.jl

@@ -366,11 +366,58 @@ function LinearAlgebra.eigen(
   V = complex(tensor(Dense(vec(VM)), Vinds))
   return D, V, spec
 end
+#
+#  Trim out zero rows of R/X within tolerance rr_cutoff. Also trim the corresponding columns
+#  of Q.  X = R or L
+#
+function trim_rows(
+  Q::AbstractMatrix, X::AbstractMatrix, rr_cutoff::Float64; rr_verbose=false, kwargs...
+)
+  #
+  #  Find and count the zero rows.  Bail out if there are none.
+  #
+  Xnr, Xnc = size(X)
+  Qnr, Qnc = size(Q)
+  @assert Xnr == Qnc #Sanity check.
+  zeros = map((r) -> (maximum(abs.(X[r, 1:Xnc])) <= rr_cutoff), 1:Xnr)
+  num_zero_rows = sum(zeros)
+  if num_zero_rows == 0
+    return Q, X
+  end
+  #
+  # Useful output for trouble shooting.
+  #
+  if rr_verbose
+    println(
+      "Rank Reveal removing $num_zero_rows rows with log10(rr_cutoff)=$(log10(rr_cutoff))"
+    )
+  end
+  #
+  # Create new Q & X martrices with reduced size.
+  #
+  X1nr = Xnr - num_zero_rows #new dim between Q & X
+  T = eltype(X)
+  X1 = Matrix{T}(undef, X1nr, Xnc)
+  Q1 = Matrix{T}(undef, Qnr, X1nr)
+  #
+  #  Transfer non-zero rows of X and corresponding columns of Q.
+  #
+  r1 = 1 #Row/col counter in new reduced Q & X
+  for r in 1:Xnr
+    if zeros[r] == false
+      X1[r1, :] = X[r, :] #transfer row
+      Q1[:, r1] = Q[:, r] #transfer column
+      r1 += 1 #next row in rank reduced matrices.
+    end #if zero
+  end #for r
+  return Q1, X1
+end
 
 function qr(T::DenseTensor{<:Any,2}; positive=false, kwargs...)
   qxf = positive ? qr_positive : qr
   return qx(qxf, T; kwargs...)
 end
+
 function ql(T::DenseTensor{<:Any,2}; positive=false, kwargs...)
   qxf = positive ? ql_positive : ql
   return qx(qxf, T; kwargs...)
@@ -380,17 +427,27 @@ end
 #  Generic function for qr and ql decomposition of dense matrix.
 #  The X tensor = R or L.
 #
-function qx(qx::Function, T::DenseTensor{<:Any,2}; kwargs...)
-  QM, XM = qx(matrix(T))
-  # Be aware that if positive==false, then typeof(QM)=LinearAlgebra.QRCompactWYQ, not Matrix
-  # It gets converted to matrix below.
-  # Make the new indices to go onto Q and R
-  q, r = inds(T)
-  q = dim(q) < dim(r) ? sim(q) : sim(r)
-  IndsT = indstype(T) #get the index type
+function qx(qx::Function, T::DenseTensor{<:Any,2}; rr_cutoff=-1.0, kwargs...)
+  QM1, XM = qx(matrix(T))
+  # When qx=qr typeof(QM1)==LinearAlgebra.QRCompactWYQ
+  # When qx=ql typeof(QM1)==Matrix and this should be a no-op
+  QM = Matrix(QM1)
+  #
+  #  Do row removal for rank revealing QR/QL.  Probably not worth it to elminate the if statement
+  #
+  if rr_cutoff >= 0.0
+    QM, XM = trim_rows(QM, XM, rr_cutoff; kwargs...)
+  end
+  #
+  # Make the new indices to go onto Q and X
+  #
+  IndsT = indstype(T) #get the indices type
+  @assert IndsT.parameters[1] == IndsT.parameters[2] #they better be the same!
+  IndexT = IndsT.parameters[1] #establish the single index type.
+  q = IndexT(size(XM)[1]) #create the Q--X link index.
   Qinds = IndsT((ind(T, 1), q))
   Xinds = IndsT((q, ind(T, 2)))
-  Q = tensor(Dense(vec(Matrix(QM))), Qinds) #Q was strided
+  Q = tensor(Dense(vec(QM)), Qinds)
   X = tensor(Dense(vec(XM)), Xinds)
   return Q, X
 end
@@ -449,7 +506,7 @@ function ql_positive(M::AbstractMatrix)
 end
 
 #
-#  Lapack replaces A with Q & R carefully packed together.  So here we just copy a
+#  Lapack replaces A with Q & L carefully packed together.  So here we just copy a
 #  before letting lapack overwirte it. 
 #
 function ql(A::AbstractMatrix; kwargs...)

NDTensors/test/linearalgebra.jl

@@ -2,6 +2,10 @@ using NDTensors
 using LinearAlgebra
 using Test
 
+# Not available on CI machine that test NDTensors.
+# using Random
+# Random.seed!(314159)
+
 @testset "random_orthog" begin
   n, m = 10, 4
   O1 = random_orthog(n, m)
@@ -20,58 +24,73 @@ end
   @test norm(U2 * U2' - Diagonal(fill(1.0, m))) < 1E-14
 end
 
-@testset "Dense $qx decomposition, elt=$elt, positve=$positive, singular=$singular" for qx in
-                                                                                        [
+@testset "Dense $qx decomposition, elt=$elt, positve=$positive, singular=$singular, rank_reveal=$rank_reveal" for qx in
+                                                                                                                  [
     qr, ql
   ],
   elt in [Float64, ComplexF64, Float32, ComplexF32],
   positive in [false, true],
-  singular in [false, true]
+  singular in [false, true],
+  rank_reveal in [false, true],
 
-  eps = Base.eps(real(elt)) * 30 #this is set rather tight, so if you increase/change m,n you may have open up the tolerance on eps.
+  eps in Base.eps(real(elt)) * 30
+  #this is set rather tight, so if you increase/change m,n you may have open up the tolerance on eps.
+  rr_cutoff = rank_reveal ? eps * 1.0 : -1.0
   n, m = 4, 8
-  Id = Diagonal(fill(1.0, min(n, m)))
   #
   # Wide matrix (more columns than rows)
   #
   A = randomTensor(elt, (n, m))
-  # We want to test 0.0 on the diagonal.  We need make all roaw equal to gaurantee this with numerical roundoff.
+  # We want to test 0.0 on the diagonal.  We need make all rows linearly dependent 
+  # gaurantee this with numerical roundoff.
   if singular
-    for i in 2:n
-      A[i, :] = A[1, :]
+    for i in 2:3
+      A[i, :] = A[1, :] * 1.05^n
     end
   end
-  Q, X = qx(A; positive=positive) #X is R or L.
+  # you can set rr_verbose=true if you want to get debug output on rank reduction.
+  Q, X = qx(A; positive=positive, rr_cutoff=rr_cutoff, rr_verbose=false) #X is R or L. 
   @test A ≈ Q * X atol = eps
-  @test array(Q)' * array(Q) ≈ Id atol = eps
-  @test array(Q) * array(Q)' ≈ Id atol = eps
+  @test array(Q)' * array(Q) ≈ Diagonal(fill(1.0, dim(Q, 2))) atol = eps
+  if dim(Q, 1) == dim(Q, 2)
+    @test array(Q) * array(Q)' ≈ Diagonal(fill(1.0, min(n, m))) atol = eps
+  end
   if positive
     nr, nc = size(X)
     dr = qx == ql ? Base.max(0, nc - nr) : 0
     diagX = diag(X[:, (1 + dr):end]) #location of diag(L) is shifted dr columns over the right.
     @test all(real(diagX) .>= 0.0)
     @test all(imag(diagX) .== 0.0)
   end
+  if rr_cutoff > 0 && singular
+    @test dim(Q, 2) == 2 #make sure the rank revealing mechanism hacked off the columns of Q (and rows of X).
+    @test dim(X, 1) == 2 #Redundant?
+  end
   #
   # Tall matrix (more rows than cols)
   #
   A = randomTensor(elt, (m, n)) #Tall array
   # We want to test 0.0 on the diagonal.  We need make all rows equal to gaurantee this with numerical roundoff.
   if singular
-    for i in 2:m
+    for i in 2:4
       A[i, :] = A[1, :]
     end
   end
-  Q, X = qx(A; positive=positive)
+  Q, X = qx(A; positive=positive, rr_cutoff=rr_cutoff, rr_verbose=false)
   @test A ≈ Q * X atol = eps
-  @test array(Q)' * array(Q) ≈ Id atol = eps
+  @test array(Q)' * array(Q) ≈ Diagonal(fill(1.0, dim(Q, 2))) atol = eps
+  #@test array(Q) * array(Q)' no such relationship for tall matrices.
   if positive
     nr, nc = size(X)
     dr = qx == ql ? Base.max(0, nc - nr) : 0
     diagX = diag(X[:, (1 + dr):end]) #location of diag(L) is shifted dr columns over the right.
     @test all(real(diagX) .>= 0.0)
     @test all(imag(diagX) .== 0.0)
   end
+  if rr_cutoff > 0 && singular
+    @test dim(Q, 2) == 4 #make sure the rank revealing mechanism hacked off the columns of Q (and rows of X).
+    @test dim(X, 1) == 4 #Redundant?
+  end
 end
 
 nothing

test/base/test_decomp.jl

@@ -1,4 +1,6 @@
-using ITensors, LinearAlgebra, Test
+using ITensors, LinearAlgebra, Test, Random
+
+Random.seed!(314159)
 
 #
 #  Decide if rank 2 tensor is upper triangular, i.e. all zeros below the diagonal.
@@ -67,6 +69,30 @@ function is_upper(A::ITensor, r::Index)::Bool
 end
 is_lower(A::ITensor, r::Index)::Bool = is_upper(r, A)
 
+#
+#  Makes all columns lineary depenedent but scaled differently.
+#
+function rank_fix(A::ITensor, Linds...)
+  Lis = commoninds(A, (Linds...))
+  Ris = uniqueinds(A, Lis)
+  #
+  #  Use combiners to render A down to a rank 2 tensor ready matrix QR routine.
+  #
+  CL, CR = combiner(Lis...), combiner(Ris...)
+  cL, cR = combinedind(CL), combinedind(CR)
+  AC = A * CR * CL
+  if inds(AC) != IndexSet(cL, cR)
+    AC = permute(AC, cL, cR)
+  end
+  At = NDTensors.tensor(AC)
+  nc = dim(At, 2)
+  @assert nc >= 2
+  for c in 2:nc
+    At[:, c] = At[:, 1] * 1.05^c
+  end
+  return itensor(At) * dag(CL) * dag(CR)
+end
+
 function diag_upper(l::Index, A::ITensor)
   At = NDTensors.tensor(A * combiner(noncommoninds(A, l)...))
   if size(At) == (1,)
@@ -362,6 +388,38 @@ end
     @test A ≈ Q * L atol = 1e-13
   end
 
+  @testset "Rank revealing QR/RQ/QL/LQ decomp on MPS dense $elt tensor" for ninds in
+                                                                            [1, 2, 3],
+    elt in [Float64, ComplexF64]
+
+    l = Index(5, "l")
+    s = Index(2, "s")
+    r = Index(5, "r")
+    A = randomITensor(elt, l, s, s', r)
+
+    Ainds = inds(A)
+    A = rank_fix(A, Ainds[1:ninds]) #make all columns linear dependent on column 1, so rank==1.
+    Q, R, q = qr(A, Ainds[1:ninds]; rr_cutoff=1e-12)
+    @test dim(q) == 1 #check that we found rank==1
+    @test A ≈ Q * R atol = 1e-13
+    @test Q * dag(prime(Q, q)) ≈ δ(Float64, q, q') atol = 1e-13
+
+    R, Q, q = rq(A, Ainds[1:ninds]; rr_cutoff=1e-12)
+    @test dim(q) == 1 #check that we found rank==1
+    @test A ≈ Q * R atol = 1e-13 #With ITensors R*Q==Q*R
+    @test Q * dag(prime(Q, q)) ≈ δ(Float64, q, q') atol = 1e-13
+
+    L, Q, q = lq(A, Ainds[1:ninds]; rr_cutoff=1e-12)
+    @test dim(q) == 1 #check that we found rank==1
+    @test A ≈ Q * L atol = 1e-13 #With ITensors L*Q==Q*L
+    @test Q * dag(prime(Q, q)) ≈ δ(Float64, q, q') atol = 1e-13
+
+    Q, L, q = ql(A, Ainds[1:ninds]; rr_cutoff=1e-12)
+    @test dim(q) == 1 #check that we found rank==1
+    @test A ≈ Q * L atol = 1e-13
+    @test Q * dag(prime(Q, q)) ≈ δ(Float64, q, q') atol = 1e-13
+  end
+
   @testset "factorize with QR" begin
     l = Index(5, "l")
     s = Index(2, "s")