# Need a rarely used binary operator for matrices

Could you suggest a good binary operator that is not commonly used for matrices in Julia Base or packages? I need to use it for some common matrix operations in my code. I don’t want to step on something else.

what’s the operation?

I don’t think there’s much problem with conflicts for everything beyond standard arithmetic. I personally like (\oplus) or '⊗ (\otimes`).

If there is a conflict Julia just informs and asks to make the statement precise, e.g MyLibrary.:⊕

multiplication of a homogeneous matrix with an inhomogeneous vector, e.g.,

hmul(A,x) = A[1:end,1:end-1]*x.+A[:,end]

This is somewhat like a affine (translation) operation.
Some possibilities
\otimesrhrim :
\oplusrhrim :

since they have a sense of directionality to them…

2 Likes

Yes, homogeneous translation would be one application of it. Thanks for the tip!

is the precedence the same as * for those operators?

Seems to be faster to do

hmul2(A,x) = A * [x; 1]

Also, check out StaticArrays for this, it will be much faster, if you are working on e.g. 3-4D vectors.

julia> using StaticArrays

julia> A = @SArray rand(4,4);

julia> x = @SArray rand(3);

julia> @btime hmul($A,$x)
110.808 ns (3 allocations: 352 bytes)
4-element SizedVector{4, Float64, Vector{Float64}} with indices SOneTo(4):
1.281701572863999
1.45284156801398
0.8569634443644434
1.3118854951752263

julia> @btime hmul2($A,$x)
2.900 ns (0 allocations: 0 bytes)
4-element SVector{4, Float64} with indices SOneTo(4):
1.281701572863999
1.45284156801398
0.8569634443644434
1.3118854951752263
3 Likes

Wow, thanks for the tip. Yes, I am using 3-4D vectors and StaticArrays. I’m impressed that the augmentation does not cost an allocation. I see that it is constructed on the stack as an SVector. Very cool!

You don’t normally get any allocations as long as the sizes are statically known. It is easy for the compiler to see that [x; 1] will have type SVector{4, Float64}, though you could also use StaticArrays.push with no difference.

If you want to get rid of the last element of the output vector (which is common), you can use StaticArrays.pop. This can be quite useful, since

julia> typeof(x[1:end-1])
Vector{Float64} (alias for Array{Float64, 1})

julia> typeof(pop(x))
SVector{2, Float64} (alias for SArray{Tuple{2}, Float64, 1, 2})
1 Like

Yes, type stability for these operations is tricky. Thanks for the tip. I have a more general question.

Suppose I have a vector of static vectors, e.g., something like,

100-element reinterpret(SVector{2, Float64}, ::Vector{Float64}

but I’d like to be able to do

y = A*x

What is the Julian and most performant way to accomplish this? E.g., this works, but seems a bit awkward,

y = Ref(A).*x

100-element Vector{SVector{2, Float64}}

Now suppose I want the output y to be a matrix. Tullio is an option, e.g.,

@tullio y[i,k] := A[i,j]*x[k][j]

2×100 Matrix{Float64}

Is there a solution using reinterpret?

This is the usual way.

That being said, if you are worried about the performance of this operation then you may be making a more fundamental mistake — the Matlab/Python style of “vectorized” programming leaves a lot of performance on the table by doing a bunch of separate “vectorized” loops that each perform one small operation (producing a bunch of temporary arrays) rather than combining the calculations into a single pass over the data whenever possible. Presumably, in your real application, you aren’t just doing Ref(A) .* x, and you should try to exploit that by writing a single loop that does all of the operations you need on each x[k] at once.

1 Like

I suppose that loop fusion should eliminate temporaries for element-wise matrix operations. (matrix multiplication clearly cant be fused). Thanks for the response. I find some tension going between vector of vectors representations and matrices (or more generally tensors).

Just using a loop should be a good solution here, but you could also have a look at GitHub - JuliaArrays/HybridArrays.jl: Arrays with both statically and dynamically sized axes in Julia which lets you mix static and dynamic dimensions, e.g. the number of rows is static, while the number of columns is dynamic. I believe it builds on top of StaticArrays.

1 Like
julia> As = @SMatrix rand(4,4);

julia> Xs = [@SVector rand(3) for _ in 1:10];

julia> reinterpret(reshape, Float64, hmul2.(Ref(As), Xs))
4×10 reinterpret(reshape, Float64, ::Vector{SVector{4, Float64}}) with eltype Float64:
1.68023  1.41434   1.33775   1.3449    2.34028  0.655653  2.2069   2.09013  1.45916   1.42397
2.12982  1.87084   1.80543   1.79198   2.79192  1.03538   2.60938  2.55969  1.67285   1.79207
1.11249  1.06255   1.02429   1.0432    1.2869   0.845418  1.24341  1.22446  1.01764   1.02749
1.09398  0.907248  0.901433  0.866234  1.47011  0.479292  1.39198  1.32744  0.951519  0.951254
1 Like

The reference operator precedence table is here:

1 Like