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.

It first generates all of the N (N-1)/2 random numbers with a single SIMD-vectorized rand! call in a contiguous temporary buffer v, and then fills this data into the matrix with a cache-oblivious algorithm that divides the array recursively into smaller and smaller blocks until ending up at tiny blocks 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)
[details="Code"]

using LinearAlgebra, Random, BenchmarkTools

function _prepare3_diag!(A, i0, n, v, k)
if n <= 16
for i = i0:i0+(n-1)
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

[/details]