Hi all, I’m new to Julia and have been using it for about a week, so apologies if the answer to this is obvious and I just haven’t been able to work it out.
The project that I’m using Julia for requires taking the logarithm of a polynomial and then optimizing it (e.g.: log(0.75*y + 125)). Initially I was doing this with SymPy, but I want to try and use some Julia packages to see how they compare. I was intending to use the DynamicPolynomials package like this:
using DynamicPolynomials
@polyvar x
f(x) = log(x)
f(x+2)
But I get an error:
ERROR: LoadError: MethodError: no method matching log(::Polynomial{true,Int64})
Closest candidates are:
log(::Float16) at math.jl:1018
log(::Complex{Float16}) at math.jl:1019
log(::Float64) at special/log.jl:254
Which I understand. Is there any sort of workaround to this? Maybe a symbolic logarithm function, or using another package would be better? I searched for a solution but I wasn’t able to find a lot of information. Thank you!
It’s a computational finance problem. I’m trying to perform utility optimization on an n-period model with one stock, where the investor buys and holds the stock. The utility function, U(x) I’m using is ln(x), where x is the capital at time n. So the process is to derive U(x), set it equal to zero, and solve for the number of shares I need to buy at time zero to maximize my expected utility at time n.
Are you just using gradient-based optimization, then? Or numerical root finding? I’m just trying to figure out first and foremost whether the optimization method you’re using inherently requires a symbolic representation of the function, or whether it’s also fine to just pass it a black-box function that given x returns the log of the polynomial evaluated at x (and possibly the gradient w.r.t. x if you’re using a gradient-based solver).
I haven’t actually gotten to the point of finding a derivative because of the error, but I was just going to use the Calculus package’s differentiate() function, set it equal to zero, and use a symbolic solver for y. I looked into gradient-based optimization, but I’m not sure if it’s the best method to use.
If I understand your question correctly, yes, it does require a symbolic representation.
You could just use the fact that (\log p)\prime = \frac{p\prime}{p}, which means that you just want p\prime to be zero. This allows you to stay in polynomial land:
using DynamicPolynomials
@polyvar y
p = 0.75 * y + 125 + y^2 + 0.1 * y^3
# let f(x) = log(x), g(y) = f(p(y)), then dg/dy = dp/dy / p;
# numerator of dg/dy:
dgdy_numerator = differentiate(p, y)
using PolynomialRoots
coeffs = [coefficient(dgdy_numerator, y^i) for i = 0 : maxdegree(dgdy_numerator)]
roots(coeffs)
Thank you so much! That works in this case, but eventually I would want to change the utility function to something like
U1(x) = -a*e^(-x*1/a)
U2(x) = 1/(-x^3)
I imagine that it would work with U2(x), but not so much with U1(x). I also couldn’t use your trick, since the derivative of an exponential is an exponential. Do you know of anything else I could try?
Well, one way to do generic U’s is to rely on automatic differentiation and generic root finding algorithms:
p(y) = 0.75 * y + 125 + y^2 + 0.1 * y^3
const a = 1000.
U1(x) = -a * exp(-x *1/a)
using ForwardDiff
dU1dy(y) = ForwardDiff.derivative(z -> U1(p(z)), y)
@show dU1dy(1)
using IntervalArithmetic, IntervalRootFinding
@show roots(dU1dy, -10..10) # *guaranteed* to find all roots in the interval [-10, 10]
Automatic differentiation (here, in the form of ForwardDiff.jl) is pretty cool if you haven’t seen it before.
Caveats:
This won’t work for U2(p(y)), which produces near-zero results in the interval [-10, 10].
IntervalRootFinding will not scale very well with ‘function complexity’; you could use e.g. Roots.jl for methods that scale better but aren’t guaranteed to find all the results in a given interval.
Thanks, you’ve been so helpful! So one issue that comes with the first method that you proposed is that my final rational polynomial ends up being the sum of a bunch of rational polynomials, so to begin with I’m not sure what the numerator and the denominator of the final rational function is. This is what my code looks like for generating the final function:
eval = 0
num = 0
for i in 0:n
tcap = (x0 - y)*(1+r)^n + u^i * d^(n-i)*y
num = DynamicPolynomials.differentiate(tcap, y)
eval += (prob)^i * (1-prob)^(n - i) * binomial(BigInt(n), BigInt(i)) * (num/tcap)
end
By putting in print statements, one function that I have is
(-0.5625*y) / (-0.5625*y^2 + 15625.0)
So now I want to find the zeros of the above expression, which happen only when the numerator is zero. Is there any way to isolate the numerator from the rational function without knowing what it is? Thank you so much for all your help!
By the way, in general you can rely on the chain rule, so you can start by finding the roots of the derivative of the polynomial; those will always be roots of the chained function as well. And then you can find the roots of the derivative of the utility function (e.g. using ForwardDiff and IntervalRootFinding) and then find where the polynomial achieves those values (more polynomial root finding problems).