Preparing an antisymmetric matrix of normal numbers as efficiently as possible

General Usage Performance
1

My goal is to prepare a fully antisymmetric matrix of normal numbers, on a single thread, as efficiently as possible. After numerous manual attempts, here is the best version I could come up with:

using BenchmarkTools
const N = 1024
const TAB = zeros(Float64, N, N)
function prepare!()
    @inbounds for i=1:N
        for j=1:(i-1)
            val = randn()
            TAB[i,j] =  val
            TAB[j,i] = -val # Antisymmetry
        end
        TAB[i,i] = 0.0 # Diagonal
    end
end 
#####
@btime prepare!() # -> 1.921 ms

I am wondering whether there could exist more efficient implementations, e.g., with better data locality, or using SIMD for the generation of the random numbers?

Looking forward to your input and all your help!

2

I think it would be faster to do two passes: first filling up the upper triangular, and second copying it to the lower triangular with a minus sign.

The reason is to replace out-of-place writes (which are expensive) with out-of-place reads (which are less expensive).

3

By any chance, would you have a MWE for this approach? Indeed, my implementation of your suggestion seems slower [on my laptop at least]

function prepare_triangular!()
    @inbounds for j=1:N
        for i=(j+1):N
            TAB[i,j] = randn()
        end
        TAB[j,j] = 0.0
    end
    @inbounds for j=1:N
        for i=1:(j-1)
            TAB[i,j] = - TAB[j,i]
        end
    end
end
#####
@btime prepare_triangular!() # -> 2.278 ms
4

You really want to do neither, in principle — it seems similar to optimizing a transpose operation, in which the best approaches involve some form of blocking/tiling in order to maximize cache-line utilization. See e.g. WIP: cache oblivious linear algebra algorithms - Pull Request #6690 - JuliaLang/julia - GitHub

(Though I’m surprised that generating anti-symmetric normal numbers is the performance-critical step in a real application. Aren’t you going to do something nontrivial with the matrix that will dominate over the benefits of micro-optimizing this? Or is this just an exercise?)

5

Here is a little recursive cache-oblivious implementation of this idea. (Recursion can be fast!)

It first generates all of the N (N-1)/2 random numbers with a single SIMD-vectorized randn! call in a contiguous temporary buffer v, and then copies this data anti-symmetrically into the matrix with a cache-oblivious algorithm that divides the array recursively into smaller and smaller blocks until ending up at tiny blocks (that should fit into L1 cache but big enough to amortize the recursion overhead) where it uses a simple loop. (I experimented a little and 16 seems like a good block size on my machine.)

julia> N = 1024; A = zeros(N, N); v = zeros((N * (N-1)) >> 1);

julia> @btime prepare!($A); # your original version, fixed to take A as an argument
  1.633 ms (0 allocations: 0 bytes)

julia> @btime prepare3!($A, $v);
  817.416 μs (0 allocations: 0 bytes)

julia> @btime randn!(v); # lower bound: just the raw random-number generation cost
  553.708 μs (0 allocations: 0 bytes)

I didn’t put many comments in my code, so I’ll leave deciphering it as an exercise :wink:. (Caveat: it’s only lightly tested.)

Code 
using LinearAlgebra, Random, BenchmarkTools

function _prepare3_diag!(A, i0, n, v, k)
    if n <= 16
        for i = i0:i0+(n-1)
            @inbounds A[i,i] = 0
            for j = i0:i-1
                @inbounds A[i,j] = -(A[j,i] = v[k += 1])
            end
        end
    else
        n2 = n >> 1
        k += _prepare3_diag!(A, i0, n2, v, k)
        k += _prepare3_diag!(A, i0+n2, n-n2, v, k)
        _prepare3_offdiag!(A, i0, i0+n2, n2, n-n2, v, k)
    end
    return (n * (n-1)) >> 1 # number of values consumed
end

function _prepare3_offdiag!(A, i0, j0, nr, nc, v, k)
    if nr * nc <= 256
        for i = i0:i0+(nr-1), j = j0:j0+(nc-1)
            @inbounds A[i,j] = -(A[j,i] = v[k += 1])
        end
    else
        nr2 = nr >> 1
        nc2 = nc >> 1
        k += _prepare3_offdiag!(A, i0, j0, nr2, nc2, v, k)
        k += _prepare3_offdiag!(A, i0+nr2, j0, nr-nr2, nc2, v, k)
        k += _prepare3_offdiag!(A, i0, j0+nc2, nr2, nc-nc2, v, k)
        _prepare3_offdiag!(A, i0+nr2, j0+nc2, nr-nr2, nc-nc2, v, k)
    end
    return nr * nc
end

function prepare3!(A::AbstractMatrix, v::AbstractVector)
    N = LinearAlgebra.checksquare(A)
    length(v) == ((N * (N-1)) >> 1) || throw(DimensionMismatchError())
    randn!(v) # SIMD-vectorized random-number generation
    _prepare3_diag!(A, 1, N, v, 0)
    return A
end

Note, by the way, that a power-of-two size like N=1024 above is usually the worst case for cache-associativity conflicts. N=1023 is substantially faster for both the original code and the cache-oblivious version:

julia> N = 1023; A = zeros(N, N); v = zeros((N * (N-1)) >> 1);

julia> @btime prepare!($A);  # original version, fixed to take A as an argument
  1.387 ms (0 allocations: 0 bytes)

julia> @btime prepare3!($A, $v);
  758.750 μs (0 allocations: 0 bytes)

(This is all on an Apple M4 from 2024.)

Of course, this doesn’t get into threading, but in principle you could spawn threads from the top-level recursive subdivisions to do them in parallel while still maintaining locality with each thread. Not sure how big N has to be for that to be beneficial, though.

