Skip to content

GitLab

Commit 83e14fc2 authored by Sebastian Wolf's avatar Sebastian Wolf
Browse files

xALS working

parent 8303563d
Pipeline #702 failed with stages
in 4 minutes and 44 seconds
Hide whitespace changes
Inline Side-by-side
Showing with 170 additions and 113 deletions
include/xerus.h
View file @ 83e14fc2
......@@ -54,6 +54,7 @@
#include "xerus/performanceData.h"
#include "xerus/measurments.h"
#include "xerus/algorithms/als.h"
#include "xerus/algorithms/xals.h"
#include "xerus/algorithms/steepestDescent.h"
#include "xerus/algorithms/cg.h"
#include "xerus/algorithms/decompositionAls.h"
......
include/xerus/basic.h
View file @ 83e14fc2
......@@ -24,6 +24,7 @@
#pragma once
#include <stddef.h>
#include "misc/allocator.h"
#include "misc/standard.h"
#include "misc/namedLogger.h"
......
include/xerus/misc/standard.h
View file @ 83e14fc2
......@@ -24,6 +24,7 @@
#pragma once
#include <stddef.h>
#include <cstdint>
#include <cstddef>
......
src/unitTests/als.cxx
View file @ 83e14fc2
......@@ -119,7 +119,6 @@ static misc::UnitTest als_tut("ALS", "tutorial", [](){
xerus::ALS_SPD(A, X, B);
// LOG(asd, frob_norm(X-B));
TEST(frob_norm(X-B) < 1e-12);
......@@ -142,10 +141,20 @@ static misc::UnitTest als_tut("ALS", "tutorial", [](){
// ALSb.useResidualForEndCriterion = true;
// std::vector<value_t> perfdata;
ALS_SPD(A, X, C, 1e-12, pd);
LOG(woopR, frob_norm(A(i/2, j/2)*X(j&0) - C(i&0))/frob_norm(C));
xerus::xals(X, A, C);
TEST(frob_norm(A(i/2, j/2)*X(j&0) - C(i&0)) < 1e-4);
LOG(woop, frob_norm(A(i/2, j/2)*X(j&0) - C(i&0))/frob_norm(C));
X = xerus::TTTensor::random(stateDims, 2);
LOG(woopR, frob_norm(A(i/2, j/2)*X(j&0) - C(i&0))/frob_norm(C));
xerus::ALS_SPD(A, X, C, 1e-12, pd);
TEST(frob_norm(A(i/2, j/2)*X(j&0) - C(i&0)) < 1e-4);
LOG(woop, frob_norm(A(i/2, j/2)*X(j&0) - C(i&0))/frob_norm(C));
// perfdata.clear();
// ALSb(A, X, B, 1e-4, &perfdata);
// TEST(!misc::approx_equal(frob_norm(A(i^d, j^d)*X(j&0) - B(i&0)), 0., 1.));
......
src/xerus/algorithms/xals.cpp
View file @ 83e14fc2
......@@ -18,14 +18,14 @@
// or contact us at contact@libXerus.org.
/**
* @file
* @brief Implementation of the ADF variants.
*/
* @file
* @brief Implementation of the ADF variants.
*/
#include <xerus/algorithms/xals.h>
#include <xerus/basic.h>
#include <xerus/misc/internal.h>
#include <xerus/indexedTensorMoveable.h>
#include <xerus/misc/basicArraySupport.h>
......@@ -34,111 +34,127 @@
#endif
namespace xerus {
class InternalSolver {
const size_t d;
// const std::vector<Tensor> opComponents;
std::vector<Tensor> leftStack;
std::vector<Tensor> rightStack;
std::vector<Tensor> leftBStack;
std::vector<Tensor> rightBStack;
TTTensor& x;
const TTOperator& A;
const TTTensor& b;
Tensor leftStackA,rightStackA;
Tensor leftBStackB,rightBStackB;
public:
static const std::vector<Tensor> get_op_comps(const TTOperator& _A) {
std::vector<Tensor> comps;
for(size_t i = 0; i < _A.degree()/2; ++i) {
comps.push_back(_A.get_component(i));
// comps.back()
}
}
InternalSolver(TTTensor& _x, const TTOperator& _A, const TTTensor& _b) : d(_x.degree()),leftStack(d), rightStack(d), x(_x), A(_A), b(_b) {
}
void calc_left_stack(const size_t _position) {
// Tensor& xi = x.get_component(_position);
// Tensor& Ai = A.get_component(_position);
// Tensor& bi = b.get_component(_position);
//
// Tensor A0 = A.get_component(_position);
//
// if(_position == 0) {
// Tensor tmp;
// xi.reinterpret_dimensions({xi.dimensions[1], xi.dimensions[2]});
// Ai.reinterpret_dimensions({A0.dimensions[1], A0.dimensions[2], A0.dimensions[3]});
// bi.reinterpret_dimensions({bi.dimensions[1], bi.dimensions[2]});
// contract(tmp, Ai, true, xi, false, 1);
// contract(leftStack[_position], tmp, true, xi, false, 1);
//
// contract(leftBStack[_position], bi, xi, 1);
// } else {
// contract(leftStack[_position], leftStack[_position-1], x.component(_position), 1);
// }
}
void calc_right_stack(const size_t _position) {
}
double calc_residual_norm() {
Index i,j;
Tensor tmp;
tmp(i&0) = A(i/2, j/2)*x(j&0);
return frob_norm(tmp);
}
void solve() {
const double solutionsNorm = frob_norm(b);
std::vector<double> residuals(10, 1000.0);
const size_t maxIterations = 1000;
class InternalSolver {
const Index i1, i2, i3, j1 , j2, j3, k1, k2;
const size_t d;
std::vector<Tensor> leftAStack;
std::vector<Tensor> rightAStack;
std::vector<Tensor> leftBStack;
std::vector<Tensor> rightBStack;
TTTensor& x;
const TTOperator& A;
const TTTensor& b;
public:
InternalSolver(TTTensor& _x, const TTOperator& _A, const TTTensor& _b) : d(_x.degree()), leftAStack(d), rightAStack(d), leftBStack(d), rightBStack(d), x(_x), A(_A), b(_b) { }
void calc_left_stack(const size_t _position) {
Tensor xi = x.get_component(_position);
Tensor Ai = A.get_component(_position);
Tensor bi = b.get_component(_position);
for(size_t iteration = 0; maxIterations == 0 || iteration < maxIterations; ++iteration) {
x.move_core(0, true);
if(_position == 0) {
xi.reinterpret_dimensions({xi.dimensions[1], xi.dimensions[2]});
Ai.reinterpret_dimensions({Ai.dimensions[1], Ai.dimensions[2], Ai.dimensions[3]});
bi.reinterpret_dimensions({bi.dimensions[1], bi.dimensions[2]});
// Rebuild right stack
for(size_t corePosition = d-1; corePosition > 0; --corePosition) {
calc_right_stack(corePosition);
leftAStack[_position](i1, i2, i3) = xi(k1, i1)*Ai(k1, k2, i2)*xi(k2, i3);
leftBStack[_position](i1, i2) = xi(k1, i1)*bi(k1, i2);
} else {
leftAStack[_position](i1, i2, i3) = leftAStack[_position-1](j1, j2, j3)*xi(j1, k1, i1)*Ai(j2, k1, k2, i2)*xi(j3, k2, i3);
leftBStack[_position](i1, i2) = leftBStack[_position-1](j1, j2)*xi(j1, k1, i1)*bi(j2, k1, i2);
}
}
void calc_right_stack(const size_t _position) {
Tensor xi = x.get_component(_position);
Tensor Ai = A.get_component(_position);
Tensor bi = b.get_component(_position);
if(_position == d-1) {
xi.reinterpret_dimensions({xi.dimensions[0], xi.dimensions[1]});
Ai.reinterpret_dimensions({Ai.dimensions[0], Ai.dimensions[1], Ai.dimensions[2]});
bi.reinterpret_dimensions({bi.dimensions[0], bi.dimensions[1]});
rightAStack[_position](i1, i2, i3) = xi(i1, k1)*Ai(i2, k1, k2)*xi(i3, k2);
rightBStack[_position](i1, i2) = xi(i1, k1)*bi(i2, k1);
} else {
rightAStack[_position](i1, i2, i3) = xi(i1, k1, j1)*Ai(i2, k1, k2, j2)*xi(i3, k2, j3)*rightAStack[_position+1](j1, j2, j3);
rightBStack[_position](i1, i2) = xi(i1, k1, j1)*bi(i2, k1, j2)*rightBStack[_position+1](j1, j2);
}
}
double calc_residual_norm() {
const Index i, j;
Tensor tmp;
tmp(i&0) = A(i/2, j/2)*x(j&0)-b(i&0);
return frob_norm(tmp);
}
void solve() {
const double solutionsNorm = frob_norm(b);
std::vector<double> residuals(10, 1000.0);
const size_t maxIterations = 10000;
// Rebuild right stack
x.move_core(0, true);
for(size_t corePosition = d-1; corePosition > 0; --corePosition) {
calc_right_stack(corePosition);
}
for(size_t iteration = 0; maxIterations == 0 || iteration < maxIterations; ++iteration) {
// Calculate residual and check end condition
residuals.push_back(calc_residual_norm()/solutionsNorm);
if(residuals.back()/residuals[residuals.size()-10] > 0.99) {
LOG(xALS, "Residual decrease from " << std::scientific << residuals[10] << " to " << std::scientific << residuals.back() << " in " << residuals.size()-10 << " iterations.");
return; // We are done!
} else {
// LOG(xALS, "Residual decrease from " << std::scientific << residuals[10] << " to " << std::scientific << residuals.back() << " in " << residuals.size()-10 << " iterations.");
}
residuals.push_back(calc_residual_norm()/solutionsNorm);
if(residuals.back()/residuals[residuals.size()-10] > 0.99) {
LOG(ALS, "Residual decrease from " << std::scientific << residuals[10] << " to " << std::scientific << residuals.back() << " in " << residuals.size()-10 << " iterations.");
return; // We are done!
}
// Sweep Right -> Left
for(size_t corePosition = 0; corePosition < x.degree(); ++corePosition) {
// Actual Work to do
for(size_t corePosition = 0; corePosition < d; ++corePosition) {
Tensor op, rhs;
if(corePosition == 0) {
leftStackA = A.get_component(corePosition);
leftStackA.reinterpret_dimensions({leftStackA.dimensions[1], leftStackA.dimensions[2], leftStackA.dimensions[3]});
leftBStackB = b.get_component(corePosition);
leftBStackB.reinterpret_dimensions({leftBStackB.dimensions[1], leftBStackB.dimensions[2]});
} else {
contract(leftStackA, leftStack[corePosition-1], A.get_component(corePosition), 1);
contract(leftBStackB, leftBStack[corePosition-1], b.get_component(corePosition), 1);
}
contract(op, leftStackA, rightStack[corePosition+1], 1);
contract(rhs, leftBStackB, rightBStack[corePosition+1], 1);
solve_least_squares(x.component(corePosition), op, rhs, 0);
Tensor Ai = A.get_component(corePosition);
Tensor bi = b.get_component(corePosition);
if(corePosition == 0) {
Ai.reinterpret_dimensions({Ai.dimensions[1], Ai.dimensions[2], Ai.dimensions[3]});
bi.reinterpret_dimensions({bi.dimensions[1], bi.dimensions[2]});
op(i2, i3, j2, j3) = Ai(i2, j2, k2)*rightAStack[corePosition+1](i3, k2, j3);
rhs(i2, i3) = bi(i2, k2)*rightBStack[corePosition+1](i3, k2);
} else if(corePosition == d-1) {
Ai.reinterpret_dimensions({Ai.dimensions[0], Ai.dimensions[1], Ai.dimensions[2]});
bi.reinterpret_dimensions({bi.dimensions[0], bi.dimensions[1]});
op(i1, i2, j1, j2) = leftAStack[corePosition-1](i1, k1, j1)*Ai(k1, i2, j2);
rhs(i1, i2) = leftBStack[corePosition-1](i1, k1)*bi(k1, i2);
} else {
op(i1, i2, i3, j1, j2, j3) = leftAStack[corePosition-1](i1, k1, j1)*Ai(k1, i2, j2, k2)*rightAStack[corePosition+1](i3, k2, j3);
rhs(i1, i2, i3) = leftBStack[corePosition-1](i1, k1)*bi(k1, i2, k2)*rightBStack[corePosition+1](i3, k2);
}
solve_least_squares(x.component(corePosition), op, rhs, 0);
if(corePosition == 0) {
x.component(corePosition).reinterpret_dimensions({1, x.component(corePosition).dimensions[0], x.component(corePosition).dimensions[1]});
} else if(corePosition == d-1) {
x.component(corePosition).reinterpret_dimensions({x.component(corePosition).dimensions[0], x.component(corePosition).dimensions[1], 1});
}
// If we have not yet reached the end of the sweep we need to take care of the core and update our stacks
if(corePosition+1 < d) {
x.move_core(corePosition+1, true);
......@@ -146,10 +162,39 @@ namespace xerus {
}
}
// Sweep Right -> Left
for(size_t corePosition = d-1; corePosition > 0; --corePosition) {
//Actual Work to do
// Tensor op, rhs;
//
// Tensor Ai = A.get_component(corePosition);
// Tensor bi = b.get_component(corePosition);
//
// if(corePosition == 0) {
// Ai.reinterpret_dimensions({Ai.dimensions[1], Ai.dimensions[2], Ai.dimensions[3]});
// bi.reinterpret_dimensions({bi.dimensions[1], bi.dimensions[2]});
//
// op(i2, i3, j2, j3) = Ai(i2, j2, k2)*rightAStack[corePosition+1](i3, k2, j3);
// rhs(i2, i3) = bi(i2, k2)*rightBStack[corePosition+1](i3, k2);
// } else if(corePosition == d-1) {
// Ai.reinterpret_dimensions({Ai.dimensions[0], Ai.dimensions[1], Ai.dimensions[2]});
// bi.reinterpret_dimensions({bi.dimensions[0], bi.dimensions[1]});
//
// op(i1, i2, j1, j2) = leftAStack[corePosition-1](i1, k1, j1)*Ai(k1, i2, j2);
// rhs(i1, i2) = leftBStack[corePosition-1](i1, k1)*bi(k1, i2);
// } else {
// op(i1, i2, i3, j1, j2, j3) = leftAStack[corePosition-1](i1, k1, j1)*Ai(k1, i2, j2, k2)*rightAStack[corePosition+1](i3, k2, j3);
// rhs(i1, i2, i3) = leftBStack[corePosition-1](i1, k1)*bi(k1, i2, k2)*rightBStack[corePosition+1](i3, k2);
// }
//
// solve_least_squares(x.component(corePosition), op, rhs, 0);
//
// if(corePosition == 0) {
// x.component(corePosition).reinterpret_dimensions({1, x.component(corePosition).dimensions[0], x.component(corePosition).dimensions[1]});
// } else if(corePosition == d-1) {
// x.component(corePosition).reinterpret_dimensions({x.component(corePosition).dimensions[0], x.component(corePosition).dimensions[1], 1});
// }
// If we have not yet reached the end of the sweep we need to take care of the core and update our stacks
if(corePosition > 0) {
x.move_core(corePosition-1, true);
......@@ -158,13 +203,13 @@ namespace xerus {
}
}
}
};
}
};
void xals(TTTensor& _x, const TTOperator& _A, const TTTensor& _b) {
InternalSolver solver(_x, _A, _b);
return solver.solve();
}
void xals(TTTensor& _x, const TTOperator& _A, const TTTensor& _b) {
InternalSolver solver(_x, _A, _b);
return solver.solve();
}
}
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