% example script %% Generate synthetic data D_true = 5; N = [1000 50 25]; % Tensor dimensions Nx = length(N); F = cell(Nx,1); for i = 1:Nx F{i} = rand(N(i),D_true); end % Diagonal identity tensor I=zeros(D_true*ones(1,Nx)); for j=1:D_true I(j,j,j)=1; end Y=tmult(I,F{1},1); % Data tensor for ip = 2:Nx Y=tmult(Y,F{ip},ip); end sig2 = 0.5; % noise level C = 5; % affine transformation to ensure non-negatitivty X = Y + sqrt(sig2)*randn(N) + C*ones(N); assert(min(X(:))>0); %% Holdout missing data p = 0.20; % holdout fraction (missing data) NE = prod(size(X)); % number of elements in tensor R = rand(NE,1)>(1-p); % holdout logical indices X(R) = nan; % missing values are treated as NaN %% Model specification D = 5; % number of latent componenents in the model Finit = cell(Nx,1); % initialization of factors (default) scale = nanstd(X(:)); % scale of data for i = 1:Nx Finit{i}=(scale.^(1/Nx))*rand(N(i),D); end % options options.maxiter = 250; % number of iterations options.mu = 0; % no multiplicative update steps are taken options.hals = 1; % hierarchical alternating least sqaures steps are taken %% Run [FACT,SSEv,CPUt]=cpNonNeg(X,D,Finit,options); % FACT gives back factors in a cell array just as Finit was initialized % SSEv is the sum of squared (reconstruction) errors in each iteration % (vector) % CPUt is the CPU time used in each iteration