Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix QR based fitting #559

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Fix QR based fitting #559

wants to merge 5 commits into from

Conversation

andreasnoack
Copy link
Member

@andreasnoack andreasnoack commented May 7, 2024

Avoid erroring out for low rank design matrices when dropcollinear=false.
Avoid unnecessary triangular solves. Avoid indexing in the Q.
Avoid slicing R matrix in a way that triggers a minimum norm solution.
Remove unnecessary temporaries.

Fixes #558

I removed the white space so please review with https://github.com/JuliaStats/GLM.jl/pull/559/files?w=1

@andreasnoack andreasnoack force-pushed the an/norank branch 2 times, most recently from 33d1584 to 35c3e94 Compare May 7, 2024 11:56
Copy link

codecov bot commented May 7, 2024

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 89.95%. Comparing base (0e145bd) to head (2b61e90).

Additional details and impacted files
@@            Coverage Diff             @@
##           master     #559      +/-   ##
==========================================
- Coverage   90.03%   89.95%   -0.09%     
==========================================
  Files           8        8              
  Lines        1124     1105      -19     
==========================================
- Hits         1012      994      -18     
+ Misses        112      111       -1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

@andreasnoack andreasnoack requested review from devmotion and removed request for mousum-github May 17, 2024 09:30
src/linpred.jl Outdated Show resolved Hide resolved
Comment on lines +89 to +88
fill!(p.delbeta, 0)
p.delbeta[1:rnk] = R \ view(p.qr.Q'r, 1:rnk)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe one could save a few allocations here?

Suggested change
fill!(p.delbeta, 0)
p.delbeta[1:rnk] = R \ view(p.qr.Q'r, 1:rnk)
mul!(view(p.delbeta, 1:rnk), view(p.qr.Q, :, 1:rnk)', r)
ldiv!(R, view(p.delbeta, 1:rnk))
fill!(@view(p.delbeta[(rnk + 1):end]), 0)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think a view of a Householder matrix would end up hitting a fallback multiplicatino which would be very costly.

src/linpred.jl Outdated Show resolved Hide resolved
src/linpred.jl Show resolved Hide resolved
src/linpred.jl Outdated Show resolved Hide resolved
src/linpred.jl Outdated Show resolved Hide resolved
src/linpred.jl Outdated Show resolved Hide resolved
src/linpred.jl Outdated Show resolved Hide resolved
src/linpred.jl Outdated Show resolved Hide resolved
src/linpred.jl Outdated
end
p.delbeta[1:rnk] = R \ (p.qr.Q'*r̃)[1:rnk]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe

Suggested change
p.delbeta[1:rnk] = R \ (p.qr.Q'*r̃)[1:rnk]
mul!(view(p.delbeta, 1:rnk), view(p.qr.Q, :, 1:rnk)', r̃)
ldiv!(R, view(p.delbeta, 1:rnk))

?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can't make a view into Q, see above, but we can make a view into the result.

Avoid erroring out for low rank design matrices when dropcollinear=false.
Avoid unnecessary triangular solves. Avoid indexing in the Q.
Avoid slicing R matrix in a way that triggers a minimum norm solution.
Remove unnecessary temporaries.
Copy link
Member Author

@andreasnoack andreasnoack left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've addressed your comments

src/linpred.jl Outdated Show resolved Hide resolved
src/linpred.jl Show resolved Hide resolved
p.delbeta = Rinv * Rinv' * p.delbeta
ỹ = sqrtW * r
p.qr = qr!(p.scratchm1)
p.delbeta = p.qr \ ỹ
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

src/linpred.jl Outdated Show resolved Hide resolved
Comment on lines +89 to +88
fill!(p.delbeta, 0)
p.delbeta[1:rnk] = R \ view(p.qr.Q'r, 1:rnk)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think a view of a Householder matrix would end up hitting a fallback multiplicatino which would be very costly.

src/linpred.jl Outdated Show resolved Hide resolved
src/linpred.jl Outdated Show resolved Hide resolved
src/linpred.jl Outdated Show resolved Hide resolved
src/linpred.jl Outdated
end
p.delbeta[1:rnk] = R \ (p.qr.Q'*r̃)[1:rnk]
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can't make a view into Q, see above, but we can make a view into the result.

@andreasnoack andreasnoack force-pushed the an/norank branch 2 times, most recently from dcc9e21 to 37249f5 Compare November 5, 2024 15:06
@andreasnoack
Copy link
Member Author

@devmotion I finally found a way to download the test dataset reliably so this one should be good to go.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

GLM fails to match certified values for the Filip dataset
2 participants