adf.h

Go to the documentation of this file.

// Xerus - A General Purpose Tensor Library
// Copyright (C) 2014-2017 Benjamin Huber and Sebastian Wolf. 
// 
// Xerus is free software: you can redistribute it and/or modify
// it under the terms of the GNU Affero General Public License as published
// by the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
// 
// Xerus is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Affero General Public License for more details.
// 
// You should have received a copy of the GNU Affero General Public License
// along with Xerus. If not, see <http://www.gnu.org/licenses/>.
//
// For further information on Xerus visit https://libXerus.org 
// or contact us at contact@libXerus.org.

#pragma once

#include "../index.h"
#include "../ttNetwork.h"
#include "../performanceData.h"
#include "../measurments.h"

namespace xerus {
    
    class ADFVariant {
        protected:
            
        template<class MeasurmentSet> 
        class InternalSolver {
            /*
             * General notation for the ADF:
             * The vector of the measured values is denoted as b
             * The measurment operator is denoted as A. It holds by definition A(x) = b (for noiseless measurments).
             * The Operator constructed from the x by removing the component at corePosition is denoted as E.
             */
            
        protected:
            const Index r1, r2, i1;
            
            TTTensor& x;
            
            const size_t degree;
            
            const std::vector<size_t> maxRanks;
            
            const MeasurmentSet& measurments;
            
            const size_t numMeasurments;
            
            const value_t normMeasuredValues;
            
            const size_t maxIterations; 
            
            const double targetResidualNorm;
            
            const double minimalResidualNormDecrease; 
            
            size_t iteration;
            
            double residualNorm;

            double lastResidualNorm;
            
            std::vector<value_t> residual;
            
            Tensor projectedGradientComponent;
            
            std::unique_ptr<Tensor*[]> forwardStackMem;
            
            Tensor* const * const forwardStack;
            
            std::unique_ptr<Tensor[]> forwardStackSaveSlots;
            
            std::vector<std::vector<size_t>> forwardUpdates;
            
            
            std::unique_ptr<Tensor*[]> backwardStackMem;
            
            Tensor* const * const backwardStack;

            std::unique_ptr<Tensor[]> backwardStackSaveSlots;
            
            std::vector<std::vector<size_t>> backwardUpdates;
            
            std::unique_ptr<double[]> measurmentNorms;
            
            PerformanceData& perfData;
            
            static double calculate_norm_of_measured_values(const MeasurmentSet& _measurments);
            
            void construct_stacks(std::unique_ptr< xerus::Tensor[] >& _stackSaveSlot, std::vector< std::vector< size_t > >& _updates, const std::unique_ptr<Tensor*[]>& _stackMem, const bool _forward);
            
            void resize_stack_tensors();
            
            std::vector<Tensor> get_fixed_components(const Tensor& _component);
            
            void update_forward_stack(const size_t _corePosition, const Tensor& _currentComponent);
            
            void update_backward_stack(const size_t _corePosition, const Tensor& _currentComponent);
            
            void calculate_residual( const size_t _corePosition );
            
            template<class PositionType>
            void perform_dyadic_product(const size_t _localLeftRank, const size_t _localRightRank, const value_t* const _leftPtr,  const value_t* const _rightPtr,  value_t* const _deltaPtr, const value_t _residual, const PositionType& _position, value_t* const _scratchSpace );
    
            
            void calculate_projected_gradient(const size_t _corePosition);
            
            std::vector<value_t> calculate_slicewise_norm_A_projGrad( const size_t _corePosition);
            
            void update_x(const std::vector<value_t>& _normAProjGrad, const size_t _corePosition);
            
            void solve_with_current_ranks();
            
        public:
            double solve();
            
            InternalSolver(TTTensor& _x, 
                            const std::vector<size_t>& _maxRanks, 
                            const MeasurmentSet& _measurments, 
                            const size_t _maxIteration, 
                            const double _targetResidualNorm, 
                            const double _minimalResidualNormDecrease, 
                            PerformanceData& _perfData ) : 
                x(_x),
                degree(_x.degree()),
                maxRanks(TTTensor::reduce_to_maximal_ranks(_maxRanks, _x.dimensions)),
                
                measurments(_measurments),
                numMeasurments(_measurments.size()),
                normMeasuredValues(calculate_norm_of_measured_values(_measurments)),
                
                maxIterations(_maxIteration),
                targetResidualNorm(_targetResidualNorm),
                minimalResidualNormDecrease(_minimalResidualNormDecrease),
                
                iteration(0),
                residualNorm(std::numeric_limits<double>::max()), 
                lastResidualNorm(std::numeric_limits<double>::max()),
                
                residual(numMeasurments),
                
                forwardStackMem(new Tensor*[numMeasurments*(degree+2)]),
                forwardStack(forwardStackMem.get()+numMeasurments),
                forwardUpdates(degree),
                    
                backwardStackMem(new Tensor*[numMeasurments*(degree+2)]),
                backwardStack(backwardStackMem.get()+numMeasurments),
                backwardUpdates(degree),
                
                measurmentNorms(new double[numMeasurments]),
                
                perfData(_perfData) 
                {
                    _x.require_correct_format();
                    XERUS_REQUIRE(numMeasurments > 0, "Need at very least one measurment.");
                    XERUS_REQUIRE(measurments.degree() == degree, "Measurment degree must coincide with x degree.");
                }
                
        };
        
    public:
        size_t maxIterations; 
        double targetResidualNorm; 
        double minimalResidualNormDecrease; // The minimal relative decrease of the residual per step  ( i.e. (lastResidual-residual)/lastResidual ). If the avg. of the last three steps is smaller than this value, the algorithm stops.
        
        ADFVariant(const size_t _maxIteration, const double _targetResidual, const double _minimalResidualDecrease)
                : maxIterations(_maxIteration), targetResidualNorm(_targetResidual), minimalResidualNormDecrease(_minimalResidualDecrease) { }
        
        template<class MeasurmentSet>
        double operator()(TTTensor& _x, const MeasurmentSet& _measurments, PerformanceData& _perfData) const {
            InternalSolver<MeasurmentSet> solver(_x, _x.ranks(), _measurments, maxIterations, targetResidualNorm, minimalResidualNormDecrease, _perfData);
            return solver.solve();
        }
        
        template<class MeasurmentSet>
        double operator()(TTTensor& _x, const MeasurmentSet& _measurments, const std::vector<size_t>& _maxRanks, PerformanceData& _perfData) const {
            InternalSolver<MeasurmentSet> solver(_x, _maxRanks, _measurments, maxIterations, targetResidualNorm, minimalResidualNormDecrease, _perfData);
            return solver.solve();
        }
    };
    
    extern const ADFVariant ADF;
}