Generalized CP Decomposition module. Provides approximate CP tensor decomposition with respect to general losses.

CPD

Tensor decomposition type for the canonical polyadic decompositions (CPD) of a tensor (i.e., a multi-dimensional array) A. This is the return type of gcp(_), the corresponding tensor decomposition function.

If M::CPD is the decomposition object, the weights λ and the factor matrices U = (U[1],...,U[N]) can be obtained via M.λ and M.U, such that A = Σ_j λ[j] U[1][:,j] ∘ ⋯ ∘ U[N][:,j].

default_algorithm(X, r, loss, constraints)

Return a default algorithm for the data tensor X, rank r, loss function loss, and tuple of constraints constraints.

See also: gcp.

default_constraints(loss)

Return a default tuple of constraints for the loss function loss.

See also: gcp.

default_init([rng=default_rng()], X, r, loss, constraints, algorithm)

Return a default initialization for the data tensor X, rank r, loss function loss, tuple of constraints constraints, and algorithm algorithm, using the random number generator rng if needed.

See also: gcp.

gcp(X::Array, r;
    loss = GCPLosses.LeastSquares(),
    constraints = default_constraints(loss),
    algorithm = default_algorithm(X, r, loss, constraints),
    init = default_init(X, r, loss, constraints, algorithm))

Compute an approximate rank-r CP decomposition of the tensor X with respect to the loss function loss and return a CPD object.

Keyword arguments:

• constraints : a Tuple of constraints on the factor matrices U = (U[1],...,U[N]).
• algorithm : algorithm to use

Conventional CP corresponds to the default GCPLosses.LeastSquares() loss with the default of no constraints (i.e., constraints = ()).

If the LossFunctions.jl package is also loaded, loss can also be a loss function from that package. Check GCPDecompositions.LossFunctionsExt.SupportedLosses to see what losses are supported.

See also: CPD, GCPLosses, GCPConstraints, GCPAlgorithms.

ncomps(M::CPD)

Return the number of components in M.

See also: ndims, size.

Algorithms for Generalized CP Decomposition.

ALS

Alternating Least Squares. Workhorse algorithm for LeastSquares loss with no constraints.

Algorithm parameters:

• maxiters::Int : max number of iterations (default: 200)

AbstractAlgorithm

Abstract type for GCP algorithms.

Concrete types ConcreteAlgorithm <: AbstractAlgorithm should implement _gcp(X, r, loss, constraints, algorithm::ConcreteAlgorithm) that returns a CPD.

FastALS

Fast Alternating Least Squares.

Efficient ALS algorithm proposed in:

Fast Alternating LS Algorithms for High Order CANDECOMP/PARAFAC Tensor Factorizations. Anh-Huy Phan, Petr Tichavský, Andrzej Cichocki. IEEE Transactions on Signal Processing, 2013. DOI: 10.1109/TSP.2013.2269903

Algorithm parameters:

• maxiters::Int : max number of iterations (default: 200)

LBFGSB

Limited-memory BFGS with Box constraints.

Brief description of algorithm parameters:

• m::Int : max number of variable metric corrections (default: 10)
• factr::Float64 : function tolerance in units of machine epsilon (default: 1e7)
• pgtol::Float64 : (projected) gradient tolerance (default: 1e-5)
• maxfun::Int : max number of function evaluations (default: 15000)
• maxiter::Int : max number of iterations (default: 15000)
• iprint::Int : verbosity (default: -1)
  • iprint < 0 means no output
  • iprint = 0 prints only one line at the last iteration
  • 0 < iprint < 99 prints f and |proj g| every iprint iterations
  • iprint = 99 prints details of every iteration except n-vectors
  • iprint = 100 also prints the changes of active set and final x
  • iprint > 100 prints details of every iteration including x and g

See documentation of LBFGSB.jl for more details.

FastALS_iter!(X, U, λ) 

Algorithm for computing MTTKRP sequences is from "Fast Alternating LS Algorithms
for High Order CANDECOMP/PARAFAC Tensor Factorizations" by Phan et al., specifically
section III-C.

_gcp(X, r, loss, constraints, algorithm)

Internal function to compute an approximate rank-r CP decomposition of the tensor X with respect to the loss function loss and the constraints constraints using the algorithm algorithm, returning a CPD object.

Loss functions for Generalized CP Decomposition.

AbstractLoss

Abstract type for GCP loss functions $f(x,m)$, where $x$ is the data entry and $m$ is the model entry.

Concrete types ConcreteLoss <: AbstractLoss should implement:

• value(loss::ConcreteLoss, x, m) that computes the value of the loss function $f(x,m)$
• deriv(loss::ConcreteLoss, x, m) that computes the value of the partial derivative $\partial_m f(x,m)$ with respect to $m$
• domain(loss::ConcreteLoss) that returns an Interval from IntervalSets.jl defining the domain for $m$

