While this does not answer the question asked (and will therefore not be accepted), I provide this response in hopes it will help others and promote worthwhile discussion.

David Baird provides a simple explanation for propagating the error through a linear least squares fit in his book Experimentation: An Introduction to Measurement Theory and Experiment Design.

1. Essentially, weight the fit by dividing $x_i$ and $y_i$ values by the corresponding $dy_i$. As commented below, when reducing this to a matrix problem as described in the OP's wiki-link, make sure to also weight the $1$, i.e., $X=[x, 1]/dy$.
2. Perform the fit to get the parameters $m$ and $b$ in $$y=mx+b$$
3. Calculate the differences of each $y_i$ value from the fit, $\delta y_i$
4. Calculate the standard deviation of the fit parameters using $$\sigma_y = \sqrt{\frac{\sum(\delta y_i)^2}{n-2}}$$ $$\sigma_m = \sigma_y \sqrt{\frac{n}{n\sum{x_i^2}-\left(\sum{x_i}\right)^2}}$$ $$\sigma_b = \sigma_y \sqrt{\frac{\sum{x_i^2}}{n\sum{x_i^2}-\left(\sum{x_i}\right)^2}}$$

I am not certain this is a complete propagation of the error but it seems reasonable.

While this does not answer the question asked (and will therefore not be accepted), I provide this response in hopes it will help others and promote worthwhile discussion.

In his book Experimentation: An Introduction to Measurement Theory and Experiment Design, David Baird provides a simple explanation for doing a linear least squares fit using the diffences between the measured and fit values to estimate the error of the fit parameters. The best fit for the parameters $m$ and $b$ in $$y=mx+b$$ will is determined using his eqn. (6.3):

$$m = \frac{n\sum x_i y_i - \sum x_i \sum y_i}{n\sum(x^2_i) - (\sum x_i)^2} $$

$$b = \frac{\sum(x_i^2)\sum y_i - \sum x_i \sum x_iy_i}{n\sum(x^2_i) - (\sum x_i)^2} $$

After obtaining $m$ and $b$, a standard deviation for the fit parameters can be obtained by calculating the differences of each $y_i$ value from the fit, $\delta y_i = y_i - (m x_i +b)$. From these differences, one calculates the standard deviation of the data from the fit line using:

$$\sigma_y = \sqrt{\frac{\sum(\delta y_i)^2}{n-2}}$$

and then the standard deviation of the parameters using $$\sigma_m = \sigma_y \sqrt{\frac{n}{n\sum{x_i^2}-\left(\sum{x_i}\right)^2}}$$ and $$\sigma_b = \sigma_y \sqrt{\frac{\sum{x_i^2}}{n\sum{x_i^2}-\left(\sum{x_i}\right)^2}}$$

To include the measurement error, $dy_i$, in the fit one would divide the initial system of equations by $dy_i$ giving $$\frac{y_i}{dy_i} = m\frac{x_i}{y_i} + b \frac{1}{dy_i}$$ then repeat Baird's derivation to get the weighted fit parameters

$$m = \frac{\sum \frac{1}{dy_i}\sum \frac{x_i y_i}{dy_i^2} - \sum \frac{x_i}{dy_i^2} \sum \frac{y_i}{dy_i}}{\sum \frac{1}{dy_i}\sum \frac{x^2_i}{dy_i^2} - \sum \frac{x_i}{dy_i}\sum \frac{x_i}{dy_i^2}} $$

$$b = \frac{\sum \frac{x_i^2}{dy_i^2}\sum \frac{y_i}{dy_i} - \sum \frac{x_i}{dy_i} \sum \frac{x_iy_i}{dy_i^2}}{\sum \frac{1}{dy_i}\sum \frac{x^2_i}{dy_i^2} - \sum \frac{x_i}{dy_i}\sum \frac{x_i}{dy_i^2}} $$

Notice that the $b\frac{1}{dy_i}$ term makes it so you cannot simply divide $x_i$ and $y_i$ by $dy_i$ (as pointed out in the comments below).

Unfortunately, this does not propagate the measurement error into an error in the fit parameters though.

While this does not answer the question asked, David Baird provides a simple explanation for propagating the error through a linear least squares fit in his book Experimentation: An Introduction to Measurement Theory and Experiment Design.

1. Essentially, weight the fit by dividing $x_i$ and $y_i$ values by the corresponding $dy_i$.
2. Perform the fit to get the parameters $m$ and $b$ in $$y=mx+b$$
3. Calculate the differences of each $y_i$ value from the fit, $\delta y_i$
4. Calculate the standard deviation of the fit parameters using $$\sigma_y = \sqrt{\frac{\sum(\delta y_i)^2}{n-2}}$$ $$\sigma_m = \sigma_y \sqrt{\frac{n}{n\sum{x_i^2}-\left(\sum{x_i}\right)^2}}$$ $$\sigma_b = \sigma_y \sqrt{\frac{\sum{x_i^2}}{n\sum{x_i^2}-\left(\sum{x_i}\right)^2}}$$

I am not certain this is a complete propagation of the error but it seems reasonable.

While this does not answer the question asked (and will therefore not be accepted), I provide this response in hopes it will help others and promote worthwhile discussion.

David Baird provides a simple explanation for propagating the error through a linear least squares fit in his book Experimentation: An Introduction to Measurement Theory and Experiment Design.

1. Essentially, weight the fit by dividing $x_i$ and $y_i$ values by the corresponding $dy_i$. As commented below, when reducing this to a matrix problem as described in the OP's wiki-link, make sure to also weight the $1$, i.e., $X=[x, 1]/dy$.
2. Perform the fit to get the parameters $m$ and $b$ in $$y=mx+b$$
3. Calculate the differences of each $y_i$ value from the fit, $\delta y_i$
4. Calculate the standard deviation of the fit parameters using $$\sigma_y = \sqrt{\frac{\sum(\delta y_i)^2}{n-2}}$$ $$\sigma_m = \sigma_y \sqrt{\frac{n}{n\sum{x_i^2}-\left(\sum{x_i}\right)^2}}$$ $$\sigma_b = \sigma_y \sqrt{\frac{\sum{x_i^2}}{n\sum{x_i^2}-\left(\sum{x_i}\right)^2}}$$

I am not certain this is a complete propagation of the error but it seems reasonable.