|
| 1 | +""" |
| 2 | + struct BPGauge |
| 3 | +
|
| 4 | +Algorithm for gauging PEPS with belief propagation fixed point messages. |
| 5 | +""" |
| 6 | +@kwdef struct BPGauge |
| 7 | + # TODO: add options |
| 8 | +end |
| 9 | + |
| 10 | +""" |
| 11 | +$(SIGNATURES) |
| 12 | +
|
| 13 | +Fix the gauge of `psi` using fixed point environment `env` of belief propagation. |
| 14 | +""" |
| 15 | +function gauge_fix(psi::InfinitePEPS, alg::BPGauge, env::BPEnv) |
| 16 | + psi′ = copy(psi) |
| 17 | + XXinv = map(eachcoordinate(psi, 1:2)) do I |
| 18 | + _, X, Xinv = _bp_gauge_fix!(CartesianIndex(I), psi′, env) |
| 19 | + return X, Xinv |
| 20 | + end |
| 21 | + return psi′, XXinv |
| 22 | +end |
| 23 | + |
| 24 | +function _sqrt_bp_messages(I::CartesianIndex{3}, env::BPEnv) |
| 25 | + dir, row, col = Tuple(I) |
| 26 | + @assert dir == NORTH || dir == EAST |
| 27 | + M12 = env[dir, dir == NORTH ? _prev(row, end) : row, dir == EAST ? _next(col, end) : col] |
| 28 | + sqrtM12, isqrtM12 = sqrt_invsqrt(twist(M12, 1)) |
| 29 | + M21 = env[dir + 2, row, col] |
| 30 | + sqrtM21, isqrtM21 = sqrt_invsqrt(M21) |
| 31 | + return sqrtM12, isqrtM12, sqrtM21, isqrtM21 |
| 32 | +end |
| 33 | + |
| 34 | +""" |
| 35 | + _bp_gauge_fix!(I, psi::InfinitePEPS, env::BPEnv) -> psi, X, X⁻¹ |
| 36 | +
|
| 37 | +For the bond at direction `I[1]` (which can be `NORTH` or `EAST`) |
| 38 | +from site `I[2], I[3]`, we identify the following gauge matrices, |
| 39 | +along the canonical direction of the PEPS arrows (`SOUTH ← NORTH` or `WEST ← EAST`): |
| 40 | +
|
| 41 | +```math |
| 42 | + I = √M₁₂⁻¹ √M₁₂ √M₂₁ √M₂₁⁻¹ |
| 43 | + = √M₁₂⁻¹ (U Λ Vᴴ) √M₂₁⁻¹ |
| 44 | + = (√M₁₂⁻¹ U √Λ) (√Λ Vᴴ √M₂₁⁻¹) |
| 45 | + = X X⁻¹ |
| 46 | +``` |
| 47 | +
|
| 48 | +Which are then used to update the gauge of `psi`. Thus, by convention `X` is attached to the `SOUTH`/`WEST` directions |
| 49 | +and `X⁻¹` is attached to the `NORTH`/`EAST` directions. |
| 50 | +""" |
| 51 | +function _bp_gauge_fix!(I::CartesianIndex{3}, psi::InfinitePEPS, env::BPEnv) |
| 52 | + dir, row, col = Tuple(I) |
| 53 | + @assert dir == NORTH || dir == EAST |
| 54 | + |
| 55 | + sqrtM12, isqrtM12, sqrtM21, isqrtM21 = _sqrt_bp_messages(I, env) |
| 56 | + U, Λ, Vᴴ = svd_compact!(sqrtM12 * sqrtM21) |
| 57 | + sqrtΛ = sdiag_pow(Λ, 1 / 2) |
| 58 | + X = isqrtM12 * U * sqrtΛ |
| 59 | + invX = sqrtΛ * Vᴴ * isqrtM21 |
| 60 | + if isdual(space(sqrtM12, 1)) |
| 61 | + X, invX = twist(flip(X, 2), 1), flip(invX, 1) |
| 62 | + end |
| 63 | + if dir == NORTH |
| 64 | + psi[row, col] = absorb_north_message(psi[row, col], X) |
| 65 | + psi[_prev(row, end), col] = absorb_south_message(psi[_prev(row, end), col], invX) |
| 66 | + elseif dir == EAST |
| 67 | + psi[row, col] = absorb_east_message(psi[row, col], X) |
| 68 | + psi[row, _next(col, end)] = absorb_west_message(psi[row, _next(col, end)], invX) |
| 69 | + end |
| 70 | + return psi, X, invX |
| 71 | +end |
| 72 | + |
| 73 | +""" |
| 74 | + SUWeight(env::BPEnv) |
| 75 | +
|
| 76 | +Construct `SUWeight` from belief propagation fixed point environment `env`. |
| 77 | +""" |
| 78 | +function SUWeight(env::BPEnv) |
| 79 | + wts = map(Iterators.product(1:2, axes(env, 2), axes(env, 3))) do (dir′, row, col) |
| 80 | + I = CartesianIndex(mod1(dir′ + 1, 2), row, col) |
| 81 | + sqrtM12, _, sqrtM21, _ = _sqrt_bp_messages(I, env) |
| 82 | + Λ = svd_vals!(sqrtM12 * sqrtM21) |
| 83 | + return isdual(space(sqrtM12, 1)) ? _fliptwist_s(Λ) : Λ |
| 84 | + end |
| 85 | + return SUWeight(wts) |
| 86 | +end |
| 87 | + |
| 88 | +""" |
| 89 | + BPEnv(wts::SUWeight) |
| 90 | +
|
| 91 | +Convert fixed point weights `wts` of trivial simple update |
| 92 | +to a belief propagation environment. |
| 93 | +""" |
| 94 | +function BPEnv(wts::SUWeight) |
| 95 | + messages = map(Iterators.product(1:4, axes(wts, 2), axes(wts, 3))) do (d, r, c) |
| 96 | + wt = if d == NORTH |
| 97 | + twist(wts[2, _next(r, end), c], 1) |
| 98 | + elseif d == EAST |
| 99 | + twist(wts[1, r, _prev(c, end)], 1) |
| 100 | + elseif d == SOUTH |
| 101 | + copy(wts[2, r, c]) |
| 102 | + else # WEST |
| 103 | + copy(wts[1, r, c]) |
| 104 | + end |
| 105 | + return TensorMap(wt) |
| 106 | + end |
| 107 | + return BPEnv(messages) |
| 108 | +end |
| 109 | + |
| 110 | +function sqrt_invsqrt(A::PEPSMessage) |
| 111 | + if isposdef(A) |
| 112 | + D, V = eigh_full(A) |
| 113 | + sqrtA = V * sdiag_pow(D, 1 / 2) * V' |
| 114 | + isqrtA = V * sdiag_pow(D, -1 / 2) * V' |
| 115 | + else |
| 116 | + D, V = eig_full(A) |
| 117 | + V⁻¹ = inv(V) |
| 118 | + sqrtA = V * sdiag_pow(D, 1 / 2) * V⁻¹ |
| 119 | + isqrtA = V * sdiag_pow(D, -1 / 2) * V⁻¹ |
| 120 | + if scalartype(A) <: Real |
| 121 | + # TODO: is this valid? |
| 122 | + sqrtA = real(sqrtA) |
| 123 | + isqrtA = real(isqrtA) |
| 124 | + end |
| 125 | + end |
| 126 | + return sqrtA, isqrtA |
| 127 | +end |
0 commit comments