From 831c4ff52c15339d39d1fb9f58f8b7fc59c797ca Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Thu, 5 Sep 2024 11:53:27 +0200 Subject: [PATCH 01/24] First pitiful attempt at affine transforms. --- src/Quantics.jl | 1 + src/affine.jl | 150 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 151 insertions(+) create mode 100644 src/affine.jl diff --git a/src/Quantics.jl b/src/Quantics.jl index 40f30ed..e21f7a4 100644 --- a/src/Quantics.jl +++ b/src/Quantics.jl @@ -31,5 +31,6 @@ include("mps.jl") include("fouriertransform.jl") include("imaginarytime.jl") include("transformer.jl") +include("affine.jl") end diff --git a/src/affine.jl b/src/affine.jl new file mode 100644 index 0000000..e643e4f --- /dev/null +++ b/src/affine.jl @@ -0,0 +1,150 @@ +using StaticArrays +using ITensors + +function affine_transform_mpo( + outsite::AbstractMatrix{<:Index}, insite::AbstractMatrix{<:Index}, + A::AbstractMatrix{<:Integer}, b::AbstractVector{<:Integer} + ) + R = size(outsite, 1) + M, N = size(A) + size(insite) == (R, N) || + throw(ArgumentError("insite is not correctly dimensioned")) + size(b) == (M,) || + throw(ArgumentError("vector is not correctly dimensioned")) + + # get the tensors so that we know how large the links must be + tensors = affine_transform_matrices( + R, SMatrix{M, N, Int}(A), SVector{M, Int}(b)) + + # Create the links + link = [Index(size(tensors[r], 1), tags="link $r") for r in 1:R-1] + + # Fill the MPO, taking care to not include auxiliary links at the edges + mpo = MPO(R) + spin_dims = ntuple(_ -> 2, M + N) + mpo[1] = ITensor(reshape(tensors[1], size(tensors[1])[1], spin_dims...), + (link[1], outsite[1,:]..., insite[1,:]...)) + for r in 2:R-1 + newshape = size(tensors[r])[1:2]..., spin_dims... + mpo[r] = ITensor(reshape(tensors[r], newshape), + (link[r], link[r-1], outsite[r,:]..., insite[r,:]...)) + end + println(size(tensors[R])) + mpo[R] = ITensor(reshape(tensors[R], size(tensors[R])[2], spin_dims...), + (link[R-1], outsite[R,:]..., insite[R,:]...)) + return mpo +end + +function affine_transform_matrices( + R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int} + ) where {M, N} + # Checks + 0 <= R <= 62 || + throw(ArgumentError("invalid value of the length R")) + + # Matrix which maps out which bits to add together. We use 2's complement, + # -x == ~x + 1, which means we use two masks: A1 operating on set bits, and + # A0 operating on unset bits. We group the +1 together with the shift. + A1 = @. clamp(A, 0, typemax(Int)) + A0 = @. clamp(-A, 0, typemax(Int)) + b = b .+ vec(sum(Bool.(A0), dims=2)) + println(A1, A0, b) + + # Amax is the maximum value that can be reached by multiplying it with set + # or unset bits. + Amax = vec(sum(A0 + A1, dims=2)) + + # The output tensors are a collection of matrices, but their first two + # dimensions (links) vary + tensors = Array{Bool, 4}[] + sizehint!(tensors, R) + + # We have to carry all but the least siginificant bit to the adjacent + # tensor. α is the index of the "outgoing" carrys, which at the beginning + # is simply all zeros. + maxcarry_α = zeros(typeof(b)) + ncarry_α = prod(maxcarry_α .+ 1) + base_α = _indexbase_from_size(maxcarry_α, 1) + for r in 1:R + # Copy the outgoing carry α of the previous tensor to the "incoming" + # carry β of the current tensor. + maxcarry_β = maxcarry_α + ncarry_β = ncarry_α + + # Figure out the current bit to add from the shift term and shift + # it out from the array + bcurr = @. b & 1 + b = @. b >> 1 + + # Determine the maximum outgoing carry. It is determined by the maximum + # previous carry plus the maximum value of the mapped indices plus the + # current shift, all divided by two. In the case of the last tensor, we + # discard all carrys. + if r == R + maxcarry_α = zeros(typeof(maxcarry_α)) + ncarry_α = 1 + base_α = _indexbase_from_size(maxcarry_α, 0) + else + maxcarry_α = @. (maxcarry_β + Amax + bcurr) >> 1 + ncarry_α = prod(maxcarry_α .+ 1) + base_α = _indexbase_from_size(maxcarry_α, 1) + end + + values = zeros(Bool, ncarry_α, ncarry_β, 1 << M, 1 << N) + println("\n r=$r $(bcurr)") + # Fill values array. The idea is the following: we iterate over all + # possible values of the carry from the previous tensor (c_β). Each of + # those corresponds to a bond index β. + allcarry_β = map(i -> 0:i, maxcarry_β) + for (β, _c_β) in enumerate(Iterators.product(allcarry_β...)) + c_β = SVector{M}(_c_β) + + # Now we iterate over all possible indices of the input legs, which + # are all combination of values of N input bits + allspin_j = ntuple(_ -> 0:1, N) + for (j, _σ_j) in enumerate(Iterators.product(allspin_j...)) + σ_j = SVector{N, Bool}(_σ_j) + + # We apply the transformation y = Ax + b to the current bits + # σ_j and adding the incoming carry c_β. The least significant + # bits are then the output σ_i while the higher-order bits + # form the outgoing carry c_α + ifull = A1 * σ_j + A0 * (.~σ_j) + bcurr + c_β + c_α = @. ifull >> 1 + σ_i = SVector{M, Bool}(@. ifull & 1) + + # Map the outgoing carry to an index and store + α = mapreduce(*, +, c_α, base_α, init=1) + i = _spins_to_number(σ_i) + 1 + println("$(c_α) * 2 + $(σ_i) <- $(c_β) * 2 + $(σ_j)") + values[α, β, i, j] = true + end + end + push!(tensors, values) + end + return tensors +end + +""" + get_int_inverse(A::AbstractMatrix{<:Integer}) + +Return inverse matrix to integer A, ensuring that it is indeed integer. +""" +function get_int_invserse(A::AbstractMatrix{<:Integer}) + Ainv = inv(A) + Ainv_int = map(eltype(A) ∘ round, Ainv) + isapprox(Ainv, Ainv_int, atol=size(A,1)*eps()) || + throw(DomainError(Ainv, "inverse is not integer")) + return typeof(A)(Ainv) +end + +_indexbase_from_size(v::SVector{M,T}, init::T) where {M,T} = + SVector{M,T}(_indexbase_from_size(Tuple(v), init)) +_indexbase_from_size(v::Tuple{}, init::T) where {T} = () +_indexbase_from_size(v::Tuple{T, Vararg{T}}, init::T) where {T} = + init, _indexbase_from_size(v[2:end], init * v[1])... + +_spins_to_number(v::SVector{<:Any, Bool}) = _spins_to_number(Tuple(v)) +_spins_to_number(v::Tuple{Bool, Vararg{Bool}}) = + _spins_to_number(v[2:end]) << 1 | v[1] +_spins_to_number(v::Tuple{}) = 0 From 1283f78c43e9882f5cf85b1b2edae4793bec665d Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Thu, 5 Sep 2024 12:46:29 +0200 Subject: [PATCH 02/24] Remove prints --- src/affine.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/affine.jl b/src/affine.jl index e643e4f..901cd95 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -48,7 +48,7 @@ function affine_transform_matrices( A1 = @. clamp(A, 0, typemax(Int)) A0 = @. clamp(-A, 0, typemax(Int)) b = b .+ vec(sum(Bool.(A0), dims=2)) - println(A1, A0, b) + #println(A1, A0, b) # Amax is the maximum value that can be reached by multiplying it with set # or unset bits. @@ -91,7 +91,7 @@ function affine_transform_matrices( end values = zeros(Bool, ncarry_α, ncarry_β, 1 << M, 1 << N) - println("\n r=$r $(bcurr)") + #println("\n r=$r $(bcurr)") # Fill values array. The idea is the following: we iterate over all # possible values of the carry from the previous tensor (c_β). Each of # those corresponds to a bond index β. @@ -116,7 +116,7 @@ function affine_transform_matrices( # Map the outgoing carry to an index and store α = mapreduce(*, +, c_α, base_α, init=1) i = _spins_to_number(σ_i) + 1 - println("$(c_α) * 2 + $(σ_i) <- $(c_β) * 2 + $(σ_j)") + #println("$(c_α) * 2 + $(σ_i) <- $(c_β) * 2 + $(σ_j)") values[α, β, i, j] = true end end From e62e99710ac07dac77bb63a04847d0fe42f5366f Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Mon, 9 Sep 2024 19:59:53 +0200 Subject: [PATCH 03/24] Now affine transforms seem to work... --- notebook/quantics.py | 98 +++++++++++++++++++++++++++++++ src/affine.jl | 136 +++++++++++++++---------------------------- 2 files changed, 146 insertions(+), 88 deletions(-) create mode 100644 notebook/quantics.py diff --git a/notebook/quantics.py b/notebook/quantics.py new file mode 100644 index 0000000..e3de7cb --- /dev/null +++ b/notebook/quantics.py @@ -0,0 +1,98 @@ +import numpy as np +import scipy.interpolate + + +def decompose_int(A, n, tol=None): + A = np.asarray(A) + Ushape = A.shape[:n] + VTshape = A.shape[n:] + + Aflat = A.reshape(np.prod(Ushape), np.prod(VTshape)) + Uflat, VTflat = decompose_int_matrix(Aflat, tol) + U = Uflat.reshape(*Ushape, -1) + VT = VTflat.reshape(-1, *VTshape) + return U, VT + + +def decompose_int_matrix(A, tol=None): + A = np.asarray(A) + if np.issubdtype(A.dtype, np.integer): + A = A.astype(np.float32) + if tol is None: + tol = max(2, np.sqrt(A.shape[0])) * np.finfo(A.dtype).eps + + # Truncated SVD + U, s, VT = np.linalg.svd(A, full_matrices=False) + s_mask = (s > tol * s[0]) + U = U[:, s_mask] + s = s[s_mask] + VT = VT[s_mask, :] + + # Truncate small elements + U_mask = np.abs(U) > tol + VT_mask = np.abs(VT) > tol + U[~U_mask] = 0 + VT[~VT_mask] = 0 + + # Undo scaling + U_count = U_mask.sum(0) + VT_count = VT_mask.sum(1) + U *= np.sqrt(U_count) + VT *= np.sqrt(VT_count)[:, None] + s /= np.sqrt(U_count * VT_count) + np.testing.assert_allclose(s, 1, atol=tol, rtol=0) + + # Round stuff + U_int = np.round(U).astype(int) + VT_int = np.round(VT).astype(int) + np.testing.assert_allclose(U_int, U, atol=tol, rtol=0) + np.testing.assert_allclose(VT_int, VT, atol=tol, rtol=0) + + # Regauge values + U_sgn = np.sign(U_int.sum(0)) + np.testing.assert_equal(np.abs(U_sgn), 1) + U_int *= U_sgn + VT_int *= U_sgn[:, None] + + # We want a bit matrix + np.testing.assert_equal(U_int, U_int.astype(bool)) + np.testing.assert_equal(VT_int, VT_int.astype(bool)) + + # Return + return U_int, VT_int + + +def interleave_bits(A, K): + A = np.asarray(A) + R = A.ndim // K + if A.shape != (2,) * (K*R): + raise ValueError("not of proper quantics shape") + + order = np.hstack([np.arange(i, K*R, R) for i in range(R)]) + return A.transpose(*order) + + +def quantize(A, interleave=False): + A = np.asarray(A) + bits = np.log2(A.shape) + if not (bits == bits.astype(int)).all(): + raise ValueError("invalid shape: " + str(A.shape)) + bits = bits.astype(int) + A = A.reshape((2,) * bits.sum()) + + if interleave: + if not (bits == bits[0]).all(): + raise ValueError("unable to interleave") + A = interleave_bits(A, len(bits)) + return A + + +def print_nonzeros(A, headers=None): + A = np.asarray(A) + N = A.ndim + if headers: + headers = tuple(headers) + print(("%2s " * N) % headers) + fmt = " ".join(("%2i",) * N) + for x in zip(*A.nonzero()): + print(fmt % x) diff --git a/src/affine.jl b/src/affine.jl index 901cd95..f352ff7 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -13,7 +13,7 @@ function affine_transform_mpo( throw(ArgumentError("vector is not correctly dimensioned")) # get the tensors so that we know how large the links must be - tensors = affine_transform_matrices( + tensors = affine_transform_tensors( R, SMatrix{M, N, Int}(A), SVector{M, Int}(b)) # Create the links @@ -35,115 +35,75 @@ function affine_transform_mpo( return mpo end -function affine_transform_matrices( +function affine_transform_tensors( R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int} ) where {M, N} # Checks - 0 <= R <= 62 || - throw(ArgumentError("invalid value of the length R")) - - # Matrix which maps out which bits to add together. We use 2's complement, - # -x == ~x + 1, which means we use two masks: A1 operating on set bits, and - # A0 operating on unset bits. We group the +1 together with the shift. - A1 = @. clamp(A, 0, typemax(Int)) - A0 = @. clamp(-A, 0, typemax(Int)) - b = b .+ vec(sum(Bool.(A0), dims=2)) - #println(A1, A0, b) - - # Amax is the maximum value that can be reached by multiplying it with set - # or unset bits. - Amax = vec(sum(A0 + A1, dims=2)) + 0 <= R || + throw(DomainError(R, "invalid value of the length R")) # The output tensors are a collection of matrices, but their first two # dimensions (links) vary - tensors = Array{Bool, 4}[] - sizehint!(tensors, R) + tensors = Vector{Array{Bool, 4}}(undef, R) - # We have to carry all but the least siginificant bit to the adjacent - # tensor. α is the index of the "outgoing" carrys, which at the beginning - # is simply all zeros. - maxcarry_α = zeros(typeof(b)) - ncarry_α = prod(maxcarry_α .+ 1) - base_α = _indexbase_from_size(maxcarry_α, 1) + # The initial carry is zero + carry = [zero(SVector{M, Int})] for r in 1:R - # Copy the outgoing carry α of the previous tensor to the "incoming" - # carry β of the current tensor. - maxcarry_β = maxcarry_α - ncarry_β = ncarry_α - # Figure out the current bit to add from the shift term and shift # it out from the array - bcurr = @. b & 1 - b = @. b >> 1 + bcurr = @. Bool(b & 1) - # Determine the maximum outgoing carry. It is determined by the maximum - # previous carry plus the maximum value of the mapped indices plus the - # current shift, all divided by two. In the case of the last tensor, we - # discard all carrys. + # Get tensor. For the last tensor, we ignore the carry and just + # or together all possible combinations. + new_carry, new_tensor = affine_transform_tensor(A, bcurr, carry) if r == R - maxcarry_α = zeros(typeof(maxcarry_α)) - ncarry_α = 1 - base_α = _indexbase_from_size(maxcarry_α, 0) + all_values = reduce(.|, eachslice(new_tensor, dims=1)) + tensors[r] = reshape(all_values, 1, size(all_values)...) else - maxcarry_α = @. (maxcarry_β + Amax + bcurr) >> 1 - ncarry_α = prod(maxcarry_α .+ 1) - base_α = _indexbase_from_size(maxcarry_α, 1) + tensors[r] = new_tensor end - values = zeros(Bool, ncarry_α, ncarry_β, 1 << M, 1 << N) - #println("\n r=$r $(bcurr)") - # Fill values array. The idea is the following: we iterate over all - # possible values of the carry from the previous tensor (c_β). Each of - # those corresponds to a bond index β. - allcarry_β = map(i -> 0:i, maxcarry_β) - for (β, _c_β) in enumerate(Iterators.product(allcarry_β...)) - c_β = SVector{M}(_c_β) - - # Now we iterate over all possible indices of the input legs, which - # are all combination of values of N input bits - allspin_j = ntuple(_ -> 0:1, N) - for (j, _σ_j) in enumerate(Iterators.product(allspin_j...)) - σ_j = SVector{N, Bool}(_σ_j) - - # We apply the transformation y = Ax + b to the current bits - # σ_j and adding the incoming carry c_β. The least significant - # bits are then the output σ_i while the higher-order bits - # form the outgoing carry c_α - ifull = A1 * σ_j + A0 * (.~σ_j) + bcurr + c_β - c_α = @. ifull >> 1 - σ_i = SVector{M, Bool}(@. ifull & 1) - - # Map the outgoing carry to an index and store - α = mapreduce(*, +, c_α, base_α, init=1) - i = _spins_to_number(σ_i) + 1 - #println("$(c_α) * 2 + $(σ_i) <- $(c_β) * 2 + $(σ_j)") - values[α, β, i, j] = true - end - end - push!(tensors, values) + # Set carry to the next value + carry = new_carry + b = @. b >> 1 end return tensors end -""" - get_int_inverse(A::AbstractMatrix{<:Integer}) +function affine_transform_tensor( + A::SMatrix{M, N, Int}, b::SVector{M, Bool}, + carry::AbstractVector{SVector{M, Int}}) where {M, N} + # The basic idea here is the following: we compute r = A*x + b + c for all + # "incoming" carrys d and all possible bit vectors, x ∈ {0,1}^N. Then we + # decompose r = 2*c + y, where y is again a bit vector, y ∈ {0,1}^M, and + # c is the "outgoing" carry, which may be negative. We then store this + # as something like out[d][c, x, y] = true. + out = Dict{SVector{M, Int}, Array{Bool, 3}}() + for (c_index, c) in enumerate(carry) + for (x_index, x) in enumerate(Iterators.product(ntuple(_ -> 0:1, N)...)) + r = A * SVector{N, Bool}(x) + b + SVector{N, Int}(c) + y = Bool.(r .& 1) + d = r .>> 1 + y_index = _spins_to_number(y) + 1 + + d_mat = get!(out, d) do + return zeros(Bool, length(carry), 1 << M, 1 << N) + end + d_mat[c_index, x_index, y_index] = true + end + end -Return inverse matrix to integer A, ensuring that it is indeed integer. -""" -function get_int_invserse(A::AbstractMatrix{<:Integer}) - Ainv = inv(A) - Ainv_int = map(eltype(A) ∘ round, Ainv) - isapprox(Ainv, Ainv_int, atol=size(A,1)*eps()) || - throw(DomainError(Ainv, "inverse is not integer")) - return typeof(A)(Ainv) + # We translate the dictionary into a vector of carrys (which we can then + # pass into the next iteration) and a 4-way tensor of output values. + carry_out = Vector{SVector{M, Int}}(undef, length(out)) + value_out = Array{Bool, 4}(undef, length(out), length(carry), 1< Date: Mon, 9 Sep 2024 22:47:11 +0200 Subject: [PATCH 04/24] Affine transforms are now correct (mostly) --- src/affine.jl | 60 ++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 45 insertions(+), 15 deletions(-) diff --git a/src/affine.jl b/src/affine.jl index f352ff7..08c6965 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -35,6 +35,20 @@ function affine_transform_mpo( return mpo end +""" + affine_transform_tensors(R, A, b) + +Compute vector of core tensors (constituent 4-way tensors) for a matrix product +operator corresponding to the affine transformation `y = A*x + b`. +""" +function affine_transform_tensors( + R::Integer, A::AbstractMatrix{<:Integer}, + b::AbstractVector{<:Integer}) + M, N = size(A) + return affine_transform_tensors( + Int(R), SMatrix{M, N, Int}(A), SVector{M, Int}(b)) +end + function affine_transform_tensors( R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int} ) where {M, N} @@ -53,15 +67,11 @@ function affine_transform_tensors( # it out from the array bcurr = @. Bool(b & 1) - # Get tensor. For the last tensor, we ignore the carry and just - # or together all possible combinations. - new_carry, new_tensor = affine_transform_tensor(A, bcurr, carry) - if r == R - all_values = reduce(.|, eachslice(new_tensor, dims=1)) - tensors[r] = reshape(all_values, 1, size(all_values)...) - else - tensors[r] = new_tensor - end + # Get tensor. For the last tensor, we discard the carry. + out_tf = (r == R) ? zero : identity + new_carry, data = affine_transform_core(A, bcurr, carry, + transform_out_carry=out_tf) + tensors[r] = data # Set carry to the next value carry = new_carry @@ -70,26 +80,46 @@ function affine_transform_tensors( return tensors end -function affine_transform_tensor( +""" + core, out_carry = affine_transform_core(A, b, in_carry; + transform_out_carry=identity) + +Construct core tensor `core` for an affine transform. The core tensor for an +affine transformation is given by: + + core[d, c, iy, ix] = + 2 * out_carry[d] + y[iy] == A * x[ix] + b + in_carry[c] + +where `A`, a matrix of integers, and `b`, a vector of bits, which define the +affine transform. `c` and `d` are indices into a set of integer vectors +`in_carry` and `out_carry`, respectively, which encode the incoming and outgoing +carry from the other core tensors. `x[ix] ∈ {0,1}^N` and `y[iy] ∈ {0,1}^M` +are the binary input and output vectors, respectively, of the affine transform. +They are indexed in a "little-endian" fashion. +""" +function affine_transform_core( A::SMatrix{M, N, Int}, b::SVector{M, Bool}, - carry::AbstractVector{SVector{M, Int}}) where {M, N} + carry::AbstractVector{SVector{M, Int}}; + transform_out_carry::Function=identity + ) where {M, N} # The basic idea here is the following: we compute r = A*x + b + c for all # "incoming" carrys d and all possible bit vectors, x ∈ {0,1}^N. Then we # decompose r = 2*c + y, where y is again a bit vector, y ∈ {0,1}^M, and # c is the "outgoing" carry, which may be negative. We then store this # as something like out[d][c, x, y] = true. out = Dict{SVector{M, Int}, Array{Bool, 3}}() + sizehint!(out, length(carry)) for (c_index, c) in enumerate(carry) for (x_index, x) in enumerate(Iterators.product(ntuple(_ -> 0:1, N)...)) r = A * SVector{N, Bool}(x) + b + SVector{N, Int}(c) - y = Bool.(r .& 1) - d = r .>> 1 + y = @. Bool(r & 1) + d::SVector{M, Int} = transform_out_carry(r .>> 1) y_index = _spins_to_number(y) + 1 d_mat = get!(out, d) do return zeros(Bool, length(carry), 1 << M, 1 << N) end - d_mat[c_index, x_index, y_index] = true + @inbounds d_mat[c_index, x_index, y_index] = true end end @@ -99,7 +129,7 @@ function affine_transform_tensor( value_out = Array{Bool, 4}(undef, length(out), length(carry), 1< Date: Tue, 10 Sep 2024 14:10:16 +0200 Subject: [PATCH 05/24] Affine transforms are now ready to be tested. --- Project.toml | 1 + notebook/quantics.py | 98 -------------------------------------------- src/affine.jl | 96 ++++++++++++++++++++++++++++++++++++++----- 3 files changed, 86 insertions(+), 109 deletions(-) delete mode 100644 notebook/quantics.py diff --git a/Project.toml b/Project.toml index 0c7ff37..34c9aac 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949" FastMPOContractions = "f6e391d2-8ffa-4d7a-98cd-7e70024481ca" ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseIR = "4fe2279e-80f0-4adb-8463-ee114ff56b7d" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/notebook/quantics.py b/notebook/quantics.py deleted file mode 100644 index e3de7cb..0000000 --- a/notebook/quantics.py +++ /dev/null @@ -1,98 +0,0 @@ -import numpy as np -import scipy.interpolate - - -def decompose_int(A, n, tol=None): - A = np.asarray(A) - Ushape = A.shape[:n] - VTshape = A.shape[n:] - - Aflat = A.reshape(np.prod(Ushape), np.prod(VTshape)) - Uflat, VTflat = decompose_int_matrix(Aflat, tol) - U = Uflat.reshape(*Ushape, -1) - VT = VTflat.reshape(-1, *VTshape) - return U, VT - - -def decompose_int_matrix(A, tol=None): - A = np.asarray(A) - if np.issubdtype(A.dtype, np.integer): - A = A.astype(np.float32) - if tol is None: - tol = max(2, np.sqrt(A.shape[0])) * np.finfo(A.dtype).eps - - # Truncated SVD - U, s, VT = np.linalg.svd(A, full_matrices=False) - s_mask = (s > tol * s[0]) - U = U[:, s_mask] - s = s[s_mask] - VT = VT[s_mask, :] - - # Truncate small elements - U_mask = np.abs(U) > tol - VT_mask = np.abs(VT) > tol - U[~U_mask] = 0 - VT[~VT_mask] = 0 - - # Undo scaling - U_count = U_mask.sum(0) - VT_count = VT_mask.sum(1) - U *= np.sqrt(U_count) - VT *= np.sqrt(VT_count)[:, None] - s /= np.sqrt(U_count * VT_count) - np.testing.assert_allclose(s, 1, atol=tol, rtol=0) - - # Round stuff - U_int = np.round(U).astype(int) - VT_int = np.round(VT).astype(int) - np.testing.assert_allclose(U_int, U, atol=tol, rtol=0) - np.testing.assert_allclose(VT_int, VT, atol=tol, rtol=0) - - # Regauge values - U_sgn = np.sign(U_int.sum(0)) - np.testing.assert_equal(np.abs(U_sgn), 1) - U_int *= U_sgn - VT_int *= U_sgn[:, None] - - # We want a bit matrix - np.testing.assert_equal(U_int, U_int.astype(bool)) - np.testing.assert_equal(VT_int, VT_int.astype(bool)) - - # Return - return U_int, VT_int - - -def interleave_bits(A, K): - A = np.asarray(A) - R = A.ndim // K - if A.shape != (2,) * (K*R): - raise ValueError("not of proper quantics shape") - - order = np.hstack([np.arange(i, K*R, R) for i in range(R)]) - return A.transpose(*order) - - -def quantize(A, interleave=False): - A = np.asarray(A) - bits = np.log2(A.shape) - if not (bits == bits.astype(int)).all(): - raise ValueError("invalid shape: " + str(A.shape)) - bits = bits.astype(int) - A = A.reshape((2,) * bits.sum()) - - if interleave: - if not (bits == bits[0]).all(): - raise ValueError("unable to interleave") - A = interleave_bits(A, len(bits)) - return A - - -def print_nonzeros(A, headers=None): - A = np.asarray(A) - N = A.ndim - if headers: - headers = tuple(headers) - print(("%2s " * N) % headers) - fmt = " ".join(("%2i",) * N) - for x in zip(*A.nonzero()): - print(fmt % x) diff --git a/src/affine.jl b/src/affine.jl index 08c6965..05c3a78 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -1,6 +1,13 @@ using StaticArrays +using SparseArrays using ITensors +""" + affine_transform_mpo(outsite, insite, A, b) + +Constructs an ITensor matrix product state for an affine transform in the +quantics representation. +""" function affine_transform_mpo( outsite::AbstractMatrix{<:Index}, insite::AbstractMatrix{<:Index}, A::AbstractMatrix{<:Integer}, b::AbstractVector{<:Integer} @@ -29,7 +36,6 @@ function affine_transform_mpo( mpo[r] = ITensor(reshape(tensors[r], newshape), (link[r], link[r-1], outsite[r,:]..., insite[r,:]...)) end - println(size(tensors[R])) mpo[R] = ITensor(reshape(tensors[R], size(tensors[R])[2], spin_dims...), (link[R-1], outsite[R,:]..., insite[R,:]...)) return mpo @@ -44,9 +50,7 @@ operator corresponding to the affine transformation `y = A*x + b`. function affine_transform_tensors( R::Integer, A::AbstractMatrix{<:Integer}, b::AbstractVector{<:Integer}) - M, N = size(A) - return affine_transform_tensors( - Int(R), SMatrix{M, N, Int}(A), SVector{M, Int}(b)) + return affine_transform_tensors(Int(R), _affine_static_args(A, b)...) end function affine_transform_tensors( @@ -111,15 +115,15 @@ function affine_transform_core( sizehint!(out, length(carry)) for (c_index, c) in enumerate(carry) for (x_index, x) in enumerate(Iterators.product(ntuple(_ -> 0:1, N)...)) - r = A * SVector{N, Bool}(x) + b + SVector{N, Int}(c) + r = A * SVector{N, Bool}(x) + b + SVector{M, Int}(c) y = @. Bool(r & 1) d::SVector{M, Int} = transform_out_carry(r .>> 1) - y_index = _spins_to_number(y) + 1 + y_index = digits_to_number(y) + 1 d_mat = get!(out, d) do return zeros(Bool, length(carry), 1 << M, 1 << N) end - @inbounds d_mat[c_index, x_index, y_index] = true + @inbounds d_mat[c_index, y_index, x_index] = true end end @@ -134,7 +138,77 @@ function affine_transform_core( return carry_out, value_out end -_spins_to_number(v::SVector{<:Any, Bool}) = _spins_to_number(Tuple(v)) -_spins_to_number(v::Tuple{Bool, Vararg{Bool}}) = - _spins_to_number(v[2:end]) << 1 | v[1] -_spins_to_number(v::Tuple{}) = 0 +""" + affine_transform_matrix(R, A, b) + +Compute full transformation matrix for `y = A*x + b`. +""" +function affine_transform_matrix( + R::Integer, A::AbstractMatrix{<:Integer}, + b::AbstractVector{<:Integer}) + return affine_transform_matrix(Int(R), _affine_static_args(A, b)...) +end + +function affine_transform_matrix( + R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int}) where {M, N} + # Checks + 0 <= R || + throw(DomainError(R, "invalid value of the length R")) + + mask = (1 << R) - 1 + y_index = Vector{Int}(undef, 1 << (R * N)) + + for (ix, x) in enumerate(Iterators.product(ntuple(_ -> 0:mask, N)...)) + y = (A * SVector{N, Int}(x) + b) .& mask + iy = digits_to_number(y, R) + 1 + @inbounds y_index[ix] = iy + end + x_index = collect(1:1 << (R * N)) + values = ones(Bool, 1 << (R * N)) + return sparse(y_index, x_index, values, 1 << (R*M), 1 << (R*N)) +end + +function affine_mpo_to_matrix( + outsite::AbstractMatrix{<:Index}, insite::AbstractMatrix{<:Index}, + mpo::ITensors.MPO) + prev_warn_order = ITensors.disable_warn_order() + try + mpo_contr = reduce(*, mpo) + tensor = Array(mpo_contr, vec(outsite)..., vec(insite)...) + matrix = reshape(tensor, 1 << length(outsite), 1 << length(insite)) + return matrix + finally + ITensors.set_warn_order!(prev_warn_order) + end +end + +function _affine_static_args( + A::AbstractMatrix{<:Integer}, b::AbstractVector{<:Integer}) + M, N = size(A) + size(b, 1) == N || + throw(ArgumentError("A and b have incompatible size")) + return SMatrix{M, N, Int}(A), SVector{M, Int}(b) +end + +""" + digits_to_number(v::AbstractVector{Bool}) + digits_to_number(v::AbstractVector{<:Integer}, bits::Integer) + +Converts a vector of digits, starting with the least significant digit, to +a number. If the digits are boolean, then they are interpreted in base-2, +otherwise, in base `2^bits`. +""" +digits_to_number(v::AbstractVector{Bool}) = _digits_to_number(Tuple(v)) +digits_to_number(v::AbstractVector{<:Integer}, bits::Integer) = + _digits_to_number(Int.(Tuple(v)), Int(bits)) + +@inline _digits_to_number(v::Tuple{}) = 0 +@inline _digits_to_number(v::Tuple{Bool, Vararg{Bool}}) = + _digits_to_number(v[2:end]) << 1 | v[1] +@inline _digits_to_number(v::Tuple{}, bits::Int) = 0 +@inline function _digits_to_number(v::Tuple{Int, Vararg{Int}}, bits::Int) + mask = (~0) << bits + iszero(v[1] & mask) || + throw(DomainError(v[1], "invalid digit in base $(1 << bits)")) + return _digits_to_number(v[2:end], bits) << bits | v[1] +end From 792fed8353dadc18d84f996c0b4666bc14b955bf Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Tue, 10 Sep 2024 14:54:47 +0200 Subject: [PATCH 06/24] Add test. --- .gitignore | 9 +++++++ src/Quantics.jl | 1 + src/affine.jl | 8 ++---- test/affine_tests.jl | 63 ++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 75 insertions(+), 6 deletions(-) create mode 100644 test/affine_tests.jl diff --git a/.gitignore b/.gitignore index 82e9f45..6c50e99 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,11 @@ +.* +*~ +\#*\# + +!.gitignore +!/.github/ +!/.gitlab* + *.jl.*.cov *.jl.cov *.jl.mem @@ -6,3 +14,4 @@ /docs/Manifest.toml /docs/build/ docs/build +__pycache__/ diff --git a/src/Quantics.jl b/src/Quantics.jl index e21f7a4..0763d3f 100644 --- a/src/Quantics.jl +++ b/src/Quantics.jl @@ -15,6 +15,7 @@ import ITensors.NDTensors: Tensor, BlockSparseTensor, blockview import SparseIR: Fermionic, Bosonic, Statistics import LinearAlgebra: I using StaticArrays +import SparseArrays: sparse import FastMPOContractions diff --git a/src/affine.jl b/src/affine.jl index 05c3a78..5e73aed 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -1,7 +1,3 @@ -using StaticArrays -using SparseArrays -using ITensors - """ affine_transform_mpo(outsite, insite, A, b) @@ -178,14 +174,14 @@ function affine_mpo_to_matrix( matrix = reshape(tensor, 1 << length(outsite), 1 << length(insite)) return matrix finally - ITensors.set_warn_order!(prev_warn_order) + ITensors.set_warn_order(prev_warn_order) end end function _affine_static_args( A::AbstractMatrix{<:Integer}, b::AbstractVector{<:Integer}) M, N = size(A) - size(b, 1) == N || + size(b, 1) == M || throw(ArgumentError("A and b have incompatible size")) return SMatrix{M, N, Int}(A), SVector{M, Int}(b) end diff --git a/test/affine_tests.jl b/test/affine_tests.jl new file mode 100644 index 0000000..8c328fb --- /dev/null +++ b/test/affine_tests.jl @@ -0,0 +1,63 @@ +@testitem "affine" begin + using Test + using ITensors + using Quantics + using LinearAlgebra + + vars = ("x", "y", "z") + insite = [Index(2; tags="$v$l") for l in 1:10, v in vars] + outsite = [Index(2; tags="$v$l")' for l in 1:10, v in vars] + + @testset "full" begin + A = [1 0; 1 1] + b = [0; 0] + + T = Quantics.affine_transform_matrix(4, A, b) + @test T' * T == I + @test T[219, 59] == 1 + + b = [4; 1] + T = Quantics.affine_transform_matrix(4, A, b) + @test T * T' == I + end + + @testset "compare_simple" begin + A = [1 0; 1 1] + b = [0; 0] + R = 3 + + T = Quantics.affine_transform_matrix(R, A, b) + mpo = Quantics.affine_transform_mpo( + outsite[1:R, 1:2], insite[1:R, 1:2], A, b) + Trec = Quantics.affine_mpo_to_matrix( + outsite[1:R, 1:2], insite[1:R, 1:2], mpo) + @test T == Trec + end + + @testset "compare_hard" begin + A = [1 0 1; 1 2 -1; 0 1 1] + b = [11; 23; -15] + R = 4 + + T = Quantics.affine_transform_matrix(R, A, b) + mpo = Quantics.affine_transform_mpo( + outsite[1:R, 1:3], insite[1:R, 1:3], A, b) + Trec = Quantics.affine_mpo_to_matrix( + outsite[1:R, 1:3], insite[1:R, 1:3], mpo) + @test T == Trec + end + + @testset "compare_rect" begin + A = [1 0 1; 1 2 0] + b = [11; -3] + R = 4 + + T = Quantics.affine_transform_matrix(R, A, b) + mpo = Quantics.affine_transform_mpo( + outsite[1:R, 1:2], insite[1:R, 1:3], A, b) + Trec = Quantics.affine_mpo_to_matrix( + outsite[1:R, 1:2], insite[1:R, 1:3], mpo) + @test T == Trec + end + +end From d21b8330520f9c729fd78b1a11e775e0b76fb54a Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Tue, 10 Sep 2024 15:16:38 +0200 Subject: [PATCH 07/24] Conform to the Quantics convention (1 ... most significant bit) --- src/affine.jl | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/src/affine.jl b/src/affine.jl index 5e73aed..23a6673 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -20,19 +20,20 @@ function affine_transform_mpo( R, SMatrix{M, N, Int}(A), SVector{M, Int}(b)) # Create the links - link = [Index(size(tensors[r], 1), tags="link $r") for r in 1:R-1] + link = [Index(size(tensors[r], 2), tags="link $r") for r in 1:R-1] # Fill the MPO, taking care to not include auxiliary links at the edges mpo = MPO(R) spin_dims = ntuple(_ -> 2, M + N) - mpo[1] = ITensor(reshape(tensors[1], size(tensors[1])[1], spin_dims...), + println(size.(tensors)) + mpo[1] = ITensor(reshape(tensors[1], size(tensors[1], 2), spin_dims...), (link[1], outsite[1,:]..., insite[1,:]...)) for r in 2:R-1 newshape = size(tensors[r])[1:2]..., spin_dims... mpo[r] = ITensor(reshape(tensors[r], newshape), - (link[r], link[r-1], outsite[r,:]..., insite[r,:]...)) + (link[r-1], link[r], outsite[r,:]..., insite[r,:]...)) end - mpo[R] = ITensor(reshape(tensors[R], size(tensors[R])[2], spin_dims...), + mpo[R] = ITensor(reshape(tensors[R], size(tensors[R], 1), spin_dims...), (link[R-1], outsite[R,:]..., insite[R,:]...)) return mpo end @@ -62,13 +63,13 @@ function affine_transform_tensors( # The initial carry is zero carry = [zero(SVector{M, Int})] - for r in 1:R + for r in R:-1:1 # Figure out the current bit to add from the shift term and shift # it out from the array bcurr = @. Bool(b & 1) - # Get tensor. For the last tensor, we discard the carry. - out_tf = (r == R) ? zero : identity + # Get tensor. For the first tensor, we discard the carry. + out_tf = (r == 1) ? zero : identity new_carry, data = affine_transform_core(A, bcurr, carry, transform_out_carry=out_tf) tensors[r] = data @@ -170,8 +171,12 @@ function affine_mpo_to_matrix( prev_warn_order = ITensors.disable_warn_order() try mpo_contr = reduce(*, mpo) - tensor = Array(mpo_contr, vec(outsite)..., vec(insite)...) - matrix = reshape(tensor, 1 << length(outsite), 1 << length(insite)) + out_indices = vec(reverse(outsite, dims=1)) + in_indices = vec(reverse(insite, dims=1)) + println(out_indices) + println(in_indices) + tensor = Array(mpo_contr, out_indices..., in_indices...) + matrix = reshape(tensor, 1 << length(out_indices), 1 << length(in_indices)) return matrix finally ITensors.set_warn_order(prev_warn_order) From d322c1582296ed5cfcf7b2e03597e2305ab69bde Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Tue, 10 Sep 2024 15:23:53 +0200 Subject: [PATCH 08/24] Remove print, add comment --- src/affine.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/affine.jl b/src/affine.jl index 23a6673..c4e7c6a 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -25,7 +25,6 @@ function affine_transform_mpo( # Fill the MPO, taking care to not include auxiliary links at the edges mpo = MPO(R) spin_dims = ntuple(_ -> 2, M + N) - println(size.(tensors)) mpo[1] = ITensor(reshape(tensors[1], size(tensors[1], 2), spin_dims...), (link[1], outsite[1,:]..., insite[1,:]...)) for r in 2:R-1 @@ -171,12 +170,16 @@ function affine_mpo_to_matrix( prev_warn_order = ITensors.disable_warn_order() try mpo_contr = reduce(*, mpo) + + # Given some variables (x, y), we have to bring the indices in the + # order (xR, ..., x1, yR, ..., y1) in order to have y = (y1 ... yR)_2 + # once we reshape a column-major array and match the order of the + # variables in the full matrix. out_indices = vec(reverse(outsite, dims=1)) in_indices = vec(reverse(insite, dims=1)) - println(out_indices) - println(in_indices) tensor = Array(mpo_contr, out_indices..., in_indices...) - matrix = reshape(tensor, 1 << length(out_indices), 1 << length(in_indices)) + matrix = reshape(tensor, + 1 << length(out_indices), 1 << length(in_indices)) return matrix finally ITensors.set_warn_order(prev_warn_order) From 02bf35a8f30783fd4c614a714615adf43f42f0a1 Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Tue, 10 Sep 2024 17:13:33 +0200 Subject: [PATCH 09/24] Add active and passive transforms. Rationals we have to think about still a bit... --- src/Quantics.jl | 2 +- src/affine.jl | 89 +++++++++++++++++++++++++++++++++++++------------ 2 files changed, 69 insertions(+), 22 deletions(-) diff --git a/src/Quantics.jl b/src/Quantics.jl index 0763d3f..2699eac 100644 --- a/src/Quantics.jl +++ b/src/Quantics.jl @@ -13,7 +13,7 @@ import ITensors import ITensors.NDTensors: Tensor, BlockSparseTensor, blockview import SparseIR: Fermionic, Bosonic, Statistics -import LinearAlgebra: I +import LinearAlgebra: I, inv using StaticArrays import SparseArrays: sparse diff --git a/src/affine.jl b/src/affine.jl index c4e7c6a..bf40635 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -6,7 +6,8 @@ quantics representation. """ function affine_transform_mpo( outsite::AbstractMatrix{<:Index}, insite::AbstractMatrix{<:Index}, - A::AbstractMatrix{<:Integer}, b::AbstractVector{<:Integer} + A::AbstractMatrix{<:Union{Integer,Rational}}, + b::AbstractVector{<:Union{Integer,Rational}} ) R = size(outsite, 1) M, N = size(A) @@ -16,8 +17,7 @@ function affine_transform_mpo( throw(ArgumentError("vector is not correctly dimensioned")) # get the tensors so that we know how large the links must be - tensors = affine_transform_tensors( - R, SMatrix{M, N, Int}(A), SVector{M, Int}(b)) + tensors = affine_transform_tensors(R, A, b) # Create the links link = [Index(size(tensors[r], 2), tags="link $r") for r in 1:R-1] @@ -44,17 +44,19 @@ Compute vector of core tensors (constituent 4-way tensors) for a matrix product operator corresponding to the affine transformation `y = A*x + b`. """ function affine_transform_tensors( - R::Integer, A::AbstractMatrix{<:Integer}, - b::AbstractVector{<:Integer}) + R::Integer, A::AbstractMatrix{<:Union{Integer,Rational}}, + b::AbstractVector{<:Union{Integer,Rational}}) return affine_transform_tensors(Int(R), _affine_static_args(A, b)...) end function affine_transform_tensors( - R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int} + R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int}, denom::Int ) where {M, N} # Checks 0 <= R || throw(DomainError(R, "invalid value of the length R")) + denom == 1 || + throw(ArgumentError("XXX only denominator one supported for now")) # The output tensors are a collection of matrices, but their first two # dimensions (links) vary @@ -102,10 +104,10 @@ function affine_transform_core( carry::AbstractVector{SVector{M, Int}}; transform_out_carry::Function=identity ) where {M, N} - # The basic idea here is the following: we compute r = A*x + b + c for all - # "incoming" carrys d and all possible bit vectors, x ∈ {0,1}^N. Then we - # decompose r = 2*c + y, where y is again a bit vector, y ∈ {0,1}^M, and - # c is the "outgoing" carry, which may be negative. We then store this + # The basic idea here is the following: we compute r = A*x + b + c + # for all "incoming" carrys d and all possible bit vectors, x ∈ {0,1}^N. + # Then we decompose r = 2*c + y, where y is again a bit vector, y ∈ {0,1}^M, + # and c is the "outgoing" carry, which may be negative. We then store this # as something like out[d][c, x, y] = true. out = Dict{SVector{M, Int}, Array{Bool, 3}}() sizehint!(out, length(carry)) @@ -140,27 +142,35 @@ end Compute full transformation matrix for `y = A*x + b`. """ function affine_transform_matrix( - R::Integer, A::AbstractMatrix{<:Integer}, - b::AbstractVector{<:Integer}) + R::Integer, A::AbstractMatrix{<:Union{Integer,Rational}}, + b::AbstractVector{<:Union{Integer,Rational}}) return affine_transform_matrix(Int(R), _affine_static_args(A, b)...) end function affine_transform_matrix( - R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int}) where {M, N} + R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int}, + denom::Int) where {M, N} # Checks 0 <= R || throw(DomainError(R, "invalid value of the length R")) mask = (1 << R) - 1 - y_index = Vector{Int}(undef, 1 << (R * N)) + y_index = Int[] + x_index = Int[] for (ix, x) in enumerate(Iterators.product(ntuple(_ -> 0:mask, N)...)) - y = (A * SVector{N, Int}(x) + b) .& mask + z = A * SVector{N, Int}(x) + b + all(iszero, z .% denom) || + continue + + # XXX there is an ambiguity here: shall we divide first, then do + # modulo, or the other way 'round? + y = @. (z ÷ denom) & mask iy = digits_to_number(y, R) + 1 - @inbounds y_index[ix] = iy + push!(y_index, iy) + push!(x_index, ix) end - x_index = collect(1:1 << (R * N)) - values = ones(Bool, 1 << (R * N)) + values = ones(Bool, size(x_index)) return sparse(y_index, x_index, values, 1 << (R*M), 1 << (R*N)) end @@ -186,12 +196,49 @@ function affine_mpo_to_matrix( end end -function _affine_static_args( - A::AbstractMatrix{<:Integer}, b::AbstractVector{<:Integer}) +function _affine_static_args(A::AbstractMatrix{<:Union{Integer,Rational}}, + b::AbstractVector{<:Union{Integer,Rational}}) M, N = size(A) size(b, 1) == M || throw(ArgumentError("A and b have incompatible size")) - return SMatrix{M, N, Int}(A), SVector{M, Int}(b) + + # Factor out common denominator and pass + denom = lcm(mapreduce(denominator, lcm, A, init=1), + mapreduce(denominator, lcm, b, init=1)) + Ai = @. Int(denom * A) + bi = @. Int(denom * b) + + # Construct static matrix + return SMatrix{M, N, Int}(Ai), SVector{M, Int}(bi), denom +end + +""" + Ainv, binv = active_to_passive(A::AbstractMatrix, b::AbstractVector) + +Change active affine transformation `y = A*x + b` to the passive (inverse) one +`x = Ainv*y + binv`. + +Note that these transformations are not strict inverses of each other once you +put them on a discrete grid: in particular, Ainv may have rational coefficients +even though A is purely integer. In this case, the inverse transformation only +maps some of the points. +""" +function active_to_passive(A::AbstractMatrix{<:Union{Rational,Integer}}, + b::AbstractVector{<:Union{Rational,Integer}}) + return active_to_passive(Rational.(A), Rational.(b)) +end + +function active_to_passive( + A::AbstractMatrix{<:Rational}, b::AbstractVector{<:Rational}) + m, n = size(A) + T = [A b; zero(b)' 1] + + # XXX: we do not support pseudo-inverses (LinearAlgebbra cannot do + # this yet over the Rationals). + Tinv = inv(T) + Ainv = Tinv[1:m, 1:n] + binv = Tinv[1:m, n+1] + return Ainv, binv end """ From 2be3d74e928c27d2ed960525f65e551f6745ca14 Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Tue, 14 Jan 2025 18:16:19 +0100 Subject: [PATCH 10/24] Correct docstring. --- Project.toml | 2 +- src/affine.jl | 50 ++++++++++++++++++++++++++++++-------------------- 2 files changed, 31 insertions(+), 21 deletions(-) diff --git a/Project.toml b/Project.toml index 3e8ca00..1ad0342 100644 --- a/Project.toml +++ b/Project.toml @@ -27,7 +27,7 @@ FastMPOContractions = "^0.2.2" ITensorMPS = "0.3.1" ITensors = "0.7" PartitionedMPSs = "0.5.2" -QuanticsTCI = "0.7.0" +QuanticsTCI = "0.7" SparseIR = "^0.96, 0.97, 1" StaticArrays = "1" TCIITensorConversion = "0.2.0" diff --git a/src/affine.jl b/src/affine.jl index dd082ac..9079944 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -1,28 +1,38 @@ """ - affine_transform_mpo(outsite, insite, A, b) - -Constructs an ITensor matrix product state for the affine transformation `y = A*x + b` in the -quantics representation. - -# Arguments -- `outsite::AbstractMatrix{<:Index}`: The output site indices for the ITensor. -- `insite::AbstractMatrix{<:Index}`: The input site indices for the ITensor. -- `A::AbstractMatrix{<:Union{Integer,Rational}}`: The matrix representing the linear transformation. -- `b::AbstractVector{<:Union{Integer,Rational}}`: The vector representing the translation in the affine transformation. -- `little_endian::Bool=false`: Whether the input and output bits are in little-endian order. - -# Returns -- `MPO`: The resulting matrix product operator representing the affine transformation. + affine_transform_mpo(y, x, A, b) + +Construct and return ITensor matrix product state for the affine transformation +`y = A*x + b` in a (fused) quantics representation: + + y[1,1] ... y[1,M] y[2,1] ... y[2,M] y[R,1] ... y[R,M] + __|__________|__ __|__________|__ __|__________|__ + | | | | | | + | T[1] |---| T[2] |--- ... ---| T[R] | + |________________| |________________| |________________| + | | | | | | + x[1,1] ... x[1,N] x[2,1] ... x[2,N] x[R,1] ... x[R,N] + + +Arguments +--------- +- `y`: An `R × M` matrix of ITensor indices, where `y[r,m]` corresponds to + the `r`-th length scale of the `m`-th output variable. +- `x`: An `R × N` matrix of ITensor indices, where `x[r,n]` corresponds to + the `r`-th length scale of the `n`-th input variable. +- `A`: An `M × N` rational matrix representing the linear transformation. +- `b`: An `M` reational vector representing the translation. """ function affine_transform_mpo( - outsite::AbstractMatrix{<:Index}, insite::AbstractMatrix{<:Index}, + y::AbstractMatrix{<:Index}, x::AbstractMatrix{<:Index}, A::AbstractMatrix{<:Union{Integer,Rational}}, b::AbstractVector{<:Union{Integer,Rational}} )::MPO - R = size(outsite, 1) + R = size(y, 1) M, N = size(A) - size(insite) == (R, N) || + size(x) == (R, N) || throw(ArgumentError("insite is not correctly dimensioned")) + size(y) == (R, M) || + throw(ArgumentError("outsite is not correctly dimensioned")) size(b) == (M,) || throw(ArgumentError("vector is not correctly dimensioned")) @@ -36,14 +46,14 @@ function affine_transform_mpo( mpo = MPO(R) spin_dims = ntuple(_ -> 2, M + N) mpo[1] = ITensor(reshape(tensors[1], size(tensors[1], 2), spin_dims...), - (link[1], outsite[1,:]..., insite[1,:]...)) + (link[1], y[1,:]..., x[1,:]...)) for r in 2:R-1 newshape = size(tensors[r])[1:2]..., spin_dims... mpo[r] = ITensor(reshape(tensors[r], newshape), - (link[r-1], link[r], outsite[r,:]..., insite[r,:]...)) + (link[r-1], link[r], y[r,:]..., x[r,:]...)) end mpo[R] = ITensor(reshape(tensors[R], size(tensors[R], 1), spin_dims...), - (link[R-1], outsite[R,:]..., insite[R,:]...)) + (link[R-1], y[R,:]..., x[R,:]...)) return mpo end From c059d7c991c37848fa715e5e96401126ca85c74b Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Wed, 15 Jan 2025 15:41:24 +0100 Subject: [PATCH 11/24] Prepare affine transforms. --- src/affine.jl | 36 ++++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/src/affine.jl b/src/affine.jl index 9079944..0e7cae6 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -61,7 +61,8 @@ end affine_transform_tensors(R, A, b) Compute vector of core tensors (constituent 4-way tensors) for a matrix product -operator corresponding to the affine transformation `y = A*x + b`. +operator corresponding to one of affine transformation `y = A*x + b` with +rational `A` and `b` """ function affine_transform_tensors( R::Integer, A::AbstractMatrix{<:Union{Integer,Rational}}, @@ -70,13 +71,13 @@ function affine_transform_tensors( end function affine_transform_tensors( - R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int}, denom::Int + R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int}, s::Int ) where {M, N} # Checks 0 <= R || throw(DomainError(R, "invalid value of the length R")) - denom == 1 || - throw(ArgumentError("XXX only denominator one supported for now")) + s == 1 || + throw(DomainError(s, "must be one for now")) # The output tensors are a collection of matrices, but their first two # dimensions (links) vary @@ -89,11 +90,16 @@ function affine_transform_tensors( # it out from the array bcurr = @. Bool(b & 1) - # Get tensor. For the first tensor, we discard the carry. - out_tf = (r == 1) ? zero : identity - new_carry, data = affine_transform_core(A, bcurr, carry, - transform_out_carry=out_tf) - tensors[r] = data + # Get tensor. + new_carry, data = affine_transform_core(A, bcurr, s, carry) + + # For the first tensor, we assume periodic boundary conditions, so + # we sum over all choices off the carry + if r == 1 + tensors[r] = sum(data, dims=1) + else + tensors[r] = data + end # Set carry to the next value carry = new_carry @@ -103,14 +109,13 @@ function affine_transform_tensors( end """ - core, out_carry = affine_transform_core(A, b, in_carry; - transform_out_carry=identity) + core, out_carry = affine_transform_core(A, b, s, in_carry) Construct core tensor `core` for an affine transform. The core tensor for an affine transformation is given by: core[d, c, iy, ix] = - 2 * out_carry[d] + y[iy] == A * x[ix] + b + in_carry[c] + 2 * out_carry[d] + s * y[iy] == A * x[ix] + b + in_carry[c] where `A`, a matrix of integers, and `b`, a vector of bits, which define the affine transform. `c` and `d` are indices into a set of integer vectors @@ -120,9 +125,8 @@ are the binary input and output vectors, respectively, of the affine transform. They are indexed in a "little-endian" fashion. """ function affine_transform_core( - A::SMatrix{M, N, Int}, b::SVector{M, Bool}, - carry::AbstractVector{SVector{M, Int}}; - transform_out_carry::Function=identity + A::SMatrix{M, N, Int}, b::SVector{M, Bool}, s::Int, + carry::AbstractVector{SVector{M, Int}} ) where {M, N} # The basic idea here is the following: we compute r = A*x + b + c # for all "incoming" carrys d and all possible bit vectors, x ∈ {0,1}^N. @@ -135,7 +139,7 @@ function affine_transform_core( for (x_index, x) in enumerate(Iterators.product(ntuple(_ -> 0:1, N)...)) r = A * SVector{N, Bool}(x) + b + SVector{M, Int}(c) y = @. Bool(r & 1) - d::SVector{M, Int} = transform_out_carry(r .>> 1) + d = r .>> 1 y_index = digits_to_number(y) + 1 d_mat = get!(out, d) do From c426819fc46024fa5d77744e380f9cb7da27650a Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Wed, 15 Jan 2025 16:59:29 +0100 Subject: [PATCH 12/24] Add modular inverse. --- src/affine.jl | 27 +++++++++++++++++++++++++++ test/affine_tests.jl | 10 ++++++++++ 2 files changed, 37 insertions(+) diff --git a/src/affine.jl b/src/affine.jl index 0e7cae6..ddfa02c 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -287,3 +287,30 @@ digits_to_number(v::AbstractVector{<:Integer}, bits::Integer) = throw(DomainError(v[1], "invalid digit in base $(1 << bits)")) return _digits_to_number(v[2:end], bits) << bits | v[1] end + +""" + modular_inverse(b, m) + +Compute the multiplicative inverse of an integer `b` modulo `2^m`, i.e., find +and return an integer `x` such that: + + x * b ≡ 1 (mod 2^m) +""" +function modular_inverse(b::Integer, m::Integer) + 0 <= m <= 8 * sizeof(b) || + throw(DomainError(m, "invalid number of bits")) + isodd(b) || + throw(DomainError(b, "only odd numbers have inverses mod power of 2")) + + # Use Dusse and Kaliski's algorithm, as reproduced in + # Arazi and Qi, IEEE Trans. Comput. 57, 10 (2008), Algorithm 1 + mask = one(b) + y = one(b) + for _ in 2:m + mask <<= 1 + if (b * y) & mask != 0 + y |= mask + end + end + return y +end diff --git a/test/affine_tests.jl b/test/affine_tests.jl index 8c328fb..a98db0b 100644 --- a/test/affine_tests.jl +++ b/test/affine_tests.jl @@ -4,6 +4,16 @@ using Quantics using LinearAlgebra + @testset "mod_inverse" begin + @test_throws DomainError Quantics.modular_inverse(423, -1) + @test_throws DomainError Quantics.modular_inverse(80, 10) + + @test Quantics.modular_inverse(1, 5) == 1 + @test Quantics.modular_inverse(77, 8) isa Signed + @test (Quantics.modular_inverse(77, 8) * 77) & 0xFF == 1 + @test (Quantics.modular_inverse(-49, 7) * -49) & 0x7F == 1 + end + vars = ("x", "y", "z") insite = [Index(2; tags="$v$l") for l in 1:10, v in vars] outsite = [Index(2; tags="$v$l")' for l in 1:10, v in vars] From b85128ff28e4c2c6632671ffb637525cb7c87a3c Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Wed, 15 Jan 2025 17:31:52 +0100 Subject: [PATCH 13/24] Allow any odd denominator in affine transform. This basically works because the reciprocal of an odd integer is again integer as long as you do in modulo arithmetic, so we can reduce it to the case we already consider. --- src/affine.jl | 33 ++++++++++++++++++++------------- test/affine_tests.jl | 13 +++++++++++++ 2 files changed, 33 insertions(+), 13 deletions(-) diff --git a/src/affine.jl b/src/affine.jl index ddfa02c..c4f8227 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -74,11 +74,19 @@ function affine_transform_tensors( R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int}, s::Int ) where {M, N} # Checks - 0 <= R || + 0 <= R <= 8 * sizeof(Int) || throw(DomainError(R, "invalid value of the length R")) - s == 1 || + isodd(s) || throw(DomainError(s, "must be one for now")) + # We are currently assuming periodic boundary conditions and s being odd. + # Then there is a multiplicative inverse such that inv_s * s ≡ 1 (mod 2^R) + # This lets us rewrite: 1/s * (A*x + b) to inv_s*(A*x + b) + base = 1 << R + inv_s = modular_inverse(s, base) + A = inv_s * A + b = inv_s * b + # The output tensors are a collection of matrices, but their first two # dimensions (links) vary tensors = Vector{Array{Bool, 4}}(undef, R) @@ -91,7 +99,7 @@ function affine_transform_tensors( bcurr = @. Bool(b & 1) # Get tensor. - new_carry, data = affine_transform_core(A, bcurr, s, carry) + new_carry, data = affine_transform_core(A, bcurr, carry) # For the first tensor, we assume periodic boundary conditions, so # we sum over all choices off the carry @@ -125,9 +133,10 @@ are the binary input and output vectors, respectively, of the affine transform. They are indexed in a "little-endian" fashion. """ function affine_transform_core( - A::SMatrix{M, N, Int}, b::SVector{M, Bool}, s::Int, + A::SMatrix{M, N, Int}, b::SVector{M, Bool}, carry::AbstractVector{SVector{M, Int}} ) where {M, N} + # The basic idea here is the following: we compute r = A*x + b + c # for all "incoming" carrys d and all possible bit vectors, x ∈ {0,1}^N. # Then we decompose r = 2*c + y, where y is again a bit vector, y ∈ {0,1}^M, @@ -173,23 +182,21 @@ end function affine_transform_matrix( R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int}, - denom::Int) where {M, N} + s::Int) where {M, N} # Checks 0 <= R || throw(DomainError(R, "invalid value of the length R")) + isodd(s) || + throw(DomainError(s, "right now we only support odd s")) - mask = (1 << R) - 1 + mask = ~(~0 << R) + inv_s = modular_inverse(s, R) y_index = Int[] x_index = Int[] for (ix, x) in enumerate(Iterators.product(ntuple(_ -> 0:mask, N)...)) - z = A * SVector{N, Int}(x) + b - all(iszero, z .% denom) || - continue - - # XXX there is an ambiguity here: shall we divide first, then do - # modulo, or the other way 'round? - y = @. (z ÷ denom) & mask + yfull = inv_s * (A * SVector{N, Int}(x) + b) + y = yfull .& mask iy = digits_to_number(y, R) + 1 push!(y_index, iy) push!(x_index, ix) diff --git a/test/affine_tests.jl b/test/affine_tests.jl index a98db0b..17d5ed8 100644 --- a/test/affine_tests.jl +++ b/test/affine_tests.jl @@ -70,4 +70,17 @@ @test T == Trec end + @testset "compare_denom_odd" begin + A = reshape([1//3], 1, 1) + b = [0] + R = 6 + + T = Quantics.affine_transform_matrix(R, A, b) + mpo = Quantics.affine_transform_mpo( + outsite[1:R, 1:1], insite[1:R, 1:1], A, b) + Trec = Quantics.affine_mpo_to_matrix( + outsite[1:R, 1:1], insite[1:R, 1:1], mpo) + @test T == Trec + end + end From c71c65eb2976d5bbdd0f16e457ee9993b1129673 Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Wed, 15 Jan 2025 18:38:49 +0100 Subject: [PATCH 14/24] Allow different boundary conditions. --- src/affine.jl | 34 +++++++++++++++++++++++++++------- 1 file changed, 27 insertions(+), 7 deletions(-) diff --git a/src/affine.jl b/src/affine.jl index c4f8227..6fd7c57 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -170,19 +170,33 @@ function affine_transform_core( end """ - affine_transform_matrix(R, A, b) + affine_transform_matrix(R, A, b; periodic=true) -Compute full transformation matrix for `y = A*x + b`. +Compute full transformation matrix for the affine transformation `y = A*x + b`, +where `y` is a `M`-vector and `x` is `N`-vector, and each component is in +`{0, 1, ..., 2^R-1}`. `A` is a rational `M × N` matrix and `b` is a rational +`N`-vector. + +Return a boolean sparse `2^(R*M) × 2^(R*N)` matrix `A`, which has true entries +whereever the condition is satisfied. The element indices `A[iy, ix]` are +mapped to `x` and `y` as follows: + + iy = 1 + y[1] + y[2] * 2^R + y[3] * 2^(2R) + ... + y[M] * 2^((M-1)*R) + ix = 1 + x[1] + x[2] * 2^R + x[3] * 2^(2R) + ... + x[N] * 2^((N-1)*R) + +If `periodic` is true, then periodic boundary conditions, `y[i] + 2^R = y[i]`, +are used. """ function affine_transform_matrix( R::Integer, A::AbstractMatrix{<:Union{Integer,Rational}}, - b::AbstractVector{<:Union{Integer,Rational}}) - return affine_transform_matrix(Int(R), _affine_static_args(A, b)...) + b::AbstractVector{<:Union{Integer,Rational}}; periodic::Bool=true + ) + return affine_transform_matrix(Int(R), _affine_static_args(A, b)..., periodic) end function affine_transform_matrix( R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int}, - s::Int) where {M, N} + s::Int, periodic::Bool) where {M, N} # Checks 0 <= R || throw(DomainError(R, "invalid value of the length R")) @@ -195,8 +209,14 @@ function affine_transform_matrix( x_index = Int[] for (ix, x) in enumerate(Iterators.product(ntuple(_ -> 0:mask, N)...)) - yfull = inv_s * (A * SVector{N, Int}(x) + b) - y = yfull .& mask + v = A * SVector{N, Int}(x) + b + if periodic + v *= inv_s + y = v .& mask + else + iszero(v .% s) || continue + y = v .÷ s + end iy = digits_to_number(y, R) + 1 push!(y_index, iy) push!(x_index, ix) From 6b02f0806408bca03f2b39c8bc19c264ef81ea4d Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Thu, 16 Jan 2025 09:58:33 +0100 Subject: [PATCH 15/24] Rename modular_inverse -> invmod_mod2. This more closely aligns to the Julia 1.11 function `invmod`. --- src/affine.jl | 11 +++++++---- test/affine_tests.jl | 14 +++++++------- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/src/affine.jl b/src/affine.jl index 6fd7c57..f407394 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -83,7 +83,7 @@ function affine_transform_tensors( # Then there is a multiplicative inverse such that inv_s * s ≡ 1 (mod 2^R) # This lets us rewrite: 1/s * (A*x + b) to inv_s*(A*x + b) base = 1 << R - inv_s = modular_inverse(s, base) + inv_s = invmod_pow2(s, base) A = inv_s * A b = inv_s * b @@ -204,7 +204,7 @@ function affine_transform_matrix( throw(DomainError(s, "right now we only support odd s")) mask = ~(~0 << R) - inv_s = modular_inverse(s, R) + inv_s = invmod_pow2(s, R) y_index = Int[] x_index = Int[] @@ -316,14 +316,17 @@ digits_to_number(v::AbstractVector{<:Integer}, bits::Integer) = end """ - modular_inverse(b, m) + invmod_pow2(b, m) Compute the multiplicative inverse of an integer `b` modulo `2^m`, i.e., find and return an integer `x` such that: x * b ≡ 1 (mod 2^m) + +For Julia 1.11 and onwards, this is equivalent to, but should be faster than, +`invmod(b, 1 << m)`. """ -function modular_inverse(b::Integer, m::Integer) +function invmod_pow2(b::Integer, m::Integer) 0 <= m <= 8 * sizeof(b) || throw(DomainError(m, "invalid number of bits")) isodd(b) || diff --git a/test/affine_tests.jl b/test/affine_tests.jl index 17d5ed8..4282129 100644 --- a/test/affine_tests.jl +++ b/test/affine_tests.jl @@ -4,14 +4,14 @@ using Quantics using LinearAlgebra - @testset "mod_inverse" begin - @test_throws DomainError Quantics.modular_inverse(423, -1) - @test_throws DomainError Quantics.modular_inverse(80, 10) + @testset "invmod_pow2" begin + @test_throws DomainError Quantics.invmod_pow2(423, -1) + @test_throws DomainError Quantics.invmod_pow2(80, 10) - @test Quantics.modular_inverse(1, 5) == 1 - @test Quantics.modular_inverse(77, 8) isa Signed - @test (Quantics.modular_inverse(77, 8) * 77) & 0xFF == 1 - @test (Quantics.modular_inverse(-49, 7) * -49) & 0x7F == 1 + @test Quantics.invmod_pow2(1, 5) == 1 + @test Quantics.invmod_pow2(77, 8) isa Signed + @test (Quantics.invmod_pow2(77, 8) * 77) & 0xFF == 1 + @test (Quantics.invmod_pow2(-49, 7) * -49) & 0x7F == 1 end vars = ("x", "y", "z") From 7f71a5430ea89e6d5d2c6079198417730ac86cf9 Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Thu, 16 Jan 2025 13:01:29 +0100 Subject: [PATCH 16/24] Boundary conditions --- src/affine.jl | 54 +++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 39 insertions(+), 15 deletions(-) diff --git a/src/affine.jl b/src/affine.jl index f407394..4362cbe 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -1,3 +1,30 @@ +""" + AbstractBoundaryConditions + +Boundary conditions for the QTT to use. Use `OpenBoundaryCondtions`` for open +boundaries and `PeriodicBoundaryConditions` for periodic ones. +""" +abstract type AbstractBoundaryConditions end + +struct PeriodicBoundaryConditions <: AbstractBoundaryConditions + R::Int +end + +@inline function divide(x, s, bc::PeriodicBoundaryConditions) + mask = ~(~0 << bc.R) + inv_s = invmod_pow2(s, bc.R) + return (x * inv_s) .& mask +end + +struct OpenBoundaryConditions <: AbstractBoundaryConditions + R::Int +end + +function divide(x, s, bc::OpenBoundaryConditions) + invmask = ~0 << bc.R + return iszero(x .& invmask) && iszero(x .% s) ? x .÷ s : nothing +end + """ affine_transform_mpo(y, x, A, b) @@ -170,7 +197,7 @@ function affine_transform_core( end """ - affine_transform_matrix(R, A, b; periodic=true) + affine_transform_matrix(R, A, b, [boundary]) Compute full transformation matrix for the affine transformation `y = A*x + b`, where `y` is a `M`-vector and `x` is `N`-vector, and each component is in @@ -184,19 +211,19 @@ mapped to `x` and `y` as follows: iy = 1 + y[1] + y[2] * 2^R + y[3] * 2^(2R) + ... + y[M] * 2^((M-1)*R) ix = 1 + x[1] + x[2] * 2^R + x[3] * 2^(2R) + ... + x[N] * 2^((N-1)*R) -If `periodic` is true, then periodic boundary conditions, `y[i] + 2^R = y[i]`, -are used. +`boundary` specifies the type of boundary conditions. """ function affine_transform_matrix( R::Integer, A::AbstractMatrix{<:Union{Integer,Rational}}, - b::AbstractVector{<:Union{Integer,Rational}}; periodic::Bool=true + b::AbstractVector{<:Union{Integer,Rational}}, + boundary::AbstractBoundaryConditions=PeriodicBoundaryConditions(R) ) - return affine_transform_matrix(Int(R), _affine_static_args(A, b)..., periodic) + return affine_transform_matrix(Int(R), _affine_static_args(A, b)..., boundary) end function affine_transform_matrix( R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int}, - s::Int, periodic::Bool) where {M, N} + s::Int, boundary::AbstractBoundaryConditions) where {M, N} # Checks 0 <= R || throw(DomainError(R, "invalid value of the length R")) @@ -204,22 +231,19 @@ function affine_transform_matrix( throw(DomainError(s, "right now we only support odd s")) mask = ~(~0 << R) - inv_s = invmod_pow2(s, R) y_index = Int[] x_index = Int[] for (ix, x) in enumerate(Iterators.product(ntuple(_ -> 0:mask, N)...)) v = A * SVector{N, Int}(x) + b - if periodic - v *= inv_s - y = v .& mask + y = divide(v, s, boundary) + if isnothing(y) + continue else - iszero(v .% s) || continue - y = v .÷ s + iy = digits_to_number(y, R) + 1 + push!(y_index, iy) + push!(x_index, ix) end - iy = digits_to_number(y, R) + 1 - push!(y_index, iy) - push!(x_index, ix) end values = ones(Bool, size(x_index)) return sparse(y_index, x_index, values, 1 << (R*M), 1 << (R*N)) From 8142d73a8c93682c34f266c8fa5de045a3f8676d Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Thu, 16 Jan 2025 16:39:44 +0100 Subject: [PATCH 17/24] Do not store R in BC. --- src/affine.jl | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/src/affine.jl b/src/affine.jl index 4362cbe..06517d9 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -6,23 +6,19 @@ boundaries and `PeriodicBoundaryConditions` for periodic ones. """ abstract type AbstractBoundaryConditions end -struct PeriodicBoundaryConditions <: AbstractBoundaryConditions - R::Int -end +struct PeriodicBoundaryConditions <: AbstractBoundaryConditions end -@inline function divide(x, s, bc::PeriodicBoundaryConditions) - mask = ~(~0 << bc.R) - inv_s = invmod_pow2(s, bc.R) - return (x * inv_s) .& mask +@inline function matching_y(v, s, R, ::PeriodicBoundaryConditions) + mask = ~(~0 << R) + inv_s = invmod_pow2(s, R) + return (v * inv_s) .& mask end -struct OpenBoundaryConditions <: AbstractBoundaryConditions - R::Int -end +struct OpenBoundaryConditions <: AbstractBoundaryConditions end -function divide(x, s, bc::OpenBoundaryConditions) - invmask = ~0 << bc.R - return iszero(x .& invmask) && iszero(x .% s) ? x .÷ s : nothing +@inline function matching_y(v, s, R, ::OpenBoundaryConditions) + invmask = ~0 << R + return iszero(v .& invmask) && iszero(v .% s) ? v .÷ s : nothing end """ @@ -216,7 +212,7 @@ mapped to `x` and `y` as follows: function affine_transform_matrix( R::Integer, A::AbstractMatrix{<:Union{Integer,Rational}}, b::AbstractVector{<:Union{Integer,Rational}}, - boundary::AbstractBoundaryConditions=PeriodicBoundaryConditions(R) + boundary::AbstractBoundaryConditions=PeriodicBoundaryConditions() ) return affine_transform_matrix(Int(R), _affine_static_args(A, b)..., boundary) end @@ -233,10 +229,13 @@ function affine_transform_matrix( mask = ~(~0 << R) y_index = Int[] x_index = Int[] + sizehint!(y_index, 1 << (R * N)) + sizehint!(x_index, 1 << (R * N)) + # XXX: we actually need to iterate whatever dimension is largest for (ix, x) in enumerate(Iterators.product(ntuple(_ -> 0:mask, N)...)) v = A * SVector{N, Int}(x) + b - y = divide(v, s, boundary) + y = matching_y(v, s, R, boundary) if isnothing(y) continue else From a408205292ca62a3cec7f2ea6a1c675d4a4c5a1a Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Thu, 16 Jan 2025 17:17:44 +0100 Subject: [PATCH 18/24] Add boundary condition argument --- src/affine.jl | 48 +++++++++++++++++++++++++++++++----------------- 1 file changed, 31 insertions(+), 17 deletions(-) diff --git a/src/affine.jl b/src/affine.jl index 06517d9..39eca8b 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -14,6 +14,8 @@ struct PeriodicBoundaryConditions <: AbstractBoundaryConditions end return (v * inv_s) .& mask end +carry_weight(c, ::PeriodicBoundaryConditions) = true + struct OpenBoundaryConditions <: AbstractBoundaryConditions end @inline function matching_y(v, s, R, ::OpenBoundaryConditions) @@ -21,8 +23,10 @@ struct OpenBoundaryConditions <: AbstractBoundaryConditions end return iszero(v .& invmask) && iszero(v .% s) ? v .÷ s : nothing end +carry_weight(c, ::OpenBoundaryConditions) = iszero(c) + """ - affine_transform_mpo(y, x, A, b) + affine_transform_mpo(y, x, A, b, [boundary]) Construct and return ITensor matrix product state for the affine transformation `y = A*x + b` in a (fused) quantics representation: @@ -44,11 +48,13 @@ Arguments the `r`-th length scale of the `n`-th input variable. - `A`: An `M × N` rational matrix representing the linear transformation. - `b`: An `M` reational vector representing the translation. +- `boundary`: boundary conditions (defaults to `PeriodicBoundaryConditions()`) """ function affine_transform_mpo( y::AbstractMatrix{<:Index}, x::AbstractMatrix{<:Index}, A::AbstractMatrix{<:Union{Integer,Rational}}, - b::AbstractVector{<:Union{Integer,Rational}} + b::AbstractVector{<:Union{Integer,Rational}}, + boundary::AbstractBoundaryConditions=PeriodicBoundaryConditions() )::MPO R = size(y, 1) M, N = size(A) @@ -60,7 +66,7 @@ function affine_transform_mpo( throw(ArgumentError("vector is not correctly dimensioned")) # get the tensors so that we know how large the links must be - tensors = affine_transform_tensors(R, A, b) + tensors = affine_transform_tensors(R, A, b, boundary) # Create the links link = [Index(size(tensors[r], 2), tags="link $r") for r in 1:R-1] @@ -81,7 +87,7 @@ function affine_transform_mpo( end """ - affine_transform_tensors(R, A, b) + affine_transform_tensors(R, A, b, [boundary]) Compute vector of core tensors (constituent 4-way tensors) for a matrix product operator corresponding to one of affine transformation `y = A*x + b` with @@ -89,26 +95,33 @@ rational `A` and `b` """ function affine_transform_tensors( R::Integer, A::AbstractMatrix{<:Union{Integer,Rational}}, - b::AbstractVector{<:Union{Integer,Rational}}) - return affine_transform_tensors(Int(R), _affine_static_args(A, b)...) + b::AbstractVector{<:Union{Integer,Rational}}, + boundary::AbstractBoundaryConditions=PeriodicBoundaryConditions()) + return affine_transform_tensors( + Int(R), _affine_static_args(A, b)..., boundary) end function affine_transform_tensors( - R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int}, s::Int - ) where {M, N} + R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int}, s::Int, + boundary::AbstractBoundaryConditions) where {M, N} # Checks 0 <= R <= 8 * sizeof(Int) || throw(DomainError(R, "invalid value of the length R")) - isodd(s) || - throw(DomainError(s, "must be one for now")) # We are currently assuming periodic boundary conditions and s being odd. # Then there is a multiplicative inverse such that inv_s * s ≡ 1 (mod 2^R) # This lets us rewrite: 1/s * (A*x + b) to inv_s*(A*x + b) - base = 1 << R - inv_s = invmod_pow2(s, base) - A = inv_s * A - b = inv_s * b + if boundary isa PeriodicBoundaryConditions + isodd(s) || + throw(DomainError(s, "must be odd for now")) + base = 1 << R + inv_s = invmod_pow2(s, base) + A = inv_s * A + b = inv_s * b + else + isone(s) || + throw(DomainError(s, "must be one for now")) + end # The output tensors are a collection of matrices, but their first two # dimensions (links) vary @@ -124,10 +137,11 @@ function affine_transform_tensors( # Get tensor. new_carry, data = affine_transform_core(A, bcurr, carry) - # For the first tensor, we assume periodic boundary conditions, so - # we sum over all choices off the carry if r == 1 - tensors[r] = sum(data, dims=1) + # For the first tensor, we examine the carry to see which elements + # contribute with which weight + weights = map(c -> carry_weight(c, boundary), new_carry) + tensors[r] = sum(data .* weights, dims=1) else tensors[r] = data end From 5a6d1213df5373c6913df9167ad47fd8dc844cc2 Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Fri, 17 Jan 2025 11:34:17 +0100 Subject: [PATCH 19/24] Support arbitrary affine transformation. But right now the QTTs are not optimal in terms of bond dimension. --- src/affine.jl | 90 ++++++++++++++++++++++++++++++++------------------- 1 file changed, 57 insertions(+), 33 deletions(-) diff --git a/src/affine.jl b/src/affine.jl index 39eca8b..de71d69 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -105,24 +105,9 @@ function affine_transform_tensors( R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int}, s::Int, boundary::AbstractBoundaryConditions) where {M, N} # Checks - 0 <= R <= 8 * sizeof(Int) || + 0 <= R * max(M, N) <= 8 * sizeof(Int) || throw(DomainError(R, "invalid value of the length R")) - # We are currently assuming periodic boundary conditions and s being odd. - # Then there is a multiplicative inverse such that inv_s * s ≡ 1 (mod 2^R) - # This lets us rewrite: 1/s * (A*x + b) to inv_s*(A*x + b) - if boundary isa PeriodicBoundaryConditions - isodd(s) || - throw(DomainError(s, "must be odd for now")) - base = 1 << R - inv_s = invmod_pow2(s, base) - A = inv_s * A - b = inv_s * b - else - isone(s) || - throw(DomainError(s, "must be one for now")) - end - # The output tensors are a collection of matrices, but their first two # dimensions (links) vary tensors = Vector{Array{Bool, 4}}(undef, R) @@ -135,7 +120,10 @@ function affine_transform_tensors( bcurr = @. Bool(b & 1) # Get tensor. - new_carry, data = affine_transform_core(A, bcurr, carry) + new_carry, data = affine_transform_core(A, bcurr, s, carry) + + # XXX do pruning: cut away carries that are dead ends in further + # tensors if r == 1 # For the first tensor, we examine the carry to see which elements @@ -160,7 +148,7 @@ Construct core tensor `core` for an affine transform. The core tensor for an affine transformation is given by: core[d, c, iy, ix] = - 2 * out_carry[d] + s * y[iy] == A * x[ix] + b + in_carry[c] + 2 * out_carry[d] == A * x[ix] + b - s * y[iy] + in_carry[c] where `A`, a matrix of integers, and `b`, a vector of bits, which define the affine transform. `c` and `d` are indices into a set of integer vectors @@ -170,28 +158,64 @@ are the binary input and output vectors, respectively, of the affine transform. They are indexed in a "little-endian" fashion. """ function affine_transform_core( - A::SMatrix{M, N, Int}, b::SVector{M, Bool}, + A::SMatrix{M, N, Int}, b::SVector{M, Bool}, s::Int, carry::AbstractVector{SVector{M, Int}} ) where {M, N} - # The basic idea here is the following: we compute r = A*x + b + c - # for all "incoming" carrys d and all possible bit vectors, x ∈ {0,1}^N. - # Then we decompose r = 2*c + y, where y is again a bit vector, y ∈ {0,1}^M, - # and c is the "outgoing" carry, which may be negative. We then store this - # as something like out[d][c, x, y] = true. + # The basic idea here is the following: we check + # + # A*x + b - s*y + c == 2*d + # + # for all "incoming" carrys c and all possible bit vectors, x ∈ {0,1}^N. + # and y ∈ {0,1}^M for some outgoing carry d, which may be negative. + # We then store this as something like out[d, c, x, y]. out = Dict{SVector{M, Int}, Array{Bool, 3}}() sizehint!(out, length(carry)) + + all_x = Iterators.product(ntuple(_ -> 0:1, N)...) + all_y = Iterators.product(ntuple(_ -> 0:1, M)...) for (c_index, c) in enumerate(carry) - for (x_index, x) in enumerate(Iterators.product(ntuple(_ -> 0:1, N)...)) - r = A * SVector{N, Bool}(x) + b + SVector{M, Int}(c) - y = @. Bool(r & 1) - d = r .>> 1 - y_index = digits_to_number(y) + 1 - - d_mat = get!(out, d) do - return zeros(Bool, length(carry), 1 << M, 1 << N) + for (x_index, x) in enumerate(all_x) + z = A * SVector{N, Bool}(x) + b + SVector{M, Int}(c) + + if isodd(s) + # if s is odd, then there is a unique y which solves satisfies + # above condition (simply the lowest bit) + y = @. Bool(z & 1) + y_index = digits_to_number(y) + 1 + + # Correct z and compute carry + z = @. z - s * y + d = z .>> 1 + + # Store this + d_mat = get!(out, d) do + return zeros(Bool, length(carry), 1 << M, 1 << N) + end + @inbounds d_mat[c_index, y_index, x_index] = true + else + # if s instead even, then the conditions for x and y decouple. + # since y cannot touch the lowest bit, z must already be even. + all(iseven, z) || + continue + + # y can take any value at this point. This may lead to + # bonds that do not contribute because they become dead ends + # in a subsequent tensor (no valid outgoing carries). We cannot + # decide this here, but must prune those "branches" from the + # right once in the driver routine (affine_transform_tensors). + for (y_index, y) in enumerate(all_y) + # Correct z and compute carry + z = @. z - s * y + d = z .>> 1 + + # Store this + d_mat = get!(out, d) do + return zeros(Bool, length(carry), 1 << M, 1 << N) + end + @inbounds d_mat[c_index, y_index, x_index] = true + end end - @inbounds d_mat[c_index, y_index, x_index] = true end end From 1a19f29747937357e8617da0e17b5f5b7cffeb4b Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Fri, 17 Jan 2025 13:20:30 +0100 Subject: [PATCH 20/24] Now we can use arbitrary divisors. --- src/affine.jl | 85 ++++++++++++-------------------------------- test/affine_tests.jl | 28 ++++++--------- 2 files changed, 33 insertions(+), 80 deletions(-) diff --git a/src/affine.jl b/src/affine.jl index de71d69..e5b1727 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -8,22 +8,21 @@ abstract type AbstractBoundaryConditions end struct PeriodicBoundaryConditions <: AbstractBoundaryConditions end -@inline function matching_y(v, s, R, ::PeriodicBoundaryConditions) +@inline function equiv(x::SVector, y::SVector, R::Int, ::PeriodicBoundaryConditions) mask = ~(~0 << R) - inv_s = invmod_pow2(s, R) - return (v * inv_s) .& mask + return iszero((x - y) .& mask) end -carry_weight(c, ::PeriodicBoundaryConditions) = true +carry_weight(c::SVector, ::PeriodicBoundaryConditions) = true struct OpenBoundaryConditions <: AbstractBoundaryConditions end -@inline function matching_y(v, s, R, ::OpenBoundaryConditions) +@inline function equiv(x::SVector, y::SVector, R::Int, ::OpenBoundaryConditions) invmask = ~0 << R - return iszero(v .& invmask) && iszero(v .% s) ? v .÷ s : nothing + return iszero(x .& invmask) && x == y end -carry_weight(c, ::OpenBoundaryConditions) = iszero(c) +carry_weight(c::SVector, ::OpenBoundaryConditions) = iszero(c) """ affine_transform_mpo(y, x, A, b, [boundary]) @@ -162,6 +161,10 @@ function affine_transform_core( carry::AbstractVector{SVector{M, Int}} ) where {M, N} + # Otherwise we have to reverse the indexing of x and y + M <= N || + throw(ArgumentError("expect wide transformation matrix")) + # The basic idea here is the following: we check # # A*x + b - s*y + c == 2*d @@ -259,27 +262,24 @@ function affine_transform_matrix( R::Int, A::SMatrix{M, N, Int}, b::SVector{M, Int}, s::Int, boundary::AbstractBoundaryConditions) where {M, N} # Checks - 0 <= R || + 0 <= R * max(M, N) <= 8 * sizeof(Int) || throw(DomainError(R, "invalid value of the length R")) - isodd(s) || - throw(DomainError(s, "right now we only support odd s")) mask = ~(~0 << R) y_index = Int[] x_index = Int[] - sizehint!(y_index, 1 << (R * N)) - sizehint!(x_index, 1 << (R * N)) + sizehint!(y_index, 1 << (R * max(M, N))) + sizehint!(x_index, 1 << (R * max(M, N))) - # XXX: we actually need to iterate whatever dimension is largest - for (ix, x) in enumerate(Iterators.product(ntuple(_ -> 0:mask, N)...)) + all_x = Iterators.product(ntuple(_ -> 0:mask, N)...) + all_y = Iterators.product(ntuple(_ -> 0:mask, M)...) + for (ix, x) in enumerate(all_x) v = A * SVector{N, Int}(x) + b - y = matching_y(v, s, R, boundary) - if isnothing(y) - continue - else - iy = digits_to_number(y, R) + 1 - push!(y_index, iy) - push!(x_index, ix) + for (iy, y) in enumerate(all_y) + if equiv(v, s * SVector{M, Int}(y), R, boundary) + push!(y_index, iy) + push!(x_index, ix) + end end end values = ones(Bool, size(x_index)) @@ -355,53 +355,12 @@ end """ digits_to_number(v::AbstractVector{Bool}) - digits_to_number(v::AbstractVector{<:Integer}, bits::Integer) Converts a vector of digits, starting with the least significant digit, to -a number. If the digits are boolean, then they are interpreted in base-2, -otherwise, in base `2^bits`. +a number. """ digits_to_number(v::AbstractVector{Bool}) = _digits_to_number(Tuple(v)) -digits_to_number(v::AbstractVector{<:Integer}, bits::Integer) = - _digits_to_number(Int.(Tuple(v)), Int(bits)) @inline _digits_to_number(v::Tuple{}) = 0 @inline _digits_to_number(v::Tuple{Bool, Vararg{Bool}}) = _digits_to_number(v[2:end]) << 1 | v[1] -@inline _digits_to_number(v::Tuple{}, bits::Int) = 0 -@inline function _digits_to_number(v::Tuple{Int, Vararg{Int}}, bits::Int) - mask = (~0) << bits - iszero(v[1] & mask) || - throw(DomainError(v[1], "invalid digit in base $(1 << bits)")) - return _digits_to_number(v[2:end], bits) << bits | v[1] -end - -""" - invmod_pow2(b, m) - -Compute the multiplicative inverse of an integer `b` modulo `2^m`, i.e., find -and return an integer `x` such that: - - x * b ≡ 1 (mod 2^m) - -For Julia 1.11 and onwards, this is equivalent to, but should be faster than, -`invmod(b, 1 << m)`. -""" -function invmod_pow2(b::Integer, m::Integer) - 0 <= m <= 8 * sizeof(b) || - throw(DomainError(m, "invalid number of bits")) - isodd(b) || - throw(DomainError(b, "only odd numbers have inverses mod power of 2")) - - # Use Dusse and Kaliski's algorithm, as reproduced in - # Arazi and Qi, IEEE Trans. Comput. 57, 10 (2008), Algorithm 1 - mask = one(b) - y = one(b) - for _ in 2:m - mask <<= 1 - if (b * y) & mask != 0 - y |= mask - end - end - return y -end diff --git a/test/affine_tests.jl b/test/affine_tests.jl index 4282129..e0d3d44 100644 --- a/test/affine_tests.jl +++ b/test/affine_tests.jl @@ -4,19 +4,10 @@ using Quantics using LinearAlgebra - @testset "invmod_pow2" begin - @test_throws DomainError Quantics.invmod_pow2(423, -1) - @test_throws DomainError Quantics.invmod_pow2(80, 10) - - @test Quantics.invmod_pow2(1, 5) == 1 - @test Quantics.invmod_pow2(77, 8) isa Signed - @test (Quantics.invmod_pow2(77, 8) * 77) & 0xFF == 1 - @test (Quantics.invmod_pow2(-49, 7) * -49) & 0x7F == 1 - end - vars = ("x", "y", "z") insite = [Index(2; tags="$v$l") for l in 1:10, v in vars] outsite = [Index(2; tags="$v$l")' for l in 1:10, v in vars] + boundaries = Quantics.OpenBoundaryConditions(), Quantics.PeriodicBoundaryConditions() @testset "full" begin A = [1 0; 1 1] @@ -73,14 +64,17 @@ @testset "compare_denom_odd" begin A = reshape([1//3], 1, 1) b = [0] - R = 6 - T = Quantics.affine_transform_matrix(R, A, b) - mpo = Quantics.affine_transform_mpo( - outsite[1:R, 1:1], insite[1:R, 1:1], A, b) - Trec = Quantics.affine_mpo_to_matrix( - outsite[1:R, 1:1], insite[1:R, 1:1], mpo) - @test T == Trec + for R in [2, 3, 6] + for bc in boundaries + T = Quantics.affine_transform_matrix(R, A, b, bc) + mpo = Quantics.affine_transform_mpo( + outsite[1:R, 1:1], insite[1:R, 1:1], A, b, bc) + Trec = Quantics.affine_mpo_to_matrix( + outsite[1:R, 1:1], insite[1:R, 1:1], mpo) + @test T == Trec + end + end end end From 969820f19351acd74749144fedee4211f9478c88 Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Fri, 17 Jan 2025 15:35:27 +0100 Subject: [PATCH 21/24] Some case of even transforms. --- .gitignore | 1 + src/affine.jl | 8 +++----- test/affine_tests.jl | 32 ++++++++++++++++++++++++++------ 3 files changed, 30 insertions(+), 11 deletions(-) diff --git a/.gitignore b/.gitignore index 6c50e99..9bc2802 100644 --- a/.gitignore +++ b/.gitignore @@ -11,6 +11,7 @@ *.jl.mem .DS_Store /Manifest.toml +/notebook*/ /docs/Manifest.toml /docs/build/ docs/build diff --git a/src/affine.jl b/src/affine.jl index e5b1727..199d1a4 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -17,10 +17,7 @@ carry_weight(c::SVector, ::PeriodicBoundaryConditions) = true struct OpenBoundaryConditions <: AbstractBoundaryConditions end -@inline function equiv(x::SVector, y::SVector, R::Int, ::OpenBoundaryConditions) - invmask = ~0 << R - return iszero(x .& invmask) && x == y -end +equiv(x::SVector, y::SVector, R::Int, ::OpenBoundaryConditions) = x == y carry_weight(c::SVector, ::OpenBoundaryConditions) = iszero(c) @@ -277,6 +274,7 @@ function affine_transform_matrix( v = A * SVector{N, Int}(x) + b for (iy, y) in enumerate(all_y) if equiv(v, s * SVector{M, Int}(y), R, boundary) + #println("$y <- $x") push!(y_index, iy) push!(x_index, ix) end @@ -345,7 +343,7 @@ function active_to_passive( m, n = size(A) T = [A b; zero(b)' 1] - # XXX: we do not support pseudo-inverses (LinearAlgebbra cannot do + # XXX: we do not support pseudo-inverses (LinearAlgebra cannot do # this yet over the Rationals). Tinv = inv(T) Ainv = Tinv[1:m, 1:n] diff --git a/test/affine_tests.jl b/test/affine_tests.jl index e0d3d44..39f83c0 100644 --- a/test/affine_tests.jl +++ b/test/affine_tests.jl @@ -41,10 +41,11 @@ R = 4 T = Quantics.affine_transform_matrix(R, A, b) + M, N = size(A) mpo = Quantics.affine_transform_mpo( - outsite[1:R, 1:3], insite[1:R, 1:3], A, b) + outsite[1:R, 1:M], insite[1:R, 1:N], A, b) Trec = Quantics.affine_mpo_to_matrix( - outsite[1:R, 1:3], insite[1:R, 1:3], mpo) + outsite[1:R, 1:M], insite[1:R, 1:N], mpo) @test T == Trec end @@ -54,10 +55,11 @@ R = 4 T = Quantics.affine_transform_matrix(R, A, b) + M, N = size(A) mpo = Quantics.affine_transform_mpo( - outsite[1:R, 1:2], insite[1:R, 1:3], A, b) + outsite[1:R, 1:M], insite[1:R, 1:N], A, b) Trec = Quantics.affine_mpo_to_matrix( - outsite[1:R, 1:2], insite[1:R, 1:3], mpo) + outsite[1:R, 1:M], insite[1:R, 1:N], mpo) @test T == Trec end @@ -68,10 +70,28 @@ for R in [2, 3, 6] for bc in boundaries T = Quantics.affine_transform_matrix(R, A, b, bc) + M, N = size(A) mpo = Quantics.affine_transform_mpo( - outsite[1:R, 1:1], insite[1:R, 1:1], A, b, bc) + outsite[1:R, 1:M], insite[1:R, 1:N], A, b, bc) Trec = Quantics.affine_mpo_to_matrix( - outsite[1:R, 1:1], insite[1:R, 1:1], mpo) + outsite[1:R, 1:M], insite[1:R, 1:N], mpo) + @test T == Trec + end + end + end + + @testset "compare_denom_even" begin + A = reshape([1//2], 1, 1) + b = [3] # b = 5 would not work :( + + for R in [3, 5] + for bc in boundaries + T = Quantics.affine_transform_matrix(R, A, b, bc) + M, N = size(A) + mpo = Quantics.affine_transform_mpo( + outsite[1:R, 1:M], insite[1:R, 1:N], A, b, bc) + Trec = Quantics.affine_mpo_to_matrix( + outsite[1:R, 1:M], insite[1:R, 1:N], mpo) @test T == Trec end end From 4ff56cad6ae195ec5937a39ceedc4d85e9caa5a8 Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Fri, 17 Jan 2025 17:09:37 +0100 Subject: [PATCH 22/24] Fix bug for even denominator - now we can do light cone coordinates. Immutability is a blessing - there was a spurious accumulation of the carry. With this change, we are now able to do the important transformation (x,x') -> 1/2 * (x + x', x - x') with periodic boundary conditions and without. --- src/affine.jl | 6 ++---- test/affine_tests.jl | 19 ++++++++++++++++++- 2 files changed, 20 insertions(+), 5 deletions(-) diff --git a/src/affine.jl b/src/affine.jl index 199d1a4..78eb903 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -185,8 +185,7 @@ function affine_transform_core( y_index = digits_to_number(y) + 1 # Correct z and compute carry - z = @. z - s * y - d = z .>> 1 + d = @. (z - s * y) .>> 1 # Store this d_mat = get!(out, d) do @@ -206,8 +205,7 @@ function affine_transform_core( # right once in the driver routine (affine_transform_tensors). for (y_index, y) in enumerate(all_y) # Correct z and compute carry - z = @. z - s * y - d = z .>> 1 + d = @. (z - s * y) >> 1 # Store this d_mat = get!(out, d) do diff --git a/test/affine_tests.jl b/test/affine_tests.jl index 39f83c0..1c39372 100644 --- a/test/affine_tests.jl +++ b/test/affine_tests.jl @@ -82,7 +82,7 @@ @testset "compare_denom_even" begin A = reshape([1//2], 1, 1) - b = [3] # b = 5 would not work :( + b = [3] # XXX b = 5 would not work :( for R in [3, 5] for bc in boundaries @@ -97,4 +97,21 @@ end end + @testset "compare_light_cone" begin + A = 1//2 * [1 1; 1 -1] + b = [2; 3] + + for R in [3, 4] + for bc in boundaries + T = Quantics.affine_transform_matrix(R, A, b, bc) + M, N = size(A) + mpo = Quantics.affine_transform_mpo( + outsite[1:R, 1:M], insite[1:R, 1:N], A, b, bc) + Trec = Quantics.affine_mpo_to_matrix( + outsite[1:R, 1:M], insite[1:R, 1:N], mpo) + @test T == Trec + end + end + end + end From 781110425f3d3bf5e69d64aa664d7c6584055f1a Mon Sep 17 00:00:00 2001 From: Markus Wallerberger Date: Fri, 17 Jan 2025 17:25:39 +0100 Subject: [PATCH 23/24] Generalize affine_transform_mpo for the cases R=0 and R=1. --- src/affine.jl | 21 +++++++++++++-------- test/affine_tests.jl | 2 +- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/src/affine.jl b/src/affine.jl index 78eb903..5830042 100644 --- a/src/affine.jl +++ b/src/affine.jl @@ -70,15 +70,20 @@ function affine_transform_mpo( # Fill the MPO, taking care to not include auxiliary links at the edges mpo = MPO(R) spin_dims = ntuple(_ -> 2, M + N) - mpo[1] = ITensor(reshape(tensors[1], size(tensors[1], 2), spin_dims...), - (link[1], y[1,:]..., x[1,:]...)) - for r in 2:R-1 - newshape = size(tensors[r])[1:2]..., spin_dims... - mpo[r] = ITensor(reshape(tensors[r], newshape), - (link[r-1], link[r], y[r,:]..., x[r,:]...)) + if R == 1 + mpo[1] = ITensor(reshape(tensors[1], spin_dims...), + (y[1,:]..., x[1,:]...)) + elseif R > 1 + mpo[1] = ITensor(reshape(tensors[1], size(tensors[1], 2), spin_dims...), + (link[1], y[1,:]..., x[1,:]...)) + for r in 2:R-1 + newshape = size(tensors[r])[1:2]..., spin_dims... + mpo[r] = ITensor(reshape(tensors[r], newshape), + (link[r-1], link[r], y[r,:]..., x[r,:]...)) + end + mpo[R] = ITensor(reshape(tensors[R], size(tensors[R], 1), spin_dims...), + (link[R-1], y[R,:]..., x[R,:]...)) end - mpo[R] = ITensor(reshape(tensors[R], size(tensors[R], 1), spin_dims...), - (link[R-1], y[R,:]..., x[R,:]...)) return mpo end diff --git a/test/affine_tests.jl b/test/affine_tests.jl index 1c39372..5848e1b 100644 --- a/test/affine_tests.jl +++ b/test/affine_tests.jl @@ -67,7 +67,7 @@ A = reshape([1//3], 1, 1) b = [0] - for R in [2, 3, 6] + for R in [1, 3, 6] for bc in boundaries T = Quantics.affine_transform_matrix(R, A, b, bc) M, N = size(A) From 0fa9a6acdce2d4519402dce13cd7d2cda19686ee Mon Sep 17 00:00:00 2001 From: Hiroshi Shinaoka Date: Tue, 21 Jan 2025 15:21:01 +0800 Subject: [PATCH 24/24] Add more tests on affine_transform_matrix() --- test/affine_tests.jl | 72 ++++++++++++++++++++++++++++++++++---------- 1 file changed, 56 insertions(+), 16 deletions(-) diff --git a/test/affine_tests.jl b/test/affine_tests.jl index 5848e1b..5a0f78a 100644 --- a/test/affine_tests.jl +++ b/test/affine_tests.jl @@ -4,6 +4,39 @@ using Quantics using LinearAlgebra + # Test results of affine_transform_matrix() + function test_affine_transform_matrix_multi_variables(R, A, b, T, boundary) + M, N = size(A) + + yranges = [0:(2^R - 1) for _ in 1:M] + xranges = [0:(2^R - 1) for _ in 1:N] + + # Iterate over all combinations of input and output variables + for x_vals in Iterators.product(xranges...), y_vals in Iterators.product(yranges...) + x = collect(x_vals) # Convert input tuple to a vector + Axb = A * x .+ b # Apply the affine transformation + if boundary == Quantics.PeriodicBoundaryConditions() + Axb = mod.(Axb, 2^R) + elseif boundary == Quantics.OpenBoundaryConditions() + Axb = map(x -> 0 <= x <= 2^R - 1 ? x : nothing, Axb) + end + ref = Axb == collect(y_vals) # Compare to the reference output + + # Calculate indices for input and output + iy = 1 + sum(y_vals[i] * (2^R)^(i - 1) for i in 1:M) + ix = 1 + sum(x_vals[i] * (2^R)^(i - 1) for i in 1:N) + + @test T[iy, ix] == ref # Verify the transformation result + end + end + + testtests = Dict( + (1, 1) => [(reshape([1], 1, 1), [1])], + (1, 2) => [([1 0], [0]), ([2 -1], [1])], + (2, 1) => [([1; 0], [0, 0]), ([2; -1], [1, -1])], + (2, 2) => [([1 0; 1 1], [0; 1]), ([2 0; 4 1], [100; -1])] + ) + vars = ("x", "y", "z") insite = [Index(2; tags="$v$l") for l in 1:10, v in vars] outsite = [Index(2; tags="$v$l")' for l in 1:10, v in vars] @@ -22,6 +55,14 @@ @test T * T' == I end + @testset "full R=$R, boundary=$(boundary), M=$M, N=$N" for R in [1, 2], boundary in boundaries, M in [1, 2], N in [1, 2] + for (A, b) in testtests[(M, N)] + A_ = reshape(A, M, N) + b_ = reshape(b, M) + T = Quantics.affine_transform_matrix(R, A_, b_, boundary) + test_affine_transform_matrix_multi_variables(R, A_, b_, T, boundary) + end + end @testset "compare_simple" begin A = [1 0; 1 1] b = [0; 0] @@ -29,9 +70,9 @@ T = Quantics.affine_transform_matrix(R, A, b) mpo = Quantics.affine_transform_mpo( - outsite[1:R, 1:2], insite[1:R, 1:2], A, b) + outsite[1:R, 1:2], insite[1:R, 1:2], A, b) Trec = Quantics.affine_mpo_to_matrix( - outsite[1:R, 1:2], insite[1:R, 1:2], mpo) + outsite[1:R, 1:2], insite[1:R, 1:2], mpo) @test T == Trec end @@ -43,9 +84,9 @@ T = Quantics.affine_transform_matrix(R, A, b) M, N = size(A) mpo = Quantics.affine_transform_mpo( - outsite[1:R, 1:M], insite[1:R, 1:N], A, b) + outsite[1:R, 1:M], insite[1:R, 1:N], A, b) Trec = Quantics.affine_mpo_to_matrix( - outsite[1:R, 1:M], insite[1:R, 1:N], mpo) + outsite[1:R, 1:M], insite[1:R, 1:N], mpo) @test T == Trec end @@ -57,14 +98,14 @@ T = Quantics.affine_transform_matrix(R, A, b) M, N = size(A) mpo = Quantics.affine_transform_mpo( - outsite[1:R, 1:M], insite[1:R, 1:N], A, b) + outsite[1:R, 1:M], insite[1:R, 1:N], A, b) Trec = Quantics.affine_mpo_to_matrix( - outsite[1:R, 1:M], insite[1:R, 1:N], mpo) + outsite[1:R, 1:M], insite[1:R, 1:N], mpo) @test T == Trec end @testset "compare_denom_odd" begin - A = reshape([1//3], 1, 1) + A = reshape([1 // 3], 1, 1) b = [0] for R in [1, 3, 6] @@ -72,16 +113,16 @@ T = Quantics.affine_transform_matrix(R, A, b, bc) M, N = size(A) mpo = Quantics.affine_transform_mpo( - outsite[1:R, 1:M], insite[1:R, 1:N], A, b, bc) + outsite[1:R, 1:M], insite[1:R, 1:N], A, b, bc) Trec = Quantics.affine_mpo_to_matrix( - outsite[1:R, 1:M], insite[1:R, 1:N], mpo) + outsite[1:R, 1:M], insite[1:R, 1:N], mpo) @test T == Trec end end end @testset "compare_denom_even" begin - A = reshape([1//2], 1, 1) + A = reshape([1 // 2], 1, 1) b = [3] # XXX b = 5 would not work :( for R in [3, 5] @@ -89,16 +130,16 @@ T = Quantics.affine_transform_matrix(R, A, b, bc) M, N = size(A) mpo = Quantics.affine_transform_mpo( - outsite[1:R, 1:M], insite[1:R, 1:N], A, b, bc) + outsite[1:R, 1:M], insite[1:R, 1:N], A, b, bc) Trec = Quantics.affine_mpo_to_matrix( - outsite[1:R, 1:M], insite[1:R, 1:N], mpo) + outsite[1:R, 1:M], insite[1:R, 1:N], mpo) @test T == Trec end end end @testset "compare_light_cone" begin - A = 1//2 * [1 1; 1 -1] + A = 1 // 2 * [1 1; 1 -1] b = [2; 3] for R in [3, 4] @@ -106,12 +147,11 @@ T = Quantics.affine_transform_matrix(R, A, b, bc) M, N = size(A) mpo = Quantics.affine_transform_mpo( - outsite[1:R, 1:M], insite[1:R, 1:N], A, b, bc) + outsite[1:R, 1:M], insite[1:R, 1:N], A, b, bc) Trec = Quantics.affine_mpo_to_matrix( - outsite[1:R, 1:M], insite[1:R, 1:N], mpo) + outsite[1:R, 1:M], insite[1:R, 1:N], mpo) @test T == Trec end end end - end