BernoulliLogit(eps::Real = 1e-10)

Loss corresponding to the statistical assumption of Bernouli data X with log odds-success rate given by the low-rank model tensor M

• Distribution: $x_i \sim \operatorname{Bernouli}(\rho_i)$
• Link function: $m_i = \log(\frac{\rho_i}{1 - \rho_i})$
• Loss function: $f(x, m) = \log(1 + e^m) - xm$
• Domain: $m \in \mathbb{R}$

BernoulliOdds(eps::Real = 1e-10)

Loss corresponding to the statistical assumption of Bernouli data X with odds-sucess rate given by the low-rank model tensor M

• Distribution: $x_i \sim \operatorname{Bernouli}(\rho_i)$
• Link function: $m_i = \frac{\rho_i}{1 - \rho_i}$
• Loss function: $f(x, m) = \log(m + 1) - x\log(m + \epsilon)$
• Domain: $m \in [0, \infty)$

BetaDivergence(β::Real, eps::Real)

BetaDivergence Loss for given β

• Loss function: $f(x, m; β) = \frac{1}{\beta}m^{\beta} - \frac{1}{\beta - 1}xm^{\beta - 1} if \beta \in \mathbb{R} \{0, 1\}, m - x\log(m) if \beta = 1, \frac{x}{m} + \log(m) if \beta = 0$
• Domain: $m \in [0, \infty)$

Gamma(eps::Real = 1e-10)

Loss corresponding to a statistical assumption of Gamma-distributed data X with scale given by the low-rank model tensor M.

• Distribution: $x_i \sim \operatorname{Gamma}(k, \sigma_i)$
• Link function: $m_i = k \sigma_i$
• Loss function: $f(x,m) = \frac{x}{m + \epsilon} + \log(m + \epsilon)$
• Domain: $m \in [0, \infty)$

Huber(Δ::Real)

Huber Loss for given Δ

• Loss function: $f(x, m) = (x - m)^2 if \abs(x - m)\leq\Delta, 2\Delta\abs(x - m) - \Delta^2 otherwise$
• Domain: $m \in \mathbb{R}$

LeastSquares()

Loss corresponding to conventional CP decomposition. Corresponds to a statistical assumption of Gaussian data X with mean given by the low-rank model tensor M.

• Distribution: $x_i \sim \mathcal{N}(\mu_i, \sigma)$
• Link function: $m_i = \mu_i$
• Loss function: $f(x,m) = (x-m)^2$
• Domain: $m \in \mathbb{R}$

NegativeBinomialOdds(r::Integer, eps::Real = 1e-10)

Loss corresponding to the statistical assumption of Negative Binomial data X with log odds failure rate given by the low-rank model tensor M

• Distribution: $x_i \sim \operatorname{NegativeBinomial}(r, \rho_i)$
• Link function: $m = \frac{\rho}{1 - \rho}$
• Loss function: $f(x, m) = (r + x) \log(1 + m) - x\log(m + \epsilon)$
• Domain: $m \in [0, \infty)$

NonnegativeLeastSquares()

Loss corresponding to nonnegative CP decomposition. Corresponds to a statistical assumption of Gaussian data X with nonnegative mean given by the low-rank model tensor M.

• Distribution: $x_i \sim \mathcal{N}(\mu_i, \sigma)$
• Link function: $m_i = \mu_i$
• Loss function: $f(x,m) = (x-m)^2$
• Domain: $m \in [0, \infty)$

Poisson(eps::Real = 1e-10)

Loss corresponding to a statistical assumption of Poisson data X with rate given by the low-rank model tensor M.

• Distribution: $x_i \sim \operatorname{Poisson}(\lambda_i)$
• Link function: $m_i = \lambda_i$
• Loss function: $f(x,m) = m - x \log(m + \epsilon)$
• Domain: $m \in [0, \infty)$

PoissonLog()

Loss corresponding to a statistical assumption of Poisson data X with log-rate given by the low-rank model tensor M.

• Distribution: $x_i \sim \operatorname{Poisson}(\lambda_i)$
• Link function: $m_i = \log \lambda_i$
• Loss function: $f(x,m) = e^m - x m$
• Domain: $m \in \mathbb{R}$

Rayleigh(eps::Real = 1e-10)

Loss corresponding to the statistical assumption of Rayleigh data X with sacle given by the low-rank model tensor M

• Distribution: $x_i \sim \operatorname{Rayleigh}(\theta_i)$
• Link function: $m_i = \sqrt{\frac{\pi}{2}\theta_i}$
• Loss function: $f(x, m) = 2\log(m + \epsilon) + \frac{\pi}{4}(\frac{x}{m + \epsilon})^2$
• Domain: $m \in [0, \infty)$

