How to replace SVD in ITensors by random SVD

ITensor Julia Questions
1

Hello everone!

I just want to ask how can I construct my own SVD algorithms (for example random singular value decomposition) to deal with itensors with QN. And then impose rsvd algorithm to TDVP or anywhere truncation is needed.

A typical random singular value decomposition was posted below

rsvd
rsvd975×408 10.9 KB

rsvd_1
rsvd_11192×649 19.1 KB

I have already realize a simpled rsvd function to fractorize itensors. It works well paritially. In such a method, one needs to make a normal Gaussian distributed random matrix with a predetermined dimension k (or the rank of the fractorized matrix to target), thus one has to define an index with dimension k to generate the random matrix.

But here is the problem: most systems have a conservative quantity, such that indices of tensors are QN indices. The storage type of tensors with QN are BlockSparse. But the random matrix are NDTensors.Dense type with one index with QN and another index without QN. The way to contract NDTensors.Dense tensors with NDTensors.BlockSparse is to dense() the BlockSparse tensors. And finally, the fractorized tensors are Dense type, and error occurs when call the orthogonalize!(psi,i) due the same reason of BlockSparse tensors and Dense tensors (they can not be contracted).

Here is my code for random SVD

using ITensors
using LinearAlgebra

function power_iteration(A::ITensor,omega::ITensor,q::Int)

    min_i=commonind(A,omega)
    Y=A*omega
    max_i=commonind(A,Y)
    l_i=commonind(Y,omega)

    for i=1:q
        Q_1,_=qr(Y,max_i)
        replaceind!(Q_1,inds(Q_1)[2],l_i)
        Q_1,_=qr(permute(A,(uniqueind(A,Q_1),commonind(A,Q_1)))*Q_1,min_i)
        replaceind!(Q_1,inds(Q_1)[2],l_i)
        Y=A*Q_1
    end

    return Y
end

function rsvd(A::ITensor,
    Linds...;kwargs...)

    maxdim=get(kwargs,:maxdim,nothing)
    q=get(kwargs,:iter,0)
    p=get(kwargs,:oversam,0)

    Lis=commoninds(A,Linds...)
    Ris=uniqueinds(A,Lis...)

    CL=combiner(Lis)
    CR=combiner(Ris)

    AC=dense(A * CR * CL)

    cL=combinedind(CL)
    cR=combinedind(CR)

    if inds(AC) != (cL, cR)
        AC = permute(AC, cL, cR)
    end

    if (maxdim==nothing || maxdim >= min(dim(cL),dim(cR)))
        maxdim=min(dim(cL),dim(cR))
    end

    k=maxdim

    l=k+p

    l_i=Index(l,"omega_l")

    if dim(cL)>=dim(cR)
        omega=ITensor(randn(dim(cR),dim(l_i)),(dag(cR),l_i))
        
        if q==0
            Y=AC*omega
        else
            Y=power_iteration(AC,omega,q)
        end

        if inds(Y)!=(cL,l_i)
            Y=permute(Y,(cL,l_i))
        end

        Q,_=qr(Y,cL)
        replaceind!(Q,inds(Q)[2],l_i)

        U_t,S,V=svd(AC*Q,l_i;maxdim=k)
        U=dag(Q*U_t)*CL
        V=dag(V)*CR

        return (U,S,V)

    elseif dim(cL)<dim(cR)
        omega=ITensor(randn(dim(l_i),dim(cL)),(l_i,dag(cL)))

        if q==0
            Y=omega*AC
        else
            Y=power_iteration(AC,omega,q)
        end

        if inds(Y)!=(l_i,cR)
            Y=permute(Y,(l_i,cR))
        end

        Q,_=qr(Y,cR)
        replaceind!(Q,inds(Q)[2],l_i)

        U,S,V_t=svd(AC*Q,cL;maxdim=k)
        Q=permute(Q,(l_i,cR))
        V=dag(V_t*Q)*CR
        U=dag(U)*CL

        return (U,S,V)
    end
end

It works well for normal itensors, but errors occur for itensors with QN because of random matrix omega=ITensor(randn(dim(cR),dim(l_i)),(dag(cR),l_i))

Any help would be appreciate.

2

Hello Wang,

I am familiar with this approach and these kinds of problem in ITensor. It looks like the problem is that you are constructing l_i as a dense index. So when the input is a BlockSparseTensor you are mixing indices

cR = (dim=2|id=315|"CMB,Link") <Out>
 1: QN() => 2
l_i = (dim=2|id=11|"omega_l")

Linear and multilinear algebra is one of my particular research interest so I am currently looking into making randomized linear algebra more available in the backend of ITensors, this will, in principle, eliminate these kinds of issues.

Karl