-
-
Notifications
You must be signed in to change notification settings - Fork 232
Commit 156a23b
committed
Refactor use of the Bareiss algorithm
This PR goes with JuliaLang/julia#42683.
When I was originally reading this code, I was having trouble
figuring out exactly what the `bareiss!` function was actually
doing. In order to hopfully make that easier to follow, this
attempts to separate the original bareiss! function into three
separate concerns:
1. The core bareiss algorithm (in LinearAlgebra)
2. The compressed eadj/cadj sparse matrix representation of the
adjancency matrix
3. MTK's pivoting heuristics
At first glance not much is gained by this refactoring, but there
are notable advantages:
1. The new sparse matrix datastructure is interconvertable with regular
sparse or dense matrices as well as displaying as a matrix.
This makes it much easier to swap in and out matrix types to
validate correctness and understand the algorithm.
2. It's much easier to gain confidence in optimizing the data structure
further. For example, currently, the inner loop has O(N^4) behavior,
rather than the expected O(N^3). This can be fixed by making the
data structure maintain coefficients in sorted order, but it's hard
to see that if everything is interleaved.
3. It's compatible with future extensions to take advantage of hierarchical
system structure.
I'm planning a number of further refactorings here, but this is the
first self-contained piece. I've kept most behavioral changes and
optimizations out of this PR, except for two:
1. I adjusted the early-out at https://github.com/SciML/ModelingToolkit.jl/blob/9bf38dd2b8d8c2fb0412d223ca95e713c666d81f/src/systems/alias_elimination.jl#L280
The base implementation of bareiss doesn't have it. Instead, I put in
a more limited early-out path that should have essentially the same
performance in MTK (where most pivots are `1` and `-1`), but isn't
mathematically objectionable.
2. I'm not doing a bariess-restart. As written, the bareiss restart
basically resets the prev_pivot to `1`. This isn't terrible for what
we need it for (it merely multiplies the rank minors by a scalar factor),
but it does technically break the coefficient-growth guarantees of the
bareiss algorithm. Instead, I simply merged the three pivot stages into
one function that simply records the boundries between the stages.
For the same reason that it shouldn't affect MTK's ultimate output to
do the restarts, removing it shouldn't either, but I felt it worth pointing out.1 parent e803bac commit 156a23bCopy full SHA for 156a23b
File tree
Expand file treeCollapse file tree
3 files changed
+449
-143
lines changedFilter options
- src/systems
- compat
Expand file treeCollapse file tree
3 files changed
+449
-143
lines changed
0 commit comments