issue and proposed fix for create_deltas

[eoin-dp-oneill]

The issue:

As the function header points out, the create_deltas function in utils.py is inefficient. This function adds a fourth leg to a 3 legged MPO by creating a new empty tensor of the desired output rank and dimension and populating the elements accordingly. For bond dims. of ~100 it takes ~3 times the time required to do the contractions in the gradient calculation.

Proposed fix:

Paul E. has written a new function create_deltas_noscramble. Rather than creating a new tensor, this increases the rank of the original tensor by 1 and populates the remaining elements accordingly. The code is as follows:

def create_delta_noscramble(
        tensor: ndarray):
    tensor_shape = tensor.shape
    ret_shape=tensor.shape+(tensor.shape[-1],)
    ret_ndarray=np.zeros(ret_shape,dtype=tensor.dtype)
    for a in range(ret_shape[-1]):
        ret_ndarray[:,:,a,a]=tensor[:,:,a]
    return ret_ndarray

The result is a computation time 10 times less than the original code. Below is a demonstration of the scaling of the computation time with bond dimension for the two versions.

[paulreastham]

I did also think about whether we could implement propagation using the 3-legged MPOs. That would speed the propagation up by another significant factor, although probably not quite as big an improvement as this. However, I think that would be specialized to the case of diagonal coupling operators (where transform_in and transform_out are None).

[gefux]

Hi @eoin-dp-oneill and @paulreastham. Thank you for spotting this and making some performence tests!

I think that the proposed solution goes into the right direction.
The general function must be able to create arbitrary deltas because that is needed on varous places in the code (other than the PT & system contraction that you are thinking of). The parameter index_scrabling is for that. The above code works for the cases in which `index_scrambling=[0,1,...,n-1,n,n]' only.
Of couse one can create special functions for special cases that are called by the general function, but we need to replace the gerneral function.

I wonder what a good solution to this could be. I guess that the biggest difference are those cases where the extension can be made inplace and those where it cannot. This would suggest that there should be two functions, one fast one which does an inplace replacement and works only for the restricted set of scramblings, and one that is more general that returns a new tensor (possibly using a copy of the fast inplace function). What do you think?

@paulreastham ad 3-legged MPOs: I think that this could be done and could increase performance significantly. As this happens on a different level of abstraction, I'd suggest that this would be a speparate Issue to discuss (feel free to open one). I briefly comment here anyways: It is currently possible to get the process tensor mpo without transformations, and it is possible that the MPO is internally stored as an MPS (which is a special case exploited by PT-TEMPO). We could adopt the BaseProcessTensor interface to have an option to return just the three legged MPS if it is available. Then the transforms could be contracted into the system propergators and the process tensor contraction made more efficient.