Merge pull request #221 from JuliaStats/dw/invwhiten_invunwhiten

Add `invwhiten` and `invunwhiten`

Author: andreasnoack

Project.toml

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

README.md

@@ -180,7 +180,7 @@ X_invA_Xt(a, x)     # compute `x * inv(a) * x'` for a matrix `x`.
 
 Xt_invA_X(a, x)     # compute `x' * inv(a) * x` for a matrix `x`.
 
-whiten(a, x)        # whitening transform. `x` can be a vector or a matrix.
+whiten(a, x)        # whitening transform defined by `a`. `x` can be a vector or a matrix.
                     #
                     # Note: If the covariance of `x` is `a`, then the
                     # covariance of the transformed result is an identity
@@ -190,8 +190,17 @@ whiten!(a, x)       # whitening transform inplace, directly updating `x`.
 
 whiten!(r, a, x)    # write the transformed result to `r`.
 
-unwhiten(a, x)      # inverse of whitening transform. `x` can be a vector or
-                    # a matrix.
+invwhiten(a, x)     # whitening transform defined by `inv(a)`. `x` can be a vector or a matrix.
+                    #
+                    # Note: If the precision of `x` is `a`, then the
+                    # covariance of the transformed result is an identity matrix.
+
+invwhiten!(a, x)    # whitening transform inplace, directly updating `x`.
+
+invwhiten!(r, a, x) # write the transformed result to `r`.
+
+unwhiten(a, x)      # inverse of whitening transform with `a`.
+                    # `x` can be a vector or a matrix.
                     #
                     # Note: If the covariance of `x` is an identity matrix,
                     # then the covariance of the transformed result is `a`.
@@ -201,6 +210,18 @@ unwhiten(a, x)      # inverse of whitening transform. `x` can be a vector or
 unwhiten!(a, x)     # un-whitening transform inplace, updating `x`.
 
 unwhiten!(r, a, x)  # write the transformed result to `r`.
+
+invunwhiten(a, x)     # inverse of whitening transform defined by `inv(a)`.
+                      # `x` can be a vector or a matrix.
+                      #
+                      # Note: If the covariance of `x` is an identity matrix,
+                      # then the precision of the transformed result is `a`.
+                      # Note: the un-whitening transform is useful for
+                      # generating Gaussian samples.
+
+invunwhiten!(a, x)    # un-whitening transform inplace, updating `x`.
+
+invunwhiten!(r, a, x) # write the transformed result to `r`.
 ```
 
 ### Fallbacks for `AbstractArray`s

src/PDMats.jl

@@ -16,8 +16,12 @@ export
     dim,
     whiten,
     whiten!,
+    invwhiten,
+    invwhiten!,
     unwhiten,
     unwhiten!,
+    invunwhiten,
+    invunwhiten!,
     pdadd,
     pdadd!,
     quad,

src/chol.jl

@@ -32,16 +32,22 @@ dim(A::Cholesky) = LinearAlgebra.checksquare(A)
 # whiten
 whiten(A::Cholesky, x::AbstractVecOrMat) = chol_lower(A) \ x
 whiten!(A::Cholesky, x::AbstractVecOrMat) = ldiv!(chol_lower(A), x)
+invwhiten(A::Cholesky, x::AbstractVecOrMat) = chol_upper(A) * x
+invwhiten!(A::Cholesky, x::AbstractVecOrMat) = lmul!(chol_upper(A), x)
 
 # unwhiten
 unwhiten(A::Cholesky, x::AbstractVecOrMat) = chol_lower(A) * x
 unwhiten!(A::Cholesky, x::AbstractVecOrMat) = lmul!(chol_lower(A), x)
+invunwhiten(A::Cholesky, x::AbstractVecOrMat) = chol_upper(A) \ x
+invunwhiten!(A::Cholesky, x::AbstractVecOrMat) = ldiv!(chol_upper(A), x)
 
 # 3-argument whiten/unwhiten
 for T in (:AbstractVector, :AbstractMatrix)
     @eval begin
         whiten!(r::$T, A::Cholesky, x::$T) = whiten!(A, copyto!(r, x))
+        invwhiten!(r::$T, A::Cholesky, x::$T) = invwhiten!(A, copyto!(r, x))
         unwhiten!(r::$T, A::Cholesky, x::$T) = unwhiten!(A, copyto!(r, x))
+        invunwhiten!(r::$T, A::Cholesky, x::$T) = invunwhiten!(A, copyto!(r, x))
     end
 end
 

src/generics.jl

@@ -38,7 +38,7 @@ LinearAlgebra.checksquare(a::AbstractPDMat) = size(a, 1)
     whiten!(a::AbstractMatrix, x::AbstractVecOrMat)
     whiten!(r::AbstractVecOrMat, a::AbstractPDMat, x::AbstractVecOrMat)
 
-Allocating and in-place versions of the `whiten`ing transform defined by `a` applied to `x`
+Allocating and in-place versions of the `whiten`ing transform defined by `a` applied to `x`.
 
 If the covariance of `x` is `a` then the covariance of the result will be `I`.
 The name `whiten` indicates that this function transforms a correlated multivariate random
@@ -66,34 +66,80 @@ julia> W * W'
  0.0  1.0
 ```
 
