SpecialMatrices overview

This page illustrates the Julia package SpecialMatrices.

This package extends the LinearAlgebra library with support for special matrices that are used in linear algebra. Every special matrix has its own type and is stored efficiently. Use Matrix(A) to access the full matrix if needed.

ToeplitzMatrices.jl supports Toeplitz, Hankel, and circulant matrices.

This page comes from a single Julia file: 1-overview.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: 1-overview.ipynb, or open it in binder here: 1-overview.ipynb.

Setup

Packages needed here.

using SpecialMatrices
using Polynomials
using LinearAlgebra: factorize, Diagonal

`Cauchy` matrix

Cauchy(1:3, 2:4)

3×3 Cauchy{Rational{Int64}, UnitRange{Int64}, UnitRange{Int64}}:
 1//3  1//4  1//5
 1//4  1//5  1//6
 1//5  1//6  1//7

Cauchy((1., 2.), 1:3)

2×3 Cauchy{Float64, Tuple{Float64, Float64}, UnitRange{Int64}}:
 0.5       0.333333  0.25
 0.333333  0.25      0.2

Cauchy(3)

3×3 Cauchy{Rational{Int64}, UnitRange{Int64}, UnitRange{Int64}}:
 1//2  1//3  1//4
 1//3  1//4  1//5
 1//4  1//5  1//6

`Companion` matrix

Companion(1:3)

3×3 Companion{Int64}:
 0  0  -1
 1  0  -2
 0  1  -3

From a polynomial

p = Polynomial(4:-1:1)

4 + 3∙x + 2∙x² + x³

Companion(p)

3×3 Companion{Float64}:
 0.0  0.0  -4.0
 1.0  0.0  -3.0
 0.0  1.0  -2.0

`Frobenius` matrix

F = Frobenius(3, 2:4) # Specify subdiagonals of column 3

6×6 Frobenius{Int64}:
 1  0  0  0  0  0
 0  1  0  0  0  0
 0  0  1  0  0  0
 0  0  2  1  0  0
 0  0  3  0  1  0
 0  0  4  0  0  1

Special form of inverse:

inv(F)

6×6 Frobenius{Int64}:
 1  0   0  0  0  0
 0  1   0  0  0  0
 0  0   1  0  0  0
 0  0  -2  1  0  0
 0  0  -3  0  1  0
 0  0  -4  0  0  1

Special form preserved if the same column has the subdiagonals

F * F

6×6 Frobenius{Int64}:
 1  0  0  0  0  0
 0  1  0  0  0  0
 0  0  1  0  0  0
 0  0  4  1  0  0
 0  0  6  0  1  0
 0  0  8  0  0  1

Otherwise it promotes to a Matrix:

F * Frobenius(2, 2:5)

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

Efficient matrix-vector multiplication:

F * (1:6)

6-element Vector{Int64}:
  1
  2
  3
 10
 14
 18

`Hilbert` matrix

H = Hilbert(5)

5×5 Hilbert{Rational{Int64}}:
  1    1//2  1//3  1//4  1//5
 1//2  1//3  1//4  1//5  1//6
 1//3  1//4  1//5  1//6  1//7
 1//4  1//5  1//6  1//7  1//8
 1//5  1//6  1//7  1//8  1//9

Inverses are also integer matrices:

inv(H)

5×5 InverseHilbert{Rational{Int64}}:
    25    -300     1050    -1400     630
  -300    4800   -18900    26880  -12600
  1050  -18900    79380  -117600   56700
 -1400   26880  -117600   179200  -88200
   630  -12600    56700   -88200   44100

`Kahan` matrix

See N. J. Higham (1987).

Kahan(5, 5, 1, 35)

5×5 Kahan{Int64, Int64}:
 1.0  -0.540302  -0.540302  -0.540302  -0.540302
 0.0   0.841471  -0.454649  -0.454649  -0.454649
 0.0   0.0        0.708073  -0.382574  -0.382574
 0.0   0.0        0.0        0.595823  -0.321925
 0.0   0.0        0.0        0.0        0.501368

Kahan(5, 3, 0.5, 0)

5×3 Kahan{Float64, Int64}:
 1.0  -0.877583  -0.877583
 0.0   0.479426  -0.420735
 0.0   0.0        0.229849
 0.0   0.0        0.0
 0.0   0.0        0.0

