Comments (7)
Using the Cholesky is typically much faster and most often the statistical error is much larger than error due to roundoff so in many cases it can preferable with the faster version.
from multivariatestats.jl.
No, Cholesky is notoriously bad. No serious library I know of (in other languages) uses it. It is very easy to come up with examples like
using MultivariateStats
function makeX(ϵ, N = 10)
X = zeros(N, 3)
k = N ÷ 2
X[:, 1] = 1.0
X[:, 2] = 1:N
X[1:k, 3] = 1:k
X[(k+1):N, 3] = ((k+1):N) + ϵ
X
end
X = makeX(1e-6)
β = Float64.(1:3)
Y = X * β
β1 = llsq(X, Y; bias = false)
β2 = X \ Y
norm(β-β1, Inf) # ouch, 0.2
norm(β-β2, Inf) # 1e-10
but the worse thing is that this phenomenon happens a lot in practice, if you have a few tens of covariates.
from multivariatestats.jl.
Well, in your example you chose zero statistical error so, of course, the errors due to rounding will dominate. If you data matrix is practically singular and your explained variable is in the range of the data matrix then sure, use QR, but if one the conditions doesn't hold then it typically doesn't matter much and the speedup can be quite significant.
No serious library I know of (in other languages) uses it.
What is the basis of this statement? SAS's proc reg
and SPARK's LinearRegression
class are both using the normal equations and both are used for quite serious work.
from multivariatestats.jl.
I am not sure I understand you argument. Are you saying that because estimation has statistical error anyway, we might as well ignore potentially significant numerical error?
I can't check SAS because it is closed source, and I don't have time to dive into SPARK's source (it would be great if you could provide a link to the relevant part, since you seem to be familiar with it). But R uses QR, of course Julia uses QR, GSL uses SVD (probably a bit overcautious).
You are of course correct about speed, MultivariateStats.llsq
is about 4-6x faster in my benchmarks (depending on matrix size). However, given that this is a general use statistical library, I think it would be better to make the default correct, and perhaps provide a fast version for users who understand the trade-off (this reminds me of a recent thread about --fast-math on Discourse).
My understanding is that this is pretty standard, numerical analysis textbooks always caution about following the naive formulas for regression, and recommend at least QR. But if your mind is set about this, feel free to close this issue, I will just write my own routines.
from multivariatestats.jl.
Are you saying that because estimation has statistical error anyway, we might as well ignore potentially significant numerical error?
I think the fear of rounding errors of least squares is exaggerated in statistical applications. If your design matrix is so ill-conditioned then your estimates will, in general, be so uncertain that the rounding errors don't matter. Do you often have data where the outcome is in the range of the regressors, i.e. where the residual is zero? I think it can happen in experimental settings but in an experimental setting, there is not really a good excuse for having an ill-conditioned design matrix.
Indeed, I should have provided some links. Here is a brief description the computational method of prog reg in SAS and if I read the source correctly then this thing here is called when there is no regularization in SPARK's regression.
But if your mind is set about this, feel free to close this issue, I will just write my own routines.
I don't have strong feelings about this package (I mainly use GLM.jl
or just \
if I don't need standard errors) but, as mentioned, I think the fear of rounding errors of least squares is exaggerated so I think you should consider if this is really a real problem for the data that your work with before rolling your own. It would probably also be fine to add a QR based version here which could be enabled with e.g. a keyword. However, this package is a bit ML focussed and in ML, the focus seems to be on selling precision for speed.
from multivariatestats.jl.
I think you are correct about the ML focus. Perhaps an explanation of favoring speed over accuracy in the README would clarify this.
from multivariatestats.jl.
GLM.jl offers both QR
and Cholesky
. pinv
uses SVD for the pseudo-inverse. Other alternatives are to use IterativeSolvers.jl. I believe Cholesky to be probably the best choice virtually always. qr
is not that great when you have few features, but large number of observations. For so ill conditioned problems, I would suggest using Ridge with a small penalization.
from multivariatestats.jl.
Related Issues (20)
- Any plan for functional PCA? HOT 1
- 1.0 Release
- `transform` and `reconstruct` HOT 1
- Request for adjoint support in MulticlassLDA HOT 4
- show function calls show(io,M), not show(io, MIME, M) HOT 1
- error while fitting PCA model HOT 2
- LDA docs clarification HOT 5
- KernelPCA reconstruct problem HOT 1
- fitting `SubspaceLDA models fails.
- `predict` for `MulticlassLDA` is throwing error HOT 1
- Mutating labels when training MulticlassLDA HOT 1
- Isotonic Regression question HOT 1
- MethodError: no method matching fit(::Type{MDS}, ::Adjoint{Float64, Matrix{Float64}} HOT 3
- Add Citation.bib and DOI with zenodo?
- Bugs in isotonic regression ? HOT 1
- CCA results do not match sklearn.cross_decomposition.CCA HOT 1
- Can I fix the number of number of factors (i.e., output dimension)?
- Inferred `outdim` for ICA sometimes fails but without an informative explanation
- LDA Printing flipped in y and x axis data on docs
- Hope to add umap
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from multivariatestats.jl.