Commit aaa6102c authored by Ben Huber's avatar Ben Huber

non-operator tt pseudo inverse as unittest

parent fe6c4a87
Pipeline #631 passed with stages
in 8 minutes and 8 seconds
......@@ -4,8 +4,8 @@
using namespace xerus;
static misc::UnitTest tt_entryprod("TT", "entrywise_product", [](){
Index i,j,k;
Index i,j,k;
TTTensor A = TTTensor::random(std::vector<size_t>(10,2), std::vector<size_t>(9,2));
TTTensor B = TTTensor::random(std::vector<size_t>(10,2), std::vector<size_t>(9,2));
......@@ -34,8 +34,8 @@ static misc::UnitTest tt_entryprod("TT", "entrywise_product", [](){
});
static misc::UnitTest tt_soft("TT", "soft_thresholding", [](){
Index i,j,k;
Index i,j,k;
TTTensor A = TTTensor::random(std::vector<size_t>(10,2), std::vector<size_t>(9,2));
TTTensor B = TTTensor::random(std::vector<size_t>(10,2), std::vector<size_t>(9,2));
......@@ -53,3 +53,30 @@ static misc::UnitTest tt_soft("TT", "soft_thresholding", [](){
TEST(frob_norm(Cf - Tensor(Co))/frob_norm(Cf) < 1e-14);
TEST(frob_norm(Cf - Tensor(C))/frob_norm(Cf) < 1e-14);
});
static misc::UnitTest tt_pseudo_inv("TT", "Non-operator Pseudo Inverse", [](){
Index i,j,k,r1,r2,r3,r4;
const size_t d = 2;
const TTTensor op = TTTensor::random(std::vector<size_t>(2*d, 10), std::vector<size_t>(2*d-1, 4));
TTTensor tmp = op;
tmp.move_core(d);
auto parts = tmp.chop(d);
Tensor U,S,V;
(U(i,r1), S(r1,r2), V(r2,j^2)) = SVD(tmp.get_component(d)(i,j^2));
S.modify_diagonal_entries([](double &_v){
if (_v>1e-10) {
_v = 1/_v;
}
});
TensorNetwork pInv;
pInv(j, i^(d-1), k^d) = parts.first(k^d,r1) * U(r1,r2) * S(r2,r3) * V(r3,j,r4) * parts.second(r4,i^(d-1));
double error = frob_norm(pInv(i^d,r1^d)*op(r1^d,r2^d)*pInv(r2^d,j^d) - pInv(i^d,j^d));
MTEST(error < 1e-10, "A^+ A A^+ != A^+ ? " << error);
error = frob_norm(op(i^d,r1^d)*pInv(r1^d,r2^d)*op(r2^d,j^d) - op(i^d,j^d));
MTEST(error < 1e-10, "A A^+ A != A ? " << error);
});
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