Using ITensorMPS to efficiently calculate all principal minors of a matrix

DMRG and Numerical Methods
,
1

Hi!

I come from statistical and quantum physics and try to calculate all 2^M principle minors (PMs) of a real M\times M matrix G using the method described by Griffin & Tsatsomeros, see also this mathoverflow posting. Below first find a direct calculation of all PMs using this method:

using ITensors, ITensorMPS, LinearAlgebra
using ITensors:ITensor # explicit import, are we redefine ITensor
function ITensor(mps::MPS)::ITensor
    T = prod(mps)
    return permute(T, sortinds(T))
end;
function sortinds(G::ITensor)::NTuple
    return (
        inds(G,tags="G",plev=1)...,
        inds(G,tags="G",plev=0)...,
        inds(G,tags="Link")...,
        inds(G,tags="Site")...
    )
end;
function myQ(M::Int)::Matrix{Float64}
    return reshape([1/(1+abs(i-j)) for i in 1:M for j in 1:M], (M,M))
end;
function myG(M::Int)::Matrix{Float64}
    Q = myQ(M)
    return inv(I + Q) * (I - Q)
end;
function PMs(G::Matrix{Float64})::Vector
    M = size(G,1)
    G = reshape(G,(M,M,1)); F = [1.0]
    [ let; if m > 1
            F = vcat(F, F .* G[1,1,:])
            G = cat(
                G[2:end,2:end,:], # left leafs
                G[2:end,2:end,:] - (G[2:end,1:1,:] .* G[1:1,2:end,:]) ./
                    G[1:1,1:1,:] # right leafs
            ; dims=3)
        end
        (G=G, F=F)
    end for m in 1:M ]
end;

function pm(G::Matrix, i::Tuple)
    return det(G[[i...],[i...]])
end;
MM = 3
G = myG(MM)
GF = PMs(G)
@show [(size(gf.G),size(gf.F)) for gf in GF]
hcat(
    vcat(1.0,(vcat([GF[m].F .* GF[m].G[1,1,:] for m in 1:MM]...))),
    [pm(G,()), pm(G,(1,)), pm(G,(2,)), pm(G,(1,2)), 
     pm(G,(3,)), pm(G,(1,3)), pm(G,(2,3)), pm(G,(1,2,3))]
)

# result:
8×2 Matrix{Float64}:
  1.0     1.0
  0.08    0.08
  0.12    0.12
 -0.048  -0.048
  0.08    0.08
 -0.008  -0.008
 -0.048  -0.048
 -0.024  -0.024

Note that PMs() calculates a list of rows of a binary tree, where in each step m the resulting 2^m matrices G (and corresponding factors F) are doubled by using either the submatrix with first row & column deleted (left leaf) or the (rank-1) Schur complement w.r.t. the (1,1)-element (right leaf). The desired results are encoded in the (1,1)-elements of all G\text s as well as in the F\text s.

Using ITensorMPS, I now try to limit this exponential growth, being similar to the Hilbert space growth in system size M for an s=1/2 chain.
Some notes:

  • States "L" and "R" stand for the left and right leafs in the btree, respectively.
  • As I could not figure out how to do the (parallel) Schur complement calculation in ITensor, I did it with usual Julia arrays and broadcasting.
  • I focus on the calculation of G.
  • I also redefined PMs() and added an option method, see below.
ITensors.disable_warn_order()
ITensors.space(::SiteType"s") = 2 # 'Ising' spin
ITensors.state(::StateName"L", ::SiteType"s") = [1.0, 0.0]; sL = 1
ITensors.state(::StateName"R", ::SiteType"s") = [0.0, 1.0]; sR = 2
ITensors.state(::StateName"+", ::SiteType"s") = [1.0, 1.0]

Mmax = 128
s    = siteinds("s", Mmax); # fixed sites for simplicity

function PMs(G::Matrix{Float64}; method=1)::Vector
    # added method flag: 1=exact Schur complement, 2=something linear in G
    M = size(G,1)
    G = reshape(G,(M,M,1)); F = [1.0]
    [ let; if m > 1
            F = vcat(F, F .* G[1,1,:])
            G = cat(
                G[2:end,2:end,:],
                if method==1
                    G[2:end,2:end,:] - (G[2:end,1:1,:] .* G[1:1,2:end,:]) ./
                                        G[1:1,1:1,:]
                elseif method==2
                    G[2:end,2:end,:] - (G[2:end,1:1,:] .* 
                    cumsum(1 .+ 0 .* G[1:1,2:end,:], dims=1))
                end
            ; dims=3)
        end
        (G=G, F=F)
    end for m in 1:M ]
