als.cpp

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.

#include <xerus/misc/math.h>

#include <xerus/algorithms/als.h>
#include <xerus/basic.h>
#include <xerus/index.h>
#include <xerus/indexedTensorList.h>
#include <xerus/tensorNetwork.h>
#include <xerus/misc/internal.h>

#include <xerus/indexedTensorMoveable.h>
#include <xerus/indexedTensor_tensor_factorisations.h>

namespace xerus {

    // -------------------------------------------------------------------------------------------------------------------------
    //                                       local solvers
    // -------------------------------------------------------------------------------------------------------------------------
    
    void ALSVariant::lapack_solver(const TensorNetwork& _A, std::vector<Tensor>& _x, const TensorNetwork& _b, const ALSAlgorithmicData& _data) {
        Tensor A(_A);
        Tensor b(_b);
        Tensor x;
        
        xerus::solve(x, A, b);
        
        Index i,j,k,l;
        if (_data.direction == Increasing) {
            for (size_t p = 0; p+1 < _data.ALS.sites; ++p) {
                Tensor U, S;
//              calculate_svd(U, S, x, x, 2, _data.targetRank[_data.currIndex+p], EPSILON); TODO
                (U(i^2,j), S(j,k), x(k,l&1)) = SVD(x(i^2,l&2), _data.targetRank[_data.currIndex+p]);
                _x[p] = std::move(U);
                x(j,l&1) = S(j,k) * x(k,l&1);
            }
            _x.back() = std::move(x);
        } else {
            // direction: decreasing index
            for (size_t p = _data.ALS.sites-1; p>0; --p) {
                Tensor S, Vt;
//              calculate_svd(x, S, Vt, x, x.degree()-1, _data.targetRank[_data.currIndex+p-1], EPSILON); TODO
                (x(i&1,j), S(j,k), Vt(k,l&1)) = SVD(x(i&2,l^2), _data.targetRank[_data.currIndex+p-1]);
                _x[p] = std::move(Vt);
                x(i&1,k) = x(i&1,j) * S(j,k);
            }
            _x[0] = std::move(x);
        }
    }
    
    void ALSVariant::ASD_solver(const TensorNetwork &_A, std::vector<Tensor> &_x, const TensorNetwork &_b, const ALSAlgorithmicData &_data) {
        // performs a single gradient step, so
        // x = x + alpha * P( A^t (b - Ax) )    or for SPD: x = x + alpha * P( b - Ax )
        // where the projection P is already part of the stacks 
        // and alpha is the exact stepsize for quadratic functionals
        REQUIRE(_data.ALS.sites == 1, "ASD only defined for single site alternation at the moment");
        Tensor grad;
        Index i,j,k;
        grad(i&0) = _b(i&0) - _A(i/2,j/2) * _x[0](j&0);
        value_t alpha;
        if (_data.ALS.assumeSPD) {
            // stepsize alpha = <y,y>/<y,Ay>
            alpha = misc::sqr(frob_norm(grad)) / value_t(grad(i&0) * _A(i/2,j/2) * grad(j&0));
        } else {
            grad(i&0) = _A(j/2,i/2) * grad(j&0);
            // stepsize alpha = <y,y>/<Ay,Ay>
            alpha = frob_norm(grad) / frob_norm(_A(i/2,j/2) * grad(j&0));
        }
        _x[0] += alpha * grad;
    }
    
    // -------------------------------------------------------------------------------------------------------------------------
    //                                       helper functions
    // -------------------------------------------------------------------------------------------------------------------------
    
