diff --git a/PyTAT/PyTAT.hpp b/PyTAT/PyTAT.hpp
index 8f6d83f4a..b6ad7b1f6 100644
--- a/PyTAT/PyTAT.hpp
+++ b/PyTAT/PyTAT.hpp
@@ -536,14 +536,16 @@ namespace TAT {
                         const DefaultName& common_name_v,
                         const DefaultName& singular_name_u,
                         const DefaultName& singular_name_v,
-                        double cut) {
+                        int cut,
+                        double relative_cut,
+                        double temperature) {
                         Cut real_cut = NoCut();
-                        if (cut > 0) {
-                           if (cut >= 1) {
-                              real_cut = RemainCut(Size(cut));
-                           } else {
-                              real_cut = RelativeCut(cut);
-                           }
+                        if (temperature > 0) { // T = 0 => RemainCut
+                           real_cut = BoltzmannCut(temperature, cut, &random_engine);
+                        } else if (relative_cut > 0) {
+                           real_cut = RelativeCut(relative_cut);
+                        } else if (cut > 0) {
+                           real_cut = RemainCut(cut);
                         }
                         auto result = tensor.svd(free_name_set_u, common_name_u, common_name_v, singular_name_u, singular_name_v, real_cut);
                         return py::make_tuple(std::move(result.U), std::move(result.S), std::move(result.V));
@@ -553,7 +555,9 @@ namespace TAT {
                      py::arg("common_name_v"),
                      py::arg("singular_name_u"),
                      py::arg("singular_name_v"),
-                     py::arg("cut") = -1,
+                     py::arg("cut") = 0,
+                     py::arg("relative_cut") = 0,
+                     py::arg("temperature") = 0,
                      "Singular value decomposition")
                .def(
                      "qr",
diff --git a/include/TAT/implement/svd.hpp b/include/TAT/implement/svd.hpp
index 426e0a7e3..8fdb49a64 100644
--- a/include/TAT/implement/svd.hpp
+++ b/include/TAT/implement/svd.hpp
@@ -25,6 +25,7 @@
 
 #include "../structure/tensor.hpp"
 #include "../utility/allocator.hpp"
+#include "../utility/common_variable.hpp"
 #include "../utility/timer.hpp"
 
 extern "C" {
@@ -448,11 +449,12 @@ namespace TAT {
                      }
                   }
                }
-               if (maximum_singular != 0) {
+               if (maximum_singular == 0) {
                   // sometimes non zero singular number is less than cut_i
-                  remain_dimension_u.at(maximum_symmetry) += 1;
-                  remain_dimension_v.at(-maximum_symmetry) += 1;
+                  break;
                }
+               remain_dimension_u.at(maximum_symmetry) += 1;
+               remain_dimension_v.at(-maximum_symmetry) += 1;
             }
 
             // delete element of tensor S
@@ -501,6 +503,70 @@ namespace TAT {
                ++it;
             }
          }
+      } else if (auto cut_value = std::get_if<BoltzmannCut>(&cut)) {
+         Size cut_i = cut_value->value;
+         if (cut_i < total_dimension) {
+            double temperature = cut_value->temperature;
+            auto engine = cut_value->engine; // this is a pointer
+            auto dist = std::uniform_real_distribution<real_scalar<ScalarType>>(0, 1);
+            for (const auto& [symmetry, vector_s] : result_s) {
+               remain_dimension_u[symmetry] = 0;
+               remain_dimension_v[-symmetry] = 0;
+            }
+            for (Size i = 0; i < cut_i; i++) {
+               // s_i' = s_i / sum of s     # for scaling invariance
+               // E_i = - ln s_i'           # s=0 should have p=0, since we can always add any number of zero singular
+               // p_i ~ exp( - E / T )      # by definition
+               // => p'_i = s_i' ^ (1 / T)
+               // => p_i = p'_i / sum of p'
+               real_scalar<ScalarType> sum_of_s = 0;
+               for (const auto& [symmetry, vector_s] : result_s) {
+                  if (auto& this_remain = remain_dimension_u.at(symmetry); this_remain != vector_s.size()) {
+                     auto this_s = vector_s[this_remain];
+                     sum_of_s += this_s;
+                  }
+               }
+               if (sum_of_s == 0) {
+                  // sometimes non zero singular number is less than cut_i
+                  break;
+               }
+               real_scalar<ScalarType> sum_of_p = 0;
+               for (const auto& [symmetry, vector_s] : result_s) {
+                  if (auto& this_remain = remain_dimension_u.at(symmetry); this_remain != vector_s.size()) {
+                     auto this_s = vector_s[this_remain];
+                     auto this_p = std::pow(this_s / sum_of_s, 1 / temperature);
+                     sum_of_p += this_p;
+                  }
+               }
+               real_scalar<ScalarType> random_number = dist(*engine) * sum_of_p;
+               real_scalar<ScalarType> resum_of_p = 0;
+               for (const auto& [symmetry, vector_s] : result_s) {
+                  if (auto& this_remain = remain_dimension_u.at(symmetry); this_remain != vector_s.size()) {
+                     auto this_s = vector_s[this_remain];
+                     auto this_p = std::pow(this_s / sum_of_s, 1 / temperature);
+                     resum_of_p += this_p;
+                     if (resum_of_p >= random_number) {
+                        // choose this symmetry
+                        remain_dimension_u[symmetry]++;
+                        remain_dimension_v[-symmetry]++;
+                        break;
+                     }
+                  }
+               }
+            }
+
+            // delete element of tensor S
+            for (auto it = result_s.begin(); it != result_s.end();) {
+               const auto& symmetry = it->first;
+               const auto& this_remain = remain_dimension_u.at(symmetry);
+               if (this_remain == 0) {
+                  it = result_s.erase(it);
+               } else {
+                  it->second.resize(this_remain);
+                  ++it;
+               }
+            }
+         }
       }
       // cut analyze done
 
diff --git a/include/TAT/structure/tensor.hpp b/include/TAT/structure/tensor.hpp
index 3d6bac347..fcd9ec9c0 100644
--- a/include/TAT/structure/tensor.hpp
+++ b/include/TAT/structure/tensor.hpp
@@ -25,6 +25,7 @@
 #include <array>
 #include <map>
 #include <memory>
+#include <random>
 #include <set>
 #include <tuple>
 #include <variant>
@@ -37,6 +38,7 @@
 #include "symmetry.hpp"
 
 namespace TAT {
+   struct NoCut {};
    struct RemainCut {
       Size value;
       explicit RemainCut(Size v) : value(v) {}
@@ -45,13 +47,18 @@ namespace TAT {
       double value;
       explicit RelativeCut(double v) : value(v) {}
    };
-   struct NoCut {};
+   struct BoltzmannCut {
+      double temperature;
+      Size value;
+      std::default_random_engine* engine;
+      BoltzmannCut(double t, Size v, std::default_random_engine* e) : temperature(t), value(v), engine(e) {}
+   };
    /**
     * Used to describle how to cut when doing svd to a tensor
     *
-    * Is one of RemainCut, RelativeCut and NoCut
+    * Is one of NoCut, RemainCut, RelativeCut and BoltzmannCut.
     */
-   using Cut = std::variant<RemainCut, RelativeCut, NoCut>;
+   using Cut = std::variant<NoCut, RemainCut, RelativeCut, BoltzmannCut>;
 
    template<typename ScalarType = double, typename Symmetry = Symmetry<>, typename Name = DefaultName>
    struct TensorShape;