ABlockedData.m

Go to the documentation of this file.

namespace data{


/* (Autoinserted by mtoc++)
 * This source code has been filtered by the mtoc++ executable,
 * which generates code that can be processed by the doxygen documentation tool.
 *
 * On the other hand, it can neither be interpreted by MATLAB, nor can it be compiled with a C++ compiler.
 * Except for the comments, the function bodies of your M-file functions are untouched.
 * Consequently, the FILTER_SOURCE_FILES doxygen switch (default in our Doxyfile.template) will produce
 * attached source files that are highly readable by humans.
 *
 * Additionally, links in the doxygen generated documentation to the source code of functions and class members refer to
 * the correct locations in the source code browser.
 * However, the line numbers most likely do not correspond to the line numbers in the original MATLAB source files.
 */

class ABlockedData
  :public handle {
    public:

        double MinRelSingularValueSize = 1e-20;
    public:

        function [matrix<double>U , matrix<double>S , matrix<double>V ] = getSVD(integer k,matrix<double> Vexclude,colvec<integer> targetdims) {
            if nargin < 4
                targetdims = " : ";
                if nargin < 3
                    Vexclude = [];
                end
            end
            /*  Get the correct dimensions if some are to be excluded */
            [n, m] = size(this);
            if ~ischar(targetdims)
                /*  Use the number of target dims for n if not all are used */
                n = length(targetdims);
            end
            /*  Assume full size if k is not set or empty */
            if nargin < 2 || isempty(k)
                k = min(n,m);
            end
            /*  Reduce target dimension if remaining dimension is smaller */
            if ~isempty(Vexclude)
                k = min(k, min(n,m)-size(Vexclude,2));
            else
                k = min(k, min(n,m));
            end

            /*  Fetch case for one block only */
            if this.getNumBlocks == 0
                U = [];
                S = [];
            elseif this.getNumBlocks == 1
                Bl = this.getBlock(1);
                Bl = Bl(targetdims,:);
                /*  Subtract exclude space if wanted */
                if ~isempty(Vexclude)
                    Bl = Bl - Vexclude*(Vexclude^t*Bl);
                end
                [U, S] = svd(Bl," econ ");
                U = U(:,1:k);
                S = S(1:k,1:k);
            else            
                /*  Fully-flavoured singular value decomp using block
                 * separation */
                opts.issym= 1;
                if KerMor.App.Verbose > 0
                    fprintf(" ABlockedData: Computing %d-partial SVD on %dx%d matrix (%d blocks)...\n ",...
                        k,n,m,this.getNumBlocks);
                end
                vb = KerMor.App.Verbose > 1 && n*m > 50000 && this.getNumBlocks > 2;
                cnt = 0;
                /*  If an exclude space is given, choose initial value outside of Vexclude */
                if ~isempty(Vexclude)
                    v0 = rand(n,1);
                    opts.v0= v0 - Vexclude*(Vexclude^t*v0);
                end
                doparallel = exist(" matlabpool "," file ") == 2 && matlabpool(" size ") > 1;
                [U,S] = eigs(@mult,n,k," la ",opts);
                if KerMor.App.Verbose > 2, fprintf(" BlockSVD: Finished after %d multiplications.\n ",cnt); end
                S = sqrt(S);
            end
            sel = sqrt(diag(S)/S(1)) >= this.MinRelSingularValueSize & diag(S)>0;
            U = U(:,sel);
            if size(U,2) < k
                warning(" KerMor:ABlockedData "," Have only %d nonzero singular values instead of %d desired ones. ",...
                    size(U,2),k);
            end
            S = S(sel,sel);
            if nargout > 2
                /*  Treat this special case here instead of cloning everything */
                if ~ischar(targetdims)
                    error(" Fixme: Not correctly implemented yet. ");
                    [nfull,~] = size(this);
                    mat = eye(nfull);
                    mat(targetdims,targetdims) = (U/S)^t;
                    V = mat*this;
                    V = V(targetdims,:);
                else
                    V = (U/S)^t*this;
                end
            end

            function w = mult(v)
                w = 0;
                nb = this.getNumBlocks;
                if vb && ~doparallel
                    pi = ProcessIndicator(" ABlockedData: Blockwise multiplication for %d-SVD on %d blocks ",...
                        nb,false,k,nb);
                end
                if doparallel
                    if vb
                        fprintf(" ABlockedData: Parallel blockwise multiplication for %d-SVD on %d blocks...\n ",k,nb);
                    end
                    parfor j = 1:nb
                        B = this.getBlock(j);/* #ok */

                        B = B(targetdims,:);
                        /*  Subtract exclude space if wanted */
                        if ~isempty(Vexclude)
                            B = B - Vexclude*(Vexclude^t*B);
                        end
                        w = w + B*(B^t*v);
                    end
                else               
                    for j = 1:nb
                        B = this.getBlock(j);
                        B = B(targetdims,:);
                        /*  Subtract exclude space if wanted */
                        if ~isempty(Vexclude)
                            B = B - Vexclude*(Vexclude^t*B);
                        end
                        w = w + B*(B^t*v);
                        if vb, pi.step; end
                    end
                end
                if vb, pi.stop; end
                cnt=cnt+1;
            end
        }
        function prod = mtimes() {
            if isa(matrix," data.ABlockedData ")
                error(" Must override in subclasses for case ABlockedData*matrix. ");
            else
                warning(" KerMor:ABlockedData ",...
                    " If you want the right singular values, you need to override the mtimes method of data.ABlockedData(matrix, this).\nReturning the same instance for V. ");
            end
            prod = this;
        }
        function A = toMemoryMatrix() {
            A = zeros(size(this));
            pos = 1;
            for i=1:this.getNumBlocks
                B = this.getBlock(i);
                s = size(B,2)-1;
                A(:,pos:pos+s) = B;
                pos = pos + s + 1;
            end
        }
    public: /* ( Abstract ) */

        virtual function varargout = size(dim) = 0;


        virtual function n = getNumBlocks() = 0;


        virtual function B = getBlock(nr) = 0;

    
    public: /* ( Static ) */

        static function res = test_BlockSVD_vs_SVD() {
            blocks = 5;
            ns = 6;
            an = [500 1000 700];
            am = [500 700 1000];
            res = true;
            for v = 1:3
                n = an(v);
                m = am(v);
                for nb=1:blocks
                    B = rand(n,m);
                    A = data.FileMatrix(n,m," BlockSize ",data.FileMatrix.blockSizeOf(B)/nb);
                    A.subsasgn(struct(" type ",[" () "]," subs ",[[" : ", " : "]]),B);
                    [ub,us] = A.getSVD(ns);
                    [u,s] = svd(B," econ ");
                    chk = norm(diag(us)-diag(s(1:ns,1:ns)));
                    res = res && chk < 1e-8;
                    /*  mess with directions - so correct */
                    ubs = sign(ub(1,:));
                    us = sign(u(1,1:ns));
                    chk2 = norm(bsxfun(@times,ub,ubs)-bsxfun(@times,u(:,1:ns),us));
                    res = res & chk2 < 1e-8;
                    fprintf(" Comparison Block-SVD vs SVD (%d singular values): Singular value vector norm diff: %g, U norm diff: %g\n ",ns,chk,chk2);
                end
            end
        }


};
}