From d97e5124ce06f2f7faaa8608e35f9d5fc46becd2 Mon Sep 17 00:00:00 2001 From: ivanslapnicar Date: Wed, 8 Oct 2014 11:46:42 -0400 Subject: [PATCH 1/3] Added Cauchy, Kahan, Riemann, Companion for polys, and inv for Hilbert. Updated README. --- README.md | 162 ++++++++++++++++++++++++++++++----------- src/SpecialMatrices.jl | 5 +- src/cauchy.jl | 29 ++++++++ src/companion.jl | 14 ++++ src/hilbert.jl | 10 +++ src/kahan.jl | 33 +++++++++ src/riemann.jl | 26 +++++++ 7 files changed, 234 insertions(+), 45 deletions(-) create mode 100644 src/cauchy.jl create mode 100644 src/kahan.jl create mode 100644 src/riemann.jl diff --git a/README.md b/README.md index 5f88be0..6b94d89 100644 --- a/README.md +++ b/README.md @@ -3,13 +3,70 @@ A [Julia](http://julialang.org) package for working with special matrix types. This package extends the `Base.LinAlg` library with support for special -matrices which are used in linear algebra. +matrices which are used in linear algebra. Every special matrix has its own type. +The full matrix is accessed by the command `full(A)`. [![Build Status](https://travis-ci.org/jiahao/SpecialMatrices.jl.svg)](https://travis-ci.org/jiahao/SpecialMatrices.jl) [![Coverage Status](https://img.shields.io/coveralls/jiahao/SpecialMatrices.jl.svg)](https://coveralls.io/r/jiahao/SpecialMatrices.jl) - ## Currently supported special matrices +## [`Cauchy`](http://en.wikipedia.org/wiki/Cauchy_matrix) matrix + +```julia +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 +julia> Circulant([1:4]) +4x4 Circulant{Int64}: + 1 4 3 2 + 2 1 4 3 + 3 2 1 4 + 4 3 2 1 +``` + +## [`Companion`](http://en.wikipedia.org/wiki/Companion_matrix) matrix + +```julia +julia> A=Companion([3,2,1]) +3x3 Companion{Int64}: + 0 0 -3 + 1 0 -2 + 0 1 -1 +``` +Also, directly from a polynomial: + +```julia +julia> using Polynomials + +julia> P=Poly([2,3,4,5]) +Poly(2 + 3x + 4x^2 + 5x^3) + +julia> C=Companion(P) +3x3 Companion{Number}: + 0 0 -0.4 + 1 0 -0.6 + 0 1 -0.8 +``` + ## [`Frobenius`](http://en.wikipedia.org/wiki/Frobenius_matrix) matrix Example @@ -63,16 +120,66 @@ julia> F*[10.0:10.0:60.0] 150.0 ``` -## [`Companion`](http://en.wikipedia.org/wiki/Companion_matrix) matrix +## [`Hankel`](http://en.wikipedia.org/wiki/Hankel_matrix) matrix + +Input is a vector of odd length. ```julia -julia> A=Companion([1,2,1]) -3x3 Companion{Int64}: - 0 0 -1 - 1 0 -2 - 0 1 -1 +julia> Hankel([-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`](http://en.wikipedia.org/wiki/Hilbert_matrix) matrix + +```julia +julia> full(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 ``` +## [`Kahan`](http://math.nist.gov/MatrixMarket/data/MMDELI/kahan/kahan.html) matrix + +```julia +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][Higham87]. + +[Higham87]:http://eprints.ma.man.ac.uk/695/01/covered/MIMS_ep2007_10.pdf"N. Higham, A +Survey of Condition Number Estimation for Triangular +Matrices, SIMAX, Vol. 29, No. 4, pp. 588, 1987" + + + ## `Strang` matrix A special `SymTridiagonal` matrix named after Gilbert Strang @@ -88,19 +195,9 @@ julia> Strang(6) 0.0 0.0 0.0 0.0 -1.0 2.0 ``` -## `Hankel` matrix - -```julia -julia> Hankel([-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 -``` +## [`Toeplitz`](http://en.wikipedia.org/wiki/Toeplitz_matrix) matrix -## `Toeplitz` matrix +Input is a vector of odd length. ```julia julia> Toeplitz([-4:4]) @@ -111,31 +208,8 @@ julia> Toeplitz([-4:4]) 3 2 1 0 -1 4 3 2 1 0 ``` +## [`Vandermonde`](http://en.wikipedia.org/wiki/Vandermonde_matrix) matrix -## `Circulant` matrix - -```julia -julia> Circulant([1:4]) -4x4 Circulant{Int64}: - 1 4 3 2 - 2 1 4 3 - 3 2 1 4 - 4 3 2 1 -``` - -## `Hilbert` matrix - -```julia -julia> full(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 -``` - -## `Vandermonde` matrix ```julia julia> Vandermonde([1:5]) diff --git a/src/SpecialMatrices.jl b/src/SpecialMatrices.jl index fb9ef1c..74fd6e6 100644 --- a/src/SpecialMatrices.jl +++ b/src/SpecialMatrices.jl @@ -2,11 +2,14 @@ module SpecialMatrices import Base: A_mul_B!, full, getindex, inv, isassigned, size, * +include("cauchy.jl") #Cauchy matrix include("companion.jl") #Companion matrix include("frobenius.jl") #Frobenius matrix -include("strang.jl") #Strang matrix include("hankel.jl") #Hankel matrix include("hilbert.jl") #Hilbert matrix +include("kahan.jl") #Kahan matrix +include("riemann.jl") #Riemann matrix +include("strang.jl") #Strang matrix include("toeplitz.jl") #Toeplitz matrix include("vandermonde.jl") #Vandermonde matrix diff --git a/src/cauchy.jl b/src/cauchy.jl new file mode 100644 index 0000000..32fa375 --- /dev/null +++ b/src/cauchy.jl @@ -0,0 +1,29 @@ +#### Cauchy matrix + +export Cauchy + +immutable Cauchy{T} <: AbstractMatrix{T} + x::Vector{T} # + y::Vector{T} # +end # immutable + +function Cauchy(k::Number) + Cauchy([1:k],[1:k]) +end + +function Cauchy(x::Vector) + Cauchy(x,x) +end + +# Define its size + +size(A::Cauchy, dim::Integer) = length(A.x) +size(A::Cauchy)= size(A,1), size(A,1) + +# Index into a Cauchy +function getindex(A::Cauchy,i::Integer,j::Integer) + return 1.0/(A.x[i]+A.y[j]) +end # getindex + +# Dense version of Cauchy +full(A::Cauchy) =[A[i,j] for i=1:size(A,1), j=1:size(A,2)] diff --git a/src/companion.jl b/src/companion.jl index c47f369..866b23c 100644 --- a/src/companion.jl +++ b/src/companion.jl @@ -3,6 +3,20 @@ export Companion immutable Companion{T} <: AbstractArray{T, 2} c :: Vector{T} end +#From polynomial + +using Polynomials +#Generate companion matrix from a polynomial + +function Companion(P::Poly) + n = length(P) + c = Array(Number,n-1) + d=P.a[n] + for i=1:n-1 + c[i]=P.a[i]/d + end + Companion(c) +end #Basic property computations size(C::Companion, r::Int) = (r==1 || r==2) ? length(C.c) : diff --git a/src/hilbert.jl b/src/hilbert.jl index fbbef00..7b45323 100644 --- a/src/hilbert.jl +++ b/src/hilbert.jl @@ -18,3 +18,13 @@ function full{T}(H::Hilbert{T}) end Hf end + +function inv{T}(H::Hilbert{T}) + invH=zeros(T,H.m,H.m) + for i=1:H.m, j=1:H.m + invH[i,j]=(-1)^(i+j)*(i+j-1)*binomial(H.m+i-1,H.m-j)* + binomial(H.m+j-1,H.m-i)*binomial(i+j-2,i-1)^2 + end + invH +end + diff --git a/src/kahan.jl b/src/kahan.jl new file mode 100644 index 0000000..680d6c6 --- /dev/null +++ b/src/kahan.jl @@ -0,0 +1,33 @@ +#### Kahan matrix + +export Kahan + +immutable Kahan{T<:Number,T1<:Number} <: AbstractMatrix{T} + m::Int # dimension + n::Int # dimension + theta::T # angle + pert::T1 # perturbation is pert*eps() +end # immutable + +# Define its size +size(A::Kahan, r::Int) = r==1 ? A.m : A.n +size(A::Kahan) = A.m, A.n + +# Index into a Kahan +function getindex(A::Kahan,i::Integer,j::Integer) + m=minimum(size(A)) + t=tan(A.theta) + c=1.0/sqrt(1.0+t^2) + s=t*c + if i>m; return 0.0 + elseif i>j; return 0.0 + elseif i==j; return s^(i-1)+A.pert*eps()*(m-i+1) + else return -c*s^(i-1) + end +end # getindex + +# Dense version of Kahan +full(A::Kahan) =[A[i,j] for i=1:size(A,1), j=1:size(A,2)] + + + diff --git a/src/riemann.jl b/src/riemann.jl new file mode 100644 index 0000000..43c4d99 --- /dev/null +++ b/src/riemann.jl @@ -0,0 +1,26 @@ +#### Riemann Matrix +import Base: size +export Riemann + +immutable Riemann{Int} <: AbstractMatrix{Int} + n::Int +end # immutable + +# Define its size + +size(A::Riemann, dim::Integer) = A.n +size(A::Riemann)= size(A,1), size(A,1) + +# Index into a Riemann +function getindex(A::Riemann,i::Integer,j::Integer) +# return (i+1)%(j+1)==0 ? i : -1 + if i<=j && (j+1)%(i+1)==0 + return i + else + return -1 + end +end # getindex + +# Dense version of Riemann +full(A::Riemann) =[A[i,j] for i=1:size(A,1), j=1:size(A,2)] + From ae24838e49b40b3d5c1c3c5ff8902f7dde8484c7 Mon Sep 17 00:00:00 2001 From: ivanslapnicar Date: Wed, 8 Oct 2014 12:23:04 -0400 Subject: [PATCH 2/3] update README.md --- README.md | 54 +++++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 49 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index 6b94d89..3a252a3 100644 --- a/README.md +++ b/README.md @@ -137,6 +137,17 @@ julia> Hankel([-4:4]) ## [`Hilbert`](http://en.wikipedia.org/wiki/Hilbert_matrix) matrix ```julia +julia> A=Hilbert(5) +Hilbert{Rational{Int64}}(5,5) + +julia> full(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> full(Hilbert(5)) 5x5 Array{Rational{Int64},2}: 1//1 1//2 1//3 1//4 1//5 @@ -145,6 +156,17 @@ julia> full(Hilbert(5)) 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 +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`](http://math.nist.gov/MatrixMarket/data/MMDELI/kahan/kahan.html) matrix @@ -162,7 +184,7 @@ julia> A=Kahan(5,3,0.5,0) 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 0.0 0.0 0.0 julia> A=Kahan(3,5,0.5,1e-3) @@ -172,11 +194,33 @@ julia> A=Kahan(3,5,0.5,1e-3) 0.0 0.0 0.229849 -0.201711 -0.201711 ``` -For more details see [N. J. Higham, 1987][Higham87]. +For more details see [N. J. Higham (1987)][Higham87]. + +[Higham87]: http://eprints.ma.man.ac.uk/695/01/covered/MIMS_ep2007_10.pdf "N. Higham, A Survey of Condition Number Estimation for Triangular Matrices, SIMAX, Vol. 29, No. 4, pp. 588, 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](http://en.wikipedia.org/wiki/Riemann_hypothesis) holds +if and only if `det(A) = O( N! N^(-1/2+epsilon))` for every `epsilon > 0`. + +```julia +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)][Roesler1986]. + +[Roesler1986]: http://www.sciencedirect.com/science/article/pii/0024379586902557 "Friedrich Roesler, Riemann's hypothesis as an eigenvalue problem, Linear Algebra and its Applications, Vol. 81, (1986)" -[Higham87]:http://eprints.ma.man.ac.uk/695/01/covered/MIMS_ep2007_10.pdf"N. Higham, A -Survey of Condition Number Estimation for Triangular -Matrices, SIMAX, Vol. 29, No. 4, pp. 588, 1987" From 84233bcd7add3520af6690ef2c62533079cdd717 Mon Sep 17 00:00:00 2001 From: ivanslapnicar Date: Wed, 8 Oct 2014 14:00:18 -0400 Subject: [PATCH 3/3] Updated REQUIRE --- REQUIRE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/REQUIRE b/REQUIRE index 1d1970f..8479dbe 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,2 +1,2 @@ julia 0.2- - +Polynomials 0.0.3 \ No newline at end of file