# Review: pth matrix roots

**URL:** <https://discourse.julialang.org/t/review-pth-matrix-roots/1866>\
**Category:** General Usage\
**Tags:** performance, review\
**Created:** [February 3, 2017, 2:01pm UTC](https://discourse.julialang.org/t/review-pth-matrix-roots/1866 "2017-02-03T14:01:43Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![felix](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/felix/32/604_2.png) [@felix](https://discourse.julialang.org/u/felix)\
**Post date:** [February 3, 2017, 2:01pm UTC](https://discourse.julialang.org/t/review-pth-matrix-roots/1866/1 "2017-02-03T14:01:43Z")

</div>

Hello,

to test the accuracy of the [Schur-Pade](https://discourse.julialang.org/t/review-schur-pade-matrix-powers-speedup/1650) approximation of `A^r` for `A` a matrix and `r` a real number, I implemented a method [`rootm`](https://github.com/felixrehren/julia/blob/fr-rootm/base/linalg/rootm.jl) to calculate `A^(1/q), q::Int`. This is a direct generalisation of `sqrtm = A^(1/2)` ([PR20214](https://github.com/JuliaLang/julia/pull/20214)) and would allow the accurate computation of `A^(p//q)` for any integers `p,q`. The algorithm solves `X^p = A` for `X` by getting a recurrence relation that directly derives from writing out the product `X^p`. **Anyway, the performance of [this code](https://github.com/felixrehren/julia/blob/fr-rootm/base/linalg/rootm.jl) sucks.** Why?

```
A = randn(127,127)
A = UpperTriangular(schurfact(A'*A)[:T])

@benchmark _sqrtm(A)
  median time: 744.535 μs (0.00% GC)
@benchmark _rootm(A,2,Val{true})
  median time: 145.209 ms (1.93% GC)

@profile _rootm(A,2,Val{true})
ProfileView.view()

```

 ![](https://global.discourse-cdn.com/julialang/original/3X/6/a/6a2445f250f548066ae7a15e961d05b9350acfe0.png)

Using `@code_warntype`, I get no red ink. The algorithm is not _that_ much more complex than `sqrtm`, but it is 200x slower. I don’t understand the profiling – where does the inference come from, and why are there several distinct calls to `_rootm` in the above profile despite the function being called once?

Performance for `p`th roots with `p>2` is obviously the more interesting point, but `p=2` is a good benchmark and I’d like to understand the performance difference.

Any comments and advice appreciated

---

<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 3, 2017, 2:38pm UTC](https://discourse.julialang.org/t/review-pth-matrix-roots/1866/2 "2017-02-03T14:38:05Z")

</div>

Two suggestions that occur to me when I glance at your code:

- Don’t call `sum` in your inner loops. Write out the loops.

- Don’t compute powers as `x^q` or similar in a loop over `q`: in each loop iteration, accumulate the product by multiplying repeatedly by `x` (essentially, use Horner’s method).

For example, replace `xij /= sum(xii^(p-1-q)*xjj^q for q in 0:(p-1))` by something like:

```julia
# compute xij /= sum(xii^(p-1-q)*xjj^q for q in 0:(p-1)):
∑ = ∏ = xii^(p-1)
ξ = xjj / xii
for q = 1:p-1
    ∏ *= ξ
    ∑ += ∏
end
xij /= ∑

```

---

<div class="post-metadata">

**Author:** ![felix](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/felix/32/604_2.png) [@felix](https://discourse.julialang.org/u/felix)\
**Post date:** [February 3, 2017, 4:28pm UTC](https://discourse.julialang.org/t/review-pth-matrix-roots/1866/3 "2017-02-03T16:28:06Z")

</div>

That worked! 🎉

I’m amazed:

```
  median time: 781.725 μs (0.00% GC)

```

That is within 10% of the `sqrtm` function! The performance gain is the same for `rootm(A,p)` with `p>2`. Thank you Steven!

Almost all of this most have come from eliminating the `sum`s, because in this case `p=2` and `1 <= q < p` means `q = 1`, the cumulative powers eliminated using Horner’s method can’t have done much for the speedup. Is there analysis anywhere why/when `sum` is slow in Julia?

---

<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 10, 2017, 1:28am UTC](https://discourse.julialang.org/t/review-pth-matrix-roots/1866/4 "2017-02-10T01:28:53Z")

</div>

sum is a perfectly good function, but you have to realize that in a sum of only three elements _any_ call of a general summation function (or any other function call) is going to have a lot of overhead compared to the trivial amount of work it is doing (just two additions). Two floating point additions are _fast_.

The cost of setting up a generator object adds even more overhead. If you were summing many numbers, this wouldn’t matter. But compared to the cost of two additions the creation of any complex object is massively costly.

Finally, computing a general exponentiation operation (albeit for an integer exponent, but which still involves a function call and various checks e.g. of the exponent sign and magnitude) for every summand is vastly more expensive than a single multiplication.

The rule of thumb is that the cheaper the operation you are performing, the more careful you have to be with high level abstractions.

---

<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 10, 2017, 2:18am UTC](https://discourse.julialang.org/t/review-pth-matrix-roots/1866/5 "2017-02-10T02:18:58Z")

</div>

(Looks like you also hit a type-inference bug: [https://github.com/JuliaLang/julia/issues/20517](https://github.com/JuliaLang/julia/issues/20517))