# Implementing the log-sum-exp function
**URL:** https://discourse.julialang.org/t/implementing-the-log-sum-exp-function/122946
**Category:** Optimization (Mathematical)
**Tags:** jump
**Created:** [November 22, 2024, 10:31am UTC](https://discourse.julialang.org/t/implementing-the-log-sum-exp-function/122946 "2024-11-22T10:31:16Z")
**Posts on this page:** 11
**Page:** 1
### Author:  [@kangqiu](https://discourse.julialang.org/u/kangqiu)
#### Post date: [November 22, 2024, 10:31am UTC](https://discourse.julialang.org/t/implementing-the-log-sum-exp-function/122946/1 "2024-11-22T10:31:16Z")
Hi everyone!
I am working on stochastic model predictive control in Julia/JuMP. I am using essentially a log-sum-exp function in the objective and encountered some issues with overflow. Now I am trying to implement a reformulation of the log-sum-exp function [1]:
I know JuMP is able to handle the max operation on x. But how is the sum without index k = argmax(x) implemented?
Any general pointers are greatly appreciated as well 🙂
[1] Pierre Blanchard, Desmond J Higham, and Nicholas J Higham. Accurately computing the log-sum-exp and softmax functions. IMA Journal of Numerical Analysis, 41(4):2311–2330, 08 2020
---
### Author:  [@tim.holy](https://discourse.julialang.org/u/tim.holy)
#### Post date: [November 22, 2024, 10:35am UTC](https://discourse.julialang.org/t/implementing-the-log-sum-exp-function/122946/2 "2024-11-22T10:35:26Z")
It’s already implemented here: [GitHub - JuliaStats/LogExpFunctions.jl: Julia package for various special functions based on `log` and `exp`.](https://github.com/JuliaStats/LogExpFunctions.jl)
I assume you know this, but your expression is analytically equivalent to \log\left(\sum\_i e^{x\_i}\right), i.e., just factoring out the largest term. You could also write it as y = x\_\mathrm{max} + \log\left(\sum\_i e^{x\_i-x\_\mathrm{max}}\right) and not worry about any index shenanigans.
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [November 22, 2024, 10:38am UTC](https://discourse.julialang.org/t/implementing-the-log-sum-exp-function/122946/3 "2024-11-22T10:38:08Z")
> [@kangqiu](#):
>
> Any general pointers are greatly appreciated as well 🙂
You might be interested in Convex.jl, whose disciplined convex programming reformulation can [handle `logsumexp`](https://jump.dev/Convex.jl/stable/manual/operations/#Exponential-Cone-Representable-Functions) provided you use an adequate solver.
---
### Author:  [@kangqiu](https://discourse.julialang.org/u/kangqiu)
#### Post date: [November 22, 2024, 10:48am UTC](https://discourse.julialang.org/t/implementing-the-log-sum-exp-function/122946/4 "2024-11-22T10:48:35Z")
Thank you for the answer! I tried that package but it since I work in the JuMP framework, it does not accept JuMP variable types as arguments. Therefore I am going about this manually…
The analytical equivalent gives overflow issues and my solver (Ipopt) also struggles when disregarding the max index stuff in the sum.
---
### Author:  [@kangqiu](https://discourse.julialang.org/u/kangqiu)
#### Post date: [November 22, 2024, 10:53am UTC](https://discourse.julialang.org/t/implementing-the-log-sum-exp-function/122946/5 "2024-11-22T10:53:56Z")
I was actually eyeing Convex.jl for its explicit handling of the logsumexp.
However I work with an SMPC with nonlinear model constraints, which is a nonconvex problem. I still use a local interior point solver (Ipopt) as we assume our initial state x\_0 at time of resolving is close enough to the optimum. I guess that is just common practice for real time optimization.
Would Convex.jl be applicable in this case?
---
### Author:  [@barucden](https://discourse.julialang.org/u/barucden)
#### Post date: [November 22, 2024, 12:11pm UTC](https://discourse.julialang.org/t/implementing-the-log-sum-exp-function/122946/6 "2024-11-22T12:11:30Z")
I don’t understand why the equivalent form x\_{\text{max}} + \log(\sum\_i \exp(x\_i - x\_{\text{max}})) suggested by Tim is not suitable.
For example, the following fails:
```julia
using JuMP
using Ipopt
function logsumexp(x)
return log(sum(exp.(x)))
end
function main()
n = 10
model = Model(Ipopt.Optimizer)
@variable(model, x[1:n])
@constraint(model, x .≥ 0)
@constraint(model, x[1] ≤ x[2] - 1000)
@objective(model, Min, logsumexp(x))
optimize!(model)
xval = value.(x)
@info xval
end
```
But it works when I change `logsumexp` as follows:
```julia
function logsumexp(x)
xmax = maximum(x)
return xmax + log(sum(exp.(x .- xmax)))
end
```
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [November 22, 2024, 1:08pm UTC](https://discourse.julialang.org/t/implementing-the-log-sum-exp-function/122946/7 "2024-11-22T13:08:22Z")
> [@kangqiu](#):
>
> However I work with an SMPC with nonlinear model constraints, which is a nonconvex problem. […]
> Would [Convex.jl](https://juliahub.com/ui/Packages/General/Convex) be applicable in this case?
I don’t think so. Convex.jl must [certify that the problem is convex](https://jump.dev/Convex.jl/stable/introduction/faq/#How-does-Convex.jl-differ-from-JuMP?), which it may even fail to do for some convex problems.
> [@barucden](#):
>
> I don’t understand why the equivalent form suggested by Tim is not suitable.
Maybe because `logsumexp` is not one of the [supported nonlinear operators](https://jump.dev/JuMP.jl/stable/manual/nonlinear/#Supported-operators)? You may want to look into [user-defined operators](https://jump.dev/JuMP.jl/stable/manual/nonlinear/#jump_user_defined_operators).
---
### Author:  [@ericphanson](https://discourse.julialang.org/u/ericphanson)
#### Post date: [November 22, 2024, 3:28pm UTC](https://discourse.julialang.org/t/implementing-the-log-sum-exp-function/122946/8 "2024-11-22T15:28:41Z")
What Convex.jl does is reformulate `logsumexp` as a conic programming problem, where you add a constraint that a particular expression belongs to an exponential cone. This is in some ways quite different to a rewriting of the formula, because with this approach, the solver does not execute the function you wrote, but rather receives an algebraic formulation of the problem (not e.g. a function pointer). To do a conic reformulation of logsumexp in JuMP, see: [How to implment `logsumexp` function in JuMP? - #2 by odow](https://discourse.julialang.org/t/how-to-implment-logsumexp-function-in-jump/84376/2)
Your solver needs to be able to handle exponential cones for this approach to work (both for Convex.jl and for JuMP.jl). I’m not sure if many nonlinear solvers can handle conic constraints, and I don’t think Ipopt can. So this approach unfortunately probably is not suitable.
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [November 22, 2024, 3:35pm UTC](https://discourse.julialang.org/t/implementing-the-log-sum-exp-function/122946/9 "2024-11-22T15:35:39Z")
> [@barucden](#):
>
> ```julia
> function logsumexp(x)
> xmax = maximum(x)
> return xmax + log(sum(exp.(x .- xmax)))
> end
>
> ```
I don’t know Jump, but I would have presumed that you would want to avoid intermediate allocations for functions used repeatedly in an optimization scenario.
So, for example
```julia
function logsumexp(x)
xmax = maximum(x)
exp_ = v -> exp(v - xmax)
return xmax + log(sum(exp_, x))
end
```
or
```julia
function logsumexp(x)
xmax = maximum(x)
return xmax + log(sum(exp(v - xmax) for v in x))
end
```
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [November 22, 2024, 3:41pm UTC](https://discourse.julialang.org/t/implementing-the-log-sum-exp-function/122946/10 "2024-11-22T15:41:50Z")
> [@DNF](#):
>
> I don’t know Jump, but I would have presumed that you would want to avoid intermediate allocations for functions used repeatedly in an optimization scenario.
AFAICT this does not matter to JuMP, because it does not actually call the function during optimization. It constructs a kind of symbolic representation internally, which is also why only a subset of operators are supported.
---
### Author:  [@kangqiu](https://discourse.julialang.org/u/kangqiu)
#### Post date: [November 23, 2024, 9:23am UTC](https://discourse.julialang.org/t/implementing-the-log-sum-exp-function/122946/12 "2024-11-23T09:23:55Z")
Your example works, thank you for clarifying! 🙂 I think I’m getting an evaluation error from somewhere else since the solver still spits out warnings in my code