ForwardDiff.jl
ForwardDiff.jl copied to clipboard
Faster Matrix{BlasFloat} * or \ VecOrMatrix{Dual}
The calculation of the gradient from a multivariate normal prior involves left-division by a triangular matrix. Currently, when the covariance matrix is fixed, and the random vector depends on model parameters (the typical use case), it is done via a fallback path in LinearAlgebra: the constant matrix gets promoted to Dual type (i.e. it is copied each time the gradient is calculated), and then the generic triangular left division implementation is called. For 100x100 and larger matrices this results in big CPU and memory overhead.
However, when the matrix is constant, we don't need to convert it to Dual. We just have to left divide the dual values vector as well as each partial by this matrix. Since it would be the operation on fixed vectors, LAPACK's trtrs() could be used for much faster division. This is what this PR does. To avoid excessive copying, it relies on a hack: the array of duals could be accessed as a vector of N+1 floats.
Codecov Report
Patch coverage: 93.02% and project coverage change: -2.55 :warning:
Comparison is base (
e3670ce) 89.65% compared to head (1d38eac) 87.11%.
:exclamation: Current head 1d38eac differs from pull request most recent head 22d8670. Consider uploading reports for the commit 22d8670 to get more accurate results
Additional details and impacted files
@@ Coverage Diff @@
## master #589 +/- ##
==========================================
- Coverage 89.65% 87.11% -2.55%
==========================================
Files 11 9 -2
Lines 967 947 -20
==========================================
- Hits 867 825 -42
- Misses 100 122 +22
| Impacted Files | Coverage Δ | |
|---|---|---|
| src/dual.jl | 81.81% <93.02%> (-0.34%) |
:arrow_down: |
:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Do you have feedback about the report comment? Let us know in this issue.
Now these functions should be fixed and properly tested. The latter would require merging JuliaDiff/DiffTests.jl#11 first and updating the deps here.
@mcabbott @fredrikekre @devmotion I've refactored the code so that in-place ldiv/mul are also supported. These changes should now be covered by tests with the updated DiffTests package. There are some linalg differences between 1.0 and the current 1.x, so some of the tests are disabled on 1.0, I guess it's the most straightforward way to handle incompatibilities.
I've rebased the PR and cleaned up commit history a bit. For tests to succeed, DiffTests 0.1.3 is required (JuliaDiff/DiffTests.jl#17).
There were concerns regarding how efficient is constant_matrix * dual_matrix case, and whether it could be reinterpreted as normal numeric matrix multiplication, but I think plain reinterpretation is not possible.