Kahan(3, 5, 0.5, 1e-3)

3×5 Kahan{Float64, Float64}:
 1.0  -0.877583  -0.877583  -0.877583  -0.877583
 0.0   0.479426  -0.420735  -0.420735  -0.420735
 0.0   0.0        0.229849  -0.201711  -0.201711

`Riemann` matrix

A Riemann matrix is defined as A = B[2:N+1, 2:N+1], where B[i,j] = i-1 if i divides j, and -1 otherwise. The Riemann hypothesis holds if and only if det(A) = O( N! N^(-1/2+ϵ)) for every ϵ > 0.

See F. Roesler (1986).

Riemann(5)

5×5 Riemann{Int64}:
  1  -1   1  -1   1
 -1   2  -1  -1   2
 -1  -1   3  -1  -1
 -1  -1  -1   4  -1
 -1  -1  -1  -1   5

`Strang` matrix

A special symmetric, tridiagonal, Toeplitz matrix named after Gilbert Strang.

S = Strang(5)

5×5 Strang{Int64}:
  2  -1   0   0   0
 -1   2  -1   0   0
  0  -1   2  -1   0
  0   0  -1   2  -1
  0   0   0  -1   2

The Strang matrix has a special $L D L'$ factorization:

F = factorize(S)

LinearAlgebra.LDLt{Float64, LinearAlgebra.SymTridiagonal{Float64, Vector{Float64}}}
L factor:
5×5 LinearAlgebra.UnitLowerTriangular{Float64, LinearAlgebra.SymTridiagonal{Float64, Vector{Float64}}}:
  1.0    ⋅          ⋅      ⋅    ⋅ 
 -0.5   1.0         ⋅      ⋅    ⋅ 
  0.0  -0.666667   1.0     ⋅    ⋅ 
  0.0   0.0       -0.75   1.0   ⋅ 
  0.0   0.0        0.0   -0.8  1.0
D factor:
5×5 LinearAlgebra.Diagonal{Float64, Vector{Float64}}:
 2.0   ⋅    ⋅        ⋅     ⋅ 
  ⋅   1.5   ⋅        ⋅     ⋅ 
  ⋅    ⋅   1.33333   ⋅     ⋅ 
  ⋅    ⋅    ⋅       1.25   ⋅ 
  ⋅    ⋅    ⋅        ⋅    1.2

Here is a verification:

F.L * Diagonal(F.D) * F.L'

5×5 Matrix{Float64}:
  2.0  -1.0   0.0   0.0   0.0
 -1.0   2.0  -1.0   0.0   0.0
  0.0  -1.0   2.0  -1.0   0.0
  0.0   0.0  -1.0   2.0  -1.0
  0.0   0.0   0.0  -1.0   2.0

`Vandermonde` matrix

a = 1:4
V = Vandermonde(a)

4×4 Vandermonde{Int64}:
 1  1   1   1
 1  2   4   8
 1  3   9  27
 1  4  16  64

Adjoint Vandermonde:

V'

4×4 adjoint(::Vandermonde{Int64}) with eltype Int64:
 1  1   1   1
 1  2   3   4
 1  4   9  16
 1  8  27  64

The backslash operator \ is overloaded to solve Vandermonde and adjoint Vandermonde systems in $O(n^2)$ time using the algorithm of Björck & Pereyra (1970).

V \ a

4-element Vector{Float64}:
 0.0
 1.0
 0.0
 0.0

V' \ V[3,:]

4-element Vector{Float64}:
 0.0
 0.0
 1.0
 0.0

Reproducibility

This page was generated with the following version of Julia:

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

11-element Vector{SubString{String}}:
 "Julia Version 1.11.1"
 "Commit 8f5b7ca12ad (2024-10-16 10:53 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"
 "  LLVM: libLLVM-16.0.6 (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/SpecialMatrices.jl/SpecialMatrices.jl/docs/Project.toml`
⌅ [e30172f5] Documenter v0.27.25
  [98b081ad] Literate v2.20.1
  [f27b6e38] Polynomials v4.0.11
  [928aab9d] SpecialMatrices v3.1.0 `~/work/SpecialMatrices.jl/SpecialMatrices.jl`
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.