# Efficient creation of power series matrix or array of arrays

**URL:** <https://discourse.julialang.org/t/efficient-creation-of-power-series-matrix-or-array-of-arrays/18988>\
**Category:** Performance\
**Tags:** question\
**Created:** [December 26, 2018, 1:32am UTC](https://discourse.julialang.org/t/efficient-creation-of-power-series-matrix-or-array-of-arrays/18988 "2018-12-26T01:32:49Z")\
**Posts on this page:** 1\
**Showing post:** 7

<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:** [December 27, 2018, 6:37pm UTC](https://discourse.julialang.org/t/efficient-creation-of-power-series-matrix-or-array-of-arrays/18988/7 "2018-12-27T18:37:11Z")

</div>

> [@Tamas\_Papp](#):
>
> Unless you just need the first few powers, coding a doubling (divide and conquer) algorithm is the best way to ensure that you get the optimal algorithm.

Not if you need _all_ of the powers in a _sequence_—then it is faster to multiply by `x` one at a time. (Note that exponentiation-by-squaring is already the algorithm used for `z^n` in Julia.)

It sounds like he just wants a [Vandermonde matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix), and the fastest way to do this is probably just two loops. For example:

```julia
function vander!(V::AbstractMatrix, x::AbstractVector, n=length(x))
    m = length(x)
    (m,n) == size(V) || throw(DimensionMismatch())
    for j = 1:m
        @inbounds V[j,1] = one(x[j])
    end
    for i = 2:n, j = 1:m
        @inbounds V[j,i] = x[j] * V[j,i-1]
    end
    return V
end
vander(x::AbstractVector, n=length(x)) = vander!(Array{eltype(x)}(undef, length(x), n), x, n)

```

Note that, as is common in Julia, I provide a `vander` routine that allocates an output matrix, and also a `vander!` routine that operates on a pre-allocated matrix (which is faster if you need to do this computation many times for different `x` with the same dimensions).

There are various ways you could modify the above to be faster (by optimizing cache-line or SIMD utilization, for example), but that kind of micro-optimization is only worth it for critical routines.

In contrast, the _simplest_ code to do the above is probably `vander(x, n) = x .^ (0:n)'`, which uses broadcasting. However, this is significantly less efficient because it computes each power of `x` independently (about 20× slower for a 1000×1000 `Float64` Vandermonde matrix on my machine).

PS. Change `one(x[j])` to `x[j]` if you want the first column to be `x^1` and not `x^0`.

---

_[View the full topic](https://discourse.julialang.org/t/efficient-creation-of-power-series-matrix-or-array-of-arrays/18988)._