Add files via upload

codes/calculateEigenValueA.m

@@ -0,0 +1,63 @@
+function opts = calculateEigenValueA(A,szy,Ny,opts)
+% Calculates the Eigenvalue Decomposition of the N dimensional tensor A,
+% asuming it is used as a design matrix in a regression problem. 
+% PROPACK package from http://sun.stanford.edu/~rmunk/PROPACK/ is used for
+% fast SVD calculation.
+% INPUTS
+% A    : N dimensional tensor to be decomposed
+% szy  : size of the observation tensor
+% Ny   : Number of dimensions of the observation tensor
+% opts : several parameters for the algorithm, check
+%        tensorRegressionCPpenSL.m for more info
+% OUTPUTS
+% opts: the resulting eigenvectors and eigenvalues are stored in the fields 
+%        A.V and A.ev respectively
+%
+% Version 1 - May 2015 
+
+% This is required to use Propack functions 
+global AA
+global Afunc;
+% Find the fastest way for matrix multiplication in Matlab
+AA = sprand(10,10,0.1);
+forwardType = findBestMultiply(AA,.2);
+Afunc       = sprintf('AforwardTen_%d',forwardType);
+
+cmoda = opts.cmoda;
+cmodx = opts.cmodx;
+
+sza = size(A); Na = ndims(A);
+
+Nx = length(cmodx) + (Ny-(Na-length(cmoda)));
+szx(cmodx) = sza(cmoda);
+szx((cmodx(end)+1):Nx) = setdiff(szy,sza(setdiff(1:Na,cmoda)));
+
+cmodat = length(cmoda)+1:Na; % contraction modes for A'A
+cmody  = 1:length(cmodat);   % contraction modes for A'Y
+
+At = permute(A,[cmoda setdiff(1:Na,cmoda)]);
+
+if prod(sza(cmody))>=prod(sza(cmoda))
+    AtA = contractTensor(At,A,cmodat,cmody,'m');
+else
+    AtA = contractTensor(A,At,cmoda,cmodx,'m');
+end
+
+AA = AtA;
+[V,ev]  = laneig(Afunc,size(AA,1),size(AA,1),'LM'); 
+
+% To use Matlab builtin function eig replace laneig with this line
+% [V,ev] = eig(AtA);
+ev = diag(ev);
+[ev,ind] = sort(ev,'descend');
+V  = V(:,ind);
+ev(ev<1e-14*ev(1))=0;
+r  = sum(ev>0);
+ev = ev(1:r);
+V  = V(:,1:r);
+
+opts.A.V = V;
+opts.A.ev= ev;
+
+
+end

codes/contractTensor.m

@@ -0,0 +1,38 @@
+function C = contractTensor(a,b,adims,bdims,varargin)
+% tensor tensor contraction
+% C = TTT(a,b,adims,bdims) computes the contracted product of
+%   tensors a and b in the dimensions specified by the row vectors
+%   adims and bdims.  The sizes of the dimensions specified by adims
+%   and bdims must match; that is, size(a,adims) must equal
+%   size(b,bdims).
+% opts: 't' result is returned as tensor  
+%       'm' result is returned as matrix
+% Modified ttt function of Tensor Toolbox of Kolda et. al.
+% MATLAB Tensor Toolbox.
+% Copyright 2010, Sandia Corporation. 
+
+    if nargin >4
+        opts = varargin{1};
+    else
+        opts = 't';
+    end
+
+    A = matricize(a,adims,'t');
+    B = matricize(b,bdims);
+    tsiz = [A.tsize(A.rindices) B.tsize(B.cindices)];
+    rindices = 1:length(A.rindices);
+    cindices = (1:length(B.cindices)) + length(A.rindices);
+    C = A.data * B.data;
+
+    if opts == 't'
+        order = [rindices,cindices];
+        data = reshape(C, [tsiz(order) 1 1]);
+        if numel(order) >= 2
+            C = ipermute(data,order);
+        else
+            C = data;
+        end
+    end
+
+end
+

codes/findBestMultiply.m

