Just use a backward stable algorithm to solve the least-squares problem $\min \|W^{1/2}(Ax-b)\|_2$ resulting from IRLS. For instance, use the QR factorization or the SVD of $W^{1/2}A$, instead of forming normal equations like you are doing here.

You are essentially using normal equations to solve the least-squares problem $\min \|W^{1/2}(Ax-b)\|_2$ resulting from IRLS. Normal equations are known not to be a backward stable algorithm. Use other standard algorithms for LS problems instead, like the QR factorization or the SVD of $W^{1/2}A$ instead. Those are guaranteed to be backward stable.