A novel machine learning based mixed integer programming model is developed for the optimal nonparametric prediction intervals (PIs) of electricity load, which minimizes interval width subject to target hit probability constraint. Binary variables are employed to analytically formulate PI coverage states and further impose calibration constraint. Considering the quantile interpretation on lower and upper bounds of PIs, an innovative binary variable reduction strategy is proposed to significantly accelerate model training process. Compared with traditional central PIs, the resultant optimal nonparametric PIs are featured with adaptive symmetric and asymmetric quantile proportion pairs for interval shortening. Numerical experiments based on actual substation data validate the remarkable superiority of the presented method with respect to both PI quality and computational efficiency.