    void ALSVariant::ALSAlgorithmicData::prepare_x_for_als() {
        const size_t d = x.degree();
        Index r1,r2,n1,cr1;
        
        size_t firstOptimizedIndex = 0;
        size_t dimensionProd = 1;
        while (firstOptimizedIndex + 1 < d) {
            const size_t localDim = x.dimensions[firstOptimizedIndex];
            size_t newDimensionProd = dimensionProd * localDim;
            if (x.rank(firstOptimizedIndex) < newDimensionProd) {
                break;
            }
            
            Tensor& curComponent = x.component(firstOptimizedIndex);
            curComponent.reinterpret_dimensions({curComponent.dimensions[0]*curComponent.dimensions[1], curComponent.dimensions[2]});
            curComponent(r1,n1,r2) = curComponent(r1,cr1) * x.get_component(firstOptimizedIndex+1)(cr1,n1,r2);
            x.set_component(firstOptimizedIndex+1, std::move(curComponent));
            
            //TODO sparse
            x.set_component(firstOptimizedIndex, Tensor(
                {dimensionProd, localDim, newDimensionProd},
                [&](const std::vector<size_t> &_idx){
                    if (_idx[0]*localDim + _idx[1] == _idx[2]) {
                        return 1.0;
                    } 
                    return 0.0;
                })
            );
            
            x.require_correct_format();
            
            firstOptimizedIndex += 1;
            dimensionProd = newDimensionProd;
        }
        
        size_t firstNotOptimizedIndex = d;
        dimensionProd = 1;
        while (firstNotOptimizedIndex > firstOptimizedIndex+ALS.sites) {
            const size_t localDim = x.dimensions[firstNotOptimizedIndex-1];
            size_t newDimensionProd = dimensionProd * localDim;
            if (x.rank(firstNotOptimizedIndex-2) < newDimensionProd) {
                break;
            }
            
            Tensor& curComponent = x.component(firstNotOptimizedIndex-1);
            curComponent.reinterpret_dimensions({curComponent.dimensions[0], curComponent.dimensions[1] * curComponent.dimensions[2]});
            curComponent(r1,n1,r2) = x.get_component(firstNotOptimizedIndex-2)(r1,n1,cr1) * curComponent(cr1,r2);
            x.set_component(firstNotOptimizedIndex-2, std::move(curComponent));
            
            //TODO sparse
            x.set_component(firstNotOptimizedIndex-1, Tensor(
                {newDimensionProd, localDim, dimensionProd},
                [&](const std::vector<size_t> &_idx){
                    if (_idx[0] == _idx[1]*dimensionProd + _idx[2]) {
                        return 1.0;
                    } 
                    return 0.0;
                })
            );

            x.require_correct_format();
            
            firstNotOptimizedIndex -= 1;
            dimensionProd = newDimensionProd;
        }
        
        if (canonicalizeAtTheEnd && corePosAtTheEnd < firstOptimizedIndex) {
            x.assume_core_position(firstOptimizedIndex);
        } else {
            if (canonicalizeAtTheEnd && corePosAtTheEnd >= firstNotOptimizedIndex) {
                x.assume_core_position(firstNotOptimizedIndex-1);
            }
            
            x.move_core(firstOptimizedIndex, true);
        }
        
        optimizedRange = std::pair<size_t, size_t>(firstOptimizedIndex, firstNotOptimizedIndex);
    }

    TensorNetwork ALSVariant::ALSAlgorithmicData::localOperatorSlice(size_t _pos) {
        //TODO optimization: create these networks without indices
        Index cr1, cr2, cr3, cr4, r1, r2, r3, r4, n1, n2, n3;
        TensorNetwork res;
        INTERNAL_CHECK(A, "ie");
        if (ALS.assumeSPD) {
            res(r1,r2,r3, cr1,cr2,cr3) = x.get_component(_pos)(r1, n1, cr1) 
                                        * A->get_component(_pos)(r2, n1, n2, cr2) 
                                        * x.get_component(_pos)(r3, n2, cr3);
        } else {
            res(r1,r2,r3,r4, cr1,cr2,cr3, cr4) = x.get_component(_pos)(r1, n1, cr1) 
                                        * A->get_component(_pos)(r2, n2, n1, cr2) 
                                        * A->get_component(_pos)(r3, n2, n3, cr3) 
                                        * x.get_component(_pos)(r4, n3, cr4);
        }
        return res;
    }
    
    TensorNetwork ALSVariant::ALSAlgorithmicData::localRhsSlice(size_t _pos) {
        //TODO optimization: create these networks without indices
        Index cr1, cr2, cr3, r1, r2, r3, n1, n2;
        TensorNetwork res;
        if (ALS.assumeSPD || (A == nullptr)) {
            res(r1,r2, cr1,cr2) = b.get_component(_pos)(r1, n1, cr1) 
                                    * x.get_component(_pos)(r2, n1, cr2);
        } else {
            res(r1,r2,r3, cr1,cr2,cr3) = b.get_component(_pos)(r1, n1, cr1) 
                                                * A->get_component(_pos)(r2, n1, n2, cr2) 
                                                * x.get_component(_pos)(r3, n2, cr3);
        }
        return res;
    }
    
