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.

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  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//1    -300//1     1050//1    -1400//1     630//1
  -300//1    4800//1   -18900//1    26880//1  -12600//1
  1050//1  -18900//1    79380//1  -117600//1   56700//1
 -1400//1   26880//1  -117600//1   179200//1  -88200//1
   630//1  -12600//1    56700//1   -88200//1   44100//1

`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

This page was generated using Literate.jl.