
    elseif size(G,1)~=size(A,2)
        disp('G has wrong dimensions.')
        return
    end
elseif (size(x,1)~=size(A,2))||(size(x,2)~=size(Y,2))
    disp('x has wrong dimensions.')
    return
end
X = x;
if ~isempty(opts.gpu_id)
    gpuDevice(opts.gpu_id);
    x = gpuArray(x);
    h = gpuArray(h);
    G= gpuArray(G);
end
rds=1:size(Y,2);

%% Compute the second term (the Gramian matrix) that does not change throughout the function
if opts.use_std
    G = @(x) G*x - (G*sum(x,2))/size(Y,2);
    h = h - A'*