# How can I use Optim/BFGS on functions I can't evaluate but whose gradients I can compute?

**URL:** <https://discourse.julialang.org/t/how-can-i-use-optim-bfgs-on-functions-i-cant-evaluate-but-whose-gradients-i-can-compute/35160>\
**Category:** Optimization (Mathematical)\
**Created:** [February 26, 2020, 10:04am UTC](https://discourse.julialang.org/t/how-can-i-use-optim-bfgs-on-functions-i-cant-evaluate-but-whose-gradients-i-can-compute/35160 "2020-02-26T10:04:05Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![marius311](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marius311/32/3953_2.png) [@marius311](https://discourse.julialang.org/u/marius311)\
**Post date:** [February 26, 2020, 10:04am UTC](https://discourse.julialang.org/t/how-can-i-use-optim-bfgs-on-functions-i-cant-evaluate-but-whose-gradients-i-can-compute/35160/1 "2020-02-26T10:04:05Z")

</div>

Suppose I have function I want to maximize with Optim’s BFGS algorithm. While I can compute gradients for this function, I cannot compute its value (this is typical if e.g. a high-dimensional determinant is involved, since gradients of a determinant reduce to a trace, but evaluating the determinant itself is N^3).

BFGS with a fixed step-size (ie no line-search) should be fine since it never needs the value, only the gradients. I can set a fixed step-size with `BFGS(; alphaguess = LineSearches.InitialStatic(alpha=0.1), linesearch = LineSearches.Static()`, but how do I tell Optim to not to try and use the target function?

I’ve tried just having it return 0 but that seems to terminate the optimization early. Any suggestions for how I can do this? Thanks.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 26, 2020, 12:09pm UTC](https://discourse.julialang.org/t/how-can-i-use-optim-bfgs-on-functions-i-cant-evaluate-but-whose-gradients-i-can-compute/35160/2 "2020-02-26T12:09:16Z")

</div>

Have you considered a root-finding method to solve ∇f=0? e.g. the analogue of BFGS for root-finding is Broyden’s method.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [February 26, 2020, 1:40pm UTC](https://discourse.julialang.org/t/how-can-i-use-optim-bfgs-on-functions-i-cant-evaluate-but-whose-gradients-i-can-compute/35160/3 "2020-02-26T13:40:21Z")

</div>

> I’ve tried just having it return 0 but that seems to terminate the optimization early. Any suggestions for how I can do this? Thanks.

Have you tried just setting a negative tolerance for the objective function?

---

<div class="post-metadata">

**Author:** ![lostella](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lostella/32/356_2.png) [@lostella](https://discourse.julialang.org/u/lostella)\
**Post date:** [February 26, 2020, 2:16pm UTC](https://discourse.julialang.org/t/how-can-i-use-optim-bfgs-on-functions-i-cant-evaluate-but-whose-gradients-i-can-compute/35160/4 "2020-02-26T14:16:19Z")

</div>

> [@stevengj](#):
>
> Have you considered a root-finding method to solve ∇f=0? e.g. the analogue of BFGS for root-finding is Broyden’s method.

I’ve implemented Broyden’s tooling in the form of an “updateable” linear operator [as part of ProximalAlgorithms](https://github.com/kul-forbes/ProximalAlgorithms.jl/blob/master/src/accel/broyden.jl): it’s actually the modified Broyden method by Powell (1970), which ensures that the operator is always invertibile. You can check [this test](https://github.com/kul-forbes/ProximalAlgorithms.jl/blob/master/test/accel/broyden.jl) for an example of how to use it, it’s very easy especially if no line search is needed.

---

<div class="post-metadata">

**Author:** ![lostella](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lostella/32/356_2.png) [@lostella](https://discourse.julialang.org/u/lostella)\
**Post date:** [February 26, 2020, 2:19pm UTC](https://discourse.julialang.org/t/how-can-i-use-optim-bfgs-on-functions-i-cant-evaluate-but-whose-gradients-i-can-compute/35160/5 "2020-02-26T14:19:09Z")

</div>

However, I guess if the gradient is C^1 then its Jacobian is symmetric and BFGS could work better

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [February 26, 2020, 2:33pm UTC](https://discourse.julialang.org/t/how-can-i-use-optim-bfgs-on-functions-i-cant-evaluate-but-whose-gradients-i-can-compute/35160/6 "2020-02-26T14:33:43Z")

</div>

> [@marius311](#):
>
> I cannot compute its value (this is typical if e.g. a high-dimensional determinant is involved, since gradients of a determinant reduce to a trace, but evaluating the determinant itself is N^3).

Also just out of curiosity, can you expand on this? To compute gradients of determinants, you need adjugates, which presumably are also N^3 ?

---

<div class="post-metadata">

**Author:** ![marius311](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marius311/32/3953_2.png) [@marius311](https://discourse.julialang.org/u/marius311)\
**Post date:** [February 27, 2020, 1:09am UTC](https://discourse.julialang.org/t/how-can-i-use-optim-bfgs-on-functions-i-cant-evaluate-but-whose-gradients-i-can-compute/35160/7 "2020-02-27T01:09:14Z")

</div>

> [@stevengj](#):
>
> Have you considered a root-finding method to solve ∇f=0? e.g. the analogue of BFGS for root-finding is Broyden’s method.

That’s a good point, thanks, I suppose the problem I’m really solving is ∇f=0. I wasn’t as familiar with Broyden’s but looking over it, it looks like the updates to the Hessian (at least as described on Wikipedia) are not exactly the same as what you’d get with BFGS on the equivalent problem. I’ll have to think about the different options (and thanks for the link @lostella). Fwiw, my problem is indeed C^1.

> [@antoine-levitt](#):
>
> Have you tried just setting a negative tolerance for the objective function?

I did but still terminated early. The only solution I’ve found so far is to artifically have the target function return a continually descreasing sequence, so it never triggers the termination criteria. Something stupid like:

```julia
f₀ = 0
optimize(
    x->(global f₀-=1), 
    x->my_actual_gradient(x),
    x₀, BFGS(...), inplace=false
)

```

> [@antoine-levitt](#):
>
> Also just out of curiosity, can you expand on this? To compute gradients of determinants, you need adjugates, which presumably are also N^3 ?

Yea sorry I needed to say more to make that make sense. You have that,

\frac{d}{dx} {\rm logdet}\, A(x) = {\rm tr}(A^{-1} \frac{dA}{dx})

In my case A is a positive definite N \times N implicit linear operator which can be applied in \mathcal{O}(N \log N), and similarly for the application of its gradient. You can do the trace with N\_{\rm MC} Monte-Carlo simulations using Hutchinson’s method, and you can do the inversion with N\_{\rm CG} conjugate gradient iterations, so the whole thing is \mathcal{O}(N\_{\rm MC} N\_{\rm CG} N \log N). Given my N\sim10^7, it works out that you can do the logdet gradient this way but the \mathcal{O}(N^3) determinant is hopeless (I’ve not found any stochastic determinant estimates that work here, although I’d be interested in hearing about suggestions!) I think this type of thing can arise often when you’re maximizing a Bayesian posterior and A is some high-dimensional covariance operator.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 27, 2020, 2:37am UTC](https://discourse.julialang.org/t/how-can-i-use-optim-bfgs-on-functions-i-cant-evaluate-but-whose-gradients-i-can-compute/35160/8 "2020-02-27T02:37:05Z")

</div>

You could use the stochastic Lanczos quadrature (SLQ) method to estimate the log determinant (= trace of log A).

Another nice thing about SLQ is that you can use the same Krylov bases to simultaneously compute your gradient trace, instead of using CG.

If you are using low-accuracy stochastic methods to estimate the function and gradient, you might want to use a stochastic optimization method like Adam. BFGS is not very robust to errors.

---

<div class="post-metadata">

**Author:** ![antoine-levitt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/antoine-levitt/32/4008_2.png) [@antoine-levitt](https://discourse.julialang.org/u/antoine-levitt)\
**Post date:** [February 27, 2020, 6:56am UTC](https://discourse.julialang.org/t/how-can-i-use-optim-bfgs-on-functions-i-cant-evaluate-but-whose-gradients-i-can-compute/35160/9 "2020-02-27T06:56:58Z")

</div>

For optim, try the allow\_f\_increases option maybe? If that doesn’t work, report it as a bug.

That’s a very fun problem you have. So you compute a sum of \<u, inv(A) dA/dx u\> = \<inv(A) u, dA/dx u\> with CG, which is independent of the size of x. Clever! Not sure if you can play the same trick with SLQ?

Conceptually you can swap out CG with the semiiterative Chebyshev method (which should even not be too inefficient if you’re in very large dimension anyway), which turns that into \<p(A) u, dA/dx u\> where p is a fixed polynomial. But even that has no nice antiderivative because the matrices don’t commute. Fun!

---

<div class="post-metadata">

**Author:** ![pkofod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pkofod/32/2179_2.png) [@pkofod](https://discourse.julialang.org/u/pkofod)\
**Post date:** [February 27, 2020, 12:29pm UTC](https://discourse.julialang.org/t/how-can-i-use-optim-bfgs-on-functions-i-cant-evaluate-but-whose-gradients-i-can-compute/35160/10 "2020-02-27T12:29:46Z")

</div>

Static is sort of an after thought, so we don’t handle this. I do get your point though… Can you post the actual error? That’ll make it easier to find. We will have to check for that line search in the lines that update the function and gradient.

---

<div class="post-metadata">

**Author:** ![marius311](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marius311/32/3953_2.png) [@marius311](https://discourse.julialang.org/u/marius311)\
**Post date:** [February 27, 2020, 11:13pm UTC](https://discourse.julialang.org/t/how-can-i-use-optim-bfgs-on-functions-i-cant-evaluate-but-whose-gradients-i-can-compute/35160/11 "2020-02-27T23:13:54Z")

</div>

Here’s what I checked, minimizing |x|^2 with a fixed step-size of 0.1. There’s no error, but if I just have my target function return e.g. 0, you can see it terminates early even with a negative `f_tol`,

```julia
using Optim, LineSearches, LinearAlgebra

bfgs = BFGS(
    alphaguess = LineSearches.InitialStatic(alpha=0.1), 
    linesearch = LineSearches.Static()
)

# do a proper minimization
optimize(
    x->dot(x,x),
    x->2x,
    ones(10),
    bfgs,
    inplace=false
).iterations # returns 181

# attempt to always return 0 from the target function,
# but with correct gradient
# this doesn't work and will terminate early
optimize(
    x->0,
    x->2x,
    ones(10),
    bfgs,
    Optim.Options(f_tol=-1, allow_f_increases=true),
    inplace=false
).iterations # returns 2

# the only workaround I've found 
f₀ = 0
optimize(
    x->(global f₀-=1),
    x->2x,
    ones(10),
    bfgs,
    inplace=false
).iterations # returns 181

```

---

<div class="post-metadata">

**Author:** ![pkofod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pkofod/32/2179_2.png) [@pkofod](https://discourse.julialang.org/u/pkofod)\
**Post date:** [February 28, 2020, 11:24am UTC](https://discourse.julialang.org/t/how-can-i-use-optim-bfgs-on-functions-i-cant-evaluate-but-whose-gradients-i-can-compute/35160/12 "2020-02-28T11:24:53Z")

</div>

Okay, try `Options(f_abstol=-1, f_reltol=-1)` perhaps.