Commit a7ade837 authored by Sebastian Wolf's avatar Sebastian Wolf

Further work on xALS

parent 2d23dabd
Pipeline #698 failed with stages
in 3 minutes and 37 seconds
......@@ -23,6 +23,8 @@
*/
#include <xerus/algorithms/xals.h>
#include <xerus/basic.h>
#include <xerus/misc/internal.h>
#include <xerus/indexedTensorMoveable.h>
#include <xerus/misc/basicArraySupport.h>
......@@ -39,19 +41,107 @@ namespace xerus {
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:
InternalSolver(TTTensor& _x, const TTOperator& _A, const TTTensor& _b) : d(_x.degree()), x(_x), A(_A), b(_b) {
InternalSolver(TTTensor& _x, const TTOperator& _A, const TTTensor& _b) : d(_x.degree()),leftStack(d), rightStack(d), x(_x), A(_A), b(_b) {
}
void solve() {
void calc_left_stack(const size_t _position) {
const Tensor& xi = x.get_component(_position);
if(_position == 0) {
Tensor tmp;
Tensor A0 = A.get_component(_position);
A0.reinterpret_dimensions({A0.dimensions[1], A0.dimensions[2], A0.dimensions[3]});
contract(tmp, A0, true, xi, false, 1);
contract(leftStack[_position], A0, true, xi, false, 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;
for(size_t iteration = 0; maxIterations == 0 || iteration < maxIterations; ++iteration) {
x.move_core(0, true);
// Rebuild right stack
for(size_t corePosition = d-1; corePosition > 0; --corePosition) {
calc_right_stack(corePosition);
}
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
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);
// 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);
calc_left_stack(corePosition);
}
}
// Sweep Right -> Left
for(size_t corePosition = d-1; corePosition > 0; --corePosition) {
//Actual Work to do
// 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);
calc_right_stack(corePosition);
}
}
}
}
};
void xals(TTTensor& _x, const TTOperator& _A, const TTTensor& _b) {
......
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