Inverse and backslash \ division

Specific Domains Statistics
,
1

Hi,

could anyone please explain to me how to use inv() and \ properly to have valid results. Alse when checking with GLM package basic OLS gives different results and sometimes it is way too far. It might be caused by the data but strangely.

Thanks a lot.

x = randn(100, 3)
y = randn(100)

julia> (transpose(x) * x) \ (transpose(x) *y) .== inv(transpose(x) * x) * (transpose(x) *y)
3-element BitArray{1}:
1
0
1

julia> (transpose(x) * x) \ (transpose(x) *y)
3-element Array{Float64,1}:
-0.13268579177064818
-0.0009498174455286444
0.06340639348585342

julia> inv(transpose(x) * x) * (transpose(x) *y)
3-element Array{Float64,1}:
-0.13268579177064818
-0.0009498174455286456
0.06340639348585342

julia> using GLM

julia> (transpose(x) * x) \ (transpose(x) *y)
3-element Array{Float64,1}:
-0.13268579177064818
-0.0009498174455286444
0.06340639348585342

julia> GLM.coef(lm(x, y))
3-element Array{Float64,1}:
-0.1326857917706482
-0.0009498174455286441
0.06340639348585342

julia> (transpose(x) * x) \ (transpose(x) *y) .== GLM.coef(lm(x, y))
3-element BitArray{1}:
0
0
1
2

Instead of checking for exact equality with ==, I would recommend using approximate equality. See isapprox.

You should be able to just replace the == with ≈.

See: Mathematics · The Julia Language

2 Likes
3

Thanks for an answer, but in my case it is more about the results, which, for data I have, lead to exploding results in OLS because of definiteness of inv(X^T X) matrix etc.

Then the difference is super significant not only on 6th+ decimal points but on 2nd, for example. And values really differ when using inv or \.

Any ideas? Thanks.

4

It sounds like you are solving a least-squares problem by the “normal-equations” approach, i.e. you are computing z = (X^T X)^{-1} X^T y to minimize the least-square error \Vert X z - y \Vert_2. If you are getting “exploding” results, it is presumably because X^T X is badly conditioned (nearly singular), in which case you have to be careful about how you solve the problem — your answer is very sensitive to roundoff errors and other errors (e.g. noise in the data).

If X^T X is badly conditioned, then inv(X'X) * (X'y) could indeed give very different answers than (X'X) \ (X'y) because of differences in roundoff errors — but the correct interpretation is that both approaches are giving you an inaccurate result.

A better approach would be to use X \ y, which is equivalent to qr(X, Val(true)) \ y — this uses a pivoted QR factorization to solve the least-squares problem, which is much more accurate than the normal-equations approach for badly conditioned X.

If you need even more control over regularization for a badly conditioned problem, you could use a pseudo-inverse (pinv) approach. See also my explanation in another thread: Efficient way of doing linear regression - #33 by stevengj

7 Likes
5

Thanks. This I do understand. And will try to use QR factorization.

Also, I will add some information/comment to the problem I’m looking at. It is an analysis where I many perform OLS on expanding window over data. And the thing I have realised is that on 1:1000 observations OLS is ok, and when window is 1:1001, the coefficients explode. This is strange since the window is expanding and 1 observation pushes the analysis into unstable solution and when I move one observation further, 1:1002, it is stable again. So the matrix XTX is probably on the edge of well/bad conditioning, or I do not have any better explanation. Does this make sense in accordance to what you’ve written?

6

I have come up with a more efficient method than the SaticArrays method for an underdetermined system:

For underdetermined systems, my new method is about ~20x faster than the SMatrix and it is numerically stable. However, for overdetermined cases, I am falling back on the traditional normal equations, which are less numerically stable.