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.

Sign up for GitHub

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

Jump to bottom

Use cholesky decomposition whenever possible #72

Draft
wants to merge 3 commits into
base: master
Choose a base branch
from
Draft

Use cholesky decomposition whenever possible #72

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a7ed326
Use cholesky decomposition whenever possible.
ayushpatnaikgit Jun 23, 2023
ddd915a
Use exception handling
ayushpatnaikgit Jun 23, 2023
2722659
Minor correction
ayushpatnaikgit Jun 24, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 11 additions & 6 deletions src/Loess.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,13 +122,18 @@ function loess(
xl *= x
end
end

if VERSION < v"1.7.0-DEV.1188"
F = qr(us, Val(true))
else
F = qr(us, ColumnNorm())

coefs = Matrix{Real}(undef, size(us, 2), 1)
try
coefs = cholesky(us' * us) \ (us' * vs)
Copy link
Contributor

@ViralBShah ViralBShah Jun 23, 2023 •
edited
Loading

Choose a reason for hiding this comment

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

With isposdef this way, you are doing cholesky twice. You really should only do it once, and check the error and if not positive definite, fall back to QR.

Copy link
Contributor

Choose a reason for hiding this comment

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

Also is this how people do it in literature? Forming the product squares the condition number, right?

Copy link
Author

@ayushpatnaikgit ayushpatnaikgit Jun 23, 2023 •
edited
Loading

Choose a reason for hiding this comment

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

With isposdef this way, you are doing cholesky twice. You really should only do it once, and check the error and if not positive definite, fall back to QR.

Do you mean like this:

coefs = Matrix{Real}(undef, size(us, 2), 1)
try
    coefs = cholesky(us' * us) \ (us' * vs)
catch PosDefException
    if VERSION < v"1.7.0-DEV.1188"
        F = qr(us, Val(true))
    else
        F = qr(us, ColumnNorm())
    end
    coefs = F \ vs
end

Copy link
Author

Choose a reason for hiding this comment

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

Also is this how people do it in literature? Forming the product squares the condition number, right?

@mousum-github what do you think?

Copy link
Contributor

Choose a reason for hiding this comment

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

Yes, I would catch the exception and trigger the QR. I believe doing QR should be numerically better than squaring and doing Cholesky, in general. @andreasnoack ?

Copy link
Author

@ayushpatnaikgit ayushpatnaikgit Jun 24, 2023 •
edited
Loading

Choose a reason for hiding this comment

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

This is indeed better.
Here are benchmarks on a different xs & ys.
Current:

@btime loess(xs, ys, span=0.5)
  150.481 μs (3037 allocations: 1.25 MiB)

If condition

@btime loess(xs, ys, span=0.5)
  125.853 μs (2544 allocations: 94.33 KiB)

Exception handling

@btime loess(xs, ys, span=0.5)
  116.999 μs (2527 allocations: 91.41 KiB)

Copy link

@mousum-github mousum-github Jun 25, 2023 •
edited
Loading

Choose a reason for hiding this comment

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

Also is this how people do it in literature? Forming the product squares the condition number, right?

@mousum-github what do you think?

@ayushpatnaikgit - Please have a look at these examples.
JuliaStats/GLM.jl#507 (comment)
In summary, up to a certain condition number of the design matrix say 1e4, both Cholesky and QR produce the same estimates, and then up to another big condition number of the design matrix say 1e7, Cholesky produces less accurate estimates whereas QR produces accurate estimates. Then up to another big condition number say 1e15, Cholesky raises PosDefException, whereas QR still produces accurate estimates. And when the condition number is above say 1.e16 both Cholesky and QR detect rank deficiency.
As per the above example, the TRY CATCH statement works well when the condition number of the design matrix is <= 1e4 or > 1e8 but it fails to produce accurate estimates when the condition is in between (1.e4, 1.e8).

Copy link
Author

Choose a reason for hiding this comment

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

@ViralBShah should we add a cholesky/qr key word argument, just as we have done in GLM.jl?

Copy link

Choose a reason for hiding this comment

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

The default will be Cholesky, and we need to have some examples showing when we should use QR.

Copy link
Author

Choose a reason for hiding this comment

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

Thanks, Mousum. I'll figure out a case where cholesky doesn't work, and then we can take the discussion from there.

catch PosDefException
if VERSION < v"1.7.0-DEV.1188"
F = qr(us, Val(true))
else
F = qr(us, ColumnNorm())
end
coefs = F \ vs
end
coefs = F\vs

predictions_and_gradients[vert] = [
us[1, :]' * coefs; # the prediction
Expand Down