We propose a subspace constrained precoding scheme that exploits the spatial channel correlation structure in massive MIMO cellular systems to fully unleash the tremendous gain provided by massive antenna array with reduced channel state information (CSI) signaling overhead. The MIMO precoder at each base station (BS) is partitioned into an inner precoder and a Transmit (Tx) subspace control matrix. The inner precoder is adaptive to the local CSI at each BS for spatial multiplexing gain. The Tx subspace control is adaptive to the channel statistics for inter-cell interference mitigation and Quality of Service (QoS) optimization. Specifically, the Tx subspace control is formulated as a QoS optimization problem which involves an SINR chance constraint where the probability of each user's SINR not satisfying a service requirement must not exceed a given outage probability. Such chance constraint cannot be handled by the existing methods due to the two stage precoding structure. To tackle this, we propose a bi-convex approximation approach, which consists of three key ingredients: random matrix theory, chance constrained optimization and semidefinite relaxation. Then we propose an efficient algorithm to find the optimal solution of the resulting bi-convex approximation problem. Simulations show that the proposed design has significant gain over various baselines.