From 7ae8adfa095e69d97a4ec24e863bc619e868185d Mon Sep 17 00:00:00 2001 From: Andreas Noack Jensen Date: Thu, 12 Feb 2015 16:54:28 -0500 Subject: [PATCH] Add qrfact(SparseMatrixCSC) by wrapping SPQR. Make \(SparseMatrixCSC) work for least squares problems. Some tests. --- base/linalg/generic.jl | 8 ++- base/sparse.jl | 1 + base/sparse/cholmod.jl | 12 +++- base/sparse/linalg.jl | 4 +- base/sparse/spqr.jl | 133 +++++++++++++++++++++++++++++++++++++++++ base/sparse/spqr_h.jl | 30 ---------- test/runtests.jl | 8 +-- test/sparse/spqr.jl | 43 +++++++++++++ test/sparsetests.jl | 4 ++ 9 files changed, 201 insertions(+), 42 deletions(-) create mode 100644 base/sparse/spqr.jl delete mode 100644 base/sparse/spqr_h.jl create mode 100644 test/sparse/spqr.jl create mode 100644 test/sparsetests.jl diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl index cd9a04e729fd7..3fbeedf7666bd 100644 --- a/base/linalg/generic.jl +++ b/base/linalg/generic.jl @@ -237,11 +237,15 @@ function inv{T}(A::AbstractMatrix{T}) A_ldiv_B!(factorize(convert(AbstractMatrix{S}, A)), eye(S, chksquare(A))) end +function \{T}(A::AbstractMatrix{T}, B::AbstractVecOrMat{T}) + size(A,1) == size(B,1) || throw(DimensionMismatch("LHS and RHS should have the same number of rows. LHS has $(size(A,1)) rows, but RHS has $(size(B,1)) rows.")) + factorize(A)\B +end function \{TA,TB}(A::AbstractMatrix{TA}, B::AbstractVecOrMat{TB}) TC = typeof(one(TA)/one(TB)) - size(A,1) == size(B,1) || throw(DimensionMismatch("LHS and RHS should have the same number of rows. LHS has $(size(A,1)) rows, but RHS has $(size(B,1)) rows.")) - \(factorize(TA == TC ? A : convert(AbstractMatrix{TC}, A)), TB == TC ? copy(B) : convert(AbstractArray{TC}, B)) + convert(AbstractMatrix{TC}, A)\convert(AbstractArray{TC}, B) end + \(a::AbstractVector, b::AbstractArray) = reshape(a, length(a), 1) \ b /(A::AbstractVecOrMat, B::AbstractVecOrMat) = (B' \ A')' # \(A::StridedMatrix,x::Number) = inv(A)*x Should be added at some point when the old elementwise version has been deprecated long enough diff --git a/base/sparse.jl b/base/sparse.jl index c7eaee447de11..51a5258b393ef 100644 --- a/base/sparse.jl +++ b/base/sparse.jl @@ -19,5 +19,6 @@ include("sparse/csparse.jl") include("sparse/linalg.jl") include("sparse/umfpack.jl") include("sparse/cholmod.jl") +include("sparse/spqr.jl") end # module SparseMatrix diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index 12908b4f6beb5..fc612e5da06af 100644 --- a/base/sparse/cholmod.jl +++ b/base/sparse/cholmod.jl @@ -13,7 +13,7 @@ export Factor, Sparse -using Base.SparseMatrix: AbstractSparseMatrix, SparseMatrixCSC, increment, increment!, indtype, decrement, decrement! +using Base.SparseMatrix: AbstractSparseMatrix, SparseMatrixCSC, increment, indtype ######### # Setup # @@ -716,7 +716,7 @@ function Sparse(filename::ASCIIString) end ## convertion back to base Julia types -function convert{T<:VTypes}(::Type{Matrix}, D::Dense{T}) +function convert{T}(::Type{Matrix{T}}, D::Dense{T}) s = unsafe_load(D.p) a = Array(T, s.nrow, s.ncol) if s.d == s.nrow @@ -730,6 +730,14 @@ function convert{T<:VTypes}(::Type{Matrix}, D::Dense{T}) end a end +convert{T}(::Type{Matrix}, D::Dense{T}) = convert(Matrix{T}, D) +function convert{T}(::Type{Vector{T}}, D::Dense{T}) + if size(D, 2) > 1 + throw(DimensionMismatch("input must be a vector but had $(size(D, 2)) columns")) + end + reshape(convert(Matrix, D), size(D, 1)) +end +convert{T}(::Type{Vector}, D::Dense{T}) = convert(Vector{T}, D) function convert{Tv,Ti}(::Type{SparseMatrixCSC{Tv,Ti}}, A::Sparse{Tv,Ti}) s = unsafe_load(A.p) diff --git a/base/sparse/linalg.jl b/base/sparse/linalg.jl index 3b5249c313699..746820a303833 100644 --- a/base/sparse/linalg.jl +++ b/base/sparse/linalg.jl @@ -686,8 +686,10 @@ function factorize(A::SparseMatrixCSC) end end return lufact(A) + elseif m > n + return qrfact(A) else - throw(ArgumentError("sparse least squares problems by QR are not handled yet")) + throw(ArgumentError("underdetermined systemes are not implemented yet")) end end diff --git a/base/sparse/spqr.jl b/base/sparse/spqr.jl new file mode 100644 index 0000000000000..83041657447f4 --- /dev/null +++ b/base/sparse/spqr.jl @@ -0,0 +1,133 @@ +module SPQR + +# ordering options */ +const ORDERING_FIXED = int32(0) +const ORDERING_NATURAL = int32(1) +const ORDERING_COLAMD = int32(2) +const ORDERING_GIVEN = int32(3) # only used for C/C++ interface +const ORDERING_CHOLMOD = int32(4) # CHOLMOD best-effort (COLAMD, METIS,...) +const ORDERING_AMD = int32(5) # AMD(A'*A) +const ORDERING_METIS = int32(6) # metis(A'*A) +const ORDERING_DEFAULT = int32(7) # SuiteSparseQR default ordering +const ORDERING_BEST = int32(8) # try COLAMD, AMD, and METIS; pick best +const ORDERING_BESTAMD = int32(9) # try COLAMD and AMD; pick best# + +# Let [m n] = size of the matrix after pruning singletons. The default +# ordering strategy is to use COLAMD if m <= 2*n. Otherwise, AMD(A'A) is +# tried. If there is a high fill-in with AMD then try METIS(A'A) and take +# the best of AMD and METIS. METIS is not tried if it isn't installed. + +# tol options +const DEFAULT_TOL = int32(-2) # if tol <= -2, the default tol is used +const NO_TOL = int32(-1) # if -2 < tol < 0, then no tol is used + +# for qmult, method can be 0,1,2,3: +const QTX = int32(0) +const QX = int32(1) +const XQT = int32(2) +const XQ = int32(3) + +# system can be 0,1,2,3: Given Q*R=A*E from SuiteSparseQR_factorize: +const RX_EQUALS_B = int32(0) # solve R*X=B or X = R\B +const RETX_EQUALS_B = int32(1) # solve R*E'*X=B or X = E*(R\B) +const RTX_EQUALS_B = int32(2) # solve R'*X=B or X = R'\B +const RTX_EQUALS_ETB = int32(3) # solve R'*X=E'*B or X = R'\(E'*B) + + +using Base.SparseMatrix: SparseMatrixCSC +using Base.SparseMatrix.CHOLMOD: C_Dense, C_Sparse, Dense, ITypes, Sparse, VTypes, common + +import Base: size +import Base.LinAlg: qrfact +import Base.SparseMatrix.CHOLMOD: convert, free! + + + +immutable C_Factorization{Tv<:VTypes,Ti<:ITypes} + xtype::Cint + factors::Ptr{Tv} +end + +type Factorization{Tv<:VTypes,Ti<:ITypes} <: Base.LinAlg.Factorization{Tv} + m::Int + n::Int + p::Ptr{C_Factorization{Tv,Ti}} +end + +size(F::Factorization) = (F.m, F.n) +function size(F::Factorization, i::Integer) + if i < 1 + throw(ArgumentError("dimension must be positive")) + end + if i == 1 + return F.m + elseif i == 2 + return F.n + end + return 1 +end + +function free!{Tv<:VTypes,Ti<:ITypes}(F::Factorization{Tv,Ti}) + ccall((:SuiteSparseQR_C_free, :libspqr), Cint, + (Ptr{Ptr{C_Factorization{Tv,Ti}}}, Ptr{Void}), + &F.p, common(Ti)) == 1 +end + +function backslash{Tv<:VTypes,Ti<:ITypes}(ordering::Integer, tol::Real, A::Sparse{Tv,Ti}, B::Dense{Tv}) + m, n = size(A) + if m != size(B, 1) + throw(DimensionMismatch("left hand side and right hand side must have same number of rows")) + end + d = Dense(ccall((:SuiteSparseQR_C_backslash, :libspqr), Ptr{C_Dense{Tv}}, + (Cint, Cdouble, Ptr{C_Sparse{Tv,Ti}}, Ptr{C_Dense{Tv}}, Ptr{Void}), + ordering, tol, A.p, B.p, common(Ti))) + finalizer(d, free!) + d +end + +function factorize{Tv<:VTypes,Ti<:ITypes}(ordering::Integer, tol::Real, A::Sparse{Tv,Ti}) + f = Factorization(size(A)..., ccall((:SuiteSparseQR_C_factorize, :libspqr), Ptr{C_Factorization{Tv,Ti}}, + (Cint, Cdouble, Ptr{Sparse{Tv,Ti}}, Ptr{Void}), + ordering, tol, A.p, common(Ti))) + finalizer(f, free!) + f +end + +function solve{Tv<:VTypes,Ti<:ITypes}(system::Integer, QR::Factorization{Tv,Ti}, B::Dense{Tv}) + m, n = size(QR) + mB = size(B, 1) + if (system == RX_EQUALS_B || system == RETX_EQUALS_B) && m != mB + throw(DimensionMismatch("number of rows in factorized matrix must equal number of rows in right hand side")) + elseif (system == RTX_EQUALS_ETB || system == RTX_EQUALS_B) && n != mB + throw(DimensionMismatch("number of columns in factorized matrix must equal number of rows in right hand side")) + end + d = Dense(ccall((:SuiteSparseQR_C_solve, :libspqr), Ptr{C_Dense{Tv}}, + (Cint, Ptr{C_Factorization{Tv,Ti}}, Ptr{C_Dense{Tv}}, Ptr{Void}), + system, QR.p, B.p, common(Ti))) + finalizer(d, free!) + d +end + +function qmult{Tv<:VTypes,Ti<:ITypes}(method::Integer, QR::Factorization{Tv,Ti}, X::Dense{Tv}) + mQR, nQR = size(QR) + mX, nX = size(X) + if (method == QTX || method == QX) && mQR != mX + throw(DimensionMismatch("Q matrix size $mQR and dense matrix has $mX rows")) + elseif (method == XQT || method == XQ) && mQR != nX + throw(DimensionMismatch("Q matrix size $mQR and dense matrix has $nX columns")) + end + d = Dense(ccall((:SuiteSparseQR_C_qmult, :libspqr), Ptr{C_Dense{Tv}}, + (Cint, Ptr{C_Factorization{Tv,Ti}}, Ptr{C_Dense{Tv}}, Ptr{Void}), + method, QR.p, X.p, common(Ti))) + finalizer(d, free!) + d +end + +qrfact(A::SparseMatrixCSC) = factorize(ORDERING_DEFAULT, DEFAULT_TOL, Sparse(A)) + +function (\){T}(F::Factorization{T}, B::StridedVecOrMat{T}) + QtB = qmult(QTX, F, Dense(B)) + convert(typeof(B), solve(RETX_EQUALS_B, F, QtB)) +end + +end # module \ No newline at end of file diff --git a/base/sparse/spqr_h.jl b/base/sparse/spqr_h.jl deleted file mode 100644 index a9f0985b4e1f1..0000000000000 --- a/base/sparse/spqr_h.jl +++ /dev/null @@ -1,30 +0,0 @@ -## SuiteSparseQR - -## ordering options -const SPQR_ORDERING_FIXED = int32(0) -const SPQR_ORDERING_NATURAL = int32(1) -const SPQR_ORDERING_COLAMD = int32(2) -const SPQR_ORDERING_GIVEN = int32(3) # only used for C/C++ interface -const SPQR_ORDERING_CHOLMOD = int32(4) # CHOLMOD best-effort (COLAMD, METIS,...) -const SPQR_ORDERING_AMD = int32(5) # AMD(A'*A) -const SPQR_ORDERING_METIS = int32(6) # metis(A'*A) -const SPQR_ORDERING_DEFAULT = int32(7) # SuiteSparseQR default ordering -const SPQR_ORDERING_BEST = int32(8) # try COLAMD, AMD, and METIS; pick best -const SPQR_ORDERING_BESTAMD = int32(9) # try COLAMD and AMD; pick best - -# Let [m n] = size of the matrix after pruning singletons. The default -# ordering strategy is to use COLAMD if m <= 2*n. Otherwise, AMD(A'A) is -# tried. If there is a high fill-in with AMD then try METIS(A'A) and take -# the best of AMD and METIS. METIS is not tried if it isn't installed. - -## Operations in qmult -const SPQR_QTX = int32(0) # Y = Q'*X -const SPQR_QX = int32(1) # Y = Q*X -const SPQR_XQT = int32(2) # Y = X*Q' -const SPQR_XQ = int32(3) # Y = X*Q - -## Types of systems to solve -const SPQR_RX_EQUALS_B = int32(0) # solve R*X=B or X = R\B -const SPQR_RETX_EQUALS_B = int32(1) # solve R*E'*X=B or X = E*(R\B) -const SPQR_RTX_EQUALS_B = int32(2) # solve R'*X=B or X = R'\B -const SPQR_RTX_EQUALS_ETB = int32(3) # solve R'*X=E'*B or X = R'\(E'*B) diff --git a/test/runtests.jl b/test/runtests.jl index f499145432a14..1f393f7148a81 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,7 +3,7 @@ testnames = [ "linalg", "core", "keywordargs", "numbers", "strings", "dates", "dict", "hashing", "remote", "iobuffer", "staged", "arrayops", "subarray", "reduce", "reducedim", "random", "intfuncs", - "simdloop", "blas", "fft", "dsp", "sparse", "bitarray", "copy", "math", + "simdloop", "blas", "fft", "dsp", "sparsetests", "bitarray", "copy", "math", "fastmath", "functional", "bigint", "sorting", "statistics", "spawn", "backtrace", "priorityqueue", "arpack", "file", "version", "resolve", "pollfd", "mpfr", "broadcast", "complex", "socket", @@ -25,12 +25,6 @@ push!(testnames, "parallel") tests = (ARGS==["all"] || isempty(ARGS)) ? testnames : ARGS -if "sparse" in tests - # specifically selected case - filter!(x -> x != "sparse", tests) - prepend!(tests, ["sparse/sparse", "sparse/cholmod", "sparse/umfpack"]) -end - if "linalg" in tests # specifically selected case filter!(x -> x != "linalg", tests) diff --git a/test/sparse/spqr.jl b/test/sparse/spqr.jl new file mode 100644 index 0000000000000..8f1861da0b9f1 --- /dev/null +++ b/test/sparse/spqr.jl @@ -0,0 +1,43 @@ +using Base.Test + +using Base.SparseMatrix.SPQR +using Base.SparseMatrix.CHOLMOD + +m, n = 1000, 10 +nn = 100 + +for eltyA in (Float64, Complex{Float64}) + for eltyB in (Float64, Complex{Float64}) + if eltyA <: Real + A = sparse(rand(1:m, nn), rand(1:n, nn), randn(nn), m, n) + else + A = sparse(rand(1:m, nn), rand(1:n, nn), complex(randn(nn), randn(nn)), m, n) + end + if eltyB <: Real + B = randn(m, 2) + else + B = complex(randn(m, 2), randn(m, 2)) + end + + @test_approx_eq A\B[:,1] full(A)\B[:,1] + @test_approx_eq A\B full(A)\B + @test_throws DimensionMismatch A\B[1:m-1,:] + + if eltyA == eltyB # promotions not defined for unexported methods + F = qrfact(A) + @test size(F) == (m,n) + @test size(F, 1) == m + @test size(F, 2) == n + @test size(F, 3) == 1 + @test_throws ArgumentError size(F, 0) + + # low level wrappers + @test_throws DimensionMismatch SPQR.solve(SPQR.RX_EQUALS_B, F, CHOLMOD.Dense(B')) + @test_throws DimensionMismatch SPQR.solve(SPQR.RTX_EQUALS_B, F, CHOLMOD.Dense(B)) + @test_throws DimensionMismatch SPQR.qmult(SPQR.QX, F, CHOLMOD.Dense(B')) + @test_throws DimensionMismatch SPQR.qmult(SPQR.XQ, F, CHOLMOD.Dense(B)) + @test_approx_eq A\B SPQR.backslash(SPQR.ORDERING_DEFAULT, SPQR.DEFAULT_TOL, CHOLMOD.Sparse(A), CHOLMOD.Dense(B)) + @test_throws DimensionMismatch SPQR.backslash(SPQR.ORDERING_DEFAULT, SPQR.DEFAULT_TOL, CHOLMOD.Sparse(A), CHOLMOD.Dense(B[1:m-1,:])) + end + end +end diff --git a/test/sparsetests.jl b/test/sparsetests.jl new file mode 100644 index 0000000000000..0b6b56472afa9 --- /dev/null +++ b/test/sparsetests.jl @@ -0,0 +1,4 @@ +include("sparse/sparse.jl") +include("sparse/umfpack.jl") +include("sparse/cholmod.jl") +include("sparse/spqr.jl") \ No newline at end of file