tensorconst_adm

 TENSORCONST_ADM - Computes the reconstruction of a partly
                   observed tensor via "Constraint" approach

 Syntax
  function [X,Z,Y,fval,gval] = tensorconst_adm(X, I, Bv, lambda, eta, tol, verbose)

 Reference
 "On the extension of trace norm to tensors"
 Ryota Tomioka, Kohei Hayashi, and Hisashi Kashima
 arXiv:1010.0789
 http://arxiv.org/abs/1010.0789
 
 Copyright(c) 2010 Ryota Tomioka
 This software is distributed under the MIT license. See license.txt

0001 % TENSORCONST_ADM - Computes the reconstruction of a partly
0002 %                   observed tensor via "Constraint" approach
0003 %
0004 % Syntax
0005 %  function [X,Z,Y,fval,gval] = tensorconst_adm(X, I, Bv, lambda, eta, tol, verbose)
0006 %
0007 % Reference
0008 % "On the extension of trace norm to tensors"
0009 % Ryota Tomioka, Kohei Hayashi, and Hisashi Kashima
0010 % arXiv:1010.0789
0011 % http://arxiv.org/abs/1010.0789
0012 %
0013 % Copyright(c) 2010 Ryota Tomioka
0014 % This software is distributed under the MIT license. See license.txt
0015 
0016 
0017 function [X,Z,Y,fval,gval] = tensorconst_adm(X, I, Bv, lambda, eta, tol, verbose)
0018 
0019 if ~exist('tol','var')
0020   tol=1e-3;
0021 end
0022 
0023 if ~exist('verbose','var')
0024   verbose=0;
0025 end
0026 
0027 sz=size(X);
0028 nd=ndims(X);
0029 m =length(I{1});
0030 
0031 if nd~=length(I)
0032   error('Number of dimensions mismatch.');
0033 end
0034 
0035 if m~=length(Bv)
0036   error('Number of samples mismatch.');
0037 end
0038 
0039 Z=cell(1,nd);
0040 Y=cell(1,nd);
0041 S=cell(1,nd);
0042 
0043 for jj=1:nd
0044   szj = [sz(jj), prod(sz)/sz(jj)];
0045   Y{jj} = zeros(szj);
0046   Z{jj} = zeros(szj);
0047 end
0048 
0049 B=zeros(sz);
0050 ind=sub2ind(sz, I{:});
0051 B(ind)=Bv;
0052 
0053 nsv=10*ones(1,nd);
0054 kk=1;
0055 while 1
0056   X1 = zeros(size(X));
0057   for jj=1:nd
0058     X1 = X1 - flatten_adj(Y{jj}-eta*Z{jj},sz,jj);
0059   end
0060   
0061   if lambda>0
0062     X1(ind) = X1(ind) + Bv/lambda;
0063     X=X1./((B~=0)/lambda + nd*eta);
0064   else
0065     X=X1/(eta*nd);
0066     X(ind)=Bv;
0067   end
0068   
0069 
0070   % Check derivative
0071 % $$$   D=zeros(size(X));
0072 % $$$   D(ind)=X(ind)-Bv;
0073 % $$$   for jj=1:nd
0074 % $$$     D=D+eta*flatten_adj(Y{jj}/eta+flatten(X,jj)-Z{jj},sz,jj);
0075 % $$$   end
0076 % $$$   fprintf('gnorm=%g\n',norm(D(:)));
0077   
0078   for jj=1:nd
0079     [Z{jj},S{jj},nsv(jj)] = softth(Y{jj}/eta+flatten(X,jj),1/eta,nsv(jj));
0080     
0081     % Check derivative
0082     % fprintf('max[%d]=%g\n',jj,max(svd(eta*(Z{jj}-flatten(X,jj)-Y{jj}/eta))));
0083   end
0084 
0085   for jj=1:nd
0086     V=flatten(X,jj)-Z{jj};
0087     Y{jj}=Y{jj}+eta*V;
0088     viol(jj)=norm(V(:));
0089   end
0090   
0091   
0092   % Compute the objective
0093   G=zeros(size(X));
0094   fval(kk)=0;
0095   for jj=1:nd
0096     fval(kk)=fval(kk)+sum(svd(flatten(X,jj)));
0097     G = G + flatten_adj(Y{jj},sz,jj);
0098   end
0099   if lambda>0
0100     fval(kk)=fval(kk)+0.5*sum((X(ind)-Bv).^2)/lambda;
0101     G(ind)=G(ind)+(X(ind)-Bv)/lambda;
0102   else
0103     G(ind)=0;
0104   end
0105   
0106   gval(kk)=norm(G(:));
0107 
0108   if verbose
0109     fprintf('k=%d fval=%g gnorm=%g viol=%s\n',...
0110             kk, fval(kk), gval(kk), printvec(viol));
0111   end
0112   
0113   if kk>1 && max(viol)<tol && gval(kk)<tol
0114     break;
0115   end
0116   
0117   if kk>1000
0118     break;
0119   end
0120   
0121   
0122   kk=kk+1;
0123 end
0124 
0125 fprintf('k=%d fval=%g gnorm=%g viol=%s\n',...
0126         kk, fval(kk), gval(kk), printvec(viol));
0127