end;
function PMs_MPS(Gmps::MPS, m::Int; cutoff=0, method=1)::MPS
    # start MPS: Gi' x l1 x Gi x l2 x s1 x l3 x s2 x l4 x ...
    # idx      : m          m+1       m+2       m+3
    # method: 1=exact Schur complement, 2=something linear in G
    @show m
    s = siteinds(Gmps); #@show s
    orthogonalize!(Gmps, m)
    Rmps = copy(Gmps) # result
    G = Gmps[m] * Gmps[m+1] * Gmps[m+2]
    i = noprime(inds(Gmps[m],tags="G")[1])
    M = dim(i) # actual size of G
    j = Index(M-1; tags = "G,i", plev = 0)
    sm = s[m+2]; # pos of site m
    S = ITensor(diagm(M-1,M, 1 => ones(M-1)), j,i)
    GL = G * onehot(sm => sL);
    iL = sortinds(GL)
    GLa = reshape(Array(GL,iL), (M,M,:))
    if method==1
        GR = ITensor(GLa - GLa[:,1:1,:] .* GLa[1:1,:,:] ./
                           GLa[1:1,1:1,:], iL)
    elseif method==2
        GR = ITensor(GLa - GLa[:,1:1,:] .* 
                     cumsum(1 .+ 0 .* GLa[1:1,:,:]; dims=1), iL)
    end
    G = S' * (GL * onehot(sm => sL) + GR * onehot(sm => sR)) * S
    # redistribute updated G with SVD: Rmps[m:m+2] = (L1*S12, L2*S23, L3)
    (L1,S12,L23) = svd(G, inds(G, tags="Link,l=$(m-1)"), sm; 
        cutoff, lefttags ="Link,left", righttags="Link,l=$(m)");
    (L2,S23,L3) = svd(L23, inds(L23, tags="Link,l=$(m)"), j'; 
        cutoff, lefttags ="Link,left", righttags="Link,l=$(m+1)");
    # put back
    Rmps[m  ] = L1 * S12
    Rmps[m+1] = L2 * S23
    Rmps[m+2] = L3
    return Rmps
end;

Starting from the left end of the tensor train GM, in each iteration (i.e. call to PMs_MPS()) the matrix G with indices i',i is moved one step right in the tensor train,

( i',i,s_1,...,s_{M-1}) \mapsto (s_1,i',i,s_2,...,s_{M-1}) \mapsto \ldots,

while the size (dimension) of indices i',i of G are reduced by one in each step.

MM = 10; mycutoff=1e-16; method=2
GF = PMs(myG(MM); method);
ii = Index(MM; tags = "G,i", plev = 0)
GT = ITensor(myG(MM), ii',ii) * ITensor(1.0, s[1:MM-1])
GM = [ MPS(GT, inds(GT), cutoff=mycutoff) ]; #@show GM
# calculate
for m in 1:MM-1
    push!(GM, PM_MPS(GM[end], m; cutoff=mycutoff, method))
end;
# show results
for m in 1:MM
    GMT = ITensor(GM[m])
    jj = inds(GMT, tags="G",plev=0)
    @show (m,
        ≈(GMT, ITensor(GF[m].G, jj',jj,s[1:m-1]) * ITensor(1.0, s[m:MM-1]); 
          rtol=1e-8),linkdims(GM[m]));
end;

Now my problem: With method=1 the exact algorithm is used, which is nonlinear in G (due to the Schur complement), and the MPS factorisation does not work after the second step, where the factored out (left-orthogonal) part becomes important:

# output for MM = 10; mycutoff=1e-16; method=1
(...) = (1, true, [10, 1, 1, 1, 1, 1, 1, 1, 1, 1])
(...) = (2, true, [2, 9, 1, 1, 1, 1, 1, 1, 1, 1])
(...) = (3, false, [2, 4, 8, 1, 1, 1, 1, 1, 1, 1])
(...) = (4, false, [2, 4, 7, 7, 1, 1, 1, 1, 1, 1])
(...) = (5, false, [2, 4, 7, 10, 6, 1, 1, 1, 1, 1])
(...) = (6, false, [2, 4, 7, 10, 13, 5, 1, 1, 1, 1])
(...) = (7, false, [2, 4, 7, 10, 13, 10, 4, 1, 1, 1])
(...) = (8, false, [2, 4, 7, 10, 13, 10, 6, 3, 1, 1])
(...) = (9, false, [2, 4, 7, 10, 13, 10, 6, 3, 2, 1])
(...) = (10, false, [2, 4, 7, 10, 13, 10, 6, 3, 1, 1])

Changing the Schur complement update to some linear function in G (method=2), the MPS dimensional reduction works as expected:

# output for MM = 10; mycutoff=1e-16; method=2
(...) = (1, true, [10, 1, 1, 1, 1, 1, 1, 1, 1, 1])
(...) = (2, true, [2, 9, 1, 1, 1, 1, 1, 1, 1, 1])
(...) = (3, true, [2, 3, 8, 1, 1, 1, 1, 1, 1, 1])
(...) = (4, true, [2, 3, 4, 7, 1, 1, 1, 1, 1, 1])
(...) = (5, true, [2, 3, 4, 5, 6, 1, 1, 1, 1, 1])
(...) = (6, true, [2, 3, 4, 5, 6, 5, 1, 1, 1, 1])
(...) = (7, true, [2, 3, 4, 5, 6, 5, 4, 1, 1, 1])
(...) = (8, true, [2, 3, 4, 5, 6, 5, 4, 3, 1, 1])
(...) = (9, true, [2, 3, 4, 5, 6, 5, 4, 3, 2, 1])
(...) = (10, true, [2, 3, 4, 5, 6, 5, 4, 3, 1, 1])

Question:

  • Any idea how to implement rank-1 updates with Schur complements using MPS?
  • How can I avoid or handle the non-linearity?

Best, Fred