Efficient creation of power series matrix or array of arrays

General Usage Performance
7

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, and the fastest way to do this is probably just two loops. For example:

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.

7 Likes