Methods · SpecialMatrices.jl

Methods list

SpecialMatrices.Cauchy
SpecialMatrices.Companion
SpecialMatrices.Frobenius
SpecialMatrices.Hilbert
SpecialMatrices.Kahan
SpecialMatrices.Riemann
SpecialMatrices.Strang
SpecialMatrices.Vandermonde
LinearAlgebra.det
SpecialMatrices.dvand!
SpecialMatrices.pvand!
SpecialMatrices.vandtype

Methods usage

SpecialMatrices.Cauchy — Type

Cauchy{T,X,Y} <: AbstractMatrix{T}
Cauchy(x [,y])

Construct lazy Cauchy matrix where

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

Both x and y can be any AbstractArray or NTuple of Numbers, but all elements of x must have the same type; likewise for y. The elements of x and of y should be distinct.

The element type T corresponds to the reciprocal of the sum of pairs of elements of x and y.

julia> Cauchy([2.0 1], (0, 1, 2))
2×3 Cauchy{Float64, Matrix{Float64}, Tuple{Int64, Int64, Int64}}:
 0.5  0.333333  0.25
 1.0  0.5       0.333333

julia> Cauchy(1: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

julia> 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

source

SpecialMatrices.Companion — Type

Companion(v::Union{AbstractVector,Polynomial})::AbstractMatrix

Construct a lazy companion matrix from the coefficients of its characteristic (monic) polynomial.

The matrix is n × n for a vector input of length n or an input Polynomial of degree n, but the storage here is only O(n). This version puts the coefficients along the last column of the matrix. Some texts put the coefficients along the first row, the transpose of the convention used here.

This type has efficient methods for mul! and inv.

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

Also, directly from a polynomial:

julia> using Polynomials

julia> P = Polynomial(2:5)
Polynomial(2 + 3*x + 4*x^2 + 5*x^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

source

SpecialMatrices.Frobenius — Type

[`Frobenius` matrix](https://en.wikipedia.org/wiki/Frobenius_matrix)

Frobenius matrices or Gaussian elimination matrices of the form

[ 1 0 ...     0 ]
[ 0 1 ...     0 ]
[ .........     ]
[ ... 1 ...     ]
[ ... c1 1 ...  ]
[ ... c2 0 1 ...]
[ ............. ]
[ ... ck ...   1]

i.e., an identity matrix with possibly nonzero subdiagonal elements along a single column.

In this implementation, the subdiagonal of the nonzero column is stored as a dense vector, so that the size can be inferred automatically as j+k, where j is the column index and k is the number of subdiagonal elements.

julia> using SpecialMatrices

julia> F = Frobenius(3, 4:6) # 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  4  1  0  0
 0  0  5  0  1  0
 0  0  6  0  0  1

julia> inv(F) # Special form of inverse
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  -5  0  1  0
 0  0  -6  0  0  1

julia> F*F # Special form preserved if the same column has the subdiagonals
6×6 Frobenius{Int64}:
 1  0   0  0  0  0
 0  1   0  0  0  0
 0  0   1  0  0  0
 0  0   8  1  0  0
 0  0  10  0  1  0
 0  0  12  0  0  1

julia> F*Frobenius(2, 2:5) # Promotes to Matrix
6×6 Matrix{Int64}:
 1   0  0  0  0  0
 0   1  0  0  0  0
 0   2  1  0  0  0
 0  11  4  1  0  0
 0  14  5  0  1  0
 0  17  6  0  0  1

julia> F * [10.0,20,30,40,50,60.0]
6-element Vector{Float64}:
  10.0
  20.0
  30.0
 160.0
 200.0
 240.0

source

SpecialMatrices.Hilbert — Type

Hilbert(m [,n])

Construct m × m or m × n Hilbert matrix from its specified dimensions, where element i,j equal to 1 / (i+j-1).

julia> A = 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

julia> Matrix(A)
5×5 Matrix{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:

julia> inv(A)
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

source

SpecialMatrices.Kahan — Type

Kahan matrix

julia> A = 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

julia> A = 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

julia> A = 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

For more details see: N. Higham, "A Survey of Condition Number Estimation for Triangular Matrices," SIMAX, Vol. 29, No. 4, pp. 588, 1987, https://doi.org/10.1137/1029112

source

SpecialMatrices.Riemann — Type

Riemann(N::Int)

Construct N × N Riemann matrix, 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.

julia> Riemann(7)
7×7 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 Friedrich Roesler, "Riemann's hypothesis as an eigenvalue problem," Linear Algebra and its Applications, Vol. 81, p.153-198, Sep. 1986. https://doi.org/10.1016/0024-3795(86)90255-7

source

SpecialMatrices.Strang — Type

Strang([T=Int,] n::Int)

Construct Strang matrix with elements of type T. This special matrix named after Gilbert Strang is symmetric, tridiagonal, and Toeplitz. See Fig. 1 of Julia paper.

julia> Strang(4)
4×4 Strang{Int64}:
  2  -1   0   0
 -1   2  -1   0
  0  -1   2  -1
  0   0  -1   2

julia> Strang(Int16, 3)
3×3 Strang{Int16}:
  2  -1   0
 -1   2  -1
  0  -1   2

source

SpecialMatrices.Vandermonde — Type

Vandermonde(c::AbstractVector)

Create a "lazy" n × n Vandermonde matrix where $A_{ij} = c_i^{j-1}$, requiring only O(n) storage for the vector c.

The transpose and adjoint operations are also lazy.

julia> a = 1:5; A = Vandermonde(a)
5×5 Vandermonde{Int64}:
 1  1   1    1    1
 1  2   4    8   16
 1  3   9   27   81
 1  4  16   64  256
 1  5  25  125  625

julia> A'
5×5 adjoint(::Vandermonde{Int64}) with eltype Int64:
 1   1   1    1    1
 1   2   3    4    5
 1   4   9   16   25
 1   8  27   64  125
 1  16  81  256  625

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), https://doi.org/10.2307/2004623.

julia> A \ a
5-element Vector{Float64}:
 0.0
 1.0
 0.0
 0.0
 0.0

julia> A' \ A[2,:]
5-element Vector{Float64}:
 0.0
 1.0
 0.0
 0.0
 0.0

source

LinearAlgebra.det — Method

det(A::InverseHilbert)

Explicit formula for the determinant. Caution: this function overflows easily.

source

SpecialMatrices.dvand! — Method

dvand!(a, b) -> b

Solve system $A*x = b$ in-place.

A is Vandermonde matrix $A_{ij} = a_i^{j-1}$.

Algorithm by Bjorck & Pereyra, Mathematics of Computation, Vol. 24, No. 112 (1970), pp. 893-903, https://doi.org/10.2307/2004623

source

SpecialMatrices.pvand! — Method

pvand!(a, b) -> b

Solve system $A^T*x = b$ in-place.

$A^T$ is transpose of Vandermonde matrix $A_{ij} = a_i^{j-1}$.

Algorithm by Bjorck & Pereyra, Mathematics of Computation, Vol. 24, No. 112 (1970), pp. 893-903, https://doi.org/10.2307/2004623

source

SpecialMatrices.vandtype — Method

vandtype(T1::Type, T2::Type)

Determine the return type of Vandermonde{T1} \ Vector{T2}.

source