The objective of this study is to incorporate the proper generalized decomposition (PGD) into the har-monic balance method (HBM) for an efficient numerical continuation of nonlinear frequency response. PGD utilizes the low-dimensional subspaces of the HBM response, approximating the solution as a low-rank separated representation of spatial and harmonic components. The progressive Galerkin approach is employed to formulate subproblems for each component, and the solution is built on the fly without prior assumption or computation of nonlinear response characteristics. During the numerical continuation, the spatial modes obtained at the previous computation are utilized as a reduced basis, and arc-length constraint is applied to follow the proper path of the frequency response curve. Numerical studies demonstrate that the proposed method allows significant computational savings compared to HBM, while accurately reflecting the nonlinear behavior.(c) 2022 Elsevier Ltd. All rights reserved.