The design of variable fractional delay (VFD) FIR filters in the weighted least squares (WLS) sense is investigated in this paper. Unlike most methods in the literature where separable weighting functions are used, we consider the more general case with an arbitrary nonnegative weighting function. To solve this problem more efficiently, the WLS design of VFD FIR filters is formulated using matrix notation, and its optimality condition is established as a matrix equation. Then, a matrix iterative algorithm is presented to solve this matrix equation, leading to a great saving in computations and memory space compared to existing methods. Finally, an example is provided to show the good performance of the proposed algorithm.