@@ -0,0 +1,78 @@
+function [forwardType,transposeType,details] = findBestMultiply(A,T)
+% [forwardType,transposeType] = findBestMultiply(A)
+%   does a speed test to determine what the fastest routine
+%   for sparse matrix multiplication is, for the matrix A.
+% ... = findBestMultiply(A,T)
+%   tries to make each test last T seconds (there are a 6 tests total)
+%
+% The fastest routine will depend on the verson of Matlab, the
+% computer OS, the cpu type (dual core? etc.), and the particular
+% size and structure of the matrix.
+%
+% Stephen Becker, srbecker@caltech.edu, 3/14/09
+
+% 5/13/09, smvp now handles complex data, so modify this accordingly
+
+if nargin < 2, T = 1; end
+
+At = A';
+x = randn(size(A,2),1);
+y = randn(size(A,1),1);
+if ~isreal(A)
+    x = x + 1i*randn(size(A,2),1);
+    y = y + 1i*randn(size(A,1),1);
+end
+
+
+% how many times to repeat the test?
+tic; A'*y; t = toc;
+
+% if we want each test to last say, T, seconds, then...
+nIter = max( 1, round(T/t) );
+
+
+% -- for the forward multiply --
+
+% -- Method 1: use standard MATLAB
+tic; for i = 1:nIter, A*x; end; t1 = toc;
+
+% -- Method 2: use transpose trick (works best in new versions)
+tic; for i = 1:nIter, At'*x; end; t2 = toc;
+
+% % -- Method 3: simple mex file
+% try
+%     tic; for i = 1:nIter, smvp(A,x); end; t3 = toc;
+% catch
+%     l = lasterror;
+%     fprintf('Error thrown: %s\n',l.message);
+%     disp('findBestMultiply: Sorry, this probably means the smvp mex file\n\t is not installed');
+%     disp('Try compiling it with:  mex -O smvp.c  in the "private" subdirectory');
+%     t3 = Inf;
+% end
+
+% [m,forwardType] = min( [t1,t2,t3] );
+[m,forwardType] = min( [t1,t2] );
+if nargout > 2, details = [t1,t2,t3]; end
+
+
+% -- for the transpose multiply --
+
+% -- Method 1: use standard MATLAB
+tic; for i = 1:nIter, At*y; end; t1 = toc;
+
+% -- Method 2: use transpose trick (works best in new versions)
+tic; for i = 1:nIter, A'*y; end; t2 = toc;
+
+% % -- Method 3: simple mex file
+% try
+%     tic; for i = 1:nIter, smvp(At,y); end; t3 = toc;
+% catch
+% %     disp('Sorry, smvp mex file not installed');
+%     t3 = Inf;
+% end
+% 
+% [m,transposeType] = min( [t1,t2,t3] );
+% if nargout > 2, details = [details,[t1,t2,t3] ]; end
+[m,transposeType] = min( [t1,t2] );
+if nargout > 2, details = [details,[t1,t2] ]; end
+

codes/initializeFactor.m

