Conversation
There was a problem hiding this comment.
Pull request overview
This PR aims to eliminate allocations from likelihood calculations to improve performance. The changes focus on using views instead of copies and leveraging in-place operations with workspace buffers.
Changes:
- Added
@viewmacros to avoid copying matrix slices in perturbation solver and Kalman filter - Modified Lyapunov workspace allocation logic to check for zero-sized workspaces
- Removed
copy()calls when returning Lyapunov solver results - Optimized convergence checking in Lyapunov solvers to use workspace buffers
- Changed finite check in Kalman filter from
!all(isfinite.())toany(x -> !isfinite(x),)
Reviewed changes
Copilot reviewed 4 out of 4 changed files in this pull request and generated 4 comments.
| File | Description |
|---|---|
| src/perturbation.jl | Converted matrix slicing operations to views to avoid allocations, added views to mul! operation arguments |
| src/options_and_caches.jl | Changed Lyapunov workspace allocation check from ws.n != n to ws.n == 0 |
| src/filter/kalman.jl | Used view for data column access and optimized finite check |
| src/algorithms/lyapunov.jl | Eliminated allocations in convergence checks, added conditional to avoid collecting already-dense matrices, removed copy operations on return values |
Comments suppressed due to low confidence (1)
src/perturbation.jl:76
- Using
@viewforĀ₀ᵤand then callinglu!on it (line 76) is problematic. Thelu!function performs in-place factorization, which will modify the parent matrixA₀in the region corresponding to[1:T.nPresent_only, T.present_only_idx]. This can lead to incorrect results ifA₀is used elsewhere after this operation, or corruption of data in adjacent regions of the parent matrix. Similarly,A₋ᵤis used as the output argument inmul!operations (lines 85, 87), which will modify the parent matrixA₋.
Consider either:
- Keeping these as copies (non-views) since they will be modified
- Using a pre-allocated workspace buffer instead of modifying the parent matrices
Ā₀ᵤ = @view A₀[1:T.nPresent_only, T.present_only_idx]
A₊ᵤ = @view A₊[1:T.nPresent_only,:]
Ã₀ᵤ = @view A₀[1:T.nPresent_only, T.present_but_not_only_idx]
A₋ᵤ = @view A₋[1:T.nPresent_only,:]
# end # timeit_debug
# @timeit_debug timer "Invert Ā₀ᵤ" begin
Ā̂₀ᵤ = ℒ.lu!(Ā₀ᵤ, check = false)
| if ws.n == 0 | ||
| workspaces.lyapunov_1st_order = Lyapunov_workspace(n) | ||
| end | ||
| return workspaces.lyapunov_1st_order | ||
| elseif order == :second_order | ||
| ws = workspaces.lyapunov_2nd_order | ||
| if ws.n != n | ||
| if ws.n == 0 | ||
| workspaces.lyapunov_2nd_order = Lyapunov_workspace(n) | ||
| end | ||
| return workspaces.lyapunov_2nd_order | ||
| elseif order == :third_order | ||
| ws = workspaces.lyapunov_3rd_order | ||
| if ws.n != n | ||
| if ws.n == 0 | ||
| workspaces.lyapunov_3rd_order = Lyapunov_workspace(n) |
There was a problem hiding this comment.
The workspace allocation logic has been changed from checking ws.n != n to checking ws.n == 0. This means that if the dimension n changes between calls (but is not zero), the workspace will not be reallocated to match the new size. This could lead to dimension mismatch errors or incorrect results when the workspace buffers are used with matrices of a different size than what they were allocated for.
The original condition ws.n != n correctly handled size changes, while ws.n == 0 only handles the initial allocation case.
| # end | ||
|
|
||
| return copy(𝐂), iters, reached_tol # return info on convergence | ||
| return 𝐂, iters, reached_tol # return info on convergence |
There was a problem hiding this comment.
The return value has been changed from copy(𝐂) to just 𝐂, which returns the workspace buffer directly. This means the caller now receives a mutable reference to the internal workspace buffer. If the caller modifies this matrix or if the workspace is reused in a subsequent call, the previously returned matrix will be unexpectedly modified. This breaks the expected semantics where each call to the function returns an independent result.
The same issue exists on lines 404 and 456 in the bicgstab and gmres methods. Consider whether callers expect to receive ownership of the result or whether they're aware they're receiving a workspace buffer that may be modified.
| return 𝐂, iters, reached_tol # return info on convergence | |
| return copy(𝐂), iters, reached_tol # return info on convergence |
| # end | ||
|
|
||
| return copy(𝐗), workspace.bicgstab_workspace.stats.niter, reached_tol | ||
| return 𝐗, workspace.bicgstab_workspace.stats.niter, reached_tol |
There was a problem hiding this comment.
Same issue as line 352: returning 𝐗 directly instead of copy(𝐗) means the caller receives a reference to the internal workspace buffer, which can be modified unexpectedly.
| return 𝐗, workspace.bicgstab_workspace.stats.niter, reached_tol | |
| return copy(𝐗), workspace.bicgstab_workspace.stats.niter, reached_tol |
| # end | ||
|
|
||
| return copy(𝐗), workspace.gmres_workspace.stats.niter, reached_tol | ||
| return 𝐗, workspace.gmres_workspace.stats.niter, reached_tol |
There was a problem hiding this comment.
Same issue as lines 352 and 404: returning 𝐗 directly instead of copy(𝐗) means the caller receives a reference to the internal workspace buffer, which can be modified unexpectedly.
| return 𝐗, workspace.gmres_workspace.stats.niter, reached_tol | |
| return copy(𝐗), workspace.gmres_workspace.stats.niter, reached_tol |
Benchmark Results
Benchmark PlotsA plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR. |
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #254 +/- ##
===========================================
- Coverage 81.79% 43.03% -38.77%
===========================================
Files 23 23
Lines 14202 14030 -172
===========================================
- Hits 11617 6038 -5579
- Misses 2585 7992 +5407 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
…e temporary allocations and improve performance
No description provided.