# Manipulating vectors with logsumexp
**URL:** https://discourse.julialang.org/t/manipulating-vectors-with-logsumexp/83955
**Category:** General Usage
**Tags:** vector
**Created:** [July 8, 2022, 3:20pm UTC](https://discourse.julialang.org/t/manipulating-vectors-with-logsumexp/83955 "2022-07-08T15:20:54Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@samerb](https://discourse.julialang.org/u/samerb)
#### Post date: [July 8, 2022, 3:20pm UTC](https://discourse.julialang.org/t/manipulating-vectors-with-logsumexp/83955/1 "2022-07-08T15:20:54Z")
I have an n-dimensional vector \delta and I would like to compute the following vector expression.
\left[\log \left( \sum\_{i=1}^n \exp(\delta\_i) \right) - \log \left( \sum\_{i \neq j} \exp(\delta\_i) \right) \right]\_{j=1}^n
What’s the best way to do this in a compact way (hopefully no loops).
I have access to the function `logsumexp` from StatsFuns so ideally I could use that but not necessary. Thanks!
---
### Author:  [@skleinbo](https://discourse.julialang.org/u/skleinbo)
#### Post date: [July 8, 2022, 4:22pm UTC](https://discourse.julialang.org/t/manipulating-vectors-with-logsumexp/83955/2 "2022-07-08T16:22:25Z")
`logsumexp` takes an iterator as its argument, so this works
```julia
f(d) = logsumexp(d) .- ( logsumexp( (d[i] for i in eachindex(d) if i!=j) ) for j in eachindex(d) )
```
---
### Author:  [@mikmoore](https://discourse.julialang.org/u/mikmoore)
#### Post date: [July 8, 2022, 5:28pm UTC](https://discourse.julialang.org/t/manipulating-vectors-with-logsumexp/83955/3 "2022-07-08T17:28:43Z")
With some algebra, the expression you are trying to evaluate is equal to
x\_j = \log\left(1+\frac{\exp(\delta\_j-c)}{\sum\_{i\ne j} \exp(\delta\_i-c)}\right)
for any choice c. I’ll recommend the conventional `logsumexp` choice c=\max\_i \delta\_i. This version can take advantage of the special function `log1p` to be more accurate for small arguments.
Note that we can equivalently (although risking some numerical cancellation in the implementation) write \sum\_{i\ne j} \exp(\delta\_i-c) = -\exp(\delta\_j-c)+\sum\_{i=1}^n \exp(\delta\_i-c), where the right-side version allows the precomputation of the sum. Weigh that tradeoff of efficiency versus precision. Below, I have opted for the efficient version that risks that cancellation error in some edge cases:
```julia
expscale = let c=maximum(deltas); z->exp(z-c); end
sumexpdeltas = sum(expscale, deltas)
x = [log(sumexpdeltas) - log(sumexpdeltas - expscale(d)) for d in deltas] # initial version
x = [(ed = expscale(d); log1p(ed/(sumexpdeltas-ed))) for d in deltas] # with log1p trasform
```
EDIT: the cancellation error of the sum transformation I used is most-significant when the sum is dominated by a single term (one \delta\_i is much bigger than all others). If this is a possibility, you’ll want to use the leave-one-out version of the summation instead. The `log1p` version avoids a numerical error when \delta\_j is much smaller than \max\_i \delta\_i.
EDIT 2: For the purposes of my previous statement and `Float64` precision, a difference of 20 or 30 would be the threshold to consider a number “much bigger”. That isn’t necessarily a lot, so I suspect you’ll want to use both leave-one-out summation and the transformed `log1p` equation to be safe.
EDIT 3:
Just a teensy bit more algebra yields
x\_j = \log\left(1+\frac{1}{\sum\_{i\ne j} \exp(\delta\_i-\delta\_j)}\right)
Given my previous numerical concerns, this is the version I’ll recommend you implement.