Taylor series native exponential implementation#243
Conversation
|
Your PR no longer requires formatting changes. Thank you for your contribution! |
lkdvos
force-pushed
the
jh/taylorexponentiate
branch
from
84e07c2 to
08d3a71
Compare
|
@Jutho I hope you don't mind but I took some time to rebase this PR now that @sanderdemeyer's PR is finished. Benchmark:
|
| n | LA (µs) | Taylor (µs) | Taylor no-balance (µs) | Tayl/LA | nobal/LA | relerr |
|---|---|---|---|---|---|---|
| 8 | 3.9 | 5.9 | 5.6 | 1.51 | 1.43 | 1.3e-15 |
| 16 | 14.6 | 13.1 | 11.8 | 0.90 | 0.81 | 5.4e-15 |
| 32 | 57.8 | 30.1 | 26.7 | 0.52 | 0.46 | 1.2e-14 |
| 64 | 263.1 | 136.1 | 134.6 | 0.52 | 0.51 | 2.1e-14 |
| 128 | 2108.7 | 1219.7 | 1188.3 | 0.58 | 0.56 | 3.5e-14 |
| 256 | 6263.0 | 4161.5 | 4050.0 | 0.66 | 0.65 | 8.4e-14 |
ComplexF64
| n | LA (µs) | Taylor (µs) | Taylor no-balance (µs) | Tayl/LA | nobal/LA | relerr |
|---|---|---|---|---|---|---|
| 8 | 8.1 | 10.0 | 9.0 | 1.23 | 1.11 | 2.2e-15 |
| 16 | 31.1 | 36.8 | 30.0 | 1.19 | 0.97 | 1.9e-15 |
| 32 | 152.1 | 133.3 | 119.1 | 0.88 | 0.78 | 5.1e-15 |
| 64 | 943.1 | 801.2 | 719.9 | 0.85 | 0.76 | 2.5e-14 |
| 128 | 3456.6 | 2833.4 | 2591.5 | 0.82 | 0.75 | 3.0e-14 |
| 256 | 13064.8 | 11885.0 | 10696.6 | 0.91 | 0.82 | 7.4e-14 |
Taylor is faster than LinearAlgebra.exp for n ≥ 16 (Float64) and competitive-to-faster for ComplexF64, only paying a penalty at the smallest sizes. Balancing costs ~10–15% on these well-scaled random matrices; it is kept on by default for robustness on badly-scaled inputs.
Arbitrary precision (256-bit)
LinearAlgebra.exp has no method for Matrix{BigFloat}.
Reported below are Taylor timings and the self-consistency residual ‖exp(A)·exp(−A) − I‖₁.
BigFloat
| n | Taylor (ms) | no-balance (ms) | ‖E·E(−A) − I‖₁ |
|---|---|---|---|
| 8 | 1.46 | 1.38 | 1.9e-74 |
| 16 | 9.67 | 9.62 | 2.7e-73 |
| 32 | 75.63 | 76.57 | 3.0e-71 |
| 64 | 769.48 | 769.52 | 1.5e-68 |
Complex{BigFloat}
| n | Taylor (ms) | no-balance (ms) | ‖E·E(−A) − I‖₁ |
|---|---|---|---|
| 8 | 5.95 | 6.00 | 4.5e-74 |
| 16 | 44.53 | 43.26 | 2.9e-72 |
| 32 | 410.79 | 418.93 | 6.3e-71 |
| 64 | 3449.70 | 3432.66 | 2.0e-69 |
lkdvos
changed the title
exponential implementation
|
I guess for the |
Codecov Report❌ Patch coverage is
🚀 New features to boost your workflow:
|
Register a pure-Julia scaling-and-squaring Taylor algorithm as `MatrixFunctionViaTaylor` and fit it to the exponential interface: the `(τ, A)` forms reuse the existing generic method, and it becomes the default algorithm for all dense matrices (`DiagonalAlgorithm` still handles `Diagonal`). Being LAPACK-free, it also covers arbitrary precision, where `LinearAlgebra.exp` has no method. The truncation order and number of squarings are chosen to minimize the Paterson-Stockmeyer matrix-multiplication count under the Taylor remainder bound, the polynomial is evaluated with incrementally built coefficients, and squaring reuses a ping-pong buffer. `balance!` is rewritten as a clean, type-generic Parlett-Reinsch scaling and wired in as an optional preprocessing step. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
lkdvos
force-pushed
the
jh/taylorexponentiate
branch
from
a7f2e42 to
0d07fe9
Compare
|
It seems like for BigFloat the decompositions are actually faster in some cases and slower in other, in particular for hermitian inputs. I'm not sure what should be the conclusion of what algorithm to choose by default, or what is causing this since that is slightly unexpected to me, but I guess the matmuls just aren't as efficient, so the more complex work needed for eigensolvers is less bad. |
lkdvos
force-pushed
the
jh/taylorexponentiate
branch
from
357d1e6 to
b111776
Compare
2 participants
This is something I started a while ago, but I need to continue working on it. Basically, for higher precision number types, the best approach for computing the exponential is just a Taylor expansion, combined with scaling and squaring (i.e. exp(A) = (exp(A/2^s))^(2^s)). But which
sin combination with which order ofmof the Taylor expansion is a nontrivial question, that has been studied in the literature, and also involves clever ways of computing higher order polynomials, i.e. polynomials ofAof degreem, thus requiring all powersA^k,k in 0:m, without actually having to dommatrix multiplications.