Physics-Informed Neural Networks (PINNs) have recently emerged as a powerful tool for solving forward and inverse problems involving partial differential equations. They offer a meshless alternative to standard discretization-based numerical solvers and enable seamless integration of available measurements for complex data assimilation tasks. Since their appearance, PINNs have sparked significant interest across the scientific machine learning community, with many authors either showing promising applications or proposing techniques to improve their predictive accuracy. However, both accuracy and robustness remain challenging: PINNs frequently fail to converge, even in seemingly simple scenarios, and even when successful, they often struggle to achieve high-precision results. Furthermore, a comprehensive mathematical theory able to explain their failure modes is lacking, and this has led many authors to propose heuristic tricks that fail to bridge the performance gap between PINNs and established numerical techniques like finite volume/element methods. In this work, we analyze the training dynamics of PINNs through the lens of the Neural Tangent Kernel (NTK) theory. From it, we derive an optimal residuals-weighting technique and demonstrate that gradient descent, when applied to the resulting modified loss function, is equivalent to the Gram-Gauss-Newton method. We make explicit the connection with the standard Gauss-Newton method, remarking that our algorithm is advantageous when training overparameterized neural networks, as computations happen in residuals space rather than parameters space. The exact application of the method would require the computation and inversion of the full NTK matrix. In order to enhance both the scalability and numerical stability of the training process, we develop an approximate procedure to compute the residuals weights based on block-Jacobi preconditioning of the NTK. The method is benchmarked on non-trivial fluid dynamics applications: Kovasznay and Beltrami flows, Taylor-Green Vortex, and Cavity flow at several Reynolds numbers. We show it consistently obtains state-of-the-art accuracies at competitive computational cost.