[WIP] [experimental] mps simulator backend

[Qottmann]

check out the demo.ipynb (or for always the latest version go "Files changed" => "demo.ipynb" => "View file")

Rough prototype, likely wont become a PennyLane feature.
A simple MPS simulator to showcase this approach's benefits.

• Ease of use: user does not need to think about optimal contraction strategies
• Well established: can rely on 20+ years of algorithms that still form the state of the art in many areas of computational quantum physics
• Allow for approximations (unlike full contraction simulators)
• Single parameter to control the degree of the approximation (bond dimension, can be used like a black box)
• Polynomial scaling in the system size ($L \chi^2 2$ coefficients to store)
  
• Faithfully simulate shallow 50 qubits circuits in ~2 minutes on laptop CPUs (no HPC knowledge necessary)
  
  [ $D$ depth of circuit, $\chi$ bond dimension, increased to see if the result saturates (as no exact result is known for such a large system size!). For $D=13$ no saturation yet, so one would need to crank up bond dimension a little further.]

The purpose of this PR is exploration and prototyping. I.e. this is a brute force prototype written in basic numpy, so dont expect clean nor optimized code. Though most theory is already well established, there are still some choices to be made. In particular how to treat long-range interactions. My current preference is writing long-distance gates as MPOs (see #3583 (comment) as well as iTensor forum discussion)

In particular things that may be improved:

• Pre-compute / hardcode 2-site decomposition of standard 2-site gates like CNOT
• Truncate MPO bond dimension (currently fixed to its max, 4)
• Combine MPOs into one bigger MPO (instead of consecutively applying MPO-MPS contractions and restoring MPS) Needs exploration if actually better
• Any room for optimizing all the seemingly pointless identity contractions in the long-range cases?
• Allow arbitrary site gates (minor improvement to current implementation, just wasnt priority so far)
• Specifically for CNOT, there might be some clever ways of just permuting some of the indices? (needs exploration)
• Expval capabilities still limited to Hamiltonians basically
• Overall code quality, thorough testing, integrating with the latest device API (this is based on a veeery old prototype)

[github-actions]

Hello. You may have forgotten to update the changelog!
Please edit doc/releases/changelog-dev.md with:

• A one-to-two sentence description of the change. You may include a small working example for new features.
• A link back to this PR.
• Your name (or GitHub username) in the contributors section.

[Qottmann]

I'm a bit confused here, do I need to do any SVD for single site updates? In single-site DMRG yes, but here this is different. Maybe @chaeyeunpark knows about this?

(tests with full state vector simulator agree atm)

[chaeyeunpark]

I will have a look in the new year. I also want to tag @ikurecic here.

[Qottmann]

I know now it is wrong, curious to know how to do it for single site (atm I am just forcing a two-site update)

[Qottmann]

This might be annoying for people to test, I needed it for some examples that wouldnt converge otherwise

[chaeyeunpark]

Few suggestions here. If this SVD algorithm (in tenpy) is not difficult, we can make our own version. Otherwise, we should make tenpy an optional dependency. It is still strange that svd in scipy is not robust enough.

[Qottmann]

Latest update lets you construct an MPO operator for a long distance 2-site operators (e.g. CNOT(1, 10)) (generalization to multi-site operators straight forward) and apply MPO-MPS products to restore an MPS (no truncation, so the new bond dimension is $4\chi$ since the MPO has bond dim 4). New update (see below) now with contraction and explicit testing

[Qottmann]

Integration of MPO-MPS contraction complete.

Added some basic explicit testing for my own sanity, see pennylane/devices/tests/test_mpssim.py

capabilities/restrictions are currently:

• Arbitrary single and two-site operations
• Only local expectation values of at most two sites

The latter point is very limiting at the moment, but I dont think I will manage to implement something more general here. A workaround and what I am currently doing is outputting the mps via measurement = qml.state() and then performing the mps methods SimpleMPS.site_expectation_value() and SimpleMPS.bond_expectation_value()

Edit: Support expval of Hamiltonians now, problem is that actually this old version of the default numpy simulator is not supporting Hamiltonian expvals.