Correct way to benchmark StaticVectors

General Usage
,
1

I am writing some code to do geometric intersection tests using StaticArrays.jl. I would like to keep some benchmarks as I am developing and would like to know if I am benchmarking static arrays correctly.

The documentation for BenchmarkTools mentions that operations that can be simplified at compile time will not be measured properly, and that you therefore have to reference and dereference the interpolated variable. I think I am doing that properly below, however the measured output is still very fast for this operation (around 6 ns), so I am suspicious that I am not measuring the right thing.

Would the below benchmark correctly measure the execution speed of the BVH.intersection function, and would anyone be able to explain / point to a source on how this interpolation and referencing / dereferencing works?

using BenchmarkTools
using Revise, StaticArrays, LinearAlgebra

function benchmark_intersection()
    for T in [Float64, Float32]

        v₁ = SVector{3, T}(1.0, 0.0, 1.0)
        v₂ = SVector{3, T}(-1.0, 1.0, 1.0)
        v₃ = SVector{3, T}(-1.0, -1.0, 1.0)

        O = SVector{3, T}(0.0, 0.0, 0.0)
        direction = SVector{3, T}(0.0, 0.0, 1.0)

        A = [-direction (v₂-v₁) (v₃-v₁)]
        b = O - v₁

        ray = Ray(O, direction)
        triangle = Triangle(v₁, v₂, v₃)
          
        println("Benchmarking $T")
        display(@benchmark intersection($(Ref(ray))[], $(Ref(triangle))[]))
    end
end

struct Ray{N, T}
    origin::SVector{N, T}
    direction::SVector{N, T}
end

struct Triangle{N, T}
    v₁::SVector{N, T}
    v₂::SVector{N, T}
    v₃::SVector{N, T}
end

function intersection(ray::Ray{N, T}, triangle::Triangle{N, T}) :: IntersectionResult{N, T} where {N, T}
    (; origin, direction) = ray
    (; v₁, v₂, v₃) = triangle

    A = [-direction (v₂-v₁) (v₃-v₁)]
    b = origin - v₁

    if det(A) ≈ 0
        return IntersectionResult(false, nothing)
    end

    (t, _, _) = A \ b
    point = origin + direction * t

    return IntersectionResult(true, point)
end
2

What is BVH?

3

The module I am writing - sorry forgot to add that code. Edited above to include.

4

Benchmarking looks fine to me. Generally, timings <1ns are suspicious and indicate a complete constant-folding of your function. But here you are using StaticArrays.jl with vectors of size 3 so each function are just a handful of FLOPS and modern CPUs are quite fast :slight_smile:

However this

is nonsense because it will always evaluate to false as just being discussed in another thread.
You’ll need to think about an absolute scale, e.g. given by the largest element of A and then set abstol accordingly. Though I am not sure what “accordingly” should mean concretely when considering det. @stevengj probably has an idea what an appropriate scale should be. I would propose abstol=sqrt(eps(eltype(A)))*maximum(A)

1 Like
5

That won’t work, because it has the wrong units. (The units of the determinant of an m \times m matrix are the matrix elements to the m'th power.) An absolute tolerance has to have the same units as the quantity you are comparing it to.

Also maximum(A) returns the largest element algebraically, whereas I think you want the largest magnitude.

The “right” way to determine whether a matrix is close to singular is to compute its condition number, which usually is computed via the SVD, though you could in principle also use norm(A)*norm(inv(A)) (except that inv might throw). A nearly as good way is to use rank-revealing QR.

But for a 3 \times 3 matrix, checking that det(A) ≤ rtol * maximum(abs ∘ float, A)^size(A,1) should be fairly reasonable, I guess (without thinking about it too hard), for some decent tolerance, e.g. rtol = sqrt(eps(float(eltype(A)))).

(Determinants can also spuriously overflow/underflow because they can grow exponentially with matrix size. But for 3x3 matrices with reasonable-magnitude entries it shouldn’t be an issue.)

2 Likes