From 25bb813f3225a8dbf56139665488b6903cb5c03a Mon Sep 17 00:00:00 2001 From: getzze Date: Fri, 2 Jun 2023 12:34:41 +0100 Subject: [PATCH 1/9] ensure same eltype for X and y --- src/robustlinearmodel.jl | 27 +++++++++++++++++++++------ 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/src/robustlinearmodel.jl b/src/robustlinearmodel.jl index 38c4e79..eb0234f 100644 --- a/src/robustlinearmodel.jl +++ b/src/robustlinearmodel.jl @@ -56,7 +56,18 @@ end ## TODO: specialize to make it faster StatsAPI.leverage(m::AbstractRobustModel) = diag(projectionmatrix(m)) -leverage_weights(m::AbstractRobustModel) = sqrt.(1 .- leverage(m)) +leverage_weights(m::AbstractRobustModel) = sqrt.(1 .- clamp.(leverage(m), 0, 1)) + +## Throw a MethodError for unspecified `args` to avoid StackOverflowError +function StatsAPI.fit( + ::Type{M}, + X::AbstractMatrix{<:AbstractFloat}, + y::AbstractVector{<:AbstractFloat}, + args...; + kwargs..., +) where {M<:AbstractRobustModel} + throw(MethodError(fit, (M, X, y, args...))) +end ## Convert to float, optionally drop rows with missing values (and convert to Non-Missing types) function StatsAPI.fit( @@ -86,7 +97,9 @@ function StatsAPI.fit( X, y, _ = missing_omit(X, y) end - return fit(M, float(X), float(y), args...; kwargs...) + # Make sure X and y have the same float eltype + T = promote_type(float(eltype(X)), float(eltype(y))) + return fit(M, convert.(T, X), convert.(T, y), args...; kwargs...) end ## Convert from formula-data to modelmatrix-response calling form @@ -156,6 +169,7 @@ function StatsModels.formula(m::RobustLinearModel) return m.formula end + """ deviance(m::RobustLinearModel) @@ -493,7 +507,7 @@ the RobustLinearModel object. """ function StatsAPI.fit( ::Type{M}, - X::Union{AbstractMatrix{T},SparseMatrixCSC{T}}, + X::AbstractMatrix{T}, y::AbstractVector{T}, est::AbstractMEstimator; method::Symbol=:chol, # :qr, :cg @@ -1350,19 +1364,20 @@ function resampling_best_estimate( ## TODO: implement something similar to DetS (not sure it could apply) ## Hubert2015 - The DetS and DetMM estimators for multivariate location and scatter ## (https://www.sciencedirect.com/science/article/abs/pii/S0167947314002175) + M = length(coef(m)) if isnothing(Nsamples) - Nsamples = resampling_minN(dof(m), 0.05, propoutliers) + Nsamples = resampling_minN(M, 0.05, propoutliers) end if isnothing(Npoints) - Npoints = dof(m) + Npoints = M end Nsubsamples = min(Nsubsamples, Nsamples) verbose && println("Start $(Nsamples) subsamples...") σis = zeros(Nsamples) - βis = zeros(dof(m), Nsamples) + βis = zeros(M, Nsamples) for i in 1:Nsamples # TODO: to parallelize, make a deepcopy of m inds = sample(rng, axes(response(m), 1), Npoints; replace=false, ordered=false) From 843a0fb8b28d21cb89880f73496b13499b9f94e0 Mon Sep 17 00:00:00 2001 From: getzze Date: Fri, 2 Jun 2023 12:52:36 +0100 Subject: [PATCH 2/9] correctly use QR --- Project.toml | 1 + src/RobustModels.jl | 16 ++++++--- src/linpred.jl | 82 +++++++++++++++++++++++++++++++++++---------- 3 files changed, 77 insertions(+), 22 deletions(-) diff --git a/Project.toml b/Project.toml index eade767..673f555 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a" IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Missings = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/src/RobustModels.jl b/src/RobustModels.jl index 7d57ca6..e22d4cb 100644 --- a/src/RobustModels.jl +++ b/src/RobustModels.jl @@ -1,5 +1,7 @@ module RobustModels +using Pkg: Pkg + include("compat.jl") # Use README as the docstring of the module and doctest README @@ -10,10 +12,10 @@ end RobustModels # Import with `using` to use the module names to prefix the methods # that are extended from these modules -using GLM -using StatsAPI -using StatsBase -using StatsModels +using GLM: GLM +using StatsAPI: StatsAPI +using StatsBase: StatsBase +using StatsModels: StatsModels ## Import to implement new methods import Base: show, broadcastable, convert, == @@ -73,6 +75,9 @@ using LinearAlgebra: inv, diag, diagm, + Diagonal, + rank, + qr, ldiv! using Random: AbstractRNG, GLOBAL_RNG @@ -83,8 +88,10 @@ using StatsBase: using StatsModels: @delegate, @formula, + formula, RegressionModel, FormulaTerm, + InterceptTerm, ModelFrame, modelcols, apply_schema, @@ -238,6 +245,7 @@ abstract type AbstractRegularizedPred{T} end Base.broadcastable(m::T) where {T<:AbstractEstimator} = Ref(m) Base.broadcastable(m::T) where {T<:LossFunction} = Ref(m) + include("tools.jl") include("losses.jl") include("estimators.jl") diff --git a/src/linpred.jl b/src/linpred.jl index 5e61e78..bdedca2 100644 --- a/src/linpred.jl +++ b/src/linpred.jl @@ -65,27 +65,73 @@ A `LinPred` type with a dense, unpivoted QR decomposition of `X` """ DensePredQR -PRED_QR_WARNING_ISSUED = false - -function qrpred(X::AbstractMatrix, pivot::Bool=false) - try - return DensePredCG(Matrix(X), pivot) - catch e - if e isa MethodError - # GLM.DensePredCG(X::AbstractMatrix, pivot::Bool) is not defined - global PRED_QR_WARNING_ISSUED - if !PRED_QR_WARNING_ISSUED - @warn( - "GLM.DensePredCG(X::AbstractMatrix, pivot::Bool) is not defined, " * - "fallback to unpivoted QR. GLM version should be >= 1.9." - ) - PRED_QR_WARNING_ISSUED = true - end - return DensePredCG(Matrix(X)) +function get_pkg_version(m::Module) + toml = Pkg.TOML.parsefile(joinpath(pkgdir(m), "Project.toml")) + return VersionNumber(toml["version"]) +end + +@static if get_pkg_version(GLM) < v"1.9" + @warn( + "GLM.DensePredQR(X::AbstractMatrix, pivot::Bool) is not defined, " * + "fallback to unpivoted QR. GLM version should be >= 1.9." + ) + + # GLM.DensePredQR(X::AbstractMatrix, pivot::Bool) is not defined + function qrpred(X::AbstractMatrix, pivot::Bool=false) + DensePredQR(Matrix(X)) + end + + # GLM.delbeta!(p::DensePredQR{T}, r::Vector{T}, wt::Vector{T}) is not defined + function delbeta!(p::DensePredQR{T}, r::Vector{T}, wt::Vector{T}) where T<:BlasReal + rnk = rank(p.qr.R) + X = p.X + W = Diagonal(wt) + sqrtW = Diagonal(sqrt.(wt)) + scratchm1 = similar(X, T) + mul!(scratchm1, sqrtW, X) + + n, m = size(X) + if n >= m + # W½ X = Q R , with Q'Q = I + # X'WX β = X'y => R'Q'QR β = X'y + # => β = R⁻¹ R⁻ᵀ X'y + qnr = qr(scratchm1) + Rinv = inv(qnr.R) + + scratchm2 = similar(X, T) + mul!(scratchm2, W, X) + mul!(p.delbeta, transpose(scratchm2), r) + + p.delbeta = Rinv * Rinv' * p.delbeta else - rethrow() + # (W½ X)' = Q R , with Q'Q = I + # W½X β = W½y => R'Q' β = y + # => β = Q . [R⁻ᵀ y; 0] + qnr = qr(scratchm1') + RTinv = inv(qnr.R)' + @assert 1 <= n <= size(p.delbeta, 1) + mul!(view(p.delbeta, 1:n), RTinv, r) + p.delbeta = zeros(size(p.delbeta)) + p.delbeta[1:n] .= RTinv * r + lmul!(qnr.Q, p.delbeta) end + return p end + + # GLM.delbeta!(p::DensePredQR{T}, r::Vector{T}) is ill-defined + function delbeta!(p::DensePredQR{T}, r::Vector{T}) where T<:BlasReal + n, m = size(p.X) + if n >= m + p.delbeta = p.qr \ r + else + qnrT = qr(p.X') + p.delbeta = qnrT' \ r + end + return p + end + +else + qrpred(X::AbstractMatrix, pivot::Bool=false) = DensePredQR(Matrix(X), pivot) end From 32f122fb282cd1998a8e1576876d6bcba10e4725 Mon Sep 17 00:00:00 2001 From: getzze Date: Fri, 2 Jun 2023 13:19:55 +0100 Subject: [PATCH 3/9] use :cholesky and dropcollinear like in GLM --- src/robustlinearmodel.jl | 25 +++++----- test/underdetermined.jl | 104 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 117 insertions(+), 12 deletions(-) create mode 100644 test/underdetermined.jl diff --git a/src/robustlinearmodel.jl b/src/robustlinearmodel.jl index eb0234f..86ee199 100644 --- a/src/robustlinearmodel.jl +++ b/src/robustlinearmodel.jl @@ -442,7 +442,7 @@ rlm(X, y, args...; kwargs...) = fit(RobustLinearModel, X, y, args...; kwargs...) ridgeG::Union{UniformScaling, AbstractArray} = I, βprior::AbstractVector = [], quantile::Union{Nothing, AbstractFloat} = nothing, - pivot::Bool = false, + dropcollinear::Bool = false, initial_scale::Union{Symbol, Real}=:mad, σ0::Union{Nothing, Symbol, Real}=initial_scale, initial_coef::AbstractVector=[], @@ -481,7 +481,8 @@ using a robust estimator. Default to `zeros(p)`; - `quantile::Union{Nothing, AbstractFloat} = nothing`: only for [`GeneralizedQuantileEstimator`](@ref), define the quantile to estimate; -- `pivot::Bool=false`: use pivoted factorization; +- `dropcollinear::Bool=false`: controls whether or not a model matrix less-than-full rank is + accepted; - `contrasts::AbstractDict{Symbol,Any} = Dict{Symbol,Any}()`: a `Dict` mapping term names (as `Symbol`s) to term types (e.g. `ContinuousTerm`) or contrasts (e.g., `HelmertCoding()`, `SeqDiffCoding(; levels=["a", "b", "c"])`, etc.). If contrasts are not provided for a variable, @@ -522,7 +523,7 @@ function StatsAPI.fit( ridgeG::Union{UniformScaling,AbstractArray{<:Real}}=I, βprior::AbstractVector{<:Real}=T[], quantile::Union{Nothing,AbstractFloat}=nothing, - pivot::Bool=false, + dropcollinear::Bool=false, contrasts::AbstractDict{Symbol,Any}=Dict{Symbol,Any}(), # placeholder __formula::Union{Nothing,FormulaTerm}=nothing, fitargs..., @@ -553,7 +554,7 @@ function StatsAPI.fit( rr = RobustLinResp(est, y, offset, wts) # Predictor object - methods = (:auto, :chol, :cg, :qr) + methods = (:auto, :chol, :cholesky, :cg, :qr) if method ∉ methods @warn("Incorrect method `:$(method)`, should be one of $(methods)") method = :auto @@ -573,22 +574,22 @@ function StatsAPI.fit( ridgeG end if method == :cg - cgpred(X, float(ridgeλ), G, βprior, pivot) + cgpred(X, float(ridgeλ), G, βprior, dropcollinear) elseif method == :qr - qrpred(X, float(ridgeλ), G, βprior, pivot) - elseif method in (:chol, :auto) - cholpred(X, float(ridgeλ), G, βprior, pivot) + qrpred(X, float(ridgeλ), G, βprior, dropcollinear) + elseif method in (:chol, :cholesky, :auto) + cholpred(X, float(ridgeλ), G, βprior, dropcollinear) else error("method :$method is not allowed, should be in: $methods") end else # No regularization if method == :cg - cgpred(X, pivot) + cgpred(X, dropcollinear) elseif method == :qr - qrpred(X, pivot) - elseif method in (:chol, :auto) - cholpred(X, pivot) + qrpred(X, dropcollinear) + elseif method in (:chol, :cholesky, :auto) + cholpred(X, dropcollinear) else error("method :$method is not allowed, should be in: $methods") end diff --git a/test/underdetermined.jl b/test/underdetermined.jl new file mode 100644 index 0000000..4367406 --- /dev/null +++ b/test/underdetermined.jl @@ -0,0 +1,104 @@ + + +seed = 1234 +rng = MersenneTwister(seed) + +function gendata( + rng::AbstractRNG, + ::Type{T}, + N::Integer, + p::Integer, + nonzero::Integer, + rho::Real, + snr::Real = 3, + alternate::Bool=true, +) where {T<:AbstractFloat} + if nonzero > 0 + beta = nonzero:-1:1 + if alternate + beta = [(nonzero - i + 1) * (i % 2 == 0 ? 1 : -1) for i in 1:nonzero] + else + beta = nonzero:-1:1 + end + beta = vcat(beta, zeros(T, p-nonzero)) + else + beta = zeros(T, p) + end + + X = randn(rng, N, p) + if rho > 0 + z = randn(N) + X .+= sqrt(rho/(1-rho)) * z + X *= sqrt(1-rho) + end + + ssd = sqrt((1-rho)*sum(abs2, beta) + rho*sum(abs, beta)) + nsd = ssd / snr + + y = X * beta + nsd * randn(N) + return X, y, beta +end + + +X3, y3, β03 = gendata(rng, Float64, 100, 1000, 30, 0.3) + +loss1 = L2Loss() +est1 = MEstimator{L2Loss}() +λ = 1.0 +pen1 = L1Penalty(λ) +pen2 = SquaredL2Penalty(λ) + +@testset "Rank deficient X: penalty ($(pen))" for (pen, methods) in zip( + (nothing, pen1, pen2), (nopen_methods, pen_methods, pen_methods) +) + def_rtol = isnothing(pen) ? 1e-6 : 1e-4 + + @testset "solver method $(method)" for method in methods + rtol = def_rtol + if method === :fista + rtol = 1e-1 + end + + ### RobustLinearModel + name = "rlm(X3, y3, L2Loss, $(pen); method=$(method))" + VERBOSE && println("\n\t\u25CF $(name)") + + if isnothing(pen) + if method in (:chol, :cholesky, :auto) + @test_throws Exception rlm(X3, y3, est1; method=method, initial_scale=1, dropcollinear=false) + end + m1 = rlm(X3, y3, est1; method=method, initial_scale=1, dropcollinear=true) + else + m1 = rlm(X3, y3, pen; method=method, initial_scale=1, maxiter=1000, dropcollinear=true) + end + + # Test printing the model + @test_nowarn show(devnull, m1) + + # interface + @testset "method: $(f)" for f in interface_methods + @test_nowarn f(m1) + end + + ### IPODRegression + name = "ipod(X3, y3, L2Loss, $(pen); method=$(method))" + VERBOSE && println("\n\t\u25CF $(name)") + + if isnothing(pen) && method in (:chol, :cholesky, :auto) + @test_throws Exception ipod(X3, y3, loss1, pen; method=method, initial_scale=1, dropcollinear=false) + end + m2 = ipod(X3, y3, loss1, pen; method=method, initial_scale=1, maxiter=1000, dropcollinear=true) + + # Test printing the model + @test_nowarn show(devnull, m2) + + # interface + @testset "method: $(f)" for f in interface_methods + if f == workingweights + # Not defined for IPODRegression + continue + end + @test_nowarn f(m2) + end + end +end From d6fd9a594e0089017e5e753079a8c6be52aa5faa Mon Sep 17 00:00:00 2001 From: getzze Date: Fri, 2 Jun 2023 12:53:27 +0100 Subject: [PATCH 4/9] test method :cholesky --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 5fed115..5794a3b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -49,7 +49,7 @@ bounded_losses = ("Geman", "Welsch", "Tukey", "YohaiZamar", "HardThreshold", "Ha losses = (convex_losses..., other_losses..., bounded_losses...) ## Solving methods -nopen_methods = (:auto, :chol, :qr, :cg) +nopen_methods = (:auto, :chol, :cholesky, :qr, :cg) ## Interface methods interface_methods = ( From 7b7bb88d7352ef747f1384e87fd399a957b827af Mon Sep 17 00:00:00 2001 From: getzze Date: Fri, 2 Jun 2023 13:22:38 +0100 Subject: [PATCH 5/9] remove unused test --- test/runtests.jl | 1 + test/underdetermined.jl | 104 ---------------------------------------- 2 files changed, 1 insertion(+), 104 deletions(-) delete mode 100644 test/underdetermined.jl diff --git a/test/runtests.jl b/test/runtests.jl index 5794a3b..60e92f4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -94,3 +94,4 @@ include("robustridge.jl") include("qreg.jl") include("weights.jl") include("univariate.jl") +include("underdetermined.jl") diff --git a/test/underdetermined.jl b/test/underdetermined.jl deleted file mode 100644 index 4367406..0000000 --- a/test/underdetermined.jl +++ /dev/null @@ -1,104 +0,0 @@ - - -seed = 1234 -rng = MersenneTwister(seed) - -function gendata( - rng::AbstractRNG, - ::Type{T}, - N::Integer, - p::Integer, - nonzero::Integer, - rho::Real, - snr::Real = 3, - alternate::Bool=true, -) where {T<:AbstractFloat} - if nonzero > 0 - beta = nonzero:-1:1 - if alternate - beta = [(nonzero - i + 1) * (i % 2 == 0 ? 1 : -1) for i in 1:nonzero] - else - beta = nonzero:-1:1 - end - beta = vcat(beta, zeros(T, p-nonzero)) - else - beta = zeros(T, p) - end - - X = randn(rng, N, p) - if rho > 0 - z = randn(N) - X .+= sqrt(rho/(1-rho)) * z - X *= sqrt(1-rho) - end - - ssd = sqrt((1-rho)*sum(abs2, beta) + rho*sum(abs, beta)) - nsd = ssd / snr - - y = X * beta + nsd * randn(N) - return X, y, beta -end - - -X3, y3, β03 = gendata(rng, Float64, 100, 1000, 30, 0.3) - -loss1 = L2Loss() -est1 = MEstimator{L2Loss}() -λ = 1.0 -pen1 = L1Penalty(λ) -pen2 = SquaredL2Penalty(λ) - -@testset "Rank deficient X: penalty ($(pen))" for (pen, methods) in zip( - (nothing, pen1, pen2), (nopen_methods, pen_methods, pen_methods) -) - def_rtol = isnothing(pen) ? 1e-6 : 1e-4 - - @testset "solver method $(method)" for method in methods - rtol = def_rtol - if method === :fista - rtol = 1e-1 - end - - ### RobustLinearModel - name = "rlm(X3, y3, L2Loss, $(pen); method=$(method))" - VERBOSE && println("\n\t\u25CF $(name)") - - if isnothing(pen) - if method in (:chol, :cholesky, :auto) - @test_throws Exception rlm(X3, y3, est1; method=method, initial_scale=1, dropcollinear=false) - end - m1 = rlm(X3, y3, est1; method=method, initial_scale=1, dropcollinear=true) - else - m1 = rlm(X3, y3, pen; method=method, initial_scale=1, maxiter=1000, dropcollinear=true) - end - - # Test printing the model - @test_nowarn show(devnull, m1) - - # interface - @testset "method: $(f)" for f in interface_methods - @test_nowarn f(m1) - end - - ### IPODRegression - name = "ipod(X3, y3, L2Loss, $(pen); method=$(method))" - VERBOSE && println("\n\t\u25CF $(name)") - - if isnothing(pen) && method in (:chol, :cholesky, :auto) - @test_throws Exception ipod(X3, y3, loss1, pen; method=method, initial_scale=1, dropcollinear=false) - end - m2 = ipod(X3, y3, loss1, pen; method=method, initial_scale=1, maxiter=1000, dropcollinear=true) - - # Test printing the model - @test_nowarn show(devnull, m2) - - # interface - @testset "method: $(f)" for f in interface_methods - if f == workingweights - # Not defined for IPODRegression - continue - end - @test_nowarn f(m2) - end - end -end From 431ba6d62700b2c18c81eb39383742f5f70315bd Mon Sep 17 00:00:00 2001 From: getzze Date: Fri, 2 Jun 2023 16:59:21 +0100 Subject: [PATCH 6/9] create own DensePredQR --- src/RobustModels.jl | 4 +- src/compat.jl | 16 ++++- src/linpred.jl | 133 ++++++++++++++++++++++++++++------------- src/regularizedpred.jl | 3 + test/robustridge.jl | 12 ++-- test/runtests.jl | 1 - 6 files changed, 117 insertions(+), 52 deletions(-) diff --git a/src/RobustModels.jl b/src/RobustModels.jl index e22d4cb..9933bb4 100644 --- a/src/RobustModels.jl +++ b/src/RobustModels.jl @@ -75,14 +75,12 @@ using LinearAlgebra: inv, diag, diagm, - Diagonal, rank, - qr, ldiv! using Random: AbstractRNG, GLOBAL_RNG using Printf: @printf, @sprintf -using GLM: FPVector, lm, SparsePredChol, DensePredChol, DensePredQR +using GLM: FPVector, lm, SparsePredChol, DensePredChol using StatsBase: AbstractWeights, CoefTable, ConvergenceException, median, mad, mad_constant, sample using StatsModels: diff --git a/src/compat.jl b/src/compat.jl index ab063f1..26619fa 100644 --- a/src/compat.jl +++ b/src/compat.jl @@ -1,4 +1,10 @@ -using LinearAlgebra: cholesky! +using LinearAlgebra: cholesky!, qr! + +function get_pkg_version(m::Module) + toml = Pkg.TOML.parsefile(joinpath(pkgdir(m), "Project.toml")) + return VersionNumber(toml["version"]) +end + ## Compatibility layers @@ -6,5 +12,13 @@ using LinearAlgebra: cholesky! @static if VERSION < v"1.8.0-DEV.1139" pivoted_cholesky!(A; kwargs...) = cholesky!(A, Val(true); kwargs...) else + using LinearAlgebra: RowMaximum pivoted_cholesky!(A; kwargs...) = cholesky!(A, RowMaximum(); kwargs...) end + +@static if VERSION < v"1.7.0" + pivoted_qr!(A; kwargs...) = qr!(A, Val(true); kwargs...) +else + using LinearAlgebra: ColumnNorm + pivoted_qr!(A; kwargs...) = qr!(A, ColumnNorm(); kwargs...) +end diff --git a/src/linpred.jl b/src/linpred.jl index bdedca2..51cde8c 100644 --- a/src/linpred.jl +++ b/src/linpred.jl @@ -50,44 +50,100 @@ leverage_weights(p::LinPred, wt::AbstractVector) = sqrt.(1 .- leverage(p, wt)) # beta0 #end -""" - DensePredQR - -A `LinPred` type with a dense, unpivoted QR decomposition of `X` - -# Members - -- `X`: Model matrix of size `n` × `p` with `n ≥ p`. Should be full column rank. -- `beta0`: base coefficient vector of length `p` -- `delbeta`: increment to coefficient vector, also of length `p` -- `scratchbeta`: scratch vector of length `p`, used in `linpred!` method -- `qr`: a `QRCompactWY` object created from `X`, with optional row weights. -""" -DensePredQR -function get_pkg_version(m::Module) - toml = Pkg.TOML.parsefile(joinpath(pkgdir(m), "Project.toml")) - return VersionNumber(toml["version"]) -end +########################################## +###### DensePredQR +########################################## -@static if get_pkg_version(GLM) < v"1.9" +#@static if get_pkg_version(GLM) < v"1.9" +@static if true @warn( - "GLM.DensePredQR(X::AbstractMatrix, pivot::Bool) is not defined, " * - "fallback to unpivoted QR. GLM version should be >= 1.9." + "GLM.DensePredQR(X::AbstractMatrix, pivot::Bool=true) is not defined, " * + "fallback to unpivoted RobustModels.DensePredQR definition. " * + "To use pivoted QR, GLM version should be greater than or equal to v1.9." ) + using LinearAlgebra: QRCompactWY, QRPivoted, Diagonal, qr!, qr + + """ + DensePredQR + + A `LinPred` type with a dense QR decomposition of `X` + + # Members + + - `X`: Model matrix of size `n` × `p` with `n ≥ p`. Should be full column rank. + - `beta0`: base coefficient vector of length `p` + - `delbeta`: increment to coefficient vector, also of length `p` + - `scratchbeta`: scratch vector of length `p`, used in `linpred!` method + - `qr`: a `QRCompactWY` object created from `X`, with optional row weights. + - `scratchm1`: scratch Matrix{T} of the same size as `X` + - `scratchm2`: scratch Matrix{T} of the same size as `X` + - `scratchR`: scratch Matrix{T} of the same size as `qr.R`, a square matrix. + """ + mutable struct DensePredQR{T<:BlasReal,Q<:Union{QRCompactWY,QRPivoted}} <: DensePred + X::Matrix{T} # model matrix + beta0::Vector{T} # base coefficient vector + delbeta::Vector{T} # coefficient increment + scratchbeta::Vector{T} + qr::Q + scratchm1::Matrix{T} + scratchm2::Matrix{T} + scratchR::Matrix{T} + + function DensePredQR(X::AbstractMatrix, pivot::Bool=false) + n, p = size(X) + T = typeof(float(zero(eltype(X)))) + + if false +# if pivot + F = pivoted_qr!(copy(X)) + else + if n >= p + F = qr(X) + else + # adjoint of X so R is square + # cannot use in-place qr! + F = qr(X) + end + end + + new{T,typeof(F)}( + Matrix{T}(X), + zeros(T, p), + zeros(T, p), + zeros(T, p), + F, + similar(X, T), + similar(X, T), + zeros(T, size(F.R)), + ) + end + end + # GLM.DensePredQR(X::AbstractMatrix, pivot::Bool) is not defined function qrpred(X::AbstractMatrix, pivot::Bool=false) DensePredQR(Matrix(X)) end + # GLM.delbeta!(p::DensePredQR{T}, r::Vector{T}) is ill-defined + function delbeta!(p::DensePredQR{T,<:QRCompactWY}, r::Vector{T}) where T<:BlasReal + n, m = size(p.X) + if n >= m + p.delbeta = p.qr \ r + else + p.delbeta = p.qr' \ r + end + return p + end + # GLM.delbeta!(p::DensePredQR{T}, r::Vector{T}, wt::Vector{T}) is not defined - function delbeta!(p::DensePredQR{T}, r::Vector{T}, wt::Vector{T}) where T<:BlasReal + function delbeta!(p::DensePredQR{T,<:QRCompactWY}, r::Vector{T}, wt::Vector{T}) where T<:BlasReal rnk = rank(p.qr.R) X = p.X W = Diagonal(wt) sqrtW = Diagonal(sqrt.(wt)) - scratchm1 = similar(X, T) + scratchm1 = p.scratchm1 = similar(X, T) mul!(scratchm1, sqrtW, X) n, m = size(X) @@ -95,10 +151,10 @@ end # W½ X = Q R , with Q'Q = I # X'WX β = X'y => R'Q'QR β = X'y # => β = R⁻¹ R⁻ᵀ X'y - qnr = qr(scratchm1) - Rinv = inv(qnr.R) + qnr = p.qr = qr(scratchm1) + Rinv = p.scratchR = inv(qnr.R) - scratchm2 = similar(X, T) + scratchm2 = p.scratchm2 = similar(X, T) mul!(scratchm2, W, X) mul!(p.delbeta, transpose(scratchm2), r) @@ -107,34 +163,29 @@ end # (W½ X)' = Q R , with Q'Q = I # W½X β = W½y => R'Q' β = y # => β = Q . [R⁻ᵀ y; 0] - qnr = qr(scratchm1') - RTinv = inv(qnr.R)' + qnrT = p.qr = qr(scratchm1') + RTinv = p.scratchR = inv(qnrT.R)' @assert 1 <= n <= size(p.delbeta, 1) mul!(view(p.delbeta, 1:n), RTinv, r) p.delbeta = zeros(size(p.delbeta)) p.delbeta[1:n] .= RTinv * r - lmul!(qnr.Q, p.delbeta) + lmul!(qnrT.Q, p.delbeta) end return p end - # GLM.delbeta!(p::DensePredQR{T}, r::Vector{T}) is ill-defined - function delbeta!(p::DensePredQR{T}, r::Vector{T}) where T<:BlasReal - n, m = size(p.X) - if n >= m - p.delbeta = p.qr \ r - else - qnrT = qr(p.X') - p.delbeta = qnrT' \ r - end - return p - end +## Use DensePredQR from GLM else - qrpred(X::AbstractMatrix, pivot::Bool=false) = DensePredQR(Matrix(X), pivot) + using GLM: DensePredQR + import GLM: qrpred end +########################################## +###### [Dense/Sparse]PredCG +########################################## + """ DensePredCG diff --git a/src/regularizedpred.jl b/src/regularizedpred.jl index 9896b39..513aac8 100644 --- a/src/regularizedpred.jl +++ b/src/regularizedpred.jl @@ -154,6 +154,9 @@ function postupdate_λ!(r::RidgePred) # Update the extended model matrix with the new value GG = r.sqrtλ * r.G @views r.pred.X[(n + 1):(n + m), :] .= GG + + # Update other fields + # TODO: update DensePredQR if isa(r.pred, DensePredChol) # Recompute the cholesky decomposition X = r.pred.X diff --git a/test/robustridge.jl b/test/robustridge.jl index 87b22a4..85a1350 100644 --- a/test/robustridge.jl +++ b/test/robustridge.jl @@ -33,12 +33,12 @@ est2 = MEstimator(loss2) kwargs = (; method=method, initial_scale=:L1) # use the dispersion from GLM to ensure that the loglikelihood is correct m2 = fit(RobustLinearModel, A, b, est; kwargs...) - m3 = fit(RobustLinearModel, A, b, est; ridgeλ=1, kwargs...) + m3 = fit(RobustLinearModel, A, b, est; ridgeλ=10, kwargs...) m4 = fit( - RobustLinearModel, A, b, est; ridgeλ=1, ridgeG=float([0 0; 0 1]), kwargs... + RobustLinearModel, A, b, est; ridgeλ=10, ridgeG=float([0 0; 0 1]), kwargs... ) m5 = fit(RobustLinearModel, A, b, est; ridgeλ=0.1, ridgeG=[0, 1], kwargs...) - m6 = rlm(A, b, est; ridgeλ=1, ridgeG=[0, 1], βprior=[0.0, 2.0], kwargs...) + m6 = rlm(A, b, est; ridgeλ=10, ridgeG=[0, 1], βprior=[0.0, 2.0], kwargs...) VERBOSE && println("\n\t\u25CF Estimator: $(name)") VERBOSE && println(" lm$(aspace) : ", coef(m1)) @@ -53,9 +53,9 @@ est2 = MEstimator(loss2) @test_nowarn println(m3) # refit - refit!(m5; ridgeλ=1, initial_scale=:L1) - @test isapprox(m5.pred.λ, 1.0; rtol=1e-5) - @test isapprox(coef(m3), coef(m5); rtol=1e-5) + refit!(m5; ridgeλ=10, initial_scale=:L1) + @test isapprox(m5.pred.λ, 10.0; rtol=1e-5) + @test isapprox(coef(m3), coef(m5); rtol=1e-6) end end diff --git a/test/runtests.jl b/test/runtests.jl index 60e92f4..5794a3b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -94,4 +94,3 @@ include("robustridge.jl") include("qreg.jl") include("weights.jl") include("univariate.jl") -include("underdetermined.jl") From ef37ed740cde04006ae24e13ea2dc2ac8669deec Mon Sep 17 00:00:00 2001 From: getzze Date: Fri, 2 Jun 2023 17:01:54 +0100 Subject: [PATCH 7/9] format --- src/linpred.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/linpred.jl b/src/linpred.jl index 51cde8c..26b0d17 100644 --- a/src/linpred.jl +++ b/src/linpred.jl @@ -55,8 +55,7 @@ leverage_weights(p::LinPred, wt::AbstractVector) = sqrt.(1 .- leverage(p, wt)) ###### DensePredQR ########################################## -#@static if get_pkg_version(GLM) < v"1.9" -@static if true +@static if get_pkg_version(GLM) < v"1.9" @warn( "GLM.DensePredQR(X::AbstractMatrix, pivot::Bool=true) is not defined, " * "fallback to unpivoted RobustModels.DensePredQR definition. " * @@ -96,7 +95,7 @@ leverage_weights(p::LinPred, wt::AbstractVector) = sqrt.(1 .- leverage(p, wt)) T = typeof(float(zero(eltype(X)))) if false -# if pivot + # if pivot F = pivoted_qr!(copy(X)) else if n >= p @@ -108,7 +107,7 @@ leverage_weights(p::LinPred, wt::AbstractVector) = sqrt.(1 .- leverage(p, wt)) end end - new{T,typeof(F)}( + return new{T,typeof(F)}( Matrix{T}(X), zeros(T, p), zeros(T, p), @@ -123,11 +122,11 @@ leverage_weights(p::LinPred, wt::AbstractVector) = sqrt.(1 .- leverage(p, wt)) # GLM.DensePredQR(X::AbstractMatrix, pivot::Bool) is not defined function qrpred(X::AbstractMatrix, pivot::Bool=false) - DensePredQR(Matrix(X)) + return DensePredQR(Matrix(X)) end # GLM.delbeta!(p::DensePredQR{T}, r::Vector{T}) is ill-defined - function delbeta!(p::DensePredQR{T,<:QRCompactWY}, r::Vector{T}) where T<:BlasReal + function delbeta!(p::DensePredQR{T,<:QRCompactWY}, r::Vector{T}) where {T<:BlasReal} n, m = size(p.X) if n >= m p.delbeta = p.qr \ r @@ -138,7 +137,9 @@ leverage_weights(p::LinPred, wt::AbstractVector) = sqrt.(1 .- leverage(p, wt)) end # GLM.delbeta!(p::DensePredQR{T}, r::Vector{T}, wt::Vector{T}) is not defined - function delbeta!(p::DensePredQR{T,<:QRCompactWY}, r::Vector{T}, wt::Vector{T}) where T<:BlasReal + function delbeta!( + p::DensePredQR{T,<:QRCompactWY}, r::Vector{T}, wt::Vector{T} + ) where {T<:BlasReal} rnk = rank(p.qr.R) X = p.X W = Diagonal(wt) @@ -175,7 +176,7 @@ leverage_weights(p::LinPred, wt::AbstractVector) = sqrt.(1 .- leverage(p, wt)) end -## Use DensePredQR from GLM + ## Use DensePredQR from GLM else using GLM: DensePredQR import GLM: qrpred From 862de71288d5bdc96736591c6f955b7999c1d979 Mon Sep 17 00:00:00 2001 From: getzze Date: Fri, 2 Jun 2023 17:42:26 +0100 Subject: [PATCH 8/9] cleanupt --- src/linpred.jl | 20 +++----------------- test/interface.jl | 4 ++-- test/linearfit.jl | 4 ++-- test/mquantile.jl | 4 ++-- test/robustridge.jl | 16 ++++++++++++---- 5 files changed, 21 insertions(+), 27 deletions(-) diff --git a/src/linpred.jl b/src/linpred.jl index 26b0d17..79c77bf 100644 --- a/src/linpred.jl +++ b/src/linpred.jl @@ -207,20 +207,8 @@ mutable struct DensePredCG{T<:BlasReal} <: DensePred scratchbeta::Vector{T} scratchm1::Matrix{T} scratchr1::Vector{T} - function DensePredCG{T}(X::Matrix{T}, beta0::Vector{T}) where {T} - n, p = size(X) - length(beta0) == p || throw(DimensionMismatch("length(β0) ≠ size(X,2)")) - return new{T}( - X, - beta0, - zeros(T, p), - zeros(T, (p, p)), - zeros(T, p), - zeros(T, (n, p)), - zeros(T, n), - ) - end - function DensePredCG{T}(X::Matrix{T}) where {T} + + function DensePredCG(X::Matrix{T}) where {T<:BlasReal} n, p = size(X) return new{T}( X, @@ -233,10 +221,8 @@ mutable struct DensePredCG{T<:BlasReal} <: DensePred ) end end -DensePredCG(X::Matrix, beta0::Vector) = DensePredCG{eltype(X)}(X, beta0) -DensePredCG(X::Matrix{T}) where {T} = DensePredCG{T}(X, zeros(T, size(X, 2))) function Base.convert(::Type{DensePredCG{T}}, X::Matrix{T}) where {T} - return DensePredCG{T}(X, zeros(T, size(X, 2))) + return DensePredCG(X) end # Compatibility with cholpred(X, pivot) diff --git a/test/interface.jl b/test/interface.jl index fe52707..ac7be31 100644 --- a/test/interface.jl +++ b/test/interface.jl @@ -37,7 +37,7 @@ est2 = MEstimator(loss2) ) β = copy(coef(m)) VERBOSE && println("rlm($name): ", β) - @test_nowarn println(m) + @test_nowarn show(devnull, m) @test isapprox(coef(m1), β; rtol=1e-5) # make sure that it is not a TableRegressionModel @@ -45,7 +45,7 @@ est2 = MEstimator(loss2) # refit refit!(m, y; verbose=false, initial_scale=:extrema) - @test_nowarn println("$m") + @test_nowarn show(devnull, "$m") @test all(coef(m) .== β) # refit diff --git a/test/linearfit.jl b/test/linearfit.jl index 7659d6a..241abe7 100644 --- a/test/linearfit.jl +++ b/test/linearfit.jl @@ -22,9 +22,9 @@ est2 = MEstimator(loss2) @test all(isfinite.(coef(m2))) if name != "L1" - @test_nowarn println(m2) + @test_nowarn show(devnull, m2) # else - # @test_warn L1_warning println(m2) + # @test_warn L1_warning show(devnull, m2) end VERBOSE && println("\n\t\u25CF Estimator: $(name)") diff --git a/test/mquantile.jl b/test/mquantile.jl index 129b514..fe2d71f 100644 --- a/test/mquantile.jl +++ b/test/mquantile.jl @@ -70,7 +70,7 @@ loss2 = RobustModels.TukeyLoss() VERBOSE && println("rlm(qr): ", coef(m4)) if name == "Expectile" - @test_nowarn println(m2) + @test_nowarn show(devnull, m2) @test isapprox(coef(m2), coef(m3); rtol=1e-4) @test isapprox(coef(m2), coef(m4); rtol=1e-4) @@ -82,7 +82,7 @@ loss2 = RobustModels.TukeyLoss() refit!(m1; quantile=τ) @test isapprox(coef(m1), coef(m2); rtol=1e-4) else - # @test_warn L1_warning println(m2) + # @test_warn L1_warning show(devnull, m2) @test isapprox(coef(m2), coef(m3); rtol=1e-2) @test isapprox(coef(m2), coef(m4); rtol=1e-2) diff --git a/test/robustridge.jl b/test/robustridge.jl index 85a1350..5fd1673 100644 --- a/test/robustridge.jl +++ b/test/robustridge.jl @@ -38,15 +38,19 @@ est2 = MEstimator(loss2) RobustLinearModel, A, b, est; ridgeλ=10, ridgeG=float([0 0; 0 1]), kwargs... ) m5 = fit(RobustLinearModel, A, b, est; ridgeλ=0.1, ridgeG=[0, 1], kwargs...) - m6 = rlm(A, b, est; ridgeλ=10, ridgeG=[0, 1], βprior=[0.0, 2.0], kwargs...) + m6 = rlm( + A, b, est; ridgeλ=0.1, ridgeG=[0, 1], dropcollinear=true, kwargs... + ) + m7 = rlm(A, b, est; ridgeλ=10, ridgeG=[0, 1], βprior=[0.0, 2.0], kwargs...) VERBOSE && println("\n\t\u25CF Estimator: $(name)") VERBOSE && println(" lm$(aspace) : ", coef(m1)) VERBOSE && println("rlm($(method)) : ", coef(m2)) - VERBOSE && println("ridge λ=1 rlm3($(method)): ", coef(m3)) - VERBOSE && println("ridge λ=1 rlm4($(method)): ", coef(m4)) + VERBOSE && println("ridge λ=10 rlm3($(method)): ", coef(m3)) + VERBOSE && println("ridge λ=10 rlm4($(method)): ", coef(m4)) VERBOSE && println("ridge λ=0.1 rlm5($(method)): ", coef(m5)) - VERBOSE && println("ridge λ=1 βprior=[0,2] rlm6($(method)): ", coef(m6)) + VERBOSE && println("ridge λ=0.1 dropcollinear=true rlm6($(method)): ", coef(m6)) + VERBOSE && println("ridge λ=10 βprior=[0,2] rlm7($(method)): ", coef(m7)) @test isapprox(coef(m3), coef(m4); rtol=1e-5) # Test printing the model @@ -56,6 +60,10 @@ est2 = MEstimator(loss2) refit!(m5; ridgeλ=10, initial_scale=:L1) @test isapprox(m5.pred.λ, 10.0; rtol=1e-5) @test isapprox(coef(m3), coef(m5); rtol=1e-6) + + refit!(m6; ridgeλ=10, initial_scale=:L1) + @test isapprox(m6.pred.λ, 10.0; rtol=1e-5) + @test isapprox(coef(m3), coef(m6); rtol=1e-6) end end From 64012b2658559b84dc7c6d0eef9061ac744a201f Mon Sep 17 00:00:00 2001 From: getzze Date: Fri, 2 Jun 2023 17:50:11 +0100 Subject: [PATCH 9/9] cleanup --- test/robustridge.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/test/robustridge.jl b/test/robustridge.jl index 5fd1673..fa9181d 100644 --- a/test/robustridge.jl +++ b/test/robustridge.jl @@ -38,9 +38,7 @@ est2 = MEstimator(loss2) RobustLinearModel, A, b, est; ridgeλ=10, ridgeG=float([0 0; 0 1]), kwargs... ) m5 = fit(RobustLinearModel, A, b, est; ridgeλ=0.1, ridgeG=[0, 1], kwargs...) - m6 = rlm( - A, b, est; ridgeλ=0.1, ridgeG=[0, 1], dropcollinear=true, kwargs... - ) + m6 = rlm(A, b, est; ridgeλ=0.1, ridgeG=[0, 1], dropcollinear=true, kwargs...) m7 = rlm(A, b, est; ridgeλ=10, ridgeG=[0, 1], βprior=[0.0, 2.0], kwargs...) VERBOSE && println("\n\t\u25CF Estimator: $(name)")