# Efficient ways to raise an array by a scalar power
**URL:** https://discourse.julialang.org/t/efficient-ways-to-raise-an-array-by-a-scalar-power/107688
**Category:** Performance
**Tags:** question
**Created:** [December 15, 2023, 10:03pm UTC](https://discourse.julialang.org/t/efficient-ways-to-raise-an-array-by-a-scalar-power/107688 "2023-12-15T22:03:57Z")
**Posts on this page:** 8
**Page:** 1
### Author:  [@valentin.m](https://discourse.julialang.org/u/valentin.m)
#### Post date: [December 15, 2023, 10:03pm UTC](https://discourse.julialang.org/t/efficient-ways-to-raise-an-array-by-a-scalar-power/107688/1 "2023-12-15T22:03:57Z")
Hello!
I’m currently working on optimizing an algorithm to match the performance of an existing Fortran implementation. A key part of the algorithm involves computing the power of large arrays of floats raised to a scalar float.
I’ve looked into the IntelVectorMath.jl and AppleAccelerate packages to optimize that. However, both packages only seem to support element-wise power operations between two vectors of floats. To use these, I would need to allocate a vector filled with the scalar value, matching the size of my base array, which seems inefficient.
My research suggests that both IntelVectorMath and AppleAccelerate could theoretically provide functions for array-to-scalar power operations, but these are not exposed in their Julia wrappers.
Given this, I am seeking advice on the most efficient way to compute the power of large float arrays by a scalar float in Julia. Are there alternative libraries or methods within Julia’s ecosystem that I should consider? Or is there a way to extend the existing wrappers to expose the functionality I need?
Any guidance or suggestions would be greatly appreciated.
Thank you!
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [December 15, 2023, 10:07pm UTC](https://discourse.julialang.org/t/efficient-ways-to-raise-an-array-by-a-scalar-power/107688/2 "2023-12-15T22:07:50Z")
> [@valentin.m](#):
>
> A key part of the algorithm involves computing the power of large arrays of floats raised to a scalar float.
Do you mean the elementwise power, i.e. `x .^ p`? ~~If `p` is a scalar and `x` is an array, then it is probably faster to compute `logp = log(p)` and then do `exp.(logp .* x)`. There are SIMD-accelerated functions for this kind of thing, too.~~ Whoops, sorry, `x .^ p` is `exp(p .* log.(x))`.
Are you sure that your speed is limited by this operation alone? Usually if you have a large array then you want to do many operations per element in a single pass, in order to utilize memory efficiently. See also [More dots: Why vectorized code is not as fast as it could be](https://julialang.org/blog/2017/01/moredots/#why_vectorized_code_is_not_as_fast_as_it_could_be). From this perspective, focusing on a single “vectorized” operation _in isolation_ is a common mistake.
---
### Author:  [@Jeff\_Emanuel](https://discourse.julialang.org/u/Jeff_Emanuel)
#### Post date: [December 15, 2023, 11:30pm UTC](https://discourse.julialang.org/t/efficient-ways-to-raise-an-array-by-a-scalar-power/107688/3 "2023-12-15T23:30:44Z")
> [@stevengj](#):
>
> `logp = log(p)` and then do `exp.(logp .* x)`.
I think you have the identity reversed
```julia
julia> x=[3, 5]
2-element Vector{Int64}:
3
5
julia> p=2
2
julia> x.^p
2-element Vector{Int64}:
9
25
julia> logp = log(p); exp.(logp .* x)
2-element Vector{Float64}:
7.999999999999998
32.0
julia> exp.(log.(x) .* p)
2-element Vector{Float64}:
9.000000000000002
24.999999999999996
```
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [December 16, 2023, 1:39am UTC](https://discourse.julialang.org/t/efficient-ways-to-raise-an-array-by-a-scalar-power/107688/4 "2023-12-16T01:39:05Z")
> [@Jeff\_Emanuel](#):
>
> I think you have the identity reversed
Oh, whoops, right, sorry!
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [December 16, 2023, 2:30am UTC](https://discourse.julialang.org/t/efficient-ways-to-raise-an-array-by-a-scalar-power/107688/5 "2023-12-16T02:30:54Z")
one other thing worth mentioning is that if you can do the computation in float32, that will likely be a good bit faster.
---
### Author:  [@Salmon](https://discourse.julialang.org/u/Salmon)
#### Post date: [December 18, 2023, 7:26am UTC](https://discourse.julialang.org/t/efficient-ways-to-raise-an-array-by-a-scalar-power/107688/6 "2023-12-18T07:26:25Z")
I am very surprised that this is expected to be faster.
Shouldn’t it be trivial to SIMD `x .^p`, possibly trivial enough for the compiler to do? Or at least reasonably quick with Loop vectorization.jl? Or is there further magic that somehow makes the approach using exp and log faster? Without benchmarking I’d expect even a naively vectorized version using powers (that might even be optimized to multiplications) to be faster than computing exp and log numerically.
What am I missing?
Using log also has the disadvantage of breaking for negative elements of x.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [December 18, 2023, 2:18pm UTC](https://discourse.julialang.org/t/efficient-ways-to-raise-an-array-by-a-scalar-power/107688/7 "2023-12-18T14:18:16Z")
> [@Salmon](#):
>
> I am very surprised that this is expected to be faster. Shouldn’t it be trivial to SIMD `x .^p`?
Computing `x ^ p` for non-integer `p` is not trivial even without SIMD, and the algorithm is complicated enough that I don’t think the compiler can SIMD-ize it.
To be clear, however, while `exp(x)` is faster than `x ^ p`, computing `exp(p * log(x))` is slower. I had at first mistakenly thought the `log` computation could be pulled out of the loop. Also, I knew there were SIMD-optimized versions of `exp` and `log`.
However, it looks like IntelVectorMath.jl has functions `pow` and `pow!` that compute `x ^ p`, and similarly for AppleAccelerate.jl, so I would try those.
> [@Salmon](#):
>
> Using log also has the disadvantage of breaking for negative elements of x.
It’s no different from `x .^ p` where `p` is not an integer, which also throws an error for negative real `x`.
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [December 18, 2023, 4:40pm UTC](https://discourse.julialang.org/t/efficient-ways-to-raise-an-array-by-a-scalar-power/107688/8 "2023-12-18T16:40:55Z")
There are various tricks to speed up computations, e.g. [Fast inverse square root - Wikipedia](https://en.wikipedia.org/wiki/Fast_inverse_square_root) was made famous in Doom by John Carmack.
If you share the distribution of values in `x` and what type of `p` to expect and when it is known, maybe it would fascilitate such approximation magic.