 # Inverse and backslash \ division

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
``````

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

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

2 Likes

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.

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

6 Likes

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?

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.