Realization of SVD compression for MPOΒΆ

From the construction of the SVD based compression scheme for MPO, we try to understand the convention used for the tensor library. This could serve as an example on how you develop your own MPO/MPS operation based on the tensor library.

In mpo_impl.h, the class mpo_impl contains the following functions

1
2
3
4
5
 std::pair<size_t, magnitude_type> svdCompressLeft(int maxbonddimension = -1, magnitude_type tol = 1.0e-12, int verbose = 0);

 std::pair<size_t, magnitude_type> svdCompressRight(int maxbonddimension = -1, magnitude_type tol = 1.0e-12, int verbose = 0);

 std::pair<size_t, magnitude_type> svdCompress(int maxbonddimension=1000000, magnitude_type tol = 1.0e-12, int verbose = 0);

Function svdCompress would call svdCompressLeft and svdCompressRight.

1
2
3
4
5
6
7
 template <typename _Tensor4>
 std::pair<size_t, typename mpo_impl<_Tensor4>::magnitude_type>
 mpo_impl<_Tensor4>::svdCompress(int maxbonddimension, magnitude_type tol, int verbose)
 {
     svdCompressLeft(-1, tol, verbose);
     return svdCompressRight(maxbonddimension, tol, verbose);
 }

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
 template <typename _Tensor4>
 std::pair<size_t, typename mpo_impl<_Tensor4>::magnitude_type>
 mpo_impl<_Tensor4>::svdCompressLeft(int maxbonddimension, magnitude_type tol, int verbose)
 {
     size_t L = super_type::size(), i;
     assert(L > 1);
     typedef std::pair<size_t, magnitude_type> bonderror_type;
     std::vector<bonderror_type> bonderror;
     std::array<size_t, 1> Is2 = {2}, Is1 = {1}, Is0 = {0};
     std::array<size_t, 4> per1 = {0, 1, 3, 2}, per2 = {1, 0, 2, 3};
     for (i = 0; i < L - 1; ++i) {
         if (verbose >= 2) {
             std::cout << "Compress mpo by balanced svd from left to right on site: " << i << '\n';
         }
         auto svdResults = gc::svdCompress(operator[](i), Is2, maxbonddimension, tol, verbose);
         bonderror.push_back(std::get<3>(svdResults));
         if (std::get<1>(svdResults).size() == 0) {
             if (verbose) {
                 std::cout << "mpo becomes zero after cutoff.\n";
             }
             for (size_t j = 0; j < L; ++j) {
                 operator[](j).clear();
             }
             break;
         }
         typename std::remove_reference<decltype(std::get<2>(svdResults))>::type sm;
         diag(std::get<1>(svdResults), sm);
         std::get<2>(svdResults) = contract(sm, Is1, std::get<2>(svdResults), Is0);
         // operator[](i) = contract(std::get<0>(svdResults), Is3, sm, Is0);
         operator[](i) = gc::permute(std::get<0>(svdResults), per1);
         operator[](i + 1) = contract(std::get<2>(svdResults), Is1, operator[](i + 1), Is1);
         operator[](i + 1) = gc::permute(operator[](i + 1), per2);
     }
     bonderror_type bet = bonderror[0];
     for (i = 1; i < bonderror.size(); ++i) {
         bet = std::make_pair(std::max(bet.first, bonderror[i].first), bet.second + bonderror[i].second);
     }
     return bet;
 }