UserDefined

Type for user-defined loss functions $f(x,m)$, where $x$ is the data entry and $m$ is the model entry.

Contains three fields:

1. func::Function : function that evaluates the loss function $f(x,m)$
2. deriv::Function : function that evaluates the partial derivative $\partial_m f(x,m)$ with respect to $m$
3. domain::Interval : Interval from IntervalSets.jl defining the domain for $m$

The constructor is UserDefined(func; deriv, domain). If not provided,

• deriv is automatically computed from func using forward-mode automatic differentiation
• domain gets a default value of Interval(-Inf, +Inf)

deriv(loss, x, m)

Compute the derivative of the (entrywise) loss function loss at the model entry m for the data entry x.

domain(loss)

Return the domain of the (entrywise) loss function loss.

grad_U!(GU, M::CPD, X::AbstractArray, loss)

Compute the GCP gradient with respect to the factor matrices U = (U[1],...,U[N]) for the model tensor M, data tensor X, and loss function loss, and store the result in GU = (GU[1],...,GU[N]).

objective(M::CPD, X::AbstractArray, loss)

Compute the GCP objective function for the model tensor M, data tensor X, and loss function loss.

value(loss, x, m)

Compute the value of the (entrywise) loss function loss for data entry x and model entry m.

Constraints for Generalized CP Decomposition.

AbstractConstraint

Abstract type for GCP constraints on the factor matrices U = (U[1],...,U[N]).

LowerBound(value::Real)

Lower-bound constraint on the entries of the factor matrices U = (U[1],...,U[N]), i.e., U[i][j,k] >= value.

Tensor kernels for Generalized CP Decomposition.

_checked_khatrirao_dims(A1, A2, ...)

Check that A1, A2, etc. have compatible dimensions for the Khatri-Rao product. If so, return a tuple of the number of rows and the shared number of columns. If not, throw an error.

_checked_mttkrp_dims(X, (U1, U2, ..., UN), n)

Check that X and U have compatible dimensions for the mode-n MTTKRP. If so, return a tuple of the number of rows and the shared number of columns for the Khatri-Rao product. If not, throw an error.

_checked_mttkrps_dims(X, (U1, U2, ..., UN))

Check that X and U have compatible dimensions for the mode-n MTTKRP. If so, return a tuple of the number of rows and the shared number of columns for the Khatri-Rao product. If not, throw an error.

create_mttkrp_buffer(X, U, n)

Create buffer to hold intermediate calculations in mttkrp!.

Always use create_mttkrp_buffer to make a buffer for mttkrp!; the internal details of buffer may change in the future and should not be relied upon.

See also: mttkrp!

khatrirao!(K, A1, A2, ...)

Compute the Khatri-Rao product (i.e., the column-wise Kronecker product) of the matrices A1, A2, etc. and store the result in K.

khatrirao(A1, A2, ...)

Compute the Khatri-Rao product (i.e., the column-wise Kronecker product) of the matrices A1, A2, etc.

mttkrp!(G, X, (U1, U2, ..., UN), n, buffer=create_mttkrp_buffer(X, U, n))

Compute the Matricized Tensor Times Khatri-Rao Product (MTTKRP) of an N-way tensor X with the matrices U1, U2, ..., UN along mode n and store the result in G.

Optionally, provide a buffer for intermediate calculations. Always use create_mttkrp_buffer to make the buffer; the internal details of buffer may change in the future and should not be relied upon.

Algorithm is based on Section III-B of the paper:

Fast Alternating LS Algorithms for High Order CANDECOMP/PARAFAC Tensor Factorizations. Anh-Huy Phan, Petr Tichavský, Andrzej Cichocki. IEEE Transactions on Signal Processing, 2013. DOI: 10.1109/TSP.2013.2269903

See also: mttkrp, create_mttkrp_buffer

mttkrp(X, (U1, U2, ..., UN), n)

Compute the Matricized Tensor Times Khatri-Rao Product (MTTKRP) of an N-way tensor X with the matrices U1, U2, ..., UN along mode n.

See also: mttkrp!

mttkrps!(G, X, (U1, U2, ..., UN))

Compute the Matricized Tensor Times Khatri-Rao Product Sequence (MTTKRPS) of an N-way tensor X with the matrices U1, U2, ..., UN and store the result in G.

See also: mttkrps

mttkrps(X, (U1, U2, ..., UN))

Compute the Matricized Tensor Times Khatri-Rao Product Sequence (MTTKRPS) of an N-way tensor X with the matrices U1, U2, ..., UN.

See also: mttkrps!