use qr, not pinv (#291)

* use qr, not pinv

* document weights usage in fit

* adjust range call for 1.0 testing

Author: jverzani

src/common.jl

@@ -69,7 +69,21 @@ fromroots(A::AbstractMatrix{T}; var::SymbolLike = :x) where {T <: Number} =
     fit(x, y, deg=length(x) - 1; [weights], var=:x)
     fit(::Type{<:AbstractPolynomial}, x, y, deg=length(x)-1; [weights], var=:x)
 
-Fit the given data as a polynomial type with the given degree. Uses linear least squares. When weights are given, as either a `Number`, `Vector` or `Matrix`, will use weighted linear least squares. The default polynomial type is [`Polynomial`](@ref). This will automatically scale your data to the [`domain`](@ref) of the polynomial type using [`mapdomain`](@ref)
+Fit the given data as a polynomial type with the given degree. Uses
+linear least squares to minimize the norm of `V⋅c - y`, where `V` is
+the Vandermonde matrix and `c` are the coefficients of the polynomial
+fit.
+
+This will automatically scale your data to the [`domain`](@ref) of the
+polynomial type using [`mapdomain`](@ref). The default polynomial type
+is [`Polynomial`](@ref).
+
+When weights are given, as either a `Number`, `Vector` or `Matrix`,
+this will use weighted linear least squares. That is, the norm of
+`W ⋅ (y - V ⋅ x)` is minimized. (As of now, the weights are specified 
+using their squares: for a number use `w^2`, for a vector `wᵢ^2`, and for a matrix
+ specify `W'*W`. This behavior may change in the future.)
+
 """
 function fit(P::Type{<:AbstractPolynomial},
              x::AbstractVector{T},
@@ -112,17 +126,21 @@ function _fit(P::Type{<:AbstractPolynomial},
     if weights !== nothing
         coeffs = _wlstsq(vand, y, weights)
     else
-        coeffs = pinv(vand) * y
+        coeffs = qr(vand) \ y
     end
     R = float(T)
     return P(R.(coeffs), var)
 end
 
 
 # Weighted linear least squares
+# TODO: Breaking change for 2.0: use non-squared weights
 _wlstsq(vand, y, W::Number) = _wlstsq(vand, y, fill!(similar(y), W))
-_wlstsq(vand, y, W::AbstractVector) = _wlstsq(vand, y, Diagonal(W))
-_wlstsq(vand, y, W::AbstractMatrix) = (vand' * W * vand) \ (vand' * W * y)
+function _wlstsq(vand, y, w::AbstractVector)
+    W = Diagonal(sqrt.(w))
+    qr(W * vand) \ (W * y)
+end
+_wlstsq(vand, y, W::AbstractMatrix) = qr(vand' * W * vand) \ (vand' * W * y)
 
 """
     roots(::AbstractPolynomial; kwargs...)

test/StandardBasis.jl

@@ -326,6 +326,16 @@ end
 
     # test default   (issue  #228)
     fit(1:3,  rand(3))
+
+    # weight test (PR #291)
+    # This should change with 2.0
+    # as for now we specify w^2.
+    x = range(0, stop=pi, length=30)
+    y = sin.(x)
+    wts = 1 ./ sqrt.(1 .+ x)
+    # cs = numpy.polynomial.polynomial.polyfit(x, y, 4, w=wts)
+    cs = [0.0006441172319036863, 0.985961582190304, 0.04999233434370933, -0.23162369757680354, 0.036864056406570644]
+    @test maximum(abs, coeffs(fit(x, y, 4, weights=wts.^2)) - cs) <= sqrt(eps())
 end
 
 @testset "Values" begin