Add `inv(::ITensor)`

[mtfishman]

Add a function [p]inv(::ITensor, Linds; kwargs...) to invert an ITensor using a Moore-Penrose pseudoinverse. We can use tolerance keyword arguments for the pseudoinverse similar to the ones used in LinearAlgebra.pinv.

A user can specify the indices to specify how the ITensor gets reshaped into a matrix, and if they aren't specified, it can use pairs of primed and unprimed indices (treating it like a square matrix from the primed to the unprimed space, like in eigen). The latter case is useful for doing non-orthogonal basis changes with gates, i.e. do a similarity transformation. This could be used in a special mode of the apply function, with an interface apply(gates, ::MPO; apply_inv = true), analogous to the current apply(gates, ::MPO; apply_dag = true) but applicable to similarity transformations of non-unitary gates.

[tipfom]

Dear @mtfishman and @kshyatt, is this still something that you want to see in the ITensor libary or should this be postponed for work on the underlying backend #1250? I need this functionality for non-Hermitian dmrg I implemented in ITensor and I'm currently using the following implementation:

function invertitensor(A::ITensor, finv, Linds...)
    Lis = commoninds(A, ITensors.indices(Linds...))
    Ris = uniqueinds(A, Lis)

    Cr = combiner(Ris...)
    Cl = combiner(Lis...)

    A = A * Cr * Cl

    Minv = nothing

    return if hasqns(A)
        Minv = deepcopy(A)
        for b in eachnzblock(A)
            Amat = matrix(A[b])
            Amatinv = finv(Amat)

            Minv[b] .= adjoint(Amatinv)
        end
        Minv * dag(Cr) * dag(Cl)
    else
        M = matrix(A, combinedind(Cr), combinedind(Cl))
        Minv = itensor(adjoint(finv(M)), combinedind(Cr), combinedind(Cl); tol=1e-8)
        Minv * dag(Cr) * dag(Cl)
    end
end

function LinearAlgebra.inv(A::ITensor, Linds...)
    return invertitensor(A, inv, Linds...)
end

function LinearAlgebra.pinv(A::ITensor, Linds...; kwargs...)
    return invertitensor(A, x -> pinv(x; kwargs...), Linds...)
end

function LinearAlgebra.inv(A::ITensor)
    Ris = filterinds(A; plev=0)
    Lis = Ris'
    return LinearAlgebra.inv(A, Lis)
end

function LinearAlgebra.pinv(A::ITensor; kwargs...)
    Ris = filterinds(A; plev=0)
    Lis = Ris'
    return LinearAlgebra.pinv(A, Lis; kwargs...)
end

It of course lacks the input validation needed for proper integration into the library, if you have any suggestions on how to proceed with this I'd be very willing to spend some time providing a more robust implementation. The current implementation returns a ITensor with the same indices. An exemplary call is, e.g.,

i = Index(20)
j = Index(20)
A = random_itensor(ComplexF64, j', i)
Ainv = pinv(deepcopy(A), j')

M = dag(replaceind(Ainv, j'=>j'')) * A

such that M is approximately the identity operator. Indices with quantum numbers are also supported.

Thanks!

[mtfishman]

I think it is better to do this as part of the rewrite.