Density functional perturbation theory (DFPT) is a crucial tool for accurately describing lattice dynamics. The adaptively compressed polarizability (ACP) method reduces the computational complexity of DFPT calculations from O(N4) to O(N3) by combining the interpolative separable density fitting (ISDF) algorithm. However, the conventional QR factorization with column pivoting (QRCP) algorithm, used for selecting the interpolation points in ISDF, not only incurs a high cubic-scaling computational cost, O(N3), but also leads to suboptimal convergence. This convergence issue is particularly pronounced when considering the complex interplay between the external potential and atomic displacement in ACP-based DFPT calculations. Here, we present a machine learning K-means clustering algorithm to select the interpolation points in ISDF, which offers a more efficient quadratic-scaling O(N2) alternative to the computationally intensive cubic-scaling O(N3) QRCP algorithm. We implement this efficient K-means-based ISDF algorithm to accelerate plane-wave DFPT calculations in KSSOLV, which is a MATLAB toolbox for performing Kohn-Sham density functional theory calculations within plane waves. We demonstrate that this K-means algorithm not only offers comparable accuracy to QRCP in ISDF but also achieves better convergence for ACP-based DFPT calculations. In particular, K-means can remarkably reduce the computational cost of selecting the interpolation points by nearly 2 orders of magnitude compared to QRCP in ISDF.