    void ALSVariant::ALSAlgorithmicData::prepare_stacks() {
        const size_t d = x.degree();
        Index r1,r2;
        
        Tensor tmpA;
        Tensor tmpB;
        if (ALS.assumeSPD || (A == nullptr)) {
            tmpA = Tensor::ones({1,1,1});
            tmpB = Tensor::ones({1,1});
        } else {
            tmpA = Tensor::ones({1,1,1,1});
            tmpB = Tensor::ones({1,1,1});
        }
        
        
        localOperatorCache.left.emplace_back(tmpA);
        localOperatorCache.right.emplace_back(tmpA);
        rhsCache.left.emplace_back(tmpB);
        rhsCache.right.emplace_back(tmpB);
        
        for (size_t i = d-1; i > optimizedRange.first + ALS.sites - 1; --i) {
            if (A != nullptr) {
                tmpA(r1&0) = localOperatorCache.right.back()(r2&0) * localOperatorSlice(i)(r1/2, r2/2);
                localOperatorCache.right.emplace_back(tmpA);
            }
            tmpB(r1&0) = rhsCache.right.back()(r2&0) * localRhsSlice(i)(r1/2, r2/2);
            rhsCache.right.emplace_back(tmpB);
        }
        for (size_t i = 0; i < optimizedRange.first; ++i) {
            if (A != nullptr) {
                tmpA(r2&0) = localOperatorCache.left.back()(r1&0) * localOperatorSlice(i)(r1/2, r2/2);
                localOperatorCache.left.emplace_back((tmpA));
            }
            tmpB(r2&0) = rhsCache.left.back()(r1&0) * localRhsSlice(i)(r1/2, r2/2);
            rhsCache.left.emplace_back(tmpB);
        }
    }
    
    void ALSVariant::ALSAlgorithmicData::choose_energy_functional() {
        if (A != nullptr) {
            if (ALS.assumeSPD) {
                residual_f = [&](){
                    Index n1, n2;
                    return frob_norm((*A)(n1/2,n2/2)*x(n2&0) - b(n1&0))/normB;
                };
                if (ALS.useResidualForEndCriterion) {
                    energy_f = residual_f;
                } else {
                    energy_f = [&](){
                        Index r1,r2;
                        Tensor res;
                        // 0.5*<x,Ax> - <x,b>
                        TensorNetwork xAx = localOperatorCache.left.back();
                        TensorNetwork bx = rhsCache.left.back();
                        for (size_t i=0; i<ALS.sites; ++i) {
                            xAx(r2&0) = xAx(r1&0) * localOperatorSlice(currIndex+i)(r1/2, r2/2);
                            bx(r2&0) = bx(r1&0) * localRhsSlice(currIndex+i)(r1/2, r2/2);
                        }
                        res() = 0.5*xAx(r1&0) * localOperatorCache.right.back()(r1&0)
                                - bx(r1&0) * rhsCache.right.back()(r1&0);
                        return res.frob_norm();
                    };
                }
            } else {
                // not Symmetric pos def
                residual_f = [&](){
                    Index r1,r2;
                    Tensor res;
                    // <Ax,Ax> - 2 * <Ax,b> + <b,b>
                    TensorNetwork xAtAx = localOperatorCache.left.back();
                    TensorNetwork bAx = rhsCache.left.back();
                    for (size_t i=0; i<ALS.sites; ++i) {
                        xAtAx(r2&0) = xAtAx(r1&0) * localOperatorSlice(currIndex+i)(r1/2, r2/2);
                        bAx(r2&0) = bAx(r1&0) * localRhsSlice(currIndex+i)(r1/2, r2/2);
                    }
                    res() = xAtAx(r1&0) * localOperatorCache.right.back()(r1&0)
                            - 2 * bAx(r1&0) * rhsCache.right.back()(r1&0);
                    return std::sqrt(res[0] + misc::sqr(normB))/normB;
                };
                energy_f = residual_f;
            }
        } else {
            // no operator A given
            residual_f = [&](){
                return frob_norm(x - b);
            };
            if (ALS.useResidualForEndCriterion) {
                energy_f = residual_f;
            } else {
                energy_f = [&](){
                    Index r1,r2;
                    Tensor res;
                    // 0.5*<x,Ax> - <x,b> = 0.5*|x_i|^2 - <x,b>
                    TensorNetwork bx = rhsCache.left.back();
                    for (size_t i=0; i<ALS.sites; ++i) {
                        bx(r2&0) = bx(r1&0) * localRhsSlice(currIndex+i)(r1/2, r2/2);
                    }
                    res() = 0.5*x.get_component(currIndex)(r1&0) * x.get_component(currIndex)(r1&0) 
                            - bx(r1&0) * rhsCache.right.back()(r1&0);
                    return res[0];
                };
            }
        }
    }
    
