• (2014?) ftp://ftp.math.ucla.edu/pub/camreport/cam12-35.pdf -- heavily cited

• (2015?) seems more general https://www-users.cs.umn.edu/~baner029/papers/15/A2DM2_kdd15.pdf

• http://proceedings.mlr.press/v37/nishihara15.pdf

@tlienart . Sorry to disturb you again. I think a more informative error message should be thrown.
I ran LogisticClassifier using the default values on a dataset having nclasses > 2. it thows the following Error

ERROR: DimensionMismatch("new dimensions (4, 4) must be consistent with array size 4")

Although setting multi_class=true works. i think there should be a more informative error message.

Here the changes since 0.1.8, which I doubt break anything in MLJLinearModels:

• Bring in is_same_except and similar gizmos from MLJBase; see JuliaAI/MLJBase.jl#178

• Let all models share same traits (#19, JuliaAI/MLJBase.jl#163)

• Add save and restore stubs for serialisation

The only issue is that would need to form X_ but can offer the option.

(see readme)

• Ridge classification
• LARS see LARS.jl
• Quantile regression see [QuantileRegression.jl] (https://github.com/pkofod/QuantileRegression.jl) see also #4 and #5
• glmnet see GLMNet wrapper
• prox Newton or L-BFGS (see also https://github.com/yuekai/PNOPT; and paper https://arxiv.org/pdf/1206.1623.pdf as well as https://papers.nips.cc/paper/5384-proximal-quasi-newton-for-computationally-intensive-l1-regularized-m-estimators.pdf)
• count -- http://hua-zhou.github.io/media/pdf/ZhangZhouZhouSun17mglm.pdf
• Ridge classifier -- https://stackoverflow.com/questions/53911663/what-does-sklearn-ridgeclassifier-do/53912015
• TheilSen -- https://en.wikipedia.org/wiki/Theil–Sen_estimator

Other stuff

• add option for size of tape in LBFGS solver(s)
• add option(s) for linesearches
• add and test options generally for solvers

It shouldn't be; this may be due to issues with the Krylov Solver; consider rewriting your own NewtonCG using the cg from IterativeSolvers & starting with a simple backtracking; maybe later with Hager-Zhang.

at https://github.com/alan-turing-institute/MLJLinearModels.jl/blob/master/src/mlj/interface.jl#L22 the function schema(X) fails. @tlienart

also use Symmetric for good measure

See this discussion for some insight on things https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/26

• Replace scratch space by allocations in the factory functions
• pass whatever is needed to apply_Xt etc

this will likely make the code a bit less readable in some places but is probably better practice and shouldn't affect performance, additionally it should help avoid non-thread safe stuff

To bring MLJLinearModels into synch with the other repos.

(see branch docs)

• the fact that, when we can, we avoid doing hcat(X, 1)
• use IterativeSolvers.cg and not IterativeSolvers.lsqr as it allows specifying the operator as a linear map which is efficient and avoids copying when having to add a column for X; anyway it should be identical apart from pathological cases
• robust regression with concomitant scale estimation, not done yet
• LAD/quantile explain weight clipping for pathological zero residual
• generally detail how algos do IWLS etc so that it's 100% clear what is being done
  • Newton
  • NewtonCG (explain briefly then link to Optim)
  • LBFGS (explain briefly then link to Optim)
  • IWLSCG

What happens if input data is entered as integer? as float32, how about the target? verify that it's consistent (at least with float64)

In case users prefer this;

that means though that things should be ScaledLosss{...} in fit methods etc, it shouldn't change much otherwise.

see https://github.com/JuliaLang/julia/pull/31853/files#diff-3a896b2398ad83022a6e197921c125c6R274-R277

Though it may be overkill to compute all of the theta for all lambdas and better to get a proxy for the GCV error and then compute just for that

Related code

• (depr) SGDOptim.jl discussed in here
• (tutorial) sgd.jl

One day there may be a refactoring to use LossFunctions, in the mean time it's irrelevant but anyway:

-( 1 - tau) * r if r < 0
tau * r otherwise

https://www.stat.berkeley.edu/~mmahoney/pubs/quantile-icml13.pdf

For 1/2 it's basically l1

Hi,
I faced this issue when using the multinomial regression (aka multiclass Logistic Regression).
In the following, I show an example with MNIST classification.

I have to shift all categories of 1 integer (0->1, 1->2 .... 9->10) otherwise it gave the error I report in the end notes

# Classify MNIST digits with a simple multi-layer-perceptron
using Flux.Data.MNIST
using MLJLinearModels
# Get MNIST dataset and transpose for (records, features)
imgs = MNIST.images()
X = Array(transpose(hcat(float.(reshape.(imgs, :))...) )
# MNIST labels: Categorical labels must be 1...c, hence add .+1 to each label
labels = MNIST.labels() .+1
# and the number of classes
n_classes = length(Set(labels))
n_features = size(X,2)
# The MNIST database does not need the intercept
intercept = false
# deploy MultinomialRegression from MLJLinearModels, λ being the strenght of the reguliser
mnr = MultinomialRegression(λ; fit_intercept=intercept)
# Fit the model
θ  = fit(mnr, X, labels)
# The model parameters are organized such we can apply X⋅θ, the following is only to clarify
params = reshape(θ, n_features +Int(intercept), n_classes)
# Get the predictions X⋅θ 
preds = MLJLinearModels.softmax(MLJLinearModels.apply_X(X,θ,n_classes))
# map each vector to its maximal element 
targets = map(x->argmax(x),eachrow(preds))
#and evaluate the model over the labels
scores = 1- sum(targets-labels)/length(preds)

Error with 0 among the categories:

DimensionMismatch("new dimensions (785, 10) must be consistent with array size 7056")
(::Base.var"#throw_dmrsa#197")(::Tuple{Int64,Int64}, ::Int64) at reshapedarray.jl:41
reshape at reshapedarray.jl:45 [inlined]
reshape at reshapedarray.jl:116 [inlined]
apply_X!(::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,1}, ::Int64) at utils.jl:66
(::MLJLinearModels.var"#102#103"{GeneralizedLinearRegression{MultinomialLoss,ScaledPenalty{LPPenalty{2}}},Array{Float64,2},Array{Int64,1},Int64,Int64,Int64,Float64})(::Float64, ::Array{Float64,1}, ::Array{Float64,1}) at d_logistic.jl:149
(::NLSolversBase.var"#61#62"{NLSolversBase.InplaceObjective{Nothing,MLJLinearModels.var"#102#103"{GeneralizedLinearRegression{MultinomialLoss,ScaledPenalty{LPPenalty{2}}},Array{Float64,2},Array{Int64,1},Int64,Int64,Int64,Float64},Nothing,Nothing,Nothing},Float64})(::Array{Float64,1}, ::Array{Float64,1}) at incomplete.jl:45
value_gradient!!(::NLSolversBase.OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}) at interface.jl:82
initial_state(::Optim.LBFGS{Nothing,LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Optim.var"#19#21"}, ::Optim.Options{Float64,Nothing}, ::NLSolversBase.OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}) at l_bfgs.jl:158
optimize(::NLSolversBase.OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}, ::Optim.LBFGS{Nothing,LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Optim.var"#19#21"}, ::Optim.Options{Float64,Nothing}) at optimize.jl:33
#optimize#93 at interface.jl:116 [inlined]
optimize(::NLSolversBase.InplaceObjective{Nothing,MLJLinearModels.var"#102#103"{GeneralizedLinearRegression{MultinomialLoss,ScaledPenalty{LPPenalty{2}}},Array{Float64,2},Array{Int64,1},Int64,Int64,Int64,Float64},Nothing,Nothing,Nothing}, ::Array{Float64,1}, ::Optim.LBFGS{Nothing,LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Optim.var"#19#21"}, ::Optim.Options{Float64,Nothing}) at interface.jl:115
optimize(::NLSolversBase.InplaceObjective{Nothing,MLJLinearModels.var"#102#103"{GeneralizedLinearRegression{MultinomialLoss,ScaledPenalty{LPPenalty{2}}},Array{Float64,2},Array{Int64,1},Int64,Int64,Int64,Float64},Nothing,Nothing,Nothing}, ::Array{Float64,1}, ::Optim.LBFGS{Nothing,LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Optim.var"#19#21"}) at interface.jl:115
_fit(::GeneralizedLinearRegression{MultinomialLoss,ScaledPenalty{LPPenalty{2}}}, ::LBFGS, ::Array{Float64,2}, ::Array{Int64,1}) at newton.jl:114
#fit#144(::LBFGS, ::typeof(fit), ::GeneralizedLinearRegression{MultinomialLoss,ScaledPenalty{LPPenalty{2}}}, ::Array{Float64,2}, ::Array{Int64,1}) at default.jl:48
fit(::GeneralizedLinearRegression{MultinomialLoss,ScaledPenalty{LPPenalty{2}}}, ::Array{Float64,2}, ::Array{Int64,1}) at default.jl:38
top-level scope at test_LR.jl:159

