Help defining Functor type

General Usage
,
1

I am almost there, but this is my first time using StaticArrays.jl and I am not sure I got the concepts right:

using StaticArrays: SMatrix

"""
    EllipsoidDistance([a₁, a₂, ...], [θ₁, θ₂, ...])

A distance defined by an ellipsoid with given semiaxes `a₁, a₂, ...`
and rotation angles `θ₁, θ₂, ...`.

The positive definite matrix A = PΛP' representing the ellipsoid is
assembled once and cached in the object as an static array. Calls to
(x-y)'*A*(x-y) are therefore very efficient.
"""
immutable EllipsoidDistance{N,T<:Real} <: AbstractDistance
  # state fields
  A::SMatrix{N}

  function EllipsoidDistance{N,T}(as, θs) where {N,T<:Real}
    @assert length(as) == length(θs) == N "number of semiaxes and rotations must match spatial dimension"

    Λ = spdiagm(one(T)./as.^2)
    P = eye(T, N) # TODO: define rotation matrix, probably using Rotations.jl
    A = P*Λ*P'

    new(A) # this is causing the error as expected since the local variable A is not static
  end
end

EllipsoidDistance(as::Vector{T}, θs::Vector{T}) where {T<:Real} = EllipsoidDistance{length(as),T}(as, θs)

(d::EllipsoidDistance)(x, y) = begin
  z = x - y
  z'*A*z
end

The error is caused by the new operator because I am passing a non-static array to it. What would be the appropriate way of caching A?

2

Just convert A to a StaticArray? Also, you probably want to have A concretely typed, aka SArray{Tuple{M, N}, T, 2, M*N} where M is number of rows and N number of columns.

1 Like
3

Thank you @kristoffer.carlsson, I replaced SMatrix{N} by the concrete type SMatrix{N,N,T}, but how to convert A to this type? Also, is this a good thing to do in general? Will the functor be much faster for 2D and 3D vectors for example?

4

Got it, just calling SMatrix{N,N,T}(A) will do it. I had issues with the inferred type T, that is why it was not working.

I am still wondering if the performance gain will be visible, starting some benchmarks.

5

What size is N? It should be significantly faster if N is smallish (assuming x and y are StaticVectors).

6

The dimension will be usually 2D or 3D (at most 4D), but I have no guarantee on the vectors x and y multiplying the matrix, they are general heap allocated arrays. Should I forget about making A static in this case?

7

Variables ending with s are Static (dim = 4):

julia> f(x, A, y) = (z = x-y; z' * A * z)
f (generic function with 1 method)

julia> @btime f($x, $A, $y);
  156.315 ns (2 allocations: 224 bytes)

julia> @btime f($x, $As, $y);
  93.785 ns (2 allocations: 224 bytes)

julia> @btime f($xs, $As, $ys);
  6.540 ns (0 allocations: 0 bytes)
3 Likes
8

I think I will stick with heap-allocated arrays for now. If this becomes a bottleneck at some point, I will hack my way with StaticArrays.

9

As a reference, for me, swapping to static arrays has lead to cleaner code with less bugs (because of array sizes in the type and immutability) as well as tremendous performance improvements. Depending on the performance characteristics of your code, the performance improvements might of course not carry over for you.

1 Like