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