Matrix-vector product

This example illustrates different ways of computing matrix-vector products using the Julia language.

This page comes from a single Julia file: mul-mat-vec.jl.

You can access the source code for such Julia documentation using the 'Edit on GitHub' link in the top right. You can view the corresponding notebook in nbviewer here: mul-mat-vec.ipynb, or open it in binder here: mul-mat-vec.ipynb.

Setup

Add the Julia packages used in this demo. Change false to true in the following code block if you are using any of the following packages for the first time.

if false
    import Pkg
    Pkg.add([
        "BenchmarkTools"
        "InteractiveUtils"
    ])
end

Tell Julia to use the following packages. Run Pkg.add() in the preceding code block first, if needed.

using BenchmarkTools: @benchmark
using InteractiveUtils: versioninfo

Overview of matrix-vector multiplication

The product between a matrix and a compatible vector is such a basic method in linear algebra that, of course, Julia has a function * built-in for it.

In practice one should simply call that method via A * x or possibly *(A, x).

This demo explores other ways of coding the product, to explore techniques for writing efficient code.

We write each method as a function because the most reliable way to benchmark different methods is to use functions.

Built-in *

Conventional high-level matrix-vector multiply function:

function mul0(A::Matrix, x::Vector)
    @boundscheck size(A,2) == length(x) || error("DimensionMismatch(A,x)")
    return A * x
end;

Double loop over m,n

This is the textbook version.

function mul_mn(A::Matrix, x::Vector)
    (M,N) = size(A)
    y = similar(x, M)
    for m in 1:M
        inprod = zero(eltype(x)) # accumulator
        for n in 1:N
            inprod += A[m,n] * x[n]
        end
        y[m] = inprod
    end
    return y
end;

Using @inbounds

function mul_mn_inbounds(A::Matrix, x::Vector)
    (M,N) = size(A)
    @assert N == length(x)
    y = similar(x, M)
    for m in 1:M
        inprod = zero(x[1]) # accumulator
        for n in 1:N
            @inbounds inprod += A[m,n] * x[n]
        end
        @inbounds y[m] = inprod
    end
    return y
end;

Double loop over n,m

We expect this way to be faster because of cache access over m.

function mul_nm(A::Matrix, x::Vector)
    (M,N) = size(A)
    y = zeros(eltype(x), M)
    for n in 1:N
        for m in 1:M
            y[m] += A[m,n] * x[n]
        end
    end
    return y
end;

With @inbounds

function mul_nm_inbounds(A::Matrix, x::Vector)
    (M,N) = size(A)
    @assert N == length(x)
    y = zeros(eltype(x), M)
    for n in 1:N
        for m in 1:M
            @inbounds y[m] += A[m,n] * x[n]
        end
    end
    return y
end;

With @inbounds and @simd

function mul_nm_inbounds_simd(A::Matrix, x::Vector)
    (M,N) = size(A)
    @assert N == length(x)
    y = zeros(eltype(x), M)
    for n in 1:N
        @simd for m in 1:M
            @inbounds y[m] += A[m,n] * x[n]
        end
    end
    return y
end;

And with @views

function mul_nm_inbounds_simd_views(A::Matrix, x::Vector)
    (M,N) = size(A)
    @assert N == length(x)
    y = zeros(eltype(x), M)
    for n in 1:N
        @simd for m in 1:M
            @inbounds @views y[m] += A[m,n] * x[n]
        end
    end
    return y
end;

Row versions

Loop over m.

function mul_row(A::Matrix, x::Vector)
    (M,N) = size(A)
    y = similar(x, M)
    for m in 1:M
        y[m] = transpose(A[m,:]) * x
    end
    return y
end;

with @inbounds

function mul_row_inbounds(A::Matrix, x::Vector)
    (M,N) = size(A)
    @assert N == length(x)
    y = similar(x, M)
    for m in 1:M
        @inbounds y[m] = transpose(A[m,:]) * x
    end
    return y
end;

with @views

