Boiling water reactor (BWR) loading pattern assessment requires solving the two-group, nodal form of the neutron diffusion equation and drift-flux form of the fluid equations simultaneously because these equation sets are strongly coupled via nonlinear feedback. To reduce the computational burden associated with the calculation of the core attributes (that is, core eigenvalue and thermal margins) of a perturbed BWR loading pattern, the analytical and numerical aspects of a higher order generalized perturbation theory (GPT) method, which correctly addresses the strong nonlinear feedbacks of two-phase flow, have been established. Inclusion of Jacobian information in the definition of the generalized flux adjoints provides for a rapidly convergent iterative method for solution of the power distribution and eigenvalue of a loading pattern perturbed from a reference state. Results show that the computational speedup of GPT compared with conventional forward solution methods demanding consistent accuracy is highly dependent on the number of spatial nodes utilized by the core simulator, varying from superior to inferior performance as the number of nodes increases.