@@ -0,0 +1,185 @@
+function U = initializeFactor(X,R,opts)
+% Initialization methods for tensor decomposition
+% INPUTS:
+% X    : N dimensional tensor
+% R    : number of components or model order of PARAFAC
+% opts : structure containing the information
+%     opts.nn   -> vector for nonnegativity constraint on the factors
+%     opts.init -> methods to initialize: random, nvecs, nn-svd  (default random)
+%                  Every factor can be initizalized with a different method
+%                  by making opts.init a cell array.
+%     opts.initmethod -> option for nn-svd function: mean, rand (default rand)
+%
+% References:
+% o MATLAB Tensor Toolbox. Copyright 2012, Sandia Corporation.
+% o Boutsidis et.al.,"SVD based initialization: A head start for 
+%   nonnegative matrix factorization", 2008
+%
+% Version 1 - May 2015 
+%
+
+% take the dimensions of the tensor
+szx = size(X);
+N   = length(szx);
+U   = cell(N,1);  
+try nonnegative = opts.nonnegative;
+catch, nonnegative = opts.nn; end    
+
+if ~isfield(opts,'init')
+    opts.init = 'random';
+end
+
+% make the opts.init as cell array if it is not
+if ~iscell(opts.init)
+    temp = cell(N,1);
+    for n = 1:N
+        temp{n} = opts.init;
+    end
+    opts.init = temp;
+end
+
+% loop over all factors
+for n = 1:N
+    switch lower(opts.init{n})
+        case 'random'
+            if nonnegative(n)
+                U{n} = posrandn(szx(n),R);
+            else
+                U{n} = randn(szx(n),R);
+            end
+        case 'nvecs'
+            if nonnegative(n)
+                fprintf('This method may give real factors \n');
+                fprintf('We will force them to be nonnegative by taking absolute value \n');
+                U{n} = nvecs(X,n,R);
+                U{n} = abs(U{n});
+            else
+                U{n} = nvecs(X,n,R);
+            end
+        case 'nn-svd'
+            fprintf('This method will give NN factor \n');
+            if ~isfield(opts,'initmethod')
+                opts.initmethod = 'mean';
+            end
+            U{n} = cp_svd_init(X,n,R,opts);
+    end
+end
+
+for n = 1:N
+    U{n} = scaleFactor(U{n});
+end
+end
+
+% -----------------------------------------------------------------------
+
+function U = cp_svd_init(X,n,R,opts)
+
+% Start nonnegative tensor factorization by using an svd-based initialization method
+% which gives nonnegative factor matrices
+
+% Ref: Boutsidis et.al.,"SVD based initialization: A head start for
+% nonnegative matrix factorization", 2008
+% If X is dense, fill the zero values in the inital factors with 
+%     opts.initmethod: 'mean' values of the tensor or 
+%                      'rand' uniform random values from [0,mean(X(:))/100]
+
+    szx = size(X);
+    N   = length(szx);
+    U   = zeros(szx(n),R); 
+    Xn  = permute(X,[n 1:n-1,n+1:N]);
+    Xn  = reshape(Xn,szx(n),prod(szx([1:n-1 n+1:N])));
+
+    [LV,D,RV] = svds(Xn,R,'L'); 
+    d = diag(D);
+
+    for j = 1:R
+        % left sv
+        x  = LV(:,j);
+        xp = (x>=0).*x;
+        xn = (x<0).*(-x);
+        xpnrm = sqrt(sum(xp.^2));
+        xnnrm = sqrt(sum(xn.^2));
+        % right sv
+        y  = RV(:,j);
+        yp = (y>=0).*y;
+        yn = (y<0).*(-y);
+        ypnrm = sqrt(sum(yp.^2));
+        ynnrm = sqrt(sum(yn.^2));
+        mp = xpnrm*ypnrm;
+        mn = xnnrm*ynnrm;
+        if mp>mn
+            u = xp/xpnrm;
+            sigma = mp;
+        else
+            u = xn/xnnrm;
+            sigma = mn;
+        end
+        U(:,j) = sqrt(d(j)*sigma)*u;
+    end
+
+    if strcmp(opts.initmethod,'mean')
+        U(U==0) = mean(Xn(:));
+    end
+
+    if strcmp(opts.initmethod,'rand')
+        U(U==0) = (mean(Xn(:))/100).*rand(Xn(:));
+    end
+
+end
+
+
+% -------------------------------------------------------------------------
+% From tensor toolbox
+function u = nvecs(X,n,r)
+% adapted from tensor toolbox
+% @tensor/nvecs.m
+%   NVECS Compute the leading mode-n vectors for a tensor.
+%
+%   U = NVECS(X,n,r) computes the r leading eigenvalues of Xn*Xn'
+%   (where Xn is the mode-n matricization of X), which provides
+%   information about the mode-n fibers. In two-dimensions, the r
+%   leading mode-1 vectors are the same as the r left singular vectors
+%   and the r leading mode-2 vectors are the same as the r right
+%   singular vectors.
+%
+%   U = NVECS(X,n,r,OPTS) specifies opts:
+%   OPTS.eigsopts: options passed to the EIGS routine [struct('disp',0)]
+%   OPTS.flipsign: make each column's largest element positive [true]
+%
+%   See also TENSOR, TENMAT, EIGS.
+%
+%MATLAB Tensor Toolbox.
+%Copyright 2012, Sandia Corporation.
+
+% This is the MATLAB Tensor Toolbox by T. Kolda, B. Bader, and others.
+% http://www.sandia.gov/~tgkolda/TensorToolbox.
+% Copyright (2012) Sandia Corporation. Under the terms of Contract
+% DE-AC04-94AL85000, there is a non-exclusive license for use of this
+% work by or on behalf of the U.S. Government. Export of this data may
+% require a license from the United States Government.
+% The full license terms can be found in the file LICENSE.txt
+
+szx = size(X);
+N   = length(szx);
+Xn = permute(X,[n 1:n-1,n+1:N]);
+Xn = reshape(Xn,szx(n),prod(szx([1:n-1 n+1:N])));
+
+Y = Xn*Xn';
+
+[u,d] = eigs(Y, r, 'LM', struct('disp',0));
+
+flipsign = true;
+
+if flipsign
+    for j = 1:r
+        negidx = u(:,j)<0;
+        isNegNormGreater = norm(u(negidx,j),'fro') > norm(u(~negidx,j),'fro');
+        if isNegNormGreater
+            u(:,j) = -u(:,j);
+        end
+    end
+end
+end
+
+%-------------------------------------------------------------------------
+

codes/matricize.m

@@ -0,0 +1,31 @@
+function Tm = matricize(T,varargin) 
+%   Tm = TENMAT(T, RDIMS, CDIMS) creates a matrix representation of
+%   tensor T.  The dimensions specified in RDIMS map to the rows of
+%   the matrix, and the dimensions specified in CDIMS map to the
+%   columns, in the order given.
+% Modified tenmat function of Tensor Toolbox of Kolda et. al.
+% MATLAB Tensor Toolbox.
+% Copyright 2010, Sandia Corporation. 
+
+    tsize = size(T);
+    n = ndims(T);
+    if nargin==2
+        rdims=varargin{1};
+        tmp = true(1,n);
+        tmp(rdims) = false;
+        cdims = find(tmp);
+    else
+        cdims = varargin{1};
+        tmp = true(1,n); 
+        tmp(cdims) = false; 
+        rdims = find(tmp);
+    end
+
+    data = reshape(permute(T,[rdims cdims]), prod(tsize(rdims)), prod(tsize(cdims)));
+    Tm.tsize = tsize;
+    Tm.rindices = rdims;
+    Tm.cindices = cdims;
+    Tm.data = data;
+end
+
+