Merge pull request #26 from ludvigak/vandsolve

Fast solver for Vandermonde systems

README.md

@@ -254,13 +254,45 @@ julia> Toeplitz(collect(-4:4))
 ```
 ## [`Vandermonde`](http://en.wikipedia.org/wiki/Vandermonde_matrix) matrix
 
+```julia
+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
+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 overator `\` 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
-julia> Vandermonde(collect(1:5))
-5x5 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\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
 ```

src/SpecialMatrices.jl

@@ -9,7 +9,8 @@ if VERSION >= v"0.7.0"
 end
 
 import Base: getindex, isassigned, size, *
-
+import Base.\
+
 include("cauchy.jl") #Cauchy matrix
 include("companion.jl") #Companion matrix
 include("frobenius.jl") #Frobenius matrix

src/vandermonde.jl

@@ -1,7 +1,7 @@
 export Vandermonde
 
 struct Vandermonde{T} <: AbstractArray{T, 2}
-	c :: Vector{T}
+    c :: Vector{T}
 end
 
 getindex(V::Vandermonde, i::Int, j::Int) = V.c[i]^(j-1)
@@ -12,10 +12,142 @@ size(V::Vandermonde, r::Int) = (r==1 || r==2) ? length(V.c) :
 size(V::Vandermonde) = length(V.c), length(V.c)
 
 function Matrix(V::Vandermonde{T}) where T
-	n=size(V, 1)
-	M=Array{T}(undef, n, n)
-	for i=1:n
-		M[:,i] = V.c.^(i-1)
-	end
-	M
+    n=size(V, 1)
+    M=Array{T}(undef, n, n)
+    M[:,1] .= 1
+    for j=2:n
+        for i=1:n
+	    M[i,j] = M[i,j-1]*V.c[i]
+        end
+    end
+    M
 end
+
+if VERSION >= v"0.7.0" # Only use Adjoint and Transpose for 0.7 and up
+
+function Matrix(V::Adjoint{T,Vandermonde{T}}) where T
+    n=size(V, 1)
+    M=Array{T}(undef, n, n)
+    M[1,:] .= 1
+    for j=1:n
+        for i=2:n
+	    M[i,j] = M[i-1,j]*adjoint(V.parent.c[j])
+        end
+    end
+    M
+end
+
+function Matrix(V::Transpose{T,Vandermonde{T}}) where T
+    n=size(V, 1)
+    M=Array{T}(undef, n, n)
+    M[1,:] .= 1
+    for j=1:n
+        for i=2:n
+	    M[i,j] = M[i-1,j]*V.parent.c[j]
+        end
+    end
+    M
+end
+
+function \(V::Adjoint{T1,Vandermonde{T1}}, y::AbstractVecOrMat{T2}) where T1 where T2
+    T = vandtype(T1,T2)
+    x = Array{T}(undef, size(y))
+    copyto!(x, y)    
+    pvand!(adjoint(V.parent.c), x)
+    return x
+end    
+
+function \(V::Transpose{T1,Vandermonde{T1}}, y::AbstractVecOrMat{T2}) where T1 where T2
+    T = vandtype(T1,T2)
+    x = Array{T}(undef, size(y))
+    copyto!(x, y)    
+    pvand!(V.parent.c, x)
+    return x
+end
+
+end # End version check
+
+function \(V::Vandermonde{T1}, y::AbstractVecOrMat{T2}) where T1 where T2
+    T = vandtype(T1,T2)
+    x = Array{T}(undef, size(y))
+    copyto!(x, y)
+    dvand!(V.c, x)
+    return x
+end
+
+
+function vandtype(T1::Type, T2::Type)
+    # Figure out the return type of Vandermonde{T1} \ Vector{T2}
+    T = promote_type(T1, T2)
+    S = typeof(oneunit(T)/oneunit(T1))
+    return S
+end
+
+"""
+    pvand!(a, b) -> b
+
+Solves 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
+"""
+function pvand!(alpha, B)
+    n = length(alpha);
+    if n != size(B,1)
+        throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows"))
+    end
+    nrhs = size(B,2)
+    @inbounds begin
+        for j=1:nrhs
+            for k=1:n-1
+                for i=n:-1:k+1
+                    B[i,j] = B[i,j]-alpha[k]*B[i-1,j]
+                end
+            end
+            for k=n-1:-1:1
+                for i=k+1:n
+                    B[i,j] = B[i,j]/(alpha[i]-alpha[i-k])
+                end
+                for i=k:n-1
+                    B[i,j] = B[i,j]-B[i+1,j]
+                end
+            end
+        end
+    end
+end
+"""
+    dvand!(a, b) -> b
+
+Solves 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
+"""
+function dvand!(alpha, B)
+    n = length(alpha)
+    if n != size(B,1)
+        throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows"))
+    end
+    nrhs = size(B,2)
+    @inbounds begin
+        for j=1:nrhs
+            for k=1:n-1
+                for i=n:-1:k+1
+                    B[i,j] = (B[i,j]-B[i-1,j])/(alpha[i]-alpha[i-k])
+                end
+            end
+            for k=n-1:-1:1
+                for i=k:n-1
+                    B[i,j] = B[i,j]-alpha[k]*B[i+1,j]
+                end
+            end
+        end
+    end
+end
+

test/vandermonde.jl

@@ -1,3 +1,52 @@
-V=Vandermonde(collect(1:5))
-@test V[4,5] == 256
+using Compat
+using Compat.Test
+using SpecialMatrices
 
+a = [1,2,3+1im,8,5]
+V=Vandermonde(a)
+
+# Test element
+@test V[4,5] == a[4]^4
+
+# Test full matrix conversions
+M = Matrix(V)
+@test Matrix(V') == M'
+@test Matrix(transpose(V)) == transpose(M)
+
+# Test solving with vector and matrix rhs
+y = [1im,1,5,0,2]
+Y = [y 2*y]
+for rhs=[y, Y]
+    # Test solution
+    @test isapprox(V\rhs, M\rhs)
+    @test isapprox(V'\rhs, M'\rhs)
+    @test isapprox(rhs'/V, rhs'/M)
+    @test isapprox(transpose(V)\rhs, transpose(M)\rhs)
+    # Check that overloading works
+    x = zero(M\rhs)
+    copyto!(x, rhs)
+    SpecialMatrices.dvand!(a, x)
+    @test V\rhs==x
+    if VERSION >= v"0.7.0"
+        copyto!(x, rhs)
+        SpecialMatrices.pvand!(a', x)
+        @test V'\rhs==x
+        @test rhs'/V==x'
+        copyto!(x, rhs)
+        SpecialMatrices.pvand!(a, x)
+        @test transpose(V)\rhs==x
+    end
+end
+
+# Test dimension errors
+rhs = zeros(2,2)
+try
+    V\rhs
+catch e
+    @test typeof(e)==DimensionMismatch
+end
+try
+    V'\rhs
+catch e
+    @test typeof(e)==DimensionMismatch
+end