function mul_row_views(A::Matrix, x::Vector)
    (M,N) = size(A)
    y = similar(x, M)
    for m in 1:M
        @views y[m] = transpose(A[m,:]) * x
    end
    return y
end;

with both

function mul_row_inbounds_views(A::Matrix, x::Vector)
    (M,N) = size(A)
    @assert N == length(x)
    y = similar(x, M)
    for m in 1:M
        @inbounds @views y[m] = transpose(A[m,:]) * x
    end
    return y
end;

Col versions

Loop over n:

function mul_col(A::Matrix, x::Vector)
    (M,N) = size(A)
    y = zeros(eltype(x), M)
    for n in 1:N
        y += A[:,n] * x[n]
    end
    return y
end;

with broadcast via @. to coalesce operations:

function mul_col_dot(A::Matrix, x::Vector)
    (M,N) = size(A)
    y = zeros(eltype(x), M)
    for n in 1:N
        @. y += A[:,n] * x[n]
    end
    return y
end;

and with @views

function mul_col_dot_views(A::Matrix, x::Vector)
    (M,N) = size(A)
    y = zeros(eltype(x), M)
    for n in 1:N
        @views @. y += A[:,n] * x[n]
    end
    return y
end;

Test and time each version

The results will depend on the computer used, of course.

M = 2^11
N = M - 4 # non-square to stress test
A = randn(Float32, M, N)
x = randn(Float32, N);

flist = (mul0,
    mul_mn, mul_mn_inbounds,
    mul_nm, mul_nm_inbounds, mul_nm_inbounds_simd, mul_nm_inbounds_simd_views,
    mul_row, mul_row_inbounds, mul_row_views, mul_row_inbounds_views,
    mul_col, mul_col_dot, mul_col_dot_views,
);

for f in flist # warm-up and test each version
    @assert A * x ≈ f(A, x)
end;

out = Vector{String}(undef, length(flist))
for (i, f) in enumerate(flist) # benchmark timing for each
    b = @benchmark $f($A,$x)
    tim = round(minimum(b.times)/10^6, digits=1) # in ms
    tim = lpad(tim, 4)
    name = rpad(f, 27)
	alloc = lpad(b.allocs, 5)
	mem = round(b.memory/2^10, digits=1)
    tmp = "$name : $tim ms $alloc alloc $mem KiB"
    out[i] = tmp
    isinteractive() && println(tmp)
end
out

14-element Vector{String}:
 "mul0                        :  0.1 ms     1 alloc 8.1 KiB"
 "mul_mn                      : 27.9 ms     1 alloc 8.1 KiB"
 "mul_mn_inbounds             : 28.5 ms     1 alloc 8.1 KiB"
 "mul_nm                      :  2.0 ms     1 alloc 8.1 KiB"
 "mul_nm_inbounds             :  0.3 ms     1 alloc 8.1 KiB"
 "mul_nm_inbounds_simd        :  0.3 ms     1 alloc 8.1 KiB"
 "mul_nm_inbounds_simd_views  :  0.3 ms     1 alloc 8.1 KiB"
 "mul_row                     : 21.0 ms  2049 alloc 16648.1 KiB"
 "mul_row_inbounds            : 23.5 ms  2049 alloc 16648.1 KiB"
 "mul_row_views               : 28.6 ms     1 alloc 8.1 KiB"
 "mul_row_inbounds_views      : 24.3 ms     1 alloc 8.1 KiB"
 "mul_col                     :  5.6 ms  6133 alloc 49830.6 KiB"
 "mul_col_dot                 :  3.2 ms  2045 alloc 16615.6 KiB"
 "mul_col_dot_views           :  0.3 ms     1 alloc 8.1 KiB"

The following results were for a 2017 iMac with 4.2GHz quad-core Intel Core i7 with macOS Mojave 10.14.6 and Julia 1.6.2.

As expected, simple A*x is the fastest, but one can come quite close to that using proper double loop order with @inbounds or using "dots" and @views to coalesce. Without @views the vector versions have huge memory overhead!