• All regressions have residuals, there could be metrics on the residuals. In most cases it would make more sense to just implement them in MLJBase and use them here
• FWD/BWD selection procedure is a bit dumb and expensive but people use it I guess.

basically just re-scale the calls by median(r). This may be what matlab does though a bit unclear.

• https://stat.ethz.ch/R-manual/R-patched/library/MASS/html/rlm.html

Expose the functions allowing to to create examples (continuous, binary etc)

• take a look at https://github.com/FugroRoames/RobustLeastSquares.jl

http://julien.mairal.org/cyanure/welcome.html

No other action required. (The MLJBase API has not changed since 0.8 as far as MLJLinearModels is concerned.)

At present, if the update is not executed, MLJLinearModels will, it seems, disappear in the next @update to the MLJ registry.

Quantile, LAD regression

• ADMM should work but if requires adapting rho then needs refactoring of H often which is wasteful, could imagine doing CG for that bit but that would also end up being expensive see also #8
• MM and other algorithms see issue #3 and #4
• IP Method with or without pre-proc ox-code
• Frisch-Newton (for L1 reg may be good) ref

NewtonCG, IWLSCG

• should add a field where the user can specify :cg or :minres or something else, unlikely it would make a big difference in perf though.

it should be possible to just pass the scratch spaces to the factory.