6

You’re assigning to the diagonal out-of-place, which clearly doesn’t help, you should have the line TAB[j, j] = 0.0 before the for loop, not after.

If that doesn’t work I have no idea what’s wrong.

stevengj is right, of course, but both ideas can be combined, see mini-optimization of transposition - Pull Request #1640 - JuliaLang/LinearAlgebra.jl - GitHub

7

I tried this with my code and it doesn’t seem to make a significant difference for my cache-oblivious implementation? (Benchmarks are about the same.)

Code 

# same as above, but exclusively in-order writes
function _prepare5_diag!(A, i0, n, v, k)
    if n <= 16
        for i = i0:i0+(n-1)
            for j = i0:i-1
                @inbounds A[j,i] = v[k += 1]
            end
            @inbounds A[i,i] = 0
        end
        for i = i0:i0+(n-1), j = i+1:i0+(n-1)
            @inbounds A[j,i] = -A[i,j]
        end
    else
        n2 = n >> 1
        k += _prepare5_diag!(A, i0, n2, v, k)
        k += _prepare5_diag!(A, i0+n2, n-n2, v, k)
        _prepare5_offdiag!(A, i0, i0+n2, n2, n-n2, v, k)
    end
    return (n * (n-1)) >> 1 # number of values consumed
end

function _prepare5_offdiag!(A, i0, j0, nr, nc, v, k)
    if nr * nc <= 256
        for i = i0:i0+(nr-1), j = j0:j0+(nc-1)
            @inbounds A[j,i] = v[k += 1]
        end
        for j = j0:j0+(nc-1), i = i0:i0+(nr-1)
            @inbounds A[i,j] = -A[j,i]
        end
    else
        nr2 = nr >> 1
        nc2 = nc >> 1
        k += _prepare5_offdiag!(A, i0, j0, nr2, nc2, v, k)
        k += _prepare5_offdiag!(A, i0+nr2, j0, nr-nr2, nc2, v, k)
        k += _prepare5_offdiag!(A, i0, j0+nc2, nr2, nc-nc2, v, k)
        _prepare5_offdiag!(A, i0+nr2, j0+nc2, nr-nr2, nc-nc2, v, k)
    end
    return nr * nc # number of values consumed
end

function prepare5!(A::AbstractMatrix, v::AbstractVector)
    N = LinearAlgebra.checksquare(A)
    length(v) == ((N * (N-1)) >> 1) || throw(DimensionMismatchError())
    randn!(v) # SIMD-vectorized random-number generation
    _prepare5_diag!(A, 1, N, v, 0)
    return A
end
8

Since the transposition uses the same block-recursive algorithm it must have the same advantage. I suppose it’s not visible here because the advantage is small and the computational cost is dominated by generating the random numbers.

9

Thank you very much for your feedback.
In practice, once this “noise” generated, it is used within a fast SIMD loop of the form (purposely simplified)

const VEC_IN = rand(Float64,N)
const MAT_IN = rand(Float64,N,N)
const VEC_OUT = zeros(Float64,N,N)
function compute!()
    @inbounds for j=1:N
        acc = 0.0
        @inbounds @simd for i=1:N
            acc += MAT_IN[i,j] + VEC_IN[i] * TAB[i,j]
        end
        VEC_OUT[j] = acc
    end
end
#####
@btime compute!() # -> 242.958 μs

So, as surprisingly as it may sound, the generation of the noise is currently the bottleneck of my code

10

Fantastic! I confirm that I could reproduce the performance gains on my side as well :slight_smile:
Following this new code, I had three remaining questions:

1/ In practice, I will be generating two such “noise” matrices. For the sake of performance, should I call the function prepare3! twice, or should I update this function so that the filling of the two matrices are done during the same pass?

2/ I noticed that, when changing the block size to 4, 8, 12, 16, 20, 24, 28, the acceleration of the code was, in essence, identical. And, it started to drop for 32. Was it obvious that such a wide range of block size would lead to similar code performance (on my laptop at least)?

3/ Finally, as I was playing around with the initial version of the code, I noticed that this simpler code [which does not draw any normal number!] had (much) worse performance

using BenchmarkTools
const N = 1024
const TAB = zeros(Float64, N, N)
function prepare_trivial!()
    @inbounds for i=1:N
        for j=1:(i-1)
            val = 1.0
            TAB[i,j] =  val
            TAB[j,i] = -val # Antisymmetry
        end
        TAB[i,i] = 0.0 # Diagonal
    end
end 
#####
@btime prepare_trivial!() # -> 2.547 ms

Unfortunately, I am not [yet ;)] strong enough in julia to dissect what could be causing the problem from @code_native and others… Is this an unexpected bug?

11

Yes, because this is a cache-oblivious algorithm (read e.g. Frigo (2012) on the theory), so the base case isn’t the same as the cache-optimized “block” size in a traditional cache-aware blocking/tiling algorithm.

The base case needs to be large enough to make the recursion overhead negligible, and small enough to fit easily inside the smallest cache, but otherwise it isn’t so critical.

The recursion divides the computation into a tree of nested blocks, and once a recursion level fits inside the cache, subsequent recursive calls incur no further cache issues. (This also means that it automatically exploits multiple levels of cache.)

For me it’s consistently about 6% slower than your original prepare!, which is only a slight slowdown. But it’s still mysterious to me that replacing val = randn() with val = 1.0 slows things at all!

12

This would be quite cool to have in SkewLinearAlgebra.jl

Im picturing an interface something like

skewhermitian(N) do i, j
    # Some i, j dependant calculation 
end
skewhermitian!(H, v) do i, j
    # ...
end