-See also [`unwhiten`](@ref).
+`whiten` is the inverse transformation of [`unwhiten`](@ref).
+
+See also [`invwhiten`](@ref) and [`invunwhiten`](ref)
 """
 whiten(a::AbstractMatrix{<:Real}, x::AbstractVecOrMat) = whiten(AbstractPDMat(a), x)
 
+"""
+    invwhiten(a::AbstractMatrix, x::AbstractVecOrMat)
+    invwhiten!(a::AbstractMatrix, x::AbstractVecOrMat)
+    invwhiten!(r::AbstractVecOrMat, a::AbstractPDMat, x::AbstractVecOrMat)
+
+Allocating and in-place versions of the `whiten`ing transform defined by `inv(a)` applied to `x`.
+
+If the precision of `x` is `a` then the covariance of the result will be `I`.
+The name `invwhiten` indicates that this function transforms a correlated multivariate random
+variable to "white noise".
+
+`invwhiten` is the inverse transformation of [`invunwhiten`](@ref).
+
+See also [`whiten`](@ref) and [`unwhiten`](@ref)
+"""
+invwhiten(a::AbstractMatrix{<:Real}, x::AbstractVecOrMat) = invwhiten(AbstractPDMat(a), x)
+
 """
     unwhiten(a::AbstractMatrix, x::AbstractVecOrMat)
     unwhiten!(a::AbstractMatrix, x::AbstractVecOrMat)
     unwhiten!(r::AbstractVecOrMat, a::AbstractPDMat, x::AbstractVecOrMat)
 
-Allocating and in-place versions of the `unwhiten`ing transform defined by `a` applied to `x`
+Allocating and in-place versions of the `unwhiten`ing transform defined by `a` applied to `x`.
 
 If the covariance of `x` is `I` then the covariance of the result will be `a`.
 The name `unwhiten` indicates that this function transforms an uncorrelated "white noise" signal
 into a signal with the given covariance.
 
 `unwhiten` is the inverse transformation of [`whiten`](@ref).