At the moment indirectly tested but they may actually be incorrect but look correct.

Before starting this, need a way to systematically:

• trace number of function calls, number of gradient calls, number of hessian calls
• have a way to stop with universal criterion OR show a plot where the objective function decreases and eventually hits the same value as that from ref package.

Against scikitlearn

expect on par or better

• ridge (in big case should see improvements from using CG)
  • analytical (should see no real diff)
  • CG
• lasso
  • FISTA
  • ISTA
• elnet
• logistic (no or l2 penalty)
  • (cannot test Newton)
  • Newton CG see issue #24
  • LBFGS
• logistic (elnet penalty)
  • FISTA
  • ISTA
• multinomial (no or l2 penalty)
  • Newton CG see issue #24
  • LBFGS
• multinomial (elnet penalty)
  • FISTA
  • ISTA

Against quantreg

expect a bit worse (quantreg is effectively in cpp)

• quantile regression

For consistency with MLJ

use JuliaNLSolvers/Optim.jl#742 in d_logistic.jl Hv!

It still lives on the branch admm-bk, should probably refactor thinking:

• if adapting rho need to refactor H
• probably should check on a simple problem first to see if works
• simple matlab code directly translated does not work

https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232

Ridge

refs

• https://arxiv.org/pdf/1509.09169.pdf
• GCV and GCVPack (probably a good place to start): http://pages.stat.wisc.edu/~wahba/stat860public/pdf1/775r.pdf
• page 22 https://elib.uni-stuttgart.de/bitstream/11682/9236/1/BA_fertig_QianKun_25072017.pdf
• slides https://www.wias-berlin.de/workshops/validation2010/workshop_slides/talk_goeman.pdf

Basically just have to minimise a criterion that can be evaluated in constant time after the first fit.

Similar stuff can be obtained for LOOCV it seems

proper refs

• http://w3.atomki.hu/~efo/hornyak/Tikhonov_references/Technometrics_Golub_Heath_Wahba

Logistic: https://arxiv.org/pdf/1711.05420.pdf

Context: alan-turing-institute/MLJ.jl#387

Fix: allow symbols to be passed around and just re-interpret them as string.

@tlienart

1. All Subset (Best K Subset):
   I've got parsimonious code for All Subset Regression which might have a nice home in this package that I'd like to share.
   While it's much faster than R's Leaps pkg, I would include a strong warning to users:
   This optimization problem is sooo not convex. Avoid for p>20.

2. Approximate All Subset (MIO):
   As you probably know,
   Bertsimas etal (2016): show how to use MIO to approximately train all subset regression for p =1000 in minutes.
   Bertsimas et al (2019): do this for p=100k in minutes.
   The corresponding Julia package is SubsetSelection.jl.
   Also see forthcoming Julia package MIP-Boost.
   It would be really awesome if these linear models were featured in MLJ.

3. Relaxed Lasso:
   Hastie et al 2017: responded that relaxed Lasso matches Lasso in low-SNR scenarios & beats All Subset in high-SNR scenarios
   glmnet: now features "relaxed Lasso"
   Would this be difficult to add to MLJLinearModels? (maybe we can work on it together?)

4. Linear Interactions:
   There are some cool new packages for interactions in linear models.
   I find these compare in out-of-sample fit w/ non-linear models (xgboost etc, no joke)
   R: sprintr & HierNet
   Python: HierScale

5. Nonconvex penalties (MCP SCAD):
   https://github.com/joshday/SparseRegression.jl/blob/master/README.md

6. Forward selection/stepwise regression:
   (stackoverflow, has Julia code)
   Textbook code

7. Ensembles of Penalized Linear Models:
   Christidis et al 2019: Paper, old code EnsembleEN, new code SplitReg