    ALSVariant::ALSAlgorithmicData::ALSAlgorithmicData(const ALSVariant &_ALS, const TTOperator *_A, TTTensor &_x, const TTTensor &_b) 
        : ALS(_ALS), A(_A), x(_x), b(_b)
        , targetRank(_x.ranks())
        , normB(frob_norm(_b))
        , canonicalizeAtTheEnd(_x.canonicalized)
        , corePosAtTheEnd(_x.corePosition)
        , lastEnergy2(1e102)
        , lastEnergy(1e101)
        , energy(1e100)
        , halfSweepCount(0)
        , direction(Increasing)
    {
        prepare_x_for_als();
        prepare_stacks();
        currIndex = optimizedRange.first;
        choose_energy_functional();
    }

    void ALSVariant::ALSAlgorithmicData::move_to_next_index() {
        Index r1,r2;
        Tensor tmpA, tmpB;
        if (direction == Increasing) {
            INTERNAL_CHECK(currIndex+ALS.sites < optimizedRange.second, "ie " << currIndex << " " << ALS.sites << " " << optimizedRange.first << " " << optimizedRange.second);
            // Move core to next position (assumed to be done by the solver if sites > 1)
            if (ALS.sites == 1) {
                x.move_core(currIndex+1, true);
            }
            
            // Move one site to the right
            if (A != nullptr) {
                localOperatorCache.right.pop_back();
                tmpA(r2&0) = localOperatorCache.left.back()(r1&0) * localOperatorSlice(currIndex)(r1/2, r2/2);
                localOperatorCache.left.emplace_back(std::move(tmpA));
            }
            
            rhsCache.right.pop_back();
            tmpB(r2&0) = rhsCache.left.back()(r1&0) * localRhsSlice(currIndex)(r1/2, r2/2);
            rhsCache.left.emplace_back(std::move(tmpB));
            currIndex++;
        } else {
            INTERNAL_CHECK(currIndex > optimizedRange.first, "ie");
            // Move core to next position (assumed to be done by the solver if sites > 1)
            if (ALS.sites == 1) {
                x.move_core(currIndex-1, true);
            }
            
            // move one site to the left
            if (A != nullptr) {
                localOperatorCache.left.pop_back();
                tmpA(r1&0) = localOperatorCache.right.back()(r2&0) * localOperatorSlice(currIndex)(r1/2, r2/2);
                localOperatorCache.right.emplace_back(std::move(tmpA));
            }
            
            rhsCache.left.pop_back();
            tmpB(r1&0) = rhsCache.right.back()(r2&0) * localRhsSlice(currIndex)(r1/2, r2/2);
            rhsCache.right.emplace_back(std::move(tmpB));
            currIndex--;
        }
    }
    
    
    TensorNetwork ALSVariant::construct_local_operator(ALSVariant::ALSAlgorithmicData& _data) const {
        INTERNAL_CHECK(_data.A, "IE");
        Index cr1, cr2, cr3, cr4, r1, r2, r3, r4, n1, n2, n3, n4, x;
        TensorNetwork ATilde = _data.localOperatorCache.left.back();
        if (assumeSPD) {
            for (size_t p=0; p<sites; ++p) {
                ATilde(n1^(p+1), n2, r2, n3^(p+1), n4) = ATilde(n1^(p+1), r1, n3^(p+1)) * _data.A->get_component(_data.currIndex+p)(r1, n2, n4, r2);
            }
            ATilde(n1^(sites+1), n2, n3^(sites+1), n4) = ATilde(n1^(sites+1), r1, n3^(sites+1)) * _data.localOperatorCache.right.back()(n2, r1, n4);
        } else {
            for (size_t p=0;  p<sites; ++p) {
                ATilde(n1^(p+1),n2, r3,r4, n3^(p+1),n4) = ATilde(n1^(p+1), r1,r2, n3^(p+1))
                        * _data.A->get_component(_data.currIndex+p)(r1, x, n2, r3)
                        * _data.A->get_component(_data.currIndex+p)(r2, x, n4, r4);
            }
            ATilde(n1^(sites+1), n2, n3^(sites+1), n4) = ATilde(n1^(sites+1), r1^2, n3^(sites+1)) * _data.localOperatorCache.right.back()(n2, r1^2, n4);
        }
        return ATilde;
    }
    
    
    TensorNetwork ALSVariant::construct_local_RHS(ALSVariant::ALSAlgorithmicData& _data) const {
        Index cr1, cr2, cr3, cr4, r1, r2, r3, r4, n1, n2, n3, n4, x;
        TensorNetwork BTilde;
        if (assumeSPD || (_data.A == nullptr)) {
            BTilde(n1,r1) = _data.rhsCache.left.back()(r1,n1);
            for (size_t p=0; p<sites; ++p) {
                BTilde(n1^(p+1), n2, cr1) = BTilde(n1^(p+1), r1) * _data.b.get_component(_data.currIndex+p)(r1, n2, cr1);
            }
            BTilde(n1^(sites+1),n2) = BTilde(n1^(sites+1), r1) * _data.rhsCache.right.back()(r1,n2);
        } else {
            BTilde(n1,r1^2) = _data.rhsCache.left.back()(r1^2,n1);
            for (size_t p=0; p<sites; ++p) {
                BTilde(n1^(p+1), n3, cr1, cr2) = BTilde(n1^(p+1), r1, r2) 
                    * _data.b.get_component(_data.currIndex+p)(r1, n2, cr1)
                    * _data.A->get_component(_data.currIndex+p)(r2, n2, n3, cr2);
            }
            BTilde(n1^(sites+1),n2) = BTilde(n1^(sites+1), r1^2) * _data.rhsCache.right.back()(r1^2,n2);
        }
        return BTilde;
    }


