Suboptimal matrix vector product

There appears to be a scaling problem with matrix vector products, I think to do with vectorization (from inspecting generated code).

Example:

julia> let N=32, M=8
           @btime mul!($(zeros(Float32, M)), $(rand(Float32, M, N)), $(rand(Float32, N)))
           @btime $(rand(SMatrix{M, N, Float32})) * $(rand(SVector{N, Float32}))
           @btime $(rand(SMatrix{M, N, Float32})) * $(rand(SMatrix{N, 1, Float32}))
       end;
  33.646 ns (0 allocations: 0 bytes)
  8.601 ns (0 allocations: 0 bytes)
  9.834 ns (0 allocations: 0 bytes)

julia> let N=64, M=8
           @btime mul!($(zeros(Float32, M)), $(rand(Float32, M, N)), $(rand(Float32, N)))
           @btime $(rand(SMatrix{M, N, Float32})) * $(rand(SVector{N, Float32}))
           @btime $(rand(SMatrix{M, N, Float32})) * $(rand(SMatrix{N, 1, Float32}))
       end;
  44.674 ns (0 allocations: 0 bytes)
  40.647 ns (0 allocations: 0 bytes)
  16.021 ns (0 allocations: 0 bytes)

julia> let N=128, M=8
           @btime mul!($(zeros(Float32, M)), $(rand(Float32, M, N)), $(rand(Float32, N)))
           @btime $(rand(SMatrix{M, N, Float32})) * $(rand(SVector{N, Float32}))
           @btime $(rand(SMatrix{M, N, Float32})) * $(rand(SMatrix{N, 1, Float32}))
       end;
  65.031 ns (0 allocations: 0 bytes)
  351.633 ns (0 allocations: 0 bytes)
  46.571 ns (0 allocations: 0 bytes)

The equivalent matrix matrix product scales fine, so my guess is the necessary functionality is there, but it needs to be wired up?

Good catch. The logic here:

StaticArrays.jl/src/matrix_multiply.jl

Lines 94 to 102 in f02280c

    
           if sa[2] != 0 
        
               retexpr = gen_by_access(wrapped_a) do access_a 
        
                   exprs = mul_smat_vec_exprs(sa, access_a) 
        
                   return :(@inbounds return similar_type(b, T, Size(sa[1]))(tuple($(exprs...)))) 
        
               end 
        
           else 
        
               exprs = [:(zero(T)) for k = 1:sa[1]] 
        
               retexpr = :(@inbounds return similar_type(b, T, Size(sa[1]))(tuple($(exprs...)))) 
        
           end

is not special-cased for such large arrays in the same way the matrix-matrix product is:

StaticArrays.jl/src/matrix_multiply.jl

Lines 135 to 158 in f02280c

    
           if a <: StaticMatrix && b <: StaticMatrix 
        
               # Julia unrolls these loops pretty well 
        
               return quote 
        
                   @_inline_meta 
        
                   return mul_loop(Sa, Sb, a, b) 
        
               end 
        
           elseif sa[1]*sa[2]*sb[2] <= 4*8*8*8*a_tri_mul*b_tri_mul*ab_tri_mul || a <: Diagonal || b <: Diagonal 
        
               return quote 
        
                   @_inline_meta 
        
                   return mul_unrolled(Sa, Sb, a, b) 
        
               end 
        
           elseif (sa[1] <= 14 && sa[2] <= 14 && sb[2] <= 14) || !(a <: StaticMatrix) || !(b <: StaticMatrix) 
        
               return quote 
        
                   @_inline_meta 
        
                   return mul_unrolled_chunks(Sa, Sb, a, b) 
        
               end 
        
           else 
        
               # we don't have any special code for handling this case so let's fall back to 
        
               # the generic implementation of matrix multiplication 
        
               return quote 
        
                   @_inline_meta 
        
                   return mul_generic(Sa, Sb, a, b) 
        
               end 
        
           end

. I would be happy to review a PR improving it 🙂

I'd take me a bit of time to get my head around that code now, but wouldn't it be possible to run everything through the matrix matrix path with reshape in and out? It should be no-op for SVector -> SMatrix -> SVector, it's the same memory.

That code is generic for other static vector types where reshaping adds a wrapper. Maybe Julia could optimize it out well enough but I'm not sure. It's certainly not an obviously simpler solution than optimizing the matrix-vector code.

Navigation Menu

Repository navigation

Suboptimal matrix vector product #1345

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

	if sa[2] != 0
	retexpr = gen_by_access(wrapped_a) do access_a
	exprs = mul_smat_vec_exprs(sa, access_a)
	return :(@inbounds return similar_type(b, T, Size(sa[1]))(tuple($(exprs...))))
	end
	else
	exprs = [:(zero(T)) for k = 1:sa[1]]
	retexpr = :(@inbounds return similar_type(b, T, Size(sa[1]))(tuple($(exprs...))))
	end

	if a <: StaticMatrix && b <: StaticMatrix
	# Julia unrolls these loops pretty well
	return quote
	@_inline_meta
	return mul_loop(Sa, Sb, a, b)
	end
	elseif sa[1]sa[2]sb[2] <= 4888a_tri_mulb_tri_mulab_tri_mul \|\| a <: Diagonal \|\| b <: Diagonal
	return quote
	@_inline_meta
	return mul_unrolled(Sa, Sb, a, b)
	end
	elseif (sa[1] <= 14 && sa[2] <= 14 && sb[2] <= 14) \|\| !(a <: StaticMatrix) \|\| !(b <: StaticMatrix)
	return quote
	@_inline_meta
	return mul_unrolled_chunks(Sa, Sb, a, b)
	end
	else
	# we don't have any special code for handling this case so let's fall back to
	# the generic implementation of matrix multiplication
	return quote
	@_inline_meta
	return mul_generic(Sa, Sb, a, b)
	end
	end

Suboptimal matrix vector product #1345

Description

Activity

mateuszbaran commented on Jun 10, 2026

ancapdev commented on Jun 10, 2026

mateuszbaran commented on Jun 10, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions