A multivariate quantile regression model with a factor structure is proposed to study data with multivariate responses with covariates. The factor structure is allowed to vary with the quantile levels, which is more flexible than the classical factor models. Assuming the number of factors is small, and the number of responses and the input variables are growing with the sample size, the model is estimated with the nuclear norm regularization. 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 loss is nonsmooth or the numerical procedure involves expensive subroutines, for example, singular value decomposition. To show that the approximate estimator is still statistically accurate, we establish a nonasymptotic bound on the Frobenius risk and prediction risk. For implementation, a numerical procedure that provably marginalizes the approximation error is proposed. The merits of our model and the proposed numerical procedures are demonstrated through the Monte Carlo simulation and an application to finance involving a large pool of asset returns.