increaseRank_mod

   Rank-1 gradient approximation to increase the rank.
   Unfortunately, not really worth the effort,
   and a random rank-1 increase works equally well.

   TTeMPS Toolbox. 
   Michael Steinlechner, 2013-2016
   Questions and contact: michael.steinlechner@epfl.ch
   BSD 2-clause license, see LICENSE.txt

0001 %   Rank-1 gradient approximation to increase the rank.
0002 %   Unfortunately, not really worth the effort,
0003 %   and a random rank-1 increase works equally well.
0004 %
0005 %   TTeMPS Toolbox.
0006 %   Michael Steinlechner, 2013-2016
0007 %   Questions and contact: michael.steinlechner@epfl.ch
0008 %   BSD 2-clause license, see LICENSE.txt
0009 function X = increaseRank_mod( X, A_Omega, Omega, idx )
0010 
0011     r = X.rank;
0012     n = X.size;
0013     d = X.order;
0014 
0015     epsilon = 1;
0016     
0017     rankInc = 2;
0018     
0019     %idx = d-1;
0020     idxL = idx;
0021     idxR = idx + 1;
0022 
0023     X = orthogonalize(X, idxR);
0024 
0025     %Omega = sub2ind( n, deal(mat2cell(Omega, size(Omega,1), ones(1,d))) );
0026     Omega_ind = sub2ind( n, Omega(:,1), Omega(:,2), Omega(:,3), Omega(:,4), Omega(:,5) );
0027     Z = zeros( n );
0028     Z(Omega_ind) =  X(Omega) - A_Omega;
0029     
0030     % right side
0031     for i = d-1:-1:idxR
0032         zz = reshape( Z, [prod(n(1:i)), n(i+1)*r(i+2)] );
0033         xx = transpose( unfold( X.U{i+1}, 'right') );
0034         Z = zz*xx;
0035     end
0036     
0037     % left side
0038     for i = 1:idxL-1
0039         Z_i = reshape( Z, [r(i)*prod(n(i)), prod(n(i+1:idxR))*r(idxR+1)] );
0040         X_i = unfold( X.U{i}, 'left'); 
0041         Z = X_i' * Z_i;
0042     end
0043     size(Z)
0044     Z = reshape( Z, [r(idxL)*n(idxL), n(idxR)*r(idxR+1)] );
0045 
0046     % truncate to rank-1
0047     
0048     norm_Z = norm(Z(:))
0049     
0050     [U,S,V] = svd( Z, 'econ');
0051     Z1 = U(:,1:rankInc)*S(1:rankInc,1:rankInc);
0052     Z2 = V(:,1:rankInc)';
0053     
0054     Z1 = reshape( Z1, [r(idxL), n(idxL), rankInc] );
0055     Z2 = reshape( Z2, [rankInc, n(idxR), r(idxR+1)]);
0056     
0057     
0058     epsilon = fminbnd( @(t) linesearch(X, t, idxL, idxR, Z1, Z2, Omega, A_Omega), -100, 1)
0059     epsilon2 = linesearch2(X, idxL, idxR, Z1, Z2, Omega, A_Omega)
0060     
0061     X.U{idxL} = cat( 3, X.U{idxL}, epsilon2*Z1 );
0062     X.U{idxR} = cat( 1, X.U{idxR}, Z2 );
0063     
0064 end
0065 
0066 
0067 function res = linesearch( X, t, idxL, idxR, Z1, Z2, Omega, A_Omega )
0068     X.U{idxL} = cat( 3, X.U{idxL}, t*Z1 );
0069     X.U{idxR} = cat( 1, X.U{idxR}, Z2 );
0070     
0071     res = 0.5*norm( A_Omega - X(Omega) )^2;
0072 end
0073 
0074 function res = linesearch2( X, idxL, idxR, Z1, Z2, Omega, A_Omega )
0075     Y = X;
0076     Y.U{idxL} = Z1;
0077     Y.U{idxR} = Z2;
0078     
0079     res = Y(Omega)'*(A_Omega - X(Omega)) / norm(Y(Omega))^2;
0080 end