Fix right division and add checks of dimensions (#153)

Author: devmotion

Project.toml

@@ -1,6 +1,6 @@
 name = "PDMats"
 uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
-version = "0.11.5"
+version = "0.11.6"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/pdiagmat.jl

@@ -39,11 +39,24 @@ function pdadd!(r::Matrix, a::Matrix, b::PDiagMat, c)
     return r
 end
 
-*(a::PDiagMat, c::T) where {T<:Real} = PDiagMat(a.diag * c)
-*(a::PDiagMat, x::AbstractVector) = a.diag .* x
-*(a::PDiagMat, x::AbstractMatrix) = a.diag .* x
-\(a::PDiagMat, x::AbstractVecOrMat) = x ./ a.diag
-/(x::AbstractVecOrMat, a::PDiagMat) = x ./ a.diag
+*(a::PDiagMat, c::Real) = PDiagMat(a.diag * c)
+function *(a::PDiagMat, x::AbstractVector)
+    @check_argdims dim(a) == length(x)
+    return a.diag .* x
+end
+function *(a::PDiagMat, x::AbstractMatrix)
+    @check_argdims dim(a) == size(x, 1)
+    return a.diag .* x
+end
+function \(a::PDiagMat, x::AbstractVecOrMat)
+    @check_argdims dim(a) == size(x, 1)
+    return x ./ a.diag
+end
+function /(x::AbstractVecOrMat, a::PDiagMat)
+    @check_argdims dim(a) == size(x, 2)
+    # return matrix for 1-element vectors `x`, consistent with LinearAlgebra
+    return reshape(x, Val(2)) ./ permutedims(a.diag) # = (a' \ x')'
+end
 Base.kron(A::PDiagMat, B::PDiagMat) = PDiagMat( vcat([A.diag[i] * B.diag for i in 1:dim(A)]...) )
 
 ### Algebra
@@ -127,21 +140,25 @@ end
 ### tri products
 
 function X_A_Xt(a::PDiagMat, x::StridedMatrix)
-    z = x .* sqrt.(reshape(a.diag, 1, dim(a)))
+    @check_argdims dim(a) == size(x, 2)
+    z = x .* sqrt.(permutedims(a.diag))
     z * transpose(z)
 end
 
 function Xt_A_X(a::PDiagMat, x::StridedMatrix)
+    @check_argdims dim(a) == size(x, 1)
     z = x .* sqrt.(a.diag)
     transpose(z) * z
 end
 
 function X_invA_Xt(a::PDiagMat, x::StridedMatrix)
-    z = x ./ sqrt.(reshape(a.diag, 1, dim(a)))
+    @check_argdims dim(a) == size(x, 2)
+    z = x ./ sqrt.(permutedims(a.diag))
     z * transpose(z)
 end
 
 function Xt_invA_X(a::PDiagMat, x::StridedMatrix)
+    @check_argdims dim(a) == size(x, 1)
     z = x ./ sqrt.(a.diag)
     transpose(z) * z
 end

src/pdmat.jl

@@ -44,11 +44,12 @@ function pdadd!(r::Matrix, a::Matrix, b::PDMat, c)
     _addscal!(r, a, b.mat, c)
 end
 
-*(a::PDMat{S}, c::T) where {S<:Real, T<:Real} = PDMat(a.mat * c)
-*(a::PDMat, x::AbstractVector{T}) where {T} = a.mat * x
-*(a::PDMat, x::AbstractMatrix{T}) where {T} = a.mat * x
+*(a::PDMat, c::Real) = PDMat(a.mat * c)
+*(a::PDMat, x::AbstractVector) = a.mat * x
+*(a::PDMat, x::AbstractMatrix) = a.mat * x
 \(a::PDMat, x::AbstractVecOrMat) = a.chol \ x
-/(x::AbstractVecOrMat, a::PDMat) = x / a.chol
+# return matrix for 1-element vectors `x`, consistent with LinearAlgebra
+/(x::AbstractVecOrMat, a::PDMat) = reshape(x, Val(2)) / a.chol
 
 ### Algebra
 
@@ -94,21 +95,25 @@ invquad!(r::AbstractArray, a::PDMat, x::StridedMatrix) = colwise_dot!(r, x, a.ma
 ### tri products
 
 function X_A_Xt(a::PDMat, x::StridedMatrix)
-    z = rmul!(copy(x), chol_lower(a.chol))
+    @check_argdims dim(a) == size(x, 2)
+    z = x * chol_lower(a.chol)
     return z * transpose(z)
 end
 
 function Xt_A_X(a::PDMat, x::StridedMatrix)
-    z = lmul!(chol_upper(a.chol), copy(x))
+    @check_argdims dim(a) == size(x, 1)
+    z = chol_upper(a.chol) * x
     return transpose(z) * z
 end
 
 function X_invA_Xt(a::PDMat, x::StridedMatrix)
-    z = rdiv!(copy(x), chol_upper(a.chol))
+    @check_argdims dim(a) == size(x, 2)
+    z = x / chol_upper(a.chol)
     return z * transpose(z)
 end
 
 function Xt_invA_X(a::PDMat, x::StridedMatrix)
-    z = ldiv!(chol_lower(a.chol), copy(x))
+    @check_argdims dim(a) == size(x, 1)
+    z = chol_lower(a.chol) \ x
     return transpose(z) * z
 end

src/pdsparsemat.jl

@@ -43,7 +43,7 @@ function pdadd!(r::Matrix, a::Matrix, b::PDSparseMat, c)
     _addscal!(r, a, b.mat, c)
 end
 
-*(a::PDSparseMat, c::T) where {T<:Real} = PDSparseMat(a.mat * c)
+*(a::PDSparseMat, c::Real) = PDSparseMat(a.mat * c)
 *(a::PDSparseMat, x::StridedVecOrMat) = a.mat * x
 \(a::PDSparseMat{T}, x::StridedVecOrMat{T}) where {T<:Real} = convert(Array{T},a.chol \ convert(Array{Float64},x)) #to avoid limitations in sparse factorization library CHOLMOD, see e.g., julia issue #14076
 /(x::StridedVecOrMat{T}, a::PDSparseMat{T}) where {T<:Real} = convert(Array{T},convert(Array{Float64},x) / a.chol )

src/scalmat.jl

@@ -37,12 +37,25 @@ function pdadd!(r::Matrix, a::Matrix, b::ScalMat, c)
     return r
 end
 
-*(a::ScalMat, c::T) where {T<:Real} = ScalMat(a.dim, a.value * c)
-/(a::ScalMat{T}, c::T) where {T<:Real} = ScalMat(a.dim, a.value / c)
-*(a::ScalMat, x::AbstractVector) = a.value * x
-*(a::ScalMat, x::AbstractMatrix) = a.value * x
-\(a::ScalMat, x::AbstractVecOrMat) = x / a.value
-/(x::AbstractVecOrMat, a::ScalMat) = x / a.value
+*(a::ScalMat, c::Real) = ScalMat(a.dim, a.value * c)
+/(a::ScalMat, c::Real) = ScalMat(a.dim, a.value / c)
+function *(a::ScalMat, x::AbstractVector)
+    @check_argdims dim(a) == length(x)
+    return a.value * x
+end
+function *(a::ScalMat, x::AbstractMatrix)
+    @check_argdims dim(a) == size(x, 1)
+    return a.value * x
+end
+function \(a::ScalMat, x::AbstractVecOrMat)
+    @check_argdims dim(a) == size(x, 1)
+    return x / a.value
+end
+function /(x::AbstractVecOrMat, a::ScalMat)
+    @check_argdims dim(a) == size(x, 2)
+    # return matrix for 1-element vectors `x`, consistent with LinearAlgebra
+    return reshape(x, Val(2)) / a.value
+end
 Base.kron(A::ScalMat, B::ScalMat) = ScalMat( dim(A) * dim(B), A.value * B.value )
 
 ### Algebra

test/pdmtypes.jl

@@ -83,4 +83,30 @@ using Test
         @test d isa PDiagMat{eltype(v),typeof(v)}
         @test d.diag === v
     end
+
+    @testset "division of vectors (dim = 1)" begin
+        A = rand(1, 1)
+        x = randn(1)
+        y = x / A
+        @assert x / A isa Matrix{Float64}
+        @assert size(y) == (1, 1)
+
+        for M in (PDiagMat(vec(A)), ScalMat(1, first(A)))
+            @test x / M isa Matrix{Float64}
+            @test x / M ≈ y
+        end
+
+        # requires https://github.com/JuliaLang/julia/pull/32594
+        if VERSION >= v"1.3.0-DEV.562"
+            @test x / PDMat(A) isa Matrix{Float64}
+            @test x / PDMat(A) ≈ y
+        end
+
+        # right division not defined for CHOLMOD:
+        # `rdiv!(::Matrix{Float64}, ::SuiteSparse.CHOLMOD.Factor{Float64})` not defined
+        if !HAVE_CHOLMOD
+            @test x / PDSparseMat(sparse(first(A), 1, 1)) isa Matrix{Float64}
+            @test x / PDSparseMat(sparse(first(A), 1, 1)) ≈ y
+        end
+    end
 end

test/testutils.jl

@@ -20,7 +20,7 @@ function test_pdmat(C::AbstractPDMat, Cmat::Matrix;
                     t_logdet::Bool=true,        # whether to test logdet method
                     t_eig::Bool=true,           # whether to test eigmax and eigmin
                     t_mul::Bool=true,           # whether to test multiplication
-                    t_rdiv::Bool=true,          # whether to test right division (solve)
+                    t_div::Bool=true,           # whether to test division
                     t_quad::Bool=true,          # whether to test quad & invquad
                     t_triprod::Bool=true,       # whether to test X_A_Xt, Xt_A_X, X_invA_Xt, and Xt_invA_X
                     t_whiten::Bool=true         # whether to test whiten and unwhiten
@@ -45,7 +45,7 @@ function test_pdmat(C::AbstractPDMat, Cmat::Matrix;
     X = rand(eltype(C),d,n) .- convert(eltype(C),0.5)
 
     t_mul && pdtest_mul(C, Cmat, X, verbose)
-    t_rdiv && pdtest_rdiv(C, Imat, X, verbose)
+    t_div && pdtest_div(C, Imat, X, verbose)
     t_quad && pdtest_quad(C, Cmat, Imat, X, verbose)
     t_triprod && pdtest_triprod(C, Cmat, Imat, X, verbose)
 
@@ -174,42 +174,60 @@ end
 
 function pdtest_mul(C::AbstractPDMat, Cmat::Matrix, verbose::Int)
     n = 5
-    X = rand(eltype(C),dim(C), n)
-
-    _pdt(verbose, "multiply")
-    @test C * X ≈ Cmat * X
-
-    for i = 1:n
-        xi = vec(copy(X[:,i]))
-        @test C * xi ≈ Cmat * xi
-    end
+    X = rand(eltype(C), dim(C), n)
+    pdtest_mul(C, Cmat, X, verbose)
 end
 
 
 function pdtest_mul(C::AbstractPDMat, Cmat::Matrix, X::Matrix, verbose::Int)
     _pdt(verbose, "multiply")
+    d, n = size(X)
+    @assert d == dim(C)
+    @assert size(Cmat) == size(C)
     @test C * X ≈ Cmat * X
 
-    y = similar(C * X, size(C, 1))
-    ymat = similar(Cmat * X, size(Cmat, 1))
-    for i = 1:size(X,2)
+    y = similar(C * X, d)
+    ymat = similar(Cmat * X, d)
+    for i = 1:n
         xi = vec(copy(X[:,i]))
         @test C * xi ≈ Cmat * xi
 
         mul!(y, C, xi)
         mul!(ymat, Cmat, xi)
         @test y ≈ ymat
     end
+
+    # Dimension mismatches
+    @test_throws DimensionMismatch C * rand(d + 1) 
+    @test_throws DimensionMismatch C * rand(d + 1, n)
 end
 
 
-function pdtest_rdiv(C::AbstractPDMat, Imat::Matrix, X::Matrix, verbose::Int)
-    _pdt(verbose, "rdivide")
+function pdtest_div(C::AbstractPDMat, Imat::Matrix, X::Matrix, verbose::Int)
+    _pdt(verbose, "divide")
+    d, n = size(X)
+    @assert d == dim(C)
+    @assert size(Imat) == size(C)
     @test C \ X ≈ Imat * X
+    # Right division with Choleskyrequires https://github.com/JuliaLang/julia/pull/32594
+    # CHOLMOD throws error since no method is found for
+    # `rdiv!(::Matrix{Float64}, ::SuiteSparse.CHOLMOD.Factor{Float64})`
+    check_rdiv = !(C isa PDMat && VERSION < v"1.3.0-DEV.562") && !(C isa PDSparseMat && HAVE_CHOLMOD)
+    check_rdiv && @test Matrix(X') / C ≈ (C \ X)'
 
-    for i = 1:size(X,2)
+    for i = 1:n
         xi = vec(copy(X[:,i]))
         @test C \ xi ≈ Imat * xi
+        check_rdiv && @test Matrix(xi') / C ≈ (C \ xi)'
+    end
+
+
+    # Dimension mismatches
+    @test_throws DimensionMismatch C \ rand(d + 1) 
+    @test_throws DimensionMismatch C \ rand(d + 1, n)
+    if check_rdiv
+        @test_throws DimensionMismatch rand(1, d + 1) / C
+        @test_throws DimensionMismatch rand(n, d + 1) / C
     end
 end
 
@@ -240,23 +258,29 @@ end
 
 
 function pdtest_triprod(C::AbstractPDMat, Cmat::Matrix, Imat::Matrix, X::Matrix, verbose::Int)
+    d, n = size(X)
+    @assert d == dim(C)
     Xt = copy(transpose(X))
 
     _pdt(verbose, "X_A_Xt")
     # default tolerance in isapprox is different on 0.4. rtol argument can be deleted
     # ≈ form used when 0.4 is no longer supported
     lhs, rhs = X_A_Xt(C, Xt), Xt * Cmat * X
     @test isapprox(lhs, rhs, rtol=sqrt(max(eps(real(float(eltype(lhs)))), eps(real(float(eltype(rhs)))))))
+    @test_throws DimensionMismatch X_A_Xt(C, rand(n, d + 1))
 
     _pdt(verbose, "Xt_A_X")
     lhs, rhs = Xt_A_X(C, X), Xt * Cmat * X
     @test isapprox(lhs, rhs, rtol=sqrt(max(eps(real(float(eltype(lhs)))), eps(real(float(eltype(rhs)))))))
+    @test_throws DimensionMismatch Xt_A_X(C, rand(d + 1, n))
 
     _pdt(verbose, "X_invA_Xt")
     @test X_invA_Xt(C, Xt) ≈ Xt * Imat * X
+    @test_throws DimensionMismatch X_invA_Xt(C, rand(n, d + 1))
 
     _pdt(verbose, "Xt_invA_X")
     @test Xt_invA_X(C, X) ≈ Xt * Imat * X
+    @test_throws DimensionMismatch Xt_invA_X(C, rand(d + 1, n))
 end