Julia Tutorial: Diagonal vs diagm

This tutorial discusses the three ways to make diagonal matrices in Julia:
Diagonal and diagm and spdiag

Of these, the Diagonal command is by far the most efficient.
Always use Diagonal instead of diagm for a pure diagonal matrix!

Jeff Fessler, University of Michigan
2019-01-15, Julia 1.0.3
2020-08-05, Julia 1.5.0
2021-08-23, Julia 1.6.2

using LinearAlgebra # this package comes with Julia, no need to install it

The Diagonal function

The Diagonal function in the LinearAlgebra package creates a diagonal "matrix" from a vector of its diagonal elements.

I put "matrix" in quotes because the data type of the output of Diagonal is not an ordinary 2D Array but rather a special Diagonal type that behaves like a matrix but is stored efficiently.

# Example of creating a diagonal matrix
# (internally Julia stores only the nonzero values)
D = Diagonal(10:10:50)

5×5 Diagonal{Int64, StepRange{Int64, Int64}}:
 10   ⋅   ⋅   ⋅   ⋅
  ⋅  20   ⋅   ⋅   ⋅
  ⋅   ⋅  30   ⋅   ⋅
  ⋅   ⋅   ⋅  40   ⋅
  ⋅   ⋅   ⋅   ⋅  50

# We can access individual elements of `D` just like one would with an ordinary matrix.
D[3,3], D[3,1], D[1,2]

(30, 0, 0)

# If we access multiple elements of `D` then we get sparse array
D[1:2,1:4]

2×4 SparseArrays.SparseMatrixCSC{Int64, Int64} with 2 stored entries:
 10   ⋅  ⋅  ⋅
  ⋅  20  ⋅  ⋅

The diagm function

The diagm function in the LinearAlgebra package creates a diagonal matrix from a vector of its diagonal elements, or a banded matrix from multiple band vectors.

Here, matrix means a standard 2D Array in Julia.

M = diagm(0 => 10:10:50)

5×5 Matrix{Int64}:
 10   0   0   0   0
  0  20   0   0   0
  0   0  30   0   0
  0   0   0  40   0
  0   0   0   0  50

The spdiagm function

The spdiagm function in the SparseArrays package creates a diagonal "matrix" from a vector of its diagonal elements, or a banded matrix from multiple band vectors.

Here, "matrix" is quoted again because it is not a standard 2D Array in Julia, but rather a special SparseMatrix type that efficiently stores only the nonzero elements and their array indexes. However, for the special case of a diagonal matrix the sparse data type is not as efficient.

using SparseArrays: spdiagm
S = spdiagm(0 => 10:10:50)

5×5 SparseArrays.SparseMatrixCSC{Int64, Int64} with 5 stored entries:
 10   ⋅   ⋅   ⋅   ⋅
  ⋅  20   ⋅   ⋅   ⋅
  ⋅   ⋅  30   ⋅   ⋅
  ⋅   ⋅   ⋅  40   ⋅
  ⋅   ⋅   ⋅   ⋅  50

Comparing matrix multiplication speeds

Next we compare the execution times of multiplying each of the above 3 diagonal matrix data types by a vector. The differences are dramatic!

using BenchmarkTools
N = 2000 # medium sized problem to illustrate the differences

D = Diagonal(1:N)
M = diagm(0 => 1:N)
S = spdiagm(0 => 1:N)

x = rand(N) # random test vector
y = zeros(N) # space for output

@btime y = D * x; # fastest
@btime y = S * x; # 3× slower
@btime y = M * x; # >100× slower!

  1.299 μs (1 allocation: 15.75 KiB)
  5.027 μs (1 allocation: 15.75 KiB)
  2.745 ms (1 allocation: 15.75 KiB)

The preceding test has to reallocate y everytime, which is inefficient. Here is a faster version that overwrites y each time:

@btime mul!(y, D, x); # fastest
@btime mul!(y, S, x); # 3× slower
@btime mul!(y, M, x); # >100× slower!

  1.175 μs (0 allocations: 0 bytes)
  4.374 μs (0 allocations: 0 bytes)
  2.738 ms (0 allocations: 0 bytes)

# The notation mul!(y, A, x) does not look much like the math y=A*x
# The following package solves this issue
using InplaceOps
@! y = D * x; # typical use
@btime @! $y = D * x; # very fast
@btime @! $y = S * x; # slower
@btime @! $y = M * x; # very slow

  1.186 μs (0 allocations: 0 bytes)
  4.370 μs (0 allocations: 0 bytes)
  2.739 ms (0 allocations: 0 bytes)

So what use is diagm then?

For making banded matrices with multiple diagonals, as follows. NOT for making an ordinary diagonal matrix!

diagm(1 => 1:5, -2 => 3:7)

7×7 Matrix{Int64}:
 0  1  0  0  0  0  0
 0  0  2  0  0  0  0
 3  0  0  3  0  0  0
 0  4  0  0  4  0  0
 0  0  5  0  0  5  0
 0  0  0  6  0  0  0
 0  0  0  0  7  0  0

Rectangular Diagonal matrices

The Diagonal function in the LinearAlgebra package can also create a diagonal matrix from the diagonal elements of a (possibly rectangular) matrix.

# Construct a diagonal matrix out of the elements of a (possibly rectangular) matrix
# When the matrix is rectangular, the size of the output diagonal matrix
# is the minimum of the dimensions of the dimensions of the input matrix
B = [1 2; 3 4; 5 6] # 3 × 2
display(B)
DB = Diagonal(B) # 2 × 2

3×2 Matrix{Int64}:
 1  2
 3  4
 5  6

2×2 Diagonal{Int64, Vector{Int64}}:
 1  ⋅
 ⋅  4

C = [1 2 3; 4 5 6] # 2 × 3
display(C)
DC = Diagonal(C) # 2 × 2

2×3 Matrix{Int64}:
 1  2  3
 4  5  6

2×2 Diagonal{Int64, Vector{Int64}}:
 1  ⋅
 ⋅  5

One must again be careful with the difference between a 1D vector and a 2D array that happens to have one dimension equal to 1, as follows.

Diagonal(ones(3)) # a 1D vector works as expected

3×3 Diagonal{Float64, Vector{Float64}}:
 1.0   ⋅    ⋅ 
  ⋅   1.0   ⋅ 
  ⋅    ⋅   1.0

colmat = ones(3,1) # This is a 2D array of size "3 × 1"
display(colmat)
Diagonal(colmat) # so this is 1 × 1

3×1 Matrix{Float64}:
 1.0
 1.0
 1.0

1×1 Diagonal{Float64, Vector{Float64}}:
 1.0

rowmat = ones(1,3) # This is a 2D array of size 1 × 3
display(rowmat)
Diagonal(rowmat) # this is 1 × 1

1×3 Matrix{Float64}:
 1.0  1.0  1.0

1×1 Diagonal{Float64, Vector{Float64}}:
 1.0