Loess.jl icon indicating copy to clipboard operation
Loess.jl copied to clipboard

Use cholesky decomposition whenever possible

Open ayushpatnaikgit opened this issue 1 year ago • 1 comments

Loess involves doing weighted least squares many times (at each vertex). The original implementation in Netlib, and the one in this package, use QR decomposition for this purpose.

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

This is the most expensive part of loess.

I propose using Cholesky factorization instead of QR decomposition.

coefs = cholesky(us' * us) \ (us' * vs)

The idea is that us' * us and us' * vs are pretty small.

Unfortunately, cholesky factorization fails when the matrix being factorized is not positive definite. Hence, I suggest using QR decomposition as a backup using an if condition. Try and catch also work, but I don't like using them.

Using my proposed method will lead to modest gains in computation time, but noticeable improvement in memory allocations. The numerical differences will between the two methods will be small. I did some tests (based on the example in the README) and the differences are < e-10. I think that's reasonable, given loess is an approximation.

Here's the performance comparison on the example in the README.

using BenchmarkTools, Loess
xs = 10 .* rand(100)
ys = sin.(xs) .+ 0.5 * rand(100)

Current implementation:

@btime loess(xs, ys, span=0.5)
164.083 μs (3573 allocations: 1.26 MiB)

Proposed:

@btime loess(xs, ys, span=0.5)
 140.478 μs (3080 allocations: 102.70 KiB)

As far as I know, this trick hasn't been tried before in the context of loess. While all the test cases pass (on Julia 1.9.0), we could be missing out on something. I've had several discussions on this with @ViralBShah and @mousum-github. Instead of opening an issue to discuss it further, I just did the PR, since it's a tiny change in terms of the number of lines.

ayushpatnaikgit avatar Jun 23 '23 13:06 ayushpatnaikgit

Codecov Report

Patch coverage: 100.00% and project coverage change: +0.14 :tada:

Comparison is base (1b34e35) 89.90% compared to head (2722659) 90.04%.

Additional details and impacted files
@@            Coverage Diff             @@
##           master      #72      +/-   ##
==========================================
+ Coverage   89.90%   90.04%   +0.14%     
==========================================
  Files           2        2              
  Lines         208      211       +3     
==========================================
+ Hits          187      190       +3     
  Misses         21       21              
Impacted Files Coverage Δ
src/Loess.jl 84.25% <100.00%> (+0.44%) :arrow_up:

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Do you have feedback about the report comment? Let us know in this issue.

codecov-commenter avatar Jun 23 '23 13:06 codecov-commenter