SpecialMatrices

A Julia package for working with special matrix types.

This package extends the LinearAlgebra library with support for special matrices which are used in linear algebra. Every special matrix has its own type. The full matrix is accessed by the command Matrix(A).

Always install the current master:

pkg> add SpecialMatrices#master

Currently supported special matrices

`Cauchy` matrix

Cauchy(x,y)[i,j]=1/(x[i]+y[j])
Cauchy(x)=Cauchy(x,x)
cauchy(k::Number)=Cauchy(collect(1:k))

julia> Cauchy([1,2,3],[3,4,5])
3x3 Cauchy{Int64}:
 0.25      0.2       0.166667
 0.2       0.166667  0.142857
 0.166667  0.142857  0.125

julia> Cauchy([1,2,3])
3x3 Cauchy{Int64}:
 0.5       0.333333  0.25
 0.333333  0.25      0.2
 0.25      0.2       0.166667

julia> Cauchy(pi)
3x3 Cauchy{Float64}:
 0.5       0.333333  0.25
 0.333333  0.25      0.2
 0.25      0.2       0.166667

`Circulant` matrix

julia> Circulant([1,2,3,4])
4x4 Circulant{Int64}:
 1  4  3  2
 2  1  4  3
 3  2  1  4
 4  3  2  1

`Companion` matrix

julia> A=Companion([3,2,1])
3x3 Companion{Int64}:
 0  0  -3
 1  0  -2
 0  1  -1

Also, directly from a polynomial:

julia> using Polynomials

julia> P=Polynomial([2.0,3,4,5])
Polynomial(2 + 3x + 4x^2 + 5x^3)

julia> C=Companion(P)
3×3 Companion{Float64}:
 0.0  0.0  -0.4
 1.0  0.0  -0.6
 0.0  1.0  -0.8

`Frobenius` matrix

Example

julia> using SpecialMatrices

julia> F=Frobenius(3, [1.0,2.0,3.0]) #Specify subdiagonals of column 3
6x6 Frobenius{Float64}:
 1.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  0.0  0.0
 0.0  0.0  2.0  0.0  1.0  0.0
 0.0  0.0  3.0  0.0  0.0  1.0

julia> inv(F) #Special form of inverse
6x6 Frobenius{Float64}:
 1.0  0.0   0.0  0.0  0.0  0.0
 0.0  1.0   0.0  0.0  0.0  0.0
 0.0  0.0   1.0  0.0  0.0  0.0
 0.0  0.0  -1.0  1.0  0.0  0.0
 0.0  0.0  -2.0  0.0  1.0  0.0
 0.0  0.0  -3.0  0.0  0.0  1.0

julia> F*F #Special form preserved if the same column has the subdiagonals
6x6 Frobenius{Float64}:
 1.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  2.0  1.0  0.0  0.0
 0.0  0.0  4.0  0.0  1.0  0.0
 0.0  0.0  6.0  0.0  0.0  1.0

julia> F*Frobenius(2, [5.0,4.0,3.0,2.0]) #Promotes to Matrix
6x6 Array{Float64,2}:
 1.0   0.0  0.0  0.0  0.0  0.0
 0.0   1.0  0.0  0.0  0.0  0.0
 0.0   5.0  1.0  0.0  0.0  0.0
 0.0   9.0  1.0  1.0  0.0  0.0
 0.0  13.0  2.0  0.0  1.0  0.0
 0.0  17.0  3.0  0.0  0.0  1.0

julia> F*[10.0,20,30,40,50,60.0]
6-element Array{Float64,1}:
  10.0
  20.0
  30.0
  70.0
 110.0
 150.0

`Hankel` matrix

Input is a vector of odd length.

julia> Hankel(collect(-4:4))
5x5 Hankel{Int64}:
 -4  -3  -2  -1  0
 -3  -2  -1   0  1
 -2  -1   0   1  2
 -1   0   1   2  3
  0   1   2   3  4

`Hilbert` matrix

julia> A=Hilbert(5)
Hilbert{Rational{Int64}}(5,5)

julia> Matrix(A)
5x5 Array{Rational{Int64},2}:
 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

julia> Matrix(Hilbert(5))
5x5 Array{Rational{Int64},2}:
 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:

julia> inv(A)
5x5 Array{Rational{Int64},2}:
    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

julia> A=Kahan(5,5,1,35)
5x5 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

julia> A=Kahan(5,3,0.5,0)
5x3 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

julia> A=Kahan(3,5,0.5,1e-3)
3x5 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

For more details see N. J. Higham (1987).

`Riemann` matrix

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. Riemann hypothesis holds if and only if det(A) = O( N! N^(-1/2+epsilon)) for every epsilon > 0.

julia> Riemann(7)
7x7 Riemann{Int64}:
  1  -1   1  -1   1  -1   1
 -1   2  -1  -1   2  -1  -1
 -1  -1   3  -1  -1  -1   3
 -1  -1  -1   4  -1  -1  -1
 -1  -1  -1  -1   5  -1  -1
 -1  -1  -1  -1  -1   6  -1
 -1  -1  -1  -1  -1  -1   7

For more details see F. Roesler (1986).

`Strang` matrix

A special SymTridiagonal matrix named after Gilbert Strang

julia> Strang(6)
6x6 Strang{Float64}:
  2.0  -1.0   0.0   0.0   0.0   0.0
 -1.0   2.0  -1.0   0.0   0.0   0.0
  0.0  -1.0   2.0  -1.0   0.0   0.0
  0.0   0.0  -1.0   2.0  -1.0   0.0
  0.0   0.0   0.0  -1.0   2.0  -1.0
  0.0   0.0   0.0   0.0  -1.0   2.0

`Toeplitz` matrix

Input is a vector of odd length.

julia> Toeplitz(collect(-4:4))
5x5 Toeplitz{Int64}:
 0  -1  -2  -3  -4
 1   0  -1  -2  -3
 2   1   0  -1  -2
 3   2   1   0  -1
 4   3   2   1   0

`Vandermonde` matrix

julia> a = collect(1.0:5.0)
julia> A = Vandermonde(a)
5×5 Vandermonde{Float64}:
 1.0  1.0   1.0    1.0    1.0
 1.0  2.0   4.0    8.0   16.0
 1.0  3.0   9.0   27.0   81.0
 1.0  4.0  16.0   64.0  256.0
 1.0  5.0  25.0  125.0  625.0

Adjoint Vandermonde:

julia> A'
5×5 LinearAlgebra.Adjoint{Float64,Vandermonde{Float64}}:
 1.0   1.0   1.0    1.0    1.0
 1.0   2.0   3.0    4.0    5.0
 1.0   4.0   9.0   16.0   25.0
 1.0   8.0  27.0   64.0  125.0
 1.0  16.0  81.0  256.0  625.0

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),

julia> A\a
5-element Array{Float64,1}:
 0.0
 1.0
 0.0
 0.0
 0.0

julia> r2 = A[2,:]
julia> A'\r2
5-element Array{Float64,1}:
 0.0
 1.0
 0.0
 0.0
 0.0