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](https://books.google.com/books?id=LicvAAAAIAAJ&dq=editions%3AD4pkF1UZ6JsC&q=least+squares+fit#v=onepage&q&f=false).  

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.