least squares implementation

[tpapp]

AFAICT the implementation of linear regression forms the X'X matrix then calculates (X'X)⁻¹(X'y) using Cholesky factorization.

I am curious why this was chosen. The the best of my knowledge, orthogonal (QR) methods are more stable (if a bit more costly), and of course SVD is the best solution for nearly rank-deficient matrices.

[andreasnoack]

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.

[tpapp]

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.

[andreasnoack]

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.

[tpapp]

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.

[andreasnoack]

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.

[tpapp]

I think you are correct about the ML focus. Perhaps an explanation of favoring speed over accuracy in the README would clarify this.

[Nosferican]

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.