PDMats.jl
Uniform interface for positive definite matrices of various structures.
Positive definite matrices are widely used in machine learning and probabilistic modeling, especially in applications related to graph analysis and Gaussian models. It is not uncommon that positive definite matrices used in practice have special structures (e.g. diagonal), which can be exploited to accelerate computation.
PDMats.jl supports efficient computation on positive definite matrices of various structures. In particular, it provides uniform interfaces to use positive definite matrices of various structures for writing generic algorithms, while ensuring that the most efficient implementation is used in actual computation.
Positive definite matrix types
This package defines an abstract type AbstractPDMat{T<:Real} as the base type for positive definite matrices with different internal representations.
Elemenent types are in princple all Real types, but in practice this is limited by the support for floating point types in Base.LinAlg.Cholesky.
Float64Fully supported from Julia 0.3.Float32Fully supported from Julia 0.4.2. Full, diagonal and scale matrix types are supported in Julia 0.3 or higher.Float16Promoted toFloat32for full, diagonal and scale matrix. Currently unsupported for sparse matrix.BigFloatSupported in Julia 0.4 for full, diagonal and scale matrix. Currently unsupported for sparse matrix.
PDMat: full covariance matrix, defined as
struct PDMat{T<:Real,S<:AbstractMatrix} <: AbstractPDMat{T}
dim::Int # matrix dimension
mat::S # input matrix
chol::Cholesky{T,S} # Cholesky factorization of mat
end
# Constructors
PDMat(mat, chol) # with both the input matrix and a Cholesky factorization
PDMat(mat) # with the input matrix, of type Matrix or Symmetric
# Remarks: the Cholesky factorization will be computed
# upon construction.
PDMat(chol) # with the Cholesky factorization
# Remarks: the full matrix will be computed upon
# construction.PDiagMat: diagonal matrix, defined as
struct PDiagMat{T<:Real,V<:AbstractVector{T}} <: AbstractPDMat{T}
dim::Int # matrix dimension
diag::V # the vector of diagonal elements
end
# Constructors
PDiagMat(v) # with the vector of diagonal elementsScalMat: uniform scaling matrix, asv * eye(d), defined as
struct ScalMat{T<:Real} <: AbstractPDMat{T}
dim::Int # matrix dimension
value::T # diagonal value (shared by all diagonal elements)
end
# Constructors
ScalMat(d, v) # with dimension d and diagonal value vPDSparseMat: sparse covariance matrix, defined as
struct PDSparseMat{T<:Real,S<:AbstractSparseMatrix} <: AbstractPDMat{T}
dim::Int # matrix dimension
mat::SparseMatrixCSC # input matrix
chol::CholTypeSparse # Cholesky factorization of mat
end
# Constructors
PDSparseMat(mat, chol) # with both the input matrix and a Cholesky factorization
PDSparseMat(mat) # with the sparse input matrix, of type SparseMatrixCSC
# Remarks: the Cholesky factorization will be computed
# upon construction.
PDSparseMat(chol) # with the Cholesky factorization
# Remarks: the sparse matrix 'mat' will be computed upon
# construction.Common interface
All subtypes of AbstractPDMat share the same API, i.e. with the same set of methods to operate on their instances. These methods are introduced below, where a is an instance of a subtype of AbstractPDMat to represent a positive definite matrix, x be a column vector or a matrix with size(x,1) == size(a, 1) == size(a, 2), and c be a positive real value.
size(a) # return the size of `a`.
size(a, i) # return the i-th dimension of `a`.
ndims(a) # the number of dimensions, which is always 2.
eltype(a) # the element type
Matrix(a) # return a copy of the matrix in full form.
diag(a) # return a vector of diagonal elements.
inv(a) # inverse of `a`, of a proper subtype of `AbstractPDMat`.
# Note: when `a` is an instance of either `PDMat`, `PDiagMat`,
# and `ScalMat`, `inv(a)` is of the same type of `a`.
# This needs not be required for customized subtypes -- the
# inverse does not always has the same pattern as `a`.
eigmax(a) # maximum eigenvalue of `a`.
eigmin(a) # minimum eigenvalue of `a`.
logdet(a) # log-determinant of `a`, computed in a numerically stable way.
a * x # multiple `a` with `x` (forward transform)
a \ x # multiply `inv(a)` with `x` (backward transform).
# The internal implementation may not explicitly instantiate
# the inverse of `a`.
a * c # scale `a` by a positive scale `c`.
# The result is in general of the same type of `a`.
c * a # equivalent to a * c.
a + b # add two positive definite matrices
pdadd(a, b, c) # add `a` with `b * c`, where both `a` and `b` are
# instances of `AbstractPDMat`.
pdadd(m, a) # add `a` to a dense matrix `m` of the same size.
pdadd(m, a, c) # add `a * c` to a dense matrix `m` of the same size.
pdadd!(m, a) # add `a` to a dense matrix `m` of the same size inplace.
pdadd!(m, a, c) # add `a * c` to a dense matrix `m` of the same size,
# inplace.
pdadd!(r, m, a) # add `a` to a dense matrix `m` of the same size, and write
# the result to `r`.
pdadd!(r, m, a, c) # add `a * c` to a dense matrix `m` of the same size, and
# write the result to `r`.
quad(a, x) # compute x' * a * x when `x` is a vector.
# perform such computation in a column-wise fashion, when
# `x` is a matrix, and return a vector of length `n`,
# where `n` is the number of columns in `x`.
quad!(r, a, x) # compute x' * a * x in a column-wise fashion, and write
# the results to `r`.
invquad(a, x) # compute x' * inv(a) * x when `x` is a vector.
# perform such computation in a column-wise fashion, when
# `x` is a matrix, and return a vector of length `n`.
invquad!(r, a, x) # compute x' * inv(a) * x in a column-wise fashion, and
# write the results to `r`.
X_A_Xt(a, x) # compute `x * a * x'` for a matrix `x`.
Xt_A_X(a, x) # compute `x' * a * x` for a matrix `x`.
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.
#
# Note: If the covariance of `x` is `a`, then the
# covariance of the transformed result is an identity
# matrix.
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.
#
# Note: If the covariance of `x` is an identity matrix,
# then the covariance of the transformed result is `a`.
# Note: the un-whitening transform is useful for
# generating Gaussian samples.
unwhiten!(a, x) # un-whitening transform inplace, updating `x`.
unwhiten!(r, a, x) # write the transformed result to `r`.Fallbacks for AbstractArrays
For ease of composability, some of these functions have generic fallbacks defined that work on AbstractArrays.
These fallbacks may not be as fast as the methods specializaed for AbstractPDMats, but they let you more easily swap out types.
While in theory all of them can be defined, at present only the following subset has:
whiten,whiten!unwhiten,unwhiten!quad,quad!invquad,invquad!
PRs to implement more generic fallbacks are welcome.
Define Customized Subtypes
In some situation, it is useful to define a customized subtype of AbstractPDMat to capture positive definite matrices with special structures. For this purpose, one has to define a subset of methods (as listed below), and other methods will be automatically provided.
# Let `M` be the name of the subtype, then the following methods need
# to be implemented for `M`:
Matrix(a::M) # return a copy of the matrix in full form, of type
# `Matrix{eltype(M)}`.
diag(a::M) # return the vector of diagonal elements, of type
# `Vector{eltype(M)}`.
pdadd!(m, a, c) # add `a * c` to a dense matrix `m` of the same size
# inplace.
* (a::M, c::Real) # return a scaled version of `a`.
* (a::M, x::DenseVecOrMat) # transform `x`, i.e. compute `a * x`.
\ (a::M, x::DenseVecOrMat) # inverse transform `x`, i.e. compute `inv(a) * x`.
inv(a::M) # compute the inverse of `a`.
logdet(a::M) # compute the log-determinant of `a`.
eigmax(a::M) # compute the maximum eigenvalue of `a`.
eigmin(a::M) # compute the minimum eigenvalue of `a`.
whiten!(r::DenseVecOrMat, a::M, x::DenseVecOrMat) # whitening transform,
# write result to `r`.
unwhiten!(r::DenseVecOrMat, a::M, x::DenseVecOrMat) # un-whitening transform,
# write result to `r`.
quad(a::M, x::DenseVector) # compute `x' * a * x`
quad!(r::AbstractArray, a::M, x::DenseMatrix) # compute `x' * a * x` in
# a column-wise manner
invquad(a::M, x::DenseVector) # compute `x' * inv(a) * x`
invquad!(r::AbstractArray, a::M, x::DenseMatrix) # compute `x' * inv(a) * x`
# in a column-wise manner
X_A_Xt(a::M, x::DenseMatrix) # compute `x * a * x'`
Xt_A_X(a::M, x::DenseMatrix) # compute `x' * a * x`
X_invA_Xt(a::M, x::DenseMatrix) # compute `x * inv(a) * x'`
Xt_invA_X(a::M, x::DenseMatrix) # compute `x' * inv(a) * x`