Currently, the fitted parameters (coefficient and intercept) for LinearRegressor does not include feature names, which makes mapping coefficients back to their input columns complicated when combined with transformations such as OneHotEncoder. E.g.

julia> fps = fitted_params(selected_model)

julia> fps.machines
3-element Array{Any,1}:
 NodalMachine{UnivariateStandardizer} @ 5…02
 NodalMachine{LinearRegressor} @ 6…58
 NodalMachine{OneHotEncoder} @ 8…55

julia> fps.fitted_params_given_machine[fps.machines[2]]
(coefs = [XXX],
 intercept = XXX,)

I think this would be a useful diagnostic feature to have something like

julia> fps.fitted_params_given_machine[fps.machines[2]]
(coefs = [XXX],
 intercept = XXX,
 names = [XXX])

Hey looks like a nice library. I would appreciate an example of a logistic regression. Also i was wondering what 0/L2 and L1/EN mean, i don't know which one to choose. Thank you

This is micro-optimization for later on but the functions that are called often like fgh! and the functions therein, could probably have a cache, would have to see if that's digestible to Optim but for instance it seems wasteful to allocate X*theta or the associated residual every time it's called.

alan-turing-institute/MLJ.jl#430

algos

• Maj Min algorithm for non-reg case (http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.206.1351&rep=rep1&type=pdf, also https://arxiv.org/pdf/1709.04126.pdf)
• ADMM for gen reg case (https://arxiv.org/pdf/1709.04126.pdf)
• L1 penalised: Frisch Newton http://www.econ.uiuc.edu/~roger/research/sparse/fn3.pdf
• non-reg also IRWLS http://users.stat.umn.edu/~sandy/courses/8053/handouts/robust.pdf

ADMM, MM, and CD approaches, while MM and CD are faster and ADMM slower than the IP algorithm available in quantreg. The results so far suggest that the MM algorithm is the best-suited for non-regularized (composite) quantile regression among the four methods tested, especially for data sets with large n and relatively small p. In regularized quantile regression, all methods perform similarly in terms of variable selection, but CD and ADMM show clear superiority in run time, particularly relative to the IP and MM methods when p is large. In the case of regularized composite quantile regression, CD and ADMM dis

(...)

Applying existing optimization algorithms to (composite) quantile regression requires a non-trivial reformulation of the problem due to the non-linearity and non-differentiability of the loss and regularization terms of the objective. The well-known quantreg package for R (Koenker, 2017) uses an interior point (IP) approach for quantile and composite quantile regression with the option of l1 (lasso) regularization for the former and no regulariza- tion options for the latter. Although advanced IP algorithms in quantreg, such as the one using prediction-correction (Mehrotra, 1992) for non-regularized quantile regression, have greatly improved upon earlier attempts using simplex methods, the time spent on matrix inversion in IP approaches (Chen and Wei, 2005) motivates us to seek faster algorithms for quantile and composite quantile regression, particularly for high-dimensional data where regularization is required. In addition, following the conjectures of Fan and Li (2001), Zou (2006) showed lasso variable selection—currently the most commonly-implemented penalty for quantile regression—to be inconsistent in certain situations and presented adaptive lasso regularization as a solution. Our work in the present paper is thus motivated by both a search for faster quantile regression algorithms as well as the lack of publicly-available meth- ods for adaptive-lasso regularized quantile and composite quantile regression, particularly for high-dimensional data.

in https://arxiv.org/pdf/1709.04126.pdf (cqreg package in R)

refs

• Hunter & Lange: MM algo for QR http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.206.1351&rep=rep1&type=pdf

References

• scad, mcp, adaptive lasso ref
• trimmed l1 ref

Could not get it to work well on LAD but maybe can make it work on something else like LASSO and can test the rest of the implementation on the go.

In the mean time maybe isolate the code in another branch as it makes code cov look bad.

A PkgEval run for a PR which changes the generated numbers for randn! indicates that the tests of this package might fail in Julia 1.5 (and on Julia current master). Apologies if this is a false positive.
cf. https://github.com/JuliaCI/NanosoldierReports/blob/7de24e455342298cbef56826b5827f0d7640d2c1/pkgeval/by_hash/b89e35c_vs_098ef24/logs/MLJLinearModels/1.5.0-DEV-71a4a114c2.log

• https://statweb.stanford.edu/~jtaylo/courses/stats203/notes/robust.pdf note discussion of scale

admm

• https://web.stanford.edu/~boyd/papers/admm/least_abs_deviations/lad.html
• https://web.stanford.edu/~boyd/papers/pdf/admm_distr_stats.pdf (page 39)

could probably also benchmark against sequence of Huber estimators with delta going to zero

check allocations and speed of elementary step such as

• apply_X!
• apply_Xt!
• get_residuals!
• ...