Remove inv_value and inv_diag (#146)

* Remove `inv_value` and `inv_diag`

* Fix compatibility with older Julia versions

* Add test

* Update README.md

Author: devmotion

README.md

@@ -50,17 +50,14 @@ PDMat(chol)         # with the Cholesky factorization
 * `PDiagMat`: diagonal matrix, defined as
 
 ```julia
-struct PDiagMat{T<:Real,V<:AbstractVector} <: AbstractPDMat{T}
+struct PDiagMat{T<:Real,V<:AbstractVector{T}} <: AbstractPDMat{T}
     dim::Int                    # matrix dimension
     diag::V                     # the vector of diagonal elements
-    inv_diag::V                 # the element-wise inverse of diag
 end
 
 # Constructors
 
-PDiagMat(v,inv_v)   # with the vector of diagonal elements and its inverse
 PDiagMat(v)         # with the vector of diagonal elements
-                    # inv_diag will be computed upon construction
 ```
 
 
@@ -70,13 +67,11 @@ PDiagMat(v)         # with the vector of diagonal elements
 struct ScalMat{T<:Real} <: AbstractPDMat{T}
     dim::Int         # matrix dimension
     value::T         # diagonal value (shared by all diagonal elements)
-    inv_value::T     # inv(value)
 end
 
 
 # Constructors
 
-ScalMat(d, v, inv_v) # with dimension d, diagonal value v and its inverse inv_v
 ScalMat(d, v)        # with dimension d and diagonal value v
 ```
 

src/deprecates.jl

@@ -11,3 +11,6 @@ using Base: @deprecate
 @deprecate full(x::AbstractPDMat) Matrix(x)
 
 @deprecate CholType Cholesky
+
+@deprecate ScalMat(d::Int, x::Real, inv_x::Real) ScalMat(d, x)
+@deprecate PDiagMat(v::AbstractVector, inv_v::AbstractVector) PDiagMat(v)

src/pdiagmat.jl

@@ -1,21 +1,12 @@
 """
 Positive definite diagonal matrix.
 """
-struct PDiagMat{T<:Real,V<:AbstractVector} <: AbstractPDMat{T}
+struct PDiagMat{T<:Real,V<:AbstractVector{T}} <: AbstractPDMat{T}
     dim::Int
     diag::V
-    inv_diag::V
-
-    PDiagMat{T,S}(d::Int,v::AbstractVector,inv_v::AbstractVector) where {T,S} =
-        new{T,S}(d,v,inv_v)
-end
-
-function PDiagMat(v::AbstractVector,inv_v::AbstractVector)
-    @check_argdims length(v) == length(inv_v)
-    PDiagMat{eltype(v),typeof(v)}(length(v), v, inv_v)
 end
 
-PDiagMat(v::AbstractVector) = PDiagMat(v, inv.(v))
+PDiagMat(v::AbstractVector{<:Real}) = PDiagMat{eltype(v),typeof(v)}(length(v), v)
 
 ### Conversion
 Base.convert(::Type{PDiagMat{T}},      a::PDiagMat) where {T<:Real} = PDiagMat(convert(AbstractArray{T}, a.diag))
@@ -51,13 +42,13 @@ 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) = a.inv_diag .* x
-/(x::AbstractVecOrMat, a::PDiagMat) = a.inv_diag .* x
+\(a::PDiagMat, x::AbstractVecOrMat) = x ./ a.diag
+/(x::AbstractVecOrMat, a::PDiagMat) = x ./ a.diag
 Base.kron(A::PDiagMat, B::PDiagMat) = PDiagMat( vcat([A.diag[i] * B.diag for i in 1:dim(A)]...) )
 
 ### Algebra
 
-Base.inv(a::PDiagMat) = PDiagMat(a.inv_diag, a.diag)
+Base.inv(a::PDiagMat) = PDiagMat(map(inv, a.diag))
 function LinearAlgebra.logdet(a::PDiagMat)
     diag = a.diag
     return isempty(diag) ? zero(log(zero(eltype(diag)))) : sum(log, diag)
@@ -71,9 +62,9 @@ LinearAlgebra.eigmin(a::PDiagMat) = minimum(a.diag)
 function whiten!(r::StridedVector, a::PDiagMat, x::StridedVector)
     n = dim(a)
     @check_argdims length(r) == length(x) == n
-    v = a.inv_diag
+    v = a.diag
     for i = 1:n
-        r[i] = x[i] * sqrt(v[i])
+        r[i] = x[i] / sqrt(v[i])
     end
     return r
 end
@@ -88,17 +79,21 @@ function unwhiten!(r::StridedVector, a::PDiagMat, x::StridedVector)
     return r
 end
 
-whiten!(r::StridedMatrix, a::PDiagMat, x::StridedMatrix) =
-    broadcast!(*, r, x, sqrt.(a.inv_diag))
+function whiten!(r::StridedMatrix, a::PDiagMat, x::StridedMatrix)
+    r .= x ./ sqrt.(a.diag)
+    return r
+end
 
-unwhiten!(r::StridedMatrix, a::PDiagMat, x::StridedMatrix) =
-    broadcast!(*, r, x, sqrt.(a.diag))
+function unwhiten!(r::StridedMatrix, a::PDiagMat, x::StridedMatrix)
+    r .= x .* sqrt.(a.diag)
+    return r
+end
 
 
 ### quadratic forms
 
 quad(a::PDiagMat, x::AbstractVector) = wsumsq(a.diag, x)
-invquad(a::PDiagMat, x::AbstractVector) = wsumsq(a.inv_diag, x)
+invquad(a::PDiagMat, x::AbstractVector) = invwsumsq(a.diag, x)
 
 function quad!(r::AbstractArray, a::PDiagMat, x::StridedMatrix)
     m, n = size(x)
@@ -116,12 +111,12 @@ end
 
 function invquad!(r::AbstractArray, a::PDiagMat, x::StridedMatrix)
     m, n = size(x)
-    ainvd = a.inv_diag
-    @check_argdims m == length(ainvd) && length(r) == n
+    ad = a.diag
+    @check_argdims m == length(ad) && length(r) == n
     @inbounds for j = 1:n
-        s = zero(promote_type(eltype(ainvd), eltype(x)))
+        s = zero(zero(eltype(x)) / zero(eltype(ad)))
         for i in 1:m
-            s += ainvd[i] * abs2(x[i,j])
+            s += abs2(x[i,j]) / ad[i]
         end
         r[j] = s
     end
@@ -132,7 +127,7 @@ end
 ### tri products
 
 function X_A_Xt(a::PDiagMat, x::StridedMatrix)
-    z = x .* reshape(sqrt.(a.diag), 1, dim(a))
+    z = x .* sqrt.(reshape(a.diag, 1, dim(a)))
     z * transpose(z)
 end
 
@@ -142,11 +137,11 @@ function Xt_A_X(a::PDiagMat, x::StridedMatrix)
 end
 
 function X_invA_Xt(a::PDiagMat, x::StridedMatrix)
-    z = x .* reshape(sqrt.(a.inv_diag), 1, dim(a))
+    z = x ./ sqrt.(reshape(a.diag, 1, dim(a)))
     z * transpose(z)
 end
 
 function Xt_invA_X(a::PDiagMat, x::StridedMatrix)
-    z = x .* sqrt.(a.inv_diag)
+    z = x ./ sqrt.(a.diag)
     transpose(z) * z
 end

src/scalmat.jl

@@ -4,13 +4,10 @@ Scaling matrix.
 struct ScalMat{T<:Real} <: AbstractPDMat{T}
     dim::Int
     value::T
-    inv_value::T
 end
 
-ScalMat(d::Int,v::Real) = ScalMat{typeof(inv(v))}(d, v, inv(v))
-
 ### Conversion
-Base.convert(::Type{ScalMat{T}}, a::ScalMat) where {T<:Real} = ScalMat(a.dim, T(a.value), T(a.inv_value))
+Base.convert(::Type{ScalMat{T}}, a::ScalMat) where {T<:Real} = ScalMat(a.dim, T(a.value))
 Base.convert(::Type{AbstractArray{T}}, a::ScalMat) where {T<:Real} = convert(ScalMat{T}, a)
 
 ### Basics
@@ -44,13 +41,13 @@ end
 /(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) = a.inv_value * x
-/(x::AbstractVecOrMat, a::ScalMat) = a.inv_value * x
+\(a::ScalMat, x::AbstractVecOrMat) = x / a.value
+/(x::AbstractVecOrMat, a::ScalMat) = x / a.value
 Base.kron(A::ScalMat, B::ScalMat) = ScalMat( dim(A) * dim(B), A.value * B.value )
 
 ### Algebra
 
-Base.inv(a::ScalMat) = ScalMat(a.dim, a.inv_value, a.value)
+Base.inv(a::ScalMat) = ScalMat(a.dim, inv(a.value))
 LinearAlgebra.logdet(a::ScalMat) = a.dim * log(a.value)
 LinearAlgebra.eigmax(a::ScalMat) = a.value
 LinearAlgebra.eigmin(a::ScalMat) = a.value
@@ -60,7 +57,7 @@ LinearAlgebra.eigmin(a::ScalMat) = a.value
 
 function whiten!(r::StridedVecOrMat, a::ScalMat, x::StridedVecOrMat)
     @check_argdims dim(a) == size(x, 1)
-    mul!(r, x, sqrt(a.inv_value))
+    _ldiv!(r, sqrt(a.value), x)
 end
 
 function unwhiten!(r::StridedVecOrMat, a::ScalMat, x::StridedVecOrMat)
@@ -72,10 +69,10 @@ end
 ### quadratic forms
 
 quad(a::ScalMat, x::AbstractVector) = sum(abs2, x) * a.value
-invquad(a::ScalMat, x::AbstractVector) = sum(abs2, x) * a.inv_value
+invquad(a::ScalMat, x::AbstractVector) = sum(abs2, x) / a.value
 
 quad!(r::AbstractArray, a::ScalMat, x::StridedMatrix) = colwise_sumsq!(r, x, a.value)
-invquad!(r::AbstractArray, a::ScalMat, x::StridedMatrix) = colwise_sumsq!(r, x, a.inv_value)
+invquad!(r::AbstractArray, a::ScalMat, x::StridedMatrix) = colwise_sumsqinv!(r, x, a.value)
 
 
 ### tri products
@@ -92,10 +89,10 @@ end
 
 function X_invA_Xt(a::ScalMat, x::StridedMatrix)
     @check_argdims dim(a) == size(x, 2)
-    lmul!(a.inv_value, x * transpose(x))
+    _rdiv!(x * transpose(x), a.value)
 end
 
 function Xt_invA_X(a::ScalMat, x::StridedMatrix)
     @check_argdims dim(a) == size(x, 1)
-    lmul!(a.inv_value, transpose(x) * x)
+    _rdiv!(transpose(x) * x, a.value)
 end

src/utils.jl

@@ -59,6 +59,15 @@ function wsumsq(w::AbstractVector, a::AbstractVector)
     return s
 end
 
+function invwsumsq(w::AbstractVector, a::AbstractVector)
+    @check_argdims(length(a) == length(w))
+    s = zero(zero(eltype(a)) / zero(eltype(w)))
+    for i = 1:length(a)
+        @inbounds s += abs2(a[i]) / w[i]
+    end
+    return s
+end
+
 function colwise_dot!(r::AbstractArray, a::AbstractMatrix, b::AbstractMatrix)
     n = length(r)
     @check_argdims n == size(a, 2) == size(b, 2) && size(a, 1) == size(b, 1)
@@ -84,3 +93,37 @@ function colwise_sumsq!(r::AbstractArray, a::AbstractMatrix, c::Real)
     end
     return r
 end
+
+function colwise_sumsqinv!(r::AbstractArray, a::AbstractMatrix, c::Real)
+    n = length(r)
+    @check_argdims n == size(a, 2)
+    for j = 1:n
+        v = zero(eltype(a))
+        @simd for i = 1:size(a, 1)
+            @inbounds v += abs2(a[i, j])
+        end
+        r[j] = v / c
+    end
+    return r
+end
+
+# `rdiv!(::AbstractArray, ::Number)` was introduced in Julia 1.2
+# https://github.com/JuliaLang/julia/pull/31179
+@static if VERSION < v"1.2.0-DEV.385"
+    function _rdiv!(X::AbstractArray, s::Number)
+        @simd for I in eachindex(X)
+            @inbounds X[I] /= s
+        end
+        X
+    end
+else
+    _rdiv!(X::AbstractArray, s::Number) = rdiv!(X, s)
+end
+
+# `ldiv!(::AbstractArray, ::Number, ::AbstractArray)` was introduced in Julia 1.4
+# https://github.com/JuliaLang/julia/pull/33806
+@static if VERSION < v"1.4.0-DEV.635"
+    _ldiv!(Y::AbstractArray, s::Number, X::AbstractArray) = Y .= s .\ X
+else
+    _ldiv!(Y::AbstractArray, s::Number, X::AbstractArray) = ldiv!(Y, s, X)
+end

test/pdmtypes.jl

@@ -7,9 +7,9 @@ using Test
             m = Matrix{T}(I, 2, 2)
             @test PDMat(m, cholesky(m)).mat == PDMat(Symmetric(m)).mat == PDMat(m).mat == PDMat(cholesky(m)).mat
             d = ones(T,2)
-            @test PDiagMat(d,d).inv_diag == PDiagMat(d).inv_diag
+            @test @test_deprecated(PDiagMat(d, d)) == PDiagMat(d)
             x = one(T)
-            @test ScalMat(2,x,x).inv_value == ScalMat(2,x).inv_value
+            @test @test_deprecated(ScalMat(2, x, x)) == ScalMat(2, x)
             s = SparseMatrixCSC{T}(I, 2, 2)
             @test PDSparseMat(s, cholesky(s)).mat == PDSparseMat(s).mat == PDSparseMat(cholesky(s)).mat
         end
@@ -75,4 +75,12 @@ using Test
     @testset "convert Matrix type to the same Cholesky type (#117)" begin
         @test PDMat([1 0; 0 1]) == [1.0 0.0; 0.0 1.0]
     end
+
+    # https://github.com/JuliaStats/PDMats.jl/pull/141
+    @testset "PDiagMat with range" begin
+        v = 0.1:0.1:0.5
+        d = PDiagMat(v)
+        @test d isa PDiagMat{eltype(v),typeof(v)}
+        @test d.diag === v
+    end
 end