  Variational Physics-Informed Neural Networks often suffer from poor
convergence when using stochastic gradient-descent-based optimizers. By
introducing a Least Squares solver for the weights of the last layer of the
neural network, we improve the convergence of the loss during training in most
practical scenarios. This work analyzes the computational cost of the resulting
hybrid Least-Squares/Gradient-Descent optimizer and explains how to implement
it efficiently. In particular, we show that a traditional implementation based
on backward-mode automatic differentiation leads to a prohibitively expensive
algorithm. To remedy this, we propose using either forward-mode automatic
differentiation or an ultraweak-type scheme that avoids the differentiation of
trial functions in the discrete weak formulation. The proposed alternatives are
up to 100 times faster than the traditional one, recovering a computational
cost-per-iteration similar to that of a conventional gradient-descent-based
optimizer alone. To support our analysis, we derive computational estimates and
conduct numerical experiments in one- and two-dimensional problems.
