Just to elaborate on John's correct answer: eye(size(Z,1)) is a
Matrix{Float64} so before the inversion, Z'Z + eye(size(Z,2)) is promoted
to a Matrix{Float64} even though Z'Z is a Matrix{Float32}. In contrast,
inv(Z'Z) is inverted as a Matrix{Float32}.
Two other observations:
1. Unless you apply the same inverse to many different right hand sides, it
is extremely wasteful to compute the inverse and then multiply. Instead,
you just want to write (Z'Z)\(Z'y). Try to make some benchmarks to see the
difference. Some people also say that \ is more precise, but the
significance of that is heavily exaggerated.
2. It is much better to use I instead eye, i.e. (Z'*Z + λ*I)\Z'*y. That
avoids a lot of allocation.
On Tuesday, March 1, 2016 at 7:33:32 PM UTC-5, John Myles White wrote:
>
> Compare and contrast the following:
>
>
> inv(map(Float32, Z'*Z))*(Z'*ys_dataset)
>
> inv(map(Float64, Z'*Z))*(Z'*ys_dataset)
>
> -- John
>
> On Tuesday, March 1, 2016 at 4:19:19 PM UTC-8, Francisco C Pereira wrote:
>>
>> Hi,
>>
>> I deeply apologize in advance if this is a repeated question... the
>> question is relatively simple. How can this happen? I was estimating a
>> basic ridge regression model but ran into strange results:
>>
>> lambda=0
>> beta_hat=inv(Z'Z+lambda*eye(size(Z,2)))*(Z'*ys_dataset)
>> beta_hat2=inv(Z'*Z)*(Z'*ys_dataset)
>> beta_hat-beta_hat2
>>
>> Output (it should be a vector with zeroes):
>> 21-element Array{Float64,1}:
>>
>> -245.183
>> -436.175
>> 2633.22
>> 73.0047
>> -3345.18
>> 1712.78
>> -354.545
>> 187.102
>> 49.8899
>> -49.4564
>> -2.91745
>> -0.253636
>> -0.186349
>> 0.457511
>> -0.031112
>> 0.0789713
>> -0.00130604
>> -0.00448352
>> -0.00324767
>> 0.000248121
>> 0.000152718
>>
>>
>>
>>
>> Just for reference, here are the other variables:
>>
>>
>> xs_dataset=[-0.8,-0.6,-0.4,-0.2,0.0,0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.2]
>>
>> ys_dataset=[-194.74997124179467,-218.7737153673673,-262.705773686837,-151.42110150479215,-147.91115080458417,-208.4047152970189,-128.36970674967745,-114.53430386917461,-158.27662488829077,-142.21198729962245,-82.80260143610138,-27.621822591723756,-25.154050675997198,-6.99434298775955,-0.06912672996766567,-5.349114073043877,32.43503138931636,-37.228578783587075,-37.22806830522013,-113.28447052959032,-71.0209779039439]
>>
>> Z=Array{Float32}([z^i for z in xs_dataset, i in 0:20])
>>
>>
>>
>> I've been arguing that Julia is the best in the world, I need help to keep
>> preaching! ;-)
>>
>>
>> Thanks!
>>
>> Francisco
>>
>>