+
+See also [`invwhiten`](@ref) and [`invunwhiten`](@ref)
 """
 unwhiten(a::AbstractMatrix{<:Real}, x::AbstractVecOrMat) = unwhiten(AbstractPDMat(a), x)
 
+"""
+    invunwhiten(a::AbstractMatrix, x::AbstractVecOrMat)
+    invunwhiten!(a::AbstractMatrix, x::AbstractVecOrMat)
+    invunwhiten!(r::AbstractVecOrMat, a::AbstractPDMat, x::AbstractVecOrMat)
+
+Allocating and in-place versions of the `unwhiten`ing transform defined by `inv(a)` applied to `x`.
+
+If the covariance of `x` is `I` then the precision of the result will be `a`.
+The name `invunwhiten` indicates that this function transforms an uncorrelated "white noise" signal
+into a signal with the given covariance.
+
+`invunwhiten` is the inverse transformation of [`invwhiten`](@ref).
+
+See also [`whiten`](@ref) and [`unwhiten`](@ref)
+"""
+invunwhiten(a::AbstractMatrix{<:Real}, x::AbstractVecOrMat) = invunwhiten(AbstractPDMat(a), x)
+
 whiten!(a::AbstractMatrix{<:Real}, x::AbstractVecOrMat) = whiten!(x, a, x)
+invwhiten!(a::AbstractMatrix{<:Real}, x::AbstractVecOrMat) = invwhiten!(x, a, x)
 unwhiten!(a::AbstractMatrix{<:Real}, x::AbstractVecOrMat) = unwhiten!(x, a, x)
+invunwhiten!(a::AbstractMatrix{<:Real}, x::AbstractVecOrMat) = invunwhiten!(x, a, x)
 
 function whiten!(r::AbstractVecOrMat, a::AbstractMatrix{<:Real}, x::AbstractVecOrMat)
     return whiten!(r, AbstractPDMat(a), x)
 end
+function invwhiten!(r::AbstractVecOrMat, a::AbstractMatrix{<:Real}, x::AbstractVecOrMat)
+    return invwhiten!(r, AbstractPDMat(a), x)
+end
 function unwhiten!(r::AbstractVecOrMat, a::AbstractMatrix{<:Real}, x::AbstractVecOrMat)
     return unwhiten!(r, AbstractPDMat(a), x)
 end
+function invunwhiten!(r::AbstractVecOrMat, a::AbstractMatrix{<:Real}, x::AbstractVecOrMat)
+    return invunwhiten!(r, AbstractPDMat(a), x)
+end
 
 ## quad
 

src/pdiagmat.jl

@@ -90,20 +90,24 @@ function whiten!(r::AbstractVecOrMat, a::PDiagMat, x::AbstractVecOrMat)
     @check_argdims a.dim == size(x, 1)
     return r .= x ./ sqrt.(a.diag)
 end
+invwhiten!(r::AbstractVecOrMat, a::PDiagMat, x::AbstractVecOrMat) = unwhiten!(r, a, x)
 function unwhiten!(r::AbstractVecOrMat, a::PDiagMat, x::AbstractVecOrMat)
     @check_argdims axes(r) == axes(x)
     @check_argdims a.dim == size(x, 1)
     return r .= x .* sqrt.(a.diag)
 end
+invunwhiten!(r::AbstractVecOrMat, a::PDiagMat, x::AbstractVecOrMat) = whiten!(r, a, x)
 
 function whiten(a::PDiagMat, x::AbstractVecOrMat)
     @check_argdims a.dim == size(x, 1)
     return x ./ sqrt.(a.diag)
 end
+invwhiten(a::PDiagMat, x::AbstractVecOrMat) = unwhiten(a, x)
 function unwhiten(a::PDiagMat, x::AbstractVecOrMat)
     @check_argdims a.dim == size(x, 1)
     return x .* sqrt.(a.diag)
 end
+invunwhiten(a::PDiagMat, x::AbstractVecOrMat) = whiten(a, x)
 
 ### quadratic forms
 

src/pdmat.jl

@@ -131,6 +131,15 @@ function whiten!(r::AbstractVecOrMat, a::PDMat, x::AbstractVecOrMat)
         return ldiv!(r, chol_lower(cholesky(a)), x)
     end
 end
+function invwhiten!(r::AbstractVecOrMat, a::PDMat, x::AbstractVecOrMat)
+    @check_argdims axes(r) == axes(x)
+    @check_argdims a.dim == size(x, 1)
+    if r === x
+        return lmul!(chol_upper(cholesky(a)), r)
+    else
+        return mul!(r, chol_upper(cholesky(a)), x)
+    end
+end
 function unwhiten!(r::AbstractVecOrMat, a::PDMat, x::AbstractVecOrMat)
     @check_argdims axes(r) == axes(x)
     @check_argdims a.dim == size(x, 1)
@@ -140,15 +149,32 @@ function unwhiten!(r::AbstractVecOrMat, a::PDMat, x::AbstractVecOrMat)
         return mul!(r, chol_lower(cholesky(a)), x)
     end
 end
+function invunwhiten!(r::AbstractVecOrMat, a::PDMat, x::AbstractVecOrMat)
+    @check_argdims axes(r) == axes(x)
+    @check_argdims a.dim == size(x, 1)
+    if r === x
+        return ldiv!(chol_upper(cholesky(a)), r)
+    else
+        return ldiv!(r, chol_upper(cholesky(a)), x)
+    end
+end
 
 function whiten(a::PDMat, x::AbstractVecOrMat)
     @check_argdims a.dim == size(x, 1)
     return chol_lower(cholesky(a)) \ x
 end
+function invwhiten(a::PDMat, x::AbstractVecOrMat)
+    @check_argdims a.dim == size(x, 1)
+    return chol_upper(cholesky(a)) * x
+end
 function unwhiten(a::PDMat, x::AbstractVecOrMat)
     @check_argdims a.dim == size(x, 1)
     return chol_lower(cholesky(a)) * x
 end
+function invunwhiten(a::PDMat, x::AbstractVecOrMat)
+    @check_argdims a.dim == size(x, 1)
+    return chol_upper(cholesky(a)) \ x
+end
 
 ## quad/invquad
 

src/pdsparsemat.jl

@@ -86,30 +86,61 @@ function whiten!(r::AbstractVecOrMat, a::PDSparseMat, x::AbstractVecOrMat)
     # Can't use `ldiv!` due to missing support in SparseArrays
     return copyto!(r, chol_lower(a.chol) \ x)
 end
+function invwhiten!(r::AbstractVecOrMat, a::PDSparseMat, x::AbstractVecOrMat)
+    @check_argdims axes(r) == axes(x)
+    @check_argdims a.dim == size(x, 1)
+    # `*` and `mul!` are not defined for `UP` factor components,
+    # so we can't use `chol_upper(C) * x`
+    # `sparse` is neither defined for `PtL` nor for `UP` nor for `U` factor components
+    C = cholesky(a)
+    PtL = sparse(C.L)[C.p, :]
+    return copyto!(r, PtL' * x)
+end
 
 function unwhiten!(r::AbstractVecOrMat, a::PDSparseMat, x::AbstractVecOrMat)
     @check_argdims axes(r) == axes(x)
     @check_argdims a.dim == size(x, 1)
     # `*` is not defined for `PtL` factor components,
-    # so we can't use `chol_lower(a.chol) * x`
-    C = a.chol
+    # so we can't use `chol_lower(C) * x`
+    # `sparse` is neither defined for `PtL` nor for `UP` nor for `U` factor components
+    C = cholesky(a)
     PtL = sparse(C.L)[C.p, :]
     return copyto!(r, PtL * x)
 end
+function invunwhiten!(r::AbstractVecOrMat, a::PDSparseMat, x::AbstractVecOrMat)
+    @check_argdims axes(r) == axes(x)
+    @check_argdims a.dim == size(x, 1)
+    # Can't use `ldiv!` due to missing support in SparseArrays
+    return copyto!(r, chol_upper(cholesky(a)) \ x)
+end
 
 function whiten(a::PDSparseMat, x::AbstractVecOrMat)
     @check_argdims a.dim == size(x, 1)
     return chol_lower(cholesky(a)) \ x
 end
+function invwhiten(a::PDSparseMat, x::AbstractVecOrMat)
+    @check_argdims a.dim == size(x, 1)
+    # `*` is not defined for `UP` factor components,
+    # so we can't use `chol_upper(C) * x`
+    # `sparse` is neither defined for `PtL` nor for `UP` nor for `U` factor components
+    C = cholesky(a)
+    PtL = sparse(C.L)[C.p, :]
+    return PtL' * x
+end
 
 function unwhiten(a::PDSparseMat, x::AbstractVecOrMat)
     @check_argdims a.dim == size(x, 1)
     # `*` is not defined for `PtL` factor components,
-    # so we can't use `chol_lower(a.chol) * x`
-    C = a.chol
+    # so we can't use `chol_lower(C) * x`
+    # `sparse` is neither defined for `PtL` nor for `UP` nor for `U` factor components
+    C = cholesky(a)
     PtL = sparse(C.L)[C.p, :]
     return PtL * x
 end
+function invunwhiten(a::PDSparseMat, x::AbstractVecOrMat)
+    @check_argdims a.dim == size(x, 1)
+    return chol_upper(cholesky(a)) \ x
+end
 
 ### quadratic forms
 

src/scalmat.jl

@@ -75,21 +75,25 @@ function whiten!(r::AbstractVecOrMat, a::ScalMat, x::AbstractVecOrMat)
     @check_argdims a.dim == size(x, 1)
     return ldiv!(r, sqrt(a.value), x)
 end
+invwhiten!(r::AbstractVecOrMat, a::ScalMat, x::AbstractVecOrMat) = unwhiten!(r, a, x)
 
 function unwhiten!(r::AbstractVecOrMat, a::ScalMat, x::AbstractVecOrMat)
     @check_argdims axes(r) == axes(x)
     @check_argdims a.dim == size(x, 1)
     return mul!(r, x, sqrt(a.value))
 end
+invunwhiten!(r::AbstractVecOrMat, a::ScalMat, x::AbstractVecOrMat) = whiten!(r, a, x)
 
 function whiten(a::ScalMat, x::AbstractVecOrMat)
     @check_argdims a.dim == size(x, 1)
     return x / sqrt(a.value)
 end
+invwhiten(a::ScalMat, x::AbstractVecOrMat) = unwhiten(a, x)
 function unwhiten(a::ScalMat, x::AbstractVecOrMat)
     @check_argdims a.dim == size(x, 1)
     return sqrt(a.value) * x
 end
+invunwhiten(a::ScalMat, x::AbstractVecOrMat) = whiten(a, x)
 
 ### quadratic forms