Commit c66920f5 authored by Ben Huber's avatar Ben Huber

changed calculate_svd to make a relative error of at most _eps (closes #215)

parent 174b5ae3
Pipeline #1087 failed with stages
in 5 minutes and 12 seconds
......@@ -73,33 +73,19 @@ static misc::UnitTest tt_round("TT", "TTTensor_Rounding", [](){
Tensor A4 = Tensor::random({2,2,2,2,2,2,2,2});
Tensor B4;
TTTensor TTA4(A4, 1e-14);
B4(j,i^7) = TTA4(j,i^7);
TEST(approx_equal(B4,A4, 1e-14));
TTA4.round(1e-14);
B4(j,i^7) = TTA4(j,i^7);
TEST(approx_equal(B4,A4, 1e-14));
TTA4.round(512);
B4(j,i^7) = TTA4(j,i^7);
TEST(approx_equal(B4,A4, 1e-14));
for (double eps = 1e-3; eps < 1.0; eps*=1.2) {
TTTensor B4(A4);
B4.round(eps);
MTEST(approx_equal(A4, B4, eps), eps << " vs " << frob_norm(A4-Tensor(B4))/frob_norm(A4));
}
Tensor A5 = Tensor::random({5,6,3,1,4,2,8,1});
Tensor B5;
TTTensor TTA5(A5, 1e-14);
B5(j,i^7) = TTA5(j,i^7);
TEST(approx_equal(B5,A5, 1e-14));
TTA5.round(1e-14);
B5(j,i^7) = TTA5(j,i^7);
TEST(approx_equal(B5,A5, 1e-14));
TTA5.round(576);
B5(j,i^7) = TTA5(j,i^7);
TEST(approx_equal(B5,A5, 1e-14));
for (double eps = 1e-3; eps < 1.0; eps*=1.2) {
TTTensor B5(A5);
B5.round(eps);
MTEST(approx_equal(A5, B5, eps), eps << " vs " << frob_norm(A5-Tensor(B5))/frob_norm(A5));
}
});
......
......@@ -1469,11 +1469,12 @@ namespace xerus {
}
// Find rank due to the Epsilon (NOTE the scaling factor can be ignored, as it does not change the ratios).
for(size_t j = 1; j < rank; ++j) {
if (tmpS[j] <= _eps*tmpS[0]) {
rank = j;
break;
}
// For the total error to be < _eps, the sum of discarded singular value squares must be smaller than _eps times the norm, squared
value_t maxErrorSqr = misc::sqr(blasWrapper::two_norm(tmpS.get(), rank) * _eps);
value_t error = 0;
while (rank > 1 && error + misc::sqr(tmpS[rank-1]) <= maxErrorSqr) {
error += misc::sqr(tmpS[rank-1]);
rank -= 1;
}
// Create tensor from diagonal values
......
......@@ -662,9 +662,11 @@ namespace xerus {
const size_t initialCorePosition = corePosition;
canonicalize_right();
auto epsPerSite = _eps / std::sqrt((double)numComponents-1);
for(size_t i = 0; i+1 < numComponents; ++i) {
round_edge(numComponents-i, numComponents-i-1, _maxRanks[numComponents-i-2], _eps, 0.0);
round_edge(numComponents-i, numComponents-i-1, _maxRanks[numComponents-i-2], epsPerSite, 0.0);
}
assume_core_position(0);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment