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:
|
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
|
That is not cool. Is that really due to the slow matrix multiplication. Would that be overhead due to the allocation of the different bigfloats. Could you try the same with something like |
|
After looking at it a bit more, one thing I noticed is that the Taylor series is truncated such that the error is To be completely fair, I would assume there is some equivalent that is true for the eigen decomposition, but I know too little of that implementation to know how that scales. I reran the benchmarks with DoubleFloat, the story is slightly different but it looks like it is hard to beat the hermitian diagonalization... |
|
I investigated a little further, and the main thing that I can come up with is to try and get a tighter bound on the operator norm estimator. Currently, we are just using There are some ways of improving on this, but things get complicated quite quickly, these methods are not free, and the actual payoff is not that large, so I think I might vote to just keep things the way they are right now. As this already is faster than LinearAlgebra, and we don't immediately have a need for hyper-efficient dense exponentials, I think this should be fine :). |
|
I don't recognize much of the original implementation anymore. I know that, since you anyway need to compute a minimum number of powers of A, it can be useful to first compute some powers of
This is based on Gelfand's theorem which states that |
|
I think I now have this implemented as well. The speedup is not gigantic, but it definitely is measurable so worth keeping. In any case, I think I'm done here now, I don't think it is worth it to try and optimize generic exponentiation much past this, since without more information about the actual input it is somewhat dangerous to not slow down the small |
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
441b4ff to
fa8a6e1
Compare
| # leaves degenerate rows/columns untouched). | ||
| @inline function _balance_factor(colnorm::R, rownorm::R, β::R, logβ::R) where {R} | ||
| (colnorm > 0 && rownorm > 0) || return one(R) | ||
| f = β^round(log(rownorm / colnorm) / (2 * logβ)) |
There was a problem hiding this comment.
What is the role of round here? It seems this yields scale factors that are powers of 2 (for the default choice of β). Is this important?
There was a problem hiding this comment.
Did some digging: this is also what LAPACK gebal does, the point is that the radix argument above is chosen so the rescaling does not lead to any floating point errors, making the similarity transform exact.
There was a problem hiding this comment.
Ok makes sense. Just wondering whether there is any point in using floating points for f in the first place then. I don't know if you can easily multiply floating point numbers by 2 as you could do with inter << 1.
There was a problem hiding this comment.
Anyway, I have to go now, and there are some other aspects of this PR that I would like to understand a bit better, so I'll try to get back to it tonight.
There was a problem hiding this comment.
I hardcoded the base to be 2 now - this is still assuming that multiplication/division by powers of 2 are exact in floating arithmetic, and I don't really foresee any element types where this is not the case for now.
I was experimenting slightly with Base's ldexp and frexp, which would be a way to use the integer powers instead, but that does not really gain any performance and hinders the generalizability since this function is for example not defined on dual numbers, which would now actually work.
| numblocks = fld(order, blocksize) | ||
| taylor_block!(result, powers, invfactorial, numblocks, blocksize, order) | ||
| for j in (numblocks - 1):-1:0 | ||
| mul!(block, result, powers[blocksize]) | ||
| taylor_block!(result, powers, invfactorial, j, blocksize, order) | ||
| result .+= block | ||
| end |
There was a problem hiding this comment.
I think the name numblocks is a bit confusing, since there is a first call to taylor_block! before the loop, and then in the loop j goes from 0 to numblocks-1 (in reverse order), so the total number of taylor_block! calls is numblocks + 1.
Indeed, my way of understanding the number of blocks is by reasoning that there are order + 1 terms in the expansion, and every block gives rise to at most blocksize terms, resulting in numblocks = cld(order + 1, blocksize). In the end this should be the same, but this gives rise to the following suggestion:
| numblocks = fld(order, blocksize) | |
| taylor_block!(result, powers, invfactorial, numblocks, blocksize, order) | |
| for j in (numblocks - 1):-1:0 | |
| mul!(block, result, powers[blocksize]) | |
| taylor_block!(result, powers, invfactorial, j, blocksize, order) | |
| result .+= block | |
| end | |
| numblocks = cld(order + 1, blocksize) | |
| taylor_block!(result, powers, invfactorial, numblocks - 1, blocksize, order) | |
| for blockindex in (numblocks - 1):-1:1 | |
| mul!(block, result, powers[blocksize]) | |
| taylor_block!(result, powers, invfactorial, blockindex - 1, blocksize, order) | |
| result .+= block | |
| end |
Rather than passing blockindex - 1, we could of course also modify taylor_block! to take blockindex as an argument, and then do base = (blockindex - 1) * blocksize (which I would probably call offset instead of base).
There was a problem hiding this comment.
In fact, cld(order + 1, blocksize) appears in the cost calculation on line 206 above. For numblocks = cld(order + 1, blocksize), there are indeed numblocks - 1 matrix multiplications.
There was a problem hiding this comment.
But that is referring to what I would call numblocks, which is one more than what is called numblocks in the current code.
| excess = log2(θ) - (log2tol + log2factorial) / (order + 1) | ||
| squarings = excess > 0 ? ceil(Int, excess) : 0 | ||
| blocksize = ceil(Int, sqrt(order + 1)) | ||
| cost = max(0, blocksize - p₀) + (cld(order + 1, blocksize) - 1) + squarings |
There was a problem hiding this comment.
I have two comments:
- I think the cost model is not exactly right, but that is also because of how
taylor_polynomialcurrently works. Currently, for a givenblocksize(let me call thisb), it builds the polynomial as
(A[0]:A[b-1] + A[b]*(A[0]:A[b-1] + A[b]*(A[0]:A[b-1] + ....)))
with A[i] indicating A^i, and A[i]:A[j] just the sum of terms of these powers. So indeed, as remarked in my earlier comment, the number of "blocks" is cld(order+1, blocksize) since there are order+1 terms, and every block contributes blocksize terms A[0]:A[b-1]. The optimal blocksize is then ceil(Int, sqrt(order+1)), e.g. for order == 15, we find blocksize=4, and then, assuming we already have the powers A[1] to A[4], we can build a order=15 polynomial with (cld(order + 1, blocksize) - 1) = 3 multiplications as
A[0]:A[3] + A[4]*(A[0]:A[3] + A[4]*(A[0]:A[3] + A[4]*(A[0]:A[3])))
However, you can easily see that for the same cost we can build an order 16 polynomial as well, namely by doing the blocking differently as
A[0] + A[1]:A[4] + A[4]*(A[1]:A[4] + A[4]*(A[1]:A[4] + A[4]*(A[1]:A[4])))
i.e. by separating the very first term in the Taylor decomposition. Alternatively, you can also modify the inner most block to do
A[0]:A[3] + A[4]*(A[0]:A[3] + A[4]*(A[0]:A[3] + A[4]*(A[0]:A[4])))
So actually I think the cost model should be:
blocksize = max(p₀, ceil(Int, sqrt(order))) # use that you anyway have `p₀` powers available for free
cost = max(0, blocksize - p₀) + (cld(order, blocksize) - 1) + squaring- The second remark is that, because the cost is not a continuous function of
order, it does actually make sense to not try all orders. If we only consider the first two terms in the cost function and compute them
orders = 1:100
costs = map(orders) do order
blocksize = max(p₀, ceil(Int, sqrt(order)))
return max(0, blocksize - p₀) + (cld(order, blocksize) - 1)
end
hcat(orders, costs)then you get the cost per order, not counting the squarings. But then you see that there are a whole bunch of orders with the same costs, and since squarings can only decrease for increasing order, the best option is always to take the highest order for a given cost. Indeed, with p₀ = 1 (assuming you have not yet computed any power, but of course you have A[1]), this code exactly reproduces the figure here:
coming from
https://www.sciencedirect.com/science/article/pii/S0024379519301454
By looping over all orders, you would of course find the optimal one, but with too much work. Of course, in order to avoid this, one needs to tabulate the optimal values of order, and they depend on p₀.
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.