Normalizing the columns of an SMatrix

General Usage Performance
, ,
1

I have an SMatrix of dimension (D, k). What I want to do is to normalize each column of this matrix, as if they were individual vectors.

What I do at the moment is always carry around with me a Vector[SVector], and do:

using StaticArrays
D = 3; k = 3
dummy = [rand(SVector{D}) for i in 1:k]
A = rand(SMatrix{D, k})
for i in 1:k 
  dummy[i] = normalize(A[:, i])
end
A = hcat(dummy...)

But that doesn’t seem very efficient…

Both D, k can be known at compile time, maybe there is a more performant way involving metaprogramming?

2

This should do the job

function normalize_impl(::Type{SMatrix{D, K, T, DK}}) where {D, K, T, DK}
    exprs = []
    for j = 1:K
        c_j = Symbol("c", j)
        push!(exprs, :($c_j = normalize(A[:, $j])))
    end

    ops = Expr[]
    for j=1:K, i=1:D
        c_j = Symbol("c", j)
        push!(ops, :($c_j[$i]))
    end


    Expr(:block,
        exprs...,
        Expr(:call, SMatrix{D, K, T, DK}, ops...)
        )
end

@generated function normalize_colums(A::SMatrix)
    normalize_impl(A)
end

You can check the generated code with

A = @SMatrix rand(4, 3)
normalize_impl(typeof(A))

which would return

begin 
    c1 = normalize(A[:, 1])
    c2 = normalize(A[:, 2])
    c3 = normalize(A[:, 3])
    (StaticArrays.SArray{Tuple{4,3},Float64,2,12})(c1[1], c1[2], c1[3], c1[4], c2[1], c2[2], c2[3], c2[4], c3[1], c3[2], c3[3], c3[4])
end
1 Like
3

You could do something like:


julia> v = @SVector [1/norm(A[:, i]) for i in 1:k];

julia> A2 = A .* v';

julia> [norm(A2[:, i]) for i in 1:k]
3-element Array{Float64,1}:
 1.0
 1.0
 1.0

Although I’m not sure exactly how the comprehensions are implemented in StaticArrays. If the comprehensions are bad you could easily write a generated function that does the equivalent.

4

Thank you so much @saschatimme and @kristoffer.carlsson ! Your help really helped me a lot!

Now I am using a modified version from Sascha, because I want to realease a new tag soon. But Kristoffer’s suggestion is definitely a good exercise for me to improve my metaprogramming skills (which are pitiful) and I will try it out at a later point!

5

This seems to be pretty efficient, though I haven’t compared the speed with the other suggestions

a = @SMatrix rand(5,4)
b = a ./ sqrt.(sum(abs2, a, Val{1}))