General affine transforms

[mwallerb]

The PR implements QTTs for arbitrary affine transforms y = A*x +b. In particular, it constructs (fused-representation) MPOs for the boolean matrices:

T[y..., x...] = iszero(y - A * x - b)                 (OpenBoundaryConditions)
T[y..., x...] = iszero((y - A * x - b) .% 2^R)        (PeriodicBoundaryConditions)

where A is an M x N matrix b is an N vector, with arbitrary rational coefficients.

The core function added is affine_transform_mpo, which returns an ITensor MPO for a given affine transform.

For PeriodicBoundaryConditions( ), I have not encountered any problems. For OpenBoundaryConditions, it also works almost all of the time. The only thing case that I found does not work for now is when b*s > 2^R.

[mwallerb]

Does this code support cases where the common denominator is not a power of 2?

Yes, it should handle these cases just fine. The test compare_denom_odd has a denominator of 3, and it works. I haven't exhaustively tested a lot of denominators, but for periodic boundary conditions I am fairly confident it should do the right thing.

For open boundary conditions, there is a strange issue that if the denominator is not one and the shift is very large that you get wrong results; I think it has to do with the way we handle the outgoing carry of the first tensor.

[shinaoka]

OK, let me examine your code. BTW, should we implement a linear interpolation in another PR?
For instance, y = x/3, x/3 is not an integer for some x.
If we want to shrink the support of the function from [0, 1] to [0, 1/3], we need a linear interpolation or mapping to the nearest neighbor integer.

[mwallerb]

BTW, should we implement a linear interpolation in another PR? For instance, y = x/3, x/3 is not an integer for some x.

I agree that this is needed. I also agree that this should be a separate PR, this code is complicated enough as it is.

In principle, once this code is merged, one can use the sum of several MPOs to emulate this. For example, one could do three MPOs: y = x/3; y = (x+1)/3; y = (x-1)/3 and use a weighted sum of them.