Alternative way to create following array raised to powers?

How to do following more efficiently (any alternative way ?)

julia> a = [-1,0,1]  # my  array
3-element Vector{Int64}:
 -1
  0
  1

julia> hcat([a.^i for i in 0:3]...)'   # array raised to power 0:n and then transposed
4×3 adjoint(::Matrix{Int64}) with eltype Int64:
  1  1  1
 -1  0  1
  1  0  1
 -1  0  1

Try this:

reduce(hcat, [a.^i for i in 0:3])' 
2 Likes

You can use broadcasting with one dimension for each input array:

julia> [-1 0 1] .^ (0:3)
4×3 Matrix{Int64}:
  1  1  1
 -1  0  1
  1  0  1
 -1  0  1
5 Likes

Wow… How I missed this I don’t understand… I was trying same but was writing array as [-1, 0,1]
It was giving error :sweat_smile:

[-1, 0, 1]' .^ (0:3)

But frankly, I’m quite surprised at the way it works, i.e. returning a matrix instead of a vector of vectors, for example…

1 Like

How so? This is how broadcasting always works. f.(row, column) always returns a matrix. And not just in Julia. Broadcasting in Python and Matlab does the same except for the dot notation.

1 Like

Possibly because [1 0 1] is a 1x3 Matrix , and not a vector… So broadcasting a matrix type is giving us a Matrix itself?

@DNF, thank you for framing this properly as I totally missed the point, but can now see this is a common pattern.

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

@stevengj @sudete
Both of your solutions are useful as well as good learning experience in julia. I wish discourse had ability to mark both solutions as Solution :sweat_smile: