Problem setting :

$ \underset{x}{\text{min}}  \|Ax-b\|$, where $A \in \mathcal{R}^{m \times n}, m\gg n $, full rank.

L1 loss is used for robust estimation using IRLS. The corresponding equation to solve turns out to be $ A^{T}WAx=A^{T}Wb$, where $W=\mathrm{diag}(d_i), d_i=1/|e_{i}|$, $e_{i}=a_{i}^{T}x-b_{i}$, $a_{i}$ is the ith row of $A$, $b_{i}$ is the ith element of $b$. For $e_{i}$ close to $0$, the value of $d_{i}$ is very large. For my specific case, the range of $d_i$ is from $10^{-3}$ to $10^5$.

To avoid high values of $d_i$, it is taken as $d_i=1/(|e_i|+\delta)$ where $\delta>0$ is a small number near $0$. Let $\delta=10^{-3}$. This brings the range of $d_i$ as $10^{-3}$ to $10^3$. The range of values of $d_i$ is still high to bring numerical stability. It makes 
$ A^{T}WA$ a near singular matrix. 

Please suggest a way to avoid numerical instability. 

Thanks in advance!