Jutho / TensorKitManifolds.jl

Useful tools for working with isometric or unitary tensors
MIT License
18 stars 6 forks source link

Optimization of unitary matrix with non-degenerate symmetry blocks #9

Open philihps opened 1 year ago

philihps commented 1 year ago

Hi Jutho,

I am trying to do an optimization of a two-mode unitary disentangling transformation, based on a cost function. Below I provide a MWE with a test cost function. Setting d = 2 works, but setting d = 1 fails. Is this expected behaviour or can it be caught?

using OptimKit
using TensorKit
using TensorKitManifolds

function fg(U::TensorMap)
    fval = real(tr(U * U'));
    gval = Stiefel.project!(2 * U', U; metric = :euclidean)
    return fval, gval;
end

# set up two-mode system
d = 2;
PL = ComplexSpace(d);
PR = ComplexSpace(d);

# initialize unitary
U = TensorMap(randhaar, Float64, PL * PR, PL * PR);

optimAlg = LBFGS(12, verbosity = 2);
res = optimize(fg, U, optimAlg;
        inner = Stiefel.inner,
        retract = Stiefel.retract,
        transport! = Stiefel.transport!,
        isometrictransport = true
    );

Optimizing a 1x1 unitary matrix is of course meaningless. In the real setting, however, the two-mode unitary has U1 quantum numbers, so it consists of different symmetry blocks. However, they are not necessarily degenerate, and it would be nice to be able to handle degenerate and non-degenerate sectors.

lkdvos commented 1 year ago

Hi Philip,

Did you get this fixed in the end? If not, can you maybe send me something with which I can reproduce the error?

philihps commented 1 year ago

Hi Lukas,

I didn't look further into this, but I can provide an example for which the error appears. Using the code

using OptimKit
using TensorKit
using TensorKitManifolds

function fg(U::TensorMap)
    fval = real(tr(U * U'));
    gval = Stiefel.project!(2 * U', U; metric = :euclidean)
    return fval, gval;
end

# set up two-mode system
PL = U1Space(0 => 1, -4 => 1, -8 => 1);
PR = U1Space(0 => 1, -3 => 1, -6 => 1);

# initialize unitary
U = TensorMap(randhaar, Float64, PL * PR, PL * PR);

# optimize unitary
optimAlg = LBFGS(12, verbosity = 2);
res = optimize(
        fg, U, optimAlg;
        inner = Stiefel.inner,
        retract = Stiefel.retract,
        transport! = Stiefel.transport!,
        isometrictransport = true
    );
optimU, optimCostFunc, normGrad, normGradHistory = res;

I get an error: "ArgumentError: matrix contains Infs or NaNs". However, this only happens for U = TensorMap(randhaar, Float64, PL * PR, PL * PR);, and it works fine for U = TensorMap(randhaar, ComplexF64, PL * PR, PL * PR);.