fzero(f, j, guess) works for me when f and j are functions. f(Af, guess)
works for me now when Af is an @anon function.
On Tuesday, July 7, 2015 at 7:34:39 PM UTC-4, j verzani wrote:
>
> Okay, this just got fixed as much as I could with v"0.1.15" (there is no
> fzero(f,j,guess) signature).
>
> On Tuesday, July 7, 2015 at 4:38:41 PM UTC-4, Andrew wrote:
>>
>> Just checked. So, Roots.fzero(f, guess) does work. However,
>> Roots.fzero(f, j, guess) doesn't work, and neither does Roots.newton(f, j,
>> guess).
>>
>> I looked at the Roots.jl source and I see ::Function annotations on the
>> methods with the jacobian, but not the regular one.
>>
>> On Tuesday, July 7, 2015 at 4:22:17 PM UTC-4, j verzani wrote:
>>>
>>> It isn't your first choice, but `Roots.fzero` can have `@anon` functions
>>> passed to it, unless I forgot to tag a new version after making that change
>>> on master not so long ago.
>>>
>>> On Tuesday, July 7, 2015 at 2:29:51 PM UTC-4, Andrew wrote:
>>>>
>>>> I'm writing this in case other people are trying to do the same thing
>>>> I've done, and also to see if anyone has any suggestions.
>>>>
>>>> Recently I have been writing some code that requires solving lots(tens
>>>> of thousands) of simple non-linear equations. The application is
>>>> economics,
>>>> I am solving an intratemporal first order condition for optimal labor
>>>> supply given the state and a savings decision. This requires solving the
>>>> same equation many times, but with different parameters.
>>>>
>>>> As far as I know, the standard ways to do this are to either define a
>>>> nested function which by the lexical scoping rules inherits the parameters
>>>> of the outer function, or use an anonymous function. Both these methods
>>>> are
>>>> slow right now because Julia can't inline those functions. However, the
>>>> FastAnonymous package lets you define an anonymous "function", which
>>>> behaves exactly like a function but isn't type ::Function, which is fast.
>>>> Crucially for me, in Julia 0.4 you can modify the parameters of the
>>>> function you get out of FastAnonymous. I rewrote some code I had which
>>>> depended on solving a lot of non-linear equations, and it's now 3 times as
>>>> fast, running in 2s instead of 6s.
>>>>
>>>> Here I'll describe a simplified version of my setup and point out a few
>>>> issues.
>>>>
>>>> 1. I store the anonymous function in a type that I will pass along to
>>>> the function which needs to solve the nonlinear equation. I use a
>>>> parametric type here since the type of an anonymous function seems to vary
>>>> with every instance. For example,
>>>>
>>>> typeof(UF.fhoursFOC)
>>>> FastAnonymous.##Closure#11431{Ptr{Void}
>>>> @0x00007f2c2eb26e30,0x10e636ff02d85766,(:h,)}
>>>>
>>>>
>>>> To construct the type,
>>>>
>>>> immutable CRRA_labor{T1, T2} <: LaborChoice # <: means "subtype of"
>>>> sigmac::Float64
>>>> sigmal::Float64
>>>> psi::Float64
>>>> hoursmax::Float64
>>>> state::State # Encodes info on how to solve itself
>>>> fhoursFOC::T1
>>>> fJACOBhoursFOC::T2
>>>> end
>>>>
>>>> To set up the anonymous functions fhoursFOC and fJACOBhoursFOC (the
>>>> jacobian), I define a constructor
>>>>
>>>> function CRRA_labor(sigmac,sigmal,psi,hoursmax,state)
>>>> fhoursFOC = @anon h -> hoursFOC(CRRA_labor(sigmac,sigmal,psi,
>>>> hoursmax,state,0., 0.) , h, state)
>>>> fJACOBhoursFOC = @anon jh -> JACOBhoursFOC(CRRA_labor(sigmac,sigmal
>>>> ,psi,hoursmax,state,0., 0.) , jh, state)
>>>> CRRA_labor(sigmac,sigmal,psi,hoursmax,state,fhoursFOC,
>>>> fJACOBhoursFOC)
>>>> end
>>>>
>>>> This looks a bit complicated because the nonlinear equation I need to
>>>> solve, hoursFOC, relies on the type CRRA_labor, as well as some aggregate
>>>> and idiosyncratic state info, to set up the problem. To encode this
>>>> information, I define a dummy instance of CRRA_labor, where I supply 0's
>>>> in
>>>> place of the anonymous functions. I tried to make a self-referential type
>>>> here as described in the documentation, but I couldn't get it to work, so
>>>> I
>>>> went with the dummy instance instead.
>>>>
>>>> @anon sets up the anonymous function. This means that code like
>>>> fhoursFOC(0.5) will return a value.
>>>>
>>>> 2. Now that I have my anonymous function taking only 1 variable, I can
>>>> use the nonlinear equation solver. Unfortunately, the existing nonlinear
>>>> equation solvers like Roots.fzero and NLsolve ask the argument to be of
>>>> type ::Function. Since anonymous functions work like functions but are
>>>> actually some different type, they wouldn't accept my argument. Instead, I
>>>> wrote my own Newton method, which is like 5 lines of code, where I don't
>>>> restrict the argument type.
>>>>
>>>> I think it would be very straightforward to make this a multivariate
>>>> Newton method.
>>>>
>>>> function myNewton(f, j, x)
>>>> for n = 1:100
>>>> fx , jx = f(x), j(x)
>>>> abs(fx) < 1e-6 && return x
>>>> d = fx/jx
>>>> x = x - d
>>>> end
>>>> println("Too many iterations")
>>>> return NaN
>>>> end
>>>>
>>>> 3. The useful thing here in 0.4 is that you can edit the parameters of
>>>> the anonymous function. The parameters are encoded in a custom type
>>>> state::State, and I update the state. Then I call my nonlinear equation
>>>> solver
>>>>
>>>> UF.fhoursFOC.state, UF.fJACOBhoursFOC.state = state, state
>>>> f = UF.fhoursFOC
>>>> j = UF.fJACOBhoursFOC
>>>> hours = myNewton(f, j, hoursguess)
>>>>
>>>> This runs much faster than my old version which used NLsolve, which
>>>> itself ran faster than a version using Roots.fzero.
>>>>
>>>> Issues:
>>>>
>>>> 1. Since the type of the anonymous function isn't ::Function, I had to
>>>> write my own solver. I'm pretty sure a 1-line edit to Roots.fzero where I
>>>> just remove the ::Function type annotation would let it work there, but
>>>> I'm
>>>> not aware of another workaround.
>>>>
>>>> 2. I would rather use NLsolve, which uses in-place updating of its
>>>> arguments ( f!(input::Array, output::Array) ), but I've tried constructing
>>>> an anonymous function that does that, and @anon didn't work. Perhaps there
>>>> is a workaround.
>>>>
>>>> 3. Since I'm using an anonymous function, I have to explicitly pass it
>>>> around. Encoding it into the type CRRA_labor wasn't really hard though.
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