    bool ALSVariant::check_for_end_of_sweep(ALSAlgorithmicData& _data, size_t _numHalfSweeps, value_t _convergenceEpsilon, PerformanceData &_perfData) const {
        if ((_data.direction == Decreasing && _data.currIndex==_data.optimizedRange.first) 
            || (_data.direction == Increasing && _data.currIndex==_data.optimizedRange.second-sites)) 
        {
            LOG(ALS, "Sweep Done");
            _data.halfSweepCount += 1;
            
            _data.lastEnergy2 = _data.lastEnergy;
            _data.lastEnergy = _data.energy;
            _data.energy = _data.energy_f();
            
            if (_perfData) {
                size_t flags = _data.direction == Increasing ? FLAG_FINISHED_HALFSWEEP : FLAG_FINISHED_FULLSWEEP;
                if (!useResidualForEndCriterion) {
                    _perfData.stop_timer();
                    value_t residual = _data.residual_f();
                    _perfData.continue_timer();
                    _perfData.add(residual, _data.x, flags);
                } else {
                    _perfData.add(_data.energy, _data.x, flags);
                }
            }
            
            // Conditions for loop termination
            if (_data.halfSweepCount == _numHalfSweeps 
                || std::abs(_data.lastEnergy-_data.energy) < _convergenceEpsilon 
                || std::abs(_data.lastEnergy2-_data.energy) < _convergenceEpsilon 
                || (_data.optimizedRange.second - _data.optimizedRange.first<=sites)) 
            {
                // we are done! yay
                LOG(ALS, "ALS done, " << _data.energy << " " << _data.lastEnergy << " " << std::abs(_data.lastEnergy2-_data.energy) << " " << std::abs(_data.lastEnergy-_data.energy) << " < " << _convergenceEpsilon);
                if (_data.canonicalizeAtTheEnd && preserveCorePosition) {
                    _data.x.move_core(_data.corePosAtTheEnd, true);
                }
                return true;
            }
                
            // change walk direction
            _data.direction = _data.direction == Increasing ? Decreasing : Increasing;
        } else {
            // we are not done with the sweep, just calculate data for the perfdata
            if (_perfData) {
                _perfData.stop_timer();
                value_t residual = _data.residual_f();
                _perfData.continue_timer();
                _perfData.add(residual, _data.x, 0);
            }
        }
        return false;
    }
    
    
    // -------------------------------------------------------------------------------------------------------------------------
    //                                   the actual algorithm
    // -------------------------------------------------------------------------------------------------------------------------
    
    
    double ALSVariant::solve(const TTOperator *_Ap, TTTensor &_x, const TTTensor &_b, size_t _numHalfSweeps, value_t _convergenceEpsilon, PerformanceData &_perfData) const {
        LOG(ALS, "ALS("<< sites <<") called");
        #ifndef XERUS_DISABLE_RUNTIME_CHECKS
            _x.require_correct_format();
            _b.require_correct_format();
            REQUIRE(_x.degree() > 0, "");
            REQUIRE(_x.dimensions == _b.dimensions, "");
            
            if (_Ap != nullptr) {
                _Ap->require_correct_format();
                REQUIRE(_Ap->dimensions.size() == _b.dimensions.size()*2, "");
                for (size_t i=0; i<_x.dimensions.size(); ++i) {
                    REQUIRE(_Ap->dimensions[i] == _x.dimensions[i], "");
                    REQUIRE(_Ap->dimensions[i+_Ap->degree()/2] == _x.dimensions[i], "");
                }
            }
        #endif
        
        if (_Ap != nullptr) {
            _perfData << "ALS for ||A*x - b||^2, x.dimensions: " << _x.dimensions << '\n'
                    << "A.ranks: " << _Ap->ranks() << '\n'
                    << "x.ranks: " << _x.ranks() << '\n'
                    << "b.ranks: " << _b.ranks() << '\n'
                    << "maximum number of half sweeps: " << _numHalfSweeps << '\n'
                    << "convergence epsilon: " << _convergenceEpsilon << '\n';
        } else {
            _perfData << "ALS for ||x - b||^2, x.dimensions: " << _x.dimensions << '\n'
                    << "x.ranks: " << _x.ranks() << '\n'
                    << "b.ranks: " << _b.ranks() << '\n'
                    << "maximum number of half sweeps: " << _numHalfSweeps << '\n'
                    << "convergence epsilon: " << _convergenceEpsilon << '\n';
        }
        _perfData.start();
        
        ALSAlgorithmicData data(*this, _Ap, _x, _b);
        
        data.energy = data.energy_f();
        
        if (_perfData) {
            _perfData.stop_timer();
            value_t residual = data.residual_f();
            _perfData.continue_timer();
            _perfData.add(residual, _x, FLAG_FINISHED_FULLSWEEP);
        }
        
        while (true) {
            LOG(ALS, "Starting to optimize index " << data.currIndex);
            
            // update current component tensor
            if (_Ap != nullptr) {
                std::vector<Tensor> tmpX;
                for (size_t p=0; p<sites; ++p) {
                    tmpX.emplace_back(_x.get_component(data.currIndex+p));
                }
                localSolver(construct_local_operator(data), tmpX, construct_local_RHS(data), data);
                for (size_t p=0; p<sites; ++p) {
                    _x.set_component(data.currIndex+p, std::move(tmpX[p]));
                }
            } else {
                //TODO?
                REQUIRE(sites==1, "approximation dmrg not implemented yet");
                _x.component(data.currIndex) = Tensor(construct_local_RHS(data));
            }
            
            if(check_for_end_of_sweep(data, _numHalfSweeps, _convergenceEpsilon, _perfData)) {
                return data.energy; // TODO residual?
            }
            
            data.move_to_next_index();
        }
    }
    
    
    const ALSVariant ALS(1, 0, ALSVariant::lapack_solver, false);
    const ALSVariant ALS_SPD(1, 0, ALSVariant::lapack_solver, true);
    
    const ALSVariant DMRG(2, 0, ALSVariant::lapack_solver, false);
    const ALSVariant DMRG_SPD(2, 0, ALSVariant::lapack_solver, true);
    
    const ALSVariant ASD(1, 0, ALSVariant::ASD_solver, false);
    const ALSVariant ASD_SPD(1, 0, ALSVariant::ASD_solver, true);
} // namespace xerus