PowerIteration.m

classdef PowerIteration < ARE.GammaCalculation.GammaCalculatorInterface
% POWERITERATION
% Calculate the (approximate) value of gamma by applying a power iteration
% to the matrix H, where H is never formed explicitly but only vector
% products are calculated H*x.
%
%  Andreas Schmidt, 2016

    properties
        % Stop the power iteration when the relative change is below this
        % value:
        pi_tolerance = 0.001;
        
        % Maximum number of iterations for the power iteration:
        pi_max_iterations = 100;
        
        % Don't accumulate summands with ||s|| <= tol_inner
        tol_inner = 1e-2;
    end

    methods
        function gamma = calculate(this, model, model_data, dsim, verbose)
            if nargin < 5, verbose=false; end
            
            [E,A,B,~,~,R] = model.assemble(model_data);
            if exist('dsim','var') && ~isempty(dsim)
                K = -inv(R + B'*dsim.Z*dsim.Z'*B)*B'*dsim.Z*dsim.Z'*A;
            else
                K = sparse(model.m, model.n);
            end
            
            gamma = this.do_power_iteration(A, E, B, K, verbose);
        end
        function gen_detailed_data(this, model, model_data)
        end
        function gen_reduced_data(this, model, detailed_data)
        end
        
        function gamma = do_power_iteration(this, A, E, B, K, verbose)
            if nargin < 3, E = speye(size(A,1)); end
            if nargin < 6, verbose=false; end
            
            tol = this.pi_tolerance;
            toli = this.tol_inner;
            maxit = this.pi_max_iterations;
            
            % Random start vector
            x0 = rand(size(A,1), 1);
            x0 = x0/norm(x0);

            [eL, eU] = lu(E);
            
            gamma_old = 0;
            
            % MAIN LOOP:
            for k = 1:maxit
                % Do the actual calculation H*x
                % Calculate E^(-1)*x0
                y = eU\(eL\x0);
                
                % hx will be the product H*x later on!
                hx = y;
                
                % Forward solve:
                i = 0;
                
                cache = y;
                while true
                    cache = eU\(eL\(A*cache + B*(K*cache)));
                    i = i + 1;
                    
                    % And the backward stuff:
                    tmp = cache;
                    for j = 1:i
                        tmp_ = eL'\(eU'\tmp);
                        tmp = A'*tmp_ + K'*(B'*tmp_);
                    end
                    
                    hx = hx + tmp;
                    if tmp < toli
                        if verbose 
                            fprintf('Needed %d summands\n', i);
                        end
                        break
                    end
                end

                hx = eL'\(eU'\hx);
                
                gamma_approx = norm(hx);
                rel_change = abs(gamma_approx-gamma_old)/abs(gamma_old);

                gamma_old = gamma_approx;
                x0 = hx/norm(hx);
                if verbose
                    fprintf('%d %d\n', [gamma_approx, rel_change]);
                end
                if rel_change < tol, break, end
            end
            gamma = gamma_approx;
        end
    end
end