I run a model (MTK) where one term of an equation does log(x^0.5). This seems to work sometimes, but also crashes for some parameter values when I do Monte Carlo simulation.
If I change the term to log(sqrt(x)), the simulation seems to always work.
Question: Is there a reason for the problem with x^0.5?
[I’d like to write the model as x^y where y happens to be y=0.5 in this particular case. There may be other formulations where y differs from 0.5.]
any reason you aren’t writing this as 0.5*log(x)? I’d expect that to be notably faster
I could, but the logarithm contains other multiplicative terms, too. But I could write it as a sum of logarithms.
in general, that will be a lot less likely to suffer from overflow and underflow (and depending on how many times you need to call log, it will likely be faster as well).
OK - I’ll try to do that. It still puzzles me that the simulation crashes when I do log(x^0.5), but not when I do log(sqrt(x)).
Will do it tomorrow… I left my job PC with the VSCode file open, and it messes up the file slightly if I open it form other computers while still open at work.
What do you mean by “crash”? Do you mean it throws an exception? What is the exact error message?
(MTK does symbolic math, so I’m guessing it has special simplification rules for sqrt that it may not apply to arbitrary powers.)
The solver seems to try to take the square root of a negative number… I do run a Monte Carlo simulation where I let some dynamic model parameters change \pm 10% from a uniform distribution. I’m pretty sure these variations are too small to make the system unstable. And the “instability”/square root of a negative number only occurs for some cases in the MC run, and only when I use log(x^0.5); not when I do log(sqrt(x)). [I had to arbitrarily truncate the error message below, because it was too long for this discussion system to handle the entire error message.]
{
"name": "DomainError",
"message": "DomainError with -327.6113408722682:\nExponentiation yielding a complex result requires a complex argument.\nReplace x^y with (x+0im)^y, Complex(x)^y, or similar.",
"stack": "DomainError with -327.6113408722682:\nExponentiation yielding a complex result requires a complex argument.\nReplace x^y with (x+0im)^y, Complex(x)^y, or similar.\n\nStacktrace:\n [1] throw_exp_domainerror(x::Float64)\n @ Base.Math .\\math.jl:41\n [2] ^\n @ .\\math.jl:1157 [inlined]\n [3] macro expansion\n @ C:\\Users\\...\\.julia\\packages\\SymbolicUtils\\jf8aQ\\src\\code.jl:387 [inlined]\n [4] macro expansion\n @ C:\\Users\\...\\.julia\\packages\\RuntimeGeneratedFunctions\\M9ZX8\\src\\RuntimeGeneratedFunctions.jl:163 [inlined]\n [5] macro expansion\n @ .\\none:0 [inlined]\n [6] generated_callfunc\n @ .\\none:0 [inlined]\n [7] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var\"#_RGF_ModTag\", ModelingToolkit.var\"#_RGF_ModTag\", (0x933cfd09, 0xe4c716be, 0xcd93182a, 0x48cad04e, 0xc7b8348d), Nothing})(::Vector{Float64}, ::Vector{Float64}, ::Vector{Float64}, ::Float64)\n @ RuntimeGeneratedFunctions C:\\Users\\...\\.julia\\packages\\RuntimeGeneratedFunctions\\M9ZX8\\src\\RuntimeGeneratedFunctions.jl:150\n [8] f\n @ C:\\Users\\...\\.julia\\packages\\ModelingToolkit\\duznJ\\src\\systems\\diffeqs\\abstractodesystem.jl:379 [inlined]\n [9] Void\n @
}
OK, so a slight summary. I have a system where a=0.5 and b=1 (but a,b could have different values). Then I originally wrote:
in a differential equation. This sometimes crashes (i.e., gets unstable??, and produces a negative x value). If I instead write:
then the code works.
Also, if I change the expression and write
the code works.
Regarding efficiency, will the first expression (\log(x^a\cdot y^b)) always be slower than the third expression (a\cdot \log(x) + b\cdot \log(y))? Maybe so, but I guess that depends on the relative time to take the logarithm vs. taking the square root/rising to power.
Generally speaking, yes. Summation and multiplication are heavily optimized for, since they are so fundamental, and the hardware itself (the computer chip) usually implements these instructions very efficiently. On the other hand, raising a number to an arbitrary power involves more computational steps. I don’t know how this is implemented in actual hardware, but even when you would do the calculations “by hand”, you’d probably have to come up with some generic algorithm that does the ^ operation as a combination of sums and multiplications, unless in some special cases, e.g. squaring a number.
But that would just be an issue about runtime and in principle the same holds true for log - so as you said, it might be that the two log operations plus one ^ take longer than one log and two ^ (here it doesn’t seem to be the case, see below). But the other issues are accuracy and stability.
As @Oscar_Smith already mentioned, if you plug in the “right” values for the base and exponent you can quickly run into situations where the result of x ^ a gets quite large\small. Of course, you take the log in the end, so the final result should be manageable again, but it might be that the intermediate value cannot be stored as a floating point number, or at least not accurately enough.
Another point regarding accuracy: Generally speaking, the more operations (adding, multiplying, …) you perform on a number, the more all the inherent tiny round-off errors
will compound. A single multiplication will have a certain maximum error, but if ^ consists of multiple operations, it might in principle have larger error (I don’t know if this is true for ^ specifically).
Even if you just sum a long list of numbers, the error can grow larger the longer the list gets. In that case there are specialized methods to account for such errors: Kahan summation algorithm - Wikipedia
When solving an ODE, you essentially also perform a lot of operations, where compounding errors could well lead to instabilities (the state at very late times could look quite different compared to early times even if there are only small errors along the way) – but what exactly happens in your case would also depend on other factors and the details of the ODE.
So yeah, I think it’s safe to assume that re-writing the log as a product and sum only has benefits. At least I can’t think of any downsides right now…
PS: Here is a quick benchmark comparing the two approaches:
julia> using BenchmarkTools
julia> f(x, y, a, b) = log( x ^ a * y ^ b )
f (generic function with 1 method)
julia> g(x, y, a, b) = a * log(x) + b * log(y)
g (generic function with 1 method)
julia> r = rand(4)
4-element Vector{Float64}:
0.30426897441393985
0.24398687655350038
0.6439939608519193
0.5045082195046893
julia> @btime f($r[1], $r[2], $r[3], $r[4])
56.441 ns (0 allocations: 0 bytes)
-1.477931723609483
julia> @btime g($r[1], $r[2], $r[3], $r[4])
13.024 ns (0 allocations: 0 bytes)
-1.4779317236094829
In general x^y is computed as exp(y*log(x)) (except you actually need to use a slower version of the logarithm that is extra accurate to prevent roundoff error. sqrt is fairly fast (about as fast as division so about 10x worse than multiplication/addition), but currently the compiler isn’t able to optimize x^0.5 into sqrt.
A post was split to a new topic: What’s the syntax to create a collapsable section on Discourse?