function AhatReg=VGnoCVversionRuben(X, Y, F, AhatIni)
    
METHOD = 'dual';%  'regression';%      % method to execute 'dual' or 'regression'
%% standardize data
        n = size(X,1);
        p = size(X,2);
        x = X;
        y = Y;

        %% normalize training data
        x = x-mean(x,2)*ones(1,p);
        dx = sqrt(1/p*sum(x.^2,2));
        x = x./(dx*ones(1,p));
        y = y-mean(y);
        
        sigmay = y*y'/p;  % sigma is output variance
    
%% other parameters (general)

    epsilon = 1e-4;         % determines gamma_min (smallest m, see Eq. 13)
    ng = F;                 % number of gammas (equivalent to step-size)
    eps_gamma = 1e-4;       % determines gamma_max
    ninits_vg = 1;          % number of random initializations
    dmmin = 1e-10;          % convergence of mean field error 1e-12
    eta0 = 1;               % initial step size of mean field update
    eta_w0 = 0.02;          % initial step size for change in w in dual.m
    maxiter = 1e3;          % maximal number of iterations 1e4
    max_sum_m = n;          % increases gamma until sum(m) == max_sum_m
    beta_max = sigmay;      % increases gamma until beta == beta_max

    
    min_total_val = inf;
    
 %% set gamma range (min, max and step size)
     b = x*y'/p;       % b is input output covariance
     chi_ii=1/p*sum(x.^2,2); % input data variance

     gamma_min = log(epsilon*sigmay/p/max(abs(b))); %gamma_min = -(p*max(abs(b))/sigmay) + log(epsilon/(1-epsilon));%
     gamma_max = eps_gamma*gamma_min;
     gamma_all =linspace(gamma_min, gamma_max, ng);


for ini = 1:ninits_vg    
    %% inits
    minklmin = inf(ng,1);
    v_mf_all = inf(ng,n);
    m_mf_all = inf(ng,n);
    iter_all = inf(ng,2);
    beta_mf_all = inf(1,ng);
    error_mf_all = inf(1,ng);
    errorv_mf_all = inf(1,ng);
    kl_all = inf(ng,2);
    v_all = inf(ng,2,n);
    m_all = inf(ng,2,n);
    beta_all = inf(ng,2);
    imin = inf(1,ng);

    %% Initialvalue
    if nargin>4
        m=AhatIni;
    else
        % random init m
        m = 0.5*ones(1,n) + 0.1*randn(1,n);
        m(m<=0.2) = 0.2;
        m(m>=0.8) = 0.8;
        eps2 = .4;
        while sum(m)>max_sum_m
            m = eps2*ones(1,n) + 0.1*randn(1,n);
            eps2 = eps2/2;
        end
    end

    %% the estimated inverse noise variance beta is initialized as the output variance
    beta = 1/sigmay;
    
        %% for gamma in gamma_min to gamma_max or when some criteria are satisfied
        i=0;
        while (beta<beta_max) && (i<ng) && (sum(m)<max_sum_m)
            i = i+1; 
            gamma = gamma_all(i);
            eval(METHOD);
            AhatReg(:,i) = v;
            m_all(i,1,:) = m;
            beta_all(i,1) = beta;
            iter_all(i,1) = iter;
            kl_all(i,1) = kl1;
        end
        
        %% for gamma is current gamma decreasing to gamma_min
        imax=i-1;
        for ii=imax:-1:1
            i = i+1;
            gamma = gamma_all(ii);
            eval(METHOD);
            AhatReg(:,i) = v;
            m_all(ii,2,:) = m;
            beta_all(ii,2) = beta;
            iter_all(ii,2) = iter;
            kl_all(ii,2) = kl1;
        end
end