A multivariate quantile regression model with a factor structure is proposed to study data with many responses of interest. The factor structure is allowed to vary with the quantile levels, which makes our framework more flexible than the classical factor models. The model is estimated with the nuclear norm regularization in order to accommodate the high dimensionality of data, but the incurred optimization problem can only be efficiently solved in an approximate manner by off-the-shelf optimization methods. Such a scenario is often seen when the empirical risk is non-smooth or the numerical procedure involves expensive subroutines such as singular value decompo- sition. To ensure that the approximate estimator accurately estimates the model, non-asymptotic bounds on error of the the approximate estimator is established. For implementation, a numerical procedure that provably marginalizes the approximate error is proposed. The merits of our model and the proposed numerical procedures are demonstrated through Monte Carlo experiments and an application to finance involving a large pool of asset returns.