function [sol, info] = proxlinblind(U, V, b, gamma, beta, init_z, maxiter, tol, ...
    early_stop, alpha_0, show_info)

% Get problem size
[m, n] = size(U);

xidx = 1:n;
yidx = (n + 1) : (2 * n);

% Get initial point
z_before = init_z;
z_after = z_before;
bestz = init_z;

% Initialize arrays for storing results
objs = zeros(m * maxiter + 1, 1);
bestobjs = zeros(m * maxiter + 1, 1);
obj = getblObj(U, V, b, z_before(xidx), z_before(yidx));
bestobj = obj;
objs = objs + obj;
bestobjs = bestobjs + bestobj;

% Number of epochs before reaching tolerance
nepochs = maxiter;
nbatchiter = maxiter * m;

% Initialize info struct
info.status = "x";

if show_info
    fprintf("%6s %10s %10s %10s\n", 'epoch', 'obj', 'bobj', 'status');
end % End if

idx = 0; 
for k = 1:maxiter
    
    if bestobj < tol && nepochs == maxiter
        nepochs = k;
        niter = k * m + idx;
        info.status = "*";
        if early_stop
            if show_info
                fprintf("%6d %10.2e %10.2e %10s\n", k, obj, bestobj, info.status);
                fprintf("Early stopped. Status: Optimal\n");
                 bestobjs((k - 1) * m + 1:end) = bestobj;
                 objs((k - 1) * m + 1:end) = obj;
            end % End if
            break;
        end
    end % End if
    
    idx = 0;
    niter = m;
    
    for i = randperm(niter) % for i = randsample(1:m, m, true)
        
        idx = idx + 1;
        
        gamma = gamma / alpha_0;
        
        % Sample from dataset
        u = U(i, :);
        v = V(i, :);
                
        % Update momentum
        w = (1 + beta) * z_after - beta * z_before;
        
        uTx = u * z_after(xidx);
        vTy = v * z_after(yidx);
        z_before = z_after;
        zeta = [vTy * u'; uTx * v'];
        ksi = uTx * vTy + zeta' * (w - z_before) - b(i);
        coef = - gamma * ksi / norm(zeta)^2;
        coef = sign(coef) * min(1, abs(coef));
        z_after = w + zeta / gamma * coef;
        obj = getblObj(U, V, b, z_after(xidx), z_after(yidx));
        
        gamma = gamma * alpha_0;
        
        if obj < bestobj
            bestobj = obj;
            bestz = z_after;
        end % End if
        
        bestobjs((k * m - m) + idx + 1) = bestobj;
        objs((k * m - m) + idx + 1) = obj;
       
        if isnan(obj) || isinf(obj)
            info.status = "Diverged";
            break;
        end % End if
        
    end % End for
    
    if isnan(obj) || isinf(obj)
        break;
    end % End if
    
    if show_info && (mod(k, 50) == 0 || k == 1)
        fprintf("%6d %10.2e %10.2e %10s\n", k, obj, bestobj, info.status);
    end % End if
    
end % End for

% Collect information

% Solution array
sol.x = z_after;
sol.bestz = bestz;

% Information array
info.nepochs = nepochs;
info.niter = niter;
info.objs = objs;
info.bestobjs = bestobjs;

% Display summary
if show_info
    if info.status == "*"
        disp("- Algorithm reaches optimal after " + nepochs + " epochs (" + ...
            niter + " iterations)");
    elseif info.status == "x"
        disp("- Algorithm fails to reach desired accuracy after " +...
            nepochs + " epochs");
        info.niter = maxiter * m;
    elseif info.status == "Diverged"
        disp("- Algorithm diverges");
        info.niter = maxiter * m;
    end % End if 
end % End if

end % End function