[
"mul0                       :  0.9 ms     1 alloc    16.1 KiB"
"mul_mn                     : 22.5 ms     1 alloc    16.1 KiB"
"mul_mn_inbounds            : 22.0 ms     1 alloc    16.1 KiB"
"mul_nm                     :  3.1 ms     1 alloc    16.1 KiB"
"mul_nm_inbounds            :  1.5 ms     1 alloc    16.1 KiB"
"mul_nm_inbounds_simd       :  1.5 ms     1 alloc    16.1 KiB"
"mul_nm_inbounds_simd_views :  1.5 ms     1 alloc    16.1 KiB"
"mul_row                    : 32.8 ms  2049 alloc 33040.1 KiB"
"mul_row_inbounds           : 32.7 ms  2049 alloc 33040.1 KiB"
"mul_row_views              : 22.4 ms     1 alloc    16.1 KiB"
"mul_row_inbounds_views     : 22.4 ms     1 alloc    16.1 KiB"
"mul_col                    : 16.0 ms  6133 alloc 98894.6 KiB"
"mul_col_dot                :  7.0 ms  2045 alloc 32975.6 KiB"
"mul_col_dot_views          :  1.5 ms     1 alloc    16.1 KiB"
];

Reproducibility

This page was generated with the following version of Julia:

using InteractiveUtils: versioninfo
io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\n')

12-element Vector{SubString{String}}:
 "Julia Version 1.10.1"
 "Commit 7790d6f0641 (2024-02-13 20:41 UTC)"
 "Build Info:"
 "  Official https://julialang.org/ release"
 "Platform Info:"
 "  OS: Linux (x86_64-linux-gnu)"
 "  CPU: 4 × AMD EPYC 7763 64-Core Processor"
 "  WORD_SIZE: 64"
 "  LIBM: libopenlibm"
 "  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)"
 "Threads: 1 default, 0 interactive, 1 GC (on 4 virtual cores)"
 ""

And with the following package versions

import Pkg; Pkg.status()

Status `~/work/book-la-demo/book-la-demo/docs/Project.toml`
  [6e4b80f9] BenchmarkTools v1.5.0
  [aaaa29a8] Clustering v0.15.7
  [35d6a980] ColorSchemes v3.24.0
  [3da002f7] ColorTypes v0.11.4
⌅ [c3611d14] ColorVectorSpace v0.9.10
  [717857b8] DSP v0.7.9
  [72c85766] Demos v0.1.0 `~/work/book-la-demo/book-la-demo`
  [e30172f5] Documenter v1.2.1
  [4f61f5a4] FFTViews v0.3.2
  [7a1cc6ca] FFTW v1.8.0
  [587475ba] Flux v0.14.12
⌅ [a09fc81d] ImageCore v0.9.4
  [71a99df6] ImagePhantoms v0.7.2
  [b964fa9f] LaTeXStrings v1.3.1
  [7031d0ef] LazyGrids v0.5.0
  [599c1a8e] LinearMapsAA v0.11.0
  [98b081ad] Literate v2.16.1
  [7035ae7a] MIRT v0.17.0
  [170b2178] MIRTjim v0.23.0
  [eb30cadb] MLDatasets v0.7.14
  [efe261a4] NFFT v0.13.3
  [6ef6ca0d] NMF v1.0.2
  [15e1cf62] NPZ v0.4.3
  [0b1bfda6] OneHotArrays v0.2.5
  [429524aa] Optim v1.9.2
  [91a5bcdd] Plots v1.40.1
  [f27b6e38] Polynomials v4.0.6
  [2913bbd2] StatsBase v0.34.2
  [d6d074c3] VideoIO v1.0.9
  [b77e0a4c] InteractiveUtils
  [37e2e46d] LinearAlgebra
  [44cfe95a] Pkg v1.10.0
  [9a3f8284] Random
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

This page was generated using Literate.jl.