Locally Weighted Regression (Loess) - Robustifying iterations

Question

As a self study exercise I am trying to understand the implementation of locally weighted regression (Loess) in python. Alexandre Gramfort (Sklearn developper) provides the following code on his github page.

def lowess(x, y, f=2. / 3., iter=3):
    '''....the number of robustifying iterations is given by iter. The
    function will run faster with a smaller number of iterations'''

    n = len(x)
    r = int(ceil(f * n))
    h = [np.sort(np.abs(x - x[i]))[r] for i in range(n)]
    w = np.clip(np.abs((x[:, None] - x[None, :]) / h), 0.0, 1.0)
    w = (1 - w ** 3) ** 3
    yest = np.zeros(n)
    delta = np.ones(n)
    for iteration in range(iter):
        for i in range(n):
            weights = delta * w[:, i]
            b = np.array([np.sum(weights * y), np.sum(weights * y * x)])
            A = np.array([[np.sum(weights), np.sum(weights * x)],
                          [np.sum(weights * x), np.sum(weights * x * x)]])
            beta = linalg.solve(A, b)
            yest[i] = beta[0] + beta[1] * x[i]

        residuals = y - yest
        s = np.median(np.abs(residuals))
        delta = np.clip(residuals / (6.0 * s), -1, 1)
        delta = (1 - delta ** 2) ** 2

    return yest

Question

What I don't understand is the part related to the 'robustifying iterations'. How does this work and how could I learn more about this ?

Vectorized implementation

For those who might be interested here is my attempt at describing the vectorized implementation mathematically.

Consider the 1D case where $\Theta = [\theta_0, \theta_1]$ and $x$ and $y$ are vectors of size $m$. The cost function $J(\theta)$ is a weighted version of the OLS regression, where the weights $w$ are defined by some kernel function

\begin{aligned} J(\theta) &= \sum_{i=1}^m w^{(i)} \left( y^{(i)} - (\theta_0 + \theta_1 x^{(i)}) \right)^2 \\ \frac{\partial J}{\partial \theta_0} &= -2 \sum_{i=1}^m w^{(i)} \left( y^{(i)} - (\theta_0 + \theta_1 x^{(i)}) \right) \\ \frac{\partial J}{\partial \theta_1} &= -2 \sum_{i=1}^m w^{(i)} \left( y^{(i)} - (\theta_0 + \theta_1 x^{(i)}) \right) x^{(i)} \end{aligned}

Cancelling the $-2$ terms, equating to zero, expanding and re-arranging the terms: \begin{aligned} & \frac{\partial J}{\partial \theta_0} = \sum_{i=1}^m w^{(i)} \left( y^{(i)} - (\theta_0 + \theta_1 x^{(i)}) \right) = 0 \\ & \sum_{i=1}^m w^{(i)} \theta_0 + \sum_{i=1}^m w^{(i)} \theta_1 x^{(i)} = \sum_{i=1}^m w^{(i)} y^{(i)} &\text{Eq. (1)} \\ \\ & \frac{\partial J}{\partial \theta_1} = \sum_{i=1}^m w^{(i)} \left( y^{(i)} - (\theta_0 + \theta_1 x^{(i)}) \right) x^{(i)} = 0 \\ & \sum_{i=1}^m w^{(i)} \theta_0 + \sum_{i=1}^m w^{(i)} \theta_1 x^{(i)} x^{(i)} = \sum_{i=1}^m w^{(i)} y^{(i)} x^{(i)} &\text{Eq. (2)} \end{aligned}

Writing Eq. (1) and Eq. (2) in matrix form $\mathbf{A \Theta = b}$ allows us to solve for $\Theta$ \begin{aligned} & \sum_{i=1}^m w^{(i)} \theta_0 + \sum_{i=1}^m w^{(i)} \theta_1 x^{(i)} = \sum_{i=1}^m w^{(i)} y^{(i)} \\ & \sum_{i=1}^m w^{(i)} \theta_0 + \sum_{i=1}^m w^{(i)} \theta_1 x^{(i)} x^{(i)} = \sum_{i=1}^m w^{(i)} y^{(i)} x^{(i)} \\ & \begin{bmatrix} \sum w^{(i)} & \sum w^{(i)} x^{(i)} \\ \sum w^{(i)} x^{(i)} & \sum w^{(i)} x^{(i)} x^{(i)} \end{bmatrix} \begin{bmatrix} \theta_0 \\ \theta_1 \end{bmatrix} = \begin{bmatrix} \sum w^{(i)} y^{(i)} \\ \sum w^{(i)} y^{(i)} x^{(i)} \end{bmatrix} \\ & \mathbf{A} \Theta = \mathbf{b} \\ & \Theta = \mathbf{A}^{-1} \mathbf{b} \end{aligned}