tensormix_adm

 TENSORMIX_ADM - Computes the reconstruction of a partly 
                 observed tensor via "Mixture" approach

 Syntax
  function [X,Z,fval,gval]=tensormix_adm(X, I, yy, 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 % TENSORMIX_ADM - Computes the reconstruction of a partly
0002 %                 observed tensor via "Mixture" approach
0003 %
0004 % Syntax
0005 %  function [X,Z,fval,gval]=tensormix_adm(X, I, yy, 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,fval,gval]=tensormix_adm(X, I, yy, lambda, eta, tol, verbose)
0018 
0019 if ~exist('tol','var')
0020   tol=1e-3;
0021 end
0022 
0023 if ~exist('eta','var')
0024   eta=1;
0025 end
0026 
0027 if ~exist('verbose','var')
0028   verbose=0;
0029 end
0030 
0031 sz=size(X);
0032 nd=ndims(X);
0033 ind=sub2ind(sz, I{:});
0034 ind0=setdiff(1:prod(sz),ind);
0035 m=size(yy,1);
0036 
0037 Z=cell(1,nd);
0038 V=cell(1,nd);
0039 for jj=1:nd
0040   szj = [sz(jj), prod(sz)/sz(jj)];
0041   Z{jj} = zeros(szj);
0042   V{jj} = zeros(szj);
0043 end
0044 
0045 nsv=10*ones(1,nd);
0046 kk=1;
0047 
0048 alpha=yy;
0049 A=zeros(sz); A(ind)=alpha;
0050 while 1
0051   for jj=1:nd
0052     Ztmp = Z{jj}+eta*flatten(A,jj);
0053     [Z{jj},S{jj},nsv(jj)]=softth(Ztmp,eta,nsv(jj));
0054     V{jj}=(Ztmp-Z{jj})/eta;
0055 
0056     viol(jj)=norm(flatten(A,jj)-V{jj},'fro');
0057   end
0058  
0059   Zsum = zeros(sz);
0060   Vsum = zeros(sz);
0061   for jj=1:nd
0062     Zsum=Zsum+flatten_adj(Z{jj},sz,jj);
0063     Vsum=Vsum+flatten_adj(V{jj},sz,jj);
0064   end
0065   
0066   alpha = (yy-Zsum(ind)+eta*Vsum(ind))/(lambda+eta*nd);
0067   A(ind)=alpha;
0068 
0069    
0070   % Compute the objective
0071   if lambda>0
0072     fval(kk)=0.5*sum((Zsum(ind)-yy).^2)/lambda;
0073   else
0074     fval(kk)=0;
0075   end
0076   
0077   gval(kk)=norm(lambda*alpha - yy + Zsum(ind));
0078 
0079   for jj=1:nd
0080     fval(kk)=fval(kk)+sum(S{jj});
0081   end
0082   
0083   % Compute the dual objective
0084 % $$$   fact=1;
0085 % $$$   for jj=1:nd
0086 % $$$     fact=min(fact, 1/norm(flatten(A,jj)));
0087 % $$$   end
0088 % $$$   aa = alpha*fact;
0089 % $$$
0090 % $$$   if kk>1
0091 % $$$     dval(kk)=max(dval(kk-1),-0.5*lambda*sum(aa.^2)+aa'*yy);
0092 % $$$   else
0093 % $$$     dval(kk)=-inf;
0094 % $$$   end
0095   
0096   if verbose 
0097     fprintf('[%d] fval=%g gval=%g viol=%s\n', kk, fval(kk), ...
0098           gval(kk), printvec(viol));
0099 %    fprintf('[%d] fval=%g dval=%g fact=%g\n', kk, fval(kk), ...
0100 %          dval(kk), fact);
0101   end
0102 %  if kk>1 && 1-dval(kk)/fval(kk)<tol
0103   if kk>1 && max(viol)<tol && gval(kk)<tol
0104    break;
0105   end
0106   
0107   if kk>1000
0108     break;
0109   end
0110 
0111   kk=kk+1;
0112 end
0113 
0114 fprintf('[%d] fval=%g gval=%g viol=%s\n', kk, fval(kk), ...
0115         gval(kk), printvec(viol));
0116 
0117 %fprintf('[%d] fval=%g dval=%g fact=%g\n', kk, fval(kk), ...
0118 %        dval(kk), fact);
0119 
0120 X=zeros(sz);
0121 for jj=1:nd
0122   X = X + flatten_adj(Z{jj},sz,jj);
0123   Z{jj}=Z{jj}*nd;
0124 end