function [ Xopt, Vopt, max_acc, max_accval] = newnew( U, Y, R, W, W_, lambda, classnum, MAX_ITER, output_period, labeled_ids, unlabeled_ids )

    M = size( U, 1 );
    N = size( U, 2 );

    V = rand( M, R );
    X = rand( R, N );
    V = V ./ repmat( sum( V ), [ M, 1 ] );
    X = X ./ repmat( sum( X ), [ R, 1 ] );

    Unonzero = sparse( U > 0 );
    
    W_org_ = W_;
    D = diag( sum( W ) );
    L = D - W;
    D_ = diag( sum( W_ ) );
    L_ = D_ - W_;
    
    max_accval = zeros( length(labeled_ids), 1 );
    max_acc = zeros( length(labeled_ids), 1 );

    iter = 0;
    objective = Inf;
    last_residual = Inf;
    min_intra_dist = Inf;
    last_intra_dist = Inf;
    consec_inc = 0;

    while ( iter < MAX_ITER )
        
        dist = ml_sqrDist( X, X );
        D = diag( sum( W ) );
        L = D - W;
        W_ = W_org_;
        W_( W_ > 0 ) = W_( W_ > 0 ) ./ ( dist( W_ > 0 ) + 0.1 );
        D_ = diag( sum( W_ ) );
        L_ = D_ - W_;
        
        % update V
        last_normV_org = norm(V, 'fro');
        last_normV = last_normV_org;
        count = 0;
        while (1)
            count = count + 1;
            if N > 3000
                for i = 1 : floor((M+999)/1000)
                    range = (i-1)*1000+1:min(i*1000,M);
                    V(range,:) = V(range,:) .* ( ( U(range,:) ./ ( V(range,:) * X ) ) * X' );
                end
            else
                V = V .* ( ( U ./ ( V * X ) ) * X' );
            end
            V = V ./ repmat( sum( V ), [ M 1 ] );
            V( V == 0 ) = eps;
            V( V == 1 ) = 1 - eps;
            break;
           
            normV = norm(V, 'fro');
            if abs( last_normV - normV ) <= 0.0001 break;
            end
            last_normV = normV;
        end
        
%         Pw_d = mex_Pw_d(U,V,X);
%         [V,XX] = mex_EMstep(U,Pw_d,V,X);
        
        residual = compute_residual( U, V, X, Unonzero );        
        [ dist_ratio intra_dist inter_dist ] = compute_distance_ratio( X, L, L_, N );
        
        if iter > 0 && mod( iter, output_period ) == 0 && output_period > 0
            dist_ratio2 = sum(X(:) .^ 2);
            fprintf( 1, '%d : %f %f %f %f %f\n', iter, residual/N, dist_ratio, intra_dist, inter_dist, dist_ratio2 );
            for k = 1 : length( labeled_ids )
                [ acc1, acc2, acc3, acc4, accval ] = compute_accuracies( X, N, R, Y, labeled_ids{k}, unlabeled_ids{k}, classnum );
                YY = Y( labeled_ids{k} );
                WW = sparse([],[],[],N,N,0);
                for i = 1 : 20
                    range = labeled_ids{k}(YY == i);
                    WW(range,range) = 1;
                end
                DD = diag( sum( WW ) );
                LL = DD - WW;
                WW_ = sparse([],[],[],N,N,0);
                WW_(labeled_ids{k},labeled_ids{k}) = 1;
                DD_ = diag( sum( WW_ ) );
                LL_ = DD_ - WW_;
                intra_dist2 = trace( X * LL * X' );
                inter_dist2 = trace( X * LL_ * X' );
                dist_ratio2 = intra_dist2 / inter_dist2;
                
                fprintf( 1, ' labeld = %d, svm: %0.2f knn: %0.2f validation: %0.2f\n', length(labeled_ids{k})/classnum, acc2, acc4, accval );
            end
        end

        % update X     
        X2 = X;
        XL_num = X2 * W + dist_ratio * X2 * D_;
        XL_dem = X2 * D + dist_ratio * X2 * W_;

        alpha = XL_num ./ XL_dem;
        toobig = X2 .* alpha >= 1;
        alpha(toobig) = 1 ./ X2(toobig);
        beta = ( 1 - alpha .* X2 ) ./ ( 1 - X2 );
        beta( beta <= 0 ) = eps;

        X2 = X2 .* alpha ./ beta;
        X2 = X2 ./ repmat( sum(X2), [R 1] );
        X2( X2 <= 0 ) = eps;
        X2( X2 >= 1 ) = 1.0 - eps;
        
        last_residual = residual;
        last_dist_ratio = dist_ratio;
        last_intra_dist = intra_dist;

        X1 = X2;
        if N > 3000
            for i = 1 : floor((N+999)/1000)
                range = (i-1)*1000+1:min(i*1000,N);
                X1(:,range) = X1(:,range) .* ( V' * ( U(:,range) ./ ( V * X1(:,range) ) ) );
            end
        else
            X1 = X1 .* ( V' * ( U ./ ( V * X1 ) ) );
        end               
        
        diff = X1 - X2;
        
        low = 0.0;
        high = 1.0;
        while( high - low > 0.0001 )
            mid = (low + high) / 2;
            XX = X2 + mid * diff;
            residual_mid = compute_residual( U, V, XX, Unonzero );
            
            if residual_mid >= last_residual
                low = mid;                
            else
                [ dist_ratio intra_dist inter_dist ] = compute_distance_ratio( XX, L, L_, N );
                slope = 2*sum(sum( (inter_dist*XX*L - intra_dist*XX*L_)/inter_dist^2 .* diff ));
                if slope > 0 
                    high = mid;
                else
                    low = mid;
                end
            end
        end

        p = high;
        XX = X2 + p * diff;
        residual = compute_residual( U, V, XX, Unonzero );
        dist_ratio = compute_distance_ratio( XX, L, L_, N );
                
        if residual < last_residual && dist_ratio < last_dist_ratio
            X = X2 + diff * p;
        end

        X( X <= 0 ) = eps;
        X( X >= 1 ) = 1.0 - eps;
        
        iter = iter + 1;
    end

Xopt = X;
Vopt = V;