Julia - Least Squares - LSMR and closed-form solution different from true `\beta`?

user_6068731 0
Viewed 86

I wanted to try using the LSMR algorithm so I generated some data and run least squares. How come the LSMR solution and the closed-form one are different from the true $\beta$ I have used to generate the data?

using Distributions: Normal
using IterativeSolvers: lsmr
# Settings
n = 500
k = 50
# Generate data
X = rand(Normal(0.0, 0.1), (n, k)) + rand(Normal(0.0, 0.2), (n, k))
β = randn(k)
y = (X * β) + rand(Normal(0.0, 0.1), n)
# Solutions
closed_form_solution = (X'X) \ (X'y)
lsmr_solution = lsmr(X, y)
# Check solutions
β ≈ closed_form_solution, β ≈ lsmr_solution  # returns false, false
julia linear-regression linear-algebra
1 Answers

I just ran your code and visualised the values, this is what I got.

using Plots
scatter(β)
scatter!(closed_form_solution)
scatter!(lsmr_solution)

enter image description here

As you can see they're actually pretty close, so the algorithms work as expected.

Now they're not the same because you've added residuals and the algorithms don't know about that. They have to adapt to them. You can find only an approximation of the underlying model if one even exists and this is basically what machine learning tries to do.

If you generate y's like this

y = X * β

and compare again you'll get much closer results (≈ might still return false though because of floating-point errors).

user_128854680
Related