Efficient creation of power series matrix or array of arrays

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.

1 Like

In my opinion, FPO count is interesting from a theoretical perspective, but not a good measure when it comes to actual performance. Things like memory accesses (and caching thereof), branching (and prediction thereof), SIMD, loop unrolling, and characteristics in the input data, make FPO count a very unreliable measure. As suggested above, I would instead just use BenchmarkTools and the various code_ methods to analyze the code.

In your example, the problem can be solved in a single FPO per element, since to create row k, you can simply multiply row k-1 with the first row. Example:

julia> x = [1,2,3,4];

julia> r = copy(x);

julia> for n = 1:3 println(r .*= x) end
[1, 4, 9, 16]
[1, 8, 27, 64]
[1, 16, 81, 256]

In fact, with SIMD, multiple elements can be handled at the same time, leading to less than a single operation per element.

To avoid the memory access (reading the previous row), and to better exploit the way Julia stores matrices (column-major order), you should probably build the matrix column by column. To take advantage of SIMD then, you could either construct a few columns at a time, or likely better, construct several rows at a time within the same column. E.g. consider the first 8 powers of 3:

julia> println(3 .^ (1:8))
[3, 9, 27, 81, 243, 729, 2187, 6561]

If you can process 4 elements at a time, the next 4 can be constructed by multiplying the previous 4 with the constant value 3^4:

julia> println([3, 9, 27, 81] .* 81)
[243, 729, 2187, 6561]

Thanks for taking time to respond and for the suggestions, I’ll look into SIMD, and appreciate the alternative you proposed.

The FPO count is a piece of the jigsaw more generally and I was hoping to explore it and Julia in this simple setting - comparing what Julia ends compiling to what was coded.

As mentioned before, ack: In some (many?) cases it is likely swamped by other differences/issues that are implementation specific.

It is a separate topic but given so much effort can be expended on different algorithms, I’m curious about whether such differences survive compilation, then whether they exert influence greater than the other issues you mention. FPO count is one piece of the (complex) jigsaw.

I’ll definitely be using BenchmarkTools.jl in the way you suggest.

Thanks again.

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

I often use this code to create a Vandermonde matrix to help introduce Julia, by the way. The two-loop version is pretty compact and readable, is completely generic to any array-like container and any number-like type supporting one and * (including user-defined types), but performs similarly to the optimized numpy.vander routine on arrays of Float64.

If you look at the numpy.vander implementation, in contrast, it is a thin Python wrapper around a rather long and complicated C routine that is itself basically a shim to dispatch to the correct inner-loop kernel; the lowest level kernels are generated from a custom C-like template syntax. The reason for all of that complexity is to handle a multiplicity of types and array formats, but even so it only handles one container type (NumPy arrays) and a small number of “built-in” numeric scalar types.

11 Likes

Nice! I posted this method above btw, but in a vector format instead of a matrix.

I find this statement a bit odd/misleading – I consider cache locality crucial in algorithms like this, and not at all a micro-optimization. To see why, reverse the two for loops like this: for j = 1:m, i = 2:n and re-run and you’ll see the performance drop 10-fold! Same algorithm, same FPO count, but it executes 10 times slower. That’s only twice as fast as the naive solution, meaning that cache locality can be more important than picking the right algorithm. (The reason for this is explained in more detail in the link in my previous post.)

SIMD can also lead to enormous improvements, there was a recent topic where we played with vectorized and branchless code and were able to make some sample code several 100 times faster. The code you posted will already be SIMD-vectorized (and loop-unrolled) automatically by the compiler btw, so it’s not something you need to enable manually.

But perhaps you knew all of this already, and meant that although your code is already high-performing, you could still optimize it further

Exactly. For example, optimizing cache lines further is not just a matter of reordering the loops, it requires blocking so that you do several outer-loop (i) iterations in every inner-loop (j) iteration. And to get the most out of SIMD the compiler typically needs a little help (I didn’t even use @simd above because it would have required me to break up the for statement).

It’s not so much that the performance gains from further optimizations of this sort are small, it’s that the tradeoff in code size/complexity may not be worth it for non-critical code.

1 Like

People interested in these issues, likely, will be interested in the Following Post by Tim Holy, updated Dec 2018 for Julia 1.1:

Multidimensional algorithms and iteration

I’m running some benchmarks that confirm the above discussions. When I have managed to work out the use of @simd in the code Steven posted I post the results here.

Thanks Steve for your pointer to Vandermonde matrix and your intro to Julia notes.

The use of SIMD is beyond my ken right now, but the improvement made by Steven’s code (and SpecifalMatrices.jl) would be enough for immediate needs.

julia> include("~/src/scratch/vandermonde.jl")
  19.904 ms (1004 allocations: 7.77 MiB)  # My naive map attempt
  2.612 ms (3 allocations: 7.63 MiB)      # SpecialMatrices.jl implementation
  12.217 ms (3 allocations: 7.63 MiB)     # SpecialMatrices.jl re-ordered for loop
  2.097 ms (2 allocations: 7.63 MiB)      # stevengj implementation
  12.770 ms (2 allocations: 7.63 MiB)     # stevengj re-rdered for loop

The code used is:

using BenchmarkTools

# The Vandermonde definitions and implementation using it comes from:
# https://github.com/JuliaMatrices/SpecialMatrices.jl
struct Vandermonde{T} <: AbstractArray{T, 2}
    c :: Vector{T}
end

getindex(V::Vandermonde, i::Int, j::Int) = V.c[i]^(j-1)
isassigned(V::Vandermonde, i::Int, j::Int) = isassigned(V.c, i)

size(V::Vandermonde, r::Int) = (r==1 || r==2) ? length(V.c) :
    throw(ArgumentError("Invalid dimension $r"))
size(V::Vandermonde) = length(V.c), length(V.c)

function vander1(V::Vandermonde{T}) where T
    n=size(V, 1)
    M=Matrix{Float64}(undef, n, n)
    M[:,1] .= 1
    for j=2:n
        for i=1:n
	    M[i,j] = M[i,j-1]*V.c[i]
        end
    end
    M
end
# This is vander1 with the order of the loop reversed
function vander2(V::Vandermonde{T}) where T
    n=size(V, 1)
    M=Matrix{Float64}(undef, n, n)
    M[:,1] .= 1
    for i=1:n
        for j=2:n
	    M[i,j] = M[i,j-1]*V.c[i]
        end
    end
    M
end

function vander3(x::AbstractVector, n=length(x))
    m = length(x)
    V=Matrix{Float64}(undef, m, n)
    V[:,1] .= 1
    for i = 2:n, j = 1:m
        @inbounds V[j,i] = x[j] * V[j,i-1]
    end
    return V
end
# This changes the vander3 loop order: result is approx 6x slower
function vander4(x::AbstractVector, n=length(x))
    m = length(x)
    V=Matrix{Float64}(undef, m, n)
    V[:,1] .= 1
    for j = 1:m, i = 2:n
        @inbounds V[j,i] = x[j] * V[j,i-1]
    end
    return V
end

function vander0(x::AbstractVector)
    m = length(x)
    map(i -> x .^i ,collect(1:m))
end

n = 1000

# @btime const M=Matrix{Float64}(undef, n, n)
# @btime M[:,1] .= 1

a = collect(1:n)

@btime V0 = vander0($a);
@btime V1 = vander1(Vandermonde($a));
@btime V2 = vander2(Vandermonde($a)); 
@btime V3 = vander3($a); 
@btime V4 = vander4($a); 

This is off-topic, but are Vandermonde matrices used in practice in numerical linear algebra applications? I thought they were very badly conditioned (except in some special cases). But since you are using them in a course, I must be mistaken.

I had the same question, tuns out: representing Discrete Fourier Transforms and fitting polynomials to points. Likely others too.

I know about DFT (one of the special cases), but for fitting general polynomials it has numerical problems, and in any case solving a linear system may not be as efficient as using eg orthogonal families, which are incidentally better conditioned.

I mainly use them for low-degree polynomial least-squares, where conditioning is not a problem. For Julia introductions, though, they are a nice example because they are easy to explain, easy to implement, and are available as a “built-in” library function in many high-level languages for comparison.

(I’ve never seen a case where the construction of the Vandermonde matrix was performance-critical, though.)

1 Like

It’s old topic, but maybe someone could help me.

I try run this function in Julia, but I have no idea what should I write and how induce it.

function matrix_vandermonde(x::AbstractVector, n=length(x))
    m = length(x)
    V = Array(m, n)
    for j = 1:m
        V[j,1] = one(x[j])
    end
    for i = 2:n
        for j = 1:m
            V[j,i] = x[j] * V[j,i-1]
        end
    end
    return V
end

and I call in Julia matrix_vandermonde(???) what should I write to have any effect execpt “matrix_vandermonde (generic function with 2 methods)”

First, the code you quoted (remember to quote your code!) isn’t correct: you can’t create an array with Array(m,n) because you have to give an element type. In this case, since you probably want to have the same element type as x, you could change the line V = Array(m, n) to

V = similar(x, m, n)

Then you call it by passing an array (for the x) argument, and it makes a Vandermonde matrix from that array, for example:

julia> matrix_vandermonde([1,2,3,4])
4×4 Array{Int64,2}:
 1  1   1   1
 1  2   4   8
 1  3   9  27
 1  4  16  64
2 Likes

Hi,

thank you for help.
Unfortunaly, when I change it and call in way that you suggested I get this:

ERROR: UndefVarError: lenght not defined
Stacktrace:
[1] vander(::Array{Int64, 1}) at .\REPL[2]:2 
[2] top-level scope at none:0   

what was wrong ?

I think it should be spelled as length .

1 Like

Oh, right, my bad.

Thank you so so so much. It works now! :slight_smile:

You can type len<TAB> to use auto-completion so that you can avoid this problem.

[I have noticed that Spanish speakers very frequently write `lenght` instead of `length`. I haven’t worked out why that might be, nor how to help them unlearn that.]

2 Likes

Thanks for tip, I’ll use it!

[I’m not Spanish speaker, but it was a typo and I didn’t see difference between my original code. And, it could because f.e different words has similar end - " thought, lightweight", For someone who is just intemediate it’s not a big difference, unfortunately :(]