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(N 4 ) to O(N 3 ) 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(N 3 ), 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(N 2 ) alternative to the computationally intensive cubic-scaling O(N 3 ) 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.