Alternative way to create following array raised to powers?

What you are constructing is a Vandermonde matrix (transposed). A much faster way to construct it is to use two nested loops, as described in another discourse thread.

The two-loop approach is about 30x faster on my machine for constructing a 100x100 Vandermonde matrix, compared to the broadcasting approach suggested above:

julia> using BenchmarkTools

julia> vander(x, n=length(x)) = .... # two-loop version from discourse

julia> bvander(x, n=length(x)) = x .^ (0:n-1)';  # broadcast-based vandermonde

julia> x = rand(100);

julia> @btime vander($x);
  5.892 μs (2 allocations: 78.17 KiB)

julia> @btime bvander($x);
  193.366 μs (2 allocations: 78.17 KiB)

The reason the two-loop version can be so much faster is that it doesn’t compute each power separately, but instead accumulates the powers one multiplication at a time from the previous powers.

(The broadcasting version is much more compact, though, so if it’s not too performance-critical I would stick with that!)

Another option is to create the matrix lazily, e.g. SpecialMatrices has methods to implicitly create such matrices and work with them. The main advantage here is that if you are doing things like solving linear systems, there are specialized algorithms for Vandermonde matrices.

7 Likes