This paper considers the problem of minimizing a differentiable function with locally Lipschitz continuous gradient on the algebraic variety of all m-by-n real matrices of rank at most r. Several definitions of stationarity exist for this nonconvex problem. Among them, Bouligand stationarity is the strongest necessary condition for local optimality. Only a handful of algorithms generate a sequence in the variety whose accumulation points are provably Bouligand stationary. Among them, the most parsimonious with (truncated) singular value decompositions (SVDs) or eigenvalue decompositions can still require a truncated SVD of a matrix whose rank can be as large as min{m, n}-r+1 if the gradient does not have low rank, which is computationally prohibitive in the typical case where r << min{m, n}. This paper proposes a first-order algorithm that generates a sequence in the variety whose accumulation points are Bouligand stationary while requiring SVDs of matrices whose smaller dimension is always at most r. A standard measure of Bouligand stationarity converges to zero along the bounded subsequences at a rate at least O(1/sqrt(i+1)), where i is the iteration counter. Furthermore, a rank-increasing scheme based on the proposed algorithm is presented, which can be of interest if the parameter r is potentially overestimated.